独立なデータ分布の対数正規分布によるベイズ分析
ベイズ分析の基本的考え方の説明を岡本（2019）、pp.169-177、岡本（2014）、pp.227-232で行っている。
確率変数
の対数
が平均
、分散
の正規分布に従うとき、確率変数は対数正規分布に従うと言い、確率密度関数は次式（１）で与えられる（Gelman et al., 2014）。
このとき、平均値、分散、モードは次式で与えられる。
対数正規分布のパラメータをと
ではなく、対数正規分布の平均値
と標準偏差
をパラメータとすると、パラメータと分布との関係が分かり易くなる。
式（２）、式（３）より
とおけば、
であり、
である。
いま、２つの独立なデータ
と
がそれぞれ対数正規分布の母集団
と
から与えられたとき、事後分布を式（７）のようにおく。
事前分布
を一様分布としたときのStanスクリプトをリスト１に示す。
リスト１のStanスクリプトを用いてデータ分析を行うスクリプト例を、リスト２とリスト３に示す。
これらのスクリプトファイルなどは圧縮ファイルbalognormalfiles.zipとしてまとめた。
次節において、これらのスクリプトによる分析例を示す。
分析例
Windows上でも実行できるように、メインのスクリプトファイルをリスト２とリスト３に分割した。リスト２のスクリプトにおいてリスト１のStanスクリプトによるMCMCサンプリングを行い、MCMCサンプリングの結果をファイルに出力している。そのファイルをリスト３のスクリプトで読み込んで、グラフ描画を行っている。リスト２のスクリプトはWindowsでは直接実行できないので、WindowsにWSL（Windows Subsystem for Linux）をインストールして、WSL上にインストールされたUbuntu上でリスト２のスクリプトを実行した。Ubuntu上での実行は、miniconda上に用意したPython3.12仮想環境で行った。すなわち、リスト２のスクリプトは、Python3.12/miniconda/Ubuntu/WSL/Windows 10で実行した。Minicondaの使用は、Anacondaと同様である。
リスト２のスクリプトを実行すると、入力データファイル名の設定が求められる。
(py312)
yasuharu@PSY-PM-PC:/mnt/d/*****/balognormalfiles$ python ba_means.py
Input data
file(*.csv) =
入力データファイルは、CSV形式のファイルとして、図１のように用意する。
・
・
・
1行目に変数名を設定する。１列目は空欄でもよい。データは２行目から設定する。１列目はデータの識別用である。２列目と３列目にデータ値を設定する。データ数が２つのグループで異なるときは、少ない方のデータは空欄にしておけばよい。空欄はスクリプトの実行時に無視される。
作成したデータファイルは、CSVファイルとして保存する。データファイルとリスト１からリスト３までのファイルは、同じフォルダ内に置いておく。
次の例では、入力データファイルとして図１のファイルsampleN47N50.csvを設定している。
(py312)
yasuharu@PSY-PM-PC:/mnt/d/*****/balognormalfiles$ python ba_means.py
Input data
file(*.csv) = sampleN47N50.csv
入力データファイル名を設定すると、計算が始まり、終了すると以下のようなメッセージが表示されてリスト２の実行が終了する。
Sampling: 100% (8000/8000), done.
Messages received during sampling:
Gradient evaluation took 2.2e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.22
seconds.
Adjust your expectations accordingly!
Gradient evaluation took 2.2e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.22
seconds.
Adjust your expectations accordingly!
Gradient evaluation took 2e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.2
seconds.
Adjust your expectations accordingly!
Gradient evaluation took 2.2e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.22
seconds.
Adjust your expectations accordingly!
(py312)
yasuharu@PSY-PM-PC:/mnt/d/*****/ba_lognormal/balognormalfiles$
リスト２の実行において読み込まれたデータとMCMCサンプルはファイルに保存されている。
保存されたファイルは、リスト３のスクリプトで読み込まれ、描画などの処理が行われる。Ubuntu/WSL/Windowsではグラフ描画ができなかった（2025.01.23）ので、リスト３はWindows上で直接実行した。
リスト２の実行後、リスト３を実行すると図２のグラフが表示される。
図２
パラメータの事後分布が描かれている。
図２のWindowを閉じると、図３のグラフが表示される。
図３
パラメータの差の事後分布が表示される。タイトル欄に確率
の値が表示されている。
の事後確率が0.5％以下であることが分かる。
および
は、データAおよびデータBに対するパラメータであることを表す（図１）。
であるので、
である。連続な確率変数が特定の値を取る確率は０である。
図３のWindowを閉じると、図４のグラフが表示される。
図４
パラメータの事後分布が描かれている。
図４のWindowを閉じると、図５のグラフが表示される。
図５
の差の事後分布が描かれている。
の事後確率は０と見なされる値であることがタイトル欄に表示されている。
図５のWindowを閉じると、図６のグラフが表示される。
図６
パラメータの事後分布のグラフである。
図６のWindowを閉じると、図７のグラフが表示される。
図７
パラメータの差の事後分布のグラフである。
の事後確率が約10％であることがタイトル欄に表示されている。パラメータの差に比べると、の差は不確かである（図３）。
図７のWindowを閉じると、図８のグラフが表示される。
図８
パラメータの事後分布である。
図８のWindowを閉じると、図９のグラフが表示される。
図９
パラメータの差の事後分布である。
の事後確率が０と見なせることがタイトル欄に表示されている。
図９のWindowを閉じると、図10のグラフが表示される。
図10
曲線は、タイトル欄に示されているパラメータ値で描かれた対数正規分布のグラフである。パラメータ値は事後分布の中央値を採用している（Gelman et al.,2021）。
図10のWindowを閉じると、図11のグラフが表示される。
図11
モード（4）の事後分布のグラフが表示されている。
図11のWindowを閉じると、図12のグラフが表示される。
図12
モードの差の事後分布のグラフである。
の事後確率が約90％であることがタイトル欄に表示されている。モードはAの方が右側にあると言えそうである。
平均値はBの方が大きかった（図10）。モードと平均値の大小関係の逆転は、標準偏差（分散）の違いによる（図４）。対数正規分布は左右対称な分布ではないので、平均だけでは分布の相対的関係を表すことができない。
Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., & Rubin, D. B. (2014). Bayesian data analysis, 3rd. Ed. CRC Press.
Gelman, A., Hill, J., & Vehtari, A. (2021). Regression and other stories. Cambridge University
岡本安晴（2019）「いまさら聞けないPythonでデータ分析」丸善出版
リスト１ 対数正規分布のStanスクリプト例（lognormal_1.stan）
data {
int n1;
array[n1] real x1;
int n2;
array[n2] real x2;
real max_v;
}
parameters {
real<lower=0.0,
upper=1.0> vmn1;
real<lower=0.0,
upper=1.0> vmn2;
real<lower=0.0,
upper=1.0> vsd1;
real<lower=0.0,
upper=1.0> vsd2;
}
transformed parameters {
real mn1;
real mn2;
real sd1;
real sd2;
real mu1;
real mu2;
real sgm1;
real sgm2;
mn1 = 0.0001 + vmn1 * max_v;
mn2 = 0.0001 + vmn2 * max_v;
sd1 = 0.0001 + vsd1 * max_v;
sd2 = 0.0001 + vsd2 * max_v;
sgm1 = (log((sd1/mn1)^2 +
1))^0.5;
mu1 = log(mn1) - 0.5 *
(sgm1^2);
sgm2 = log((sd2/mn2)^2 +
1)^0.5;
mu2 = log(mn2) - 0.5 *
(sgm2^2);
}
model {
vmn1 ~ beta(1, 1);
vmn2 ~ beta(1, 1);
vsd1 ~ beta(1, 1);
vsd2 ~ beta(1, 1);
x1 ~ lognormal(mu1, sgm1);
x2 ~ lognormal(mu2, sgm2);
}
generated quantities {
real diff_mns;
real diff_sds;
real diff_mus;
real diff_sgms;
real mode1;
real mode2;
real diff_modes;
diff_mns = mn1 - mn2;
diff_sds = sd1 - sd2;
diff_mus = mu1 - mu2;
diff_sgms = sgm1 - sgm2;
mode1 = exp(mu1 - sgm1^2);
mode2 = exp(mu2 - sgm2^2);
diff_modes = mode1 - mode2;
}
リスト２ リスト１のStanスクリプトを用いてMCMCサンプリングを行うスクリプト（ba_means.py）
import numpy as np
import pandas as pd
import stan
import pickle
inf_nm = input('Input data file(*.csv)
= ')
dframe = pd.read_csv(inf_nm)
with open('in_dframe.pkl', 'wb') as f:
pickle.dump(dframe, f)
# Save the input data
df_keys = dframe.keys()
x_name = df_keys[1]
y_name = df_keys[2]
print(x_name, y_name)
x = dframe[x_name].values
x = x[np.isnan(x)==False]
print('X =\n', x)
y = dframe[y_name].values
y = y[np.isnan(y)==False]
print('Y =\n', y)
mean1 = np.mean(x)
sd1 = np.std(x)
print('mean1 =', mean1, ' sd1 =', sd1)
mean2 = np.mean(y)
sd2 = np.std(y)
print('mean2 =', mean2, ' sd2 =', sd2)
max_v = np.max(x) if np.max(x) >
np.max(y) else np.max(y);
with open('lognormal_1.stan', 'r') as
f_stan:
sm = stan.build(f_stan.read(),
data = {'n1':len(x), 'x1':x,
'n2':len(y), 'x2':y,
'max_v':max_v})
fit = sm.sample()
fit_dframe = fit.to_frame()
with open('fit_dframe.pkl', 'wb') as
f:
pickle.dump(fit_dframe,
f)
# Save MCMC sample
リスト３ リスト２によるMCMCサンプルを分析するスクリプト（draw_ba.py）
import numpy as np
import pickle
import matplotlib.pyplot as plt
import seaborn as sb
def lognormal_pdf(x, mu, sgm):
v1 = sgm * x *
((2*np.pi)**0.5)
v2 = np.exp(-((np.log(x) - mu)**2) / (2 *
(sgm**2)))
return v2 / v1
with open('in_dframe.pkl', 'rb') as f:
in_dframe =
pickle.load(f)
# Read out the input
data
df_keys = in_dframe.keys()
x_name = df_keys[1]
y_name = df_keys[2]
print(x_name, y_name)
x = in_dframe[x_name].values
x = x[np.isnan(x)==False]
print('X =\n', x)
y = in_dframe[y_name].values
y = y[np.isnan(y)==False]
print('Y =\n', y)
with open('fit_dframe.pkl', 'rb') as
f:
fit_dframe =
pickle.load(f)
# Read out MCMC sample
print(fit_dframe.keys())
plt_mns = {x_name:fit_dframe['mn1'],
y_name:fit_dframe['mn2']}
sb.kdeplot(data=plt_mns)
plt.xlabel('$mn$', fontsize=14)
plt.title('Posterior distrbutions of
$mn$s')
plt.show()
diff_means = fit_dframe['diff_mns']
p_pos = np.mean(diff_means > 0)
sb.kdeplot(diff_means)
plt.xlabel(f'$mn$[{x_name}] -
$mn$[{y_name}]', fontsize=14)
plt.title(f'P($mn$[{x_name}] -
$mn$[{y_name}] > 0) = {p_pos}' + '\n' +
fr'Posterior distribution of $mn$[{x_name}] - $mn$[{y_name}]')
plt.show()
plt_sds = {x_name:fit_dframe['sd1'],
y_name:fit_dframe['sd2']}
sb.kdeplot(data=plt_sds)
plt.xlabel('$sd$', fontsize=14)
plt.title('Posterior distrbutions of
$sd$s')
plt.show()
diff_sds = fit_dframe['diff_sds']
p_pos = np.mean(diff_sds > 0)
sb.kdeplot(diff_sds)
plt.xlabel(fr'$sd$[{x_name}] -
$sd$[{y_name}]', fontsize=14)
plt.title(fr'P($sd$[{x_name}] -
$sd$[{y_name}] > 0) = {p_pos}' + '\n' +
fr'Posterior distribution of $sd$[{x_name}] - $sd$[{y_name}]')
plt.show()
plt_mus = {x_name:fit_dframe['mu1'],
y_name:fit_dframe['mu2']}
sb.kdeplot(data=plt_mus)
plt.xlabel(r'$\mu$', fontsize=14)
plt.title(r'Posterior distrbutions of
$\mu$s')
plt.show()
diff_mus = fit_dframe['diff_mus']
p_pos = np.mean(diff_mus > 0)
sb.kdeplot(diff_mus)
plt.xlabel(fr'$\mu$[{x_name}] -
$\mu$[{y_name}]', fontsize=14)
plt.title(fr'P($\mu$[{x_name}] -
$\mu$[{y_name}] > 0) = {p_pos}' + '\n' +
fr'Posterior distribution of $\mu$[{x_name}] - $\mu$[{y_name}]')
plt.show()
plt_sgms = {x_name:fit_dframe['sgm1'],
y_name:fit_dframe['sgm2']}
sb.kdeplot(data=plt_mns)
plt.xlabel(r'$\sigma$', fontsize=14)
plt.title(r'Posterior distrbutions of
$\sigma$s')
plt.show()
diff_sgms = fit_dframe['diff_sgms']
p_pos = np.mean(diff_sgms > 0)
sb.kdeplot(diff_sgms)
plt.xlabel(fr'$\sigma$[{x_name}] -
$\sigma$[{y_name}]', fontsize=14)
plt.title(fr'P($\sigma$[{x_name}] -
$\sigma$[{y_name}] > 0) = {p_pos}' + '\n' +
fr'Posterior distribution of $\sigma$[{x_name}] - $\sigma$[{y_name}]')
plt.show()
mu1 = np.median(fit_dframe['mu1'])
sgm1 = np.median(fit_dframe['sgm1'])
max_x = np.max(x) if np.max(x) >
np.max(y) else np.max(y)
min_x = np.min(x) if np.min(x) <
np.min(y) else np.min(y)
rng_x = np.linspace(min_x, max_x,
1000)
d1 = lognormal_pdf(rng_x, mu1, sgm1)
plt.plot(rng_x, d1, label=x_name)
mu2 = np.median(fit_dframe['mu2'])
sgm2 = np.median(fit_dframe['sgm2'])
d2 = lognormal_pdf(rng_x, mu2, sgm2)
plt.plot(rng_x, d2, label=y_name)
mn1 = np.median(fit_dframe['mn1'])
sd1 = np.median(fit_dframe['sd1'])
mn2 = np.median(fit_dframe['mn2'])
sd2 = np.median(fit_dframe['sd2'])
mode1 = np.median(fit_dframe['mode1'])
mode2 = np.median(fit_dframe['mode2'])
plt.hist([x, y], density=True,
label=[x_name, y_name])
plt.title('Histograms of data and
Probability densities\n' +
fr'{x_name}: $\mu={mu1:.1f}, \sigma={sgm1:.1f}$' +
fr', $mn$={mn1:.1f}, $sd$={sd1:.1f}, $mode$={mode1:.1f}' + '\n' +
fr' {y_name}:
$\mu={mu2:.1f}, \sigma={sgm2:.1f}$' +
fr', $mn$={mn2:.1f}, $sd$={sd2:.1f}, $mode$={mode2:.1f}')
plt.ylabel('Density')
plt.legend()
plt.tight_layout()
plt.show()
mode_samples =
{x_name:fit_dframe['mode1'], y_name:fit_dframe['mode2']}
sb.kdeplot(data=mode_samples)
plt.xlabel('mode', fontsize=14)
plt.title('Posterior distribution of
modes')
plt.show()
diff_modes = fit_dframe['diff_modes']
p_pos = np.mean(diff_modes > 0)
sb.kdeplot(diff_modes)
plt.xlabel(fr'$mode$[{x_name}] -
$mode$[{y_name}]', fontsize=14)
plt.title(fr'P($mode$[{x_name}] -
$mode$[{y_name}] > 0) = {p_pos}' + '\n' +
fr'Posterior distribution of $mode$[{x_name}] - $mode$[{y_name}]')
plt.show()