Home

拙著「Pythonでデータ分析」の第11章「Pythonによる単回帰モデル分析」の更新

Gelmanら（2021）は、回帰モデルのデフォルト事前分布としてsoft constraintを説明している（p.124）。変数が標準化されている場合、係数の事前分布は平均が０，標準偏差が2.5の正規分布、誤差の事前分布は平均が１の指数分布が勧められている。

変数が標準化されているときのStanスクリプトをリスト１のように用意した。

リスト１　単回帰分析のStanスクリプト（PyStan 3.7.0対応）

data{

int N;

vector[N] Y;

vector[N] X;

}

parameters{

real b0;

real b1;

real<lower = 0> sgm;

}

transformed parameters{

vector[N] mu;

mu = b0 + b1 * X;

}

model{

b0 ~ normal(0.0, 2.5);

b1 ~ normal(0.0, 2.5);

sgm ~ exponential(1/1.0);

Y ~ normal(mu, sgm);

}

リスト１のStanスクリプトを用いて単回帰分析を行うPythonスクリプトをリスト２のように用意した。これらのスクリプトファイルなどは、圧縮ファイルc11files.zipにまとめている。ダウンロード・展開して自由に利用されたい。

リスト２のスクリプトでは、入力データはExcelファイルの形式のものを読み込む。入力ファイル例を図１に示す。

図１

リスト２の実行は、Windows上のUbuntu、すなわちWSL (Windows Subsystem for Linux)で行った。WSL/WindowsのUbuntuにインストールされたAnacondaのPython 3.11仮想環境においてPyStan 3をコマンド

pip install pystan

によりインストールした。その後、Python 3.11の仮想環境においてリスト２を実行した。

リスト２のスクリプトファイルs_regress.pyを実行すると、入力ファイル名を聞いてくる。

(py311) ・・・/c11files$ python s_regress.py

入力データファイル(*.xlsx) = 気温と緯度.xlsx

上の例では、図１のExcelファイルを指定している。

入力データファイルの指定を行うと、スクリプトの実行が進み、Stanスクリプトのビルド、MCMCサンプリングが行われる。MCMCサンプリング後、まずパラメータb0の事後分布が表示される（図２）。

図２

事後分布の代表値として中央値（med）、標準偏差に相当するものとしてmad sdが表示されている。Gelmanら（2021）は、平均値より安定している中央値とそれに基づく統計量を薦めている（ｐ.73）。モードは、データと事前分布に総合的に適合すると論じている（p.107）が、ここではGelmanらで実際に採用されている中央値に基づく統計量を採用した。。

図２のグラフは、保存ボタンをクリックすると保存されるが、保存先のフォルダ名はWSLから見た保存先である。Windowsから見た保存先は、/mntディレクトリに移ってから選ぶことに注意。

図２のフォームを閉じると、b1の事後分布のグラフが表示される（図３）。

図３

図３のフォームを閉じると、σの事後分布のグラフが表示される（図４）。

図４

図４のフォームを閉じると、パラメータb0とb1の散布図が表示される（図５）。

図５

図５のフォームを閉じると、データの散布図と回帰直線のグラフが表示される（図６）。

図６

回帰直線は、中央値を推定値として求めたものを用いて計算されている。

図６のフォームを閉じると、出力ファイルResults.txtが保存されたというメッセージが出力されて、実行終了となる。

Results.txt was saved.

(py311) ・・・/c11files$

出力ファイルを開くと以下のようになっている。

X: mean = 35.13 sd = 4.252

Y: mean = 21.10 sd = 3.330

mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat

b0 -0.000 0.206 -0.390 0.384 0.004 0.004 3039.0 2153.0 1.0

b1 -0.865 0.204 -1.236 -0.465 0.004 0.003 2903.0 2103.0 1.0

sgm 0.623 0.175 0.359 0.949 0.004 0.003 2158.0 2485.0 1.0

mu[0] -1.622 0.432 -2.476 -0.820 0.008 0.006 3037.0 2023.0 1.0

mu[1] -0.645 0.255 -1.135 -0.181 0.005 0.003 3115.0 2222.0 1.0

mu[2] -0.564 0.244 -1.044 -0.124 0.004 0.003 3113.0 2189.0 1.0

mu[3] -0.299 0.217 -0.699 0.119 0.004 0.003 3094.0 2193.0 1.0

mu[4] -0.116 0.208 -0.503 0.279 0.004 0.003 3059.0 2132.0 1.0

mu[5] 0.087 0.208 -0.310 0.472 0.004 0.004 3020.0 2151.0 1.0

mu[6] 0.311 0.220 -0.108 0.726 0.004 0.003 2976.0 2139.0 1.0

mu[7] 0.311 0.220 -0.108 0.726 0.004 0.003 2976.0 2139.0 1.0

mu[8] 0.718 0.269 0.221 1.229 0.005 0.004 2925.0 2245.0 1.0

mu[9] 1.816 0.478 0.932 2.732 0.009 0.007 2893.0 2221.0 1.0

b0: med.= 0.001, mad sd = 0.191

b1: med.= -0.865, mad sd = 0.190

sgm: med.= 0.592, mad sd = 0.148

引用文献

Gelman, A., Hill, J., & Vehtari, A. (2021). Regression and other stories. Cambridge University Press.

リスト２　リスト１のStanスクリプト（Regression.stan）を用いる単回帰分析Pythonスクリプト（s_regress.py）

import stan

import pickle

import pandas as pd

import numpy as np

import arviz as az

import matplotlib.pyplot as plt

import scipy.stats as scst

import seaborn as sb

import japanize_matplotlib

"""

Excelファイルからの読み込み

"""

in_file_name = input('入力データファイル(*.xlsx) = ')

xlsx = pd.ExcelFile(in_file_name)

data_xlsx = pd.read_excel(xlsx)

ID = data_xlsx.values[:,0]

Y = data_xlsx.values[:,1]

X = data_xlsx.values[:,2]

N = len(Y)

for v in zip(ID, Y, X):

print(v)

f = open('Results.txt', 'w') # 出力用ファイル

mean_X = np.mean(X)

sd_X = np.std(X)

print(f'X: mean = {mean_X:.2f} sd = {sd_X:.3f}')

f.write(f'X: mean = {mean_X:.2f} sd = {sd_X:.3f}\n')

X = (X - mean_X) / sd_X

mean_Y = np.mean(Y)

sd_Y = np.std(Y)

print(f'Y: mean = {mean_Y:.2f} sd = {sd_Y:.3f}')

f.write(f'\nY: mean = {mean_Y:.2f} sd = {sd_Y:.3f}\n')

Y = (Y - mean_Y) / sd_Y

Data = {'N': N, 'Y': Y, 'X': X}

#

# ビルドとサンプリング

#

with open('Regression.stan', 'r') as g:

sm = stan.build(g.read(), data = Data)

fit = sm.sample()

d_frm = fit.to_frame()

i_data = az.from_pystan(posterior = fit, posterior_model = sm)

print(az.summary(i_data))

f.write('\n')

f.write(az.summary(i_data).__str__())

def med_mad_sd(x):

""" 中央値とmad sdの計算

Gelman et al. (2021), p.73

"""

med = np.median(x)

mad = np.median(np.abs(x - med))

mad_sd = 1.483 * mad

return med, mad_sd

med_b0, mad_sd_b0 = med_mad_sd(d_frm['b0'])

med_b1, mad_sd_b1 = med_mad_sd(d_frm['b1'])

med_sgm, mad_sd_sgm = med_mad_sd(d_frm['sgm'])

f.write('\n\n')

f.write(f'b0: med.= {med_b0:.3f}, mad sd = {mad_sd_b0:.3f}\n')

f.write(f'\nb1: med.= {med_b1:.3f}, mad sd = {mad_sd_b1:.3f}\n')

f.write(f'\nsgm: med.= {med_sgm:.3f}, mad sd = {mad_sd_sgm:.3f}\n')

sb.kdeplot(d_frm['b0'])

plt.xlabel('b0', fontsize = 16)

plt.title(f'b0: med. = {med_b0:.3f}, mad sd = {mad_sd_b0:.3f}', fontsize = 16)

plt.show()

sb.kdeplot(d_frm['b1'])

plt.xlabel('b1', fontsize = 16)

plt.title(f'b1: med. = {med_b1:.3f}, mad sd = {mad_sd_b1:.3f}', fontsize = 16)

plt.show()

sb.kdeplot(d_frm['sgm'])

plt.xlabel('$\sigma$', fontsize = 14)

plt.title(f'$\sigma$: med. = {med_sgm:.3f}, mad sd = {mad_sd_sgm:.3f}',

fontsize = 16)

plt.show()

plt.plot(d_frm['b0'], d_frm['b1'], 'o', alpha = 0.3)

plt.xlabel('b0', fontsize = 16)

plt.ylabel('b1', fontsize = 16)

plt.tight_layout()

plt.show()

#

# 回帰直線と散布図

#

for i in range(len(ID)):

plt.plot(X[i], Y[i], 'o', c='b')

plt.text(X[i], Y[i], ID[i])

min_x = np.min(X)

max_x = np.max(X)

plt.plot([min_x, max_x], [med_b0 + med_b1*min_x, med_b0 + med_b1*max_x],

label = f'Y = {med_b0:.3f} + {med_b1:.3f} * X')

plt.xlabel('X', fontsize = 14)

plt.ylabel('Y', fontsize = 14)

plt.legend()

plt.tight_layout()

plt.show()

f.close()

print('Results.txt was saved.')

Home