简单易学的机器学习算法——马尔可夫链蒙特卡罗方法MCMC

版权声明:本文为博主原创文章,未经博主允许不得转载。 https://blog.csdn.net/google19890102/article/details/51755242

对于一般的分布的采样,在很多的编程语言中都有实现,如最基本的满足均匀分布的随机数,但是对于复杂的分布,要想对其采样,却没有实现好的函数,在这里,可以使用马尔可夫链蒙特卡罗(Markov Chain Monte Carlo, MCMC)方法,其中Metropolis-Hastings采样和Gibbs采样是MCMC中使用较为广泛的两种形式。

MCMC的基础理论为马尔可夫过程,在MCMC算法中,为了在一个指定的分布上采样,根据马尔可夫过程,首先从任一状态出发,模拟马尔可夫过程,不断进行状态转移,最终收敛到平稳分布。

一、马尔可夫链

1、马尔可夫链

Xt表示随机变量X在离散时间t时刻的取值。若该变量随时间变化的转移概率仅仅依赖于它的当前取值,即

P(Xt+1=sjX0=s0,X1=s1,,Xt=si)=P(Xt+1=sjXt=si)

也就是说状态转移的概率只依赖于前一个状态。称这个变量为马尔可夫变量,其中,s0,s1,,si,sjΩ为随机变量X可能的状态。这个性质称为马尔可夫性质,具有马尔可夫性质的随机过程称为马尔可夫过程。

马尔可夫链指的是在一段时间内随机变量X的取值序列(X0,X1,,Xm),它们满足如上的马尔可夫性质。

2、转移概率

马尔可夫链是通过对应的转移概率定义的,转移概率指的是随机变量从一个时刻到下一个时刻,从状态si转移到另一个状态sj的概率,即:

P(ij):=Pi,j=P(Xt+1=sjXt=si)

π(t)k表示随机变量X在时刻t的取值为sk的概率,则随机变量X在时刻t+1的取值为si的概率为:

π(t+1)i=P(Xt+1=si)=kP(Xt+1=siXt=sk)P(Xt=sk)=kPk,iπ(t)k

假设状态的数目为n,则有:

(π(t+1)1,,π(t+1)n)=(π(t)1,,π(t)n)P1,1P2,1Pn,1P1,2P2,2Pn,2P1,nP2,nPn,n

3、马尔可夫链的平稳分布

对于马尔可夫链,需要注意以下的两点:

  • 1、周期性:即经过有限次的状态转移,又回到了自身;
  • 2、不可约:即两个状态之间相互转移;

如果一个马尔可夫过程既没有周期性,又不可约,则称为各态遍历的。

对于一个各态遍历的马尔可夫过程,无论初始值π(0)取何值,随着转移次数的增多,随机变量的取值分布最终都会收敛到唯一的平稳分布π,即:

limtπ(0)Pt=π

且这个平稳分布π满足:

πP=π

其中,P=(pi,j)n×n为转移概率矩阵。

二、马尔可夫链蒙特卡罗方法

1、基本思想

对于一个给定的概率分布P(X),若是要得到其样本,通过上述的马尔可夫链的概念,我们可以构造一个转移矩阵为P的马尔可夫链,使得该马尔可夫链的平稳分布为P(X),这样,无论其初始状态为何值,假设记为x0,那么随着马尔科夫过程的转移,得到了一系列的状态值,如:x0,x1,x2,,xn,xn+1,,,如果这个马尔可夫过程在第n步时已经收敛,那么分布P(X)的样本即为xn,xn+1,

2、细致平稳条件

对于一个各态遍历的马尔可夫过程,若其转移矩阵为P,分布为π(x),若满足:

π(i)Pi,j=π(j)Pj,i

π(x)是马尔可夫链的平稳分布,上式称为细致平稳条件。

3、Metropolis采样算法

Metropolis采样算法是最基本的基于MCMC的采样算法。

3.1、Metropolis采样算法的基本原理

假设需要从目标概率密度函数p(θ)中进行采样,同时,θ满足<θ<。Metropolis采样算法根据马尔可夫链去生成一个序列:

θ(1)θ(2)θ(t)

其中,θ(t)表示的是马尔可夫链在第t代时的状态。

在Metropolis采样算法的过程中,首先初始化状态值θ(1),然后利用一个已知的分布q(θθ(t1))生成一个新的候选状态θ(),随后根据一定的概率选择接受这个新值,或者拒绝这个新值,在Metropolis采样算法中,概率为:

α=min1,p(θ())p(θ(t1))

这样的过程一直持续到采样过程的收敛,当收敛以后,样本θ(t)即为目标分布p(θ)中的样本。

3.2、Metropolis采样算法的流程

基于以上的分析,可以总结出如下的Metropolis采样算法的流程:

  • 初始化时间t=1
  • 设置u的值,并初始化初始状态θ(t)=u
  • 重复一下的过程:
    • t=t+1
    • 从已知分布q(θθ(t1))中生成一个候选状态θ()
    • 计算接受的概率:α=min(1,p(θ())p(θ(t1)))
    • 从均匀分布Uniform(0,1)生成一个随机值a
    • 如果aα,接受新生成的值:θ(t)=θ();否则:θ(t)=θ(t1)
  • 直到t=T

3.3、Metropolis算法的解释

要证明Metropolis采样算法的正确性,最重要的是要证明构造的马尔可夫过程满足如上的细致平稳条件,即:

π(i)Pi,j=π(j)Pj,i

对于上面所述的过程,分布为p(θ),从状态i转移到状态j的转移概率为:

Pi,j=αi,jQi,j

其中,Qi,j为上述已知的分布。

对于选择该已知的分布,在Metropolis采样算法中,要求该已知的分布必须是对称的,即Qi,j=Qj,i,即

q(θ=θ(t)θ(t1))=q(θ=θ(t1)θ(t))

常用的符合对称的分布主要有:正态分布,柯西分布以及均匀分布等。

接下来,需要证明在Metropolis采样算法中构造的马尔可夫链满足细致平稳条件。

p(θ(i))Pi,j=p(θ(i))αi,jQi,j=p(θ(i))min1,p(θ(j))p(θ(i))Qi,j=min{p(θ(i))Qi,j,p(θ(j))Qi,j}=p(θ(j))minp(θ(i))p(θ(j)),1Qj,i=p(θ(j))αj,iQj,i=p(θ(j))Pj,i

因此,通过以上的方法构造出来的马尔可夫链是满足细致平稳条件的。

3.4、实验

假设需要从柯西分布中采样数据,我们利用Metropolis采样算法来生成样本,其中,柯西分布的概率密度函数为:

f(θ)=1π(1+θ2)

那么,根据上述的Metropolis采样算法的流程,接受概率α的值为:

α=min1,1+[θ(t)]21+[θ()]2

代码如下:

'''
Date:20160629
@author: zhaozhiyong
'''
import random
from scipy.stats import norm
import matplotlib.pyplot as plt

def cauchy(theta):
    y = 1.0 / (1.0 + theta ** 2)
    return y

T = 5000
sigma = 1
thetamin = -30
thetamax = 30
theta = [0.0] * (T+1)
theta[0] = random.uniform(thetamin, thetamax)

t = 0
while t < T:
    t = t + 1
    theta_star = norm.rvs(loc=theta[t - 1], scale=sigma, size=1, random_state=None)
    #print theta_star
    alpha = min(1, (cauchy(theta_star[0]) / cauchy(theta[t - 1])))

    u = random.uniform(0, 1)
    if u <= alpha:
        theta[t] = theta_star[0]
    else:
        theta[t] = theta[t - 1]

ax1 = plt.subplot(211)
ax2 = plt.subplot(212) 
plt.sca(ax1)
plt.ylim(thetamin, thetamax)
plt.plot(range(T+1), theta, 'g-')
plt.sca(ax2)
num_bins = 50
plt.hist(theta, num_bins, normed=1, facecolor='red', alpha=0.5)
plt.show()

实验的结果:

这里写图片描述

对于Metropolis采样算法,其要求选定的分布必须是对称的,为了弥补这样的一个缺陷,在下一篇中,介绍一下Metropolis-Hastings采样算法,其是Metropolis采样算法的推广形式。

参考文献

没有更多推荐了,返回首页