抽样方法(2)--拒绝抽样和高斯分布

mklf

2016-10-03 11:39

在上篇博客抽样方法(1)--随机数和简单随机抽样中，我介绍了如何生成随机数和均匀分布。进而通过均匀分布我们可以用Inverse Transform Method对一些简单分布进行抽样，但在博客的最后,我提到到由于高斯分布的累积分布函数是不可积的，因此用 Inverse Transform Method 不能直接对高斯分布进行抽样。这篇博客，我们就来为高斯分布寻找一个抽样方法。在继续后面的内容之前，注意到由于$N(\mu,\delta^2)=\mu+\delta N(0,1)$,因此要生成均值为$\mu$，方差为$\delta^2$的高斯分布，只需要我们能生成标准的高斯分布即可。针对前面高斯分布抽样的存在问题，我们可以从两个方向入手：

既然高斯分布的CDF: $F(x) = \int_{-\infty}^x \frac{1}{\sqrt{2\pi}} exp(-\frac{1}{2} x^2) dx$是不可积的，那么我们可以通过一些手段使得上式可积或者得到这个积分的一个近似。
采用其他不需要计算累计分布函数的抽样方法。

解决了第一个方向的问题，我们就可以继续用Inverse Transform Method 来对高斯分布进行抽样，著名的Box-Muller变换就是从第一种方法入手，它每次抽样两个独立的高斯分布样本，这样累计分布函数就变为一个二重积分:

$$\begin{aligned}F(x,y)= \int_{-\infty}^y \int_{-\infty}^x \frac{1}{2\pi} exp(-\frac{1}{2} (x^2 + y^2) ) dxdy\quad \end{aligned}$$

上面的这个积分在很多微积分教材中都有提及(比如我学习高数时用的同济第六版)，把x,y变换到极坐标下就可以求得积分。这样，在求得极坐标下的样本后，再通过直角坐标到极坐标的逆变换就可以到到x,y的值,即高斯分布的样本。除此之外，另一种和Box-Muller变换类似的方法叫Marsaglia polar method，这种方法相比 Box-Muller变换，由于不需要在逆变换过程计算sin,cos这样的超越函数，节省了不少运算量，从而使得抽样的速度比Box-Muller变换更快。

上面的两种方法已经能很好的解决高斯分布的抽样问题了，不过由于这个系列的博客的目的是介绍抽样方法，所以我们主要从第二个方向入手来解决这个问题，也就是用这篇博客的主题——拒绝采样来对高斯分布进行抽样。

1.Rejection sampling:¶

Rejection sampling 的想法是，假设我们能生成概率分布为$p(x)$的简单随机变量(这个分布也叫envelope分布)，现在要生成概率分布为$q(x)$的随机变量，我们可以从p中生成样本Y，再按照一定的概率接受这个样本。设常量c满足:

$$\begin{aligned}\forall y,c\geq\frac{q(y)}{p(y)} \end{aligned}\label{1}\tag{公式1}$$

拒绝采样的步骤如下:

生成具有概率分布p的随机变量的样本Y
生成均匀分布的样本U
如果 $ U\leq \frac{q(Y)}{cp(Y)}$ 则接受Y，令X=Y，否则回到第1步

通过以上步骤得到的样本X满足概率分布为为q(x)

我们注意到在一次迭代过程中，一个样本值为y且被接受的概率为: $$P(接受,x=y) \\{}= P\{Y=y\}P\{接受|Y=y\} \\{}= p(y)*\frac{q(y)}{cp(y)}=\frac{q(y)}{c}$$

通过对y求和，我们求得算法接受任意样本Y的概率为: $$P(接受) = \sum_{y}{\frac{q(y)}{c}}=\frac{1}{c}$$

由于每一次迭代是独立进行的，上面的过程看做是一个几何概型：整个过程以$P(接受)=\frac{1}{c}$的概率接受，如果不接受，则继续进行,直到第一次接受为止。我们知道参数为$p$的几何分布的期望为$E(x)=\frac{1}{p}$,因此接受率为$\frac{1}{c}$的拒绝采样中，平均每接受一个样本需要进行$c$次迭代。我们通常希望算法的迭代次数尽量少，,因此c应取为$\frac{q(y)}{p(y)}$的上确界。

注意到拒绝采样是不需要对目标分布进行积分的，这样我们就可以愉快的对高斯分布进行采样了。根据拒绝采样算法，首先我们要选择一个我们能进行采样的分布来作为envelope分布，上篇博客中我们实现对指数分布进行采样，这里正好能用上，那我们就以均值为1的指数分布作为 envelope 分布。但由于指数分布是非负的，我们把高斯分布也限制在正半轴,这样正半轴的高斯分布密度函数为$q(x) = \frac{2}{\sqrt{2\pi}}e^{-x^2/2},0<x<\infty$,而作为envelope分布的均值为1的指数分布密度函数为$p(x)=e^{-x},0<x<\infty$

根据拒绝采样算法，第二步是计算$c$的取值，$c$ 应该为$\frac{q(x)}{p(x)}$的上确界:

$$\begin{equation}\frac{q(x)}{p(x)}=\sqrt{2/\pi}e^{x-x^2/2}\\{}c=\max\frac{q(x)}{p(x)}\end{equation}$$

通过求导可以得到当$x=1$时,$\frac{q(x)}{p(x)}$的最大值，此时$c = max \frac{q(x)}{p(x)} =\sqrt{2e/\pi}$,带入c的值，可以得到 $$\frac{q(x)}{c\cdot p(x)} =\exp(-\frac{(x-1)^2}{2})$$

最后，以均值为1的指数分布为envelope分布,我们的拒绝采样的高斯分布抽样算法为：

生成一个$\lambda =1$ 的指数分布样本y
生成一个均匀分布样本U
如果$u\leq exp\{-(Y-1)^2/2\}$,则令x=y,否则返回步骤1

上面的拒绝采样算法生成了正半轴上的高斯分布数据，要想获得整个实数域上的高斯分布，利用高斯分布的对称性，我们知道对每一个样本点Y，它属于正半轴和负半轴的概率都是0.5，因此我们再生成一个均匀分布样本U，如果U>0.5,X = Y, 否则X=-Y，代表它在负半轴。

下面用代码来实现这个算法,然后检验它所产生的样本是否符合均值为0,方差为1的高斯分布。首先定义指数分布的采样函数(方法和上篇博客相同)

In [1]:

%pylab inline
from functools import partial
def exponential_sampling(param_lambda,N=10000):
    """
        生成指数分布样本
        param_lambda :float, 参数\lambda 
        N: int, 样本点数目
        return : 产生的样本
    """
    assert(param_lambda>0)
    x = random.rand(N)
    return (-1/param_lambda)*log(x)

Populating the interactive namespace from numpy and matplotlib

然后实现拒绝采样的通用函数以及产生高斯分布的拒绝采样

In [3]:

def rejection_sampling(envelope_sampler,cond,N):
    """
        拒绝采样
        envelope_sampler: envelope分布样本生成函数,
        返回样本Y
        cond:cond(Y),计算q(y)/(c*p(y)) 
        N: int,迭代次数
        return :numpy.ndarray 样本点
    """
    Y = envelope_sampler(N)
    U = random.rand(N)
    T = cond(Y)
    accept = U<T
    return np.extract(accept,Y)
exp_sampler = partial(exponential_sampling,1)
unary_op = lambda y: exp(-(y-1)**2/2)

def gaussian_sampler(exp_sampler,unary_op,N):
    pos = rejection_sampling(exp_sampler,unary_op,N)
    U = random.rand(pos.shape[0])>0.5
    pos[U] = -pos[U]
    return pos
gaussian_sampler = partial(gaussian_sampler,
                           exp_sampler,unary_op)

我们来验证一下算法的参数的分布，首先看看抽样分布的均值和方差

In [4]:

N = 1000000
samples = gaussian_sampler(N)
s_mean,s_var = mean(samples),var(samples)
print("sample mean:{:.6f},variance:{:.6f}".format(s_mean,s_var))

sample mean:-0.000462,variance:1.001364

进一步用下面的代码画出分布函数,图中红色的区现对应的是我们算法得到的样本密度函数，绿线是numpy的高斯样本的抽样分布。

In [5]:

hist(samples,bins=200,histtype='step',normed=True,color='r');
hist(random.randn(100000),bins=200,histtype='step',normed=True,color='g');

最后我们来检查一下样本的接受率，我们在前面推算出拒绝采样的平均采样次数为$c$,在这个例子中$c=\sqrt{2e/\pi}$,算法接受样本的概率大约应该为$\frac{1}{c}\simeq 0.760$,我们可以用一行代码来验证计算这个值:

In [6]:

print('accpet probability: {}'.format(len(samples)/float(N)))

accpet probability: 0.759996

2.Adaptive Rejection Sampling¶

拒绝采样的 envelope 分布极大的影响了采样效率，如果envelope分布和目标分布形状类似，此时$c=\max \frac{q(x)}{p(x)}$小，平均采样次数就少，相反，如果在目标分布概率大的区域，envelope分布的概率反而小，根据$\ref{1}$, c的值将很大，从而导致采样相同数目的样本需要更多的采样次数。基于这个原因，我们总是希望envelope分布和目标分布越接近越好。一种特殊的情况是，如果目标分布是log concave的分布，由于concave函数的函数值总是小于或等于切线上的函数值，我们在函数上寻找多个点，分布做出这些点上的切线，然后用各个切线段得到的分段线性函数把目标分布的分布函数包裹起来,而且随着切线的数目增多，这个分段线性函数和目标分布函数会越来越接近。注意这个分段线性函数在log scale下是线性的，因此原分布是一个分段的指数分布，而指数分布是易于抽样的。我们用这个和目标分布非常接近的分布来进行采样，可以得到很高的样本接受率。

还是以高斯分布为例，它的log 密度函数$\log P(x) = -\frac{1}{2}\cdot x^2+constant$,因此高斯分布是 log concave 的，我们可以用分段线性函数来近似它。下面这段代码分别用4,7,10个线性函数来近似高斯分布的log密度函数。可以看到随着线性函数增多，这个分段线性函数将高斯分布包裹得越来越紧,从而在采样时接受率越来越高。

In [6]:

# 高斯核以及用分段线性函数包裹高斯核函数
def origin_concave(x):
    return -1/2.0*(x**2-sqrt(2*pi))
def grad_of_origin(x):
    return -x
def pisewise_liear_of(origin,grad,sample_points):
    y = origin(sample_points)
    gradient = grad(sample_points)
    intercept = y-gradient*sample_points
    def func(x):
        y_approx = np.outer(x,gradient)+intercept
        return np.min(y_approx,axis=1)
    return func

figsize(16,4)
for i,n in enumerate([4,7,10]):
    subplot(1,3,i+1)
    # 随机选取n个点
    samples = np.linspace(-5,5,n)+np.random.rand(n)-0.5
    #对这n个点求分段线性函数
    pisewise_linear = pisewise_liear_of(origin_concave,grad_of_origin,samples)
    x = np.arange(-5,5,0.01)
    linear = pisewise_linear(x)
    plot(x,origin_concave(x),x,linear)
    xlabel("{} linears".format(n))

无论是Inverse transform method 还是 Rejection sampling，当它们面对高维分布的时候都存在各种问题。不幸的是，在机器学习中我们面对的问题很多都是高维的概率分布，因此我们需要一种在高维下相对高效的抽样方法。下篇博客我将介绍MCMC这一系列更加强大的抽样算法。

1.Rejection sampling:¶

2.Adaptive Rejection Sampling¶

评论