Trying PyStan 3
I tried PyStan 3 (2022.11.22) in Ubuntu. I installed PyStan 3 in a Python 3.10 virtual environment in Anaconda in Ubuntu.
Install PyStan 3 in the Python 3.10 virtual environment by running the following pip command.
pip install
pystan
When using PyStan 3 in a Python script, import the module as follows:
import stan
In the following, examples of a binomial model, a simple regression model, and a bivariate normal model are shown.
In Python 3.10 virtual environment in Anaconda in Ubuntu of WSL in Windows 10, I had trouble that matplotlib caused (2022.12.04). To avoid this trouble, I did MCMC sampling using PyStan 3 in the Python 3.10 virtual environment of Ubuntu onWSL, and saved the samples in files. Then, in Python 3.11 in Windows 10, I loaded the data of sampling from the files, and analyzed them.
For details, check this.
But, now (2023.06), matplotlib causes no
trouble in WSL of Windows. So, we can do MCMC sampling by PyStan 3 and analyzing
the sample by matplotlib both in Ubuntu on WSL of Windows.
I tried a script in Listing 1.
Listing 1. A sample script of a binomial model.
import stan
import matplotlib.pyplot as plt
import seaborn as sb
import arviz as az
binomial_model = """
data {
int
N;
int
k;
}
parameters {
real<lower = 0.0, upper = 1.0> p;
}
model {
k ~
binomial(N, p);
}
"""
binomial_data = {'N':10, 'k':7}
sm = stan.build(binomial_model, data =
binomial_data)
fit = sm.sample()
print('fit =\n', fit, '\n')
p_smpls = fit['p']
print('p_smpls =\n', p_smpls)
sb.kdeplot(p_smpls[0])
plt.savefig('FigKDE.png')
plt.show()
plt.close()
df = fit.to_frame()
sb.kdeplot(df['p'])
plt.savefig('FigPD.png')
plt.show()
plt.close()
inference_data =
az.from_pystan(posterior = fit, posterior_model = sm)
az.plot_trace(inference_data)
plt.savefig('FigTrace.png')
plt.show()
plt.close()
Execute the script in Listing 1, the following messages are displayed on the terminal.
fit =
<stan.Fit>
Parameters:
p: ()
Draws: 4000
p_smpls =
[[0.62546345 0.68693505 0.64059199 ...
0.76143154 0.69859012 0.78014042]]
The following graphs are dislayed and saved in files.
FigKDE.png
FigPD.png
FigTrace.png
A script of a simple regression model is shown in Listing 2.
Listing 2. A script of a simple regression model
import stan
import scipy.stats as ss
import matplotlib.pyplot as plt
import seaborn as sb
import arviz as az
Xs = ss.uniform.rvs(loc = 10, scale =
20, size = 100)
Es = ss.norm.rvs(loc = 0, scale = 1,
size = 100)
Ys = 10 + 0.5 * Xs + Es
plt.scatter(Xs, Ys)
plt.savefig('FigXY.png')
plt.show()
plt.close()
simple_reg_model = """
data {
int
N;
vector[N] X;
vector[N] Y;
}
parameters {
real
a;
real
b;
real<lower = 0> sgm;
}
model {
a ~
normal(0, 10);
b ~
normal(0, 10);
sgm
~ exponential(0.1);
Y ~
normal(a + b*X, sgm);
}
"""
simple_reg_data = {'N':len(Xs),
'X':Xs, 'Y':Ys}
sm = stan.build(simple_reg_model, data
= simple_reg_data)
fit = sm.sample()
print('fit =\n', fit, '\n')
infer_data =
az.from_pystan(posterior=fit, posterior_model=sm)
smry = az.summary(infer_data)
print(smry)
az.plot_trace(infer_data)
plt.savefig('FigTrace.png')
print('FigTrace.png was saved.')
plt.show()
plt.close()
a_smpls = fit['a']
b_smpls = fit['b']
plt.scatter(a_smpls, b_smpls)
plt.xlabel('a')
plt.ylabel('b')
plt.savefig('Fig_ab.png')
plt.show()
plt.close()
Run the script in Listing 2, the following messages are displayed on the terminal.
fit =
<stan.Fit>
Parameters:
a: ()
b: ()
sgm: ()
Draws: 4000
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
a 9.913 0.362 9.211 10.556 0.010 0.007 1272.0 1507.0 1.0
b 0.502 0.018 0.469 0.535 0.000 0.000 1273.0 1509.0 1.0
sgm 0.996 0.074 0.860 1.128 0.002 0.001 1856.0 1865.0 1.0
The following graphs are displayed and saved in files.
FigXY.png
FigTrace.png
Fig_ab.png
A sample script of a bivariate normal model is shown in Listing 3.
Listing 3. A sample script of a bivariate normal model
import stan
import scipy.stats as ss
import matplotlib.pyplot as plt
import seaborn as sb
import arviz as az
X =
ss.multivariate_normal.rvs([0,0],[[1, 0.8],[0.8,1]], size = 100)
print(X)
plt.scatter(X.T[0], X.T[1])
plt.savefig('FigXY.png')
plt.show()
plt.close()
binormal_model = """
data {
int
N;
vector[2] X[N];
}
parameters {
vector[2] mu;
vector[2] sgm;
real<lower = -1, upper = 1> rho;
}
transformed parameters {
cov_matrix[2] cov;
cov[1][1] = sgm[1]*sgm[1];
cov[1][2] = rho*sgm[1]*sgm[2];
cov[2][1] = rho*sgm[1]*sgm[2];
cov[2][2]
= sgm[2]*sgm[2];
}
model {
mu[1] ~ normal(0, 10);
mu[2] ~ normal(0, 10);
sgm
~ exponential(0.1);
rho
~ uniform(-1,1);
X ~
multi_normal(mu, cov);
}
"""
binormal_data = {'N':len(X), 'X':X}
sm = stan.build(binormal_model, data =
binormal_data)
fit = sm.sample()
print('fit =\n', fit, '\n')
inference_data =
az.from_pystan(posterior = fit, posterior_model = sm)
smry = az.summary(inference_data)
print('Summary =\n',smry)
az.plot_trace(inference_data)
plt.savefig('FigTrace.png')
print('FigTrace.png was saved.')
plt.show()
plt.close()
mu_smpls = fit['mu']
print('mu_smpls =\n', mu_smpls)
df = fit.to_frame()
print('df =\n', df, '\n\n')
for v in zip(df['mu.1'][:5],
df['mu.2'][:5], df['cov.2.2'][:5]):
print(v)
sb.kdeplot(df['rho'])
plt.savefig('FigRho.png')
plt.show()
plt.close()
Run the script in Listing 3, the following messages are shown on the terminal.
fit =
<stan.Fit>
Parameters:
mu: (2,)
sgm: (2,)
rho: ()
cov: (2, 2)
Draws: 4000
Summary =
mean
sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
mu[0] 0.187 0.105 -0.007 0.383 0.002 0.001 2472.0 2695.0 1.01
mu[1] 0.210 0.097 0.024 0.388 0.002 0.001 2498.0 2694.0 1.00
sgm[0] 1.065 0.076 0.928 1.210 0.002 0.001 2180.0 2205.0 1.00
sgm[1] 0.984 0.069 0.853 1.112 0.001 0.001 2366.0 2355.0 1.00
rho 0.770 0.040 0.691 0.839 0.001 0.001 2295.0 2620.0 1.00
cov[0,0] 1.140 0.164 0.862 1.463 0.004 0.003 2180.0 2205.0 1.00
cov[0,1] 0.813 0.132 0.582 1.062 0.003 0.002 1880.0 2135.0 1.00
cov[1,0] 0.813 0.132 0.582 1.062 0.003 0.002 1880.0 2135.0 1.00
cov[1,1] 0.974 0.138 0.721 1.229 0.003 0.002 2366.0 2355.0 1.00
FigTrace.png was saved.
mu_smpls =
[[0.18040335 0.16903765 0.14978146 ...
0.43856426 0.01501792 0.24142373]
[0.20285295 0.15172237 0.22550795 ...
0.43964007 0.14393881 0.22252268]]
df =
parameters lp__ accept_stat__ stepsize__ treedepth__ ... cov.1.1 cov.2.1 cov.1.2 cov.2.2
draws
...
0
-58.263544
0.827678
0.458740
3.0 ... 1.133158 0.922741 0.922741 1.122094
1
-64.044667
0.562658
0.457346
2.0 ... 1.353705 0.763930 0.763930 0.725858
2
-57.713119
0.988581 0.430489
2.0 ... 1.037119 0.775723 0.775723 0.910998
3
-58.287036
0.957960
0.466302
3.0 ... 1.052202 0.752236 0.752236 0.959562
4
-57.652823
0.958025
0.458740
3.0 ... 1.041470 0.715984 0.715984 0.926422
...
...
...
...
... ... ... ... ... ...
3995
-59.680630
0.975303
0.466302
3.0 ... 1.104269 0.755304 0.755304 1.030205
3996 -57.830879
0.998800
0.458740
4.0 ... 0.920630 0.625547 0.625547 0.809076
3997
-60.108280
0.971335
0.457346
3.0 ... 1.202069 0.903426 0.903426 1.077758
3998
-58.572426
0.983354
0.430489
3.0 ...
1.125427 0.783323 0.783323 0.908921
3999
-58.904657
0.998236
0.466302
3.0 ... 1.055764 0.672293 0.672293 0.739689
[4000 rows x 16 columns]
(0.18040335026726861,
0.20285294781077376, 1.122094106032836)
(0.16903765146625016,
0.15172236853563245, 0.7258576891080232)
(0.1497814603342748,
0.22550794533864094, 0.9109977685363938)
(0.22654457223130214,
0.3351371542279542, 0.959561947061149)
(0.2556725357232705,
0.2796258006706391, 0.9264216126600003)
The following graphs are displayed and saved in files.
FigXY.png
FigTrace.png
FigRho.png
Trying
PyStan 3 in a Python 3.10 virtual environment of Anaconda in Ubuntu in WSL in
Windows 10
I could not use matplotlib in a Python 3.10 virtual environment in Anaconda in WSL in Windows 10 (2022.12.04). So, I did analysis in two steps: Step 1 and Step 2.
But, now (2023.06), matplotlib causes no trouble, so sampling by PyStan 3 and analyses by matplotlib can be done both in Ubuntu/WSL/Windows.
Step 1: MCMC sampling using PyStan 3 in a Python 3.10 virtual Anaconda in Ubuntu in WSL in Windows 10.
I did MCMC sampling by PyStan 3 in WSL, and saved the sampling data into files using the script in Listing W1.
Listing W1 Script for MCMC sampling using PyStan 3
import stan
import arviz as az
import pickle
binomial_model = """
data {
int
N;
int
k;
}
parameters {
real<lower = 0.0, upper = 1.0> p;
}
model {
k ~
binomial(N, p);
}
"""
binomial_data = {'N':10, 'k':7}
sm = stan.build(binomial_model, data =
binomial_data)
fit = sm.sample()
with open('fit.pkl', 'wb') as f:
pickle.dump(fit['p'], f)
print('fit.pkl was saved.')
#
# Convert data of
fit into an InferenceData object
#
inference_data =
az.from_pystan(posterior = fit, posterior_model = sm)
smry = az.summary(inference_data)
print('Summary =\n',smry)
with open('i_data.pkl', 'wb') as f:
pickle.dump(inference_data,
f)
print('i_data.pkl was saved.')
#
# Convert data of
fit into a Pandas DataFRame
#
df = fit.to_frame()
with open('d_frame.pkl', 'wb') as f:
pickle.dump(df, f)
print('d_frame.pkl was saved.')
print('\nStep-1 ended.')
Run the script in Listing W1, the following output will be shown after MCMC sampling
fit.pkl was saved.
Summary =
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
p
0.67 0.125 0.433 0.89 0.003 0.002 1630.0 1936.0 1.0
i_data.pkl was saved.
d_frame.pkl was saved.
Step-1 ended.
After execution of the script completed, the files fit.pkl, i_data.pkl, d_frame.pkl of the sampling data were saved.
In the next step, Step 2, the sampling data in the files are loaded, and they are processed in Python 3.11 in Windows 10.
Step 2: Drawing the posterior distribution in Python 3.11 in Windows 10
The script in Listing W2 loads the files saved by the script in Listing W1, and processes them in Python 3.11 in Windows 10.
Listing W2 Script drawing the graph.
import matplotlib.pyplot as plt
import seaborn as sb
import arviz as az
import pickle
with open('fit.pkl', 'rb') as f:
p_smpls =
pickle.load(f) # An ndarray of the samples
print('p_smpls =\n', p_smpls)
sb.kdeplot(p_smpls[0], label = 'p')
plt.xlabel('p')
plt.legend()
plt.show()
with open('i_data.pkl', 'rb') as f:
inference_data =
pickle.load(f) # An InferenceData object
smry = az.summary(inference_data)
print('Summary =\n',smry)
az.plot_trace(inference_data)
plt.tight_layout()
plt.show()
with open('d_frame.pkl', 'rb') as f:
df = pickle.load(f)
# A Pandas DataFrame
print('df =\n', df, '\n\n')
sb.kdeplot(df['p'], label = 'p')
plt.xlabel('p')
plt.legend()
plt.show()
Run the script in Listing W2, the graph of the posterior distribution of parameter p will be shown (Figure W1).
Figure W1
Close the window shown in Figure W1, then the next graph will be shown (Figure W2).
Figure W2
Close the window shown in Figure W2, then the next graph will be shown (Figure W3).
Figure W3
Close the window in Figure W3, then script completes.
In the terminal window, the following ourput are displayed.
p_smpls =
[[0.63121086 0.58010096 0.6628399 ... 0.65729322 0.69610772 0.79043097]]
Summary =
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
p
0.67 0.125 0.433 0.89 0.003 0.002 1630.0 1936.0 1.0
df =
parameters lp__ accept_stat__ ...
energy__
p
draws
...
0
-7.671043
1.000000 ... 7.779637 0.631211
1
-7.827389
1.000000 ... 7.932186 0.580101
2
-7.638564 0.995635 ...
7.703345 0.662840
3
-8.332925
0.800349 ... 8.926897 0.810378
4
-8.184562
0.895592 ... 8.185370 0.517570
...
...
... ... ... ...
3995
-7.951949 0.882511 ...
8.476777 0.767434
3996
-7.638595
0.995531 ... 7.688506 0.662693
3997
-7.640521
0.999555 ... 7.652110 0.657293
3998
-7.662335
1.000000 ... 7.714833 0.696108
3999
-8.132224
0.962709 ... 8.167014 0.790431
[4000 rows x 8 columns]