r - 模拟来自(非标准)密度函数的数据

Question

我想模拟来自非标准密度函数的数据。我已经找到以下链接( How do I best simulate an arbitrary univariate random variate using its probability function? )。然而，这会产生奇怪的结果。不知何故，这个累积密度函数( cdf() )效果不佳。从某些值来看，它给出了非常奇怪的结果。例如，看一下下面的代码:

density=function(x)(25*200.7341^25/x^26*exp(-(200.7341/x)^25))
cdf<-function(x) integrate(density,1,x)[[1]]

cdf(9701)
[1] 1

cdf(9702)
[1] 6.33897e-05

所以我的问题是，如何创建一个“好的”CDF 函数？或者更直接地说，如何模拟 PDF 中的数据？

Answer 1

正如 @pjs 所指出的，我们可以使用 Rejection sampling (查看 wiki 了解详细信息)。

这是此方法的一种实现。

最重要的一步是找到一个分布 g，我们可以从中采样，并从中存在 M，使得所有点的 M * g > f

f <- function(x) (25 * 200.7341^25 / x^26 * exp(-(200.7341/x)^25))
g <- function(x) dnorm(x, mean = 200.7341, sd = 40)
M <- 5
curve(f, 0, 500)
curve(M * g(x), 0, 500, add = TRUE, lty = "dashed")

现在，我们可以执行算法

set.seed(42)
k <- 1
count <- 0
res <- vector(mode = "numeric", length = 1000)
while(k < 1001) {
          z <- rnorm(n = 1, mean = 200.7341, sd = 40)
          R <- f(z) / (M * g(z))
          if (R > runif(1)) {
              res[k] <- z
              k  <- k + 1
          }
          count <- count + 1
    }

(accept_rate <- (k / count) * 100)
## [1] 19.7086

require(MASS) ## for truehist
truehist(res)
curve(f, 0, 250, add = TRUE)

录取率不是很高。您可以尝试寻找更好的包络函数或使用 Metropolis Hasting 算法。