Python Script for WAIC

A Simple Example of Python Script for WAIC

WAIC is a statistic, which assesses model fit. WAIC can be calculated by the following formula (Gelman et al., 2014, p. 174).

is given by the following formula (Gelman et al. 2014, p. 169).

Gelman et al., (2014, p. 173) explain two formulas for , and recommend the following one (p. 174).

A simple Python script to calculate WAIC in Listing 1 was developed, using the example in Gelman et al. 2014, p. 177. The data in the script was obtained by deleting the last row from the data file downloaded from the website explained in Gelman et al. 2021, p. XIII).

WAIC is calculated by the following code.

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))

WAIC = -2 * (lppd - pWAIC2)

When you run the script in Listing 1, you will get ouput like this:

Inference for Stan model: anon_model_105ce2197a3e94689c866d927437a11f.

4 chains, each with iter=2000; warmup=1000; thin=1;

post-warmup draws per chain=1000, total post-warmup draws=4000.

mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat

a 45.88 0.05 1.89 42.22 44.68 45.85 47.06 49.78 1433 1.0

b 3.15 0.02 0.79 1.54 2.67 3.16 3.67 4.74 1394 1.0

log_sgm 1.38 4.9e-3 0.2 1.02 1.24 1.36 1.5 1.81 1714 1.0

sgm 4.05 0.02 0.85 2.76 3.44 3.91 4.5 6.1 1663 1.0

lp__ -28.18 0.04 1.37 -31.79 -28.8 -27.8 -27.18 -26.63 1153 1.0

Samples were drawn using NUTS at Tue Sep 14 10:04:00 2021.

For each parameter, n_eff is a crude measure of effective sample size,

and Rhat is the potential scale reduction factor on split chains (at

convergence, Rhat=1).

WAIC = 87.06883542066456

Actual output depends on the sample by MCMC sampling.

References

Gelman, A., Carlin, J. B., Stern, H. S., Dunson, J. B., Vehtari, A., & Rubin, D. B. (2014). Bayesian Data Analysis, 3^rd ed. CRC Press.

Gelman, A., Hill, J., & Vehtari, A. (2021). Regression and Other Stories. Cambridge University Press.

Listing 1.

"""

Yasuharu Okamoto

"""

import numpy as np

import pystan

import scipy.stats as ss

# Simple Regression Model. See Gelman et al., 2014, p. 177

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);

}

"""

# Data. See Gelman et al., 2014. p. 165

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)