I've been able to model the probability of multiple successes where the probability of each success is the same, but I am unable to figure out how to model the combined probability where the probabilities of individual successes are different.

By example:

If I roll 10 six-sided dice, the probability to roll five fours or more is a result of a cumulative binomial distribution:

Number of trials = n = 10 Probability of success = P = 0.5 Number of successes = x >= 5

Binomial Distribution: b(x; n, P) = nCx * P^x * (1 - P)^(n - x)

and the Cumulative Binomial Distribution:

b(x >= 5; 10, 0.5) = b(x = 0; 10, 0.5) + b(x = 1; 10, 0.5) + .. + b(x = 5; 10, 0.5) = 0.623046875

(source: http://stattrek.com/probability-distributions/binomial.aspx)

So far, so good. My problem now is figuring out how to calculate the probability if the target values are, for example: two fours or greater, two fives or greater and one six out of 10 trials.

Gut feel: The probability for each set of equivalent successes (sixes, fives and fours) should be evaluated in isolation and, once evaluated, the number of trials should be reduced by the size of the set. This also implies that the order makes a difference and for that purpose, I'm assuming that obtaining a six is a higher priority then obtaining a five, and so on. This would then mean that the probability of obtaining a six is calculated first using 10 trials, followed by the fives using 9 trials and the fours using 7 trials.

b(x >= 1; 10, 0.166667) = 0.838495063131323

b(x >= 2; 9, 0.333333) = 0.856932373512168

b(x >= 2; 7, 0.5) = 0.9375

What I don't know is how to put the individual probabilities back together again for a cumulative probability. This method also doesn't feel right, since it seems like you are conduction a total of 26 trials to achieve these results.

I sense I may be over-complicating the calculation since it may just be a normal cumulative binomial function with changing probabilities:

b(x >= 5; 10, {0.166667;0.333333;0.333333;0.5;0.5}) = b(x = 0; 10, 0.166667) + b(x = 1; 10, 0.333333) + b(x = 2; 10, 0.333333) + b(x = 3; 10, 0.5) + b(x = 4; 10, 0.5) + b(x = 5; 10, ???)

As you can see, I have no idea what the probability should be for the last term where x = 5... I'm assuming that the formula should iterate 6 times from 0 to 5, since that is the number of iterations for b(x >= 5; 10, 0.5).

I have not been able to find an equivalent question here and would appreciate any assistance!

  • $\begingroup$ "If I roll 10 six-sided dice, the probability to roll five fours is a result of a cumulative binomial distribution: Number of trials = n = 10 Probability of success = P = 0.5 Number of successes = x >= 5" If the dice are fair, the probability of rolling a 4 is $\frac{1}{6}$, not $\frac{1}{2}$. $\endgroup$
    – Sycorax
    Commented Nov 24, 2017 at 15:34
  • $\begingroup$ Hi, yes, the probability of rolling any single value on the die is 1/6; however the probability of rolling 4 OR MORE on a single die is 3/6 (4, 5 and 6 out of 6). $\endgroup$
    – Justus
    Commented Nov 24, 2017 at 18:49
  • $\begingroup$ That's not what you wrote. $\endgroup$
    – Sycorax
    Commented Nov 24, 2017 at 22:06
  • $\begingroup$ Maybe you misread? But just in case, I'll clarify the problem statement: "Number of trials = n = 10 Probability of success = P = 0.5 Number of successes = x >= 5" - The trials are the number of rolls with a fair die (or for simplicity a number of dice with a single roll). The probability of P = 0.5 refers to die results with a value of 4, 5 and 6 (3/6), which translates to a 0.5 probability. The number of successes x >= 5 refers to the number of dice out of the 10 rolls that must reflect a value of 4, 5 or 6. $\endgroup$
    – Justus
    Commented Nov 26, 2017 at 15:42
  • $\begingroup$ So instead of "the probability to roll five fours" you mean "the probability of five or more rolls of four or more." The text "the probability to roll five fours" means that exactly five dice show 4 (but not 5 or 6). You should clarify that in your question. $\endgroup$
    – Sycorax
    Commented Nov 26, 2017 at 16:07

1 Answer 1


In general, you have to compute these multinomial probabilities individually and add them.

Consider a die with $d$ sides, each appearing with probability $p_i, i=1, 2, \ldots, d$. The experiment consists of rolling this die $n$ times independently and recording the number of times $x_i$ that each face $i$ appears.

Suppose you want to find the chance that the $x_i$ simultaneously are one of the counts in a specified set $S_i$. For instance, the problem asks for the chance that $x_4\in S_4=\{2,3,\ldots,10\}$ (two or more fours), $x_5\in S_5=\{2,3,\ldots,10\}$ (two or more fives), $x_6\in S_6=\{1,2,\ldots,10\}$ (one or more sixes), and $x_1,x_2,x_3$ may have any value in $S_1=S_2=S_3=\{0,1,\ldots, 10\}$.

Let's abstractly name this chance $f(n, \mathbf{p}, \mathbf{S})$ (where $\mathbf p = (p_1,p_2, \ldots, p_d)$ and $\mathbf{S} = (S_1, S_2, \ldots, S_d)$).

The outcomes have a multinomial distribution. Thus

$$f(n, \mathbf{p}, \mathbf{S}) = \sum_{a_i\in S_i; a_1+\cdots+a_d=n} \binom{n}{a_1,a_2,\ldots,a_d}p_1^{a_1} p_2^{a_2} \cdots p_d^{a_d}$$

where the multinomial coefficient can be computed as

$$\binom{n}{a_1,a_2,\ldots,a_d} = \frac{n!}{a_1!a_2!\cdots a_d!}.$$

Evaluate this sum one index at a time. For instance, beginning with the last index $a_d$ we obtain

$$\eqalign{ f(n, \mathbf{p}, \mathbf{S}) &= \sum_{a_d\in S_d} \binom{n}{a_d} p_d^{a_d}\sum_{a_i\in S_i; a_1+\cdots+a_{d-1}=n-a_d} \binom{n-a_d}{a_1,a_2,\ldots,a_{d-1}}p_1^{a_1} p_2^{a_2} \cdots p_{d-1}^{a_{d-1}}\\ &= \sum_{a_d\in S_d} \binom{n}{a_d} p_d^{a_d} f(n-a_d, \mathbf{p}^\prime, \mathbf{S}^\prime) }$$

where the primes indicate the last element has been deleted from each vector. This is a recursive program. It gets started with the case $n=1$, where

$$f(n, p, S) = \left\{\matrix{p^n & \text{if } \in S\\ 0 & \text{otherwise.}}\right.$$

This is potentially a time-consuming calculation, because when the $S_i$ are numerous, there may be a great many cases to examine. For small $d$ and (depending on $d$) small $n$, it's practicable--the example in the question takes no measurable time. For large $d$, $n$, and numerous $S_i$, a cleverer algorithm is needed to screen out cases that contribute so little to the answer that their computation can be ignored.

Because the coefficients $\binom{n}{a_d}$ can be very large and the powers $p_d^{a_d}$ very small, it's a good idea to sum logarithms rather than perform these multiplications directly. That is the method implemented in the code below.

For the example described in the question, the result is $10\,593\,330/6^{10} \approx 0.1751943$. The following R code computed it. It includes a simulation afterwards to run the experiment a million times, and then compares the frequency with which $\mathbf x$ was found in $\mathbf S$ to the value previously computed. A Z-score is also returned: typically, a Z-score less than $2.5$ or so in size is significant confirmation, especially when the simulation is large enough to have observed this event many times. Here is its output:

Theory Simulated         Z 
0.1752    0.1752    0.0255 

The simulation agreed almost exactly with the computed value.

# Return the logarithm of the chance that a die described by probability
# vector `p`, after `n` independent rolls, will exhibit face counts given
# by the numbers in the (nonempty) list `s`.
f <- function(n, p, s) {
  if (length(s)==1) return(ifelse(n %in% s[[1]], n * log(p[1]), -Inf)) # Base case

  a <- s[[1]]
  a <- a[0 <= a & a <= n]
  if (length(a)==0) return(-Inf) # There's nothing satisfying the criteria

  z <-  a * log(p[1]) + lchoose(n, a) + sapply(a, function(a) f(n-a, p[-1], s[-1]))
# Examples and testing.
n <- 10
d <- 6
p <- rep(1/d, d)
s <- list(0:n, 0:n, 0:n, 2:n, 2:n, 1:n)
# s <- lapply(p, function(i) sample(0:n, ceiling(0.9*n))) # For testing
theory <- exp(f(n, p, s))
# Simulation.
n.sample <- 1e6
x <- rmultinom(n.sample, n, p)
e <- paste(sapply(1:length(s), function(i) paste0("x[", i, ", ] %in% s[[", i, "]]")), 
           collapse=" & ")
s <- eval(parse(text=paste0("mean(", e, ")")))
round(c(Theory=theory, Simulated=s, 
        Z=(s-theory) / sqrt(theory*(1-theory)/n.sample)), ceiling(log10(n.sample)/2+1))

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.