Probability to roll a sum of less than 3600 with 1000 dice throws

Question

I try to estimate the probability that the sum of the numbers when rolling a die 1000 times is less than 3600. I am confident that I can approximate the probability with a normal distribution with expectation value $\mu=3500$ and the correct standard deviation. Is $\sigma=\sqrt{\frac{35n}{12}} (n=1000)$ the correct value? If I use this $\sigma$ then $p=96.8\%$ which looks pretty high to me. How can I verify my result?

You can repeatedly simulate $1000$ dice rolls, it shouldn't be too bad. — Qiaochu Yuan, Jul 25 '22 at 09:18
96.8% looks about right to me, you are expecting to get 3,500 over 1,000 trials, so 3,600 is pretty unlikely. — Suzu Hirose, Jul 25 '22 at 09:20

score 2 · Accepted Answer · answered Jul 25 '22 at 09:35

Pretty easy to simulate in python.

import numpy as np
num_rolls = 1000
num_trials = 100000
cutoff = 3600
samples = np.random.choice(range(1,7),size=(num_trials, num_rolls)).sum(axis=1)
print((samples<cutoff).mean())

My sim agrees with your calculation.

score 2 · Answer 2 · answered Jul 25 '22 at 09:59

You should have a continuity correction in your normal approximation. For example in R, depending whether "less than $3600$" is intended to be strict or not:

maxrolls <- 1000
faces <- 6
pnorm(3600-0.5, (faces+1)*maxrolls/2, sqrt(maxrolls*(faces^2-1)/12))
# 0.9672904
pnorm(3600+0.5, (faces+1)*maxrolls/2, sqrt(maxrolls*(faces^2-1)/12)) 
# 0.9686207

and it is possible to do a more exact calculation of the probabilities

maxrolls <- 1000
faces <- 6
probs <- numeric(faces*maxrolls+1) # offset by 1
probs[1] <- 1
for (roll in 1:maxrolls){
  matprobs <- matrix(
    rep(c(rep(0,faces+1), probs), faces)[-(1:(faces))], 
    ncol=faces)
  probs <- (rowSums(matprobs) / faces)[1:(faces*maxrolls+1)]
}
sum(probs) # should be 1
# 1
sum(probs * (0:(faces*maxrolls))) - (faces+1)*maxrolls/2 # should be 0
# 1.818989e-12
sum(probs[1:3600]) / sum(probs) # prob strictly less than 3600
# 0.9672951
sum(probs[1:3601]) / sum(probs) # prob 3600 or less
# 0.9686258

and this shows that the normal approximation with a continuity correction is very close in this case, with both calculations within $0.000005$ of the correct values

score 1 · Answer 3 · answered Jul 25 '22 at 09:46

The mean value $\mu$ for throwing dice $\mu=\frac16(1+2+...+6)=3.5$, and the standard deviation is $$ \begin{align} \sigma&=\sqrt{\frac16\{(6-\mu)^2+...+(1-\mu)^2\}}\\ &=\sqrt{35\over 12} \end{align} $$ so your numbers for the normal approximation for the distribution after $n$ trials seem to be OK.

Numerically,

package main
import (
    "fmt"
    "math/rand"
    "time"
)
func main() {
    nTrials := 100000
    nRolls := 1000
    over := 0
    rand.Seed(time.Now().UnixNano())
    for t := 0; t < nTrials; t++ {
        total := 0
        for r := 0; r < nRolls; r++ {
            total += rand.Intn(6) + 1
        }
        if total >= 3600 {
            over++
        }
    }
    fmt.Printf("%d %d %g\n", over, nTrials,
        float64(over)/float64(nTrials))
}

gives a value of

3246 100000 0.03246

which is the same as your answer

Run online here.

so it looks good.

Probability to roll a sum of less than 3600 with 1000 dice throws

3 Answers3