Rolling a dice until the average of the outcome hits a specific value.

Question

I roll the dice until the average of the outcomes is $2.5$. Once the average is $2.5$, I will stop rolling. For example, if I get a sequence of 2, 6, 4, 1, 1, 1 outcomes, I will stop rolling the dice at the sixth trial.

Let $E$ be the event of stopping rolling the dice and $E_n$ be the event of getting the average $2.5$ after rolling the dice exactly $n$ times. Then, $E = E_1 \cup E_2 \cup E_3 \cup E_4 \cup \cdots$. Since it will not end with odd trials, I can rewrite it as $E = E_2 \cup E_4 \cup \cdots$.

Since ${E_{2n}}^\complement \supseteq \bigcup_{k=2n}^{\infty} E_{2k}$, $E_{2k}$'s are disjoint. Thus, $Pr(E) = \sum_{k=1}^{\infty} Pr(E_{2k})$. However, I am stuck with calculating the exact probability of each event $E_{2k}$ for large $k$. How can I calculate the probability of each event $E_{2k}$ ($Pr(E_{2k})$ for $k>1$) and finally the probability of stopping rolling the dice $Pr(E)$?

It seems that $Pr(E_{2k})$ goes to zero as $k\rightarrow \infty$. Then, would the $Pr(E)<1$? Or, is the probability of stopping rolling the dice eventually one?

Also, I wonder if there is any related theorem that I can use to calculate such events.

Answer 1

In this answer, we

1. prove $\mathbf{P}(E) < 1$ in OP' scenario, and
2. provide a theoretical way of approximating the value of $\mathbf{P}(E)$ to any prescribed accuracy. For example, this will allow us to obtain the first $20$ digits in the decimal expansion of $\mathbf{P}(E)$:

$$\mathbf{P}(E) \approx 0.21257819505882691800$$

---

Setting.

Let $X_n$ denote the outcome of the $n$th roll, and let $S_n = \sum_{k=1}^{n} X_k$ be the partial sum. Also, fix $\alpha$ and define $E_n$ be the event that $n$ is the first time $S_n = \alpha n$ holds, excluding $n=0$. That is,

$$ E_n = \{ \text{$S_n = \alpha n$ and $S_k \neq \alpha k$ for all $k \in [n-1]$} \},$$

where $[n-1] = \{1, 2, \ldots, n-1\}$. Then we are interested in the value of $\mathbf{P}(E)$, where

$$ E = \bigcup_{n=1}^{\infty} E_n = \{\text{$S_n = \alpha n$ for some $n\geq 1$} \}. $$

---

Finding a representation of $\mathbf{P}(E)$.

To simplify this, we write $A_n = \{S_{n} = \alpha n\}$. Then by the inclusion-exclusion principle,

\begin{align*} \mathbf{P}(E_n) = \mathbf{P}(A_n \cap A_{n-1}^c \cap \cdots \cap A_1^c) = \sum_{J \subseteq [n-1]} (-1)^{|J|} \mathbf{P}\biggl(A_n \cap \bigcap_{j \in J} A_j \biggr). \end{align*}

To further simplify the last sum, we write $l = |J|+1$ and $J = \{j_1 < \ldots < j_{l-1}\}$, together with $j_l = n$. Then

\begin{align*} (-1)^{|J|} \mathbf{P}\biggl(A_n \cap \bigcap_{j \in J} A_j \biggr) &= (-1)^{l-1} \mathbf{P}(A_{j_1} \cap A_{j_2} \cap \cdots \cap A_{j_l} ) \\ &= (-1)^{l-1} \mathbf{P}(A_{j_1})\mathbf{P}(A_{j_2-j_1}) \cdots \mathbf{P}(A_{j_l-j_{l-1}}) \\ &= (-1)^{l-1} \mathbf{P}(A_{d_1})\mathbf{P}(A_{d_2}) \cdots \mathbf{P}(A_{d_l}), \end{align*}

where $d_1 = j_1$, $d_2 = j_2 - j_1$, $\ldots$, and $d_l = j_l - j_{l-1}$. (In the second step, we utilized the fact that $S_n$ is a random walk, having stationary and independent increments.) As $J$ runs over the $(l-1)$-subsets of $[n-1]$, the tuple $(d_1, \ldots, d_l)$ runs over all the positive values satisfying $d_1, \ldots, d_l \geq 1$ and $d_1 + \cdots + d_l = n$. Hence,

\begin{align*} \mathbf{P}(E_n) &= \sum_{l=1}^{n} (-1)^{l-1} \sum_{\substack{d_1, \ldots, d_l \geq 1 \\ d_1 + \cdots + d_l = n}} \mathbf{P}(A_{d_1})\mathbf{P}(A_{d_2}) \cdots \mathbf{P}(A_{d_l}). \end{align*}

Now we want to plug this to $\mathbf{P}(E) = \sum_{n=1}^{\infty} \mathbf{P}(E_n)$ and rearrange this sum, but this is not always possible. To overcome this issue, we consider regularized sums

\begin{align*} f(z) &= \sum_{n=1}^{\infty} \mathbf{P}(E_n) z^n, \\ g(z) &= \sum_{n=1}^{\infty} \mathbf{P}(A_n) z^n. \end{align*}

Note that each of $f(z)$ and $g(z)$ defines an analytic function on the region $|z| < 1$. Moreover, $\mathbf{P}(E) = \lim_{z\uparrow 1} f(z)$. Now let $\delta > 0$ be small enough so that $g(\delta) < 1$. Then for $z \in [0, \delta]$, Fubini's theorem is applicable, yielding

\begin{align*} f(z) &= \sum_{n=1}^{\infty} \biggl[ \sum_{l=1}^{n} (-1)^{l-1} \sum_{\substack{d_1, \ldots, d_l \geq 1 \\ d_1 + \cdots + d_l = n}} \mathbf{P}(A_{d_1})\mathbf{P}(A_{d_2}) \cdots \mathbf{P}(A_{d_l}) \biggr] z^n \\ &= \sum_{l=1}^{\infty} (-1)^{l-1} \sum_{d_1, \ldots, d_l \geq 1} \mathbf{P}(A_{d_1})\mathbf{P}(A_{d_2}) \cdots \mathbf{P}(A_{d_l}) z^{d_1 + \cdots + d_l} \\ &= \sum_{l=1}^{\infty} (-1)^{l-1} g(z)^l \\ &= \frac{g(z)}{1 + g(z)}. \end{align*}

Although this identity has been proved a priori for $z \in [0, \delta]$, both sides are analytic function on $|z| < 1$ and hence the validity of the equality remains true on all of $|z| < 1$ by the principle of analytic continuation. In particular,

$$ \bbox[color:navy; border:1px dotted navy; padding:8px;]{ \mathbf{P}(E) = \lim_{z\uparrow 1} f(z) = \lim_{z\uparrow 1} \frac{g(z)}{1 + g(z)} = \frac{g(1)}{1+g(1)}, } $$

where $g(1) = \sum_{n=1}^{\infty} \mathbf{P}(A_n)$ and we adopt the convention that $f(1) = 1$ if $g(1) = \infty$. In particular, $\mathbf{P}(E) < 1$ if and only if $g(1) < \infty$.

---

Estimating $\mathbf{P}(E)$ for $\alpha = 5/2$. Now we turn to OP's problem. Let $\alpha = 5/2$. Then $\mathbf{P}(A_n) = 0$ for odd $n$'s, hence

$$ g(1) = \sum_{d=1}^{\infty} \mathbf{P}(A_{2d}) = \sum_{d=1}^{\infty} \mathbf{P}(S_{2d} = 5d). $$

By a simple combinatorial argument,

\begin{align*} \mathbf{P}(S_{2d} = 5d) &= \frac{1}{6^{2d}} \#\{ (x_1, \ldots, x_{2d}) \in [6]^{2d} : x_1 + \cdots + x_{2d} = 5d \} \\ &\leq \frac{1}{6^{2d}} \#\{ (x_1, \ldots, x_{2d}) \in \mathbb{N}_1^{2d} : x_1 + \cdots + x_{2d} = 5d \} \\ &= \frac{1}{6^{2d}} \binom{5d-1}{2d-1}. \tag{*} \end{align*}

Also, it is not hard to prove that $\binom{5(d+1)-1}{2(d+1)-1}\big/\binom{5d-1}{2d-1} \leq \frac{5^5}{2^2 \cdot 3^3}$ for all $d \geq 1$. Hence

$$ \mathbf{P}(S_{2d} = 5d) \leq \frac{1}{6^{2d}} \binom{5-1}{2-1} \left(\frac{5^5}{2^2 \cdot 3^3}\right)^{d-1} = \frac{432}{3125} \left( \frac{3125}{3888} \right)^d. $$

This shows that $g(1) < \infty$, hence $\mathbf{P}(E) < 1$ as required.

This estimate can also be used to estimate $\mathbf{P}(E)$. To this end, we first derivate an exact formula for $\mathbf{P}(S_{2n} = 5n)$. Invoking an inclusion-exclusion type computation, We can give an explicit formula for $\mathbf{P}(S_{2n} = 5n)$:

\begin{align*} \mathbf{P}(S_{2n} = 5n) &= \frac{1}{6^{2n}} \sum_{x_1, \ldots, x_{2n} \geq 1} \mathbf{1}[x_1 + \cdots + x_{2n} = 5n] \prod_{j=1}^{2n} \mathbf{1}[x_j \leq 6] \\ &= \frac{1}{6^{2n}} \sum_{x_1, \ldots, x_{2n} \geq 1} \mathbf{1}[x_1 + \cdots + x_{2n} = 5n] \prod_{j=1}^{2n} (1 - \mathbf{1}[x_j > 6]) \\ &= \frac{1}{6^{2n}} \sum_{J \subseteq [2n]} (-1)^{|J|} \sum_{x_1, \ldots, x_{2n} \geq 1} \mathbf{1}[x_1 + \cdots + x_{2n} = 5n] \prod_{j \in J} \mathbf{1}[x_j > 6] \\ &= \frac{1}{6^{2n}} \sum_{J \subseteq [2n]} (-1)^{|J|} \sum_{x_1, \ldots, x_{2n} \geq 1} \mathbf{1}[x_1 + \cdots + x_{2n} = 5n - 6|J|] \\ &= \frac{1}{6^{2n}} \sum_{J \subseteq [2n]} (-1)^{|J|} \binom{5n-6|J|-1}{2n-1} \cdot \mathbf{1}[5n-6|J|\geq 2n] \\ &= \frac{1}{6^{2n}} \sum_{l=0}^{\lfloor n/2 \rfloor} (-1)^{l} \binom{5n-6l-1}{2n-1} \binom{2n}{l} \end{align*}

Moreover, using the bound on $\mathbf{P}(S_{2n} = 2n)$ above, we get

$$ \epsilon_N := \sum_{n=N+1}^{\infty} \mathbf{P}(S_{2n} = 5n) \leq \frac{432}{763} \left(\frac{3125}{3888}\right)^N. $$

Altogether, when $\alpha = 5/2$ we can decompose $g(1)$ into

$$ \bbox[color:navy; border:1px dotted navy; padding:8px;]{ g(1) = \sum_{n=1}^{N} \frac{1}{6^{2n}} \sum_{l=0}^{\lfloor n/2 \rfloor} (-1)^{l} \binom{5n-6l-1}{2n-1} \binom{2n}{l} + \epsilon_N, } $$

which then allows us to produce rigorous bounds on the value of $\mathbf{P}(E) = \frac{g(1)}{1+g(1)}$ up to any desired accuracy. For example, this allows us to prove that

$$\mathbf{P}(E) \approx 0.21257819505882691800$$

by choosing $N = 208$ and using $\epsilon_{208} < 10^{-20}$.

Answer 2

Let the probability we achieve an average score of $2.5$ after rolling the die $n$ times be $p_n$ (NB: this ignores the possibility that the target average has been achieved earlier, but we'll come to that later). If $T_n$ is the cumulative total score after roll $n$, then $p_n$ is the probability that $2T_n=5n$.

Clearly we need $n$ to be even for this to be possible, so put $n=2m$; then $T_n=5m$ and we're looking to count the number of ways we can make this total with $2m$ ordered die rolls.

A trick we can use is that this number is given by the coefficient of $x^{5m}$ in the expansion of $$F_m(x)=\left(x+x^2+x^3+x^4+x^5+x^6\right)^{2m}$$

Now, $$\begin{align}F_m(x)&=x^{2m}\left(1+x+x^2+x^3+x^4+x^5\right)^{2m} \\ &=x^{2m}\left[\frac{1-x^6}{1-x}\right]^{2m} \\ &=x^{2m}\left(1-x^6\right)^{2m}(1-x)^{-2m} \end{align}$$

The two bracketed terms are binomials, which we can expand; we have $$\left(1-x^6\right)^{2m} = \sum_{k=0}^n (-1)^k \binom{2m}{k} x^{6k}$$ and $$(1-x)^{-2m} = \sum_{l=0}^\infty \binom{2m-1+l}{2m-1} x^l$$

Combining all this, the coefficient of $x^{5m}$ in $F_m(x)$ is $$C_m=\sum_{k=0}^\left\lfloor \frac{m}{2} \right\rfloor (-1)^k \binom{2m}{k} \binom{5m-6k-1}{2m-1}$$

---

To recap, we've now got a formula for $C_m$, the number of ways we can achieve a score of $5m$ by rolling a die $2m$ times. The probability we have an average score of $2.5$ after $2m$ rolls is $p_{2m}=\frac{C_m}{6^{2m}}$.

This is good, but not quite what we need. We've not factored in the fact that once the target average score is achieved, we stop rolling.

(I had a mistake here earlier - thanks to Daniel Mathias for pointing it out. The question has now been answered better elsewhere; if I get anything useful to add I'll edit this again)

Answer 3

The answers already given contain great approaches to the problem. I want to provide another one, but only schematically show how it would work on a simpler example. The method should also be adaptable to the example given in OPs question.

My easier example is a two-sided die, i.e. it can only roll a $1$ or a $2$ (with equal probabilities). Actually, to make it easier I shift and rescale the variable so it takes the values $-1$ or $1$ (so the mean is $0$). This can be reinterpreted as a random walk on $\mathbb{Z}$. At each step of the walk the walker randomly chooses one of the two neighbouring points. This random choice corresponds to rolling either a $1$ or a $-1$ on the dice.

The first thing I want to calculate is the probability that at some point the average is equal to $0$. I will later generalize that to the probability that at some point the average is $\alpha$.

In the language of random walks the probability that eventually our steps sum up to $0$ again after we started rolling (or walking) is called the return probability. It is well known that this is 1 in our case (due to a result of the name "Pólya’s recurrence theorem"), but I want to rederive that now. To do that call $p(x)$ the probability that we eventually reach zero when we start at point $x$. Trivially $p(0)=1$ because we already are at $0$. Also it is clear that $p$ fulfills the difference equation $p(x-1)=\frac{1}{2}(p(x)+p(x-2))$. This is because with a probability of $\frac{1}{2}$ we go to $x+2$ in the next step/roll and with the same probability we go to $x$. Defining the difference operator $\Delta f(x)=f(x)-f(x-1)$ this equation can be rewritten as $\Delta^2 u(x)=0$. This should remind you of the differential equation $u''(x)=0$ and it has the same general solution (as you can easily show): A linear polynomial, i.e. $p(x)=Ax+B$.

We would like to wish for one final thing: $p$ denotes a probability, thus should never exceed $1$. But the general solution is unbounded unless $A=0$. So we must take $A=0$. Also we already know $B=p(0)=1$. Combining everything we have $p(x)=1$. So we know that from every possible starting point on our random walk we will eventually get to $0$ with probability $1$.

The thing we want to know (i.e. the probability that we eventually return to $0$) can now be easily computed. At our first roll we either go to $1$ or $-1$. So our return probability will be $\frac{1}{2}(p(1)+p(-1))$ which is $1$.

A little warning: Something I didn't really pay attention to is that we can't necessarily take one solution of the difference equation derived above for $p(x)$ for all $x$. Instead we need on solution for positiv $x$ and one for negative $x$. Why? Well, we to get from negative to positive we have to cross $0$, but $0$ is a special point where we would "stop rolling" because we reached our goal. But the argument with the unboundedness will of course work in the positive and negative regions seperately and nothing changes. But we will have to keep that in mind for the generalization.

Now we want to calculate the probability that tha average of our rolls eventually is $\alpha$ (for some $\alpha$). We can encompass that in our discussion from above by just shifting $x$ by $-\alpha$ each roll. Why? Well, $x$ is the sum of all rolls. So if the average is supposed to be equal to $\alpha$, then the sum needs to be equal to $n\alpha$ where $n$ is the number of rolls. The latter will change by $\alpha$ each roll. We can instead shift each roll $x$ by $-\alpha$ to account for that. This basically means that at the $n$-th roll $x$ would be shifted to $x'=x-n\alpha$ in total. So if $x=n\alpha$ before the shift, after the shift we have $x'=0$. This means we reformulated our problem to look for another return probability, just looking at $x'$ instead of $x$. So the new difference equation would be $p(x'-1)=\frac{1}{2}(p(x'-\alpha)+p(x'-2-\alpha))$.

The procedure would now be similar to before, but significantly more complicated. I will show the case $\alpha=0.5$ but any rational $\alpha$ should work. Now the variables are half integers, so I express everything in terms of $y=2x'$ instead (which is an integer). Plugging that into our difference equation and taking $q(y)=p(x')=q(\frac{y}{2})$, we have $q(y-2)=\frac{1}{2}(q(y-1)+q(y-5))$ or (to make it look nicer) $q(y-1)=\frac{1}{2}(q(y)+q(y-4))$. This can be written as $(\Delta^3 +4\Delta^2+2\Delta-2)\Delta q(x)=0$. Again, we solve this like a differential equation: First compute the roots of the characteristic polynomial $(\lambda^3 +4\lambda^2+2\lambda-2)\lambda$. One root is obviously $0$ and corresponds to a constant solution. To get the others we coul use Cardano's formula for cubic equations, but generally we will have to solve it numerically (because for different $\alpha$ or more faces on our die we will get polynomials of higher degree for which no general formulas exist). The three nonzero roots (let's call them $\lambda_{1,2,3}$ correspond to exponential solutions. But be careful: We won't have $\lambda^y$ as a solution! To see that, write the eigenvalue equation $\Delta f(x)=\lambda f(x)$ and insert $f(x)=a^x$. This gives $a^x-a^{x-1}=\lambda a^x$ which is only true for $1-\frac{1}{a}=\lambda$. Great, we can now obtain $a_{1,2,3}$ from $\lambda_{1,2,3}$ and now our general solution has the form $q(y)=A+Ba_1^y+Ca_2^y+Da_3^y$.

Now we have to be careful again. Naively we might just use the same argument as before: This function is unbounded or constant, but our probabilities need to be bounded, so it has to be constant. But this would again lead to $q(y)=1$ for all $y$ which we know is not true (e.g. from simulation). What went wrong? Well, as mentioned before, we have to take seperate solutions for positive and negative $y$.The roots give us $a_1>1$ and $0<a_{2,3}<1$ (I have approximately $a_1\approx 1.92$, $a_2\approx 0.43$ and $a_3\approx 0.24$). So the general bounded solution would be of the form $q(y)=A+Ba_1^y$ for $y\le 0$ and $q(y)=C+Da_2^y+Ea_3^y$ for $y\ge 0$. Our "boundary conditions" connecting the two regions are $q(0)=1$, $q(1)=\frac{1}{2}(q(-2)+q(2))$ and $q(2)=\frac{1}{2}(q(-1)+q(3))$. These are not enough equations, so we are missing something. I think we can assume $A=\lim_{y\to-\infty}q(y)$ to be $0$. This also implies $B=1$

With this we have $3$ unknowns $C,D,E$ and $3$ linear equations. Solving that system I got $D\approx 0.49$, $E\approx 0.19$ and $C\approx 0.32$.

Now to compute the return probability: The first step will take us with probability $\frac{1}{2}$ to $y=1$ and with probability $\frac{1}{2}$ to $y=-3$. So the return probability is $\frac{1}{2}(q(1)+q(-3))=\frac{1}{2}(C+Da_2+Ea_3+a_1^{-3})\approx 0.36$. I tried writing a program to simulate that, and I got values close to that but typically below it (which makes sense because I can only have it roll finitely many times and convergence might be slow).

I hope I didn't make any significant mistakes here and this gives OP a different perspective on the problem as the solution of a difference equation.