waic example in manual

[bob-carpenter]

Write up Andrew's WAIC example for the manual. See 12/28/13 e-mail to stan-dev and stan-users for the examples.

waic_example.R:

# R script for Waic example.  Also needed in this directory:  data file hibbs.dat and Stan file lm_waic.stan

# Little function to calculate posterior variances from simulation
colVars <- function (a){
  diff <- a - matrix (colMeans(a), nrow(a), ncol(a), byrow=TRUE)
  vars <- colMeans (diff^2)*nrow(a)/(nrow(a)-1)
  return (vars)
}

# The calculation of Waic!  Returns lppd, p_waic_1, p_waic_2, and waic, which we define
# as 2*(lppd - p_waic_2), as recommmended in BDA
waic <- function (stanfit){
  log_lik <- extract (stanfit, "log_lik")$log_lik
  lppd <- sum (log (colMeans(exp(log_lik))))
  p_waic_1 <- 2*sum (log(colMeans(exp(log_lik))) - colMeans(log_lik))
  p_waic_2 <- sum (colVars(log_lik))
  waic_2 <- -2*lppd + 2*p_waic_2
  return (list (waic=waic_2, p_waic=p_waic_2, lppd=lppd, p_waic_1=p_waic_1))
}

# Read in and prepare the data 
hibbs <- read.table ("hibbs.dat", header=TRUE)
year <- hibbs[,1]
growth <- hibbs[,2]
vote <- hibbs[,3]
inc <- hibbs[,4]
other <- hibbs[,5]
y <- vote
N <- length (y)
X <- cbind (rep(1,N), growth)
K <- ncol (X)

# Fit the model in Stan!
library ("rstan")
lm_waic <- stan (file="lm_waic.stan", data=c("N","K","X","y"), iter = 2000, chains = 4)
print (lm_waic)

# Calculate and print Waic
print (waic (lm_waic))

lm_waic.stan:

# This is the simple linear regression model from p.165 of Bayesian Data Analysis (3rd edition)
# with Waic as calculated on p.177.

# This code illustrates the calculation of Waic in Stan:  we define the log likelihood as a vector
# in the transformed parameters block, then sum it up in the model block.  The output for the N # # # components of the vector "loglik" are combined in the waic() function in R.

data {
  int N;
  int K;
  vector[N] y;
  matrix[N,K] X;
}
parameters {
  vector[K] b;
  real<lower=0> sigma;
}
transformed parameters {
  vector[N] log_lik;
  for (n in 1:N){
    log_lik[n] <- normal_log (y[n], X[n]*b, sigma);
  }
}
model {
  increment_log_prob (-log(sigma));    # log prior for p(sigma) proportional to 1/sigma
  increment_log_prob (sum (log_lik));   # log likelihood
}

hibbs.dat:

year growth inc.party.vote incumbent.party other
1952 2.4 44.6 "Stevenson" "Eisenhower"
1956 2.89 57.76 "Eisenhower" "Stevenson"
1960 .85 49.91 "Nixon" "Kennedy"
1964 4.21 61.34 "Johnson" "Goldwater"
1968 3.02 49.60 "Humphrey" "Nixon"
1972 3.62 61.79 "Nixon" "McGovern"
1976 1.08 48.95 "Ford" "Carter"
1980 -.39 44.70 "Carter" "Reagan"
1984 3.86 59.17 "Reagan" "Mondale"
1988 2.27 53.94 "Bush, Sr." "Dukakis"
1992 .38 46.55 "Bush, Sr." "Clinton"
1996 1.04 54.74 "Clinton" "Dole"
2000 2.36 50.27 "Gore" "Bush, Jr."
2004 1.72 51.24 "Bush, Jr." "Kerry"
2008 .1 46.32 "McCain" "Obama"

[andrewgelman]

Maybe it's not worth writing it up in the manual, if we're planning to incorporate it more fully into Stan and Rstan?
A

On Dec 29, 2013, at 8:56 PM, Bob Carpenter wrote:

Write up Andrew's WAIC example for the manual. See 12/28/13 e-mail to stan-dev and stan-users for the examples.

—
Reply to this email directly or view it on GitHub.

[bob-carpenter]

I think it should go in the manual until we get
it incorporated into Stan.

We're moving to put everything in Stan so that
the command-line, RStan and PyStan interfaces can
share implementations.

On 12/29/13, 4:13 PM, Andrew Gelman wrote:

Maybe it's not worth writing it up in the manual, if we're planning to incorporate it more fully into Stan and Rstan?
A

On Dec 29, 2013, at 8:56 PM, Bob Carpenter wrote:

Write up Andrew's WAIC example for the manual. See 12/28/13 e-mail to stan-dev and stan-users for the examples.

—
Reply to this email directly or view it on GitHub.

—
Reply to this email directly or view it on GitHub #473 (comment).

[rtrangucci]

repeat of #438 ? We're now moving towards PSIS-LOO, and away from WAIC.

[bob-carpenter]

Yes.