ディリクレ分布に関するメモ (3) ではガンマ乱数を正規化することでディリクレ分布に従う乱数を発生させました。
今回は別の方法を試してみます。
ポリアの壺
ポリアの壺(Polya’s urn)というのは以下のような手続きのことです。
- n種類の色のボールが入った壺を考える
- 壺から一つボールを取り出す
- 取り出したボールを壺に戻し、同じ色のボールを壺に追加する
- ステップ2に戻る
この手続きによって、壺の中のボールは増加し続けます。その色のシェアはランダムに変動しますが、やがて収束します。この収束した時のシェアを \(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)) } }
結果
以前のガンマ乱数を正規化して発生させた時の結果と比べてみます。
ほとんど同じ分布が得られました。
まとめ
ということで証明はしていませんが、ポリアの壺の極限がディリクレ分布になることがわかったのですが、この乱数生成の方法は実際上役に立つことはほとんどないでしょう。
ただ、このイメージは CRP (Chinese Restaurant Process) を理解するときに重要かもしれません。おしまい。