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.$
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))
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))
