受限玻尔兹曼机与MCMC-Gibbs采样计算

神经网络模型：

  1. 向前的神经网络DNN和CNN，
  2. 有反馈的神经网络RNN和LSTM。
  3. 玻尔兹曼机，此处主要关注受限玻尔兹曼机(Restricted Boltzmann Machine)，玻尔兹曼机主要应用领域在于推荐系统。

受限玻尔兹曼机

RBM模型结构

玻尔兹曼机是一大类的神经网络模型，但是在实际应用中主要使用受限玻尔兹曼机。RBM本身模型很简单，只是一个两层的神经网络，因此其不能算作深度学习。不过深度玻尔兹曼机(Deep Boltzmann Machine)可以看作是RBM的推广。理解了RBM再去理解DBM就比较简单了。

RBM的结构，是一个双层的神经网络，如下图所示：
image
上面一层神经元组成隐藏层(hidden layer)，我们使用\(h\)向量代表隐藏层神经元的值。下面一层的神经元组成可见层(visible layer)，用\(v\)向量代表可见层神经元的值。其中隐藏层和可见层是全连接的，这点和DNN类似，隐藏神经元之间是独立的，可见层神经元之间也是独立的。连接权重用矩阵W表示。和DNN的区别是，RBM不区分前向和反向，可见层的状态可以作用域隐藏层，而隐藏层的状态也可作用于可见层。隐藏层的偏执系数是向量\(b\)，而可见层的偏执系数是向量\(a\)。

常用的RBM一般是二值的，即不管是隐藏层还是可见层，他们的神经元取值只为\(0\)或\(1\)。本文只讨论而止RBM。

RBM模型结构的结构：主要是权重矩阵\(W\)，偏执系数向量a和b，隐藏层神经元状态向量h和可见层神经元状态向量v。

RBM概率分布

RBM是基于能量的概率分布。其分为两部分，第一部分是能量函数，第二部分是基于能量函数的概率分布函数 。

对于给定的状态向量h和v，则RBM当前的呢改良函数可以表示为；

\[E(v,h)=-a^Tv-b^Th-h^TWv \]

上面公式得到的结果是一个值，组成为可见层值*偏执，隐藏层值*偏执，可见层值*权重*隐藏层值。

有了关于\(v,h\)的能量函数\(E(v,h)\)，则我们可以定义RBM的状态为给定\(v,h\)的概率分布为：

\[P(v,h)=\frac{e^{-E(v,h)}}{\sum_{v,h}e^{-E(v,h)}} \]

有了概率分布之后，我们看一下条件分布

\[P(h|v)=\frac{P(h,v)}{P(v)}\\ =\frac{1}{P(V)}\frac{1}{Z}exp\{a^Tv+b^Th+h^TWv\}\\ =\frac{1}{Z'}exp\{b^Th+h^TWv\}\\ =\frac{1}{Z'}exp\{\sum_{j=1}^{nh}(b_j^Th_j+h_j^TW_j,:v)\}\\ =\frac{1}{Z'}\prod_{j=1}^{nh}exp\{\sum_{j=1}^{nh}(b_j^Th_j+h_j^TW_j,:v)\} \]

其中\(Z'\)为新的归一化系数，表达式为:

\[\frac{1}{Z'}=\frac{1}{P(v)}\frac{1}{Z}exp(a^Tv) \]

和上面同样的方式，我们也可以求出\(P(v|h)\)，这里不再列出。

有了条件概率分布，现在我们可以看一下RBM的激活函数，提到神经网络，一定逃不开激活函数，但是上面我们没有提到。由于使用的是能量概率模型，RBM的基于条件分布的激活函数是很容易推到出来的，我们以\(P(h_j=1|v)\)为例推导如下：

\[p(h_j=1|v)=\frac{P(h_j=1|v)}{P(h_j=1|v)+P(h_j=0|v)}\\ =\frac{exp\{b_j+W_j,:v\}}{exp\{0\}+exp\{b_j+W_j,:v\}}\\ =\frac{1}{1+exp\{-(b_j+W_j,:v)\}}\\ =sigmoid(b_j+W_j,:v) \]

从上面可看出，RBM里从可见层到隐藏层用的起始就是sigmoid激活函数。同样的方法，我们可以得到隐藏层到可见层用的也是sigmoid激活函数，即：

\[P(v_j=1|h)=sigmoid(a_j+W_j,:v) \]

有了激活函数，我们就可以从可见层参数推导出隐藏层神经元的取值概率了。对于\(0,1\)的取值情况，则大于0.5即取值为\(1\)。从隐藏层参数推导出可见的神经元取值方法也是一样的。

RBM模型的损失函数与优化

RBM模型的关键就是求出我们模型中参数\(W,a,b\)。如何求出呢？对于训练集的\(m\)个样本，RBM一般采用对数损失函数，即期望最小化下式：

\[L(W,a,b)=-\sum_{i=1}^mln(P(V^{(i)})) \]

对于优化过程，我们首先想到的当然是梯度下降法来迭代出\(W,a,b\)。我们首先来看单个样本的梯度计算，单个样本的损失函数为：\(-ln(P(V))\)具体的内容如下：

\[-ln(P(V))=-ln(\frac{1}{Z}\sum_he^{-E(V,h)})\\ =lnZ-ln(\sum_he^{-E(V,h)})\\ =ln(\sum_{v,h}e^{-E(v,h)})-ln(\sum_he^{-E(V,h)}) \]

注意这里面\(V\)表示的是某个特定的训练样本，而\(v\)指的是任意一个样本。
我们以\(a_i\)的梯度计算为例：

\[\frac{\partial(-ln(P(V)))}{\partial a_i}=\frac{1}{\partial a_i}\partial ln(\sum_{v,h}e^{-E(v,h)})-\frac{1}{\partial a_i}\partial ln(\sum_{h}e^{-E(V,h)})\\ =-\frac{1}{\sum_{v,h}e^{-E(v,h)}}(\sum_{v,h}e^{-E(v,h)}\frac{\partial E(v,h)}{\partial a_i})+\frac{1}{\sum_{v,h}e^{-E(V,h)}}(\sum_{V,h}e^{-E(V,h)}\frac{\partial E(V,h)}{\partial a_i})\\ =\sum_hP(h|V)\frac{\partial E(V,h)}{\partial a_i}-\sum_{v,h}P(h,v)\frac{\partial E(v,h)}{\partial a_i}\\ =-\sum_hP(h|V)V_i+\sum_{v,h}P(h,v)v_i\\ =-\sum_hP(h|V)V_i+\sum_vP(v)\sum_hP(h|v)v_i\\ =\sum_vP(v)v_i-V_i \]

其中用到了：

\[\sum_hP(h|v)=1 \:, \sum_hP(h|V)=1 \]

同样的话，可以得到\(W,b\)的梯度。这里就不求导了，直接给出结果：

\[\frac{\partial(-ln(P(V)))}{\partial b_i}=\sum_vP(V)P(h_i=1|v)-P(h_i=1|V)\\ \frac{\partial(-ln(P(V)))}{\partial W_{ij}}=\sum_vP(v)P(h_i=1|v)v_j-P(h_i=1|V)V_j \]

虽然在理论上梯度下降法可以解决RBM的优化问题，但是在实际应用中，由于概率分布的计算量很大，因为概率分布有\(2^{n_v+n_h}\)种情况，所以往往不直接按照上面的梯度求和公式去求所有样本的梯度和，而是基于MCMC的方法来模拟计算求解每个样本的梯度损失再求损失和，常用的方法是基于Gibbs采样的对比三度方法来求解，对于对比散度方法，需要MCMC的知识具体看这里

MCMC Gibbs采样计算方法

RBM推广到DBM

RBM很容易推广到深层的RBM，即DBM。推广的方法就是加入更多的隐藏层，比如一个三层的DBM如下：
image

当然隐藏层的层数可以实任意的，随着层数越来越复杂，那模型怎么表示呢？其实DBM也可以看作是一个非全连接的RBM，比如下图的一个4层的DBM，稍微加以改变就可以看作是一个RBM。
image
将可见层和偶数层放到左边，将奇数层放到右边，我们就得到了RBM，和之前的RBM的区别就是此处的RBM不是全连接的，起始可以看作部分权重为0的全连接RBM。RBM的算法思想可以在DBM上使用，只是此时模型参数更多，而且迭代求解参数也更加复杂。

RBM小结

RBM所在的玻尔兹曼机流派是深度学习中三大流派之一，也是目前比较热门的创新区域之一，目前在实际应用中的比较成功的是推荐系统。

蒙特卡洛方法

蒙特卡洛是一个赌场的名称，用它作为名字大概是因为蒙特卡罗方法是一种随机模拟的方法，这很像赌博场里面的扔骰子的过程。最早的蒙特卡罗方法都是为了求解一些不太好求解的求和或者积分问题。比如积分：

\[\theta=\int_a^bf(x)dx \]

一个简单的近似求解方法是在\([a,b]\)之间随机的采样一个点。比如\(x_0\),然后用\(f(x_0)\)代表在\([a,b]\)区间上所有的\(f(x)\)的值。那么上面的定积分的近似求解为:

\[(b-a)f(x_0) \]

用一个值代表\([a,b]\)区间上所有的\(f(x)\)的值，这个假设太粗糙。那么我们可以采样\([a,b]\)区间的\(n\)个值：\(x_0,x_1,\dots x_{n−1}\),用它们的均值来代表\([a,b]\)区间上所有的\(f(x)\)的值。这样我们上面的定积分的近似求解为:

\[\int_a^bf(x)dx \approx \frac{b-a}{n}\sum_{i=1}^{n-1}f(x_i) \]

上面的方法虽然可以求解，但是它隐含了一个假定，即\(x\)在\([a,b]\)之间是均匀分布的，然而在大部分情况下刚好相反，为了解决这个问题，我们可以假定我们得到了\(x\)在\([a,b]\)之间的概率分布函数\(p(x)\)，那么我们的定积分可以这样进行:

\[\theta=\int_a^bf(x)dx=\int_a^b\frac{f(x)}{p(x)}p(x)dx \approx \color{red}{\frac{1}{n}\sum_{i=0}^{n-1}\frac{f(x_i)}{p(x_i)}} \]

上面公式最右边的这个形式就是蒙特卡洛方法的一般形式。

可以看出的是，最上面我们假设\(x\)在\([a,b]\)之间是均匀分布的时候，\(p(x_i)=\frac{1}{b-a}\)，带入我们上面的带有概率分布的蒙特卡洛积分公式可以得到：

\[\frac{1}{n}\sum_{i=0}^{n-1}\frac{f(x_i)}{\frac{1}{b-a}}=\frac{b-a}{n}\sum_{i=0}^{n-1}f(x_i) \]

也就是说均匀分布可以看作是蒙特卡洛积分积分概率为均匀概率情况下的特例。

\(\color{red}{我们的目的是求积分，但是根据上面的蒙特卡洛积分公式得到的结果，可以确定，我们目前的问题转换为了，如何求出x的分布p(x)对应的若干个样本上。\\只要采样足够多，我们就可以近似出来结果。}\)

概率分布采样

话接上回，我们讲到蒙特卡洛方法求解的关键是得到\(x\)的概率分布，如果求出了\(x\)的概率分布，我们可以基于概率分布去采样基于这个概率分布的\(n\)个\(x\)样本。

对于常见的均匀分布uniform(0,1)是很容易采样的，一般通过线性同余发生器可以很方便的生成\((0,1)\)之间的伪随机数样本。而其他常见的概率分布，无论是连续的还是离散的分布，它们的样本都可以通过uniform(0,1)的样本转换得到，比如二维正态分布的样本\((Z_1,Z_2)\)可以通过独立采样得到的uniform(0,1)样本对\((X_1,X_2)\)通过如下式子的转换而得:

\[Z_1=\sqrt{-2lnX_1cos(2\pi X_2)}\\ Z_2=\sqrt{-2lnX_1sin(2\pi X_2)} \]

对于一些常见的连续分布，比如\(t\)分布,\(F\)分布,\(Beta\)分布等，都可以通过类似的方式从uniform(0,1)的样本转换而得，比如二维正态分布样本\((Z_1,Z_2)\)可以通过独立采样转换得到，在python的numpy,scikit-learn等类库中，都有生成这些常用分布样本的函数可以使用。

不过在大多数情况下，\(x\)的分布不是常见的分布，这意味着我们没办法方便的得到这些常见的概率分布的采样样本集，这个问题怎么解决呢？

接受-拒绝采样

对于概率分布不是常见的分布，一个可行的办法是采用接受-拒绝采样来得到该分布的样本。既然p(x)太过复杂没有办法采样得到，那么我们可以设定一个程序可采样的分布\(q(x)\)比如高斯分布，然后按照一定的方法拒绝某些样本，以达到接近\(p(x)\)分布的目的，其中\(q(x)\)叫做proposal distribution
image

具体采样过程如下，设定一个方便采样的的常用概率分布函数q(x)以及一个常量k，使得p(x)总在kq(x)的下方，如上图。

首先，采样得到q(x)的一个样本\(z_0\)，采样方法如上一节，然后，从均匀分布\((0,kq(z_0))\)中采样得到一个值\(u\)。如果\(u\)落在了上图中的灰色区域，则拒绝这次抽样，否则接受这个样本\(z_0\)。重复以上过程得到\(n\)个接受的样本\(z_0,z_1,\dots z_{n−1}\),则最后的蒙特卡罗方法求解结果为：

\[\frac{1}{n}\sum_{i=0}^{n-1}\frac{f(z_i)}{p(z_i)} \]

整个过程中，我们通过一系列的接受拒绝决策来达到用q(x)模拟p(x)概率分布的目的。

蒙特卡洛小结

使用接受-拒绝采样，我们可以解决一些概率分布不是常见的分布的时候，得到其采样集并使用蒙特卡洛方法求和的目的，但是接受拒绝采样也只能部分满足我们的需求，在很多时候我们还是很难的到我们的概率分布的样本集，比如：

  1. 对于一些二维分布\(p(x',y)\)，有时候我们只能得到条件分布\(p(x|y)\)和\(p(y|x)\)，却很难得到二维分布\(P(x,y)\)一般形式，这是我们无法用接受-拒绝采样得到其样本集。
  2. 对于一些\(n\)维复杂非常见样本分布\(p(x_1,x_2,\dots,x_n)\)，我们要找到一个合适的q(x)和k非常困难。

由上可得，要想将蒙特卡洛方法作为一个通用的采样模拟求和的方法，必须解决如何渐变的得到各种复杂概率分布的而对应的采样样本集的问题，下面的马尔科夫链就是帮助我们找到这些复杂概率分布的对应的采样样本集合的方法。

马尔科夫链

概述

  • 使用蒙特卡洛方法去随机模拟一些复杂的连续积分或者离散求和的方法，但是这个方法需要得到\(\color{red}{对应的概率分布的样本集}\)，而想要得到这样的样本集非常困难，因此我们需要马尔科夫链去帮忙。

马尔科夫链定义本身比较简单，它假设某一时刻状态转移的概率只依赖于他的前一个状态，举个形象的比喻，假如每天的天气是一个状态的话，那么今天的天气情况只依赖于昨天的天气，而与前天的天气无关，以此去简化模型的复杂度，因此马尔科夫链在很多时间序列模型中得到了广泛的应用，比如\(\color{red}{RNN，HNN，MCMC}\)

如果使用精确的数学定义来描述，则假设我们的序列状态是\(\dots X_{t-2},X_{t-1},X_t,X_{t+1},\dots\)，那么我们的时刻在\(X_{t+1}\)的状态的条件概率仅仅依赖于时刻\(X_t\)，即：

\[P(X_{t+1}|\dots X_{t-2},X_{t-1},X_t)=P(X_{t+1}|X_t) \]

\(\color{green}{称为在X_t的条件下X_{t+1}发生的概率}\)

既然在某一时刻状态转移的概率只依赖于它的前一个状态，那么我们只要能求出系统中任意两个状态之间的转移概率，这个马尔科夫链的模型就定了，我们来看看下图这个马尔可夫模型的具体的例子。

image

这个马尔科夫链是表示股市模型的，共有三种状态：牛市(bull market)，熊市(bear market)，横盘(stagnant market)，每一个状态都以一定的概率转换到下一个状态，比如牛市以\(0.025\)的概率转换到横盘的状态。这个状态转化图可以用矩阵的形式表示，如果我们定义矩阵\(P\)某一位置\(P(i,j)\)的值为\(P(j|i)\)，即从状态\(i\)转换到状态\(j\)的概率，并定义牛市的状态为\(0\)，熊市为状态\(1\)，横盘状态\(2\)，这样我们就得到了状态转移矩阵为：

\[P = \begin{pmatrix} 0.9&0.075&0.025\\ 0.15&0.8&0.05\\ 0.25&0.25&0.5\\ \end{pmatrix}\]

在基本了解了马尔科夫链模型的状态转移矩阵之后，我们需要将其和蒙特卡洛方法联系起来，这时候我们需要从马尔科夫链模型的状态转移矩阵的性质讲起。

马尔科夫链模型状态转移矩阵的性质

我们仍然上面的这个状态转移矩阵为例，假设我们当前故事的概率分布为：\([0.3,0.4,0.3]\)，即\(30\%\)的的牛市，\(40\%\)概率的熊市和\(30\%\)的横盘，这个状态作为初始状态\(t_0\)，将其带入这个状态转移矩阵计算\(t_1,t_2,t_3,\dots\)的状态，代码如下：

点击查看代码 
import numpy as np
matrix = np.matrix([[0.9,0.075,0.025],[0.15,0.8,0.05],[0.25,0.25,0.5]], dtype=float)
vector1 = np.matrix([[0.3,0.4,0.3]], dtype=float)
for i in range(100):
    vector1 = vector1*matrix
    print ("Current round:" , i+1)
    print (vector1)
输出结果如下：
点击查看代码 
Current round: 1
[[0.405  0.4175 0.1775]]
Current round: 2
[[0.4715  0.40875 0.11975]]
Current round: 3
[[0.5156 0.3923 0.0921]]
Current round: 4
[[0.54591  0.375535 0.078555]]
Current round: 5
[[0.567288 0.36101  0.071702]]
Current round: 6
[[0.5826362 0.3492801 0.0680837]]
Current round: 7
[[0.59378552 0.34014272 0.06607176]]
Current round: 8
[[0.60194632 0.33316603 0.06488765]]
Current round: 9
[[0.6079485  0.32790071 0.06415079]]
Current round: 10
[[0.61237646 0.3239544  0.06366914]]
Current round: 11
[[0.61564926 0.32100904 0.0633417 ]]
Current round: 12
[[0.61807111 0.31881635 0.06311253]]
Current round: 13
[[0.61986459 0.31718655 0.06294886]]
Current round: 14
[[0.62119333 0.3159763  0.06283037]]
Current round: 15
[[0.62217803 0.31507813 0.06274383]]
Current round: 16
[[0.62290791 0.31441182 0.06268027]]
Current round: 17
[[0.62344896 0.31391762 0.06263343]]
Current round: 18
[[0.62385006 0.31355112 0.06259882]]
Current round: 19
[[0.62414743 0.31327936 0.06257322]]
Current round: 20
[[0.62436789 0.31307785 0.06255426]]
Current round: 21
[[0.62453135 0.31292843 0.06254022]]
Current round: 22
[[0.62465253 0.31281765 0.06252982]]
Current round: 23
[[0.62474238 0.31273552 0.0625221 ]]
Current round: 24
[[0.624809   0.31267462 0.06251639]]
Current round: 25
[[0.62485839 0.31262947 0.06251215]]
Current round: 26
[[0.624895   0.31259599 0.06250901]]
Current round: 27
[[0.62492215 0.31257117 0.06250668]]
Current round: 28
[[0.62494228 0.31255277 0.06250495]]
Current round: 29
[[0.62495721 0.31253912 0.06250367]]
Current round: 30
[[0.62496827 0.31252901 0.06250272]]
Current round: 31
[[0.62497648 0.31252151 0.06250202]]
Current round: 32
[[0.62498256 0.31251594 0.0625015 ]]
Current round: 33
[[0.62498707 0.31251182 0.06250111]]
Current round: 34
[[0.62499041 0.31250876 0.06250082]]
Current round: 35
[[0.62499289 0.3125065  0.06250061]]
Current round: 36
[[0.62499473 0.31250482 0.06250045]]
Current round: 37
[[0.62499609 0.31250357 0.06250034]]
Current round: 38
[[0.6249971  0.31250265 0.06250025]]
Current round: 39
[[0.62499785 0.31250196 0.06250018]]
Current round: 40
[[0.62499841 0.31250146 0.06250014]]
Current round: 41
[[0.62499882 0.31250108 0.0625001 ]]
Current round: 42
[[0.62499912 0.3125008  0.06250008]]
Current round: 43
[[0.62499935 0.31250059 0.06250006]]
Current round: 44
[[0.62499952 0.31250044 0.06250004]]
Current round: 45
[[0.62499964 0.31250033 0.06250003]]
Current round: 46
[[0.62499974 0.31250024 0.06250002]]
Current round: 47
[[0.6249998  0.31250018 0.06250002]]
Current round: 48
[[0.62499985 0.31250013 0.06250001]]
Current round: 49
[[0.62499989 0.3125001  0.06250001]]
Current round: 50
[[0.62499992 0.31250007 0.06250001]]
Current round: 51
[[0.62499994 0.31250005 0.06250001]]
Current round: 52
[[0.62499996 0.31250004 0.0625    ]]
Current round: 53
[[0.62499997 0.31250003 0.0625    ]]
Current round: 54
[[0.62499998 0.31250002 0.0625    ]]
Current round: 55
[[0.62499998 0.31250002 0.0625    ]]
Current round: 56
[[0.62499999 0.31250001 0.0625    ]]
Current round: 57
[[0.62499999 0.31250001 0.0625    ]]
Current round: 58
[[0.62499999 0.31250001 0.0625    ]]
Current round: 59
[[0.62499999 0.3125     0.0625    ]]
Current round: 60
[[0.625  0.3125 0.0625]]
Current round: 61
[[0.625  0.3125 0.0625]]
Current round: 62
[[0.625  0.3125 0.0625]]
Current round: 63
[[0.625  0.3125 0.0625]]
Current round: 64
[[0.625  0.3125 0.0625]]
Current round: 65
[[0.625  0.3125 0.0625]]
Current round: 66
[[0.625  0.3125 0.0625]]
Current round: 67
[[0.625  0.3125 0.0625]]
Current round: 68
[[0.625  0.3125 0.0625]]
Current round: 69
[[0.625  0.3125 0.0625]]
Current round: 70
[[0.625  0.3125 0.0625]]
Current round: 71
[[0.625  0.3125 0.0625]]
Current round: 72
[[0.625  0.3125 0.0625]]
Current round: 73
[[0.625  0.3125 0.0625]]
Current round: 74
[[0.625  0.3125 0.0625]]
Current round: 75
[[0.625  0.3125 0.0625]]
Current round: 76
[[0.625  0.3125 0.0625]]
Current round: 77
[[0.625  0.3125 0.0625]]
Current round: 78
[[0.625  0.3125 0.0625]]
Current round: 79
[[0.625  0.3125 0.0625]]
Current round: 80
[[0.625  0.3125 0.0625]]
Current round: 81
[[0.625  0.3125 0.0625]]
Current round: 82
[[0.625  0.3125 0.0625]]
Current round: 83
[[0.625  0.3125 0.0625]]
Current round: 84
[[0.625  0.3125 0.0625]]
Current round: 85
[[0.625  0.3125 0.0625]]
Current round: 86
[[0.625  0.3125 0.0625]]
Current round: 87
[[0.625  0.3125 0.0625]]
Current round: 88
[[0.625  0.3125 0.0625]]
Current round: 89
[[0.625  0.3125 0.0625]]
Current round: 90
[[0.625  0.3125 0.0625]]
Current round: 91
[[0.625  0.3125 0.0625]]
Current round: 92
[[0.625  0.3125 0.0625]]
Current round: 93
[[0.625  0.3125 0.0625]]
Current round: 94
[[0.625  0.3125 0.0625]]
Current round: 95
[[0.625  0.3125 0.0625]]
Current round: 96
[[0.625  0.3125 0.0625]]
Current round: 97
[[0.625  0.3125 0.0625]]
Current round: 98
[[0.625  0.3125 0.0625]]
Current round: 99
[[0.625  0.3125 0.0625]]
Current round: 100
[[0.625  0.3125 0.0625]]

可以发现，从第\(60\)轮开始，我们的状态概率就不变了，一直保持在\([0.625,0.3125,0.0625]\)，即\(0.625\)的牛市，\(0.3125\)的熊市，\(0.0625\)的横盘。我们现在换一个初始概率分布试一下\([0.7,0.1,0.2]\)重新带入如下：

点击查看代码 
Current round: 1
[[0.695  0.1825 0.1225]]
Current round: 2
[[0.6835  0.22875 0.08775]]
Current round: 3
[[0.6714 0.2562 0.0724]]
Current round: 4
[[0.66079  0.273415 0.065795]]
Current round: 5
[[0.652172 0.28474  0.063088]]
Current round: 6
[[0.6454378 0.2924769 0.0620853]]
Current round: 7
[[0.64028688 0.29791068 0.06180244]]
Current round: 8
[[0.6363954  0.30180067 0.06180393]]
Current round: 9
[[0.63347695 0.30462117 0.06190188]]
Current round: 10
[[0.6312979  0.30668318 0.06201892]]
Current round: 11
[[0.62967532 0.30819862 0.06212607]]
Current round: 12
[[0.62846909 0.30931606 0.06221485]]
Current round: 13
[[0.6275733  0.31014174 0.06228495]]
Current round: 14
[[0.62690847 0.31075263 0.0623389 ]]
Current round: 15
[[0.62641525 0.31120496 0.06237979]]
Current round: 16
[[0.62604941 0.31154006 0.06241053]]
Current round: 17
[[0.62577811 0.31178839 0.0624335 ]]
Current round: 18
[[0.62557693 0.31197244 0.06245062]]
Current round: 19
[[0.62542776 0.31210888 0.06246336]]
Current round: 20
[[0.62531716 0.31221003 0.06247282]]
Current round: 21
[[0.62523515 0.31228501 0.06247984]]
Current round: 22
[[0.62517435 0.31234061 0.06248505]]
Current round: 23
[[0.62512926 0.31238182 0.06248891]]
Current round: 24
[[0.62509584 0.31241238 0.06249178]]
Current round: 25
[[0.62507106 0.31243504 0.0624939 ]]
Current round: 26
[[0.62505268 0.31245184 0.06249548]]
Current round: 27
[[0.62503906 0.31246429 0.06249665]]
Current round: 28
[[0.62502896 0.31247352 0.06249752]]
Current round: 29
[[0.62502147 0.31248037 0.06249816]]
Current round: 30
[[0.62501592 0.31248545 0.06249863]]
Current round: 31
[[0.6250118  0.31248921 0.06249899]]
Current round: 32
[[0.62500875 0.312492   0.06249925]]
Current round: 33
[[0.62500649 0.31249407 0.06249944]]
Current round: 34
[[0.62500481 0.3124956  0.06249959]]
Current round: 35
[[0.62500357 0.31249674 0.06249969]]
Current round: 36
[[0.62500264 0.31249758 0.06249977]]
Current round: 37
[[0.62500196 0.31249821 0.06249983]]
Current round: 38
[[0.62500145 0.31249867 0.06249988]]
Current round: 39
[[0.62500108 0.31249901 0.06249991]]
Current round: 40
[[0.6250008  0.31249927 0.06249993]]
Current round: 41
[[0.62500059 0.31249946 0.06249995]]
Current round: 42
[[0.62500044 0.3124996  0.06249996]]
Current round: 43
[[0.62500033 0.3124997  0.06249997]]
Current round: 44
[[0.62500024 0.31249978 0.06249998]]
Current round: 45
[[0.62500018 0.31249984 0.06249998]]
Current round: 46
[[0.62500013 0.31249988 0.06249999]]
Current round: 47
[[0.6250001  0.31249991 0.06249999]]
Current round: 48
[[0.62500007 0.31249993 0.06249999]]
Current round: 49
[[0.62500005 0.31249995 0.0625    ]]
Current round: 50
[[0.62500004 0.31249996 0.0625    ]]
Current round: 51
[[0.62500003 0.31249997 0.0625    ]]
Current round: 52
[[0.62500002 0.31249998 0.0625    ]]
Current round: 53
[[0.62500002 0.31249999 0.0625    ]]
Current round: 54
[[0.62500001 0.31249999 0.0625    ]]
Current round: 55
[[0.62500001 0.31249999 0.0625    ]]
Current round: 56
[[0.62500001 0.31249999 0.0625    ]]
Current round: 57
[[0.625  0.3125 0.0625]]
Current round: 58
[[0.625  0.3125 0.0625]]
Current round: 59
[[0.625  0.3125 0.0625]]
Current round: 60
[[0.625  0.3125 0.0625]]
Current round: 61
[[0.625  0.3125 0.0625]]
Current round: 62
[[0.625  0.3125 0.0625]]
Current round: 63
[[0.625  0.3125 0.0625]]
Current round: 64
[[0.625  0.3125 0.0625]]
Current round: 65
[[0.625  0.3125 0.0625]]
Current round: 66
[[0.625  0.3125 0.0625]]
Current round: 67
[[0.625  0.3125 0.0625]]
Current round: 68
[[0.625  0.3125 0.0625]]
Current round: 69
[[0.625  0.3125 0.0625]]
Current round: 70
[[0.625  0.3125 0.0625]]
Current round: 71
[[0.625  0.3125 0.0625]]
Current round: 72
[[0.625  0.3125 0.0625]]
Current round: 73
[[0.625  0.3125 0.0625]]
Current round: 74
[[0.625  0.3125 0.0625]]
Current round: 75
[[0.625  0.3125 0.0625]]
Current round: 76
[[0.625  0.3125 0.0625]]
Current round: 77
[[0.625  0.3125 0.0625]]
Current round: 78
[[0.625  0.3125 0.0625]]
Current round: 79
[[0.625  0.3125 0.0625]]
Current round: 80
[[0.625  0.3125 0.0625]]
Current round: 81
[[0.625  0.3125 0.0625]]
Current round: 82
[[0.625  0.3125 0.0625]]
Current round: 83
[[0.625  0.3125 0.0625]]
Current round: 84
[[0.625  0.3125 0.0625]]
Current round: 85
[[0.625  0.3125 0.0625]]
Current round: 86
[[0.625  0.3125 0.0625]]
Current round: 87
[[0.625  0.3125 0.0625]]
Current round: 88
[[0.625  0.3125 0.0625]]
Current round: 89
[[0.625  0.3125 0.0625]]
Current round: 90
[[0.625  0.3125 0.0625]]
Current round: 91
[[0.625  0.3125 0.0625]]
Current round: 92
[[0.625  0.3125 0.0625]]
Current round: 93
[[0.625  0.3125 0.0625]]
Current round: 94
[[0.625  0.3125 0.0625]]
Current round: 95
[[0.625  0.3125 0.0625]]
Current round: 96
[[0.625  0.3125 0.0625]]
Current round: 97
[[0.625  0.3125 0.0625]]
Current round: 98
[[0.625  0.3125 0.0625]]
Current round: 99
[[0.625  0.3125 0.0625]]
Current round: 100
[[0.625  0.3125 0.0625]]

可以看到尽管上一个是在第\(60\)轮达到稳定状态，目前这个在\(57\)轮达到稳定状态，但是最终的稳定状态是一致的。
也就是说我们的马尔科夫链模型的状态转移矩阵是否收敛到稳定概率分布与初始状态无关。这是一个非常好的性质，也就是说，如果我们得到了这个稳定概率分布对应的马尔科夫链模型的状态转移矩阵，则我们可以用任意的概率分布样本开始，带入马尔科夫链模型的状态转移矩阵，这样经过一系列的转换，最终就可以得到符合对于稳定概率分布的样本。

这个性质不单单对于状态转移矩阵有效，对于绝大多数的马尔科夫链模型的状态转移矩阵也有效，同时不光是离散状态，\(\color{blue}{连续状态是也成立}\)。

同时对于一个确定的状态转移矩阵\(P\)，他的\(n\)次幂\(P^n\)在当\(n\)大于一定值的时候也可以发现是确定的。计算代码如下：

点击查看代码 
matrix = np.matrix([[0.9,0.075,0.025],[0.15,0.8,0.05],[0.25,0.25,0.5]], dtype=float)
for i in range(10):
    matrix = matrix*matrix
    print ("Current round:" , i+1)
    print (matrix)
结果如下：
点击查看代码 
Current round: 1
[[0.8275  0.13375 0.03875]
 [0.2675  0.66375 0.06875]
 [0.3875  0.34375 0.26875]]
Current round: 2
[[0.73555  0.212775 0.051675]
 [0.42555  0.499975 0.074475]
 [0.51675  0.372375 0.110875]]
Current round: 3
[[0.65828326 0.28213131 0.05958543]
 [0.56426262 0.36825403 0.06748335]
 [0.5958543  0.33741675 0.06672895]]
Current round: 4
[[0.62803724 0.30972343 0.06223933]
 [0.61944687 0.3175772  0.06297594]
 [0.6223933  0.3148797  0.062727  ]]
Current round: 5
[[0.62502532 0.31247685 0.06249783]
 [0.6249537  0.31254233 0.06250397]
 [0.62497828 0.31251986 0.06250186]]
Current round: 6
[[0.625  0.3125 0.0625]
 [0.625  0.3125 0.0625]
 [0.625  0.3125 0.0625]]
Current round: 7
[[0.625  0.3125 0.0625]
 [0.625  0.3125 0.0625]
 [0.625  0.3125 0.0625]]
Current round: 8
[[0.625  0.3125 0.0625]
 [0.625  0.3125 0.0625]
 [0.625  0.3125 0.0625]]
Current round: 9
[[0.625  0.3125 0.0625]
 [0.625  0.3125 0.0625]
 [0.625  0.3125 0.0625]]
Current round: 10
[[0.625  0.3125 0.0625]
 [0.625  0.3125 0.0625]
 [0.625  0.3125 0.0625]]

我们可以发现，在\(n\geq 6\)以后，\(P^n\)的值稳定不再变化，而且每一行都为\([0.625,0.3125,0.0625]\)，这和前面的稳定分布是一致的，这个性质不单单是离散状态，\(\color{blue}{连续状态是也成立}\)。

现在我们可以使用数学语言总结马尔科夫链的收敛性之：

如果一个非周期性的马尔科夫链有状态转移矩阵\(P\)，并且他的任何两个状态是连通的，那么\(\lim\limits_{x \to \infty} P_{ij}^n\)与\(i\)无关，我们有：

\[\lim\limits_{x \to \infty} P_{ij}^n = \pi(j)\tag{1} \]

\[\lim\limits_{x \to \infty} P^n= \begin{pmatrix} \pi(1) & \pi(2) & \dots & \pi(j) & \dots\\ \pi(1) & \pi(2) & \dots & \pi(j) & \dots\\ \dots & \dots & \dots & \dots & \dots\\ \pi(1) & \pi(2) & \dots & \pi(j) & \dots\\ \dots & \dots & \dots & \dots & \dots\\ \end{pmatrix}\tag{2} \]

\[\pi(j)=\sum_{i=0}^\infty \pi(i)P_{ij} \tag{3} \]

\(\pi\)是方程\(\pi P=\pi\)的唯一非负解，其中：

\[\pi=[\pi(1),\pi(2),\dots,\pi(j),\dots] \quad \sum_{j=0}^\infty \pi(j)=1 \]

上面的性质需要解释的有：

  1. 非周期性的马尔科夫链：主要指马尔科夫链的状态转化不是循环的，如果是循环的则永远不会收敛。幸运的是我们遇到的马尔科夫链一般都是非周期性的。
  2. 任何两个状态之间是连通的：这个指的是从任意一个状态可以通过有限步到达其他的任意一个状态，不会出现条件概率率一直为0导致不可到达的情况发生。
  3. 马尔科夫链的状态数是有限的，也可以是无限的。\(\color{blue}{因此可以用于连续概率分布和离散概率分布}\)。
  4. \(\pi\)通常称为马尔可夫链的平稳分布。

基于马尔科夫链采样

蒙特卡洛方法已经讲的很清楚了，关键在于求出\(x\)的概率分布，本节我们又补充了马尔科夫链的基础知识，下面就讲讲如何基于马尔可夫链进行采样。

如果我们得到了某个平稳分布所对应的马尔科夫链状态转移矩阵，我们就可以很容易采样出这个平稳分布的样本集。

\(\color{red}{通过平稳分布获取到马尔可夫链转台转移矩阵矩阵，容易采样该平稳分布的样本集}\)

假设我们任意初始的概率分布是\(\pi_0(x)\)，经过第一轮马尔可夫链状态转移矩阵计算之后得到的概率分布是\(\pi_1(x)\),第\(i\)轮得到的是\(\pi_i(x)\)。假设经过\(n\)轮之后马尔可夫链收敛到我们的平稳分布\(\pi(x)\)即：

\[\pi_n(x)=\pi_{n+1}(x)=\pi_{n+2}(x)=\dots=\pi(x) \]

对于每个分布\(\pi_i(x)\)，我们有：

\[\pi_i(x)=\pi_{i-1}(x)P=\pi_{i-2}(x)P^2=\pi_{0}(x)P^i \]

\(\color{brown}{现在我们可以进行采样工作}\)，首先,基于初始任意简单概率分布比如高斯分布\(\pi_0(x)\)采样得到状态值\(x_0\)，基于条件概率分布\(P(x|x_0)\)采样状态值\(x_1\)，一直进行下去，当状态转移到一定次数之后，比如到\(n\)此时，我们认为此时的采样集和\((x_n,x_{n+1},x_{n+2},\dots)\)即是符合我们平稳分布对应的样本集，可以用来做蒙特卡洛模拟求和了。

总结一下马尔可夫链的采样过程：

  1. 输入马尔科夫链状态转移矩阵\(P\)，设定状态转移次数阈值\(n_1\)，需要的样本数量\(n_2\)。
  2. 从任意简单概率分布采样得到初始状态值\(x_0\)
  3. \(for\quad t=0\quad to\quad n_1+n_2-1\)：从概率分布\(P(x|x_t)\)中采样得到样本\(x\)。样本集\(x_{n_1},x_{n_1+1},x_{n_1+2,\dots,x_{n_1+n_2-1}}\)就是我们需要的平稳分布对应的样本集。

???状态转移矩阵是怎么获得的？ 前面讲的是如果我们得到了某个平稳分布所对应的马尔科夫链状态转移矩阵，然后有了后面的东西， 现在还没有讲，状态转移矩阵是怎么得到的。

马尔可夫链小结

如果假定我们可以得到我们需要采样样本的平稳分布所对应的马尔科夫链状态转移矩阵，那么我们就可以用马尔可夫链采样得到我们需要的样本集，进行蒙特卡洛模拟，但是一个重要的问题是，随意给定一个平稳分布\(\pi\)如何得到它所对应的马尔可夫链状态转移矩阵\(P\)呢？这是一个大问题，我们绕了一圈把解决蒙特卡洛的问题转到--求采样集和-转到马尔科夫链状态转移矩阵上。

马尔可夫链的细致平稳条件

  • 解决积分不可求解问题。
  1. 使用蒙特卡洛方法进行近似求解，但是关于\(x\)的样本集输入我们无法解决。
  2. 使用马尔可夫链采样解决\(x\)的样本集问题，但是马尔可夫状态转移矩阵无法求解。
  3. 此处解决马尔科夫链状态转移矩阵问题。

在上一节马尔可夫链采样中我们讲到给定一个概率平稳分布\(\pi\)，很难直接找到对应的马尔科夫链状态转移矩阵\(P\)，而只要解决这个问题，我们就可以找到一种\(\color{red}{通用}\)的概率分布采样方法，进而用于蒙特卡洛模拟，这里我们就介绍解决这个问题的方法：MCMC采样和他的建议版本M-H采样。

对于任意给定的目标分布\(\pi(x)\)，我们如何找到以它为唯一平稳分布的马尔可夫链，并且基于马尔可夫链采样的方法对其进行近似采样呢？

找到这样一个马尔可夫链，本质就是要寻找他的转移概率矩阵\(P\)，那么首先确立一个思考路径：有没有什么条件，使得只要我们的转移矩阵\(P\)满足了，就意味着目标分布\(\pi(x)\)就是转移矩阵\(P\)对应的马尔可夫链的平稳分布呢？

在解决从平稳分布\(\pi\)，找到对应的马尔可夫链状态转移矩阵\(P\)之前，我们需要看一下，\(\color{red}{马尔可夫链的细致平稳条件}\)。定义如下：
如果非周期性马尔可夫链的状态转移矩阵\(P\)和概率分布\(\pi(x)\)对于所有的\(i,j\)满足:

\[\pi(i)P(i,j)=\pi(j)P(j,i)\quad \forall i,j \tag{4} \]

则称概率分布\(\pi(x)\)是状态转移矩阵\(P\)的平稳分布。
证明如下，我们同时对等式4两边进行求和，由细致平稳条件得出：

\[\sum_{i=1}^\infty \pi(i)P(i,j)=\sum_{i=1}^\infty \pi(j)P(j,i)\tag{5} \]

观察等式6右边部分，发现\(\pi(j)\)的取值是独立于参数\(i\)，可以将其提出，并且由于马尔科夫链状态转移矩阵的性质可得：

\[\sum_{i=1}^\infty \pi(i)P(i,j)=\sum_{i=1}^\infty \pi(j)P(j,i)=\pi_j\sum_i^\infty P(j,i)=\pi_j\tag{6} \]

将公式6用矩阵可表示为：\(\pi P=\pi\)
即满足马尔可夫链的收敛性质，也就是说，只要我们找到了可以使概率分布\(\pi(x)\)满足细致平稳条件分布的矩阵\(P\)即可。这给了我们从平稳分布\(\pi\)找到对应的马尔可夫链状态转移矩阵\(P\)的新思路。

不过不幸的是，仅仅从细致平稳条件还是很难找到合适的矩阵\(P\)，比如我们的目标平稳分布是\(\pi(x)\)，随即找一个马尔科夫链状态转移矩阵\(Q\)，他是很难满足细致平稳条件的，即：

\[\pi(i)P(i,j) \neq \pi(j)P(j,i)\quad \forall i,j \]

那么如何是\(P\)满足上面的等式呢？ 下面使用MCMC采样去解决这个问题。

MCMC采样

由于一般情况下，目标平稳分布\(\pi(x)\)和某一个马尔可夫链状态转移矩阵\(Q\)不满足细致平稳条件即：

\[\pi(i)P(i,j) \neq \pi(j)P(j,i)\quad \forall i,j \]

我们可以对上面的式子做一个改造，使细致平稳条件成立，方法是引入一个\(\alpha(i,j)\)使上式可以取等号，即：

\[\pi(i)P(i,j)\alpha(i,j) = \pi(j)P(j,i)\alpha(j,i)\quad \forall i,j \tag{7} \]

问题是什么样的\(\alpha\)矩阵可以使等式7成立呢？

\[\alpha(i,j)=\pi(j)Q(j,i) \]

\[\alpha(j,i)=\pi(i)Q(i,j) \]

这样我们就得到了满足我们分布\(\pi(x)\)的对应的马尔可夫链状态转移矩阵\(P\)，满足：

\[P(i,j)=Q(i,j)\alpha(i,j)\tag{8} \]

这里的马尔科夫链状态转移矩阵\(Q\)是随意的，\(\pi\)是已经给出的内容，$\alpha可以从上面的到，所所以上面的公式8，可以由此推出。

也就是书，我们的目标矩阵即马尔科夫链状态转移矩阵\(P\)，可以通过任意一个马尔科夫链状态转移矩阵\(Q\)乘以\(\alpha(i,j)\)获得，\(\alpha(i,j)\)我们一般称之为接受率，取值在\([0,1]\)，可以理解为一个概率值，即目标矩阵\(P\)可以通过任意一个马尔科夫链状态转移矩阵\(Q\)以一定的接受率获得，这个很想我们在之前看的那个接受-拒绝采样哪里是一个常见分布函数通过一定的接受-拒绝概率得到一个非常见分布的概率函数，这里是一个常见的马尔科夫链状态转移矩阵\(Q\)通过一定的接受-拒绝概率得到目标转移矩阵\(P\)，两者的解决问题的思路是一致的。

下面描述以下MCMC采样的过程。

  1. 输入我们任意选定的马尔可夫状态转移矩阵\(Q\)，平稳分布\(\pi(x)\)，设定状态转移阈值\(n_1\)，需要的样本个数\(n_2\)
  2. 从任意简单概率分布采样得到初始状态值\(x_0\)
  3. \(for t=0 to n_1+n_2-1:\)
    1. 从条件分布\(Q(x|x_t)\)中采样得到样本\(x_*\)
    2. 从均匀分布采样\(u\sim uniform[0,1]\)
    3. 如果\(u<\alpha(x_t,x_*)=\pi(x_*)Q(x_*,x_t)\)，则接受转移\(x_t\Rightarrow x_*\)，则\(x_{t+1}=x_*\)
    4. 否则不接受转移，即\(x_{t+1}=x_t\)

样本集\((x_{n_1},x_{n_1+1},x_{n_1+2,\dots,x_{n_1+n_2-1}})\)就是我们需要的平稳分布对应的样本集。

上面这个过程基本上就是MCMC采样的的完整采样理论了，但是这个采样算法还是比较难 在实际中应用，为什么呢？问题在上面的3.3步，在接受率这里，由于\(\alpha(x_t,x_*)\)可能非常的小，比如\(0.1\)，导致我们大部分采样值都被拒绝转移了，采样效率芬迪，有可能我们采样了很多次的马尔可夫链还没有收敛，就是说上面这个\(n_1\)要非常非常大，从而导致拟合的效果，这一点是无法接受的，这时候就需要M-H采样出厂了。

M-H采样

M-H采样是Metropolis-Hastings采样的简称，这个算法首先由Metropolis提出，被Hastings改进，因此被称之为Metropolis-Hastings采样或M-H采样

M-H采样解决了我们上一节MCMC采样接受率过低的问题。

我们回到MCMC采样的细致平稳条件：

\[\pi(i)Q(i,j)\alpha(i,j) = \pi(j)Q(j,i)\alpha(j,i) \]

我们采样效率低的原因是\(\alpha(i,j)\)太小了，比如为0.1，而\(\alpha(j,i)\)为0.2。即：

\[\pi(i)Q(i,j)\times 0.1 = \pi(j)Q(j,i)\times 0.2 \]

这时我们可以看到，如果两边同时扩大五倍，接受率提高到了0.5，但是细致平稳条件却仍然是满足的，即：

\[\pi(i)Q(i,j)\times 0.5 = \pi(j)Q(j,i)\times 1 \]

我们发现，依据正态分布曲线关于均值\(\mu\)对称，因此对于任意两个\(x\)和\(x_*\)，均值为\(x\)，方差为1的正态分布在\(x_*\)处的概率值和均值为\(x_*\)，方差为1的正态分布在\(x\)处的概率值显然是相等的，所以得到\(Q(i,j)=Q(j,i)\)
这样我们的接受率可以做如下改进，即：

\[\alpha(i,j) = min\{ \frac{\pi(j)Q(j,i)}{\pi(i)Q(i,j)},1\} \]

\[\frac{\pi(j)Q(j,i)}{\pi(i)Q(i,j)}=\frac{\alpha(i,j)}{\alpha(j,i)}=\frac{\pi(j)}{\pi(i)} \]

通过这个微小的改造，我们就得到了可以在实际应用中使用的M-H采样算法过程如下：

对于目标采样分布\(\pi(x)\):

  1. 随机选定一个起始点\(x\)，指定燃烧期\(m\)和稳定期\(M\)。
  2. 开始采样，每一轮采样都以上一轮的采样值\(x\)为均值，方差为1，生成一个正态分布，然后在这个正态分布中依据概率随机选取一个值\(x_*\)
  3. 在\([0,1]\)的均匀分布中随机生成一个数\(U\)，并指定接受概率\(\alpha=min\{1,\frac{\pi(x_*)}{\pi(x)}\}\)，如果\(U<\alpha\)，则本轮新的采样值为\(x=x_*\)，否则本轮新的采样值仍为上一轮的\(x\)。

重复上述第二步~第三步的采样过程\(m+N\)次之后，保留后\(N\)此采样结果作为目标分布的近似采样。我们还是对之前使用过的目标分布:

\[\pi(x)=\frac{1}{1.2113}(0.3exp(-(z-0.3)^2)+0.7exp(\frac{-(z-2)^2}{0.3})) \]

，进行采样，燃烧采样个数为\(m\)，最终保留有效采样数量为\(N\)。

点击查看代码 
import random
from scipy.stats import norm
import matplotlib.pyplot as plt
import numpy as np
import seaborn
seaborn.set()
# 目标采样分布pi
def pi(x):
    return (0.3 * np.exp(-(x - 0.3) ** 2) + 0.7 * np.exp(-(x - 2.) ** 2 / 0.3)) / 1.2113
m = 1 # 燃烧期样本数
N = 100000 # 实际保留的有效样本数
sample = [0 for i in range(m + N)] # 采样数组
sample[0] = 0 # 任选一个起始点，选择默认的0也可以，效果一样
# 基于接受概率，在建议马尔科夫链上随机游走采样
for t in range(1, m + N):
    x = sample[t - 1]
    x_star = norm.rvs(loc=x, scale=1, size=1)[0] # 生成下一时刻随机游走的点x*
    alpha = min(1, (pi(x_star) / pi(x))) # 接受概率
    u = random.uniform(0, 1) # 生成满足0~1之间均匀分布的随机数
    if u < alpha: # 接受-拒绝的过程
        sample[t] = x_star
    else:
        sample[t] = x
x = np.arange(-2, 4, 0.01)
plt.plot(x, pi(x), color='r') # 实际目标分布
plt.hist(sample[m:], bins=100, density=1, color='b', edgecolor='k', alpha=0.6) # 实际
plt.show()

Gibbs采样

在上一篇中，我们讲到了细致平稳条件：如果非周期马尔科夫链的状态转移矩阵\(P\)和概率分布\(π(x)\)对于所有的\(i,j\)满足：

\[\pi(i)P(i,j) = \pi(j)P(j,i) \]

则称概率分布\(π(x)\)是状态转移矩阵\(P\)的平稳分布。
在M-H采样中我们通过引入接受率使细致平稳条件满足。现在我们换一个思路。
\(\color{red}{注意此处具体就是三维数据了}\)
从二维数据分布开始，假设\(\pi(x_1,x_2)\)是一个二维联合数据分布，观察第一个特征维度相同的两个点\(A(x_1^{(1)},x_2^{(1)})\)和\(B(x_1^{(1)},x_2^{(2)})\)，根据条件概率公式，可以发现下面两式成立：

\[\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)}) \]

也就是：

\[\pi(A) \pi(x_2^{(2)} | x_1^{(1)}) = \pi(B) \pi(x_2^{(1)} | x_1^{(1)}) \]

观察上式再观察细致平稳条件的公式，我们发现在\(x_1=x^{(1)}_1\)这条直线上，如果用条件概率分布\(\pi(x_2| x_1^{(1)})\)作为马尔科夫链的状态转移概率，则任意两个点之间的转移满足细致平稳条件！这真是一个开心的发现，同样的道理，在在\(x_2 = x_2^{(1)}\)这条直线上，如果用条件概率分布\(\pi(x_1| x_2^{(1)})\)作为马尔科夫链的状态转移概率，则任意两个点之间的转移也满足细致平稳条件。那是因为假如有一点\(C(x_1^{(2)},x_2^{(1)})\),我们可以得到：

\[\pi(A) \pi(x_1^{(2)} | x_2^{(1)}) = \pi(C) \pi(x_1^{(1)} | x_2^{(1)}) \]

基于上面的发现，我们可以这样构造分布\(\pi(x_1,x_2)\)的马尔可夫链对应的状态转移矩阵\(P\)：

\[P(A \to B) = \pi(x_2^{(B)}|x_1^{(1)})\;\; if\; x_1^{(A)} = x_1^{(B)} =x_1^{(1)} \]

\[P(A \to C) = \pi(x_1^{(C)}|x_2^{(1)})\;\; if\; x_2^{(A)} = x_2^{(C)} =x_2^{(1)} \]

\[P(A \to D) = 0\;\; else \]

有了上面这个状态转移矩阵，我们很容易验证二维平面上的任意两点E,F，满足细致平稳条件时：

\[\pi(E)P(E \to F) = \pi(F)P(F \to E) \]

于是这个二维空间上的马氏链将收敛到平稳分布 \(π(x,y)\)

二维Gibbs采样

利用上一节找到的状态转移矩阵，我们就得到了二维Gibbs采样，这个采样需要两个维度之间的条件概率。具体过程如下：

  1. 输入平稳分布\(\pi(x_1,x_2)\)，设定状态转移次数阈值\(n_1\)，需要的样本个数\(n_2\)
  2. 随机初始化初始状态值\(x_1^{(0)}\)和\(x_2^{(0)}\)
  3. \(for \quad t=0\quad to\quad n_1+n_2-1:\)
    1. 从条件概率分布\(P(x_2|x_1^{(t)})\)中采样得到样本\(x_2^{t+1}\)
    2. 从条件概率分布\(P(x_1|x_2^{(t+1)})\)中采样得到样本\(x_1^{t+1}\)

样本集\(\{(x_1^{(n_1)}, x_2^{(n_1)}), (x_1^{(n_1+1)}, x_2^{(n_1+1)}), ..., (x_1^{(n_1+n_2-1)}, x_2^{(n_1+n_2-1)})\}\)即为我们需要的平稳分布对应的样本集。

整个采样过程中，我们通过轮换坐标轴，采样的过程为：

\[(x_1^{(1)}, x_2^{(1)}) \to (x_1^{(1)}, x_2^{(2)}) \to (x_1^{(2)}, x_2^{(2)}) \to ... \to (x_1^{(n_1+n_2-1)}, x_2^{(n_1+n_2-1)}) \]

用下图可以很直观的看出，采样是在两个坐标轴上不停的轮换的。当然，坐标轴轮换不是必须的，我们也可以每次随机选择一个坐标轴进行采样。不过常用的Gibbs采样的实现都是基于坐标轴轮换的。

image

采样实例

假设我们要采样的是一个二维正态分布\(Norm(\mu,\sum)\)，其中

\[\mu=(\mu_1,\mu_2)=(5,-1) \]

\[\sum= \begin{pmatrix} \sigma_1^2 & p\sigma_1\sigma_2 \\ p\sigma_1\sigma_2 & \sigma_2^2 \end{pmatrix}=\begin{pmatrix} 1 & 1 \\ 1 & 4 \end{pmatrix}\]

而采样过程中的需要的状态转移条件分布为：

\[P(x_1|x_2) = Norm\left ( \mu _1+\rho \sigma_1/\sigma_2 \left ( x _2-\mu _2 \right ), (1-\rho ^2)\sigma_1^2 \right ) \]

\[P(x_2|x_1) = Norm\left ( \mu _2+\rho \sigma_2/\sigma_1 \left ( x _1-\mu _1 \right ), (1-\rho ^2)\sigma_2^2 \right ) \]

代码如下：

点击查看代码 
import random
import math
from scipy.stats import norm
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.stats import multivariate_normal
samplesource = multivariate_normal(mean=[5,-1], cov=[[1,1],[1,4]])

def p_ygivenx(x, m1, m2, s1, s2):
    return (random.normalvariate(m2 + rho * s2 / s1 * (x - m1), math.sqrt((1 - rho ** 2) * (s2**2))))

def p_xgiveny(y, m1, m2, s1, s2):
    return (random.normalvariate(m1 + rho * s1 / s2 * (y - m2), math.sqrt((1 - rho ** 2) * (s1**2))))

N = 5000
K = 20
x_res = []
y_res = []
z_res = []
m1 = 5
m2 = -1
s1 = 1
s2 = 2

rho = 0.5
y = m2

for i in range(N):
    for j in range(K):
        x = p_xgiveny(y, m1, m2, s1, s2)
        y = p_ygivenx(x, m1, m2, s1, s2)
        z = samplesource.pdf([x,y])
        x_res.append(x)
        y_res.append(y)
        z_res.append(z)

num_bins = 50
plt.hist(x_res, num_bins, density=1, facecolor='green', alpha=0.5)
plt.hist(y_res, num_bins, density=1, facecolor='red', alpha=0.5)
plt.title('Histogram')
plt.show()

输出的两个特征各自的分布如下：

image

然后我们看看样本集生成的二维正态分布，代码如下：

点击查看代码 
fig = plt.figure()
ax = Axes3D(fig, rect=[0, 0, 1, 1], elev=30, azim=20)
ax.scatter(x_res, y_res, z_res,marker='o')
plt.show()

image
image

Gibbs采样小结

由于Gibbs采样在高维特征时的优势，目前我们通常意义上的MCMC采样都是用的Gibbs采样。当然Gibbs采样是从M-H采样的基础上的进化而来的，同时Gibbs采样要求数据至少有两个维度，一维概率分布的采样是没法用Gibbs采样的,这时M-H采样仍然成立。

有了Gibbs采样来获取概率分布的样本集，有了蒙特卡罗方法来用样本集模拟求和，他们一起就奠定了MCMC算法在大数据时代高维数据模拟求和时的作用。MCMC系列就在这里结束吧。

参考来源：

刘建平Pinard

posted @ 2021-10-29 20:10  X-POWER  阅读(335)  评论(0)  编辑  收藏  举报