参考文献は有名なこちらです。
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でソースを貼付ける良い方法、何かないですか？
