単回帰分析
下図のデータdata.csvの単回帰分析を行う。
図１
すなわち、データ
において
と
の関係を次式のモデル（１）で表すことを考える。
ここで、
は誤差（残差）を表し、平均が０の正規分布に従うものとする。
は、平均０，分散
の正規分布を表し、記号
は、左辺の確率変数が右辺の分布に従うことを表す。
単回帰モデルによるベイズ分析を行うためのStanスクリプトをリスト１に示す。
リスト１ 単回帰分析のStanスクリプト（simple_regression.stan）
data {
int
N;
real X[N];
real Y[N];
}
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）
import numpy
as np
import pandas as pd
import matplotlib.pyplot
as plt
import pystan
import seaborn
as sb
df = pd.read_csv('data.csv')
#print(df)
data = df.values
#print('data =\n', data)
X = data.T[1]
Y = data.T[2]
for x, y in zip(X, Y):
print(x, y)
sm = pystan.StanModel(file
= 'simple_regression.stan')
print('Compiling end.')
Data = {'N':len(X),
'X':X, 'Y':Y}
fit = sm.sampling(data
= Data, n_jobs = 1)
print(fit)
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()
実行すると、以下の出力を得る。
図２ aの事後分布
図３ bの事後分布
図４ (b,a)の散布図
図５ 回帰直線
中心化（Centering）
モデル（１）の定数が独立変数の平均値における値を示すようにするために、独立変数の中心化を行う。
この中心化されたZを独立変数として回帰分析を行うPythonスクリプトをリスト３に示す。
リスト３ 中心化された独立変数による単回帰分析（SRegression_centered.py）
import numpy
as np
import pandas as pd
import matplotlib.pyplot
as plt
import pystan
import seaborn
as sb
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()
sm = pystan.StanModel(file
= 'simple_regression.stan')
print('Compiling end.')
Data = {'N':len(Z),
'X':Z, 'Y':Y}
fit = sm.sampling(data
= Data, n_jobs = 1)
print(fit)
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, $\hat{a}$)')
+ ('\n$\hata$={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()
リスト３のスクリプトを実行すると、以下の出力を得る。
図６ aの事後分布
図７ bの事後分布
図８ (b,a)の散布図
図９ 回帰直線