Watanabe Akaike information criterion (WAIC)

Waic examples from Vehtari & Gelman 2014.

Forecasting elections from the economy:

y = [44.6 57.76 49.91 61.34 49.60 61.79 48.95 44.70 59.17 53.94...
     46.55 54.74 50.27 51.24 46.32]';
growth = [2.4 2.89 .85 4.21 3.02 3.62 1.08 -.39 3.86 2.27 .38 ...
     1.04 2.36 1.72 .1]';
X = [ones(size(y)) , growth];

model = {
   'data {'
   '   int N;'
   '   int J;'
   '   vector[N] y;'
   '   matrix[N,J] X;'
   '}'
   'parameters {'
   '   vector[J] b;'
   '   real<lower=0> sigma;'
   '}'
   'model {'
   '   y ~ normal(X*b, sigma); // data model'
   '   increment_log_prob(-log(sigma)); // log prior for p(sigma) propto 1/sigma'
   '}'
   'generated quantities {'
   '   vector[N] log_lik;'
   '   for (n in 1:N){'
   '      log_lik[n] <- normal_log(y[n], X[n]*b, sigma);'
   '   }'
   '}'
};

data = struct('N',numel(y),'J',size(X,2),'X',X,'y',y);

hibbs_fit = stan('model_code',model,'data',data,'verbose',true);

log_lik = hibbs_fit.extract.log_lik;
mstan.waic(log_lik)

Scaled eight schools:

schools_code = {
   'data {'
   '   int J;'
   '   vector[J] y;vector<lower=0>[J] sigma;'
   '}'
   'parameters {'
   '   real mu;'
   '   real<lower=0> tau;'
   '   vector[J] eta;'
   '}'
   'transformed parameters {'
   '   vector[J] theta;'
   '   theta <- mu + tau*eta;'
   '}'
   'model {'
   '   eta ~ normal(0, 1);'
   '   y ~ normal(theta, sigma);'
   '}'
   'generated quantities {'
   '   vector[J] log_lik;'
   '   for (j in 1:J){'
   '      log_lik[j] <- normal_log(y[j], theta[j], sigma[j]);'
   '   }'
   '}'
};

schools_dat = struct('J',8,...
                     'y',[28 8 -3 7 -1 1 18 12],...
                     'sigma',[15 10 16 11 9 11 10 18]);

schools_fit = stan('model_code',schools_code,'data',schools_dat,'verbose',true);

log_lik = schools_fit.extract.log_lik;
mstan.waic(log_lik)

% Scale data
schools_dat = struct('J',8,...
                     'y',5*[28 8 -3 7 -1 1 18 12],...
                     'sigma',[15 10 16 11 9 11 10 18]);

schools_fit2 = stan('fit',schools_fit,'data',schools_dat,'verbose',true);

log_lik = schools_fit2.extract.log_lik;
mstan.waic(log_lik)