ベータ分布の最尤推定（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でソースを貼付ける良い方法、何かないですか？