基于接受拒绝法的随机抽样

基于接受拒绝法的随机抽样

在通常面对的一些实际问题中的概率分布可能不是常见的哪几种概率分布，如何能够生成服从特殊分布的随机数呢？生成服从指定分布的随机数序列的过程我们称之为抽样过程，本文将介绍基于接受拒绝发的随机抽样过程。

接受拒绝法的思想可以形象地理解为制作沙雕，经历由粗到细的雕琢过程。该方法首先要像沙雕堆砌出一个简单雏形那样，设置一个形状简单的概率密度曲线，使得抬高它到一定的高度就能够完全罩住所要的抽样的概率密度曲线。前者的概率密度函数称为建议概率密度函数，而后者称为目标概率密度函数，当然两者的定义域要相同。因为人们都知道，做沙雕时，最初堆出来的那堆沙堆要比最终的雕像大。接着就好似将多出去部分削去的雕刻那样，我们将从建议概率密度函数生成的随机数按照一定的判定概率予以接受或者拒绝。用这样的方式来产生的随机数就服从概率密度函数为f(x)的分布了。

接受拒绝法的具体算法步骤如下：
（1）首先选择一个容易抽样的某个分布作为建议分布，设它为$g(x)$，接着要确定一个常数$M>1$，使得在x的定义域上均有$f(x)\le Mg(x)$成立（这相当于$Mg(x)$的图形完全罩住了$f(x)$），如下图所示

（2）生成服从概率密度函数为$g(x)$的建议随机数y；
（3）再生成一个服从均匀分布$U(0,1)$的随机数u；
（4）计算接收准则概率函数$h(x)=\frac{f(x)}{Mg(x)}$，如果$u<h(y)$，则接受所生成的随机数y；反之，则丢弃该随机数y，并转到第（2）步。
由此所生成的随机数就是我们需要的服从概率密度函数为$f(x)$的随机数。

举例：
假设我们要抽样的目标概率密度函数是：
$$f(x)=\frac{6(x-0.5{)}^2}{7},x\in (0,2)$$
这个密度函数如下图所示，

生成服从它的随机数的具体算法如下：
（1）我们选用在相同定义域上的均匀分布U(0, 2)作为建议分布，即$g(x)=0.5$，由于要求$f(x)<Mg(x)$成立，不难求得f(x)在(0, 2)的最大值为1.9286，则我们可以选择M=3.858。因为M选得越大，意味着我们堆砌的原始雏形就越大，需要削去的部分也越多，效率会越低，所以我们要选择尽可能小的M，上图中就绘出了接受-拒绝法所要求的几条曲线的位置
（2）生成建议的随机数y~U(0, 2);
（3）然后在省城一个随机数u~U(0, 1)
（4）计算接收准则函数$h(y)=\frac{f(y)}{Mg(y)}$，如果u<h(y)，则生成的y就是我们需要的随机数；反之则丢弃生成的y，并转至第（2）步重新生成；
不断重复（2）到（4）步骤就可以生成出一组服从目标概率密度函数f(x)的随机数序列

N = 500000;
M = 3.858;
y = unifrnd(0, 2, 1, N);
u = unifrnd(0, 1, 1, N);
gy = 0.5;
fy = 6 * (y - 0.5).^2 / 7;
x = y(u < fy ./ gy / M);
sample = length(x);
[xNum, xCen] = hist(x, 50);
bar(xCen, xNum / sample);
title('f(x) = 6(x-0.5)^2/7')
hold on
t = 0:0.04:2;
z = 6 * (t-0.5).^2 / 7;
scale = 2 / 50;
plot(t,  scale * z, 'r')
hold off

接受-拒绝法的原理：
接受拒绝法的原理很直观，证明也简单。设F和G分别是对应目标概率密度函数f和建议概率密度函数g的分布函数，我们只需要证明
$$P(Y\le y|U\le \frac{f(Y)}{Mg(Y)})=F(y)$$

记事件$A=Y\le y$和$B=U\le \frac{f(Y)}{Mg(Y)}$，因为我们知道，条件概率
$$P(U\le \frac{f(Y)}{Mg(Y)}|Y=y)=\frac{f(y)}{Mg(y)}$$

那么，它的无条件概率是
$$\begin{array}{rl}P(U\le \frac{f(Y)}{Mg(y)})& ={\int }_{-\infty }^{\infty }P(U\le \frac{f(Y)}{Mg(y)}|Y=y)g(y)dy\\ & ={\int }_{-\infty }^{\infty }\frac{f(y)}{Mg(y)}g(y)dy\\ & =\frac{1}{M}{\int }_{-\infty }^{\infty }f(y)dy=\frac{1}{M}\end{array}$$

这个概率反映了算法的效率，M越大，则效率越低。
由条件概率$P(A|B)=P(B|A)P(A)/P(B)$可得
$$\begin{array}{rl}P(Y\le y|U\le \frac{f(Y)}{Mg(Y)})& =P(A|B)=\frac{P(B|A)P(A)}{P(B)}\\ & =\frac{P(U\le \frac{f(Y)}{Mg(Y)}|Y\le y)P(Y\le y)}{1/M}\\ & =M{\int }_{-\infty }^{y}P(U\le \frac{f(Y)}{Mg(y)}|Y=\omega <y)g(\omega )d\omega \\ & =M{\int }_{-\infty }^y\frac{f(\omega )}{Mg(\omega )}g(\omega )d\omega \\ & =M{\int }_{-\infty }^{y}f(\omega )d\omega \\ & =F(y)\end{array}$$

这就证明了接受拒绝法；接受拒绝法似乎更简单，不需要求分布函数的反函数，但它也有如下缺陷：
（1）由于算法要随机地拒绝掉许多建议的随机数y，根据算法效率，我们估计迭代N次后最终会得到随机数的数量大约是N/M；
（2）选择合适的建议概率密度函数g(x)是算法的关键地方，在完全使Mg(x)罩住f(x)的前提下，g(x)的选择原则是：
（a）M要尽可能小；
（b）建议概率密度函数g(x)要容易被抽样，即要求它是很常用的分布或者可以应用其他方法构造出来抽样
（c）在满足前两个要求的基础上，Mg(x)尽可能与f(x)形似，需要削去的部分就越少，这种算法的效率就会越高

需要注意的是：很多时候，人们不采用这种方法的原因几乎都在于它的效率过低。

接受拒绝法抽样举例：

1. 设函数定义如下
   $$g(x)=\frac{8}{7}+\frac{118}{63}{x}^2-\frac{74}{63}{x}^4+\frac{10}{63}{x}^6,x\in [-2,2]$$

试求常数c，使得f(x)= cg(x)是一个概率密度函数，并说明如何从这个密度抽样。

要使得f(x)为一个概率密度函数，那么f(x)在[-2, 2]上的积分必须为1，于是可以求出c的大小，c = 0.1876

g = @(x)(8 / 7 + 118 / 63 * x.^2 -74 / 63 * x.^4 + 10 / 63 * x.^6);
int_g = integral(g, -2, 2);
c = 1 / int_g;
fprintf('c = %f\n', c);

x = -2:0.001:2;
f = c * g(x);
f_max = max(f);

M = f_max / (1 / (2 - (-2)));
fprintf('c1 = %f\n', c1);

N = 100000;

y = unifrnd(-2, 2, 1, N);
u = unifrnd(0, 1, 1, N);
f1 = c * g(y);
h = f1 / M / 0.25;
x = y(u < h);
sample = length(x);
[xNum, xCen] = hist(x, 50);
bar(xCen, xNum / sample);
title('f(x)=8/7 + 118/63*x^2-74/63*x^4+10/63*x^6')
t = -2:0.001:2;
f1 = c * g(t);
scale = 4 / 50;
hold on
plot(t, scale * f1)