重回帰モデル／ベイズ分析／Python

重回帰分析

2025.11：CmdStanPy対応に改訂

ベイズ分析のためのソフトStanへのPythonインターフェースCmdStanPy対応に改訂した。CmdStanPyの簡単な説明を、このウェブサイトに用意した。

重回帰分析として独立変数が２個の場合について考える。

独立変数がとの２個である仮想データ例を図１に示す。

図１

図１に対する重回帰モデルを以下のように設定する。

パラメータの事前分布は、平らな事前分布として標準偏差の大きい正規分布と指数分布を以下のように設定する。

　上記の確率モデルによるベイズ分析を行うためにStanスクリプトをリスト１のように用意した。

スクリプトファイルとデータファイルは、圧縮ファイルfiles_mregress.zipにまとめた。ダウンロード・展開すれば自由に利用できる。

リスト１　重回帰分析のStanスクリプト（mv_regress.stan）

data {

int N;

array[N] real Y;

array[N] real X1;

array[N] real X2;

}

parameters {

real b0;

real b1;

real b2;

real<lower = 0.0> sgm;

}

model {

b0 ~ normal(0, 1000);

b1 ~ normal(0, 1000);

b2 ~ normal(0, 1000);

sgm ~ exponential(0.01);

for (i in 1:N)

Y[i] ~ normal(b0 + b1 * X1[i] + b2 * X2[i], sgm);

}

リスト１のStanスクリプトを使って図１のデータの重回帰分析を行うPythonスクリプトをリスト２のように用意した。

リスト２　重回帰分析を行うPythonスクリプト（MRegress.py）

from cmdstanpy import CmdStanModel

import numpy as np

import pandas as pd

import matplotlib.pyplot as plt

import arviz as az

from matplotlib import cm

import seaborn as sb

df = pd.read_csv('MVData.csv')

data = df.values

print('data =')

Y = data.T[1]

X1 = data.T[2]

X2 = data.T[3]

for y, x1, x2 in zip(Y, X1, X2):

print(y, x1, x2)

Data = {'N':len(Y), 'Y':Y, 'X1':X1, 'X2':X2}

model = CmdStanModel(stan_file='mv_regress.stan')

fit = model.sample(data=Data)

print(fit.summary())

inf_data = az.from_cmdstanpy(fit) # Arviz DataFrame型に変換

az.plot_trace(inf_data) # トレースプロット

plt.tight_layout()

plt.show()

fit = fit.draws_pd() # Pandas DataFrame型に変換

b0_hat = np.median(fit['b0'])

b1_hat = np.median(fit['b1'])

b2_hat = np.median(fit['b2'])

sb.kdeplot(fit['b0'])

plt.title('b0: median = {0:.2f}'.format(b0_hat))

plt.savefig('Figb0.png')

plt.show()

sb.kdeplot(fit['b1'])

plt.title('b1: median = {0:.2f}'.format(b1_hat))

plt.savefig('Figb1.png')

plt.show()

sb.kdeplot(fit['b2'])

plt.title('b2: median = {0:.2f}'.format(b2_hat))

plt.savefig('Figb2.png')

plt.show()

fig = plt.figure()

ax = fig.add_subplot(projection='3d')

ax.scatter(X1, X2, Y, color='blue')

Xs = np.linspace(min(X1), max(X1), 100)

Ys = np.linspace(min(X2), max(X2), 100)

Xg, Yg = np.meshgrid(Xs, Ys)

Z = np.empty((100, 100))

for iy in range(100):

for ix in range(100):

Z[iy][ix] = b0_hat + b1_hat * Xs[ix] + b2_hat * Ys[iy]

ax.plot_surface(Xg, Yg, Z, cmap = cm.viridis, alpha = 0.3)

ax.set_xlabel('X1')

ax.set_ylabel('X2')

ax.set_zlabel('Y')

ax.set_title('Multivariate Regression')

for i in range(len(Y)):

plt.plot([X1[i],X1[i]], [X2[i],X2[i]],

[b0_hat + b1_hat*X1[i] + b2_hat*X2[i], Y[i]], c = 'm')

plt.plot([0], [0], [b0_hat], 'D', color = 'r',

label = '(0, 0, b0)' + '\n' + f'b0={b0_hat:.2f}')

plt.legend()

plt.title('You can rotate this graph by dragging')

plt.savefig('FigRegress.png')

plt.show()

リスト１のスクリプトファイル（mv_regress.stan）、リスト２のスクリプトファイル（MRegress.py）、および入力デーらファイル（例えば図１のファイルMVData.csv）を同じフォルダに置く。そのフォルダにカレントディレクトリ（フォルダ）を移し、CmdStanPyのインストールされた環境で次のコマンドを実行する。

(stan) *****/files_mregress$ python MRegress.py

上のコマンドを実行すると、ファイルが読み込まれ、Stanスクリプトのコンパイルが始まる。コンパイルに多少時間が掛かる。コンパイルが終了すると、MCMCサンプリングが行われ、サンプリング後、トレース図（図２）が表示される。

図２

左側に事後分布のグラフ、右側にサンプリングのトレース図が表示されている。

図２のWindowを閉じると、図３のグラフが表示される。

図３　b0の事後分布

b0の事後分布のグラフである。点推定値として中央値が示されている。Gelmanら（2021）は、代表値として、平均値より安定した中央値を勧めている。

図３のWindowを閉じると、図４のグラフが表示される。

図４　b1の事後分布

b1の事後分布のグラフである。

図４のWindowを閉じると、図５のグラフが表示される。

図５　b2の事後分布

b2の事後分布のグラフである。

図５のWindowを閉じると、図６の３Dグラフが表示される。

図６　回帰平面

回帰平面のグラフは、グラフ上をドラッグすれば回転できる。

図６では、重回帰モデルの定数項b0は点における値であるので、データ点の集まりから離れた位置にが赤色の小菱形としてプロットされている。

独立変数を中心化すると、定数項は独立変数の平均値における重回帰平面の値となる。中心化については、次に扱う。

図６のWindowを閉じると、スクリプトの実行終了である。

中心化

偏回帰係数b1、b2の値は、それぞれの独立変数の値が１だけ異なるときの回帰平面の値の違いを表す。これを独立変数の値が標準偏差だけ異なるときの回帰平面の値の違いを表すようにするときは、それぞれの独立変数の値をその標準偏差で割ったものを用いればよい。

変数の中心化と標準偏差を１に揃えることは変量の標準化と呼ばれているが、標準化後の変量をで表すと、次式で与えられる。

ここで、は平均値、は標準偏差を表す。

独立変数の標準化を行い、重回帰分析を行うPythonスクリプトをリスト３に示す。

リスト３　独立変数の標準化による重回帰分析（MRegress_Z.py）

from cmdstanpy import CmdStanModel

import numpy as np

import pandas as pd

import matplotlib.pyplot as plt

import arviz as az

import seaborn as sb

from matplotlib import cm

df = pd.read_csv('MVData.csv')

data = df.values

Y = data.T[1]

X1 = data.T[2]

X2 = data.T[3]

for y, x1, x2 in zip(Y, X1, X2):

print(y, x1, x2)

# 標準化

Z1 = (X1 - np.mean(X1)) / np.std(X1)

Z2 = (X2 - np.mean(X2)) / np.std(X2)

print('X1: mean = {0:.3f} std = {1:.3f}'.format(np.mean(X1), np.std(X1)))

print('X2: mean = {0:.3f} std = {1:.3f}'.format(np.mean(X2), np.std(X2)))

Data = {'N':len(Y), 'Y':Y, 'X1':Z1, 'X2':Z2}

model = CmdStanModel(stan_file='mv_regress.stan')

fit = model.sample(data=Data)

print(fit.summary())

inf_data = az.from_cmdstanpy(fit) # Arviz DataFrame型に変換

az.plot_trace(inf_data) # トレースプロット

plt.tight_layout()

plt.show()

fit = fit.draws_pd() # Pandas DataFrame型に変換

b0_hat = np.median(fit['b0'])

b1_hat = np.median(fit['b1'])

b2_hat = np.median(fit['b2'])

sb.kdeplot(fit['b0'])

plt.title('b0: median = {0:.2f}'.format(b0_hat))

plt.savefig('Figb0_c.png')

plt.show()

sb.kdeplot(fit['b1'])

plt.title('b1: median = {0:.2f}'.format(b1_hat))

plt.savefig('Figb1_c.png')

plt.show()

sb.kdeplot(fit['b2'])

plt.title('b2: median = {0:.2f}'.format(b2_hat))

plt.savefig('Figb2_c.png')

plt.show()

fig = plt.figure()

ax = fig.add_subplot(projection='3d')

ax.scatter(Z1, Z2, Y, color='blue')

Xs = np.linspace(min(Z1), max(Z1), 100)

Ys = np.linspace(min(Z2), max(Z2), 100)

Xg, Yg = np.meshgrid(Xs, Ys)

Z = np.empty((100, 100))

for iy in range(100):

for ix in range(100):

Z[iy][ix] = b0_hat + b1_hat * Xs[ix] + b2_hat * Ys[iy]

ax.plot_surface(Xg, Yg, Z, cmap = cm.viridis, alpha = 0.3)

ax.set_xlabel('Z1')

ax.set_ylabel('Z2')

ax.set_zlabel('Y')

ax.set_title('Multivariate Regression')

for i in range(len(Y)):

plt.plot([Z1[i],Z1[i]], [Z2[i],Z2[i]],

[b0_hat + b1_hat*Z1[i] + b2_hat*Z2[i], Y[i]], c = 'm')

plt.plot([0], [0], [b0_hat], 'D', color = 'r',

label = '(0, 0, b0)' + '\n' + f'b0={b0_hat:.2f}')

plt.legend()

plt.title('You can rotate this graph by dragging')

plt.savefig('FigRegress_c.png')

plt.show()

リスト３のPythonスクリプトファイル（MRegress_Z.py）、リスト１のStanスクリプトファイル（mv_regress.stan）、および入力データファイル（図１：MVData.csv）を同じフォルダに置く。そのフォルダにカレントディレクトリ（フォルダ）を移し、CmdStanPyのインストールされた環境において、次のコマンドを実行する。

(stan) *****/files_mregress$ python MRegress_Z.py

上のコマンドを実行すると、ファイルが読み込まれ、Stanスクリプトのコンパイルが始まる。コンパイル後、MCMCサンプリングが行われ、サンプリング後、トレース図が表示される（図７）。