matlab泊松分布随机数和图像_R——概率统计与模拟(三) 变换均匀分布对特定分布进行抽样...

原创：hxj7

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

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

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

... 等等

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

指数分布抽样

下面的例子是对一个指数分布进行抽样，并绘制了三条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 Havingna 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 实现模拟指数分布抽样是可行的。实际上：

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

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

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 Havingna 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 Havingna 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")

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

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

上面只是对一个均匀分布变量进行变换，我们也可以对多个均匀分布变量进行变换。比如，如果要对标准正态分布抽样，我们可以[0,1]上的两个均匀分布变量X和Y做如下变换：

即可以得到一个标准正态分布的抽样样本，该方法被称作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 Havingna 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.曲线的结果如下：

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

Probability Integral Transformation

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

补充证明

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

(公众号：生信了)