Home

重回帰モデルのベイズ分析

 

拙著「Pythonでデータ分析」の「第５章　重回帰分析」のベイズ分析用スクリプトを用意した。

重回帰分析を、変数をすべて標準化して行い、事前分布として、係数の事前分布は平均が０，標準偏差が2.5の正規分布、誤差の事前分布は平均が１の指数分布を採用した（soft constraint; Gelmanら、2021、p.124）。Stanスクリプトをリスト１に示す。スクリプトファイルなどは圧縮ファイルmlregfiles.zipにまとめた。ダウンロード・展開して自由に利用されたい。

 

 

リスト１　重回帰モデルのベイズ分析Stanスクリプト（標準化変数：MLReg.stan）

 

data {

    int N;

    int M;

    vector[N] y;

    array[N] vector[M] X;

}

parameters {

    real b0;

    vector[M] b;

    real<lower = 0.0> sgm;

}

model {

    b0 ~ normal(0.0, 2.5);

    b ~ normal(0.0, 2.5);

    sgm ~ exponential(1/1.0);

    for (i in 1:N) {

        y[i] ~ normal(b0 + dot_product(b,X[i]), sgm);

    }

}

 

 

リスト１のStanスクリプト（PyStan 3.7.0対応）を用いるベイズ分析Pythonスクリプトをリスト２に示す。リスト２のスクリプトを、WindowsのWSL (Windows Subsystem for Linux)のUbuntuにインストールしたAnacondaのPython 3.11仮想環境上で実行した。PyStan 3は、Python 3.11仮想環境において、コマンド

 

pip install pystan

 

を実行してインストールした。Pythonスクリプトでのimportは

 

import stan

 

で行う。

WindowsでのPyStan 3の利用については、このウェブサイトで説明した。

リスト２のスクリプトを実行すると、入力データファイル名を聞いてくるので、設定する。入力データファイルは、図１の形式でExcelファイルとして用意する。

図１

 

1行目に変数名を設定して、データ値は２行目以降に書く。１列目は、データの識別文字列、２列目に従属変数、３列目以降に独立変数を置く。変数値は、スクリプトの実行時に標準得点化される。

次の例では、図１のファイル予想最高気温xlsxが設定されている。

 

(py311) ・・・/mlregfiles$ python MLRegr.py

入力データファイル(*.xlsx) = 予想最高気温.xlsx

 

入力データファイル名を設定すると、計算が始まる。

MCMCサンプリングが終了すると、図１のトレース図が表示される。

図２

 

表示グラフは、保存ボタンをクリックすると保存できる。Ubuntu/WSLで保存ボタンをクリックするとUbuntuでのフォルダが表示される。Windowsのフォルダに保存するときは、/mntフォルダ（ディレクトリ）の下にWindowsのディレクトリ（フォルダ）があるので、/mntディレクトリに移れば、Windowsのディレクトリが表示される。Windowsのディレクトリは、例えばC:ドライブは小文字１文字cで表されるので注意。

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

図３

 

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

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

図４

 

標準化回帰係数（standardized regression coefficient）は、標準化独立変数の値が１異なるときの回帰モデルの値の違いを表す（Gelmanら、2021、p.133；岡本、2021、p.55）。

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

図５

 

図５のフォームを閉じると、図６の散布図が表示される。

図６

 

図６のフォームを閉じると、スクリプトの実行終了である。端末には次のメッセージが表示される。

 

Results.txt was saved.

 

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

 

 

入力データファイル = 予想最高気温.xlsx

 

y_mean = 20.750,   y_sd = 6.193

 

b1(緯度):   mean = 37.458,   sd = 5.784

 

b2(西経):   mean = 95.342,   sd = 17.235

 

 

       mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_bulk  ess_tail  r_hat

b0    0.002  0.221  -0.402    0.434      0.004    0.004    3770.0    2294.0    1.0

b[0] -0.805  0.233  -1.238   -0.388      0.004    0.003    3404.0    2184.0    1.0

b[1]  0.268  0.221  -0.121    0.692      0.004    0.003    3547.0    2278.0    1.0

sgm   0.743  0.197   0.444    1.116      0.005    0.003    1987.0    1892.0    1.0

 

b0_med = 0.003,    b0_mad_sd = 0.203

 

b1(緯度):   med = -0.802,  mad_sd = 0.209

 

b2(西経):   med = 0.271,  mad_sd = 0.204

 

 

引用文献

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

岡本安晴（2022）．データを使いこなすための統計入門．22世紀アート（電子書籍）

 

 

 

リスト２　ベイズ分析Pythonスクリプト（MLRegr.py）

 

import pandas as pd

import numpy as np

import scipy.stats as ss

import matplotlib.pyplot as plt

import japanize_matplotlib

import seaborn as sb

import stan

import arviz as az

 

"""

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

"""

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

xlsx = pd.ExcelFile(in_file_name)

data_xlsx = pd.read_excel(xlsx)

 

"""

            出力ファイルの設定

"""

f = open('Results.txt', 'w')

f.write('\n入力データファイル = ' + in_file_name)

 

keys = data_xlsx.keys()         #     １行目の変数名の読み込み

print('Keys...\n', keys)

case_name = keys[0]

 

y_name = keys[1]

x_names = keys[2:]

print(y_name)

print(x_names)

Data = data_xlsx.values

print(Data)

case_id = Data[:,0]

y = Data[:,1]

X = Data[:,2:]

 

print(case_id)

print(y)

print(X)

y_mean = np.mean(y)

y_sd = np.std(y)

f.write(f'\n\ny_mean = {y_mean:.3f},   y_sd = {y_sd:.3f}\n')

X = np.array(X, dtype='float')

x_means = np.mean(X, axis=0)

print(x_means)

x_sds = np.std(X, axis=0)

print(x_sds)

for j in range(len(x_names)):

    f.write(f'\nb{j+1}({x_names[j]}):   ')

    f.write(f'mean = {x_means[j]:.3f},   sd = {x_sds[j]:.3f}\n')

 

y = ss.zscore(np.array(y, dtype ='float'))

X = ss.zscore(np.array(X, dtype = 'float'))

print(y)

print(X)

 

N, M = np.shape(X)

print(N, M)

 

with open('MLReg.stan', 'r') as f_stan:

    sm = stan.build(f_stan.read(), data = {'N':N, 'M':M, 'y':y, 'X':X})

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\n')

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

 

az.plot_trace(i_data)

plt.show()

 

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

 

b0_med, b0_mad_sd = med_mad_sd(d_frm['b0'])

 

print(f'\nb0_med = {b0_med:.3f},    b0_mad_sd = {b0_mad_sd:.3f}')

f.write(f'\n\nb0_med = {b0_med:.3f},    b0_mad_sd = {b0_mad_sd:.3f}\n')

sb.kdeplot(d_frm['b0'])

plt.title('Posterior distribution for b0')

plt.show()

 

b_med = np.empty(M)

b_mad_sd = np.empty(M)

for j in range(M):

    b_med[j], b_mad_sd[j] = med_mad_sd(d_frm[f'b.{j+1}'])

    print(x_names[j])

    print(b_med[j], b_mad_sd[j])

    f.write(f'\nb{j+1}({x_names[j]}):' +

            f'   med = {b_med[j]:.3f},  mad_sd = {b_mad_sd[j]:.3f}\n')

    sb.kdeplot(d_frm[f'b.{j+1}'])

    plt.title(f'Posterir distribution for b{j+1}: ' + x_names[j])

    plt.show()

 

y_hat = np.empty(N)

for i in range(N):

    y_hat[i] = b0_med + np.dot(b_med, X[i])

 

for i in range(N):

    plt.plot(y_hat[i], y[i], 'o', c='b')

    plt.text(y_hat[i], y[i], case_id[i])

plt.xlabel('$\hat{y}$(model)')

plt.ylabel('$y$(data)')

x_min = np.min(y_hat)

x_max = np.max(y_hat)

plt.plot([x_min, x_max], [x_min, x_max], label ='y(data) = $\hat{y}$(model)')

plt.legend()

plt.show()

 

print('Results.txt was saved.')

 

 

 

Home