蒙特卡洛准则，蒙特卡洛求pi近似值

MCMC是马尔可夫链蒙特卡罗的缩写，但传统的模拟做了一个非常重要的假设，即样本是独立的，这是贝叶斯统计很难实现的，尤其是高纬度模型。因此，MCMC的目的是使用蒙特卡罗模拟马尔可夫链。

随机建模现在很流行，但当我第一次了解它时，有一件事总是困扰着我。许多贝叶斯建模技术都需要计算积分，但到目前为止我见过的实际示例似乎使用高斯或伯努利分布。原因很简单，如果您尝试使用比这更复杂的东西，分析起来就会变成一场噩梦。将贝叶斯模型限制为“运行良好”的分布的一小部分子集可能会严重妨碍它们对问题进行建模的能力，因此我们需要找到克服这一限制的方法。

蒙特卡罗近似如果您不想通过解析方式计算繁琐的积分怎么办？您可以使用蒙特卡罗近似。

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

对连续随机变量的期望。相同的过程适用于离散情况，只是改变总和的积分。

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

使用蒙特卡罗样本计算积分很棒，但是如何从目标分布中采样？采样高斯或均匀样本很容易，但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)，然后仅通过接受概率为(i) 的“跳跃”得到p(i ,j)。j).i,j)。如果“跳转”被拒绝，则状态保持j=i。这种“接受”的想法并不是Metropolis 算法所独有的，并且存在于MCMC 的大多数变体中。

跳跃概率的推导

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

大都会-黑斯廷斯跳跃概率

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

否则，您可以保留其完整格式，称为Metropolis-Hasting MCMC。

现在我们可以保证细节平衡，我们可以让马尔可夫链来接管。如果马尔可夫链是遍历的（所有状态都是不可约的），则在某个点该链将达到平稳分布，我们可以从目标分布中采样。

您可能还注意到alpha 是(j)/(i) 的函数。这样的结果意味着不需要对目标分布进行归一化。当使用Metropolis 进行贝叶斯后验估计（其中证据项难以计算）时，这特别有用。

实现说明Metropolis 算法的常见版本称为“Random Walk 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 table=[] while len(samples) n_samples: theta=theta_ + np.random .randn(dims)*varratio=pi(theta)/pi(theta_) if np.random.rand(1)rate: 样本=theta theta_=thetasample.append(sample) else: 样本=theta_sample.append(sample) 样本=np.array(samples) return table[int(samples*burn_in):] 您可以看到它如何与两个高斯函数的总和一起工作（请注意，这是一个非正态分布。请）。

从scipy.stats 导入multivariate_normaldef make_pdf(mean1,means2,cov1,cov2): pdf1=multivariate_normal(mean1,cov1) pdf2=multivariate_normal(mean2,cov2) def pdf(x): return pdf1.pdf(x) * pdf2.pdf (x) 返回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),)样本=Metropolis(pdf, 2, 1_000, 0.1)

上面的gif 显示了算法如何在分布中移动，有时在两种不同的分布模式之间跳转。请注意，这也凸显了Metropolis 算法的弱点，该算法处理多模型分布的能力相对较差。这是因为探索新模式需要足够大的步骤才能从一种模式过渡到另一种模式。这需要大的步长或紧密的模式分布。哈密顿MCMC 等修改可以帮助解决这个问题，但总的来说，这是大多数MCMC 技术的问题。

论文地址：arxiv/1504.01896.pdf

作者： 亚历山大·贝利

Deephub翻译团队