はじめに

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

　この記事では、SEMの最尤推定について、数式を使って解説します。

【前の内容】

www.anarchive-beta.com

【他の内容】

www.anarchive-beta.com

【今回の内容】

はじめに
0.2.2 空間誤差モデル(SEM)の最尤法の導出
- モデルの確認
- パラメータの推定
参考文献
おわりに

0.2.2 空間誤差モデル(SEM)の最尤法の導出

　空間誤差モデル(SEM・Spatial Error Model)に対する最尤推定(MLE・Maximum Likelihood Estimation)を導出します。
　SEMの定義式や対数尤度関数については「【Python】0.1.3：空間誤差モデル(SEM)の定義式【はじめての地理空間DSのノート】 - からっぽのしょこ」を参照してください。

モデルの確認

　まずは、SEMの定義(仮定)を数式で確認します。
　空間重み行列については「空間重み行列の定義式」、多変量正規分布については「多次元ガウス分布の定義式 - からっぽのしょこ」を参照してください。

定義式

　SEMの定義式については「SEMの定義式」を参照してください。

　SEMは、次の式で定義されます。

$\begin{align} \mathbf{y} &= \mathbf{X} \boldsymbol{\beta} + \mathbf{u} \tag{0.36.a}\\ \mathbf{u} &= \lambda \mathbf{W} \mathbf{u} + \boldsymbol{\epsilon} \tag{0.36.b}\\ y_n &= \mathbf{x}_n^{\top} \boldsymbol{\beta} + u_n \\ &= \beta_0 + \sum_{k=1}^K x_{nk} \beta_k + u_n \\ u_n &= \lambda \mathbf{w}_n^{\top} \boldsymbol{u} + \epsilon_n \\ &= \lambda \sum_{j=1}^N w_{nj} u_j + \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$ 次元ベクトル)、 $u_n$ は地域 $n$ の空間誤差項(スカラ)、 $\mathbf{u}$ は $N$ 個の地域の空間誤差項( $N$ 次元ベクトル)、 $\mathbf{w}_n$ は地域 $n$ に関する重み( $N$ 次元ベクトル)、 $\mathbf{W}$ は空間重み行列( $N \times N$ の行列)、 $\lambda$ は空間パラメータ(スカラ)です。

　独立誤差項 $\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{0.12.b}\\ \epsilon_n &\sim \mathcal{N}(0, \sigma^2) \end{align}$

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

計算式

　SEMに関する計算式については「SEMの定義式」を参照してください。

　空間自己回帰に関する項について、次のようにおきます。

$\displaystyle \mathbf{A} = \mathbf{I} - \lambda \mathbf{W} \tag{1}$

　定義式(0.36)を独立誤差項 $\boldsymbol{\epsilon}$ について整理すると、次の式となります。

$\begin{align} \boldsymbol{\epsilon} &= (\mathbf{I} - \lambda \mathbf{W}) \mathbf{y} - (\mathbf{I} - \lambda \mathbf{W}) \mathbf{X} \boldsymbol{\beta} \\ &= \mathbf{A} \mathbf{y} - \mathbf{A} \mathbf{X} \boldsymbol{\beta} \\ &= \mathbf{A} (\mathbf{y} - \mathbf{X} \boldsymbol{\beta}) \tag{2} \end{align}$

尤度関数

　SEMの尤度関数については「SEMの尤度関数の導出」を参照してください。

　SEMのパラメータ $\boldsymbol{\beta}, \sigma^2, \lambda$ をまとめて、パラメータベクトル $\boldsymbol{\theta} = (\boldsymbol{\beta}, \sigma^2, \lambda)^{\top}$ とします。

　尤度関数 $L(\boldsymbol{\theta})$ は、次の式となります。

$\begin{align} L(\boldsymbol{\theta}) &= p(\mathbf{y} \mid \mathbf{X}, \boldsymbol{\theta}) \\ &= (2 \pi)^{-\frac{N}{2}} (\sigma^2)^{-\frac{N}{2}} \exp \Bigl( - \frac{1}{2} \frac{1}{\sigma^2} \boldsymbol{\epsilon}^{\top} \boldsymbol{\epsilon} \Bigr) |\mathbf{A}| \tag{3} \end{align}$

　ここで、 $|\mathbf{X}|$ は行列式です。

　対数尤度関数は、次の式となります。

$\begin{align} \log L(\boldsymbol{\theta}) &= \log p(\mathbf{y} \mid \mathbf{X}, \boldsymbol{\theta}) \\ &= - \frac{N \log (2 \pi)}{2} - \frac{N}{2} \log \sigma^2 - \frac{1}{2} \frac{1}{\sigma^2} \boldsymbol{\epsilon}^{\top} \boldsymbol{\epsilon} + \log |\mathbf{A}| \tag{0.37} \end{align}$

対数尤度関数の微分

　SEMの対数尤度関数の微分については「SEMのヘッセ行列の導出」を参照してください。

　対数尤度関数の各パラメータ $\boldsymbol{\beta}, \sigma^2, \lambda$ に関する微分は、それぞれ次の式になります。

$\begin{align} \frac{\partial \log L(\boldsymbol{\theta})}{\partial \boldsymbol{\beta}} &= \frac{1}{\sigma^2} \mathbf{X}^{\top} \mathbf{A}^{\top} \boldsymbol{\epsilon} \tag{4}\\ \frac{\partial \log L(\boldsymbol{\theta})}{\partial \sigma^2} &= - \frac{N}{2} \frac{1}{\sigma^2} + \frac{1}{2} \frac{1}{(\sigma^2)^2} \boldsymbol{\epsilon}^{\top} \boldsymbol{\epsilon} \tag{5}\\ \frac{\partial \log L(\boldsymbol{\theta})}{\partial \lambda} &= - \mathrm{Tr}(\mathbf{A}^{-1} \mathbf{W}) + \frac{1}{\sigma^2} \boldsymbol{\epsilon}^{\top} \mathbf{W} (\mathbf{y} - \mathbf{X} \boldsymbol{\beta}) \tag{0.40} \end{align}$

　以上の式を用いてパラメータの計算式を求めます。

スポンサードリンク

パラメータの推定

　次は、最尤推定によるパラメータの計算式(最尤推定量・最尤解)を導出します。

回帰パラメータの最尤推定量

　回帰パラメータの最尤推定量を求めます。

　対数尤度関数の回帰パラメータ $\boldsymbol{\beta}$ に関する微分を変形します。

$\begin{align} \frac{\partial \log L(\boldsymbol{\theta})}{\partial \boldsymbol{\beta}} &= \frac{1}{\sigma^2} \mathbf{X}^{\top} \mathbf{A}^{\top} \boldsymbol{\epsilon} \tag{4}\\ &= \frac{1}{\sigma^2} (\mathbf{A} \mathbf{X})^{\top} \mathbf{A} (\mathbf{y} - \mathbf{X} \boldsymbol{\beta}) \\ &= - \frac{1}{\sigma^2} \Bigl( - (\mathbf{A} \mathbf{X})^{\top} \mathbf{A} \mathbf{y} + (\mathbf{A} \mathbf{X})^{\top} \mathbf{A} \mathbf{X} \boldsymbol{\beta} \Bigr) \tag{4'} \end{align}$

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

1: 式(4)を再掲しています。
2: 転置の性質 $(\mathbf{A} \mathbf{B})^{\top} = \mathbf{B}^{\top} \mathbf{A}^{\top}$ より、 $\mathbf{A}, \mathbf{X}$ をまとめます。
2: $\boldsymbol{\epsilon}$ に、式(2)を代入します。
3: $- (\mathbf{A} \mathbf{X})^{\top} \mathbf{A}$ を括弧の中に入れます。

　 $\frac{\partial \log L(\boldsymbol{\theta})}{\partial \boldsymbol{\beta}} = \mathbf{0}$ とおき、 $\boldsymbol{\beta}$ について解きます。

$\begin{align} && - \frac{1}{\sigma^2} \Bigl( - (\mathbf{A} \mathbf{X})^{\top} \mathbf{A} \mathbf{y} + (\mathbf{A} \mathbf{X})^{\top} \mathbf{A} \mathbf{X} \boldsymbol{\beta} \Bigr) &= \mathbf{0} \\ \Rightarrow && - (\mathbf{A} \mathbf{X})^{\top} \mathbf{A} \mathbf{y} + (\mathbf{A} \mathbf{X})^{\top} \mathbf{A} \mathbf{X} \boldsymbol{\beta} &= \mathbf{0} \\ \Rightarrow && (\mathbf{A} \mathbf{X})^{\top} \mathbf{A} \mathbf{X} \boldsymbol{\beta} &= (\mathbf{A} \mathbf{X})^{\top} \mathbf{A} \mathbf{y} \\ \Rightarrow && \boldsymbol{\beta} &= \Bigl( (\mathbf{A} \mathbf{X})^{\top} \mathbf{A} \mathbf{X} \Bigr)^{-1} (\mathbf{A} \mathbf{X})^{\top} \mathbf{A} \mathbf{y} \equiv \hat{\boldsymbol{\beta}}_{\mathrm{ML}}(\lambda) \tag{6}\\ && &= \Bigl( \mathbf{X}^{\top} \mathbf{A}^{\top} \mathbf{A} \mathbf{X} \Bigr)^{-1} \mathbf{X}^{\top} \mathbf{A}^{\top} \mathbf{A} \mathbf{y} \end{align}$

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

1: 式(4')を0ベクトルとおきます。
2: 両辺に $-\sigma^2$ を掛けます。
3: $\boldsymbol{\beta}$ 以外の項を右辺に移します。
4: 両辺に $(\mathbf{A} \mathbf{X})^{\top} \mathbf{A} \mathbf{X}$ の逆行列を左から掛けます。
5: 転置の性質 $(\mathbf{A} \mathbf{B})^{\top} = \mathbf{B}^{\top} \mathbf{A}^{\top}$ より、転置を戻します。(確認用)

　空間パラメータ $\lambda$ を固定すると、対数尤度関数 $\log L(\boldsymbol{\theta})$ を最大化する回帰パラメータ $\hat{\boldsymbol{\beta}}$ の値が定まることから、 $\boldsymbol{\beta}$ の最尤推定量を $\lambda$ の関数 $\hat{\boldsymbol{\beta}}_{\mathrm{ML}}(\lambda)$ として扱います。

　式(0.38)について、次のようにおきます。

$\begin{align} \mathbf{X}(\lambda) &= \mathbf{A} \mathbf{X} \tag{7}\\ &= (\mathbf{I} - \lambda \mathbf{W}) \mathbf{X} \\ \mathbf{y}(\lambda) &= \mathbf{A} \mathbf{y} \tag{8}\\ &= (\mathbf{I} - \lambda \mathbf{W}) \mathbf{y} \end{align}$

　回帰パラメータ $\boldsymbol{\beta}$ の最尤推定量について、式(7)(8)で置き換えます。

$\begin{align} \hat{\boldsymbol{\beta}}_{\mathrm{ML}}(\lambda) &= \Bigl( (\mathbf{A} \mathbf{X})^{\top} \mathbf{A} \mathbf{X} \Bigr)^{-1} (\mathbf{A} \mathbf{X})^{\top} \mathbf{A} \mathbf{y} \tag{6}\\ &= \Bigl( \mathbf{X}(\lambda)^{\top} \mathbf{X}(\lambda) \Bigr)^{-1} \mathbf{X}(\lambda)^{\top} \mathbf{y}(\lambda) \tag{0.38} \end{align}$

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

1: 式(6)を再掲しています。
2: $\mathbf{A} \mathbf{X}, \mathbf{A} \mathbf{y}$ を、式(7)(8)で置き換えます。

　以上で、回帰パラメータの最尤推定量の式が得られました。

分散パラメータの最尤推定量

　分散パラメータの最尤推定量を求めます。

　対数尤度関数の分散パラメータ $\sigma^2$ に関する微分を変形します。

$\begin{align} \frac{\partial \log L(\boldsymbol{\theta})}{\partial \sigma^2} &= - \frac{N}{2} \frac{1}{\sigma^2} + \frac{1}{2} \frac{1}{(\sigma^2)^2} \boldsymbol{\epsilon}^{\top} \boldsymbol{\epsilon} \tag{5}\\ &= - \frac{1}{2 \sigma^4} \Bigl( N \sigma^2 - \boldsymbol{\epsilon}^{\top} \boldsymbol{\epsilon} \Bigr) \tag{5'} \end{align}$

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

1: 式(5)を再掲しています。
2: $- \frac{1}{2 \sigma^4}$ を括り出します。

　 $\frac{\partial \log L(\boldsymbol{\theta})}{\partial \sigma^2} = 0$ とおき、 $\sigma^2$ について解きます。

$\begin{align} && - \frac{1}{2 \sigma^4} \Bigl( N \sigma^2 - \boldsymbol{\epsilon}^{\top} \boldsymbol{\epsilon} \Bigr) &= 0 \\ \Rightarrow && N \sigma^2 - \boldsymbol{\epsilon}^{\top} \boldsymbol{\epsilon} &= 0 \\ \Rightarrow && N \sigma^2 &= \boldsymbol{\epsilon}^{\top} \boldsymbol{\epsilon} \\ \Rightarrow && \sigma^2 &= \frac{1}{N} \boldsymbol{\epsilon}^{\top} \boldsymbol{\epsilon} \equiv \hat{\sigma}_{\mathrm{ML}}^2(\lambda) \tag{0.39}\\ && &= \frac{1}{N} (\mathbf{y} - \mathbf{X} \boldsymbol{\beta})^{\top} \mathbf{A}^{\top} \mathbf{A} (\mathbf{y} - \mathbf{X} \boldsymbol{\beta}) \end{align}$

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

1: 式(0.5')を0とおきます。
2: 両辺に $-2 \sigma^4$ を掛けます。
3: $\sigma^2$ 以外の項を右辺に移します。
4: 両辺に $\frac{1}{N}$ を掛けます。
5: $\boldsymbol{\epsilon}$ を、式(2)で戻します。(確認用)

　空間パラメータ $\lambda$ を固定すると、対数尤度関数 $\log L(\boldsymbol{\theta})$ を最大化する分散パラメータ $\hat{\sigma}^2$ の値が定まることから、 $\sigma^2$ の最尤推定量を $\lambda$ の関数 $\hat{\sigma}_{\mathrm{ML}}^2(\lambda)$ として扱います。

　独立誤差項 $\boldsymbol{\epsilon}$ について、 $\boldsymbol{\beta}$ の最尤推定量を用いて、 $\hat{\boldsymbol{\beta}}_{\mathrm{ML}}(\lambda), \lambda$ の関数 $\boldsymbol{\epsilon}(\hat{\boldsymbol{\beta}}_{\mathrm{ML}}, \lambda)$ とおきます。

$\begin{align} \boldsymbol{\epsilon} \Bigl( \hat{\boldsymbol{\beta}}_{\mathrm{ML}}(\lambda), \lambda \Bigr) &= \mathbf{A} \Bigl( \mathbf{y} - \mathbf{X} \hat{\boldsymbol{\beta}}_{\mathrm{ML}}(\lambda) \Bigr) \tag{2'}\\ &= \mathbf{A} \mathbf{y} - \mathbf{A} \mathbf{X} \hat{\boldsymbol{\beta}}_{\mathrm{ML}}(\lambda) \\ &= \mathbf{y}(\lambda) - \mathbf{X}(\lambda) \hat{\boldsymbol{\beta}}_{\mathrm{ML}}(\lambda) \tag{9}\\ &= \mathbf{y}(\lambda) - \mathbf{X}(\lambda) \Bigl( \mathbf{X}(\lambda)^{\top} \mathbf{X}(\lambda) \Bigr)^{-1} \mathbf{X}(\lambda)^{\top} \mathbf{y}(\lambda) \\ &= \mathbf{A} \mathbf{y} - \mathbf{A} \mathbf{X} \Bigl( \mathbf{X}^{\top} \mathbf{A}^{\top} \mathbf{A} \mathbf{X} \Bigr)^{-1} \mathbf{X}^{\top} \mathbf{A}^{\top} \mathbf{A} \mathbf{y} \end{align}$

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

1: 式(2)について、 $\boldsymbol{\beta}$ を式(0.38)で置き換えた式を立てます。
2: 括弧を展開します。
3: $\mathbf{A} \mathbf{X}, \mathbf{A} \mathbf{y}$ を、式(7)(8)で置き換えます。
4: $\hat{\boldsymbol{\beta}}_{\mathrm{ML}}(\lambda)$ に、式(0.38)を代入します。(確認用)
5: $\mathbf{X}(\lambda), \mathbf{y}(\lambda)$ を、式(7)(8)で戻します。(確認用)

　分散パラメータ $\sigma^2$ の最尤推定量について、 $\boldsymbol{\beta}$ の最尤推定量を用いて、 $\hat{\boldsymbol{\beta}}_{\mathrm{ML}}(\lambda), \lambda$ の関数 $\hat{\sigma}_{\mathrm{ML}}^2(\hat{\boldsymbol{\beta}}_{\mathrm{ML}}, \lambda)$ とおきます。

$\begin{align} \hat{\sigma}_{\mathrm{ML}}^2 \Bigl( \hat{\boldsymbol{\beta}}_{\mathrm{ML}}(\lambda), \lambda \Bigr) &= \frac{1}{N} \boldsymbol{\epsilon} ( \hat{\boldsymbol{\beta}}_{\mathrm{ML}}, \lambda )^{\top} \boldsymbol{\epsilon} ( \hat{\boldsymbol{\beta}}_{\mathrm{ML}}, \lambda ) \tag{0.39'}\\ &= \frac{1}{N} \Bigl( \mathbf{y}(\lambda) - \mathbf{X}(\lambda) \hat{\boldsymbol{\beta}}_{\mathrm{ML}}(\lambda) \Bigr)^{\top} \Bigl( \mathbf{y}(\lambda) - \mathbf{X}(\lambda) \hat{\boldsymbol{\beta}}_{\mathrm{ML}}(\lambda) \Bigr) \\ &= \frac{1}{N} \Biggl( \mathbf{y}(\lambda)^{\top} \mathbf{y}(\lambda) - \Bigl( \mathbf{X}(\lambda) \hat{\boldsymbol{\beta}}_{\mathrm{ML}}(\lambda) \Bigr)^{\top} \mathbf{y}(\lambda) \Biggr. \\ &\qquad \Biggl. - \mathbf{y}(\lambda)^{\top} \mathbf{X}(\lambda) \hat{\boldsymbol{\beta}}_{\mathrm{ML}}(\lambda) + \Bigl( \mathbf{X}(\lambda) \hat{\boldsymbol{\beta}}_{\mathrm{ML}}(\lambda) \Bigr)^{\top} \mathbf{X}(\lambda) \hat{\boldsymbol{\beta}}_{\mathrm{ML}}(\lambda) \Biggr) \\ &= \frac{1}{N} \Bigl( \mathbf{y}(\lambda)^{\top} \mathbf{y}(\lambda) - \mathbf{y}(\lambda)^{\top} \mathbf{X}(\lambda) \hat{\boldsymbol{\beta}}_{\mathrm{ML}}(\lambda) - \mathbf{y}(\lambda)^{\top} \mathbf{X}(\lambda) \hat{\boldsymbol{\beta}}_{\mathrm{ML}}(\lambda) + \mathbf{y}(\lambda)^{\top} \mathbf{X}(\lambda) \hat{\boldsymbol{\beta}}_{\mathrm{ML}}(\lambda) \Bigr) \\ &= \frac{1}{N} \mathbf{y}(\lambda)^{\top} \Bigl( \mathbf{y}(\lambda) - \mathbf{X}(\lambda) \hat{\boldsymbol{\beta}}_{\mathrm{ML}}(\lambda) \Bigr) \\ &= \frac{1}{N} \mathbf{y}(\lambda)^{\top} \boldsymbol{\epsilon} ( \hat{\boldsymbol{\beta}}_{\mathrm{ML}}, \lambda ) \tag{10} \end{align}$

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

1: 式(0.39)について、 $\boldsymbol{\epsilon}$ を式(9)で置き換えた式を立てます。
2: $\boldsymbol{\epsilon}(\hat{\boldsymbol{\beta}}_{\mathrm{ML}}, \lambda)$ に、式(9)を代入します。
3: 転置の性質 $(\mathbf{A} + \mathbf{B})^{\top} = \mathbf{A}^{\top} + \mathbf{B}^{\top}$ より、括弧を展開します。
4: 転置の性質 $(\mathbf{A} \mathbf{B})^{\top} = \mathbf{B}^{\top} \mathbf{A}^{\top}$ より、括弧を展開します。
4: $\hat{\boldsymbol{\beta}}_{\mathrm{ML}}(\lambda)$ に、式(0.38)を代入します。

　4番目の項は、次のように変形できます。

$\displaystyle \begin{aligned} \Bigl( \mathbf{X}(\lambda) \hat{\boldsymbol{\beta}}_{\mathrm{ML}}(\lambda) \Bigr)^{\top} \mathbf{X}(\lambda) \hat{\boldsymbol{\beta}}_{\mathrm{ML}}(\lambda) &= \hat{\boldsymbol{\beta}}_{\mathrm{ML}}(\lambda)^{\top} \mathbf{X}(\lambda)^{\top} \mathbf{X}(\lambda) \hat{\boldsymbol{\beta}}_{\mathrm{ML}}(\lambda) \\ &= \hat{\boldsymbol{\beta}}_{\mathrm{ML}}(\lambda)^{\top} \mathbf{X}(\lambda)^{\top} \mathbf{X}(\lambda) \Bigl( \mathbf{X}(\lambda)^{\top} \mathbf{X}(\lambda) \Bigr)^{-1} \mathbf{X}(\lambda)^{\top} \mathbf{y}(\lambda) \\ &= \hat{\boldsymbol{\beta}}_{\mathrm{ML}}(\lambda)^{\top} \mathbf{X}(\lambda)^{\top} \mathbf{y}(\lambda) \end{aligned}$

4: 2から4番目の項は二次形式(スカラ)なので、転置できます。

　後2つの項が打ち消し合います。

$\displaystyle \begin{aligned} \Bigl( \mathbf{X}(\lambda) \hat{\boldsymbol{\beta}}_{\mathrm{ML}}(\lambda) \Bigr)^{\top} \mathbf{y}(\lambda) &= \hat{\boldsymbol{\beta}}_{\mathrm{ML}}(\lambda)^{\top} \mathbf{X}(\lambda)^{\top} \mathbf{y}(\lambda) \\ &= \Bigl( \hat{\boldsymbol{\beta}}_{\mathrm{ML}}(\lambda)^{\top} \mathbf{X}(\lambda)^{\top} \mathbf{y}(\lambda) \Bigr)^{\top} \\ &= \mathbf{y}(\lambda)^{\top} \mathbf{X}(\lambda) \hat{\boldsymbol{\beta}}_{\mathrm{ML}}(\lambda) \end{aligned}$

5: $\mathbf{y}(\lambda)$ を括り出します。
6: 括弧を、式(9)で置き換えます。

　以上で、分散パラメータの最尤推定量の式が得られました。

集約対数尤度関数

　対数尤度関数の式において空間パラメータが行列式や逆行列の中に含まれるため、空間パラメータの最尤推定量は、解析的に求められません。そこで、集約対数尤度関数(集中対数尤度関数)が最大となる値を最尤推定量とします。

　独立な変数としてのパラメータ $\boldsymbol{\beta}, \sigma^2$ を $\lambda$ の関数である最尤推定量 $\hat{\boldsymbol{\beta}}_{\mathrm{ML}}(\lambda), \hat{\sigma}_{\mathrm{ML}}^2(\lambda)$ に置き換えることで、 $\boldsymbol{\beta}, \sigma^2, \lambda$ の関数である対数尤度関数 $\log L(\boldsymbol{\theta})$ を $\lambda$ の関数である集約対数尤度関数 $\log L_{\mathrm{c}}(\hat{\boldsymbol{\beta}}_{\mathrm{ML}}, \hat{\sigma}_{\mathrm{ML}}^2, \lambda)$ として扱います。

$\begin{align} \log L_{\mathrm{c}} \Bigl( \hat{\boldsymbol{\beta}}_{\mathrm{ML}}(\lambda), \hat{\sigma}_{\mathrm{ML}}^2(\lambda), \lambda \Bigr) &= - \frac{N \log (2 \pi)}{2} - \frac{N}{2} \log \hat{\sigma}_{\mathrm{ML}}^2 ( \hat{\boldsymbol{\beta}}_{\mathrm{ML}}, \lambda ) - \frac{1}{2} \frac{ \boldsymbol{\epsilon} ( \hat{\boldsymbol{\beta}}_{\mathrm{ML}}, \lambda )^{\top} \boldsymbol{\epsilon} ( \hat{\boldsymbol{\beta}}_{\mathrm{ML}}, \lambda ) }{ \hat{\sigma}_{\mathrm{ML}}^2 ( \hat{\boldsymbol{\beta}}_{\mathrm{ML}}, \lambda ) } + \log |\mathbf{A}| \tag{0.37'}\\ &= - \frac{N \log (2 \pi)}{2} - \frac{N}{2} \log \Bigl( \frac{1}{N} \boldsymbol{\epsilon} ( \hat{\boldsymbol{\beta}}_{\mathrm{ML}}, \lambda )^{\top} \boldsymbol{\epsilon} ( \hat{\boldsymbol{\beta}}_{\mathrm{ML}}, \lambda ) \Bigr) - \frac{1}{2} \frac{ N \hat{\sigma}_{\mathrm{ML}}^2 ( \hat{\boldsymbol{\beta}}_{\mathrm{ML}}, \lambda ) }{ \hat{\sigma}_{\mathrm{ML}}^2 ( \hat{\boldsymbol{\beta}}_{\mathrm{ML}}, \lambda ) } + \log |\mathbf{A}| \\ &= - \frac{N}{2} - \frac{N}{2} \log (2 \pi) - \frac{N}{2} \log \frac{1}{N} - \frac{N}{2} \log \Bigl( \boldsymbol{\epsilon} ( \hat{\boldsymbol{\beta}}_{\mathrm{ML}}, \lambda )^{\top} \boldsymbol{\epsilon} ( \hat{\boldsymbol{\beta}}_{\mathrm{ML}}, \lambda ) \Bigr) + \log |\mathbf{A}| \\ &= - \frac{N}{2} \Biggl( 1 + \log \Bigl( \frac{2 \pi}{N} \Bigr) \Biggr) - \frac{N}{2} \log \Bigl( \boldsymbol{\epsilon} ( \hat{\boldsymbol{\beta}}_{\mathrm{ML}}, \lambda )^{\top} \boldsymbol{\epsilon} ( \hat{\boldsymbol{\beta}}_{\mathrm{ML}}, \lambda ) \Bigr) + \log |\mathbf{A}| \tag{0.41}\\ &= - \frac{N}{2} \Biggl( 1 + \log \Bigl( \frac{2 \pi}{N} \Bigr) \Biggr) - \frac{N}{2} \log \Bigl( \mathbf{y}(\lambda)^{\top} \boldsymbol{\epsilon} \Bigl( \hat{\boldsymbol{\beta}}_{\mathrm{ML}}, \lambda \Bigr) \Bigr) + \log |\mathbf{A}| \end{align}$

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

1: 式(0.37)について、 $\sigma^2$ を式(0.39')(10)、 $\boldsymbol{\epsilon}$ を式(2')で置き換えた式を立てます。
2: $\hat{\sigma}_{\mathrm{ML}}^2(\hat{\boldsymbol{\beta}}_{\mathrm{ML}}, \lambda)$ に、式(10)を代入します。
2: 式(0.39')より、 $N \hat{\sigma}_{\mathrm{ML}}^2(\hat{\boldsymbol{\beta}}_{\mathrm{ML}}, \lambda) = \boldsymbol{\epsilon}(\hat{\boldsymbol{\beta}}_{\mathrm{ML}}, \lambda)^{\top} \boldsymbol{\epsilon}(\hat{\boldsymbol{\beta}}_{\mathrm{ML}}, \lambda)$ で置き換えます。
3: 対数の性質 $\log(x y) = \log x + \log y$ より、対数の和に分割します。
4: 対数の性質 $\log(x y) = \log x + \log y$ より、積の対数にまとめます。
5: $\boldsymbol{\epsilon}(\hat{\boldsymbol{\beta}}_{\mathrm{ML}}, \lambda)^{\top} \boldsymbol{\epsilon}(\hat{\boldsymbol{\beta}}_{\mathrm{ML}}, \lambda)$ を、式(0.39')(10)で置き換えます。(確認用)

　集約対数尤度関数の式が得られました。
　3つのパラメータ $\boldsymbol{\beta}, \sigma^2, \lambda$ の関数 $\log L(\boldsymbol{\theta})$ から、1つのパラメータ $\lambda$ の関数 $\log L_{\mathrm{c}}(\hat{\boldsymbol{\beta}}_{\mathrm{ML}}, \hat{\sigma}_{\mathrm{ML}}^2, \lambda)$ に変数(次元)を集約することで、 $\lambda$ の最適化問題に帰着させることができます。

　集約対数尤度関数(0.41)を最大化する $\lambda$ を最尤推定量 $\hat{\lambda}_{\mathrm{ML}}$ とします。

$\displaystyle \begin{aligned} \hat{\lambda}_{\mathrm{ML}} &= \mathop{\mathrm{argmax}}\limits_{\lambda}\ L_{\mathrm{c}} \Bigl( \hat{\boldsymbol{\beta}}_{\mathrm{ML}}, \hat{\sigma}_{\mathrm{ML}}^2, \lambda \Bigr) \\ &= \mathop{\mathrm{argmax}}\limits_{\lambda}\ \log L_{\mathrm{c}} \Bigl( \hat{\boldsymbol{\beta}}_{\mathrm{ML}}, \hat{\sigma}_{\mathrm{ML}}^2, \lambda \Bigr) \end{aligned}$

　対数関数は単調増加するので、尤度関数を最大化する値と対数尤度関数を最大化する値は一致します。
　集約対数尤度関数を最大化する値は、解析的に求められない(解析的最適化できない)ので、探索アルゴリズムなど(数値最適化)により求めます。

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

　この記事では、SEMの最尤法を数式で確認しました。

参考文献

Pythonによるはじめての地理空間データサイエンス

作者:植杉大
森北出版

Amazon

おわりに

　10章あるいは付録の導出編が完了です。おつかれさまでした～。

　久し振りにまともに予約投稿を仕込んで更新を続けられたので心地いい日々を過ごせております。
　投稿の予約時には気付きませんでしたがこの記事が投稿された直後のはてなのダッシュボードによると、この記事で777記事目だそうです。目次ページのような内容の無い記事も含めてですけど。
　ここまでになると、記事を書くよりも管理する方が大変なんです。それらが、最近のはてなの仕様変更(あるいはそれに伴うバグ)により一部の数式がこれまで通りに表示されなくなっており、その対応が本当にもう……

　最後に、Juice=Juiceのライブ映像をどうぞ♪

【次の内容】

からっぽのしょこ

はてなブログの仕様変更の影響で、数式の表示が崩壊しています。読んだら書く！書いたら読む！同じ事は二度調べ(たく)ない

0.2.2：空間誤差モデル(SEM)の最尤法の導出【はじめての地理空間DSのノート】

はじめに

0.2.2 空間誤差モデル(SEM)の最尤法の導出

モデルの確認

定義式

計算式

尤度関数

対数尤度関数の微分

パラメータの推定

回帰パラメータの最尤推定量

分散パラメータの最尤推定量

集約対数尤度関数

参考文献

おわりに