# Using R for dice probabilities

Using R, what function(s) would I use to obtain the following probabilities?

• Roll at least one 1 when rolling 2 six-sided dice (2d6) = 11/36
• Roll at least one 1 when rolling 3 six-sided dice (3d6) = 91/216
• Roll at least one 1 when rolling 1d4, 1d6, 1d8, and 1d8 = 801/1536

First I hope my answers above are correct! I did these pretty much manually.

I think I need to use binomial distributions and/or probability-generating functions, but not sure if I'm over-complicating things. I've tried using R's *binom() functions but can't seem to arrive at the answers I need.

• The last answer of $654/1536$ is too low. For a brute force (but illuminating) calculation, try dice <- expand.grid(1:4, 1:6, 1:8, 1:8); dim(subset(dice, subset=(Var1==1 | Var2==1 | Var3==1 | Var4==1)))[1] and inspect the contents of dice afterwards: it shows all $1536$ possible outcomes.
– whuber
Commented Mar 24, 2013 at 22:29
• Made the edit, 801 / 1536 is the correct probability for the last question. Commented Mar 25, 2013 at 9:26
• Thanks for fixing that last one. My main goal here is to arrive at these probabilities without having to load a potentially very large table into memory. Also I want to use standard statistical methods instead of a package that abstracts it all for you, e.g. the 'dice' package. But I do thank @djhurio for providing that info.
– Tim
Commented Mar 25, 2013 at 16:14
• You can do dice$containsOne <- apply(dice, 1, function(roll) 1 %in% roll), then use dplyr %>% pipe or data.table := directly into a summarize. You could save memory by making each row-index a unique string summarizing the rolls e.g. "3287" or "3:2:8:7". (for dice with <- 10 sides, you could directly use a decimal integer like 3287 as row-index, save more memory) – smci Commented Sep 17, 2018 at 0:05 ## 3 Answers Solution 1 There is a package dice that could help you. I got the same probabilities for the first two examples, but not for the last one. See the example code: require(dice) # Roll at least one 1 when rolling 2 six-sided dice (2d6) = 11/36 getEventProb(nrolls = 2, ndicePerRoll = 1, nsidesPerDie = 6, eventList = list(1)) 11/36 # Roll at least one 1 when rolling 3 six-sided dice (3d6) = 91/216 getEventProb(nrolls = 3, ndicePerRoll = 1, nsidesPerDie = 6, eventList = list(1)) 91/216 # Roll at least one 1 when rolling 1d4, 1d6, 1d8, and 1d8 = 654/1536 p1 <- getEventProb(nrolls = 1, ndicePerRoll = 1, nsidesPerDie = 4, eventList = list(1)) p2 <- getEventProb(nrolls = 1, ndicePerRoll = 1, nsidesPerDie = 6, eventList = list(1)) p3 <- getEventProb(nrolls = 2, ndicePerRoll = 1, nsidesPerDie = 8, eventList = list(1)) p1 + p2 + p3 - p1 * p2 - p2 * p3 - p3 * p1 + p1 * p2 * p3 801/1536 ### The same as 1 - (1-p1) * (1-p2) * (1-p3) 801/1536  Solution 2 This is a solution with out usage of any package. You can compute the probability to draw at least one$1$by this formula (mentioned by @whuber): $$p = 1 - \prod_{i=1}^n \left( 1 - \frac{1}{d_i} \right)$$ where$n$is the number of dices and$d_i$is the number of sides of dice$i$. Then you can define a function in R with one argument dices, where dices is a vector of sides. See the example code: dice.prob <- function(dices) 1 - prod(1 - 1 / dices) dice.prob(c(6, 6)) 11/36 dice.prob(c(6, 6, 6)) 91/216 dice.prob(c(4, 6, 8, 8)) 801/1536  • What is the basis of your last calculation? – whuber Commented Mar 24, 2013 at 22:25 • @whuber, it is just another way hot to calculate the probability of three independent events. Added alternative to the answer. Commented Mar 25, 2013 at 6:50 • It looks better now, but how do you wind up with an incorrect answer? The full formula is$1 - (1-1/4)(1-1/6)(1-1/8)(1-1/8)$=$801/1536\$.
– whuber
Commented Mar 25, 2013 at 13:13
• @whuber, both solutions are equivalent. The wrong answer was taken from the question. Made a correction to the answer by deleting the wrong answer. Commented Mar 25, 2013 at 14:38
• @djhurio thanks but I'm more interested in knowing the methods behind the scenes. So I want to know how the dice package arrived at these answers.
– Tim
Commented Mar 25, 2013 at 16:54

These are arithmetically very simple problems. Therefore, it's not clear in what way you intend to use R as a calculator to solve your problems for you. Focusing on example 3 as the hardest example: For U=d4, X=1d6, Y=1d8, and Z=1d8:

Pr(U=1 or X=1 or Y=1 or Z=1) = 1 - Pr(U!=1 and X!=1 and Y != 1 and Z!=1)

This is 1- 3/4 * 5/6 * 7/8 * 7/8 = 801/1536

You can of course plug in 1-3/4 * 5/6 * 7/8 * 7/8  to your terminal and obtain the probability. tada.

You could simulate:

set.seed(123)
n <- 1e8
u = sample(1:4, n, replace=T)
x = sample(1:6, n, replace=T)
y = sample(1:8, n, replace=T)
z = sample(1:8, n, replace=T)

mean(u==1 | x==1 | y==1 | z==1)


Gives 0.5215

Other than that, the existence of and the need for a package to do thinking on a probability problem like this implies a lack of thought.

nothrows <- seq(6, by = 6, length.out = 200)

probsix <- 0
for(i in 1: length(nothrows)) {
samp <- sample.int(6, nothrows[i], replace = TRUE)
probsix[i] <- sum(samp == 6)/nothrows[i]
}
probsix


I would do something like this and just change the values manually as required.