３パラメータ対数正規分布による独立な２つのデータの平均値のベイズ分析
CmdStanPy対応
2026.03改訂
確率変数
の対数
が平均
、分散
の正規分布に従うとき、確率変数は３パラメータの対数正規分布に従うという（Johnson et al., 1994）。確率密度関数は次式（１）のように書くことができる。
このとき、平均値、分散、モードは次式で与えられる。
対数正規分布のパラメータをと
ではなく、対数正規分布の平均値
と標準偏差
をパラメータとすると、パラメータと分布との関係が分かり易くなる。
式（２）、式（３）より
とおけば、
であり、
である。
いま、２つの独立なデータ
と
がそれぞれ３パラメータ対数正規分布
と
の母集団から与えられたとき、事後分布を式（７）のようにおく。
事前分布
を一様分布としたときのStanによるMCMCサンプリングのスクリプトをリスト１に示す。ファイルはlognormal3pcspfiles.zipにまとめた。
リスト１ MCMCサンプリングのStanスクリプト（lognormal3p.stan）
data {
int n1;
array[n1] real x1;
int n2;
array[n2] real x2;
}
transformed data {
real min_x1;
real min_x2;
min_x1 = min(x1);
min_x2 = min(x2);
}
parameters {
real mu1;
real mu2;
real<lower=0.0001,
upper=10000.0> sgm1;
real<lower=0.0001,
upper=10000.0> sgm2;
real<lower=0.0,
upper=min_x1> theta1;
real<lower=0.0,
upper=min_x2> theta2;
}
model {
mu1 ~ uniform(-1000, 1000);
mu2 ~ uniform(-1000, 1000);
sgm1 ~ uniform(0.0001,
10000);
sgm2 ~ uniform(0.0001,
10000);
theta1 ~ uniform(0, min_x1);
theta2 ~ uniform(0, min_x2);
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 mode1;
real sd1;
real mean2;
real mode2;
real sd2;
real diff_means;
real diff_sds;
real diff_mus;
real diff_sgms;
real diff_modes;
mean1 = theta1 + exp(mu1 +
0.5*(sgm1^2));
mode1 = theta1 + exp(mu1 -
(sgm1^2));
sd1 = (exp(2*mu1 + (sgm1^2))
* (exp(sgm1^2) - 1.0))^0.5;
mean2 = theta2 + exp(mu2 +
0.5*(sgm2^2));
mode2 = theta2 + exp(mu2 -
(sgm2^2));
sd2 = (exp(2*mu2 + (sgm2^2))
* (exp(sgm2^2) - 1.0))^0.5;
diff_means
= mean1 - mean2;
diff_sds
= sd1 - sd2;
diff_mus
= mu1 - mu2;
diff_sgms
= sgm1 - sgm2;
diff_modes
= mode1 - mode2;
}
Generated quantities部において、式（２）、（３）、（４）により、平均値、標準偏差（分散の平方根）、モードを算出している。
リスト１のStanスクリプトを用いて、ベイズ分析を行うPythonスクリプトをリスト２のように用意した。ファイルはlognormal3pcspfiles.zipにまとめている。
リスト２ リスト１のStanスクリプトを用いてベイズ分析を行うPythonスクリプト（ba_lognormal3p_csp.py）
from cmdstanpy
import CmdStanModel
import numpy
as np
import pandas as pd
import matplotlib.pyplot
as plt
import seaborn as sb
import arviz
as az
inf_nm = input('Input data file(*.xlsx) = ')
dframe = pd.read_excel(inf_nm)
# データの読み込み
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)
#
Stanスクリプトのコンパイル
model = CmdStanModel(stan_file="lognormal3p.stan")
#
MCMCサンプリング
fit = model.sample(data={'n1':len(x),
'x1':x, 'n2':len(y), 'x2':y},
adapt_delta=0.99)
#'max_v':max_v, 'min_v1':min_v1, 'min_v2':min_v2},)
print(fit.diagnose())
print(fit.summary())
# MCMCサンプリングの統計
inf_data = az.from_cmdstanpy(fit)
# Arviz(az)用のデータに変換
# inf_dataを用いてのトレースプロット
az.plot_trace(inf_data, var_names=['mean1','sd1','mean2','sd2',
'theta1', 'theta2'], figsize=(12,8))
plt.tight_layout()
plt.show()
fit_dframe = fit.draws_pd() # PandasのDataFrame型に変換
print(fit_dframe.keys())
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
"""
パラメータの点推定値として中央値をとる
"""
mu1 = np.median(fit_dframe['mu1'])
sgm1 = np.median(fit_dframe['sgm1'])
theta1 = np.median(fit_dframe['theta1'])
mu2 = np.median(fit_dframe['mu2'])
sgm2 = np.median(fit_dframe['sgm2'])
theta2 = np.median(fit_dframe['theta2'])
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'])
"""
thetaの事後分布
"""
plt.figure(figsize=(10,5))
plt.subplot(1,2,1)
plt.title(r'Posterior
distributions of $\theta$', fontsize=16)
sb.kdeplot(fit_dframe['theta1'],
label=x_name+'(theta1)')
sb.kdeplot(fit_dframe['theta2'],
label=y_name+'(theta2)')
plt.xlabel(r'$\theta$', fontsize=16)
plt.legend()
plt.tight_layout()
plt.subplot(1,2,2)
p_pos = np.mean(fit_dframe['theta1'] - fit_dframe['theta2']
> 0)
plt.title(r'Posterior
disgtribution of $\theta1-\theta2$' + '\n' +
fr'$P(\theta1-\theta2>0)$={p_pos:.3f}', fontsize=16)
sb.kdeplot(fit_dframe['theta1']
- fit_dframe['theta2'])
plt.xlabel(r'$\theta1-\theta2$', fontsize=16)
plt.legend()
plt.tight_layout()
plt.show()
"""
平均値パラメータmean1とmean2の事後分布
"""
plt.figure(figsize=(10,5))
plt.subplot(1,2,1)
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.subplot(1,2,2)
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()
"""
モードmode1とmode2の事後分布
"""
plt.figure(figsize=(10,5))
plt.subplot(1,2,1)
plt_mds = {x_name:fit_dframe['mode1'],
y_name:fit_dframe['mode2']}
sb.kdeplot(data=plt_mds)
plt.xlabel('$mode$', fontsize=14)
plt.title('Posterior distrbutions
of $mode$s')
plt.subplot(1,2,2)
diff_modes = fit_dframe['diff_modes']
p_pos = np.mean(diff_modes > 0)
sb.kdeplot(diff_modes)
plt.xlabel(f'$mode$[{x_name}] - $mode$[{y_name}]', fontsize=14)
plt.title(f'P($mode$[{x_name}] - $mode$[{y_name}] >
0) = {p_pos}' + '\n' +
fr'Posterior distribution of $mode$[{x_name}] - $mode$[{y_name}]')
plt.show()
"""
標準偏差パラメータsd1とsd2の事後分布
"""
plt.figure(figsize=(10,5))
plt.subplot(1,2,1)
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.subplot(1,2,2)
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()
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, c='g')
#
データ２の対数正規分布
d2 = []
for vx in rng_x:
d2.append(lognormal_pdf(vx, mu2, sgm2,
theta2))
plt.plot(rng_x, d2,
label=y_name, c='r')
#plt.hist([x, y], density=True,
label=[x_name, y_name])
plt.hist(x, density=True, label=x_name, color='g', alpha=0.5)
plt.hist(y, density=True, label=y_name, color='r', alpha=0.5)
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()
リスト１のファイル、リスト２のファイル、および入力データファイルを同じフォルダにおいて、CmdStanPyのインストールされた環境においてリスト２のスクリプトを実行する。CmdStanPyのインストールの簡単な説明を、このウェブサイトなどに用意した。
リスト２のスクリプトを実行すると、入力データファイル名を聞いてくる。
(stan) ****/lognormal3pcspfiles$
python ba_lognormal3p_csp.py
Input data file(*.xlsx) =
sampledata.xlsx
上の例では、sampledata.xlsxが設定されている。
入力データファイルは、Excelファイルとして用意する（図１A、図１B）。
図１A
図１B
1行目に変数名を設定する。１列目は空欄でもよい。データは２行目から設定する。１列目はデータの識別用である。２列目と３列目にデータ値を設定する。データ数が２つのグループで異なるときは、少ない方のデータは空欄にしておけばよい。空欄はスクリプトの実行時に無視される。
x = x[np.isnan(x)==False] # 空セルデータの削除
y = y[np.isnan(y)==False] # 空セルデータの削除
入力データファイル名を設定して、Enterキーを押すと、入力データファイルが読み込まれて、リスト１のStanスクリプトのコンパイル・ビルドが始まる。このコンパイル・ビルドには多少の時間が掛かる。コンパイル後、MCMCサンプリングが始まり、サンプリング後、トレース図が表示される（図２）。
図２
図２のWindowを閉じると、図３のグラフが表示される。
図３
パラメータ
の事後分布と、差
の事後分布、および差が正である事後確率が表示されている。正である事後確率が
であるので２つのパラメータ
と
に差があるとは言い難い。なお、本ウェブサイトにおいては、変数ｘの添字はAあるいは１、変数ｙの添字はBあるいは２が用いられている。
図３のWindowを閉じると、図４のグラフが表示される。
図４
平均値パラメータmeanの事後分布、および差
の事後分布と差が正である事後確率が表示されている。
であるので、パラメータmean[VarA]はパラメータmean[VarB]より確かに大きい値であると言える。事後確率が、1.0と小数点１位までの表記であるのは、スクリプト言語の仕様である。
図４のWindowを閉じると、図５のグラフが表示される。
図５
モードmodeの事後分布と、モードの差の事後分布、および差が正である事後確率が表示されている。差が正である事後確率は
であるので、差はないと考えられる。平均値には明瞭な差が認められたが（図４）、モードには差が認められない（図５）。
図５のWindowを閉じると、図６のグラフが表示される。
図６
標準偏差の事後分布と、標準偏差の差の事後分布、および差が正である事後確率が表示されている。差が正である事後確率は、
であるので、VarAの標準偏差の方が確かに大きいと言えそうである。確率が1.0と小数点以下１位で切れているのは、スクリプト言語の仕様である。
図６のWindowを閉じると、図７のグラフが表示される。
図７
データのヒストグラムと、それぞれのデータに対する３パラメータ対数正規分布の確率密度関数が表示されている。確率密度関数のパラメータ値は、事後分布の中央値を点推定値としている。Gelmanら（2021）は、平均値より安定しているとして、中央値を代表値とすることを勧めている。
3パラメータ対数正規分布がデータに良く当てはまっていること、2つの３パラメータ対数正規分布のモードがほぼ同じであること、平均値と標準偏差が異なることがよく表されていると言える。
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.