Using Stan (PyStan) 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 in this website.
The following sample programs use PyStan for Bayesian analysis of univariate and bivariate (the simplest case of multivariate) data. The files of sample scripts are archived in files.zip , which can be downloaded and the extracted files can be used freely.
The sample programs are written in PyStan 2. For PyStan 3, check this website.
Stan script for univariate normal distribution model is shown in Listing 1.
Listing
1 Stan script for univariate normal
distribution (OneVar.stan)
data {
int<lower=0> J;
real
y[J];
}
parameters {
real<lower=-1000, upper=1000> mu;
real<lower=0, upper=1000> sigma;
}
transformed parameters {
}
model {
y ~
normal(mu, sigma);
}
Listing 2 shows Python script, which uses the Stan script of Listing 1 to analyze univariate data.
Listing
2 Python script (OneVar.py), which
uses Stan script (OneVar.stan; Listing 1)
import pystan
import matplotlib.pyplot as plt
import numpy as np
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 =
pystan.StanModel(file='OneVar.stan')
fit = sm.sampling(data=OneVar_dat,
iter=5000, chains=4, n_jobs=1)
print(fit)
LResults = fit.extract(permuted=True)
Lmu = LResults['mu']
meanMu = np.mean(Lmu)
plt.figure()
plt.title("Posterior Distribution
of mu Mean =
{0:.7}".format(meanMu))
plt.hist(Lmu, bins = 50)
plt.show()
Lsigma = LResults['sigma']
meanSigma = np.mean(Lsigma)
plt.figure()
plt.title("Posterior Distribution
of sigma Mean =
{0:.7}".format(meanSigma))
plt.hist(Lsigma, bins = 50)
plt.show()
The Python script OneVar.py can be executed in a virtual environment in which the module pystan is installed, as shown in Figure 1. Why the pystan module is installed in a virtual environment is explained in this website.
Figure 1
When sampling by Stan finished, first the posterior distribution of mu is shown (Figure 2).
Figure 2
Close the window in Figure 2, then the posterior distribution of sigma is shown (Figure 3).
Figure 3
Close the window in Figure 3, the execution of the script ends (Figure 4).
Figure 4
The Stan script in Listing 1 can be included in the Python script in Listing 2 as string code. The combined script is shown in Listing 3. The Stan script is coded as string OneVar_code.
Listing
3 Stan script in Python script
(OneVar.py)
import pystan
import matplotlib.pyplot as plt
import numpy as np
OneVar_code = """
data {
int<lower=0> J;
real
y[J];
}
parameters {
real<lower=-1000,
upper=1000> mu;
real<lower=0, upper=1000> sigma;
}
transformed parameters {
}
model {
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 =
pystan.StanModel(model_code=OneVar_code)
fit = sm.sampling(data=OneVar_dat,
iter=5000, chains=4, n_jobs=1)
print(fit)
LResults = fit.extract(permuted=True)
Lmu = LResults['mu']
meanMu = np.mean(Lmu)
plt.figure()
plt.title("Posterior Distribution
of mu Mean =
{0:.7}".format(meanMu))
plt.hist(Lmu, bins = 50)
plt.show()
Lsigma = LResults['sigma']
meanSigma = np.mean(Lsigma)
plt.figure()
plt.title("Posterior Distribution
of sigma Mean = {0:.7}".format(meanSigma))
plt.hist(Lsigma, bins = 50)
plt.show()
Run the script by python command, we see the same execution as in the case, where the Stan code is prepared as a separate file (Listing 1).
Figure 5
Stan script for bivariate normal model is shown in Listing 4.
Listing
4 Stan script for bivariate normal
model (TwoVarNormal.stan)
data {
int
N;
vector[2]
Data[N];
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 {
mu1
~ uniform(l1, u1);
mu2
~ uniform(l2, u2);
sgm1
~ uniform(sgm_l1, sgm_u1);
sgm2
~ uniform(sgm_l2, sgm_u2);
r
~ uniform(-1.0, 1.0);
for
(i in 1:N){
Data[i]
~ multi_normal(mu, cov);
}
}
The Python script in Listing 5 uses the Stan script in Listing 4, and analyze the data given by a list object Data.
Listing
5 Python script which uses the Stan
script in Listing 4.
import numpy as np
import pystan
import matplotlib.pyplot as plt
import seaborn as sb
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]
def f_init():
return
{'mu1': mu[0], 'mu2': mu[1], 'sgm1': sgms[0], 'sgm2': sgms[1],
'r':
0.0}
sm = pystan.StanModel(file =
'TwoVarNormal.stan')
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}
fit = sm.sampling(data = stan_data,
init = f_init, n_jobs = 1)
print(fit)
samples_mu1 = fit['mu1']
samples_mu2 = fit['mu2']
plt.title('Posterior Distributions',
fontsize = 20)
sb.kdeplot(samples_mu1, label = 'mu1')
sb.kdeplot(samples_mu2, label = 'mu2')
plt.legend(fontsize = 20)
plt.show()
samples_sgm1 = fit['sgm1']
samples_sgm2 = fit['sgm2']
plt.title('Posterior Distributions',
fontsize = 20)
sb.kdeplot(samples_sgm1, label =
'$\sigma_1$')
sb.kdeplot(samples_sgm2, label =
'$\sigma_2$')
plt.legend(fontsize = 20)
plt.show()
samples_r = fit['r']
plt.title('Posterior Distributions',
fontsize = 20)
sb.kdeplot(samples_r, label = 'r')
plt.xlabel('r', fontsize = 20)
plt.legend(fontsize = 20)
plt.show()
Run the script in Listing 5 (Figure 6).
Figure 6
After sampling by MCMC, the posterior distribution of mu1 and mu2 are displayed (Figure 7).
Figure 7
Close the window in Figure 7, then the posterior distributions of sgm1 and sgm2 are shown (Figure 8).
Figure 8
Close the window in Figure 8, then the posterior distribution of r is displayed (Figure 9).
Figure 9
Close the window in Figure 9, then the program ends (Figure 10).
Figure 10