はじめに

　『Pythonで学ぶ空間データサイエンス入門』の独学ノートです。本の内容から寄り道してアレコレ考えます。
　本を読んだ上で補助的に読んでください。　

　この記事では、線形回帰モデルの最小二乗法について、図を使って解説します。

【前の内容】

www.anarchive-beta.com

【他の内容】

www.anarchive-beta.com

【今回の内容】

はじめに
3.2.1-3 線形回帰モデルの最小二乗法の可視化：1次元の場合：
- 最小二乗法の図
- データ数の影響
参考文献
おわりに

3.2.1-3 線形回帰モデルの最小二乗法の可視化：1次元の場合：

　線形回帰モデル(linear regression model)に対する最小二乗法(OLS・ordinary least square)を図で確認します。この記事では、説明変数(係数パラメータ)が1次元(定数項を含めると2次元)の場合を扱います。
　線形回帰モデルの定義や計算式については「3.2.1-3：線形回帰モデルの最小二乗法の導出【空間データサイエンス入門のノート】 - からっぽのしょこ」、推定処理については「【Python】3.2.1-3：線形回帰モデルの最小二乗法の実装【空間データサイエンス入門のノート】 - からっぽのしょこ」、2次元の場合については「【Python】3.2.1-3：線形回帰モデルの最小二乗法の可視化：2次元の場合【空間データサイエンス入門のノート】 - からっぽのしょこ」を参照してください。

　利用するライブラリを読み込みます。

# ライブラリを読込
import numpy as np
from scipy.stats import norm # 確率密度の計算用
import matplotlib.pyplot as plt

最小二乗法の図

　まずは、線形回帰モデルに対するOLSの定義や性質を可視化します。
　この例では、2次元の図で可視化するため、説明変数の次元(係数の数)が1の場合のみ対応します。

真のモデルの設定

　データ数を指定して、説明変数を作成します。

# データ数を指定
N = 10

# 係数パラメータ数を指定:(固定)
K = 1

# 説明変数の値を生成
dat = np.random.uniform(low=-5.0, high=5.0, size=(N, K)) # 範囲を指定
print(dat[:5].round(2))

# 説明変数を作成
X = np.hstack([np.ones(shape=(N, 1)), dat])
print(X[:5].round(2))
print(X.shape)

[[4.4 ]
 [2.03]
 [3.23]
 [2.46]
 [2.6 ]]
[[1.   4.4 ]
 [1.   2.03]
 [1.   3.23]
 [1.   2.46]
 [1.   2.6 ]]
(10, 2)

　観測データ数(サンプルサイズ) $N$ 、係数パラメータ数 $K=1$ を指定します。
　データごとに1番目の次元の値 $x_{i1}$ を作成して、0番目の次元(定数項用)の値 $x_{i0} = 1$ と列方向に結合して、 $N$ 個の説明変数 $\mathbf{x}_i = (x_{i0}, x_{i1})$ 、 $\mathbf{X} = (\mathbf{x}_1, \cdots, \mathbf{x}_N)$ を作成します。
　(0番目を含めるので)説明変数や係数パラメータの要素(列)数は $K+1$ です。

　説明変数のヒストグラムを作成します。

作図コード(クリックで展開)

# プロット位置の調整用
offset_y = 0.1

# 説明変数を作図
fig, ax = plt.subplots(figsize=(8, 6), dpi=100, facecolor='white', constrained_layout=True)
ax.hist(X[:, 1], color='gray') # 説明変数の度数
ax.scatter(X[:, 1], np.zeros(N)-offset_y, 
           facecolor='none', edgecolor='C0', s=50, label='$x_{i1}$') # 説明変数
ax.set_xlabel('$x_{i1}$')
ax.set_ylabel('count')
ax.set_title(f'$N = {N}, K = {K}$', loc='left')
fig.suptitle('explanatory variables', fontsize=20)
ax.grid()
ax.legend(loc='upper right')
plt.show()

　横軸は各説明変数 $x_{i1}$ の値、縦軸は $N$ 個の説明変数 $\mathbf{X}$ の度数です。
　目安として、データ点を描画しています。

　係数パラメータを作成します。

# 真の係数パラメータを指定
beta_true = np.array([2.0, 0.5])
print(beta_true)

[2.  0.5]

　真の係数パラメータ $\boldsymbol{\beta} = (\beta_0, \beta_1)$ を指定します。

　真のモデルのグラフを作成します。

作図コード(クリックで展開)

# グラフサイズを設定
x_min = np.floor(X[:, 1].min())
x_max = np.ceil(X[:, 1].max())

# 真のモデルの座標を作成
line_x      = np.linspace(start=x_min, stop=x_max, num=1001)
line_y_true = beta_true[0] + beta_true[1] * line_x

# グラフサイズを設定
y_min = np.floor(line_y_true.min())
y_max = np.ceil(line_y_true.max())

# ラベル用の文字列を作成
param_label  = f'$N = {N}, K = {K}, '
param_label += f'\\beta = ({beta_true[0]:.2f}, {beta_true[1]:.2f})$'

# 真のモデルを作図
fig, ax = plt.subplots(figsize=(8, 6), dpi=100, facecolor='white', constrained_layout=True)
ax.quiver(x_min, 0.0, x_max-x_min, 0.0, 
          angles='xy', scale_units='xy', scale=1.0, 
          units='dots', width=2.5, headwidth=5.0, headlength=7.5, headaxislength=7.5) # x軸線
ax.quiver(0.0, y_min, 0.0, y_max-y_min, 
          angles='xy', scale_units='xy', scale=1.0, 
          units='dots', width=2.5, headwidth=5.0, headlength=7.5, headaxislength=7.5) # y軸線
ax.plot(line_x, line_y_true, 
        color='red', label='$y = \\beta_0 + \\beta_1 x_1$') # 真のモデル
ax.plot([0.0, 1.0], 
        [beta_true[0], beta_true[0]], 
        color='black') # 長さ1の線分
ax.plot([0.0, 0.0], [0.0, beta_true[0]], 
        color='black') # 定数の線分
ax.plot([1.0, 1.0], [beta_true[0], beta_true[0]+beta_true[1]], 
        color='black') # 係数の線分
ax.text(x=0.5, y=beta_true[0], s='$1$', 
        size=15, ha='center', va='top') # 1ラベル
ax.text(x=0.0, y=0.5*beta_true[0], s='$\\beta_0$', 
        size=15, ha='left', va='center') # 定数ラベル
ax.text(x=1.0, y=beta_true[0]+0.5*beta_true[1], s='$\\beta_1$', 
        size=15, ha='left', va='center') # 係数ラベル
ax.set_xticks(ticks=np.arange(x_min, x_max+1))
ax.set_yticks(ticks=np.arange(y_min, y_max+1))
ax.set_xlim(xmin=x_min, xmax=x_max)
ax.set_ylim(ymin=y_min, ymax=y_max)
ax.set_xlabel('$x_1$')
ax.set_ylabel('$y$')
ax.set_title(param_label, loc='left')
fig.suptitle('true model', fontsize=20)
ax.grid()
ax.legend()
plt.show()

　横軸は説明変数 $x_1$ の値、縦軸は被説明変数 $y$ の値です。
　 $K=1$ のとき、真のモデル $y = \mathbf{x}^{\top} \boldsymbol{\beta}$ は $y = \beta_0 + \beta_1 x_1$ であり、直線で表わせます。
　 $\beta_0$ は定数項(直線の切片)、 $\beta_1$ は説明変数の係数(直線の傾き・x軸方向に1変化したときのy軸方向の変化量)に対応します。

　分散パラメータを指定して、誤差項を作成します。

# 誤差項の真の標準偏差パラメータを指定
sigma_true = 0.5

# 誤差項の真の分散パラメータを計算
sigma2_true = sigma_true**2
print(sigma2_true)

# 誤差項を生成
epsilon = np.random.normal(loc=0.0, scale=sigma_true, size=N)
print(epsilon[:5].round(2))
print(epsilon.shape)

0.25
[-0.6   0.06  0.75  0.05 -0.12]
(10,)

　誤差項の標準偏差パラメータ $\sigma$ を指定して、分散パラメータ $\sigma^2$ を計算します。
　平均 $\mu = 0$ 、分散 $\sigma^2$ の正規分布に従って、 $N$ 個の誤差項 $\boldsymbol{\epsilon} = (\epsilon_1, \cdots, \epsilon_N)$ を作成します。

　誤差項のヒストグラムを作成します。

作図コード(クリックで展開)

# 正規分布の確率密度を計算
eps_size = max(epsilon.max(), abs(epsilon.min()), 2.0*sigma_true)
line_eps = np.linspace(start=-eps_size, stop=eps_size, num=1000)
dens_eps_true = norm.pdf(line_eps, loc=0.0, scale=sigma_true)

# プロット位置の調整用
offset_y = 0.1

# ラベル用の文字列を作成
param_label = f'$N = {N}, \\mu = 0, \\sigma = {sigma_true:.2f}, \\sigma^2 = {sigma2_true:.2f}$'

# 誤差項を作図
fig, ax = plt.subplots(figsize=(8, 6), dpi=100, facecolor='white', constrained_layout=True)
ax.plot(line_eps, dens_eps_true, 
        color='red', label='$N(0, \\sigma^2)$') # 真の確率密度
ax.hist(epsilon, density=True, color='gray') # 誤差項の密度
ax.scatter(epsilon, np.zeros(N)-offset_y, 
           facecolor='none', edgecolor='deeppink', s=50, label='$\\epsilon_i \\sim N(0, \\sigma^2)$') # 誤差項
ax.set_xlabel('$\\epsilon_i$')
ax.set_ylabel('density')
ax.set_title(param_label, loc='left')
fig.suptitle('error term', fontsize=20)
ax.grid()
ax.legend(loc='upper right')
plt.show()

　横軸は各誤差 $\epsilon_i$ の値、縦軸は $N$ 個の誤差 $\boldsymbol{\epsilon}$ の密度です。
　目安としてデータ点と確率密度曲線を描画しています。

　被説明変数(の実現値)を作成します。

# 被説明変数を計算
y = X @ beta_true + epsilon
print(y[:5].round(2))
print(y.shape)

[3.6  3.08 4.37 3.28 3.18]
(10,)

　被説明変数 $\mathbf{y} = (y_1, \cdots, y_N)$ を $\mathbf{y} = \mathbf{X} \boldsymbol{\beta} + \boldsymbol{\epsilon}$ で計算します。

　説明変数と被説明変数の関係のグラフを作成します。

作図コード(クリックで展開)

# グラフサイズを設定
y_min = np.floor(min(y.min(), line_y_true.min()))
y_max = np.ceil(max(y.max(), line_y_true.max()))

# ラベル表示用のデータ番号を設定
i = np.argmax(abs(epsilon))

# 真のモデル上の点の座標を計算
y_true = X @ beta_true

# ラベル用の文字列を作成
param_label  = f'$N = {N}, K = {K}, '
param_label += f'\\beta = ({beta_true[0]:.2f}, {beta_true[1]:.2f})$'
rv_label  = '$(x_{i1}, y_i), '
rv_label += 'y_i = \\beta_0 + \\beta_1 x_{i1} + \\epsilon_i$'

# 観測データを作図
fig, ax = plt.subplots(figsize=(8, 6), dpi=100, facecolor='white', constrained_layout=True)
ax.quiver(x_min, 0.0, x_max-x_min, 0.0, 
          angles='xy', scale_units='xy', scale=1.0, 
          units='dots', width=2.5, headwidth=5.0, headlength=7.5, headaxislength=7.5) # x軸線
ax.quiver(0.0, y_min, 0.0, y_max-y_min, 
          angles='xy', scale_units='xy', scale=1.0, 
          units='dots', width=2.5, headwidth=5.0, headlength=7.5, headaxislength=7.5) # y軸線
ax.plot(line_x, line_y_true, 
        color='red', label='$y = \\beta_0 + \\beta_1 x_1$') # 真のモデル
ax.scatter(X[:, 1], y, 
           color='C0', s=50, label=rv_label) # 実現値
ax.plot([X[:, 1], X[:, 1]], [y_true, y], 
        color='black', linestyle='dotted') # 誤差
ax.text(x=X[i, 1], y=y[i]-0.5*epsilon[i], s='$\\epsilon_i$', 
        size=15, ha='left', va='center') # 誤差ラベル
ax.set_xlim(xmin=x_min, xmax=x_max)
ax.set_ylim(ymin=y_min, ymax=y_max)
ax.set_xlabel('$x_{i1}$')
ax.set_ylabel('$y_i$')
ax.set_title(param_label, loc='left')
fig.suptitle('observation data', fontsize=20)
ax.grid()
ax.legend()
plt.show()

　横軸は説明変数 $x_{i1}$ の値、縦軸は被説明変数 $y_i$ の値です。
　「観測データの点 $(x_{i1}, y_i)$ 」と「真のモデル上の点 $(x_{i1}, \mathbf{x}_i^{\top} \boldsymbol{\beta})$ 」のy座標の差が誤差 $\epsilon_i = y_i - \mathbf{x}_i^{\top} \boldsymbol{\beta}$ に対応します。

係数パラメータの推定

　係数パラメータの推定値を計算します。

# 係数パラメータの推定値を計算
beta_hat = np.linalg.inv(X.T @ X) @ X.T @ y
print(beta_hat.round(2))
print(beta_hat.shape)

[2.12 0.45]
(2,)

　係数パラメータの推定値 $\hat{\boldsymbol{\beta}} = (\hat{\beta}_0, \hat{\beta}_1)$ を $\hat{\boldsymbol{\beta}} = (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^{\top} \mathbf{y}$ で計算します。

　回帰直線のグラフを作成します。

作図コード(クリックで展開)

# 回帰直線の座標を作成
line_y_hat = beta_hat[0] + beta_hat[1] * line_x

# グラフサイズを設定
y_min = np.floor(min(y.min(), line_y_true.min(), line_y_hat.min()))
y_max = np.ceil(max(y.max(), line_y_true.max(), line_y_hat.min()))

# ラベル用の文字列を作成
param_label  = f'$N = {N}, K = {K}, '
param_label += f'\\beta = ({beta_true[0]:.2f}, {beta_true[1]:.2f}), '
param_label += f'\\hat{{\\beta}} = ({beta_hat[0]:.2f}, {beta_hat[1]:.2f})$'

# 回帰直線を作図
fig, ax = plt.subplots(figsize=(8, 6), dpi=100, facecolor='white', constrained_layout=True)
ax.quiver(x_min, 0.0, x_max-x_min, 0.0, 
          angles='xy', scale_units='xy', scale=1.0, 
          units='dots', width=2.5, headwidth=5.0, headlength=7.5, headaxislength=7.5) # x軸線
ax.quiver(0.0, y_min, 0.0, y_max-y_min, 
          angles='xy', scale_units='xy', scale=1.0, 
          units='dots', width=2.5, headwidth=5.0, headlength=7.5, headaxislength=7.5) # y軸線
ax.plot(line_x, line_y_true, 
        color='red', linestyle='dashed', label='$y = \\beta_0 + \\beta_1 x_1$') # 真のモデル
ax.plot(line_x, line_y_hat, 
        color='orange', linewidth=2.0, label='$\\hat{y} = \\hat{\\beta}_0 + \\hat{\\beta}_1 x_1$') # 回帰直線
ax.scatter(X[:, 1], y, 
           color='C0', s=50, label=rv_label) # 実現値
ax.plot([0.0, 1.0], 
        [beta_hat[0], beta_hat[0]], 
        color='black') # 長さ1の線分
ax.plot([0.0, 0.0], [0.0, beta_hat[0]], 
        color='black') # 定数の線分
ax.plot([1.0, 1.0], [beta_hat[0], beta_hat[0]+beta_hat[1]], 
        color='black') # 係数の線分
ax.text(x=0.5, y=beta_hat[0], s='$1$', 
        size=15, ha='center', va='top') # 1ラベル
ax.text(x=0.0, y=0.5*beta_hat[0], s='$\\hat{\\beta}_0$', 
        size=15, ha='left', va='center') # 定数ラベル
ax.text(x=1.0, y=beta_hat[0]+0.5*beta_hat[1], s='$\\hat{\\beta}_1$', 
        size=15, ha='left', va='center') # 係数ラベル
ax.set_xticks(ticks=np.arange(x_min, x_max+1))
ax.set_yticks(ticks=np.arange(y_min, y_max+1))
ax.set_xlim(xmin=x_min, xmax=x_max)
ax.set_ylim(ymin=y_min, ymax=y_max)
ax.set_xlabel('$x_1$')
ax.set_ylabel('$\\hat{y}$')
ax.set_title(param_label, loc='left')
fig.suptitle('estimated coefficient parameters', fontsize=20)
ax.grid()
ax.legend()
plt.show()

　真の回帰直線(真のモデル)を赤色の破線、回帰直線(推定したモデル)をオレンジ色の実線で示します。

　真のモデルと同様に、推定したモデル $\hat{y} = \mathbf{x}^{\top} \hat{\boldsymbol{\beta}}$ は $\hat{y} = \beta_0 + \hat{\beta}_1 x_1$ であり、直線で表わせます。
　 $\hat{\beta}_0$ は定数項(直線の切片)、 $\hat{\beta}_1$ は説明変数の係数(直線の傾き)に対応します。

　被説明変数の推定値(理論値)を計算します。

# 被説明変数の推定値を計算
y_hat = X @ beta_hat
print(y_hat[:5].round(2))
print(y_hat.shape)

[4.12 3.05 3.59 3.24 3.3 ]
(10,)

　被説明変数の推定値 $\hat{\mathbf{y}} = (\hat{y}_1, \cdots, \hat{y}_N)$ を $\hat{\mathbf{y}} = \mathbf{X} \hat{\boldsymbol{\beta}}$ で計算します。

　説明変数と被説明変数の推定値の関係のグラフを作成します。

作図コード(クリックで展開)

# ラベル用の文字列を作成
tv_label  = '$(x_{i1}, \\hat{y}_i), '
tv_label += '\\hat{y}_i = \\hat{\\beta}_0 + \\hat{\\beta}_1 x_{i1}$'

# 理論値を作図
fig, ax = plt.subplots(figsize=(8, 6), dpi=100, facecolor='white', constrained_layout=True)
ax.quiver(x_min, 0.0, x_max-x_min, 0.0, 
          angles='xy', scale_units='xy', scale=1.0, 
          units='dots', width=2.5, headwidth=5.0, headlength=7.5, headaxislength=7.5) # x軸線
ax.quiver(0.0, y_min, 0.0, y_max-y_min, 
          angles='xy', scale_units='xy', scale=1.0, 
          units='dots', width=2.5, headwidth=5.0, headlength=7.5, headaxislength=7.5) # y軸線
ax.plot(line_x, line_y_true, 
        color='red', linestyle='dashed', label='$y = \\beta_0 + \\beta_1 x_1$') # 真のモデル
ax.plot(line_x, line_y_hat, 
        color='orange', linewidth=2.0, label='$\\hat{y} = \\hat{\\beta}_0 + \\hat{\\beta}_1 x_1$') # 回帰直線
ax.scatter(X[:, 1], y, 
           color='C0', s=50, label=rv_label) # 実現値
ax.scatter(X[:, 1], y_hat, 
           color='C1', s=50, label=tv_label) # 理論値
ax.set_xticks(ticks=np.arange(x_min, x_max+1))
ax.set_yticks(ticks=np.arange(y_min, y_max+1))
ax.set_xlim(xmin=x_min, xmax=x_max)
ax.set_ylim(ymin=y_min, ymax=y_max)
ax.set_xlabel('$x_{i1}$')
ax.set_ylabel('$y_i, \\hat{y}_i$')
ax.set_title(param_label, loc='left')
fig.suptitle('explained variables', fontsize=20)
ax.grid()
ax.legend()
plt.show()

　横軸は説明変数 $x_{i1}$ の値、縦軸は被説明変数 $y_i, \hat{y}_i$ の値です。
　 $\mathbf{y}$ の推定値 $\hat{\mathbf{y}}$ は、回帰直線 $\hat{y} = \mathbf{x}^{\top} \hat{\boldsymbol{\beta}}$ 上の点であり、理論値とも呼ばれます。

　残差を作成します。

# 残差を計算
e = y - y_hat
print(e[:5].round(2))
print(e.shape)

[-0.52  0.03  0.78  0.04 -0.12]
(10,)

　残差 $\mathbf{e} = (e_1, \cdots, e_N)$ を $\mathbf{e} = \mathbf{y} - \hat{\mathbf{y}}$ で計算します。

　実現値と理論値の関係のグラフを作成します。

作図コード(クリックで展開)

# ラベル表示用のデータ番号を設定
i = np.argmax(abs(e))

# 残差を作図
fig, ax = plt.subplots(figsize=(8, 6), dpi=100, facecolor='white', constrained_layout=True)
ax.quiver(x_min, 0.0, x_max-x_min, 0.0, 
          angles='xy', scale_units='xy', scale=1.0, 
          units='dots', width=2.5, headwidth=5.0, headlength=7.5, headaxislength=7.5) # x軸線
ax.quiver(0.0, y_min, 0.0, y_max-y_min, 
          angles='xy', scale_units='xy', scale=1.0, 
          units='dots', width=2.5, headwidth=5.0, headlength=7.5, headaxislength=7.5) # y軸線
ax.plot(line_x, line_y_true, 
        color='red', linestyle='dashed', label='$y = \\beta_0 + \\beta_1 x_1$') # 真のモデル
ax.plot(line_x, line_y_hat, 
        color='orange', linewidth=2.0, label='$\\hat{y} = \\hat{\\beta}_0 + \\hat{\\beta}_1 x$') # 回帰直線
ax.scatter(X[:, 1], y, 
           color='C0', s=50, label=rv_label) # 実現値
ax.scatter(X[:, 1], y_hat, 
           color='C1', s=50, label=tv_label) # 理論値
ax.plot([X[:, 1], X[:, 1]], [y_hat, y], 
        color='black', linestyle='dotted') # 残差
ax.text(x=X[i, 1], y=y[i]-0.5*e[i], s='$e_i$', 
        size=15, ha='left', va='center') # 残差ラベル
ax.set_xlim(xmin=x_min, xmax=x_max)
ax.set_ylim(ymin=y_min, ymax=y_max)
ax.set_xlabel('$x_{i1}$')
ax.set_ylabel('$y_i, \\hat{y}_i$')
ax.set_title(param_label, loc='left')
fig.suptitle('residual error', fontsize=20)
ax.grid()
ax.legend()
plt.show()

　「実現値の点 $(x_{i1}, y_i)$ 」と「理論値(回帰直線上)の点 $(x_{i1}, \hat{y}_i)$ 」のy座標の差が残差 $e_i = y_i - \hat{y}_i$ に対応します。

分散パラメータの推定

　分散パラメータの推定値を作成します。

# 残差平方和を計算
rss = e.T @ e
print(rss)

# 分散パラメータの推定値を計算
sigma2_hat = rss / (N - K - 1)
print(sigma2_hat)

# 標準偏差パラメータの推定値を計算
sigma_hat = np.sqrt(sigma2_hat)
print(sigma_hat)

1.6642842633765853
0.20803553292207316
0.45610912391890734

　分散パラメータの推定値 $\hat{\sigma}^2 = \frac{\mathbf{e}^{\top} \mathbf{e}}{N - (K+1)}$ と標準偏差パラメータの推定値 $\hat{\sigma} = \sqrt{\hat{\sigma}^2}$ を計算します。

　誤差項の真の分布と推定した分布のグラフを作成します。

作図コード(クリックで展開)

# 正規分布の確率密度を計算
dens_eps_hat = norm.pdf(line_eps, loc=0.0, scale=sigma_hat)

# ラベル用の文字列を作成
param_label  = f'$N = {N}, \\mu = 0, '
param_label += f'\\sigma = {sigma_true:.2f}, \\sigma^2 = {sigma2_true:.2f}, '
param_label += f'\\hat{{\\sigma}} = {sigma_hat:.2f}, \\hat{{\\sigma}}^2 = {sigma2_hat:.2f}$'

# 誤差項を作図
fig, ax = plt.subplots(figsize=(8, 6), dpi=100, facecolor='white', constrained_layout=True)
ax.plot(line_eps, dens_eps_true, 
        color='red', linestyle='dashed', label='$N(0, \\sigma^2)$') # 真の確率密度
ax.plot(line_eps, dens_eps_hat, 
        color='black', label='$N(0, \\hat{\\sigma}^2)$') # 確率密度の推定値
ax.scatter(epsilon, np.zeros(N), 
           facecolor='none', edgecolor='deeppink', s=50, label='$\\epsilon_i \\sim N(0, \\sigma^2)$') # 誤差項
ax.set_xlabel('$\\epsilon$')
ax.set_ylabel('density')
ax.set_title(param_label, loc='left')
fig.suptitle('error term', fontsize=20)
ax.grid()
ax.legend(loc='upper right')
plt.show()

　横軸は誤差 $\epsilon$ の値、縦軸は誤差の生成分布(正規分布)の確率密度です。
　平均 $\mu = 0$ 、分散 $\sigma^2$ の分布(真の分布)を赤色の破線、分散 $\hat{\sigma}^2$ の分布(推定分布)を黒色の実線で示します。

係数パラメータの分散共分散行列の推定

　係数パラメータの推定値の分散共分散行列を作成します。

# 係数パラメータの推定値の分散共分散行列を計算
Sigma_beta_hat = sigma2_hat * np.linalg.inv(X.T @ X)
print(Sigma_beta_hat.round(3))
print(Sigma_beta_hat.shape)

[[ 0.023 -0.002]
 [-0.002  0.002]]
(2, 2)

　係数パラメータの推定値 $\hat{\boldsymbol{\beta}}$ の分散共分散行列 $\boldsymbol{\Sigma}_{\hat{\boldsymbol{\beta}}} = \sigma^2 (\mathbf{X}^{\top} \mathbf{X})^{-1}$ を計算します。

　係数パラメータの推定値のマハラノビス距離を作成します。

# マハラノビス距離を計算
delta_val = np.sqrt(
    (beta_hat - beta_true).T @ np.linalg.inv(Sigma_beta_hat) @ (beta_hat - beta_true)
)
print(delta_val)

1.1471505715707588

　 $\hat{\boldsymbol{\beta}}$ の期待値 $\boldsymbol{\beta}$ と $\boldsymbol{\Sigma}_{\hat{\boldsymbol{\beta}}}$ を用いた $\hat{\boldsymbol{\beta}}$ のマハラノビス距離 $\Delta_{\hat{\boldsymbol{\beta}}} = \sqrt{(\hat{\boldsymbol{\beta}} - \boldsymbol{\beta})^{\top} \boldsymbol{\Sigma}_{\hat{\boldsymbol{\beta}}}^{-1} (\hat{\boldsymbol{\beta}} - \boldsymbol{\beta})}$ を計算します。

処理コード(クリックで展開)

　等高線の描画用のマハラノビス距離を作成します。

# 各軸の範囲の調整値を指定
sgm_num = 3.0

# 各軸の範囲を設定
beta0_min = beta_true[0] - sgm_num * np.sqrt(Sigma_beta_hat[0, 0])
beta0_max = beta_true[0] + sgm_num * np.sqrt(Sigma_beta_hat[0, 0])
beta1_min = beta_true[1] - sgm_num * np.sqrt(Sigma_beta_hat[1, 1])
beta1_max = beta_true[1] + sgm_num * np.sqrt(Sigma_beta_hat[1, 1])
print(beta0_min.round(2), beta0_max.round(2))
print(beta1_min.round(2), beta1_max.round(2))

# 各軸の値を作成
beta0_vec = np.linspace(start=beta0_min, stop=beta0_max, num=50)
beta1_vec = np.linspace(start=beta1_min, stop=beta1_max, num=50)
print(beta0_vec[:5].round(2))
print(beta1_vec[:5].round(2))

# 格子点を作成
beta0_grid, beta1_grid = np.meshgrid(beta0_vec, beta1_vec)

# 座標を作成
beta_arr = np.vstack([beta0_grid.flatten(), beta1_grid.flatten()]).T
print(beta_arr[:5].round(2))
print(beta_arr.shape)

# マハラノビス距離を計算
inv_Sigma_beta_hat = np.linalg.inv(Sigma_beta_hat) # 精度行列
delta_grid = np.array(
    [np.sqrt((beta - beta_true).T @ inv_Sigma_beta_hat @ (beta - beta_true)) for beta in beta_arr]
).reshape(beta0_grid.shape)
print(delta_grid[:5, :5].round(2))
print(delta_grid.shape)

1.55 2.45
0.37 0.63
[1.55 1.57 1.59 1.6  1.62]
[0.37 0.37 0.38 0.38 0.39]
[[1.55 0.37]
 [1.57 0.37]
 [1.59 0.37]
 [1.6  0.37]
 [1.62 0.37]]
(2500, 2)
[[4.99 4.89 4.79 4.69 4.59]
 [4.89 4.79 4.69 4.59 4.49]
 [4.79 4.69 4.58 4.48 4.38]
 [4.69 4.59 4.48 4.38 4.28]
 [4.59 4.49 4.38 4.28 4.18]]
(50, 50)

　 $\hat{\beta}_0, \hat{\beta}_1$ がとり得る値をそれぞれ作成し、格子点 $(\hat{\beta}_0, \hat{\beta}_1)$ に変換して、マハラノビス距離を計算します。この例では、真の値 $\beta_0, \beta_1$ を中心として標準偏差の sgm_num 倍を範囲とします。

　真のパラメータと推定パラメータの関係のグラフを作成します。

作図コード(クリックで展開)

# ラベル用の文字列を作成
sigma_str = '(' + ', '.join('(' + ', '.join(str(val.round(2)) for val in vec) + ')' for vec in Sigma_beta_hat) + ')'
param_label  = f'$N = {N}, K = {K}, ' + '\\Sigma_{\\hat{\\beta}} = ' + sigma_str
param_label += f', \\Delta_{{\\hat{{\\beta}}}} = {delta_val:.2f}$'

# 係数パラメータのマハラノビス距離を作図
fig, ax = plt.subplots(figsize=(8, 6), dpi=100, facecolor='white', constrained_layout=True)
cs = ax.contour(beta0_grid, beta1_grid, delta_grid) # マハラノビス距離
ax.scatter(*beta_true, 
           color='red', s=100, 
           label='$\\beta = ({:.2f}, {:.2f})$'.format(*beta_true)) # 真のパラメータ
ax.scatter(*beta_hat, 
           color='orange', s=100, 
           label='$\\hat{{\\beta}} = ({:.2f}, {:.2f})$'.format(*beta_hat)) # 推定パラメータ
ax.plot([beta_true[0], beta_hat[0]], [beta_true[1], beta_hat[1]], 
        color='black', linestyle='dotted') # 距離
ax.text(x=beta_true[0]+0.5*(beta_hat[0]-beta_true[0]), 
        y=beta_true[1]+0.5*(beta_hat[1]-beta_true[1]), 
        s='$\\Delta_{\\hat{\\beta}}$', size=15, ha='left', va='center') # 距離ラベル
ax.set_xlabel('$\\hat{\\beta}_0$')
ax.set_ylabel('$\\hat{\\beta}_1$')
ax.set_title(param_label, loc='left')
fig.suptitle("Mahalanobis' distance", fontsize=20)
fig.colorbar(cs, label='distance')
ax.grid()
ax.legend(loc='upper left')
#ax.set_aspect('equal', adjustable='box')
plt.show()

# 距離軸の範囲を設定
delta_min = 0.0
delta_max = np.ceil(delta_grid.max())

# 軸サイズを設定
beta0_size = beta0_max - beta0_min
beta1_size = beta1_max - beta1_min

# 係数パラメータのマハラノビス距離を作図
fig, ax = plt.subplots(figsize=(9, 8), dpi=100, facecolor='white', constrained_layout=True, 
                       subplot_kw={'projection': '3d'})
ax.contour(beta0_grid, beta1_grid, delta_grid, offset=delta_min) # マハラノビス距離
ax.scatter(*beta_true, delta_min, 
           color='red', s=100, 
           label='$\\beta = ({:.2f}, {:.2f})$'.format(*beta_true)) # 真のパラメータ
ax.scatter(*beta_hat, delta_min, 
           facecolor='none', edgecolor='orange', s=100) # 推定パラメータ(座標平面上)
ax.scatter(*beta_hat, delta_val, 
           color='orange', s=100, 
           label='$\\hat{{\\beta}} = ({:.2f}, {:.2f})$'.format(*beta_hat)) # 推定パラメータ(曲面上)
ax.plot([beta_true[0], beta_hat[0]], [beta_true[1], beta_hat[1]], [delta_min, delta_min], 
        color='black', linestyle='dotted', label='Euclid') # ユークリッド距離
ax.plot([beta_hat[0], beta_hat[0]], [beta_hat[1], beta_hat[1]], [delta_min, delta_val], 
        color='black', linestyle='dashed', label='Mahalanobis') # マハラノビス距離
ax.plot_surface(beta0_grid, beta1_grid, delta_grid, 
                cmap='viridis', vmin=delta_min, vmax=delta_max, alpha=0.8) # 曲面
ax.set_xlabel('$\\hat{\\beta}_0$')
ax.set_ylabel('$\\hat{\\beta}_1$')
ax.set_zlabel('$\\Delta$', labelpad=0.0)
ax.set_title(param_label, loc='left')
fig.suptitle("Mahalanobis' distance", fontsize=20)
ax.legend(loc='upper left')
#ax.set_box_aspect([1.0, beta1_size/beta0_size, delta_max/beta0_size])
plt.show()

　横軸は定数 $\hat{\beta}_0$ の値、縦軸は係数 $\hat{\beta}_1$ の値、高さ軸はマハラノビス距離 $\Delta_{\hat{\boldsymbol{\beta}}}$ です。
　真の係数パラメータの点 $(\beta_0, \beta_1)$ を赤色、推定した係数パラメータの点 $(\hat{\beta}_0, \hat{\beta}_1)$ をオレンジ色で示します。
　推定パラメータの期待値(真のパラメータ) $\mathbb{E}[\hat{\boldsymbol{\beta}}] = \boldsymbol{\beta}$ を中心とする $\boldsymbol{\Sigma}_{\hat{\boldsymbol{\beta}}}$ を用いたマハラノビス距離を等高線または曲面で示します。
　マハラノビス距離については「【Python】2次元マハラノビス距離の作図 - からっぽのしょこ」を参照してください。

　以上で、線形回帰モデルに対するOLSを可視化できました。続いて、パラメータ類の影響を見ていきます。

データ数の影響

　続いて、データ数と推定結果の関係を図で確認します。
　アニメーションの作図コードについては「Introduction-to-SpatialDataScience/code/ch3/1d_ols.py at main · anemptyarchive/Introduction-to-SpatialDataScience · GitHub」を参照してください。

係数パラメータの図

　データ数と回帰直線(係数パラメータ)の関係をアニメーションで確認します。

　真の回帰直線(真のモデル)を赤色の破線、回帰直線(推定したモデル)をオレンジ色の実線で示します。

　観測データが増えると回帰直線(推定パラメータ $\hat{\boldsymbol{\beta}}$ )が真のモデル(真のパラメータ $\boldsymbol{\beta}$ )に近付くのを確認できます。

分散パラメータの図

　データ数と誤差項の生成分布(分散パラメータ)の関係をアニメーションで確認します。

　真の生成分布(真のモデル)を赤色の破線、推定した生成分布(推定したモデル)を黒色の実線で示します。

　観測データが増えると推定したパラメータによる確率密度曲線(推定パラメータ $\hat{\sigma}^2$ の正規分布)が真のモデル(真のパラメータ $\sigma^2$ の正規分布)に近付くのを確認できます。
　ただし、分散パラメータ $\hat{\sigma}^2 = \frac{\mathbf{e}^{\top} \mathbf{e}}{N - (K+1)}$ に関して、0以上の値の条件 $\hat{\sigma}^2 \geq 0$ を満たす必要があるため、分子(残差平方和)は常に0以上の値 $\mathbf{e}^{\top} \mathbf{e} \geq 0$ なので分母(データ数と次元数の大小関係)が0以上の値 $N - (K+1) \gt 0$ である必要があります。つまり、データ数( $\mathbf{X}$ の行数)が次元数( $\mathbf{X}$ の列数)よりも大きい $N \gt (K + 1)$ 必要があります。(計算できないデータ数の範囲のグラフもアニメーションに含めています。)

係数パラメータの分散共分散行列の図

　データ数とマハラノビス距離(分散共分散行列)の関係をアニメーションで確認します。

　真の係数パラメータの点(真のモデル)を赤色、推定した係数パラメータの点(推定したモデル)をオレンジ色で示します。
　また、マハラノビス距離の等高線の配色は共通ですが間隔はフレームごとに異なります。

　観測データが増えると推定パラメータ $\hat{\boldsymbol{\beta}}$ が真のパラメータ $\boldsymbol{\beta}$ に近付く(マハラノビス距離が小さくなる)のを確認できます。
　また、分散共分散行列 $\boldsymbol{\Sigma}_{\hat{\boldsymbol{\beta}}}$ が小さく(精度行列が大きく)なり、同じ位置(座標)のマハラノビス距離が大きくなります。
　ただし、分散共分散行列 $\boldsymbol{\Sigma}_{\hat{\boldsymbol{\beta}}}$ の計算において分散パラメータ $\hat{\sigma}^2$ が含まれるので、こちらもデータ数( $\mathbf{X}$ の行数)が次元数( $\mathbf{X}$ の列数)よりも大きい $N \gt (K + 1)$ 必要があります。(計算できないデータ数の範囲のグラフもアニメーションに含めています。)

　この記事では、線形回帰モデルの最小二乗法を1次元の場合の図で確認しました。次の記事では、2次元の場合を確認します。

参考文献

Pythonで学ぶ空間データサイエンス入門: 地域の特徴を発見する方法

作者:植杉　大
Independently published

Amazon

おわりに

　回帰直線の図を作るだけだし実装編の記事の中で扱うつもりだったのですが、3つの推定パラメータの図はそれぞれで作った方がいいか、誤差と残差もややこしいから比較する図がいるな、いや一緒に描画すると見分けにくいからそれぞれ図解した方がいいな、あぁもう面倒だから計算する度に全部作るか、という感じで盛り沢山な記事になりました。
　と言いつつ共分散行列をどう表現したもの悩んで記事から落とす直前でした。が、そういやなんか書いたことあったなとマハラノビス距離を閃いて(思い出して)、自分の記事をググってなんとか形にできました。いやまぁ普通にt分布の図にすればよかったんですが。

　次節からのモデル評価も何か書くつもりだったのですが、ここまでの内容を読めたのなら本を読めば十分だと思いますし、深入りするにはそれぞれが重そうですし、一旦飛ばして4章以降で必要になったら書こうと思います。

　と上まで書き予約投稿を設定して4章を読んでたのですが、2次元版の図も作ろうと思っていたのを思い出したので書き始めました。まだ3章の記事が続きます。

【次の内容】

からっぽのしょこ

読んだら書く！書いたら読む！同じ事は二度調べ(たく)ない

【Python】3.2.1-3：線形回帰モデルの最小二乗法の可視化：1次元の場合【空間データサイエンス入門のノート】

はじめに

3.2.1-3 線形回帰モデルの最小二乗法の可視化：1次元の場合：

最小二乗法の図

真のモデルの設定

係数パラメータの推定

分散パラメータの推定

係数パラメータの分散共分散行列の推定

データ数の影響

係数パラメータの図

分散パラメータの図

係数パラメータの分散共分散行列の図

参考文献

おわりに