Roll a 6-sided die until the total $\geq M$. Mean amount by which $M$ is exceeded?

Question

Here is the question:

You roll a fair 6-sided dice iteratively until the sum of the dice rolls is greater than or equal to M. What is the mean and standard deviation of the sum minus M when M=300?

Should I write a code to answer these kind of questions?

Please give me some hints on that. thanks!

Answer 1

You can certainly use code, but I wouldn't simulate.

I'm going to ignore the "minus M" part (you can do that easily enough at the end).

You can compute the probabilities recursively very easily, but the actual answer (to a very high degree of accuracy) can be calculated from simple reasoning.

Let the rolls be $X_1, X_2, ...$. Let $S_t=\sum_{i=1}^t X_i$.

Let $\tau$ be the smallest index where $S_\tau\geq M$.

$P(S_\tau=M)=P(\text{got to }M-6\text{ at }\tau-1\text{ and rolled a }6) \\\qquad\qquad\qquad+ P(\text{got to }M-5\text{ at }\tau-1\text{ and rolled a }5)\\\qquad\qquad\qquad +\:\: \vdots\\\qquad\qquad\qquad+\:P(\text{got to }M-1\text{ at }\tau-1\text{ and rolled a }1)\\ \qquad\qquad\quad=\frac{1}{6}\sum_{j=1}^6 P(S_{\tau-1}=M-j)$

similarly

$P(S_\tau=M+1)=\frac{1}{6}\sum_{j=1}^5 P(S_{\tau-1}=M-j)$

$P(S_\tau=M+2)=\frac{1}{6}\sum_{j=1}^4 P(S_{\tau-1}=M-j)$

$P(S_\tau=M+3)=\frac{1}{6}\sum_{j=1}^3 P(S_{\tau-1}=M-j)$

$P(S_\tau=M+4)=\frac{1}{6}\sum_{j=1}^2 P(S_{\tau-1}=M-j)$

$P(S_\tau=M+5)=\frac{1}{6} P(S_{\tau-1}=M-1)$

Equations similar to the first one above could then (at least in principle) be run back until you hit any of the initial conditions to get an algebraic relationship between the initial conditions and the probabilities we want (which would be tedious and not especially enlightening), or you can construct the corresponding forward equations and run them forward from the initial conditions, which is easy to do numerically (and which is how I checked my answer). However, we can avoid all that.

The points' probabilities are running weighted averages of previous probabilities; these will (geometrically quickly) smooth out any variation in probability from the initial distribution (all probability at point zero in the case of our problem). The

To an approximation (a very accurate one) we can say that $M-6$ to $M-1$ should be almost equally probable at time $\tau-1$ (really close to it), and so from the above we can write down that the probabilities will be very close to being in simple ratios, and since they must be normalized, we can just write down probabilities.

Which is to say, we can see that if the probabilities of starting from $M-6$ to $M-1$ were exactly equal, there are 6 equally likely ways of getting to $M$, 5 of getting to $M+1$, and so on down to 1 way of getting to $M+5$.

That is, the probabilities are in the ratio 6:5:4:3:2:1, and sum to 1, so they're trivial to write down.

Computing it exactly (up to accumulated numerical round off errors) by running the probability recursions forward from zero (I did it in R) gives differences on the order of .Machine$double.eps ($\approx$2.22e-16 on my machine) from the above approximation (which is to say, simple reasoning along the above lines gives effectively exact answers, since they're as close to the answers computed from recursion as we'd expect the exact answers should be).

Here's my code for that (most of it's just initializing the variables, the work is all in one line). The code starts after the first roll (to save me putting in a cell 0, which is a small nuisance to deal with in R); at each step it takes the lowest cell which could be occupied and moves forward by a die roll (spreading the probability of that cell over the next 6 cells):

 p = array(data = 0, dim = 305)
 d6 = rep(1/6,6)
 i6 = 1:6
 p[i6] = d6
 for (i in 1:299) p[i+i6] = p[i+i6] + p[i]*d6

(we could use rollapply (from zoo) to do this more efficiently - or a number of other such functions - but it will be easier to translate if I keep it explicit)

Note that d6 is a discrete probability function over 1 to 6, so the code inside the loop in the last line is constructing running weighted averages of earlier values. It's this relationship that makes the probabilities smooth out (until the last few values we're interested in).

So here's the first 50-odd values (the first 25 values marked with circles). At each $t$, the value on the y-axis represents the probability that accumulated in the hindmost cell before we rolled it forward into the next 6 cells.

As you see it smooths out (to $1/\mu$, the reciprocal of the mean of the number of steps each die roll takes you) quite quickly and stays constant.

And once we hit $M$, those probabilities drop away (because we're not putting the probability for values at $M$ and beyond forward in turn)

So the idea that the values at $M-1$ to $M-6$ should be equally likely because the fluctuations from the initial conditions will get smoothed out is clearly seen to be the case.

Since the reasoning doesn't depend on anything but that $M$ is large enough that the initial conditions wash out so that $M-1$ to $M-6$ are nearly equally probable at time $\tau-1$, the distribution will be essentially the same for any large $M$, as Henry suggested in comments.

In retrospect, Henry's hint (which is also in your question) to work with the sum minus M would save a little effort, but the argument would follow very similar lines. You could proceed by letting $R_t=S_t-M$ and writing similar equations relating $R_0$ to the preceding values, and so on.

From the probability distribution, the mean and the variance of the probabilities are then simple.

Edit: I suppose I should give the asymptotic mean and standard deviation of the final position minus $M$:

The asymptotic mean excess is $\frac{5}{3}$ and the standard deviation is $\frac{2\sqrt 5}{3}$. At $M=300$ this is accurate to a much greater degree than you're likely to care about.

Answer 2

Let $\Omega$ be the set of sequences of partial sums of the dice rolls (with each sequence beginning at $0$). For any integer $n$, let $E_n$ be the event that $n$ appears in a sequence; that is,

$$E_n = \{\omega\in\Omega\,|\, n \in \omega\}.$$

Define $X_M(\omega)$ to be the first value in $\omega$ that equals or exceeds $M$. The question asks for properties of $X_M-M$. We can obtain the exact distribution of $X_M$, and from that everything follows.

First, notice that $X_M(\omega) - M \in \{0, 1, 2, 3, 4, 5\}$. By partitioning the event $X_M - M = k$ according to the immediately preceding value in $\omega$, and letting $p(i) = 1/6$ be the probability of observing face $i$ on one roll of the die ($i=1, 2, 3, 4, 5, 6$), it follows that

$$\Pr(X_M - M = k) = \sum_{j=k}^6 \Pr(E_{M+k-j}) p(j) = \frac{1}{6} \sum_{j=k}^6 \Pr(E_{M+k-j}).$$

At this point we could argue heuristically that, to a very good approximation for all but the smallest $M$, $$\Pr(E_i) \approx 2/7.$$ This is because the expected value of a roll is $(1+2+3+4+5+6)/6 = 7/2$ and its reciprocal should be the limiting, stable long-run frequency of any particular value in $\omega$.

A rigorous way to demonstrate this considers how $E_i$ could occur. Either $E_{i-1}$ occurs and the subsequent roll was a $1$; or $E_{i-2}$ occurs and the subsequent roll was a $2$; or ... or $E_{i-6}$ occurs and the subsequent roll was a $6$. This is an exhaustive partition of the possibilities, whence

$$\Pr(E_i) = \sum_{j=1}^6 \Pr(E_{i-j}) p(j) = \frac{1}{6} \sum_{j=1}^6 \Pr(E_{i-j}).$$

The initial values of this sequence are

$$\Pr(E_0) = 1;\quad \Pr(E_{-i}) = 0, i = 1, 2, 3, \ldots .$$

This plot of $\Pr(E_i)$ against $i$ shows how rapidly the chances settle down to a constant $2/7$, shown by the horizontal dotted line.

There is a standard theory of such recursive sequences. It can be developed by means of generating functions, Markov chains, or even algebraic manipulation. The general result is that a closed-form formula for $\Pr(E_i)$ exists. It will be a linear combination of a constant and the $i^\text{th}$ powers of roots of the polynomial

$$x^6 - p(1) x^5 - p(2) x^4 - p(3) x^3 \cdots - p(6) = x^6 - (x^5 + x^4+x^3+x^2+x+1)/6.$$

The largest magnitude of these roots is approximately $\exp(-0.314368)$. In a double precision floating point representation, $\exp(-36.05)$ is essentially zero. Therefore, for $i \gg -36.05 / -0.314368 = 115$, we may completely ignore all but the constant. This constant is $2/7$.

Consequently, for $M = 300 \gg 115$, for all practical purposes we may take $E_{M+k-j} = 2/7$, whence

$$\Pr(X_M - M = (0,1,2,3,4,5)) = \left(\frac{2}{7}\right)\left(\frac{1}{6}\right)(6,5,4,3,2,1).$$

Computing the mean and variance of this distribution is straightforward and easy.

---

Here is an R simulation to confirm these conclusions. It generates almost 100,000 sequences through $M+5 = 305$, tabulates the values of $X_{300} - 300$, and applies a $\chi^2$ test to assess whether the results are consistent with the foregoing. The p-value (in this case) of $0.1367$ is large enough to indicate they are consistent.

M <- 300
n.iter <- 1e5
set.seed(17)
n <- ceiling((2/7) * (M + 3*sqrt(M)))
dice <- matrix(ceiling(6*runif(n*n.iter)), n, n.iter)
omega <- apply(dice, 2, cumsum)
omega <- omega[, apply(omega, 2, max) >= M+5]
omega[omega < M] <- NA
x <- apply(omega, 2, min, na.rm=TRUE)
count <- tabulate(x)[0:5+M]
(cbind(count, expected=round((2/7) * (6:1)/6 * length(x), 1)))
chisq.test(count, p=(2/7) * (6:1)/6)