3

我想在已知 CDF 的情况下快速生成离散随机数。本质上，该算法是：

  1. 构造 CDF 向量（从 0 开始到 1 结束的递增向量）cdf
  2. 生成一个 uniform(0, 1) 随机数u
    • 如果u < cdf[1]选1
    • 否则如果u < cdf[2]选择2
    • 否则，如果u < cdf[3]选择 3 *...

例子

首先生成一个cdf：

cdf = cumsum(runif(10000, 0, 0.1))
cdf = cdf/max(cdf)

接下来生成N统一的随机数：

N = 1000
u = runif(N)

现在对值进行采样：

##With some experimenting this seemed to be very quick
##However, with N = 100000 we run out of memory
##N = 10^6 would be a reasonable maximum to cope with
colSums(sapply(u, ">", cdf))
4

3 回答 3

4

如果您知道概率质量函数（如果您知道累积分布函数，您会这样做），您可以使用 R 的内置sample函数，您可以在其中使用参数定义离散事件的概率prob。

cdf = cumsum(runif(10000, 0, 0.1))
cdf = cdf/max(cdf)

system.time(sample(size=1e6,x=1:10000,prob=c(cdf[1],diff(cdf)),replace=TRUE))
   user  system elapsed 
   0.01    0.00    0.02
于 2013-02-28T14:19:52.927 回答
3

如何使用cut：

N <- 1e6
u <- runif(N)
system.time(as.numeric(cut(u,cdf)))
   user  system elapsed 
   1.03    0.03    1.07 

head(table(as.numeric(cut(u,cdf))))

  1   2   3   4   5   6 
 51  95 165 172 148  75
于 2013-02-28T14:42:11.630 回答
2

如果您有有限数量的可能值，那么您可以使用@Hemmo 提到的findIntervalorcut或更好的值。sample

但是，如果您想从理论上达到无穷大的分布（如几何、负二项式、泊松等）生成数据，那么这里有一个可行的算法（这也适用于有限数量的值，如果通缉）：

从您的统一值向量开始，循环遍历从统一向量中减去它们的分布值，随机值是值变为负数的迭代。这是一个更容易看到的例子。这会从平均值为 5 的 Poisson 生成值（将dpois调用替换为您的计算值），并将其与使用逆 CDF 进行比较（在这种情况下它存在的情况下更有效）。

i <- 0
tmp <- tmp2 <- runif(10000)
randvals <- rep(0, length(tmp) )

while( any(tmp > 0) ) {
    tmp <- tmp - dpois(i, 5)
    randvals <- randvals + (tmp > 0)
    i <- i + 1
}

randvals2 <- qpois( tmp2, 5 )

all.equal(randvals, randvals2)
于 2013-03-01T01:34:25.740 回答