からっぽのしょこ

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

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

はじめに

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

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

【前の内容】

www.anarchive-beta.com

【他の内容】

www.anarchive-beta.com

【今回の内容】

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

 線形回帰モデル(linear regression model)に対する最小二乗法(OLS・ordinary least square)におけるパラメータの計算式を導出します。

モデルの設定

 まずは、線形回帰モデルの定義を数式で確認します。

  N 個のデータについて、 i 番目の説明変数(入力データ)を  K+1 次元ベクトル  \mathbf{x}_i とします。

 \displaystyle
\mathbf{x}_i
    = \begin{pmatrix}
          x_{i0} \\ x_{i1} \\ x_{i2} \\ \vdots \\ x_{iK}
      \end{pmatrix}

  N 個の説明変数をまとめて、 N \times (K+1) の行列  \mathbf{X} とします。

 \displaystyle
\mathbf{X}
    = \begin{pmatrix}
          \mathbf{x}_1^{\top} \\
          \mathbf{x}_2^{\top} \\
          \vdots \\
          \mathbf{x}_N^{\top}
      \end{pmatrix}
    = \begin{pmatrix}
          x_{10} & x_{11} & x_{12} & \cdots & x_{1K} \\
          x_{20} & x_{21} & x_{22} & \cdots & x_{2K} \\
          \vdots & \vdots & \vdots & \ddots & \vdots \\
          x_{N0} & x_{N1} & x_{N2} & \cdots & x_{NK}
      \end{pmatrix}
    = \begin{pmatrix}
          1      & x_{11} & x_{12} & \cdots & x_{1K} \\
          1      & x_{21} & x_{22} & \cdots & x_{2K} \\
          \vdots & \vdots & \vdots & \ddots & \vdots \\
          1      & x_{N1} & x_{N2} & \cdots & x_{NK}
      \end{pmatrix}

 各変数に関して、 0 番(列)目の要素を  x_{i0} = 1 とします。(インデックスを0から割り当てているので)  K+1 次元(列)です。

  k 番目の要素  x_{ik} に対する真の係数パラメータを  \beta_k として、 K+1 個の真の係数パラメータをまとめて、 K+1 次元ベクトル  \boldsymbol{\beta} とします。

 \displaystyle
\boldsymbol{\beta}
    = \begin{pmatrix}
          \beta_0 \\ \beta_1 \\ \beta_2 \\ \vdots \\ \beta_K
      \end{pmatrix}

  0 番目の要素  \beta_0 は定数項に対応します。

  i 番目のデータの誤差項を  \epsilon_i として、 N 個のデータの誤差項をまとめて、 N 次元ベクトル  \boldsymbol{\epsilon} とします。

 \displaystyle
\boldsymbol{\epsilon}
    = \begin{pmatrix}
          \epsilon_1 \\ \epsilon_2 \\ \vdots \\ \epsilon_N
      \end{pmatrix}

 各誤差  \epsilon_i は、平均が  \mu = 0 で分散が  \sigma^2 (標準偏差が  \sigma )の正規分布に従って、独立に生成されると仮定します。

 \displaystyle
\epsilon_i
    \sim
      \mathcal{N}(0, \sigma^2)

  i 番目と  j 番目の誤差項  \epsilon_i, \epsilon_j の共分散を  \sigma_{ij} とすると、誤差項間は無相関(無関係に生成される)なので、 \sigma_{ij} = 0 です。各誤差項自身との共分散は分散  \sigma_{ii} = \sigma_{jj} = \sigma^2 です。

  i 番目の被説明変数(出力データ)を  y_i として、 N 個の被説明変数をまとめて、 N 次元ベクトル  \mathbf{y} とします。

 \displaystyle
\mathbf{y}
    = \begin{pmatrix}
          y_1 \\ y_2 \\ \vdots \\ y_N
      \end{pmatrix}


 (1つのデータにおける)線形回帰モデルは、次の式で定義されます。

 \displaystyle
\begin{aligned}
y_i
   &= \mathbf{x}_i^{\top}
      \boldsymbol{\beta}
      + \epsilon_i
\\
   &= \begin{pmatrix}
          x_{i0} & x_{i1} & x_{i2} & \cdots & x_{iK}
      \end{pmatrix}
      \begin{pmatrix}
          \beta_0 \\ \beta_1 \\ \beta_2 \\ \vdots \\ \beta_K
      \end{pmatrix}
      + \epsilon_i
\\
   &= x_{i0} \beta_0
      + x_{i1} \beta_1
      + x_{i2} \beta_2
      + \cdots
      + x_{iK} \beta_K
      + \epsilon_i
\\
   &= \sum_{k=0}^K
          \beta_k x_{ik}
      + \epsilon_i
\\
   &= \beta_0
      + \sum_{k=1}^K
          \beta_k x_{ik}
      + \epsilon_i
\end{aligned}

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


  • 1: 説明変数  \mathbf{x}_i とパラメータ  \boldsymbol{\beta} の内積と誤差項  \epsilon_i の和の式を立てます。
  • 2: ベクトル  \mathbf{x}_i, \boldsymbol{\beta} の要素を明示します。
  • 3: ベクトルの内積(要素ごとの積の和)を計算します。
  • 4:  K+1 個の項の和を  \sum_{k=0}^K でまとめます。
  • 5: (項の形が異なる)  0 番目の項  \beta_0 x_{i0} = \beta_0 を取り出し(  K 個の要素の和を  \sum_{k=1}^K でまとめ)ます。

  N 個のデータをまとめると、次の式で計算できます。

 \displaystyle
\begin{align}
&&
\mathbf{y}
   &= \mathbf{X}
      \boldsymbol{\beta}
      + \boldsymbol{\epsilon}
\tag{1}\\
\Leftrightarrow &&
\begin{pmatrix}
    y_1 \\ y_2 \\ \vdots \\ y_N
\end{pmatrix}
   &= \begin{pmatrix}
          1      & x_{11} & x_{12} & \cdots & x_{1K} \\
          1      & x_{21} & x_{22} & \cdots & x_{2K} \\
          \vdots & \vdots & \vdots & \ddots & \vdots \\
          1      & x_{N1} & x_{N2} & \cdots & x_{NK}
      \end{pmatrix}
      \begin{pmatrix}
          \beta_0 \\ \beta_1 \\ \beta_2 \\ \vdots \\ \beta_K
      \end{pmatrix}
      + \begin{pmatrix}
          \epsilon_1 \\ \epsilon_2 \\ \vdots \\ \epsilon_N
        \end{pmatrix}
\\
&&
   &= \begin{pmatrix}
          \beta_0 + x_{11} \beta_1 + x_{12} \beta_2 + \cdots + x_{1K} \beta_K \\
          \beta_0 + x_{21} \beta_1 + x_{22} \beta_2 + \cdots + x_{2K} \beta_K \\
          \vdots \\
          \beta_0 + x_{N1} \beta_1 + x_{N2} \beta_2 + \cdots + x_{NK} \beta_K
      \end{pmatrix}
      + \begin{pmatrix}
          \epsilon_1 \\ \epsilon_2 \\ \vdots \\ \epsilon_N
        \end{pmatrix}
\\
&&
   &= \begin{pmatrix}
          \beta_0 + \sum_{k=1}^K \beta_k x_{1k} + \epsilon_1 \\
          \beta_0 + \sum_{k=1}^K \beta_k x_{2k} + \epsilon_2 \\
          \vdots \\
          \beta_0 + \sum_{k=1}^K \beta_k x_{Nk} + \epsilon_N
      \end{pmatrix}
\end{align}

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


  • 1: 説明変数  \mathbf{X} とパラメータ  \boldsymbol{\beta} の積と誤差項  \boldsymbol{\epsilon} の和の式を立てます。
  • 2: 行列とベクトル  \mathbf{X}, \boldsymbol{\beta} の要素を明示します。
  • 3: 行列とベクトルの積(行とベクトルの内積)を計算します。
  • 4: ベクトルの和(要素ごとの和)を計算します。

 説明変数の定義  x_{i0} = 1 より、 0 番(列)目の項は  \beta_0 x_{i0} = \beta_0 となり、 \beta_0 は(説明変数の係数ではなく)定数項になります。

 ここまでは、線形回帰モデルを確認しました。次からは、観測データとして  \mathbf{X}, \mathbf{y} が得られたときの線形回帰モデルのパラメータを推定します。

係数パラメータの計算式

 線形回帰モデルに対するOLSによる係数パラメータの推定値の計算式を導出します。

 説明変数の  k 番目の要素  x_{ik} に対する係数パラメータの推定値(最小二乗推定量)を  \hat{\beta}_k として、 K+1 個の係数パラメータの推定値をまとめて、 K+1 次元ベクトル  \hat{\boldsymbol{\beta}} とします。

 \displaystyle
\hat{\boldsymbol{\beta}}
    = \begin{pmatrix}
          \hat{\beta}_0 \\ \hat{\beta}_1 \\ \hat{\beta}_2 \\ \vdots \\ \hat{\beta}_K
      \end{pmatrix}

  i 番目の被説明変数の推定値を  \hat{y}_i として、 N 個の被説明変数の推定値をまとめて、 N 次元ベクトル  \hat{\mathbf{y}} とします。

 \displaystyle
\hat{\mathbf{y}}
    = \begin{pmatrix}
          \hat{y}_1 \\ \hat{y}_2 \\ \vdots \\ \hat{y}_N
      \end{pmatrix}


 各被説明変数  y_i の推定値  \hat{y}_i は、推定パラメータ  \hat{\boldsymbol{\beta}} を用いて求めた値です。

 \displaystyle
\begin{aligned}
\hat{y}_i
   &= \mathbf{x}_i^{\top} \hat{\boldsymbol{\beta}}
\\
   &= \begin{pmatrix}
          x_{i0} & x_{i1} & x_{i2} & \cdots & x_{iK}
      \end{pmatrix}
      \begin{pmatrix}
          \hat{\beta}_0 \\ \hat{\beta}_1 \\ \hat{\beta}_2 \\ \vdots \\ \hat{\beta}_K
      \end{pmatrix}
\\
   &= x_{i0} \hat{\beta}_0
      + x_{i1} \hat{\beta}_1
      + x_{i2} \hat{\beta}_2
      + \cdots
      + x_{iK} \hat{\beta}_K
\\
   &= \sum_{k=0}^K
          \hat{\beta}_k x_{ik}
\\
   &= \hat{\beta}_0
      + \sum_{k=1}^K
          \hat{\beta}_k x_{ik}
\end{aligned}

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


  • 1: 説明変数  \mathbf{x}_i とパラメータ  \hat{\boldsymbol{\beta}} の内積の式を立てます。
  • 2: ベクトル  \mathbf{x}_i, \hat{\boldsymbol{\beta}} の要素を明示します。
  • 3: ベクトルの内積(要素ごとの積の和)を計算します。
  • 4:  K+1 個の項の和を  \sum_{k=0}^K でまとめます。
  • 5: (項の形が異なる)  0 番目の項  \hat{\beta}_0 x_{i0} = \hat{\beta}_0 を取り出し(  K 個の要素の和を  \sum_{k=1}^K でまとめ)ます。

  N 個のデータをまとめると、次の式で計算できます。

 \displaystyle
\begin{align}
&&
\hat{\mathbf{y}}
   &= \mathbf{X} \hat{\boldsymbol{\beta}}
\tag{2}\\
\Leftrightarrow &&
\begin{pmatrix}
    \hat{y}_1 \\
    \hat{y}_2 \\
    \vdots \\
    \hat{y}_N
\end{pmatrix}
   &= \begin{pmatrix}
          1      & x_{11} & x_{12} & \cdots & x_{1K} \\
          1      & x_{21} & x_{22} & \cdots & x_{2K} \\
          \vdots & \vdots & \vdots & \ddots & \vdots \\
          1      & x_{N1} & x_{N2} & \cdots & x_{NK}
      \end{pmatrix}
      \begin{pmatrix}
          \hat{\beta}_0 \\
          \hat{\beta}_1 \\
          \hat{\beta}_2 \\
          \vdots \\
          \hat{\beta}_K
      \end{pmatrix}
\\
&&
   &= \begin{pmatrix}
          \hat{\beta}_0 + \sum_{k=1}^K \hat{\beta}_k x_{1k} \\
          \hat{\beta}_0 + \sum_{k=1}^K \hat{\beta}_k x_{2k} \\
          \vdots \\
          \hat{\beta}_0 + \sum_{k=1}^K \hat{\beta}_k x_{Nk}
      \end{pmatrix}
\end{align}

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


  • 1: 説明変数  \mathbf{X} とパラメータ  \hat{\boldsymbol{\beta}} の積の式を立てます。
  • 2: 行列とベクトル  \mathbf{X}, \hat{\boldsymbol{\beta}} の要素を明示します。
  • 3: 行列とベクトルの積(行とベクトルの内積)を計算します。

 観測データとして得られた  \mathbf{y} を実現値、推定パラメータから求めた  \hat{\mathbf{y}} を理論値とも呼びます。

 各データに関して、実現値  y_i と理論値  \hat{y}_i の差を残差と言い、 e_i とします。

 \displaystyle
e_i = y_i - \hat{y}_i

  N 個のデータの残差をまとめて、 N 次元ベクトル  \mathbf{e} とします。

 \displaystyle
\begin{aligned}
&&
\mathbf{e}
   &= \mathbf{y} - \hat{\mathbf{y}}
\\
\Leftrightarrow &&
\begin{pmatrix}
    e_1 \\ e_2 \\ \vdots \\ e_N
\end{pmatrix}
   &= \begin{pmatrix}
          y_1 \\ y_2 \\ \vdots \\ y_N
      \end{pmatrix}
      - \begin{pmatrix}
          \hat{y}_1 \\ \hat{y}_2 \\ \vdots \\ \hat{y}_N
        \end{pmatrix}
\\
&&
   &= \begin{pmatrix}
          y_1 - \hat{y}_1 \\ y_2 - \hat{y}_2 \\ \vdots \\ y_N - \hat{y}_N
      \end{pmatrix}
\end{aligned}

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


  • 1: 実現値  \mathbf{y} と理論値  \hat{\mathbf{y}} の差の式を立てます。
  • 2: ベクトル  \mathbf{y}, \hat{\mathbf{y}} の要素を明示します。
  • 3: ベクトルの差(要素ごとの差)を計算します。

 残差  \mathbf{e} の定義式に関して、理論値  \hat{\mathbf{y}} の式(2)で置き換えて使います。

 \displaystyle
\begin{align}
\mathbf{e}
   &= \mathbf{y} - \hat{\mathbf{y}}
\\
   &= \mathbf{y}
      - \mathbf{X} \hat{\boldsymbol{\beta}}
\tag{3}
\end{align}


 残差  \mathbf{e} の2乗和(残差平方和・ \mathbf{e} の内積)を考えます。

 \displaystyle
\begin{aligned}
\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}
      - \mathbf{y}^{\top}
        \mathbf{X} \hat{\boldsymbol{\beta}}
      - (\mathbf{X} \hat{\boldsymbol{\beta}})^{\top}
        \mathbf{y}
      + (\mathbf{X} \hat{\boldsymbol{\beta}})^{\top}
        \mathbf{X} \hat{\boldsymbol{\beta}}
\\
   &= \mathbf{y}^{\top} \mathbf{y}
      - \hat{\boldsymbol{\beta}}^{\top} \mathbf{X}^{\top}
        \mathbf{y}
      - \hat{\boldsymbol{\beta}}^{\top} \mathbf{X}^{\top}
        \mathbf{y}
      + \hat{\boldsymbol{\beta}}^{\top} \mathbf{X}^{\top}
        \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}}
\end{aligned}

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


  • 1:  \mathbf{e} の式(3)より、 \mathbf{e} の内積の式を立てます。
  • 2: 括弧を展開します。
  • 3: スカラは転置しても変化しないので、2つ目の項は  \mathbf{y}^{\top} \mathbf{X} \hat{\boldsymbol{\beta}} = (\mathbf{y}^{\top} \mathbf{X} \hat{\boldsymbol{\beta}})^{\top} で置き換えられ、転置行列の性質  (\mathbf{A} \mathbf{B} \mathbf{C})^{\top} = \mathbf{C}^{\top} \mathbf{B}^{\top} \mathbf{A}^{\top} より、 (\mathbf{y}^{\top} \mathbf{X} \hat{\boldsymbol{\beta}})^{\top} = \hat{\boldsymbol{\beta}}^{\top} \mathbf{X}^{\top} \mathbf{y} となります。
  • 3: 転置行列の性質  (\mathbf{A} \mathbf{B})^{\top} = \mathbf{B}^{\top} \mathbf{A}^{\top} より、3つ目の項は  (\mathbf{X} \hat{\boldsymbol{\beta}})^{\top} = \hat{\boldsymbol{\beta}}^{\top} \mathbf{X}^{\top} となります。
  • 4: 式を整理します。

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

 \displaystyle
\begin{aligned}
\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\}
      + \frac{\partial}{\partial \hat{\boldsymbol{\beta}}} \Bigl\{
          - 2 \hat{\boldsymbol{\beta}}^{\top} \mathbf{X}^{\top}
            \mathbf{y}
      \Bigr\}
      + \frac{\partial}{\partial \hat{\boldsymbol{\beta}}} \Bigl\{
          \hat{\boldsymbol{\beta}}^{\top}
          \mathbf{X}^{\top} \mathbf{X}
          \hat{\boldsymbol{\beta}}
      \Bigr\}
\\
   &= \mathbf{0}
      - 2
        \frac{\partial \hat{\boldsymbol{\beta}}^{\top}}{\partial \hat{\boldsymbol{\beta}}}
        \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}}
\end{aligned}

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


  • 1:  \mathbf{e}^{\top} \mathbf{e} に対する  \hat{\boldsymbol{\beta}} の偏微分の式を立てます。 \hat{\boldsymbol{\beta}} に関する微分なので、 \mathbf{X}, \mathbf{y} は定数として扱います。
  • 2: 和の微分  {f(x) + g(x)}' = f'(x) + g'(x) より、項ごとの微分の和に分割します。
  • 3: 定数の微分  {a}' = 0 より、1つ目の項は0ベクトルになり消えます。
  • 3: 定数倍の微分  {a f(x)}' = a f'(x) より、 \hat{\boldsymbol{\beta}} の係数を  \frac{\partial}{\partial \hat{\boldsymbol{\beta}}} の外に出します。
  • 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: 変数の微分  {x}' = 1 より、2つ目の項の微分が単位行列になり消えます。

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

 \displaystyle
\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}
\tag{4}
\end{align}

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


  • 1:  \frac{\partial \mathbf{e}^{\top} \mathbf{e}}{\partial \hat{\boldsymbol{\beta}}} = \mathbf{0} とおきます。
  • 2-3:  \hat{\boldsymbol{\beta}} について式を整理します。
  • 3: 両辺に左から  \mathbf{X}^{\top} \mathbf{X} の逆行列を掛けます。

  \hat{\boldsymbol{\beta}} の計算式が得られました。

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

分散パラメータの計算式

 線形回帰モデルに対するOLSによる誤差項の生成分布の分散パラメータの推定値の計算式を導出します。

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

 \displaystyle
\begin{align}
\mathbf{e}
   &= \mathbf{y}
      - \mathbf{X} \hat{\boldsymbol{\beta}}
\tag{3}\\
   &= \mathbf{y}
      - \mathbf{X} \Bigl(
            (\mathbf{X}^{\top} \mathbf{X})^{-1}
            \mathbf{X}^{\top}
            \mathbf{y}
        \Bigr)
\\
   &= (\mathbf{X} \boldsymbol{\beta} + \boldsymbol{\epsilon})
      - \mathbf{X} \Bigl(
            (\mathbf{X}^{\top} \mathbf{X})^{-1}
            \mathbf{X}^{\top} (
                \mathbf{X} \boldsymbol{\beta}
                + \boldsymbol{\epsilon}
            )
        \Bigr)
\\
   &= \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} \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{5}
\end{align}

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


  • 2:  \hat{\boldsymbol{\beta}} を式(4)で置き換えます。
  • 3:  \mathbf{y} を式(1)で置き換えます。
  • 4: 括弧を展開します。
  • 5: 逆行列の性質  \mathbf{A}^{-1} \mathbf{A} = \mathbf{I} より、 (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^{\top} \mathbf{X} が単位行列になり消えます。
  • 6:  \boldsymbol{\epsilon} を括り出して、式を整理します。
  • 7: 冪等行列の定義  \mathbf{M} = \mathbf{I} - \mathbf{X} (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^{\top} より、置き換えます。

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

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

 \displaystyle
\begin{align}
\mathbf{e}
   &= \mathbf{y}
      - \mathbf{X} \hat{\boldsymbol{\beta}}
\tag{3}\\
   &= \mathbf{y}
      - \mathbf{X}
        (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^{\top}
        \mathbf{y}
\\
   &= \mathbf{y}
      - \mathbf{S} \mathbf{y}
\\
   &= (\mathbf{I} - \mathbf{S})
      \mathbf{y}
\end{align}

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


  • 2:  \hat{\boldsymbol{\beta}} を式(4)で置き換えます。
  • 3: ハット行列の定義  \mathbf{S} = \mathbf{X} (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^{\top} より、置き換えます。
  • 4:  \mathbf{y} を括り出します。

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

 \displaystyle
\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{5}
\end{align}

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


  • 1: ハット行列を用いた冪等行列の定義  \mathbf{M} = \mathbf{I} - \mathbf{S} より、置き換えます。
  • 2:  \mathbf{y} を式(1)で置き換えます。
  • 3: 括弧を展開します。
  • 4-5: 冪等行列の性質  \mathbf{M} \mathbf{X} = \mathbf{0} より、0ベクトルになり消えます。

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

 残差(観測値と推定値の差)  \mathbf{e} と誤差(観測値と真値の差)  \boldsymbol{\epsilon} の関係式が得られました。

 残差  \mathbf{e} の式(5)を用いて、残差平方和を考えます。

 \displaystyle
\begin{aligned}
\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}
\end{aligned}

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


  • 1:  \mathbf{e} の式(5)より、 \mathbf{e} の内積の式を立てます。
  • 2: 転置行列の性質  (\mathbf{A} \mathbf{B})^{\top} = \mathbf{B}^{\top} \mathbf{A}^{\top} より、 (\mathbf{M} \boldsymbol{\epsilon})^{\top} = \boldsymbol{\epsilon}^{\top} \mathbf{M}^{\top} となります。
  • 3: 冪等行列の性質  \mathbf{M}^{\top} = \mathbf{M} \mathbf{M} \mathbf{M} = \mathbf{M} より、式を整理します。

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

 \displaystyle
\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{6}
\end{align}

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


  • 2: 横ベクトルと行列  \boldsymbol{\epsilon}^{\top}, \mathbf{M} の積を計算します。
  • 3: ベクトル  \boldsymbol{\epsilon}^{\top} \mathbf{M}, \boldsymbol{\epsilon} の内積を計算します。

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

 \displaystyle
\begin{align}
\mathbb{E}[\mathbf{e}^{\top} \mathbf{e}]
   &= \mathbb{E} \left[
          \sum_{i=1}^N \sum_{j=1}^N
              \epsilon_i m_{ij} \epsilon_j
      \right]
\\
   &= \sum_{i=1}^N \sum_{j=1}^N
          m_{ij}
          \mathbb{E}[\epsilon_i \epsilon_j]
\tag{7}
\end{align}

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


  • 1:  \mathbf{e}^{\top} \mathbf{e} の期待値の式を立てます(  \mathbf{e}^{\top} \mathbf{e} の式(6)の両辺で期待値をとります)。
  • 2:  \mathbf{X} によって定義される  \mathbf{M} \boldsymbol{\epsilon} と無関係なので、期待値の性質  \mathbb{E}[a x] = a \mathbb{E}[x] より、 \mathbf{M} の各要素  m_{ij} \mathbb{E}[\cdot] の外に出せます。

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

 \displaystyle
\mathbb{E}[\epsilon_i]
    = \mu = 0

 よって、 \epsilon_i, \epsilon_j の積の期待値は  \epsilon_i, \epsilon_j の共分散  \sigma_{ij} になります。

 \displaystyle
\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}
\tag{8}
\end{align}

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


  • 1: 左辺  \mathbb{E}[\epsilon_i \epsilon_j] に関して、 \epsilon_i = \epsilon_i - 0 に変形します。
  • 2:  \mathbb{E}[\epsilon_i] = 0 で置き換えます。
  • 3: 共分散の定義  \sigma_{xy} = \mathbb{E}[(x - \mathbb{E}[x]) (y - \mathbb{E}[y])] より、置き換えます。

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

 \displaystyle
\sigma_{ij}
    = \begin{cases}
          \sigma^2
             &\quad (i = j) \\
          0
             &\quad (i \neq j)
      \end{cases}
\tag{9}

 この関係を用いて、残差平方和の期待値  \mathbb{E}[\mathbf{e}^{\top} \mathbf{e}] の式(6)を整理します。

 \displaystyle
\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\}
\\
   &= \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:  \mathbb{E}[\mathbf{e}^{\top} \mathbf{e}] の式(7)の  j に関する項の和  \sum_{j=1}^N を展開します。
  • 2: 各組み合わせの積の期待値  \mathbb{E}[\epsilon_i \epsilon_j] を式(8)で置き換えます。
  • 3: 各組み合わせの共分散  \sigma_{ij} を式(9)で置き換えます。
  • 4: 同じデータインデックスの組み合わせ  i = j の項のみ残ります。

 さらに変形します。

 \displaystyle
\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)
\end{align}

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


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

 この式を分散パラメータ  \sigma^2 について解きます。

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

  \sigma^2 の計算式が得られました。

 残差平方和の期待値  \mathbb{E}[\mathbf{e}^{\top} \mathbf{e}] を残差平方和  \mathbf{e}^{\top} \mathbf{e} で近似

 \displaystyle
\mathbb{E}[\mathbf{e}^{\top} \mathbf{e}]
    \approx
      \mathbf{e}^{\top} \mathbf{e}

して、置き換えます(なぜこんな操作をしてもいいのか?)。

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


 以上で、線形回帰モデルの誤差項の分散パラメータの計算式が求まりました。

係数パラメータの分散共分散行列の計算式

 線形回帰モデルに対するOLSによる係数パラメータの推定値の分散共分散行列の計算式を導出します。

 誤差項  \boldsymbol{\epsilon} の期待値は、0ベクトルになります。

 \displaystyle
\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{10}
\end{align}

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


  • 1:  \boldsymbol{\epsilon} の期待値の式を立てて、ベクトル  \boldsymbol{\epsilon} の要素を明示します。
  • 2: ベクトルの期待値を要素ごとの期待値に変形します。
  • 3: 誤差項の仮定より、各誤差項の期待値は  \mathbb{E}[\epsilon_i] = 0 です。

 誤差項  \boldsymbol{\epsilon} の分散共分散行列は、 \boldsymbol{\epsilon} の2乗和(内積)の期待値になります。

 \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-2:  \mathbb{E}[\boldsymbol{\epsilon}] を式(10)で置き換えると、0ベクトルになり消えます。

 さらに、要素を目指して計算すると、分散パラメータ  \sigma^2 と単位行列  \mathbf{I} の積になります。

 \displaystyle
\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_{11} & \sigma_{12} & \cdots & \sigma_{1N} \\
          \sigma_{21} & \sigma_{22} & \cdots & \sigma_{2N} \\
          \vdots & \vdots & \ddots & \vdots \\
          \sigma_{N1} & \sigma_{N2} & \cdots & \sigma_{NN}
      \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{11}
\end{align}

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


  • 1: ベクトル  \boldsymbol{\epsilon} の要素を明示します。
  • 2: ベクトル  \boldsymbol{\epsilon}, \boldsymbol{\epsilon}^{\top} の積を計算します。
  • 3: 行列の期待値を要素ごとの期待値に変形します。
  • 4: 各組み合わせの積の期待値  \mathbb{E}[\epsilon_i \epsilon_j] を式(8)で置き換えます。
  • 5: 各組み合わせの共分散  \sigma_{ij} を式(9)で置き換えます。
  • 6-7: スカラと行列の積に分解します。


 係数パラメータの推定値  \hat{\boldsymbol{\beta}} の式(4)を変形します。

 \displaystyle
\begin{align}
\hat{\boldsymbol{\beta}}
   &= (\mathbf{X}^{\top} \mathbf{X})^{-1}
      \mathbf{X}^{\top} \mathbf{y}
\tag{4}\\
   &= (\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{12}
\end{align}

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


  • 2:  \mathbf{y} を式(1)で置き換えます。
  • 4: 括弧を展開します。
  • 5: 逆行列の性質  \mathbf{A}^{-1} \mathbf{A} = \mathbf{I} より、 (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^{\top} \mathbf{X} が単位行列になり消えます。

 推定パラメータ  \hat{\boldsymbol{\beta}} の期待値は、真のパラメータ  \boldsymbol{\beta} になります。

 \displaystyle
\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{13}
\end{align}

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


  • 1:  \hat{\boldsymbol{\beta}} の期待値の式を立てます(  \hat{\boldsymbol{\beta}} の式(12)の両辺で期待値をとります)。
  • 2: 和の期待値  \mathbb{E}[x + y] = \mathbb{E}[x] + \mathbb{E}[y] より、項ごとの期待値に分割します。
  • 3: 観測値  \mathbf{X} は定数として扱うので、 \mathbb{E}[\cdot] の外に出せます。
  • 4-5:  \mathbb{E}[\boldsymbol{\epsilon}] を式(10)で置き換えると、後の項が0ベクトルになり消えます。

 よって、推定パラメータ  \hat{\boldsymbol{\beta}} と推定パラメータの期待値  \mathbb{E}[\hat{\boldsymbol{\beta}}] は、次の式になります。

 \displaystyle
\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{14}
\end{align}

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


  • 1:  \hat{\boldsymbol{\beta}}, \mathbb{E}[\hat{\boldsymbol{\beta}}] の差の式を立てて、 \mathbb{E}[\hat{\boldsymbol{\beta}}] を式(13)で置き換えます。
  • 2:  \hat{\boldsymbol{\beta}} を式(12)で置き換えます。
  • 3: 式を整理します。

 推定パラメータ  \hat{\boldsymbol{\beta}} の分散共分散行列を考えます。

 \displaystyle
\begin{aligned}
\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}
\end{aligned}

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


  • 1:  \hat{\boldsymbol{\beta}} の分散共分散行列の式を立てます。
  • 1-2:  \hat{\boldsymbol{\beta}}, \mathbb{E}[\hat{\boldsymbol{\beta}}] の差を式(14)で置き換えます。
  • 3: 転置行列の性質  (\mathbf{A} \mathbf{B} \mathbf{C})^{\top} = \mathbf{C}^{\top} \mathbf{B}^{\top} \mathbf{A}^{\top}、転置行列と逆行列の関係  (\mathbf{A}^{-1})^{\top} = (\mathbf{A}^{\top})^{-1} より、 ((\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^{\top}\boldsymbol{\epsilon})^{\top} = \boldsymbol{\epsilon}^{\top} \mathbf{X} (\mathbf{X}^{\top} \mathbf{X})^{-1} となります。
  • 4: 観測値  \mathbf{X} は定数として扱うので、 \mathbb{E}[\cdot] の外に出せます。
  • 5:  \mathbb{E}[\boldsymbol{\epsilon} \boldsymbol{\epsilon}^{\top}] を式(11)で置き換えます。
  • 6: スカラ  \sigma^2 を前に出します。
  • 7: 逆行列の性質より、 (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^{\top} \mathbf{X} または  \mathbf{X}^{\top} \mathbf{X} (\mathbf{X}^{\top} \mathbf{X})^{-1} が単位行列になり消えます。

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

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

参考文献

おわりに

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

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

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

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

【次の内容】

www.anarchive-beta.com