第3章采样方法
上一章介绍了采样的概念。例如，想知道一所大学里所有男生的平均身高。但是因为学校里的男生可能有上万人之多，所以为每个人都测量一下身高存在很大困难，于是从每个学院随机挑选出100名男生来作为样本，这个过程就是采样。然而，本章将要讨论的采样则有另外一层含义。现实中的很多问题可能求解起来是相当困难的。这时就可能会想到利用计算机模拟的方法来帮助求解。在使用计算机进行模拟时，所说的采样，是指从一个概率分布中生成观察值的方法。而这个分布通常是由其概率密度函数来表示的。但即使在已知概率密度函数的情况下，让计算机自动生成观测值也不是一件容易的事情。
3.1蒙特卡洛法求定积分
蒙特卡洛(Monte Carlo)法是一类随机算法的统称。它是20世纪40年代中期由于科学技术的发展，尤其是电子计算机的发明，而被提出并发扬光大的一种以概率统计理论为基础的数值计算方法。它的核心思想就是使用随机数(或更准确地说是伪随机数)来解决一些复杂的计算问题。现今，蒙特卡洛法已经在诸多领域展现出了超强的能力。本节，我们将通过蒙特卡洛法最为常见的一种应用——求解定积分，来演示这类算法的核心思想。

3.1.1无意识统计学家法则

作为一个预备知识，先来介绍一下无意识统计学家法则(Law of the Unconscious Statistician，LOTUS)。在概率论与统计学中，如果知道随机变量X的概率分布，但是并不显式地知道函数g(X)的分布，那么LOTUS就是一个可以用来计算关于随机变量X的函数g(X)之期望的定理。该法则的具体形式依赖于随机变量X之概率分布的描述形式。
如果随机变量X的分布是离散的，而且我们知道它的PMF是fX，但不知道fg(X)，那么g(X)的期望是
E［g(X)］=∑xg(x)fX(x)
其中和式是在取遍X的所有可能之值x后求得。
如果随机变量X的分布是连续的，而且我们知道它的PDF是fX，但不知道fg(X)，那么g(X)的期望是
E［g(X)］=∫∞-∞g(x)fX(x)

简而言之，已知随机变量X的概率分布，但不知道g(X)的分布，此时用LOTUS公式能计算出函数g(X)的数学期望。其实就是在计算期望时，用已知的X的PDF(或PMF)代替未知的g(X)的PDF(或PMF)。

3.1.2投点法


图31投点法求定积分


投点法是讲解蒙特卡洛法基本思想的一个最基础也最直观的实例。这个方法也常常被用来求圆周率π。现在我们用它来求函数的定积分。如图31所示，有一个函数f(x)，若要求它从a到b的定积分，其实就是求曲线下方的面积。
可以用一个比较容易算得面积的矩型罩在函数的积分区间上(假设其面积为Area)。然后随机地向这个矩形框里面投点，其中落在函数f(x)下方的点为菱形，其他点为三角形。然后统计菱形点的数量占所有点(菱形+三角形)数量的比例为r，那么就可以据此估算出函数f(x)从a到b的定积分为Area×r。
注意由蒙特卡洛法得出的值并不是一个精确值，而是一个近似值。而且当投点的数量越来越大时，这个近似值也越接近真实值。

3.1.3期望法
下面来重点介绍利用蒙特卡洛法求定积分的第二种方法——期望法，有时也称为平均值法。
任取一组相互独立、同分布的随机变量{Xi}，Xi在［a，b］上服从分布律fX，也就是说fX是随机变量X的PDF(或PMF)。令g*(x)=g(x)fX(x)，则g*(Xi)也是一组独立同分布的随机变量，而且因为g*(x)是关于x的函数，所以根据LOTUS可得
E［g*(Xi)］=∫bag*(x)fX(x)dx=∫bag(x)dx=I
由强大数定理
PrlimN→∞1N∑Ni=1g*(Xi)=I=1
若选
=1N∑Ni=1g*(Xi)
则依概率1收敛到I。平均值法就用作为I的近似值。
假设要计算的积分有如下形式
I=∫bag(x)dx
其中，被积函数g(x)在区间［a，b］上可积。任意选择一个有简便办法可以进行抽样的概率密度函数fX(x)，使其满足下列条件： 
(1) 当g(x)≠0时，fX(x)≠0，a≤x≤b； 
(2) ∫bafX(x)dx=1。
如果记
g*(x)=g(x)fX(x)
，fX(x)≠0

0，fX(x)=0

那么原积分式可以写成
I=∫bag*(x)fX(x)dx
因而求积分的步骤是： 
(1) 产生服从分布律fX的随机变量Xi，i=1，2，…，N； 
(2) 计算均值
=1N∑Ni=1g*(Xi)
并用它作为I的近似值，即I≈。
如果a，b为有限值，那么fX可取作为均匀分布
fX(x)=1b-a，a≤x≤b

0，其他
此时原来的积分式变为
I=(b-a)∫bag(x)1b-adx
因而求积分的步骤是： 
(1) 产生［a，b］上的均匀分布随机变量Xi，i=1，2，…，N； 
(2) 计算均值
=b-aN
∑Ni=1g(Xi)

并用它作为I的近似值，即I≈。

最后来看一下平均值法的直观解释。注意积分的几何意义就是［a，b］区间曲线下方的面积，如图32所示。

当在［a，b］随机取一点x时，它对应的函数值就是f(x)，然后便可以用f(x)·(b-a)来粗略估计曲线下方的面积(也就是积分)，如图33所示，当然这种估计(或近似)是非常粗略的。


图32积分的几何意义




图33对积分值进行粗略估计


于是我们想到在［a，b］随机取一系列点xi时(xi满足均匀分布)，然后把估算出来的面积取平均来作为积分估计的一个更好的近似值，如图34所示。可以想象，如果这样的采样点越来越多，那么对于这个积分的估计也就越来越接近。


图34对积分值进行估计


按照上面这个思路，得到积分公式为
=(b-a)1N
∑N-1i=0f(Xi)=1N
∑N-1i=0f(Xi)1b-a

其中，1b-a 就是均匀分布的PMF。这跟之前推导出来的蒙特卡洛积分公式是一致的。

3.2蒙特卡洛采样
在3.1节中，通过求解定积分这个具体的例子来演示了蒙特卡洛方法的基本思想，而其中的核心就是使用随机数。当所求解问题可以转化为某种随机分布的特征数(如随机事件出现的概率，或者随机变量的期望值等)时，往往就可以考虑使用蒙特卡洛方法。通过随机抽样的方法，以随机事件出现的频率估计其概率，或者以抽样的数字特征估算随机变量的数字特征，并将其作为问题的解。这种方法多用于求解复杂的高维积分问题。
实际应用中，所要面对的第一个问题就是如何抽样？注意，在计算机模拟时，这里所说的抽样其实是指从一个概率分布中生成观察值(observations)的方法。而这个分布通常是由其概率密度函数来表示的。前面曾经提过，即使在已知PDF的情况下，让计算机自动生成观测值也不是一件容易的事情。从本质上来说，计算机只能实现对均匀分布的采样。幸运的是，仍然可以在此基础上对更为复杂的分布进行采样。

3.2.1逆采样
比较简单的一种情况是，可以通过PDF与CDF之间的关系，求出相应的CDF。或者根本就不知道PDF，但是知道CDF。此时就可以使用CDF反函数(及分位数函数)的方法来进行采样。这种方法又称为逆变换采样(Inverse Transform Sampling)。
假设已经得到了CDF的反函数F-1(u)，如果想得到m个观察值，则重复下面的步骤m次： 
(1) 从U(0，1)中随机生成一个值(计算机可以实现从均匀分布中采样)，用u表示； 
(2) 计算F-1(u)的值x，则x就是从目标分布f(x)中得出的一个采样点。
面对一个具有复杂表达式的函数，逆变换采样法真有效吗？来看一个例子，假设现在希望从具有下面这个PDF的分布中采样
f(x)=8x，0≤x<0.25

83-8x3，0.25≤x≤1

0，其他
可以算得相应的CDF为
F(x)=0，x<0

4x2，0≤x<0.25

8x3-4x23-13，0.25≤x≤1

1，x>1
对于u∈［0，1］，上述CDF的反函数为
F-1(u)=u2，0≤x<0.25

1-3(1-u)2，0.25≤x≤1
下面在R中利用上面的方法采样10000个点，并以此来演示抽样的效果。



m <- 10000

u <- runif(m， 0， 1)

invcdf.func <- function(u) {

+ if (u >= 0 && u < 0.25)

+sqrt(u)/2

+else if (u >= 0.25 && u <= 1)

+1 - sqrt(3 * (1 - u))/2

+ }

x <- unlist(lapply(u， invcdf.func))



curve(8*x， from = 0， to = 0.25， xlim=c(0，1)， ylim=c(0，2)， 

+ col = "red"，xlab = ""， ylab="")

par(new=TRUE)

curve((8/3)-(8/3)*x， from = 0.25， to = 1， xlim=c(0，1)， ylim=c(0，2)， 

+ col = "red"，xlab = ""， ylab="")

par(new=TRUE)

plot(density(x)， xlim=c(0，1)， ylim=c(0，2)， 

+ col = "blue"， xlab = "x"， ylab="density")






图35逆变换采样举例


将所得结果与真实的PDF函数图形进行对照，如图35所示。可见由逆变换采样法得到的点所呈现处理的分布与目标分布非常吻合。

下面再举一个稍微复杂一点的例子，已知分布的PDF如下
h(x)=2m2(1-m2)x3
，x∈［m，1］

可以算得相应的CDF为
H(x)=∫x-∞h(t)dt0，x<m

11-m2
-m2(1-m2)x2
，x∈［m，1］

1，x>1

对于u∈［0，1］，它的反函数为
H-1(u)=m21-(1-m2)u


同样，给出R中的示例代码如下。



invcdf <- function(u， m) {

return(sqrt(m^2/(1 - (1 - m^2) * u)))

}



sample1 <- sapply(runif(1000)， invcdf， m = .5)




下面这段代码利用R中提供的一些内置函数实现了已知PDF时基于逆变换方法的采样，将新定义的函数命名为samplepdf()。当然，对于那些过于复杂的PDF函数(例如很难积分的)，samplepdf()确实有力所不及的情况。但是对于标准的常规PDF，该函数的效果还是不错的。



endsign <- function(f， sign = 1) {

b <- sign

while (sign * f(b) < 0) b <- 10 * b

return(b)

}



samplepdf <- function(n， pdf， ...， spdf.lower = -Inf， spdf.upper = Inf) {

vpdf <- function(v) sapply(v， pdf， ...)  # vectorize

cdf <- function(x) integrate(vpdf， spdf.lower， x)$value

invcdf <- function(u) {

subcdf <- function(t) cdf(t) - u

if (spdf.lower == -Inf) 

spdf.lower <- endsign(subcdf， -1)

if (spdf.upper == Inf) 

spdf.upper <- endsign(subcdf)

return(uniroot(subcdf， c(spdf.lower， spdf.upper))$root)

}

sapply(runif(n)， invcdf)

}




下面就用samplepdf()函数来对上面给定的h(x)进行采样，然后再跟之前所得结果进行对比。



h <- function(t， m) {

if (t >= m & t <= 1) 

return(2 * m^2/(1 - m^2)/t^3)

return(0)







}



sample2 <- samplepdf(1000， h， m = .5)



plot(density(sample1)， xlim=c(0.4， 1.1)， ylim=c(0， 4)， 

+ col = "red"， xlab = ""， ylab=""， main="")

par(new=TRUE)

plot(density(sample2)， xlim=c(0.4， 1.1)， ylim=c(0， 4)， 

+ col ="blue"， xlab = "x， N=1000"， ylab="density"， main="")

text.legend = c("my_invcdf"，"samplepdf")

legend("topright"， legend = text.legend， lty = c(1，1)， col = c( "red"， "blue"))




代码执行结果如图36所示。


图36逆采样代码执行结果



3.2.2博克斯穆勒变换
博克斯穆勒变换(BoxMuller Transform)最初由乔治·博克斯(George Box)与默文·穆勒(Mervin Muller)在1958年共同提出。博克斯是统计学的一代大师，统计学中的很多名词术语都以其名字命名。博克斯与统计学的家学渊源相当深厚，他的导师是统计学开山鼻祖皮尔逊的儿子，英国统计学家埃贡·皮尔逊(Egon Pearson)，博克斯还是统计学的另外一位巨擘级奠基人费希尔的女婿。统计学中的名言“所有模型都是错的，但其中一些是有用的”也出自博克斯之口。
本质上来说，计算机只能生产符合均匀分布的采样。如果要生成其他分布的采样，就需要借助一些技巧性的方法。而在众多的“其他分布”中，正态分布无疑占据着相当重要的地位。下面这个定理，就为生成符合正态分布的采样(随机数)提供了一种方法，而且这也是很多软件或者编程语言的库函数中生成正态分布随机数时所采样的方法。
定理(BoxMuller变换)： 如果随机变量U1和U2是独立同分布的，且U1，U2~U［0，1］，则
Z0=-2lnU1cos(2πU2)

Z1=-2lnU1sin(2πU2)

其中，Z0和Z1独立且服从标准正态分布。
如何来证明这个定理呢？这需要用到一些微积分中的知识，首先回忆一下二重积分化为极坐标下累次积分的方法
Df(x，y)dxdy=∫βαdθ∫ρ2(θ)ρ1(θ)f(ρcosθ，ρsinθ)ρdρ
假设现在有两个独立的标准正态分布X~N(0，1)和Y~N(0，1)，由于两者相互独立，则联合概率密度函数为
p(x，y)=p(x)·p(y)=12π
e-x22
·12π
e-y22
=12π
e-x2+y22

做极坐标变换，则x=Rcosθ，y=Rsinθ，则有
12πe-x2+y22=12πe-R22
这个结果可以看成是两个概率分布的密度函数的乘积，其中一个可以看成是［0，2π］上均匀分布，将其转换为标准均匀分布则有θ~U(0，2π)=2πU2。
另外一个的密度函数为
P(R)=e-R22

则其累计分布函数CDF为
P(R≤r)=∫r0e-ρ22ρdρ=-e-ρ22r0=-e-r22+1
这个CDF函数的反函数可以写成
F-1(u)=-2log(1-u)
根据逆变换采样的原理，如果有个PDF为P(R)的分布，那么对齐CDF的反函数进行均匀采样所得的样本分布将符合P(R)的分布，而如果u是均匀分布的，那么U1=1-u也将是均匀分布的，于是用U1替换1-u，最后可得
X=R·cosθ=-2logU1·cos(2πU2)

Y=R·sinθ=-2logU1·sin(2πU2)
结论得证。最后来总结一下利用BoxMuller变换生成符合高斯分布的随机数的方法： 
(1) 产生两个随机数U1，U2~U［0，1］； 
(2) 用它们来创造半径R=-2log(U1)和和夹角θ=2πU2； 
(3) 将(R，θ)从极坐标转换到笛卡儿坐标： (Rcosθ，Rsinθ)。
3.2.3拒绝采样与自适应拒绝采样
读者已经看到逆变换采样的方法确实有效。但其实它的缺点也是很明显的，那就是有些分布的CDF可能很难通过对PDF的积分得到，再或者CDF的反函数也很不容易求。这时可能需要用到另外一种采样方法，这就是下面即将要介绍的拒绝采样(Reject Sampling)。


图37拒绝采样的原理


图37阐释了拒绝采样的基本思想。假设对PDF为p(x)的函数进行采样，但是由于种种原因(例如这个函数很复杂)，对其进行采样是相对困难的。但是另外一个PDF为q(x)的函数则相对容易采样，例如采用逆变换方法可以很容易对对它进行采样，甚至q(x)就是一个均匀分布(别忘了计算机可以直接进行采样的分布就只有均匀分布)。那么，当我们将q(x)与一个常数M相乘之后，可以实现如图37所示关系，即M·q(x)将p(x)完全“罩住”。

然后重复如下步骤，直到获得m个被接受的采样点： 
(1) 从q(x)中获得一个随机采样点xi； 
(2) 对于xi计算接受概率(acceptanceprobability)
α=p(xi)Mq(xi)
(3) 从U(0，1)中随机生成一个值，用u表示； 
(4) 如果α≥u，则接受xi作为一个来自p(x)的采样值，否则就拒绝xi并回到第一步。

当然可以采用严密的数学推导来证明拒绝采样的可行性。但它的原理从直观上来解释也是相当容易理解的。可以想象一下在图27的例子中，从哪些位置抽出的点会比较容易


图38拒绝采样示例


被接受。显然，红色曲线(位于上方)和绿色曲线(位于下方)所示之函数更加接近的地方接受概率较高，即是更容易被接受，所以在这样的地方采到的点就会比较多，而在接受概率较低(即两个函数差距较大)的地方采到的点会比较少，这也就保证了这个方法的有效性。

还是以本章前面给出的那个分段函数f(x)为例来演示拒绝采样方法。如图38所示，所选择的参考分布是均匀分布(当然也可以选择其他的分布，但采用均匀分布显然是此处最简单的一种处理方式)。而且令常数M=3。

下面给出R中的示例代码，易见此处采样点数目为10000。



f.x <- function(x) {

+ if (x >= 0 && x < 0.25)

+ 8 * x

+ else if (x >= 0.25 && x <= 1)

+ 8/3 - 8 * x/3

+ else 0

+ }



g.x <- function(x) {

+ if (x >= 0 && x <= 1)

+ 1

+ else 0

+ }



M <- 3

m <- 10000

n.draws <- 0

draws <- c()

x.grid <- seq(0， 1， by = 0.01)







while (n.draws < m) {

+ x.c <- runif(1， 0， 1)

+ accept.prob <- f.x(x.c)/(M * g.x(x.c))

+ u <- runif(1， 0， 1)

+ if (accept.prob >= u) {

+draws <- c(draws， x.c)

+n.draws <- n.draws + 1

+}

+ }




上述代码的执行结果如图39所示，可见采样结果是非常理想的。
下面的例子演示了对(表达式非常复杂的)beta(3，6)分布进行拒绝采样的效果。这里采用均匀分布作为参考分布。而且这里的Mq(x)所取之值就是beta(3，6)分布的极大值，它的函数图形应该是与beta(3，6)的极值点相切的一条水平直线。



sampled <- data.frame(proposal = runif(50000，0，1))

sampled$targetDensity <- dbeta(sampled$proposal， 3，6)



maxDens = max(sampled$targetDensity， na.rm = T)

sampled$accepted = ifelse(runif(50000，0，1) 

+ < sampled$targetDensity / maxDens， TRUE， FALSE)



hist(sampled$proposal［sampled$accepted］， freq = F， 

+ col = "grey"， breaks = 100， main="")

curve(dbeta(x， 3，6)，0，1， add =T， col = "red"， main="")




图310给出了采样50000个点后的密度分布情况，可见采样分布与目标分布beta(3，6)非常吻合。


图39程序执行结果




图310拒绝采样举例


拒绝采样的方法确实可以解决我们的问题。但是它的一个不足涉及其采样效率的问题。针对上面给出的例子而言，我们选择了离目标函数最近的参考函数，就均匀分布而言，已经不能有更进一步的方法了。但即使这种，在这个类似钟形的图形两侧其实仍然会拒绝掉很多很多采样点，这种开销相当浪费。最理想的情况下，参考分布应该跟目标分布越接近越好，从图形上来看就是包裹的越紧实越好。但是这种情况的参考分布往往又不那么容易得到。在满足某些条件的时候也确实可以采用所谓的改进方法，即自适应的拒绝采样(Adaptive Rejection Sampling)。


拒绝采样的弱点在于当被拒绝的点很多时，采样的效率会非常不理想。同时我们也知道，如果能够找到一个跟目标分布函数非常接近的参考函数，那么就可以保证被接受的点占大多数(被拒绝的点很少)。这样一来便克服了拒绝采样效率不高的弱点。如果函数是logconcave的话，那么就可以采用自适应的拒绝采样方法。什么是logconcave呢？还是回到之前介绍过的贝塔分布，用下面的代码来绘制beta(2，3)的概率密度函数图像，以及将beta(2，3)的函数取对数之后的图形。



> integrand <- function(x) {(x^1)*((1-x)^2)}

> integrate(integrand， lower = 0， upper = 1)

0.08333333 with absolute error < 9.3e-16

> f<-function(x，a，b){log(1/0.08333)+(a-1)*log(x)+(b-1)*log(1-x)}

> curve(f(x， 2， 3))

> curve(dbeta(x， 2， 3))




上述代码的执行结果如图311所示，其中(a)是beta(2，3)的概率密度函数图形，(b)是将beta(2，3)的函数取对数之后的图形，你可以发现结果是一个凹函数(concave)。那么beta(2，3)就满足logconcave的要求。


图311贝塔函数与其取对数后的函数图形


然后在对数图像上找一些点做图像的切线，如图312所示。因为取对数后的函数是凹函数，所以每个切线都相当于一个超平面，而且对数图像只会位于超平面的一侧。

同时给出用以绘制图312的代码，要知道R语言的一个强项就是绘图。



log_f <- function(x，a，b){log(1/0.08333)+(a-1)*log(x)+(b-1)*log(1-x)}

g <- function(x，a，b){(a-1)/x-(b-1)/(1-x)}



log_f_y1 <- log_f(0.18， 2， 3)

log_f_y2 <- log_f(0.40， 2， 3)

log_f_y3 <- log_f(0.65， 2， 3)

log_f_y4 <- log_f(0.95， 2， 3)



g1 <- g(0.18， 2， 3)







b1 <- log_f_y1 - g1*0.18

y1 <- function(x) {g1*x+b1}



g2 <- g(0.40， 2， 3)

b2 <- log_f_y2 - g2*0.40

y2 <- function(x) {g2*x+b2}



g3 <- g(0.65， 2， 3)

b3 <- log_f_y3 - g3*0.65

y3 <- function(x) {g3*x+b3}



g4 <- g(0.95， 2， 3)

b4 <- log_f_y4 - g4*0.95

y4 <- function(x) {g4*x+b4}



curve(log_f(x， 2， 3)， col = "blue"， xlim = c(0，1)， ylim = c(-7， 1))

curve(y1， add = T， lty = 2， col = "red"， to = 0.38)

curve(y2， add = T， lty = 2， col = "red"， from = 0.15， to=0.78)

curve(y3， add = T， lty = 2， col = "red"， from = 0.42)

curve(y4， add = T， lty = 2， col = "red"， from = 0.86)



par(new=TRUE)



xs = c(0.18， 0.40， 0.65， 0.95)

ys = c(log_f_y1， log_f_y2， log_f_y3， log_f_y4)

plot(xs， ys， col="green"， xlim=c(0，1)， ylim=c(-7， 1)， xlab=""， ylab="")




再把这些切线转换回原始的beta(2，3)图像中，显然原来的线性函数会变成指数函数，它们将对应图313中的一些曲线，这些曲线会被原函数的图形紧紧包裹住。特别是当这些的指数函数变得很多很稠密时，以彼此的交点作为分界线，其实相当于得到了一个分段函数。这个分段函数是原函数的一个逼近。用这个分段函数来作为参考函数再执行拒绝采样，自然就完美地解决了之前的问题。


图312做取对数后的图形的切线





图313切线取指数函数后的变换结果


下面是用来画出图313的R语言代码。



e_y1 <- function(x) {exp(g1*x+b1)}

e_y2 <- function(x) {exp(g2*x+b2)}

e_y3 <- function(x) {exp(g3*x+b3)}

e_y4 <- function(x) {exp(g4*x+b4)}



curve(dbeta(x， 2， 3)， col="blue"， ylim=c(0， 2.0))

curve(e_y1(x)， add=T， lty=3， col = "red")

curve(e_y2(x)， add=T， lty=3， col = "red"， from = 0.2， to = 0.75)

curve(e_y3(x)， add=T， lty=3， col = "red"， from=0.48)

curve(e_y4(x)， add=T， lty=3， col = "red"， from=0.86)




这无疑是一种绝妙的想法。而且这种想法，在前面其实已经暗示过。在上一部分最后一个例子中，我们其实就是选择了一个与原函数相切的均匀分布函数来作为参考函数。我们当然会想去选择更多与原函数相切的函数，然后用这个函数的集合来作为新的参考函数。只是由于原函数的凹凸性无法保证，所以直线并不是一种好的选择。而自适应拒绝采样(Adaptive Rejection Sampling，ARS)所采用的策略则非常巧妙地解决了我们的问题。当然函数是logconcave的条件必须满足，否则就不能使用ARS。
下面给出一个在R中进行自适应拒绝采样的例子。显然，该例子要比之前的代码简单许多。因为R中ars包已经提供了一个现成的用于执行自适应拒绝采样的函数，即ars()。关于这个函数在用法上的一些细节，读者还可以进一步参阅R的帮助文档，这里不再赘言。此次我们需要指出： ars()函数中两个重要参数，一个是对原分布的PDF取对数，另外一个则是对PDF的对数形式再进行求导(在求导时我们忽略了前面的系数项)，其实也就是为了确定切线。



f<-function(x，a，b){(a-1)*log(x)+(b-1)*log(1-x)}

fprima<-function(x，a，b){(a-1)/x-(b-1)/(1-x)}

mysample<-ars(20000， f， fprima， x=c(0.3，0.6)， m=2， lb=TRUE， xlb=0， 

+ ub=TRUE， xub=1， a=1.3， b=2.7)

hist(mysample， freq=F)

curve(dbeta(x，1.3，2.7)， add=T， col="red")




上述代码的执行结果如图314所示。


图314自适应采样代码执行结果



3.3矩阵的极限与马尔科夫链
先来看一个例子。社会学家经常把人按其经济状况分成3类： 下层(lowerclass)、中层(middleclass)、上层(upperclass)，用1、2、3分别代表这三个阶层。社会学家们发现决定一个人的收入阶层的最重要的因素就是其父母的收入阶层。如果一个人的收入属于下层类别，那么他的孩子属于下层收入的概率是0.65，属于中层收入的概率是0.28，属于上层收入的概率是0.07。从父代到子代，收入阶层的变化的转移概率如图315所示。


图315阶层变化的转移概率


使用矩阵的表示方式，转移概率矩阵记为
P=0.650.280.07

0.150.670.18

0.120.360.52
假设当前这一代人处在下、中、上层的人的比例是概率分布向量π0=［π0(1)，π0(2)，π0(3)］，那么他们的子女的分布比例将是π1=π0P，孙子代的分布比例将是π2=π1P=π1P2=…，第n代子孙的收入分布比例将是πn=πn-1P=…=π0Pn。假设初始概率分布为π0=［0.21，0.68，0.11］，则可以计算前n代人的分布状况如下： 



第n代人下层中层上层
00.2100.6800.110
10.2520.5540.194
20.2700.5120.218
30.2780.4970.225
40.2820.4900.226
50.2850.4890.225
60.2860.4890.225
70.2860.4890.225
80.2890.4880.225
90.2860.4890.225
100.2860.4890.225
…………

发现从第7代人开始，这个分布就稳定不变了，这个是偶然的吗？我们换一个初始概率分布π0=［0.75，0.15，0.1］试试看，继续计算前n代人的分布状况如下



第n代人下层中层上层
00.750.150.1
10.5220.3470.132
20.4070.4260.167
30.3490.4590.192
40.3180.4750.207
50.3030.4820.215
60.2950.4850.220
70.2910.4870.222
80.2890.4880.225
90.2860.4890.225
100.2860.4890.225
…………

结果，到第9代人的时候，分布又收敛了。最为奇特的是，两次给定不同的初始概率分布，最终都收敛到概率分布π=［0.286，0.489，0.225］，也就是说收敛的行为和初始概率分布π0无关。这说明这个收敛行为主要是由概率转移矩阵P决定的。计算一下Pn，

P20=P21=…=P100=…=0.2860.4890.225

0.2860.4890.225

0.2860.4890.225

我们发现，当n足够大的时候，这个Pn矩阵的每一行都是稳定地收敛到π=［0.286，0.489，0.225］这个概率分布。自然地，这个收敛现象并非是这个例子所独有，而是绝大多数“马尔科夫链”的共同行为。
为了求得一个理论上的结果，再来看一个更小规模的例子(这将方便后续的计算演示)，假设在一个区域内，人们要么是住在城市，要么是住在乡村。下面的矩阵表示人口迁移的一些规律(或倾向)。例如，第1行第1列就表示，当前住在城市的人口，明年将会有90%的人选择继续住在城市。

当前住在

城市
当前住在

乡村

下一年住在城市

下一年住在乡村0.900.02

0.100.98=A



图316两年后人口的转移


作为一个简单的开始，来计算一下现今住在城市的人中2年后会住在乡村的概率有多大。分析可知，当前住在城市的人，1年后，会有90%继续选择住在城市，另外10%的人则会选择搬去乡村居住。然后又过了1年，上一年中选择留在城市的人中又会有10%的人搬去乡村。而上一年中搬到乡村的人中将会有98%的选择留在乡村。这个分析过程如图316所示，最终可以算出现今住在城市的人中2年后会住在乡村的概率=0.90×0.10+0.10×0.98。

其实上述计算过程就是在做矩阵的平方。在下面给出的矩阵乘法中，你会发现结果矩阵中第2行第1列的计算就是在执行上面的操作。在此基础，我们还可以继续计算n年后的情况，也就是计算矩阵A自乘n次后的结果。
A2=0.900.02
0.100.980.900.02
0.100.98=(A2)11(A2)12

(A2)21(A2)22
如果假设最开始的时候，城乡人口的比例为7∶3，可以用一个列向量来表示它，即P=［0.7，0.3］T，想知道最终城乡人口的比例为何，则就是要计算。如果最初城乡人口的比例为9∶1，结果又如何呢？这些都要借助矩阵的极限，对角化操作以及马尔科夫链等概念来辅助计算。
来辨析三个概念： 随机过程、马尔科夫过程、马尔科夫链。这三个概念中，都涉及对象和它们对应的状态这两个要素。在刚刚给出的例子中，研究的对象就是人，人的状态分为住在城市或者住在乡村两种。
三者之中，最宽泛的概念就是随机过程，限制最多的就是马尔科夫链。对于马尔科夫链，必须满足两个条件： ①当前状态仅跟上一个状态有关； ②总共的状态数是有限的。如果状态数可以是无限多个，那样的问题就称为马尔科夫过程。但在马尔科夫过程中仍然要求，时间是离散的，例如前面的例子是以“年”为单位的。如果时间允许是连续的，那样的问题就称为随机过程。本书仅仅讨论马尔科夫链。
在某个时间点上，对象的状态为s1，下一个时刻，它的状态以某种概率转换到其他状态(也包含原状态s1)，这里所说的“以某种概率转换”最终是通过状态转移矩阵(或称随机矩阵)的形式来给出的。转移矩阵的定义如下： 
令A∈Mn×n(R)，假设： 
（1） Aij≥0； 
（2） 对于所有的1≤j≤n，∑ni=1Aij=1。

那么，A称为一个转移矩阵(或随机矩阵)。矩阵A的列向量被称为概率向量。
此外，如果矩阵的某个次幂仅包含正数项，该转移矩阵称为是正则的。这里，“某个”的意思就是存在一个整数n使得对于所有的i，j，有(An)ij>0。
从状态转移矩阵中，(结合之前的例子)可以看出Aij元素给出的信息就是(在一个单位时间间隔内)对象从状态j转移到状态i的概率。令P=(p0，p1，…，pn)T是一个向量，如果对于所有的i，有pi≥0以及∑pi=1，那么P就称为一个概率向量(probability vector)。所以可以看出，任意一个转移矩阵中的某一列都是一个概率向量。
定理： 令A是一个n×n的正则转移矩阵。那么： 
(1) 1一定是矩阵A的一个特征值，并且1的几何重数等于1，除此之外，所有其他的特征值之绝对值都小于1； 
(2) A可以对角化，并且limm→∞Am存在； 
(3) L=limm→∞Am是一个转移矩阵； 
(4) AL=LA=L； 
(5) 矩阵L=［�瘙經，�瘙經，…，�瘙經］，即L的每一列都一样，都是�瘙經。而且�瘙經就是矩阵A相对于特征值1的特征向量； 
(6) 对于任意概率向量w，都有limm→∞(Amw)=�瘙經。
这个定理非常非常奇妙的地方就是它解答了之前那个令人困扰的问题！原问题是： 如果假设最开始的时候，城乡人口的比例为7∶3，可以用一个列向量来表示它P=［0.7，0.3］T，欲知道最终城乡人口的比例为何，则就是要计算AnP，如果最开始的城乡人口的比例为9∶1，结果又如何。上述定理中的最后一条就表明，当n趋近于无穷大的时候，AnP就等于�瘙經，而且与P是无关的。更精妙的地方还在于，这个定理还告诉我们v是一个概率向量，而且它就是特征值1所对应的特征向量。
这个定理的证明已经超出了本书的范围，但可以用之前给出的例子来验证一下它。注意到如果想要计算AnP，其实就是要先设法计算矩阵A自乘n次的结果，这时为了计算方便应该先将矩阵A对角化。为此，先求矩阵A的特征多项式，通过其特征多项式便知道矩阵A有两个特征值，一个是1，一个是0.88。根据定理，1必然是该矩阵的特征值。更进一步，特征值1对应的特征向量是［1，5］T，特征值0.88对应的特征向量是［1，-1］T。所以知道矩阵对角化时所用的Q和Q-1分别为

Q=11
5-1，Q-1=-16-1-1
-51

于是可知矩阵A的对角化结果如下： 
10
00.88=D=Q-1AQ
所以有
A=QDQ-1Am=QDmQ-1

limm→∞Am=limm→∞QDmQ-1=Qlimm→∞DmQ-1
然后把值带进去就能算出最终结果如下
limm→∞Am=1/61/6

5/65/6

之前计算过特征值1对应的一个特征向量是［1，5］T，特征向量乘以一个系数仍然是特征向量(注意要求最后的特征向量同时是一个概率向量，所以会得到［1/6，5/6］T，可见上述计算与定理所揭示的结果是完全一致的。
矩阵的极限，其实就是在讨论limm→∞Am的存在性，也就是把矩阵limm→∞Am自乘m次后，如果结果矩阵中每个元素的极限都存在，就说这个矩阵的极限是存在的。而矩阵极限是否存在可以由下面的定理保证。
定理矩阵极限存在的充要条件： 
(1) |λ|<1或者λ=1，其中λ是A的任意特征值。
(2) 如果λ=1，那么它对应的几何重数等于代数重数。
定理矩阵极限存在的充分条件： 
(1) |λ|<1或者λ=1，其中λ是A的任意特征值。
(2) 矩阵A是可对角化的。

3.4查普曼柯尔莫哥洛夫等式
还是先从一个例子来谈起。图315中左侧是状态转移矩阵，其中每个位置都表示从一个状态转移到另外一个状态的概率。例如，从状态1转移到状态2的概率是0.28，所以在第一行第二列的位置，给出的数值是0.28。
在马尔科夫链中，随机变量在一个按时间排序的数组T1，T2，…，Tn中，根据状态转移矩阵让我们可以非常直观地得知当前时刻某一状态在下一时刻变到任意状态的概率，如图317所示。
现在的问题是能不能做更进一步的预测。例如，当前时刻T1时的状态为a，能否知道在下下时刻T3时状态为b的概率是多少呢？这其实也相当容易做到，如图318所示。


图317按时间排序的状态转移情况





图318计算下下时刻某状态的概率


这个过程用转移矩阵的自乘来表示是非常直观且方便的。矩阵P中第1行就表示，当前时刻状态a在下一时刻变到状态a、b、c的概率，矩阵P中第2列就表示，由状态a、b、c，在下一时刻转移到状态b的概率。那么根据矩阵的乘法公式，P2中第1行第2列的位置就表示“当前时刻T1时的状态为a，在下下时刻T3时状态为b的概率”。也就是说，如果想得到跨越2个时刻的转移矩阵，就把跨越1个时刻的转移矩阵乘以跨越1个时刻的转移矩阵即可。同理，如果想跨越3个时刻的转移矩阵，就把跨越2个时刻的转移矩阵乘以跨越1个时刻的转移矩阵即可。更普遍地有Pm+n=PmPn，而这个关系就被称为查普曼柯尔莫哥洛夫等式(ChapmanKolmogorov Equation)。
下面是查普曼柯尔莫哥洛夫等式的一个表述： 令T是一个离散状态空间中的n步马尔科夫链，其状态转移矩阵为
Pn=［pn(j，k)］j，k∈S
其中，pn(j，k)=P(Tm+n=k|Tm=j)=pnj，k是n步状态转移概率。那么，Pm+n=PmPn，或等价地有
pn+mi，j=∑k∈Spni，kpmk，j
最后给出查普曼柯尔莫哥洛夫等式的证明。
首先考虑在n时刻，系统正处在什么状态。给定T0=i，pni，k=P(Tn=k|T0=i)就表示(在时刻n)系统处于k状态的概率。但是，此后再给定Tn=k，在第n时刻之后的情况就与过去的历史彼此独立了。于是，再过m个时间单位后，处在状态j的概率是pmk，j将满足乘积
pni，kpmk，j=P(Tn+m=j，Tn=k|T0=i)
将所有k的情况进行加总就会得到要证明之结论。一个更加严格的证明如下 
pn+mi，j=P(Tn+m=j|T0=i)

=∑k∈SP(Tn+m=j，Tn=k|T0=i)

=∑k∈SP(Tn+m=j，Tn=k，T0=i)P(T0=i)

=∑k∈SP(Tn+m=j|Tn=k，T0=i)P(Tn=k，T0=i)P(T0=i)

=∑k∈Spmk，jP(Tn=k，T0=i)P(T0=i)

=∑k∈Spni，kpmk，j
上述证明中运用了马尔科夫链的性质来得到
P(Tn+m=j|Tn=k，T0=i)=P(Tn+m=j|Tn=k)=P(Tm=j|T0=k)=pmk，j
最终，结论得证。

3.5马尔科夫链蒙特卡洛
在以贝叶斯方法为基础的各种机器学习技术中，常常需要计算后验概率，再通过最大后验概率(MAP)的方法进行参数推断和决策。然而，在很多时候，后验分布的形式可能非常复杂，这个时候寻找其中的最大后验估计或者对后验概率进行积分等计算往往非常困难，此时可以通过采样的方法来求解。这其中非常重要的一个基础就是马尔科夫链蒙特卡洛(Markov chain Monte Carlo，MCMC)。

3.5.1重要性采样
前面已经详细介绍了基于随机算法进行定积分求解的技术。这里主要用到其中的平均值法。这里仅做简单回顾。
在计算定积分∫bag(x)dx时，会把g(x)拆解成两项的乘积，即g(x)=f(x)p(x)，其f(x)是某种形式的函数，而p(x)是关于随机变量X的概率分布(也就是PDF或PMF)。如此一来，上述定积分就可以表示为求f(x)的期望，即
∫bag(x)dx=∫baf(x)p(x)dx=E［f(x)］
当然，g(x)的分布函数可能具有很复杂的形式，仍然无法直接求解，这时就可以用采样的方法去近似。这时积分可以表示为
∫bag(x)dx=E［f(x)］=1n∑ni=1f(xi)
在贝叶斯分析中，蒙特卡洛积分运算常常被用来在对后验概率进行积分时做近似估计。比如要计算
I(y)=∫baf(y|x)p(x)dx
便可以使用下面这个近似形式
I^(y)=1n∑ni=1f(y|xi)

不难发现，在利用蒙特卡洛法进行积分求解时，非常重要的一个环节就是从特定的分布中进行采样。这里的“采样”意思也就是生成满足某种分布的观测值。之前已经介绍过“逆采样”和“拒绝采样”等方法。
采样的目的很多时候都是为了做近似积分运算，前面的采样方法(逆采样和拒绝采样)都是先对分布进行采样，然后再用采样的结果近似计算积分。下面要介绍的另外一种方法“重要性采样”(Importance sampling)则两步并做一步，直接做近似计算积分。
现在的目标是计算下面的积分
E［f(x)］=∫f(x)p(x)dx
按照蒙特卡洛求定积分的方法，将从满足p(x)的概率分布中独立地采样出一系列随机变量xi，然后便有
E［f(x)］≈1n∑ni=1f(xi)
但是现在的困难是对满足p(x)的概率分布进行采样非常困难，毕竟实际中很多p(x)的形式都相当复杂。这时该怎么做呢？于是想到做等量变换，于是有
∫f(x)p(x)dx=∫f(x)p(x)q(x)q(x)dx
如果把其中的f(x)p(x)q(x)看成是一个新的函数h(x)，则有
∫f(x)p(x)dx=∫h(x)q(x)dx≈1n∑ni=1f(xi)p(xi)q(xi)
其中p(xi)q(xi)被称为是xi的权重或重要性权重(Importance Weights)。所以这种采样的方法就被称为是重要性采样(Importance Sampling)。


图319重要性采样原理


如图319所示，在使用重要性采样时，并不会拒绝掉某些采样点，这与在使用拒绝采样时不同。此时，所有的采样点都将为我们所用，但是它们的权重是不同的。因为权重为p(xi)q(xi)，所以在图319中，当p(xi)>q(xi)时，采样点xi的权重就大(深色线在浅色线上方时)，反之就小。重要性采样就是通过这种方式来从一个“参考分布”q(x)中获得“目标分布”p(x)的。


3.5.2马尔科夫链蒙特卡洛的基本概念
马尔科夫链蒙特卡洛方法是一类用于从一个概率分布中进行采样的算法，该类方法以构造马尔科夫链为基础，而被构造的马尔科夫链分布应该同需要之分布等价。
MCMC构造马尔科夫链，使其稳态分布等于要采样的分布，这样就可以通过马尔科夫链来采样。这种等价如何理解是深入探讨具体操作方法之前，需要先攻克的一个问题。在此之前，希望读者对马尔科夫链已经有了一个比较清晰的认识。
现在，用下面的式子来表示每一步(时刻推进)中从状态si到状态sj的转移概率
p(i，j)=p(i→j)=P(Xt+1=sj|Xt=si)
这里的一步是指时间从时刻t过渡到下一时刻t+1。
马尔科夫链在时刻t处于状态j的概率可以表示为πj(t)=P(Xt=sj)。这里用向量π(t)来表示在时刻t各个状态的概率向量。在初始时刻，需要选取一个特定的π(0)，通常情况下可以使向量中一个元素为1，其他元素均为0，即从某个特定的状态开始。随着时间的推移，状态会通过转移矩阵，向各个状态扩散。
马尔科夫链在t+1时刻处于状态si的概率可以通过t时刻的状态概率和转移概率来求得，并可通过查普曼柯尔莫哥洛夫等式得到
πi(t+1)=P(Xt+1=si)=∑kP(Xt+1=si|Xt=sk)P(Xt=sk)

=∑kp(k→i)πk(t)=∑kp(i|k)πk(t)
用转移矩阵写成矩阵的形式如下
π(t+1)=π(t)P

其中，转移矩中的元素p(i，j)表示p(i→j)。因此，π(t)=π(0)Pt。此外，pni，j表示矩阵Pn中第ij个元素，即pni，j=P(Xt+n=sj|Xt=si)。
一条马尔科夫链有一个平稳分布π*，是指给定任意一个初始分布π(0)，经过有限步之后，最终都会收敛到平稳分布π*。平稳分布具有性质π*=π*P。可以结合本章前面谈到的矩阵极限的概念来理解马尔科夫链的平稳分布。
一条马尔科夫链拥有平稳分布的一个充分条件是对于任意两个状态i和j，其满足细致平衡(Detailed Balance)： p(j→i)πj=p(i→j)πi。可以看到，此时
(πP)j=∑iπip(i→j)=∑iπjp(j→i)=πj∑ip(j→i)=πj
所以π=πP，明显满足平稳分布的条件。如果一条马尔科夫链满足细致平衡，就说它是可逆的。


图320MCMC的基本原理



最后来总结一下MCMC的基本思想。在拒绝采样和重要性采样中，当前生成的样本点与之前生成的样本点之间是没有关系的，它的采样都是独立进行的。然而，MCMC是基于马尔科夫链进行的采样。这就表明，当前的样本点生成与上一时刻的样本点是有关的。如图320所示，假设当前时刻生成的样本点是x，下一次时刻采样到x′的(转移)概率就是p(x|x′)，或者写成p(x→x′)。我们希望的是这个过程如此继续下去，最终的概率分布收敛到π(x)，也就是要采样的分布。

易见，转移概率p(x|x′)与采样分布π(x)之间必然存在某种联系，因为希望依照p(x|x′)来进行转移采样最终能收敛到π(x)。那这个关系到底是怎么样的呢？首先，给定所有可能的x，然后求从x转移到x′的概率，所得之结果就是x′出现的概率，如果用公式表示即
π(x′)=∫π(x)p(x′|x)dx
这其实也就是查普曼柯尔莫哥洛夫等式所揭示的关系。如果马尔科夫链是平稳的，那么基于该马尔科夫链的采样过程最终将收敛到平稳状态
π*(x′)=∫π*(x′)p(x|x′)dx′
而平稳状态的充分条件又是满足细致平衡
π(x)p(x′|x)dx=π(x)p(x|x′)
可见细致平衡的意思是说从x转移到x′的概率应该等于从x′转移到x的概率，在这样的情况下，采样过程将最终收敛到目标分布。也就是说，设计的MCMC算法，只要验证其满足细致平衡的条件，那么用这样的方法做基于马尔科夫链的采样，最终就能收敛到一个稳定状态，即目标分布。
实际应用中，有两个马尔科夫链蒙特卡洛采样算法十分常用，它们是MetropolisHastings(梅特罗波利斯黑斯廷斯)算法与Gibbs Sampling采样（吉布斯），可以证明吉布斯采样是梅特罗波利斯黑斯廷斯算法的一种特殊情况，而且两者都满足细致平衡的条件。
3.5.3MetropolisHastings算法
对给定的概率分布，如果从中直接采样比较困难，那么MetropolisHastings算法就是一个从中获取一系列随机样本的MCMC方法。这里所说的MetropolisHastings算法和模拟退火方法中所使用的MetropolisHastings算法本质上是一回事。用于模拟退火和MCMC的原始算法最初是由Metropolis提出的，后来Hastings对其进行了推广。Hastings提出没有必要使用一个对称的参考分布(proposal distribution)。当然达到均衡分布的速度与参考分布的选择有关。
MetropolisHastings算法的执行步骤如下。 



1.初始化x0

2.对于i=0到N-1

u~U(0，1)

x*~q(x*|x(i))

如果u<α(x*)=min1，π(x*)p(x|x*)π(x)p(x*|x)

x(i+1)=x*

否则

x(i+1)=x(i)


MetropolisHastings算法的执行步骤是先随便指定一个初始的样本点x0，u是来自均匀分布的一个随机数，x*是来自p分布的样本点，样本点是否从前一个位置跳转到x*，由α(x*)和u的大小关系决定。如果接受，则跳转，即新的采样点为x*，否则就拒绝，即新的采用点仍为上一轮的采样点。这个算法看起来非常简单，难免让人疑问照此步骤来执行，最终采样分布能收敛到目标分布吗？根据之前的讲解，可知要想验证这一点，就需要检查细致平衡是否满足。下面是具体的证明过程，不再赘言。

π(x)p(x*|x)α(x*)=π(x)p(x*|x)min1，π(x*)p(x|x*)π(x)p(x*|x)

=minπ(x)p(x*|x)，π(x*)p(x|x*)

=π(x*)p(x|x*)min1，π(x)p(x*|x)π(x*)p(x|x*)

=π(x*)p(x|x*)α(x)

既然已经从理论上证明MetropolisHastings算法的确实可行，下面举一个简单的例子来看看实际中MetropolisHastings算法的效果。稍微有点不一样的地方是，这里并未出现π(x)p(x′|x)这样的后验概率形式，作为一个简单的示例，下面只演示从柯西分布中进行采样的做法。
柯西分布也叫作柯西洛伦兹分布，它是以奥古斯丁·路易·柯西与亨德里克·洛伦兹名字命名的连续概率分布，其概率密度函数为
f(x； x0，γ)=1πγ1+x-x0γ2=1πγ(x-x0)2+γ2
其中，x0是定义分布峰值位置的位置参数，γ是最大值一半处的一半宽度的尺度参数。x0=0且γ=1的特例称为标准柯西分布，其概率密度函数为
f(x； 0，1)=1π［1+x2］
下面是在R中进行基于MetropolisHastings算法对柯西分布进行采样的示例代码。



rm(list=ls()) ## 清除全部对象

set.seed(201912)

reps=40000



# target density: the cauchy distribution

cauchy<-function(x， x0=0， gamma=1){

out<-1/(pi*gamma*(1+((x-x0)/gamma)^2))

return(out)

}



chain<-c(0)

for(i in 1:reps){

proposal<-chain［i］+runif(1， min=-1， max=1)

accept<-runif(1)<cauchy(proposal)/cauchy(chain［i］)

chain［i+1］<-ifelse(accept==T， proposal， chain［i］)

}



plot(chain， type="l")

plot(density(chain［1000:reps］)， xlim=c(-5，5)， ylim=c(0，0.4)， col="red")

den<-cauchy(seq(from=-5， to=5， by=0.1)， x=0， gamma=1)

lines(den~seq(from=-5， to=5， by=0.1)， lty=2， col="blue")



图321所示为每次采样点的数值分布情况。
为了更清晰明确地看到采样结果符合预期分布，这里也绘制了分布的密度图并与实际的分布密度进行对照，如图322所示，深色实线是采样分布的密度图，而浅色虚线则是实际柯西分布的密度图，可见两者吻合地相当好。


图321采样点数值分布情况




图322密度图对比


当然，作为一个简单例子的开始，上面的例子中并没有涉及p(y|x)，接下来的这个例子则会演示涉及后验概率的基于MetropolisHastings算法的MCMC。这里将用MetropolisHastings算法从瑞利分布(Rayleigh Distribution)中采样，瑞利分布的密度为
f(x)=xσ2e-x22σ2，x≥0，σ>0
然后取自由度为Xt的卡方分布为参考分布。实现代码如下。 



rm(list=ls()) ## 清除全部对象



f <- function(x， sigma) {

if (any(x < 0)) return (0)

stopifnot(sigma > 0)

return((x / sigma^2) * exp(-x^2 / (2*sigma^2)))

}



m <- 40000

sigma <- 4

x <- numeric(m)

x［1］ <- rchisq(1， df=1)

k <- 0

u <- runif(m)



for (i in 2:m) {

xt <- x［i-1］

y <- rchisq(1， df = xt)

num <- f(y， sigma) * dchisq(xt， df = y)

den <- f(xt， sigma) * dchisq(y， df = xt)

if (u［i］ <= num/den) x［i］ <- y else {

x［i］ <- xt







k <- k+1     #y is rejected

}

}




然后要验证一下，生产的采样数据是否真的符合瑞利分布。注意在R中使用瑞利分布的相关函数，需要加装VGAM包。下面的代码可以让读者直观地感受到采样结果的分布情况。



> curve(drayleigh(x， scale = 4， log = FALSE)， from = -1， to = 15， xlim=c(-1，15)， ylim=c(0，0.2)， 

+ lty = 2， col="blue"，xlab = ""， ylab="") 

> par(new=TRUE) 

> plot(density(x［1000:m］) ， xlim=c(-1，15)， ylim=c(0，0.2)，col="red")




从图323中不难看出，采样点分布确实符合预期。当然，这仅仅是一个演示用的小例子。显然，它并不高效。因为采用自由度为Xt的卡方分布来作为参考函数，大约有40%的采样点都被拒绝掉了。如果换做其他更加合适的参考函数，可以大大提升采样的效率。


图323采样点分布



3.5.4Gibbs采样
Gibbs采样是一种用以获取一系列观察值的MCMC算法，这些观察值近似于来自一个指定的多变量概率分布，但从该分布中直接采样较为困难。Gibbs采样可以被看成是MetropolisHastings算法的一个特例。而这个特殊之处就在于，使用Gibbs采样时，通常是对多变量分布进行采样的。比如说，现在有一个分布p(z1，z2，z3)，可见它是定义在三个变量上的分布。这种分布在实际中还有很多，比如一维的正态分布也是关于其均值和方差这两个变量的分布。在算法的第t步，选出了三个值z(t)1，z(t)2，z(t)3。这时，既然z(t)2和z(t)3是已知的，那么可以由此获得一个新的z(t+1)1，并用这个新的值替代原始的z(t)1，而新的z(t+1)1可以由下面这个条件分布来获得
p(z1|z(t)2，z(t)3)
接下来同样的道理再用一个新的z(t+1)2来替代原始的z(t)2，而这个新的z(t+1)2可以由下面这个条件分布来获得
p(z2|z(t+1)1，z(t)3)
可见刚刚获得的新值z(t+1)1已经被直接用上了。同理，最后用z(t+1)3来更新z(t)3，而新的值由下面这个条件分布来获得
p(z3|z(t+1)1，z(t+1)2)
按照这个顺序如此往复即可。下面给出更为完整的吉布斯采样的算法描述。



1. 初始化{zi： i=1，…，M}

2. 对于τ=1，…，T： 

采样 z(τ+1)1~p(z1|z(τ)2，z(τ)3，…，z(τ)M)

采样 z(τ+1)2~p(z2|z(τ+1)1，z(τ)3，…，z(τ)M)

︙

采样 z(τ+1)j~p(zj|z(τ+1)1，…，z(τ+1)j-1，z(τ)j+1…，z(τ)M)

︙

采样 z(τ+1)M~p(zM|z(τ+1)1，z(τ+1)2，…，z(τ+1)M-1)


当然如果要从理论上证明Gibbs采样确实可以得到预期的分布，就应该考察它是否满足细致平衡。但是前面也讲过，Gibbs采样是MetropolisHastings算法的一个特例。所以其实甚至无需费力去考察其是否满足细致平衡，如果能够证明吉布斯采样就是MetropolisHastings算法，而MetropolisHastings算法是满足细致平衡，其实也就得到了我们想要的。下面是证明的过程，其中x-k表示一个向量(x1，…，xk-1，xk+1，…，xD)，也就是从中剔除了xk元素。
令x=x1，…，xD。
当采样第k个元素，pk(x*|x)=π(x*k|x-k)
当采样第k个元素，x*-k=x-k
π(x*)p(x|x*)π(x)p(x*|x)=π(x*)π(xk|x*-k)π(x)π(x*k|x-k)=π(x*k|x*-k)π(xk|x*-k)π(xk|x-k)π(x*k|x-k)=1
以Gibbs采样为基础实现的MCMC在自然语言处理中的LDA里有重要应用。如果读者正在或者后续准备研究LDA的话，这些内容将是非常重要的基础。