1. 最尤推定の基本概念
最尤推定（Maximum Likelihood Estimation, MLE）は、観測データが与えられた時に、そのデータが得られる確率を最大にするようなパラメータを推定する手法です。
対数尤度関数を $L(\theta)$ とすると、最尤推定は以下の最適化問題として定式化されます：
\text{maximize} \quad L(\theta), \quad \theta \in \mathbb{R}^n
2. ニュートン法による最尤推定
2.1 基本的なアプローチ
ニュートン法による更新式：
\Delta \theta = -[\nabla^2 L(\theta)]^{-1} \nabla L(\theta)
ここで：
- $\nabla L(\theta)$ は対数尤度関数の勾配（スコア関数）
- $\nabla^2 L(\theta)$ は対数尤度関数のヘッセ行列
2.2 フィッシャー情報量との関係
対数尤度関数の負のヘッセ行列の期待値は、フィッシャー情報量行列 $\mathcal{I}(\theta)$ として知られています：
\mathcal{I}(\theta) = -\mathbb{E}[\nabla^2 L(\theta)]
これを用いると、ニュートン法の更新式は以下のように書き換えられます：
\Delta \theta \approx \mathcal{I}(\theta)^{-1} \nabla L(\theta)
3. 具体例：多変量正規分布の最尤推定
3.1 パラメータ化と対数尤度関数
多変量正規分布のパラメータを $\theta = (\boldsymbol{\mu}, \text{vec}(\Sigma))$ とベクトル化して表現できます。
ここで $\text{vec}(\Sigma)$ は共分散行列を列ベクトルに変換する演算子です。
サンプル集合 $\mathbf{X} = {\mathbf{x}_1, \mathbf{x}_2, \dots, \mathbf{x}_m}$ に対する対数尤度関数：
L(\theta) = -\frac{mn}{2}\log(2\pi) - \frac{m}{2}\log |\Sigma| - \frac{1}{2}\sum_{i=1}^m (\mathbf{x}_i - \boldsymbol{\mu})^T \Sigma^{-1} (\mathbf{x}_i - \boldsymbol{\mu})
3.2 勾配とヘッセ行列
パラメータ $\theta$ に関する勾配ベクトルは以下のように表されます：
\nabla L(\theta) = \begin{bmatrix}
\nabla_{\boldsymbol{\mu}} L \\
\nabla_{\text{vec}(\Sigma)} L
\end{bmatrix}
ここで：
\nabla_{\boldsymbol{\mu}} L = \sum_{i=1}^m \Sigma^{-1}(\mathbf{x}_i - \boldsymbol{\mu})
\nabla_{\text{vec}(\Sigma)} L = \text{vec}\left(-\frac{m}{2} \Sigma^{-1} + \frac{1}{2}\sum_{i=1}^m \Sigma^{-1}(\mathbf{x}_i - \boldsymbol{\mu})(\mathbf{x}_i - \boldsymbol{\mu})^T\Sigma^{-1}\right)
ヘッセ行列は以下のブロック行列として表されます：
\nabla^2 L(\theta) = \begin{bmatrix}
\nabla^2_{\boldsymbol{\mu}\boldsymbol{\mu}} L & \nabla^2_{\boldsymbol{\mu}\Sigma} L \\
\nabla^2_{\Sigma\boldsymbol{\mu}} L & \nabla^2_{\Sigma\Sigma} L
\end{bmatrix}
3.3 ニュートン法による直接的な導出
ニュートン法の更新式 $\Delta \theta = -[\nabla^2 L(\theta)]^{-1} \nabla L(\theta)$ を用いて、最適解を直接導出することができます。
Step 1: ブロック行列による表現
更新式をブロック行列で表現すると：
\begin{bmatrix}
\Delta \boldsymbol{\mu} \\
\Delta \text{vec}(\Sigma)
\end{bmatrix} = -
\begin{bmatrix}
\nabla^2_{\boldsymbol{\mu}\boldsymbol{\mu}} L & \nabla^2_{\boldsymbol{\mu}\Sigma} L \\
\nabla^2_{\Sigma\boldsymbol{\mu}} L & \nabla^2_{\Sigma\Sigma} L
\end{bmatrix}^{-1}
\begin{bmatrix}
\nabla_{\boldsymbol{\mu}} L \\
\nabla_{\text{vec}(\Sigma)} L
\end{bmatrix}
具体的な値を代入すると：
\begin{bmatrix}
\Delta \boldsymbol{\mu} \\
\Delta \text{vec}(\Sigma)
\end{bmatrix} = -
\begin{bmatrix}
-m\Sigma^{-1} & * \\
* & *
\end{bmatrix}^{-1}
\begin{bmatrix}
\sum_{i=1}^m \Sigma^{-1}(\mathbf{x}_i - \boldsymbol{\mu}) \\
\text{vec}\left(-\frac{m}{2} \Sigma^{-1} + \frac{1}{2}\sum_{i=1}^m \Sigma^{-1}(\mathbf{x}_i - \boldsymbol{\mu})(\mathbf{x}_i - \boldsymbol{\mu})^T\Sigma^{-1}\right)
\end{bmatrix}
ここで、$*$ は計算に直接影響しない項を表します。この式を解くと：
\begin{bmatrix}
\Delta \boldsymbol{\mu} \\
\Delta \text{vec}(\Sigma)
\end{bmatrix} =
\begin{bmatrix}
\frac{1}{m}\sum_{i=1}^m \mathbf{x}_i - \boldsymbol{\mu} \\
\text{vec}\left(\frac{1}{m}\sum_{i=1}^m (\mathbf{x}_i - \hat{\boldsymbol{\mu}})(\mathbf{x}_i - \hat{\boldsymbol{\mu}})^T - \Sigma\right)
\end{bmatrix}
したがって、最適解は一回の更新で以下のように得られます：
\begin{bmatrix}
\hat{\boldsymbol{\mu}} \\
\text{vec}(\hat{\Sigma})
\end{bmatrix} =
\begin{bmatrix}
\frac{1}{m}\sum_{i=1}^m \mathbf{x}_i \\
\text{vec}\left(\frac{1}{m}\sum_{i=1}^m (\mathbf{x}_i - \hat{\boldsymbol{\mu}})(\mathbf{x}_i - \hat{\boldsymbol{\mu}})^T\right)
\end{bmatrix}
この結果は、多変量正規分布の最尤推定量としてよく知られた以下の式と一致します：
\hat{\boldsymbol{\mu}} = \frac{1}{m}\sum_{i=1}^m \mathbf{x}_i
\hat{\Sigma} = \frac{1}{m}\sum_{i=1}^m (\mathbf{x}_i - \hat{\boldsymbol{\mu}})(\mathbf{x}_i - \hat{\boldsymbol{\mu}})^T