Using Stan (PyStan 3) in Python
Stan is a famous library for Bayesian analysis, and PyStan is the Python interface for Stan. PyStan can be easily installed using Anaconda. How to install Anaconda and PyStan is explained at this website and at that one.
The sample scripts shown below were executed with PyStan 3.9.1 in Python 3.11 virtual environment in Anaconda on Ubuntu installed on Windows.
Sample scripts were prepared for a univariate normal distribution model and a bivariate normal one. Files used at the website were compressed into the file filesPS3.zip, which can be downloaded and from which the files can be extracted, freely.
Univariate Normal Distribution Model
Stan script is written in a string form, and used as the first argument of the Stan build function. But, the string can be given as a return value of a file read function.
In the following, in the first example, Stan script is given as a string in a Python script, and in the second example, Stan script is given in a separate file from the Python script file.
Stan script in the Python script
Listing 1 shows an example of PyStan 3 script for a univariate normal distribution model.
Listing 1. Stan script for univariate normal
distribution (OneVar.stan)
import stan
import matplotlib.pyplot as plt
import numpy as np
import arviz as az
OneVar_code = """
data {
int<lower=0> J;
array[J] real y;
}
parameters {
real
vmu;
real<lower=0> vsigma;
}
transformed parameters {
// Rescaling
of vmu and vsigma
real
mu;
real
sigma;
mu =
vmu * 1000.0;
sigma = vsigma * 1000.0 + 0.001;
}
model {
vmu
~ normal(0.0, 1.0);
vsigma ~ exponential(1.0);
y ~
normal(mu, sigma);
}
"""
LData = [ 44, 44, 56, 59, 44, 44, 63, 60, 34, 58,
55, 50, 32, 61, 47, 25, 38, 47, 30, 42]
OneVar_dat = { 'J' : len(LData), 'y':
LData }
sm = stan.build(OneVar_code,
data=OneVar_dat)
fit = sm.sample()
#
# InferenceData class
#
infdata = az.from_pystan(posterior =
fit, posterior_model = sm)
summary = az.summary(infdata,
var_names=['mu', 'sigma'])
print('\nSummary...\n', summary)
az.plot_trace(infdata,
var_names=['mu', 'sigma'])
plt.tight_layout()
plt.show()
#
# Pandas DataFrame
#
dframe = fit.to_frame()
#
# Calculation of
median based-summaries,
# which are
recommended by Gelman et al.(2021)
#
mu_med = np.median(dframe['mu'])
mu_mad_sd = 1.483 *
np.median(np.abs(dframe['mu'] - mu_med))
print(f'\nmu_med = {mu_med:.3f}, mu_mad_sd = {mu_mad_sd:.3f}')
sigma_med = np.median(dframe['sigma'])
sigma_mad_sd = 1.483 *
np.median(np.abs(dframe['sigma'] - sigma_med))
print(f'sigma_med =
{sigma_med:.3f}, sigma_mad_sd
= {sigma_mad_sd:.3f}')
The Stan script is given as a string OneVar_code.
To set weak prior distributios, the parameters are rescaled as follows:
real mu;
real sigma;
mu = vmu * 1000.0;
sigma = vsigma * 1000.0 + 0.001;
The Stan script is used in the following build function:
sm = stan.build(OneVar_code, data=OneVar_dat)
MCMC sampling is run by the following code:
fit = sm.sample()
To analyze the MCMC sample, the result fit by sm.sample() is transformed to InferenceData and DataFrame class objects as follows:
infdata = az.from_pystan(posterior =
fit, posterior_model = sm)
and
dframe = fit.to_frame()
Using infdata, summary statistics are calculated and shown on the terminal by the following scripts:
summary = az.summary(infdata,
var_names=['mu', 'sigma'])
print('\nSummary...\n', summary)
Summary statistics are shown as follows:
Summary...
mean
sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
mu 46.627 2.638 41.400 51.518 0.045 0.032 3542.0 2778.0 1.0
sigma 11.949 2.124 8.274 15.682 0.053 0.038 1677.0 2042.0 1.0
Trace plot is shown by the following script:
az.plot_trace(infdata,
var_names=['mu', 'sigma'])
plt.tight_layout()
plt.show()
We get trace plot as shown in Figure 1.
Figure 1
Since dframe is a Pandas DataFrame, a sample from the posterior distribution for parameter mu is obtained by dframe['mu'].
From this, we can calculate the median and the median absolute deviation(MAD). Gelman et al. (2021, p.73) recommend median-based summaries, instead of the mean. They use “mad sd”, which replaces sd. mad sd is given by
Note that
Medians and mad sds for parameters mu and sigma are calculated by the following script:
mu_med = np.median(dframe['mu'])
mu_mad_sd = 1.483 *
np.median(np.abs(dframe['mu'] - mu_med))
print(f'\nmu_med = {mu_med:.3f}, mu_mad_sd = {mu_mad_sd:.3f}')
sigma_med = np.median(dframe['sigma'])
sigma_mad_sd = 1.483 *
np.median(np.abs(dframe['sigma'] - sigma_med))
print(f'sigma_med =
{sigma_med:.3f}, sigma_mad_sd
= {sigma_mad_sd:.3f}')
The following messages are displayed on the terminal.
mu_med = 46.677, mu_mad_sd = 2.733
sigma_med = 11.697, sigma_mad_sd = 1.995
Stan script and Python script are in separate files
The Stan script can be written in a file other than in the Python script file. The Stan script of the Python script in Listing 1 is extracted and constructed as the file named OneVar.stan, which is shown in Listing 2.
Listing 2. The Stan script extracted from Listing 1
(OneVar.stan)
data {
int<lower=0> J;
array[J] real y;
}
parameters {
real vmu;
real<lower=0> vsigma;
}
transformed parameters {
//
Rescaling of vmu and vsigma
real
mu;
real
sigma;
mu =
vmu * 1000.0;
sigma = vsigma * 1000.0 + 0.001;
}
model {
vmu
~ normal(0.0, 1.0);
vsigma ~ exponential(1.0);
y ~
normal(mu, sigma);
}
To use the file OneVar.stan, the python script in Listing 1 was rewritten as shown in Listing3.
Listing 3. The Python script, which uses the Stan
script file OneVar.stan of Listing 2.
import stan
import matplotlib.pyplot as plt
import numpy as np
import arviz as az
LData = [ 44, 44, 56, 59, 44, 44, 63, 60, 34,
58,
55, 50, 32, 61, 47, 25, 38, 47, 30, 42]
OneVar_dat = { 'J' : len(LData), 'y':
LData }
with open('OneVar.stan', 'r') as
f_stan:
sm =
stan.build(f_stan.read(), data=OneVar_dat)
fit = sm.sample()
#
#
InferenceData class
#
infdata = az.from_pystan(posterior =
fit, posterior_model = sm)
summary = az.summary(infdata,
var_names=['mu', 'sigma'])
print('\nSummary...\n', summary)
az.plot_trace(infdata,
var_names=['mu', 'sigma'])
plt.tight_layout()
plt.show()
#
# Pandas DataFrame
#
dframe = fit.to_frame()
#
# Calculation of
median based-summaries,
# which are
recommended by Gelman et al.(2021)
#
mu_med = np.median(dframe['mu'])
mu_mad_sd = 1.483 *
np.median(np.abs(dframe['mu'] - mu_med))
print(f'\nmu_med = {mu_med:.3f}, mu_mad_sd = {mu_mad_sd:.3f}')
sigma_med = np.median(dframe['sigma'])
sigma_mad_sd = 1.483 *
np.median(np.abs(dframe['sigma'] - sigma_med))
print(f'sigma_med =
{sigma_med:.3f}, sigma_mad_sd
= {sigma_mad_sd:.3f}')
The script in Listing 3 is the same as in Listing 1, except it uses the Stan script file OneVar.stan.
In the script in Listing 3, the build function is called as follows:
with open('OneVar.stan', 'r') as
f_stan:
sm =
stan.build(f_stan.read(), data=OneVar_dat)
The Stan script is read in from the file, so the stan script code OneVar_code in Listing 1 is removed. The other parts in Listing 3 are the same as ones in Listing 1.
Bivariate Normal Distribution Model
Stan script for a bivariate normal model is shown in Listing 4.
Listing 4. Stan script for a bivariate normal model
(TwoVarNormal.stan)
data {
int
N;
array[N]
vector[2] Data;
real
l1;
real
u1;
real
l2;
real
u2;
real
sgm_l1;
real
sgm_u1;
real
sgm_l2;
real
sgm_u2;
}
parameters {
real<lower
= l1, upper = u1> mu1;
real<lower
= l2, upper = u2> mu2;
real<lower
= sgm_l1, upper = sgm_u1> sgm1;
real<lower
= sgm_l2, upper = sgm_u2> sgm2;
real<lower
= -1.0, upper = 1.0> r;
}
transformed parameters {
vector[2]
mu;
matrix[2,
2] cov;
mu[1]
= mu1;
mu[2]
= mu2;
cov[1][1]
= sgm1 * sgm1;
cov[1][2]
= sgm1 * sgm2 * r;
cov[2][1]
= sgm1 * sgm2 * r;
cov[2][2]
= sgm2 * sgm2;
}
model {
for
(i in 1:N){
Data[i]
~ multi_normal(mu, cov);
}
}
Prior distributions are given as weakly informative ones (Gelman et al. 2014, p. 55). They are set as uniform distributions, whose ranges depend on statistics of the data.
The stan script in Listing 4 is used by the Python script in Listing 5.
Listing 5. Python script which uses the Stan script in Listing 4
import numpy as np
import stan
import matplotlib.pyplot as plt
import arviz as az
Data = [
[179, 121], [185, 89], [181, 106], [210,
100], [202, 104],
[197, 91], [227, 107], [208, 101], [188,
99], [189, 112],
[240, 111], [203, 103], [202, 85], [190,
95], [221, 96],
[150, 83], [205, 103], [177, 90], [221,
92], [172, 86],
[190, 87], [218, 101], [205, 100], [208,
120], [224, 98],
[199, 116], [203, 96], [213, 91], [200,
109], [190, 100],
[211, 115], [204, 85], [185, 92], [178,
86], [231, 122],
[206, 84], [197, 93], [236, 93], [213,
106], [205, 103],
[212, 102], [177, 86], [193, 96], [183,
93], [204, 87],
[203, 108], [203, 103], [180, 100], [196,
101], [167, 75],
[202, 106], [182, 98], [225, 109], [219,
115], [205, 102],
[200, 99], [219, 94], [189, 102], [220,
86], [186, 102],
[203, 101], [210, 124], [254, 103], [208,
95], [192, 106],
[176, 91], [173, 97], [195, 100], [189,
98], [220, 115],
[202, 112], [257, 108], [207, 104], [199,
96], [173, 98],
[230, 101], [230, 111], [212, 105], [196,
82], [183, 82],
[230, 103], [207, 112], [225, 101], [192,
104], [175, 83],
[223, 105], [192, 106], [202, 87], [204,
87], [184, 97],
[193, 94], [204, 109], [196, 92], [216,
109], [219, 113],
[256, 117], [210, 86], [217, 103], [206,
87], [213, 105] ]
mu = np.mean(Data, axis = 0)
print('mu =\n', mu)
cov = np.cov(np.array(Data).T)
print('cov =\n', cov)
sgms = [cov[0][0]**0.5,
cov[1][1]**0.5]
print('sgm =\n', sgms)
L1 = mu[0] - 5 * sgms[0]
U1 = mu[0] + 5 * sgms[0]
L2 = mu[1] - 5 * sgms[1]
U2 = mu[1] + 5 * sgms[1]
sgm_L1 = 0.0
sgm_U1 = 10 * sgms[0]
sgm_L2 = 0.0
sgm_U2 = 10 * sgms[1]
stan_data = {'N': len(Data), 'Data':
Data, 'l1': L1, 'u1': U1, 'l2': L2, 'u2': U2,
'sgm_l1':
sgm_L1, 'sgm_u1': sgm_U1,
'sgm_l2':
sgm_L2, 'sgm_u2': sgm_U2}
with open('TwoVarNormal.stan', 'r') as
f_stan:
sm =
stan.build(f_stan.read(), data=stan_data)
fit = sm.sample()
#
#
InferenceData class
#
infdata = az.from_pystan(posterior =
fit, posterior_model = sm)
summary =
az.summary(infdata,var_names=['mu1', 'mu2', 'sgm1', 'sgm2', 'r'])
print('\nSummary...\n', summary)
az.plot_trace(infdata,
var_names=['mu1', 'mu2', 'sgm1',
'sgm2', 'r'],
figsize=(10,6))
plt.tight_layout()
plt.show()
#
# Pandas DataFrame
#
dframe = fit.to_frame()
#
# Calculation of
median based-summaries,
# which are
recommended by Gelman et al.(2021)
mu1_med = np.median(dframe['mu1'])
mu1_mad_sd = 1.483 *
np.median(np.abs(dframe['mu1'] - mu1_med))
print(f'\nmu1_med =
{mu1_med:.3f}, mu1_mad_sd =
{mu1_mad_sd:.3f}')
mu2_med = np.median(dframe['mu2'])
mu2_mad_sd = 1.483 *
np.median(np.abs(dframe['mu2'] - mu2_med))
print(f'mu2_med = {mu2_med:.3f}, mu2_mad_sd = {mu2_mad_sd:.3f}')
sgm1_med = np.median(dframe['sgm1'])
sgm1_mad_sd = 1.483 *
np.median(np.abs(dframe['sgm1'] - sgm1_med))
print(f'sgm1_med =
{sgm1_med:.3f}, sgm1_mad_sd =
{sgm1_mad_sd:.3f}')
sgm2_med = np.median(dframe['sgm2'])
sgm2_mad_sd = 1.483 *
np.median(np.abs(dframe['sgm2'] - sgm2_med))
print(f'sgm2_med =
{sgm2_med:.3f}, sgm2_mad_sd =
{sgm2_mad_sd:.3f}')
r_med = np.median(dframe['r'])
r_mad_sd = 1.483 *
np.median(np.abs(dframe['r'] - r_med))
print(f'r_med = {r_med:.3f}, r_mad_sd = {r_mad_sd:.3f}')
The Stan script is read in and compiled by the following code:
with open('TwoVarNormal.stan', 'r') as
f_stan:
sm =
stan.build(f_stan.read(), data=stan_data)
MCMC sampling is run by the code:
fit =
sm.sample()
To analyze the MCMC sample, the result fit by sm.sample() is transformed to InferenceData and DataFrame class objects as follows:
infdata = az.from_pystan(posterior =
fit, posterior_model = sm)
and
dframe = fit.to_frame()
Statistics of MCMC sample are calculated and shown by the following script:
summary =
az.summary(infdata,var_names=['mu1', 'mu2', 'sgm1', 'sgm2', 'r'])
print('\nSummary...\n', summary)
We get the following output on the terminal display.
Summary...
mean
sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
mu1 203.002 1.969 199.381 206.723 0.031 0.022 3930.0 2854.0 1.0
mu2 99.618 1.033 97.774 101.614 0.017 0.012 3864.0 2748.0 1.0
sgm1 19.435 1.390 16.920 22.125 0.023 0.016 3702.0 3079.0 1.0
sgm2 10.352 0.771 8.980 11.775 0.012 0.009 3975.0 3005.0 1.0
r 0.423 0.082 0.271 0.579 0.001 0.001 4363.0 3016.0 1.0
Trace plot is shown by the following script:
az.plot_trace(infdata,
var_names=['mu1', 'mu2', 'sgm1',
'sgm2', 'r'],
figsize=(10,6))
plt.tight_layout()
plt.show()
We get trace plot like Figure 2:
Figure 2
As statistics which describe posterior distributions, Gelman et al. (2021, p.73) recommend median-based summaries, instead of means and sd’s.
We calculate the median and the median absolute deviation(MAD), from which mad sd is calculated.
mad sd is given by
Note that
Medians and mad sd’s are calculated and shown by the following script:
mu1_med = np.median(dframe['mu1'])
mu1_mad_sd = 1.483 *
np.median(np.abs(dframe['mu1'] - mu1_med))
print(f'\nmu1_med =
{mu1_med:.3f}, mu1_mad_sd =
{mu1_mad_sd:.3f}')
mu2_med = np.median(dframe['mu2'])
mu2_mad_sd = 1.483 *
np.median(np.abs(dframe['mu2'] - mu2_med))
print(f'mu2_med = {mu2_med:.3f}, mu2_mad_sd = {mu2_mad_sd:.3f}')
sgm1_med = np.median(dframe['sgm1'])
sgm1_mad_sd = 1.483 *
np.median(np.abs(dframe['sgm1'] - sgm1_med))
print(f'sgm1_med =
{sgm1_med:.3f}, sgm1_mad_sd =
{sgm1_mad_sd:.3f}')
sgm2_med = np.median(dframe['sgm2'])
sgm2_mad_sd = 1.483 *
np.median(np.abs(dframe['sgm2'] - sgm2_med))
print(f'sgm2_med =
{sgm2_med:.3f}, sgm2_mad_sd =
{sgm2_mad_sd:.3f}')
r_med = np.median(dframe['r'])
r_mad_sd = 1.483 *
np.median(np.abs(dframe['r'] - r_med))
print(f'r_med = {r_med:.3f}, r_mad_sd = {r_mad_sd:.3f}')
We get these statistics as follows:
mu1_med = 203.015, mu1_mad_sd = 1.961
mu2_med = 99.664, mu2_mad_sd = 1.025
sgm1_med = 19.318, sgm1_mad_sd = 1.330
sgm2_med = 10.289, sgm2_mad_sd = 0.733
r_med = 0.426, r_mad_sd = 0.082
Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., & Rubin, D. B. (2014). Bayesian data analysis, third edition. CRC Press.
Gelman, A., Hill, J., & Vehtari A. (2021). Regression and other stories. Cambridge University Press.