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でソースを貼付ける良い方法、何かないですか?

2009年7月21日火曜日

夏休み突入

上半期最後のゼミ発表、毎年恒例のB4ゼミ合宿と、ゼミ関連のイベントが目白押しで忙しく過ごしていました。
僕が指導を担当した後輩達の発表では、自分が発表するときよりも緊張。
きっと僕の上司もそう思っているのでしょう(いつもありがとうございます)。

指導担当といっても実質的には助手の方に任せっきりでありました。
より一層精進せねば。


ところで、今日でひとしきりイベントに片がつきました。
学生として最後の夏休みの始まりです。
同期に夏の予定を聞いたところ、「バイトと研究」だそうです。
僕はとりあえず「論文執筆と研究」です。
普段と変わらない...(普段より研究するかも...)

2009年7月8日水曜日

Mozila Open Web Tools Directory

こちらの記事より

FireFoxで有名なMozilaがWeb上のツールを集めてくれています。


ほとんどシラナイ・・・。
開発者とは到底呼べない僕ですが、いくつかを自在に操り面白いものを作れたらいいなと思ってます。

2009年7月1日水曜日

イタリア所感

忘れないうちにイタリアの思い出を。
強行軍の合間を縫って、1泊2日でローマ観光。
街中に歴史的・芸術的建造物が点在する、贅沢な街でした。

ローマの街並

サン・ピエトロ寺院

コロッセオ

滞在中に感じたのは

 ・所謂イタリア人のイメージ(陽気、好色 etc...)どおり
 ・ローマもこの時期暑いよ...
 ・ローマにも三越あるよ...
 ・あぁ、1日かけてバチカン美術館に浸りたい
 ・「インターネット使えます」って書いてあったじゃんかよ!
 ・スーパーはなんでも格安だよね
 ・ローマ人よりもトリノ人の方が日本人に優しいなぁ

なんてところ。
ところで、イタリアはなんといっても飯がうまい!
新鮮な魚介をメインに、ピザ、ラザニア、パスタ、ジェラート etc...を堪能してきました。
10日近くの滞在で、はずれは1回だけ!


ムール貝のトマト蒸し

もう、色々


でもやっぱりメインは学会だったかな。
家具なんて到底買いに行く暇もなく、手に入れたのはこいつだけ

A di ALESSI ボトルオープナー

ALESSIのボトルオープナーです。
1800円くらいで購入できたので、2割引くらい?
機会を作って必ずもう一回行こうと思います。