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

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

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

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

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

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

CRPのテーブルの数の分布

Chinese Restaurant Process (→ 以前の記事) でデータ数(レストランに来る人数/壺からボールを取り出す回数)や \(\alpha\) が変化した時に利用されるテーブルの数の分布がどうなるか実験してみた(下の図をクリックで拡大)。


“CRPのテーブルの数の分布”の続きを読む

正規分布/Normal-inverse-Wishart が事後分布に収束していく様子

正規分布のパラメータ \(\mu, \Sigma\) の共役事前分布は Normal-inverse-wishart (NIW) 分布。データ数が増加した時に真のパラメータに収束していく様子を図示してみた (クリックで拡大する、かも)。


絵のせつめい。

  • データ数 n を変化させた各グラフの中で、事後分布から20組のパラメータをサンプリング。
  • 一組のパラメータを90%信頼楕円として表している
  • 点線は真の分布。黒丸は事前分布のモード(最頻値)。
  • データが少ないうちは事後分布がばらつく、つまりパラメータの不確実性が大きいが、データ数が256を超えたあたりからほぼ真の分布に収束する

下のコードのコメントにも書いたが、ハイパーパラメータを学習しない場合のレシピとして、NIW分布の平均パラメータ \(\mu_0\)、分散パラメータ \(V_0\) に自信がないときはそれぞれ \(k_0, \nu_0\) を小さく設定すればいい、はず。事前分布の影響を小さくできる。

参考その1:前回のエントリ → 逆Wishart分布を図示してみる
参考その2:事後分布のパラメータの求め方 → Wikipedia/Conjugate prior

お絵かきスクリプト in R。ellipse, MCMCpack, mvtnorm パッケージは入っていなければインストールする必要あり。

library(ellipse)   # conffidence ellipse
library(MCMCpack)  # wishart
library(mvtnorm)

#----------------------------------------------------------------------------
# Baysian estimation of 2 dimensional normal distribution
#----------------------------------------------------------------------------

d <- 2     # dimension

#----------------------------------------------------------------------------
# Hyperparameters for NIW (Normal inverse wishart)
k0 <- 0.1        # see below
mu0 <- rep(0,d)  # hyper mean
v0 <- 3.5        # see below
V0 <- diag(rep(10,d)) # hyper variance
Phi0 <- V0*(v0-d-1)
prior.hyper.par <- list(k=k0,mu=mu0,v=v0,Phi=Phi0)

#------------------------------------------------------------------------------------
# [Recipe for determining hyperparameters v0 and k0]
#  1. Less confidence for the hyper mean(mu0), take smaller k0; typically, 0<k0<<1
#  2. Less confidence for the hyper variance(V0), take smaller v0; typically, v0~1+p
#------------------------------------------------------------------------------------


#-----------------------------------------
# Bayesian estimation ( bayesian update)
bayes.update <- function(X,hyper.par){

  k0 <- hyper.par$k
  mu0 <- hyper.par$mu
  v0 <- hyper.par$v
  Phi0 <- hyper.par$Phi

  if(!is.matrix(X)){
    # This is n=1 case that we have to deal with as special because of R issue
    n <- 1
    EX <- X
    X <- t(X)
    C <- matrix(rep(0,d*d),nc=d)
  } else {
    n <- dim(X)[1]
    EX <- colMeans(X)
    C <- (n-1)*cov(X)
  }
  mu <- (k0*mu0 + n*EX)/(k0+n)
  k <- k0+n
  v <- v0+n
  Phi <- Phi0 + C + k0*n/(k0+n)*((EX-mu0) %*% t(EX-mu0))

  list(k=k,mu=mu,v=v,Phi=Phi)  # This new parameters list is a bayesian update result!!

}

#------------------------------------------------------------------
# Sampling from posterior
sample.posterior <- function(n,bayes.fit){
  k <- bayes.fit$k
  v <- bayes.fit$v
  Phi <- bayes.fit$Phi
  mu <- bayes.fit$mu
  result <- list()
  for(i in 1:n){
    # Sampling mu and V from NIW
    V <- riwish(v,Phi)  # 1. sampling covariance matrix by inverse wishart
    result[[i]] <- list()
    result[[i]]$V <- V
    result[[i]]$mu <- rmvnorm(1,mu,V/k) # 2. sampling mu by normal dist.
  }
  result
}

#--------------------------------------------------------------
# Test case
#  --- true distribution
V.true <- matrix(c(3,-1.5,-1.5,1.8),nc=d)
mu.true <- rep(10,d)
#  --- generating random variables
X <- rmvnorm(20000,mu.true,V.true)


par(mfrow=c(3,4))  # dividing graphic device
j <- 1
ids <- 2^seq(0,11)
hc <- rainbow(12)
for(i in ids){
  Y <- X[1:i,]
  bayes.fit <- bayes.update(Y,prior.hyper.par)  # estimating posterior
  post.samp <- sample.posterior(20,bayes.fit)  # sampling from posterior
  plot(Y,pch=".",col="gray",xlim=c(-10,20),ylim=c(-10,20),xlab="",ylab="") # just plotting raw data

  # plotting posterior samples
  for( theta in post.samp ){
    mu <- theta$mu
    V <- theta$V
    elp <- t(apply(ellipse(V,level=0.9),1,function(x) x+mu))
    lines(elp,col=hc[j],lw=2)
  }

  # plotting true distribution
  true.dist <- t(apply(ellipse(V.true,level=0.9),1,function(x) x+mu.true))
  lines(true.dist,type="l",lw=2,lt=4)

  # plotting mode of the prior
  prior <- ellipse(Phi0/(v0+3), level=0.9) + mu0
  lines(prior,lw=2)

  text(20,-7,paste("n =",i),pos=2,cex=2)
  j <- j+1
}

dev2bitmap(file="normal_posterior.jpg",taa=4,gaa=4,width=15,height=12)

逆Wishart分布を図示してみる

正規分布のベイズ推定で共分散行列の(共役)事前分布として使われる逆Wishart分布(→ wikipedia:en)というものがありますが、あまりイメージしにくかったので図にしてパラメータの変化がどのような違いをもたらすか調べてみました。

といっても逆Wishart分布は n x n 行列を値にとる確率分布なので、2×2 行列を考えても3次元なので(※対称行列)そのままプロットするとわかりにくい。そこで ellipse パッケージを使って可視化してみます。
“逆Wishart分布を図示してみる”の続きを読む

事後分布に収束していく様子をグラフ化してみる(ポアソン/ガンマ分布)

ポアソン分布のパラメータの事前分布としてガンマ分布を仮定すると事後分布もガンマ分布になる(過去の記事 → 共役事前分布について)。

平均3のポアソン分布から100個サンプルをとってきて、オンライン学習的に次々にデータを食わせて事後分布がどのように変化していくかを図示してみた。事前分布としては \(\alpha=0.1, \beta=0.1\) のガンマ分布を使用。

データを増やしていくとなんとなくいいかんじに分布が収束していくことがわかる。こういう図をうまく見せるのって難しい。以下、Rスクリプト。

library(RColorBrewer)
x <- rpois(100,3)
lam <- seq(0,8,0.01)
hc=brewer.pal(11,"RdYlGn")
plot(lam,dgamma(lam,0.1,0.1),type="n",
     ylim=c(0,3),xlab="Lambda",ylab="Density",col="gray")
lines(lam,dgamma(lam,0.1,0.1),col="gray",lw=2)
i <- 1
for( n in c(seq(1,10,2),seq(20,100,20)) ){
  lines(lam,dgamma(lam,0.1+sum(x[1:n]),0.1+n),
        col=hc[i],lw=2)
  i=i+1
}
title("Posterior of Poisson distribution parameter with Gamma(0.1,0.1) prior")
lines(c(3,3),c(0,5),lw=2,lt=2,col="darkblue")
legend("topright",
       legend=c("Prior (Gamma(0.1,0.1))",
                paste("n =",c(seq(1,10,2),seq(20,100,20)))),
       lw=2,col=c("gray",hc))
dev2bitmap(file="PoissonPosterior.png",
           gaa=4,taa=4,width=9,height=6)

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の結果を可視化するときのほぼ事実上の標準形式

Chinese Restaurant Process を発生させてみた

タイトルの通りです。単なる遊びです。Blackwell-McQueen Polya’s urn scheme とか表記されたりもするみたいです。

Chinese Restaurant Process (CRP)

よく(名前にもあるように)中国料理屋の例が出ますが、個人的には「初期の壺に連続的なスペクトルの色(虹色)を持つボールが入っているバージョンのポリアの壺」と考えるほうがわかりやすかったです。つまり、以前に記事にしたポリアの壺の色数が(可算とは限らない)無限集合に拡張されたもの、ということです。兎に角、CRP では以下のスキームで無限の長さの確率変数列を発生させます。

  1. \(G_0\) を適当な確率分布、\(\alpha\) を正の実数とします
  2. \(k\leftarrow1\)とする
  3. \(G_0\) から乱数(色)をひとつサンプリング
  4. その乱数を記録し、同じ色のボールを空の壺にいれる
  5. \(k\leftarrow k+1\)とする
  6. \(\alpha/(\alpha+k-1)\) の確率で \(G_0\) からサンプリングし、それ以外のときは壺からボールを一つ取り出してそれをサンプルとする
  7. サンプルと同じ色のボールを壺に入れる
  8. ステップ5に戻る

実際には無限に数列をつくることはできないので、10000個サンプルした段階でやめることにします。長さ10000の系列を15回生成して、その分布をプロットしてみます。\(G_0\) は標準正規分布、\(\alpha=3\) としました。

# Generating n random variables via CRP(alpha, G0)
rcrp<-function(n,alpha,G0){
  res<-numeric(n); res[1]<-G0(1);
  for(i in 2:n){
    p<-alpha/(alpha+i-1)
    res[i]<-ifelse(runif(1)<p,G0(1),sample(res[1:(i-1)],1))
  }
  class(res)<-"crp"
  res
}

# Simple simulation
par(mfrow=c(5,3)); par(mar=rep(1.5,4)); par(oma=c(1,1,1,1))

N <- 10000

for(i in 1:15){
  plot(table(rcrp(N,3,rnorm))/N,xlim=c(-3,3),ylim=c(0,0.4),ylab="",xaxt="n")
  axis(side=1,at=-2:2)
  curve(dnorm,add=T)
}

ここからわかること:

  • \(G_0\) は連続値をとる密度関数だが、得られる分布は離散的(実はこれは10000個で生成をやめずに、無限個生成しても離散的)
  • 10000個のサンプルのほとんどは数~十数のポイントに落ちる

いろいろ調べると(証明すると)わかること:

  • 実は得られた15個の離散分布はディリクレ過程(←ここでは説明していない)に従う「ランダム確率分布」である
  • (これは、もちろん有限で止めているという意味では正確な表現ではない)
  • 本当のディリクレ過程からのサンプルされる分布は(離散だけども)サポートは無限集合になる
  • CRPによって得られる確率変数列は過去に依存しているので、iidには見えない。しかし、実際にはこれらは得られた離散分布からのiid とみなすことができる

共役事前分布について

定義など少し混乱気味になって、調べてみたので書いてみました。

定義

パラメトリックな確率分布の族 \(\{p(\theta); \theta\in\Theta\}\) を考える。確率変数 \(X\) の従う分布が \(p(x|\theta)\) のとき(もしくはそのように仮定するとき)、ベイズの定理より事後分布は

\begin{align*}
p(\theta|x)\propto p(x|\theta)p(\theta)
\end{align*}

となるが、このとき事前分布 \(p(\theta)\) と事後分布 \(p(\theta|x)\) が同じ分布族に属するときに \(p(\theta)\) は \(p(x|\theta)\) の共役事前分布であるという。同じ分布族に属するので、事前パラメータ \(\theta_0\) とデータ \(x\) が決定されると事後パラメータ \(\tilde\theta\) は \(\theta_0, x\) の関数として\(\tilde\theta(\theta_0, x)\) と書ける。
“共役事前分布について”の続きを読む