R-概率统计与模拟（三）变换均匀分布对特定分布进行抽样

本文介绍了如何变换均匀分布以便对特定分布进行抽样。

如果你要进行随机抽样，R语言提供了诸多现成的函数供你使用，比如：

runif: 均匀分布抽样
rbinom: 二项分布抽样
rpois: 泊松分布抽样
rnorm: 正态分布抽样
rexp: 指数分布抽样
rgamma: 伽马分布抽样

…\ldots… 等等

那么，如果不用现成的函数，我们能自己实现抽样功能吗？比如，我们是否可以不用 rexp 函数而实现指数分布抽样？答案是肯定的，只需对均匀分布进行适当地变换即可。

指数分布抽样

下面的例子是对一个指数分布进行抽样，并绘制了三条 c.d.f.\text{c.d.f.}c.d.f. 曲线，分别是：

用 runif 函数模拟指数分布抽样的 c.d.f.\text{c.d.f.}c.d.f. 曲线；
用R自带的 rexp 实现指数分布抽样的 c.d.f.\text{c.d.f.}c.d.f. 曲线；
指数分布的理论 c.d.f.\text{c.d.f.}c.d.f. 曲线。

代码如下：

# Q1
N <- 100000
lambda <- 1   # 指数分布参数
set.seed(123)
x <- runif(N, 0, 1)
y <- log(1 - x) / -lambda   # 用runif模拟指数分布抽样
ey <- ecdf(y)
set.seed(123)
r <- rexp(N, lambda)     # R自带的rexp抽样函数
er <- ecdf(r)
plot(ey, xlim = c(0,3), col="red", main="Generating Pseudo-Random Numbers Having\na Exponential Distribution with lambda=1",ylab="Cum prob F(x)")   # runif模拟抽样的c.d.f.曲线
lines(er, xlim=c(0,3), col="green")    # rexp模拟的c.d.f.曲线
curve(1 - exp(-lambda * x), from=0, to=3, add=T, col="blue")  # 指数分布的理论c.d.f.曲线
legend("bottomright", legend=c("simulative (using runif)", "simulative (using rexp)", "theoretical"),col=c("red", "green", "blue"), lty=1, bty="n")

效果如下：

可以看出，三条曲线几乎完全重合，说明用 runif 实现模拟指数分布抽样是可行的。实际上：

假设我们想要根据某个特定的概率分布（我们已经知道该分布的 c.d.f.\text{c.d.f.}c.d.f. 是FFF）进行抽样，可以这样：首先我们对一个在[0,1][0,1][0,1]上服从均匀分布的变量XXX进行随机抽样，得到x1,x2,…,xnx_1, x_2, \ldots, x_nx1,x2,…,xn；然后我们构造一个变量Y=F−1(X)Y=F^{-1}(X)Y=F−1(X)，并据此计算得到 y1,y2,…,yny_1, y_2, \ldots, y_ny1,y2,…,yn。那么，y1,y2,…,yny_1, y_2, \ldots, y_ny1,y2,…,yn就是我们想要的抽样样本。因为实际上可以证明 YYY 的 c.d.f.\text{c.d.f.}c.d.f. 就是FFF（证明过程在文末）。
上面F−1F^{-1}F−1其实就是 quantile 函数，值得注意的是当有很多个值F(x1)=F(x2)=⋯=F(xn)=yF(x_1)=F(x_2)=\cdots=F(x_n)=yF(x1)=F(x2)=⋯=F(xn)=y时，F−1(y)=min⁡{x1,x2,…,xn}F^{-1}(y)=\min\{x_1,x_2,\ldots,x_n\}F−1(y)=min{x1,x2,…,xn}。

也就是说，我们可以通过一种间接的方法，即变换均匀分布来对特定分布进行抽样。为什么要这样拐弯抹角呢？因为有时候我们碰到的分布不是很常见，R语言并没有提供相应的函数，这时候我们就可以用上述间接的方法自己写函数实现抽样。

变换均匀分布对一种特定分布抽样

比如，如果我们想要对下面这一个分布进行抽样，其p.d.f.\text{p.d.f.}p.d.f.是：
f(x)={2x0≤x≤1,0otherwise.f(x) = \begin{cases} 2x \quad & 0 \le x \le 1, \\ 0 & \text{otherwise.} \end{cases}f(x)={2x00≤x≤1,otherwise.

R语言并没有提供一个现成的函数可以对上面的分布进行抽样。但是我们可以自行编写代码实现：

首先我们知道上面分布的 c.d.f.\text{c.d.f.}c.d.f.是：
F(x)={0x<0,x20≤x≤1,1x>1.F(x) = \begin{cases} 0 \quad & x < 0, \\ x^2 \quad & 0 \le x \le 1, \\ 1 & x > 1. \end{cases}F(x)=⎩⎪⎨⎪⎧0x21x<0,0≤x≤1,x>1.
那么：
F−1(x)=x12for 0≤x≤1.F^{-1}(x) = x^{\frac{1}{2}} \quad \text{for $0 \le x \le 1.$} F−1(x)=x21for 0≤x≤1.
按照上面的间接方法，我们首先对[0,1][0, 1][0,1]上的均匀分布进行抽样，得到一系列 xxx 值：

N <- 100000     # 抽样的样本大小
set.seed(123)
x <- runif(N, 0, 1)

然后根据Y=F−1(X)Y=F^{-1}(X)Y=F−1(X)计算出相应的 yyy 值：

y <- sqrt(x)    # 用runif模拟实现特定分布抽样

计算出的一系列 yyy 值就构成了我们想要的抽样样本。我们可以比较一下用这种间接方法得出的模拟 c.d.f.\text{c.d.f.}c.d.f. 曲线和理论 c.d.f.\text{c.d.f.}c.d.f. 曲线相差多少：

ey <- ecdf(y)
plot(ey, xlim = c(0,1), col="red", main="Generating Pseudo-Random Numbers Having\na Specified Distribution",ylab="Cum prob F(x)")   # 模拟抽样的c.d.f.曲线
curve(x ^ 2, from=0, to=1, add=T, col="blue")  # 上述特定分布的理论曲线
legend("bottomright", legend=c("simulative", "theoretical"),col=c("red", "blue"), lty=1, bty="n")

结果如下：

从图中可以看出，模拟曲线和理论曲线几乎完全重合，也就是说我们的间接方法的确很好地模拟了特定分布抽样。

上述特定分布的完整代码如下：

N <- 100000     # 抽样的样本大小
set.seed(123)
x <- runif(N, 0, 1)
y <- sqrt(x)    # 用runif模拟实现特定分布抽样
ey <- ecdf(y)
plot(ey, xlim = c(0,1), col="red", main="Generating Pseudo-Random Numbers Having\na Specified Distribution",ylab="Cum prob F(x)")   # 模拟抽样的c.d.f.曲线
curve(x ^ 2, from=0, to=1, add=T, col="blue")  # 上述特定分布的理论曲线
legend("bottomright", legend=c("simulative", "theoretical"),col=c("red", "blue"), lty=1, bty="n")

我们用两个例子说明了可以用一种间接抽样的方法对特定分布进行抽样，不过这种方法是有前提的，即我们要知道FFF的解析表达式，以及FFF要存在一个反函数。由于这些限制，该方法在实际应用中有诸多限制。

Box-Muller方法：变换均匀分布对标准正态分布抽样

上面只是对一个均匀分布变量进行变换，我们也可以对多个均匀分布变量进行变换。比如，如果要对标准正态分布抽样，我们可以对[0,1][0,1][0,1]上的两个均匀分布变量XXX和YYY做如下变换：
Z=cos⁡(2πX)log⁡(1Y2)Z = \cos(2\pi X)\sqrt{\log(\frac{1}{Y^2})}Z=cos(2πX)log(Y21)

即可以得到一个标准正态分布的抽样样本，该方法被称作Box-Muller方法。

具体代码如下：

N <- 100000     # 抽样的样本大小
set.seed(123)
x <- runif(N, 0, 1)
y <- runif(N, 0, 1)
z <- cos(2 * pi * x) * sqrt(log(1 / y ^ 2))    # 用runif模拟实现标准正态分布抽样
ez <- ecdf(z)
set.seed(123)
r <- rnorm(N, 0, 1)     # 用rnorm模拟实现标准正态分布抽样
er <- ecdf(r)
plot(ez, col="red", main="Generating Pseudo-Random Numbers Having\na Standard Normal Distribution",ylab="Cum prob F(x)")   # runif模拟抽样的c.d.f.曲线
lines(er, col="blue")   # 用rnorm模拟抽样的c.d.f.曲线
legend("bottomright", legend=c("simulative (using runif)", "simulative (using rnorm)"),col=c("red", "blue"), lty=1, bty="n")

c.d.f.\text{c.d.f.}c.d.f. 曲线的结果如下：

可以看出，两条曲线几乎完全重叠，说明变换的效果是成功的。

Probability Integral Transformation

最后，我们简单介绍一下Probability Integral Transformation，就是将上述间接方法的过程反过来，以实现一种伪随机数生成工具。具体来说，就是：

如果一个连续型变量XXX的c.d.f.\text{c.d.f.}c.d.f.是FFF，那么变量Y=F(x)Y=F(x)Y=F(x)在[0,1][0,1][0,1]上服从均匀分布。

补充证明

最后，我们给出上文的一个结论的证明：

假设我们想要根据某个特定的概率分布（我们已经知道该分布的 c.d.f.\text{c.d.f.}c.d.f. 是FFF）进行抽样，可以这样：首先我们对一个在[0,1][0,1][0,1]上服从均匀分布的变量XXX进行随机抽样，得到x1,x2,…,xnx_1, x_2, \ldots, x_nx1,x2,…,xn；然后我们构造一个变量Y=F−1(X)Y=F^{-1}(X)Y=F−1(X)，并据此计算得到 y1,y2,…,yny_1, y_2, \ldots, y_ny1,y2,…,yn。那么，y1,y2,…,yny_1, y_2, \ldots, y_ny1,y2,…,yn就是我们想要的抽样样本。因为实际上可以证明 YYY 的 c.d.f.\text{c.d.f.}c.d.f. 就是FFF。
上面F−1F^{-1}F−1其实就是 quantile 函数，值得注意的是当有很多个值F(x1)=F(x2)=⋯=F(xn)=yF(x_1)=F(x_2)=\cdots=F(x_n)=yF(x1)=F(x2)=⋯=F(xn)=y时，F−1(y)=min⁡{x1,x2,…,xn}F^{-1}(y)=\min\{x_1,x_2,\ldots,x_n\}F−1(y)=min{x1,x2,…,xn}。

证明：
因为XXX在[0,1][0,1][0,1]上服从均匀分布，所以其c.d.f.\text{c.d.f.}c.d.f. 是：
G(x)=Pr⁡(X≤x)=xG(x) = \Pr(X \le x) = x G(x)=Pr(X≤x)=x

因此，YYY 的 c.d.f.\text{c.d.f.}c.d.f. 是：
H(y)=Pr⁡(Y≤y)=Pr⁡[F−1(X)≤y]=Pr⁡[F(F−1(X))≤F(y)]=Pr⁡[X≤F(y)]=G(F(y))=F(y)\begin{aligned} H(y) & = \Pr(Y \le y) \\ & = \Pr[F^{-1}(X) \le y] \\ & = \Pr[F(F^{-1}(X)) \le F(y)] \\ & = \Pr[X \le F(y)] \\ & = G(F(y)) \\ & = F(y) \end{aligned}H(y)=Pr(Y≤y)=Pr[F−1(X)≤y]=Pr[F(F−1(X))≤F(y)]=Pr[X≤F(y)]=G(F(y))=F(y)

注意，从第二个等式到第三个等式成立是因为FFF是一个c.d.f.\text{c.d.f.}c.d.f.，所以FFF是一个单调不减函数，再加上F−1F^{-1}F−1 quantile 函数的特点，因此两个等式是等价的。

也就是说，YYY 的 c.d.f.\text{c.d.f.}c.d.f. 就是FFF。证毕。

（公众号：生信了）