CODEWITHSUNDEEP
r
AKUGIZIBWE EDWIN
AKUGIZIBWE EDWIN

Reputation: 11

r code for rolling n dice with elimination

I am trying to write an r code for this challenge: Suppose we roll n dice, remove all the dice that come up 1, and roll the rest again. If we repeat this process, eventually all the dice will be eliminated. How many rolls, on average, will we make? Here is what I have tried so far but it is not working. The text book theoretical answer for 5 dice is 13.02

CODE ATTEMPT

N=10000
myfun <- function(...) {
  for(i in list(...)){
  num=1
  S=sample(1:6,i,replace = TRUE)
  i=i-length(which(S==1))
  while(i!=0){
    i=i-length(which(S==1))
    num=num+1
  }
  result[i]=num
}

}

replicate(N,myfun(1:100))

Upvotes: 0

Views: 312

Answers (3)

AKUGIZIBWE EDWIN
AKUGIZIBWE EDWIN

Reputation: 11

Here is the updated code for the above problem

N=10000 # number of simulation
    set.seed(1873)
    NoOfRolls<-function(d){
      num=0
      while(d!=0){
        S=sample(1:6,d,replace = TRUE)
        d=d-length(which(S==1)) #Gives difference between number of dice and number of times 1 appears. 
        num=num+1
      }
      return(num)
    }

    Nrolls=replicate(N,NoOfRolls(5))
    hist(Nrolls, breaks=0:60, prob=T)
    mean(Nrolls)#Average changes depending on no of dice thrown. This is the  average for 5 dice.
13.03

Upvotes: 0

Maurits Evers
Maurits Evers

Reputation: 50678

To compare with @TimBiegeleisen's answer, here is an alternative approach

We define a function that simulates rolling a 6-sided die, and that returns the minimum number of rolls necessary to get all sides at least once.

myfun <- function(Nmax = 10^5) {
    smpl <- sample(1:6, Nmax, replace = T)
    for (n in 1:Nmax) if (length(table(smpl[1:n])) == 6) break;
    return(n)
}

We now repeat the process 1000 times

set.seed(2018)
x <- replicate(1000, myfun())

summary(x)
#Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
#6.00   10.00   13.00   14.45   17.00   51.00

We plot the distribution

ggplot(data.frame(x = x), aes(x)) +
    geom_density()

enter image description here

Note that the value is in good agreement with the theoretical value

6 * sum(1 / (1:6))
[1] 14.7

So the result is close, with an absolute error percentage of

abs(mean(x) - 6 * sum(1 / (1:6))) / (6 * sum(1 / (1:6))) * 100
#[1] 1.727891

Upvotes: 2

Tim Biegeleisen
Tim Biegeleisen

Reputation: 521502

Here is a working script which counts how many times a die must be rolled to generate each of the six values:

numRolls <- function() {
    cnt <- 0
    x <- c(1:6)
    while (length(x) > 0) {
        rand <- sample(1:6,1,replace = TRUE)   # generate random value 1 to 6
        x <- x[which(x!=rand)]                 # remove this value if not yet seen
        cnt <- cnt + 1                         # increment number of rolls
    }

    return(cnt)
}

totalRolls <- 0

for (i in 1:1000) {
    totalRolls <- totalRolls + numRolls()
}

totalRolls / 1000
[1] 14.819

I ran 1000 tests, and got an average of 14.819 rolls needed to cover every value on a die.

Upvotes: 3

Related Questions