単回帰分析／Python

単回帰分析

2025.11：CmdStanPy対応に改訂

StanのPythonインターフェースCmdStanPy対応とした。CmdStanPyの簡単な説明を、このウェブサイトに用意した。

下図のデータdata.csvの単回帰分析を行う。

図１

すなわち、データにおいてとの関係を次式のモデル（１）で表すことを考える。

ここで、は誤差（残差）を表し、平均が０の正規分布に従うものとする。

は、平均０，分散の正規分布を表し、記号は、左辺の確率変数が右辺の分布に従うことを表す。

単回帰モデルによるベイズ分析を行うためのStanスクリプトをリスト１に示す。

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

data {

int N;

array[N] real X;

array[N] real Y;

}

parameters {

real a;

real b;

real<lower = 0> sgm;

}

model {

a ~ normal(0, 100);

b ~ normal(0, 100);

sgm ~ exponential(0.001);

for (i in 1:N)

Y[i] ~ normal(a + b*X[i], sgm);

}

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

Stanについては、岡本安晴（2019）で説明している。

リスト１のStanスクリプトファイルを使って、図１のデータを分析するPythonスクリプトをリスト２に示す。

リスト２　単回帰分析のPythonスクリプト例（SimpleRegression.py）

from cmdstanpy import CmdStanModel

import numpy as np

import pandas as pd

import matplotlib.pyplot as plt

import seaborn as sb

import arviz as az

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

data = df.values

X = data.T[1]

Y = data.T[2]

for x, y in zip(X, Y):

print(x, y)

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

model = CmdStanModel(stan_file='simple_regression.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型に変換

a_hat = np.median(fit['a'])

b_hat = np.median(fit['b'])

sb.kdeplot(fit['a'])

plt.title('a: median = {0:.2f}'.format(a_hat), fontsize = 18)

plt.savefig('Fig_a.png')

plt.show()

sb.kdeplot(fit['b'])

plt.title('b: median = {0:.2f}'.format(b_hat), fontsize = 18)

plt.savefig('Fig_b.png')

plt.show()

plt.plot(fit['b'], fit['a'], 'o')

plt.xlabel('b', fontsize = 20)

plt.ylabel('a', fontsize = 20)

plt.savefig('Fig_scatter_ab.png')

plt.show()

plt.plot(X, Y, 'o', label = 'Data')

plt.plot([min(X), max(X)], [a_hat + b_hat * min(X), a_hat + b_hat * max(X)],

label = 'Regression Line')

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

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

plt.legend(fontsize = 14)

plt.savefig('FigSRegress.png')

plt.show()

リスト１のStanスクリプトファイル（simple_regression.stan）、リスト２のPythonスクリプトファイル（SimpleRegression.py）、入力データファイル（data.csv）を同じフォルダに置く。そのフォルダにカレントディレクトリ（フォルダ）を移して、CmdStanPyのインストールされている環境で、次のコマンドを実行する。CmdStanPyの簡単な解説は、このウェブサイトに用意した。

(stan) *****/filesSRegression$ python SimpleRegression.py

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

図２

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

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

図３　aの事後分布

パラメータaの事後分布のグラフである。

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

図４　bの事後分布

パラメータbの事後分布のグラフである。

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

図５　(b,a)の散布図

パラメータaとbのMCMCサンプル(b、a)の散布図である。強い相関が表れている。パラメータ間の相関を避ける方法として、中心化がある。

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

図６　回帰直線

データの散布図と回帰直線が描かれている。回帰直線の傾きが小さくなると、ｙ切片の値が大きなるという図５で表されていることが納得できる。

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

傾きbとy切片aの相関が無いようにする方法として、中心化がある。次に、中心化について説明する。

中心化（Centering）

モデル（１）の定数aが,、独立変数の平均値における値を示すようにするために、独立変数の中心化を行う。

この中心化されたZを独立変数として回帰分析を行うPythonスクリプトをリスト３に示す。

リスト３　中心化された独立変数による単回帰分析（SRegression_centered.py）

from cmdstanpy import CmdStanModel

import numpy as np

import pandas as pd

import matplotlib.pyplot as plt

import seaborn as sb

import arviz as az

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

data = df.values

X = data.T[1]

Y = data.T[2]

for x, y in zip(X, Y):

print(x, y)

Z = X - np.mean(X)

plt.plot(Z, Y, 'o')

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

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

plt.show()

Data = {'N':len(Z), 'X':Z, 'Y':Y}

model = CmdStanModel(stan_file='simple_regression.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型に変換

a_hat = np.median(fit['a'])

b_hat = np.median(fit['b'])

sb.kdeplot(fit['a'])

plt.title('a: median = {0:.2f}'.format(a_hat))

plt.savefig('FigZ_a.png')

plt.show()

sb.kdeplot(fit['b'])

plt.title('b: median = {0:.2f}'.format(b_hat))

plt.savefig('FigZ_b.png')

plt.show()

plt.plot(fit['b'], fit['a'], 'o')

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

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

plt.savefig('FigZ_ab.png')

plt.show()

plt.plot(Z, Y, 'o')

plt.plot([min(Z), max(Z)], [a_hat + b_hat * min(Z), a_hat + b_hat * max(Z)],

label = 'Regression Line')

plt.plot(0, a_hat, 's', c = 'r',

label = (r'(0, a_hat)') + ('\na_hat={0:.2f}'.format(a_hat)))

plt.xlabel(r'Z = X - $\bar{X}$', fontsize = 16)

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

plt.legend(fontsize = 14)

plt.savefig('FigZ_regress.png')

plt.show()

リスト３のPythonスクリプトファイル（SRegression_centered.y）、リスト１のStanスクリプトファイル（simple_regression.stan）、入力データファイル（data.csv：図１）を同じフォルダに置く。そのフォルダに、カレントディレクトリ（フォルダ）を移し、CmdStanPyのインストールされた環境において、次のコマンドを実行する。CmdStanPyの簡単な説明を、このウェブサイトに用意した。