Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
PyStan version of two compartment model from "Stan: A probabilistic programming language for Bayesian inference and optimization" Gelman, Lee, Guo (2015)
import pystan
import numpy as np
# Two compartment model from
# "Stan: A probabilistic programming language for
# Bayesian inference and optimization" Gelman, Lee, Guo (2015)
# http://www.stat.columbia.edu/~gelman/research/published/stan_jebs_2.pdf
a = np.array([0.8, 1.0])
b = np.array([2, 0.1])
sigma = 0.2
x = np.arange(0, 1000, dtype='float')/100
N = len(x)
# The two compartment model we are attempting to fit
y_pred = a[0]*np.exp(-b[0]*x) + a[1]*np.exp(-b[1]*x)
# Include multiplicative noise
y = y_pred * np.exp(np.random.normal(0, sigma, N))
# Compile and sample model
fit = pystan.stan(file='two_compartment.stan',
data={'N': N, 'x': x, 'y': y},
iter=1000, chains=4)
# Plot parameter estimates of interest
fit.plot(pars=['a', 'b', 'sigma'])
# Print all parameter estimates (limitation of PyStan 2.0)
print(fit)
data {
int N;
vector[N] x;
vector[N] y;
}
parameters {
vector[2] log_a;
ordered[2] log_b;
real<lower=0> sigma;
}
transformed parameters {
vector<lower=0>[2] a;
vector<lower=0>[2] b;
a = exp(log_a);
b = exp(log_b);
}
model {
vector[N] y_pred;
y_pred = a[1]*exp(-b[1]*x) + a[2]*exp(-b[2]*x);
y ~ lognormal(log(y_pred), sigma);
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment