R - Probability of a number on multiple variable number of die

Question

I am trying to write an r script of modeling rolling the top number on a die (in this case a 6) on at least one die for a cumulative set of die (1,2,3,4,5,6 dice) on an arbitrary die (in this case a d6; the code should be easily adjustable to other die like d8 d10 d12). The result I am looking for is roughly this for a d6:

Dice (p)All 6s 1 0.16667 2 0.2334 3 0.2927 ...

z <- 1
b <- c(numeric)
while(z < 6) {
  b[z] <- sapply(1:1, function(z) mean(rbinom(10000, z, 1/6) > 0))
  z <- z + 1
}
b

I end up with this:

[1] 0.1631 0.1637 0.1716 0.1623 0.1659 0.1648

I can't get the while loop to step correctly. Any suggestions?

I am able to get 20 terms when I modify the second line of your code to `b <- numeric()`. However, all terms are very close to 1/6 instead of increasing in your example. — Ben Norris, Mar 15 '21 at 22:08
I don't really understand what you're trying to do. Maybe try describing your actual goal in better detail instead of the code because your attempts at describing the code weren't successful in my mind. — Dason, Mar 15 '21 at 22:10
What I mean by that is describe what you care about and give some actual examples. — Dason, Mar 15 '21 at 22:11
I still don't understand what the desired output represents. If you are trying to determine the probability of rolling a 1 on 1d6, 2d6, 3d6, etc., then the results should be 0.166, 0, 0, ... With any other number, the probability of rolling that exact number will decrease with more dice rather than increase. Are you trying to model the probability of rolling *at least* one of a certain number with increasing numbers of dice? — Ben Norris, Mar 15 '21 at 22:32
@EnglishmanBob - Thanks for the clarification. However, that probability should decrease with increasing dice number. For example, the probability of rolling a 6 on 1d6 is 1/6. The probability of rolling 6 and 6 on 2d6 is 1/36. The probability of rolling three sixes on 3d6 is 1/216. — Ben Norris, Mar 15 '21 at 22:52

score 2 · Accepted Answer · answered Mar 15 '21 at 23:33

What I think you want is an approach to systematically simulating the roll of a die an increasing number of times. Let's try a function built on sample:

dice <- function(sides, times){
  sample(1:sides, times, replace = TRUE)
}

If you want to roll 2d6, then you do

set.seed(9) # for reproducibility
dice(6, 2)  
[1] 3 5

Let's say you want to repeat this from 1 to 6 d6s. Now we need sapply. You will get a list with all of the outputs.

set.seed(9)
sapply(1:6, function(z) dice(6, z))

[[1]]
[1] 3

[[2]]
[1] 5 6

[[3]]
[1] 3 3 3

[[4]]
[1] 4 3 6 4

[[5]]
[1] 5 2 5 2 4

[[6]]
[1] 3 4 6 1 6 1

Now you want to check to see if they are all equal to some other value (say 6). You can compare your output and use all.

set.seed(9)
all(dice(6, 2) == 6)
[1] FALSE

Couple this with sapply, and you get a vector for each number of throws.

sapply(1:6, function(z) all(dice(6, z) == 6))
[1] FALSE FALSE FALSE FALSE FALSE FALSE

However, you want to repeat this a large number of times and estimate the number probability of TRUE. Using sapply inside sapply returns a matrix, and we can convert to a data.frame and do some dplyr action to it.

library(dplyr)
set.seed(9)
sapply(1:1000, function(i) sapply(1:6, function(z) all(dice(6, z) == 6))) %>% 
  t %>% data.frame %>% 
  summarise(across(everything(), mean))
     X1    X2    X3    X4 X5 X6
1 0.156 0.017 0.005 0.002  0  0

If you want to wrap this up in a single function so you can choose the sides, number of rolls, target number, and repetitions, you can.

my_fun <- function(sides, times, target, reps, seed = NULL) {
  set.seed = seed
  sapply(1:reps, function(i) sapply(1:times, function(z) all(dice(sides, z) == target))) %>% 
    t %>% data.frame %>% 
    summarise(across(everything(), mean))
}
my_fun(8, 20, 2, 100000) 

       X1      X2      X3      X4    X5 X6 X7 X8 X9 X10 X11 X12 X13 X14 X15 X16 X17 X18 X19 X20
1 0.12556 0.01551 0.00194 0.00023 1e-05  0  0  0  0   0   0   0   0   0   0   0   0   0   0   0

However, without more information from you (e.g. a mathematical formula that generates your sequence), I cannot reproduce your desired output. The probabilities you say you want just do not increase as the number of dice increase.

They are also easy to generate by the sequence of 1/(sides^n).

R - Probability of a number on multiple variable number of die

1 Answers1