はじめに

　『Pythonで学ぶ空間データサイエンス入門』の独学時のまとめノートです。「定義編・導出編」「作図編・実装編」「可視化編」の三部構成でモデルやアルゴリズムの理解を目指します。
　本の内容から寄り道・回り道しながら進めます。本を読んだ上で補助的に読んでください。

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

【前の内容】

www.anarchive-beta.com

【他の内容】

www.anarchive-beta.com

【今回の内容】

はじめに
3.2.1-3 線形回帰モデルの最小二乗法(OLS)の導出
参考文献
おわりに

3.2.1-3 線形回帰モデルの最小二乗法(OLS)の導出

　線形回帰モデル(linear regression model)に対する最小二乗法(OLS・ordinary least square)を導出します。
　線形回帰モデルについては「3.2.0：線形回帰モデルの定義式【空間DS入門のノート】 - からっぽのしょこ」を参照してください。

モデルの設定

　まずは、線形回帰モデルの定義(仮定)を数式で確認します。
　多変量正規分布については「多次元ガウス分布の定義式 - からっぽのしょこ」を参照してください。

定義式

　線形回帰モデルの定義式については「線形回帰モデルの定義式」を参照してください。

　線形回帰モデルは、次の式で定義されます。

$\begin{align} \mathbf{y} &= \mathbf{X} \boldsymbol{\beta} + \boldsymbol{\epsilon} \tag{1}\\ y_n &= \mathbf{x}_n^{\top} \boldsymbol{\beta} + \epsilon_n \\ &= \beta_0 + \sum_{k=1}^K \beta_k x_{nk} + \epsilon_n \end{align}$

　ここで、 $\mathbf{x}_n$ はデータ $n$ の説明変数( $N$ 次元ベクトル)、 $\mathbf{X}$ は $N$ 個のデータの説明変数( $N \times (K+1)$ の行列)、 $y_n$ はデータ $n$ の被説明変数(スカラ)、 $\mathbf{y}$ は $N$ 個のデータの被説明変数( $N$ 次元ベクトル)、 $\boldsymbol{\beta}$ は回帰パラメータ( $K+1$ 次元ベクトル)、 $\epsilon_n$ はデータ $n$ の誤差項(スカラ)、 $\boldsymbol{\epsilon}$ は $N$ 個のデータの誤差項( $N$ 次元ベクトル)です。

　誤差項 $\boldsymbol{\epsilon}$ は、平均ベクトル $\boldsymbol{\mu} = \mathbf{0}$ ・分散共分散行列 $\boldsymbol{\Sigma} = \sigma^2 \mathbf{I}$ の多変量正規分布に従うと仮定します。

$\begin{align} \boldsymbol{\epsilon} &\sim \mathcal{N}(\mathbf{0}, \sigma^2 \mathbf{I}) \tag{2}\\ \epsilon_n &\sim \mathcal{N}(0, \sigma^2) \end{align}$

　ここで、 $\sigma^2$ は分散パラメータ(スカラ)です。

推定値

　線形回帰モデルに関する推定値については「線形回帰モデルの定義式」を参照してください。

　回帰パラメータの推定値 $\hat{\boldsymbol{\beta}}$ を用いて、被説明変数の推定値(予測値) $\hat{\mathbf{y}}$ を求めます。

$\begin{align} \hat{\mathbf{y}} &= \mathbf{X} \hat{\boldsymbol{\beta}} \tag{3}\\ \hat{y}_n &= \mathbf{x}_n^{\top} \hat{\boldsymbol{\beta}} \\ &= \hat{\beta}_0 + \sum_{k=1}^K \hat{\beta}_k x_{nk} \end{align}$

　観測値(実現値) $\mathbf{y}$ と予測値 $\hat{\mathbf{y}}$ の差を残差 $\mathbf{e}$ とします。

$\begin{align} \mathbf{e} &= \mathbf{y} - \hat{\mathbf{y}} \tag{4}\\ e_n &= y_n - \hat{y}_n \end{align}$

計算式

　線形回帰モデルに関する計算式については「線形回帰モデルの定義式」を参照してください。

　残差 $\mathbf{e}$ の式(4)の予測値 $\hat{\mathbf{y}}$ を式(3)で置き換えます。

$\begin{align} \mathbf{e} &= \mathbf{y} - \hat{\mathbf{y}} \tag{4}\\ &= \mathbf{y} - \mathbf{X} \hat{\boldsymbol{\beta}} \tag{4'} \end{align}$

　観測値 $\mathbf{X}, \mathbf{y}$ を用いた式に変形しました。

　残差 $\mathbf{e}$ の内積(残差平方和)は、次の式となります。

$\begin{align} \mathbf{e}^{\top} \mathbf{e} &= (\mathbf{y} - \mathbf{X} \hat{\boldsymbol{\beta}})^{\top} (\mathbf{y} - \mathbf{X} \hat{\boldsymbol{\beta}}) \\ &= \mathbf{y}^{\top} \mathbf{y} - 2 \hat{\boldsymbol{\beta}}^{\top} \mathbf{X}^{\top} \mathbf{y} + \hat{\boldsymbol{\beta}}^{\top} \mathbf{X}^{\top} \mathbf{X} \hat{\boldsymbol{\beta}} \tag{5} \end{align}$

　ハット行列 $\mathbf{S}$ は、次の式で定義されます。

$\displaystyle \mathbf{S} = \mathbf{X} (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^{\top} \tag{6}$

　ハット行列(射影行列)については「3.2.2：ハット行列の性質の導出【空間データサイエンス入門のノート】 - からっぽのしょこ」を参照してください。

　冪等行列 $\mathbf{M}$ は、次の式で定義されます。

$\begin{align} \mathbf{M} &= \mathbf{I} - \mathbf{S} \\ &= \mathbf{I} - \mathbf{X} (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^{\top} \tag{7} \end{align}$

　冪等行列(べき等行列)については「3.2.2：冪等行列の性質の導出【空間データサイエンス入門のノート】 - からっぽのしょこ」を参照してください。

パラメータの推定

　次は、線形回帰モデルに対する最小二乗法によるパラメータの計算式(最小二乗推定量・最小二乗解)を導出します。

回帰パラメータの最小二乗推定量

　回帰パラメータの最小二乗推定量(残差平方和が最小となる推定値)を求めます。

　残差平方和 $\mathbf{e}^{\top} \mathbf{e}$ を回帰パラメータの推定値 $\hat{\boldsymbol{\beta}}$ に関して微分します。

$\begin{align} \frac{ \partial \mathbf{e}^{\top} \mathbf{e} }{ \partial \hat{\boldsymbol{\beta}} } &= \frac{\partial}{\partial \hat{\boldsymbol{\beta}}} \Bigl\{ \mathbf{y}^{\top} \mathbf{y} - 2 \hat{\boldsymbol{\beta}}^{\top} \mathbf{X}^{\top} \mathbf{y} + \hat{\boldsymbol{\beta}}^{\top} \mathbf{X}^{\top} \mathbf{X} \hat{\boldsymbol{\beta}} \Bigr\} \\ &= \frac{\partial}{\partial \hat{\boldsymbol{\beta}}} \Bigl\{ \mathbf{y}^{\top} \mathbf{y} \Bigr\} - 2 \frac{\partial \hat{\boldsymbol{\beta}}^{\top}}{\partial \hat{\boldsymbol{\beta}}} \mathbf{X}^{\top} \mathbf{y} + \frac{ \partial \hat{\boldsymbol{\beta}}^{\top} \mathbf{X}^{\top} \mathbf{X} \hat{\boldsymbol{\beta}} }{ \partial \hat{\boldsymbol{\beta}} } \\ &= \mathbf{0} - 2 \mathbf{I} \mathbf{X}^{\top} \mathbf{y} + 2 \mathbf{X}^{\top} \mathbf{X} \hat{\boldsymbol{\beta}} \\ &= - 2 \mathbf{X}^{\top} \mathbf{y} + 2 \mathbf{X}^{\top} \mathbf{X} \hat{\boldsymbol{\beta}} \tag{8}\\ &= - 2 \mathbf{X}^{\top} \Bigl( \mathbf{y} - \mathbf{X} \hat{\boldsymbol{\beta}} \Bigr) \\ &= - 2 \mathbf{X}^{\top} \mathbf{e} \end{align}$

途中式の途中式(クリックで展開)

1: 式(5)より、 $\mathbf{e}$ の内積の偏微分の式を立てます。 $\hat{\boldsymbol{\beta}}$ に関する微分なので、 $\hat{\boldsymbol{\beta}}$ 以外の項は定数として扱います。
2: 和の微分 $\frac{d \{f(x) + g(x)\}}{d x} = \frac{d f(x)}{d x} + \frac{d g(x)}{d x}$ より、微分の和に分割します。
2: 定数倍の微分 $\frac{d \{a f(x)\}}{d x} = a \frac{f(x)}{d x}$ より、係数を微分の外に出します。
3: 式(9)より、定数スカラの微分は0ベクトルになり消えます。
3: 式(10)より、変数ベクトルの微分は単位行列になり消えます。
3: 対称行列の二次形式の微分 $\frac{\partial \mathbf{x}^{\top} \mathbf{A} \mathbf{x}}{\partial \mathbf{x}} = 2 \mathbf{A} \mathbf{x}$ より、 $\mathbf{X}^{\top} \mathbf{X}$ を1つの行列とみなして変形します。
4: 式を整理します。
5: 共通の項と $-1$ を括り出します。
6: 式(4')より、丸括弧を残差に置き換えます。

　被説明変数 $\mathbf{y}$ の内積を回帰パラメータの推定値 $\hat{\boldsymbol{\beta}}$ に関して微分します。

$\begin{align} \frac{\partial \mathbf{y}^{\top} \mathbf{y}}{\partial \hat{\boldsymbol{\beta}}} &= \begin{pmatrix} \frac{\partial \mathbf{y}^{\top} \mathbf{y}}{\partial \hat{\beta}_1} \\ \frac{\partial \mathbf{y}^{\top} \mathbf{y}}{\partial \hat{\beta}_2} \\ \vdots \\ \frac{\partial \mathbf{y}^{\top} \mathbf{y}}{\partial \hat{\beta}_N} \end{pmatrix} \\ &= \begin{pmatrix} 0 \\ 0 \\ \vdots \\ 0 \end{pmatrix} \\ &= \mathbf{0} \tag{9} \end{align}$

途中式の途中式(クリックで展開)

1: 縦ベクトルによるスカラの微分の要素を明示します。
2: 定数の微分 $\frac{d a}{d x} = 0$ より、定数ベクトルの微分は0ベクトルになります。

　回帰パラメータの推定値 $\hat{\boldsymbol{\beta}}$ の転置を回帰パラメータの推定値 $\hat{\boldsymbol{\beta}}$ に関して微分します。

$\begin{align} \frac{\partial \hat{\boldsymbol{\beta}}^{\top}}{\partial \hat{\boldsymbol{\beta}}} &= \begin{pmatrix} \frac{\partial \hat{\beta}_1}{\partial \hat{\beta}_1} & \frac{\partial \hat{\beta}_2}{\partial \hat{\beta}_1} & \cdots & \frac{\partial \hat{\beta}_N}{\partial \hat{\beta}_1} \\ \frac{\partial \hat{\beta}_1}{\partial \hat{\beta}_2} & \frac{\partial \hat{\beta}_2}{\partial \hat{\beta}_2} & \cdots & \frac{\partial \hat{\beta}_N}{\partial \hat{\beta}_2} \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial \hat{\beta}_1}{\partial \hat{\beta}_N} & \frac{\partial \hat{\beta}_2}{\partial \hat{\beta}_N} & \cdots & \frac{\partial \hat{\beta}_N}{\partial \hat{\beta}_N} \end{pmatrix} \\ &= \begin{pmatrix} 1 & 0 & \cdots & 0 \\ 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 1 \end{pmatrix} \\ &= \mathbf{I} \tag{10} \end{align}$

途中式の途中式(クリックで展開)

1: 縦ベクトルによる横ベクトルの微分の要素を明示します。
2: 定数の微分 $\frac{d a}{d x} = 0$ または変数の微分 $\frac{d x}{d x} = 1$ より、変数ベクトルの微分は単位行列になります。

　対角要素は $\frac{\partial y_n}{\partial y_n} = 1\ (n = 1, 2, \dots, N)$ 、非対角要素は $\frac{\partial y_i}{\partial y_j} = 0\ (i \neq j)$ となり、対角要素のみが残ります。

　残差平方和の微分 $\frac{\partial \mathbf{e}^{\top} \mathbf{e}}{\partial \hat{\boldsymbol{\beta}}}$ が0ベクトル $\mathbf{0}$ となる回帰パラメータの推定値 $\hat{\boldsymbol{\beta}}$ を求めます。

$\begin{align} && - 2 \mathbf{X}^{\top} \mathbf{y} + 2 \mathbf{X}^{\top} \mathbf{X} \hat{\boldsymbol{\beta}} &= \mathbf{0} \\ \Rightarrow && 2 \mathbf{X}^{\top} \mathbf{X} \hat{\boldsymbol{\beta}} &= 2 \mathbf{X}^{\top} \mathbf{y} \\ \Rightarrow && \hat{\boldsymbol{\beta}} &= (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^{\top} \mathbf{y} \equiv \hat{\boldsymbol{\beta}}_{\mathrm{OLS}} \tag{11} \end{align}$

途中式の途中式(クリックで展開)

1: 式(8)より、 $\frac{\partial \mathbf{e}^{\top} \mathbf{e}}{\partial \hat{\boldsymbol{\beta}}} = \mathbf{0}$ とおいて、 $\hat{\boldsymbol{\beta}}$ について解きます。
2: $\hat{\boldsymbol{\beta}}$ と無関係な項を右辺に移します。
3: 両辺に $\mathbf{X}^{\top} \mathbf{X}$ の逆行列を左から掛けます。

　回帰パラメータ $\boldsymbol{\beta}$ の最小二乗推定量 $\hat{\boldsymbol{\beta}}_{\mathrm{OLS}}$ の計算式が得られました。

分散パラメータの不偏推定量

　続いて、残差平方和を用いて、分散パラメータの不偏推定量を求めます。

　冪等行列 $\mathbf{M}$ を用いて、残差 $\mathbf{e}$ の式(4')を変形します。

$\begin{align} \mathbf{e} &= \mathbf{y} - \mathbf{X} \hat{\boldsymbol{\beta}} \tag{4'}\\ &= \mathbf{y} - \mathbf{X} (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^{\top} \mathbf{y} \\ &= (\mathbf{X} \boldsymbol{\beta} + \boldsymbol{\epsilon}) - \mathbf{X} (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^{\top} ( \mathbf{X} \boldsymbol{\beta} + \boldsymbol{\epsilon} ) \\ &= \mathbf{X} \boldsymbol{\beta} + \boldsymbol{\epsilon} - \mathbf{X} (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^{\top} \mathbf{X} \boldsymbol{\beta} - \mathbf{X} (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^{\top} \boldsymbol{\epsilon} \\ &= \mathbf{X} \boldsymbol{\beta} + \boldsymbol{\epsilon} - \mathbf{X} \mathbf{I} \boldsymbol{\beta} - \mathbf{X} (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^{\top} \boldsymbol{\epsilon} \\ &= \Bigl( \mathbf{I} - \mathbf{X} (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^{\top} \Bigr) \boldsymbol{\epsilon} \\ &= \mathbf{M} \boldsymbol{\epsilon} \tag{12} \end{align}$

途中式の途中式(クリックで展開)

1: 式(4')を再掲しています。
2: 式(11)より、 $\hat{\boldsymbol{\beta}}$ に式を代入します。
3: 式(1)より、 $\mathbf{y}$ に式を代入します。
4: 括弧を展開します。
5: 逆行列の性質 $\mathbf{A}^{-1} \mathbf{A} = \mathbf{A} \mathbf{A}^{-1} = \mathbf{I}$ より、 $\mathbf{X}^{\top} \mathbf{X}$ を1つの行列とみなして積を計算すると、単位行列 $\mathbf{I}$ になり消えます。
6: $\boldsymbol{\epsilon}$ を括り出します。
7: 式(7)より、 $\boldsymbol{\epsilon}$ の係数の項を冪等行列 $\mathbf{M}$ に置き換えます。

　ちなみに、ハット行列 $\mathbf{S}$ を用いて、残差 $\mathbf{e}$ の式(4')を変形すると

途中式の途中式(クリックで展開)

1: 式(4')を再掲しています。
2: 式(11)より、 $\hat{\boldsymbol{\beta}}$ に式を代入します。
3: 式(6)より、 $\mathbf{y}$ の係数の項をハット行列 $\mathbf{S}$ に置き換えます。
4: $\mathbf{y}$ を括り出します。

となり、さらに冪等行列 $\mathbf{M}$ を用いて変形すると、同じ式が得られます。

$\begin{align} \mathbf{e} &= \mathbf{M} \mathbf{y} \\ &= \mathbf{M} ( \mathbf{X} \boldsymbol{\beta} + \boldsymbol{\epsilon} ) \\ &= \mathbf{M} \mathbf{X} \boldsymbol{\beta} + \mathbf{M} \boldsymbol{\epsilon} \\ &= \mathbf{0} \boldsymbol{\beta} + \mathbf{M} \boldsymbol{\epsilon} \\ &= \mathbf{M} \boldsymbol{\epsilon} \tag{12} \end{align}$

途中式の途中式(クリックで展開)

1: 式(7)より、 $\mathbf{y}$ の係数の項を冪等行列 $\mathbf{M}$ に置き換えます。
2: 式(1)より、 $\mathbf{y}$ に式を代入します。
3: 括弧を展開します。
4: 冪等行列の性質 $\mathbf{M} \mathbf{X} = \mathbf{0}$ より、変数との積は0行列になり消えます。

　残差 $\mathbf{e}$ の計算式が得られました。
　この式(12)は、残差(観測値と予測値の差) $\mathbf{e}$ と誤差(観測値と理論値の差) $\boldsymbol{\epsilon}$ の関係式と言えます。

　冪等行列 $\mathbf{M}$ を用いて、残差平方和 $\mathbf{e}^{\top} \mathbf{e}$ の式(5)を変形します。

$\begin{align} \mathbf{e}^{\top} \mathbf{e} &= (\mathbf{M} \boldsymbol{\epsilon})^{\top} (\mathbf{M} \boldsymbol{\epsilon}) \\ &= \boldsymbol{\epsilon}^{\top} \mathbf{M}^{\top} \mathbf{M} \boldsymbol{\epsilon} \\ &= \boldsymbol{\epsilon}^{\top} \mathbf{M} \boldsymbol{\epsilon} \tag{13} \end{align}$

途中式の途中式(クリックで展開)

1: 式(11)より、 $\mathbf{e}$ の内積の式を立てます。
2: 転置の性質 $(\mathbf{A} \mathbf{B})^{\top} = \mathbf{B}^{\top} \mathbf{A}^{\top}$ より、括弧を展開します。
3: 冪等行列の性質 $\mathbf{M}^{\top} = \mathbf{M}$ 、 $\mathbf{M} \mathbf{M} = \mathbf{M}$ より、自身との積は1つの $\mathbf{M}$ になります。

　さらに、ベクトル・行列の要素を明示して計算します。

$\begin{align} \mathbf{e}^{\top} \mathbf{e} &= \begin{pmatrix} \epsilon_1 & \epsilon_2 & \cdots & \epsilon_N \end{pmatrix} \begin{pmatrix} m_{11} & m_{12} & \cdots & m_{1N} \\ m_{21} & m_{22} & \cdots & m_{2N} \\ \vdots & \vdots & \ddots & \vdots \\ m_{N1} & m_{N2} & \cdots & m_{NN} \end{pmatrix} \begin{pmatrix} \epsilon_1 \\ \epsilon_2 \\ \vdots \\ \epsilon_N \end{pmatrix} \\ &= \begin{pmatrix} \sum_{i=1}^N \epsilon_i m_{i1} & \sum_{i=1}^N \epsilon_i m_{i2} & \cdots & \sum_{i=1}^N \epsilon_i m_{iN} \end{pmatrix} \begin{pmatrix} \epsilon_1 \\ \epsilon_2 \\ \vdots \\ \epsilon_N \end{pmatrix} \\ &= \sum_{i=1}^N \sum_{j=1}^N \epsilon_i m_{ij} \epsilon_j \tag{13'} \end{align}$

途中式の途中式(クリックで展開)

1: 内積の要素を明示します。
2: 横ベクトルと行列の積を計算します。
3: ベクトルの内積を計算します。

　残差平方和 $\mathbf{e}^{\top} \mathbf{e}$ の計算式が得られました。

　残差平方和 $\mathbf{e}^{\top} \mathbf{e}$ の期待値を求めます。

$\begin{align} \mathbb{E}[\mathbf{e}^{\top} \mathbf{e}] &= \mathbb{E} \Biggl[ \sum_{i=1}^N \sum_{j=1}^N m_{ij} \epsilon_i \epsilon_j \Biggr] \\ &= \sum_{i=1}^N \sum_{j=1}^N \mathbb{E}[m_{ij} \epsilon_i \epsilon_j] \\ &= \sum_{i=1}^N \sum_{j=1}^N m_{ij} \mathbb{E}[\epsilon_i \epsilon_j] \tag{14} \end{align}$

途中式の途中式(クリックで展開)

1: 式(12')より、 $\mathbf{e}$ の内積の期待値の式を立てます。
2: 和の期待値 $\mathbb{E}[x + y] = \mathbb{E}[x] + \mathbb{E}[y]$ より、期待値の和に分割します。

　期待値の総和を、総和の期待値に変形します。

$\displaystyle \begin{aligned} \mathbb{E} \Biggl[ \sum_{i=1}^N \sum_{j=1}^N m_{ij} \epsilon_i \epsilon_j \Biggr] &= \mathbb{E} \Bigl[ m_{11} \epsilon_1 \epsilon_1 + \cdots + m_{ij} \epsilon_i \epsilon_j + \cdots + m_{NN} \epsilon_N \epsilon_N \Bigr] \\ &= \mathbb{E}[m_{1} \epsilon_1 \epsilon_1] + \cdots + \mathbb{E}[m_{ij} \epsilon_i \epsilon_j] + \cdots + \mathbb{E}[m_{NN} \epsilon_N \epsilon_N] \\ &= \sum_{i=1}^N \sum_{j=1}^N \mathbb{E}[m_{ij} \epsilon_i \epsilon_j] \end{aligned}$

3: 定数倍の期待値 $\mathbb{E}[a x] = a \mathbb{E}[x]$ より、 $\boldsymbol{\epsilon}$ と無関係な $\mathbf{M}$ の要素を期待値の外に出します。

　各誤差項 $\epsilon_n$ は平均 $\mu = 0$ の正規分布に従い生成されると仮定しているので、 $\epsilon_n$ の期待値は

$\displaystyle \mathbb{E}[\epsilon_n] = \mu = 0$

であり、 $\epsilon_i, \epsilon_j$ の積の期待値は $\epsilon_i, \epsilon_j$ の共分散 $\sigma_{ij}$ に変形できます。

$\begin{align} \mathbb{E}[\epsilon_i \epsilon_j] &= \mathbb{E}[(\epsilon_i - 0) (\epsilon_j - 0)] \\ &= \mathbb{E}[(\epsilon_i - \mathbb{E}[\epsilon_i]) (\epsilon_j - \mathbb{E}[\epsilon_j])] \\ &= \sigma_{ij} \end{align}$

途中式の途中式(クリックで展開)

1: $\mathbb{E}[\epsilon_i \epsilon_j]$ に関して、 $\epsilon_n = \epsilon_n - 0$ とします。
2: モデルの定義(2)より、 $\epsilon_n$ の期待値は0です。
3: 共分散の定義 $\sigma_{xy} = \mathbb{E}[(x - \mathbb{E}[x]) (y - \mathbb{E}[y])]$ より、式全体を誤差項の共分散に置き換えます。

　各誤差項 $\epsilon_n$ は独立に生成されると仮定しているので、異なるインデックスの要素 $\epsilon_i, \epsilon_j$ の共分散は $\sigma_{ij} = 0$ です。また、各誤差項 $\epsilon_n$ は分散 $\sigma^2$ の正規分布に従い生成されると仮定しているので、同じインデックスの要素 $\epsilon_n, \epsilon_n$ の共分散( $\epsilon_n$ の分散)は $\sigma_{nn} = \sigma^2$ です。

$\displaystyle \mathbb{E}[\epsilon_i \epsilon_j] = \sigma_{ij} = \begin{cases} \sigma^2 &\quad (i = j) \\ 0 &\quad (i \neq j) \end{cases} \tag{15}$

　 $\epsilon_n, \epsilon_n$ の期待値と分散・共分散の関係(15)を用いて、残差平方和 $\mathbf{e}^{\top} \mathbf{e}$ の期待値の式(14)を変形します。

$\begin{align} \mathbb{E}[\mathbf{e}^{\top} \mathbf{e}] &= \sum_{i=1}^N \Bigl\{ m_{i1} \mathbb{E}[\epsilon_i \epsilon_1] + \cdots + m_{ii} \mathbb{E}[\epsilon_i \epsilon_i] + \cdots + m_{iN} \mathbb{E}[\epsilon_i \epsilon_N] \Bigr\} \tag{14}\\ &= \sum_{i=1}^N \Bigl\{ m_{i1} \sigma_{i1} + \cdots + m_{ii} \sigma_{ii} + \cdots + m_{iN} \sigma_{iN} \Bigr\} \\ &= \sum_{i=1}^N \Bigl\{ m_{i1} 0 + \cdots + m_{ii} \sigma^2 + \cdots + m_{iN} 0 \Bigr\} \\ &= \sum_{i=1}^N m_{ii} \sigma^2 \end{align}$

途中式の途中式(クリックで展開)

1: 式(14)の $j$ に関する項の和 $\sum_j$ を展開します。
2: 式(15)より、 $\epsilon_n, \epsilon_n$ の積の期待値は共分散です。
3: 式(15)より、同じインデックス $i, i$ の共分散は分散です。
3: 式(15)より、異なるインデックス $i, j$ の共分散は0です。

　さらに、冪等行列 $\mathbf{M}$ を用いて、式を変形します。

$\begin{align} \mathbb{E}[\mathbf{e}^{\top} \mathbf{e}] &= \sigma^2 \sum_{i=1}^N m_{ii} \\ &= \sigma^2 \mathrm{trace}(\mathbf{M}) \\ &= \sigma^2 \Bigl( N - (K+1) \Bigr) \tag{16} \end{align}$

途中式の途中式(クリックで展開)

1: 各データと無関係な $\sigma^2$ を総和の外に出します。
2: 行列のトレースの定義(対角要素の和) $\mathrm{tr(\mathbf{A})} = \sum_{n=1}^N a_{nn}$ より、 $\mathbf{M}$ の対角要素の和 $\sum_{i=1}^N m_{ii}$ を $\mathbf{M}$ のトレース $\mathrm{tr(\mathbf{M})}$ に置き換えます。
3: 冪等行列の性質 $\mathrm{tr(\mathbf{M})} = N - M$ より、 $\mathbf{M}$ のトレースは $\mathbf{X}$ の行数と列数の差になります。

　残差平方和 $\mathbf{e}^{\top} \mathbf{e}$ の期待値の計算式が得られました。

　残差平方和 $\mathbf{e}^{\top} \mathbf{e}$ の期待値の式(16)を分散パラメータ $\sigma^2$ について解きます。

$\displaystyle \sigma^2 = \frac{ \mathbb{E}[\mathbf{e}^{\top} \mathbf{e}] }{ N - (K+1) } \tag{17}$

　分散パラメータ $\sigma^2$ の計算式が得られました。

　不偏推定量の条件(期待値が真の値と一致する)

$\displaystyle \mathbb{E}[\hat{\sigma}^2] = \sigma^2$

を満たす式として、分散パラメータの推定値 $\hat{\sigma}^2$ を次のようにおきます。

$\displaystyle \hat{\sigma}^2 = \frac{ \mathbf{e}^{\top} \mathbf{e} }{ N - (K+1) } \tag{18}$

　この式の両辺の期待値をとると

$\displaystyle \begin{aligned} \mathbb{E}[\hat{\sigma}^2] &= \mathbb{E} \Biggl[ \frac{ \mathbf{e}^{\top} \mathbf{e} }{ N - (K+1) } \Biggr] \\ &= \frac{ \mathbb{E}[\mathbf{e}^{\top} \mathbf{e}] }{ N - (K+1) } \\ &= \sigma^2 \end{aligned}$

途中式の途中式(クリックで展開)

1: 式(18)より、 $\hat{\sigma}^2$ の期待値の式を立てます。
2: 定数倍の期待値 $\mathbb{E}[a x] = a \mathbb{E}[x]$ より、分母を期待値の外に出します。
3: 式(17)より、式全体を分散パラメータに置き換えます。

不偏推定量の条件を満たすことが分かります。
　分散パラメータの不偏推定量 $\hat{\sigma}^2$ の計算式(18)が得られました。

　以上で、線形回帰モデルのパラメータの推定値の計算式が求まりました。

パラメータの統計量

　最後に、分散パラメータの推定値を評価するための統計量を導出します。

回帰パラメータの分散共分散行列

　回帰パラメータの最小二乗推定量の分散共分散行列を求めます。

　誤差項 $\boldsymbol{\epsilon}$ の期待値を計算します。

$\begin{align} \mathbb{E}[\boldsymbol{\epsilon}] &= \mathbb{E} \left[ \begin{pmatrix} \epsilon_1 \\ \epsilon_2 \\ \vdots \\ \epsilon_N \end{pmatrix} \right] \\ &= \begin{pmatrix} \mathbb{E}[\epsilon_1] \\ \mathbb{E}[\epsilon_2] \\ \vdots \\ \mathbb{E}[\epsilon_N] \end{pmatrix} \\ &= \begin{pmatrix} 0 \\ 0 \\ \vdots \\ 0 \end{pmatrix} \\ &= \mathbf{0} \tag{19} \end{align}$

途中式の途中式(クリックで展開)

1: ベクトルの要素を明示します。
2: ベクトルの期待値を、期待値のベクトルに変形します。
3: モデルの定義(2)より、誤差項の期待値は $\mathbb{E}[\epsilon_n] = \mu = 0$ です。

　誤差項 $\boldsymbol{\epsilon}$ の期待値の計算式が得られました。0ベクトル $\mathbf{0}$ になります。

　誤差項 $\boldsymbol{\epsilon}$ の分散共分散行列を計算します。

$\displaystyle \begin{aligned} \mathbb{E} \Bigl[ (\boldsymbol{\epsilon} - \mathbb{E}[\boldsymbol{\epsilon}]) (\boldsymbol{\epsilon} - \mathbb{E}[\boldsymbol{\epsilon}])^{\top} \Bigr] &= \mathbb{E} \Bigl[ (\boldsymbol{\epsilon} - \mathbf{0}) (\boldsymbol{\epsilon} - \mathbf{0})^{\top} \Bigr] \\ &= \mathbb{E}[\boldsymbol{\epsilon} \boldsymbol{\epsilon}^{\top}] \end{aligned}$

途中式の途中式(クリックで展開)

1: 分散共分散行列の定義 $\boldsymbol{\Sigma}_{\mathbf{x}\mathbf{y}} = \mathbb{E}[(\mathbf{x} - \mathbb{E}[\mathbf{x}]) (\mathbf{y} - \mathbb{E}[\mathbf{y}])^{\top}]$ より、 $\boldsymbol{\epsilon}$ の分散共分散行列の式を立てます。
1: 式(19)より、 $\boldsymbol{\epsilon}$ の期待値は0ベクトルになり消えます。

　さらに、行列の要素を明示して計算します。

$\begin{align} \mathbb{E} \Bigl[ (\boldsymbol{\epsilon} - \mathbb{E}[\boldsymbol{\epsilon}]) (\boldsymbol{\epsilon} - \mathbb{E}[\boldsymbol{\epsilon}])^{\top} \Bigr] &= \mathbb{E} \left[ \begin{pmatrix} \epsilon_1 \\ \epsilon_2 \\ \vdots \\ \epsilon_N \end{pmatrix} \begin{pmatrix} \epsilon_1 & \epsilon_2 & \cdots & \epsilon_N \end{pmatrix} \right] \\ &= \mathbb{E} \left[ \begin{pmatrix} \epsilon_1 \epsilon_1 & \epsilon_1 \epsilon_2 & \cdots & \epsilon_1 \epsilon_N \\ \epsilon_2 \epsilon_1 & \epsilon_2 \epsilon_2 & \cdots & \epsilon_2 \epsilon_N \\ \vdots & \vdots & \ddots & \vdots \\ \epsilon_N \epsilon_1 & \epsilon_N \epsilon_2 & \cdots & \epsilon_N \epsilon_N \end{pmatrix} \right] \\ &= \begin{pmatrix} \mathbb{E}[\epsilon_1 \epsilon_1] & \mathbb{E}[\epsilon_1 \epsilon_2] & \cdots & \mathbb{E}[\epsilon_1 \epsilon_N] \\ \mathbb{E}[\epsilon_2 \epsilon_1] & \mathbb{E}[\epsilon_2 \epsilon_2] & \cdots & \mathbb{E}[\epsilon_2 \epsilon_N] \\ \vdots & \vdots & \ddots & \vdots \\ \mathbb{E}[\epsilon_N \epsilon_1] & \mathbb{E}[\epsilon_N \epsilon_2] & \cdots & \mathbb{E}[\epsilon_N \epsilon_N] \end{pmatrix} \\ &= \begin{pmatrix} \sigma^2 & 0 & \cdots & 0 \\ 0 & \sigma^2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \sigma^2 \end{pmatrix} \\ &= \sigma^2 \begin{pmatrix} 1 & 0 & \cdots & 0 \\ 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 1 \end{pmatrix} \\ &= \sigma^2 \mathbf{I} \tag{20} \end{align}$

途中式の途中式(クリックで展開)

1: ベクトルの要素を明示します。
2: ベクトルの外積を計算します。
3: 行列の期待値を、期待値の行列に変形します。
4: 式(15)より、誤差項の分散は $\mathbb{E}[\epsilon_n \epsilon_n] = \sigma^2$ 、共分散は $\mathbb{E}[\epsilon_i \epsilon_j] = \sigma_{ij} = 0$ です。
5: 行列の定数倍を、スカラと行列の積に変形します。

　誤差項 $\boldsymbol{\epsilon}$ の分散共分散行列の計算式が得られました。分散パラメータ $\sigma^2$ と単位行列 $\mathbf{I}$ の積になります。

　回帰パラメータの推定値 $\hat{\boldsymbol{\beta}}$ の式(11)を変形します。

$\begin{align} \hat{\boldsymbol{\beta}} &= (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^{\top} \mathbf{y} \tag{11}\\ &= (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^{\top} ( \mathbf{X} \boldsymbol{\beta} + \boldsymbol{\epsilon} ) \\ &= (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^{\top} \mathbf{X} \boldsymbol{\beta} + (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^{\top} \boldsymbol{\epsilon} \\ &= \boldsymbol{\beta} + (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^{\top} \boldsymbol{\epsilon} \tag{21} \end{align}$

途中式の途中式(クリックで展開)

1: 式(11)を再掲しています。
2: 式(1)より、 $\mathbf{y}$ に式を代入します。
3: 括弧を展開します。
4: 逆行列の性質 $\mathbf{A}^{-1} \mathbf{A} = \mathbf{A} \mathbf{A}^{-1} = \mathbf{I}$ より、 $\mathbf{X}^{\top} \mathbf{X}$ を1つの行列とみなして積を計算すると、単位行列 $\mathbf{I}$ になり消えます。

　回帰パラメータの推定値 $\hat{\boldsymbol{\beta}}$ の期待値をとります。

$\begin{align} \mathbb{E}[\hat{\boldsymbol{\beta}}] &= \mathbb{E} \Bigl[ \boldsymbol{\beta} + (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^{\top} \boldsymbol{\epsilon} \Bigr] \\ &= \mathbb{E}[\boldsymbol{\beta}] + \mathbb{E} \Bigl[ (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^{\top} \boldsymbol{\epsilon} \Bigr] \\ &= \boldsymbol{\beta} + (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^{\top} \mathbb{E}[\boldsymbol{\epsilon}] \\ &= \boldsymbol{\beta} + (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^{\top} \mathbf{0} \\ &= \boldsymbol{\beta} \tag{22} \end{align}$

途中式の途中式(クリックで展開)

1: 式(21)より、 $\hat{\boldsymbol{\beta}}$ の期待値の式を立てます。
2: 和の期待値 $\mathbb{E}[x + y] = \mathbb{E}[x] + \mathbb{E}[y]$ より、期待値の和に分割します。
3: 定数倍の期待値 $\mathbb{E}[a x] = a \mathbb{E}[x]$ より、 $\boldsymbol{\epsilon}$ と無関係な $\mathbf{X}$ を期待値の外に出します。
4: 式(19)より、 $\boldsymbol{\epsilon}$ の期待値は0ベクトルになり消えます。

　回帰パラメータの推定値の期待値 $\mathbb{E}[\hat{\boldsymbol{\beta}}]$ と真の値 $\boldsymbol{\beta}$ が一致することから、不偏推定量の条件を満たしていることが分かります。

　回帰パラメータの推定値 $\hat{\boldsymbol{\beta}}$ の偏差を計算します。

$\begin{align} \hat{\boldsymbol{\beta}} - \mathbb{E}[\hat{\boldsymbol{\beta}}] &= \hat{\boldsymbol{\beta}} - \boldsymbol{\beta} \\ &= \boldsymbol{\beta} + (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^{\top} \boldsymbol{\epsilon} - \boldsymbol{\beta} \\ &= (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^{\top} \boldsymbol{\epsilon} \tag{23} \end{align}$

途中式の途中式(クリックで展開)

1: 偏差の定義 $x - \mathbb{E}[x]$ より、 $\hat{\boldsymbol{\beta}}$ の偏差の式を立てます。
1: 式(22)より、 $\hat{\boldsymbol{\beta}}$ の期待値は真の値になります。
2: 式(21)より、 $\hat{\boldsymbol{\beta}}$ に式を代入します。

　回帰パラメータの推定値 $\hat{\boldsymbol{\beta}}$ の分散共分散行列を求めます。

$\begin{align} \mathbb{E} \Bigl[ (\hat{\boldsymbol{\beta}} - \mathbb{E}[\hat{\boldsymbol{\beta}}]) (\hat{\boldsymbol{\beta}} - \mathbb{E}[\hat{\boldsymbol{\beta}}])^{\top} \Bigr] &= \mathbb{E} \Bigl[ (\hat{\boldsymbol{\beta}} - \boldsymbol{\beta}) (\hat{\boldsymbol{\beta}} - \boldsymbol{\beta})^{\top} \Bigr] \\ &= \mathbb{E} \Bigl[ \Bigl( (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^{\top} \boldsymbol{\epsilon} \Bigr) \Bigl( (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^{\top} \boldsymbol{\epsilon} \Bigr)^{\top} \Bigr] \\ &= \mathbb{E} \Bigl[ (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^{\top} \boldsymbol{\epsilon} \boldsymbol{\epsilon}^{\top} \mathbf{X} (\mathbf{X}^{\top} \mathbf{X})^{-1} \Bigr] \\ &= (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^{\top} \mathbb{E}[\boldsymbol{\epsilon} \boldsymbol{\epsilon}^{\top}] \mathbf{X} (\mathbf{X}^{\top} \mathbf{X})^{-1} \\ &= (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^{\top} \sigma^2 \mathbf{I} \mathbf{X} (\mathbf{X}^{\top} \mathbf{X})^{-1} \\ &= \sigma^2 (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^{\top} \mathbf{X} (\mathbf{X}^{\top} \mathbf{X})^{-1} \\ &= \sigma^2 (\mathbf{X}^{\top} \mathbf{X})^{-1} \tag{24} \end{align}$

途中式の途中式(クリックで展開)

1: 分散共分散行列の定義 $\boldsymbol{\Sigma}_{\mathbf{x}\mathbf{y}} = \mathbb{E}[(\mathbf{x} - \mathbb{E}[\mathbf{x}]) (\mathbf{y} - \mathbb{E}[\mathbf{y}])^{\top}]$ より、 $\hat{\boldsymbol{\beta}}$ の分散共分散行列の式を立てます。
1-2: 式(23)より、 $\hat{\boldsymbol{\beta}}$ の偏差に式を代入します。
3: 転置の性質 $(\mathbf{A} \mathbf{B} \mathbf{C})^{\top} = \mathbf{C}^{\top} \mathbf{B}^{\top} \mathbf{A}^{\top}$ 、 $(\mathbf{A}^{\top})^{\top} = \mathbf{A}$ 、転置行列と逆行列の関係 $(\mathbf{A}^{-1})^{\top} = (\mathbf{A}^{\top})^{-1}$ より、括弧を展開します。

　転置により、次のように変形します。

$\displaystyle \begin{aligned} \Bigl( (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^{\top} \boldsymbol{\epsilon} \Bigr)^{\top} &= \boldsymbol{\epsilon}^{\top} (\mathbf{X}^{\top})^{\top} \Bigl( (\mathbf{X}^{\top} \mathbf{X})^{-1} \Bigr)^{\top} \\ &= \boldsymbol{\epsilon}^{\top} \mathbf{X} \Bigl( (\mathbf{X}^{\top} \mathbf{X})^{\top} \Bigr)^{-1} \\ &= \boldsymbol{\epsilon}^{\top} \mathbf{X} \Bigl( \mathbf{X}^{\top} (\mathbf{X}^{\top})^{\top} \Bigr)^{-1} \\ &= \boldsymbol{\epsilon}^{\top} \mathbf{X} (\mathbf{X}^{\top} \mathbf{X})^{-1} \end{aligned}$

4: 定数倍の期待値 $\mathbb{E}[a x] = a \mathbb{E}[x]$ より、 $\boldsymbol{\epsilon}$ と無関係な $\mathbf{X}$ を期待値の外に出します。
5: 式(20)より、 $\mathbb{E}[\boldsymbol{\epsilon} \boldsymbol{\epsilon}^{\top}]$ に式を代入します。
6: スカラ $\sigma^2$ を前に出します。
7: 逆行列の性質 $\mathbf{A}^{-1} \mathbf{A} = \mathbf{A} \mathbf{A}^{-1} = \mathbf{I}$ より、 $\mathbf{X}^{\top} \mathbf{X}$ を1つの行列とみなして積を計算すると、単位行列 $\mathbf{I}$ になり消えます。

　 $\hat{\boldsymbol{\beta}}$ の分散共分散行列の計算式が得られました。

　以上で、最小二乗法によるパラメータ推定の式が得られました。

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

参考文献

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

作者:植杉　大
Independently published

Amazon

線形回帰分析 (統計ライブラリー)

作者:蓑谷千凰彦
朝倉書店

Amazon

おわりに

　これまでにいくつかの回帰モデルを扱ってきましたが、ナントカ回帰じゃないシンプルなモデルはやってなかったかもしれません。いやこのモデルは重回帰なので、よりシンプルな単回帰をまだやってないのですが。(忘れてるだけで扱ってるかもしれませんが…)
　シンプルなら簡単かというとそんなことはなく、行列計算が普通に重かったです。1年以上ぶりにthe matrix cookbookを引くことになりました。二次形式の微分を公式的に書いてしまったのが心残りです。行列の性質シリーズを導出してみたいなぁ(と思うこと3年ほど？)、ネタ切れしたらやるかも。

　このシリーズの想定読者が定まっていなかったのですが、軽い行列計算はできるレベルということでいいですかね。

　2024年8月4日は、私立恵比寿中学の結成15周年の記念日です。

　15年前に自分が何をしていたかぱっと思い出せませんが、結成当初(かは分かりませんがももクロの妹分と呼ばれてた(よね？)頃)には少し聞いたことがあったはずなんですよね。当時の私はアイドル自体に興味がなかったのですが、10年以上の時を経てまさかライブに行くほどハマることになるとは、人生何が起こるか分かりません。今後どうなるかは分かりませんが、とりあえず年度末のライブのチケットを申し込んだので今はそれが楽しみです(次こそエニエニをお願いします！)。
　おめでとうございます！これからも永遠に中学生♪

2026.8.05：加筆修正の際に記事から「線形回帰モデルの定義式」を分割しました。

　訳あって想定読者をいつもよりも低めにするか迷いながら初稿を書き始めていたのですが、どの道この記事辺りから付いてこれなくなるんでした。そんなことは忘れて、結局いつも通りのレベルで書き直しています。訳とやらも不要そうですし。
　あとがきを読んで気付きましたが、気まぐれに書いたGMLの導出で単回帰を扱ったことになるんですね。素のモデルではなく変則的で、思っていた形と違いますが。
　幸いにもネタ切れしていないので行列の公式シリーズは未だに取り組めていません。やりたい気持ちはまだありますよ。

　ところでもう1日早く更新しておけば丁度2年でした。半月前には書き終わっていたのですが、これだけ更新するのを忘れていただけなんですが。

　その後も何度かえびちゅうのライブに行ったもののまだエニエニは聴けていません。迷った末に行かなかった地元でのライブで日替わり曲として歌われるというニアミスはありました。無念。
　でも、夢にも思わなかったソラシド～ねえねえ～のカバーを聴けるという僥倖もありました。何があるか分からないものですね。

【次の内容】

　線形回帰モデルのOLSをプログラムで確認します。

www.anarchive-beta.com

　GWRのOLSを数式で確認します。

www.anarchive-beta.com

　SLMの最尤推定を数式で確認します。

からっぽのしょこ

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

3.2.1-3：線形回帰モデルの最小二乗法(OLS)の導出【空間DS入門のノート】

はじめに

3.2.1-3 線形回帰モデルの最小二乗法(OLS)の導出

モデルの設定

定義式

推定値

計算式

パラメータの推定

回帰パラメータの最小二乗推定量

分散パラメータの不偏推定量

パラメータの統計量

回帰パラメータの分散共分散行列

参考文献

おわりに