data {
int<lower=0> N; // N is a non-negative integer
int y[N]; // y is a length-N vector of integers
}
parameters {
real<lower=0,upper=1> theta; // theta is between 0 and 1
}
model {
theta ~ beta (1,1);
y ~ bernoulli(theta);
}
The chunk header for the chunk above looks like this:
{stan output.var = "simple_stan", cache = TRUE}
* `stan`: This is a Stan model, not R code.
* `output.var`: Name of R object for storing the stan DSO.
* `cache = TRUE`: Don't recompile this unless something changes.
Identify as many differences as you can between this Stan code and the corresponding JAGS code.
What is this model? Draw a model diagram for this model. What sort of data must you provide?
simple_stanfit <-
sampling(
simple_stan,
data = list(
N = 50,
y = c(rep(1, 15), rep(0, 35))
),
chains = 3, # default is 4
iter = 1000, # default is 2000
warmup = 200 # default is half of iter
)
This looks pretty similar to jags()
, but it isn’t identical. List as many differences as you can.
Running sampling()
produces a lot of output (hidden here with results = "hide"
in the R chunk header). Here’s just a portion of what is produced.
SAMPLING FOR MODEL 'stan-126b674f2d49f' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 1.7e-05 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.17 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1: Elapsed Time: 0.006157 seconds (Warm-up)
Chain 1: 0.021385 seconds (Sampling)
Chain 1: 0.027542 seconds (Total)
simple_stanfit
## Inference for Stan model: 83dbe9f99dbf55ff04494fdddf566a3d.
## 3 chains, each with iter=1000; warmup=200; thin=1;
## post-warmup draws per chain=800, total post-warmup draws=2400.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## theta 0.31 0.00 0.06 0.19 0.26 0.31 0.35 0.44 978 1
## lp__ -32.63 0.02 0.75 -34.74 -32.81 -32.35 -32.15 -32.10 1091 1
##
## Samples were drawn using NUTS(diag_e) at Fri Mar 22 11:53:55 2019.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).
Use CalvinBayes::posterior()
to create a dataframe with posterior samples. These can be plotted or explored using ggformula
or other familiar tools.
Use as.matrix()
or as.mcmc.list()
to create an object that can be used with bayesplot
just as we did when we used JAGS.
gf_dens(~theta, data = posterior(simple_stanfit))
mcmc_areas(as.mcmc.list(simple_stanfit), prob = 0.9, pars = "theta")
mcmc_trace(as.matrix(simple_stanfit))