# Dice probability - amount of “hits” in a pool

What I want is to calculate the probability of a minimum number of rolls with a minimum value for at given number of dice, e.g. "if I roll 20 dice, what is the probability of getting at least 3 dice with a value of 5 or more?"

I have a brute force setup, that creates the Cartesian product and just enumerates through that, but it completely explodes at 9 dice. So what I'm really looking for is some kind of function that does not require me to expand all possibilities, so that I can calculate for 20-30 dice.

"if I roll 20 dice, what is the probability of getting at least 3 dice with a value of 5 or more".

Break it down into two stages:

a) if I roll 1 die, what is the probability of getting a value of 5 or more?

I assume these are 6-sided dice numbered 1-6. In which case $p=\frac{1}{3}$

You seem to call getting that result 'a hit'. I'll refer to it as "a success".

b) In 20 trials, what is the probability of at least 3 successes?

This is a binomial distribution question -

Given $X\sim \text{Bin}(n=20, p=\frac{1}{3})$, what is $P(X\geq 3)$?

By hand, that probability is relatively easily computed as $1-[P(X=0)+P(X=1)+P(X=2)]$

$=1-[(1-1/3)^{20}+ {{20}\choose{1}}(1/3)^1(1-1/3)^{19}+{{20}\choose{2}}(1/3)^2(1-1/3)^{18}]$

With bigger numbers of values to compute, this rapidly gets tedious.

These values can be obtained by tables or by computer programs (like R or Excel, for example). So let's look at some examples:

In R:

Calculated as if by hand:

> 1-((2/3)^20+20*(1/3)*(2/3)^19+20*19/2*(1/3)^2*(2/3)^18)
[1] 0.9824074


calculated using the inbuilt cdf:

> 1-pbinom(2,20,1/3)
[1] 0.9824074


So the probability of 3 or more successes in 20 dice is 0.9824

In Excel:

=1-BINOM.DIST(2,20,1/3,TRUE)


gives:

0.982407


(in older versions of Excel use BINOMDIST - with the same arguments - instead)

So how do we write it just as a single function?

In R you could do this:

poolprob <- function(target, dice, minhits) 1-pbinom(minhits-1,dice,1-(target-1)/6)


then:

> poolprob(5, 20, 3)
[1] 0.9824074


(there are one or two improvements that could be made, but this is just for illustration)

You could similarly set up an Excel spreadsheet that did the same kind of calculation.

So here it is for 3-30 dice:

> n=3:30
> data.frame(n,p=poolprob(5, n, 3))


which gives:

    n          p
1   3 0.03703704
2   4 0.11111111
3   5 0.20987654
4   6 0.31961591
5   7 0.42935528
6   8 0.53177869
7   9 0.62282172
8  10 0.70085861
9  11 0.76588935
10 12 0.81887735
11 13 0.86126776
12 14 0.89466626
13 15 0.92064288
14 16 0.94062489
15 17 0.95584927
16 18 0.96735214
17 19 0.97597930
18 20 0.98240737
19 21 0.98716891
20 22 0.99067741
21 23 0.99325031
22 24 0.99512894
23 25 0.99649521
24 26 0.99748526
25 27 0.99820030
26 28 0.99871513
27 29 0.99908475
28 30 0.99934942

• That looks great! I'll try and see if i can implement it :) – Brestfloda Jun 12 '14 at 5:25
• You can compare the results with complete enumeration for the cases you have both; it would be a handy check you have them both implemented correctly. – Glen_b Jun 12 '14 at 8:21