采样方法（一）

https://blog.csdn.net/dark_scope/article/details/70992266?utm_source&＃61;itdadao&utm_medium&＃61;referral

引子
最近开始拾起来看一些NLP相关的东西&＃xff0c;特别是深度学习在NLP上的应用&＃xff0c;发现采样方法在很多模型中应用得很多&＃xff0c;因为训练的时候如果预测目标是一个词&＃xff0c;直接的softmax计算量会根据单词数量的增长而增长。恰好想到最开始深度学习在DBN的时候采样也发挥了关键的作用&＃xff0c;而自己对采样相关的方法了解不算太多&＃xff0c;所以去学习记录一下&＃xff0c;经典的统计的方法确实巧妙&＃xff0c;看起来非常有收获。

本篇文章先主要介绍一下经典的采样方法如Inverse Sampling、Rejective Sampling以及Importance Sampling和它在NLP上的应用&＃xff0c;后面还会有一篇来尝试介绍MCMC这一组狂炫酷拽的算法。才疏学浅&＃xff0c;行文若有误望指正。

Why Sampling
采样是生活和机器学习算法中都会经常用到的技术&＃xff0c;一般来说采样的目的是评估一个函数在某个分布上的期望值&＃xff0c;也就是 
E[f(x)],x∼p,p is a distribution.
E[f(x)],x∼p,p is a distribution.

比如我们都学过的抛硬币&＃xff0c;期望它的结果是符合一个伯努利分布的&＃xff0c;定义正面的概率为pp,反面概率为1−p1−p。最简单地使f(x)&＃61;xf(x)&＃61;x&＃xff0c;在现实中我们就会通过不断地进行抛硬币这个动作&＃xff0c;来评估这个概率p。 
E[f(x)]≈1m∑i&＃61;1mf(xi).  xi∼p
E[f(x)]≈1m∑i&＃61;1mf(xi).  xi∼p

这个方法也叫做蒙特卡洛法&＃xff08;Monte Carlo Method&＃xff09;&＃xff0c;常用于计算一些非常复杂无法直接求解的函数期望。 
对于抛硬币这个例子来说: 
E[f(x)]&＃61;p≈1m∑i&＃61;1mxi&＃61;cntum
E[f(x)]&＃61;p≈1m∑i&＃61;1mxi&＃61;cntum

其期望就是抛到正面的计数cntucntu除以总次数mm。 
而我们抛硬币的这个过程其实就是采样&＃xff0c;如果要用程序模拟上面这个过程也很简单&＃xff0c;因为伯努利分布的样本很容易生成&＃xff1a; 
yi∼Uniform(0,1)&＃xff0c;so xi&＃61;I(y yi∼Uniform(0,1)&＃xff0c;so xi&＃61;I(y

而在计算机中的随机函数一般就是生成0到1的均匀分布随机数。
Sampling Method
可以看到蒙特卡洛法其实就是按一定的概率分布中获取大量样本&＃xff0c;用于计算函数在样本的概率分布上的期望。其中最关键的一个步骤就是如何按照指定的概率分布p进行样本采样&＃xff0c;抛硬币这个case里伯努利分布是一个离散的概率分布&＃xff0c;它的概率分布一般用概率质量函数&＃xff08;pmf&＃xff09;表示&＃xff0c;相对来说比较简单&＃xff0c;而对于连续概率分布我们需要考虑它的概率密度函数&＃xff08;pdf&＃xff09;&＃xff1a; 

比如上图示例分别是标准正态分布概率密度函数&＃xff0c;它们的面积都是1&＃xff08;这是概率的定义&＃xff09;&＃xff0c;如果我们可以按照相应概率分布生成很多样本&＃xff0c;那这些样本绘制出来的直方图应该跟概率密度函数是一致的。 
而在实际的问题中&＃xff0c;p的概率密度函数可能会比较复杂&＃xff0c;我们由浅入深&＃xff0c;看看如何采样方法如何获得服从指定概率分布的样本。

Inverse Sampling
对于一些特殊的概率分布函数&＃xff0c;比如指数分布&＃xff1a; 
pexp(x)&＃61;{λexp(−λx)0,x≥0,x<0
pexp(x)&＃61;{λexp(−λx),x≥00,x<0

我们可以定义它的概率累积函数(Cumulative distribution function)&＃xff0c;也就是(ps.这个’F’和前面的’f’函数并没有关系) 
F(x)&＃61;∫x−∞p(x)dx
F(x)&＃61;∫−∞xp(x)dx

从图像上看就是概率密度函数小于x部分的面积。这个函数在x≥0x≥0的部分是一个单调递增的函数(在定义域上单调非减)&＃xff0c;定义域和值域是[0,&＃43;∞)→[0,1)[0,&＃43;∞)→[0,1)&＃xff0c;画出来大概是这样子的一个函数&＃xff0c;在p(x)p(x)大的地方它增长快&＃xff08;梯度大&＃xff09;&＃xff0c;反之亦然&＃xff1a;

因为它是唯一映射的&＃xff08;在>0的部分&＃xff0c;接下来我们只考虑这一部分&＃xff09;&＃xff0c;所以它的反函数可以表示为F−1(a)&＃xff0c;a∈[0,1),值域为[0,&＃43;∞)F−1(a)&＃xff0c;a∈[0,1),值域为[0,&＃43;∞)
因为F单调递增&＃xff0c;所以F−1F−1也是单调递增的&＃xff1a; 
x≤ya≤b⇒F(x)≤F(y)⇒F−1(a)≤F−1(b)
x≤y⇒F(x)≤F(y)a≤b⇒F−1(a)≤F−1(b)

利用反函数的定义&＃xff0c;我们有&＃xff1a; 
F−1(a) F−1(a)

我们定义一下[0,1]均匀分布的CDF,这个很好理解&＃xff1a; 
P(a≤x)&＃61;H(x)&＃61;⎧⎩⎨1x0,x≥1,0≤x≤1,x<0
P(a≤x)&＃61;H(x)&＃61;{1,x≥1x,0≤x≤10,x<0

所以 
P(F−1(a)≤x)&＃61;P(a≤F(x))&＃61;H(F(x))因为F(x)的值域[0,1)&＃xff0c;所以上式P(F−1(a)≤x)&＃61;H(F(x))&＃61;F(x)
P(F−1(a)≤x)&＃61;P(a≤F(x))&＃61;H(F(x))因为F(x)的值域[0,1)&＃xff0c;所以上式P(F−1(a)≤x)&＃61;H(F(x))&＃61;F(x)

根据F(x)F(x)的定义&＃xff0c;它是exp分布的概率累积函数&＃xff0c;所以上面这个公式的意思是F−1(a)F−1(a)符合exp分布,我们通过F的反函数将一个0到1均匀分布的随机数转换成了符合exp分布的随机数&＃xff0c;注意&＃xff0c;以上推导对于cdf可逆的分布都是一样的&＃xff0c;对于exp来说&＃xff0c;它的反函数的形式是&＃xff1a; 
F−1exp(a)&＃61;−1λ∗log(1−a)
Fexp−1(a)&＃61;−1λ∗log(1−a)

具体的映射关系可以看下图(a)&＃xff0c;我们从y轴0-1的均匀分布样本&＃xff08;绿色&＃xff09;映射得到了服从指数分布的样本&＃xff08;红色&＃xff09;。 

我们写一点代码来看看效果,最后绘制出来的直方图可以看出来就是exp分布的图&＃xff0c;见上图(b)&＃xff0c;可以看到随着采样数量的变多&＃xff0c;概率直方图和真实的CDF就越接近&＃xff1a;
def sampleExp(Lambda &＃61; 2,maxCnt &＃61; 50000):
    ys &＃61; []
    standardXaxis &＃61; []
    standardExp &＃61; []
    for i in range(maxCnt):
        u &＃61; np.random.random()
        y &＃61; -1/Lambda*np.log(1-u) #F-1(X)
        ys.append(y)
    for i in range(1000):
        t &＃61; Lambda * np.exp(-Lambda*i/100)
        standardXaxis.append(i/100)
        standardExp.append(t)
    plt.plot(standardXaxis,standardExp,&＃39;r&＃39;)
    plt.hist(ys,1000,normed&＃61;True)
    plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
Rejective Sampling
我们在学习随机模拟的时候通常会讲到用采样的方法来计算ππ值&＃xff0c;也就是在一个1×1的范围内随机采样一个点&＃xff0c;如果它到原点的距离小于1,则说明它在1/4圆内&＃xff0c;则接受它&＃xff0c;最后通过接受的占比来计算1/4圆形的面积&＃xff0c;从而根据公式反算出预估的ππ值&＃xff0c;随着采样点的增多&＃xff0c;最后的结果π^π^会越精准。 

上面这个例子里说明一个问题&＃xff0c;我们想求一个空间里均匀分布的集合面积&＃xff0c;可以尝试在更大范围内按照均匀分布随机采样&＃xff0c;如果采样点在集合中&＃xff0c;则接受&＃xff0c;否则拒绝。最后的接受概率就是集合在‘更大范围’的面积占比。

当我们重新回过头来看想要sample出来的样本服从某一个分布p&＃xff0c;其实就是希望样本在其概率密度函数p(x)p(x)高的地方出现得更多&＃xff0c;所以一个直觉的想法&＃xff0c;我们从均匀分布随机生成一个样本xixi,按照一个正比于p(xi)p(xi)的概率接受这个样本&＃xff0c;也就是说虽然是从均匀分布随机采样&＃xff0c;但留下的样本更有可能是p(x)p(x)高的样本。

这样的思路很自然&＃xff0c;但是否是对的呢。其实这就是Rejective Sampling的基本思想&＃xff0c;我们先看一个很intuitive的图 

假设目标分布的pdf最高点是1.5,有三个点它们的pdf值分别是 
p(x1)p(x2)p(x3)&＃61;1.5;&＃61;0.5;&＃61;0.3;
p(x1)&＃61;1.5;p(x2)&＃61;0.5;p(x3)&＃61;0.3;

因为我们从x轴上是按均匀分布随机采样的&＃xff0c;所以采样到三个点的概率都一样&＃xff0c;也就是 
q(x1)&＃61;q(x2)&＃61;q(x3)
q(x1)&＃61;q(x2)&＃61;q(x3)

接下来需要决定每个点的接受概率acc(xi)acc(xi)&＃xff0c;它应该正比于p(xi)p(xi)&＃xff0c;当然因为是概率值也需要小于等于1. 
我们可以画一根y&＃61;2y&＃61;2的直线&＃xff0c;因为整个概率密度函数都在这根直线下&＃xff0c;我们设定 
acc(xi)&＃61;p(xi)2;所以acc(x1)&＃61;0.75;acc(x2)&＃61;0.25;acc(x3)&＃61;0.15;
acc(xi)&＃61;p(xi)2;所以acc(x1)&＃61;0.75;acc(x2)&＃61;0.25;acc(x3)&＃61;0.15;

我们要做的就是生成一个0-1的随机数xixi&＃xff0c;如果它小于接受概率acc(xi)acc(xi)&＃xff0c;则留下这个样本。因为acc∝pacc∝p&＃xff0c;所以可以看到因为p(x1)p(x1)是p(x2)p(x2)的3倍&＃xff0c;所以acc(x1)&＃61;3acc(x2)acc(x1)&＃61;3acc(x2)。同样采集100次&＃xff0c;最后留下来的样本数期望也是3倍。这根本就是概率分布的定义!
我们将这个过程更加形式化一点&＃xff0c;我们我们又需要采样的概率密度函数p(x)p(x)&＃xff0c;但实际情况我们很有可能只能计算出p~(x)p~(x),有p(x)&＃61;p~(x)Zpp(x)&＃61;p~(x)Zp。我们需要找一个可以很方便进行采样的分布函数q(x)q(x)并使 
cq(x)≥p~(x)
cq(x)≥p~(x)

其中c是需要选择的一个常数。然后我们从qq分布中随机采样一个样本xixi,并以 
acc(xi)&＃61;p~(xi)cq(xi)
acc(xi)&＃61;p~(xi)cq(xi)

的概率决定是否接受这个样本。重复这个过程就是「拒绝采样」算法了。
在上面的例子我们选择的q分布是均匀分布&＃xff0c;所以从图像上看其pdf是直线&＃xff0c;但实际上cq(x)cq(x)和p~(x)p~(x)越接近&＃xff0c;采样效率越高&＃xff0c;因为其接受概率也越高&＃xff1a; 

Importance Sampling
上面描述了两种从另一个分布获取指定分布的采样样本的算法&＃xff0c;对于1.在实际工作中&＃xff0c;一般来说我们需要sample的分布都及其复杂&＃xff0c;不太可能求解出它的反函数&＃xff0c;但p(x)p(x)的值也许还是可以计算的。对于2.找到一个合适的cq(x)cq(x)往往很困难&＃xff0c;接受概率有可能会很低。 
那我们回过头来看我们sample的目的&＃xff1a;其实是想求得E[f(x)],x∼pE[f(x)],x∼p,也就是 
E[f(x)]&＃61;∫xf(x)p(x)dx
E[f(x)]&＃61;∫xf(x)p(x)dx

如果符合p(x)分布的样本不太好生成&＃xff0c;我们可以引入另一个分布q(x)q(x)&＃xff0c;可以很方便地生成样本。使得 
∫xf(x)p(x)dxwhere  g(x)&＃61;∫xf(x)p(x)q(x)q(x)dx&＃61;∫xg(x)q(x)dx&＃61;f(x)p(x)q(x)&＃61;f(x)w(x)
∫xf(x)p(x)dx&＃61;∫xf(x)p(x)q(x)q(x)dx&＃61;∫xg(x)q(x)dxwhere  g(x)&＃61;f(x)p(x)q(x)&＃61;f(x)w(x)

我们将问题转化为了求g(x)g(x)在q(x)q(x)分布下的期望&＃xff01;&＃xff01;&＃xff01; 
我们称其中的w(x)&＃61;p(x)q(x)w(x)&＃61;p(x)q(x) 叫做Importance Weight.
Importance Sample 解决的问题
首先当然是我们本来没办法sample from p&＃xff0c;这个是我们看到的&＃xff0c;IS将之转化为了从q分布进行采样&＃xff1b;同时IS有时候还可以改进原来的sample&＃xff0c;比如说: 

可以看到如果我们直接从p进行采样&＃xff0c;而实际上这些样本对应的f(x)f(x)都很小&＃xff0c;采样数量有限的情况下很有可能都无法获得f(x)f(x)值较大的样本&＃xff0c;这样评估出来的期望偏差会较大&＃xff1b;

而如果我们找到一个qq分布&＃xff0c;使得它能在f(x)∗p(x)f(x)∗p(x)较大的地方采集到样本&＃xff0c;则能更好地逼近[Ef(x)][Ef(x)]&＃xff0c;因为有Importance Weight控制其比重&＃xff0c;所以也不会导致结果出现过大偏差。

所以选择一个好的p也能帮助你sample出来的效率更高&＃xff0c;要使得f(x)p(x)f(x)p(x)较大的地方能被sample出来。

无法直接求得p(x)的情况
上面我们假设g(x)g(x)和q(x)q(x)都可以比较方便地计算&＃xff0c;但有些时候我们这个其实是很困难的&＃xff0c;更常见的情况市我们能够比较方便地计算p~(x)p~(x)和q~(x)q~(x) 
p(x)&＃61;p~(x)Zpq(x)&＃61;q~(x)Zq
p(x)&＃61;p~(x)Zpq(x)&＃61;q~(x)Zq

其中Zp/qZp/q是一个标准化项&＃xff08;常数&＃xff09;&＃xff0c;使得p~(x)p~(x)或者p~(x)p~(x)等比例变化为一个概率分布&＃xff0c;你可以理解为softmax里面那个除数。也就是说 
Zp&＃61;∫xp~(x)dxZq&＃61;∫xq~(x)dx
Zp&＃61;∫xp~(x)dxZq&＃61;∫xq~(x)dx

这种情况下我们的importance sampling是否还能应用呢&＃xff1f; 
∫xf(x)p(x)dxwhere  g~(x)&＃61;∫xf(x)p(x)q(x)q(x)dx&＃61;∫xf(x)p~(x)/Zpq~(x)/Zqq(x)dx&＃61;ZqZp∫xf(x)p~(x)q~(x)q(x)dx&＃61;ZqZp∫xg~(x)q(x)dx&＃61;f(x)p~(x)q~(x)&＃61;f(x)w~(x)
∫xf(x)p(x)dx&＃61;∫xf(x)p(x)q(x)q(x)dx&＃61;∫xf(x)p~(x)/Zpq~(x)/Zqq(x)dx&＃61;ZqZp∫xf(x)p~(x)q~(x)q(x)dx&＃61;ZqZp∫xg~(x)q(x)dxwhere  g~(x)&＃61;f(x)p~(x)q~(x)&＃61;f(x)w~(x)
而ZqZpZqZp我们直接计算并不太好计算&＃xff0c;而它的倒数&＃xff1a; 
ZpZq&＃61;1Zq∫xp~(x)dx而Zq&＃61;q~(x)/q(x)所以ZpZq&＃61;∫xp~(x)q~(x)q(x)dx&＃61;∫xw~(x)q(x)dx
ZpZq&＃61;1Zq∫xp~(x)dx而Zq&＃61;q~(x)/q(x)所以ZpZq&＃61;∫xp~(x)q~(x)q(x)dx&＃61;∫xw~(x)q(x)dx

因为我们家设能很方便地从q采样&＃xff0c;所以上式其实又被转化成了一个蒙特卡洛可解的问题&＃xff0c;也就是 
ZpZq&＃61;1m∑i&＃61;1mw~(xi).   xi∼q
ZpZq&＃61;1m∑i&＃61;1mw~(xi).   xi∼q

最终最终&＃xff0c;原来的蒙特卡洛问题变成了&＃xff1a; 
Ef(x)&＃61;∑i&＃61;1mw^(xi)f(xi).xi∼q其中w^(xi)&＃61;w~(xi)∑mj&＃61;1w~(xj)
Ef(x)&＃61;∑i&＃61;1mw^(xi)f(xi).xi∼q其中w^(xi)&＃61;w~(xi)∑j&＃61;1mw~(xj)

所以我们完全不用知道q(x)确切的计算值&＃xff0c;就可以近似地从中得到在q分布下f(x)f(x)的取值&＃xff01;&＃xff01;amazing&＃xff01;
Importance Sampling在深度学习里面的应用
在深度学习特别是NLP的Language Model中&＃xff0c;训练的时候最后一层往往会使用softmax函数并计算相应的梯度。 

而我们知道softmax函数的表达式是&＃xff1a; 
P(yi)&＃61;softmax(xTwi)&＃61;exTwiZZ&＃61;∑j&＃61;1|V|exTwj
P(yi)&＃61;softmax(xTwi)&＃61;exTwiZZ&＃61;∑j&＃61;1|V|exTwj

要知道在LM中m的大小是词汇的数量决定的&＃xff0c;在一些巨大的模型里可能有几十万个词&＃xff0c;也就意味着计算Z的代价十分巨大。
而我们在训练的时候无非是想对softmax的结果进行求导&＃xff0c;也就是说 
Δθ(logP(yi))&＃61;Δθ(xTwi)−∑y′∈VP(y′)Δθ(xTw′)
Δθ(logP(yi))&＃61;Δθ(xTwi)−∑y′∈VP(y′)Δθ(xTw′)
后面那一块&＃xff0c;我们好像看到了熟悉的东西&＃xff0c;没错这个形式就是为采样量身定做似的。 
∑y′∈VP(y′)Δθ(xTw′)&＃61;EP[Δθ(xTw′)]
∑y′∈VP(y′)Δθ(xTw′)&＃61;EP[Δθ(xTw′)]

经典的蒙特卡洛方法就可以派上用途了&＃xff0c;与其枚举所有的词&＃xff0c;我们只需要从V里sample出一些样本词&＃xff0c;就可以近似地逼近结果了。
同时直接从PP中sample也不可取的&＃xff0c;而且计算PP是非常耗时的事情&＃xff08;因为需要计算Z&＃xff09;&＃xff0c;我们一般只能计算P~(y)P~(y)&＃xff0c;而且直接从PP中sample也不可取&＃xff0c;所以我们选择另一个分布QQ进行Importance Sample即可。

一般来说可能选择的QQ分布是简单一些的n−gramn−gram模型。下面是论文中的算法伪代码&＃xff0c;基本上是比较标准的流程&＃xff08;论文图片的符号和上面的描述稍有出入&＃xff0c;理解一下过程即可&＃xff09;&＃xff1a; 

References
【1】mathematicalmonk’s machine learning course on y2b. machine learing 
【2】Pattern Recognition And Machine Learning 
【3】Adaptive Importance Sampling to Accelerate Training 
of a Neural Probabilistic Language Model.Yoshua Bengio and Jean-Sébastien Senécal.
--------------------- 
作者&＃xff1a;Dark_Scope 
来源&＃xff1a;CSDN 
原文&＃xff1a;https://blog.csdn.net/dark_scope/article/details/70992266 
版权声明&＃xff1a;本文为博主原创文章&＃xff0c;转载请附上博文链接&＃xff01;