单分量MH算法,简单来说就是对高维随机变量遍历其每一个维度,并对每个维度单独进行采样,这样就将一个复杂的采样问题分拆成了一个个更加简单的采样问题,聪明童鞋可能已经想到了,要怎么操作才能保证细致平衡方程成立勒,这就需要用到满条件分布咯,它将会为我们的单个维度采样打开大门,下面先说明一下
MCMC的目标分布通常是多元的联合概率分布$p(x)=p(x_1,x_2,...,x_k)$,其中$x=(x_1,x_2,...,x_k)^T$为$k$维随机变量。如果条件概率分布$p(x_I\mid x_{-I})$中所有$k$个变量全部出现,其中$x_I=\{x_i,i\in I\},x_{-I}=\{x_i,i\notin I\},I\subset K=\{1,2,...,k\}$,那么称这种条件概率分布为满条件分布。
满条件分布有一个非常有用的性质,那就是对于随机变量的任意取值$x,x'$,和任意的$I\subset K$,以下等式成立:
$$ \frac{p(x'_I\mid x'_{-I})}{p(x_I\mid x_{-I})}=\frac{p(x')}{p(x)} $$惊不惊喜,意不意外,有了上面的等式那就意味着对$\frac{p(x')}{p(x)}$的计算,只需落到对单个维度上的条件概率进行计算即可!这在方便联合概率之比的计算的同时也为单分量MH算法提供了理论支撑,下面先简单推导一下上面的等式为什么成立:
$$ p(x_I\mid x_{-I})=\frac{p(x)}{p(x_{-I})}=\frac{p(x)}{\int_{x_I}p(x)dx_I}\propto p(x) $$在写具体算法之前,先对相应的符号表达做个解释:
(1)设在第$i$次迭代结束时分量$x_j$的取值为$x_j^{(i)}$,其中$j=1,2,...,k,i=0,1,2,...$;
(2)$x_{-j}^{(i)}$表示在第$i$轮,更新到第$j-1$维时,除去第$j$维后的所有值,即前$j-1$维已经更新为第$i$轮的值,而$j+1$维直至最后维都还是$i-1$轮的值,即:
$$ x_{-j}^{(i)}=(x_1^{i},...x_{j-1}^i,x_{j+1}^{(i-1)},...x_k^{(i-1)})^T $$好了,可以描述算法流程了
输入:目标概率分布的密度函数$p(x)$,正整数$m,n,m<n$;
输出:$p(x)$的随机样本$x_{m+1},x_{m+2},...,x_n$
(1)初始化样本$x^{(0)}=(x_1^{(0)},x_2^{(0)},...,x_k^{(0)})$
(2)对$i=1,2,...,m,m+1,...n$:
对$j=1,2,...,k$:
(2.1)根据建议分布$q(x_j^{(i-1)}\rightarrow x_j^{'(i)}\mid x_{-j}^{(i)})$抽样第$j$维的候选值$x_j^{'(i)}$
$$ \alpha(x_j^{(i-1)}\rightarrow x_j^{'(i)}\mid x_{-j}^{(i)})=min\{1,\frac{p(x_j^{'(i)}\mid x_{-j}^{(i)})q(x_j^{'(i)}\rightarrow x_j^{(i-1)}\mid x_{-j}^{(i)})}{p(x_j^{(i-1)}\mid x_{-j}^{(i)})q(x_j^{(i-1)}\rightarrow x_j^{'(i)}\mid x_{-j}^{(i)})}\} $$(2.2)计算第$j$维的接收率:
(2.3)从$(0,1)$均匀分布中随机采样$u$,如果$u<\alpha(x_j^{(i-1)}\rightarrow x_j^{'(i)}\mid x_{-j}^{(i)})$,则令$x_j^{(i)}=x_j^{'(i)}$,否则$x_j^{(i)}=x_j^{(i-1)}$
(3)返回样本集$\{x_{m+1},x_{m+2},...,x_n\}$
完事儿了,要证明单分量MH算法满足细致平衡方程也很简单,下面简单说明一下:
(1)首先条件概率之比$\frac{p(x_j^{'(i)}\mid x_{-j}^{(i)})}{p(x_j^{(i-1)}\mid x_{-j}^{(i)})}$根据满条件分布中的等式可以替换为联合概率之比;
(2)其次,任意建议分布$q(x_j^{(i-1)}\rightarrow x_j^{'(i)}\mid x_{-j}^{(i)})$均不会影响到细致平衡方程(见前一个note的证明);
综合这两条我们可以推导出(下面的$\sim$表示分布相同):
$$ p(x_1^{(i-1)},x_2^{(i-1)},...,x_k^{(i-1)})\\ \sim p(x_1^{(i)},x_2^{(i-1)},...,x_k^{(i-1)})\\ \cdots\\ \sim p(x_1^{(i)},x_2^{(i)},...x_j^{(i)},x_{j+1}^{(i-1)},...,x_k^{(i-1)})\\ \cdots\\ \sim p(x_1^{(i)},x_2^{(i)},...,x_k^{(i)})\\ $$所以:$p(x_1^{(i-1)},x_2^{(i-1)},...,x_k^{(i-1)})\sim p(x_1^{(i)},x_2^{(i)},...,x_k^{(i)}),i=1,2,...$
由于单分量MH是用于高维空间的采样,所以这里我造一个二元高斯分布$N(u,\Sigma)$,其对应的均值、协方差如下:
$$ u=[0,0]^T\\ \Sigma=\begin{bmatrix} 1 &\rho \\ \rho & 1 \end{bmatrix} $$那么,它对应的条件概率分布为一元高斯分布,如下:
$$ p(x_1\mid x_2)=N(\rho x_2,(1-\rho^2))\\ p(x_2\mid x_1)=N(\rho x_1,(1-\rho^2)) $$我们不妨取$\rho=0.5$
import os
os.chdir('../')
from ml_models import utils
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
#定义均值,协方差
u=np.asarray([0,0])
sigma=np.asarray([[1,0.5],
[0.5,1]])
#看看目标分布的效果
X=np.random.multivariate_normal(mean=u, cov=sigma, size=200)
utils.plot_contourf(X,lambda x:utils.gaussian_nd(x,u,sigma),8)
D:\ProgramData\Anaconda3\lib\site-packages\matplotlib\contour.py:960: UserWarning: The following kwargs were not used by contour: 'linewidth' s)
建议分布使用以当前维度的取值为均值,方差1为一元高斯分布,这样它们在该维的建议分布取值对称:
$$ q(x_j^{'(i)}\rightarrow x_j^{(i-1)}\mid x_{-j}^{(i)})=q(x_j^{(i-1)}\rightarrow x_j^{'(i)} \mid x_{-j}^{(i)}) $$对于计算接受概率就更方便了:
$$ \alpha(x_j^{(i-1)}\rightarrow x_j^{'(i)}\mid x_{-j}^{(i)})=min\{1,\frac{p(x_j^{'(i)}\mid x_{-j}^{(i)})}{p(x_j^{(i-1)}\mid x_{-j}^{(i)})}\} $$import copy
#采样的样本量
nums=1000
count=0
points=[]
#采样x0
point=[np.random.randn(),np.random.randn()]
points.append(point)
while count<nums:
for k in (0,1):
new_point=copy.deepcopy(point)
#按照q(x,x')
new_point[k]=np.random.randn()+point[k]
#alpha(x,x')
alpha=min(1.,utils.gaussian_1d(new_point[k],0.5*new_point[(k+1)%2],0.75)/utils.gaussian_1d(point[k],0.5*point[(k+1)%2],0.75))
#从(0,1)均匀采样一个u
u=np.random.random()
#判断是否接收新点还是旧点
if u<alpha:
point=new_point
points.append(point)
count+=1
utils.plot_contourf(np.asarray(points[-400:]),lambda x:utils.gaussian_nd(x,u,sigma),8)
D:\ProgramData\Anaconda3\lib\site-packages\matplotlib\contour.py:960: UserWarning: The following kwargs were not used by contour: 'linewidth' s)
效果看起来还不错,然后我们继续想一种情况,如果目标分布$p(x)$的采样困难是由于随机变量维度太高造成的,而它在单个维度(满条件概率分布)上采样是相对比较容易的,那么我们可以用单维度上的满条件概率分布替换建议分布来采样,即:
$$ q(x_j^{(i-1)}\rightarrow x_j^{'(i)} \mid x_{-j}^{(i)})=p(x_j^{'(i)} \mid x_{-j}^{(i)}) $$这种情况下又会发生什么神奇的事情呢?请看下一节Gibbs采样