r - 在自定义密度下模拟值

Question

我有一个与密度和模拟值有关的理论和编码问题。

我正在通过 density(x) 命令构建自定义密度。但是我希望从这个密度中生成 1000-10000 个模拟值。总体目标是采用密度（x$y）形式构建的两个密度并运行模拟，并说这个密度 A 超过密度 B x% 的时间。我只是取每个模拟值，看看哪个更高，然后编写代码来计算 A 比 B 高多少倍。

有没有办法做到这一点？或者有没有办法用这些密度完成类似的事情？谢谢！

Answer 1

该sample函数可以取样本密度区间的中点，然后使用密度作为概率参数。

mysamp <- sample(x= dens$x, size=1000  , prob=dens$y, repl=TRUE)

这样做的缺点是您可能需要抖动结果以避免大量重复。

 mysamp <- jitter(mysamp)

另一种方法是使用approxfunand ecdf。您可能需要反转函数（x 和 y 的相反角色）以便使用runif(1000)结果中的输入进行采样。我很确定在 SO 中有这样的工作示例，而且我很确定我是过去将此类代码发布到 R-help 的众多人之一。（如果您的搜索未能找到，则发布搜索策略，其他人可以尝试改进它们。）

Answer 2

按照@DWin 的技巧来反转ecdf，这里是如何实现这种方法，使用样条曲线来拟合反转的阶跃函数：

给定

z <- c(rnorm(40), runif(40))
plot(density(z))

定义

spl <- with(environment(ecdf(z)), splinefun(y, x))

sampler <- function(n)spl(runif(n))

现在你可以sampler()用你想要的大小调用：

plot(density(sampler(1000)))

最后说明：这将永远不会生成原始数据范围之外的值，但重复将非常罕见：

> anyDuplicated(sampler(1e4))
[1] 0