In [3]:

    

from matplotlib.pyplot import plot
from ssa import SSA
from pystan import stan

%matplotlib inline


    

In [5]:

    

N = 100
alpha = 1
mu = 10

x, t = SSA(100,N,a=alpha,mu=mu)

x = x.astype(int) # path data supposed to be integers.

path_data = {
    'N' : N,
    't' : t,
    'x' : x
}


    

In [7]:

    

# Setup STAN :
model_description = """
data{
      int<lower=0> N; ## number of time steps
      vector[N] t; ## time value at each time step
      int<lower=0> x[N]; ## population value at each time step
    }

    transformed data{
      int<lower=0> x0; ## starting population value
      x0 <- x[1];
    }

    parameters{
      real<lower=0> alpha; ## birth rate parameter 
      real<lower=0> mu; ## death rate parameter
    }

    model {
          x ~ poisson(alpha/mu+(x0-alpha/mu)*exp(-mu*t)); 
    }
"""


    

In [8]:

    

fit = stan(model_code=model_description,
           data=path_data, chains=1,iter=1000)

print fit


    

    




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

        mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
alpha   5.71    0.22   5.01   0.17   1.83    4.4   8.05  19.84  500.0    1.0
mu      8.32  8.6e-3   0.19   7.99   8.18   8.31   8.44   8.72  500.0    1.0
lp__   1.6e4    0.04   0.95  1.6e4  1.6e4  1.6e4  1.6e4  1.6e4  500.0    1.0

Samples were drawn using NUTS(diag_e) at Thu Apr  7 19:54:06 2016.
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).

In [ ]: