stan_ep1

6_1_1_catalog.py

import pystan
import numpy as np

stan_code = """
data {
    int<lower=0> N;
    real<lower=0> x[N];
}
parameters {
    real mu;
    real<lower=0> sigma;
}
model {
    for(n in 1:N) {
        x[n] ~ normal(mu, sigma);
    }
}
generated quantities {
    real<lower=0, upper=1> mu_over;
    real<lower=0, upper=1> mu_over2;
    real es;
    real<lower=0, upper=1> es_over;

    mu_over <- step(mu - 2500);
    mu_over2 <- step(mu - 3000);
    es <- (mu - 2500) / sigma;
    es_over <- step(es-0.8);
}
"""

sample = [3060, 2840, 1780, 3280, 3550, 2450, 2200, 3070, 2100, 4100,
          3630, 3060, 3280, 1870, 2980, 3120, 2150, 3830, 4300, 1880]

stan_model = pystan.StanModel(model_code=stan_code)
print np.mean(sample)
op = stan_model.sampling(data=dict(x=sample, N=len(sample)), chains=1, iter=10000)
print op