K&M model in Python using Stan

From: Bayesian Models for Astrophysical Data, Cambridge Univ. Press

you are kindly asked to include the complete citation if you used this material in a publication

Code 9.4 K and M model in Python using Stan

===================================

import numpy as np
import pystan

from scipy.linalg import toeplitz
from scipy.stats import norm, uniform, nbinom, multivariate_normal, bernoulli

# Data
np.random.seed(1056)

nobs = 500
nvar = 15
rho = 0.6
p = bernoulli.rvs(0.2, size=nvar)

# use these to generate new values of beta
# beta1 = p * norm.rvs(loc=0, scale=5.0, size=nvar)
# beta = np.array([round(item, 2) for item in beta1])

# use same beta values as Code 9.1
beta = np.array([-2.43, 1.60, 0.00, 5.01, 4.12, 0.00, 0.00, 0.00,
-0.89, 0.00, -2.31, 0.00, 0.00, 0.00, 0.00])

# covariance matrix
d = beta.shape[0]
Sigma = toeplitz(np.insert(np.repeat(rho, d-1), 0, 1))

# multivariate sampling - default mean is zero
M = multivariate_normal.rvs(cov=Sigma, size=nobs)
xb = np.dot(M, beta)

y = [norm.rvs(loc=xb[i], scale=2.0) for i in range(xb.shape[0])]

# fit
mydata = {}
mydata['X'] = M - 1.0
mydata['K'] = mydata['X'].shape[1]
mydata['Y'] = y
mydata['N'] = nobs
mydata['Ind'] = bernoulli.rvs(0.2, size=nvar)

stan_model = '''
data{
int<lower=0> N;
int<lower=0> K;
matrix[N, K] X;
vector[N] Y;
int Ind[K];
}
parameters{
vector[K] beta;
real<lower=0> sigma;
real<lower=0> sdBeta;
}
transformed parameters{
vector[N] mu;
real<lower=0> tau;
real<lower=0> tauBeta;

mu = X * beta;
tau = pow(sigma, 2);
tauBeta = pow(sdBeta, 2);
}
model{

sdBeta ~ gamma(0.01, 0.01);

for (i in 1:K){
if (Ind[i] > 0) beta[i] ~ normal(0, tauBeta);
}

sigma ~ gamma(0.01, 0.01);

Y ~ normal(mu, tau);
}
'''

fit = pystan.stan(model_code=stan_model, data=mydata, iter=5000, chains=3, thin=1,
warmup=2500, n_jobs=3)

# Output
nlines = 21 # number of lines in screen output

output = str(fit).split('\n')

for item in output[:nlines]:
print(item)

===================================

GET SOURCE

Output on screen:

Inference for Stan model: anon_model_06a21ab591fac250162f7b0f276b42d7.
3 chains, each with iter=5000; warmup=2500; thin=1;
post-warmup draws per chain=2500, total post-warmup draws=7500.

mean se_mean sd 2.5% 25% 50% 75% 97.5%     n_eff Rhat
beta[0] -2.52 3.1e-3 0.27 -3.05 -2.7   -2.52 -2.34    -1.99    7500 1.0
beta[1]    1.11 2.9e-3 0.26    0.61    0.93   1.11    1.28       1.6    7500    1.0
beta[2]    -0.44 3.1e-3 0.27   -0.97   -0.63    -0.44    -0.27    0.09    7500 1.0
beta[3] 4.75 3.0e-3    0.26    4.25    4.58       4.75    4.93        5.27    7500 1.0
beta[4]   4.23    2.9e-3    0.25    3.75    4.06       4.23 4.4        4.72    7500 1.0
beta[5] -0.32    3.0e-3 0.26   -0.83    -0.5      -0.32    -0.15        0.17    7500 1.0
beta[6]    -0.59    2.9e-3    0.25    -1.09   -0.76      -0.59 -0.42        -0.09    7500 1.0
beta[7]    -0.27    2.9e-3    0.25    -0.76   -0.44    -0.27    -0.1        0.24        7500 1.0
beta[8]    -1.27    3.0e-3 0.26   -1.79    -1.45      -1.28 -1.1     -0.76        7500 1.0
beta[9]    0.28 2.8e-3 0.24 -0.2    0.12       0.28 0.44           0.76      7500 1.0
beta[10]   -2.41    2.8e-3    0.24   -2.88 -2.57      -2.41 -2.25        -1.93        7500 1.0
beta[11]   -0.21 2.9e-3    0.25    -0.7 -0.38      -0.21    -0.04     0.28       7500 1.0
beta[12]   -0.11    2.9e-3 0.25 -0.58   -0.28    -0.11 0.06        0.38    7500 1.0
beta[13] -1.4e-3   2.8e-3 0.24    -0.48   -0.16   2.3e-3    0.16        0.48    7500 1.0
beta[14]   -0.25    2.9e-3    0.26   -0.74    -0.42    -0.25 -0.08     0.25        7500 1.0
sigma 1.91    3.5e-4 0.03    1.86    1.89     1.91    1.93     1.98    7500 1.0

HSI

HSI