PyStan 2.19.1.1のWindowsでの使用例
2026.06改訂
PyStan 2.19.1.1をWindowsで使ってみた。
PyStan 2.19.1.1を使うために、WSL上にPython 3.8の仮想環境をMiniconda/Ubuntu/WSL/Windowsに作成した。Python 3.8仮想環境上でリスト１のスクリプトを実行した。リスト１のスクリプトは、岡本（2019）「いまさら聞けないPythonでデータ分析」丸善出版の第11章の改訂版である。スクリプトファイルなどは圧縮ファイルRegressWSL.zipにまとめた。
WindowsにLinux（Ubuntuなど）をインストールする方法は、このウェブサイトなどで説明している。現時点（2026.05.29）では、UnbuntuはUbuntu 24がよいようである。Ubuntu 26は、まだ不安定で使い難いことがある。Ubuntuのインストール後、PyStanはC++コンパイラが必要なので、g++をインストールすればよい。インストールは、次の子マンでできる。
sudo apt install g++
PyStan 2を用いるためにPythonのバージョンを選ぶが、これは仮想環境を利用すればよい。仮想環境は、Minicondaを利用すると簡単にできるが、UbuntuへのMinicondaのインストールは、このウェブサイトで説明している。
Minicondaをインストール後、Ubuntuを再起動した後、Python3.8の仮想環境を作成する。コマンドは次のようである。
(base) yasuh@pcyo:~$ conda create -n
py38 python=3.8
Python 3.8の仮想環境を作成したら、その仮想環境をactiveにして、pystanをインストールする。コマンドは、以下のようである。
(base) yasuh@pcyo:~$ conda activate py38
(py38) yasuh@pcyo:~$ conda install pystan
Activeな仮想環境名が、先頭に（py38）のように表示されていることに注意。
他のライブラリもインストールする。次のようなコマンドでインストールしておけば、たいていの場合、間に合うと思われる。
(py38) yasuh@pcyo:~$ conda install numpy matplotlib scipy seaborn
pandas openpyxl arviz
他のライブラリも必要になれば、その都度インストールすればよい。condaコマンドでインストールできない（エラーになる）ときは、pipコマンドを試してみればよい。例えば、次のようである。
(py38) ****/RegressWSL$ pip install japanize_matplotlib
リスト１のスクリプトファイル（Regression.py）の置かれているフォルダにカレントディレクトリを移す。Pythonスクリプトファイルをコマンドpythonにより実行するときは、カレントディレクトリを実行するスクリプトファイルの置かれているフォルダに移してからpythonコマンドでPythonスクリプトを実行する。これに対して、conda コマンドあるいはpipコマンドでライブラリのインストールを行うときは、カレントディレクトリがどこであるかは関係ない。そのときのカレントディレクトリでcondaコマンドなどを実行すればよい。
次のコマンドでは、lsコマンドでカレントディレクトリ内のファイルを確認した後、pythonコマンドでRegression.pyを実行している。Regression.pyはリスト１のスクリプトのファイル名である。
(py38) ****/RegressWSL$
ls
Regression.py Regression.stan
(py38) ****/RegressWSL$
python Regression.py
リスト１では、配列
RawData = [ ['City', 'Temperature',
'Latitude'],
['Sapporo',
15,
43.1],
・
・
・
['Naha',
26,
26.2]
]
により、データが設定されている。Temperatureは５日後の予想最高気温であり、https://tenki.jp/から読み込み編集したものである。Latitudeは、https://www.geocoding.jp/から読み込み編集した。
データRawDataに設定されている内容を、次のスクリプトでYに気温を、Xに緯度を、IDに都市名を、Nにデータ数を設定している。
RawData = np.array(RawData)
# 数値データは、文字列に変換されている
print(RawData)
ID = RawData[1:,0]
Y = [float(v) for v in RawData[1:,1]]
X = [float(v) for v in RawData[1:,2]]
N = len(ID)
リスト２では、次の単回帰モデルを設定している。
上式を、次の２つの式に分解してStanスクリプトが書かれている。
リスト２のStanスクリプトによるMCMCサンプリングが、リスト１の次のスクリプトにより行われている。
Data = {'N': N, 'Y': Y, 'X': X}
sm = pystan.StanModel(file='Regression.stan')
fit = sm.sampling(data=Data,
control={'adapt_delta':0.999})
サンプリングの結果fitのトレース図が、次のスクリプトによりArvizの関数plot_traceで描画される。
i_data = az.from_pystan(posterior
= fit, posterior_model = sm)
az.plot_trace(i_data)
plt.show()
図１のようなトレース図が得られる。
図１
図１のWindowを閉じると、次のスクリプトによりパラメータb0の事後分布のヒストグラムが表示される（図２）。
plt.title('Samples of b0')
plt.xlabel('b0')
plt.ylabel('frequency')
plt.hist(fit['b0'])
plt.show()
図２
図２のWindowを閉じると、次のスクリプトにより、パラメータb1の事後分布のヒストグラムが表示される（図３）。
ppos = np.mean(fit['b1']
> 0)
plt.title(('Samples of b1\n') +
('Prop(b1 > 0) = {0:.3f}'.format(ppos))) # npos/nb1)))
plt.xlabel('b1')
plt.ylabel('frequency')
plt.hist(fit['b1'])
plt.show()
図３
「b1>0」の頻度が
であるので、単回帰係数b1は統計的に０より大きいと判断できる。
図３のWindowを閉じると、パラメータsgmの事後分布のヒストグラムが表示される（図４）。
plt.title('Samples of sgm\n')
plt.xlabel('sgm')
plt.ylabel('frequency')
plt.hist(fit['sgm'])
plt.show()
図４
図４のWindowを閉じると、b0とb1の事後分布の散布図が、次のスクリプトにより表示される（図５）。
r, p = scst.pearsonr(fit['b0'],
fit['b1'])
print('r = ', r)
plt.xlabel('b0')
plt.ylabel('b1')
plt.title('r = {0:.3f}'.format(r))
plt.plot(fit['b0'], fit['b1'], 'b.' )
plt.show()
図５
強い負の相関
が見られる。
図５のWindowを閉じると、単回帰直線のグラフが表示される（図６）。
図６
橙色の細い実線は、事後分布からランダムにb0とb1を50組選び描画したものであり、次のスクリプトにより描いている。
b0s = fit['b0']
b1s = fit['b1']
rng = np.random.default_rng()
indx_ch = rng.choice(np.arange(len(b0s)), size=50,
replace=False)
min_x = np.min(X)
max_x = np.max(X)
for i in indx_ch:
plt.plot([min_x, max_x], [b0s[i] + b1s[i]*min_x,
b0s[i] + b1s[i]*max_x],
c='orange', alpha=0.3)
傾きの小さい回帰直線のｙ切片が小さい傾向にあることが認められるが、これは図５のb0とb1の散布図における相関係数が負であることに対応している。
青の実線は、事後分布の中央値をパラメータの推定値として回帰直線を描いたもので、次のスクリプトにより描いている。
b0_med = np.median(b0s)
b1_med = np.median(b1s)
plt.plot([min_x, max_x], [b0_med + b1_med*min_x,
b0_med + b1_med*max_x],
label=f'y={b0_med:.2f}-{-b1_med:.2f}*x')
Gelman, A., Hill, J., & Vehtari, A. (2021). Regression and other stories. では、事後分布の代表値として中央値が薦められている。
都市名とそのデータ値のプロットは、次のスクリプトにより描いている。
for nm, vx, vy in zip(ID, X, Y):
plt.text(vx + 0.1, vy + 0.1, nm,
rotation=25)
plt.plot(vx, vy, 'o', c='g')
図６のWindowを閉じると、スクリプトの実行終了である。
リスト１ MCMCサンプリングのPythonスクリプト（Regression.py；リスト２のStanスクリプトRegression.stanを使用）
import pystan
import arviz
as az
import pandas as pd
import matplotlib.pyplot
as plt
import scipy.stats
as scst
import numpy
as np
RawData = [ ['City', 'Temperature',
'Latitude'],
['Sapporo',
15,
43.1],
['Sendai',
20,
38.3],
['Niigata',
17,
37.9],
['Kanazawa',
18,
36.6],
['Tokyo', 24,
35.7],
['Osaka', 22,
34.7],
['Fukuoka',
22,
33.6],
['Kochi', 24,
33.6],
['Kagoshima', 23,
31.6],
['Naha',
26,
26.2]
]
RawData = np.array(RawData)
# 数値データは、文字列に変換されている
print(RawData)
ID = RawData[1:,0]
Y = [float(v) for v in RawData[1:,1]]
X = [float(v) for v in RawData[1:,2]]
N = len(ID)
for v in zip(ID, Y, X):
print(v)
Data = {'N': N, 'Y': Y, 'X': X}
sm = pystan.StanModel(file='Regression.stan')
fit = sm.sampling(data=Data,
control={'adapt_delta':0.999})
print(fit)
"""
Arvizによるトレースプロットの描画
"""
i_data = az.from_pystan(posterior
= fit, posterior_model = sm)
az.plot_trace(i_data)
plt.show()
"""
PndasのDataFrameによる統計量の計算
"""
d_frame = {'b0':fit['b0'], 'b1': fit['b1'], 'sgm': fit['sgm']}
for i in
range(10):
d_frame[f'mu.{i+1}'] = fit['mu'][:,i]
d_frame = pd.DataFrame(d_frame)
pd.set_option('display.max_columns',
None)
print('Statistics...')
print(d_frame.describe(percentiles
= [0.025, 0.25, 0.5, 0.75, 0.975]).T)
"""
事後分布の描画
"""
plt.title('Samples of b0')
plt.xlabel('b0')
plt.ylabel('frequency')
plt.hist(fit['b0'])
plt.show()
ppos = np.mean(fit['b1']
> 0)
plt.title(('Samples of b1\n') +
('Prop(b1 > 0) = {0:.3f}'.format(ppos))) # npos/nb1)))
plt.xlabel('b1')
plt.ylabel('frequency')
plt.hist(fit['b1'])
plt.show()
plt.title('Samples of sgm\n')
plt.xlabel('sgm')
plt.ylabel('frequency')
plt.hist(fit['sgm'])
plt.show()
r, p = scst.pearsonr(fit['b0'],
fit['b1'])
print('r = ', r)
plt.xlabel('b0')
plt.ylabel('b1')
plt.title('r = {0:.3f}'.format(r))
plt.plot(fit['b0'], fit['b1'], 'b.' )
plt.show()
"""
回帰直線の描画
"""
b0s = fit['b0']
b1s = fit['b1']
rng = np.random.default_rng()
indx_ch = rng.choice(np.arange(len(b0s)), size=50,
replace=False)
print(indx_ch)
"""
MCMCアンプルからランダムに50組を選んで描画
"""
min_x = np.min(X)
max_x = np.max(X)
for i in indx_ch:
plt.plot([min_x, max_x], [b0s[i] + b1s[i]*min_x,
b0s[i] + b1s[i]*max_x],
c='orange', alpha=0.3)
"""
回帰係数の事後分布の中央値を推定値とする単回帰直線
"""
b0_med = np.median(b0s)
b1_med = np.median(b1s)
plt.plot([min_x, max_x], [b0_med + b1_med*min_x,
b0_med + b1_med*max_x],
label=f'y={b0_med:.2f}-{-b1_med:.2f}*x')
"""
データのプロット
"""
for nm, vx, vy in zip(ID, X, Y):
plt.text(vx + 0.1, vy + 0.1, nm,
rotation=25)
plt.plot(vx, vy, 'o', c='g')
plt.xlabel('Latitude')
plt.ylabel('Temperature')
plt.legend()
plt.show()
リスト２ リスト１で用いられるStanスクリプト（Regression.stan）
data{
int N;
real Y[N];
real X[N];
}
parameters{
real b0;
real b1;
real<lower = 0> sgm;
}
transformed parameters{
real mu[N];
for (i
in 1:N)
mu[i] = b0 + b1 * X[i];
}
model{
Y ~ normal(mu, sgm);
}
PyStan 2の説明を、岡本安晴（2019）「いまさら聞けないPythonでデータ分析」丸善出版で行っている。