wfolta

library (coda)

gelman.diag (mcmc.list (mcmc (extract (stan.model, pars="alpha", permuted=F)[,1,1]),
                             mcmc (extract (stan.model, pars="alpha", permuted=F)[,2,1]),
                             mcmc (extract (stan.model, pars="alpha", permuted=F)[,3,1]),
                             mcmc (extract (stan.model, pars="alpha", permuted=F)[,4,1])))

heidel.diag (mcmc.list (mcmc (extract (stan.model, pars="alpha", permuted=F)[,1,1]),
                        mcmc (extract (stan.model, pars="alpha", permuted=F)[,2,1]),
                        mcmc (extract (stan.model, pars="alpha", permuted=F)[,3,1]),

wfolta

stan.model.s <- summary (stan.model)

a1 <- rnorm (3000, stan.model.s$summary["alpha[1]","mean"], stan.model.s$summary["alpha[1]","sd"])   # low-cooler intercepts
b1 <- rnorm (3000, stan.model.s$summary["beta[1]", "mean"], stan.model.s$summary["beta[1]", "sd"])   # low-cooler slopes
s1 <- rnorm (3000, stan.model.s$summary["sigma[1]","mean"], stan.model.s$summary["sigma[1]","sd"])   # low-cooler noise
g1 <- rnorm (3000, stan.model.s$summary["gamma[5]","mean"], stan.model.s$summary["gamma[5]","sd"])   # low-cooler regime R5

ci1 <- foreach (x=30:55, .combine="cbind") %do% c(x, quantile (a1 + s1 + g1 + b1 * (x - 55), c(0.025, 0.5, 0.975)))

a2 <- rnorm (3000, stan.model.s$summary["alpha[2]","mean"], stan.model.s$summary["alpha[2]","sd"])   # low-warmer intercepts

wfolta

Inference for Stan model: Electricity 55.
4 chains, each with iter=3e+05; warmup=150000; thin=1; 
post-warmup draws per chain=150000, total post-warmup draws=6e+05.

              mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
alpha[1]     446.4     1.1 31.7  385.0  427.1  445.6  464.9  509.7   810    1
alpha[2]     486.5     0.3 31.3  425.6  467.5  486.3  504.3  550.5  9341    1
alpha[3]     588.3     0.3 41.3  507.6  562.0  588.2  614.1  671.4 23192    1
alpha_alpha  513.7     0.4 59.9  400.6  477.7  510.0  547.8  644.4 22189    1
alpha_sigma   95.1     1.1 64.5   33.7   57.0   78.8  111.9  255.4  3251    1

wfolta

stan.code <-
"data
   {
   int cut ;
   int nu ;
   int N ;
   real dollars[N] ;
   real dca[N] ;
   int rate[N] ;
   int regime[N] ;