MCMC
蒙特卡罗方法马尔可夫链
蒙特卡罗方法
方法引入
若f(x)非均匀分布,假设x的概率密度函数为p(x),则
∫baf(x)dx=∫baf(x)p(x)p(x)dx≈1nn−1∑i=0f(xi)p(xi)
也就是说先根据x的分布函数p(x)来采样x,获取n个数值f(xi)p(xi),再取n个值的平均值即可得到积分的近似值。
若\(f(x)\)均匀分布,则
\[\begin{equation*}
\int_a^bf(x)\,dx
\end{equation*} \approx
\frac{b-a}{n}\sum\limits_{i=0}^{n-1}f(x_i)
\]
接受拒绝采样
采样过程
利用常见的分布函数如高斯分布\(g(x)\)和常数\(c\),使得任意的 \(x\)都有\(f(x)<=cg(x)\)
在\(g(x)\)中采样获取 \(x\) 变量的一个值 \(z\)
- 在均匀分布(0, 1)中采样一个值 \(u\),若\(u<=\frac{f(z)}{cg(z)}\),则接受这次采样,否则重新采样。
注意:
- \(c=sup_x \{\frac{f(x)}{g(x)}\}\), \(c\)等于\(\frac{f(x)}{g(x)}\)的最小上界
- 采样方法有效性证明见:
缺点:
- 对于一些二维分布𝑝(𝑥,𝑦),有时候我们只能得到条件分布𝑝(𝑥|𝑦)和𝑝(𝑦|𝑥), 却很难得到二维分布𝑝(𝑥,𝑦)一般形式,这时我们无法用接受-拒绝采样得到其样本集
- 对于一些高维的复杂非常见分布𝑝(𝑥1,𝑥2,...,𝑥𝑛),我们要找到一个合适的g(𝑥)和c非常困难
马尔可夫链
定义:
某一时刻状态转移的概率只依赖于它的前一个状态。
假如状态序列为 \(x_1 , \,x_2, \dots , \, x_t, \, x_{t+1} \),则\(p(x_{t+1}|x_1, x_2, \dots, x_t)=p(x_{t+1}|x_t)\)
状态转移矩阵:
假设现在共有三种状态:牛市(Bull market), 熊市(Bear market)和横盘(Stagnant market)。每一个状态都以一定的概率转化到下一个状态。如果我们定义矩阵阵𝑃某一位置𝑃(𝑖,𝑗)的值为𝑃(𝑗|𝑖),即从状态i转化到状态j的概率,并定义牛市为状态0, 熊市为状态1, 横盘为状态2. 牛市转为3种状态的概率为[0.9, 0.075, 0.025],熊市转为3种状态的概率为[0.15, 0.8, 0.05], 横盘转为3种状态的概率为[0.25, 0.25, 0.25]这样我们得到了马尔科夫链模型的状态转移矩阵为:
\[
\begin{bmatrix}
0.9 & 0.075 & 0.025\\
0.15 & 0.8 & 0.05\\
0.25 & 0.25 & 0.5
\end{bmatrix}
\]
状态转移矩阵性质:
如果一个非周期的马尔科夫链有状态转移矩阵𝑃, 并且它的任何两个状态是连通的,那么 \( \lim_ {𝑛→∞}𝑃_{𝑖𝑗}^n\) 与\(i\)无关,我们有:
- \[\lim\limits_{𝑛→∞}P_{ij}^n=\pi(j)\]
- \[
\lim\limits_{𝑛→∞}P^n=
\begin{bmatrix}
\pi(1) & \pi(2) & \dots \pi(j) &\dots\\
\pi(1) & \pi(2) & \dots \pi(j) &\dots\\
\dots & \dots & \dots & \dots\\
\pi(1) & \pi(2) & \dots \pi(j) &\dots\\
\dots & \dots & \dots & \dots
\end{bmatrix}
\] - \[\pi(j)=\sum\limits_{i=0}^∞\pi(i)𝑃_{𝑖𝑗}\]
- \(\pi\) 是 \(\pi P=\pi\)的唯一非负解:
\[\pi=[\pi(1), \,\pi(2), \, \dots, \, \pi(j), \, \dots], \sum\limits_{i=0}^∞\pi(i)=1\]
马尔可夫链采样过程:
- 已知稳态分布的状态转移矩阵 \(P\), 状态转移至稳态的次数阈值 \(n_1\), 抽样数据集大小 \(n_2\)
- 从简单概率分布函数(通常是均匀分布)中采样出初始值 \(x_0\)
- for t in 0 to n1+ n2 - 1: 从 \(P(x_{t+1}|x_t)\) 中采样得到样本 \(x_{t+1}\)
- 样本集 \((x_{n1}, x_{n1+1}, \dots, x_{n1+n2-1}\,) \) 就是我们需要的样本集
问题:
如何获取任一稳态分布 \(\pi\) 的状态转移矩阵 \(P\) 呢?
MCMC和M-H采样
马尔可夫链的细致平稳条件:
\[\pi(i)P(i, j)=\pi(j)P(j, i)\]
MCMC采样:
对于任一稳态分布,很难通过细致平稳条件找到对应的状态转移矩阵 \(P\)。但是我们可以通过引入 \(\alpha\) 矩阵来达到条件。
\[\pi(i)Q(i, j)\alpha(i, j)=\pi(j)Q(j, i)\alpha(j, i)\]
其中 \(Q(i, j)\) 表示任一状态转移矩阵,\(\alpha(i, j)=\pi(j)Q(j, i), \alpha(j, i)=\pi(i)Q(i, j)\)
MCMC采样过程:
- 给定任一状态转移矩阵 \(Q\) 和稳态分布 \(\pi\), 状态转移次数阈值 \(n1\), 采样样本大小 \(n2\)
- 从简单概率分布函数(通常是均匀分布)中采样出初始值 \(x_0\)
- for t in 0 to n1+ n2 - 1:
- 从 \(Q(x_{t+1}|x_t)\) 中采样得到样本 \(x_{*}\)
- 在均匀分布(0, 1)中采样一个值 \(u\)
- 若\(u<=\alpha(x_t, x_{*})=\pi(x_*)Q(x_*, x_t)\), 则接受转移 \(x_{t+1}=x_*\),否则不接受 \(x_{t+1}=x_t\)
MCMC采样缺点:
\(\alpha(i, j)\) (又称接收率)通常较小,导致大部分的采样值都被拒绝转移,采样效率很低。
M-H采样:
MCMC的改进版,通过扩大 \(\alpha(i, j)\) 解决采样效率低下的问题。
\[\alpha(i, j) = \min\{\frac{\pi(j)Q(j, i)}{\pi(i)Q(i, j)}, 1\}\]
因为 \(Q\) 一般是对称的,即 \(Q(i, j)=Q(j, i)\),则
\[\alpha(i, j) = \min\{\frac{\pi(j)}{\pi(i)}, 1\}\]
大数据时代M-H采样缺点:
- 数据维度多,计算 \(\alpha(i, j)\)较为复杂,计算效率低下;\(\alpha(i, j)<1\), 采样效率还是低下,能不能不拒绝采样呢?
- 由于特征维度大,很多时候我们甚至很难求出目标的各特征维度联合分布,但是可以方便求出各个特征之间的条件概率分布。这时候我们能不能只有各维度之间条件概率分布的情况下方便的采样呢?
Gibbs采样(吉布斯采样)
多维数据的细致平稳条件:
以二维数据为例:
\[
\begin{aligned}
\pi(x_1^{1}, x_2^{1})\pi(x_2^{2}|x_1^{1})&=\pi(x_1^1)\pi(x_2^1|x_1^1)\pi(x_2^{2}|x_1^{1})\\
\pi(x_1^{1}, x_2^{2})\pi(x_2^{1}|x_1^{1})&=\pi(x_1^1)\pi(x_2^2|x_1^1)\pi(x_2^{1}|x_1^{1})\\
\pi(x_1^{1}, x_2^{1})\pi(x_2^{2}|x_1^{1})&=\pi(x_1^{1}, x_2^{2})\pi(x_2^{1}|x_1^{1})
\end{aligned}
\]
观察上式,我们发现在 \(x_1=x_1^1\) 这条直线上,如果用条件概率分布 \(\pi(x_2|x_1^1)\) 作为马尔科夫链的状态转移概率,则任意两个点之间的转移满足细致平稳条件.
多维数据吉布斯采样过程:
- 输入平稳分布𝜋(𝑥1,𝑥2,...,𝑥𝑛)或者对应的所有特征的条件概率分布,设定状态转移次数阈值𝑛1,需要的样本个数𝑛2
- 随机初始化初始状态值(𝑥(0)1,𝑥(0)2,...,𝑥(0)𝑛
- for 𝑡=0 to 𝑛1+𝑛2−1:
- 从条件概率分布𝑃(𝑥1|𝑥(𝑡)2,𝑥(𝑡)3,...,𝑥(𝑡)𝑛)中采样得到样本𝑥𝑡+11
- 从条件概率分布𝑃(𝑥2|𝑥(𝑡+1)1,𝑥(𝑡)3,𝑥(𝑡)4,...,𝑥(𝑡)𝑛)中采样得到样本𝑥𝑡+12
- ...
- 从条件概率分布𝑃(𝑥𝑗|𝑥(𝑡+1)1,𝑥(𝑡+1)2,...,𝑥(𝑡+1)𝑗−1,𝑥(𝑡)𝑗+1...,𝑥(𝑡)𝑛)中采样得到样本𝑥𝑡+1𝑗
- ...
- 从条件概率分布𝑃(𝑥𝑛|𝑥(𝑡+1)1,𝑥(𝑡+1)2,...,𝑥(𝑡+1)𝑛−1)中采样得到样本𝑥𝑡+1𝑛
样本集 \(\{(𝑥_1^{𝑛1},𝑥_2^{𝑛1},...,𝑥_n^{𝑛1}),...,(𝑥_1^{𝑛1+𝑛2−1},𝑥_2^{𝑛1+𝑛2−1},...,𝑥_n^{𝑛1+𝑛2−1})\}\) 即为我们需要的平稳分布对应的样本集