CODEWITHSUNDEEP
simulationprobabilitymodelingmontecarlo
raffy cee
raffy cee

Reputation: 51

Monte Carlo Simulations: Probability that 5 consecutive rolls of seven will occur

I just need your help about Monte Carlo Simulations. Here's the problem.

On a roll of two dice, a total of seven occurs with probability 1/6. In 100 rolls of the dice, what is the probability that five consecutive rolls of seven will occur?

(This is an exercise from Mathematical Modeling and Simulations by M. Meerschaert)

I made a pseudo-code using R:

dice.sims1 = function(trials){
occur
for (j in 1:trials){
x = 0
y = 0
n = 0
result = NULL
while (n<100){
n = n+1
result[n] = sample(1:6,1) + sample(1:6,1)
if(result[n] == 7){
x = x+1
}else {x = 0}
if(x>=5){
y = y+1
}
}
occur[j] = y 
}
print(mean(occur)/100)
}

Is this correct? If so, how do I interpret my result? If it's not correct, can you help me correct it?

Thanks and have a good day.

Upvotes: 1

Views: 333

Answers (1)

Severin Pappadeux
Severin Pappadeux

Reputation: 20130

Apart from bugs in the code, there is an unclear condition - if you have 6 rolls of 7 in a row, is it counted as one or as two?

You count it as two here 

if(x>=5){
    y = y+1
}

Then 7 rolls in a row would make it three etc.

UPDATE

dice.sims1 = function(trials) {
    occur = rep(0, times=trials)
    for (j in 1:trials) {
        x = 0
        y = 0
        n = 0
        while (n<100) {
            n = n + 1
            result = sample(1:6,1) + sample(1:6,1)
            if(result == 7){
                x = x+1
            } else {
                x = 0
            }
            if(x>=5) {
                y = y+1
                x = 0 # !!! unclear problem setup
            }
        }
        occur[j] = y
    }
    mean(occur)
}

set.seed(12345)

print(dice.sims1(10000))

will produce 0.0109 with marked line and 0.0131 with marked line removed

UPDATE II

Version without an array, direct summation

dice.sims <- function(trials) {
    s <- 0L
    for (i in 1:trials) {

        ncons <- 0L
        for(k in 1:100) {
            roll <- sample.int(6,1) + sample.int(6,1)
            if (roll == 7L) {
                ncons <- ncons + 1L # good roll
            } else {
                ncons <- 0L # bad one, back to square one
            }
            if (ncons >= 5L) {
                s <- s + 1L
                # ncons <- 0L # start from the beginning
            }
        }
    }
    as.numeric(s)/as.numeric(trials)
}

Upvotes: 1

Related Questions