随机过程--Metropolis-Hastings算法

随机过程–Metropolis-Hastings算法

蒙特卡罗方法

  蒙特卡罗(Monte Carlo)方法又称随机抽样或统计试验方法，简单地理解就是利用随机数去解决许多计算问题，通过实验去求解一些数学问题。通常是通过一些随机模拟实验去求解一些概率或者期望的问题。

生成随机数

  我们知道，其实计算机只能产生均匀分布的伪随机数，但很多时候，我们希望得到一连串服从一定分布（如高斯分布，泊松分布）的随机数。我们就要想办法把均匀分布的数映射到服从一定分布的数，这样就通过模拟生成一些服从一定分布的随机数来解决一些复杂的数学问题。

求解概率和期望问题

1. 求解概率问题 
     贝叶斯先验转后验： 
   
   p(θ|D)=p(D|θ)p(θ)Dp(θ|D)=p(D|θ)p(θ)D
   
      概率归一化问题： 
   
   p(θ)=p˜(θ)Zp(θ)=p~(θ)Z
   
      通常，我们需要从贝叶斯后验和归一化后的概率中抽样，但是通常这两个概率很难求，主要是因为分母很难求，分母可能涉及到复杂的积分计算。因此，我们无法通过简单的数学公式将均匀分布的随机数映射到服从该分布的随机数。
2. 求解期望问题 
     如上所述，如果一个概率的表达式没法求出，我们无法求得期望。但是如果我们能生成许多连续的服从该分布的随机数，根据切比雪夫大数定理就能通过简单地加和近似地求得概率。 
   切比雪夫大数定理： 
   
   limn→∞p(|1S∑s=1Sf(xs)−1S∑s=1SE[f（X）]|<ϵ)=1limn→∞p(|1S∑s=1Sf(xs)−1S∑s=1SE[f（X）]|<ϵ)=1
   
   计算期望： 
   
   E[f（X）]=∫f（x）p（x）dx≃1S∑s=1Sf(xs)E[f（X）]=∫f（x）p（x）dx≃1S∑s=1Sf(xs)
   
   
   xs∼p(X)xs∼p(X)

栗子

估计圆周率：我们可以在一个包含圆的正方形中随机生成数据点，然后统计一下落在圆内的点的个数与总个数的比例，就可以大约估算出圆周率的值。

还有一个比较出名的求圆周率的实验就是蒲丰投针，这是1777年法国科学家布丰提出的一种计算圆周率的方法：

1. 取一张白纸，在上面画上许多条间距为a的平行线。
2. 取一根长度为l（l=a/2） 的针，随机地向画有平行直线的纸上掷n次，观察针与直线相交的次数，记为m。
3. 计算针与直线相交的概率：p=2lπap=2lπa。
4. 可以计算得到π=2lapπ=2lap 
   

马尔可夫链

  马尔可夫链(Markov Chain )简单地来说就是状态机，能实现状态之间的转移，下一个状态只与当前的状态有关，与之前的状态没有关系。 

p(xN|x1,…,xN−1)=p(xN|xN−1)p(xN|x1,…,xN−1)=p(xN|xN−1)

现实生活中也有很多事情的变化规律是形成一个马尔可夫链的： 

1. 蛋白质分子的变化过程：蛋白质的下一个状态只与它现在的结构有关。
2. 股市行情：下一段时间股市的行情只与当前的股价有关。
3. 赌徒赌资：一个赌徒经过下一个赌局后的资金只与他现在所拥有的资金有关。
4. 交通拥堵状态：一个路口下一个时间点的拥堵状况只与当前的交通状态有关。

现在让我们来举一个栗子： 
  假设现在社会根据人民所拥有的资产分为上层，中层，下层，且它们之间可以转换，也就是每一层的人民都有一定的几率变成其他层，它们之间的跳转概率为：

-	下层	中层	上层
下层	0.65	0.28	0.17
中层	0.15	0.67	0.18
下层	0.12	0.36	0.52

它们之间的状态转移图为： 

现在，我们要计算第N代后，整个社会中各阶层人数的比例：

1. 计算概率转移矩阵
   P=⎡⎣⎢0.650.150.120.280.670.360.170.180.52⎤⎦⎥P=[0.650.280.170.150.670.180.120.360.52]
2. 初始化第0代各阶层的人数比例： 
   
   π0=[a,b,c],a+b+c=1π0=[a,b,c],a+b+c=1
3. 不断地迭代更新知道收敛： 
   
   πN=πN−1∗P=πN−2∗P2=...π0∗PNπN=πN−1∗P=πN−2∗P2=...π0∗PN

我们随机初始化，且做了A与B两组实验，实验结果分别如下： 
A组

代数	下层	中层	上层
0	0.210	0.680	0.110
1	0.252	0.554	0.194
2	0.270	0.512	0.218
3	0.278	0.497	0.225
4	0.282	0.490	0.226
5	0.285	0.489	0.225
6	0.286	0.489	0.225
7	0.289	0.489	0.225
8	0.286	0.488	0.225
9	0.286	0.489	0.225

B组

代数	下层	中层	上层
0	0.750	0.150	0.100
1	0.522	0.347	0.132
2	0.407	0.426	0.167
3	0.349	0.459	0.192
4	0.318	0.475	0.207
5	0.303	0.482	0.215
6	0.295	0.485	0.220
7	0.291	0.487	0.222
8	0.289	0.488	0.225
9	0.286	0.489	0.225

  根据实验结果，我们可以发现，就算给予不同的初始化，最终收敛后的结果都是相同的。当马尔科夫链满足细致平稳性就会有这样的结果。

细致平稳性（detail balance condition）：当非周期性的马尔可夫链的概率转移矩阵和每一个状态的概率满足： 

π(i)p(j|i)=π(j)p(i|j)π(i)p(j|i)=π(j)p(i|j)

最终得到的状态 ππ 是该马尔可夫链的平稳分布。

证明: 

∑i=1∞π(i)p(j|i)=∑i=1∞π(j)p(i|j)=π(j)∑i=1∞p(i|j)=π(j)∑i=1∞π(i)p(j|i)=∑i=1∞π(j)p(i|j)=π(j)∑i=1∞p(i|j)=π(j)

Metropolis算法

我们先介绍提议分布，其实提议分布可以随意给，通常我们可以假设提议分布服从高斯分布： 

q(j|i)=N(j|i,σ)q(j|i)=N(j|i,σ)

但是这样的话，通常不能满足细致平稳性，即： 

π(i)q(j|i)≠π(j)q(i|j)π(i)q(j|i)≠π(j)q(i|j)

  为了满足细致平稳性，我们要让上面的不等式的不等号变成等号，方法其实很简单，就是让左边乘上右边，右边乘上左边，乘上的表 
达式被称为接受概率（acceptance probability） αα ： 

α(j|i)=π(j)q(i|j)α(j|i)=π(j)q(i|j)

π(i)q(j|i)α(j|i)=π(j)q(i|j)α(i|j)π(i)q(j|i)α(j|i)=π(j)q(i|j)α(i|j)

  现在我们就构造出了满足马尔科夫细致平稳性地马尔可夫链，我们可以把q(j|i)α(j|i)q(j|i)α(j|i)视为转移概率。 

Metropolis-Hastings算法

  但是Metropolis算法构造出的接受概率可能会很小，这样造成算法要经过很多的迭代才能到达平稳分布。为了满足细致平稳，不等式的两边都乘以了一个小的接受概率，那我们可以把其中一个接受概率乘以一个数变为1，另外一边的接受概率也乘上相同的倍数，具体的做法如下： 

π(i)q(j|i)α(j|i)=π(j)q(i|j)α(i|j)π(i)q(j|i)π(j)q(i|j)=π(j)q(i|j)π(i)q(j|i)π(i)q(j|i)π(j)q(i|j)π(i)q(j|i)=π(j)q(i|j)orπ(i)q(j|i)=π(j)q(i|j)π(i)q(j|i)π(j)q(i|j)π(i)q(j|i)α(j|i)=π(j)q(i|j)α(i|j)π(i)q(j|i)π(j)q(i|j)=π(j)q(i|j)π(i)q(j|i)π(i)q(j|i)π(j)q(i|j)π(i)q(j|i)=π(j)q(i|j)orπ(i)q(j|i)=π(j)q(i|j)π(i)q(j|i)π(j)q(i|j)

  整合上面等式可以得到： 

α(j|i)=min{π(j)q(i|j)π(i)q(j|i),1}α(j|i)=min{π(j)q(i|j)π(i)q(j|i),1}

具体的算法流程如下：

1. 当前的状态i
2. 从提议分布q(j|i)q(j|i)中产生新的状态j（如果是简单的分布，如高斯分布，可以通过将均匀分布的值通过高斯分布的反函数映射得到服从高斯分布的值。或者使用其他采样方法采样得到，如重采样与拒绝采样。）
3. 计算接受概率： 
   
   α(j|i)=min{π(j)q(i|j)π(i)q(j|i),1}α(j|i)=min{π(j)q(i|j)π(i)q(j|i),1}
4. 从[0，1]区间中抽样出一个数据点（均匀分布）u∼U[0,1]u∼U[0,1]

5. 将随机生成的数和接受率对比，如果u<αu<α,状态就从i跳转到j，否则继续保留在状态i。

     让我们举个栗子，我们抛两枚筛子，计算出两个筛子正面的值的和为x的概率分布是多少，虽然这个很简单，我们理论计算都能得到，为了说明MH算法嘛。在这个过程中，我们利用MH方法不断抽样，最终再统计各个结果的个数。并且在这个过程中，我们没有使用每个结果出现的概率（知道的话我们就没必要抽样计算了），而是知道每个结果出现的组合数。当然某个结果出现的组合数除以总的组合数就是概率，这个例子比较简单，在很多复杂情况下，这个总数，也就是分母，并不是这么好求，很多要涉及到无限积分。 
   

    
   实验代码： 
    
   实验的结果：随着模拟次数额增大，抽样结果越来越接近理论结果。 
   

     让我们再来举个栗子，我们要从一个未归一化（归一化涉及到无限积分）的概率分布f(x)=0.5x2e−xf(x)=0.5x2e−x中抽样出满足这个分布的数据点。 
   实验的代码是： 
   

    
   实验结果为如下图所示，其中左上角的图片是真实的概率分布，当抽样的点多的时候，抽样的点的分布与真实的分布相似。