Gibbs 采样器无法收敛

Question

一段时间以来，我一直在尝试理解 Gibbs 抽样。最近，我看到一个视频，很有意义。

https://www.youtube.com/watch?v=a_08GKWHFWo

作者使用 Gibbs 采样收敛到一个二元正态分布的均值（theta_1 和 theta_2），使用过程如下：

init：将 theta_2 初始化为一个随机值。

循环：

1. 以 theta_2 为条件的样本 theta_1 为 N~(p(theta_2), [1-p**2])
2. 以 theta_1 为条件的样本 theta_2 为 N~(p(theta_1), [1-p**2])

（重复直到收敛。）

我自己尝试了这个并遇到了一个问题：

import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal

rv = multivariate_normal(mean=[0.5, -0.2], cov=[[1, 0.9], [0.9, 1]])

rv.mean
>>> 
array([ 0.5, -0.2])

rv.cov
>>>
array([[1. , 0.9],
       [0.9, 1. ]])

import numpy as np
samples = []

curr_t2 = np.random.rand()
def gibbs(iterations=5000):
    theta_1 = np.random.normal(curr_t2, (1-0.9**2), None)
    theta_2 = np.random.normal(theta_1, (1-0.9**2), None)
    samples.append((theta_1,theta_2))
    for i in range(iterations-1):
        theta_1 = np.random.normal(theta_2, (1-0.9**2), None)
        theta_2 = np.random.normal(theta_1, (1-0.9**2), None)
        samples.append((theta_1,theta_2))
gibbs()

sum([a for a,b in samples])/len(samples)
>>>
4.745736136676516

sum([b for a,b in samples])/len(samples)
>>>
4.746816908769834

现在，我知道我哪里搞砸了。我发现 theta_1 取决于 theta_2 的实际值，而不是它的概率。同样，我发现 theta_2 取决于 theta_1 的实际值，而不是它的概率。

我遇到的问题是，如何评估任一 theta 取任何给定观察值的概率？

我看到了两个选项：概率密度（基于正态曲线上的位置）和 p 值（从无穷大（和/或负无穷大）到观察值的积分）。这些解决方案听起来都不“正确”。

我应该如何进行？

Answer 1

也许我的视频不够清晰。该算法不会收敛于“平均值”，而是收敛于分布中的样本。尽管如此，来自分布的样本的平均值将收敛到它们各自的平均值。

问题在于您的条件方法。在视频中，我选择了零的边际均值来减少符号。如果您有非零边际均值，conditional expectation for a bivariate normal 涉及边际均值、相关性和标准差（在您的双变量正态中为 1）。更新后的代码是

import numpy as np
from scipy.stats import multivariate_normal

mu1 = 0.5
mu2 = -0.2
rv = multivariate_normal(mean=[mu1, mu2], cov=[[1, 0.9], [0.9, 1]])

samples = []

curr_t2 = np.random.rand()
def gibbs(iterations=5000):
    theta_1 = np.random.normal(mu1 + 0.9 * (curr_t2-mu2), (1-0.9**2), None)
    theta_2 = np.random.normal(mu2 + 0.9 * (theta_1-mu1), (1-0.9**2), None)
    samples.append((theta_1,theta_2))
    for i in range(iterations-1):
        theta_1 = np.random.normal(mu1 + 0.9 * (theta_2-mu2), (1-0.9**2), None)
        theta_2 = np.random.normal(mu2 + 0.9 * (theta_1-mu1), (1-0.9**2), None)
        samples.append((theta_1,theta_2))

gibbs()

sum([a for a,b in samples])/len(samples)
sum([b for a,b in samples])/len(samples)