Simple Calculation of LOO-CV and WAIC in Python

Given a MCMC sample by PyStan, LOO-CV and WAIC can be calculated by ArviZ. A simple example is shown in Listing 1 below.

The script in Listing 1 is modeled for the 15 elections (Gelman et al. 2014, pp.165, 176-177). The data in the script was obtained from the data file downloaded from the website explained in Gelman et al. 2021, p. XIII).

LOO-CV can be obtained by the function loo as follows:

results_loo = az.loo(inference_data, scale='deviance')

WAIC can be obtained by the fuction waic as follows:

results_waic = az.waic(inference_data, scale='deviance')

The object inference_data in the above codes is an InferenceData object constructed from the MCMC sample as follows:

sm = stan.build(stan_code, data = Data)

fit = sm.sample(num_samples = 5000, num_warmup = 5000)

inference_data = az.from_pystan(posterior = fit, log_likelihood = {'LogL':'LL'},

posterior_model = sm)

The variable for the log likelihood, ‘LL’, in the above code, was generated in the MCMC sampling as follows:

generated quantities {

array[N] real LL; // Log-Likelihood

for (i in 1:N)

LL[i] = normal_lpdf(y[i] | a + b*x[i], sgm);

}

Log likelihood data given in the dictionary as below is used in calculation by the functions loo and waic.

log_likelihood = {'LogL':'LL'}

The script in Listing 1 was run in a Python 3.12 virtual environment of Anaconda on Ubuntu in WSL (Windows Subsystem for Linux) of Windows.

The following output was displayed on the terminal.

mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat

a 45.807 1.896 42.230 49.373 0.023 0.016 7119.0 8380.0 1.0

b 3.175 0.787 1.686 4.639 0.009 0.006 7422.0 8555.0 1.0

sgm 4.153 0.903 2.664 5.846 0.009 0.007 9384.0 10504.0 1.0

LOO-CV = 87.388

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(

WAIC = 87.014

Statistics

LOOCV:

Computed from 20000 posterior samples and 15 observations log-likelihood matrix.

Estimate SE

deviance_loo 87.39 6.63

p_loo 2.74 -

------

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%

WAIC:

Computed from 20000 posterior samples and 15 observations log-likelihood matrix.

Estimate SE

deviance_waic 87.01 6.37

p_waic 2.55 -

For comparison, check Gelman et al. (2014), p.177.

The script file loowaic.py shown in Listing 1 is archived in the file scriptfile.zip , which can be freely downloaded.

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.

import numpy as np

import stan

import arviz as az

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

stan_code = """

data {

int N;

array[N] real y;

array[N] real x;

}

parameters {

real p_a;

real p_b;

real<lower=0.0> p_sgm;

}

transformed parameters {

real a;

real b;

real sgm;

a = p_a * 5.0;

b = p_b * 1.0;

sgm = p_sgm * 1.0 + 0.001;

}

model {

p_a ~ normal(0.0, 10.0);

p_b ~ normal(0.0, 10.0);

p_sgm ~ exponential(0.1);

for (i in 1:N) {

y[i] ~ normal(a + b*x[i], sgm);

}

generated quantities {

array[N] real LL; // Log-Likelihood

for (i in 1:N)

LL[i] = normal_lpdf(y[i] | a + b*x[i], sgm);

}

"""

# Data. See Gelman et al., 2021.

# from 15 ElectionsEconomy data (with the last record deleted)

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(num_samples = 5000, num_warmup = 5000)

inference_data = az.from_pystan(posterior = fit, log_likelihood = {'LogL':'LL'},