Генерация полиномиальных случайных чисел с максимальным размером в каждом поле

Question

В R я хотел бы сгенерировать мультиномиально распределенный вектор случайных чисел заданного размера N, например, используя rmultinom, но с максимальным размером для каждого из K блоков. Например:

set.seed(1)
draw = rmultinom(n = 1, size = 1000, prob = c(67,211,264,166,144,52,2,175))

В этом случае размер равен 1000, указывая общее количество объектов, которые помещаются в восемь полей (length из prob), и prob = c(67,211,264,166,144,52,2,175) вектор вероятностей для восьми ящиков (который внутренне нормализован до суммы 1). Кроме того, я хотел бы, чтобы c(67,211,264,166,144,52,2,175) было вектором максимального размера для каждого из восьми блоков.

Однако в этом случае можно сгенерировать числа, превышающие c(67,211,264,166,144,52,2,175) (для в приведенном выше примере draw[7,]=4 больше, чем 2), тогда как я хотел бы, чтобы каждое число было меньше или равно максимальному размеру каждого поля, указанного в prob, в дополнение к draw, суммируя к size = 1000 .

Знаете ли вы какую-нибудь функцию или какой-нибудь простой способ сделать это? Я не нашел ответа.

Answer 1

Из wikipedia : «Для n независимых испытаний, каждое из которых приводит к успеху ровно для одной из k категорий, при этом каждая категория имеет заданную фиксированную вероятность успеха, полиномиальное распределение дает вероятность любой конкретной комбинации количества успехов для различных категорий ».

Ключевое слово здесь независимый . Ваше ограничение на количество раз, когда каждая категория может быть нарисована, означает, что выборка не является независимой. Если бы ваша проблема была полиномиальной, было бы возможно - хотя и очень маловероятно - что все числа можно было бы взять из поля 7. Это не то, что вам нужно, поэтому вы не можете использовать rmultinom.

Вот другой подход:

# vector of item counts
m <- c(67,211,264,166,144,52,2,175)
# expand the item counts in to a single vector with i repeated m[i] times
d <- unlist(lapply(1:length(m), function(x) rep(x, m[x])))
# sample from d without replacement
s <- sample(d, size=1000, replace=FALSE)
# count the number of items of each type were sampled
table(factor(s))

  1   2   3   4   5   6   7   8 
 63 197 242 153 135  48   2 160