马尔科夫链蒙特卡洛法学习笔记

马尔科夫链蒙特卡洛法

蒙特卡洛法

抽样方法

蒙特卡洛法的抽样方法本质就是一个间接抽样方法,比如目标分布p(x)不可以直接抽样,但是q(x)容易抽样,那么可以设定$a(x) = \frac{p(x)}{q(x)}$为拒绝分布,堆q(x)的分布进行抽样,然后再以a(x)的概率接受这次抽样,1-a(x)的概率拒绝抽样,这么最后的结果就是分布p(x).

也叫接受-拒绝抽样法

数学期望估计

$$E_{p(x)}[f(x)] =\underset{n \rightarrow \infty }{\lim} \frac{1}{n}\displaystyle \sum_{i=1}^n f(x_i)$$

p(x)是概率密度函数,f(x)是定义在X上的函数,抽取n个样本$x_1,…,x_n$,计算出样本均值$\bar f_n =\frac{1}{n}\displaystyle \sum_{i=1}^n f(x_i)$

积分计算

用样本均值的计算代替积分运算.如果$h(x) = f(x)p(x)$,那么$\int_xh(x) d(x) = \int_xf(x)p(x) d(x) = E_{p(x)}[f(x)]$

马尔科夫链

基本概念

满足下列性质的随机序列称为n阶马尔科夫链.
$$P(X_t|X_0,…,X_{t-1}) = P(X_t|X_{t-1}…X_{t-n})$$
其中n=1称为一阶马尔可夫

离散状态马尔科夫链

离散状态马尔科夫链$X={X_0,…X_n}$,随机变量$X_t$定义在离散空间S中,转移概率分布可以由矩阵表示.
$$P = (p_{ij}) = P(X_t=i|X_{t-1} = j)$$
且$\displaystyle \sum p_{ij} = 1$
考虑马尔科夫链在时刻t的概率分布称为时刻t的状态分布,记作
$$\pi(t) = \begin{bmatrix}
\pi_1(t)\
\vdots\
\pi_n(t)
\end{bmatrix}$$
其中$\pi_i(t) = P(X_t = i)$
如果在状态空间S上存在这样一个分布:
$$\pi = P\pi$$
则称$\pi$为马尔科夫链X的平稳分布,其意义是以平稳分布作为起始分布,那么之后任何一个时刻都是该平稳分布.

连续状态的马尔科夫链

连续状态马尔科夫链$X = {X_0,…,X_n}$随机变量$X_t$定义在连续空间S,转移概率分布由概率转移或转移核表示.
设S是连续空间,对任意$x \in S, A \subset S$,转移核P(x,A)表示$P(X_t \subset A | X_{t-1} = x )$定义为
$$P(x,A) = \int_A p(x,y)dy$$
同理,其也有平稳分布.
$$\pi(y) = \int p(x,y)\pi(x)dx$$

马尔可夫链的性质

不可约

对于任意$i,j \in S$,存在t, $s.tP(X_t = i | X_{t-1} = j) \gt 0$

非周期

对任意$i \in S$,如果时刻0从i出发,t时刻返回i状态的所有市场的最大公约数为1,则称其为非周期的,否则称其为周期的.

正常返

对于任意$i,j \in S$,,定义概率$p_{ij}^t = P(X_t = i,X_s \ne i,s = 1,…,t-1|X_0 = j)$,若对所有状态i,j满足$\underset{t \rightarrow \infty}{lim}p_{ij}^t \gt 0$,则称马尔科夫链X是正常返的

遍历性质

如果马尔可夫链X是不可约的,非周期的,正常返的,则该马尔科夫链由唯一平稳分布$\pi$,并且转移概率的极限分布是马尔科夫链的平稳分布.
$$\underset{t \rightarrow \infty }{lim}P(X_t = i |X_0 =j) = \pi_i$$
####　可逆马尔可夫链
如果对于马尔科夫链X的任意一个时刻皆有
$$P(X_t |X_{t-1} =j)\pi_j = P(X_{t-1} = j|X_t = i)\pi_i(细致平衡方程)$$
则称起为可逆的.

马尔科夫链蒙特卡洛法

基本思想

假设多元随机变量x,满足$x \in X$,其概率分布为p(x),f(x)是定义在$x \in X$上的函数,求概率分布p(x)的样本集合,以及函数f(x)的数学期望.
马尔科夫链蒙特卡洛法的基本思路是:

• 定义一个满足遍历定理的马尔科夫链,使其平稳分布为抽样的目标分布
• 在该马尔科夫链上随机游走,由遍历性质可得在足够长时间之后,样本的分布趋近平稳分布,样本的函数均值接近函数的数学期望.

算法的难点是找到这个马尔科夫链.

这个算法满足遍历定理,所以起始位置对结果没有影响.
对于最后的收敛,可以多次游走之后判断收敛值是否相同,或者多个点同时多次游走,判断是否相同.

Metropois-Hastings算法

基本原理

该原理类似接受拒绝抽样法.
假设要抽样的概率分布为p(x),该算法采用转移核为$p(x,x^\prime)$的马尔科夫链:
$$p(x,x^\prime) = q(x,x^\prime)\alpha(x,x^\prime)$$
其中$q(x,x^\prime)$和$\alpha(x,x^\prime)$分别称为建议分布和接受分布.

建议分布$q(x,x^\prime)$是另一个马尔科夫链的转移核,并且其是不可约的,即其概率值恒不为0,同时,这个分布是一个容易抽样的分布.接受分布则是
$$\alpha(x,x^\prime) = \text min(1,\frac{p(x^\prime)q(x^\prime,x)}{p(x)q(x,x^\prime)})$$
若$p(x^\prime)q(x^\prime,x) \gt p(x)q(x,x^\prime)$,则
$$ \begin{aligned} p(x)p(x,x^\prime) &= p(x)q(x,x^\prime)= q(x,x^\prime)p(x^\prime)\frac{p(x)q(x^\prime,x)}{p(x^\prime)q(x^\prime,x)} \
&= p(x^\prime)q(x^\prime,x)\alpha(x^\prime,x)\
&= p(x^\prime)p(x^\prime,x)
\end{aligned}$$

若$p(x^\prime)q(x^\prime,x) \lt p(x)q(x,x^\prime)$,则
$$
\begin{aligned}
p(x)p(x,x^\prime) &= p(x)q(x,x^\prime)\frac{p(x^\prime)q(x^\prime,x)}{p(x)q(x,x^\prime)}\
&=p(x^\prime)q(x^\prime,x)\
&=p(x^\prime)q(x^\prime,x)\alpha((x^\prime,x)\
&=p(x^\prime)p(x^\prime,x)
\end{aligned}
$$
综上,$p(x,x^\prime)$是可逆的.

建议分布

建议分布的常用形式:

• 第一种形式:建议分布是对称的,有$q(x,x^\prime) = q(x^\prime,x)$.
• • 特例之一是取条件概率分布$q(x,x^\prime) = p(x^\prime|x)$
• • 特例之二是令$q(x,x^\prime) = q(|x-x^\prime|)$,称为随机游走的Metropolis算法
• 第二种形式:独立抽样,假设$q(x,x^\prime)$与当前状态无关,建议分布按照$q(x^\prime)$独立抽样进行.

满条件分布

对于多元联合概率分布$p(x) = p(x_1,…,x_k)$,其中x为k维随机变量.如果条件概率分布$p(x_I|x_{-I})$中所有k个变量全部出现,其中$x_I = {x_i,i\in I},x_{-I} = {x_i,i\not\in I},I \subset K = {1,…,k}$,那么称这总条件概率分布维满条件分布.

它有一下性质:
$$p(x_I|x_{-I}) = \frac{p(x)}{\int p(x)dx_i} \propto p(x)$$
另有
$$\frac{p(x_I^\prime|x_{-I}^\prime)}{p(x_I|x_{-I})} = \frac{p(x^\prime)}{p(x)}$$
其意义是通过计算满条件概率的bi来计算联合概率的比,而前者更容易计算.

算法

输入:抽样的目标分布的密度函数p(x),函数f(x)
输出:p(x)的随机抽样样本$x_{m+1},…,x_ {n}$,函数样本均值$f_{mn}$
参数:收敛步数m,迭代步数n.

1. 任意选择一个$x_0$作为起始状态.
2. 对i=1,…,n进行循环执行

• 设状态$x_{i-1} = x$,按照建议分布$q(x,x^\prime)$随机抽样一个候选状态$x^\prime$
• 计算接受概率$\alpha(x,x^\prime)$
• 从区间[0,1]中按照均匀分布随机抽取一个数,若$u \leq \alpha(x,x^\prime)$则接受,反之则拒绝

3. 得到样本集合$x_{m+1},…,x_{n}$

单分量的M-H算法

对于多元变量分布的抽样是困难的,所以对多元变量的每一变量的条件分布一次分别进行抽样,从而实现对整个多元变量的一次抽样,这就是但分量M-H算法.

假设马尔科夫链的状态有k维随机变量表示
$$x = (x_1,…,x_k)^T$$
其中$x_j$表示变量x的第j个分量,$x^i$表示在时刻i的状态.

为了生成容量为n的样本集合${x^1,…,x^n}$,单分量M-H算法由下面的k步迭代实现一次迭代.

设在第(i-1)次迭代结束时分量$x_j$的取值时$x_j^{i-1}$,在第i次迭代的第j步,对分量$x_j$更加M-H算法更新,得到其新的取值$x_j^i$.

• 由建议分布$q(x_j^{i-1},x_j^i|x_{-j}^i)$抽样产生分量的候选值$x_j^i$,这里的$x_{-j}^i$表示第i次迭代的第(j-1)步后的$x_i$出去$x_j^{i-1}$的所有值,即
  $$x_{-j}^i = (x_1^i,…,x_{j-1}^i,x_{j+1}^{i-1},…,x_k^{i-1})$$
• 根据接受概率$\alpha(x_j^{i-1},x_j^{i}|x_{-j}^i) = \text min(1,\frac{p(x_j^i|x_{-j}^i)q(x_j^i,x_j^{i-1}|x_{-j}^i)}{p(x_j^{i-1}|x_{-j}^i)q(x_j^{i-1},x_j^{i}|x_{-j}^i)})$ 抽样决定是否接受