対数正規分布による独立な２つのデータの平均値のベイズ分析
CmdStanPy対応、2026.03改訂
確率変数
の対数
が平均
、分散
の正規分布に従うとき、確率変数は対数正規分布に従うと言い、確率密度関数は次式（１）で与えられる（Gelman et al., 2014）。
このとき、平均値、分散、モードは次式で与えられる。
対数正規分布では、平均値（２）とモード（４）が異なる。データ全体に対する代表値としての平均値と、観察頻度の高い領域の位置を示すモード（最頻値）の２つの値が観察値の代表値として押さえておく必要がある。
対数正規分布のパラメータを対数変換した値の正規分布のパラメータと
ではなく、対数正規分布の平均値
と標準偏差
をパラメータとすると、パラメータと分布との関係が直接的になる。
式（２）、式（３）より
とおけば、
であり、
である。
対数正規分布による分析結果は、パラメータ平均値（２）、標準偏差（分散（３）の平方根）、モード（４）を示すと２つの分布の比較ができる。
いま、２つの独立なデータ
と
がそれぞれ対数正規分布
と
の母集団から与えられたとき、事後分布を式（７）のようにおく。
事前分布
を一様分布としたときのStanスクリプトをリスト１に示す。ファイルはlognormalcspfiles.zipにまとめた。
リスト１ 対数正規分布のStanスクリプト（ba_lognormal.stan）
data {
int n1;
array[n1] real x1;
int n2;
array[n2] real x2;
}
transformed data {
real log_max_x1;
real log_min_x1;
real log_max_x2;
real log_min_x2;
log_max_x1 = log(max(x1));
log_min_x1 = log(min(x1));
log_max_x2 = log(max(x2));
log_min_x2 = log(min(x2));
}
parameters {
real<lower=log_min_x1,
upper=log_max_x2> mu1;
real<lower=log_min_x2,
upper=log_max_x2> mu2;
real<lower=0.0001,
upper=log_max_x1 - log_min_x1> sgm1;
real<lower=0.0001,
upper=log_max_x2 - log_min_x2> sgm2;
}
model {
mu1 ~ uniform(log_min_x1,
log_max_x1);
sgm1 ~ uniform(0.0001,
log_max_x1 - log_min_x1);
mu2 ~ uniform(log_min_x2,
log_max_x2);
sgm2 ~ uniform(0.0001,
log_max_x2 - log_min_x2);
x1 ~ lognormal(mu1, sgm1);
x2 ~ lognormal(mu2, sgm2);
}
generated quantities {
real mn1;
real mn2;
real sd1;
real sd2;
real diff_mns;
real diff_sds;
real diff_mus;
real diff_sgms;
real mode1;
real mode2;
real diff_modes;
mn1 = exp(mu1 +
0.5*(sgm1^2));
mode1 = exp(mu1 - (sgm1^2));
sd1 = (exp(2*mu1 + (sgm1^2))
* (exp(sgm1^2) - 1.0))^0.5;
mn2 = exp(mu2 +
0.5*(sgm2^2));
mode2 = exp(mu2 - (sgm2^2));
sd2 = (exp(2*mu2 + (sgm2^2))
* (exp(sgm2^2) - 1.0))^0.5;
diff_mns = mn1 - mn2;
diff_sds = sd1 - sd2;
diff_modes = mode1 - mode2;
diff_sds = sd1 - sd2;
diff_mus = mu1 - mu2;
diff_sgms = sgm1 - sgm2;
diff_modes = mode1 - mode2;
}
上のStanで記述されたモデルは、リスト２のスクリプトで用いられている。
リスト２ リスト１のStanスクリプトを利用するPythonスクリプト（ba_lognormal.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)
model =
CmdStanModel(stan_file="ba_lognormal.stan") # コンパイル
fit = model.sample(data={'n1':len(x),
'x1':x,'n2':len(y), 'x2':y})
print(fit.summary())
# MCMCサンプリングの統計
inf_data = az.from_cmdstanpy(fit)
# Arviz(az)用のデータに変換
# トレースプロット
az.plot_trace(inf_data,
var_names=['mn1','sd1','mn2','sd2'])
plt.tight_layout()
plt.show()
fit_dframe = fit.draws_pd() # PandasのDataFrame型に変換
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
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)
print(fit_dframe.keys())
"""
平均値パラメータmn1とmn2の事後分布
"""
plt.figure(figsize=(10,7))
plt.subplot(2,2,1)
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.tight_layout()
plt.subplot(2,2,2)
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.tight_layout()
"""
標準偏差パラメータsd1とsd2の事後分布
"""
plt.subplot(2,2,3)
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.tight_layout()
plt.subplot(2,2,4)
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.tight_layout()
plt.show()
"""
mode1とmode2の事後分布
"""
plt.figure(figsize=(10,4))
plt.subplot(121)
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.tight_layout()
plt.subplot(122)
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.tight_layout()
plt.show()
plt.hist(x, density=True,
label=x_name, color='r', alpha=0.5)
plt.hist(y, density=True,
label=y_name, color='g', alpha=0.5)
"""
パラメータの点推定値として中央値をとる
"""
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,
color='r',)
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,
color='g')
"""
事後分布の代表値として中央値をとる
"""
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.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()
リスト１のスクリプトファイル、リスト２のスクリプトファイル、および入力データファイルを同じフォルダに置き、CmdStanPyのインストールされた環境で実行する。CmdStanPyのインストールの簡単な説明をこのウェブサイトに用意した。
リスト２のスクリプトを実行すると、入力データファイル名を聞いてくる。
(stan) ****/lognormalcspfiles$ python
ba_lognormal.py
Input data file(*.xlsx) = data.xlsx
(mystan)
*****/lognormalcspfiles$ python ba_lognormal.py
Input data
file(*.xlsx) = data.xlsx
上の例では、入力データファイル名としてdata.xlsxが設定されている。入力データファイルは、図１a、図１bのExcelファイルの形式で用意する。
図１a
図１b
Excelファイルは、1行目に変数名を設定する。１列目は空欄でもよい。データは２行目から設定する。１列目はデータの識別用である。２列目と３列目にデータ値を設定する。データ数が２つのグループで異なるときは、少ない方のデータは空欄にしておけばよい。空欄はスクリプトの実行時に次のスクリプトで削除される。
x = x[np.isnan(x)==False] # 空セルデータの削除
y = y[np.isnan(y)==False] # 空セルデータの削除
入力データファイル名を設定してEnterキーを押すと、データファイルが読み込まれ、続いてStanスクリプトのコンパイルが始まる。コンパイルには多少時間が掛かるが、コンパイルが終了すると、MCMCサンプリングが行われる。MCMCサンプリングが終了すると、トレース図が表示される（図２）。
図２
図２のWindowの右上角のX印をクリックして閉じると、図３のグラフが表示される。
図３
上の段の左には、パラメータmn1とmn2の事後分布がmn[A]とmn[B]として表示されている。１および２が、AおよびBとなっているのは、スクリプトの作成の手順でそのようになったもので、深い意味はない。上段右側には、差mn[A]-mn[B]の事後分布と、差の事後分布の正の領域の事後確率が表示されている。
であるので、差mn[A]-mn[B]の事後分布の正である領域は１％以下であり、差mn[A]-mn[B]は負であると考えられる。
下の段には、標準偏差sd[A]とsd[B]の事後分布が表示されている。下段右側には差sd[A]-sd[B]の事後分布と差sd[A]-sd[B]が正である事後確率が表示されている。
であるので（小数点以下１位までの表記は、スクリプト言語の仕様による）、差sd[A]-sd[B]はほぼ確実に負であると考えられる。
図３のWindowを閉じると、図４のグラフが表示される。
図４
モードの事後分布が描かれている。右側には、差mode[A]-mode[B]の事後分布と、差が正である事後確率が表示されている。差が正である確率は
とかなり高いので、mode[A]はmode[B]より大きい傾向にあると考えられる。
図４のWindowを閉じると、図５のグラフが表示される。
図５
データのヒストグラムと２つのデータに対する対数正規分布のグラフが描かれている。対数正規分布のパラメータは、との事後分布の中央値を点推定値とした。Gelmanら（2021）は、点推定値として、平均値より安定した中央値を勧めている。グラフのタイトル部に表示されている点推定値としての中央値を見ると、データA（すなわち、データ1あるいはデータX）の方がデータB（すなわち、データ２あるいはY）より平均値が低い（図３上段）が、モードはデータAの方がデータBより大きい（図４）。これは、標準偏差（分散）の違い（図３下段）による。
図５のWindowを閉じると、スクリプトの実行終了である。
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