2009年7月24日金曜日

ベータ分布の最尤推定(Rで実装)

研究の息抜きに実装してみました。
参考文献は有名なこちらです。

 Estimating a Dirichlet distribution (T.Minka, 2003)

二項分布の共役事前分布であるベータ分布から、n個の確率が生成されたとします。
この確率ベクトルからベータ分布のパラメータaとbを最尤推定にて求める関数が以下です。


newton_beta<-function(x,shape1=1.0,shape2=1.0,eps=0.0001){

 euqlid<-function(y){
  return(sqrt(sum(y*y)))
 }

 N <- length(x)
 alpha <- c(shape1,shape2)
 pn <- 0

 while(TRUE){
  palpha <- alpha
  #gradient
  g <- c(N * digamma(sum(alpha)) - N * digamma(alpha[1]) + sum(log(x)) , N * digamma(sum(alpha)) - N * digamma(alpha[2]) + sum(log(1 - x)))
  #Parameters for Hessian
  q <- c(- N * trigamma(alpha[1]) , - N * trigamma(alpha[2]))
  z <- N * trigamma(sum(alpha))

  #Update
  b <- sum(g/q) / ( 1/z + sum(1/q))
  alpha <- palpha - (g - b)/q

  n <- euqlid(alpha)
  print(alpha)
  print(n)

  if( abs(pn - n)/pn < eps ){
   break
  }
  pn <- n
 }

 return(list("shape1"=alpha[1],"shpae2"=alpha[2]))

}


厳密には閉じた形では解けないので、ニュートン法を用いて求めます。
ベータ分布(ディリクレ分布)で効率的に計算を行うための工夫をしているところがポイントのようです。
百聞は一見にしかずということで、ためしにやってみましょう。
まず、適当なベータ分布から確率ベクトルを作ります。


a <- 3
b <- 5
x <- rbeta(100,a,b)


真のパラメータがのベータ分布は以下のようになります。



発生させた確率ベクトルよりパラメータを推定するとこんな感じで計算します。


newton_beta(x)
[1] 1.377578 1.730856
[1] 2.212144
[1] 1.916810 2.800353
[1] 3.393544
[1] 2.487916 3.960042
[1] 4.676715
[1] 2.833155 4.670363
[1] 5.462514
[1] 2.907663 4.824839
[1] 5.633256
[1] 2.910250 4.830231
[1] 5.63921
[1] 2.910253 4.830237
[1] 5.639217
$shape1
[1] 2.910253

$shpae2
[1] 4.830237



これは推定されたパラメータがであることを表しています。概ねよさげですが、ちょっと精度が悪い?
このパラメータを持つベータ分布は以下のようになります。


初期値を変えると結果も変わるかもしれません。ベータ分布のように、二つのパラメータを初期値にとることができます。


newton_beta(x,初期値1,初期値2)



間違い等ございましたらご指摘ください。
ところで、Bloggerでソースを貼付ける良い方法、何かないですか?

2 件のコメント:

TAKA さんのコメント...

RSSフィードを整理していたら久しぶりに見つけたのでついでに書き込みします。
といってもベータ分布なんぞわかるはずもなく…
Macintoshについて何か書いてもらえるとうれしいです。


追伸:卒論のやる気が薄くなってきて停滞気味です。

Unknown さんのコメント...

そろそろ授業開始?
Macではないけど更新しておいたので是非。
メール、確認したら返信よろしく。

コメントを投稿