ディリクレ分布に関するメモ (6)

ディリクレ分布に関するメモ (3) ではガンマ乱数を正規化することでディリクレ分布に従う乱数を発生させました。

今回は別の方法を試してみます。

ポリアの壺

ポリアの壺（Polya’s urn）というのは以下のような手続きのことです。

n種類の色のボールが入った壺を考える
壺から一つボールを取り出す
取り出したボールを壺に戻し、同じ色のボールを壺に追加する
ステップ２に戻る

この手続きによって、壺の中のボールは増加し続けます。その色のシェアはランダムに変動しますが、やがて収束します。この収束した時のシェアを $p^*$ とします。この極限のシェア $p^*$ はボールを取り出した経路に依存し、これ自体が確率変数です。これは、例えば壺に赤と青のボールが一つずつ入っていたとして、最初の数回でたまたま赤が連続して出ると、その後も赤が出やすくなり、極限では赤のシェアが90%ということも起こりえるわけです。

さて、n色のポリアの壺を考えます。最初に $\alpha=(\alpha_1,\cdots,\alpha_n)$ 個のボールが壺に入っていたとすると、このポリアの壺の極限分布が $\text{Dirichlet}(\alpha)$ に従います（ちなみに、これは $\alpha$ が自然数でなくても成立します）。

R でポリアの壺を実装してみる

以下のコードの rdirichlet.polya がポリアの壺の実装になります。極限までは計算できないので、napprox で打ち切ります。ディリクレ分布に関するメモ (3) で計算したときと同じパラメータ $\alpha$ で乱数を発生させてみます。コードは以下のとおりです（非常に重いです）。

rdirichlet.polya <- function(n, alpha, napprox=1000){
  res <- matrix(0,n,length(alpha))
  for(k in 1:n){
    p <- alpha
    idx <- 1:length(alpha)
    for(i in 1:napprox){
      id <- sample(idx,1,F,p)
      p[id] <- p[id] + 1
    }
    res[k,] <- p
  }
  res <- res/rowSums(res)
  class(res) <- "dirichlet"
  attr(res,"alpha") <- alpha
  res
}

plot.dirichlet <- function(x){
  # x must be Nx3 matrix!!
  X <- x[,2]-x[,1]
  Y <- x[,3]
  par(xaxt="n")
  par(yaxt="n")
  par(bty="n")
  par(mar=rep(0,4))
  plot(X,Y,xlim=c(-1.2,1.2),ylim=c(-0.2,1.2),xlab="",ylab="",pch=20,cex=0.7)
  lines(c(-1,1,0,-1),c(0,0,1,0),col="gray")
  text(-1,0,"x1"); text(1,0,"x2"); text(0,1,"x3")
  alpha <- attr(x,"alpha")
  stralpha <- paste("a1=",alpha[1],"\na2=",alpha[2],"\na3=",alpha[3],sep="")
  text(-0.8,0.8,stralpha)
  print(alpha)
  symbols(-1,0,circles=0.08*sqrt(alpha[1]),inches=F,add=T)
  symbols(1,0,circles=0.08*sqrt(alpha[2]),inches=F,add=T)
  symbols(0,1,circles=0.08*sqrt(alpha[3]),inches=F,add=T)
}

alphas <- list(c(1,1,1),c(3,2,1),c(6,2,1))
scale <- c(0.1,1,8)

par(mfrow=c(3,3))
for(alpha in alphas){
  for(s in scale){
    plot(rdirichlet.polya(1000,s*alpha))
  }
}