著名な倫理学者ピーター・シンガー氏は、日本のIWC脱退をどう見るか

Peter Singer(ピーター・シンガー)氏は、応用倫理学の世界的権威として知られており、2005年にタイム誌によって世界の最も影響力のある100人にも選ばれた著名な倫理学者です。

彼自身は動物の権利論や菜食主義思想の新たな枠組みとして独自の道徳規範を提唱し、過去40年以上にわたる論争によって、生命倫理学に関する知見を飛躍的に進歩させてきました。

彼の理論はほぼ全て利益に対する平等な配慮という唯一つの道徳規範に基づいており、今までに

など数多くの著書を執筆しながら、多くの社会問題に関して応用倫理学的見地から疑問を投じています。

また彼は一部の動物種に対して人格(Personhood)が認められると主張し、従来のパーソン論と動物倫理学との橋渡しを行うなど、これまでの倫理学の枠組みに大きな変化をもたらしてきました。

そんな彼が先日、カナダのThe Globe and Mail誌に社説を掲載しました。内容は日本のIWC脱退に対するもので、現代倫理学的観点から見た捕鯨のあり方と問題点、そしてIWCが果たすべき役割について執筆なされています。

今回の記事は、この社説の全文和訳です。捕鯨問題の説明資料として是非ご活用頂ければと思います。

稚拙な和訳ですがご了承ください。

「海には多くのクジラがいるが、人間が獲るための『資源』ではない」

12月26日、日本は国際捕鯨委員会IWC)から脱退することを発表した。菅義偉内閣官房長官は沿岸地域における捕鯨事業の文化的重要性を強調したうえで、IWCが鯨類の保護ばかりに注目し、持続的な捕鯨事業の発展という本来の目的を見失っていると主張した。

IWCが本来の目的に沿っていないというのは否定しがたい事実だ。IWCはもともと1946年に採択された国際捕鯨取締条約に基づいて設立された機関である。同条約の前文では鯨類を「大きな天然資源」と捉えており、また条約の目的が「鯨類資源の適正な保全を行い、それにより捕鯨事業の秩序ある発展を実現する」ということが明記されている。

最初の25年間、IWCはその方針に従ってきた。しかし1970年代以降、鯨類に対する情勢は変化し始める。捕鯨を続行するために資源の非持続的な乱獲を防止する目的で加盟していた国々が、市民の意見を積極的に取り入れ始めたのだ。その結果、1986年に商業捕鯨モラトリアムが採択された。すべての鯨種に対していかなる捕鯨によっても持続可能性が危ぶまれるとは言い難い現在においても、このモラトリアムは解除されていない。

日本は公然とモラトリアムに違反していたわけではない。しかし、科学的研究のために鯨を殺すことができるという条約の抜け穴をくぐり抜けてきた。日本は科学的研究という名目で毎年約300頭の鯨を捕獲していたのである。鯨肉を食べたいという人は年々減少しているにもかかわらず、鯨類の死骸は捕鯨船により回収され、その肉を鯨肉として販売していた。

2010年、オーストラリアは日本を国際司法裁判所に提訴し、その結果日本の捕鯨は実質的には商業捕鯨であり、条約に違反するという判決がなされた。しかし日本は研究計画にわずかな変更を加えただけで、その後も以前とほぼ同じ数のクジラを捕殺し続けた。

昨年9月、ブラジルのフロリアノポリスで開かれた会合で、商業捕鯨モラトリアムを続行するというブラジルの提案(フロリアノポリス宣言)が賛成40、反対27で可決され、IWCの目的に変化が芽生えた。捕鯨はもはや不必要な経済活動とみなされるようになったのである。国家の意地にかけても捕鯨を存続したかった日本にとって、この採決は無益なIWCへ加盟し続けるための最後の藁となった。

しかし我々が無視できない事実は、IWCにおける情勢の変化を招いた鯨類への新たな姿勢は、神聖な大型動物を殺すことに対する感情的な嫌悪によるものでも、西洋の価値観の押し付けによるものでもないということだ。その根底にあるのは、鯨類に関する科学的知見の発展と、人間を超えて他の種にまで道徳的配慮の輪を広げようとする人類の倫理的進歩である。こういった道徳的配慮は、あらゆる生命に対して畏敬の念を抱くという日本の仏教における戒律に非常に類似している。

1946年以降鯨類に関する知見は飛躍的に進歩し、彼らが巨大な脳を持つ社会的な動物であること、多様な鳴音を用いて他の個体とコミュニケーションをとることが分かってきた。彼らは自分の子供や社会的な群れと強く結ばれている。クジラは非常に長生きする動物で、特にホッキョククジラは他のどの哺乳類よりも寿命が長く、ある個体からは200年も前の牙の先端が肉の中から発見されている。他のクジラも40年以上という長い寿命を持っている。また、彼らは感情と痛覚の両方を併せ持つと言われている。それも物理的苦痛を感じるのみならず、群れの仲間や自分の子供を失った際に精神的苦痛を伴う可能性が非常に高い。

それ故に鯨類は、国々が石炭を資源として蓄えるという意味での「資源」ではない。また、小麦畑のように収穫されるような「資源」でもない。彼らは一つの道徳的存在であり、彼らの生き方に従って彼らの行く末が決定される。

現代の商業捕鯨では、移動する標的に向けて移動する船舶から爆発銛を発射することでクジラを捕殺している。この方法を用いてクジラの急所に命中させ、相手の意識を一瞬で消失させることは非常に困難である。また一瞬で絶命させるほどの爆発を起こそうとするとクジラの体が吹っ飛んで粉々になってしまうが、捕鯨船の乗組員の目的はなるべく無傷で鯨体を回収することであるため、あまり多くの爆薬を使用することを望まない。もし仮に我々がクジラを食べないと死んでしまうという状況であれば、感受性を持つ社会的動物に対してそのような捕殺方法を取ることも必要悪と認められるかもしれない。しかし日本含む飽食の国にとっては、それは正当性を持たない。

捕鯨事業が古来からの遺産であるような日本の一部の地域のおいても、捕鯨は十分な正当性があるとは言えない。中国における纏足は古来からの文化的遺産であったが、女性にとっては多大な負担であった。今となっては、この文化が衰退し過去の出来事となったことは好ましいことだ。捕鯨も同じ道を辿るべきであろう。

いや、実際そうなるかもしれない。日本がIWCを脱退すれば、「科学的研究」の名目のもとで南大洋で行っている捕鯨は中止となる。この事実を認識したうえで、日本は、自国の領海または排他的経済水域EEZ)内、つまり、国土の周囲約450万平方キロメートルの海域でのみ捕鯨を行うことを宣言した。面積だけで言えば十分広いのだが、南大洋と比較するとクジラの生息数は遥かに少ない。また日本は持続可能な捕鯨を計画しているため、捕鯨船による捕獲枠は大きく制限される。

我々は日本がIWCを脱退したという事実にうろたえるのではなく、むしろ1994年にIWCで採択された南大洋サンクチュアリが、今では粗暴な「科学的研究」が二度と行われることのない、本当の意味での「聖域」になるということを喜ぶべきなのかもしれない。

一部の日本国民を含む多くの国の人々が持つ正当な道徳的関心を無視し、日本はIWCからの脱退を決定した。これにより日本は国際社会からの孤立への道を一歩進むことになる。しかし次世代の日本で政権を握る人々は、間違いなくこれを誤ったステップと捉え、そして再び逆転させようと考えるだろう。

ちなみにPeter Singer氏は以前にも捕鯨問題に関する記事を執筆なさっています。こちらは別の方が既に全訳されているようなのでリンクを貼っておきます。

davitrice.hatenadiary.jp

2018年度IWC総会 「フロリアノポリス宣言」の和訳

2018年度IWC総会にて採択された、"THE FLORIANÓPOLIS DECLARATION ON THE ROLE OF THE INTERNATIONAL WHALING COMMISSION IN THE CONSERVATION AND MANAGEMENT OF WHALES IN THE 21st CENTURY"の和訳です。

ソース:

 

商業捕鯨モラトリアムの必要性や、21世紀の国際社会におけるIWCの役割などが記載されています。

稚拙な和訳ですがご了承ください。誤訳などありましたらコメント欄にてご指摘いただけると幸いです。

フロリアノポリス宣言(和訳: Kogia_sima)

国際捕鯨委員会(IWC)が鯨類保護や捕鯨管理の役割を担う中心的な国際機関として広く認識されている一方、

鯨類研究や管理の代替的手法又は鯨類資源の持続的な利用の発展、及び国際捕鯨取締条約(ICRW, 1946〜)の制定以来の数多くの国際法の制定によって、委員会は、鯨類保護に向けた多くの解決策及び鯨類資源の非致死的管理を含む多様な附表修正、並びに海洋生態系における重要な生態的役割や炭素循環への貢献を満たすための健全な鯨類資源の維持に関する役割を担うことを認識し、

生存的文化的目的により鯨類に依存している先住民の需要を支援することの重要性を(加盟国が)認めているにも関わらず、鯨類やその生息域の保護に関する幅広い人々の関心を満たすことにおいては、IWCの権利の使用方法について加盟国間で意見が分かれていることを認め、

鯨類の非致死的利用に関する決議(Resolution 2007-3)を再考し、鯨類が生態系において重要な役割を担い、また自然環境や人々に利益をもたらす存在であって、さらに持続可能な非致死的・非抽出的な利用が急速に高まりつつあり、これらが発展途上国を中心として世界中の沿岸地域共同体に大きな社会経済的利益をもたらしていることを認め、

1986年に施行された商業捕鯨モラトリアムがいくつかの鯨類の個体群の回復に貢献していることを再確認し、漁網への絡まり(entanglement)、混獲、騒音、船舶との衝突、海洋ゴミ、及び気候変動など、鯨類に対する過去及び現在の脅威が累積していることを認識し、

沿岸地域共同体に対して科学的知見、事業、収益をもたらす非致死的な活動が行われている海域において、締約国政府の大多数の支持のもと、ICRW第5条に従って加盟国から頻繁に鯨類サンクチュアリが提案されていることを確認し、

さらに国際捕鯨委員会の独立審査に対する回答に関する決議(Resolution [2018-X])を確認する。

これらの前提のもと、委員会は、

鯨類の個体数を近代捕鯨以前の水準まで回復させることの責任を含む21世紀の国際捕鯨委員会の役割に同意し、この背景のもと商業捕鯨モラトリアムの重要性を再確認し、

現代において鯨類の非致死的調査法が多数存在することを認め、またそれによって致死的な調査法の適用が不必要であることに同意し、

漁師の安全や鯨類に対する動物福祉のため、先住民に利益をもたらす先住民生存捕鯨IWCにより管理・保護されるように努め、

Resolution [2018-X] に従ってOEWG(Working Group on Operational Effectiveness)により提案された計画を実行する際、鯨類の保護資金の必要性や非致死的管理に係る問題を考慮に入れるよう各小委員会等に指示し、

(IWCにより提案される南大西洋サンクチュアリ案の内容を確認するとともに、)移動性野生動物種の保全に関する条約委員会の2017年の第12回総会において採択された、南大西洋における鯨類及びその生息域の管理と保護に関する決議の内容を確認し、南大西洋沿岸各国に対してサンクチュアリの適正な施行のために協力を依頼し、

持続可能な非致死的利用の推進などの鯨類保護活動を連携させるため、生物の多様性に関する条約、移動性野生動物種の保全に関する条約、南極の海洋生物資源の保存に関する条約、および世界観光機関など、関連する他の条約や組織とのさらなる協力を求めるよう事務局に依頼し、さらに

この宣言の内容が、国際連合事務総長、国際連合総会、国連環境計画、移動性野生動物種の保全に関する条約、生物の多様性に関する条約、南極の海洋生物資源の保存に関する条約、絶滅のおそれのある野生動植物の種の国際取引に関する条約海洋法に関する国際連合条約、その他委員会で一般的な議論や協力が行われているような国際条約に取り入れられるよう、事務局に依頼する。

本文を引用する際は、なるべく原文の該当する部分も併せて載せるようにお願いします。

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

さて、$ \hat{w} $、$ n $、$ \bar{s} $の分散をそれぞれ考えてきました。また一番最初の記事で紹介した

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

という式を用いれば、生息数密度Dの分散が求まります。これをもとに、いよいよDの信頼区間について考えていきたいと思います。

 

$ \hat{D} $は通常、対数正規分布に従うと仮定されます。なので、$ \log(\hat{D}) $は正規分布に従います。正規分布に従う確率密度関数の95%信頼区間は、統計の教科書に載っているはずなので省略します。信頼区間は以下のようになります。

\log(\hat{D}) - 1.96\sqrt{Var[\log(\hat{D})]} \leqq \log(D) \leqq \log(\hat{D}) + 1.96\sqrt{Var[\log(\hat{D})]} \tag{1}

となります。

ここで、前回の行った計算と同様にすると、$ \mathcal{N}(μ, σ) $に従う確率変数を$ X $としたとき、$ \exp(X) $の期待値と分散はそれぞれ$ \exp(μ+σ/2), \exp(2μ+σ) \cdot (\exp(σ)-1) $となります。これとここから逆に考えると、期待値$ m $,分散$ v $の対数正規分布に従う確率変数を$ Y $とした時、$ \log(Y) $の期待値と分散はそれぞれ、

 

Var[\log(Y)] = \log \left(\frac{v}{m^2}+1 \right)

となります。2つ目の式より、

Var[\log(D)] = \log \left(\frac{Var[ \hat{D} ]}{\hat{D}^2} + 1\right)

となります。これを(1)式に代入すると、

\log(\hat{D}) - 1.96\sqrt{\log \left(\frac{Var[\hat{D}]}{\hat{D}^2} + 1\right)} \leqq \log(D) \leqq \log(\hat{D}) + 1.96\sqrt{\log \left(\frac{Var[\hat{D}]}{\hat{D}^2} + 1\right)}

となります。$ D $の95%信頼区間を求めるには、各辺にexpをつけてあげればよいだけです。結果は以下のようになります。

\frac{\hat{D}}{C} \leqq D \leqq C\hat{D}

ただし、

C = \exp \left(1.96\sqrt{\log \left(\frac{Var[\hat{D}]}{\hat{D}^2} + 1\right)} \right)

です。

 

以上で説明は終わりですが、ここまでの記事は面倒がっていたせいでかなり適当になっているので、今後時間があるときに修正していこうと思っています。またもし分からない部分、間違いなどありましたら、コメント欄に書いていただければと思います。それではここまでこの記事をお読みいただき、ありがとうございました。

 

ここまでの参考資料

Buckland S.T., Anderson D.R., Burnham K.P., Laake J.L.(1993), Distance Sampling:  Estimating Abundance of Biological Populations, Chapman and Hall, London.

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

今回は$ \bar{s} $の分散について説明します。

 

 まず、何も調整を行わずに、単純に得られたデータの群れサイズを平均する場合です。この場合は、

\bar{s} = \frac{1}{n}\sum_{i=1}^{n} s_i

とするだけです。siはi番目の群れサイズを表します。分散は、

\begin{aligned} Var[\bar{s}] = \frac{1}{n^2}\sum_{i=1}^{n} Var[s_i] & = \frac{1}{n^2}・n・\frac{1}{n-1} \sum_{i=1}^{n} \left( s_i - \bar{s} \right)^2 \\ & = \frac{1}{n(n-1)} \sum_{i=1}^n \left(s_i - \bar{s} \right)^2 \end{aligned}

となります。

 

しかしこれには問題があります。通常ライントランセクト調査では、トラックラインから遠い場所だと、群れサイズが大きい群ればかり発見され、群れサイズが小さい群れのほうが見逃しやすいという事態が発生するのです。

このような場合、単純にデータの群れサイズを平均しただけでは群れサイズを過小推定する恐れがあるため、群れサイズの調整を行うことが多いです。

ではどのように調整すればよいのでしょうか。実際には、線形回帰分析により平均群れサイズを推定する手法がしばしば用いられます(Buckland et al., 1993)。モデル式は以下の通りです。

z = \log(s) = a+bg(x)+\epsilon, \hspace{8mm} \epsilon ~ \mathcal{N}(0,\sigma^2)

つまり、データ上では$ z=\log(s)とg(x) $は一次関数の関係になっていて、誤差項は正規分布に従う、ということを仮定するのです。

ここからどのように群れサイズ平均を推定するかというと、すべての個体が発見される場合、つまりg(x)=1の地点におけるsの確率分布が、真の群れサイズの確率分布と同じになると考えるのです。したがって、$ g(x)=1 $におけるsの確率分布を考えれば良いということになります。

大半の人は回帰分析を大学1年で習うので基礎的な部分は省略しますが、パラメータaとbの推定値は以下のようになります。

\hat{a} = \bar{z} - b \bar{g}(x) \tag{1}

\hat{b} = \frac{S_{xy}}{S_{xx}}

ただし、

S_{xy} = \frac{1}{n} \sum_{i=1}^{n} (g(x_i) - \bar{g}(x))(z_i - \bar{z})

S_{xx} = \frac{1}{n} \sum_{i=1}^{n} (g(x_i) - \bar{g}(x))^2

です。

モデル式より、g(x)=1の時のときのzの$ z_0 $とすると、

z_0 = a+b+\epsilon

(1)のaの推定値を代入すると、

z_0 = \bar{z} - b \bar{g}(x) +b+\epsilon

\hspace{8mm} = \bar{z}+b(1-\bar{g}(x)) +\epsilon

よって、

\begin{aligned}Var[z_0] = & Var[\bar{z}] + (1-\bar{g}(x))^{2}Var[b] + Var[\epsilon] + \\ & 2(1-\bar{g}(x))Cov[\bar{z},b] + 2(1-\bar{g}(x))Cov[b,\epsilon] + 2Cov[\bar{z},\epsilon] \end{aligned} \tag{2}

となります。ただしここでは、$ \epsilon $はデータと独立なため、$ C ov[b, \epsilon], C ov[\bar{z}, \epsilon] $はともに0になります。

 

回帰分析に関する本はたくさんあるので詳しい説明はそちらに譲りますが、各分散、共分散の値は以下のようになります。

Var[\bar{z}] = \frac{\sigma^2}{n}

Var[b] = \frac{\sigma^2}{S_{xx}}

Var[\epsilon] = \sigma^2

Cov[\bar{z},b] = 0

となります。ただし、$ σ^2 $は残差分散で、その不偏推定量

\hat{\sigma^2} = \frac{1}{n-2} \sum_{i=1}^{n} (z_i - a+bg(x_i))^2

と求められます。

これらの値を式(2)に代入して整理すると、

Var[z_0] = \sigma^2 \left(1 + \frac{1}{n} + \frac{(1-\bar{g}(x))^{2}}{S_{xx}} \right)

となります。

 

以上より、zは期待値$ a+b $、分散$ Var[z_0] $の正規分布に従うことが分かりました。しかし今求めたいのは、$ s_0=\exp(z_0) $の確率分布です。

正規分布に従う確率変数を指数変換した時、その分布は対数正規分布と呼ばれるものになります。)

期待値と分散を考えてみます。これらは、正規分布確率密度関数まで戻り、期待値と分散の定義に従って積分で求めることが出来ます。

E[s_0] = E[\exp(z_0)] = \int_{-∞}^{∞}\frac{\exp(z)}{\sqrt{2 \pi Var[z_0]}}\exp \left(-\frac{(z-(a+b))^2}{2Var[z_0]} \right) dx

途中計算は面倒なので省略します。exp()の中身を平方完成して整理していくと、以下の式が導けます。

E[s_0] = \exp \left(a+b+\frac{Var[z_0]}{2} \right)

同様に$ E[s_0^2] $は

E[{s_0}^2] = \exp(2(a+b)+2Var[z_0])

となります。故に、

Var[s_0] = E[{s_0}^2] = E[s_0]^2 = \exp(2(a+b)+Var[z_0])・(\exp(Var[z_0])-1)

となります。これが、$ g(x)=1 $となる点における群れサイズsの確率分布の期待値と分散です。

これをもとに群れサイズの平均について考えれば良いわけです。平均の期待値は$ E[s_0] $に一致します。分散は、上で求めた$ s_0 $の分散をnで割れば良いだけです。

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

さて今回は発見群数$ n $の分散について考えようと思います。

 

まずは皆さんに問題です。ある海域において、1kmのライントランセクト調査を行った翌日、同じ海域で2kmの調査を行ったとします。ただし、天候や有効探索幅などの条件はすべて同じであると仮定します。この時、2つの調査における発見群数の期待値や分散の関係はどうなるでしょうか。

 

発見群数の期待値は、もちろん調査距離に比例します。では分散はどうでしょうか。2kmの調査の記録を、半分までの1kmと残りの1kmに分けて考えます。双方の調査は完全に独立です。なので、総発見群数の分散は1kmの時に比べて2倍になります。

3kmのときはどうでしょうか。同じように考えると、分散は3倍になりそうですね。

 

このように考えていくと、発見群数の分散も調査距離に比例することが分かるかと思います。つまり、発見群数の期待値と分散は常に比例することが分かります。

 

 

この性質を使って、発見群数nの分散を求める方法を説明します。ライントランセクト調査ではこれを求めるために、データを分割して部分トランセクトというものに分けます。この分け方は通常、日付によって分けたり、あるいは調査コースをいくつかに分割することによって分けます。

各部分トランセクトにおける調査距離を$ L_i $とし、全体の調査距離$ L_i $の総和を$ L $とします。同様にして、各部分トランセクトにおける発見群数を$ n_i $とし、全体の発見群数($ n_i $の総和)をnとします。さらに$ n $の期待値($ E[n] $)を$ n_0 $とします。

先ほど見たように、発見群数の期待値と調査距離は比例するので、i番目の部分トランセクトにおける発見群数の期待値は、

E[n_i] = \frac{L_i}{L} n_0 \tag{1}

となります。同様に分散も調査距離に比例するので、

Var[n_i] = \frac{L_i}{L} Var[n]

この式を以下のように変形してみます。

Var[n] = \frac{L}{L_i} Var[n_i] \tag{2}

ここで、$ n_i $の分散の定義より、

\begin{aligned} Var[n_i] &= E[(n_i - E[n_i])^2] \\ &= E\left[\left(n_i - \frac{L_i}{L} n_0 \right)^2\right] \end{aligned} \tag{3}

(3)を(2)に代入すると、

Var[n] = E\left[\frac{L}{L_i} \left(n_i - \frac{L_i}{L} n_0 \right)^2\right]

故に、$ n $の分散の不偏推定量は、

\hat{Var[n]} = \frac{L}{L_i} \left(n_i - \frac{L_i}{L} n_0 \right)^2

これはすべてのiについて成り立つので、さらに精度を上げるために、これらをiについて平均します。

\hat{Var[n]} = \frac{1}{s} \sum_{i=1}^{s} \frac{L}{L_i} \left(n_i - \frac{L_i}{L} n_0 \right)^2

ここで、sは部分トランセクトの数を表します。

実際には、$ n_0 $はデータから求めることが出来ません。そこで、$ n $で代用することになるのですが、統計学を学んでいる人は、不偏分散の式の分母が(標本数)-1になることを知っていると思います。これは、期待値を平均で代用した時に不偏推定量ではなくなる、という性質によるものでした。今回も同じです。そのまま$ n $で代用してしまうと、不偏推定量ではなくなるという事態が発生します。従って、$ n $の分散の不偏推定量は以下のようになります。

\hat{Var[n]} = \frac{1}{s-1} \sum_{i=1}^{s} \frac{L}{L_i} \left(n_i - \frac{L_i}{L} n \right)^2

暇ができれば厳密な証明をこちらに書くかもしれません。

 

ちなみにですが、n/Lの値のことを遭遇率と呼ぶことがあります。遭遇率の分散の不偏推定量は、上で求めた式と、分散の基本的性質を使えば、下のようになります。

\begin{aligned} \hat{Var[\frac{n}{L}]} &= \frac{1}{L^2} \hat{Var[n]} \\ &= \frac{1}{s-1} \sum_{i=1}^{s} \frac{1}{{L_i} L} \left(n_i - \frac{L_i}{L} n \right)^2 \\ &= \frac{1}{s-1} \sum_{i=1}^{s} \frac{L_i}{L} \left(\frac{n_i}{L_i} - \frac{n}{L} \right)^2 \end{aligned}

 

関連資料:

岡村寛(2004), 「海産哺乳類を中心とした生態系モデリングのための数理統計学的研究」, 水研センター研報 No.10, pp.18-100

ライントランセクト調査による個体数の区間推定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}) $をそれぞれ求めてあげればよいことが分かります。これらの値の求め方を次回以降解説していくことにします。