対数正規分布による独立な２つのデータの平均値のベイズ分析
ベイズ分析の基本的考え方の説明を岡本（2019）、pp.169-177、岡本（2014）、pp.227-232で行っている。
確率変数
の対数
が平均
、分散
の正規分布に従うとき、確率変数は対数正規分布に従うと言い、確率密度関数は次式（１）で与えられる（Gelman et al., 2014）。
このとき、平均値、分散、モードは次式で与えられる。
対数正規分布のパラメータをと
ではなく、対数正規分布の平均値
と標準偏差
をパラメータとすると、パラメータと分布との関係が分かり易くなる。
式（２）、式（３）より
とおけば、
であり、
である。
いま、２つの独立なデータ
と
がそれぞれ対数正規分布
と
の母集団から与えられたとき、事後分布を式（７）のようにおく。
事前分布
を一様分布としたときのStanによるMCMCサンプリングのスクリプトをリスト１に示す。
リスト１によるMCMCサンプルを用いてデータ分析を行うスクリプト例を、リスト２に示す。
これらのスクリプトファイルなどは圧縮ファイル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.02.09）ので、リスト２はWindows上で直接実行した。
リスト１の実行後、リスト２を実行すると図２のグラフが表示される。
図２
パラメータの事後分布が描かれている。
図２のWindowを閉じると、図３のグラフが表示される。
図３
パラメータの差の事後分布が表示される。タイトル欄に確率
の値が表示されている。
の事後確率が0.5％以下であることが分かる。
および
は、データAおよびデータBに対するパラメータであることを表す（図１）。
であるので、
である。連続な確率変数が特定の値を取る確率は０である。
図３のWindowを閉じると、図４のグラフが表示される。
図４
パラメータの事後分布が描かれている。
図４のWindowを閉じると、図５のグラフが表示される。
図５
の差の事後分布が描かれている。
の事後確率は０と見なされる値であることがタイトル欄に表示されている。
図５のWindowを閉じると、図６のグラフが表示される。
図６
パラメータの事後分布のグラフである。
図６のWindowを閉じると、図７のグラフが表示される。
図７
パラメータの差の事後分布のグラフである。
の事後確率が約10％であることがタイトル欄に表示されている。パラメータの差に比べると、の差は不確かである（図３）。
図７のWindowを閉じると、図８のグラフが表示される。
図８
パラメータの事後分布である。
図８のWindowを閉じると、図９のグラフが表示される。
図９
パラメータの差の事後分布である。
の事後確率が０と見なせることがタイトル欄に表示されている。
図９のWindowを閉じると、図10のグラフが表示される。
図10
モード（4）の事後分布のグラフが表示されている。
図10のWindowを閉じると、図１１のグラフが表示される。
図12
モードの差の事後分布のグラフである。
の事後確率が約90％であることがタイトル欄に表示されている。モードはAの方が右側にあると言えそうである。
平均値はBの方が大きかった（図６）。モードと平均値の大小関係の逆転は、標準偏差（分散）の違いによる（図４）。対数正規分布は左右対称な分布ではないので、平均だけでは分布の相対的関係を表すことができない。
図12のWindowを閉じると、図13のグラフが表示される。
図10
曲線は、タイトル欄に示されているパラメータ値で描かれた対数正規分布のグラフである。パラメータ値は事後分布の中央値を採用している（Gelman et al.,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
岡本安晴（2019）「いまさら聞けないPythonでデータ分析」丸善出版
リスト１ MCMCサンプリングスクリプト（ba_means.py）
import numpy
as np
import pandas as pd
import stan
import pickle
stan_code = """
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;
/*
パラメータのスケーリング
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);
}
model {
/*
一様分布を宣言するといろいろ警告メッセージが
出されることがあるが、同じ分布をbeta分布で
与えると警告を避けることができるようである
*/
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;
}
"""
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);
sm = stan.build(stan_code, 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)
# MCMCサンプルの保存
リスト２ リスト１の実行で得られたMCMCサンプルの分析（draw_means.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) # データファイルの読み込み
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)
# MCMCサンプルの読み込み
print(fit_dframe.keys())
"""
平均値パラメータmn1とmn2の事後分布
"""
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()
"""
標準偏差パラメータ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()
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()
"""
パラメータの点推定値として中央値をとる
"""
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()