一日一膳(当社比)

RとJavaと時々数学

ブラウン運動の最尤推定

 \require{AMScd,amsmath,amssymb,amsfonts,amscd}

 \require{color}

keywords :時系列解析 , 最尤法 , Kullback-Leibuler情報量 , ブラウン運動

今回のテーマは、前回に引き続き時系列解析の話題です。

前回の記事ではワイトノイズ、ブラウン運動、AR(1)過程について 

その平均、自己共分散がどのように変化するかについて書きました。 

sonnamonyaro.hatenablog.com

 

今回の記事では、これらのモデルを実際に応用する局面で最も基本的な問題

 

問題:与えられた時系列データを最も適切に再現するように時系列モデルを決定せよ 

 

について、最尤法という手法により上の問題を扱いたいとおもいます

 

目次

 

最尤法

最尤法は統計学で広く使われる考え方で、パラメーターを持った確率変数

に従うサンプルから、適当なパラメーターを決定する方法 と言えます。

例えば、ここに表が出る確率がp\ (0<p<1)のコインがあったとして、

pを求めろと上司に言われたとします。そうしたら、あなたはおそ

らく(上司に逆らわずに)10回、100回、・・・とコインをひたすら投げ、

表の出た割合をの値として報告するでしょう。

実はそれが、最尤法なのです。

それでは、最尤法の説明にはいります

尤度関数

\theta \in \mathbb{R}をパラメータに持つ確率変数

X_\thetaを考えます。またはX_\theta連続な確率密度f(x,\theta)

を持つと仮定します。先の例でいうと、pがパラメータ\theta

コインを1回投げた時、表が出たら1裏が出たら0とするカウントがX_\theta 

になります。

{\bf definition}

尤度関数(Likelihood function)とは、上のように用意されたパラメーター

(\theta)付き確率変数と、観測されたN個のデータ\mathsf{D}=(x_1,x_2,\cdots,x_N)に対して定まる、\thetaの関数でつぎで定義される

{\bf (LF)}\ \mathcal{L}(\theta|\mathsf{ D})=\mathsf{D}\text{が観測される確率}

もし、個のデータが一回ずつ採取され、それらが独立と仮定する場合には

(先ほどのコインの例)積の法則により

{\bf (LF')}\ \mathcal{L}(\theta|\mathsf{ D})=\prod_{i=1}^{N}f(x_i,\theta)

 この式が、よく皆さんが目にする式ではないかと思います。なお、section2で

は(LF')の式ではなく(LF)の式を使うので注意してください

最尤法は次のような考え方です

 尤度関数を最大化するパラメーター\theta^\astは"適切"な選び方である

ここで導入された、尤度関数を最大化するパラメーターを最尤推定量 

({\bf Maximum \ Likelihood \ Estimator/\  M.L.E.})といいます。覚えておきましょう。

また、確率密度が微分可能ならば、最尤推定量での尤度関数の微分は消えることから

次の尤度方程式が成り立つことがわかります

 \frac{d\ \text{ln}\ \mathcal{L}(\theta^\ast|\mathsf{ D})}{d\ \theta}=0
 

 

 

ここで、適切の意味は次のような事実をさします

{\bf Fact}

最尤推定量は、サンプルを死ぬほど採りまくった場合に、正しい値を教えてくれる

ことが保証されている(強一致性,漸近正規性,...)

これらの事実は統計学のまともな本に譲るとして、うえの事実を納得するための

考え方を紹介するにとどめます 。

Kullback-Leibuler情報量

Kullback-Leibuler 情報量は、2つの確率分布が等しいか否かをはかる量です。

今回の記事では、Kullback-Leibuler情報量と最尤法のかかわりをみてみます。

{\bf definition}
\mathsf{A},\mathsf{B}を確率分布、それらの確率密度f(x),g(x)\ (x \in \mathbb{R})
は連続とする。このときKullback-Leibuler 情報量は次の式で定義される
I(\mathsf{A}||\mathsf{B})=\int_{\mathbb{R}}dx\ f(x){\rm ln} \Bigl(\frac{f(x)}{g(x)}\Bigr)

 

次のことが知られています。

重要性質
\mathsf{A},\mathsf{B}が確率分布であるならば、
I(\mathsf{A}||\mathsf{B})\geq 0\ ,\text{等号成立は}\ \mathsf{A}=\mathsf{B}\ \text{のときに限る}

 

うえの性質より、2つの確率分布が等しいかどうかをみるには、それらのKullback-Leibuler情報量が0かどうかに注目すればよいのです。

 

さてこのことを次のパラメーター推定問題

{\bf Problem}

確率分布{\sf A}(確率密度f(x))と、パラメーター付き確率分布
{\sf A}(\theta)(確率分布f(x,\theta))が与えられたとき、
\mathsf{A}_\thetaの分布が\mathsf{A}に最も分布が近くなる\thetaをあたえよ

に応用してみます。
0\leq I(\mathsf{A}||\mathsf{A}_\theta)
=\int_{\mathbb{R}}dx\ f(x){\rm ln}\ f(x)-\int_{\mathbb{R}}dx\ f(x){\rm ln}\ f(x,\theta)
なので分布\mathsf{A}_\theta\mathsf{A}に近づけるためには、上の式の第二項
\int_{\mathbb{R}}dx\ f(x){\rm ln}\ f(x,\theta)
を最大にすればよいことがわかります。

 

すこしずつ最尤法がにちかずいてきました。もうすこし頑張りましょう。
結局上の積分直接計算するには、f(x)を知らなければなりません。しかし、\mathsf{A}から採ったサンプル\mathsf{D}=\{x_1,x_2,\cdots,x_N\}]があれば、

そこからモンテカルロ法を利用して間接的積分を計算できます:
\int_{\mathbb{R}}dx\ f(x){\rm ln}\ f(x,\theta)\sim \frac{1}{N}\sum_{i=1}^N{\rm ln}\ f(x_i,\theta)

これって対数尤度ですね。

つまり、最尤法の原理「尤度関数を最大化せよ」というのは、真の分布と
モデル分布のKullback-Leibuler情報量が小さいことが期待されるという点で、適当なパラメーターの選び方であると言えるわけです。

ブラウン運動最尤推定

実際に最尤法のコツをつかむために一次元ブラウン運動最尤推定について
解説します。

今回考えるブラウン運動は、2つのパラメーター\alpha,\sigma^2
を持つ時系列モデルX_t\ (t=0,1,2,\cdots)で次の式で与えられるものです

 ({\bf BM1})X_0=0

 ({\bf BM2})X_{t}=\alpha +X_{t-1}+w_t
 ({\bf BM3})w_t \underset{\text{i.i.d.}}{\sim} N(0,\sigma^2)
 
さて、このデータに従う時系列データが得られたとします。

簡単のため、時系列データは時刻1から$T \geq 1$まで連続的に取ったものとし、

それを\mathsf{D}=\{x_0=0, x_1,\cdots,x_T\}
とします。 

尤度関数の計算

今回は、パラメーターは\alpha,\sigma^2なので尤度関数は2変数に
なります。それを、 \mathcal{L}( \alpha,\sigma^2|\mathsf{D})
とおきます。

\mathcal{L}( \alpha,\sigma^2|\mathsf{D})
=\mathbb{P}[ X_1=x_1\  \text{and}\  X_2=x_2 \ \text{and}\  \cdots \ \text{and}\  X_T=x_T]
=\mathbb{P}[ X_1=x_1 \ \text{and}\  X_2-X_1=x_2-x_1 \ \text{and}\  \cdots \ \text{and} \ X_T-X_{T-1}=x_T-x_{T-1}]
ここで式({\bf BM2})よりX_i-X_{i-1}=w_i+\alphaで、これらが({\bf BM3})

より独立であったことから、尤度関数が
\mathcal{L}( \alpha,\sigma^2|\mathsf{D})
\mathbb{P}[ w_1+ \alpha =x_1\  \text{and}\  w_2 +\alpha=x_2-x_1 \ \text{and}\  \cdots \ \text{and}\  w_T +\alpha=x_T-x_{T-1}]
=\prod_{i=1}^T \Bigl(\frac{1}{\sqrt{2 \pi \sigma^2}}\text{exp}\bigl( \frac{(x_i-x_{i-1}-\alpha)^2}{2\sigma^2}\bigr) \Bigr)
ともとまりましたね。目がちかちかするので、対数を取りましょう
\text{ln}\ \mathcal{L}( \alpha,\sigma^2|\mathsf{D})
=-\frac{T}{2}\text{ln}(2\pi \sigma^2)-\frac{1}{2\sigma^2}\sum_{i=1}^T(\alpha-(x_i-x_{i-1}))^2

尤度方程式
ここまでこれば、あとは尤度方程式を解けばオッケーです
0=\frac{\partial \ \text{ln}\ \mathcal{L} }{\partial \alpha}=-\frac{1}{ \sigma^2}\sum_{i=1}^T(\alpha-(x_i-x_{i-1}))
0=\frac{\partial \ \text{ln}\ \mathcal{L} }{\partial \sigma^2}=-\frac{T}{2\sigma^2}+\frac{1}{2(\sigma^2)^2}\sum_{i=1}^T(\alpha-(x_i-x_{i-1}))^2

これを連立して解けば次のように最尤推定量がもとまります
\alpha^\ast=\frac{x_T}{T}
\sigma^{\ast 2}=\frac{1}{T}\sum_{i=1}^T(x_i-x_{i-1})^2- \frac{x_T^2}{T^2}

とめでたく求まりました。

式の形を見れば気づくかもしれませんが、この最尤推定は暗算でやる方法

があるので、考えてみるとよいでしょう。