草薙の研究ログ

英語の先生をやってます。

幾何分布関係メモ: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関数と同じ使い方,