机器学习算法——马尔可夫链蒙特卡罗（MCMC）

什么是MCMC，什么时候使用它

MCMC只是一种从分布中抽样的算法。这个术语代表“马尔可夫链蒙特卡罗”，因为它是一种使用了“马尔可夫链”的“蒙特卡罗”（即随机）方法。MCMC只是一种蒙特卡罗方法。

为什么我要从分布抽样呢?

从分布抽样是解决一些问题的最简单的方法。

也许在贝叶斯推断中最常见的方式是使用MCMC来从某些模型的后验概率分布中抽取样本。有了这些示例，你就可以问这样的问题：“什么是参数的平均值和可信区间？”。

例如，假设你有合适的参数模型的后验概率密度是某个函数f(x,y)。然后，计算参数 x的平均值，你可以这样计算

你可以简单地读作“x乘以参数(x,y)的概率，积分在x和y可能采用的所有可能值上。”

另一种计算这个值的方法是模拟观察值K ，，从f(x,y)计算样本平均值为

其中x(j)为第j个样本的x的值。

如果这些样品是分布的独立样本，那么随着k→∞，x的估计平均数会收敛于真平均数。

假设我们的目标分布是一个正态分布，平均数m和标准偏差S。显然，这种分布的平均值是m，让我们试着通过从分布中抽取样本来显示。

举一个例子，估计一个平均值m和标准偏差S的正态分布的值（在这里，我只会使用与普遍的常态相对应的参数）：

我们可以使用 rnorm 函数很容易地从这个分布中抽样：

样本的平均值与真平均数（零）非常接近：

事实上，在这种情况下，n个样本的估计的期望方差是1/n，所以我们预期最值位于真平均数为10000点的

中。

此函数计算累积平均数（即元素1,2，…，K的和除以K）。

这里是收敛到接近于真平均值（红线为0）。

将X轴变换到对数刻度上，并显示另外30个随机方法：

如何得出结果？考虑积分

如果这可以分解成函数f（x）和概率密度函数P（x）的乘积，则

注意，右边仅仅是期望值E[f(x)]。根据大数定律，当样本量增长到无穷大时，期望值是样本均值的极限。因此，我们可以近似E[f(x)]为

你可以用这种方法做很多类似的事情。例如，如果你想在估计值附近画一个95%的可信区间，你可以通过求解a来估计它的底部分量。或者，你可以从你的一系列采样点中取样本分位数。

下面是分析地计算概率密度为2.5%的点：

在这种情况下，我们可以直接使用积分来估计（使用上述参数）

并且用Monte Carlo积分估计点：

请注意，这有一个错误：

但当样本大小趋于无穷大时，这将收敛。如果我们重复采样100次，我们会得到一系列估计，其具有与平均值周围的误差大约相同量级的误差：

仍然不相信？

在贝叶斯框架中，你将计算出在所有其他参数中你感兴趣的参数的临界分布。如果你有其他50个参数，这是一个很难的积分！

为什么这是很难的，考虑参数空间区域中包含“有趣”的参数值：即参数值都明显大于零的概率。

为了说明这个问题，

l.假设一个半径为R的圆的和一个边长为2R的正方形；空间的“有趣的”区域为，所以我们有一个很好的机会，随机选择点在圆内。

2.在一个边长为2R的立方体内有一个半径为R的球体，球体的体积是和立方体的体积，所以的体积是“有趣的”。

3.d作为问题的维度，增加（使用超立方体中的超球面）这个比值为

因此，假设一个有零协方差项，在原点的平均值和单位方差的多元正态分布。这些在原点有具有不同的模式（最大值），并且点和模式处的概率比为

关于这个函数，任何大的值会导致概率低。所以，随着问题的维数增加，有趣的空间变得非常小。考虑在区域内采样，并且计数10,000个采样点中有多少具有大于1/1000的相对概率，其像超球面情况一样下降。

为什么“正常统计”不使用蒙特卡罗方法？

对于传统讲授的统计中的许多问题，你可以最大化或最大化一个函数而不是从分布中取样。所以我们使用一些函数来描述可能性，并最大限度地提高它（最大似然推理），或一些函数计算平方和和最小化。

为了避免必须从分布中抽样，在概率估计的误差估计倾向于是渐近的大数据估计或者可能基于Bootstrap估计。

然而，蒙特卡罗方法在贝叶斯推断中所起的作用与频率统计中的优化程序相同;它只是做的推理算法。

马尔可夫链蒙特卡罗

定义

设表示在t时间一些随机变量的值。马尔可夫链起始于产生一系列样点，接着一系列随机步骤。

马尔可夫链满足马尔可夫属性。马尔可夫属性是已知现在状态的条件下，将来所处的状态与过去状态无关，即“遗忘性”；基本上不管你如何达到某个状态x，x的转移概率不变：

从一个步骤到下一个步骤的转换是由转移核描述，它可以由状态i到状态j的转变的概率（或连续变量的概率密度）描述为

表示状态 j 在 t 时间（步骤）的链的概率，并且定义为可能的状态的概率的向量。然后，给出，我们利用 Chapman-Kolmogorov 方程计算。这个概率是我们在状态 k 乘以从 k 到 i 的转变的概率，在所有可能的状态 k上求和。设 P是概率转移矩阵 ——矩阵的第i,j元素是P(i→j)，并且将上述等式重写为

我们可以轻松地迭代这个方程：

平稳分布

如果有一些矢量满足

那么$\vec\pi^$是这个马尔可夫链的平稳分布。直观地，认为这个系统将设置的状态的最终特征集；运行足够长的时间，该系统会“遗忘”它的初始状态，则这个向量的第i个元素是这个系统将在状态i的概率。

如果这个过程是不可约的和非周期性的，马尔可夫链将会有一个平稳分布。

在数学上，是特征值为1的左特征向量。

这里有一个快速的定义，使事情更具体（但请注意，这无关MCMC本身–这只是对马尔可夫链的观点！）假设我们有一个三态马尔可夫过程。设P为链的概率转移矩阵：

请注意，P行的共计为一：

入口P[i,j] 给出了从状态i到状态j的概率（所以这就是上述的 P(i→j)）。

注意，与行不同，列不一定共计为1：

这个函数取一个状态向量x（其中x[i]是处于状态i的概率），并通过乘以转移矩阵P来迭代它，使系统前进n个步骤。

从状态1中的系统开始（所以x是向量[ 1,0,0 ]，指示存在处于状态1的100％概率，并且没有处于任何其他状态的机会），并且迭代10个步骤：

同样，对于其他两个可能的起始状态：

这表明在平稳分布上的收敛性。

这意味着不管起始分布如何，不管它在哪里开始，在大约10次或更多次迭代之后，链处于状态1的概率为32％。

我们可以用R的特征函数来提取系统的主导特征向量（这里的t()转置矩阵得到左特征向量）。

然后添加点到之前显示我们接近收敛的图中：

根据特征向量的定义，乘以特征向量的转移矩阵返回特征向量本身：

（严格地说，这应该是由V相乘的特征值，但这些矩阵的主特征值总是1）。

在这里运行的函数用一个状态（这一次，只是一个整数代表其中系统的状态1,2,3），如上相同的转移矩阵，和一些运行的步骤。每个步骤，查看可能转换到的地方，并选择1（这使用R的样本函数）。

这是100步左右的链：

绘制我们在每个状态随着时间的推移的时间分数。

运行这步再长一点（5000步）

存在稳定分布的充分（但不是必要的）条件是细致平衡，其为：

这意味着链是可逆的。这种情况意味着一个平稳分布存在的原因是

将状态j的细致平衡方程的两侧求和

左边的项等于的第k个元素，右边的项可以作为因子

然后，因为因为P是一个转移概率函数，由总概率定律决定，概率为1）所以右边是

，所以我们有

其适用于所有k，所以

马尔可夫链具有平稳分布，如果我们运行它们足够长的时间，我们可以看看链在哪里花费时间，并得到该平稳分布的合理估计。

Metropolis算法

这是最简单的MCMC算法。我们要做的是具有一些我们想要采样的分布，并且我们将能够评估与目标分布的概率密度成比例的函数f(x)（也就是说，如果 p(x)是概率密度函数本身，，即f(x)= p(x)/ Z，其中Z =∫f(x)dx。注意，x可以是向量或标量。

我们还需要一个概率密度函数P，我们可以从中抽取样本。对于最简单的算法，建议分布是对称的，即P(x→x’)= P(x'→x)。

该算法如下进行。

1.在一些状态下启动t。

2.创建新状态x'

3.计算“接受概率”α= min[1,f(x’)f(x)]

4.从[0,1]中绘制一些均匀分布的随机数u; 如果u<α，接受该随机数，设置$x{t+1} = x^\prime。否则拒绝它并设置x{t+1}=x_t$。

注意，在上面的步骤3中，未知的归一化常数丢失，因为

这将生成一系列样本$ {x0，x1，\ldots}$。在所提出的样本被拒绝的情况下，相同的值将存在于连续样本中。

这些不是来自目标分布的独立样本，它们是依赖样本。也就是说，样本取决于，等等。然而，因为链接近稳定分布，所以只要我们采样足够的点，这种依赖就没有关系。

MCMC采样在1d（单参数）问题

这是从中抽样的目标分布，它是两个正态分布的加权和。概率密度函数是

这是一个人为的例子，但是这样的分布不是完全不可能的，并且可能出现在从混合物中采样东西时。

这里有一些参数和目标密度的定义。

这里是概率密度绘制在“重要”领域的一部分

让我们定义一个非常简单的建议的算法，从一个正态分布为中心的电流点的标准偏差为4

这实现了核心算法，如上所述：

这只是运行MCMC的一些步骤。它将在点x处开始返回一个具有nsteps行的矩阵，并且具有与x相同的列数。如果在标量x上运行，它将返回一个向量。

我们选择一个地方开始（-10好了）

这里是马尔可夫链的前1000步，目标密度在右边：

即使只有一千个（非独立的）样本，我们开始类似的目标分布。

运行更长时间会看起来更好一些：

现在，运行不同的建议机制——一个具有非常宽的标准偏差（33个单位），另一个具有非常小的标准偏差（3个单位）。

这里是和上面相同的图——注意三个轨迹移动的不同的方式。

灰线轨迹相当自由地弹跳。

相反，红色轨迹意味着它倾向于长期呆在一起的时间。

蓝色轨迹倾向于小移动，它对于大多数轨迹随机行走而移动。它需要数百次迭代甚至达到大部分的概率密度。

你可以在后续参数之间的自相关中看到不同建议步骤的效果——这些图表示不同滞后的步长之间的自相关系数的衰减，蓝色线表示统计独立性。

由此，可以计算独立样本的有效数量：

两个“混合”的链比第一个更差。

这更清楚地显示了链运行更长时间会发生什么：

显示100，1,000，10,000和100,000步：

MCMC在二维中

这是一个给定装置（分布的中心）和方差 - 协方差矩阵的向量的多变量法向密度的函数。

如上所述，将目标密度定义为两个mvns的总和（此时未加权）：

这里有一系列不同的策略——我们建议同时在两个维度上的移动，或者我们可以独立地沿着每个轴进行采样。这两种策略都会奏效，但是它们的组合速度会有所不同。

假设我们实际上不知道如何从mvn中抽样，让我们做一个方案，在两维分布均匀，从宽度为d的正方形采样。

绘制样品

将抽样分布与已知分布进行比较：

然后，我们可以轻松地用难以直接做的样本做事情。例如，参数1的边际分布是什么：

（这是第一参数采用的分布，对第二参数可能采用的所有可能的值求平均，由它们的概率加权）。

正确地计算这一点是很复杂的，我们需要将第一个参数的第二个参数的所有可能的值合并到第一个值中。然后，因为目标函数本身不是归一化的，所以我们必须将它除以第一维上的积分值（这是分布下的总面积）。

文章原标题《Markov Chain Monte Carlo》，作者：Rich FitzJohn