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

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

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

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

$\dots$ 等等

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

指数分布抽样

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

• 用 runif 函数模拟指数分布抽样的 $\text{c.d.f.}$ 曲线；
• 用R自带的 rexp 实现指数分布抽样的 $\text{c.d.f.}$ 曲线；
• 指数分布的理论 $\text{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 实现模拟指数分布抽样是可行的。实际上：

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