In [1]:

    

import pystan


    

In [2]:

    

stan_code = """
## http://openbugs.net/Examples/Surgical.html
## random effects model 
data {
  int<lower=0> N;
  int r[N];
  int n[N];
}
parameters {
   real mu;
   real<lower=0> sigmasq;
   real b[N];
}
transformed parameters {
  real<lower=0> sigma;
  real<lower=0,upper=1> p[N];
  sigma <- sqrt(sigmasq); 
  for (i in 1:N)
    p[i] <- inv_logit(b[i]);
}
model {
  mu ~ normal(0.0, 1000.0); 
  sigmasq ~ inv_gamma(0.001, 0.001);
  b ~ normal(mu, sigma);
  r ~ binomial_logit(n, b);
}
generated quantities {
  real pop_mean;
  pop_mean <- inv_logit(mu); 
}
"""


    

In [3]:

    

stan_model = pystan.StanModel(model_code=stan_code)


    

    




INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_b75445567ed4638c70765f73ec25f477 NOW.

In [4]:

    

init = {
    'b': [0.1,   0.1,   0.1,   0.1,   0.1,   0.1,   0.1,   0.1,   0.1,   0.1,   0.1,   0.1],
    'sigmasq': 1,
    'mu': 0
}


    

In [13]:

    

data = {
    'n': [47, 148, 119, 810, 211, 196, 148, 215, 207, 97, 256, 360],
    'r': [0, 18, 8, 46, 8, 13, 9, 31, 14, 8, 29, 24],
    'N': 12
}


    

In [15]:

    

init = [init, init, init, init] # one for each chain


    

In [16]:

    

fit = stan_model.sampling(init=init, data=data, iter=1000, chains=4)


    

In [17]:

    

fit


    

    
Out[17]:





Inference for Stan model: anon_model_b75445567ed4638c70765f73ec25f477.
4 chains, each with iter=1000; warmup=500; thin=1; 
post-warmup draws per chain=500, total post-warmup draws=2000.

           mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
mu        -2.56  3.9e-3   0.15  -2.86  -2.65  -2.55  -2.47  -2.28   1498    1.0
sigmasq    0.19  4.9e-3   0.15   0.03   0.09   0.15   0.22   0.61    997    1.0
b[0]      -2.95    0.01   0.42  -3.94  -3.19   -2.9  -2.66  -2.27   1131    1.0
b[1]      -2.19  5.4e-3   0.24  -2.64  -2.36  -2.19  -2.02  -1.71   2000    1.0
b[2]      -2.62  6.2e-3   0.28   -3.2  -2.79  -2.61  -2.44  -2.11   2000    1.0
b[3]      -2.77  3.3e-3   0.15  -3.07  -2.87  -2.77  -2.67  -2.49   2000    1.0
b[4]      -2.94  6.2e-3   0.28  -3.53  -3.12  -2.93  -2.75  -2.46   2000    1.0
b[5]      -2.62  5.2e-3   0.23  -3.13  -2.76  -2.61  -2.47  -2.19   2000    1.0
b[6]      -2.67  5.8e-3   0.26  -3.23  -2.84  -2.67   -2.5  -2.19   2000    1.0
b[7]      -1.98  6.2e-3   0.21   -2.4  -2.11  -1.97  -1.84  -1.59   1102    1.0
b[8]      -2.61  4.8e-3   0.22  -3.05  -2.75   -2.6  -2.46  -2.21   2000    1.0
b[9]       -2.5  6.2e-3   0.28   -3.1  -2.67  -2.49  -2.33  -1.95   2000    1.0
b[10]     -2.19  4.2e-3   0.19  -2.57  -2.32  -2.19  -2.07  -1.83   2000    1.0
b[11]     -2.62  4.1e-3   0.18  -2.96  -2.74  -2.62  -2.51  -2.27   2000    1.0
sigma       0.4  5.6e-3   0.15   0.16    0.3   0.38   0.47   0.78    736    1.0
p[0]       0.05  4.9e-4   0.02   0.02   0.04   0.05   0.07   0.09   1484    1.0
p[1]        0.1  5.0e-4   0.02   0.07   0.09    0.1   0.12   0.15   2000    1.0
p[2]       0.07  3.9e-4   0.02   0.04   0.06   0.07   0.08   0.11   2000    1.0
p[3]       0.06  1.8e-4 8.2e-3   0.04   0.05   0.06   0.06   0.08   2000    1.0
p[4]       0.05  2.9e-4   0.01   0.03   0.04   0.05   0.06   0.08   2000    1.0
p[5]       0.07  3.3e-4   0.01   0.04   0.06   0.07   0.08    0.1   2000    1.0
p[6]       0.07  3.5e-4   0.02   0.04   0.06   0.07   0.08    0.1   2000    1.0
p[7]       0.12  5.0e-4   0.02   0.08   0.11   0.12   0.14   0.17   2000    1.0
p[8]       0.07  3.1e-4   0.01   0.05   0.06   0.07   0.08    0.1   2000    1.0
p[9]       0.08  4.4e-4   0.02   0.04   0.06   0.08   0.09   0.12   2000    1.0
p[10]       0.1  3.8e-4   0.02   0.07   0.09    0.1   0.11   0.14   2000    1.0
p[11]      0.07  2.6e-4   0.01   0.05   0.06   0.07   0.08   0.09   2000    1.0
pop_mean   0.07  2.5e-410.0e-3   0.05   0.07   0.07   0.08   0.09   1553    1.0
lp__     -724.9    0.13   3.17 -732.0 -726.8 -724.6 -722.6 -719.5    584    1.0

Samples were drawn using NUTS at Thu Dec  7 18:16:55 2017.
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).