指数型分布族のメモ

モチベーション

推定や検定の理論を扱うときに、正規分布の場合~、ポアソン分布の場合~、みたいに分布ごとに性質を調べてたんじゃやってられない、ということで、さまざまな密度関数をある種のクラスとしてまとめて考えることで、一般論を展開したい。そのための枠組みの一つが指数型分布族。

実際、指数型分布族に該当する分布は推定や検定などで「極めていいかんじの」性質を持っている。一方で実用上よく使う分布の多くが、指数型分布族というクラスを「継承」している(ほとんどが、と言ってしまってもいいかもしれないが、機械学習系の確率モデルには当てはまらないことが多い)。

要するに、指数型分布族の理論について知っていれば、実用上よく使う分布たちが極めてよい性質を持っていることが統一的に理解できる。ちなみに、GLM(一般化線形モデル)も基本的には指数型分布族前提で定式化されている。

適当に書いて見ましたが、こんなところでしょうか。。。

定義

n次元の確率変数 \(X\) が d-パラメータの指数型分布族に従うとは、あるパラメトリゼーション \(\boldsymbol{\theta}\) と d 個の関数 \(T_i:\mathbb{R^n}\to\mathbb{R}\hspace{5pt}(i=1\cdots d)\) が存在して、\(X\) の密度関数が
\begin{align*}
f(x|\boldsymbol{\theta}):=h(x)\cdot\exp\big(\boldsymbol{\theta}\cdot\boldsymbol{T}(x)-A(\boldsymbol\theta)\big)
\end{align*}
と書けること。\(h, A\) はそれぞれスカラー値をとる関数。\(T_i(\cdot)\) は十分統計量と呼ばれる。

指数型分布族の分布の例

正規分布、ポアソン分布、指数分布、ガンマ分布、ディリクレ分布、負の二項分布などなど。数え出したらきりがない。

指数型分布族でない分布の例

一様分布、コーシー分布など。混合ガウス分布とか隠れマルコフモデルとか、その手のいわゆる「特異モデル」は該当しない。

Natural exponential family

分布族を
\begin{align*}
f(x|\boldsymbol{\theta}):=h(x)\cdot\exp\big(\boldsymbol{\theta}\cdot x-A(\boldsymbol\theta)\big)
\end{align*}
に限定したものが natural exponential family。

平均と共分散の計算

指数分布族に従う d 次元確率変数 \(T(X)=(T_1(X),\cdots,T_d(X))\) の平均と共分散は
\begin{align*}
E(T_i(X))&=\frac{\partial A}{\partial\theta_i}\\
Cov(T_i(X),T_j(X))&=\frac{\partial^2 A}{\partial\theta_i\partial\theta_j}
\end{align*}
となる。なぜなら、正規化条件
\begin{align*}
1=\int h(x)\cdot\exp\big(\boldsymbol{\theta}\cdot\boldsymbol{T}(x)-A(\boldsymbol\theta)\big)dx
\end{align*}
の両辺を \(\theta_i\) で微分すれば
\begin{align*}
0=E_\theta(T_i(X))-\frac{\partial A}{\partial\theta_i}\cdot1
\end{align*}
となり、もう一度、今度は \(\theta_j\) で微分すれば
\begin{align*}
0=&\int T_i(x)T_j(x)\cdot h(x)\cdot\exp\big(\boldsymbol{\theta}\cdot\boldsymbol{T}(x)-A(\boldsymbol\theta)\big)dx\\
&-\frac{\partial A}{\partial \theta_j}\int T_i(x)\cdot h(x)\cdot\exp\big(\boldsymbol{\theta}\cdot\boldsymbol{T}(x)-A(\boldsymbol\theta)\big)dx
-\frac{\partial^2 A}{\partial\theta_i\partial\theta_j}\\
=&E(T_i(X)T_j(X))-E(T_i(X))E(T_j(X))-\frac{\partial^2 A}{\partial\theta_i\partial\theta_j}\\
=&Cov(T_i(X),T_j(X))-\frac{\partial^2 A}{\partial\theta_i\partial\theta_j}
\end{align*}
ということでOK(積分と極限の交換の条件とかは無視してしまいましたが)。

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()

JAGSを使ってギブスサンプリングを試してみた (1)

前から気になってた JAGS を試してみました。Just Another Gibbs Sampler を略して JAGS。Bugs ファミリーの中では現在、最もアクティブに開発がすすめられているようです。

ここでは、とにかく何も知らない状態から R+JAGS で簡単なモデル(ガウスノイズな線形回帰)の推定までできるように最低限のことをメモしました。

準備

JAGS本体のインストール

JAGS は R とは分離された独立なソフトウェアなので R の CRAN 経由で一発インストールというわけにはいきません。ですが、Ubuntu の場合は apt で一発です。他のOSは試していないので、インストールしていない人は別途調べてください。

sudo apt-get install jags

rjagsのインストール

JAGSをRから呼び出すためには rjags というパッケージをインストールする必要があります。apt 使うなら

sudo apt-get install r-cran-rjags

もしくはRのプロンプトで

install.pacakges("rjags")

とすればOK。

テスト

テスト問題として

\begin{align*}
y_i&\sim \text{Normal}(\mu_i,\sigma^2)\\
\mu_i&=\alpha+\beta x_i
\end{align*}
というデータの発生メカニズムを仮定します。これに対応した次のような「モデルファイル」をつくって、適当に名前をつけておきます(lm.jagsとする)。

model {
  for(i in 1:N){
    mu[i] <- alpha + beta*x[i]
    y[i] ~ dnorm(mu[i],1/(tau^2))
  }
  alpha ~ dnorm(0,0.0001)
  beta ~ dnorm(0,0.0001)
  tau ~ dgamma(0.1,0.1)
}

見た感じ、Rのコードにかなり似ているけど、実際には細かい点で異なります。dnormが正規分布なのは同じだけど、第一引数は存在しないし、JAGSでは分散ではなくて精度(precision;分散の逆数)を指定することになっているので注意。Rを知っていればなんとなくわかってしまう構文は素晴らしいですが、上で挙げたような注意の他に、分布の名前自体がちがっていたりすることもあるのでハマる可能性もあります。ちなみに、for ブロックの後ろの3行は未知パラメータの事前分布を設定しています。

データの生成

\(\alpha=2, \beta=3\) としてデータを発生させます。

# determine the number of samples
N <- 50
# data generation
x <- rnorm(N) 
y <- 2 + 3*x + rnorm(N,0.1)

JAGSをRから呼び出す

モデルファイルを保存したら、Rで以下のスクリプトを走らせます。すると下のようなグラフが表示されるはずです。coda[1]でプロットさせているのでMCMCglmmなど他のRのベイズパッケージを使った時と同じ図になります。

library(rjags)
library(coda)

# compiling and initializing
lm.mcmc <- jags.model("lm.jags",
                       data=list('x'=x, 'y'=y, 'N'=N),
                       n.chains=4,  # the number of parallel chains
                       n.adapt=100)

# mcmc sampling
update(lm.mcmc,1000) # burn-in
posterior <- coda.samples(lm.mcmc, c('alpha','beta'), 2000)
# plotting the result
plot(posterior) # plot.coda is called


こんな感じで推定出来ました。すこしスクリプトについて書いておくと

  • jags.model 関数の第一引数では先程のモデルファイルを指定する。
  • 第二引数 data では JAGS に渡したいRオブジェクトをリスト形式で与える。クオーテーションする必要があるので注意。
  • n.chains で並列パスの数を指定
  • update で burn-in している
  • coda.samples で指定した数のサンプルを coda 形式で取得する。

まとめといろいろと蛇足

階層ベイズでも何でもないつまらない例を使って、ほんとに最低限のJAGSの使い方について書きました。たぶん続編があるはずです。次は混合ポアソン回帰でも試してみようかな。

わりとどうでもいいことですが、’J’AGS という名前(←JRubyみたいな意味で。)とクロスプラットフォームを謳っていることから Java で書かれていると思っていましたが、c++で実装されているようです。この勘違いのせいでずいぶん長い間敬遠してきた、という経緯があります。

ちなみに、これまで自分でMCMC実装するときは

  • cで書く
  • Rでものすごく遅いやつを書く(ちかいうちにこのブログにも登場予定)
  • MCMCglmmみたいなRのパッケージを使う

という選択肢の中から選んでいたのですが、これにJAGSを追加して比較するとこんな感じでしょうか?(いろいろと異論はあるかもしれませんが)

自由度 実行スピード 実装の手軽さ
c言語 ×
R ×
R-packages × ○ (?)
JAGS ○ (?)

おしまい。

  1. [1] RでMCMCの結果を可視化するときのほぼ事実上の標準形式