はじめに

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

　この記事では、SEMの漸近分散共分散行列について、数式を使って解説します。

【前の内容】

www.anarchive-beta.com

【他の内容】

www.anarchive-beta.com

【今回の内容】

はじめに
0.2.2 空間誤差モデル(SEM)の漸近分散共分散行列の導出
- モデルの確認
- 漸近分散共分散行列の導出
参考文献
おわりに

0.2.2 空間誤差モデル(SEM)の漸近分散共分散行列の導出

　空間誤差モデル(SEM・Spatial Error Model)における漸近分散共分散行列(asymptotic variance-covariance matrix)を導出します。
　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の定義式」を参照してください。

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

$\begin{align} \mathbf{A} &= \mathbf{I} - \lambda \mathbf{W} \tag{1}\\ \mathbf{B} &= \mathbf{W} \mathbf{A}^{-1} \tag{2} \end{align}$

　定義式(0.36)を被説明変数 $\mathbf{y}$ について整理すると、次の式になります。

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

　定義式(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{4} \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{5} \end{align}$

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

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

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

ヘッセ行列

　SEMのヘッセ行列については「SEMのヘッセ行列の導出」を参照してください。

　ヘッセ行列は、次の式となります。

$\begin{align} \mathbf{H}(\boldsymbol{\theta}) &= \frac{\partial^2 \log L(\boldsymbol{\theta})}{\partial \boldsymbol{\theta} \partial \boldsymbol{\theta}^{\top}} \\ &= \begin{pmatrix} \frac{\partial^2 \log L(\boldsymbol{\theta})}{\partial \boldsymbol{\beta} \partial \boldsymbol{\beta}^{\top}} & \frac{\partial^2 \log L(\boldsymbol{\theta})}{\partial \boldsymbol{\beta} \partial \sigma^2} & \frac{\partial^2 \log L(\boldsymbol{\theta})}{\partial \boldsymbol{\beta} \partial \lambda} \\ \frac{\partial^2 \log L(\boldsymbol{\theta})}{\partial \sigma^2 \partial \boldsymbol{\beta}^{\top}} & \frac{\partial^2 \log L(\boldsymbol{\theta})}{\partial \sigma^2 \partial \sigma^2} & \frac{\partial^2 \log L(\boldsymbol{\theta})}{\partial \sigma^2 \partial \lambda} \\ \frac{\partial^2 \log L(\boldsymbol{\theta})}{\partial \lambda \partial \boldsymbol{\beta}^{\top}} & \frac{\partial^2 \log L(\boldsymbol{\theta})}{\partial \lambda \partial \sigma^2} & \frac{\partial^2 \log L(\boldsymbol{\theta})}{\partial \lambda \partial \lambda} \end{pmatrix} \\ &= \begin{pmatrix} - \frac{1}{\sigma^2} \mathbf{X}^{\top} \mathbf{A}^{\top} \mathbf{A} \mathbf{X} & - \frac{1}{(\sigma^2)^2} \mathbf{X}^{\top} \mathbf{A}^{\top} \boldsymbol{\epsilon} & - \frac{1}{\sigma^2} \Bigl( \mathbf{X}^{\top} (\mathbf{A}^{\top} \mathbf{W} + \mathbf{W}^{\top} \mathbf{A}) (\mathbf{y} - \mathbf{X} \boldsymbol{\beta}) \Bigr) \\ \cdot & \frac{N}{2} \frac{1}{(\sigma^2)^2} - \frac{1}{(\sigma^2)^3} \boldsymbol{\epsilon}^{\top} \boldsymbol{\epsilon} & - \frac{1}{(\sigma^2)^2} \boldsymbol{\epsilon}^{\top} \mathbf{W} (\mathbf{y} - \mathbf{X} \boldsymbol{\beta}) \\ \cdot & \cdot & - \mathrm{Tr}(\mathbf{A}^{-1} \mathbf{W} \mathbf{A}^{-1} \mathbf{W}) - \frac{1}{\sigma^2} (\mathbf{y} - \mathbf{X} \boldsymbol{\beta})^{\top} \mathbf{W}^{\top} \mathbf{W} (\mathbf{y} - \mathbf{X} \boldsymbol{\beta}) \end{pmatrix} \tag{6} \end{align}$

　以上の式を用いて漸近分散共分散行列を求めます。

スポンサードリンク

漸近分散共分散行列の導出

　次は、SEMの対数尤度関数の漸近分散共分散行列を導出します。

フィッシャー情報行列の設定

　対数尤度関数のフィッシャー情報行列を確認します。

　対数尤度関数 $\log L(\boldsymbol{\theta})$ のヘッセ行列の負の期待値を求めます。ヘッセ行列 $\mathbf{H}(\boldsymbol{\theta})$ の負の期待値をフィッシャー情報行列 $\mathbb{I}(\boldsymbol{\theta})$ と呼びます。

$\begin{align} \mathbb{I}(\boldsymbol{\theta}) &\equiv - \mathbb{E}[\mathbf{H}(\boldsymbol{\theta})] \\ &= - \mathbb{E} \left[ \begin{pmatrix} \frac{\partial^2 \log L(\boldsymbol{\theta})}{\partial \boldsymbol{\beta} \partial \boldsymbol{\beta}^{\top}} & \frac{\partial^2 \log L(\boldsymbol{\theta})}{\partial \boldsymbol{\beta} \partial \sigma^2} & \frac{\partial^2 \log L(\boldsymbol{\theta})}{\partial \boldsymbol{\beta} \partial \lambda} \\ \frac{\partial^2 \log L(\boldsymbol{\theta})}{\partial \sigma^2 \partial \boldsymbol{\beta}^{\top}} & \frac{\partial^2 \log L(\boldsymbol{\theta})}{\partial \sigma^2 \partial \sigma^2} & \frac{\partial^2 \log L(\boldsymbol{\theta})}{\partial \sigma^2 \partial \lambda} \\ \frac{\partial^2 \log L(\boldsymbol{\theta})}{\partial \lambda \partial \boldsymbol{\beta}^{\top}} & \frac{\partial^2 \log L(\boldsymbol{\theta})}{\partial \lambda \partial \sigma^2} & \frac{\partial^2 \log L(\boldsymbol{\theta})}{\partial \lambda \partial \lambda} \end{pmatrix} \right] \\ &= - \begin{pmatrix} \mathbb{E}[\frac{\partial^2 \log L(\boldsymbol{\theta})}{\partial \boldsymbol{\beta} \partial \boldsymbol{\beta}^{\top}}] & \mathbb{E}[\frac{\partial^2 \log L(\boldsymbol{\theta})}{\partial \boldsymbol{\beta} \partial \sigma^2}] & \mathbb{E}[\frac{\partial^2 \log L(\boldsymbol{\theta})}{\partial \boldsymbol{\beta} \partial \lambda}] \\ \mathbb{E}[\frac{\partial^2 \log L(\boldsymbol{\theta})}{\partial \sigma^2 \partial \boldsymbol{\beta}^{\top}}] & \mathbb{E}[\frac{\partial^2 \log L(\boldsymbol{\theta})}{\partial \sigma^2 \partial \sigma^2}] & \mathbb{E}[\frac{\partial^2 \log L(\boldsymbol{\theta})}{\partial \sigma^2 \partial \lambda}] \\ \mathbb{E}[\frac{\partial^2 \log L(\boldsymbol{\theta})}{\partial \lambda \partial \boldsymbol{\beta}^{\top}}] & \mathbb{E}[\frac{\partial^2 \log L(\boldsymbol{\theta})}{\partial \lambda \partial \sigma^2}] & \mathbb{E}[\frac{\partial^2 \log L(\boldsymbol{\theta})}{\partial \lambda \partial \lambda}] \end{pmatrix} \tag{7} \end{align}$

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

1: 式(6)の符号を反転した期待値の式を立てます。
2: 行列の要素を明示します。
3: 行列の期待値を、期待値の行列に変形します。

　フィッシャー情報行列の各要素は、ヘッセ行列の各要素(各パラメータ $\boldsymbol{\beta}, \sigma^2, \lambda$ の組み合わせによる偏微分)の期待値で求まるのが分かります。

　フィッシャー情報行列(6)の各要素の式を求めていきます。

回帰・回帰パラメータによる偏微分の期待値

　対数尤度関数の回帰パラメータに関する2階微分の期待値を求めます。

　フィッシャー情報行列(7)の1行1列目の要素は、対数尤度関数の空間パラメータ $\boldsymbol{\beta}$ に関する2階微分の期待値をとって求まります。

$\begin{align} \mathbb{E} \Biggl[ \frac{\partial^2 \log L(\boldsymbol{\theta})}{\partial \boldsymbol{\beta} \partial \boldsymbol{\beta}^{\top}} \Biggr] &= \mathbb{E} \Biggl[ - \frac{1}{\sigma^2} \mathbf{X}^{\top} \mathbf{A}^{\top} \mathbf{A} \mathbf{X} \Biggr] \tag{0.42'}\\ &= - \frac{1}{\sigma^2} \mathbf{X}^{\top} \mathbf{A}^{\top} \mathbf{A} \mathbf{X} \tag{0.48} \end{align}$

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

1: 式(0.42)より、式(7)の $(1, 1)$ 成分の期待値の式を立てます。 $\boldsymbol{\epsilon}$ についての期待値なので、 $\boldsymbol{\epsilon}$ 以外の項は定数として扱います。
2: 期待値の性質 $\mathbb{E}[a] = a$ より、期待値を外します。

　フィッシャー情報行列(7)の1行1列目の要素の式が得られました。

回帰・分散パラメータによる偏微分の期待値

　対数尤度関数の回帰パラメータと分散パラメータに関する微分の期待値を求めます。

　フィッシャー情報行列(7)の1行2列目の要素は、対数尤度関数の回帰パラメータ $\boldsymbol{\beta}$ と分散パラメータ $\sigma^2$ に関する微分の期待値をとって求まります。

$\begin{align} \mathbb{E} \Biggl[ \frac{\partial^2 \log L(\boldsymbol{\theta})}{\partial \boldsymbol{\beta} \partial \sigma^2} \Biggr] &= \mathbb{E} \Biggl[ - \frac{1}{(\sigma^2)^2} \mathbf{X}^{\top} \mathbf{A}^{\top} \boldsymbol{\epsilon} \Biggr] \tag{0.43'}\\ &= - \frac{1}{(\sigma^2)^2} \mathbf{X}^{\top} \mathbf{A}^{\top} \mathbb{E}[\boldsymbol{\epsilon}] \\ &= - \frac{1}{(\sigma^2)^2} \mathbf{X}^{\top} \mathbf{A}^{\top} \mathbf{0} \\ &= \mathbf{0} \tag{0.49} \end{align}$

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

1: 式(0.43)より、式(7)の $(1, 2)$ 成分の期待値の式を立てます。 $\boldsymbol{\epsilon}$ についての期待値なので、 $\boldsymbol{\epsilon}$ 以外の項は定数として扱います。
2: 期待値の性質 $\mathbb{E}[a x] = a \mathbb{E}[x]$ より、係数を期待値の外に出します。
3: $\boldsymbol{\epsilon}$ の期待値に、式(8)を代入します。

　誤差項 $\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{8} \end{align}$

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

1: ベクトルの要素を明示します。
2: ベクトルの期待値を、期待値のベクトルに変形します。
3: SEMの定義(0.12.b)より、誤差項の平均 $\mathbb{E}[\epsilon_n] = \mu = 0$ で置き換えます。
3-4: 0ベクトルとなります。

　以上で、フィッシャー情報行列(7)の1行2列目の要素の式が得られました。

回帰・空間パラメータによる偏微分の期待値

　対数尤度関数の回帰パラメータと空間パラメータに関する微分の期待値を求めます。

　フィッシャー情報行列(7)の1行3列目の要素は、対数尤度関数の回帰パラメータ $\boldsymbol{\beta}$ と空間パラメータ $\lambda$ に関する微分の期待値をとって求まります。

$\begin{align} \mathbb{E} \Biggl[ \frac{\partial^2 \log L(\boldsymbol{\theta})}{\partial \boldsymbol{\beta} \partial \lambda} \Biggr] &= \mathbb{E} \Biggl[ - \frac{1}{\sigma^2} \Bigl( \mathbf{X}^{\top} (\mathbf{A}^{\top} \mathbf{W} + \mathbf{W}^{\top} \mathbf{A}) (\mathbf{y} - \mathbf{X} \boldsymbol{\beta}) \Bigr) \Biggr] \tag{0.44'}\\ &= \mathbb{E} \Biggl[ - \frac{1}{\sigma^2} \Bigl( \mathbf{X}^{\top} (\mathbf{A}^{\top} \mathbf{W} + \mathbf{W}^{\top} \mathbf{A}) \mathbf{y} - \mathbf{X}^{\top} (\mathbf{A}^{\top} \mathbf{W} + \mathbf{W}^{\top} \mathbf{A}) \mathbf{X} \boldsymbol{\beta} \Bigr) \Biggr] \\ &= - \frac{1}{\sigma^2} \Bigl( \mathbf{X}^{\top} (\mathbf{A}^{\top} \mathbf{W} + \mathbf{W}^{\top} \mathbf{A}) \mathbb{E}[\mathbf{y}] - \mathbb{E} \Bigl[ \mathbf{X}^{\top} (\mathbf{A}^{\top} \mathbf{W} + \mathbf{W}^{\top} \mathbf{A}) \mathbf{X} \boldsymbol{\beta} \Bigl] \Bigr) \\ &= - \frac{1}{\sigma^2} \Bigl( \mathbf{X}^{\top} (\mathbf{A}^{\top} \mathbf{W} + \mathbf{W}^{\top} \mathbf{A}) \mathbf{X} \boldsymbol{\beta} - \mathbf{X}^{\top} (\mathbf{A}^{\top} \mathbf{W} + \mathbf{W}^{\top} \mathbf{A}) \mathbf{X} \boldsymbol{\beta} \Bigr) \\ &= - \frac{1}{\sigma^2} \mathbf{0} \\ &= \mathbf{0} \tag{0.50} \end{align}$

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

1: 式(0.44)より、式(7)の $(1, 3)$ 成分の期待値の式を立てます。 $\boldsymbol{\epsilon}$ についての期待値なので、 $\boldsymbol{\epsilon}$ 以外の項は定数として扱います。
2: 括弧を展開します。
3: 期待値の性質 $\mathbb{E}[a x] = a \mathbb{E}[x]$ より、係数を期待値の外に出します。
3: 期待値の性質 $\mathbb{E}[x + y] = \mathbb{E}[x] + \mathbb{E}[y]$ より、期待値の和に分割します。
4: $\mathbf{y}$ の期待値に、式(9)を代入します。
4: 期待値の性質 $\mathbb{E}[a] = a$ より、期待値を外します。

　被説明変数 $\mathbf{y}$ の期待値をとります。

$\begin{align} \mathbb{E}[\mathbf{y}] &= \mathbb{E} \Bigl[ \mathbf{X} \boldsymbol{\beta} + \mathbf{A}^{-1} \boldsymbol{\epsilon} \Bigr] \\ &= \mathbb{E}[\mathbf{X} \boldsymbol{\beta}] + \mathbf{A}^{-1} \mathbb{E}[\boldsymbol{\epsilon}] \\ &= \mathbf{X} \boldsymbol{\beta} + \mathbf{A}^{-1} \mathbf{0} \\ &= \mathbf{X} \boldsymbol{\beta} \tag{9} \end{align}$

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

1: $\mathbf{y}$ に、式(3)を代入します。
2: 期待値の性質 $\mathbb{E}[x + y] = \mathbb{E}[x] + \mathbb{E}[y]$ より、期待値の和に分割します。
2: 期待値の性質 $\mathbb{E}[a x] = a \mathbb{E}[x]$ より、係数を期待値の外に出します。
3: 期待値の性質 $\mathbb{E}[a] = a$ より、期待値を外します。
3: $\boldsymbol{\epsilon}$ の期待値に、式(8)を代入します。

　フィッシャー情報行列(7)の1行3列目の要素の式が得られました。

分散・分散パラメータによる偏微分の期待値

　対数尤度関数の分散パラメータに関する2階微分の期待値を求めます。

　フィッシャー情報行列(7)の2行2列目の要素は、対数尤度関数の分散パラメータ $\sigma^2$ に関する2階微分の期待値をとって求まります。

$\begin{align} \mathbb{E} \Biggl[ \frac{\partial \log L(\boldsymbol{\theta})}{\partial \sigma^2 \partial \sigma^2} \Biggr] &= \mathbb{E} \Biggl[ \frac{N}{2} \frac{1}{(\sigma^2)^2} - \frac{1}{(\sigma^2)^3} \boldsymbol{\epsilon}^{\top} \boldsymbol{\epsilon} \Biggr] \tag{0.45'}\\ &= \mathbb{E} \Biggl[ \frac{N}{2} \frac{1}{(\sigma^2)^2} \Biggr] - \frac{1}{(\sigma^2)^3} \mathbb{E}[\boldsymbol{\epsilon}^{\top} \boldsymbol{\epsilon}] \\ &= \frac{N}{2} \frac{1}{(\sigma^2)^2} - \frac{1}{(\sigma^2)^3} N \sigma^2 \\ &= \frac{N}{2} \frac{1}{(\sigma^2)^2} - N \frac{1}{(\sigma^2)^2} \\ &= - \frac{N}{2} \frac{1}{(\sigma^2)^2} \tag{0.51} \end{align}$

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

1: 式(0.45)より、式(7)の $(2, 2)$ 成分の期待値の式を立てます。 $\boldsymbol{\epsilon}$ についての期待値なので、 $\boldsymbol{\epsilon}$ 以外の項は定数として扱います。
2: 期待値の性質 $\mathbb{E}[x + y] = \mathbb{E}[x] + \mathbb{E}[y]$ より、期待値の和に分割します。
2: 期待値の性質 $\mathbb{E}[a x] = a \mathbb{E}[x]$ より、係数を期待値の外に出します。
3: 期待値の性質 $\mathbb{E}[a] = a$ より、期待値を外します。
3: $\boldsymbol{\epsilon}$ の内積の期待値に、式(10)を代入します。

　 $\boldsymbol{\epsilon}$ の内積の期待値をとります。

$\begin{align} \mathbb{E}[\boldsymbol{\epsilon}^{\top} \boldsymbol{\epsilon}] &= \mathbb{E} \Bigl[ \mathrm{Tr}( \boldsymbol{\epsilon} \boldsymbol{\epsilon}^{\top} ) \Bigr] \\ &= \mathrm{Tr} \Bigl( \mathbb{E}[ \boldsymbol{\epsilon} \boldsymbol{\epsilon}^{\top} ] \Bigr) \\ &= \mathrm{Tr}(\sigma^2 \mathbf{I}) \\ &= \sigma^2 \mathrm{Tr}(\mathbf{I}) \\ &= N \sigma^2 \tag{10} \end{align}$

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

1: トレースの性質 $\mathbf{x}^{\top} \mathbf{x} = \mathrm{Tr}(\mathbf{x} \mathbf{x}^{\top})$ より、内積をトレースに置き換えます。
2: トレースの性質 $\mathrm{Tr}(\mathbb{E}[\mathbf{X}]) = \mathbb{E}[\mathrm{Tr}(\mathbf{X})]$ より、トレースと期待値の順番を入れ換えます。
3: $\boldsymbol{\epsilon}$ の積の期待値に、式(11)を代入します。
4: トレースの性質 $\mathrm{Tr}(a \mathbf{X}) = a \mathrm{Tr}(\mathbf{X})$ より、係数をトレースの外に出します。
5: トレースの性質 $\mathrm{Tr}(\mathbf{I}_n) = n$ より、トレースが対角要素数になります。

　誤差項 $\boldsymbol{\epsilon}$ の積の期待値をとります。

$\begin{align} \mathbb{E}[\boldsymbol{\epsilon} \boldsymbol{\epsilon}^{\top}] &= \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{11} \end{align}$

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

1: ベクトルの要素を明示します。
2: ベクトルの積を計算します。
3: 行列の期待値を、期待値の行列に変形します。
4: SEMの定義(0.12.b)より、誤差項の分散 $\mathbb{E}[\epsilon_n \epsilon_n] = \sigma^2$ 、共分散 $\mathbb{E}[\epsilon_i \epsilon_j] = \sigma_{ij} = 0$ で置き換えます。
5-6: $\sigma^2$ と単位行列の積になります。

　以上で、フィッシャー情報行列(6)の2行2列目の要素の式が得られました。

分散・空間パラメータによる偏微分の期待値

　対数尤度関数の分散パラメータと空間パラメータに関する微分の期待値を求めます。

　フィッシャー情報行列(7)の2行3列目の要素は、対数尤度関数の分散パラメータ $\sigma^2$ と空間パラメータ $\lambda$ に関する微分の期待値をとって求まります。

$\begin{align} \mathbb{E} \Biggl[ \frac{\partial \log L(\boldsymbol{\theta})}{\partial \sigma^2 \partial \lambda} \Biggr] &= \mathbb{E} \Biggl[ - \frac{1}{(\sigma^2)^2} \boldsymbol{\epsilon}^{\top} \mathbf{W} (\mathbf{y} - \mathbf{X} \boldsymbol{\beta}) \Biggr] \tag{0.46'}\\ &= - \frac{1}{(\sigma^2)^2} \mathbb{E} \Bigl[ \boldsymbol{\epsilon}^{\top} \mathbf{W} \Bigl( \mathbf{X} \boldsymbol{\beta} + \mathbf{A}^{-1} \boldsymbol{\epsilon} - \mathbf{X} \boldsymbol{\beta} \Bigr) \Bigr] \\ &= - \frac{1}{(\sigma^2)^2} \mathbb{E} \Bigl[ \boldsymbol{\epsilon}^{\top} \mathbf{W} \mathbf{A}^{-1} \boldsymbol{\epsilon} \Bigr] \\ &= - \frac{1}{(\sigma^2)^2} \mathbb{E} \Bigl[ \boldsymbol{\epsilon}^{\top} \mathbf{B} \boldsymbol{\epsilon} \Bigr] \\ &= - \frac{1}{(\sigma^2)^2} \sigma^2 \mathrm{Tr}(\mathbf{B}) \\ &= - \frac{1}{\sigma^2} \mathrm{Tr}(\mathbf{B}) \tag{0.52} \end{align}$

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

1: 式(0.46)より、式(7)の $(2, 3)$ 成分の期待値の式を立てます。 $\boldsymbol{\epsilon}$ についての期待値なので、 $\boldsymbol{\epsilon}$ 以外の項は定数として扱います。
2: 期待値の性質 $\mathbb{E}[a x] = a \mathbb{E}[x]$ より、係数を期待値の外に出します。
2: $\mathbf{y}$ に、式(3)を代入します。
3: 括弧を展開します。
4: $\mathbf{W} \mathbf{A}^{-1}$ を、式(2)で置き換えます。
5: $\boldsymbol{\epsilon}, \mathbf{B}$ の二次形式の期待値に、式(12)を代入します。

　 $\boldsymbol{\epsilon}, \mathbf{B}$ の二次形式の期待値をとります。

$\begin{align} \mathbb{E}[ \boldsymbol{\epsilon}^{\top} \mathbf{B} \boldsymbol{\epsilon} \Bigr] &= \mathbb{E} \Bigl[ \mathrm{Tr} ( \mathbf{B} \boldsymbol{\epsilon} \boldsymbol{\epsilon}^{\top} ) \Bigr] \\ &= \mathrm{Tr} \Bigl( \mathbb{E} [ \mathbf{B} \boldsymbol{\epsilon} \boldsymbol{\epsilon}^{\top} ] \Bigr) \\ &= \mathrm{Tr} \Bigl( \mathbf{B} \mathbb{E}[\boldsymbol{\epsilon} \boldsymbol{\epsilon}^{\top}] \Bigr) \\ &= \mathrm{Tr} \Bigl( \mathbf{B} \sigma^2 \mathbf{I} \Bigr) \\ &= \sigma^2 \mathrm{Tr}(\mathbf{B}) \tag{12} \end{align}$

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

1: トレースの性質 $\mathbf{x}^{\top} \mathbf{A} \mathbf{x} = \mathrm{Tr}(\mathbf{A} \mathbf{x} \mathbf{x}^{\top})$ より、二次形式をトレースに置き換えます。
2: トレースの性質 $\mathrm{Tr}(\mathbb{E}[\mathbf{X}]) = \mathbb{E}[\mathrm{Tr}(\mathbf{X})]$ より、トレースと期待値の順番を入れ換えます。
3: 期待値の性質 $\mathbb{E}[a x] = a \mathbb{E}[x]$ より、係数を期待値の外に出します。
4: $\boldsymbol{\epsilon}$ の積の期待値に、式(11)を代入します。
5: トレースの性質 $\mathrm{Tr}(a \mathbf{X}) = a \mathrm{Tr}(\mathbf{X})$ より、係数をトレースの外に出します。

　以上で、フィッシャー情報行列(7)の2行3列目の要素の式が得られました。

空間・空間パラメータによる偏微分の期待値

　対数尤度関数の空間パラメータに関する2階微分の期待値を求めます。

　フィッシャー情報行列(7)の3行3列目の要素は、対数尤度関数の空間パラメータ $\lambda$ に関する2階微分の期待値をとって求まります。

$\begin{align} \mathbb{E} \Biggl[ \frac{\partial \log L(\boldsymbol{\theta})}{\partial \lambda \partial \lambda} \Biggr] &= \mathbb{E} \Biggl[ - \mathrm{Tr}(\mathbf{A}^{-1} \mathbf{W} \mathbf{A}^{-1} \mathbf{W}) - \frac{1}{\sigma^2} (\mathbf{y} - \mathbf{X} \boldsymbol{\beta})^{\top} \mathbf{W}^{\top} \mathbf{W} (\mathbf{y} - \mathbf{X} \boldsymbol{\beta}) \Biggr] \tag{0.47'}\\ &= - \mathbb{E} \Bigl[ \mathrm{Tr}(\mathbf{A}^{-1} \mathbf{W} \mathbf{A}^{-1} \mathbf{W}) \Bigr] - \frac{1}{\sigma^2} \mathbb{E} \Bigl[ (\mathbf{y} - \mathbf{X} \boldsymbol{\beta})^{\top} \mathbf{W}^{\top} \mathbf{W} (\mathbf{y} - \mathbf{X} \boldsymbol{\beta}) \Bigr] \\ &= - \mathrm{Tr}(\mathbf{A}^{-1} \mathbf{W} \mathbf{A}^{-1} \mathbf{W}) - \frac{1}{\sigma^2} \sigma^2 \mathrm{Tr}(\mathbf{B}^{\top} \mathbf{B}) \\ &= - \mathrm{Tr}(\mathbf{B} \mathbf{B}) - \mathrm{Tr}(\mathbf{B}^{\top} \mathbf{B}) \\ &= - \mathrm{Tr} \Bigl( \mathbf{B} \mathbf{B} + \mathbf{B}^{\top} \mathbf{B} \Bigr) \\ &= - \mathrm{Tr} \Bigl( (\mathbf{B} + \mathbf{B}^{\top}) \mathbf{B} \Bigr) \tag{0.53} \end{align}$

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

1: 式(0.47)より、式(7)の $(3, 3)$ 成分の期待値の式を立てます。 $\boldsymbol{\epsilon}$ についての期待値なので、 $\boldsymbol{\epsilon}$ 以外の項は定数として扱います。
2: 期待値の性質 $\mathbb{E}[x + y] = \mathbb{E}[x] + \mathbb{E}[y]$ より、期待値の和に分割します。
2: 期待値の性質 $\mathbb{E}[a x] = a \mathbb{E}[x]$ より、係数を期待値の外に出します。
3: 期待値の性質 $\mathbb{E}[a] = a$ より、期待値を外します。
3: $\mathbf{X}, \mathbf{y}, \mathbf{W}$ などの二次形式の期待値に、式(13)を代入します。
4: トレースの性質 $\mathrm{Tr}(\mathbf{A} \mathbf{B} \mathbf{C}) = \mathrm{Tr}(\mathbf{B} \mathbf{C} \mathbf{A})$ より、トレース内の行列の積の順番を入れ換えます。
4: $\mathbf{W} \mathbf{A}^{-1}$ を、式(2)で置き換えます。

　トレースの項は、次のように変形できます。

$\displaystyle \begin{aligned} \mathrm{Tr}(\mathbf{A}^{-1} \mathbf{W} \mathbf{A}^{-1} \mathbf{W}) &= \mathrm{Tr}(\mathbf{W} \mathbf{A}^{-1} \mathbf{W} \mathbf{A}^{-1}) \\ &= \mathrm{Tr}(\mathbf{B} \mathbf{B}) \end{aligned}$

5: トレースの性質 $\mathrm{Tr}(\mathbf{A} + \mathbf{B}) = \mathrm{Tr}(\mathbf{A}) + \mathrm{Tr}(\mathbf{B})$ より、トレースをまとめます。
6: $\mathbf{B}$ を括り出します。

　 $\mathbf{X}, \mathbf{y}, \mathbf{W}$ などの二次形式の期待値をとります。

$\begin{align} \mathbb{E} \Bigl[ (\mathbf{y} - \mathbf{X} \boldsymbol{\beta})^{\top} \mathbf{W}^{\top} \mathbf{W} (\mathbf{y} - \mathbf{X} \boldsymbol{\beta}) \Bigr] &= \mathbb{E} \Bigl[ (\mathbf{y} - \mathbf{X} \boldsymbol{\beta})^{\top} \mathbf{A}^{\top} (\mathbf{A}^{-1})^{\top} \mathbf{W}^{\top} \mathbf{W} \mathbf{A}^{-1} \mathbf{A} (\mathbf{y} - \mathbf{X} \boldsymbol{\beta}) \Bigr] \\ &= \mathbb{E} \Bigl[ \boldsymbol{\epsilon}^{\top} \mathbf{B}^{\top} \mathbf{B} \boldsymbol{\epsilon} \Bigr] \\ &= \sigma^2 \mathrm{Tr}(\mathbf{B}^{\top} \mathbf{B}) \tag{13} \end{align}$

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

1: $\mathbf{I} = \mathbf{A}^{-1} \mathbf{A}$ を掛けます。
2: $\mathbf{A}$ と括弧の積を、式(4)で置き換えます。
2: $\mathbf{W} \mathbf{A}^{-1}$ を、式(2)で置き換えます。

　ただし、次のように変形しています。

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

　転置の性質 $(\mathbf{A} \mathbf{B})^{\top} = \mathbf{B}^{\top} \mathbf{A}^{\top}$ より、変形しています。

3: $\boldsymbol{\epsilon}, \mathbf{B}$ の二次形式の期待値に、式(14)を代入します。

　 $\boldsymbol{\epsilon}, \mathbf{B}$ の二次形式の期待値をとります。

$\begin{align} \mathbb{E}[ \boldsymbol{\epsilon}^{\top} \mathbf{B}^{\top} \mathbf{B} \boldsymbol{\epsilon} \Bigr] &= \mathbb{E} \Bigl[ \mathrm{Tr} ( \mathbf{B}^{\top} \mathbf{B} \boldsymbol{\epsilon} \boldsymbol{\epsilon}^{\top} ) \Bigr] \\ &= \mathrm{Tr} \Bigl( \mathbb{E} [ \mathbf{B}^{\top} \mathbf{B} \boldsymbol{\epsilon} \boldsymbol{\epsilon}^{\top} ] \Bigr) \\ &= \mathrm{Tr} \Bigl( \mathbf{B}^{\top} \mathbf{B} \mathbb{E}[\boldsymbol{\epsilon} \boldsymbol{\epsilon}^{\top}] \Bigr) \\ &= \mathrm{Tr} \Bigl( \mathbf{B}^{\top} \mathbf{B} \sigma^2 \mathbf{I} \Bigr) \\ &= \sigma^2 \mathrm{Tr}(\mathbf{B}^{\top} \mathbf{B}) \tag{14} \end{align}$

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

$\mathbf{B}^{\top} \mathbf{B}$ を1つの行列とみなして、式(11)と同様に変形します。

　以上で、フィッシャー情報行列(7)の3行3列目の要素の式が得られました。

フィッシャー情報行列

　対数尤度関数のフィッシャー情報行列を求めます。

　フィッシャー情報行列(7)の各要素をそれぞれの式(0.48)(0.49)(0.50)(0.51)(0.52)(0.53)と置き換えます。

$\displaystyle \mathbb{I}(\boldsymbol{\theta}) = \begin{pmatrix} - \frac{1}{\sigma^2} \mathbf{X}^{\top} \mathbf{A}^{\top} \mathbf{A} \mathbf{X} & \mathbf{0} & \mathbf{0} \\ \cdot & - \frac{N}{2} \frac{1}{(\sigma^2)^2} & - \frac{1}{\sigma^2} \mathrm{Tr}(\mathbf{B}) \\ \cdot & \cdot & - \mathrm{Tr} \Bigl( (\mathbf{B} + \mathbf{B}^{\top}) \mathbf{B} \Bigr) \end{pmatrix} \tag{7'}$

　フィッシャー情報行列の式が得られました。

漸近分散共分散行列

　漸近分散共分散行列を求めます。

　対数尤度関数 $\log L(\boldsymbol{\theta})$ のフィッシャー情報行列の逆行列を求めます。フィッシャー情報行列 $\mathbb{I}(\boldsymbol{\theta})$ の逆行列(ヘッセ行列 $\mathbf{H}(\boldsymbol{\theta})$ の負の期待値の逆行列)を漸近分散共分散行列 $\mathbb{Avar}(\boldsymbol{\theta})$ と呼びます。

$\begin{align} \mathbb{Avar}(\boldsymbol{\theta}) &\equiv \mathbb{I}(\boldsymbol{\theta})^{-1} \\ &= - \mathbb{E}[\mathbf{H}(\boldsymbol{\theta})]^{-1} \\ &= \begin{pmatrix} \frac{1}{\sigma^2} \mathbf{X}^{\top} \mathbf{A}^{\top} \mathbf{A} \mathbf{X} & \mathbf{0} & \mathbf{0} \\ \cdot & \frac{N}{2} \frac{1}{(\sigma^2)^2} & \frac{1}{\sigma^2} \mathrm{Tr}(\mathbf{B}) \\ \cdot & \cdot & \mathrm{Tr} \Bigl( (\mathbf{B} + \mathbf{B}^{\top}) \mathbf{B} \Bigr) \end{pmatrix}^{-1} \tag{15} \end{align}$

　以上で、対数尤度関数の漸近分散共分散行列の式が得られました。

　この記事では、SEMの漸近分散共分散行列を数式で確認しました。次の記事では最尤法を数式で確認します。

参考文献

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

作者:植杉大
森北出版

Amazon

おわりに

　SLMと比べて、ヘッセ行列の式はSEMの方が複雑でしたが、漸近分散共分散行列の式はSEMの方がシンプルになって面白かったです。未だに何に使うのか全然分からな行列ですが、やっといてよかったです。
　行列AまたはAの逆行列が、yの式でεとXβの両方に掛かるのか、εの式でyとXβの両方に掛かるのか、この2つの違いがSLMとSEMの異なる点であり、どの式が複雑になるのかに影響してるようですね。

　最後に、Juice=Juiceの楽曲をどうぞ♪

　夜明けに夢を見た日⚽

【次の内容】

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

からっぽのしょこ

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

0.2.2：空間誤差モデル(SEM)の漸近分散共分散行列の導出【はじめての地理空間DSのノート】

はじめに

0.2.2 空間誤差モデル(SEM)の漸近分散共分散行列の導出

モデルの確認

定義式

計算式

尤度関数

ヘッセ行列

漸近分散共分散行列の導出

フィッシャー情報行列の設定

回帰・回帰パラメータによる偏微分の期待値

回帰・分散パラメータによる偏微分の期待値

回帰・空間パラメータによる偏微分の期待値

分散・分散パラメータによる偏微分の期待値

分散・空間パラメータによる偏微分の期待値

空間・空間パラメータによる偏微分の期待値

フィッシャー情報行列

漸近分散共分散行列

参考文献

おわりに