r语言 - 从密度对象生成随机数(或者更广泛地说从一组数字)

Question

假设我有一组数字，我怀疑它们来自同一个分布。

set.seed(20130613)
x <- rcauchy(10)

我想要一个函数，从相同的未知分布随机生成一个数字。我想到的一种方法是创建一个density对象，然后从中获得CDF，并取随机均匀变量的逆CDF(参见Wikipedia)。

den <- density(x)
#' Generate n random numbers from density() object
#' 
#' @param n The total random numbers to generate
#' @param den The density object from which to generate random numbers
rden <- function(n, den)
{
        diffs <- diff(den$x)
        # Making sure we have equal increments
        stopifnot(all(abs(diff(den$x) - mean(diff(den$x))) < 1e-9))
        total <- sum(den$y)
        den$y <- den$y / total
        ydistr <- cumsum(den$y)
        yunif <- runif(n)
        indices <- sapply(yunif, function(y) min(which(ydistr > y)))
        x <- den$x[indices]
        return(x)
}
rden(1, den)
## [1] -0.1854121

我的问题如下:

1. 是否有更好的(或内置R)方式来从密度对象生成随机数?
2. 关于如何从一组数字生成随机数(除了sample)还有其他想法吗?

Answer 1

要从密度估计中生成数据，只需随机选择一个原始数据点，并根据密度估计的内核添加一个随机的"误差"块，对于默认的"高斯"，这只是意味着从原始向量中选择一个随机元素，并添加一个随机正态，平均值为0,sd等于所使用的带宽:

den <- density(x)
N <- 1000
newx <- sample(x, N, replace=TRUE) + rnorm(N, 0, den$bw)

另一种选择是使用logspline包中的logspline函数拟合密度(使用不同的估计密度的方法)，然后使用该包中的rlogspline函数从估计的密度生成新数据。

Answer 2

如果您所需要的只是从现有的数字池中绘制值，那么sample就是您的选择。
如果你想从假定的底层分布中绘制，那么使用density，并将其拟合到假定的分布中以获得必要的系数(mean, sd等)，并使用适当的R分布函数。

除此之外，我会看一看C中数值食谱的第7.3章("拒绝方法")，以了解根据任何分布"选择性"抽样的方法。代码非常简单，可以很容易地翻译成R。我敢打赌已经有人这样做了，并且会发布一个比这个更好的答案。

Answer 3

Greg Snow的回答对我很有帮助，我意识到密度函数的输出包含了从输入分布中创建随机数所需的所有数据。基于他的示例，您可以执行以下操作来使用密度输出获得随机值。

x <- rnorm(100)   # or any numeric starting vector you desire
dens <- density(x)
N <- 1000
newx <- sample(x = dens$x, N, prob = dens$y, replace=TRUE) + rnorm(N, 0, dens$bw)

甚至可以创建一个简单的随机数生成函数

rdensity <- function(n, dens) {
    return(sample(x = dens$x, n, prob = dens$y, replace=TRUE) + rnorm(n, 0, dens$bw))
}