Probability that the histogram for the sum of two dice has an expected shape 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$?
 A: 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.
