# What is the limiting probability distribution of the output of a particular die-rolling process as the number of rolls goes to $\infty$?

Suppose I have a die with six arithmetic operations--- $${-}2, {-}1, \times 0, +1, +2, +3$$ ---and that each roll of the die is uniformly distributed.

To any finite sequence of rolls of the die, assign the value given by successively applying the operations to a starting value of $$0$$. So, for example, for the $$5$$-term sequence $$(+1,+3,{-}2,\times 0,-1)$$, the value would be: $$((((0 + 1) + 3) - 2) \times 0) - 1 = -1 .$$ For any number $$n$$ of rolls, we can ask for the fraction $$P_n(k)$$ of $$n$$-roll sequences that have value $$k$$.

How can I find an explicit expression for the function $$f(k) := \lim_{n \to \infty} P_n(k)?$$

• Not sure I follow the question. The sequence $0333$ is thrown with probability $\frac 1{6^4}$ so if you throw the die in blocks of four often enough you will throw this sequence with probability $1$. Thus in any sufficiently long string of tosses, this sequence will occur with probability $1$. Does that answer your question? Of course there are plenty of other ways in which you can land on $9$, but just looking at this particular one suffices.
– lulu
Jun 13, 2020 at 0:48
• I am also confused. Suppose the die is rolled precisely $n$ times, where $n$ is some finite number. Are you asking for the probability that the result will be $9$ after $n$ rolls, or are you asking that, as $n$ tends to infinity, what is the probability that $9$ appears at least once in our running count after some number of rolls, even if we continue to roll after that $9$ appeared? Jun 13, 2020 at 0:53
• Please look at my answer and Travis Willse's comment. We read the question quite differently. Please clarify. Jun 13, 2020 at 2:57
• @TravisWillse It's perfectly worded. Jun 13, 2020 at 20:39

You can model this as a Markov chain and there are known techniques for how to solve these problems. I'll explain how we can solve this example.

Let $$p_n$$ be the probability of having a total of $$n$$ after a large number of rolls. If we've gone on long enough then we should expect these probabilities to not change after our next roll. So, $$p_n = \frac{1}{6}(p_{n-3} + p_{n-2} + p_{n-1} + p_{n+1} + p_{n+2})$$ except when $$n=0$$, in which case $$p_0 = \frac{1}{6}(p_{-3} + p_{-2} + p_{-1} + p_{1} + p_{2} + 1)$$. That extra one represents the chance that we could come back to 0 from any number.

Note that if we didn't have the set 0 roll, then this technique wouldn't work because the solution would be that all the $$p_n$$s are equal, but this is impossible because there are infinitely many of them and they sum to 1. In that case we would have to limit ourselves to questions like what happens after $$t$$ throws instead of what happens asymptotically. I believe in this case we can solve this system, but I'm not exactly sure how off the top of my head.

Once we do solve this system, the solution is called the stationary distribution because it stays stationary after the next roll. There is a handy theorem that for any Markov chain that has a stationary distribution, it will approach the stationary distribution given enough time. I'm not sure the exact statement, but I believe it holds in this case. So all you have to do is solve that infinite system of equations.

• I think this will just give you the fact that on average you hit every number twice because the upward bias is $\frac 12$ per roll and does not answer the chance you hit a given number even if it is far beyond the startup transient. Jun 13, 2020 at 1:59
• We can deduce the asymptotic behavior of $p_n$ without solving the infinite system. Since for $n > 2$ we have $p_n + p_{n - 1} - 6 p_{n - 2} + p_{n - 3} + p_{n - 4} + p_{n - 5}$, $p_n$ (for $n > 2$) is some linear combination $\sum a_i r_i^n$, where the $r_i$ are the roots of the characteristic polynomial $x^5 + x^4 - 6 x^3 + x^2 + x + 1$. (Alas, its Galois group is the nonsolvable group $S_5$ and so the roots cannot be written in terms of radicals.) (cont.) Jun 13, 2020 at 3:40
• Our polynomial has two real roots $r_i$ with $|r_i| > 1$, and since each $p_n$ is bounded above by $1$, we know that the coefficients $a_i$ of those roots must be zero. The remaining roots are $\alpha = 0.82140\ldots$ and the conjugate nonreal roots $\beta = -0.27496\ldots + i 0.38561\ldots$ and $\bar \beta$. It follows that (for $n > 2$) we have $p_n = A \alpha^n + B (\beta^n + \bar\beta^n)$ for some constants $A, B$, and thus $p_n \sim A \alpha^n$ for large $n$. Jun 13, 2020 at 3:51
• @TravisWillse, awesome, I guess that's how we go about solving it. My only concern is what happens below $0$, because then the $\alpha^n$ term explodes and the $r_i^n$ terms do not. I guess then $A$ and $B$ are $0$ and the coefficients of $r_i^n$ are non-zero? Jun 13, 2020 at 10:42
• @A.Kriegman Yes, the rule must be different for (sufficiently negative) $n$, and like you say, the coefficients $A', B'$ of $\alpha, \beta$ must be zero. We don't know a priori that the coefficients (say, $C'$ and $D'$) of the other terms are both nonzero (at least one of them must be, since any negative value is possible). A quick Monte Carlo simulation suggests that the coefficient of the root of smaller absolute value, $\gamma = 1.77912 \ldots$, is nonzero, and so $p_{-n} \sim C' \gamma^{-n}$ for large $n$, for some constant $C'$. Jun 13, 2020 at 20:11

This answer expands on A. Kriegman's and folds in some of my comments thereunder.

Let $$P_n(k)$$ denote the fraction of values of $$n$$-term sequences with value $$k$$, which we can interpret as the probability that the value of a uniformly randomly selected $$n$$-term sequence has value $$k$$.

The limiting probabilities $$p_k := \lim_{n \to \infty} P_n(k)$$ are stable under the application of a uniformly selected die roll, giving an infinite set of equalities: $$\begin{array}{rcll} p_k &=& \frac{1}{6}(p_{k - 3} + p_{k - 2} + p_{k - 1} + p_{k + 1} + p_{k + 2}), & k \neq 0 \\ p_0 &=& \frac{1}{6}(p_{- 3} + p_{- 2} + p_{- 1} + p_{1} + p_{2} + 1) . \\ \end{array}\qquad (\ast)$$

The first equation defines a linear recurrence with characteristic polynomial $$p(r) = r^5 + r^4 - 6 r^3 + r^2 + r + 1,$$ and so the half-infinite sequences $$\{p_k\}_{k \leq 0}$$ and $$\{p_k\}_{k \geq 0}$$ can be given as linear combinations of powers $$\alpha_i^k$$ of the roots $$\alpha_i$$ of $$p$$ (possibly with different coefficients for $$k > 0$$ and $$k < 0$$).

The roots of $$p$$ are : $$\alpha = 0.82140\ldots, \quad \beta = -0.27496\ldots+i 0.38561 \ldots, \quad \bar\beta, \quad \gamma = 1.77912\ldots, \quad \delta = -3.05060\ldots .$$ Since $$0 \leq p_k \leq 1$$ for all $$k$$, the coefficients of $$\gamma, \delta$$ (whose real parts have absolute value $$> 1$$) must be zero for the sequence $$\{p_k\}_{k \geq 0}$$, and the coefficients of $$\alpha, \beta, \bar\beta$$ (whose real parts have absolute value $$< 1$$) must be zero for $$\{p_k\}_{k \leq 0}$$, and so $$\boxed{\begin{array}{rcll} p_k &=& A \alpha^k + B (\beta^k + \bar\beta^k), &k \geq 0 \\ p_k &=& C \gamma^k + D \delta^k , &k \leq 0 \end{array}\qquad (\ast\ast)}$$ for some constants $$A, B, C, D$$. (NB we can rewrite $$\beta^k + \bar\beta^k$$ as a manifestly real expression, namely, as $$2 e^{\operatorname{Re}(\beta) k} \cos (\operatorname{Im}(\beta) k)$$.) We can find those constants by producing an independent linear system in those variables and solving; one option is to substitute the expressions $$(\ast\ast)$$, $$k = -1,0,1$$ in $$(\ast)$$. We get one equation each from substituting the first and second equations in $$(\ast\ast)$$ in $$(\ast)$$, or we can replace one of those two equations with the condition $$A + 2 B = C + D$$ given by substituting $$k = 0$$ in both of the equations in $$(\ast\ast)$$.

Appealing to a C.A.S. produces explicit formulae for $$A, B, C, D$$ as rational polynomials in $$\alpha, \beta, \gamma, \delta$$, but the expressions are unwieldy (hundreds of thousands of characters among them), and it is not evident that they can be simplified further. Their numerical values are: \boxed{\begin{align*} A &= 0.13210\ldots\\ B &= 0.04359\ldots\\ C &= 0.15602\ldots\\ D &= 0.06328\ldots . \end{align*}} In particular, $$p_0 = 0.21930\ldots$$.

Since $$A, C \neq 0$$, the limiting behaviors of $$p_k$$ are \begin{align*} p_k \sim A \alpha^k ,&\quad k \to \phantom{-}\infty \\ p_k \sim C \gamma^k ,&\quad k \to -\infty . \end{align*}

Remark One might ask whether we can produce exact expressions for the roots $$\alpha, \beta, \ldots$$ of the (quintic) polynomial $$p$$. If we restrict ourselves to algebraic expressions, we cannot: By reducing modulo $$2$$ we can efficiently deduce that $$p$$ is irreducible over $$\Bbb Q$$, so its Galois group contains a $$5$$-cycle. On the other hand, we've seen that $$p$$ has exactly $$2$$ nonreal roots, and hence the complex conjugation map is a transposition in the Galois group of $$p$$. But a transposition and a $$5$$-cycle generate all of $$S_5$$, which is hence the Galois group. In particular, it is not solvable, so the roots $$\alpha, \beta, \ldots$$ are not expressible in terms of radicals.

The limit that you hit any number $$k$$, positive or negative, goes to $$1$$ as $$n \to \infty$$. Say we want the chance of hitting $$k=27$$. This is higher than the chance we get a $$0$$ and then $$9\ +3$$'s in a row because there are other ways to get to $$27$$, but the chance you get that string in $$n$$ throws is $$1-6^{9-n}$$. This goes to $$1$$ as $$n \to \infty$$. The reset to $$0$$ allows us to overcome the upward bias of $$\frac 12$$ per throw.

• Perhaps I misunderstand OP's question---which, to be honest, I found not entirely clear---but as I understand it, each time we roll a zero on the die, we reset our total $T$ to zero. So, for example, if we denote our total after $n$ rolls by $T_n$, then $\Bbb P(T_n = 0) > \frac{1}{6}$ for all $n$. Jun 13, 2020 at 2:44
• @TravisWillse: I think that is an excellent point. Yes, it makes my approach badly wrong. In that case it is easy. For any $n$, as the number of rolls goes to infinity the chance you hit $n$ goes to $1$. For example, the chance we eventually hit $27$ is greater than the chance we have had a run of $0$ followed by $9\ +3$s because there are other ways to get to $27$. For $n\gt 9$ throws the chance we have seen that string is $1-6^{9-n}$ which goes to $1$. Jun 13, 2020 at 2:56
• Again, I might well be interpreting OP incorrectly, but I read OP as asking not for the probability for a given $k \in \Bbb Z$ that $T_n = k$ for some $n$, but rather for a formula for $\lim_{n \to \infty} \Bbb P(T_n = k)$ as a function of $k$. It would anyway help if OP clarified their intent. Jun 13, 2020 at 3:31
• How do you get the fancy T_n and stuff that might help in the future? I'm new at this. Jun 13, 2020 at 14:16
• @BriggyT: It is called MathJax and is similar to $\LaTeX$. You put it between dollar signs. A tutorial is here. You can then right click on any piece and choose Show Math As -> TeX Commands to see how it was done. Jun 13, 2020 at 14:23