美文网首页
马尔科夫链蒙特卡罗法

马尔科夫链蒙特卡罗法

作者: 单调不减 | 来源:发表于2019-06-22 21:50 被阅读0次

蒙特卡罗法(Monte Carlo Method)也称为统计模拟方法,是通过概率模型的随机抽样进行近似数值计算的方法。

马尔科夫链蒙特卡罗法(Markov Chain Monte Carlo,MCMC)则是以马尔科夫链为概率模型的蒙特卡罗法。MCMC构建一个马尔科夫链,使其平稳分布就是要进行抽样的分布,首先基于该马尔科夫链进行随机游走,产生样本的序列,之后使用该平稳分布的样本进行近似数值计算。

1、蒙特卡罗法

在我浅显的认知里,蒙特卡罗法就是在正方形里随机取点用概率近似求内接圆面积的近似方法。看过《统计学习方法》后的简单介绍后发现自己还是naive,蒙特卡罗法还是很强大的。

蒙特卡罗法要解决的问题是:假设概率分布定义已知,通过抽样获得概率分布的随机样本,并通过得到的随机样本对概率分布的特征进行分析。蒙特卡罗法的核心是随机抽样。

一般的蒙特卡罗法有直接抽样法、接受-拒绝抽样法、重要性抽样法等。后两者适用于概率密度复杂、不能直接抽样的情况。

接受-拒绝法的基本思想如下:假设p(x)不可直接抽样,找到一个可以直接抽样的分布,称为建议分布。假设q(x)为建议分布的概率密度函数,并且有cq(x)\geq p(x),其中c>0。按照q(x)进行抽样,假设得到结果是x^*,再按照\frac{p(x^*)}{cq(x^*)}的比例决定是否接受x^*。如下图所示,落到p(x^*)范围内就接受(绿色),否则拒绝(红色)。

具体算法流程如下:

(1)选择概率密度函数为q(x)的概率函数作为建议分布,使其对任一x满足cq(x)\geq p(x),其中c>0

(2)按照建议分布q(x)随机抽样得到样本x^*,再按照均匀分布在(0,1)范围内抽样得到u

(3)若u\leq\frac{p(x^*)}{cq(x^*)},则将x*作为抽样结果,否则返回(2)。

(4)直至得到n个随机样本,结束。

接受-拒绝法实际是按照p(x)的涵盖面积占cq(x)的涵盖面积的比例进行抽样。其优点是容易实现,缺点是效率可能不高。若p(x)涵盖面积占cq(x)的涵盖面积的比例很低,就会导致拒绝比例很高,抽样效率很低。

一般蒙特卡罗法可以用于数学期望估计。按概率分布p(x)独立地抽取n个样本x_1,x_2,\dots,x_n,之后计算函数f(x)的样本均值\hat{f}_n

\hat{f}_n=\frac{1}{n}\sum_{i=1}^n f(x_i)

作为数学期望E_{p(x)}[f(x)]的近似值。这里可以通过离散的情况简单的理解,比如x=x_i被抽到m次,把m\frac{1}{n}相乘就得到x=x_i的频率,可用于估计概率p(x_i),从而上式可视为\sum_x f(x)p(x)即期望E_{p(x)}[f(x)]的近似值。

此外,一般蒙特卡罗法还可以用于积分计算。假设有一个函数h(x),目标是计算该函数积分:

\int_{x}h(x)dx

h(x)分解成一个函数f(x)和一个概率密度函数p(x)的乘积,则有:

\int_x h(x)dx=\int_x f(x)p(x)dx=E_{p(x)}[f(x)]

于是函数h(x)的积分就可以表示为函数f(x)关于p(x)的期望。实际上给定一个概率密度函数p(x),只要取f(x)=\frac{h(x)}{p(x)},即可得到上述形式。接下来,我们用期望的估计即可计算积分:

\int_x h(x)dx=E_{p(x)}[f(x)]\approx \frac{1}{n}\sum_{i=1}^n f(x_i)

2、马尔科夫链

马尔科夫链定义:

考虑一个随机变量的序列X=\{X_0,X_1,\dots,X_t,\dots\},这里X_t表示时刻t的随机变量,t=0,1,2,\dots。每个随机变量X_t(t=0,1,2,\dots)的取值集合相同,称为状态空间,表示为S。随机变量可以是离散的,也可以是连续的。以上随机变量的序列构成随机过程

假设在时刻0的随机变量X_0遵循概率分布P(X_0)=\pi_0,称为初始状态分布。在某时刻t\geq 1的随机变量X_t与前一时刻的随机变量X_{t-1}有条件分布P(X_t|X_{t-1})。若X_t只依赖于X_{t-1}而不依赖于过去的随机变量\{X_0,X_1,\dots,X_{t-2} \},这一性质称为马尔可夫性,即:

P(X_t|X_0,X_1,\dots,X_{t-1})=P(X_t|X_{t-1}),\quad t=1,2,\dots

具有马尔科夫性的随机序列X=\{X_0,X_1,\dots,X_t,\dots\}称为马尔科夫链马尔科夫过程。条件概率分布P(X_t|X_{t-1})称为马尔科夫链的转移概率分布。

若马尔科夫链在时刻(t-1)处于状态j,在时刻t移动到状态i,将转移概率记作:

p_{ij}=(X_t=i|X_{t-1}=j),\quad i=1,2,\dots,\quad j=1,2,\dots

满足:

p_{ij}\geq 0,\quad \sum_i p_{ij}=1

马尔科夫链的转移概率p_{ij}可用矩阵表示:

P= \begin{bmatrix} p_{11} & p_{12} & p_{13} & \dots\\ p_{21} & p_{22} & p_{23} & \dots\\ p_{31} & p_{32} & p_{33} & \dots\\ \dots & \dots & \dots & \dots\\ \end{bmatrix} \qquad

考虑马尔科夫链X=\{X_0,X_1,\dots,X_t,\dots\}在时刻t的概率分布,称为时刻t的状态分布:

\pi(t)=\begin{bmatrix} \pi_1(t)\\ \pi_2(t)\\ \dots\\ \dots\\ \end{bmatrix}

其中\pi_i(t)表示时刻t状态为i的概率P(X_t=i)

马尔科夫链X在时刻t的状态分布可以由在时刻(t-1)的状态分布以及转移概率分布决定:

\pi(t)=P\pi(t-1)

递推得到:

\pi(t)=P^t \pi(0)

若状态空间S上存在一个分布\pi使得:

\pi=P\pi

则称\pi为马尔科夫链X=\{X_0,X_1,\dots,X_t,\dots\}平稳分布

下面给出几个定义及定理:

(1)不可约:对于任意状态i,j\in S,存在时刻t,使得时刻0从状态j出发,时刻t到达状态i的概率大于0,则称此马尔科夫链X是不可约的。否则称马尔科夫链是可约的。

(2)非周期:对任意状态i\in S,若时刻0从状态i出发,t时刻返回状态的所有时间长\{ t:P(X_t=i|X_0=i)>0 \}的最大公约数为1,则称此马尔科夫链X是非周期的,否则称马尔科夫链是周期的。

(3)正常返:对于任意状态i,j\in S,定义概率p_{ij}^t为时刻0从状态j出发,时刻t首次转移到状态i的概率,若对任意状态i,j\in S都满足:\lim_{t\rightarrow \infty}p_{ij}^t>0,则称马尔科夫链X是正常返的。

(4)可逆:若有一个状态分布\pi=(\pi_1,\pi_2,\dots)^T,对任意状态i,j\in S,对任意时刻t,满足:

p_{ji}\pi_i=p_{ij}\pi_j

则称此马尔科夫链X可逆马尔科夫链。上式称为细致平衡方程

定理1:不可约且非周期的有限状态马尔科夫链有唯一平稳分布存在。

定理2:不可约、非周期且正常返的马尔科夫链有唯一平稳分布存在。

定理3(遍历定理):不可约、非周期且正常返的马尔科夫链有唯一平稳分布\pi=(\pi_1,\pi_2,\dots)^T,且转移概率的极限分布式马尔科夫链的平稳分布:

\lim_{t\rightarrow \infty}P(X_t=i|X_0=j)=\pi_i,\quad i=1,2,\dots,\quad j=1,2,\dots

f(X)是定义在状态空间上的函数,E_{\pi}[|f(X)|]<\infty,则:

P \{\hat{f}_t\rightarrow E_{\pi}[|f(X)|] \}=1

这里:

\hat{f}_t=\frac{1}{t}\sum_{s=1}^t f(x_s)

遍历定理是说,满足相应条件的马尔科夫链,当时间趋于无穷时,马尔科夫链的状态分布趋近于平稳分布,随机变量的函数的样本均值以概率1收敛于该函数的数学期望。

定理4:满足细致平衡方程的状态分布\pi就是该马尔科夫链的平稳分布,即:

P\pi=\pi

证明:

(P\pi)_i=\sum_j p_{ij}\pi_j=\sum_j p_{ji}\pi_i=\pi_i\sum_j p_{ji}=\pi_i

3、马尔科夫链蒙特卡罗法

假设多元随机变量x概率密度函数为p(x)f(x)x的函数,目标是获得概率分布p(x)的样本集合,以及函数f(x)的数学期望E_{p(x)}[f(x)]

应用马尔科夫链蒙特卡罗法解决这个问题的基本想法是:在随机变量x的状态空间S上定义一个满足遍历定理的马尔科夫链X=\{X_0,X_1,\dots,X_t,\dots \},使其平稳分布就是抽样的目标分布p(x)。然后在这个马尔科夫链上进行随机游走,每个时刻得到一个样本。根据遍历定理,当时间趋于无穷时,样本的分布趋于平稳分布,样本的函数均值趋近于函数的数学期望。

所以,当时间足够长时(时刻大于某个正整数m0\sim m的时间段称为燃烧期),在之后的时间(时刻小于某个给定正整数nn>m)里随机游走得到的样本集合\{x_{m+1},x_{m+2},\dots,x_n\}就是目标概率分布的抽样结果,得到的函数均值就是要计算的期望值:

\hat{E}f=\frac{1}{n-m}\sum_{i=m+1}^n f(x_i)

马尔科夫链蒙特卡罗法可概括为以下三步:

(1)在随机变量x的状态空间S上定义一个满足遍历定理的马尔科夫链X=\{X_0,X_1,\dots,X_t,\dots \},使其平稳分布就是抽样的目标分布p(x)

(2)从状态空间的某一点x_0出发,用构造的马尔科夫链进行随机游走,产生样本序列x_0,x_1,\dots,x_t,\dots

(3)应用马尔科夫链的遍历定理,确定整数mn,得到样本集合\{x_{m+1},x_{m+2},\dots,x_n\},求得函数f(x)的均值:

\hat{E}f=\frac{1}{n-m}\sum_{i=m+1}^n f(x_i)

相关文章

网友评论

      本文标题:马尔科夫链蒙特卡罗法

      本文链接:https://www.haomeiwen.com/subject/kdveqctx.html