top of page

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

(c) 2017,  Joseph M. Hilbe, Rafael S. de Souza and Emille E. O. Ishida

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

â€‹

â€‹

â€‹

Code 10.2 Normal linear model, in Python using Stan, for assessing the relationship between central black hole mass and bulge velocity dispersion

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

import numpy as np
import pandas as pd
import pystan

â€‹

â€‹

â€‹

# prepare data for Stan
data = {}
data['obsx'] = np.array(data_frame['obsx'])
data['errx'] = np.array(data_frame['errx'])
data['obsy'] = np.array(data_frame['obsy'])
data['erry'] = np.array(data_frame['erry'])
data['N'] = len(data['obsx'])

â€‹

# Stan Gaussian model with errors
stan_code="""
data{
int<lower=0> N;                                                          # number of data points
vector[N] obsx;                                                            # obs velocity dispersion
vector<lower=0>[N] errx;                                           # errors in velocity dispersion measurements
vector[N] obsy;                                                            # obs black hole mass
vector<lower=0>[N] erry;                                           # errors in black hole mass measurements
}
parameters{
real alpha;                                                                    # intercept
real beta;                                                                      # angular coefficient
real<lower=0> epsilon;                                               # scatter around true black hole mass
vector[N] x;                                                                  # true velocity dispersion
vector[N] y;                                                                  # true black hole mass
}
model{

# likelihood
obsx ~ normal(x, errx);
y ~ normal(alpha + beta * x, epsilon);
obsy ~ normal(y, erry);
}
"""

â€‹

# Run mcmc
fit = pystan.stan(model_code=stan_code, data=data, iter=15000, chains=3,
warmup=5000, thin=10, n_jobs=3)

â€‹

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

â€‹

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

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

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

Output on screen:

â€‹

Inference for Stan model: anon_model_17e10ffd8db7fefeb22cee9be2011e11.
3 chains, each with iter=15000; warmup=5000; thin=10;
post-warmup draws per chain=1000, total post-warmup draws=3000.

â€‹

mean     se_mean         sd       2.5%        25%        50%        75%       97.5%        n_eff        Rhat
alpha         8.35         1.1e-3      0.06       8.23         8.31        8.35        8.39           8.46         2941          1.0
beta           4.45         6.2e-3      0.33         3.8         4.22        4.45        4.67           5.11         2864          1.0
epsilon      0.29         1.2e-3      0.06       0.18         0.25        0.29        0.32           0.41         2564          1.0

bottom of page