はじめに

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

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

【前の内容】

www.anarchive-beta.com

【他の内容】

www.anarchive-beta.com

【今回の内容】

はじめに
3.2.1-3 線形回帰モデルの最小二乗法の実装
参考文献
おわりに

3.2.1-3 線形回帰モデルの最小二乗法の実装

　線形回帰モデル(linear regression model)に対する最小二乗法(OLS・ordinary least square)におけるパラメータの推定処理を実装します。
　線形回帰モデルの定義や計算式については「3.2.1-3：線形回帰モデルの最小二乗法の導出【空間データサイエンス入門のノート】 - からっぽのしょこ」を参照してください。

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

# ライブラリを読込
import numpy as np

疑似データの生成

　まずは、線形回帰モデルにおける観測データを人工的に作成します。

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

# データ数を指定
N = 10

# 係数パラメータ数を指定
K = 5

# 説明変数の値を生成
dat = np.random.uniform(low=-10.0, high=10.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)

[[ 3.75  8.4   3.7  -9.44  3.99]
 [-6.14  6.82  8.14  3.18 -0.98]
 [-8.29  0.42  4.45 -7.92  7.2 ]
 [ 9.84  5.7   1.37  3.96  5.95]
 [-5.56  3.54  6.37 -9.78  7.27]]
[[ 1.    3.75  8.4   3.7  -9.44  3.99]
 [ 1.   -6.14  6.82  8.14  3.18 -0.98]
 [ 1.   -8.29  0.42  4.45 -7.92  7.2 ]
 [ 1.    9.84  5.7   1.37  3.96  5.95]
 [ 1.   -5.56  3.54  6.37 -9.78  7.27]]
(10, 6)

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

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

# 真の係数パラメータを指定:(範囲を指定)
#beta_true = np.random.uniform(low=-1.0, high=1.0, size=K+1)      # (横ベクトルで扱う場合)
beta_true = np.random.uniform(low=-1.0, high=1.0, size=(K+1, 1)) # (縦ベクトルで扱う場合)
print(beta_true.round(2))

[[-0.04]
 [-0.4 ]
 [ 0.42]
 [ 0.17]
 [ 0.67]
 [-0.87]]

　真の係数パラメータ $\boldsymbol{\beta} = (\beta_0, \beta_1, \cdots, \beta_K)$ を指定します。
　この例では、連続一様分布に従う乱数を使います。random モジュールの uniform() のサンプルサイズの引数 size にパラメータの要素数 $K+1$ を指定します。

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

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

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

# 誤差項を生成
#epsilon = np.random.normal(loc=0.0, scale=sigma_true, size=N)      # (横ベクトルで扱う場合)
epsilon = np.random.normal(loc=0.0, scale=sigma_true, size=(N, 1)) # (縦ベクトルで扱う場合)
print(epsilon[:5].round(2))
print(epsilon.shape)

0.25
[[ 0.88]
 [ 0.87]
 [ 0.62]
 [-0.14]
 [ 0.08]]
(10, 1)

　誤差項の真の標準偏差パラメータ $\sigma$ を指定して、真の分散パラメータ $\sigma^2$ を計算します。
　平均 $\mu = 0$ 、分散 $\sigma^2$ の正規分布に従って、 $N$ 個の誤差項 $\boldsymbol{\epsilon} = (\epsilon_1, \cdots, \epsilon_N)$ を作成します。
　正規分布に従う乱数は random モジュールの normal() で生成できます。平均の引数 loc に $0$ 、標準偏差の引数 scale に標準偏差パラメータ $\sigma$ 、サンプルサイズの引数 size にデータ数 $N$ を指定します。

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

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

[[-6.28]
 [10.59]
 [-6.75]
 [-4.  ]
 [-8.04]]
(10, 1)

　被説明変数 $\mathbf{y} = (y_1, \cdots, y_N)$ を次の式で計算します。

$\displaystyle \mathbf{y} = \mathbf{X} \boldsymbol{\beta} + \boldsymbol{\epsilon}$

　係数パラメータ beta_true と誤差項 epsilon の形状によって、被説明変数 y の形状が (N,) の1次元配列(横ベクトル)または (N, 1) の2次元配列(縦ベクトル)になります。
　どちらの形状でも処理できるように実装します。

　以上で、計算に利用するデータを用意できました。

最小二乗法の計算

　次は、線形回帰モデルに対するOLSによるパラメータ推定の処理を確認します。

係数パラメータの推定

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

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

[[ 0.06]
 [-0.42]
 [ 0.47]
 [ 0.21]
 [ 0.65]
 [-0.91]]
(6, 1)

　係数パラメータの推定値 $\hat{\boldsymbol{\beta}} = (\hat{\beta}_0, \hat{\beta}_1, \cdots, \hat{\beta}_K)$ を次の式で計算します。

$\displaystyle \hat{\boldsymbol{\beta}} = (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^{\top} \mathbf{y}$

　逆行列は linalg モジュールの inv() で計算できます。

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

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

[[-6.59]
 [10.55]
 [-7.  ]
 [-3.92]
 [-7.57]]
(10, 1)

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

$\displaystyle \hat{\mathbf{y}} = \mathbf{X} \hat{\boldsymbol{\beta}}$

　残差を作成します。

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

[[ 0.3 ]
 [ 0.05]
 [ 0.25]
 [-0.08]
 [-0.47]]
(10, 1)

　残差 $\mathbf{e} = (e_1, \cdots, e_N)$ を次の式で計算します。

$\displaystyle \mathbf{e} = \mathbf{y} - \hat{\mathbf{y}}$

分散パラメータの推定

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

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

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

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

[[0.70734406]]
[[0.17683601]]
[[0.42051874]]

　分散パラメータの推定値 $\hat{\sigma}^2$ と標準偏差パラメータの推定値 $\hat{\sigma}$ を次の式で計算します。

$\displaystyle \begin{aligned} \hat{\sigma}^2 &= \frac{ \mathbf{e}^{\top} \mathbf{e} }{ N - (K+1) } \\ \hat{\sigma} &= \sqrt{\hat{\sigma}^2} \end{aligned}$

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

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

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

[[ 0.047 -0.005  0.001 -0.006  0.004  0.   ]
 [-0.005  0.001 -0.     0.001 -0.001 -0.   ]
 [ 0.001 -0.     0.001 -0.     0.     0.   ]
 [-0.006  0.001 -0.     0.002 -0.001 -0.   ]
 [ 0.004 -0.001  0.    -0.001  0.001  0.001]
 [ 0.    -0.     0.    -0.     0.001  0.001]]
(6, 6)

　係数パラメータの推定値 $\hat{\boldsymbol{\beta}}$ の分散共分散行列を次の式で計算します。

$\displaystyle \mathbb{E} \Bigl[ (\hat{\boldsymbol{\beta}} - \mathbb{E}[\hat{\boldsymbol{\beta}}]) (\hat{\boldsymbol{\beta}} - \mathbb{E}[\hat{\boldsymbol{\beta}}])^{\top} \Bigr] = \sigma^2 (\mathbf{X}^{\top} \mathbf{X})^{-1}$

　以上の処理を実装します。

最小二乗法の実装

　続いて、線形回帰モデルに対するOLSによるパラメータ推定を行うクラスを作成します。

　OLSによるパラメータ推定を自作クラスとして定義します。

# 線形回帰モデルに対するOLSを実装
class OLS:

    # データの設定
    def __init__(self, X, y):
        
        # 観測データを格納
        self.X = X # 説明変数
        self.y = y # 被説明変数
        self.N = X.shape[0] # 観測データ数
        self.K = X.shape[1] - 1 # 係数パラメータ数

    # 係数パラメータの推定
    def get_beta_hat(self):
        
        # 観測データを取得
        X = self.X
        y = self.y.flatten() # (横ベクトルで扱う場合)
        
        # 係数パラメータの推定値を計算
        beta_hat = np.linalg.inv(X.T @ X) @ X.T @ y
        
        return beta_hat
    
    # 誤差項の分散パラメータの推定
    def get_sigma2_hat(self):
        
        # 観測データを取得
        X = self.X
        y = self.y.flatten() # (横ベクトルで扱う場合)

        # 係数パラメータの推定値を計算
        beta_hat = self.get_beta_hat()

        # 被説明変数の推定値を計算
        y_hat = X @ beta_hat
        
        # 残差を計算
        e = y - y_hat
        
        # 残差平方和を計算
        rss = e.T @ e
        
        # 分散パラメータの推定値を計算
        sigma2_hat = rss / (self.N - self.K - 1)
        
        return sigma2_hat
    
    # 係数パラメータの推定値の分散共分散行列の推定
    def get_sigma_mat_hat(self):

        # 観測データを取得
        X = self.X

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

        return Sigma_beta

　出力変数 y が (N, 1) の2次元配列(縦ベクトル)の場合、.flatten() で (N,) の1次元配列(横ベクトル)に変形して処理します。y が (N,) の1次元配列の場合はそのまま処理されます。
　y を縦ベクトルで扱う( .flatten() を使わない)場合は、(残差平方和 rss や)分散パラメータの推定値 sigma2_hat が (1, 1) の2次元配列になるので、.item() でスカラに変換して出力する処理を追加する必要があります。

　実装したクラスを使ってOLSによるパラメータ推定の計算をします。

# OLSのインスタンスを作成
ols = OLS(X, y)

# 係数パラメータの推定値を計算
beta_hat = ols.get_beta_hat()
print(beta_hat.round(2))
print(beta_hat.shape)

# 分散パラメータの推定値を計算
sigma2_hat = ols.get_sigma2_hat()
print(sigma2_hat)

# 係数パラメータの推定値の分散共分散行列を計算
Sigma_beta_hat = ols.get_sigma_mat_hat()
print(Sigma_beta_hat.round(3))
print(Sigma_beta_hat.shape)

[ 0.06 -0.42  0.47  0.21  0.65 -0.91]
(6,)
0.1768360142975876
[[ 0.047 -0.005  0.001 -0.006  0.004  0.   ]
 [-0.005  0.001 -0.     0.001 -0.001 -0.   ]
 [ 0.001 -0.     0.001 -0.     0.     0.   ]
 [-0.006  0.001 -0.     0.002 -0.001 -0.   ]
 [ 0.004 -0.001  0.    -0.001  0.001  0.001]
 [ 0.    -0.     0.    -0.     0.001  0.001]]
(6, 6)

　以上で、OLSによるパラメータ推定を実装できました。

　この記事では、線形回帰モデルの最小二乗法をプログラムで確認しました。次の記事では、図で確認します。

参考文献

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

作者:植杉　大
Independently published

Amazon

おわりに

　一応この記事は準備で、次の記事がメインです。この記事だけで理解できる人はわざわざこのブログを読んでないと思うので、次の記事の図もセットで読んでください。

【次の内容】

からっぽのしょこ

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

【Python】3.2.1-3：線形回帰モデルの最小二乗法の実装【空間データサイエンス入門のノート】