Gibbs Sampler Algorithmによって多変量正規分布からのサンプル抽出を行なう

もちろんRで。

2変数正規分布

2変数の正規分布からのサンプリングはここのをそのまままねすると

B <- 2000
x1 <- rep(NA,B+1)
x2 <- rep(NA,B+1) 
x1[1] <- -2
x2[1] <- 1 
for (i in 1:B){ 
  x1[i+1] <- rnorm(1,1+0.7*(x2[i]-2),sqrt(1-0.7^2)) 
  x2[i+1] <- rnorm(1,2+0.7*(x1[i+1]-1),sqrt(1-0.7^2)) 
}

相関係数もこんな感じになって、

> cor(x1,x2)
[1] 0.712666

散布図もこんな感じ。うまくいってますね。

面倒なのでbuild-inとか考えてない(ry。

3変数正規分布

Gibbs Sampler Algorithmを使っているので、full conditioningした分布が分かってないといけないです。が、すでにPRMLで勉強したので、どうなるかは知っています。P85の付近ですね。どうなっているかと言えば


となっています。これさえあれば多次元正規分布からのサンプリングもやれちゃうぜ!ということでやりました。サンプリングしてきたい3変量正規分布はこんなのです。

mu1 <- 1
mu2 <- 3
mu3 <- 5

cor12 <- 0.3
cor13 <- 0.5
cor23 <- 0.7

cor21 <- cor12
cor31 <- cor13
cor32 <- cor23

B <-2000

x1 <- rep(NA,B+1)
x2 <- rep(NA,B+1)
x3 <- rep(NA,B+1) 

x1[1] <- 1
x2[1] <- 3
x3[1] <- 5

for (i in 1:B){
  x1[i+1] <- rnorm(1,mu1+c(cor12,cor13) %*% solve(matrix(c(1,cor23,cor32,1),2,2,byrow=TRUE)) %*% c(x2[i]- mu2,x3[i]-mu3),sqrt(1 - c(cor12,cor13) %*% solve(matrix(c(1,cor23,cor32,1),2,2,byrow=TRUE)) %*% c(cor21,cor31)))
  x2[i+1] <- rnorm(1,mu2+c(cor23,cor21) %*% solve(matrix(c(1,cor31,cor13,1),2,2,byrow=TRUE)) %*% c(x3[i]- mu3,x1[i+1]-mu1),sqrt(1 - c(cor23,cor21) %*% solve(matrix(c(1,cor31,cor13,1),2,2,byrow=TRUE)) %*% c(cor23,cor21)))
  x3[i+1] <- rnorm(1,mu3+c(cor31,cor32) %*% solve(matrix(c(1,cor12,cor21,1),2,2,byrow=TRUE)) %*% c(x1[i+1]- mu1,x2[i+1]-mu2),sqrt(1 - c(cor31,cor32) %*% solve(matrix(c(1,cor12,cor21,1),2,2,byrow=TRUE)) %*% c(cor13,cor23)))
}

平均ベクトル、共分散行列も大体よい感じですね。

> apply(data.frame(x1,x2,x3),2,mean)
       x1        x2        x3 
0.9911853 2.9345036 4.9296849 
> cor(data.frame(x1,x2,x3))
          x1        x2        x3
x1 1.0000000 0.3081783 0.5300602
x2 0.3081783 1.0000000 0.6915105
x3 0.5300602 0.6915105 1.0000000

一応iterationの経過も載せておく。

plot.iteration <- function(x){
  plot(x,type="l",main=paste("trace of ",deparse(substitute(x)),sep=""),xlab="iteration")
}

par(mfrow=c(3,1))
plot.iteration(x1)
plot.iteration(x2)
plot.iteration(x3)

問題なさそう。