幾何分布関係メモ:Rで乱数発生,最尤推定,負の二項回帰
幾何分布ってなに
幾何分布(geometric distribution)は,離散分布のひとつ。
表になるまでコイン投げを繰り返すとして,はじめて表になった回数(または,表になるまでにかかった施行回数)などにあてはまる。もちろん,無記憶だということで,各回の施行は独立。はじめて表になった回数と失敗した回数は違うので,明確に区別しなければならない。ここでは前者を基本的な前提とする。
基本的に,コイン投げの回数をk = 1, 2, 3,...とし,(イカサマ)コインで表になる確率をpとしたとき,累積分布関数は1-(1-p)^kであたえられる。はじめにひとつのケースが,k回目に成功する確率は,(1-p)^(k-1)*pというわけ。
ちなみに,平均は,1/p,分散は(1-p)/p^2。最頻値は常に1。
基本:Rの関数
累積分布関数は,そのまま
1-(1-p)^k
xは回数なので,1:10のようにする。seq(0,100,1)みたいなのでもいい。実数はとれない。pは,(独立した)単独回の成功確率ね。pgeomという関数もあるけど,rのpgeom関数は,失敗した回数の分布を考えているので。
そっちでいいのだったら,
pgeom(k,p)
がもちろん,楽。
累積でない分布関数は,dgeomを使っちゃおう。
dgeom(k, p)/p
でいいね。/pで調整している。
乱数を作るのなら,
rgeom(n, p)+1
1を足すのも,rの関数が失敗した回数の分布によるから調整。
得られたデータからpを推定するのには,簡単で,方程式を解くとよい。
length(dat)/sum(dat)
としたら最尤推定値になる。しかしこれは成功回数の分布。
いろんな分布をサポートしている最尤推定の関数がパッケージMASSにあるので,使ってもいい。
library(MASS)
fitdistr(dat,densfun="geometric")
ただし,fitdistr関数も失敗した回数の分布を想定しているので,
fitdistr(dat-1,densfun="geometric")
としたらよい。fitdistrは対数尤度を返してくれる。
fitdistr(dat-1,densfun="geometric")$loglik
これを,AICに渡すこともできる。
AIC(fitdistr(dat-1,densfun="geometric"))
ここまでの例をあげる。
rgeom(30, .2)+1->d
> d
[1] 2 12 11 7 2 8 10 2 10 2 4 2 4 3 3 4 4 2 6 2 14 1 3 1 5 13 4 4 1> length(d)/sum(d)
[1] 0.2040816> fitdistr(d-1,densfun="geometric")
prob
0.20408163
(0.03324127)> fitdistr(d-1,densfun="geometric")$loglik
[1] -71.55461> AIC(fitdistr(d-1,densfun="geometric"))
[1] 145.1092
これは,たまたまいい推定値だね。
負の二項回帰
負の二項分布(negative binomial distributioin)は,幾何分布を一般化したもので,初回の成功回数ではなく,r回の成功回数が得られるまでの分布。なので,r = 1の負の二項分布が幾何分布というように考える。
幾何分布回帰ってのはなくて,負の二項回帰(negative binomial regression)ってのはある。もちろん,一般化線形モデルの一種。
これをするには,MASSパッケージのglm.nb関数を使う。
基本的には,
library(MASS)
m<-glm.nb(d1~d2)
summary(m)
これでいい。基本的には,lmやglm関数と同じ使い方,