回帰直線の分布
ブートストラッピング
Galtonのデータの母と娘の身長の回帰直線の分布をブートストラッピングによって求めることがSpiegelhalter(2019)で行われている。Spiegelhalter(2019), p. xviに記載されているウェブサイトからダウンロードしたファイルを参考にして、Pythonスクリプトを書いてみた。入力データファイル07-4-galton-x.csvは、ダウンロードしたファイルに含まれていたものである。
作成したPythonスクリプトをリスト１に示す。リスト１のスクリプトを実行すると、図１のグラフが表示される。
図１
データの散布図と回帰直線（黒の実線）が描かれている。散布図は、データの重なりを避けるため、ジッターが加えられている。
回帰直線の係数は、次のコードで求めている。
y = daughter_data[:,4] # Daughter's height
x = daughter_data[:,2] # Mother's height
n = len(y)
X = np.concatenate((np.ones((n,1)), np.reshape(x,
(n,1))), axis=1)
X = np.array(X,
dtype=float)
y = np.array(y,
dtype=float)
B0 = np.linalg.pinv(X)
@ np.reshape(y, (n,1)) # Regression coefficients
上のコードにおいて、pinvはムーア・ペンローズ逆行列を与える関数であるが、回帰分析の係数がムーア・ペンローズ逆行列で与えられることは、岡本（2009）、p.26で説明している。
ブートストラッピングは、ブートストラッピング用のデータXyを次のコード
Xy = np.concatenate((X,
np.reshape(y, (n,1))), axis=1)
で作成して、サンプリングを次のコードで行うという方法である。
smpl = rng.choice(Xy, n)
# サンプリング
各サンプリングに対して、次のコードで回帰直線の計算と描画、および回帰直線の値の保存を行っている。
B = np.linalg.pinv(smpl[:,:2]) @ smpl[:,2]
plt.plot([x_min, x_max], [reg(x_min,
B), reg(x_max, B)],
lw=1, c='orange', alpha=0.3)
# 回帰直線のxp100における値をregp100に格納
regp100.append(reg(xp100,
B))
保存された回帰直線の値の2.5％および97.5％値を次のコード
pcntp = np.percentile(regp100,
[2.5, 97.5], axis=0)
で求め、次のコードで描画している。
plt.plot(xp100, pcntp[1],
c='aqua', ls='--', label='97.5%')
plt.plot(xp100, pcntp[0],
c='lime', ls='-.', label='2.5%')
最後に、入力データに対する回帰直線と散布図を次のコードで描いている。
#
# 入力データに対する回帰直線の描画
#
plt.plot([x_min, x_max], [reg(x_min,
B0), reg(x_max, B0)],
lw=3, c='k', label='Regression line')
#
# データの散布図
# 重なりを避けるためにジッターを加えてある
#
plt.plot(x + rng.uniform(-1,
1, n), y + rng.uniform(-1, 1, n),
'o', c='k', lw=0, fillstyle='none')
回帰直線と散布図の描画を最後に行うのは、他の描画で覆われて隠されるのを避けるためである。
import pandas as pd
import matplotlib.pyplot
as plt
import seaborn
as sb
import numpy
as np
# ダウンロードしたファイルの読み込み
d_csv = pd.read_csv('07-4-galton-x.csv')
v_data = d_csv.values
daughter_data = v_data[v_data[:,3]=='F']
y = daughter_data[:,4] # Daughter's height
x = daughter_data[:,2] # Mother's height
n = len(y)
X = np.concatenate((np.ones((n,1)), np.reshape(x,
(n,1))), axis=1)
X = np.array(X,
dtype=float)
y = np.array(y,
dtype=float)
B0 = np.linalg.pinv(X)
@ np.reshape(y, (n,1)) # Regression coefficients
x_min = np.min(x)
x_max = np.max(x)
def reg(v,
B):
# regression
line
return B[0] + B[1] * v
rng = np.random.default_rng()
#
# Xとyをまとめてブートストラッピングに於けるサンプリングの対象とする
#
Xy = np.concatenate((X,
np.reshape(y, (n,1))), axis=1)
# x座標上に100個のポイントをとる
xp100 = np.linspace(x_min, x_max, 100)
#
# 各ブートストラッピングにおけるサンプルに対する
# 回帰直線のxp100における値をリストとしてregp100に
# 格納する
#
regp100 = []
#
ブートストラッピング
for i in
range(1000):
smpl
= rng.choice(Xy, n) # サンプリング
#
# サンプルsmplに対する回帰直線と、直線の描画
#
B = np.linalg.pinv(smpl[:,:2]) @ smpl[:,2]
plt.plot([x_min, x_max], [reg(x_min, B), reg(x_max, B)],
lw=1, c='orange', alpha=0.3)
# 回帰直線のxp100における値をregp100に格納
regp100.append(reg(xp100, B))
#
# ブートストラッピングにおける回帰直線の値の
# 2.5%点、97.5%点の値を取り出す。
#
pcntp = np.percentile(regp100,
[2.5, 97.5], axis=0)
#
# 2.5%点、97.5%点を結ぶ曲線の描画
#
plt.plot(xp100, pcntp[1],
c='aqua', ls='--', label='97.5%')
plt.plot(xp100, pcntp[0],
c='lime', ls='-.', label='2.5%')
#
# 入力データに対する回帰直線の描画
#
plt.plot([x_min, x_max], [reg(x_min,
B0), reg(x_max, B0)],
lw=3, c='k', label='Regression line')
#
# データの散布図
# 重なりを避けるためにジッターを加えてある
#
plt.plot(x + rng.uniform(-1,
1, n), y + rng.uniform(-1, 1, n),
'o', c='k', lw=0, fillstyle='none')
plt.xlabel("Mother's height")
plt.ylabel("Daughter's height")
plt.legend(loc='lower
right')
plt.title('Regression Lines by Bootstrapping')
plt.tight_layout()
plt.show()
岡本安晴（2009） 統計学を学ぶための数学入門［下］、培風館
岡本安晴（2019） いまさら聞けないPythonでデータ分析．丸善出版
Spiegelhalter, D. (2019). The art of
statistics: Learning from data. Pelican