简单易学!一步步带你理解机器学习算法——马尔可夫链蒙特卡罗(MCMC)

  1. 云栖社区>
  2. 博客>
  3. 正文

简单易学!一步步带你理解机器学习算法——马尔可夫链蒙特卡罗(MCMC)

【方向】 2016-12-22 00:16:56 浏览29334
展开阅读全文

本文由北邮爱可可-爱生活老师推荐,阿里云云栖社区组织翻译。

以下为译文:

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

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

为什么我要从分布抽样?

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

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

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

92efeb2b155e2454e5f8dc717079edbc8303128b

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

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

e8f8d0d3e5837addd0febac259a9cc3f0c4835ae

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

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

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

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

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

ae60e37c500e51eedd5e6a98497b8bc107aa58f1

样本的平均值与真平均数(零)非常接近:

2dc5683223b0959f7647d410d044b990b7afd645

1e92db7973f65383aadfb2bb58445a0f72e1be9d

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


中。

5ae3f4d1dce738eb691c98d4f62ee2e55e518049

4ece26f3474cffd61b2bbdfcbb4e21cd19065f41

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

146d6d88ea1f9deaf0c2cc06dbac2ac44b5bd7fc

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

67c65775d088a8c0356b10e2d134a69cd0d8728a

60da4e7164d8ba0b43198c9162d61ed54e258398

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

027dcdfa2792c66b42fddf6945bd7f9d879a2250

775bff0b3c01b397ced7e5f4c1750ef4dc0e23ea

如何得出结果?考虑积分

1ac554f2961db9e9a1021582ca4cd09a7c42aec8

如果这可以分解成函数f(x)和概率密度函数P(x)的乘积,则

2bc4e265939383a8b669a7ae30912de9a186b0df

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

fb7d0df5319cd615f14a78b1358cbe40dcd8efce

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

32648b0497939f63e27c470da91df6752d259f30

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

2b86073a41ef569ad4aa5d119331371b30327c24

086f46ba0d38449d4c321e9100991724c4384b18

在这种情况下,我们可以直接使用积分来估计(使用上述参数)

0ebef837c16cf40c8f91f03e50f0ea5b275c6eda

ab3631f26721faec31f12585b2a1736e8a695dbc

并且用Monte Carlo积分估计点:

a18a1c6b2154fc7794d81c4fdac0b448178393b3

abc38d62fd5eb5893addd1f060e7399b35fb5543

请注意,这有一个错误:

34cbb9136a305bb343443c9275671f0cb32b11fb

3b6bece0364d06f6026a3f8368da50a33a4fae93

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

fb77093d3040f56e8ab6a5988613aa23405bf95f

6e7424df2cdfc5dd00b1f5dabbb2af07a192e3b8

仍然不相信?

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

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

为了说明这个问题,

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

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

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

fd22d41f8eddaef85ccc21ccf6b63fa3e44a698e

700e23946991493ac5cc6a4e4e6174939012d70784c8ff3e84039a87c22ad5fa0bea87241ace6899

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

0f8463e50208497117dbdf00484d02a365d36412

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

e7db32afb566ccd770ce35875ea323c168302d61

2ee73b6e5973477c5a578ced46c819ea7e277aac

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

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

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

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

马尔可夫链蒙特卡罗

定义

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

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

381b38d158c5068fbe23f0218dc82c6651e7f12a

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

d49a03c424b6428f453d40d693ab0041fa7e25de

76034340b4b4e683b187614edf815efc209088a6表示jt时间(步骤的链概率并且定义07068ec9d8d9e47f051b6741b1e1288707d5774a可能状态的概率的向量。然后,给出07068ec9d8d9e47f051b6741b1e1288707d5774a我们利用Chapman-Kolmogorov方程计算fd4f7cfb832cb15d324e27743f1e85285d269e62

41720fed885ea1d1235a04b62bf5bbf9f05577bf个概率是我们在状态k乘以从ki转变的概率,在所有可能的状态k上求和。P是概率转移矩阵——矩阵的第i,j元素是P(i→j),并且将上述等式重写为 b1549831d0eb5b62f25e8ce233d04e3de7631ee6

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

ecf0ffc7bb02109fcea859355f22b1bd00b568f4

平稳分布

如果有一些矢量67d5b47fd468c6e84e14df3830e4baf71696eefe满足

a8d6f009876aa6ea37d194468bf29af84bb29c0f

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

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

在数学上,67d5b47fd468c6e84e14df3830e4baf71696eefe特征值为1的左特征向量。

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

dc52e3991da20534e7b7ad147f2ce365673fb367

3c25576d712c98602ab44c04008fe3ab8e2da29d

请注意,P行的共计为一:

83ba3747a99f3d41fe19dc650db9174718f7a9af

43921ebd7eb771b0f391b061db1a5fba65a7433d

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

注意,与行不同,列不一定共计为1:

1a169ca63c2e213969fc3c37953152204bf0ef1d

e44ef537f093c4ec6fe2e93dd7bd4ad0b1535dc1

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

371a2b484436c5f230d92bde763b450fb036ea6f

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

8ca1cb7194800de768e54ed8f0c2454e13f31fb2

同样,对于其他两个可能的起始状态:

587dc1a380cb6f48740ddc8b8b365d85debf4073

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

7afc4c12245f6d2ea43d02b7d221d07991ece2bd11fdaafd182588f887532bf2a2d372de54a4977d

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

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

6df7b87feb41537472154bffe8abe4238874f181

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

05fcb68304e9089ee7fd4a482e33c84c3847e7ea

5a023d6ae98e9eb9ae72d63e108d8d2db5f2bb15

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

48b860fe2baa8e5fcaf25ac0f21cc862e0bc6c88

a64c0464e5a5000ee00e08ec1e85f2c3056f4805

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

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

7b5f41c5dc8c3ce560220f3240a33cc8fca466cd

这是100步左右的链:

f1efe6e276586deb97401df7977c67eed4f4c50c

6240eb59da9b5cc4d45f3b4493e3c13ab5c766da

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

aeb79753f120c5aaaf63a4f7a310335dd41c97ac

3a39ad0f38d87acd3e8accecf91834288e3cffc8

运行这步再长一点(5000步)

c56f498964ee8a478815e0d27c6890d4aca61c2d

aa6ad36ad506b257a53a2df1d33440390c4d7016

存在稳定分布的充分(但不是必要的)条件是细致平衡,其

dd8d0935af75ef6de80be93d2acab8e9ec1d80d0

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

a8d6f009876aa6ea37d194468bf29af84bb29c0f

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

6dc39dfac306b31bf3280a596927cb8b7a3a87a8

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

b929e58145a44bc8eee3cf7b9f8ea3a370838e64

然后,因为5c9470bbd762076fdeca941956dc670dacbe5bc0因为P是一个转移概率函数,由总概率定律决定,概率为1)所以右边是ad1b30780c6047e172271813de74b444f8687275

,所以我们有

9d58224bb33278916409bd702cc398d9da0a232d

其适用于所有k,所以

4da740cf6521d69160f7a5e8ba01d06d8c8f2d85

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

Metropolis算法

这是最简单的MCMC算法。我们要做的是具有一些我们想要采样的分布,并且我们将能够评估与目标分布的概率密度成比例的函数f(x)(也就是说,如果 p(x)是概率密度函数本身,bac366b7391b7ed2f4496000a1aa366627df79c3,即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中,未知的归一化常数丢失,因为

1ecd9a30cc07d58b130954db23d0e5ccdb84c58b

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

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

MCMC采样在1d(单参数)问题

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

3def1e4a168db3375ba064975abee63cea8d2c0f

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

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

28c3ec6f3421a4c618f782eac5e2792857f6345d

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

6d506fd9b49073939f622f119118b832180088cb283caac2b7ae6eddc1cc95170b3e3a523cd2d606

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

d56a24858dcc7531d770acc6344e34e44f531a6c

这实现了核心算法,如上所述:

9ba8ddb558eaea7dbd6dc68b54efb0c0e665fe47

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

e19769120c9275cec17bcc4a85166165721e8cbc

我们选择一个地方开始(-10好了)

c89bbfc8e63b17e9cabc2fd4cdb0a48c31f2ad14

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

b04f6e2a360a3848387b3cc1b5074550640f86f7

5c2aadf0ac33f26fbbc7e7d5644fd41d7323f656

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

c1b898e6472c66ca9e402de8c453f261199ece36

98952b6bd1fdfa8b00fdba58979d9bb044b32ab2

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

9c6933c4c59cfc53bdf4b9090312e2af253a91d8

a62a836088274dabc24022540c37c6b24cabef51

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

e2579269aad4633d143487a418fb31857f95dc67

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

a689eed9fc7c0f2a03a1b7850cabc4fb60290fb6

d26b870f61f4d9f08298bc92a187d00f5a0553b7

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

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

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

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

e72cd09f5d103088e4d062f07b7b801324565007

16c1a3f6bbb6af20ecc26202a320280d857b88c5

由此,可以计算独立样本的有效数量:

039007998c8460bcfaf5562754e4bb42f6585479

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

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

08cf83bf3cc3987f4e28cd6869b6693c58f0f4a8

显示100,1,000,10,000和100,000步:

3436583a2c6071ac5430385564352c8eb7b5f20e

3f48400c2e453fac8c8875e030c7bda6aaa49313

MCMC在二维中

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

593bb90bd2d74fd75d049287dc65e14abc0ce47f

如上所述,将目标密度定义为两个mvns的总和(此时未加权):

6bb1269aee427d8cc6ed0e2624fb8c89f32e79b3

fd7eac84c7458083a4eb74c6f1a78dc769a4772e

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

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

fab8058dbf162401d276e740f52eea985df0640f

a0fc969a5666325678c4ffd86ff042ac2476c9e8

绘制样品

2793d30553cce29ede2487a93ac19df319e72ab0

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

08328da64d6c9da73856d7092389a59cd8c19d69

f801dc5458e841f08b4e030262f3890b36adb218

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

279effc94ca90ff40fe7a62978cebac737e7b885

8879fda695e51942c6951f84b0917abf814b526a

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

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

44f24e5aebfac3613813bc8a294bb7740f693f44

aaa423e37f5647583a19f3b74f695e292802553c

数十款阿里云产品限时折扣中,赶紧点击领劵开始云上实践吧!

文章原标题Markov Chain Monte Carlo》,作者:Rich FitzJohn 文章为简译,更为详细的内容,请查看原文

网友评论

登录后评论
0/500
评论
【方向】
+ 关注