吉布斯采样器无法收敛

Question

一段时间以来，我一直在努力理解 Gibbs 抽样。最近，我看到一个视频，很有道理。

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

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

init: 将 theta_2 初始化为运行dom 值。

循环：

样本 theta_1 以 theta_2 作为 N~(p(theta_2), [1-p**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 的实际值，而不是它的概率。

我被卡住的地方是，我如何评估任一 theta 取任何给定观测值的概率？

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

我该如何进行？

Answer 1

可能是我的视频不够清晰。该算法不收敛 "on the mean values" 而是收敛于分布中的样本。尽管如此，来自分布的样本的平均值将收敛到它们各自的平均值。

问题出在你的条件手段上。在视频中，我选择了零的边际均值来减少符号。如果您的边际均值非零，则 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)

吉布斯采样器无法收敛

Gibbs sampler fails to converge

python

random

normal-distribution

sampling

mcmc