这篇博客是在自学李航老师的《统计学习方法》19章MCMC的过程中,进行简要的复述和记录,之后会在二刷三刷博客的时候增加自己对方法对口语化描述及理解,关注我一起加油吧~~
文章目录
蒙特卡罗法(也称为统计模拟方法)是通过从概率模型的随机抽样进行近似数据计算的方法。 马尔可夫链蒙特卡罗(MCMC)法是以马尔可夫链为概率模型的蒙特卡罗法。
MCMC方法的基本思想是:通过蒙特卡罗法构建一个马尔可夫链,使其平稳分布就是要进行抽样的分布,首先基于该马尔可夫链进行随机游走,产生样本序列,之后使用该平稳分布的样本进行近似的数值计算。
19.1 蒙特卡罗法
19.1.1 随机抽样
统计学习和机器学习的目的是基于数据对概率分布的特征进行判断。蒙特卡罗法要解决的问题是,假设概率分布的定义已知,通过抽样获得概率分布的随机样本,通过得到的样本对概率分布的特征进行分析,蒙特卡罗法的核心是随机抽样。
一般的蒙特卡罗法有直接抽样法、接受-拒绝抽样法、重要性抽样法等。后两种方法适合于概率密度函数复杂(如密度函数有多个变量,各变量相互不独立,密度函数形式复杂),不能直接抽样的方法。
- 接受-拒绝抽样法
假设有随机变量 x x x,取值 x ∈ χ x\in\chi x∈χ,其概率密度函数为 p ( x ) p(x) p(x).目标是得到该概率分布的随机样本,进而对这个概率分布进行分析。基本思想如下:假设 p ( x ) p(x) p(x)不可以直接抽样,找一个可以直接抽样的分布,称为建议分布。假设 q ( x ) q(x) q(x)是建议分布的概率密度函数,并且有 c q ( x ) ≥ p ( x ) , 且 c ≥ 0 cq(x)\geq p(x),且c\geq0 cq(x)≥p(x),且c≥0。对 q ( x ) q(x) q(x)进行抽样,假设得到的结果是 x ∗ x^* x∗,再按照 p ( x ) c q ( x ) \frac {p(x)}{cq(x)} cq(x)p(x)的例随机决定是否接受 x ∗ x^* x∗。接受拒绝法实际就是按照 p ( x ) p(x) p(x)的涵盖面积(涵盖体积)占 c q ( x ) cq(x) cq(x)的涵盖面积的比例进行抽样。 - 接受-拒绝法
输入:抽样的目标概率分布的概率密度函数 p ( x ) p(x) p(x);
输出:概率分布的随机样本 x 1 , x 2 , ⋅ ⋅ ⋅ , x n x_1,x_2,\cdot\cdot\cdot,x_n x1,x2,⋅⋅⋅,xn.
参数:样本 n n n.
(1)选择概率密度函数为 q ( x ) q(x) q(x)的概率分布作为建议分布,使其对任一 x x x满足 c q ( x ) ≥ p ( x ) cq(x)\geq p(x) cq(x)≥p(x),其中 c ≥ 0 c\geq0 c≥0.
(2)按照建议分布 q ( x ) q(x) q(x)随机抽样得到样本 x ∗ x^* x∗,再按照均匀分布在 ( 0 , 1 ) (0,1) (0,1)范围内进行抽样得到 u u u.
(3)如果 u ≤ p ( x ∗ ) c q ( x ∗ ) u\leq\frac {p(x^*)}{cq(x^*)} u≤cq(x∗)p(x∗),则将 x ∗ x^* x∗作为抽样结果;否则,回到步骤(2).
(4)直至得到 n n n个随机样本,结束。 - 接受-拒绝法的优点是容易实现,缺点是效率可能不高。(如: p ( x ) p(x) p(x)的涵盖体积占 c q ( x ) cq(x) cq(x)的涵盖体积的比例较低,就会导致拒绝的比例很高,抽样效率很低。)
19.2 数学期望估计
一般的蒙特卡罗法,如直接抽样法、接受-拒绝法、重要性抽样法等,也可以直接用于数学期望估计。假设随机变量
x
x
x,取值
x
∈
χ
x\in\chi
x∈χ,其概率密度函数为
p
(
x
)
p(x)
p(x),
f
(
x
)
f(x)
f(x)为定义在
χ
\chi
χ上的函数,目标是求函数
f
(
x
)
f(x)
f(x)关于密度函数
p
(
x
)
p(x)
p(x)的数学期望
E
p
(
x
)
[
f
(
x
)
]
E_{p(x)}[f(x)]
Ep(x)[f(x)].
蒙特卡罗法按照概率分布
p
(
x
)
p(x)
p(x)独立地选取
n
n
n个样本
x
1
,
x
2
,
⋅
⋅
⋅
,
x
n
x_1,x_2,\cdot\cdot\cdot,x_n
x1,x2,⋅⋅⋅,xn,比如用以上的抽样方法,计算函数
f
(
x
)
f(x)
f(x)的样本均值
f
^
n
\hat f_n
f^n,
f
^
n
=
1
n
∑
i
=
1
n
f
(
x
i
)
\hat f_n=\frac {1}{n}\sum_{i=1}^nf(x_i)
f^n=n1i=1∑nf(xi)作为数学期望
E
p
(
x
)
[
f
(
x
)
]
E_{p(x)}[f(x)]
Ep(x)[f(x)]的近似值。根据大数定律可知,当样本容量增大时,样本均以概率1收敛于数学期望,因此得到了数学期望的近似计算方法:
E
p
(
x
)
[
f
(
x
)
]
≈
1
x
=
n
∑
i
=
1
n
f
(
x
i
)
E_{p(x)}[f(x)]\approx \frac{1}{x=n}\sum_{i=1}^nf(x_i)
Ep(x)[f(x)]≈x=n1i=1∑nf(xi)
19.3 积分计算
一般的蒙特卡罗法也可以用于定积分的近似计算,称为蒙特卡罗积分。假设有一个函数
h
(
x
)
h(x)
h(x),目标是计算该函数的积分
∫
χ
h
(
x
)
d
x
\int_\chi h(x)dx
∫χh(x)dx
如果能够将函数
h
(
x
)
h(x)
h(x)分解成一个函数
f
(
x
)
f(x)
f(x)和一个概率密度函数
p
(
x
)
p(x)
p(x)的乘积的形式,那么就有
∫
χ
h
(
x
)
d
x
=
∫
χ
f
(
x
)
p
(
x
)
d
x
=
E
p
(
x
)
[
f
(
x
)
]
\int_\chi h(x)dx=\int_\chi f(x)p(x)dx=E_{p(x)}[f(x)]
∫χh(x)dx=∫χf(x)p(x)dx=Ep(x)[f(x)]
于是函数
h
(
x
)
h(x)
h(x)的积分可以表示为函数
f
(
x
)
f(x)
f(x)关于概率密度函数
p
(
x
)
p(x)
p(x)的数学期望。任何一个函数的积分都可以表示为某一个函数的数学期望的形式:给定一个概率密度函数
p
(
x
)
p(x)
p(x),只要取
f
(
x
)
=
h
(
x
)
p
(
x
)
f(x)=\frac{h(x)}{p(x)}
f(x)=p(x)h(x) .函数的数学期望又可以通过函数的样本均值估计,于是就可以采用样本均值来近似计算积分。这就是蒙特卡罗积分的基本思想。
∫
χ
h
(
x
)
d
x
=
E
p
(
x
)
[
f
(
x
)
]
≈
1
n
∑
i
=
1
n
f
(
x
i
)
\int_\chi h(x)dx=E_{p(x)}[f(x)]\approx \frac {1}{n}\sum_{i=1}^nf(x_i)
∫χh(x)dx=Ep(x)[f(x)]≈n1i=1∑nf(xi)
马尔可夫链蒙特卡罗法也适用于概率密度函数复杂,不能直接抽样的情况,旨在解决一般的蒙特卡罗法抽样效率不高的问题。一般的蒙特卡罗法中的抽样样本是独立的,而马尔可夫链蒙特卡罗法中的抽样样本不是独立的,样本序列形成马尔可夫链。
19.2 马尔可夫链
19.2.1 基本定义
19.2.3 连续状态马尔可夫链
连续状态的马尔可夫链
X
=
X
0
,
X
1
,
⋅
⋅
⋅
,
X
t
X={X_0,X_1,\cdot\cdot\cdot,X_t}
X=X0,X1,⋅⋅⋅,Xt,随机变量
X
t
(
t
=
0
,
1
,
2
,
⋅
⋅
⋅
)
X_t(t=0,1,2,\cdot\cdot\cdot)
Xt(t=0,1,2,⋅⋅⋅)定义在连续状态空间
S
S
S,转移概率分布由概率转移核或转移核表示。
设
S
S
S是连续状态空间,对任意的
x
∈
S
,
A
∈
S
x \in S,A \in S
x∈S,A∈S,转移核
P
(
x
,
A
)
P(x,A)
P(x,A)定义为
P
(
x
,
A
)
=
∫
A
p
(
x
,
y
)
d
y
P(x,A)=\int_Ap(x,y)dy
P(x,A)=∫Ap(x,y)dy 其中
p
(
x
,
⋅
)
p(x,\cdot)
p(x,⋅)表示概率密度函数,满足
p
(
x
,
⋅
)
≥
0
,
P
(
x
,
S
)
=
∫
S
p
(
x
,
y
)
d
y
=
1
p(x,\cdot)\geq 0,P(x,S)=\int_Sp(x,y)dy=1
p(x,⋅)≥0,P(x,S)=∫Sp(x,y)dy=1。转移核
P
(
x
,
A
)
P(x,A)
P(x,A)表示从
x
~
A
x~A
x~A的转移概率。
P
(
X
t
=
A
∣
X
t
−
1
=
x
)
=
P
(
x
,
A
)
P(X_t=A|X_{t-1}=x)=P(x,A)
P(Xt=A∣Xt−1=x)=P(x,A)有时也将概率密度函数
p
(
x
,
⋅
)
p(x,\cdot)
p(x,⋅)称为转移核。
若马尔可夫链的状态空间
S
S
S上的概率分布
π
(
x
)
\pi(x)
π(x)满足条件
π
(
y
)
=
∫
p
(
x
,
y
)
π
(
x
)
d
x
,
∀
y
∈
S
\pi(y)=\int p(x,y)\pi(x)dx,\forall y\in S
π(y)=∫p(x,y)π(x)dx,∀y∈S则称分布
π
(
x
)
\pi(x)
π(x)为该马尔可夫链的平稳分布。等价地,
π
(
A
)
=
∫
P
(
x
,
A
)
π
(
x
)
d
x
,
∀
A
⊂
S
\pi(A)=\int P(x,A)\pi(x)dx,\forall A\subset S
π(A)=∫P(x,A)π(x)dx,∀A⊂S或简写为
π
=
P
π
\pi=P\pi
π=Pπ
19.2.4 马尔可夫链的性质
-
不可约
-
非周期
-
正常返
-
遍历定理
-
可逆马尔可夫链
19.3 马尔可夫链蒙特卡罗法
19.4 Metropolis-Hastings算法
19.5 吉布斯抽样(Gibbs Sampling)
吉布斯抽样用于多元变量联合分布的抽样和估计。其基本做法是:从联合概率分布定义满条件概率分布,依次对满条件概率分布进行抽样,得到样本的序列。抽样的过程就是在马尔可夫链上的随机游走,每一个样本对应着马尔可夫链的状态,平稳分布目标的联合分布。