ライントランセクト調査による個体数の区間推定2 - 迷子のオガワコマッコウの放浪日記

前回の記事の続きです。今回は$ \hat{w} $の分散を考えます。そのために、前回扱ったデルタ法を一般の二変数に拡張してみます。

確率変数$ X_1,X_2 $があり、それぞれの期待値を$ \mu_1,\mu_2 $とします。それらを説明変数とした任意の関数$ f(X_1,X_2) $の一次近似式は

f(X_1,X_2) \cong f(\mu_1,\mu_2) + (X_1 - \mu_1) \frac{\partial f}{\partial X_1}(\mu_1,\mu_2) + (X_2 - \mu_2) \frac{\partial f}{\partial X_2}(\mu_1,\mu_2) \tag{1}

となります。

前回は2つの説明変数が独立だったので簡単でしたが、今回は独立でない場合も含めて一般化してみます。そこで、分散と共分散に関する以下の性質を利用します。

Var[aX+bY] = a^{2}Var[X] + b^{2}Var[Y] + 2ab・Cov[X,Y] \tag{2}

これは、分散と共分散に関する以下の基本的性質を用いることで簡単に証明できます。

\displaystyle Var[X] = E[X^2] - (E[X])^2

Cov[X,Y] = E[XY] - E[X]E[Y]

さて、(1)式と(2)式を用いると、$ f(X_1,X_2) $の分散は、

\begin{aligned} Var[f(X_1,X_2)] \cong & Var[X_1] \left( \frac{\partial f}{\partial X_1}(\mu_1,\mu_2) \right)^2 + Var[X_2] \left(\frac{\partial f}{\partial X_2}(\mu_1,\mu_2) \right)^2 \\ & + 2\frac{\partial f}{\partial X_1}(\mu_1,\mu_2)・\frac{\partial f}{\partial X_2}(\mu_1,\mu_2)・Cov[X,Y] \end{aligned}

となります。

ここで、$ \frac{\partial f}{\partial X_1}(\mu_1,\mu_2),\frac{\partial f}{\partial X_2}(\mu_1,\mu_2) $に注目してみると、線形代数に詳しい人なら分かるとおり、二次形式になっています。なので、$ \displaystyle f(X_1,X_2) $の分散は下のように行列形式で表すことが出来ます。

Var[f(X_1,X_2)] \cong \left( \frac{\partial f}{\partial X_1}(\mu_1, \mu_2) \hspace{4mm} \frac{\partial f}{\partial X_2}(\mu_1,\mu_2) \right) \left( \begin{matrix} Var[X_1] & Cov[X_1,X_2] \\[6pt] Cov[X_1,X_2] & Var[X_2] \end{matrix} \right) \left( \begin{matrix} \frac{\partial f}{\partial X_1}(\mu_1, \mu_2) \\[6pt] \frac{\partial f}{\partial X_2}(\mu_1,\mu_2) \end{matrix} \right)

線形代数が苦手な人でも、下の式を展開すると上の式に一致することを確かめてください。

右辺の真ん中の2×2行列のように、対角成分が各変数の分散、それ以外の成分が共分散になっているような行列を分散共分散行列と言います。正確な定義はここでは不要なので割愛します。

また、分散共分散行列の両側にある偏微分で表されたベクトルは勾配ベクトルと呼ばれ、$ \nabla f(\mu_1,\mu_2) $と表記することもあります。

デルタ法を一般の2変数関数に拡張したこの段階で、いよいよ$ \hat{w} $の分散を考えます。$ \hat{w} $は発見関数のパラメータに依存し、またそのパラメータはデータから得られる確率変数です。例えば発見関数がhazard-rate関数の時は以下のように表されるのでした。

\hat{w} = \int_0^{x_M} 1 - \exp \left( - \left( \frac{x}{a} \right)^{-b} \right) dx

$ x_M $は予め設定される定数なので、右辺の中の確率変数はaとbの二つです。つまり、$ \hat{w} $は確率変数a,bを説明変数とした二変数関数であることがわかります。a,bの期待値を最尤推定量$ \hat{a},\hat{b} $で評価すると、デルタ法より、

Var[\hat{w}(a,b)] \cong \left( \frac{\partial \hat{w}}{\partial a}(\hat{a}, \hat{b}) \hspace{5mm} \frac{\partial \hat{w}}{\partial b}(\hat{a},\hat{b}) \right) \left( \begin{matrix} Var[a] & Cov[a,b] \\[6pt] Cov[a,b] & Var[b] \end{matrix} \right) \left( \begin{matrix} \frac{\partial \hat{w}}{\partial a}(\hat{a}, \hat{b}) \\[6pt] \frac{\partial \hat{w}}{\partial b}(\hat{a},\hat{b}) \end{matrix} \right)

となります。偏微分の定義を考えれば、勾配ベクトルの近似値は数値計算で簡単に求まりそうですね。ところが、aとbの分散共分散行列は定義を考えてもうまくいきそうにありません。求めるのは不可能なのでしょうか。

心配ご無用。実は、最尤推定量の分散を考えるのにすごーく便利な法則があるのです。それが下の定理です。

Hessian行列の逆行列は漸近的に(つまりサンプル数を無限大に近づけた時に)最尤推定量の分散共分散行列に一致する。

尤度関数のパラメータが二つの場合は、Hessian行列は以下のように定義されます。

H(\hat{ \theta_1 },\hat{ \theta_2 }) = - \left( \begin{matrix} \frac{\partial^2}{\partial {\theta_1}^2}\log L(X|\hat{\theta_1},\hat{\theta_2}) & \frac{\partial^2}{\partial \theta_1 \partial \theta_2}\log L(X|\hat{\theta_1},\hat{\theta_2}) \\[6pt] \frac{\partial^2}{\partial \theta_2 \partial \theta_1}\log L(X|\hat{\theta_1},\hat{\theta_2}) & \frac{\partial^2}{\partial {\theta_2}^2}\log L(X|\hat{\theta_1},\hat{\theta_2}) \end{matrix} \right)

ここで、$ \hat{ \theta_1},\hat{ \theta_2 } $は各パラメータの最尤推定量を、$ L(X|\theta_1,\theta_2) $は尤度関数を表します。

この定理を使って、aとbの分散共分散行列を求めてみます。aとbの尤度関数は

f(X|a,b) = \prod_{i=1}^{n} \frac{g(x_i)}{w}

と表されるのでした。この関数fを用いると、aとbのHessian行列は

H(\hat{a},\hat{b}) = - \left( \begin{matrix} \frac{\partial^2}{\partial a^2}\log f(X|\hat{a},\hat{b}) & \frac{\partial^2}{\partial a \partial b}\log f(X|\hat{a},\hat{b}) \\[6pt] \frac{\partial^2}{\partial b \partial a}\log f(X|\hat{a},\hat{b}) & \frac{\partial^2}{\partial b^2}\log f(X|\hat{a},\hat{b}) \end{matrix} \right)

となります。

以上の内容をまとめると、$ \hat{w} $の分散は以下のようになります。

Var[\hat{w}(a,b)] \cong \nabla \hat{w}(\hat{a},\hat{b})^\mathrm{T} \cdot H(\hat{a},\hat{b})^{-1} \cdot \nabla \hat{w}(\hat{a},\hat{b})

(Tは転置を表す)

ただし、

\nabla \hat{w}(\hat{a},\hat{b}) = \left( \begin{matrix} \frac{\partial \hat{w}}{\partial a}(\hat{a}, \hat{b}) \\[6pt] \frac{\partial \hat{w}}{\partial b}(\hat{a},\hat{b}) \end{matrix} \right)

次回は$ n $の分散について考えます。ポアソン過程の知識があった方が分かりやすいかもしれません。