R, использовать биномиальное распределение с более чем двумя возможностями

Question

Я знаю, что это, вероятно, элементарно, но у меня, кажется, есть умственный блок. Допустим, вы хотите вычислить вероятность броска 4, 5 или 6 на броске одного кубика. В R это достаточно просто:

sum(1/6, 1/6, 1/6)

Это дает 1/2, который является правильным ответом. Тем не менее, я в глубине души (где это возможно должно остаться), что я должен быть в состоянии использовать биномиальное распределение для этого. Я пробовал различные комбинации аргументов для pbinom и dbinom, но не могу получить правильный ответ.

С бросками монет это работает отлично. Это совершенно неуместно в ситуациях, когда есть более двух возможных результатов? (Я программист, а не статистик, поэтому я ожидаю, что меня убьют статисты.)

Вопрос: Как я могу использовать pbinom() или dbinom() для вычисления вероятности броска 4, 5 или 6 с одним броском кубика? Я знаком с пакетами prob и dice, но я действительно хочу использовать один из встроенных дистрибутивов.

Благодарю.

Answer 1

Как упоминалось выше @Alex, бросание игральных костей может быть представлено в терминах полиномиальных вероятностей. Например, вероятность броска 4

dmultinom(c(0, 0, 0, 1, 0, 0), size = 1, prob = rep(1/6, 6)) 
# [1] 0.1666667

и вероятность броска 4, 5 или 6 равна

X <- cbind(matrix(rep(0, 9), nc = 3), diag(1, 3))
X
#      [,1] [,2] [,3] [,4] [,5] [,6]
# [1,]    0    0    0    1    0    0
# [2,]    0    0    0    0    1    0
# [3,]    0    0    0    0    0    1
sum(apply(X, MAR = 1, dmultinom, size = 1, prob = rep(1/6, 6)))
# [1] 0.5

Answer 2

Хотя это не совсем очевидно, это можно сделать с помощью pmultinom, реализованного либо в моем пакете pmultinom на CRAN, либо на этом другом пакете pmultinom на Github.

Вы понимаете это как событие, которое не является 1, 2 или 3. Затем вы пишете эту вероятность как

P(X_1 ≤ 0, X_2 ≤ 0, X_3 ≤ 0, X_4 ≤ ∞, X_5 ≤ ∞, X_6 ≤ ∞)

где X_i - количество вхождений стороны i. Все X вместе имеют полиномиальное распределение с параметром размера 1, и все вероятности равны 1/6. Эта вероятность выше может быть рассчитана (используя мой пакет) как

pmultinom(upper=c(0, 0, 0, Inf, Inf, Inf), size=1,
          probs=c(1/6, 1/6, 1/6, 1/6, 1/6, 1/6), method="exact")
# [1] 0.5

Хотя это немного неловко, но мне нравится, потому что я предпочитаю использовать функцию "p", а не сумму "d".