2013-02-28 179 views
3

我想快速生成离散随机数,我有一个已知的CDF。本质上,该算法是:高效地生成离散随机数

  1. 构建CDF矢量(0,1)随机数u
    • 如果u < cdf[1]选择(从0开始以1增加矢量和结束)cdf
    • 产生均匀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)) 

回答

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 
4

如果你知道概率密度函数(你做什么,如果你知道的累积分布函数),您均可以使用内置的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 
+0

而作为“如果替换为真,则使用沃克的别名法(里普利,1987年)时,有超过250个合理可能的值”,它是有效的时间复杂度是O(n)的 – colinfang 2013-11-21 14:52:21

2

如果有可能的值的数量有限,那么你可以使用findIntervalcut或更好sample由@Hemmo提及。然而,如果你想从理论上走向无穷大(如几何,负二项式,泊松等)的分布生成数据,那么这里是一个算法,它将起作用(这也将与有限的如果需要值的数量):

从您的统一值向量开始,循环遍历分布值,然后从统一向量中减去它们,随机值是值变为负值的迭代。这是一个更容易看到的例子。这将生成平均值为5的泊松(将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) 
+0

大约分布好点无限的支持,不知何故我忘了那些。 – 2013-03-01 04:10:26

+0

这正是我的问题。但是,如果写入的算法在R中会有可怕的性能。目前,我使用大量的“i”步骤,我想我会使用'cut'来生成随机数。 – csgillespie 2013-03-03 22:35:55