ライントランセクト調査による個体数の区間推定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)

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)

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

ライントランセクト調査による個体数の区間推定1

ライントランセクト調査において、個体数密度の点推定値は以下のように与えられます。

\hat{D} = \frac{n \bar{s}}{2L \hat{w}} \tag{1}

ここで、$ n $は発見群数、$ \bar{s} $は調整された群れサイズ平均、$ L $は調査コースの全長、$ \hat{w} $は有効探索幅の推定値を表します。この式の導き方についてはわかりやすい資料や本などがたくさんあると思うので、省略します。

個体数密度Dの点推定値は求まりましたが、それがどの程度信頼性を持つのかをはっきりさせなければいけません。例えば下のIWCのHPでは、個体数の点推定値とともに95%信頼区間が示されています。

https://iwc.int/estimate#table

このブログでも、同様に個体数の95%信頼区間を求めることを最終目標にしたいと思います。なるべく高度な数学は使わないようにはしますが、それでも学部1年レベルの統計学や微分積分学の知識がどうしても必要になると思います。

さて、まずは$ \hat{D} $の分散について考えてみます。式(1)の右辺に注目してください。このうち、Lはあらかじめ設定される値ですが、$ n, \bar{s}, \hat{w} $は調査の結果得られた値で、これらは確率的に変動します。つまり、これらの値はすべて確率変数です。つまり、$ \hat{D} $の分散を求めるためには、確率変数の積や商の分散を考えなければなりません。そのときに便利なのが、デルタ法と呼ばれる方法です。

<デルタ法>

互いに独立な確率変数X、Yがあり、それぞれの期待値をu、vとします。この時、積XYは一次近似式を用いて次のように変形できます。

\begin{aligned} XY  &\cong  uv + (X-u)\frac{\partial XY}{\partial X}(u,v) + (Y-v)\frac{\partial XY}{\partial Y}(u,v) \\ &= uv + (X-u)v + (Y-v)u \\ &= vX + uY - uv \end{aligned}

XとYは互いに独立なので、分散の性質($ Var(aX+bY)=a^2 Var(X)+b^2 Var(Y) $)より、

Var(XY)  \cong  v^2 Var(X) + u^2 Var(Y)

となります。商X/Yについても同様にすると、

\begin{aligned} \frac{X}{Y} &\cong \frac{u}{v} + (X-u)\frac{\partial \frac{X}{Y}}{\partial X}(u,v) + (Y-v)\frac{\partial \frac{X}{Y}}{\partial Y}(u,v) \\ &= \frac{u}{v} + \frac{X-u}{v} - (Y-v)\frac{u}{v^2} \\ &= \frac{1}{v}X - \frac{u}{v^2}Y + \frac{u}{v} \\ &= \frac{u}{v} \left(\frac{X}{u} - \frac{Y}{v} + 1 \right) \end{aligned}

Var\left(\frac{X}{Y}\right)  \cong  \left(\frac{u}{v}\right)^2 \left(\frac{Var(X)}{u^2} + \frac{Var(Y)}{v^2}\right)

式(1)にデルタ法を用いることにより、個体数密度の分散の近似値を得ることができます。結果だけ示すと、以下のようになります。

Var(\hat{D})  \cong  \left( \frac{n \bar{s}}{2L \hat{w}} \right)^2 \left(\frac{Var(n)}{n^2}+\frac{Var(\bar{s})}{\bar{s}^2}+\frac{Var(\hat{w})}{\hat{w}^2} \right)

以上より、$ Var(\hat{D}) $の値を求めるには、$ Var(n), Var(\bar{s}), Var(\hat{w}) $をそれぞれ求めてあげればよいことが分かります。これらの値の求め方を次回以降解説していくことにします。