拙著「Pythonでデータ分析」の第14章「PyMCによる単回帰モデル分析」の更新
スクリプトをリスト１に示すように更新した。スクリプトファイルは、圧縮ファイルc14files.zipとして用意した。ダウンロード展開すれば自由に利用できる。
更新スクリプトでは、独立変数は中心化（平均値を引いて、平均値が０になるようにする）してから分析した。実行は、Python 3.10/Anaconda/Ubuntu/WSL/Windowsで行った。すなわち、WindowsにWSL (Windows Subsystem for Linux)をインストールするとUbuntuもインストールされるので、そのUbuntuにAnacondaをインストールしてPython 3.10仮想環境を用意した。Python 3.10仮想環境においてコマンド
conda install
pymc3
によりPyMC3をインストールした。詳しくは、このウェブサイトを参照されたい。
リスト１のスクリプトをPython 3.10仮想環境において実行した。
(py310) ・・・/c14files$ python
RegressionPyMC3.py
上記コマンドを実行するとスクリプトが実行され、MCMCサンプリング後、トレース図が表示される。
図１
図１のフォームを閉じると、図２のフォームが表示される。
図２
切片b0は、独立変数が中心化されているので、独立変数の平均値における回帰モデルの値である。
図２のフォームを閉じると、図３のフォームが表示される。
図３
図３のフォームを閉じると、図４のフォームが表示される。
図４
独立変数が中心化されているので、b0とb1のサンプルに相関はほとんど見られない。
図４のフォームを閉じると、図５のフォームが表示される。
図５
散布図に回帰直線が描かれている。
図５のフォームを閉じると、次のメッセージが表示されて、実行終了となる。
Results.txt
was saved.
出力ファイルResults.txtを開くと、以下のような内容である。
Mean of X = 35.130
Summary...
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
b0 21.119 0.726 19.742 22.552 0.013 0.009 3262.0 2326.0 1.0
b1 -0.677 0.174 -1.004 -0.342 0.003 0.002 2924.0 1987.0 1.0
sgm 2.196 0.701 1.170 3.494 0.017 0.013 1854.0 2026.0 1.0
b0: med = 21.109, (Q1, Q3) = (20.667, 21.557)
b1: med = -0.676, (Q1, Q3) = (-0.783, -0.571)
リスト１．単下位モデルのベイズ分析（RegressionPyMC3.py）
import matplotlib.pyplot as plt
import japanize_matplotlib
import seaborn as sb
import numpy as np
import scipy.stats as scst
import pymc3 as pm
import arviz as az
RawData = [ ['都市', '気温', '緯度'],
['札幌', 15, 43.1],
['仙台', 20, 38.3],
['新潟', 17, 37.9],
['金沢', 18, 36.6],
['東京', 24, 35.7],
['大阪', 22, 34.7],
['福岡', 22, 33.6],
['高知', 24, 33.6],
['鹿児島', 23, 31.6],
['那覇', 26, 26.2]
]
f = open('Results.txt', 'w')
RawData = np.array(RawData)
var_names = RawData[0]
ID = RawData[1:,0]
PreData = np.array(RawData[1:,1:],
dtype='float')
Y = PreData[:, 0]
X = PreData[:, 1]
for v in zip(ID, Y, X):
print(v)
N = len(X)
Xc = X - np.mean(X)
f.write(f'\nMean of X =
{np.mean(X):.3f}\n')
with pm.Model() as RegModel:
b0 = pm.Uniform('b0', lower
= -100.0, upper = 100.0)
b1 = pm.Uniform('b1', lower
= -100.0, upper = 100.0)
sgm = pm.Uniform('sgm',
lower = 0.0, upper = 100.0)
muY = b0 + b1 * Xc
Y_obs = pm.Normal('Y_obs',
mu = muY, sd = sgm, observed = Y)
trace = pm.sample()
az.plot_trace(trace)
plt.show()
summary_trace =
pm.summary(trace)
print(summary_trace)
f.write('\nSummary...\n{}'.format(summary_trace))
b0_trace = trace['b0']
b0_q1, b0_med, b0_q3 =
np.percentile(b0_trace, [25, 50, 75])
sb.kdeplot(b0_trace)
plt.title(f'b0: med = {b0_med:.3f}, (Q1, Q3) = ({b0_q1:.3f},
{b0_q3:.3f})')
plt.xlabel('b0')
plt.show()
f.write(f'\n\nb0: med = {b0_med:.3f}, (Q1, Q3) = ({b0_q1:.3f},
{b0_q3:.3f})\n')
b1_trace = trace['b1']
b1_q1, b1_med, b1_q3 =
np.percentile(b1_trace, [25, 50, 75])
sb.kdeplot(b1_trace)
plt.title(f'b1: med = {b1_med:.3f}, (Q1, Q3) = ({b1_q1:.3f},
{b1_q3:.3f})')
plt.xlabel('b1')
plt.show()
f.write(f'\nb1: med = {b1_med:.3f}, (Q1, Q3) = ({b1_q1:.3f},
{b1_q3:.3f})')
r, p =
scst.pearsonr(b0_trace, b1_trace)
plt.title('r =
{0:.3f}'.format(r))
plt.xlabel('b0')
plt.ylabel('b1')
plt.plot(b0_trace, b1_trace,
'b.')
plt.show()
for i in range(N):
plt.plot(Xc[i], Y[i], 'o', c='b')
plt.text(Xc[i], Y[i], ID[i])
min_Xc = np.min(Xc)
max_Xc = np.max(Xc)
plt.plot([min_Xc, max_Xc],
[b0_med + b1_med*min_Xc, b0_med + b1_med*max_Xc],
label = f'Y = {b0_med:.3f} + {b1_med:.3f} * Xc')
plt.legend()
plt.xlabel('Xc: ' +
var_names[2] + '(centered)')
plt.ylabel('Y: ' +
var_names[1])
plt.show()
f.close()
print('Results.txt was
saved.')