はじめに

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

　この記事では、空間誤差モデルのヘッセ行列について、数式を使って解説します。

【前の内容】

www.anarchive-beta.com

【他の内容】

www.anarchive-beta.com

【今回の内容】

はじめに
0.2.2 空間誤差モデル(SEM)のヘッセ行列の導出
参考文献
おわりに

0.2.2 空間誤差モデル(SEM)のヘッセ行列の導出

　空間誤差モデル(SEM・Spatial Error Model)におけるヘッセ行列(Hessian 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の定義式」を参照してください。

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

$\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}$

　独立誤差項 $\boldsymbol{\epsilon}$ の内積(2乗和)は、次の式となります。

$\begin{align} \boldsymbol{\epsilon}^{\top} \boldsymbol{\epsilon} &= \Bigl( \mathbf{A} (\mathbf{y} - \mathbf{X} \boldsymbol{\beta}) \Bigr)^{\top} \Bigl( \mathbf{A} (\mathbf{y} - \mathbf{X} \boldsymbol{\beta}) \Bigr) \\ &= (\mathbf{y} - \mathbf{X} \boldsymbol{\beta})^{\top} \mathbf{A}^{\top} \mathbf{A} (\mathbf{y} - \mathbf{X} \boldsymbol{\beta}) \\ &= \mathbf{y}^{\top} \mathbf{A}^{\top} \mathbf{A} \mathbf{y} - 2 \boldsymbol{\beta}^{\top} \mathbf{X}^{\top} \mathbf{A}^{\top} \mathbf{A} \mathbf{y} + \boldsymbol{\beta}^{\top} \mathbf{X}^{\top} \mathbf{A}^{\top} \mathbf{A} \mathbf{X} \boldsymbol{\beta} \tag{3} \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{4} \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}$

　以上の式を用いてヘッセ行列を求めます。

スポンサードリンク

対数尤度関数の1階微分の導出

　次は、SEMの対数尤度関数の1階微分を導出します。

対数尤度の微分の設定

　対数尤度関数の微分を確認します。

　対数尤度関数 $\log L(\boldsymbol{\theta})$ の1階微分を求めます。

$\displaystyle \frac{\partial \log L(\boldsymbol{\theta})}{\partial \boldsymbol{\theta}} = \begin{pmatrix} \frac{\partial \log L(\boldsymbol{\theta})}{\partial \boldsymbol{\beta}} \\ \frac{\partial \log L(\boldsymbol{\theta})}{\partial \sigma^2} \\ \frac{\partial \log L(\boldsymbol{\theta})}{\partial \lambda} \end{pmatrix} \tag{5}$

　各要素は、各パラメータ $\boldsymbol{\beta}, \sigma^2, \lambda$ による偏微分で求まるのが分かります。

　対数尤度関数の1階微分(5)の各要素の式を求めていきます。

回帰パラメータによる偏微分

　各種の変数の回帰パラメータに関する微分を求めます。

　対数尤度関数の1階微分(5)の1番目の要素は、対数尤度関数(0.37)を回帰パラメータ $\boldsymbol{\beta}$ に関して微分して求まります。

$\begin{align} \frac{\partial \log L(\boldsymbol{\theta})}{\partial \boldsymbol{\beta}} &= \frac{\partial}{\partial \boldsymbol{\beta}} \Biggl\{ \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} \Biggr\} \\ &= \frac{\partial}{\partial \boldsymbol{\beta}} \Bigl\{ \log |\mathbf{A}| \Bigr\} + \frac{\partial}{\partial \boldsymbol{\beta}} \Biggl\{ - \frac{N \log (2 \pi)}{2} \Biggr\} + \frac{\partial}{\partial \boldsymbol{\beta}} \Biggl\{ - \frac{N}{2} \log \sigma^2 \Biggr\} - \frac{1}{2} \frac{1}{\sigma^2} \frac{ \partial \boldsymbol{\epsilon}^{\top} \boldsymbol{\epsilon} }{ \partial \boldsymbol{\beta} } \\ &= 0 + 0 + 0 - \frac{1}{2} \frac{1}{\sigma^2} \Bigl( - 2 \mathbf{X}^{\top} \mathbf{A}^{\top} \boldsymbol{\epsilon} \Bigr) \\ &= \frac{1}{\sigma^2} \mathbf{X}^{\top} \mathbf{A}^{\top} \boldsymbol{\epsilon} \tag{6} \end{align}$

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

1: 式(0.37)の偏微分の式を立てます。 $\boldsymbol{\beta}$ に関する微分なので、 $\boldsymbol{\beta}$ 以外の項は定数として扱います。
2: 和の微分 $\frac{d \{f(x) + g(x)\}}{d x} = \frac{d f(x)}{d x} + \frac{d g(x)}{d x}$ より、微分の和に分割します。
3: $\boldsymbol{\epsilon}$ の内積の微分に、式(7)を代入します。

　独立誤差項 $\boldsymbol{\epsilon}$ の内積を $\boldsymbol{\beta}$ に関して微分します。

$\begin{align} \frac{\partial \boldsymbol{\epsilon}^{\top} \boldsymbol{\epsilon}}{\partial \boldsymbol{\beta}} &= \frac{\partial}{\partial \boldsymbol{\beta}} \Bigl\{ \Bigl( \mathbf{A} (\mathbf{y} - \mathbf{X} \boldsymbol{\beta}) \Bigr)^{\top} \Bigl( \mathbf{A} (\mathbf{y} - \mathbf{X} \boldsymbol{\beta}) \Bigr) \Bigr\} \\ &= \frac{\partial}{\partial \boldsymbol{\beta}} \Bigl\{ \mathbf{y}^{\top} \mathbf{A}^{\top} \mathbf{A} \mathbf{y} - 2 \boldsymbol{\beta}^{\top} \mathbf{X}^{\top} \mathbf{A}^{\top} \mathbf{A} \mathbf{y} + \boldsymbol{\beta}^{\top} \mathbf{X}^{\top} \mathbf{A}^{\top} \mathbf{A} \mathbf{X} \boldsymbol{\beta} \Bigr\} \\ &= \frac{\partial}{\partial \boldsymbol{\beta}} \Bigl\{ \mathbf{y}^{\top} \mathbf{A}^{\top} \mathbf{A} \mathbf{y} \Bigr\} - 2 \frac{\partial \boldsymbol{\beta}^{\top}}{\partial \boldsymbol{\beta}} \mathbf{X}^{\top} \mathbf{A}^{\top} \mathbf{A} \mathbf{y} + \frac{\partial}{\partial \boldsymbol{\beta}} \Bigl\{ \boldsymbol{\beta}^{\top} \mathbf{X}^{\top} \mathbf{A}^{\top} \mathbf{A} \mathbf{X} \boldsymbol{\beta} \Bigr\} \\ &= 0 - 2 \mathbf{I} \mathbf{X}^{\top} \mathbf{A}^{\top} \mathbf{A} \mathbf{y} + 2 \mathbf{X}^{\top} \mathbf{A}^{\top} \mathbf{A} \mathbf{X} \boldsymbol{\beta} \\ &= - 2 \mathbf{X}^{\top} \mathbf{A}^{\top} \mathbf{A} \mathbf{y} + 2 \mathbf{X}^{\top} \mathbf{A}^{\top} \mathbf{A} \mathbf{X} \boldsymbol{\beta} \\ &= - 2 \mathbf{X}^{\top} \mathbf{A}^{\top} \mathbf{A} (\mathbf{y} - \mathbf{X} \boldsymbol{\beta}) \\ &= - 2 \mathbf{X}^{\top} \mathbf{A}^{\top} \boldsymbol{\epsilon} \tag{7} \end{align}$

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

1-2: 式(3)の偏微分の式を立てます。
3: 和の微分 $\frac{d \{f(x) + g(x)\}}{d x} = \frac{d f(x)}{d x} + \frac{d g(x)}{d x}$ より、微分の和に分割します。
4: $\mathbf{X}, \boldsymbol{\beta}$ などの二次形式の微分に、式(8)を代入します。
6: 共通の項と $-1$ を括り出します。
7: $\mathbf{A}$ と括弧の積を式(2)で置き換えます。

　説明変数・回帰パラメータ $\mathbf{X}, \boldsymbol{\beta}$ などの二次形式を $\boldsymbol{\beta}$ に関して微分します。

$\begin{align} \frac{\partial}{\partial \boldsymbol{\beta}} \Bigl\{ \boldsymbol{\beta}^{\top} \mathbf{X}^{\top} \mathbf{A}^{\top} \mathbf{A} \mathbf{X} \boldsymbol{\beta} \Bigr\} &= \Bigl( \mathbf{X}^{\top} \mathbf{A}^{\top} \mathbf{A} \mathbf{X} + ( \mathbf{X}^{\top} \mathbf{A}^{\top} \mathbf{A} \mathbf{X} )^{\top} \Bigr) \boldsymbol{\beta} \\ &= \Bigl( \mathbf{X}^{\top} \mathbf{A}^{\top} \mathbf{A} \mathbf{X} + \mathbf{X}^{\top} \mathbf{A}^{\top} \mathbf{A} \mathbf{X} \Bigr) \boldsymbol{\beta} \\ &= 2 \mathbf{X}^{\top} \mathbf{A}^{\top} \mathbf{A} \mathbf{X} \boldsymbol{\beta} \tag{8} \end{align}$

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

1: 二次形式の微分 $\frac{\partial \mathbf{x}^{\top} \mathbf{A} \mathbf{x}}{\partial \mathbf{x}} = (\mathbf{A} + \mathbf{A}^{\top}) \mathbf{x}$ より、 $\mathbf{X}^{\top} \mathbf{A}^{\top} \mathbf{A} \mathbf{X}$ を1つの行列とみなして微分を計算します。
2: 転置の性質 $(\mathbf{A} \mathbf{B})^{\top} = \mathbf{B}^{\top} \mathbf{A}^{\top}$ より、括弧を展開します。

　以上で、対数尤度関数の1階微分(5)の1番目の要素の式が得られました。

分散パラメータによる偏微分

　各種の変数の分散パラメータに関する微分を求めます。

　対数尤度関数の1階微分(5)の2番目の要素は、対数尤度関数(0.37)を分散パラメータ $\sigma^2$ に関して微分して求まります。

$\begin{align} \frac{\partial \log L(\boldsymbol{\theta})}{\partial \sigma^2} &= \frac{\partial}{\partial \sigma^2} \Biggl\{ \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} \Biggr\} \\ &= \frac{\partial}{\partial \sigma^2} \Bigl\{ \log |\mathbf{A}| \Bigr\} + \frac{\partial}{\partial \sigma^2} \Biggl\{ - \frac{N \log (2 \pi)}{2} \Biggr\} - \frac{N}{2} \frac{\partial \log \sigma^2}{\partial \sigma^2} - \frac{1}{2} \frac{\partial (\sigma^2)^{-1}}{\partial \sigma^2} \boldsymbol{\epsilon}^{\top} \boldsymbol{\epsilon} \\ &= 0 + 0 - \frac{N}{2} \frac{1}{\sigma^2} - \frac{1}{2} (- 1) (\sigma^2)^{-2} \boldsymbol{\epsilon}^{\top} \boldsymbol{\epsilon} \\ &= - \frac{N}{2} \frac{1}{\sigma^2} + \frac{1}{2} \frac{1}{(\sigma^2)^2} \boldsymbol{\epsilon}^{\top} \boldsymbol{\epsilon} \tag{9} \end{align}$

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

1: 式(0.37)の偏微分の式を立てます。 $\sigma^2$ に関する微分なので、 $\sigma^2$ 以外の項は定数として扱います。
2: 和の微分 $\frac{d \{f(x) + g(x)\}}{d x} = \frac{d f(x)}{d x} + \frac{d g(x)}{d x}$ より、微分の和に分割します。
2: 微分の計算が分かりやすいように、分数の指数表記 $\frac{1}{x} = x^{-1}$ に書き換えています。
3: 自然対数の微分 $\frac{d \log x}{d x} = \frac{1}{x}$ より、微分を計算します。
3: べき乗の微分 $\frac{d x^n}{d x} = n x^{n-1}$ より、 $\sigma^2$ を1つの変数とみなして微分を計算します。
4: 分数の指数表記 $\frac{1}{x^n} = x^{-n}$ を戻します。

　以上で、対数尤度関数の1階微分(5)の2番目の要素の式が得られました。

空間パラメータによる偏微分

　対数尤度関数の1階微分(5)の3番目の要素は、対数尤度関数(0.37)を空間パラメータ $\lambda$ に関して微分して求まります。

$\begin{align} \frac{\partial \log L(\boldsymbol{\theta})}{\partial \lambda} &= \frac{\partial}{\partial \lambda} \Biggl\{ \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} \Biggr\} \\ &= \frac{\partial \log |\mathbf{A}|}{\partial \lambda} + \frac{\partial}{\partial \lambda} \Biggl\{ - \frac{N \log (2 \pi)}{2} \Biggr\} + \frac{\partial}{\partial \lambda} \Biggl\{ - \frac{N}{2} \log \sigma^2 \Biggr\} - \frac{1}{2} \frac{1}{\sigma^2} \frac{\partial \boldsymbol{\epsilon}^{\top} \boldsymbol{\epsilon}}{\partial \lambda} \\ &= - \mathrm{Tr}(\mathbf{A}^{-1} \mathbf{W}) + 0 + 0 - \frac{1}{2} \frac{1}{\sigma^2} \Bigl( - 2 \boldsymbol{\epsilon}^{\top} \mathbf{W} (\mathbf{y} - \mathbf{X} \boldsymbol{\beta}) \Bigr) \\ &= - \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}$

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

1: 式(0.37)の偏微分の式を立てます。 $\lambda$ に関する微分なので、 $\lambda$ 以外の項は定数として扱います。
2: 和の微分 $\frac{d \{f(x) + g(x)\}}{d x} = \frac{d f(x)}{d x} + \frac{d g(x)}{d x}$ より、微分の和に分割します。
3: $\mathbf{A}$ の対数行列式の微分に、式(0.15)を代入します。
3: $\boldsymbol{\epsilon}$ の内積の微分に、式(11)を代入します。

　 $\mathbf{A}$ の対数行列式を $\lambda$ に関して微分します。

$\begin{align} \frac{\partial \log |\mathbf{A}|}{\partial \lambda} &= |\mathbf{A}|^{-1} \frac{\partial |\mathbf{A}|}{\partial \lambda} \\ &= |\mathbf{A}|^{-1} |\mathbf{A}| \mathrm{Tr} \Bigl( \mathbf{A}^{-1} \frac{\partial \mathbf{A}}{\partial \lambda} \Bigr) \\ &= \mathrm{Tr} \Bigl( \mathbf{A}^{-1} (- \mathbf{W}) \Bigr) \\ &= - \mathrm{Tr}(\mathbf{A}^{-1} \mathbf{W}) \tag{0.15} \end{align}$

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

1: 合成関数の微分(連鎖率) $\frac{d y}{d x} = \frac{d y}{d u} \frac{d u}{d x}$ 、自然対数の微分 $\frac{d \log x}{d x} = x^{-1}$ より、 $|\mathbf{A}|$ を中間変数とみなして微分を計算します。
2: $x$ に依存する行列 $\mathbf{A}$ における行列式の微分(ヤコビの公式) $\frac{\partial |\mathbf{A}|}{\partial x} = |\mathbf{A}| \mathrm{Tr}(\mathbf{A}^{-1} \frac{\partial \mathbf{A}}{\partial x})$ より、トレースを用いた微分に変形します。
3: $\mathbf{A}$ の微分に、式(10)を代入します。
4: トレースの性質 $\mathrm{Tr}(- \mathbf{X}) = - \mathrm{Tr}(\mathbf{X})$ より、符号をトレースの外に出します。

　 $\mathbf{A}$ を $\lambda$ に関して微分します。

$\begin{align} \frac{\partial \mathbf{A}}{\partial \lambda} &= \frac{\partial}{\partial \lambda} \Bigl\{ \mathbf{I} - \lambda \mathbf{W} \Bigr\} \\ &= \frac{\partial}{\partial \lambda} \Bigl\{ \mathbf{I} \Bigr\} - \frac{\partial \lambda}{\partial \lambda} \mathbf{W} \\ &= \mathbf{0} - 1 \mathbf{W} \\ &= - \mathbf{W} \tag{10} \end{align}$

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

1: 式(1)の偏微分の式を立てます。
2: 和の微分 $\frac{d \{f(x) + g(x)\}}{d x} = \frac{d f(x)}{d x} + \frac{d g(x)}{d x}$ より、微分の和に分割します。

　独立誤差項 $\boldsymbol{\epsilon}$ の内積を $\lambda$ に関して微分します。

$\begin{align} \frac{\partial \boldsymbol{\epsilon}^{\top} \boldsymbol{\epsilon}}{\partial \lambda} &= \frac{\partial \boldsymbol{\epsilon}^{\top}}{\partial \lambda} \boldsymbol{\epsilon} + \boldsymbol{\epsilon}^{\top} \frac{\partial \boldsymbol{\epsilon}}{\partial \lambda} \\ &= \Bigl( - \mathbf{W} (\mathbf{y} - \mathbf{X} \boldsymbol{\beta}) \Bigr)^{\top} \boldsymbol{\epsilon} + \boldsymbol{\epsilon}^{\top} \Bigl( - \mathbf{W} (\mathbf{y} - \mathbf{X} \boldsymbol{\beta}) \Bigr) \\ &= - (\mathbf{y} - \mathbf{X} \boldsymbol{\beta})^{\top} \mathbf{W}^{\top} \boldsymbol{\epsilon} - \boldsymbol{\epsilon}^{\top} \mathbf{W} (\mathbf{y} - \mathbf{X} \boldsymbol{\beta}) \\ &= - 2 \boldsymbol{\epsilon}^{\top} \mathbf{W} (\mathbf{y} - \mathbf{X} \boldsymbol{\beta}) \tag{11} \end{align}$

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

1: 積の微分 $\frac{d \{f(x) g(x)\}}{d x} = \frac{d f(x)}{d x} g(x) + f(x) \frac{d g(x)}{d x}$ より、2つの微分の和に変形します。
2: $\boldsymbol{\epsilon}$ の微分に、式(12)を代入します。
3: 転置の性質 $(\mathbf{A} \mathbf{B})^{\top} = \mathbf{B}^{\top} \mathbf{A}^{\top}$ より、括弧を展開します。
4: 二次形式はスカラなので、転置できます。

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

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

　独立誤差項 $\boldsymbol{\epsilon}$ を $\lambda$ に関して微分します。

$\begin{align} \frac{\partial \boldsymbol{\epsilon}}{\partial \lambda} &= \frac{\partial}{\partial \lambda} \Bigl\{ \mathbf{A} (\mathbf{y} - \mathbf{X} \boldsymbol{\beta}) \Bigr\} \\ &= \frac{\partial}{\partial \lambda} \Bigl\{ (\mathbf{I} - \lambda \mathbf{W}) (\mathbf{y} - \mathbf{X} \boldsymbol{\beta}) \Bigr\} \\ &= \frac{\partial}{\partial \lambda} \Bigl\{ \mathbf{y} - \lambda \mathbf{W} \mathbf{y} - \mathbf{X} \boldsymbol{\beta} + \lambda \mathbf{W} \mathbf{X} \boldsymbol{\beta} \Bigr\} \\ &= \frac{\partial}{\partial \lambda} \Bigl\{ \mathbf{y} \Bigr\} - \frac{\partial \lambda}{\partial \lambda} \mathbf{W} \mathbf{y} + \frac{\partial}{\partial \lambda} \Bigl\{ - \mathbf{X} \boldsymbol{\beta} \Bigr\} + \frac{\partial \lambda}{\partial \lambda} \mathbf{W} \mathbf{X} \boldsymbol{\beta} \\ &= \mathbf{0} - 1 \mathbf{W} \mathbf{y} + \mathbf{0} + 1 \mathbf{W} \mathbf{X} \boldsymbol{\beta} \\ &= - \mathbf{W} \mathbf{y} + \mathbf{W} \mathbf{X} \boldsymbol{\beta} \\ &= - \mathbf{W} (\mathbf{y} - \mathbf{X} \boldsymbol{\beta}) \tag{12} \end{align}$

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

1: 式(2)の偏微分の式を立てます。
2: $\mathbf{A}$ に、式(1)を代入します。
3: 括弧を展開します。
4: 和の微分 $\frac{d \{f(x) + g(x)\}}{d x} = \frac{d f(x)}{d x} + \frac{d g(x)}{d x}$ より、微分の和に分割します。
7: $\mathbf{W}$ を括り出します。

　以上で、対数尤度関数の1階微分(5)の3番目の要素の式が得られました。

対数尤度の微分

　対数尤度関数の微分を求めます。

　対数尤度関数の1階微分(5)の各要素をそれぞれの式(6)(9)(0.40)と置き換えます。

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

　以上で、対数尤度関数の1階微分の式が得られました。

スポンサードリンク

対数尤度関数の2階微分の導出

　続いて、SEMの対数尤度関数の2階微分を導出します。

ヘッセ行列の設定

　対数尤度関数のヘッセ行列を確認します。

　対数尤度関数 $\log L(\boldsymbol{\theta})$ の2階微分を求めます。関数の2階偏微分を並べた行列をヘッセ行列 $\mathbf{H}(\boldsymbol{\theta})$ と呼びます。

$\displaystyle \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} \equiv \mathbf{H}(\boldsymbol{\theta}) \tag{13}$

　各要素は、各パラメータ $\boldsymbol{\beta}, \sigma^2, \lambda$ の組み合わせによる偏微分で求まるのが分かります。

　ヘッセ行列(13)の各要素の式を求めていきます。

回帰・回帰パラメータによる偏微分

　各種の変数の回帰パラメータに関する2階微分を求めます。

　ヘッセ行列(13)の1行1列目の要素は、対数尤度関数(0.37)を空間パラメータ $\boldsymbol{\beta}$ に関して2階微分して求まります。

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

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

1: 式(0.37)の偏微分の式を立てます。 $\boldsymbol{\beta}$ に関する微分なので、 $\boldsymbol{\beta}$ 以外の項は定数として扱います。
2: $\boldsymbol{\beta}, \boldsymbol{\beta}^{\top}$ による微分の順番を入れ換えます。
3: $\log L(\boldsymbol{\theta})$ の微分に、式(6)を代入します。
4: $\boldsymbol{\beta}$ と無関係な項を微分の外に出します。
5: $\boldsymbol{\beta}$ の微分に、式(14)を代入します。

　独立誤差項 $\boldsymbol{\epsilon}$ を $\boldsymbol{\beta}$ の転置に関して微分します。

$\begin{align} \frac{\partial \boldsymbol{\epsilon}}{\partial \boldsymbol{\beta}^{\top}} &= \frac{\partial}{\partial \boldsymbol{\beta}^{\top}} \Big\{ \mathbf{A} (\mathbf{y} - \mathbf{X} \boldsymbol{\beta}) \Bigr\} \\ &= \frac{\partial}{\partial \boldsymbol{\beta}^{\top}} \Big\{ \mathbf{A} \mathbf{y} - \mathbf{A} \mathbf{X} \boldsymbol{\beta} \Bigr\} \\ &= \frac{\partial}{\partial \boldsymbol{\beta}^{\top}} \Big\{ \mathbf{A} \mathbf{y} \Bigr\} - \mathbf{A} \mathbf{X} \frac{\partial \boldsymbol{\beta}}{\partial \boldsymbol{\beta}^{\top}} \\ &= \mathbf{0} - \mathbf{A} \mathbf{X} \mathbf{I} \\ &= - \mathbf{A} \mathbf{X} \tag{14} \end{align}$

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

1: 式(2)の偏微分の式を立てます。
2: 括弧を展開します。
3: 和の微分 $\frac{d \{f(x) + g(x)\}}{d x} = \frac{d f(x)}{d x} + \frac{d g(x)}{d x}$ より、微分の和に分割します。

　以上で、ヘッセ行列(13)の1行1列目の要素の式が得られました。

回帰・分散パラメータによる偏微分

　各種の変数の回帰パラメータと分散パラメータに関する微分を求めます。

　ヘッセ行列(13)の1行2列目の要素は、対数尤度関数(0.37)を回帰パラメータ $\boldsymbol{\beta}$ と分散パラメータ $\sigma^2$ に関して微分して求まります。

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

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

1: 式(0.37)の偏微分の式を立てます。 $\boldsymbol{\beta}, \sigma^2$ に関する微分なので、それぞれ他の項は定数として扱います。
2: $\log L(\boldsymbol{\theta})$ の微分に、式(9)を代入します。
3: 和の微分 $\frac{d \{f(x) + g(x)\}}{d x} = \frac{d f(x)}{d x} + \frac{d g(x)}{d x}$ より、微分の和に分割します。
4: $\boldsymbol{\epsilon}$ の内積の微分に、式(7)を代入します。

　以上で、ヘッセ行列(13)の1行2列・2行1列目の要素の式が得られました。

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

　各種の変数の回帰パラメータと空間パラメータに関する微分を求めます。

　ヘッセ行列(13)の1行3列目の要素は、対数尤度関数(0.37)を回帰パラメータ $\boldsymbol{\beta}$ と空間パラメータ $\lambda$ に関して微分して求まります。

$\begin{align} \frac{\partial^2 \log L(\boldsymbol{\theta})}{\partial \boldsymbol{\beta} \partial \lambda} &= \frac{\partial}{\partial \boldsymbol{\beta}} \frac{\partial}{\partial \lambda} \log L(\boldsymbol{\theta}) \\ &= \frac{\partial}{\partial \boldsymbol{\beta}} \Biggl\{ - \mathrm{Tr}(\mathbf{A}^{-1} \mathbf{W}) + \frac{1}{\sigma^2} \boldsymbol{\epsilon}^{\top} \mathbf{W} (\mathbf{y} - \mathbf{X} \boldsymbol{\beta}) \Biggr\} \\ &= \frac{\partial}{\partial \boldsymbol{\beta}} \Biggl\{ - \mathrm{Tr}(\mathbf{A}^{-1} \mathbf{W}) + \frac{1}{\sigma^2} \Bigl( \boldsymbol{\epsilon}^{\top} \mathbf{W} \mathbf{y} - \boldsymbol{\epsilon}^{\top} \mathbf{W} \mathbf{X} \boldsymbol{\beta} \Bigr) \Biggr\} \\ &= \frac{\partial}{\partial \boldsymbol{\beta}} \Biggl\{ - \mathrm{Tr}(\mathbf{A}^{-1} \mathbf{W}) + \frac{1}{\sigma^2} \Biggl( \boldsymbol{\epsilon}^{\top} \mathbf{W} \mathbf{y} - \Bigl( \boldsymbol{\beta}^{\top} \mathbf{X}^{\top} \mathbf{W}^{\top} \mathbf{A} \mathbf{y} - \boldsymbol{\beta}^{\top} \mathbf{X}^{\top} \mathbf{A}^{\top} \mathbf{W} \mathbf{X} \boldsymbol{\beta} \Bigr) \Biggr) \Biggr\} \\ &= \frac{\partial}{\partial \boldsymbol{\beta}} \Bigr\{ - \mathrm{Tr}(\mathbf{A}^{-1} \mathbf{W}) \Bigr\} \\ &\qquad + \frac{1}{\sigma^2} \Biggl( \frac{\partial \boldsymbol{\epsilon}^{\top}}{\partial \boldsymbol{\beta}} \mathbf{W} \mathbf{y} - \frac{\partial \boldsymbol{\beta}^{\top}}{\partial \boldsymbol{\beta}} \mathbf{X}^{\top} \mathbf{W}^{\top} \mathbf{A} \mathbf{y} + \frac{\partial}{\partial \boldsymbol{\beta}} \Bigl\{ \boldsymbol{\beta}^{\top} \mathbf{X}^{\top} \mathbf{A}^{\top} \mathbf{W} \mathbf{X} \boldsymbol{\beta} \Bigr\} \Biggr) \\ &= 0 + \frac{1}{\sigma^2} \Bigl( (- \mathbf{X}^{\top} \mathbf{A}^{\top}) \mathbf{W} \mathbf{y} - \mathbf{I} \mathbf{X}^{\top} \mathbf{W}^{\top} \mathbf{A} \mathbf{y} + \mathbf{X}^{\top} \mathbf{A}^{\top} \mathbf{W} \mathbf{X} \boldsymbol{\beta} + \mathbf{X}^{\top} \mathbf{W}^{\top} \mathbf{A} \mathbf{X} \boldsymbol{\beta} \Bigr) \\ &= - \frac{1}{\sigma^2} \Bigl( \mathbf{X}^{\top} \mathbf{A}^{\top} \mathbf{W} \mathbf{y} - \mathbf{X}^{\top} \mathbf{A}^{\top} \mathbf{W} \mathbf{X} \boldsymbol{\beta} + \mathbf{X}^{\top} \mathbf{W}^{\top} \mathbf{A} \mathbf{y} - \mathbf{X}^{\top} \mathbf{W}^{\top} \mathbf{A} \mathbf{X} \boldsymbol{\beta} \Bigr) \\ &= - \frac{1}{\sigma^2} \Bigl( \mathbf{X}^{\top} \mathbf{A}^{\top} \mathbf{W} (\mathbf{y} - \mathbf{X} \boldsymbol{\beta}) + \mathbf{X}^{\top} \mathbf{W}^{\top} \mathbf{A} (\mathbf{y} - \mathbf{X} \boldsymbol{\beta}) \Bigr) \\ &= - \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) \tag{0.44} \end{align}$

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

1: 式(0.37)の偏微分の式を立てます。 $\boldsymbol{\beta}, \lambda$ に関する微分なので、それぞれ他の項は定数として扱います。
2: $\log L(\boldsymbol{\theta})$ の微分に、式(0.40)を代入します。
3: 括弧を展開します。
4: $\boldsymbol{\epsilon}$ に、式(2)を代入します。

　式を変形して置き換えます。

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

5: 和の微分 $\frac{d \{f(x) + g(x)\}}{d x} = \frac{d f(x)}{d x} + \frac{d g(x)}{d x}$ より、微分の和に分割します。
6: $\boldsymbol{\epsilon}$ の内積の微分に、式(15)を代入します。
6: $\mathbf{X}, \boldsymbol{\beta}$ などの二次形式の微分に、式(16)を代入します。
7: $-1$ を括り出します。
8: 共通の項を括り出します。
9: $\mathbf{X}$ を括り出します。

　独立誤差項 $\boldsymbol{\epsilon}$ の転置を $\boldsymbol{\beta}$ に関して微分します。

$\begin{align} \frac{\partial \boldsymbol{\epsilon}^{\top}}{\partial \boldsymbol{\beta}} &= \frac{\partial}{\partial \boldsymbol{\beta}} \Biggl\{ \Bigl( \mathbf{A} (\mathbf{y} - \mathbf{X} \boldsymbol{\beta}) \Bigr)^{\top} \Biggr\} \\ &= \frac{\partial}{\partial \boldsymbol{\beta}} \Bigl\{ (\mathbf{A} \mathbf{y})^{\top} - (\mathbf{A} \mathbf{X} \boldsymbol{\beta})^{\top} \Bigr\} \\ &= \frac{\partial}{\partial \boldsymbol{\beta}} \Bigl\{ \mathbf{y}^{\top} \mathbf{A}^{\top} - \boldsymbol{\beta}^{\top} \mathbf{X}^{\top} \mathbf{A}^{\top} \Bigr\} \\ &= \frac{\partial}{\partial \boldsymbol{\beta}} \Bigl\{ \mathbf{y}^{\top} \mathbf{A}^{\top} \Bigr\} - \frac{\partial \boldsymbol{\beta}^{\top}}{\partial \boldsymbol{\beta}} \mathbf{X}^{\top} \mathbf{A}^{\top} \\ &= \mathbf{0} - \mathbf{I} \mathbf{X}^{\top} \mathbf{A}^{\top} \\ &= - \mathbf{X}^{\top} \mathbf{A}^{\top} \tag{15} \end{align}$

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

1: 式(2)の偏微分の式を立てます。
2: 転置の性質 $(\mathbf{A} + \mathbf{B})^{\top} = \mathbf{A}^{\top} + \mathbf{B}^{\top}$ より、括弧を展開します。
3: 転置の性質 $(\mathbf{A} \mathbf{B})^{\top} = \mathbf{B}^{\top} \mathbf{A}^{\top}$ より、括弧を展開します。
4: 和の微分 $\frac{d \{f(x) + g(x)\}}{d x} = \frac{d f(x)}{d x} + \frac{d g(x)}{d x}$ より、微分の和に分割します。

　説明変数・回帰パラメータ $\mathbf{X}, \boldsymbol{\beta}$ などの二次形式を $\boldsymbol{\beta}$ に関して微分します。

$\begin{align} \frac{\partial}{\partial \boldsymbol{\beta}} \Bigl\{ \boldsymbol{\beta}^{\top} \mathbf{X}^{\top} \mathbf{A}^{\top} \mathbf{W} \mathbf{X} \boldsymbol{\beta} \Bigr\} &= \Bigl( \mathbf{X}^{\top} \mathbf{A}^{\top} \mathbf{W} \mathbf{X} + ( \mathbf{X}^{\top} \mathbf{A}^{\top} \mathbf{W} \mathbf{X} )^{\top} \Bigr) \boldsymbol{\beta} \\ &= \Bigl( \mathbf{X}^{\top} \mathbf{A}^{\top} \mathbf{W} \mathbf{X} + \mathbf{X}^{\top} \mathbf{W}^{\top} \mathbf{A} \mathbf{X} \Bigr) \boldsymbol{\beta} \\ &= \mathbf{X}^{\top} \mathbf{A}^{\top} \mathbf{W} \mathbf{X} \boldsymbol{\beta} + \mathbf{X}^{\top} \mathbf{W}^{\top} \mathbf{A} \mathbf{X} \boldsymbol{\beta} \tag{16} \end{align}$

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

1: 二次形式の微分 $\frac{\partial \mathbf{x}^{\top} \mathbf{A} \mathbf{x}}{\partial \mathbf{x}} = (\mathbf{A} + \mathbf{A}^{\top}) \mathbf{x}$ より、 $\mathbf{X}^{\top} \mathbf{A}^{\top} \mathbf{W} \mathbf{X}$ を1つの行列とみなして微分を計算します。
2: 転置の性質 $(\mathbf{A} \mathbf{B})^{\top} = \mathbf{B}^{\top} \mathbf{A}^{\top}$ より、括弧を展開します。

　対数尤度関数の回帰パラメータと空間パラメータに関する微分(0.44)について、 $\mathbf{A}$ の逆行列を用いて、次のようにも変形できます。

$\begin{align} \frac{\partial^2 \log L(\boldsymbol{\theta})}{\partial \boldsymbol{\beta} \partial \lambda} &= - \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) \tag{0.44}\\ &= - \frac{1}{\sigma^2} \Bigl( \mathbf{X}^{\top} (\mathbf{A}^{\top} \mathbf{W} + \mathbf{W}^{\top} \mathbf{A}) \mathbf{A}^{-1} \mathbf{A} (\mathbf{y} - \mathbf{X} \boldsymbol{\beta}) \Bigr) \\ &= - \frac{1}{\sigma^2} \Bigl( \mathbf{X}^{\top} ( \mathbf{A}^{\top} \mathbf{W} \mathbf{A}^{-1} + \mathbf{W}^{\top} \mathbf{A} \mathbf{A}^{-1} ) \boldsymbol{\epsilon} \Bigr) \\ &= - \frac{1}{\sigma^2} \Bigl( \mathbf{X}^{\top} (\mathbf{A}^{\top} \mathbf{W} \mathbf{A}^{-1} + \mathbf{W}^{\top}) \boldsymbol{\epsilon} \Bigr) \tag{0.44'} \end{align}$

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

1: 式(0.44)を再掲しています。
2: 中の括弧に右から $\mathbf{I} = \mathbf{A}^{-1} \mathbf{A}$ を掛けます。
3: $\mathbf{A}^{-1}$ を中の括弧に掛けます。
3: $\mathbf{A}$ と後の括弧の積を式(2)で置き換えます。

　さらに、対数尤度関数の回帰パラメータと空間パラメータに関する微分(0.44)について考えます。
　空間重み行列 $\mathbf{W}$ が対称行列(2つの地域 $i, j$ 間の重みが等しい $w_{ij} = w_{ji}$ )のとき、 $\mathbf{W}^{\top} = \mathbf{W}$ なので、式(1)より

$\displaystyle \begin{aligned} \mathbf{A}^{\top} &= (\mathbf{I} - \lambda \mathbf{W})^{\top} \\ &= \mathbf{I}^{\top} - \lambda \mathbf{W}^{\top} \\ &= \mathbf{I} - \lambda \mathbf{W} \\ &= \mathbf{A} \end{aligned}$

が成り立ち、 $\mathbf{A}$ も対称行列になります。
　また、 $\mathbf{W}, \mathbf{A}$ の積は

$\displaystyle \begin{aligned} \mathbf{A} \mathbf{W} &= (\mathbf{I} - \lambda \mathbf{W}) \mathbf{W} \\ &= \mathbf{W} - \lambda \mathbf{W} \mathbf{W} \\ &= \mathbf{W} (\mathbf{I} - \lambda \mathbf{W}) \\ &= \mathbf{W} \mathbf{A} \end{aligned}$

が成り立ち、可換(積の入れ換えが可能)になります。
　よって、式(0.44)は、次のように変形できます。

$\begin{align} \frac{\partial^2 \log L(\boldsymbol{\theta})}{\partial \boldsymbol{\beta} \partial \lambda} &= - \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) \tag{0.44}\\ &= - \frac{1}{\sigma^2} \Bigl( \mathbf{X}^{\top} (\mathbf{A} \mathbf{W} + \mathbf{W} \mathbf{A}) (\mathbf{y} - \mathbf{X} \boldsymbol{\beta}) \Bigr) \\ &= - \frac{1}{\sigma^2} \Bigl( \mathbf{X}^{\top} (2 \mathbf{W} \mathbf{A}) (\mathbf{y} - \mathbf{X} \boldsymbol{\beta}) \Bigr) \\ &= - 2 \frac{1}{\sigma^2} \mathbf{X}^{\top} \mathbf{W} \boldsymbol{\epsilon} \tag{0.44''} \end{align}$

　以上で、ヘッセ行列(13)の1行3列・3行1列目の要素の式が得られました。

分散・分散パラメータによる偏微分

　各種の変数の分散パラメータに関する2階微分を求めます。

　ヘッセ行列(13)の2行2列目の要素は、対数尤度関数(0.37)を分散パラメータ $\sigma^2$ に関して2階微分して求まります。

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

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

1: 式(0.37)の偏微分の式を立てます。 $\sigma^2$ に関する微分なので、 $\sigma^2$ 以外の項は定数として扱います。
2: $\log L(\boldsymbol{\theta})$ の微分に、式(9)を代入します。
3: 和の微分 $\frac{d \{f(x) + g(x)\}}{d x} = \frac{d f(x)}{d x} + \frac{d g(x)}{d x}$ より、微分の和に分割します。
3: 微分の計算が分かりやすいように、分数の指数表記 $\frac{1}{x} = x^{-1}$ に書き換えています。
4: べき乗の微分 $\frac{d x^n}{d x} = n x^{n-1}$ より、 $\sigma^2$ を1つの変数とみなして微分を計算します。
5: 分数の指数表記 $\frac{1}{x^n} = x^{-n}$ を戻します。

　以上で、ヘッセ行列(13)の2行2列目の要素の式が得られました。

分散・空間パラメータによる偏微分

　各種の変数の分散パラメータと空間パラメータに関する微分を求めます。

　ヘッセ行列(13)の2行3列目の要素は、対数尤度関数(0.37)を分散パラメータ $\sigma^2$ と空間パラメータ $\lambda$ に関して微分して求まります。

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

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

1: 式(0.37)の偏微分の式を立てます。 $\sigma^2, \lambda$ に関する微分なので、それぞれ他の項は定数として扱います。
2: $\log L(\boldsymbol{\theta})$ の微分に、式(0.40)を代入します。
3: 和の微分 $\frac{d \{f(x) + g(x)\}}{d x} = \frac{d f(x)}{d x} + \frac{d g(x)}{d x}$ より、微分の和に分割します。
3: 微分の計算が分かりやすいように、分数の指数表記 $\frac{1}{x} = x^{-1}$ に書き換えています。
4: べき乗の微分 $\frac{d x^n}{d x} = n x^{n-1}$ より、 $\sigma^2$ を1つの変数とみなして微分を計算します。
5: 分数の指数表記 $\frac{1}{x^n} = x^{-n}$ を戻します。

　以上で、ヘッセ行列(13)の2行3列・3行2列目の要素の式が得られました。

空間・空間パラメータによる偏微分

　各種の変数の空間パラメータに関する2階微分を求めます。

　ヘッセ行列(13)の3行3列目の要素は、対数尤度関数(0.37)を空間パラメータ $\lambda$ に関して2階微分して求まります。

$\begin{align} \frac{\partial \log L(\boldsymbol{\theta})}{\partial \lambda \partial \lambda} &= \frac{\partial}{\partial \lambda} \frac{\partial}{\partial \lambda} \log L(\boldsymbol{\theta}) \\ &= \frac{\partial}{\partial \lambda} \Biggl\{ - \mathrm{Tr}(\mathbf{A}^{-1} \mathbf{W}) + \frac{1}{\sigma^2} \boldsymbol{\epsilon}^{\top} \mathbf{W} (\mathbf{y} - \mathbf{X} \boldsymbol{\beta}) \Biggr\} \\ &= - \frac{\partial \mathrm{Tr}(\mathbf{A}^{-1} \mathbf{W})}{\partial \lambda} + \frac{1}{\sigma^2} \frac{\partial \boldsymbol{\epsilon}^{\top}}{\partial \lambda} \mathbf{W} (\mathbf{y} - \mathbf{X} \boldsymbol{\beta}) \\ &= - \mathrm{Tr}(\mathbf{A}^{-1} \mathbf{W} \mathbf{A}^{-1} \mathbf{W}) + \frac{1}{\sigma^2} \Bigl( - \mathbf{W} (\mathbf{y} - \mathbf{X} \boldsymbol{\beta}) \Bigr)^{\top} \mathbf{W} (\mathbf{y} - \mathbf{X} \boldsymbol{\beta}) \\ &= - \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}) \tag{0.47} \end{align}$

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

1: 式(0.37)の偏微分の式を立てます。 $\lambda$ に関する微分なので、 $\lambda$ 以外の項は定数として扱います。
2: $\log L(\boldsymbol{\theta})$ の微分に、式(0.40)を代入します。
3: 和の微分 $\frac{d \{f(x) + g(x)\}}{d x} = \frac{d f(x)}{d x} + \frac{d g(x)}{d x}$ より、微分の和に分割します。
4: $\mathbf{A}^{-1} \mathbf{W}$ のトレースの微分に、式(17)を代入します。
4: $\boldsymbol{\epsilon}$ の微分に、式(12)を代入します。
5: 転置の性質 $(\mathbf{A} \mathbf{B})^{\top} = \mathbf{B}^{\top} \mathbf{A}^{\top}$ より、括弧を展開します。

　 $\mathbf{A}^{-1} \mathbf{W}$ のトレースを $\lambda$ に関して微分します。

$\begin{align} \frac{\partial \mathrm{Tr}(\mathbf{A}^{-1} \mathbf{W})}{\partial \lambda} &= \mathrm{Tr} \Bigl( \frac{\partial}{\partial \lambda} \Bigl\{ \mathbf{A}^{-1} \mathbf{W} \Bigr\} \Bigr) \\ &= \mathrm{Tr} \Bigl( \frac{\partial \mathbf{A}^{-1}}{\partial \lambda} \mathbf{W} \Bigr) \\ &= \mathrm{Tr} \Bigl( - \mathbf{A}^{-1} \frac{\partial \mathbf{A}}{\partial \lambda} \mathbf{A}^{-1} \mathbf{W} \Bigr) \\ &= \mathrm{Tr} \Bigl( - \mathbf{A}^{-1} (- \mathbf{W}) \mathbf{A}^{-1} \mathbf{W} \Bigr) \\ &= \mathrm{Tr}( \mathbf{A}^{-1} \mathbf{W} \mathbf{A}^{-1} \mathbf{W} ) \tag{17} \end{align}$

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

1: $x$ に依存する行列 $\mathbf{A}(x)$ におけるトレースの微分 $\frac{\partial \mathrm{Tr}(\mathbf{A})}{\partial x} = \mathrm{Tr}(\frac{\partial \mathbf{A}}{\partial x})$ より、微分のトレースに変形します。
2: $\mathbf{A}$ と無関係な項を微分の外に出します。
3: $x$ に依存する行列 $\mathbf{A}(x)$ における逆行列の微分 $\frac{\partial \mathbf{A}^{-1}}{\partial x} = - \mathbf{A}^{-1} \frac{\partial \mathbf{A}}{\partial x} \mathbf{A}^{-1}$ より、元の行列の微分と逆行列の積に変形します。
4: $\mathbf{A}$ の微分に、式(10)を代入します。

　以上で、ヘッセ行列(13)の3行3列目の要素の式が得られました。

ヘッセ行列

　対数尤度関数のヘッセ行列を求めます。

　ヘッセ行列(13)の各要素をそれぞれの式(0.42)(0.43)(0.44)(0.45)(0.46)(0.47)と置き換えます。

$\displaystyle \mathbf{H}(\boldsymbol{\theta}) = \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{13'}$

　ヘッセ行列の式が得られました。

　以上で、対数尤度関数の2階微分の式が得られました。

　この記事では、SEMのヘッセ行列を数式で確認しました。次の記事では最尤法を数式で確認します。

参考文献

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

作者:植杉大
森北出版

Amazon

おわりに

　SEMの記事から読み始めて「いったい何を言っているんだ」となった人は、回り道になりますが一旦OLSの記事に立ち返って読んだ後に、SLMの記事も読んでからSEMの記事(この記事)に戻って読み比べると、仮定や式の違いが対比されて理解しやすくなると思います。試してみてください。そんなことをしなくても理解できる人はなんでもいいです。

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

　さぁ、決戦の日だ⚽

【次の内容】

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

からっぽのしょこ

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

0.2.2：空間誤差モデル(SEM)のヘッセ行列の導出【はじめての地理空間DSのノート】

はじめに

0.2.2 空間誤差モデル(SEM)のヘッセ行列の導出

モデルの確認

定義式

計算式

尤度関数

対数尤度関数の1階微分の導出

対数尤度の微分の設定

回帰パラメータによる偏微分

分散パラメータによる偏微分

空間パラメータによる偏微分

対数尤度の微分

対数尤度関数の2階微分の導出

ヘッセ行列の設定

回帰・回帰パラメータによる偏微分

回帰・分散パラメータによる偏微分

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

分散・分散パラメータによる偏微分

分散・空間パラメータによる偏微分

空間・空間パラメータによる偏微分

ヘッセ行列

参考文献

おわりに