高效地生成離散隨機數

Question

我想快速生成離散隨機數，我有一個已知的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)) 

Answer 1

如何使用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 

Answer 2

如果你知道概率密度函數（你做什麼，如果你知道的累積分佈函數），您均可以使用內置的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 

Answer 3

如果有可能的值的數量有限，那麼你可以使用findInterval或cut或更好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)