2009年11月28日土曜日

Probabilistic Latent Semantic Analysis : PLSA (Rで実装)

前回のエントリからはや一ヶ月。月日が立つのは早いものです。
修論に向け、bag-of-featuresの実装をもくろんでおりますが、その一環としてPLSAを試してみました。

参考文献はこちら(リンク先pdf)。

T. Hofmann. Probabilistic latent semantic analysis. In Proceedings of the 15th Conference on Uncertainty in AI, 1999.

ちょうど10年前に提案されたモデルですが、LDAの元となったり、現在でも多くの論文が発表されたりと、良い言語モデルのようです。
これをRで素直に実装したのがこちら。
plsi <- function(x, K=10, eps=0.9, max_itr=200,...){
 #logsumexp
 logsumexp<-function(x,y,flg){
  if(flg){
   return(y)
  }
  if(x == y){
   return(x + 0.69314718055)
  }

  vmin = min(x,y)
  vmax = max(x,y)

  if(vmax > (vmin + 50.0)){
   return(vmax)
  }else{
   return(vmax + log(exp(vmin - vmax) + 1.0))
  }
 }

 #random variable from Dirichlet distribution
 rdirichlet<-function(alpha,prec=0,...){
  if(prec != 0){
   x<-rgamma(length(alpha),alpha*prec,1.0)
  }else{
   x<-rgamma(length(alpha),alpha,1.0)
  }

  return (x/sum(x))
 }

 #initialize
 D <- dim(x)[1]
 W <- dim(x)[2]
 T <- 1 #temperature
 pllik <- 0

 #parameters
 pz_dw <- list()
 for(k in 1:K){
  pz_dw[[k]] <- matrix(0, D, W)
 }
 pz <- rep(1/K, length=K)
 pd_z <- matrix(0, K, D)
 pw_z <- matrix(0, K, W)
 unigram <- 1/W * apply(x,2,sum)
 for(k in 1:K){
  pd_z[k,] <- rdirichlet(rep(1/D, length=D))
  pw_z[k,] <- rdirichlet(unigram)
 }

 ##############################

 for(t in 1:max_itr){
  
  cat("Iteration : ",t," Temperature : ",T,"  ")

  #E-step 
  #initialize expected count
  cz <- rep(0, length=K)
  cd_z <- matrix(0, K, D)
  cw_z <- matrix(0, K, W)

  for(d in 1:D){
   for(w in 1:W){
    z <- 0
    for(k in 1:K){
     pz_dw[[k]][d,w] <- T * ( log(pz[k]) + log(pd_z[k,d]) + log(pw_z[k,w]) )
     z <- logsumexp(z, pz_dw[[k]][d,w], (k==1))
    }
    
    for(k in 1:K){
     pz_dw[[k]][d,w] <- exp(pz_dw[[k]][d,w] - z)

     if(x[d,w] != 0){
      cz[k] <- cz[k] + x[d,w] * pz_dw[[k]][d,w]
      cd_z[k,d] <- cd_z[k,d] + x[d,w] * pz_dw[[k]][d,w]
      cw_z[k,w] <- cw_z[k,w] + x[d,w] * pz_dw[[k]][d,w]
     }
    }

   } 
  }

  #M-step
  pz <- cz / sum(cz)
  for(k in 1:K){
   pd_z[k, ] <- cd_z[k, ] / sum(cd_z[k, ])
   pw_z[k, ] <- cw_z[k, ] / sum(cw_z[k, ])
  }

  #converged?
  llik <- 0
  for(d in 1:D){
   for(w in 1:W){
    if(x[d,w] != 0){
     l <- 0
     pdw <- 0
     for(k in 1:K){
      l <- log(pz[k]) + log(pd_z[k,d]) + log(pw_z[k,w])
      pdw <- logsumexp(pdw, l, (k==1))
     }

     llik <- llik + x[d,w] * pdw
    }
   }
  }
  cat(" log-likelihood : ",llik, " PPL : ",- 1/W * llik,"\n")

  if(t > 1){
   if( abs((llik - pllik) / llik) < 1e-5 || pllik > llik){
    cat("Converged.\n")
    break
   }
  }
  pllik <- llik

  T <- eps * T
 }

 return(list("pz_dw"=pz_dw, "pz"=pz, "pd_z"=pd_z, "pw_z"=pw_z))
}

素直に参考文献の式(3)~(6)を実装したものになります。推定アルゴリズムはTemperature EMアルゴリズムで、温度スケジューリングは簡単のため
  1. 温度T=1.0からスタート
  2. EMアルゴリズムの各反復において, T = eps × T (eps < 1.0)で更新
としています。論より証拠。サンプルで試してみましょう。
用いたデータは以下のような12×9の行列データです。
1 2 2 0 0 0 0 0 0
2 1 3 0 0 0 0 0 0
2 0 2 0 0 0 0 0 0
1 1 1 0 0 0 0 0 0
0 0 0 1 1 2 0 0 0
0 0 0 2 1 1 0 0 0
0 0 0 1 3 1 0 0 0
0 0 0 1 2 2 0 0 0
0 0 0 0 0 0 2 3 3
0 0 0 0 0 0 3 2 2
0 0 0 0 0 0 2 1 3
0 0 0 0 0 0 1 1 4
これがテキストデータであるとすると、データ数が12、単語数が9、行列要素はデータ内で観測された頻度に相当します。
目視で確認できるように、左側・中央・右側に出現頻度が偏った、いわば3クラスが含まれるデータになります。

これを、以下のようにしてplsiを実行します。
> opt <- plsi(x, K=3, eps=0.95)
オプションKはクラス数、epsはTEMの温度に関するパラメータです。これにより、以下のような出力がされます。
Iteration :  1  Temperature :  1    log-likelihood :  -277.8714 
Iteration :  2  Temperature :  0.9    log-likelihood :  -275.7771 
Iteration :  3  Temperature :  0.81    log-likelihood :  -274.7316 
Iteration :  4  Temperature :  0.729    log-likelihood :  -272.8197 
Iteration :  5  Temperature :  0.6561    log-likelihood :  -268.9356 
Iteration :  6  Temperature :  0.59049    log-likelihood :  -265.5543 
Iteration :  7  Temperature :  0.531441    log-likelihood :  -264.1783 
Iteration :  8  Temperature :  0.4782969    log-likelihood :  -264.3499 
Converged.
最後の反復で若干対数尤度が悪化しているのは玉に傷です。

各クラスからの文章の出現確率、すなわちp(d|z)は、以下によって取得できます(見やすさのため小数点以下二桁で丸めています)。
> round(opt$pd_z,2)
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
[1,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00  0.3  0.26  0.22  0.22
[2,] 0.28 0.33 0.22 0.17 0.00 0.00 0.00 0.00  0.0  0.00  0.00  0.00
[3,] 0.00 0.01 0.00 0.00 0.22 0.22 0.27 0.27  0.0  0.00  0.00  0.00
行がクラス、列がデータを表しています。概ね正しく3つのクラスに分類できていることが伺えます。
同様にして、各クラスからの単語の出現確率p(w|z)は、以下によって取得できます。
> round(opt$pw_z,2)
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
[1,] 0.00 0.00 0.00 0.00 0.00 0.00  0.3 0.26 0.44
[2,] 0.33 0.22 0.44 0.00 0.00 0.00  0.0 0.00 0.00
[3,] 0.01 0.00 0.01 0.27 0.38 0.33  0.0 0.00 0.00
行がクラス、列が単語を表しています。こちらも概ね正しく各クラスからの単語の出現確率を推定できていることが伺えます。


初期値やK/epsの値によってかなり推定精度にバラつきがありますが、シンプルな理論/実装で動くモデルなので良いですね。例によって、間違い、またより良い実装方法がございましたらご指摘ください。

2009年10月29日木曜日

Tweetup with @syou6162

本日、素敵なTweetup@syou6162さんとご一緒しました。
Web上でも、現実でも、精力的に活動されている姿は尊敬するばかりです。

ところで、日本国籍でない友人に日本人の悪い点を聞くと、すぐにtoo shyと返ってきました。
国民性や個人の性格的なものは無理に変える必要はないと考えていますが、それでもコミュニティに対して何らかのコミットをするのは、コミュニティにとっても、自身にとっても有益なものとなる場合が多いように思います。

しかしながら、そういった行動をおこすには、またそういった方々とお近づきになるには、それなりの労力が必要です。
Twitterに限らず、SNSやWebは、それらの敷居を下げる良い道具ではないかと思います。


抽象的ではありましたが、非常にポジティブな気持ちを想起させる、そんなTweetupでした。
より詳細な会話の内容が知りたい方、また、R, C++, Emacs etc...の良質なtipsをご覧になりたい方は、こちらのサイトをどうぞ。

 Seeking for my unique color.
 Seeking for my unique color.:@TakaakiTalkさんとお食事

Thanks a lot @syou6162!

2009年10月27日火曜日

RとImageMagickで混合分布推定アニメーション

EMアルゴリズムや変分ベイズ、MCMCでは最適化やサンプリングを繰り返すことで確率分布を求めます。
その過程を可視化する際のtipsをまとめてみます。

Rの他に必要なのはこちら。


MacPortでしたら

sudo port install ImageMagick +lcms +jpeg2
sudo port install ffmpeg

でインストール可能です。
動画作成までの道のりは以下のとおりです。

 1. Rプログラムの反復演算の都度画像を生成する
 2. ImageMagickで連結・mpeg変換を行う

いやぁ、簡単。このままではtipsにならないので、サンプルとその生成方法を記載してみます。
まずRにおいて、以下のような画像生成項目を記載します。

name = paste(t, ".png" , sep="")
png(file = paste("img/", name, sep=""), sep="")
 #プロットしたい図のplot
dev.off()

尚、tは反復回数t回目を表します。これによりimgフォルダにt=1から(任意の)収束回数までの連番pngファイルが生成されます。
続いて、imgフォルダ内に以下のようなシェルスクリプトを保存します。

#!/bin/bash
num=${1:?"Usage: $./conv.sh npics"}
for((i=1;i<$(($1+1));++i))
do
 name+=$i.png" "
done
#echo $name
`convert -delay 15 -quality 100% -antialias -compress None $name opt.mpg`

名前をconv.shとしましょう。また、t=100で収束したとします。その際、実行権限を与えた後、以下のように実行します。

chmod u+x conv.sh
./conv.sh 100

これにより、opt.mpgができあがります。表示間隔はdelayの引数を変えることで調整できます。
混合ガウス分布の推定を行った結果がこちら

ImageMagick自体は検索すれば多くのtipsが出てきますが、Rのtipsとして参考にしたのはこちら。

 RjpWiki:グラフィックス参考実例集:自作グラフィックス投稿欄:GIFアニメーション

一枚一枚画像を生成するので非常に時間がかかります(MCMCでやることは考えたくないです)が、視覚的に確認できると楽しいですね。是非お試しあれ。

2009年10月15日木曜日

統計のための線形代数 in C++ (Boost uBLASの参考サイトとtipsまとめ)

最近の僕はより大規模な実験を行うべく、これまで作ってきたtoy programをC++に移植する毎日です。
その過程でどうしてもベクトル・行列演算が必要になったので、時期標準の呼び声高いBoost内のベクトル・行列ライブラリ、uBLASを用いてみました。
今回の記事は、特に僕が必要であった演算を実現すべく、参考にしたサイト及びそれらをまとめたテンプレート関数群をご紹介します。

2009年9月24日木曜日

Snow Leopard導入記録その2 その後

つい先日、高校の部活仲間2人が婚約したとメールをくれました。
本当におめでとう。
前回のエントリ然り、ここへきて僕の周りでは色々変化があるようです。

さて、GoogleAnalyticsでこのblogの動向を見ていると、Snow LeopardにおけるTex環境整備について書いたエントリの注目度が高いです。
おおまかに書いたため対応しきれていない可能性が高いですが、お気づきの点等あればコメントいただければ幸いです。


では、iMac(Mid 07)にインストールしたSnow Leopardのその後を少し。

 ・MacPortで入れたlftpが動かない
   →upgradeを試みるとncurseswとncursesに"+darwin_10"が
    matchしないよと警告される
   →ncurseswとncursesを再インストール
   ⇒lftpを再インストールでfix
  
 ・e-mobileが動かないらしい
  参考:理系学生日記:E-mobile のモデムが Snow Leopard で
     使用できなくなったあなたへ
   ⇒32bitモードで起動とのこと

 ・Cyberduckが起動しない
   ⇒対応バージョンを再インストール

 ・Preview.appでTexが壊れる
  参考:Snow LeopardのプレビューだとLaTeXで書いた数式が
     こんな風に見えます。
   ⇒とりあえずAdobe Acrobat Readerで開く
 
 ・TimeMachineがLeopard時に比べ非常に早い!
 ・劇的ではないが以前より動作が軽い


僕の以前の環境にも依存するのでなんとも言えませんが、メインで研究等に使用しているMacへのインストールは未だ不安かもしれません。
個人的にはiPhone OS 3.1へのアップグレードの方が恩恵を大きく感じました。
これは速度改善に依るものが大きかったため、ラップトップでこそSnow Leopardが生きてくるのかもしれません。
環境整備に時間がかかるやもしれませんが、自身で問題なく対応できる、もしくは十分に時間のとれる方は検討してみては。

2009年9月8日火曜日

帰省をしておりました

一週間ばかりの実家帰省から舞い戻ってきました。
今回の帰省のベストショットはこれ。

父親デビュー!

なんて使い古された冗談さておき、同い歳の友人の娘です。
僕の実家は地方であることもあり、就職、結婚、出産が早いのですが、今でも連絡を取り合う友人では初の結婚&出産。
それはもう、感動一塩でした。

学生として最後の夏休みなので、他にも多数の幼なじみや友人達と再開。
話の大半が思い出話であることに年月を感じつつも、笑いのツボは当時のままで、こちらも感動でした。

多いにリフレッシュできたので、今日からまた頑張ります。

2009年8月30日日曜日

Snow Leopard導入記録その1 MacPort関連簡易メモ

Pre-orderの甲斐あって、きっかり発売日(8月28日)に到着致しました。
MacOS 10.6 Snow Leopardです。

Snow Leopard has com! on Twitpic
早速メインのiMacに

 ・バックアップはTimeCapsule任せ
 ・上書き(通常)インストール

を慣行。時間がなかったので強行軍でした。
で、Apple純正ソフト以外で通常どおり動いたもの(○)、動きそうなもの(△)、動かなかったもの(×)はこんな感じ

 ・Chrome 4.*** ◎
 ・Firefox 3.*** ◎
 ・R 2.8-2.9.2 ◎
 ・Adobe CS3 Illustrator Photoshop ○
 ・TexShop *** ×(pTexに問題あり)
 ・LatexI 1.16*** ×(pTexに問題あり)
 ・MacPort *** ×
 ・Parallels Desktop for Mac 3.*** ×

MacPort、Parallelsはともかく、Tex関連はそのまま動いてて欲しかった...

Tex関連のMacPortによるSnow Leopard環境構築(簡易表記版)は以下

 1. MacPortの再インストール
  こちらを参考に

 2. Tex環境の再構築(E:error, S:solution)

  E1 とりあえずsudo port upgrade outdated
    →fontconfigのアップデートでXMLがおかしいよと怒られる

  S1 XML関連ライブラリ(libxml?)のアップデート
    →sudo port upgrade fontconfigでfix

  E2 再度sudo port upgrade outdated
    →libpngのアップデートでzlibが見つからないよ!と怒られる

  S2 こちらを参考にzlibを+universalで再インストール
    さらにsudo port edit libpngでconfigureの引数に--libdir=
    /opt/local/libを加える(不要?)
    ⇒sudo port upgrade libpngでfix

  E3 再度sudo port upgrade outdated
    →pTexのアップデートでkpathsea.aがないよと怒られる

  S3 sudo port edit pTexでconfigureの引数に--libdir=
    /opt/local/libを加える
    ⇒sudo port upgrade pTexでfix

 3. こちらからTexShopをダウンロード・アップデート
 
でうまくいきまいした。要点をまとめると

  • インストールしてあるライブラリがおかしいと怒られたら
  • とりあえず怒られたライブラリ側をアップデート/再インストール
  • それでもダメなら
  • +universalをつけて怒られたライブラリをアップデート/再インストール
  • インストールしてあるライブラリがないよと言われたら
  • sudo port edit パッケージ名 でライブラリの場所を追加して(configure.argsに書いて)やる
の3点です。パッケージソフトの利点は何処へ...
詳細なエラー内容、アップデート前後のバージョン、対応等、実機から離れてしまっているため不備が多いです。参考程度として下さい。
何かあればtwitterでお気軽に。

2009年8月26日水曜日

真の夏休み突入

すっかりご無沙汰になりました。
特に山梨に潜伏しているわけではございません。

7月〜8月は通常時より大学(研究室)に顔を出しておりました。
細かい調整はさておき、協力させていただいていた論文がやっと一段落です。
結局以前のエントリのとおりに事が進んでいるということですね。
良かった良かった。

今後は

 ・9月〜10月にTechnical Reportを報告
 ・12月迄に後輩と学会
 ・3月までに論文投稿

なんてできたらいいなと思っています。
積み重ねているものがあるわけではないので、このうち1〜2つ実践できれば万々歳でしょう。


さて、Blogを書かないうちに色々面白げなガジェットの噂/詳細が出ていました。


Leaked Mac Tablet?



Sony Ericsson Rachael UI



Nokia Booklet 3G first video




最初の動画はfakeですが、よくできているなぁ。
iPhoneを持ち始めてから感じたのですが、良いUIを持つデバイスは生活を変えますね。
twitterとyoutubeはiPhoneから利用する方が便利かもしれません。

そういった意味ではNokiaのNetbookは特にUIにこだわっているわけではありませんが、世界で実績のあるメーカがNetbookに進出ということで、期待を込めて。12hバッテリっていいな!Ovi servicesってどうなの?


今週末からは息抜きとして、興味のあることだけをする1週間にします。

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割引くらい?
機会を作って必ずもう一回行こうと思います。

2009年6月28日日曜日

無事帰国

イタリアより本日無事帰国。
これに行ってました。

 http://bnpworkshop.carloalberto.org/

研究分野の大物達がひしめく学会で、芸能人ばかりに遭遇した感覚。
なんて幸せ!
ただインターネット環境が貧弱、かつタイトなスケジュールであったことにより、現地からBlogを更新できなかったことが心残り。

学会会場

頑張る後輩

頑張る後輩?


今回の研究成果は先人達の偉業があってこそのもの。
次は完全にオリジナルで出せるよう精進しよう。
イタリア所感はまた今度!

2009年6月14日日曜日

出発前に色々

いよいよ来週末からイタリアへ。
ポスターは問題ないが、論文の再投稿に一枚噛むことになり、これから急いで実験。
理系学生らしい忙しさにアップアップ。


前回の海外は中国だった。昨年の10月頃。
4泊5日の日程だったが、デニム・Tシャツ2枚・パーカ+洗濯で乗り切った。
今回は8泊10日ということもあり、もう少しボリュームを増やしたいところ。
現地調達でもいいかななんて考えている。

ところで今回、Tシャツの他にも現地で調達したいものがある。


DANESE PinUp
 


イタリアが誇る家具ブランド、DANESEの壁掛けボード。
良く行くカフェに置いてあって一目惚れ。
日本で購入するのに比べ、3〜4割は安く購入できると睨んでいる。
問題は現地での知名度(取り扱いの多さ)とそのサイズ。
ミラノを本拠地としているが、残念ながら行く時間はない。
経由地のローマや、トリノにあるかどうか。
また、例え購入できたとしても持って帰るのに一苦労。
うーん。

2009年6月8日月曜日

変分ベイズでディリクレ過程混合モデル(Blei論文を解く)

DPMに対する変分ベイズ法を書いたBleiの論文、Variational inference for Dirichlet Process mixturesを解いてみました。
ざっと書き上げたまま確認をしておりませんので、間違いがあればご指摘ください。
式展開だけ妙に丁寧にしていますが、全体的に分かりにくい文だなぁ...

 https://www.box.net/shared/2e6k7abce5

Google DocumentsがPDF Web共有できればいいのに!

2009年6月1日月曜日

線形代数事始め

統計をたしなむ以上避けて通れない線形代数。
高校依頼避け続けてきたけどそろそろ限界だなぁ。

 統計のための行列代数(参考書)
 統計のための行列代数 上 (1) 統計のための行列代数 上 (1)

 Invertible matrix:Wikipedia
 http://en.wikipedia.org/wiki/Invertible_matrix

 Woodbury matrix identity:Wikipedia
 http://en.wikipedia.org/wiki/Woodbury_matrix_identity

米国の学生は上記参考書を授業で使うらしい。
まずいぞ日本の学生!

2009年5月31日日曜日

Japanimation

僕は夢中になると疲れを感じない体質なのですが、その分遊びに集中すると仕事は多いに後回しです。
このエントリも後回しにした学会ポスター作成の後に書いています。

本日の遊びはこちら。

FLCL Fooly Cooly

手塚治虫,火の鳥


Japanimation万歳!!
片やエヴァで一世を風靡したGAINAX制作アニメ、片や巨匠手塚治虫のLifeworkと、新旧の名作を堪能(?)致しました。
FLCLは動画サイト、火の鳥は漫画を借りるとエコな遊びでしたが、決めました。どちらも買います。


どちらも古いため新古・中古がよさげ。

2009年5月29日金曜日

Native Google Chrome for Mac OSX

Windows版の軽快な動作から、Mac版にも期待しているGoogle Chromeですが、Nativeの開発版が既にあるようですね。

 DOWNLOAD UPDATED NATIVE GOOGLE CHROME FOR MAC OS X

早速ダウンロードしてインストール...

 おぉ 動いてる

やっぱり軽い・速いように思います。
他にもMac上の仮想環境(?)のひとつ、Cross Overで動くCross Over Chromiumなんてのがあるようです。

お試しあれ。

2009年5月27日水曜日

海外学会ポスター作成メモ

来月のWorkshopに向けてポスター作成を開始。
以下やったこと。

 1. MS PowerPointで作成開始
  →A0サイズがうまくpdf化されない
  →Macじゃ見れないじゃないか
  ⇒Illustrator作成に変更

 2. Illustratorで作成開始
  →フォントに悩む

 3. Illustrator CS3にTexのフォントを入れてみる
  →Texの数式をDTPソフトに
  →BoldならBoldフォント、ItalicならItalicフォントを利用する…
  →なんだかスマートじゃない
  ⇒却下

 4. 一般的なフォントを利用する
  →What are the best fonts to use for a presentation?
  →Wikipedia:Sans-serif
  →Arialは好きじゃない
  ⇒Sans-selfじゃないけどTimes New Romanでいっかー

 5. ポスター下書き.pdfを先生に送る
  →いやpptで送ってよ
  ⇒1に戻る

作成の過程でPowerPointにepsを貼付けたくなりました。
OSX CS3からは

 別名で保存
  →フォーマット:Illustrator eps
  →EPSオプション:
    バージョン:Illustrator 10 eps
    プレビュー:TIFF 8ビットカラー

でEPSを保存しておいてやれば問題なく貼付けできました。
海外ではMacユーザは少ないからPowerPointで書いた方が良いなんて記事を読んだけど、本当でしょうか?

2009年5月25日月曜日

BloggerはFireFoxが便利?

最近は専らSafariを使っておりますが、Bloggerの編集だけはFirefoxが良いようで。
いかんせんSafariだと画像のアップロードができなかったりします。

また、かの強力アドオンGreasemonkeyでLatexまで!

 クリボウのBlogger Tips
 Blogger に LaTeX 数式を挿入するユーザースクリプト
 http://www.kuribo.info/2009/05/blogger-latex-latex-for-blogger.html

試しに...

 

ベイズの定理だってこのとおり!

昨夜は実家の両親とskypeでビデオ通話をしてみたり。
オンライン環境の進歩は目を見張るものがありますね。

2009年5月24日日曜日

ガウス過程事始め

夏のゼミに向けて、ガウス過程(Gaussian Process)の勉強を開始。
資料のメモ。

 Mackay, Introduction to Gaussian Process
 http://www.inference.phy.cam.ac.uk/mackay/gpB.pdf

 正田備也, ガウス過程に関するメモ(1)
 http://www.iris.dti.ne.jp/~tmasada/2007071101.pdf

 CE Rasmussen, Gaussian Processes for Machine Learning
 (Adaptive Computation and Machine Learning)
 


普段はひたすらGoogle Scholarに頼っているのだけれど、やっぱり本はいいですね。
研究室にもあるのだけれど、自費購入を検討中。

2009年5月20日水曜日

nkf拡張シェルスクリプト

普段はMac、外や研究室ではWindows、実験ではLinuxと環境を使い分けていると、文字コードに悩まされることがある。
そんなときの強い味方、Network Kanji Filter : nkfであるが、複数のファイルを一度に変換するにはいちいち面倒くさい。
どうせならフォルダ内の対応したファイルを全て変換したいっ!

そんな煩わしさを解消するシェルスクリプトを書いてみた。
ソースは続きから。


「あるディレクトリのC++ファイル(.*\.hや.*\.cpp)を全てunicodeに変換したい!」場合の使い方は以下。

 nkf.sh -w -d -r ディレクトリ名

オプション-dでディレクトリ指定を。-rで元のファイルの削除を。-w/-sでutf-8/sjisの変換を行う。
もちろん単一ファイルも対応可能で、その際は-dオプションを外せばよい。

 nkf.sh -w -r ファイル名

.*\.h、.*\.cpp、.*\.txtファイルを対象としているが、適宜ソースを編集すればOK。
nkfのオプションも適宜編集のこと。
ソースはこんな感じ。もっと短く書けるだろうなぁ。



#!/bin/sh
#extended nkf
#Author:Takaaki
#Last update:09/05/20

is_rm=0 # whether remove original files or not defalt:false
is_dir=0 # whether $1 is a directory or not defalt:false
nkf_opt="" # format for nkf
while getopts "rdsw" opt; do

case $opt in
r ) is_rm=1 ;;
d ) is_dir=1 ;;
s ) nkf_opt="-s" ;;
w ) nkf_opt="-w" ;;
\? ) echo 'usage: nkf -format [-r] [-d] args'
exit 1
esac

done

shift $(($OPTIND - 1))

if [ $nkf_opt = ""];then
 echo 'usage: -format must be given (-w / -s)'
 exit 1
fi

if [ $is_dir = 1 ];then
 for file in $(ls $1);do
  ext=`echo $file | sed -e "s/.*\.\([^\.]*\)\$/\1/g"`
  if [ "$ext" = "cpp" -o "$ext" = "h" -o "$ext" = "txt" ];then
   if [ $is_rm = 1 ];then
    nkf $nkf_opt $1/$file > $1/$file~
    mv $1/$file~ $1/$file
   else
    filename=`echo $file | sed -e "s/\.$ext//g"`
    nkf $nkf_opt $1/$file > $1/$filename~.$ext
   fi
  fi
 done
else
 ext=`echo $1 | sed -e "s/.*\.\([^\.]*\)\$/\1/g"`
 if [ "$ext" = "cpp" -o "$ext" = "h" -o "$ext" = "txt" ];then
  if [ $is_rm = 1 ];then
   nkf $nkf_opt $1 > $1~
   mv $1~ $1
  else
   filename=`echo $1 | sed -e "s/\.$ext//g"`
   nkf $nkf_opt $1 > $filename~.$ext
  fi
 fi
fi



2009年5月18日月曜日

ドクターフィッシュ

週末を利用して母校の生徒の卒論指導のため帰省。
父と出かけた日帰り温泉で初体験。


旅ぃゆけぇばぁ・・・



ひぃあぁぁぁ 喰われとるやないか



魚による不規則な吸引は、くすぐったさを余り感じない僕をもってしてでもなかなかのもの。
気持ちいいというより、面白い体験でした。

ところで、たっぷり15分間の吸引の最中、役目を終えた後の魚達のことばかりを考えていました。
食すわけにもいかないしなぁ。

2009年5月12日火曜日

変分ベイズで混合正規分布推定(Rで実装)

以前のエントリの先駆けとして簡単な実装を。
初期値と推定結果はこんな感じになる。

初期値と

推定結果

参考文献は上田さんの学会誌

 上田修功,ベイズ学習[IV・完] : 変分ベイズ学習の応用例

基本的にこのとおりに実装すれば動く。素晴らしすぎる。
他に参考文献はこの辺り

 荒木佑季,変分ベイズ学習による混合正規分布推定

ソースは続きから。

Rでのソースはこんな感じ。
2変量正規分布での実装である。
基本的に参考文献通りに計算をしているだけ。
参考文献のnotationでΣが分散共分散行列でない点に注意。

ちなみに、plot.gaussは2次元プロット上に正規分布をplotするための関数。
パッケージにもあった気がするが使い勝手が悪く実装した。


vb<-function(x,K=10,maxItr=100,cirCol="red",...){
source("plot.gauss.r")

#initialize
N <- dim(x)[1]
d <- dim(x)[2]
phi_0 <- N / K
xi_0 <- 1.0
eta_0 <- d + 2
nu_0 <- c(mean(x[,1]),mean(x[,2]))
SS_0 <- matrix(c(var(x[,1]),0,0,var(x[,2])),2,2)
B_0 <- SS_0

x_h<-cbind(rep(0,K),rep(0,K)) #R^K*2
N_h <- rep(N/K,length=K)
phi <- rep(phi_0,length=K)
eta <- rep(eta_0,length=K)
p<-rep(1/K,length=K)
f <- eta_0 + N_h + 1 - d #R^K
mu <- list()
B <- list()
sigma <- list()
C_h <- list()
S <- list() #precision matrix
SS <- list() #variance covariance matrix
for(k in 1:K){
mu[[k]]<-c(rnorm(1,nu_0[1],2),rnorm(1,nu_0[2],2))
B[[k]]<-B_0
sigma[[k]] <- B[[k]]/(f[k] * (N_h[k] + xi_0))
C_h[[k]]<-matrix(0,2,2)
S[[k]] <- (eta_0 + N_h[k]) * solve(B[[k]])
SS[[k]] <- solve(S[[k]])
}
#print(list("N_h"=N_h,"phi"=phi,"mu"=mu,"eta"=eta,"f"=f,"B"=B,"sigma"=sigma,"S"=S,"SS"=SS))

#plot input data
#x11()
plot(x,main="Training data & initial points",...)
plot.gauss(p,mu,SS,add=T,col=cirCol,xlab="",ylab="",xaxt="n",yaxt="n",...)



#########################################################

for(t in 1:maxItr){
cat("Iteration ",t,"\n")
#VB-Estep
A0 <- digamma(phi_0 + N_h) - digamma(K*phi_0 + sum(N_h))
A1 <-0
for(j in 1:d){
A1 <- A1 + digamma((eta_0+N_h+1-j)/2)
}

gamma<-matrix(0,N,K)
for(k in 1:K){
for(n in 1:N){
gamma[n,k] <- A0[k] + 1/2 * A1[k] - 1/2 * log(det(B[[k]])) - 1/2 * sum(diag( (eta_0 + N_h[k]) * solve(B[[k]]) %*% ( f[k]/(f[k]-2) * sigma[[k]] + matrix(x[n,]-mu[[k]]) %*% t(matrix(x[n,]-mu[[k]])) ) ))
}
}
z0<-exp(gamma)
z1<-rep(0,N)
for(n in 1:N){z1[n]<-sum(z0[n,])}
z_h<-z0/z1 #R^N*K

#VB-Mstep
for(k in 1:K){
#N_h and x_h must calc at first
N_h[k] <- sum(z_h[,k]) #R^K
x_h[k,] <- c((sum(z_h[,k]*x[,1]) / N_h[k]) , (sum(z_h[,k]*x[,2]) / N_h[k])) #R^K*2

C_h[[k]]<-matrix(0,2,2)
for(n in 1:N){
C_h[[k]] <- C_h[[k]] + z_h[n,k] * matrix(x[n,] - x_h[k,]) %*% t(matrix(x[n,] - x_h[k,]))
}

mu[[k]] <- (N_h[k] * x_h[k,] + xi_0 * nu_0) / (N_h[k] + xi_0)
B[[k]] <- B_0 + C_h[[k]] + (N_h[k] * xi_0)/(N_h[k] + xi_0) * matrix(x_h[k,] - nu_0) %*% t(matrix(x_h[k,] - nu_0))
}

phi <- phi_0 + N_h
eta <- eta_0 + N_h
f <- eta + 1 - d

for(k in 1:K){
sigma[[k]] <- B[[k]]/(f[k] * (N_h[k] + xi_0))
S[[k]] <- (eta_0 + N_h[k]) * solve(B[[k]])
SS[[k]] <- solve(S[[k]])
}
#print(list("N_h"=N_h,"x_h"=x_h,"phi"=phi,"mu"=mu,"eta"=eta,"f"=f,"C_h"=C_h,"B"=B,"sigma"=sigma,"S"=S,"SS"=SS))
}

#plot estimates
pz<-K*phi_0 + sum(N_h)
for(k in 1:K){
p[k]<-(phi_0+N_h[k])/pz
}
print(list("Mixing parameter"=p))

#quartz()
#x11()
quartz()
plot(x,main="Estimated value (VB)",...)
plot.gauss(p,mu,SS,add=T,col=cirCol,xlab="",ylab="",xaxt="n",yaxt="n",...)

return(list("N_h"=N_h,"x_h"=x_h,"phi"=phi,"mu"=mu,"eta"=eta,"f"=f,"C_h"=C_h,"B"=B,"sigma"=sigma,"S"=S,"SS"=SS,"weights"=p))
}



ちなみに、このエントリを記載中にVB for GMMのパッケージを発見

 vabayelMix

マジかよ、、、
パフォーマンスチェックをした上で、分があればちゃんと公開しよう。

2009年5月10日日曜日

物件探し その1

ゼミ発表も無事終わり、弟と物件探しへ不動産屋へ。
親戚のアパートに転がり込んで早5年。
不動産屋を訪ねて物件を探すのはこれが初めてである。

引越先では弟と完全に同居の予定。
当初の条件は2DK、2部屋は完全セパレート、バス・トイレ別、都心アクセスそこそこ良好。
Webで目星をつけた物件は残念ながら既に入居済みであったので、店舗にていくつかピックアップの後にいざ内見へ。


最初の物件は一軒家、11万5千円の破格物件。
逐年数は古く間取りも独特。風呂追い炊きの操作装置はあるが肝心の本体がない等どこかチグハグ。
しかし都心に近い一軒家物件などそうない。住めば都にしてしまえばいい。
そう心で自身に言い聞かせていた矢先、担当者が怪訝な顔で

「すみません。僕も知らなかったのですが、窓の向こうはとんでもないビューが・・・。」
「!((初回から隣が墓地物件))」

引きの強さ(?)をひしひしと感じる。
これからの季節にはもってこいだ。いつだって涼しく快眠だろう。
もしかしたら永久に目を覚まさないやもしれん。
遊びにも事欠かないだろう。卒塔婆でチャンバラする機会なんて地元でもなかった。
極めつけはその物件担当者の言葉。

「いやぁ、でもまずい話は聞かないですけどねぇ。」

聞いてたら困る。


気を取り直して2件目、3件目へ。
いずれも申し分ない物件であったが、「完全セパレート」の実現が難しい。
引越は夏までにするとして、気長に探す事にしよう。

2009年5月7日木曜日

連休明け

6月末にイタリアで開かれるNonparametric Bayes Workshopの手続きのため、1時間睡眠で家を出てつい先程帰宅。

昼過ぎでも慌ただしい銀行員を眺めつつ今日から仕事始まりであることに気づく。
この辺りは学生の特権か。

仮眠をとって明日のゼミ準備。




2009年5月6日水曜日

変分ベイズでディリクレ過程混合モデル(現状)

Bleiの論文が2004年、持橋さんが2006年だからもう最先端ではないのかな?
理由あって後2日以内にフォローしなければなりません。

現状:
実装済み・テスト未済
VBフォロー済み・DPVBフォロー未済


読んでいる論文:
Variational Inference for Dirichlet Process Mixtures

無限混合ディリクレ文書モデル


全て終わったら色々UPします。