3

既知の CDF がある離散乱数をすばやく生成したいと考えています。基本的に、アルゴリズムは次のとおりです。

  1. CDF ベクトルを作成します (0 から始まり 1 で終わる増加するベクトル)cdf
  2. 一様 (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 で述べられているように、findIntervalまたはcutそれ以上を使用できます。sample

ただし、理論的に無限大になる分布 (幾何学、負の二項、ポアソンなど) からデータを生成する場合は、次のアルゴリズムが機能します (これは、次の場合に有限数の値でも機能します)。希望):

均一値のベクトルから始めて、分散値をループして均一値のベクトルから減算します。ランダム値は、値が負になる反復です。これは例で見るのが簡単です。これは、平均 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)
于 2013-03-01T01:34:25.740 に答える