感情也存在着可以隐瞒的倾斜的程度。

如果有一个随机变量XX,服从概率密度函数为π(x)\pi(x) 的分布。

我们想对这个随机变量抽一组样本出来时,若π(x)\pi(x) 比较简单(如均匀分布,正态分布),则显然比较好模拟。但如果π(x)\pi(x) 比较复杂,不好模拟怎么办。

这时候,需要构造一个马尔可夫过程,即构造一个马尔可夫转移矩阵PPPyxP_{yx} 表示从 x 转移到 y 的概率),使得:

y,π(y)=xPyxπ(x)\forall y,\pi(y)=\sum_xP_{yx}*\pi(x)

矩阵化地,即有πn×1=Pn×nπn×1\pi_{n\times 1}=P_{n\times n}*\pi_{n\times 1}。一个遍历的马尔科夫链有一个很好的性质:对于任意一个分布列pRn×1p\in\mathbb{R}^{n\times 1} 满足p=1,pi0||p||=1,p_i\geq 0,都有limnPnp=π\lim_{n\rightarrow\infty}P^np=\pi。即常说的,“以任意一个状态开始,不断进行马尔可夫转移,最后都会趋于稳定态”。

需要注意的是,马尔科夫链我认为实际上是随机变量间的转移。譬如一系列随机变量为X1,X2,...,XnX_1,X_2,...,X_n,它们对应的样本值是x1,x2,...,xnx_1,x_2,...,x_n

  • X1U[0,1]X_1\sim U[0,1],我们很容易地随机取了个x1x_1 出来。
  • 然后再根据转移矩阵,以概率分布P1x1,P2x1,P3x1,...,Pnx1P_{1x_1},P_{2x_1},P_{3x_1},...,P_{nx_1} 随机抽取得x2x_2,则显然x2x_2 对应的随机变量X2X_2 的概率分布就不是U[0,1]U[0,1] 了。
  • 而重复这个过程,当趋于无穷时,我们就会发现limnXnπ\lim_{n\rightarrow\infty}X_n\sim\pi

那么,其实X1,X2,...,XnX_1,X_2,...,X_n 就服从分布U[0,1],PU[0,1],...,π,π,π,...U[0,1],P*U[0,1],...,\pi,\pi,\pi,...,因此我们考虑足够大的一些样本值,它们都是同分布于π\pi 的,即做到了抽取一组服从分布π\pi 的样本值。


然而如何得到PP 就是一个很大的问题。我发现很多地方都写的PP 是依赖于π\pi 的,那既然都可以模拟PP 了为啥不直接模拟π\pi 呢?

但其实,其简化模拟的核心在于,抽样π\pi 需要知道π\pi 的连续信息,而马尔科夫链的转移只需要知道π\pi 的点信息

考虑 Metropolis-Hastings 算法:

  • QyxQ_{yx} 是一个好模拟的,简单的分布。譬如Qyx=U[0,1]Q_{yx}=U[0,1] 与 x 无关,或Qyx=G(x,0.01)Q_{yx}=G(x,0.01) 即以 x 为均值的正态分布。
  • α(y,x)=π(x)Qxy\alpha(y,x)=\pi(x)*Q_{xy},则可以构造Pyx=Qyxα(y,x)P_{yx}=Q_{yx}*\alpha(y,x)

这样做有什么意义呢?注意到在已知样本x1x_1,想要得到下一个样本x2x_2 时,我们的操作并不是把P1x1,P2x2,...P_{1x_1},P_{2x_2},... 都算出来,然后按照分布抽取。而应该先按照Qyx1Q_{yx_1} 抽取出来一个yy,再计算α(y,x1)\alpha(y,x_1),然后再以α(y,x1)\alpha(y,x_1) 的概率转移到x2=yx_2=y。以1α(y,x)1-\alpha(y,x) 的概率停留在x2=x1x_2=x_1这就是一个简单分布 + 两点分布解决了转移的过程。简单分布负责给出候选值,再计算α\alpha,再根据两点分布决定是否转移。而最终转移的概率就是PyxP_{yx}

于是这样做,显然就只需要知道π\pi 的点信息。

然而这么做为什么是对的?

首先,x2x_2 停留在x1x_1 的概率显然是yQyx(1α(y,x))=1yQyxα(y,x)=1yPyx\sum_yQ_{yx}(1-\alpha(y,x))=1-\sum_{y}Q_{yx}\alpha(y,x)=1-\sum_y P_{yx} 这很合理。

其次,这样构造的PP 满足了Pyxπ(x)Pxyπ(y)P_{yx}*\pi(x)\equiv P_{xy}*\pi(y)。而一个重要的结论就是,满足这个恒等式(detailed balance)的马尔科夫链是遍历的,即其存在一个稳定的输出π(x)\pi(x),即我们需要抽样的随机分布。


所以说,我认为 MCMC 之所以可以模拟抽样非常复杂的分布,一个重要内容就是把复杂分布的连续信息,变成了单点 check 的信息。我通过一个简单的分布给出一个候选值,再根据π(x)\pi(x) 来决定是否采用,而不是根据π(x)\pi(x) 来得到候选值。

# 一个例子 —— 指数分布

譬如我想抽取满足概率密度函数f(x)f(x) 的随机变量XX 的一组样本:

f(x)={exx00elsef(x)=\begin{cases}e^{-x}&x\geq 0\\ 0&else\end{cases}

很显然它是连续的,无穷的,想进行抽样非常麻烦。

但我们可以令x1=10x_1=10,然后在U[103,10+3]U[10-3,10+3] 中随机抽取一个数,譬如10.310.3

那么计算α(10.3,10)=f(10)16=e106\alpha(10.3,10)=f(10)*\frac{1}{6}=\frac{e^{-10}}{6}。然后以这个概率进行二项分布抽样,譬如抽出来选择不接受。

x2=x1=10x_2=x_1=10

简单分析一下就会发现,这个过程如果重复足够多次,xx 是会大概率不断减小的。这就是指数分布在xx 较小时概率较大的原因。由于我们初始值x1=10x_1=10 这个取得不是很好,因此如果想要逼近稳定输出f(x)f(x),需要重复很多很多次。

譬如重复了 10000 次后,我们再取出x10000,...,x10010x_{10000},...,x_{10010} 这 11 个样本值,就可以当作指数分布的一组样本了。

Edited on Views times