すこし予定を変更して、今回はディリクレ分布を可視化してみます。
次元が高いとプロットできないので、３次元のディリクレ分布を可視化することにします。
三次元なのでこのような三角形の上に確率変数が分布します。ディリクレ分布のパラメータをいろいろ変えながら、この三角形上にランダムにサンプルされたディリクレ分布をプロットしていきます。そのためにはディリクレ分布に従う乱数を発生させるアルゴリズムが必要になります。
今回は理由は省略しますが、パラメータ \(\alpha\) をもつディリクレ分布からのサンプリング方法は以下のようになります。【追記】理由を書きました。
- \(i=1,\cdots,d\) に対してシェープパラメータが \(\alpha_i\) 、スケールパラメータが１のガンマ分布を発生させる。これを \(X_i\) とする。
- \(X_0=X_1+\cdots+X_d\) として、\(Y_i=X_i/X_0\) としてこの \(Y\) を記録する
パラメータをいろいろと変化させたグラフを以下に示します。三角形の頂点にある円の大きさは \(\alpha_i\) の大きさに対応しています。サンプリングおよび図を作るための R のプログラムはエントリーの最後に書いておきます。
ここからおおよそ次のようなことがわかります。
- \(\alpha\) のノルムが小さいと三角形の周辺に張り付くように分布し、ノルムが大きいと一点に集中していく。
- \(\alpha_i\) の大きな頂点に分布が集中していく。
これを前回計算した、平均値と分散の数式と見比べると直感的に理解できると思います。
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)
}
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
}
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(1000,s*alpha))
}
}
