人工知能してみる

人工知能の中の人が機械学習とか統計とかAI的なことを書き連ねます

Rを使ったポアソン分布における最尤推定

本格的に統計を学ぶためにデータ解析のための統計モデリング入門、いわゆる緑本を読み始めました。
ナナメ読みで概観を捉えながら、追ってRで手を動かして勉強してます。
統計はツマミ食いでしか学んでいなかったので、こうやって体系的に学ぶと、いままでの知識が繋がって面白いですね。

今回は緑本の最初にあるポアソン分布について実際にRを使って勉強します。
内容は本に書いてあるコードを追っているだけです。
※コード中の"data.RData"は緑本の参考データです。

load("data.RData")
hist(data,breaks = seq(-0.5,9.5,1))

これでこんな感じのヒストグラムが得られます。
f:id:Grahamian:20160912221454p:plain
さて、この得られたヒストグラムについて最も適したパラメータを求めます。
今回は得られたデータがポアソン分布に基づいていると仮定したとこから始まります。
ポアソン分布は平均値λのみをパラメータとして持ちます。
そのため、今回はλのみ考えます。

最も適したλを求めるために最尤推定法を使います。
最尤推定法とは尤度という値の当てはまりのよさを示す統計量です。
これは尤度関数L(λ)で求まります。
L(λ)は各応答変数の出現確率の積です。
このままでは使いにくいので、対数をとったlogLを対数尤度関数として使います。
こうすることで対数尤度logL最も大きい(ゼロに近い)とき分布への値の当てはまりが良い、ということになります。

では、実際にRで求めてみましょう。

logL <- function(m) sum(dpois(data,m,log=TRUE))
lambda <- seq(2,5,0.1)
plot(lambda,sapply(lambda, logL),type="l")

これにより次のようなλvsLogLのグラフが得られます。
f:id:Grahamian:20160912221701p:plain
このグラフを微分したときゼロになるのが最尤推定値です。
ね、簡単でしょ?

ということで、Rを使うことで非常に簡単に最尤推定を行うことができました。
まだまだ分からないことだらけですが、手を動かすのは大事ですね。