Gibbs采样器无法收敛

Question

bounty到期在7天内。这个问题的答案有资格获得+ 50声望赏金。 jbuddy_13正在寻找来自信誉良好的来源的答案：

有资格获得赏金奖励的回答者将根据需要调整代码，以创建功能齐全的Gibbs采样器，该采样器可从双变量分布中进行采样。根据需要参考youtube链接以获取更多信息。

我一直在尝试了解Gibbs采样。最近，我看了一段很有意义的视频。

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

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

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

循环：

1. 以then_2为条件，以n〜（p（theta_2），[1-p ** 2]）采样theta_1
2. 以then_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值（从无穷大（和/或负无穷大）到观测值的积分）。这些解决方案都不听起来“正确”。

我应该如何进行？

我一直在尝试了解Gibbs采样。最近，我看了一段很有意义的视频。 https://www.youtube.com/watch?v=a_08GKWHFWo作者使用Gibbs采样来......>

Answer 1

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

问题在于您有条件的手段。在视频中，我选择了零的边际均值以减少表示法。如果您的边际均值非零，则conditional expectation for a bivariate normal包含边际均值，相关性和标准偏差（在双变量法线中为1）。更新的代码是