STAN Demo

STAN Overview

• Before class on Friday, either make sure your computer is equipped to run stan https://github.com/stan-dev/rstan/wiki/RStan-Getting-Started or that you have access to the R studio cloud https://andrewgelman.com/2018/10/12/stan-on-the-web-for-free-thanks-to-rstudio/.

1. Normal Data

Simulate Normal Data

## Simulate Data
set.seed(11122017)
N = 5000
mu = 0
sigmasq = 10
y = rnorm(N, mu, sqrt(sigmasq))

## set up priors
mu.0 <- 0
tausq.0 <- 100
nu.0 <- .01
sigmasq.0 <- 1
rate.param <- nu.0 * sigmasq.0 / 2
shape.param <- nu.0 / 2

Gibbs Sampler

library(invgamma)
### initialize vectors and set starting values 
num.sims <- 10000
mu.samples <- rep(0, num.sims)
sigmasq.samples <- rep(1, num.sims)

mean.y <- mean(y)
nu.n <- nu.0 + N

for (iter in 2:num.sims){
  # sample theta from full conditional
  mu.n <- (mu.0 / tausq.0 + N * mean.y / sigmasq.samples[iter - 1]) / (1 / tausq.0 + N / sigmasq.samples[iter - 1] )
  tausq.n <- 1 / (1/tausq.0 + N / sigmasq.samples[iter - 1])
  mu.samples[iter - 1] <- rnorm(1,mu.n,sqrt(tausq.n))
  
  # sample (1/sigma.sq) from full conditional
  sigmasq.n.theta <- 1/nu.n*(nu.0*sigmasq.0 + sum((y - mu.samples[iter])^2))
  sigmasq.samples[iter] <- rinvgamma(1,shape = nu.n/2, rate = nu.n*sigmasq.n.theta/2)
}


# plot marginal posterior of mu
hist(mu.samples,xlab=expression(mu),main=expression('Marginal Posterior of ' ~ mu),probability=T)
abline(v=mu,col='red',lwd=2)
# plot marginal posterior of sigmasq
hist(sigmasq.samples,xlab=expression(sigma[2]),main=expression('Marginal Posterior of ' ~ sigma[2]),probability=T)
abline(v=sigmasq,col='red',lwd=2)

# plot trace plots
plot(mu.samples,type='l',ylab=expression(mu), main=expression('Trace plot for ' ~ mu))
abline(h=mu,lwd=2,col='red')
plot(sigmasq.samples,type='l',ylab=expression(sigma[2]), main=expression('Trace plot for ' ~ sigma[2]))
abline(h=sigmasq,lwd=2,col='red')

# compute posterior mean and quantiles
mean(mu.samples)
quantile(mu.samples, probs = c(.025, .975))

mean(sigmasq.samples)
quantile(sigmasq.samples, probs = c(.025, .975))

Metropolis-Hastings Sampler

# Run Metropolis Algorithm
num.mcmc <- 10000
stepsize.mu <- .1
stepsize.sigmasq <- .75
acceptratio.mu <- acceptratio.sigmasq <- rep(0,num.mcmc)
mu.samplesmh <- rep(0,num.mcmc)
sigmasq.samplesmh <- sigmasq.star <- rep(3, num.mcmc)

for (iter in 2:num.mcmc){
  # mu
  mu.star <- mu.samplesmh[iter] + rnorm(1,0,stepsize.mu)

  log.p.current <- sum(dnorm(y, mean = mu.samplesmh[iter - 1],
                             sd = sqrt(sigmasq.samplesmh[iter - 1]), log=T)) +
  dnorm(mu.samplesmh[iter - 1],mean = mu.0, sd = sqrt(tausq.0), log=T)

  log.p.star <- sum(dnorm(y, mean = mu.star, sd = sqrt(sigmasq.samplesmh[iter - 1]), log=T)) +
    dnorm(mu.star, mean = mu.0, sd = sqrt(tausq.0), log=T)

  log.r <- log.p.star - log.p.current

  if (log(runif(1)) < log.r){
    mu.samplesmh[iter] <- mu.star
    acceptratio.mu[iter] <- 1
  } else{
    mu.samplesmh[iter] <- mu.samplesmh[iter - 1]
  }

  # sigma
  sigmasq.star[iter] <- sigmasq.samplesmh[iter-1] + rnorm(1,0,stepsize.sigmasq)
   
  log.p.current <- sum(dnorm(y, mean = mu.samplesmh[iter], 
                             sd = sqrt(sigmasq.samplesmh[iter - 1]), log=T)) + 
     dinvgamma(sigmasq.samplesmh[iter - 1], rate = rate.param, shape = shape.param, log=T)
   
  log.p.star <- sum(dnorm(y, mean = mu.samplesmh[iter], sd = sqrt(sigmasq.star[iter]), log=T)) +
    dinvgamma(sigmasq.star[iter], rate = rate.param, shape = shape.param, log=T)

  log.r <- log.p.star - log.p.current
  if (log(runif(1)) < log.r){
    sigmasq.samplesmh[iter] <- sigmasq.star[iter]
    acceptratio.sigmasq[iter] <- 1
  } else{
    sigmasq.samplesmh[iter] <- sigmasq.samplesmh[iter - 1]
  }
}

mean(acceptratio.mu)
mean(acceptratio.sigmasq)


plot(mu.samplesmh,type='l')
abline(h=mu,lwd=2,col='red')

plot(sigmasq.samplesmh,type='l')
abline(h=sigmasq,lwd=2,col='red')

# compute posterior mean and quantiles
mean(mu.samplesmh)
quantile(mu.samplesmh, probs = c(.025, .975))

mean(sigmasq.samplesmh)
quantile(sigmasq.samplesmh, probs = c(.025, .975))

JAGS

library(rjags)
# Define the model:
modelString = "
model {
    for (i in 1:N) {
        y[i] ~ dnorm(mu, tau.sq)
    }
    mu ~ dnorm(0, 1/ 100)
    tau.sq  ~ dgamma(.005, .005)
  sigmasq <- 1 / tau.sq
}
" 
writeLines( modelString , con="TEMPmodel.txt" )

jags <- jags.model('TEMPmodel.txt',
                   data = list('y' = y,
                               'N' = N),
                   n.chains = 4,
                   n.adapt = 100)
 


codaSamples = coda.samples( jags , variable.names=c("mu","tau.sq", 'sigmasq') ,
                            n.iter=1000)

summary(codaSamples)
plot(codaSamples)

Stan

library(rstan)
set.seed(11122017)
N = 5000
mu = 0
sigmasq = 10
y = rnorm(N, mu, sqrt(sigmasq))
num.mcmc <- 10000

stancode <- 'data {
  int<lower=0> N; 
  vector[N] y;
}
parameters {
  real<lower=0> sigmasq;
  real mu; 
} 
model {
  mu ~ normal(0, sqrt(100));
  sigmasq ~ inv_gamma(.005, .005)
  y ~ normal(mu, sqrt(sigmasq));
}' 

stan_data = list(N=N, y=y) # data passed to stan 
# set up the model
stan_model = stan(model_code = stancode, data = stan_data, chains = 0)
stanfit = stan(fit = stan_model, data = stan_data,
               iter=num.mcmc) # run the model print(stanfit,digits=2)
print(stanfit)

2. Logistic Regression

Select another model type, such as logistic regression, and work on samplers using Metropolis-Hastings, JAGS, and/or stan.

Metropolis-Hastings Sampler

JAGS

Stan