马尔科夫链蒙特卡洛方法 metropolis算法运用

MCMC 是Markov Chain Monte Carlo 的简称，但在传统模拟中有一个很重要的假设是样本是独立的（independent samples），这一点在贝叶斯统计尤其是高纬度的模型中很难做到。所以MCMC的目的就是运用蒙特卡洛模拟出一个马可链（Markov chain)。

如今，概率建模风靡一时，但是当我第一次了解它时，总有一件事情困扰我。 许多贝叶斯建模方法都需要计算积分，而我看到的任何工作示例似乎都使用高斯或伯努利分布，原因很简单如果您尝试使用比这更复杂的方法，它将成为分析的噩梦。 将贝叶斯模型限制在”表现良好”的分布的小子集中，可能会极大地阻碍你对问题建模的能力，所以我们必须找到克服这一限制的方法。

蒙特卡洛近似

如果我不想分析计算某个讨厌的积分怎么办#可以使用蒙特卡洛近似。

我们知道，我们可以通过使用目标分布的样本值计算期望通过使用目标分布的样本值计算样本均值。为什么重要#那么，期望是什么呢#

连续随机变量的期望。同样的过程也适用于离散的情况，只要改变求和的积分。

这种估计积分的方法由中心极限定理提供了一些很好的保证。首先，这是期望的无偏估计，其次，我们可以计算估计的方差。

使用蒙特卡罗样本计算积分是非常好的，但是我们如何从目标分布中抽取样本呢#绘制高斯或均匀样本很容易，但np.random会让你失望。

画样本最简单的方法是使用逆CDF方法但这依赖于获得逆CDF函数它通常没有一个很好的解析形式只对一维随机变量有意义。

Metropolis算法是许多马尔可夫链蒙特卡洛（MCMC）采样方法的组成部分之一。 当您可以访问的只是目标分布的pdf时，它使我们能够绘制样本。

MCMC方法需要注意的是我们不再采用独立样本所以我们不能保证估计的方差如何随着样本数量的增加而减少。如果样本是独立的，中心极限定理告诉我们估计值的方差将与样本数量(N)成反比地减少。对于MCMC，我们可以通过将样本数量从N调整为Neff来对其贴现。Neff(几乎)总是小于N，与链中样本的相关性有关。

Metropolis采样

Metropolis算法的步骤如下：

1.从目标分布域或先前分布的域中均匀采样起点。 2.在那时pdf。 3.根据一些状态转换函数，从当前位置走一步，为新样本提出建议。 4.计算新的pdf值。 5.计算新pdf /旧pdf的值。 6.如果比率大于1，请接受该步骤。 7.如果比率小于1： 1.采样一个统一的随机变量。 2.如果比率大于随机样本，请接受该步骤。 8.否则拒绝该建议，将当前职位作为示例添加并采取新的步骤。

请注意，在5–8中描述的过程等效于以概率为min（1，p（new）/ p（old））的伯努利概率接受样本，记住这个后面会用到……

Metropolis为何起作用？

对于任何MCMC方法，我们都希望确保一种称为详细平衡或可逆性的属性。

详细平衡

如果π满足，则π是马尔可夫链（1）的平稳分布。 我们可以通过对等式两边求和来证明这一点。 如果我们可以保证详细的平衡，那么我们也知道我们正在从马尔可夫链的固定分布中取样，我们将其作为目标分布。

详细平衡的思想是，由于每次转移的概率”质量”从状态i到状态j的转移与从状态j到状态i的转移是相同的，因此在链的每次转移之后，我们保持在平稳分布 。

现在，让我们展示一下Metropolis算法如何满足这一条件……

为了找到满足详细平衡的p（i，j），我们首先提出任意的”跳跃”概率q（i，j），然后仅通过接受概率为α（的”跳跃”来获得p（i，j）。 i，j）。 当”跳转”被拒绝时，状态保持j = i。 这种”接受”的想法并不是Metropolis算法独有的，它存在于MCMC的大多数变体中。

跳跃概率推导

这取决于α是有效概率分布。 因此，α的有效形式为：

Metropolis-Hastings 跳跃概率

如果跳跃概率是对称的，则可以简化为：

否则，它可以保留为完整形式，称为Metropolis-Hasting MCMC。

现在我们可以保证详细的平衡，我们可以让马尔可夫链式接管。 如果马尔可夫链是遍历的（所有状态都是不可约的），那么在某个时候，该链将到达平稳分布，并且我们能够从目标分布中获取样本。

你可能还注意到，由于alpha是π（j）/π（i）的函数。 这样的结果意味着目标分布无需标准化。 当将Metropolis用于难以计算证据项的贝叶斯后验估计时，这特别有用。

实现的注意事项

Metropolis算法的通用版本称为”随机行走Metropolis”，其中建议的状态为当前状态，再加上均值为零且协方差矩阵为σ²I的多元高斯。σ应选择为足够使得足够多的样本被拒绝大。 这是为了确保对目标分布进行良好的探索。

要注意的第二件事是老化的概念。 在马尔可夫链到达平稳分布之前采集的样本应删除，因为它们在链收敛之前不能代表目标分布。 确定要删除的样本数量有些困难，但通常，要删除的样本为10–20％（或有效的样本为10–100）。

Numpy代码实现

def metropolis(pi, dims, n_samples, burn_in=0.1, var=1):
    theta_ = np.random.randn(dims)*var    samples = []
    while len(samples) < n_samples:
        theta = theta_ + np.random.randn(dims)*varratio = pi(theta)/pi(theta_)
        if np.random.rand(1) < ratio:
            sample = theta
            theta_ = theta
            samples.append(sample)
        else:
            sample = theta_
            samples.append(sample)    samples = np.array(samples)
    return samples[int(samples*burn_in):,:]

我们可以看到这在两个高斯函数的和上的表现(注意这是一个非正态分布)。

from scipy.stats import multivariate_normal
def make_pdf(mean1, mean2, cov1, cov2):    pdf1 = multivariate_normal(mean1, cov1)
    pdf2 = multivariate_normal(mean2, cov2)    def pdf(x):
        return pdf1.pdf(x) * pdf2.pdf(x)    return pdfpdf = make_pdf(
    [3, 3],
    [-1, -1],
    np.array([[1,0.1],[0.1,1]], dtype=float),
    np.array([[1,0.1],[0.1,1]], dtype=float),
)
samples = metropolis(pdf, 2, 1_000, 0.1)

上面的gif显示了算法是如何遍历分布的，偶尔会在分布的两种不同模式之间跳转。注意，这也突出了metropolis算法的一个弱点，它处理相对较差的多模型分布。这是由于要探索一种新模式，步骤必须足够大，以便从一种模式希望到另一种模式。这要么需要一个大的步长，要么需要一个模态紧密分布。修改如哈密顿量MCMC可以帮助解决这一问题，但一般来说，这是大多数MCMC方法的一个问题。