吉布斯抽样未产生预期结果

Question

我正在尝试学习 Gibbs sampling 的机制。我有 2 个变量，我试图从中进行推断。此示例仅假设高斯分布。我在 R 中的代码如下所示。

library(condMVNorm)

rm(list=ls())

means <- c(0, 25)
cov <- matrix(c(1.09, 1.95, 1.95, 4.52), 2, 2)

k <- 10
initSample <- c(0, 0)
traceSamples <- matrix(, k, 2)

for (i in 1:k) {
    X <- initSample[1]
    c1 <- rcmvnorm(n=1, mean=means, sigma=cov, dep=2, given=1, X=X)

    X <- c1
    c2 <- rcmvnorm(n=1, mean=means, sigma=cov, dep=1, given=2, X=X)

    currentSample <- c(c1, c2)
    traceSamples[i, ] <- currentSample
    initSample <- currentSample
}

colMeans(traceSamples)

我得到的输出如下。

[1] 2220.7619  947.3168

我原以为第一个变量会非常接近 25，第二个变量会非常接近 0。

我不知道我对 Gibbs 抽样的理解是否有误，因为文献总是说你从条件分布中抽样 p(X1=x1|X2=x2)。对我来说，p(X1=x1|X2=x2) 是给定 X2=x2 的 X1=x1 的密度估计，人们会将其映射到 dcmvnorm 而不是 rcmvnorm。

打印 traceSamples 矩阵，我得到以下内容。

           [,1]         [,2]
 [1,]   22.0574   -0.7827272
 [2,]   63.6865   16.3375931
 [3,]  138.7078   49.2994688
 [4,]  272.0850  107.3952335
 [5,]  510.2272  208.3522406
 [6,]  940.7504  395.4438929
 [7,] 1708.2603  725.3048137
 [8,] 3080.5096 1317.7650679
 [9,] 5538.0734 2378.8674730
[10,] 9933.2615 4275.1848015

值似乎在增加（所以这表明我的 R 代码有问题）。此外，我还做了一个没有for循环的非常简单的采样。

means <- c(0, 25)
cov <- matrix(c(1.09, 1.95, 1.95, 4.52), 2, 2)
x1 <- rcmvnorm(n=1, mean = means, sigma = cov, dep=2, given=1, X=c(0))
x2 <- rcmvnorm(n=1, mean = means, sigma = cov, dep=1, given=2, X=c(x1))

x1 <- rcmvnorm(n=1, mean = means, sigma = cov, dep=2, given=1, X=c(x2))
x2 <- rcmvnorm(n=1, mean = means, sigma = cov, dep=1, given=2, X=c(x1))

我的 x1 和 x2 值如下。

23.40496 -0.01044726
22.67643 -0.6836546

对我做错了什么有什么想法吗？

请注意，我能够使用以下代码获得更好的预期结果。

means <- c(0, 25)
cov <- matrix(c(1.09, 1.95, 1.95, 4.52), 2, 2)

k <- 9000
x1 <- 0
x2 <- 0
traceSamples <- matrix(, k, 2)

for (i in 1:k) {
    x1 <- rcmvnorm(n=1, mean=means, sigma=cov, dep=2, given=1, X=x2)
    x2 <- rcmvnorm(n=1, mean=means, sigma=cov, dep=1, given=2, X=x1)

    traceSamples[i, ] <- c(x1, x2)
}

colMeans(traceSamples)

有人能告诉我重用和重新分配有什么问题吗initSample？

Answer 1

我在这里解决了为什么 Gibbs 在模拟中提供错误值的问题，但我认为以这种方式编写代码时会变得越来越复杂，我认为可以删除一些行以构建代码一种更有效的方法，也更快。但是，请注意我在 x <-initSample 和 X = X[1] 以及 X = X[2].

中所做的更改

library(condMVNorm)
rm(list=ls())

means <- c(0, 25)
cov <- matrix(c(1.09, 1.95, 1.95, 4.52), 2, 2)
k <- 9000
initSample <- c(0,0)
traceSamples <- matrix(, k, 2)
for (i in 1:k){
  X <- initSample
  c1 <- rcmvnorm(n=1, mean=means, sigma=cov, dep=2, given=1, X=X[2])
  X <- c1
  c2 <- rcmvnorm(n=1, mean=means, sigma=cov, dep=1, given=2, X=X[1])
  currentSample <- c(c1, c2)
  traceSamples[i, ] <- currentSample
  initSample <- currentSample
}
> head(traceSamples,10)
                  [,1]                [,2]
 [1,] 23.8233821520619 -0.9169596237697860
 [2,] 22.8293033255339 -1.6287517329781345
 [3,] 21.3923155517845 -1.9104909272586084
 [4,] 20.5331401021848 -2.3320921649401360
 [5,] 21.4287399563041 -1.1376683051591154
 [6,] 23.4335659872032 -0.4379604108831421
 [7,] 25.4074041761893 -0.0613743089436460
 [8,] 24.2471298284230  0.0764901351102767
 [9,] 24.7450703427834 -1.2443499508519478
[10,] 24.2193799579308 -0.4995919725966815
> cov.wt(traceSamples)
$cov
                 [,1]             [,2]
[1,] 4.54864368811939 1.96444834328156
[2,] 1.96444834328156 1.09723665614730

$center
[1] 24.9626145462517535 -0.0163323659130855

$n.obs
[1] 9000

吉布斯采样器是一种马尔可夫链Monte Carlo (MCMC) 算法。因此，您应该检查链的收敛性。 coda 包提供了一些非常有用的测试。

library(coda)
MC <- mcmc(traceSamples)
plot(MC)
heidel.diag(MC)
Stationarity start     p-value
     test         iteration        
var1 passed       1         0.231  
var2 passed       1         0.193  

     Halfwidth Mean    Halfwidth
     test                       
var1 passed    24.9626 0.1228   
var2 failed    -0.0163 0.0598

其中接受马尔可夫链来自平稳分布的原假设。

吉布斯抽样未产生预期结果

Gibbs sampling is not producing expected results

r

random

simulation

markov-chains