If I were that good at weighting dice, I'm not sure why I'd choose the uniform distribution

When you roll a pair of fair dice, the most likely outcome is $7$ (which occurs $1/6$ of the time) and the least likely outcomes are $2$ and $12$ (which each occur $1/36$ of the time). Annoyed by the variance of these probabilities, I set out to create a pair of “uniform dice.” These dice still have sides that are uniquely numbered from $1$ to $6,$ and they are identical to each other. However, they are weighted so that their sum is more uniformly distributed between $2$ and $12$ than that of fair dice. Unfortunately, it is impossible to create a pair of such dice so that the probabilities of all 11 sums from $2$ to $12$ are identical (i.e., they are all $1/11$), but I bet we can get pretty close.

So how should I make my dice as uniform as possible? In other words, which specific weighting of the dice minimizes the variance among the $11$ probabilities?

The canonical setup assumes that I am making two identical dice with weighted distributions such that the face $i$ will land with probability $p_i,$ for $i = 1, \dots, 6.$ However, I will set up the generic problem where I can have weights $p = (p_1, \dots, p_6)$ for one die and $q = (q_1, \dots, q_6)$ for the other.

Given $p$ and $q$, the probability of rolling $s$ for some $s \in \{ 2, \dots, 12 \}$ is $$\mathbb{P} \{ S = s \mid p, q \} = \sum_{i= \max\{s-6,1\}}^{\min\{s-1, 6\}} p_i q_{s-i}.$$ So the variance of the probabilities is given by $$V(p, q) = \frac{1}{11} \sum_{s=2}^12 \left( \mathbb{P} \{S = s \mid p, q \} - \frac{1}{11} \right)^2.$$ Obviously the choice of weights $p$ and $q$ must satisfy $p^T \mathbb{1} = q^T \mathbb{1} = 1$ and $p, q \geq 0$ in order for them to be true probabhility weights.

Therefore, the optimal weightings to obtain as uniform as possible a sum of the two dice rolls is therefore given by the optimization problem \begin{align*} \text{minimize} \,\,\, & \,\, V(p, q) \\ \text{such that} & \,\, p^T \mathbb{1} = q^T \mathbb{1} = 1 \\ & \,\, p, q \in \mathbb{R}^6_+.\end{align*} In the canonical version of this problem, the additional linear constraint that $p = q$ is added.

I set up a small Python script (see below) to do the optimization using the SLSQP solver method for scipy.optimize.minimize. For the canonical problem, the solution is given by $$\hat{p} = \hat{q} = (0.2439, 0.1375, 0.1186, 0.1186, 0.1375, 0.2439),$$ which gives minimal probability variance of $V(\hat{p}, \hat{p}) = 0.00121758.$ In the general problem, however, we can do quite a bit better, reducing the probability variance by about a factor of 6. In the general case, the solution is \begin{align*}\hat{p} &= (0.1250, 0.1875, 0.1875, 0.1875, 0.1875, 0.1250) \,\,\text{and}\\ \hat{q} &= (0.5000, 0.0000, 0.0000, 0.0000, 0.0000, 0.5000),\end{align*} (or vice versa), with minimal probability variance of $V(\hat{p}, \hat{q}) = 0.00025826.$

In [1]:

import numpy as np
from scipy.optimize import Bounds, minimize, LinearConstraint

def F(x):
    p = x[:6].copy()
    q = x[6:].copy()
    prob = np.zeros((11,))
    jac = np.zeros((12,))
    for s in range(1,12):
        innerj = np.zeros((12,))
        lower = max(0,s-6)
        upper = min(6,s)
        for i in range(lower, upper):
            prob[s-1] += p[i] * q[s-i-1]
            innerj[i] += q[s-i-1]
            innerj[6+s-i-1] += p[i]
        jac += 2.0 / 11.0 * (prob[s-1] - 1.0/11.0) * innerj
    var = (prob - 1.0/11.0).dot(prob - 1.0/11.0) / 11.0
    return var, jac, prob

def objective(x):
    o, _, _ = F(x)
    return o

def jac(x):
    _, j, _ = F(x)
    return j

def probs(x):
    _, _, p = F(x)
    return p

A = np.concatenate([np.concatenate([np.ones((1,6)), np.zeros((1,6))], axis=1), 
                    np.concatenate([np.zeros((1,6)), np.ones((1,6))], axis=1)], axis=0)

bounds = Bounds(np.zeros((12,)), np.ones((12,)))
r = np.random.uniform(size=(12,))
x0 = np.concatenate([r[:6]/r[:6].sum(), r[6:]/r[6:].sum()])
soln = minimize(objective, x0, jac=jac, method='SLSQP', bounds=bounds, options={'ftol': 1e-16, 'disp': True},
        constraints=[{'type' : 'eq', 'fun' : lambda x: A.dot(x) - np.ones((2,)), 'jac' : lambda x : A}], 
)
print('The solution for the general problem is given by p=(%0.4f, %0.4f, %0.4f, %0.4f, %0.4f, %0.4f) and q=(%0.4f, %0.4f, %0.4f, %0.4f, %0.4f, %0.4f) (or vice versa!).\nThese optimizers yield optimal value of V(p,q)=%0.8f.' % (*soln.x, soln.fun))

Optimization terminated successfully.    (Exit mode 0)
            Current function value: 0.0002582644628099179
            Iterations: 47
            Function evaluations: 47
            Gradient evaluations: 47
The solution for the general problem is given by p=(0.1250, 0.1875, 0.1875, 0.1875, 0.1875, 0.1250) and q=(0.5000, 0.0000, 0.0000, 0.0000, 0.0000, 0.5000) (or vice versa!).
These optimizers yield optimal value of V(p,q)=0.00025826.

In [2]:

def Fc(x):
    xx = np.concatenate([x,x])
    return F(xx)

def objectivec(x):
    oc, _, _ = Fc(x)
    return oc

def jacc(x):
    _, jc, _ = Fc(x)
    return jc[:6]

def probsc(x):
    _, _, pc = Fc(x)
    return pc

x0c = np.ones((6,))/6.0 
boundsc = Bounds(np.zeros((6,)), np.ones((6,)))
solnc = minimize(objectivec, x0c, jac=jacc, method='SLSQP', bounds=boundsc, options={'ftol': 1e-16, 'disp': True}, 
        constraints=[{'type' : 'eq', 'fun' : lambda x : np.ones((1,6)).dot(x) - 1.0, 'jac' : lambda x : np.ones((1,6))}],
)
print('The solution for the canonical problem is given by p=q=(%0.4f, %0.4f, %0.4f, %0.4f, %0.4f, %0.4f).\nThis optimizer yields optimal value of  V(p,p)=%0.8f.' % (*solnc.x, solnc.fun))

Optimization terminated successfully.    (Exit mode 0)
            Current function value: 0.0012175833292183136
            Iterations: 12
            Function evaluations: 12
            Gradient evaluations: 12
The solution for the canonical problem is given by p=q=(0.2439, 0.1375, 0.1186, 0.1186, 0.1375, 0.2439).
This optimizer yields optimal value of  V(p,p)=0.00121758.