一个简单的解决方案是rejection sampling（尽管在Severin Pappadeux's answer below 上查看我的cmets）。一个通用的拒绝采样算法很容易在 R 中实现（如果它没有在 几个 包中实现，我会感到惊讶，我只是没有查找哪些包具有这样的功能）：
#' Rejection sampling
#'
#' @param f A function calculating the density of interest
#' @param g A function giving the proposal density
#' @param rg A function providing random samples from g
#' @param M A numeric vector of length one giving a bound on the the ratio
#' f(x) / g(x). M must be > 1 and greater than or equal to f(x) / g(x) over
#' the whole support of X.
#' @param n An integer vector of length one giving the number of samples to
#' draw; the default is ten thousand.
#' @param ... Further arguments to be passed to g and rg
rejection_sampling <- function(f, g, rg, M, n = 10000, ...) {
result <- numeric(n)
for ( i in 1:n ) {
reject <- TRUE
while ( reject ) {
y <- rg(n = 1, ...)
u <- runif(1)
if ( u < ( f(y) / (M * g(y, ...)) ) ) {
result[i] <- y
reject <- FALSE
}
}
}
return(result)
}
然后我们可以使用它从您的分布中获取样本，并将样本的密度与真实概率密度一起绘制，看看它的效果如何：
x <- seq(0.01, 15, 0.01)
f <- function(x) 4 * x * exp(-x)
y <- f(x)
set.seed(123)
z <- rejection_sampling(f = f, g = dexp, rg = rexp, M = 10, n = 1e3, rate = 1/4)
dens <- density(z, from = 0.01, to = 15)
scaling_constant <- max(y) / max(dens$y)
plot(x, y, type = "l", xlab = "x", ylab = "f(x)", lty = 2, col = "blue")
lines(dens$x, dens$y * scaling_constant, col = "red", lty = 3)
legend("topright", bty = "n", lty = 2:3, col = c("blue", "red"),
legend = c("True f(x)", "(Re-scaled) density of sample"))
拒绝抽样的工作原理是从提案分布中抽取样本，如果随机均匀偏差大于比率f(x) / M g(x)，则拒绝它们，其中g(x) 是您的提案密度，M 是f(x) / g(x) 的界限如上面 Roxygen 文档中的简要描述。
我在上面使用了速率参数为 1/4 的指数建议分布。你可以使用其他人。
这个 p.d.f.如Severin Pappadeux's answer 所述，与形状为 2、比例为 1 的 Gamma 分布成正比。 （也就是说，如果你将这些参数插入到 Gamma p.d.f. 中，你会发现它与你的只是一个缩放常数不同）。根据您这样做的目的，这可能是更好的方式，或者这种方式。我不确定您的目标是从任意分布中生成样本例如作为示例，或者您需要来自此分布本身的样本等。通常最好识别是否你的任意pdf实际上是一个已实现分布的示例，但如果您不在这种情况下，则需要诸如拒绝抽样之类的东西。