Newton–Pepys problem

The Newton–Pepys problem is a probability problem concerning the probability of throwing sixes from a certain number of dice.

In 1693 Samuel Pepys and Isaac Newton corresponded over a problem posed to Pepys by a school teacher named John Smith. The problem was:

Pepys initially thought that outcome C had the highest probability, but Newton correctly concluded that outcome A actually has the highest probability.

The probabilities of outcomes A, B and C are:

${\displaystyle P(A)=1-\left({\frac {5}{6}}\right)^{6}={\frac {31031}{46656}}\approx 0.6651\,,}$

${\displaystyle P(B)=1-\sum _{x=0}^{1}{\binom {12}{x}}\left({\frac {1}{6}}\right)^{x}\left({\frac {5}{6}}\right)^{12-x}={\frac {1346704211}{2176782336}}\approx 0.6187\,,}$

${\displaystyle P(C)=1-\sum _{x=0}^{2}{\binom {18}{x}}\left({\frac {1}{6}}\right)^{x}\left({\frac {5}{6}}\right)^{18-x}={\frac {15166600495229}{25389989167104}}\approx 0.5973\,.}$

These results may be obtained by applying the binomial distribution (although Newton obtained them from first principles). In general, if P(n) is the probability of throwing at least n sixes with 6n dice, then:

${\displaystyle P(n)=1-\sum _{x=0}^{n-1}{\binom {6n}{x}}\left({\frac {1}{6}}\right)^{x}\left({\frac {5}{6}}\right)^{6n-x}\,.}$

As n grows, P(n) decreases monotonically towards an asymptotic limit of 1/2.

The solution outlined above can be implemented in R as follows:

for (s in 1:3) {          # looking for s = 1, 2 or 3 sixes
  n = 6*s                 # ... in n = 6, 12 or 18 dice
  q = pbinom(s-1, n, 1/6) # q = Prob( <s sixes in n dice )
  cat("Probability of at least", s, "six in", n, "fair dice:", 1-q, "\n")
}

Although Newton correctly calculated the odds of each bet, he provided a separate intuitive explanation to Pepys. He imagined that B and C toss their dice in groups of six, and said that A was most favorable because it required a 6 in only one toss, while B and C required a 6 in each of their tosses. This explanation assumes that a group does not produce more than one 6, so it does not actually correspond to the original problem.

A natural generalization of the problem is to consider n non-necessarily fair dice, with p the probability that each die will select the 6 face when thrown (notice that actually the number of faces of the dice and which face should be selected are irrelevant). If r is the total number of dice selecting the 6 face, then ${\displaystyle P(r\geq k;n,p)}$is the probability of having at least k correct selections when throwing exactly n dice. Then the original Newton–Pepys problem can be generalized as follows:

Let ${\displaystyle \nu _{1},\nu _{2}}$be natural positive numbers s.t. ${\displaystyle \nu _{1}\leq \nu _{2}}$. Is then ${\displaystyle P(r\geq \nu _{1}k;\nu _{1}n,p)}$not smaller than ${\displaystyle P(r\geq \nu _{2}k;\nu _{2}n,p)}$for all n, p, k?

Notice that, with this notation, the original Newton–Pepys problem reads as: is ${\displaystyle P(r\geq 1;6,1/6)\geq P(r\geq 2;12,1/6)\geq P(r\geq 3;18,1/6)}$?

As noticed in Rubin and Evans (1961), there are no uniform answers to the generalized Newton–Pepys problem since answers depend on k, n and p. There are nonetheless some variations of the previous questions that admit uniform answers:

(from Chaundy and Bullard (1960)):

If ${\displaystyle k_{1},k_{2},n}$are positive natural numbers, and ${\displaystyle k_{1}<k_{2}}$, then ${\displaystyle P(r\geq k_{1};k_{1}n,{\frac {1}{n}})>P(r\geq k_{2};k_{2}n,{\frac {1}{n}})}$.

If ${\displaystyle k,n_{1},n_{2}}$are positive natural numbers, and ${\displaystyle n_{1}<n_{2}}$, then ${\displaystyle P(r\geq k;kn_{1},{\frac {1}{n_{1}}})>P(r\geq k;kn_{2},{\frac {1}{n_{2}}})}$.

(from Varagnolo, Pillonetto and Schenato (2013)):

If ${\displaystyle \nu _{1},\nu _{2},n,k}$are positive natural numbers, and ${\displaystyle \nu _{1}\leq \nu _{2},k\leq n,p\in [0,1]}$then ${\displaystyle P(r=\nu _{1}k;\nu _{1}n,p)\geq P(r=\nu _{2}k;\nu _{2}n,p)}$.