ベイズ分析による正規分布、対数正規分布、３パラメータ対数正規分布の比較例
2026.06.15改訂
モデルの予測精度を表すベイズ的指標として、LOO-CV(leave-one-out cross-validation)があげられる（Gelman et al., 2014; Lambert, 2018）。予測精度（Expected log pointwise predictive density for a new dataset）は、以下のように与えられる（Vehtari et al., 2016）。
LOO-CVの場合：
である。
が大きいことは、
が全体として大きいことを意味するので、モデルのデータへの当てはまりが良いことを示していると考えられる。
確率は、データ
を除いた
について求められる（Leave One Out）ので、
にわたって計算することになり時間が掛かる。Pareto smoothed importance sampling (PSIS) では、次式により計算時間の短縮が図られている。
ここで、
はvector weight of Pareto smoothed importance samplingである。
これにより、
の計算値は次式で与えられる。
PSIS-LOO-CV（以後、LOOと記す場合がある）は、Arvizの関数looにより求めることができる。Arvizの関数を用いるときは、Stanスクリプトにおいて変数log_likを用意しておいて
の値を設定しておくとよい。例えば、リスト１では以下のようになっている。
generated quantities {
array[n] real log_lik;
// log likelihood for
LOO
for (i
in 1:n) {
log_lik[i] = normal_lpdf(x[i] | mu, sgm);
}
}
上のコードを含むStanスクリプトの実行とLOO-CVの計算は、以下のようである（リスト４より）。
def init_f():
return {'mu':np.mean(x), 'sgm':np.std(x)}
sm_n = CmdStanModel(stan_file="normalmodel.stan") # コンパイル
fit_n = sm_n.sample(data={'n':len(x), 'x':x},
chains=4, inits=[init_f()]*4)
i_data_n = az.from_cmdstanpy(fit_n, log_likelihood=True) # Arviz(az)用のデータに変換
results = az.loo(i_data_n)
print(results.elpd)
Arviz用のデータ型にMCMCサンプルを変換するときに、関数from_cmdstanpyに引数log_likelihood=Trueを設定しておく。
LOO-CVの値は、関数looの戻り値のアトリビュートelpdで得ることができる。
デモンストレーション用に３つのデータ、データ１（正規分布）、データ２（対数正規分布）、データ３（３パラメータ対数正規分布）を用意した。以下に分析例を示す。
データdata_normal.txtについて考える。ファイルは、comp_modelsfile.zipにまとめた。
データdata_normal.txtのヒストグラムは、図１のようである。
図１
正規分布モデル（リスト１）、対数正規分布モデル（リスト２）、３パラメータ対数正規分布モデル（リスト３）のStanスクリプトを用意して、PSIS-LOO-CVを求めるPythonスクリプトを作成した（リスト４）。
リスト４のスクリプトを実行すると、入力データファイル名の設定が求められるので、図１のファイル名を設定する。
(stan) ****/comp_modelsfile$
python comp_models.py
Input File Name(*.txt) =
data_normal.txt
ファイル名を設定すると、計算が始まる。
計算が終了すると、図２のグラフが表示される。
図２
入力データのヒストグラム（ビンが異なるので、図１のヒストグラムとは異なった形状である）と、事後分布の中央値をパラメータの推定値とする３つのモデルの確率密度関数のグラフが描かれている。中央値を推定値とするのは、Gelmanら（2014）に従ったものである。
正規分布（赤一点鎖線）のLOOが最大である。図１のデータには、正規分布が最もよく当てはまっているといえる。
対数正規分布と３パラメータ対数正規分布のグラフは、データの分布より少し左に寄っているようであるが、これは両者の分布が右に裾が伸びた形状によるものと思われる。
図２のWindowを閉じると、スクリプトの実行終了となり、以下のようなメッセージが端末に出力される。
LOO.txt was saved.
(stan) ****/comp_modelsfile$
「ファイルが保存された」とあるので、内容を確認してみる。
LOO.txt was saved.
(stan) ****/comp_modelsfile$
cat LOO.txt
LOO-CV
Normal
-4390.240
LogNormal
-4395.643
3P LogNormal
-4396.455
(stan) ****/comp_modelsfile$
３つのモデルに対するLOO-CVが出力されている。
データdata_lognormal.txtについて考える。ファイルは、comp_modelsfile.zipにまとめた。
データdata_lognormal.txtのヒストグラムは、図３のようである。
図３
正規分布モデル（リスト１）、対数正規分布モデル（リスト２）、３パラメータ対数正規分布モデル（リスト３）のStanスクリプトを用意して、PSIS-LOO-CVを求めるPythonスクリプトを作成した（リスト４）。
リスト４のスクリプトを実行すると、入力データファイル名の設定が求められるので、図３のファイル名を設定する。
(stan) ****/comp_modelsfile$
python comp_models.py
Input File Name(*.txt) =
data_lognormal.txt
ファイル名を設定すると、計算が始まる。
計算が終了すると、図４のグラフが表示される。
図４
入力データのヒストグラムと３つのモデルの確率密度関数のグラフである。ヒストグラムは、図３のものとはビンが異なるので形状が少し異なっている。
３つの分布のグラフとそれぞれのLOOの値が表示されている。対数正規分布のLOOの値が最も大きいので、データに最もよく適合していると言える。３パラメータ対数正規分布は、対数正規分布にシフトパラメータｔを追加したものであるので、対数正規分布が適合しているときは３パラメータ対数正規分布も適合しているが余分なパラメータをもつモデルとなるので、ＬＯＯの値は対数正規分布より小さくなる。
図４のグラフを見ると、データのヒストグラムは、正規分布より右側の裾においてもデータが認められ、これは正規分布の観点からは外れ値として扱われることがある。また、左側の裾では、正規分布の観点からはデータが存在することが期待されるが欠損値があるかのような様子である。外れ値のときは、正規分布を基準にして分析対象から外すのであれば、欠損値があると見なされるときは欠損値に対応するデータを追加する必要があると思われるが、データに除外とか追加の操作を加えることは不自然と思われる。対数正規分布を想定すれば、データに手を加えずにそのまま分析することができる。
図４のWindowを閉じると、スクリプトの実行終了となり、以下のようなメッセージが端末に出力される。
LOO.txt was saved.
(stan) ****/comp_modelsfile$
「ファイルが保存された」とあるので、内容を確認してみる。
LOO.txt was saved.
(stan) ****/comp_modelsfile$
cat LOO.txt
LOO-CV
Normal
-5837.142
LogNormal
-5685.951
3P LogNormal
-5686.795
(stan) ****/comp_modelsfile$
３つのモデルに対するLOO-CVが出力されている。
データdata_3plognormal.txtについて考える。ファイルは、comp_modelsfile.zipにまとめた。
データdata_3plognormal.txtのヒストグラムは、図５のようである。
図５
正規分布モデル（リスト１）、対数正規分布モデル（リスト２）、３パラメータ対数正規分布モデル（リスト３）のStanスクリプトを用意して、PSIS-LOO-CVを求めるPythonスクリプトを作成した（リスト４）。
リスト４のスクリプトを実行すると、入力データファイル名の設定が求められるので、図５のファイル名を設定する。
(stan) ****/comp_modelsfile$
python comp_models.py
Input File Name(*.txt) =
data_3plognormal.txt
ファイル名を設定すると、計算が始まる。
計算が終了すると、図６のグラフが表示される。
図６
データのヒストグラムと、３つの確率密度関数のグラフである。ヒストグラムは、図５のものとはビンが異なるので形状が少し異なっている。
３つの確率密度関数は、パラメータの事後分布の中央値を推定値として描かれたものである。中央値を事後分布の代表値とすることは、Gelmanら（2014）に従ったものである。
各分布のLOOの値を見ると、３パラメータ正規分布が最も大きい。グラフを見ると、３パラメータ対数正規分布（水色点線）がよくデータの分布に合っていることが分かる。これに対して、正規分布と対数正規分布のグラフは、分布のボリュームが右にずれている様に思われる。
図６のWindowを閉じると、スクリプトの実行終了となり、以下のようなメッセージが端末に出力される。
LOO.txt was saved.
(stan) ****/comp_modelsfile$
「ファイルが保存された」とあるので、内容を確認してみる。
LOO.txt was saved.
(stan) ****/comp_modelsfile$
cat LOO.txt
LOO-CV
Normal
-4605.273
LogNormal
-4585.819
3P LogNormal
-4503.611
(stan) ****/comp_modelsfile$
３つのモデルに対するLOO-CVが出力されている。
Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., & Rubin, D. B. (2014). Bayesian data analysis, 3rd. Ed. CRC Press.
Lambert, B. (2018). A student’s guide to Bayesian statistics. SAGE.
Vehtari,A., Gelman,A., & Gabry,J. (2017). Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. Statistics and Computing, 27,1413-1432.
付録
リスト１ 正規分布のStanスクリプト（normalmodel.stan）
data {
int n;
array[n] real x;
}
transformed data {
real x_min;
real x_max;
x_min
= min(x);
x_max
= max(x);
}
parameters {
real mu;
real<lower=0.001> sgm;
}
model {
mu ~ uniform(x_min, x_max);
sgm
~ uniform(0.001, x_max-x_min);
x ~ normal(mu, sgm);
}
generated quantities {
array[n] real log_lik;
// log likelihood for
LOO
for (i
in 1:n) {
log_lik[i] = normal_lpdf(x[i] | mu, sgm);
}
}
リスト２ 対数正規分布のStanスクリプト（lognormalmodel.stan）
data {
int n;
array[n] real x;
}
transformed data {
real log_max_x;
real log_min_x;
log_max_x
= log(max(x));
log_min_x
= log(min(x));
}
parameters {
real<lower=log_min_x, upper=log_max_x>
mu;
real<lower=0.0001, upper=log_max_x - log_min_x> sgm;
}
model {
mu ~ uniform(log_min_x, log_max_x);
sgm
~ uniform(0.0001, log_max_x - log_min_x);
x ~ lognormal(mu, sgm);
}
generated quantities {
real mn;
real sd;
real mode;
array[n] real log_lik; // log likelihood for LOO
mn
= exp(mu + 0.5*(sgm^2));
mode = exp(mu - (sgm^2));
sd
= (exp(2*mu + (sgm^2)) * (exp(sgm^2) - 1.0))^0.5;
for (i
in 1:n) {
log_lik[i] = lognormal_lpdf(x[i] | mu, sgm);
}
}
リスト３ ３パラメータ対数正規分布のStanスクリプト（lognormal3pmodel.stan）
data {
int n;
array[n] real x;
}
transformed data {
real max_t;
max_t
= min(x);
}
parameters {
real<lower=0, upper=max_t> t;
real mu;
real<lower=0.0001, upper=1000>
sgm;
}
model {
t ~ uniform(0, max_t);
mu ~ normal(0, 1000);
sgm
~ uniform(0.0001, 1000);
for (i
in 1:n) {
x[i] - t ~ lognormal(mu, sgm);
}
}
generated quantities {
real mn;
real sd;
real mode;
array[n] real log_lik;
mn
= t + exp(mu + 0.5*(sgm^2));
mode = t + exp(mu -
(sgm^2));
sd
= (exp(2*mu + (sgm^2)) * (exp(sgm^2) - 1.0))^0.5;
for (i
in 1:n) {
log_lik[i] = lognormal_lpdf(x[i]-t | mu, sgm);
}
}
import matplotlib.pyplot
as plt
import numpy
as np
from cmdstanpy
import CmdStanModel
import arviz
as az
import pandas as pd
import scipy.stats
as ss
"""
対数正規分布確率密度関数
"""
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
"""
３パラメータ対数正規分布確率密度関数
"""
def lognormal_shift_pdf(x,
mu, sgm, t):
if x > t:
v1 =
sgm * (x - t) * ((2*np.pi)**0.5)
v2 =
np.exp(-((np.log(x - t) - mu)**2) / (2 * (sgm**2)))
return v2 / v1
else:
return 0.0
infl_nm = input('Input File Name(*.txt) = ')
x = []
with open(infl_nm,
'r') as f:
for v in f:
x.append(float(v))
print(x[:5], x[-5:])
def init_f():
return {'mu':np.mean(x), 'sgm':np.std(x)}
sm_n = CmdStanModel(stan_file="normalmodel.stan") # コンパイル
fit_n = sm_n.sample(data={'n':len(x), 'x':x},
chains=4, inits=[init_f()]*4)
print(fit_n.summary())
# MCMCサンプリングの統計
pd_fit_n = fit_n.draws_pd()
mu_med_n = np.median(pd_fit_n['mu'])
sgm_med_n = np.median(pd_fit_n['sgm'])
print('mu_med_n
', mu_med_n)
print('sgm_med_n
=', sgm_med_n)
i_data_n = az.from_cmdstanpy(fit_n, log_likelihood=True) # Arviz(az)用のデータに変換
#print('LOO-CV')
results = az.loo(i_data_n)
print(results)
print(results.elpd)
loos = [results.elpd]
sm_ln = CmdStanModel(stan_file="lognormalmodel.stan") # コンパイル
fit_ln = sm_ln.sample(data={'n':len(x), 'x':x})
print(fit_ln.summary())
# MCMCサンプリングの統計
pd_fit_ln = fit_ln.draws_pd()
mu_med_ln = np.median(pd_fit_ln['mu'])
sgm_med_ln = np.median(pd_fit_ln['sgm'])
i_data_ln = az.from_cmdstanpy(fit_ln, log_likelihood=True)
# Arviz(az)用のデータに変換
#print('LOO-CV')
results_ln = az.loo(i_data_ln)
print(results_ln)
print(results_ln.elpd)
loos.append(results_ln.elpd)
sm_3p = CmdStanModel(stan_file="lognormal3pmodel.stan") # コンパイル
fit_3p = sm_3p.sample(data={'n':len(x), 'x':x}, adapt_delta=0.999,
max_treedepth=17)
print(fit_3p.summary())
# MCMCサンプリングの統計
pd_fit_3p = fit_3p.draws_pd()
mu_med_3p = np.median(pd_fit_3p['mu'])
sgm_med_3p = np.median(pd_fit_3p['sgm'])
t_med_3p = np.median(pd_fit_3p['t'])
i_data_3p = az.from_cmdstanpy(fit_3p,
log_likelihood=True)
# Arviz(az)用のデータに変換
#print('LOO-CV')
results_3p = az.loo(i_data_3p)
print(results_3p)
print(results_3p.elpd)
loos.append(results_3p.elpd)
print('loo:\n', loos)
m_labels = ['Normal', 'LogNormal',
'3P LogNormal']
print(f"\n{'
':20s}{'LOO-CV':>20s}")
for i in
range(3):
print(f"{m_labels[i]:>20s}{loos[i]:>20.3f}")
outf_nm = 'LOO.txt'
with open(outf_nm,
'w') as f:
f.write(f"{'
':20s}{'LOO-CV':>20s}\n")
for i
in range(3):
f.write(f"{m_labels[i]:>20s}{loos[i]:>20.3f}\n")
print(f'\n{outf_nm}
was saved.\n')
plt.figure(figsize=(12,5))
min_x = np.min(x)
max_x = np.max(x)
xcoord = np.linspace(min_x, max_x, 1000)
Bins = np.linspace(min_x, max_x, 20)
plt.hist(x, bins=Bins, density=True,
alpha=0.7)
y_normal = ss.Normal(mu=mu_med_n, sigma= sgm_med_n).pdf(xcoord)
y_lognormal = lognormal_pdf(xcoord, mu_med_ln, sgm_med_ln)
y_lognormal_3p = []
for vx in xcoord:
y_lognormal_3p.append(lognormal_shift_pdf(vx,
mu_med_3p, sgm_med_3p,
t_med_3p))
plt.plot(xcoord, y_normal, lw=2, ls='-.', c='m',
label=f"Normal: LOO={loos[0]:.2f}")
plt.plot(xcoord, y_lognormal, lw=3, ls='--',
c='g',
label=f'Lognormal: LOO={loos[1]:.2f}')
plt.plot(xcoord,
y_lognormal_3p, lw=4, ls=':', c='aqua',
label=f'3P Lognormal: LOO={loos[2]:.2f}')
plt.legend(fontsize=12)
plt.tight_layout()
plt.show()