Probability that the histogram for the sum of two dice has an expected shape

Question

Suppose that two 6-sided dice are thrown $n$ times and that the sum after each throw is plotted on a histogram. Let $s_i$ be the frequency of the sum $i \in \{2, 3, \dots, 12 \}$. As $n \rightarrow \infty$, the probability that

$s_i \le s_{i+1}$ for $i \in \{2, 3, \dots, 6\}$ and $s_{i} \ge s_{i+1}$ for $i \in \{7, 8, \dots, 11\}$

approaches $1$ (because the further the sum is from $7$, the less likely it is to occur). But after only a finite number of throws, what is the probability that the above expression is true? Can this be expressed as a closed-form function of $n$?

Answer 1

As long as you can express a "dice pool" evaluation like this as a finite state machine with transitions determined by tuples (outcome, number of throws rolling that outcome) and keep the number of states in check, it's possible to compute the result in polynomial time of reasonable order. Here's a crude explanation of the algorithm; I'm working on a more refined explanation but it's not ready yet.

I implemented this approach as part of my hdroller Python library. You can try this computation in your browser using this JupyterLite notebook.

State transition function

By default, hdroller presents outcomes to the transition function in ascending order. Therefore, we can define the state transition like this:

import hdroller

class IsTriangleIsh(hdroller.EvalPool):
    def __init__(self, rise_end, fall_start):
        self.rise_end = rise_end
        self.fall_start = fall_start
    
    def next_state(self, state, outcome, count):
        # State consists of the previous count,
        # or -1 if the rise-fall pattern was violated.
        if state is None:
            return count
        if state == -1:
            return -1
        if outcome <= self.rise_end and count < state:
            return -1
        if outcome > self.fall_start and count > state:
            return -1
        return count
    
    def final_outcome(self, final_state, _):
        return final_state >= 0

evaluator = IsTriangleIsh(7, 7)
# 2d6
die = 2 @ hdroller.d6

Manually checking the result for 1 and 2 throws:

print(evaluator.eval(die.pool(1)))
print(evaluator.eval(die.pool(2)))

Denominator: 36

Outcome	Weight	Probability
False	30	83.333333%
True	6	16.666667%

With a single throw, the only accepting result is if the roll is exactly 7. There are 6 ways for this to happen.

Denominator: 1296

Outcome	Weight	Probability
False	1140	87.962963%
True	156	12.037037%

With two throws, there are three accepting possibilities:

• Both throws roll 7. (6 * 6 = 36 ways)
• The first throw rolls 7, the second throw rolls 6 or 8. (6 * 10 = 60 ways)
• The reverse of the above. (60 ways)

Total = 156 ways as computed.

Plotting

The exact chance for 500 throws takes about 4 minutes on my computer when running in JupyterLite. Results are memoized so we are not starting from scratch at each number of throws.

max_n = 500
x = [n for n in range(1, max_n+1)]
y = []

import time

start_ns = time.perf_counter_ns()
for n in x:
    y.append(evaluator.eval(die.pool(n)).mean() * 100.0)
end_ns = time.perf_counter_ns()
elapsed_s = round((end_ns - start_ns) * 1e-9)
print(f'Computation time: {elapsed_s} s')
    
import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(16, 9))
ax.plot(x, y)
ax.grid()
ax.set_xlabel('Number of throws')
ax.set_ylabel('Chance of fitting shape (%)')
ax.set_xlim(left=0, right=max_n)
ax.set_ylim(bottom=0)
plt.show()

Unless I've made an error somewhere, the convergence to 100% chance is pretty slow by the standards of throwing physical dice.

By the time you reach 500 throws you're probably better off using statistical methods rather than exact calculation, but I'll leave that to someone more versed in that approach.