対数正規分布による独立な２つのデータの平均値のベイズ分析
CmdStanPy対応
確率変数
の対数
が平均
、分散
の正規分布に従うとき、確率変数は対数正規分布に従うと言い、確率密度関数は次式（１）で与えられる（Gelman et al., 2014）。
このとき、平均値、分散、モードは次式で与えられる。
対数正規分布では、平均値（２）とモード（４）が異なる。データ全体に対する代表値としての平均値と、観察頻度の高い領域の位置を示すモード（最頻値）の２つの値が観察値の代表値として押さえておく必要がある。
対数正規分布のパラメータを対数変換した値の正規分布のパラメータと
ではなく、対数正規分布の平均値
と標準偏差
をパラメータとすると、パラメータと分布との関係が直接的になる。
式（２）、式（３）より
とおけば、
であり、
である。
対数正規分布による分析結果は、パラメータとしてはとを報告すれば十分であろう。
いま、２つの独立なデータ
と
がそれぞれ対数正規分布
と
の母集団から与えられたとき、事後分布を式（７）のようにおく。
事前分布
を一様分布としたときのStanスクリプトをリスト１に示す。ファイルはlognormalcspfiles.zipにまとめた。
リスト１ 対数正規分布のStanスクリプト（ba_lognormal.stan）
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宣言部で変域の広いパラメータに
出される警告を、スケーリングで避ける
*/
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);
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;
mode1 = exp(mu1 - sgm1^2);
mode2 = exp(mu2 - sgm2^2);
diff_modes = mode1 - mode2;
}
乱数として変域を
に制限したものをmodelブロックでまず生成している。
vmn1 ~ beta(1, 1);
vmn2 ~ beta(1, 1);
vsd1 ~ beta(1, 1);
vsd2 ~ beta(1, 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);
これらのパラメータにより、データが対数正規分布で表されるモデルが設定されている。
x1 ~ lognormal(mu1, sgm1);
x2 ~ lognormal(mu2, sgm2);
上の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)
max_v = np.max(x) if np.max(x) >
np.max(y) else np.max(y);
model =
CmdStanModel(stan_file="ba_lognormal.stan") # コンパイル
fit = model.sample(data={'n1':len(x),
'x1':x,'n2':len(y), 'x2':y,
'max_v':max_v})
# サンプリング
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,7))
plt.subplot(2,4,(1,2))
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(2,4,(3,4))
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.subplot(2,4,(6,7))
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()
リスト１のStanスクリプトは、次のコードでコンパイルしている。
model = CmdStanModel(stan_file="ba_lognormal.stan")
コンパイルしたmodelの実行によるMCMCサンプリングは、次のコードで行っている。
fit = model.sample(data={'n1':len(x),
'x1':x,'n2':len(y), 'x2':y,
'max_v':max_v})
MCMCサンプリングの統計は、次のコードで表示している。
print(fit.summary())
MCMCサンプリングのトレース図を描くために、MCMCサンプリングの結果を、次のコードでArviz（az）のInferenceData型に変換している。
inf_data = az.from_cmdstanpy(fit)
変換したデータinf_dataを用いて、次のスクリプトでトレース図を描いている。
az.plot_trace(inf_data,
var_names=['mn1','sd1','mn2','sd2'])
plt.tight_layout()
plt.show()
MCMCサンプルを取り出すために、MCMCサンプルをPandasのDataFrame型に次のコードで変換している。
fit_dframe = fit.draws_pd()
パラメータのMCMCサンプルの値は、パラメータ名を文字列で表してfit_dframeのキーとして用いることによりできる。
例えば、パラメータmn1のサンプルは、次のコードで得ることができる。
fit_dframe['mn1']
次のスクリプトは、パラメータ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.tight_layout()
リスト２のスクリプトは、CmdStaPyのインストールされた環境において、次のコマンドで実行できる。CmdStanPyのインストールの簡単な説明をこのウェブサイトに用意した。
(mystan) *****/lognormalcspfiles$
python ba_lognormal.py
スクリプトが実行されると、入力データファイル名の設定が求められる。次の例では、data.xlsxが設定されている。
(mystan) *****/lognormalcspfiles$
python ba_lognormal.py
Input data
file(*.xlsx) = data.xlsx
入力データファイルは、図１のExcelファイルの形式で用意する。
図１a
図１b
Excelファイルは、1行目に変数名を設定する。１列目は空欄でもよい。データは２行目から設定する。１列目はデータの識別用である。２列目と３列目にデータ値を設定する。データ数が２つのグループで異なるときは、少ない方のデータは空欄にしておけばよい。空欄はスクリプトの実行時に無視される。
データファイルとリスト１およびリスト２のファイルは、同じフォルダ内に置いておく。
入力データファイル名を設定してEnterキーを押すと、データファイルが読み込まれ、続いてStanスクリプトのコンパイルが始まる。コンパイルには多少時間が掛かるが、コンパイルが終了すると、MCMCサンプリングが行われる。MCMCサンプリングが終了すると、トレース図が表示される（図２）。
図２
図２のWindowを閉じると、図３のグラフが表示される。
図３
上の段の左には、パラメータ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]より大きい傾向にあると考えられる。
下段には、データのヒストグラムと２つのデータに対する対数正規分布のグラフが描かれている。対数正規分布のパラメータは、との事後分布の中央値を点推定値とした。Gelmanら（2021）は、点推定値として、平均値より安定した中央値を勧めている。グラフのタイトル部に表示されている点推定値としての中央値を見ると、データA（すなわち、データ1あるいはデータX）の方がデータB（すなわち、データ２あるいはY）より平均値が低いが、モードはデータAの方がデータBより大きい。これは、標準偏差（分散）の違いによる。
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