WAICの計算例
確率モデルのデータへの適合度の指標の１つとしてWAICがある。WAICは次式で与えられる（Gelman et al., 2014, p. 174）。
ここで、
は次式で与えられる（Gelman et al. 2014, p. 169）。
WAICは２通りの計算法が示されているが、pWAIC2として次式が推奨されている（Gelman et al., 2014, p. 173）。
上の計算方法によるWAICの計算プログラムをリスト１のように書いてみた。データと分析モデルは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, ddof = 1))
print('\nWAIC =', -2 * (lppd -
pWAIC2))
リスト１を実行すると、WAICが算出されて以下のように出力される。
WAIC =
86.88745443786172
MCMCサンプリングに基づく計算であるので、サンプリングごとに値は多少変化することに注意。
PyhtonにおけるStan (PyStan) については、岡本（2019）で説明している。
import scipy.stats as ss
import pickle
stan_code = """
data {
int
N;
real
y[N];
real
x[N];
}
parameters {
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);
}
}
"""
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]
sm = pystan.StanModel(model_code =
stan_code)
Data = {'N':len(y), 'y':y, 'x':x}
fit = sm.sampling(data = Data, n_jobs
= 1)
print(fit)
n = len(y)
a = fit['a']
b = fit['b']
sgm = fit['sgm']
S = len(a)
p_y_theta = np.empty((n, S))
for i in range(n):
for s in range(S):
p_y_theta[i][s] = ss.norm.pdf(y[i], a[s] + b[s] * x[i], sgm[s])
lppd =
np.sum(np.log(np.mean(p_y_theta, axis = 1)))
pWAIC2 =
np.sum(np.var(np.log(p_y_theta), axis = 1, ddof = 1))
print('\nWAIC =', -2 * (lppd -
pWAIC2))
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でデータ分析．丸善出版