３パラメータ対数正規分布による独立な２つのデータの平均値のベイズ分析
CmdStanPy対応
確率変数
の対数
が平均
、分散
の正規分布に従うとき、確率変数は３パラメータの対数正規分布に従うという（Johnson et al., 1994）。確率密度関数は次式（１）のように書くことができる。
このとき、平均値、分散、モードは次式で与えられる。
対数正規分布のパラメータをと
ではなく、対数正規分布の平均値
と標準偏差
をパラメータとすると、パラメータと分布との関係が分かり易くなる。
式（２）、式（３）より
とおけば、
であり、
である。
いま、２つの独立なデータ
と
がそれぞれ対数正規分布
と
の母集団から与えられたとき、事後分布を式（７）のようにおく。
事前分布
を一様分布としたときのStanによるMCMCサンプリングのスクリプトをリスト１に示す。ファイルはlognormal3pcspfiles.zipにまとめた。
リスト１ MCMCサンプリングのStanスクリプト（lognormal3p.stan）
/*
Yasuharu
Okamoto, 2025.07,2025.10
The
script can be used freely
*/
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 {
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;
}
まず、変域が
である一様乱数を次のスクリプトで生成している。
vmn1 ~ uniform(0, 1);
vmn2 ~ uniform(0, 1);
vsd1 ~ uniform(0, 1);
vsd2 ~ uniform(0, 1);
vtheta1 ~ uniform(0, 1);
vtheta2 ~ uniform(0, 1);
これらの乱数を、次のスクリプトでデータに合わせている。
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;
これらの調整されたパラメータにより、次のスクリプトでモデルがデータに適用されている。
for (i in 1:n1) {
x1[i] - theta1 ~ lognormal(mu1, sgm1);
}
for (i in 1:n2) {
x2[i] - theta2 ~ lognormal(mu2, sgm2);
リスト１のStanスクリプトを用いて、ベイズ分析を行うPythonスクリプトをリスト２のように用意した。ファイルはlognormal3pcspfiles.zipにまとめている。
リスト２ リスト１のStanスクリプトを用いてベイズ分析を行うPythonスクリプト（ba_lognormal3p_csp.py）
"""
Yasuharu
Okamoto, 2025.07
The script can
be used freely
"""
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,
'max_v':max_v, 'min_v1':min_v1, 'min_v2':min_v2})
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)
sb.kdeplot(fit_dframe['theta2'],
label=y_name)
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)
#
データ２の対数正規分布
d2 = []
for vx in rng_x:
d2.append(lognormal_pdf(vx,
mu2, sgm2, theta2))
plt.plot(rng_x, d2, label=y_name)
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()
リスト１のStanファイル（lognormal3p.stan）を、次のコードでコンパイルしている。
model =
CmdStanModel(stan_file="lognormal3p.stan")
コンパイルされたmodelにより、次のコードでMCMCサンプリングを行っている。
fit = model.sample(data={'n1':len(x),
'x1':x, 'n2':len(y), 'x2':y,
'max_v':max_v, 'min_v1':min_v1, 'min_v2':min_v2})
MCMCサンプリングの統計は、次のコードで出力されている。
print(fit.summary())
Arviz（az）を用いてトレース図を描くために、fitをArvizのInferenceData型に次のコードで変換している。
inf_data = az.from_cmdstanpy(fit)
このinf_dataを用いて、次のコードでトレース図を描いている。
az.plot_trace(inf_data,
var_names=['mean1','sd1','mean2','sd2',
'theta1', 'theta2'], figsize=(12,8))
plt.tight_layout()
plt.show()
パラメータのMCMCサンプルを取り出すために、fitをPandasのDataFrame型に次のコードで変換している。
fit_dframe = fit.draws_pd()
パラメータのMCMCサンプル値は、変数名を文字列のキーとして得ることができる。
次のコードは、パラメータmu1のMCMCサンプルを取り出して、その中央値を求めるものである。
mu1 = np.median(fit_dframe['mu1'])
リスト２のPythonスクリプトは、CmdStanPy対応である。CmdStanPyのインストールの簡単な説明を、このウェブサイトに用意した。
リスト２のファイル（ba_lognormal3p_csp.py）とリスト１のファイルを、入力データファイルとともに同じフォルダにおいて、次のコマンドで実行する。
(stan) *****/lognormal3pcspfiles$
python ba_lognormal3p_csp.py
スクリプトが実行されると、入力データファイル名の設定が求められる。
(stan) *****/lognormal3pcspfiles$
python ba_lognormal3p_csp.py
Input data
file(*.xlsx) = sampledata.xlsx
上の例では、sampledata.xlsxが設定されている。
入力データファイルは、Excelファイルとして用意する（図１）。
図１A
図１B
1行目に変数名を設定する。１列目は空欄でもよい。データは２行目から設定する。１列目はデータの識別用である。２列目と３列目にデータ値を設定する。データ数が２つのグループで異なるときは、少ない方のデータは空欄にしておけばよい。空欄はスクリプトの実行時に無視される。
入力データファイル名を設定して、Enterキーを押すと、入力データファイルが読み込まれて、リスト１のStanスクリプトのコンパイルが始まる。コンパイルには多少の時間が掛かる。コンパイル後、MCMCサンプリングが始まり、サンプリング後、トレース図が表示される（図２）。
図２
図２のWindowを閉じると、図３のグラフが表示される。
図３
パラメータ
の事後分布と、差
の事後分布、および差が正である事後確率が表示されている。正である事後確率が
であるので２つのパラメータ
と
に差はないと考えられる。なお、本ウェブサイトにおいては、変数ｘの添字はAあるいは１、変数ｙの添字はBあるいは２が用いられている。
図３のWindowを閉じると、図４のグラフが表示される。
図４
平均値パラメータmeanの事後分布、および差
の事後分布と差が正である事後確率が表示されている。
であるので、パラメータmean[VarA]はパラメータmean[VarB]より確かに高いと言える。事後確率が、1.0と小数点１位までの表記であるのは、スクリプト言語の仕様である。
図４のWindowを閉じると、図５のグラフが表示される。
図５
モードmodeの事後分布と、モードの差の事後分布、および差が正である事後確率が表示されている。差が正である事後確率は
であるので、差はないと考えられる。平均値には明瞭な差が認められたが（図４）、モードには差が認められない（図５）。
図５のWindowを閉じると、図６のグラフが表示される。
図６
標準偏差の事後分布と、標準偏差の差の事後分布、および差が正である事後確率が表示されている。差が正である事後確率は、
であるので、VarAの標準偏差の方が確かに大きいと言えそうである。確率が1.0と小数点以下１位で切れているのは、スクリプト言語の仕様である。
図６のWindowを閉じると、図７のグラフが表示される。
図７
データのヒストグラムと、それぞれのデータに対する３パラメータ対数正規分布の確率密度関数が表示されている。確率密度関数のパラメータ値は、事後分布の中央値を点推定値としている。Gelmanら（2021）は、平均値より安定しているとして、中央値を代表値とすることを勧めている。
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.