Search code examples
rmontecarlodice

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


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?


Solution

  • 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).