# ライブラリを読込
import numpy as np
from scipy.stats import norm # 確率密度の計算用
import matplotlib.pyplot as plt
import matplotlib.cm as cm # 配色の統一用
from matplotlib.animation import FuncAnimation # アニメーションの作図用

最小二乗法の図

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

真のモデルの設定

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

# データ数を指定
N = 10

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

# 説明変数の値を生成
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)

[[ 1.73 -1.74]
 [-1.02 -1.12]
 [ 4.43 -3.93]
 [-1.   -4.16]
 [ 1.79  3.26]]
[[ 1.    1.73 -1.74]
 [ 1.   -1.02 -1.12]
 [ 1.    4.43 -3.93]
 [ 1.   -1.   -4.16]
 [ 1.    1.79  3.26]]
(10, 3)

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

　説明変数の散布図を作成します。

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

# グラフサイズを設定
x_size = np.ceil(max(X.max(), abs(X.min())))

# 説明変数を作図
fig, ax = plt.subplots(figsize=(6, 6), dpi=100, facecolor='white', constrained_layout=True)
ax.scatter(X[:, 1], X[:, 2], 
           facecolor='none', edgecolor='C0', s=50, label='$(x_{i1}, x_{i2})$') # 説明変数
ax.set_xlim(xmin=-x_size, xmax=x_size)
ax.set_ylim(ymin=-x_size, ymax=x_size)
ax.set_xlabel('$x_{i1}$')
ax.set_ylabel('$x_{i2}$')
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()

　横軸は各説明変数の $1$ 番目の次元 $x_{i1}$ の値、縦軸は各説明変数の $2$ 番目の次元 $x_{i2}$ の値として、 $N$ 個の説明変数 $\mathbf{X}$ をデータ点 $(x_{i1}, x_{i2})$ の散布図で示します。

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

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

[ 2.   1.  -0.5]

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

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

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

# グラフサイズを設定
x1_min = np.floor(X[:, 1].min())
x1_max = np.ceil(X[:, 1].max())
x2_min = np.floor(X[:, 2].min())
x2_max = np.ceil(X[:, 2].max())

# 平面用の値を作成
x1_vec = np.linspace(start=x1_min, stop=x1_max, num=11)
x2_vec = np.linspace(start=x2_min, stop=x2_max, num=11)

# 格子点を作成
x1_grid, x2_grid = np.meshgrid(x1_vec, x2_vec)

# 真のモデルの平面の座標を作成
y_grid = beta_true[0] + beta_true[1] * x1_grid + beta_true[2] * x2_grid

# 真のモデルの直線の座標を作成
y_true_x1_vec = beta_true[0] + beta_true[1] * x1_vec
y_true_x2_vec = beta_true[0] + beta_true[2] * x2_vec

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

# 矢サイズを指定
alr = 0.5

# ラベル用の文字列を作成
param_label  = f'$N = {N}, K = {K}, '
param_label += f'\\beta = (' + ', '.join(str(val.round(2)) for val in beta_true) +')$'

# 真のモデルを作図
fig, ax = plt.subplots(figsize=(8, 8), dpi=100, facecolor='white', constrained_layout=True, 
                       subplot_kw={'projection': '3d'})

# x1・x2・y空間
ax.quiver(x1_min, 0, 0, x1_max-x1_min, 0, 0, 
          arrow_length_ratio=alr/(x1_max-x1_min), color='black', linewidth=1) # x1軸線
ax.quiver(0, x2_min, 0, 0, x2_max-x2_min, 0, 
          arrow_length_ratio=alr/(x2_max-x2_min), color='black', linewidth=1) # x2軸線
ax.quiver(0, 0, y_min, 0, 0, y_max-y_min, 
          arrow_length_ratio=alr/(y_max-y_min), color='black', linewidth=1) # y軸線
ax.plot_surface(x1_grid, x2_grid, y_grid, 
                color='red', alpha=0.5, 
                label='$y = \\beta_0 + \\beta_1 x_1 + \\beta_2 x_2$') # 真のモデル
ax.plot(x1_vec, np.zeros(len(x1_vec)), y_true_x1_vec, 
        color='red', linewidth=1) # 真のモデル:x1軸方向
ax.plot(np.zeros(len(x2_vec)), x2_vec, y_true_x2_vec, 
        color='red', linewidth=1) # 真のモデル:x2軸方向

# x2・y平面
ax.quiver(x1_min, x2_max, 0, x1_max-x1_min, 0, 0, 
          arrow_length_ratio=alr/(x1_max-x1_min), 
          color='black', linewidth=1, linestyle='dashed') # x1軸線
ax.quiver(0, x2_max, y_min, 0, 0, y_max-y_min, 
          arrow_length_ratio=alr/(y_max-y_min), 
          color='black', linewidth=1, linestyle='dashed') # y軸線
ax.plot(x1_vec, x2_max.repeat(len(x1_vec)), y_true_x1_vec, 
        color='red', linewidth=1, linestyle='dashed', 
        label='$y = \\beta_0 + \\beta_k x_k$') # 真のモデル:x1軸方向
ax.plot([0, 1], [x2_max, x2_max], [beta_true[0], beta_true[0]], 
        color='black', linewidth=1) # 長さ1の線分
ax.plot([0, 0], [x2_max, x2_max], [0, beta_true[0]], 
        color='black', linewidth=1) # 定数の線分
ax.plot([1, 1], [x2_max, x2_max], [beta_true[0], beta_true[0]+beta_true[1]], 
        color='black', linewidth=1) # 係数の線分
ax.text(0.5, x2_max, beta_true[0], s='$1$', 
        ha='center', va='top') # 1ラベル
ax.text(0, x2_max, 0.5*beta_true[0], s='$\\beta_0$', 
        ha='right', va='center') # 定数ラベル
ax.text(1, x2_max, beta_true[0]+0.5*beta_true[1], s='$\\beta_1$', 
        ha='left', va='center') # 係数ラベル

# x1・y平面
ax.quiver(x1_min, x2_min, 0, 0, x2_max-x2_min, 0, 
          arrow_length_ratio=alr/(x2_max-x2_min), 
          color='black', linewidth=1, linestyle='dashed') # x2軸線
ax.quiver(x1_min, 0, y_min, 0, 0, y_max-y_min, 
          arrow_length_ratio=alr/(y_max-y_min), 
          color='black', linewidth=1, linestyle='dashed') # y軸線
ax.plot(x1_min.repeat(len(x2_vec)), x2_vec, y_true_x2_vec, 
        color='red', linewidth=1, linestyle='dashed') # 真のモデル:x2軸方向
ax.plot([x1_min, x1_min], [0, 1], [beta_true[0], beta_true[0]], 
        color='black', linewidth=1) # 長さ1の線分
ax.plot([x1_min, x1_min], [0, 0], [0, beta_true[0]], 
        color='black', linewidth=1) # 定数の線分
ax.plot([x1_min, x1_min], [1, 1], [beta_true[0], beta_true[0]+beta_true[2]], 
        color='black', linewidth=1) # 係数の線分
ax.text(x1_min, 0.5, beta_true[0], s='$1$', 
        ha='center', va='bottom') # 1ラベル
ax.text(x1_min, 0, 0.5*beta_true[0], s='$\\beta_0$', 
        ha='right', va='center') # 定数ラベル
ax.text(x1_min, 1, beta_true[0]+0.5*beta_true[2], s='$\\beta_2$', 
        ha='left', va='center') # 係数ラベル

ax.set_xticks(ticks=np.arange(x1_min, x1_max+1))
ax.set_yticks(ticks=np.arange(x2_min, x2_max+1))
ax.set_zticks(ticks=np.arange(y_min, y_max+1))
ax.set_xlim(xmin=x1_min, xmax=x1_max)
ax.set_ylim(ymin=x2_min, ymax=x2_max)
ax.set_zlim(zmin=y_min, zmax=y_max)
ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$')
ax.set_zlabel('$y$', labelpad=0.0)
ax.set_title(param_label, loc='left')
fig.suptitle('true model', fontsize=20)
ax.legend()
#ax.view_init(elev=90, azim=270) # 0・1軸平面
#ax.view_init(elev=0, azim=270) # 0・2軸平面
#ax.view_init(elev=0, azim=0) # 1・2軸平面
plt.show()

　横軸は説明変数の $1$ 番目の次元 $x_1$ の値、縦軸は説明変数の $2$ 番目の次元 $x_2$ の値、高さ軸は被説明変数 $y$ の値です。
　 $K=2$ のとき、真のモデル $y = \mathbf{x}^{\top} \boldsymbol{\beta}$ は $y = \beta_0 + \beta_1 x_1 + \beta_2 x_2$ であり、平面で表わせます。
　 $\beta_0$ は定数項(切片)、 $\beta_1$ は説明変数の $1$ 番目の次元の係数( $x_1$ 軸方向の傾き)、 $\beta_2$ は説明変数の $2$ 番目の次元の係数( $x_2$ 軸方向の傾き)に対応します。
　目安としてグラフ領域の側面に、1つの軸の値が $x_k = 0$ のときのもう1つの軸の変化を直線として描画しています。平面上の実線と対応します。

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

# 誤差項の真の標準偏差パラメータを指定
sigma_true = 2.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)

6.25
[-0.01  1.28 -3.29  0.47  2.76]
(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.05

# ラベル用の文字列を作成
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)

[4.59 2.82 5.1  3.55 4.92]
(10,)

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

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

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

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

# 真のモデルの平面上の点の座標を作成
y_true = X @ beta_true

# 矢サイズを指定
alr = 0.5

# ラベル用の文字列を作成
param_label  = f'$N = {N}, K = {K}, '
param_label += f'\\beta = (' + ', '.join(str(val.round(2)) for val in beta_true) +')$'
rv_label  = '$(x_{i1}, y_i), '
rv_label += 'y_i = \\beta_0 + \\beta_1 x_{i1} + \\epsilon_i$'

# 真のモデルを作図
fig, ax = plt.subplots(figsize=(8, 8), dpi=100, facecolor='white', constrained_layout=True, 
                       subplot_kw={'projection': '3d'})

# x1・x2・y空間
ax.plot_surface(x1_grid, x2_grid, y_grid, 
                color='red', alpha=0.5, 
                label='$y = \\beta_0 + \\beta_1 x_1 + \\beta_2 x_2$') # 真のモデル
ax.scatter(X[:, 1], X[:, 2], y, color='C0', s=50, label=rv_label) # 実現値
for i in range(N):
    ax.plot([X[i, 1], X[i, 1]], [X[i, 2], X[i, 2]], [y_true[i], y[i]], 
            color='black', linestyle='dotted') # 誤差
i = np.argmax(abs(epsilon)) # ラベル表示用のデータ番号を設定
ax.text(x=X[i, 1], y=X[i, 2], z=y[i]-0.5*epsilon[i], s='$\\epsilon_i$', 
        size=15, ha='left', va='center') # 誤差ラベル

# x2・y平面
ax.quiver(x1_min, x2_max, 0, x1_max-x1_min, 0, 0, 
          arrow_length_ratio=alr/(x1_max-x1_min), 
          color='black', linewidth=1) # x1軸線
ax.quiver(0, x2_max, y_min, 0, 0, y_max-y_min, 
          arrow_length_ratio=alr/(y_max-y_min), 
          color='black', linewidth=1) # y軸線
ax.plot(x1_vec, x2_max.repeat(len(x1_vec)), y_true_x1_vec, 
        color='red', linewidth=1, linestyle='dashed', 
        label='$y = \\beta_0 + \\beta_k x_k$') # 真のモデル:x1軸方向

# x1・y平面
ax.quiver(x1_min, x2_min, 0, 0, x2_max-x2_min, 0, 
          arrow_length_ratio=alr/(x2_max-x2_min), 
          color='black', linewidth=1) # x2軸線
ax.quiver(x1_min, 0, y_min, 0, 0, y_max-y_min, 
          arrow_length_ratio=alr/(y_max-y_min), 
          color='black', linewidth=1) # y軸線
ax.plot(x1_min.repeat(len(x2_vec)), x2_vec, y_true_x2_vec, 
        color='red', linewidth=1, linestyle='dashed') # 真のモデル:x2軸方向

ax.set_xlim(xmin=x1_min, xmax=x1_max)
ax.set_ylim(ymin=x2_min, ymax=x2_max)
ax.set_zlim(zmin=y_min, zmax=y_max)
ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$')
ax.set_zlabel('$y$')
ax.set_title(param_label, loc='left')
fig.suptitle('observation data', fontsize=20)
ax.legend()
#ax.view_init(elev=90, azim=270) # 0・1軸平面
#ax.view_init(elev=0, azim=270) # 0・2軸平面
#ax.view_init(elev=0, azim=0) # 1・2軸平面
plt.show()

　横軸は説明変数の $1$ 番目の次元 $x_{i1}$ の値、縦軸は説明変数の $2$ 番目の次元 $x_{i2}$ の値、高さ軸は被説明変数 $y_i$ の値です。
　「観測データの点 $(x_{i1}, x_{i2}, y_i)$ 」と「真のモデル上の点 $(x_{i1}, x_{i2}, \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)

[ 1.78  0.75 -0.3 ]
(3,)

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

　回帰平面のグラフを作成します。

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

# 回帰平面の座標を作成
y_hat_grid = beta_hat[0] + beta_hat[1] * x1_grid + beta_hat[2] * x2_grid

# 回帰直線の座標を作成
y_hat_x1_vec = beta_hat[0] + beta_hat[1] * x1_vec
y_hat_x2_vec = beta_hat[0] + beta_hat[2] * x2_vec

# 矢サイズを指定
alr = 0.5

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

# 回帰平面を作図
fig, ax = plt.subplots(figsize=(8, 8), dpi=100, facecolor='white', constrained_layout=True, 
                       subplot_kw={'projection': '3d'})

# x1・x2・y空間
ax.quiver(x1_min, 0, 0, x1_max-x1_min, 0, 0, 
          arrow_length_ratio=alr/(x1_max-x1_min), color='black', linewidth=1) # x1軸線
ax.quiver(0, x2_min, 0, 0, x2_max-x2_min, 0, 
          arrow_length_ratio=alr/(x2_max-x2_min), color='black', linewidth=1) # x2軸線
ax.quiver(0, 0, y_min, 0, 0, y_max-y_min, 
          arrow_length_ratio=alr/(y_max-y_min), color='black', linewidth=1) # y軸線
ax.plot_wireframe(x1_grid, x2_grid, y_grid, 
                  color='red', alpha=0.5, linewidth=1, 
                  label='$y = \\beta_0 + \\beta_1 x_1 + \\beta_2 x_2$') # 真のモデル
ax.plot_surface(x1_grid, x2_grid, y_hat_grid, 
                color='orange', alpha=0.5, 
                label='$\\hat{y} = \\hat{\\beta}_0 + \\hat{\\beta}_1 x_1 + \\hat{\\beta}_2 x_2$') # 回帰平面
ax.plot(x1_vec, np.zeros(len(x1_vec)), y_hat_x1_vec, 
        color='orange', linewidth=1) # 回帰直線:x1軸方向
ax.plot(np.zeros(len(x2_vec)), x2_vec, y_hat_x2_vec, 
        color='orange', linewidth=1) # 回帰直線:x2軸方向
ax.scatter(X[:, 1], X[:, 2], y, color='C0', s=50, label=rv_label) # 実現値

# x2・y平面
ax.quiver(x1_min, x2_max, 0, x1_max-x1_min, 0, 0, 
          arrow_length_ratio=alr/(x1_max-x1_min), 
          color='black', linewidth=1, linestyle='dashed') # x1軸線
ax.quiver(0, x2_max, y_min, 0, 0, y_max-y_min, 
          arrow_length_ratio=alr/(y_max-y_min), 
          color='black', linewidth=1, linestyle='dashed') # y軸線
ax.plot(x1_vec, x2_max.repeat(len(x1_vec)), y_true_x1_vec, 
        color='red', linewidth=1, linestyle='dashed', 
        label='$y = \\beta_0 + \\beta_k x_k$') # 真のモデル:x1軸方向
ax.plot(x1_vec, x2_max.repeat(len(x1_vec)), y_hat_x1_vec, 
        color='orange', linewidth=1, linestyle='dashed', 
        label='$\\hat{y} = \\hat{\\beta}_0 + \\hat{\\beta}_k x_k$') # 回帰直線
ax.plot([0, 1], [x2_max, x2_max], [beta_hat[0], beta_hat[0]], 
        color='black', linewidth=1) # 長さ1の線分
ax.plot([0, 0], [x2_max, x2_max], [0, beta_hat[0]], 
        color='black', linewidth=1) # 定数の線分
ax.plot([1, 1], [x2_max, x2_max], [beta_hat[0], beta_hat[0]+beta_hat[1]], 
        color='black', linewidth=1) # 係数の線分
ax.text(0.5, x2_max, beta_hat[0], s='$1$', 
        ha='center', va='top') # 1ラベル
ax.text(0, x2_max, 0.5*beta_hat[0], s='$\\hat{\\beta}_0$', 
        ha='right', va='center') # 定数ラベル
ax.text(1, x2_max, beta_hat[0]+0.5*beta_hat[1], s='$\\hat{\\beta}_1$', 
        ha='left', va='center') # 係数ラベル

# x1・y平面
ax.quiver(x1_min, x2_min, 0, 0, x2_max-x2_min, 0, 
          arrow_length_ratio=alr/(x2_max-x2_min), 
          color='black', linewidth=1, linestyle='dashed') # x2軸線
ax.quiver(x1_min, 0, y_min, 0, 0, y_max-y_min, 
          arrow_length_ratio=alr/(y_max-y_min), 
          color='black', linewidth=1, linestyle='dashed') # y軸線
ax.plot(x1_min.repeat(len(x2_vec)), x2_vec, y_true_x2_vec, 
        color='red', linewidth=1, linestyle='dashed') # 真のモデル:x2軸方向
ax.plot(x1_min.repeat(len(x2_vec)), x2_vec, y_hat_x2_vec, 
        color='orange', linewidth=1, linestyle='dashed') # 回帰直線
ax.plot([x1_min, x1_min], [0, 1], [beta_hat[0], beta_hat[0]], 
        color='black', linewidth=1) # 長さ1の線分
ax.plot([x1_min, x1_min], [0, 0], [0, beta_hat[0]], 
        color='black', linewidth=1) # 定数の線分
ax.plot([x1_min, x1_min], [1, 1], [beta_hat[0], beta_hat[0]+beta_hat[2]], 
        color='black', linewidth=1) # 係数の線分
ax.text(x1_min, 0.5, beta_hat[0], s='$1$', 
        ha='center', va='bottom') # 1ラベル
ax.text(x1_min, 0, 0.5*beta_hat[0], s='$\\hat{\\beta}_0$', 
        ha='right', va='center') # 定数ラベル
ax.text(x1_min, 1, beta_hat[0]+0.5*beta_hat[2], s='$\\hat{\\beta}_2$', 
        ha='left', va='center') # 係数ラベル

ax.set_xticks(ticks=np.arange(x1_min, x1_max+1))
ax.set_yticks(ticks=np.arange(x2_min, x2_max+1))
ax.set_zticks(ticks=np.arange(y_min, y_max+1))
ax.set_xlim(xmin=x1_min, xmax=x1_max)
ax.set_ylim(ymin=x2_min, ymax=x2_max)
ax.set_zlim(zmin=y_min, zmax=y_max)
ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$')
ax.set_zlabel('$y$')
ax.set_title(param_label, loc='left')
fig.suptitle('estimated coefficient parameters', fontsize=20)
ax.legend()
#ax.view_init(elev=90, azim=270) # 0・1軸平面
#ax.view_init(elev=0, azim=270) # 0・2軸平面
#ax.view_init(elev=0, azim=0) # 1・2軸平面
plt.show()

　真の回帰平面(真のモデル)を赤色の格子面、回帰平面(推定したモデル)をオレンジ色の塗りつぶし面で示します。

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

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

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

[3.61 1.36 6.3  2.3  2.13]
(10,)

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

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

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

# 矢サイズを指定
alr = 0.5

# ラベル用の文字列を作成
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, 8), dpi=100, facecolor='white', constrained_layout=True, 
                       subplot_kw={'projection': '3d'})

# x1・x2・y空間
ax.plot_wireframe(x1_grid, x2_grid, y_grid, 
                  color='red', alpha=0.5, linewidth=1, 
                  label='$y = \\beta_0 + \\beta_1 x_1 + \\beta_2 x_2$') # 真のモデル
ax.plot_wireframe(x1_grid, x2_grid, y_hat_grid, 
                  color='orange', alpha=0.5, linewidth=1, 
                  label='$\\hat{y} = \\hat{\\beta}_0 + \\hat{\\beta}_1 x_1 + \\hat{\\beta}_2 x_2$') # 回帰平面
ax.scatter(X[:, 1], X[:, 2], y, color='C0', s=50, label=rv_label) # 実現値
ax.scatter(X[:, 1], X[:, 2], y_hat, color='C1', s=50, label=tv_label) # 理論値

# x2・y平面
ax.quiver(x1_min, x2_max, 0, x1_max-x1_min, 0, 0, 
          arrow_length_ratio=alr/(x1_max-x1_min), 
          color='black', linewidth=1) # x1軸線
ax.quiver(0, x2_max, y_min, 0, 0, y_max-y_min, 
          arrow_length_ratio=alr/(y_max-y_min), 
          color='black', linewidth=1) # y軸線
ax.plot(x1_vec, x2_max.repeat(len(x1_vec)), y_true_x1_vec, 
        color='red', linewidth=1, linestyle='dashed', 
        label='$y = \\beta_0 + \\beta_k x_k$') # 真のモデル:x1軸方向
ax.plot(x1_vec, x2_max.repeat(len(x1_vec)), y_hat_x1_vec, 
        color='orange', linewidth=1, linestyle='dashed', 
        label='$\\hat{y} = \\hat{\\beta}_0 + \\hat{\\beta}_k x_k$') # 回帰直線

# x1・y平面
ax.quiver(x1_min, x2_min, 0, 0, x2_max-x2_min, 0, 
          arrow_length_ratio=alr/(x2_max-x2_min), 
          color='black', linewidth=1) # x2軸線
ax.quiver(x1_min, 0, y_min, 0, 0, y_max-y_min, 
          arrow_length_ratio=alr/(y_max-y_min), 
          color='black', linewidth=1) # y軸線
ax.plot(x1_min.repeat(len(x2_vec)), x2_vec, y_true_x2_vec, 
        color='red', linewidth=1, linestyle='dashed') # 真のモデル:x2軸方向
ax.plot(x1_min.repeat(len(x2_vec)), x2_vec, y_hat_x2_vec, 
        color='orange', linewidth=1, linestyle='dashed') # 回帰直線

ax.set_xlim(xmin=x1_min, xmax=x1_max)
ax.set_ylim(ymin=x2_min, ymax=x2_max)
ax.set_zlim(zmin=y_min, zmax=y_max)
ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$')
ax.set_zlabel('$y$')
ax.set_title(param_label, loc='left')
fig.suptitle('explained variables', fontsize=20)
ax.legend()
#ax.view_init(elev=90, azim=270) # 0・1軸平面
#ax.view_init(elev=0, azim=270) # 0・2軸平面
#ax.view_init(elev=0, azim=0) # 1・2軸平面
plt.show()

　横軸は説明変数の $1$ 番目の次元 $x_{i1}$ の値、縦軸は説明変数の $2$ 番目の次元 $x_{i2}$ の値、高さ軸は被説明変数 $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.98  1.47 -1.19  1.26  2.79]
(10,)

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

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

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

# 矢サイズを指定
alr = 0.5

# 残差を作図
fig, ax = plt.subplots(figsize=(8, 8), dpi=100, facecolor='white', constrained_layout=True, 
                       subplot_kw={'projection': '3d'})

# x1・x2・y空間
ax.plot_wireframe(x1_grid, x2_grid, y_grid, 
                  color='red', alpha=0.5, linewidth=1, 
                  label='$y = \\beta_0 + \\beta_1 x_1 + \\beta_2 x_2$') # 真のモデル
ax.plot_surface(x1_grid, x2_grid, y_hat_grid, 
                color='orange', alpha=0.5, 
                label='$\\hat{y} = \\hat{\\beta}_0 + \\hat{\\beta}_1 x_1 + \\hat{\\beta}_2 x_2$') # 回帰平面
ax.scatter(X[:, 1], X[:, 2], y, color='C0', s=50, label=rv_label) # 実現値
for i in range(N):
    ax.plot([X[i, 1], X[i, 1]], [X[i, 2], X[i, 2]], [y_hat[i], y[i]], 
            color='black', linestyle='dotted') # 残差
i = np.argmax(abs(e)) # ラベル表示用のデータ番号を設定
ax.text(x=X[i, 1], y=X[i, 2], z=y[i]-0.5*e[i], s='$e_i$', 
        size=15, ha='left', va='center') # 残差ラベル

# x2・y平面
ax.quiver(x1_min, x2_max, 0, x1_max-x1_min, 0, 0, 
          arrow_length_ratio=alr/(x1_max-x1_min), 
          color='black', linewidth=1) # x1軸線
ax.quiver(0, x2_max, y_min, 0, 0, y_max-y_min, 
          arrow_length_ratio=alr/(y_max-y_min), 
          color='black', linewidth=1) # y軸線
ax.plot(x1_vec, x2_max.repeat(len(x1_vec)), y_true_x1_vec, 
        color='red', linewidth=1, linestyle='dashed', 
        label='$y = \\beta_0 + \\beta_k x_k$') # 真のモデル:x1軸方向
ax.plot(x1_vec, x2_max.repeat(len(x1_vec)), y_hat_x1_vec, 
        color='orange', linewidth=1, linestyle='dashed', 
        label='$\\hat{y} = \\hat{\\beta}_0 + \\hat{\\beta}_k x_k$') # 回帰直線

# x1・y平面
ax.quiver(x1_min, x2_min, 0, 0, x2_max-x2_min, 0, 
          arrow_length_ratio=alr/(x2_max-x2_min), 
          color='black', linewidth=1) # x2軸線
ax.quiver(x1_min, 0, y_min, 0, 0, y_max-y_min, 
          arrow_length_ratio=alr/(y_max-y_min), 
          color='black', linewidth=1) # y軸線
ax.plot(x1_min.repeat(len(x2_vec)), x2_vec, y_true_x2_vec, 
        color='red', linewidth=1, linestyle='dashed') # 真のモデル:x2軸方向
ax.plot(x1_min.repeat(len(x2_vec)), x2_vec, y_hat_x2_vec, 
        color='orange', linewidth=1, linestyle='dashed') # 回帰直線

ax.set_xlim(xmin=x1_min, xmax=x1_max)
ax.set_ylim(ymin=x2_min, ymax=x2_max)
ax.set_zlim(zmin=y_min, zmax=y_max)
ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$')
ax.set_zlabel('$y$')
ax.set_title(param_label, loc='left')
fig.suptitle('residual error', fontsize=20)
ax.legend()
#ax.view_init(elev=90, azim=270) # 0・1軸平面
#ax.view_init(elev=0, azim=270) # 0・2軸平面
#ax.view_init(elev=0, azim=0) # 1・2軸平面
plt.show()

　「実現値の点 $(x_{i1}, x_{i2}, y_i)$ 」と「理論値(回帰平面上)の点 $(x_{i1}, x_{i2}, \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)

27.413373787171995
3.916196255310285
1.9789381635893237

　分散パラメータの推定値 $\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)

# ラベル用の文字列を作成
dist_label  = f'$N = {N}, \\mu = 0, '
dist_label += f'\\sigma = {sigma_true:.2f}, \\sigma^2 = {sigma2_true:.2f}, '
dist_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(dist_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.44  0.021 0.047]
 [0.021 0.05  0.008]
 [0.047 0.008 0.048]]
(3, 3)

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

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

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

1.6502594844859126

　 $\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])
beta2_min = beta_true[2] - sgm_num * np.sqrt(Sigma_beta_hat[2, 2])
beta2_max = beta_true[2] + sgm_num * np.sqrt(Sigma_beta_hat[2, 2])
print(beta0_min.round(2), beta0_max.round(2))
print(beta1_min.round(2), beta1_max.round(2))
print(beta2_min.round(2), beta2_max.round(2))

# 2軸方向の線の数を指定
line_num = 100

# 各軸の値を作成
beta0_vec = np.linspace(start=beta0_min, stop=beta0_max, num=50) # 0軸方向の点の数を指定
beta1_vec = np.linspace(start=beta1_min, stop=beta1_max, num=50) # 1軸方向の点の数を指定
beta2_vec = np.linspace(start=beta2_min, stop=beta2_max, num=line_num)
print(beta0_vec[:5].round(2))
print(beta1_vec[:5].round(2))
print(beta2_vec[:5].round(2))

# 0・1軸の格子点を作成
beta0_grid, beta1_grid = np.meshgrid(beta0_vec, beta1_vec)

# 0・1軸の形状を設定
grid_shape = beta0_grid.shape
grid_size  = beta0_grid.size
print(grid_shape)
print(grid_size)

# 高さごとに処理
inv_Sigma_beta_hat = np.linalg.inv(Sigma_beta_hat) # 精度行列
delta_arr = np.tile(np.nan, reps=(line_num, grid_size)) # 受け皿
for l in range(line_num):
    
    # 座標を作成
    beta_arr = np.stack(
        [beta0_grid.flatten(), beta1_grid.flatten(), beta2_vec[l].repeat(grid_size)], 
        axis=1
    )
    
    # マハラノビス距離を計算
    delta_vec = np.array(
        [np.sqrt((beta - beta_true).T @ inv_Sigma_beta_hat @ (beta - beta_true)) for beta in beta_arr]
    )
    
    # 距離を格納
    delta_arr[l] = delta_vec.copy()
print(delta_arr[:5, :5].round(2))
print(delta_arr.shape)

0.01 3.99
0.33 1.67
-1.16 0.16
[0.01 0.09 0.17 0.25 0.33]
[0.33 0.36 0.38 0.41 0.44]
[-1.16 -1.14 -1.13 -1.12 -1.1 ]
(50, 50)
2500
[[4.38 4.32 4.27 4.22 4.17]
 [4.35 4.29 4.24 4.19 4.14]
 [4.32 4.27 4.21 4.16 4.11]
 [4.3  4.24 4.18 4.13 4.08]
 [4.27 4.21 4.16 4.11 4.06]]
(100, 2500)

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

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

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

# グラデーションの範囲を設定
delta_min = 0.0
delta_max = np.ceil(delta_arr.max())

# 等高線の位置を指定
delta_levels = np.linspace(start=delta_min, stop=delta_max, num=8) # 線の数を指定

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

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

# 係数パラメータのマハラノビス距離を作図
fig, ax = plt.subplots(figsize=(9, 8), dpi=100, facecolor='white', constrained_layout=True, 
                       subplot_kw={'projection': '3d'})
cs = ax.contour(beta0_grid, beta1_grid, np.linspace(delta_min, delta_max, num=grid_size).reshape(grid_shape), offset=0.0, 
                cmap='viridis', vmin=delta_min, vmax=delta_max, levels=delta_levels) # カラーバー表示用のダミー
fig.colorbar(cs, ax=ax, shrink=0.6, label='$\\Delta$')
ax.cla() # カラーバー表示用
for l in range(line_num):
    ax.contour(beta0_grid, beta1_grid, delta_arr[l].reshape(grid_shape), offset=beta2_vec[l], 
               cmap='viridis', vmin=delta_min, vmax=delta_max, levels=delta_levels, alpha=0.5, 
               linewidths=1) # マハラノビス距離
ax.scatter(*beta_true, 
           c='red', s=100, 
           label='$\\beta = ('+', '.join([f'{val:.2f}' for val in beta_true])+')$') # 真のパラメータ
ax.scatter(*beta_hat, 
           c='orange', s=100, 
           label='$\\hat{\\beta} = ('+', '.join([f'{val:.2f}' for val in beta_hat])+')$') # 推定パラメータ
ax.plot([beta_true[0], beta_hat[0]], 
        [beta_true[1], beta_hat[1]], 
        [beta_true[2], beta_hat[2]], 
        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]), 
        z=beta_true[2]+0.5*(beta_hat[2]-beta_true[2]), 
        s='$\\Delta_{\\hat{\\beta}}$', size=15, ha='left', va='bottom') # 距離ラベル
ax.set_xlabel('$\\hat{\\beta}_0$')
ax.set_ylabel('$\\hat{\\beta}_1$')
ax.set_zlabel('$\\hat{\\beta}_2$')
ax.set_title(delta_label, loc='left')
fig.suptitle("Mahalanobis' distance", fontsize=20)
ax.legend(loc='upper left')
ax.set_zlim(zmin=beta2_min, zmax=beta2_max)
#ax.set_box_aspect([1.0, beta1_size/beta0_size, beta2_size/beta0_size])
#ax.view_init(elev=90, azim=270) # 0・1軸平面
#ax.view_init(elev=0, azim=270) # 0・2軸平面
#ax.view_init(elev=0, azim=0) # 1・2軸平面
plt.show()

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

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

データ数の影響

　続いて、データ数と推定結果の関係を図で確認します。
　アニメーションの作図コードについては「Introduction-to-SpatialDataScience/code/ch3/2d_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)$ 必要があります。(計算できないデータ数の範囲のグラフもアニメーションに含めています。)

　ここまでの記事では、線形回帰モデルの最小二乗法を確認しました。次からの記事では、モデルの評価を確認します。

参考文献

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

作者:植杉　大
Independently published

Amazon

おわりに

　回帰平面ってどんな感じだろう回帰直線の図と並べてみたいなとちょっと思ってしまったばっかりにとても手間のかかることになりました。まぁでも満足に可視化できたので良かった。
　特に3次元のマハラノビス距離のグラフを作れたのが良かったです。色々な苦労の上のいくつかの小細工の末の図なので忘れる前に記事として書き残しておきたいです。それとつまり同じ方法で3次元のガウス分布も表現できるってことですよね。それも書きたい。
　~~しかしアニメーションの作図コードがこんなに増えるならGitHubで管理すればよかった。~~マハラノビス距離の計算をミスって(ルートをとり忘れて)て作り直すことになったので、その際(半月と経たず)にGitHubに上げてそっちを見てもらうことにしました。

　さて1つ前の記事にも書きましたが、3章の内容はこれで一旦終わりにします。4章の内容は重くなるんですかね、それとも係数パラメータと別に重みパラメータが増えるだけで難易度はそれほど変わらないんですかね。
　2章の途中からは数記事が書き溜まってから予約投稿をしてたのですが、この記事を追加で書いている間にストックが尽きたので、次の記事まで空きます。その前に重み行列の図を作り直したくなってます。その他にもほぼ書けてるのに更新が止まってる記事(本)もあって、、

　2024年8月7日は、BEYOOOOONDSのメジャーデビュー5周年の記念日です。

　5年が経って体制も変化しつつある今日この頃ですが、これからもビヨらしい活動が続いてほしいです。おめでとうございます！！！！！

【次の内容】

つづく

からっぽのしょこ

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

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

はじめに

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

最小二乗法の図

真のモデルの設定

係数パラメータの推定

分散パラメータの推定

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

データ数の影響

係数パラメータの図

分散パラメータの図

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

参考文献

おわりに