Нарисуйте огромный образец (1e09) из полиномиального распределения с образцом - PullRequest
0 голосов
/ 27 октября 2011

Я бы хотел сделать выборку из полиномиального распределения. Я бы сделал это, используя образец и указав некоторые вероятности. Например: у меня есть 3 категории, и я хочу сделать выборку 10 раз.

> my_prob = c(0.2, 0.3, 0.5)
> x = sample(c(0:2), 100, replace = T, prob = my_prob)
> head(x)
[1] 2 0 2 1 1 2

Мои настройки теперь отличаются только в следующем аспекте: я хочу сэмплировать много (например, 1e09) чисел. И на самом деле меня интересует только частота каждой категории. Таким образом, в приведенном выше примере это будет означать:

> table(x)
x
 0  1  2 
27 29 44 

У кого-нибудь есть идеи, как рассчитать это максимально эффективно?

спасибо, Штеффи

Ответы [ 2 ]

6 голосов
/ 27 октября 2011

Вам нужно rmultinom.

my_prob <- c(0.2,0.3,0.5)
number_of_experiments <- 10
number_of_samples <- 100
experiments <- rmultinom(n=number_of_experiments, size=number_of_samples, prob=my_prob)
experiments

     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
     [1,]   14   18   15   19   14   17   23   18   24    15
     [2,]   33   34   36   30   40   30   27   38   24    30
     [3,]   53   48   49   51   46   53   50   44   52    55
1 голос
/ 27 октября 2011

Если проблема в том, что вы не можете поместить вектор длины 1e9 в ОЗУ, то вы можете многократно рассчитать таблицу для меньшего числа выборок и сложить итоговые значения.

n_total <- 1e9
n_chunk <- 1e6
n_iter <- n_total / n_chunk
my_prob = c(0.2, 0.3, 0.5)
totals <- numeric(3)
for(i in seq_len(n_iter))
{
  totals <- totals + table(sample(0:2, n_chunk, replace = TRUE, prob = my_prob))
}
totals
stopifnot(sum(totals) == n_total)

Как и Макс сказал , вы можете предпочесть rmultinom сэмплу. Возьмите rowSums его experiments переменной.

...