吉布斯采样器无法收敛

Question

我一直在尝试理解吉布斯采样一段时间。最近，我看到了一个很有意义的视频。

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

作者使用 Gibbs 抽样收敛于二元正态分布的平均值(theta_1 和 theta_2(，使用的过程如下：

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

圈：

1. 样品theta_1以theta_2为条件为 N~(p(theta_2(， [1-p**2](
2. 样品theta_2以theta_1为条件为 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的实际价值为条件，而不是它的概率。

我陷入困境的是，如何评估θ具有任何给定观测值的概率？

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

我应该怎么做？

Answer 1

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

问题在于您的条件均值。在视频中，我选择边际均值为零以减少符号。如果您有非零边际均值，则二元正态的条件期望值涉及边际均值、相关性和标准差(在二元正态中为 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)