３パラメータ対数正規分布による独立な２つのデータの平均値のベイズ分析
確率変数
の対数
が平均
、分散
の正規分布に従うとき、確率変数は３パラメータの対数正規分布に従うという（Johnson et al., 1994）。確率密度関数は次式（１）のように書くことができる。
このとき、平均値、分散、モードは次式で与えられる。
対数正規分布のパラメータをと
ではなく、対数正規分布の平均値
と標準偏差
をパラメータとすると、パラメータと分布との関係が分かり易くなる。
式（２）、式（３）より
とおけば、
であり、
である。
いま、２つの独立なデータ
と
がそれぞれ対数正規分布
と
の母集団から与えられたとき、事後分布を式（７）のようにおく。
事前分布
を一様分布としたときのStanによるMCMCサンプリングのスクリプトをリスト１に示す。
リスト１によるMCMCサンプルを用いてデータ分析を行うスクリプト例を、リスト２に示す。
これらのスクリプトファイルなどは圧縮ファイルbalognormal3pfiles.zipとしてまとめた。
次節において、これらのスクリプトによる分析例を示す。
分析例
Windows上でも実行できるように、メインのスクリプトファイルをリスト１とリスト２に分割した。リスト１のスクリプトにおいてStanスクリプトによるMCMCサンプリングを行い、MCMCサンプリングの結果をファイルに出力している。そのファイルをリスト２のスクリプトで読み込んで、グラフ描画を行っている。リスト１のスクリプトはWindowsでは直接実行できないので、WindowsにWSL（Windows Subsystem for Linux）をインストールして、WSL上にインストールされたUbuntu上でリスト１のスクリプトを実行した。Ubuntu上での実行は、miniconda上に用意したPython3.12仮想環境で行った。すなわち、リスト１のスクリプトは、Python3.12/miniconda/Ubuntu22/WSL/Windows 11で実行した。リスト２は、Python3.13/Windows 11で実行した。
リスト１のスクリプトを実行すると、入力データファイル名の設定が求められる。
(py312)
*****/balognormal3pfiles$ python ba_means_lognormal3p.py
Input data
file(*.csv) =
入力データファイルは、CSV形式のファイルとして、図１のように用意する。
・
・
・
1行目に変数名を設定する。１列目は空欄でもよい。データは２行目から設定する。１列目はデータの識別用である。２列目と３列目にデータ値を設定する。データ数が２つのグループで異なるときは、少ない方のデータは空欄にしておけばよい。空欄はスクリプトの実行時に無視される。
作成したデータファイルは、CSVファイルとして保存する。データファイルとリスト１およびリスト２のファイルは、同じフォルダ内に置いておく。
次の例では、入力データファイルとして図１のファイルsampledata.csvを設定している。ファイルampledata.csvは圧縮ファイルbalognormal3pfiles.zipに含まれている。
(py312)
*****/balognormal3pfiles$ python ba_means_lognormal3p.py
Input data
file(*.csv) = sampledata.csv
入力データファイル名を設定すると、計算が始まり、終了すると以下のようなメッセージが表示されてリスト１の実行が終了する。
theta1
93.358 18.027 59.080 119.618 ... 0.550 1221.0 1194.0 1.0
theta2
83.577 39.052 12.247 142.386 ... 0.512 1173.0 1672.0 1.0
mean1
197.778 5.214 188.082 207.643 ... 0.086 2443.0 2124.0 1.0
mean2
170.033 1.070 168.139 172.175 ... 0.016 3428.0 2983.0 1.0
diff_means 27.745 5.325 18.277 38.206 ... 0.087 2462.0 2251.0 1.0
diff_sds 41.147 6.871 29.460 53.817 ... 0.178 1609.0 1530.0 1.0
diff_mus 0.179 0.531 -0.721 1.135 ... 0.007 1184.0 1624.0 1.0
diff_sgms 0.325 0.126 0.082 0.547 ... 0.002 1148.0 1348.0 1.0
mode1
168.462 6.015 157.702 180.008 ... 0.095 1764.0 1982.0 1.0
mode2
167.657 1.543 164.459 170.226 ... 0.024 1903.0 2243.0 1.0
diff_modes 0.805 6.184 -10.194 12.457 ... 0.099 1758.0 1969.0 1.0
[25 rows x 9 columns]
(py312) *****/balognormal3pfiles$
リスト１の実行において読み込まれたデータとMCMCサンプルはファイルに保存されている。
保存されたファイルは、リスト２のスクリプトで読み込まれ、描画などの処理が行われる。リスト２はWindows上で直接実行した。
リスト１の実行後、リスト２を実行すると図２のグラフが表示される。
図２
パラメータの事後分布の中央値を推定値とする対数正規分布のグラフが描かれている。２つの分布のモード(mode)はほぼ等しいが、平均値は異なっている。
図２のWindowを閉じると、図３のグラフが表示される。
図３
パラメータ
の事後分布である。
図３のWindowを閉じると、図４のグラフが表示される。
図４
平均値の事後分布である。２つの事後分布は、明瞭に分離されている。
図４のWindowを閉じると、図５のグラフが表示される。
図５
平均値の差
の事後分布と、事後確率
が示されている。平均値は異なると考えられる。
図５のWindowを閉じると、図６のグラフが表示される。
図６
パラメータの事後分布である。
図６のWindowを閉じると、図７のグラフが表示される。
図７
パラメータの差
の事後分布と確率
が表示されている。
図７のWindowを閉じると、図８のグラフが表示される。
図８
パラメータの事後分布である。
図８のWindowを閉じると、図９のグラフが表示される。
図９
パラメータの差
の事後分布と事後確率
が表示されている。２つのパラメータには、明瞭な差は認められない。
図９のWindowが閉じられると、図10のグラフが表示される。
図10
パラメータの差の事後分布である。明瞭に分離されている。
図10のWindowを閉じると、図11のグラフが表示される。
図11
パラメータの差
の事後分布と、事後確率
が表示されている。２つのパラメータの差は、ほぼ確実にあるといえそうである。
図11のWindowを閉じると、図12のグラフが表示される。
図12
最頻値（mode）の事後分布である。
図12のWindowを閉じると、図13のグラフが表示される。
図13
モードの差
の事後分布と事後確率
が表示されている。最頻値に差はないと考えてよさそうである。
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
Johnson, N. L., Kotz, S., & Barakrishnan, N. (1994) Continuous univariate distributions, Vol.1, 2nd Ed. John Wiley & Sons, Inc.
リスト１ MCMCサンプリングスクリプト（ba_means_lognormal3p.py）
"""
Yasuharu Okamoto, 2025.07
The script can
be used freely
"""
import numpy
as np
import pandas as pd
import stan
import arviz
as az
import pickle
stan_code = """
data {
int
n1;
array[n1]
real x1;
int
n2;
array[n2] real x2;
real
max_v;
// パラメータの変域設定用
real
min_v1;
// パラメータの変域設定用
real
min_v2;
// パラメータの変域設定用
}
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;
real<lower=0.0, upper=1.0> vtheta1;
real<lower=0.0, upper=1.0> vtheta2;
}
transformed parameters {
real
mn1;
real
mn2;
real
sd1;
real
sd2;
real
mu1;
real
mu2;
real
sgm1;
real
sgm2;
real
theta1;
real
theta2;
/*
パラメータのスケーリング
model宣言部で変域の広いパラメータに
出される警告を、スケーリングで避ける
*/
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);
theta1 = vtheta1 * min_v1;
theta2 = vtheta2 * min_v2;
}
model {
vmn1
~ uniform(0, 1);
vmn2
~ uniform(0, 1);
vsd1
~ uniform(0, 1);
vsd2
~ uniform(0, 1);
vtheta1 ~ uniform(0, 1);
vtheta2 ~ uniform(0, 1);
for
(i in 1:n1) {
x1[i] - theta1 ~ lognormal(mu1, sgm1);
}
for
(i in 1:n2) {
x2[i] - theta2 ~ lognormal(mu2, sgm2);
}
}
generated quantities {
real
mean1;
real
mean2;
real
diff_means;
real
diff_sds;
real
diff_mus;
real
diff_sgms;
real
mode1;
real
mode2;
real
diff_modes;
mean1 = theta1 + mn1;
mean2 = theta2 + mn2;
diff_means = mean1 - mean2;
diff_sds = sd1 - sd2;
diff_mus = mu1 - mu2;
diff_sgms = sgm1 - sgm2;
mode1 = theta1 + exp(mu1 - sgm1^2);
mode2 = theta2 + exp(mu2 - sgm2^2);
diff_modes = mode1 - mode2;
}
"""
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)
# 読み込みデータの保存
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);
min_v1 = np.min(x)
min_v2 = np.min(y)
sm = stan.build(stan_code, data = {'n1':len(x), 'x1':x,
'n2':len(y), 'x2':y,
'max_v':max_v,
'min_v1':min_v1,
'min_v2':min_v2})
fit = sm.sample()
i_data = az.from_pystan(posterior=fit,
posterior_model=sm)
print('Summary...\n', az.summary(i_data))
fit_dframe = fit.to_frame()
with open('fit_dframe.pkl',
'wb') as f:
pickle.dump(fit_dframe, f)
# MCMCサンプルの保存
リスト２ リスト１の実行で得られたMCMCサンプルの分析（draw_ba_lognormal3p.py）
"""
Yasuharu Okamoto, 2025.07
The script can
be used freely
"""
import numpy
as np
import scipy.stats
as ss
import pickle
import matplotlib.pyplot
as plt
import seaborn as sb
def lognormal_pdf(x,
mu, sgm, theta):
x -= theta
if x <= 0.0:
return 0.0
else:
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)
# データファイルの読み込み
with open('fit_dframe.pkl',
'rb') as f:
fit_dframe
= pickle.load(f)
# MCMCサンプルの読み込み
print(fit_dframe.keys())
df_keys = in_dframe.keys()
print('df_keys
=\n', df_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)
"""
パラメータの点推定値として中央値をとる
"""
mu1 = np.median(fit_dframe['mu1'])
sgm1 = np.median(fit_dframe['sgm1'])
theta1 = np.median(fit_dframe['theta1'])
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 = []
for vx in rng_x:
d1.append(lognormal_pdf(vx, mu1, sgm1,
theta1))
plt.plot(rng_x, d1,
label=x_name)
mu2 = np.median(fit_dframe['mu2'])
sgm2 = np.median(fit_dframe['sgm2'])
theta2 = np.median(fit_dframe['theta2'])
#
データ２の対数正規分布
d2 = []
for vx in rng_x:
d2.append(lognormal_pdf(vx, mu2, sgm2,
theta2))
plt.plot(rng_x, d2,
label=y_name)
"""
事後分布の代表値として中央値をとる
"""
mean1 = np.median(fit_dframe['mean1'])
sd1 = np.median(fit_dframe['sd1'])
mean2 = np.median(fit_dframe['mean2'])
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}:
$\theta$={theta1:.1f}, $\mu={mu1:.1f}, \sigma={sgm1:.1f}$' +
fr', $mean$={mean1:.1f}, $sd$={sd1:.1f},
$mode$={mode1:.1f}' + '\n' +
fr'
{y_name}: $\theta$={theta2:.1f},
$\mu={mu2:.1f}, \sigma={sgm2:.1f}$' +
fr', $mean$={mean2:.1f}, $sd$={sd2:.1f},
$mode$={mode2:.1f}')
plt.ylabel('Density')
plt.legend()
plt.tight_layout()
plt.show()
plt.title(r'Posterior
disgtributions of $\theta$', fontsize=16)
sb.kdeplot(fit_dframe['theta1'],
label=x_name)
sb.kdeplot(fit_dframe['theta2'],
label=y_name)
plt.xlabel(r'$\theta$', fontsize=16)
plt.legend()
plt.show()
"""
平均値パラメータmean1とmean2の事後分布
"""
plt_mns = {x_name:fit_dframe['mean1'],
y_name:fit_dframe['mean2']}
sb.kdeplot(data=plt_mns)
plt.xlabel('$mean$', fontsize=14)
plt.title('Posterior distrbutions
of $mean$s')
plt.show()
diff_means = fit_dframe['diff_means']
p_pos = np.mean(diff_means > 0)
sb.kdeplot(diff_means)
plt.xlabel(f'$mean$[{x_name}] - $mean$[{y_name}]', fontsize=14)
plt.title(f'P($mean$[{x_name}] - $mean$[{y_name}] >
0) = {p_pos}' + '\n' +
fr'Posterior distribution of $mean$[{x_name}] - $mean$[{y_name}]')
plt.show()
"""
標準偏差パラメータsd1とsd2の事後分布
"""
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()
"""
変換パラメータmu1とmu2の事後分布
"""
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()
"""
パラメータsgm1とsgm2の事後分布
"""
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()
"""
パラメータmode1とmode2の事後分布
"""
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()
"""
パラメータの差
mode1-mode2 の事後分布
"""
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()