Adding levels of randomness

Peter Ralph

16 November 2020 – Advanced Biological Statistics

Hierarchical coins

Motivating problem: more coins

Suppose now we have data from \(n\) different coins from the same source. We don’t assume they have the same \(\theta\), but don’t know what its distribution is, so try to learn it.

\[\begin{aligned} Z_i &\sim \Binom(N_i, \theta_i) \\ \theta_i &\sim \Beta(\alpha, \beta) \\ \alpha &\sim \Unif(0, 100) \\ \beta &\sim \Unif(0, 100) \end{aligned}\]

Goal: find the posterior distribution of \(\alpha\), \(\beta\).

Problem: we don’t have a nice mathematical expression for this posterior distribution.

MC Stan

“Quick and easy” MCMC: Stan

The skeletal Stan program

data {
    // stuff you input
}
transformed data {
    // stuff that's calculated from the data (just once, at the start)
}
parameters {
    // stuff you want to learn the posterior distribution of
}
transformed parameters {
    // stuff that's calculated from the parameters (at every step)
}
model {
    // the action!
}
generated quantities {
    // stuff you want computed also along the way
}

How to do everything: see the user’s manual.

Beta-Binomial with Stan

From before:

We’ve flipped a coin 10 times and got 6 Heads. We think the coin is close to fair, so put a \(\Beta(20,20)\) prior on it’s probability of heads, and want the posterior distribution.

\[\begin{aligned} Z &\sim \Binom(10, \theta) \\ \theta &\sim \Beta(20, 20) \end{aligned}\] Sample from \[\theta \given Z = 6\]

\[\begin{aligned} Z &\sim \Binom(10, \theta) \\ \theta &\sim \Beta(20, 20) \end{aligned}\]

data {
    // stuff you input
}
parameters {
    // stuff you want to learn 
    // the posterior distribution of
}
model {
    // the action!
}

\[\begin{aligned} Z &\sim \Binom(10, \theta) \\ \theta &\sim \Beta(20, 20) \end{aligned}\]

data {
    int N;   // number of flips
    int Z;   // number of heads
}
parameters {
    // stuff you want to learn 
    // the posterior distribution of
}
model {
    // the action!
}

\[\begin{aligned} Z &\sim \Binom(10, \theta) \\ \theta &\sim \Beta(20, 20) \end{aligned}\]

data {
    int N;   // number of flips
    int Z;   // number of heads
}
parameters {
    // probability of heads
    real<lower=0,upper=1> theta;  
}
model {
    // the action!
}

\[\begin{aligned} Z &\sim \Binom(10, \theta) \\ \theta &\sim \Beta(20, 20) \end{aligned}\]

data {
    int N;   // number of flips
    int Z;   // number of heads
}
parameters {
    // probability of heads
    real<lower=0,upper=1> theta;
}
model {
    Z ~ binomial(N, theta);
    theta ~ beta(20, 20);
}

Compiling Stan model, in R

library(rstan)
stan_block <- "
data {
    int N;   // number of flips
    int Z;   // number of heads
}
parameters {
    // probability of heads
    real<lower=0,upper=1> theta;
}
model {
    Z ~ binomial(N, theta);
    theta ~ beta(20, 20);
}
"
bb_model <- stan_model(
               model_code=stan_block)

Optimization: maximum likelihood

With a uniform prior, the “maximum posterior” parameter values are also the maximum likelihood values.

> help(optimizing)

optimizing                package:rstan                R Documentation

Obtain a point estimate by maximizing the joint posterior

Description:

     Obtain a point estimate by maximizing the joint posterior from the
     model defined by class 'stanmodel'.

Usage:

     ## S4 method for signature 'stanmodel'
     optimizing(object, data = list(),

Maximum posterior

library(rstan)
(fit <- optimizing(bb_model,  # stan model from above
                   data=list(N=10, Z=6)))

## $par
##     theta 
## 0.5208333 
## 
## $value
## [1] -33.22939
## 
## $return_code
## [1] 0
## 
## $theta_tilde
##          theta
## [1,] 0.5208333

The answer!

The maximum a posteriori probability estimate (MAP) is 0.5208333.

post_fun <- function (p) {
    n <- 10; z <- 6; a <- b <- 20
    lik <- dbinom(z, size=n, prob=p)
    prior <- dbeta(p, a, b)
    return( prior * lik )
}
curve(post_fun, 0, 1, xlab=expression(theta), ylab='posterior prob')
points(fit$par, post_fun(fit$par), col='red', pch=20)

Note: Since the prior was uniform, then this is the maximum likelihood estimate (MLE) also.

plot of chunk r stan_check

Sampling from the posterior distribution

library(rstan)
fit <- sampling(bb_model,  # compiled stan model from above
                data=list(N=10, Z=6),
                chains=3, iter=10000)

lp__ is the log posterior density. Note n_eff and Rhat.

print(fit)

## Inference for Stan model: 5f3115ac9593cded393e67d2d0e3ce84.
## 3 chains, each with iter=10000; warmup=5000; thin=1; 
## post-warmup draws per chain=5000, total post-warmup draws=15000.
## 
##         mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
## theta   0.52    0.00 0.07   0.38   0.47   0.52   0.57   0.66  6136    1
## lp__  -35.12    0.01 0.71 -37.17 -35.28 -34.85 -34.67 -34.62  5608    1
## 
## Samples were drawn using NUTS(diag_e) at Mon Nov 15 21:07:21 2021.
## 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).

Fuzzy caterpillars are good.

stan_trace(fit)

plot of chunk r trace_rstan

Stan uses ggplot2.

stan_hist(fit, bins=20) + xlim(0,1)

## Warning: Removed 2 rows containing missing values (geom_bar).

plot of chunk r plot_rstan

The samples:

sampled_theta <- extract(fit, permuted=FALSE, inc_warmup=TRUE)

plot of chunk r the_samples

What’s the posterior probability that \(\theta < 0.5\)?

samples <- extract(fit)
mean(samples$theta < 0.5)

## [1] 0.384

# compare to analytic solution
pbeta(0.5, shape1=20+6, shape2=20+4)

## [1] 0.3877248

Hierarchical Coins

\[\begin{aligned} Z_i &\sim \Binom(N_i, \theta_i) \\ \theta_i &\sim \Beta(\alpha, \beta) \\ \alpha &\sim \Unif(0, 100) \\ \beta &\sim \Unif(0, 100) \end{aligned}\]

data {
    // the data (input)
}
parameters {
    // the parameters (output)
}
model {
    // how they are related
}

\[\begin{aligned} Z_i &\sim \Binom(N_i, \theta_i) \\ \theta_i &\sim \Beta(\alpha, \beta) \\ \alpha &\sim \Unif(0, 100) \\ \beta &\sim \Unif(0, 100) \end{aligned}\]

data {
    int n;      // number of coins
    int N[n];   // number of flips
    int Z[n];   // number of heads
}
parameters {
    // the parameters (output)
}
model {
    // how they are related
}

\[\begin{aligned} Z_i &\sim \Binom(N_i, \theta_i) \\ \theta_i &\sim \Beta(\alpha, \beta) \\ \alpha &\sim \Unif(0, 100) \\ \beta &\sim \Unif(0, 100) \end{aligned}\]

data {
    int n;      // number of coins
    int N[n];   // number of flips
    int Z[n];   // number of heads
}
parameters {
    // probability of heads
    real<lower=0,upper=1> theta[n];
    real<lower=0,upper=100> alpha;
    real<lower=0,upper=100> beta;
}
model {
    // how they are related
}

\[\begin{aligned} Z_i &\sim \Binom(N_i, \theta_i) \\ \theta_i &\sim \Beta(\alpha, \beta) \\ \alpha &\sim \Unif(0, 100) \\ \beta &\sim \Unif(0, 100) \end{aligned}\]

data {
    int n;      // number of coins
    int N[n];   // number of flips
    int Z[n];   // number of heads
}
parameters {
    // probability of heads
    real<lower=0,upper=1> theta[n];
    real<lower=0, upper=100> alpha;
    real<lower=0, upper=100> beta;
}
model {
    Z ~ binomial(N, theta);
    theta ~ beta(alpha, beta);
    // uniform priors "go without saying"
    // alpha ~ uniform(0, 100);
    // beta ~ uniform(0, 100);
}

Exercise

Data:

set.seed(23)
ncoins <- 100
true_theta <- rbeta(ncoins, 20, 50)
N <- rep(50, ncoins)
Z <- rbinom(ncoins, size=N, prob=true_theta)

Find the posterior distribution on alpha and beta: check convergence with print() and stan_trace(), then plot using stan_hist(), pairs(), and/or stan_scat() (with e.g., pars=c("alpha", "beta", "theta[1]", "theta[2]")).

data {
    int n;      // number of coins
    int N[n];   // number of flips
    int Z[n];   // number of heads
}
parameters {
    // probability of heads
    real<lower=0,upper=1> theta[n];
    real<lower=0,upper=100> alpha;
    real<lower=0,upper=100> beta;
}
model {
    Z ~ binomial(N, theta);
    theta ~ beta(alpha, beta);
    // uniform priors 'go without saying'
    // alpha ~ uniform(0, 100);
    // beta ~ uniform(0, 100);
}

IN CLASS

model_code <- "
data {
    int n;      // number of coins
    int N[n];   // number of flips
    int Z[n];   // number of heads
}
parameters {
    // probability of heads
    real<lower=0,upper=1> theta[n];
    real<lower=0,upper=100> alpha;
    real<lower=0,upper=100> beta;
}
model {
    Z ~ binomial(N, theta);
    theta ~ beta(alpha, beta);
    // uniform priors 'go without saying'
    // alpha ~ uniform(0, 100);
    // beta ~ uniform(0, 100);
}
"
model <- stan_model(model_code=model_code)

data:

set.seed(23)
ncoins <- 100
true_theta <- rbeta(ncoins, 20, 50)
N <- rep(50, ncoins)
Z <- rbinom(ncoins, size=N, prob=true_theta)

fit mmmodel

fit1 <- sampling(model,
         data=list(N=N, Z=Z, n=ncoins),
         chains=3, iter=1000)

## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## http://mc-stan.org/misc/warnings.html#bulk-ess

fit2 <- sampling(model,
         data=list(N=N, Z=Z, n=ncoins),
         chains=3, iter=5000)

plot trace

stan_trace(fit2, pars=c("alpha", "beta"))

plot of chunk r plottrace

look at results

Note the posterior intervals for alpha and beta:

rstan::summary(fit2)$summary

##                     mean      se_mean          sd          2.5%           25%           50%           75%         97.5%     n_eff      Rhat
## theta[1]       0.2917988 0.0004125683  0.03779479     0.2211084     0.2661861     0.2907644     0.3160459     0.3700863 8392.1310 0.9998365
## theta[2]       0.2581785 0.0004120545  0.03704161     0.1879117     0.2325430     0.2578913     0.2825079     0.3332164 8081.1009 1.0000029
## theta[3]       0.3065083 0.0004536602  0.03886474     0.2344475     0.2795073     0.3055743     0.3318121     0.3863964 7339.2248 1.0000731
## theta[4]       0.2851268 0.0004214931  0.03821356     0.2129591     0.2585586     0.2844941     0.3102555     0.3603522 8219.6637 1.0002799
## theta[5]       0.2571079 0.0004446583  0.03672410     0.1867113     0.2321856     0.2561545     0.2816435     0.3295946 6821.0259 0.9997497
## theta[6]       0.3270950 0.0004603358  0.04010376     0.2523229     0.2995214     0.3257214     0.3533930     0.4091024 7589.6286 0.9998880
## theta[7]       0.2992070 0.0004153123  0.03810265     0.2275179     0.2728440     0.2980592     0.3243058     0.3769895 8417.0662 0.9998509
## theta[8]       0.2922524 0.0003978970  0.03824006     0.2214836     0.2660296     0.2904321     0.3173091     0.3702462 9236.2518 0.9999761
## theta[9]       0.2639424 0.0004116463  0.03707327     0.1939021     0.2378778     0.2632794     0.2887504     0.3395545 8110.9801 1.0002949
## theta[10]      0.2713837 0.0004104860  0.03684887     0.1995005     0.2465886     0.2712412     0.2959662     0.3437765 8058.4545 0.9997515
## theta[11]      0.2919100 0.0004198478  0.03805388     0.2201065     0.2663409     0.2912990     0.3165783     0.3683227 8215.1243 1.0000383
## theta[12]      0.2643035 0.0003818231  0.03686586     0.1954648     0.2391527     0.2628249     0.2890451     0.3397814 9322.3309 0.9997950
## theta[13]      0.2915575 0.0004187071  0.03726929     0.2208711     0.2663372     0.2901505     0.3172342     0.3663520 7922.8532 0.9997019
## theta[14]      0.2501941 0.0004474661  0.03684892     0.1797143     0.2254120     0.2492363     0.2752540     0.3255106 6781.5558 0.9996698
## theta[15]      0.2511710 0.0004257149  0.03696623     0.1800006     0.2267502     0.2502266     0.2751937     0.3277070 7540.0210 0.9997739
## theta[16]      0.2724192 0.0004118264  0.03802260     0.2008984     0.2454932     0.2714542     0.2983664     0.3492893 8524.2280 0.9996837
## theta[17]      0.3687138 0.0006017706  0.04295469     0.2897383     0.3395782     0.3674180     0.3960810     0.4571714 5095.1771 0.9998620
## theta[18]      0.3349516 0.0005211066  0.04090971     0.2593766     0.3062662     0.3337778     0.3619035     0.4199811 6163.1053 1.0000127
## theta[19]      0.2785679 0.0004030213  0.03800727     0.2075628     0.2523441     0.2780688     0.3031291     0.3573856 8893.5955 1.0004885
## theta[20]      0.2923364 0.0003881004  0.03765584     0.2221160     0.2669556     0.2909912     0.3171463     0.3688302 9414.0539 0.9996694
## theta[21]      0.2923944 0.0004105208  0.03752207     0.2209383     0.2665856     0.2912131     0.3164970     0.3680485 8354.1681 0.9999173
## theta[22]      0.2788560 0.0003996906  0.03733491     0.2090631     0.2538924     0.2776593     0.3026852     0.3549537 8725.3379 0.9998152
## theta[23]      0.2928908 0.0004114481  0.03824679     0.2222196     0.2662670     0.2912846     0.3179989     0.3712505 8640.9188 1.0000555
## theta[24]      0.2644473 0.0004153146  0.03671163     0.1949200     0.2390569     0.2632830     0.2891419     0.3378905 7813.6310 0.9999130
## theta[25]      0.3128055 0.0004385096  0.03895709     0.2388283     0.2861345     0.3116477     0.3384651     0.3926516 7892.5041 0.9998247
## theta[26]      0.2986913 0.0004200183  0.03929545     0.2240903     0.2718427     0.2974430     0.3238565     0.3797093 8752.8221 1.0001100
## theta[27]      0.3060174 0.0004120814  0.03868151     0.2333721     0.2798390     0.3048613     0.3318317     0.3845329 8811.3142 1.0003280
## theta[28]      0.2989248 0.0004290114  0.03957001     0.2251945     0.2717036     0.2974334     0.3247360     0.3792510 8507.3569 1.0000547
## theta[29]      0.3061401 0.0004472804  0.03953055     0.2324903     0.2787192     0.3048072     0.3319467     0.3873361 7810.9899 0.9999754
## theta[30]      0.2928218 0.0004119275  0.03830807     0.2202530     0.2666469     0.2913593     0.3180223     0.3715707 8648.4633 0.9998400
## theta[31]      0.2642508 0.0004251426  0.03727759     0.1938331     0.2391509     0.2632630     0.2885809     0.3386117 7688.2297 0.9998156
## theta[32]      0.3058192 0.0004292515  0.03879178     0.2339007     0.2789955     0.3045386     0.3312443     0.3849721 8166.8740 0.9998484
## theta[33]      0.3274622 0.0004201846  0.03946135     0.2542896     0.3001349     0.3263858     0.3531745     0.4092664 8819.9021 0.9997113
## theta[34]      0.2848898 0.0003833603  0.03692793     0.2153058     0.2600348     0.2837112     0.3089870     0.3591429 9278.8845 0.9996602
## theta[35]      0.3414920 0.0005061411  0.04010711     0.2670948     0.3139972     0.3396154     0.3670868     0.4244887 6279.1325 0.9998830
## theta[36]      0.2781216 0.0004056042  0.03820720     0.2070069     0.2519129     0.2764561     0.3038706     0.3565034 8873.3065 0.9997304
## theta[37]      0.3343937 0.0005101920  0.04053305     0.2597412     0.3062192     0.3332731     0.3618454     0.4162061 6311.7724 1.0002280
## theta[38]      0.2710310 0.0004135586  0.03713213     0.2005977     0.2458790     0.2698862     0.2954022     0.3473106 8061.6825 1.0000961
## theta[39]      0.2852881 0.0003816384  0.03794432     0.2131113     0.2597463     0.2837482     0.3100398     0.3630662 9885.2901 0.9998564
## theta[40]      0.2297447 0.0004762846  0.03635624     0.1603494     0.2046664     0.2291501     0.2541400     0.3016015 5826.7267 1.0001311
## theta[41]      0.2090910 0.0005760593  0.03632870     0.1405706     0.1842842     0.2086480     0.2334374     0.2826001 3977.0885 1.0011642
## theta[42]      0.2995069 0.0004164803  0.03896647     0.2270683     0.2731376     0.2981535     0.3250133     0.3794263 8753.7339 0.9997469
## theta[43]      0.2920939 0.0004072035  0.03788036     0.2198289     0.2666680     0.2913138     0.3161774     0.3709978 8653.7661 0.9999513
## theta[44]      0.2715980 0.0004048654  0.03722735     0.2007269     0.2460824     0.2710163     0.2957196     0.3475937 8454.7916 0.9997736
## theta[45]      0.2922163 0.0004003830  0.03841135     0.2189823     0.2662122     0.2910298     0.3168344     0.3698336 9203.8171 1.0000835
## theta[46]      0.3546832 0.0005101880  0.04171936     0.2779532     0.3254227     0.3528418     0.3820520     0.4415040 6686.7446 0.9997996
## theta[47]      0.3059269 0.0004506704  0.03934227     0.2326989     0.2791032     0.3050701     0.3319493     0.3848617 7620.8021 0.9999296
## theta[48]      0.2784598 0.0003887562  0.03692207     0.2080397     0.2535158     0.2777081     0.3026433     0.3527193 9020.2268 0.9997928
## theta[49]      0.3204225 0.0004595129  0.03942231     0.2468644     0.2930574     0.3189228     0.3459725     0.4016114 7360.1856 1.0000674
## theta[50]      0.2781997 0.0003887698  0.03814484     0.2060448     0.2518953     0.2772301     0.3030379     0.3559719 9626.8983 1.0005048
## theta[51]      0.2918973 0.0004146772  0.03905966     0.2188104     0.2649327     0.2907285     0.3174387     0.3725411 8872.3058 0.9996323
## theta[52]      0.2783465 0.0003987032  0.03702295     0.2087075     0.2527973     0.2779787     0.3024503     0.3535757 8622.6850 0.9998797
## theta[53]      0.2984610 0.0004325383  0.03838253     0.2272434     0.2721811     0.2970202     0.3233230     0.3767968 7874.4117 0.9998195
## theta[54]      0.2573666 0.0004639541  0.03777003     0.1863292     0.2318089     0.2565776     0.2824067     0.3324148 6627.4231 1.0002692
## theta[55]      0.2920894 0.0003936185  0.03749328     0.2203673     0.2666705     0.2912075     0.3172593     0.3677279 9073.1062 0.9999331
## theta[56]      0.3128227 0.0004406603  0.03893309     0.2396030     0.2857130     0.3120786     0.3384010     0.3938547 7806.0218 1.0000059
## theta[57]      0.2712220 0.0004232062  0.03689594     0.2016821     0.2464009     0.2701305     0.2951135     0.3451827 7600.6948 0.9996794
## theta[58]      0.2921724 0.0004076824  0.03771271     0.2202968     0.2661902     0.2911625     0.3167772     0.3694597 8557.1962 0.9998622
## theta[59]      0.3264725 0.0004729908  0.03998185     0.2527415     0.2992321     0.3255940     0.3518972     0.4094379 7145.2961 1.0000292
## theta[60]      0.2857274 0.0003970438  0.03797955     0.2138729     0.2601644     0.2843445     0.3112081     0.3615038 9150.0336 0.9999082
## theta[61]      0.2515131 0.0004455979  0.03727520     0.1801835     0.2266858     0.2507712     0.2753264     0.3270156 6997.6746 0.9998969
## theta[62]      0.2298571 0.0005030731  0.03608710     0.1593349     0.2057095     0.2292898     0.2534322     0.3022709 5145.6672 1.0001168
## theta[63]      0.3410529 0.0004986582  0.04115278     0.2645825     0.3129934     0.3399292     0.3663483     0.4268985 6810.7116 1.0000850
## theta[64]      0.2849522 0.0003886903  0.03759836     0.2138198     0.2593137     0.2844731     0.3094068     0.3614122 9356.8648 0.9998152
## theta[65]      0.3336771 0.0004723182  0.04023587     0.2591633     0.3059794     0.3323851     0.3596622     0.4152885 7257.0031 0.9998823
## theta[66]      0.2993343 0.0003916411  0.03752453     0.2287496     0.2740452     0.2980983     0.3236076     0.3753244 9180.2393 0.9996591
## theta[67]      0.2431904 0.0004649919  0.03716742     0.1716023     0.2186350     0.2425173     0.2675705     0.3186644 6389.0192 1.0005350
## theta[68]      0.2573065 0.0004386811  0.03687291     0.1868474     0.2318513     0.2569924     0.2817758     0.3299359 7065.0773 0.9999714
## theta[69]      0.2714580 0.0004005096  0.03769085     0.2008549     0.2454348     0.2711711     0.2963628     0.3475691 8856.1730 0.9998527
## theta[70]      0.2855112 0.0004342834  0.03848544     0.2139243     0.2586820     0.2842544     0.3105883     0.3633542 7853.1957 0.9997942
## theta[71]      0.2579963 0.0004127618  0.03745715     0.1872082     0.2323497     0.2569180     0.2834602     0.3350401 8235.1276 1.0002755
## theta[72]      0.2505337 0.0004100033  0.03683692     0.1819202     0.2254840     0.2494439     0.2753326     0.3247486 8072.1981 0.9996681
## theta[73]      0.3199342 0.0004500445  0.03935861     0.2468163     0.2928441     0.3184078     0.3455522     0.4022636 7648.3622 1.0000962
## theta[74]      0.2573790 0.0004417862  0.03725100     0.1872251     0.2312967     0.2568115     0.2821153     0.3325960 7109.7044 0.9999936
## theta[75]      0.2992437 0.0004111894  0.03817072     0.2273318     0.2736403     0.2980117     0.3240436     0.3765505 8617.4118 0.9998924
## theta[76]      0.3059712 0.0004151197  0.03822134     0.2341072     0.2799243     0.3048010     0.3312102     0.3832158 8477.4482 0.9997668
## theta[77]      0.2576736 0.0004198204  0.03710522     0.1867541     0.2323357     0.2566088     0.2819400     0.3329811 7811.6530 1.0002743
## theta[78]      0.2924022 0.0003978701  0.03830415     0.2210830     0.2667495     0.2914186     0.3175172     0.3697365 9268.4957 0.9997812
## theta[79]      0.2790787 0.0004143416  0.03820962     0.2073831     0.2526505     0.2781731     0.3047033     0.3564319 8504.1004 0.9997270
## theta[80]      0.2849260 0.0003948108  0.03776196     0.2150635     0.2585816     0.2840913     0.3103482     0.3615194 9148.1049 1.0000432
## theta[81]      0.2506442 0.0004118227  0.03584663     0.1821977     0.2254052     0.2499782     0.2747080     0.3227381 7576.6298 0.9999518
## theta[82]      0.2925154 0.0004165494  0.03844362     0.2195283     0.2667300     0.2909902     0.3174840     0.3727409 8517.5677 1.0003595
## theta[83]      0.2508219 0.0004447594  0.03643722     0.1823824     0.2262477     0.2507050     0.2741846     0.3238242 6711.8165 1.0003658
## theta[84]      0.2848214 0.0003959930  0.03769013     0.2136895     0.2587300     0.2841507     0.3093756     0.3618961 9058.9978 0.9996221
## theta[85]      0.3068041 0.0004346726  0.03874220     0.2338116     0.2801897     0.3052356     0.3322418     0.3855132 7944.0876 0.9999517
## theta[86]      0.2995141 0.0003999847  0.03787046     0.2290823     0.2730097     0.2982571     0.3242332     0.3772465 8964.2609 0.9999162
## theta[87]      0.2993881 0.0004029828  0.03822569     0.2270593     0.2737849     0.2984900     0.3234363     0.3786745 8997.8252 0.9997244
## theta[88]      0.2783920 0.0004154374  0.03837951     0.2060811     0.2520944     0.2778191     0.3034364     0.3583292 8534.6878 0.9996714
## theta[89]      0.2434345 0.0004715800  0.03662163     0.1702311     0.2192862     0.2438964     0.2674153     0.3149767 6030.6559 1.0005370
## theta[90]      0.2925462 0.0004028093  0.03829529     0.2189288     0.2665436     0.2914624     0.3176834     0.3709523 9038.4033 0.9997619
## theta[91]      0.2923922 0.0004198644  0.03826316     0.2192218     0.2666656     0.2916989     0.3169135     0.3709579 8305.0770 0.9998412
## theta[92]      0.3059846 0.0004215965  0.03815562     0.2330506     0.2799007     0.3048453     0.3306180     0.3846887 8190.7416 0.9998313
## theta[93]      0.2923159 0.0004056526  0.03722206     0.2232239     0.2664263     0.2916418     0.3167751     0.3683684 8419.6185 0.9998251
## theta[94]      0.2501512 0.0004232692  0.03661459     0.1828189     0.2244148     0.2486887     0.2746619     0.3236030 7482.9873 0.9999446
## theta[95]      0.3054844 0.0004244169  0.03884440     0.2327900     0.2787842     0.3041714     0.3305114     0.3876175 8376.6752 0.9998497
## theta[96]      0.3194484 0.0004464505  0.03899804     0.2459608     0.2930759     0.3178136     0.3444551     0.3993687 7630.2527 1.0000965
## theta[97]      0.3059958 0.0004338551  0.03850390     0.2324243     0.2806907     0.3051084     0.3311212     0.3845098 7876.2629 1.0004275
## theta[98]      0.3196896 0.0004572148  0.03847958     0.2476291     0.2931237     0.3186093     0.3442675     0.3993554 7083.0447 1.0000261
## theta[99]      0.2851845 0.0004202724  0.03835971     0.2128151     0.2587722     0.2840756     0.3102396     0.3636934 8330.8431 0.9999119
## theta[100]     0.2999776 0.0003926419  0.03800322     0.2297143     0.2741687     0.2986800     0.3240882     0.3780506 9368.0143 0.9998246
## alpha         28.4027881 0.2187025384  6.91769972    15.5058941    23.0760503    28.5251951    33.9733535    40.0037681 1000.4969 1.0026045
## beta          70.1630980 0.5322864068 16.99793449    38.4191475    57.1308486    70.5970520    83.9629617    98.2313826 1019.7684 1.0025123
## lp__       -2962.9709479 0.3547401671 10.93604299 -2986.3872794 -2970.2070950 -2962.1079018 -2955.0642400 -2943.9111804  950.3852 1.0030760

Some notes on Stan

Variable types in Stan:

int x;       // an integer
int y[10];   // ten integers
real z;      // a number
real z[2,5]; // a 2x5 array of numbers

vector[10] u;      // length 10 vector
matrix[10,10] v;   // 10x10 matrix
vector[10] w[10];  // ten length 10 vectors

don’t forget the ;
make sure R types match
read the error messages

Steps in fitting a model

Write model as a text block.
Compile the model with stan_model( ).
Sample from the model with sampling( ).
Read warnings and errors.
Check convergence diagnostics (Rhat, stan_trace( ), pairs( ), etc).
Summarize result (credible intervals, stan_plot( ), stan_hist( ), etc).