STAT 812: Computational Statistics

STAN Sampling

library (rstan)
code <- '
  data {
    int<lower=0> T;
    vector[T] y;
  }

  parameters { 
    real mu;
    real<lower=0,upper=1> phis;
    real<lower=0.01> tausq;
    real h0;
    vector[T] h;
  }

  transformed parameters {
    real phi;
    real<lower=0> tau;
    phi <- 2*phis - 1;
    tau <- sqrt(tausq);
  }

  model {
    
    mu ~ normal(-10,5);
    phis ~ beta(20, 1.5);
    tausq ~ inv_gamma(2.5, 0.025);
    h0 ~ normal(mu, tau);
    h[1] ~ normal(mu + phi*(h0-mu), tau);
    for (t in 2:T)
      h[t] ~ normal(mu + phi*(h[t-1]-mu), tau);
    for (t in 1:T)
      y[t] ~ normal(0,exp(h[t]/2)); 
  }

  generated quantities {
    vector[T] log_lik;
    for (t in 1:T){
      log_lik[t] <- normal_log (y[t], 0 ,exp (h[t]/2)  );
    }
  }
'

N <- 500
mu <- 3
phi <- 0.3
tau <- 1
h0 <- 2
y<- h <- rep (0, N)
h[1] <- rnorm (1, mu + phi * (h0 - mu), tau)
for (i in 2:N) h[i] <- rnorm (1,mu + phi* (h[i-1] - mu), tau)
for (t in 1:N) y[t] <-  rnorm (1, 0, exp (h[t]/2) )

plot(y)

## fit stan model
fit <- stan(model_code = code, 
            data = list(y = y, T = N), 
            iter = 5000, chains = 2)

## DIAGNOSTIC(S) FROM PARSER:
## Info (non-fatal): assignment operator <- deprecated in the Stan language; use = instead.
## Info (non-fatal): assignment operator <- deprecated in the Stan language; use = instead.
## Info (non-fatal): assignment operator <- deprecated in the Stan language; use = instead.
## Info: Deprecated function 'normal_log'; please replace suffix '_log' with '_lpdf' for density functions or '_lpmf' for mass functions
## 
## 
## SAMPLING FOR MODEL '020e0bda59b3614430b9105cefc2bbd4' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0.000118 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.18 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 5000 [  0%]  (Warmup)
## Chain 1: Iteration:  500 / 5000 [ 10%]  (Warmup)
## Chain 1: Iteration: 1000 / 5000 [ 20%]  (Warmup)
## Chain 1: Iteration: 1500 / 5000 [ 30%]  (Warmup)
## Chain 1: Iteration: 2000 / 5000 [ 40%]  (Warmup)
## Chain 1: Iteration: 2500 / 5000 [ 50%]  (Warmup)
## Chain 1: Iteration: 2501 / 5000 [ 50%]  (Sampling)
## Chain 1: Iteration: 3000 / 5000 [ 60%]  (Sampling)
## Chain 1: Iteration: 3500 / 5000 [ 70%]  (Sampling)
## Chain 1: Iteration: 4000 / 5000 [ 80%]  (Sampling)
## Chain 1: Iteration: 4500 / 5000 [ 90%]  (Sampling)
## Chain 1: Iteration: 5000 / 5000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 4.97468 seconds (Warm-up)
## Chain 1:                4.63997 seconds (Sampling)
## Chain 1:                9.61465 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL '020e0bda59b3614430b9105cefc2bbd4' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 6.3e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.63 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 5000 [  0%]  (Warmup)
## Chain 2: Iteration:  500 / 5000 [ 10%]  (Warmup)
## Chain 2: Iteration: 1000 / 5000 [ 20%]  (Warmup)
## Chain 2: Iteration: 1500 / 5000 [ 30%]  (Warmup)
## Chain 2: Iteration: 2000 / 5000 [ 40%]  (Warmup)
## Chain 2: Iteration: 2500 / 5000 [ 50%]  (Warmup)
## Chain 2: Iteration: 2501 / 5000 [ 50%]  (Sampling)
## Chain 2: Iteration: 3000 / 5000 [ 60%]  (Sampling)
## Chain 2: Iteration: 3500 / 5000 [ 70%]  (Sampling)
## Chain 2: Iteration: 4000 / 5000 [ 80%]  (Sampling)
## Chain 2: Iteration: 4500 / 5000 [ 90%]  (Sampling)
## Chain 2: Iteration: 5000 / 5000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 5.93471 seconds (Warm-up)
## Chain 2:                5.7974 seconds (Sampling)
## Chain 2:                11.7321 seconds (Total)
## Chain 2:

## Warning: There were 2 chains where the estimated Bayesian Fraction of Missing Information was low. See
## http://mc-stan.org/misc/warnings.html#bfmi-low

## Warning: Examine the pairs() plot to diagnose sampling problems

## visualize fitting results
pars <- c("mu", "phi", "tau", "log_lik[1]", "h[1]")
plot(fit, pars = pars, plotfun="stan_plot")

## ci_level: 0.8 (80% intervals)

## outer_level: 0.95 (95% intervals)

plot(fit, pars = pars, plotfun="stan_trace")

## extract MCMC from stan objects
extract (fit, pars = c("mu", "phis", "tau")) -> fit.list

str(fit.list)

## List of 3
##  $ mu  : num [1:5000(1d)] 2.98 3.12 3.2 2.96 3.11 ...
##   ..- attr(*, "dimnames")=List of 1
##   .. ..$ iterations: NULL
##  $ phis: num [1:5000(1d)] 0.833 0.787 0.815 0.765 0.751 ...
##   ..- attr(*, "dimnames")=List of 1
##   .. ..$ iterations: NULL
##  $ tau : num [1:5000(1d)] 0.855 0.655 0.844 0.705 0.647 ...
##   ..- attr(*, "dimnames")=List of 1
##   .. ..$ iterations: NULL