はじめに

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

　この記事では、空間ラグモデルの漸近分散共分散行列について、数式を使ってを使って解説します。

【前の内容】

www.anarchive-beta.com

【他の内容】

www.anarchive-beta.com

【今回の内容】

はじめに
0.2.1 空間ラグモデル(SLM)の漸近分散共分散行列の導出
- モデルの確認
- 漸近分散共分散行列の導出
参考文献
おわりに

0.2.1 空間ラグモデル(SLM)の漸近分散共分散行列の導出

　空間ラグモデル(SLM・Spatial Lag Model・空間自己回帰モデル・SARモデル・Spatial Autoregressive Model)における漸近分散共分散行列(asymptotic variance-covariance matrix)を導出します。
　SLMの定義式については「【Python】0.1.2：空間ラグモデル(SLM)の定義式【はじめての地理空間DSのノート】 - からっぽのしょこ」を参照してください。

モデルの確認

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

定義式

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

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

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

　誤差項 $\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$ は分散パラメータ(スカラ)です。

計算式

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

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

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

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

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

尤度関数

　SLMの尤度関数については「【Python】0.2.1：空間ラグモデル(SLM)の尤度関数の導出【はじめての地理空間DSのノート】 - からっぽのしょこ」を参照してください。

　SLMのパラメータ $\boldsymbol{\beta}, \sigma^2, \rho$ をまとめて、パラメータベクトル $\boldsymbol{\theta} = (\boldsymbol{\beta}, \sigma^2, \rho)^{\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{4} \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.13} \end{align}$

ヘッセ行列

　SLMのヘッセ行列については「【Python】0.2.1：空間ラグモデル(SLM)のヘッセ行列の導出【はじめての地理空間DSのノート】 - からっぽのしょこ」を参照してください。

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

$\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 \rho} \\ \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 \rho} \\ \frac{\partial^2 \log L(\boldsymbol{\theta})}{\partial \rho \partial \boldsymbol{\beta}^{\top}} & \frac{\partial^2 \log L(\boldsymbol{\theta})}{\partial \rho \partial \sigma^2} & \frac{\partial^2 \log L(\boldsymbol{\theta})}{\partial \rho \partial \rho} \end{pmatrix} \\ &= \begin{pmatrix} - \frac{1}{\sigma^2} \mathbf{X}^{\top} \mathbf{X} & - \frac{1}{(\sigma^2)^2} \mathbf{X}^{\top} \boldsymbol{\epsilon} & - \frac{1}{\sigma^2} \mathbf{X}^{\top} \mathbf{W} \mathbf{y} \\ \cdot & \frac{N}{2 (\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} \\ \cdot & \cdot & - \mathrm{Tr}(\mathbf{A}^{-1} \mathbf{W} \mathbf{A}^{-1} \mathbf{W}) - \frac{1}{\sigma^2} \mathbf{y}^{\top} \mathbf{W}^{\top} \mathbf{W} \mathbf{y} \end{pmatrix} \tag{5} \end{align}$

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

スポンサードリンク

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

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

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

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

　対数尤度関数 $\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 \rho} \\ \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 \rho} \\ \frac{\partial^2 \log L(\boldsymbol{\theta})}{\partial \rho \partial \boldsymbol{\beta}^{\top}} & \frac{\partial^2 \log L(\boldsymbol{\theta})}{\partial \rho \partial \sigma^2} & \frac{\partial^2 \log L(\boldsymbol{\theta})}{\partial \rho \partial \rho} \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 \rho}] \\ \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 \rho}] \\ \mathbb{E}[\frac{\partial^2 \log L(\boldsymbol{\theta})}{\partial \rho \partial \boldsymbol{\beta}^{\top}}] & \mathbb{E}[\frac{\partial^2 \log L(\boldsymbol{\theta})}{\partial \rho \partial \sigma^2}] & \mathbb{E}[\frac{\partial^2 \log L(\boldsymbol{\theta})}{\partial \rho \partial \rho}] \end{pmatrix} \tag{6} \end{align}$

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

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

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

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

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

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

　フィッシャー情報行列(6)の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{X} \Biggr] \tag{0.24'}\\ &= - \frac{1}{\sigma^2} \mathbf{X}^{\top} \mathbf{X} \tag{0.30} \end{align}$

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

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

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

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

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

　フィッシャー情報行列(6)の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} \boldsymbol{\epsilon} \Biggr] \tag{0.25'}\\ &= - \frac{1}{(\sigma^2)^2} \mathbf{X}^{\top} \mathbb{E}[\boldsymbol{\epsilon}] \\ &= - \frac{1}{(\sigma^2)^2} \mathbf{X}^{\top} \mathbf{0} \\ &= \mathbf{0} \tag{0.31} \end{align}$

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

　フィッシャー情報行列(6)の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.27'}\\ &= \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.33} \end{align}$

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

1: 式(0.27)より、式(6)の $(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}$ の内積の期待値に、式(9)を代入します。

　 $\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{9} \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}$ の積の期待値に、式(10)を代入します。
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{10} \end{align}$

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

1: ベクトルの要素を明示します。
2: ベクトルの積を計算します。
3: 行列の期待値を、期待値の行列に変形します。
4: SLMの定義(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列目の要素の式が得られました。

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

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

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

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

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

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

　 $\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{11} \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}$ の積の期待値に、式(10)を代入します。
5: トレースの性質 $\mathrm{Tr}(a \mathbf{X}) = a \mathrm{Tr}(\mathbf{X})$ より、係数をトレースの外に出します。

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

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

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

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

$\begin{align} \mathbb{E} \Biggl[ \frac{\partial \log L(\boldsymbol{\theta})}{\partial \rho \partial \rho} \Biggr] &= \mathbb{E} \Biggl[ - \mathrm{Tr}(\mathbf{A}^{-1} \mathbf{W} \mathbf{A}^{-1} \mathbf{W}) - \frac{1}{\sigma^2} \mathbf{y}^{\top} \mathbf{W}^{\top} \mathbf{W} \mathbf{y} \Biggr] \tag{0.29'}\\ &= - \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}^{\top} \mathbf{W}^{\top} \mathbf{W} \mathbf{y} \Bigr] \\ &= - \mathrm{Tr}(\mathbf{A}^{-1} \mathbf{W} \mathbf{A}^{-1} \mathbf{W}) - \frac{1}{\sigma^2} \Bigl( (\mathbf{B} \mathbf{X} \boldsymbol{\beta})^{\top} \mathbf{B} \mathbf{X} \boldsymbol{\beta} + \sigma^2 \mathrm{Tr}(\mathbf{B}^{\top} \mathbf{B}) \Bigr) \\ &= - \mathrm{Tr}(\mathbf{B} \mathbf{B}) - \mathrm{Tr}(\mathbf{B}^{\top} \mathbf{B}) - \frac{1}{\sigma^2} (\mathbf{B} \mathbf{X} \boldsymbol{\beta})^{\top} \mathbf{B} \mathbf{X} \boldsymbol{\beta} \\ &= - \mathrm{Tr} \Bigl( \mathbf{B} \mathbf{B} + \mathbf{B}^{\top} \mathbf{B} \Bigr) - \frac{1}{\sigma^2} (\mathbf{B} \mathbf{X} \boldsymbol{\beta})^{\top} \mathbf{B} \mathbf{X} \boldsymbol{\beta} \\ &= - \mathrm{Tr} \Bigl( (\mathbf{B} + \mathbf{B}^{\top}) \mathbf{B} \Bigr) - \frac{1}{\sigma^2} (\mathbf{B} \mathbf{X} \boldsymbol{\beta})^{\top} \mathbf{B} \mathbf{X} \boldsymbol{\beta} \tag{0.35} \end{align}$

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

1: 式(0.29)より、式(6)の $(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{y}, \mathbf{W}$ の二次形式の期待値に、式(12)を代入します。
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}$

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

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

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

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

1: $\mathbf{y}$ に、式(3)を代入します。
2: $\mathbf{W} \mathbf{A}^{-1}$ を、式(2)で置き換えます。

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

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

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: 中2つの項は二次形式(スカラ)なので、転置できます。

　2つの項が一致するのでまとめます。

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

5: 期待値の性質 $\mathbb{E}[x + y] = \mathbb{E}[x] + \mathbb{E}[y]$ より、期待値の和に分割します。
5: 期待値の性質 $\mathbb{E}[a x] = a \mathbb{E}[x]$ より、係数を期待値の外に出します。
6: 期待値の性質 $\mathbb{E}[a] = a$ より、期待値を外します。
6: $\boldsymbol{\epsilon}$ の期待値に、式(7)を代入します。
6: $\boldsymbol{\epsilon}, \mathbf{B}$ の二次形式の期待値に、式(13)を代入します。

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

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

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

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

フィッシャー情報行列

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

　フィッシャー情報行列(6)の各要素をそれぞれの式(0.30)(0.31)(0.32)(0.33)(0.34)(0.35)と置き換えます。

$\displaystyle \mathbb{I}(\boldsymbol{\theta}) = \begin{pmatrix} - \frac{1}{\sigma^2} \mathbf{X}^{\top} \mathbf{X} & \mathbf{0} & - \frac{1}{\sigma^2} \mathbf{X}^{\top} \mathbf{B} \mathbf{X} \boldsymbol{\beta} \\ \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) - \frac{1}{\sigma^2} (\mathbf{B} \mathbf{X} \boldsymbol{\beta})^{\top} \mathbf{B} \mathbf{X} \boldsymbol{\beta} \end{pmatrix} \tag{6'}$

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

漸近分散共分散行列

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

　対数尤度関数 $\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{X} & \mathbf{0} & \frac{1}{\sigma^2} \mathbf{X}^{\top} \mathbf{B} \mathbf{X} \boldsymbol{\beta} \\ \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) + \frac{1}{\sigma^2} (\mathbf{B} \mathbf{X} \boldsymbol{\beta})^{\top} \mathbf{B} \mathbf{X} \boldsymbol{\beta} \end{pmatrix}^{-1} \tag{14} \end{align}$

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

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

参考文献

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

作者:植杉大
森北出版

Amazon

おわりに

　漸近分散共分散行列とは何なのかよく分からないまま最後まで式を追いました。最尤法を実装することになれば使うかもしれませんが、最尤法の導出には使わないようですね。最尤法の方を先に書いたのですが、そこで最後まで登場しなかったことで気付きました。
　どちらかというとフィッシャー情報行列を求めたという気持ちです。その逆行列ですって最後に投げてる気分です。これは、書く内容の目次(小見出し)を作ってるときに気付きました。
　では、次がSLMの導出編のラストです。

　最後に、いぎなり東北産のライブ映像をどうぞ♪

【次の内容】

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

www.anarchive-beta.com

　SEMの漸近分散共分散行列を数式で確認します。

https://www.anarchive-beta.com/entry/2026/06/30/180000www.anarchive-beta.com

からっぽのしょこ

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

0.2.1：空間ラグモデル(SLM)の漸近分散共分散行列の導出【はじめての地理空間DSのノート】

はじめに

0.2.1 空間ラグモデル(SLM)の漸近分散共分散行列の導出

モデルの確認

定義式

計算式

尤度関数

ヘッセ行列

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

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

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

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

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

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

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

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

フィッシャー情報行列

漸近分散共分散行列

参考文献

おわりに