重点サンプリング (2)

前回は正規分布の裾の積分をなんとなく決めた提案分布による重点サンプリングで求めた。今回は提案分布の違いがどのような誤差の違いを生むのかについて実験した。ただし、今回の積分範囲は [4,∞) とした。ようするに \(f\) を標準正規分布のpdfとして

\begin{align}
\int_4^\infty f(x)dx
\end{align}

を求める問題。
“重点サンプリング (2)”の続きを読む

ハイブリッドモンテカルロの実験

相関の強い二変量正規分布に対してハイブリッドモンテカルロを使ってみた。上から順に、サンプリング結果、x1の自己相関、x2の自己相関。

自己相関ほぼ完全になし、という結果になった。ギブスサンプラーだとこうはいかない。ただし、

  • 微分方程式を解く時間のスケールが小さすぎると自己相関が出たので良い感じのスケールをちょっとだけ探索した。
  • ステップ幅を固定にしたら怪しげな自己相関の挙動がでた。
  • 計算時間でギブスサンプラーと比較してどちらが有利かは今回は検討してません。

という点は追記しておきます。

“ハイブリッドモンテカルロの実験”の続きを読む

チューリングマシンと限定合理性 : 「行動ゲーム理論入門」を読んだ

この本「行動ゲーム理論入門」はたまたま本屋で見かけてパラパラと見ていたら、経済学の本にもかかわらず「チューリングマシン」だとか「強化学習」だとかいう一見経済学とは関連の薄そうな単語があったので、興味深いな、と思って脊髄反射的に購入した。僕はこの分野は全く知らない状態でこの本を読み始めたのだけど、非常に刺激的な本だったので記憶が鮮明なうちに書いておくことにする。

“チューリングマシンと限定合理性 : 「行動ゲーム理論入門」を読んだ”の続きを読む

LDA (Latent Dirichlet Allocation) 実装してみた

LDA の理解が正しいかどうかをチェックするために R で実装してみました。ついでにトピック数を真の値からずらしたときに何が起こるか観察しました。例題は論文

Finding scientific topics, Thomas L. Griffiths, and Mark Steyvers, PNAS April 6, 2004 vol. 101 no. Suppl 1 5228-5235

から。手法は同論文で提案された collapsed Gibbs sampler。

例題データ


5×5=25の単語 (w) で10のトピック (z) をもつドキュメント (d) の集合を仮定する。図の正方形は 5×5 になっていて 25 の単語の頻度(もしくは確率)に対応している。一番上の行は10のトピックそれぞれに対する単語の発生確率 \(P(w|z)\)(z=1, 2, …, 10)。その下三行は生成した各ドキュメントの単語の頻度。2000個生成したうちの最初の30個だけを示している。ハイパーパラメータはそれぞれ、\(\alpha=1, \beta=0\) (\(\beta\) は \(P(w|z)\) の事前分布(ディリクレ)のパラメータだけれども、これに関してはゼロというか、この例ではそもそも事前分布を考えていない。)

結果

▼ 真のトピック数を既知とした場合。

上から下へ向かって計算が進んでいきます。左側の縦長のグラフは尤度が反復と共にどのように変化するか。縦軸が反復回数で横軸が尤度。右側のマトリクスは左側のグラフで○があるところでの反復回数の時 (1, 2, 5, 10, 20, 30, 40, 50, 75, 100) の \(P(w|z)\) をそれぞれのトピックについて(横方向)、反復回数と共に(縦方向)示している。初期状態のランダムなトピック割り当てからだんだん真の解に収束。

▼ トピック数を多め(12個)にとった場合。

大体収束しますが、余った分のトピックはぼんやりとした感じになっている(左から三番目と一番右)。もう少し反復を増やすと変化があるかもしれない。

▼ トピック数を少なめ(8個)に取った場合。

一番左と一番右が2つのトピックがかぶったような感じになってしまっている。それでも残り6つはだいたい正しい感じになっている?こちらも反復数が足りないかもしれない。

使ったスクリプト (in R)

とにかく実装しただけで、信じられないくらい激遅なので使わないほうがいいですが、一応残しておきます。(社会的責任という意味では公開しないほうがいいかもしれませんが。まあ、それくらいのものです。)ちなみに、上の実験結果の反復数が100になっているのは実行速度が遅いためです。

■ データの生成部分。

# Settings
ndoc <- 2000   # The number of documents to be generated
nrow <- 5      # The number of rows for each document
ncol <- 5      # The number of columns for each document
nTopics <- nrow+ncol  # The number of topics
W <- nrow*ncol # The number of vocabralies
words.par.doc <- 100
alpha <- 1     # Hyper parameter for theta (topic prior probability for each document)
beta <- 0      # Hyper parameter for phi

# Helper methods for data generator
#  generator for phi
generate.phi <- function(i,is.row,nrow,ncol){
  res <- matrix(0,nrow,ncol)
  if(is.row){
    res[i,] <- 1
  } else {
    res[,i] <- 1
  }
  res/sum(res)
}

# Dirichlet random variables generator
rdirichlet <- function(N,alpha){
  d <- length(alpha)
  res <- matrix(0,N,d)
  for( i in 1:d ){
    res[,i] <- rgamma(N,alpha[i])
  }
  res <- res / rowSums(res)
  class(res) <- "dirichlet"
  attr(res,"alpha") <- alpha
  res
}

# Generating p(w|z) for each topic
phi <- list()
for(i in 1:nrow){
  phi[[i]] <- generate.phi(i,T,nrow,ncol)
}
for(i in 1:ncol){
  phi[[nrow+i]] <- generate.phi(i,F,nrow,ncol)
}
theta <- rdirichlet(ndoc,rep(1,nTopics))

# Generating p(z|d) for each document
documents <- matrix(0,ndoc,W)
for(d in 1:ndoc){
  prob <- matrix(0,nrow,ncol)
  for(i in 1:nTopics){
    prob <- prob + theta[d,i] * phi[[i]]
  }
  the.document <- factor(sample(1:W,words.par.doc,T,prob=as.numeric(prob)),levels=1:W)
  documents[d,] <- as.numeric(table(the.document))
}

# Writing out the generated word distribution
# over each document to CSV file
write.table(documents,file="random_document.csv",sep=",",row.names=F,col.names=F)

# Visualization setting
grays <- gray.colors(10,0.05,0.95,gamma=0.3)
par(mfrow=c(5,nTopics)) # deviding plot area
par(mar=rep(1,4))  # defining plot margin
# Visualize topics (phi)
lapply(phi,function(x) image(x,col=grays,xaxt="n",yaxt="n"))
# Visualize documents
par(mfg=c(3,1)) # identify plot region
apply(documents[1:(nTopics*3),],1,function(x) image(matrix(x,nr=nrow,nc=ncol),col=grays,xaxt="n",yaxt="n"))

■ ギブスサンプラー。メモリも for も使いすぎ。

ndoc <- 2000   # The number of documents to be generated
nrow <- 5      # The number of rows for each document
ncol <- 5      # The number of columns for each document
nTopics <- nrow+ncol  # The number of topics
W <- nrow*ncol # The number of vocabralies
words.par.doc <- 100
alpha <- 1     # Hyper parameter for theta
               # (topic prior probability for each document)
beta <- 0      # Hyper parameter for phi


library(plyr)

X <- read.table("random_document.csv",sep=",",header=F)

# initialize topic
random.topic <- function(doc,nTopics){
  nword <- sum(doc)
  topics <- factor(sample(1:nTopics,nword,T),levels=1:nTopics)
  words <- factor(rep(1:W,doc),levels=1:W)
  docstr <- data.frame(word=words,topic=topics)
  list(wordCount=table(docstr),topicCount=table(topics),docstr=docstr)
}

# initialize topic for each document
docstr <- list()
n.wj <- matrix(0,W,nTopics)
n.dj <- matrix(0,ndoc,nTopics)
for(d in 1:ndoc){
  nn <- random.topic(X[d,],nTopics)
  n.wj <- n.wj + as.matrix(nn$wordCount)
  n.dj[d,] <- as.numeric(nn$topicCount)
  docstr[[d]] <- nn$docstr
}

# log likelihood function
log.lik.fun <- function(nwj){
  Pwj=t(apply(nwj,1,function(x){ifelse(x==0,0,log(x/sum(x)))}))
  sum(nwj*Pwj)
}

# Iteration
nrep <- 100
log.lik <- rep(0,nrep)
nwj.list <- list()
record.index <- c(1,2,5,10,20,30,40,50,75,100)
record.count <- 1
for(r in 1:nrep){

  # Recording transition
  log.lik[r] <- log.lik.fun(n.wj)
  if( r %in% record.index ){
    nwj.list[[ record.count ]] <- n.wj
    record.count <- record.count + 1
  }

  # LDA Iteration
  for(d in 1:ndoc){
    nword <- dim(docstr[[d]])[1]
    for(w in 1:nword){
      word <- docstr[[d]][w,1]
      oldtopic <- docstr[[d]][w,2]
      n.wj[word,oldtopic] <- n.wj[word,oldtopic] - 1
      n.dj[d,oldtopic] <- n.dj[d,oldtopic] - 1
      p1 <- (n.wj[word,]+beta)/(colSums(n.wj)+W*beta)
      # p2 <- (n.jd[d,]+alpha)/(sum(n.jd[d,])+nTopics*alpha)
      p2 <- n.dj[d,]+alpha
      p <- p1*p2
      newtopic <- sample(1:nTopics,1,prob=p)
      docstr[[d]][w,2] <- newtopic
      n.wj[word,newtopic] <- n.wj[word,newtopic] + 1
      n.dj[d,newtopic] <- n.dj[d,newtopic] + 1
    }
  }
}

jpeg(file="LDA_result.jpg",width=1200,height=1000)
tmp <- t(matrix(1:(nTopics*length(nwj.list))+1,nr=nTopics))
layout.matrix <- rbind(cbind(matrix(rep(1,2*length(nwj.list)),nc=2),tmp),
                       c(1,1,rep(max(tmp)+1,nTopics)))
layout(layout.matrix,height=c(rep(1,length(nwj.list)),0.35))
par(mar=c(4.5,4,1,1))
par(oma=rep(0,4))
plot(log.lik,1:length(log.lik),ylim=rev(range(1:length(log.lik))),
     xlab="LogLik",ylab="Iteration",type="l",cex.lab=1.5,cex.axis=1.5)
points(log.lik[record.index],record.index,cex=2)
par(mar=rep(1,4))
grays <- gray.colors(10,0.05,0.95,gamma=0.3)
lapply(nwj.list,function(x){
  apply(x,2,function(y){
    image(matrix(y,nr=nrow),col=grays,xaxt="n",yaxt="n")
    box()
  })
})
dev.off()

[Ruby] 10行で書ける Dijkstra 法

Dijkstra法はグラフ上の最短経路を求めるアルゴリズムです。このアルゴリズムの戻り値はある基点となるノードからその他のすべてのノードまでの距離の配列です。

解説はウェブ上にも豊富にあるので省略することにして、ここではどのくらい短いコード量で Dijkstra法を実装できるのか、という遊びをしてみます。仕様としては

  • 隣接リストではなく隣接行列を受け取る
  • ある二点間に辺が存在しない ⇔ 隣接行列の要素が負の値
  • ある辺の距離は隣接行列の正の要素で表現
  • 出力は指定した頂点からすべての頂点までの距離
  • 外部ライブラリを使用しない

くらいのものを考えましょう。

こんな感じになりました。

“[Ruby] 10行で書ける Dijkstra 法”の続きを読む