PyStan 3 試用記
2022.10, 2024.09
WinodwsとUbuntuとCentOS Stream 9でPyStan 3を使ってみた（2024.09）。Windowsは、WindowsにインストールされたUbuntu22（以下、Ubuntu22/Windowsと略記）での使用である。WindowsへのUbuntuのインストールは、このウェブサイトなどを参照されたい。PCに直接インストールされたUbuntu22（以下、Ubuntu22/PCと略記）でのPyStan3の使用は、CentOS Stream 9の場合と同様であった。
Ubuntu 22をインストールするときは、PowerShellにおいてコマンド
wsl.exe
--install Ubuntu-22.04
を実行すればよい。
Ubuntu22/WindowsにおいてPyStan 3が使い難くなっていた。しばらくすれば、Ubuntu22/Windowsにおいても、CentOS Stream 9の場合と同様に使い易くなると思われるが、まず、CentOS Stream 9での使用を報告し、その後、Ubuntu22/Windowsの場合について説明する。いずれの場合も、Anacondaをインストールして、Python 3.12の仮想環境を用意し、PyStan 3をコマンド
pip install
pystan
でインストールした。Pythonのスクリプト内でPyStan 3をimportするときは
import stan
である。
PyStan 3をインストールするとき、あるいはPyStan 3実行するとき、インストールできない、あるいは実行できないときは、C++コンパイラがインストールされていない可能性がある。コマンド「g++」を実行してみて、「g++」が見つかりません、コマンド「*****」でインストールして下さいというようなメッセージが表示されたら、その通り実行してみる。
「g++」のインストールが済めば、PyStan 3のインストールができて、PyStan 3スクリプトの実行ができるはずである。
まず、CentOS Stream 9およびUbuntu22/PCでの利用について説明する。Ubuntu22/Windowsの場合については、このウェブページの後半で説明する。
スクリプトファイルは、まとめて圧縮ファイルsample_script.zipとした。
PyStanについては、＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞で説明している。
CentOS Stream 9およびPCに直接インストールされたUbuntu22/PCの場合
CentOS Stream 9の場合も、Ubuntu22/PCの場合も同じであるので、CentOS Stream 9の場合について説明する。Anacondaをインストールして、AnacondaのPython仮想環境においてPyStan 3を利用する。Anacondaのインストールについては、別のウェブサイトで説明している。
Anacondaのインストールの最後において、初期化を行うかどうか聞いてくるが、「yes」と答えておくと楽である。インストール後、日本語入力などで問題が生じたときは、コマンド「conda config –set
auto_activate_base false」を実行すればよい。実行後、再起動すれば、端末の起動時に自動的に仮想環境がアクティブになっているということは無くなるが、コマンド「conda activate」を実行すればアクティブになる。
実行環境は、CentOS
Stream 9 にインストールしたAnacondaの仮想環境としてPython 3.12の仮想環境を用意した。PyStan 3は、用意した仮想環境においてコマンド
pip
install pystan
によりインストールした。
必要に応じてライブラリをインストールするが、scipyはバージョン1.12のものをコマンド
conda
install scipy=1.12
でインストールする。すでに他のバージョンのものがインストールされているときは、それと置き換えられる。
以下、２項分布、単回帰モデル、２変量正規分布について説明する。
リスト１にスクリプト例を示す。
リスト１のスクリプトを実行すると、FigKDE.pngが表示される。
FigKDE.png
これは、スクリプト
sb.kdeplot(p_smpls[0])
により表示されるものである。
FigKDE.pngの表示されているWindowを閉じると、FigPD.pngが表示される。
FigPD.png
これは、コード
sb.kdeplot(df['p'])
により表示されるものである。いずれも、pの事後分布で、同じものである。
図FigPD.pngのWindowを閉じると、図FigTrace.pngが表示される。
FigTrace.png
図FigTrace.pngのWindowを閉じると、スクリプトの実行終了である。
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()
単回帰モデルの例として、リスト２を用意した。
リスト２を実行すると、下図の散布図FigXY.pngが表示される。
FigXY.png
これは、次のスクリプトにより作成されたデータXsとYsの散布図である。
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)
このデータに、次の単回帰モデルを当てはめている。
Y ~
normal(a + b*X, sgm);
図FigXY.pngのWindowを閉じると、トレース図FigTrace.pngが表示される。
FigTrace.png
トレース図FigTrace.pngのWindowを閉じると、係数aとbの事後分布の散布図Fig_ab.pngが表示される。
Fig_ab.png
独立変数の中心化（centering）を行っていないので、強い相関が認められる。
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()
2変量正規分布を用いたスクリプト例をリスト３のように用意した。
リスト３を実行すると、下図の散布図FigXY.pngが表示される。
FigXY.png
これは、次のコードで生成されたデータの散布図である。
X =
ss.multivariate_normal.rvs([0,0],[[1, 0.8],[0.8,1]], size = 100)
このデータに対して、次の多変量正規分布のモデルを適用している。
X ~ multi_normal(mu, cov);
散布図FigXY.pngのWindowを閉じると、トレース図FigTrace.pngが表示される。
FigTrace.png
トレース図FigTrace.pngのWindowを閉じると、相関係数rhoの事後分布が表示される。
FigRho.png
事後分布のグラフFigRho.pngのWindowを閉じると、スクリプトの実行終了である。
import stan
import scipy.stats as ss
import numpy as np
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;
array[N] vector[2] X;
}
parameters {
vector[2] mu;
vector<lower=0.0>[2] v_sgm;
real<lower = -1, upper = 1> rho;
}
transformed parameters {
cov_matrix[2] cov;
vector[2] sgm;
sgm[1] = 0.001 + v_sgm[1];
sgm[2] = 0.001 + v_sgm[2];
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 ~
normal(0, 10);
v_sgm ~ exponential(0.1);
rho
~ uniform(-1,1);
X ~
multi_normal(mu, cov);
}
"""
X = np.array(X)
print(len(X))
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)
#plt.figure(figsize=(10,7))
az.plot_trace(inference_data,
figsize=(10,6))
plt.tight_layout()
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()
WindowsにインストールされたUbuntu22/Windowsの場合におけるPyStan 3の利用
WindowsにインストールされたUbuntu22/Windowsの場合は、StanによるMCMCサンプリングのみをUbuntu22/Windowsで行い、MCMCサンプルをファイルに保存して、Windowsにおいて保存したMCMCサンプルを読み出して分析する。
UbuntuはWindowsのWSLにインストールしたUbuntu 22を利用した。Ubuntu 22/Windowsにおいて、上のCentOS Stream 9と同様の手続きで、Python 3.12の仮想環境をAnaconda/Ubuntu22/WSL/Windowsに用意してPyStan 3をインストールした。Scipyはバージョン1.12のscipy1.12をインストールした。
以上の環境の下で、2項分布のMCMCサンプリングを行うリストU1をPython3.12/Anaconda/Ubuntu22/WSL/Windowsにおいて実行した。
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()
df = fit.to_frame()
with open('df.pkl', 'wb') as f:
pickle.dump(df, f)
# DataFrame型を保存
inference_data = az.from_pystan(posterior
= fit, posterior_model = sm)
with open('inf_d.pkl', 'wb') as f:
pickle.dump(inference_data,
f) # InferenceData型を保存
smry = az.summary(inference_data)
print('Summary =\n',smry)
# MCMCサンプリングの統計量の出力
MCMCサンプルは次のコードで、DataFrame型のものとInferenceData型のものの2つを保存している。
df = fit.to_frame()
with open('df.pkl', 'wb') as f:
pickle.dump(df, f)
# DataFrame型を保存
inference_data =
az.from_pystan(posterior = fit, posterior_model = sm)
with open('inf_d.pkl', 'wb') as f:
pickle.dump(inference_data,
f) # InferenceData型を保存
上の保存されたデータdf.pklとinf_d.pklを、Windows上において読み込み、処理を行うスクリプトを、リストU２のように用意した。
リストU2 MCMCサンプルのファイルを読み込み処理するスクリプト
import matplotlib.pyplot as plt
import seaborn as sb
import arviz as az
import pickle
with open('df.pkl', 'rb') as f:
df = pickle.load(f)
sb.kdeplot(df['p'])
plt.savefig('FigPD.png')
plt.show()
plt.close()
with open('inf_d.pkl', 'rb') as f:
inference_data =
pickle.load(f)
az.plot_trace(inference_data)
plt.savefig('FigTrace.png')
plt.show()
plt.close()
print(az.summary(inference_data))
リストU2のスクリプトを、リストU1で保存されたファイルの置かれているフォルダで実行すると、まず、ｐの事後分布のグラフ図Ｕ１が表示される。
図U1
図Ｕ１のWindowを閉じると、トレース図（図U2）が表示される。
図U2
図U2のWindowsを閉じると、スクリプトの実行は終了して、端末には次の出力が表示されている。
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
p
0.666 0.135 0.394 0.884 0.004 0.003 1298.0 1717.0 1.0
これは、次のコード
print(az.summary(inference_data))
によるものである。