MyEnigma

とある自律移動システムエンジニアのブログです。#Robotics #Programing #C++ #Python #MATLAB #Vim #Mathematics #Book #Movie #Traveling #Mac #iPhone

最尤推定法による直線近似



前回のメモの続きです.

尤度とは何者なのか? - MY ENIGMA



では,尤度の正体がわかったので,

尤度を使用した最もポピュラーな方法である最尤推定法で

直線近似パラメータを推定してみたいと思います.



ちなみに前回と同様に参考書はこちら.

これなら分かる最適化数学―基礎原理から計算手法まで
これなら分かる最適化数学―基礎原理から計算手法まで



平面上の同一直線上にない(つまり誤差を含む)

N個の点(x_1,y_1),(x_2,y_2),....,(x_N,y_N)がある時,

これらの点を直線近似してみます.



この問題は,各点(x_i,y_i)は直線上の真の位置(x_true-i,y_true-i)から

誤差によってズレてしまっていると考えることができます.

そこで,求めたい直線をy=ax+bとし,

入力x_iに対して,出力y_iにランダムな誤差e_iが加わっているものとすると,

(実は,この誤差のモデル化が最尤推定法の最も重要な部分の一つです.

 ここでのモデル化によって,最終的な答えがかなり変わります)

このモデルは下記の式で表現できます.

 y_i=ax_i+b+e_i\quad i=1,2,....,N

ちなみに,e_iは平均0,分散σ^2の正規分布に従うとします.

(ここで正規分布だと仮定したのは,

 解析的に解きやすいかつ,実際の誤差の分布に近いと思われるからです.)



では,前回の記事の内容で述べた通り,

このモデルの尤度は,

に変形できます.


一行目の左辺は尤度の定義そのものですね.

右辺は,誤差分布が正規分布に従うとしているので,

正規分布の確率密度関数の積和で表されます.

そして,今回の問題では先ほどの式を

x=y_i-axi-bと変形でき,

また平均0の誤差だとしているので,μ=0とできます.

それらを一行目に代入したものが二行目です.

そして,積和の部分を無くすために,

指数の法則を使って綺麗にしたのが三行目ですね.



最尤推定法ではこの尤度を最大化すればいいわけです.

しかし,先ほどの三行目の式はいささか複雑なので,

L=-log(p(y_1)p(y_2).....p(y_N))と置いて,

このLを最小化するような値を決めればいいということになります.

logをとったのは自然対数を消すためで,

マイナスを付けたのは最大化より,最小化のほうが少しだけわかりやすいからです.



すると,

式変形できます.

今回は,本当に高校で習った対数の公式を使って変形しただけです.



ここで,σやNは定数なので,

(今回の問題ではσの値はわかりませんが,ある一特定の値だとしています.

 これがデータによって変わると話はもっと複雑になります)

第二項は定数です.

すると,第一項の値が最小になるようなaとbを求めれば

最尤推定法によるパラメータ設定になります.



実はこの第一項からσを抜くと,

最小二乗法による線形近似とまったく同じ値になります.

つまり,最小二乗法は,

出力値に正規分布の誤差が生じていると仮定した場合の最尤推定法と同じなのです.

最小二乗法は最尤推定法の特殊な場合ともいえます.



これで,最尤推定法による線形近似ができました.

めでたしめでたし.




ちなみに今回つかった数式のTexコードはこちら.

p(y_1)p(y_2),...,p(y_N)&=& \prod^N_{i=1} \frac{e^{-(x-\mu)^2/2\sigma^2}}{\sqrt{2\pi\sigma^2}}\\
&=& \prod^N_{i=1} \frac{e^{-(y_{i}-ax_{i}-b)^2/2\sigma^2}}{\sqrt{2\pi\sigma^2}}\\
&=& \frac{e^{-\sum^N_{i=1}(y_{i}-ax_{i}-b)^2/2\sigma^2}}{(2\pi\sigma^2)^{\frac{N}{2}}}\\


L&=&-log(\frac{e^{-\sum^N_{i=1}(y_{i}-ax_{i}-b)^2/2\sigma^2}}{(2\pi\sigma^2)^{\frac{N}{2}}})\\
&=&-log(e^{-\sum^N_{i=1}(y_{i}-ax_{i}-b)^2/2\sigma^2})+log{(2\pi\sigma^2)^{\frac{N}{2}}})\\
&=&\frac{1}{2\sigma^2}\sum^N_{i=1}(y_{i}-ax_{i}-b)^2+\frac{N}{2}log2\pi\sigma^2