How to Use PyStan 3
In Windows, PyStan 3 can be used in Ubuntu installed on Windows by WSL (Windows Subsystem for Linux). For some information about WSL, check this website.
PyStan 3 can be installed on Python 3.11 virtual environment of Anaconda in Ubuntu.
But, before installing PyStan 3, install g++ on Ubuntu by the following command:
sudo apt
install g++
Then install PyStan 3 on Python 3.11 virtual environment of Anaconda by running the following command:
pip install
pystan
To import pystan 3, run the following code
import stan
Two sample scripts, which use PyStan 3, are shown below. The script in Listing 1 analyzes data for a binomial model, and that in Listing 2 for a multivariate normal model. Sample script files are archived in the file scriptfiles.zip .
A sample script for a binomial model is shown in Listing 1.
Listing 1. Binomial model(binomial.py).
import stan
import arviz as az
import scipy.stats as ss
import matplotlib.pyplot as plt
model_script = """
data {
int
N;
int
k;
}
parameters {
real<lower = 0.0, upper = 1.0> p;
}
model {
k ~
binomial(N, p);
}
"""
Stan_Data = {'N': 10, 'k': 7}
sm = stan.build(model_script, data =
Stan_Data)
fit = sm.sample()
#
# Convert sampling
data into an InferenceData class
#
i_data = az.from_pystan(posterior=fit,
posterior_model=sm)
smry = az.summary(i_data)
print('Summary...\n', smry)
az.plot_trace(i_data)
plt.tight_layout()
plt.savefig('FigTrace.png')
plt.show()
az.plot_posterior(i_data,
point_estimate = 'median')
plt.savefig('FigRho.png')
plt.show()
Execute the script , binomial.py, in Listing 1 by running the following command:
(py311) ******/scriptfiles$
python binomial.py
After MCMC sampling, a trace plot is displayed (Figure 1).
Figure 1
Close the window of Figure 1, then posterior distribution of p is displayed.
Figure 2
Close the window of Figure 2, then the script ends.
The following message is displayed on the display.
Summary...
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
p
0.669 0.128 0.44 0.909 0.003 0.002 1476.0 1697.0 1.0
(py311) ******/scriptfiles$
A sample script for a multivariate normal model is shown in Listing 2.
Listing 2. Multivariate normal model (multivarnormal.py).
import stan
import arviz as az
import scipy.stats as ss
import matplotlib.pyplot as plt
Data = ss.multivariate_normal.rvs(mean
= [0,0], cov = [[1,0.7],[0.7,1]],
size = 100)
model_script = """
data {
int
N;
array[N] vector[2] XY;
}
parameters {
real
mu1;
real
mu2;
real<lower = 0> v_sgm1;
real<lower = 0> v_sgm2;
real<lower=-1, upper=1> rho;
}
transformed parameters {
vector[2] mus;
matrix[2,2] cov;
real
sgm1;
real
sgm2;
sgm1
= v_sgm1 * 10.0;
sgm2 = v_sgm2 * 10.0;
mus[1] = mu1;
mus[2] = mu2;
cov[1][1] = sgm1^2;
cov[2][2] = sgm2^2;
cov[1][2] = sgm1*sgm2*rho;
cov[2][1] = cov[1][2];
}
model {
mu1
~ normal(0,10);
mu2
~ normal(0,10);
v_sgm1 ~ exponential(0.1);
v_sgm2 ~ exponential(0.1);
rho
~ uniform(-1, 1);
XY ~
multi_normal(mus, cov);
}
"""
Stan_Data = {'N':len(Data), 'XY':Data}
sm = stan.build(model_script, data =
Stan_Data)
fit = sm.sample()
#
# Convert sampling
data into an InferenceData class
#
i_data = az.from_pystan(posterior=fit,
posterior_model=sm)
smry = az.summary(i_data)
print('Summary...\n', smry)
az.plot_trace(i_data, var_names =
['mu1', 'mu2', 'sgm1', 'sgm2', 'rho'])
plt.tight_layout()
plt.savefig('FigTrace.png')
plt.show()
az.plot_posterior(i_data,
point_estimate = 'median', var_names = ['rho'])
plt.savefig('FigRho.png')
plt.show()
Run the script in Listing 2 by the following command:
(py311) ******/scriptfiles$ python multivarnormal.py
After MCMC sampling, a trace plot graph is displayed (Figure 3).
Figure 3
Close the window of Figure 3, the posterior distribution of rho is displayed (Figure 4).
Figure 4
Close the window of Figure 4, then the script ends.
The following message is displayed on the display.
Summary...
mean
sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
mu1
-0.112 0.109 -0.326 0.078 0.002 0.002 2204.0 2390.0 1.0
mu2
-0.053 0.110 -0.260 0.153 0.002 0.002 2226.0 2835.0 1.0
v_sgm1 0.106 0.008 0.092 0.120 0.000 0.000 2441.0 2601.0 1.0
v_sgm2 0.108 0.008 0.094 0.123 0.000 0.000 2198.0 2115.0 1.0
rho
0.755 0.044 0.673 0.839 0.001 0.001 2265.0 2174.0 1.0
mus[0] -0.112 0.109 -0.326 0.078 0.002 0.002 2204.0 2390.0 1.0
mus[1] -0.053 0.110 -0.260 0.153 0.002 0.002 2226.0 2835.0 1.0
cov[0, 0] 1.119 0.163 0.850 1.442 0.003 0.002 2441.0 2601.0 1.0
cov[0, 1] 0.870 0.150 0.609 1.158 0.003 0.003 1904.0 1940.0 1.0
cov[1, 0] 0.870 0.150 0.609 1.158 0.003 0.003 1904.0 1940.0 1.0
cov[1, 1] 1.180 0.176 0.865 1.504 0.004 0.003 2198.0 2115.0 1.0
sgm1 1.055 0.076 0.922 1.201 0.002 0.001 2441.0 2601.0 1.0
sgm2 1.083 0.080 0.939 1.234 0.002 0.001 2198.0 2115.0 1.0
(py311) ******/scriptfiles$