原創：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，就是將上述間接方法的過程反過來，以實現一種偽隨機數生成工具。具體來說，就是：

補充證明

最后，我們給出上文的一個結論的證明：

(公眾號：生信了)

總結

以上是生活随笔為你收集整理的曲线均匀分布_R——概率统计与模拟(三) 变换均匀分布对特定分布进行抽样的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。