WAICとPSIS-LOO-CVの計算例
確率モデルのデータへの適合度の指標として、ベイズ分析ではWAIC（Watanabe-Akaike or widely available information criterion）とLooCV（Leave-one-out cross-validation）などがある。
WAICは次式で与えられる（Gelman et al., 2014, p. 174）。
ここで、
は次式で与えられる（Gelman et al. 2014, p. 169）。
WAICは２通りの計算法が示されているが、pWAIC2として次式が推奨されている（Gelman et al., 2014, p. 173）。
LooCVは、次式が説明されている（Gelman et al. 2014, pp. 175-177）。
ここで、leave-one-outにおける事後分布をPareto-smoothed importance samplingにより求めたときのLooCVとして次式が考えられる。
WAICとPSIS-LOO-CVの計算プログラムをリスト１のように書いてみた。データと分析モデルはGelman et al. 2014, p. 177で扱われているもので、データはGelman et al. 2021, p. XIIIに説明されているウェブサイトからダウンロードしたファイルから最後の１行のデータを削除して得ることができる。
WAICの計算部分は、以下の簡単なスクリプトで表されている。
lppd =
np.sum(np.log(np.mean(p_y_theta, axis = 1)))
pWAIC2 =
np.sum(np.var(np.log(p_y_theta), axis = 1))
print('\nWAIC =', -2 * (lppd -
pWAIC2))
ここで、p_y_thetaはStanの実行時に計算される尤度probsから以下のコードで設定されている。
data_f = fit.to_frame()
p_list = []
for i in range(len(y)):
p_list.append(data_f[f'probs.{i+1}'])
p_y_theta = np.array(p_list)
尤度probsは次のスクリプトで与えられている。
generated quantities {
array[N]
real probs; // for calculation of WAIC by
Formulas
array[N] real log_lik;
// for arviz.waic and
loo
for
(i in 1:N) {
log_lik[i] = normal_lpdf(y[i] | a + b*x[i], sgm);
probs[i] = exp(log_lik[i]);
}
}
リスト１を実行すると、WAICが算出されて以下のように出力される。
WAIC =
87.18113720366594
MCMCサンプリングに基づく計算であるので、サンプリングごとに値は多少変化する。
リスト１のスクリプトの実行は、Ubuntu 24/WSL/Windows 10においてPython 3.12/Anacondaの環境にPyStan 3.10.0をインストールして行った（2024.07.06）。Arvizは0.18.0を用いた。
WAICは、ライブラリArvizの関数arviz.waicによっても求めることができる。
Arviz.waicを利用するとき、対数尤度を設定した変数をStanの実行時にlog_likの名前で用意しておくと、Arviz.waicの呼び出しにおいて自動的にlog_likに対数尤度が設定されているとして処理される。リスト１では以下のようにlog_likを用意している（再掲）。
generated quantities {
array[N] real
probs; // for calculation of WAIC by
Formulas
array[N] real log_lik;
// for arviz.waic and
loo
for
(i in 1:N) {
log_lik[i] = normal_lpdf(y[i] | a + b*x[i], sgm);
probs[i] = exp(log_lik[i]);
}
}
変数probsは、直接WAICを計算する次のスクリプトで用いている（一部再掲）。
data_f = fit.to_frame()
p_list = []
for i in range(len(y)):
p_list.append(data_f[f'probs.{i+1}'])
p_y_theta = np.array(p_list)
#
#
Calculation of WAIC
#
lppd =
np.sum(np.log(np.mean(p_y_theta, axis = 1)))
pWAIC2 =
np.sum(np.var(np.log(p_y_theta), axis = 1))
print('\nWAIC =', -2 * (lppd -
pWAIC2))
Arviz.waicによりWAICを算出するときは、MCMCサンプルのデータをarviz.waicで扱えるようにarviz.from_pystanで変換してからarviz.waicを呼び出す。
i_data = az.from_pystan(posterior=fit,
posterior_model=sm)
vWAIC = az.waic(i_data,
scale='deviance')
print('\narviz.WAIC =',
vWAIC.elpd_waic)
上のスクリプトの実行により、以下の出力が得られる。
/home/yasuharu/anaconda3/envs/py312/lib/python3.12/site-packages/arviz/stats/stats.py:1632:
UserWarning: For one or more samples the posterior variance of the log
predictive densities exceeds 0.4. This could be indication of WAIC starting to
fail.
See http://arxiv.org/abs/1507.04544
for details
warnings.warn(
arviz.WAIC = 87.18113720366594
使用したarvizはarviz 0.18.0である。計算において、0.4を超える分散が得られたものがあるという警告メッセージも出力されている。
次のスクリプトで分散を出力してみる。
print('\nVariances =\n',
np.var(np.log(p_y_theta), axis = 1))
以下の出力を得る。
Variances =
[1.15644306 0.06719878 0.05396204
0.13037081 0.31930752 0.25458442
0.04780136 0.07868891 0.0706631 0.04486204 0.05653934 0.26759746
0.05697368 0.04567997 0.06178181]
0.4を超える値が含まれていることが確認できる。
次に、PSIS-LOO-CVを以下のスクリプトで求める。
vLOOCV = az.loo(i_data,
scale='deviance')
print('\nPSIS-LOO-CV =',
vLOOCV.elpd_loo)
MCMCサンプルのデータi_dataに含まれる対数尤度の変数log_likが自動的に認識されて計算が行われ、次の出力を得る。
PSIS-LOO-CV =
87.59311722621717
Arviz.waicおよびarviz.looの実行により得られた戻り値を、そのまま次のスクリプトで出力してみる。
print('\nWAIC:\n', vWAIC)
print('\nPSIS-LOO-CV:\n', vLOOCV)
以下の出力が得られた。
WAIC:
Computed from 4000 posterior samples and
15 observations log-likelihood matrix.
Estimate SE
deviance_waic 87.18 6.73
p_waic
2.71 -
There has been a warning during the
calculation. Please check the results.
PSIS-LOO-CV:
Computed from 4000 posterior samples and 15
observations log-likelihood matrix.
Estimate SE
deviance_loo 87.59 7.01
p_loo
2.92 -
------
Pareto k diagnostic values:
Count Pct.
(-Inf, 0.5] (good) 14 93.3%
(0.5, 0.7] (ok)
1 6.7%
(0.7, 1] (bad)
0 0.0%
(1, Inf) (very bad) 0 0.0%
なお、PyhtonにおけるStanについては、岡本（2019）で説明している。
リスト１ （スクリプトファイルcalcwaicfile.zip）
"""
Gelman et al. Bayesian Data Analysis, 3rd ed, 2014
"""
import numpy as np
import stan
import arviz as az
stan_code = """
data {
int
N;
array[N] real y;
array[N] real x;
}
parameters {
// A
noninformative prior distribution, BDA 3, p.165.
real
a;
real b;
real
log_sgm;
}
transformed parameters {
real
sgm;
sgm
= exp(log_sgm);
}
model {
for
(i in 1:N) {
y[i] ~ normal(a + b*x[i], sgm);
}
}
generated quantities {
array[N]
real probs; // for calculation of WAIC by
Formulas
array[N] real log_lik;
// for arviz.waic and
loo
for
(i in 1:N) {
log_lik[i] = normal_lpdf(y[i] | a + b*x[i], sgm);
probs[i] = exp(log_lik[i]);
}
}
"""
xy = [
[44.6, 2.4],
[57.76, 2.89],
[49.91, 0.85],
[61.34, 4.21],
[49.6, 3.02],
[61.79, 3.62],
[48.95, 1.08],
[44.7, -0.39],
[59.17, 3.86],
[53.94, 2.27],
[46.55, 0.38],
[54.74, 1.04],
[50.27, 2.36],
[51.24, 1.72],
[46.32, 0.1]]
x = np.array(xy).T[1]
y = np.array(xy).T[0]
Data = {'N':len(y), 'y':y, 'x':x}
sm = stan.build(stan_code,
data=Data)
fit = sm.sample()
data_f = fit.to_frame()
p_list = []
for i in range(len(y)):
p_list.append(data_f[f'probs.{i+1}'])
p_y_theta = np.array(p_list)
#
#
Calculation of WAIC
#
lppd =
np.sum(np.log(np.mean(p_y_theta, axis = 1)))
pWAIC2 =
np.sum(np.var(np.log(p_y_theta), axis = 1))
print('\nWAIC =', -2 * (lppd -
pWAIC2))
#
# For
arviz.waic and arviz.loo
#
i_data = az.from_pystan(posterior=fit,
posterior_model=sm)
vWAIC = az.waic(i_data,
scale='deviance')
print('\narviz.WAIC =',
vWAIC.elpd_waic)
print('\nVariances =\n',
np.var(np.log(p_y_theta), axis = 1, ddof = 1))
vLOOCV = az.loo(i_data,
scale='deviance')
print('\nPSIS-LOO-CV =', vLOOCV.elpd_loo)
print('\nWAIC:\n', vWAIC)
print('\nPSIS-LOO-CV:\n', vLOOCV)
Gelman, A., Carlin, J. B., Stern, H. S., Dunson, J. 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 Press.
岡本安晴（2019）いまさら聞けないPythonでデータ分析．丸善出版