You write “In the case of the Newton-Pepys problem we need to find the probability $P(k)$ of throwing at least $k$ sixes with $6n$ dice,” and you go on to say: \[ P(X=k)=\binom{6n}{k} \left( \frac{1}{6} \right)^k \left(\frac{5}{6}\right)^{6n−k} \] But this is the probability $P(k)$ of throwing _exactly_ $k$ sixes with $6n$ dice, whereas what _we_ want is $P(X \geq k)$ for $6k$ dice. To sum these for all $k ~.. 6k$ would be tedious in the extreme and the solution given prior to this comment would be computationally preferable, for small $k$. @StevePowel Thanks for the heads up! Good point. I had a small typo in my equation. Just revised my previous annotation. @Micael Thanks for the corrections. There is still a small slip in the equation, because it should say $P(X \geq n) = \cdots $ and not $P(X = k) = \cdots$. It would be good to just mention that you employ the complement strategy, too. Samuel Pepys was a gambling man, and the Newton-Pepys problem is related to a wager he planned to make (he planned to stake 10 pounds, the equivalent to about $1500 today on a similar bet). In the correspondence he exchanged with Isaac Newton he wanted to know which of the following three scenarios was the most probable: * A. Throwing 6 dice and getting at least one six. * B. Throwing 12 dice and getting at least two sixes. * C. Throwing 18 dice and getting at least three sixes. Initially **Pepys thought that C was the most probable** but Newton showed him that ** A is the most probable scenario**. We can easily show that Newton was right and A is in fact the most probable scenario. Today this seems like a trivial problem but the field of [probability theory](https://en.wikipedia.org/wiki/Probability_theory) was in it's infancy. Probability table of obtaining **n** or more 6 when **n dice** are thrown. The **probability decreases with increasing n**.  Newton gave Pepys the correct answer to the problem, although his explanation was wrong. He imagined that A, B and C were tossing their dice in groups of six. For Newton 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 wrongly assumes that a group of dice does not produce more than one 6, so it does not actually correspond to the original problem posed by Pepys. If Newton's explanation was correct then A would always be more probable independently of p (the probability of a dice rolling to 6). If you start playing with biased dice (changing $p$) Newton's explanation no longer works. The figure below shows the probability of A, B and C winning the game. The ranking is initially unchanged ($P(A)>P(B)>P(C)$). As you can see as the die bias increases, the ranking of the games inverts, and with a highly biased die the games where more dice are rolled are more likely to win ($P(C)>P(B)>P(A)$), thus Newton's explanation no longer holds true.  Am I correct in saying that De Moivre's approximation here is just a special case of the normal approximation to the binomial distribution? This would explain why $P(X=Np) = \frac{1}{2\pi\sigma}$. @RobertAndrewMartin You are correct. In this specific case de Moivre–Laplace theorem shows that the probability function of the number of "successes" observed in a series of n independent throws converges to the probability density function of the normal distribution. You can learn more about it here: [De Moivre–Laplace theorem](https://en.wikipedia.org/wiki/De_Moivre–Laplace_theorem) Harvard Professor Joe Blitzstein introducing the Newton-Pepys problem (Harvard class Stats110): [](https://www.youtube.com/watch?v=P7NE4WF8j-Q&feature=youtu.be&t=17m47s) There is a problem with this part of Newton's argument. He assumes that you can only get one 6 per throw, which does not correspond to this problem. Assume that Peter and James throw groups of six dice. James does not need to throw the dice twice to win. He can win with the first throw by getting two or more 6. Here Newton makes the **wrong argument that you cannot get more than one 6 with one group of dice.** [Samuel Pepys](https://en.wikipedia.org/wiki/Samuel_Pepys "samuel pepys"), was a famous diarist of 17th Century London. He kept a **his private diary from 1660 until 1669 which is one of the most important sources for the English Restoration period**. Pepys was born in modest circumstances and educated as a scholarship student at Cambridge. He was a naturally curious person and undertook at the age of 29 to learn arithmetic. He found the multiplication tables particularly challenging, and used to wake up early and stay up late to study them. During the restoration of the Stuart monarchy, **Pepys was one of the most influential men in England**. As Secretary of the Affairs of the Admiralty he presided over the recovery of the British Navy and helped make it dominant over those of France and the Netherlands, which had dealt Britain an ignominious naval defeat in 1667. **Despite his lack of science training as a student, he had natural curiosity and became quite interested in science.**  A binomial distribution with parameters $n$ and $p$ is the discrete probability distribution of the number of successes in **a sequence of n independent experiments**, with **success probability p**: $$P(X=k)={\binom {n}{k}}p^{k}(1-p)^{n-k}$$ for k = 0, 1, 2, ..., n, where: $${\binom {n}{k}}={\frac {n!}{k!(n-k)!}}$$ is the binomial coefficient. In the case of the Newton-Pepys problem we need to find the probability **P(k) of throwing at least n sixes with 6n dice**. By applying the binomial distribution to this problem we can write: $$P(X=k)=1 - \sum_{k=0}^{n-1} {\binom {6n}{k}}\left(\frac{1}{6}\right)^{k}\left(\frac{5}{6}\right)^{6n-k}$$  Figure: Probability of rolling at least *k* six when we roll *6k* dice. The probability decreases as k increases. Newton and Pepys exchanged a series of 6 letters in late 1693. Pepys in London and Newton in Cambridge discussed a problem about gambling odds. Here are the links to the letters Pepys and Newton exchanged: You can read the **transcripts of their correspondence** here (chronological order): [Pepys to Newton - 11/22/1693](https://webspace.yale.edu/chem125/125/history99/2Pre1800/SPepysINewton/PepyNewtonPDFs/431P2N_112293.pdf) [Newton to Pepys - 11/26/1693](https://webspace.yale.edu/chem125/125/history99/2Pre1800/SPepysINewton/PepyNewtonPDFs/432N2P_112693.pdf) [Pepys to Newton - 12/9/1693](https://webspace.yale.edu/chem125/125/history99/2Pre1800/SPepysINewton/PepyNewtonPDFs/433P2N_12993.pdf) [Newton to Pepys - 12/16/1693](https://webspace.yale.edu/chem125/125/history99/2Pre1800/SPepysINewton/PepyNewtonPDFs/434N2P_121693.pdf) [Pepys to Newton - 12/21/1693](https://webspace.yale.edu/chem125/125/history99/2Pre1800/SPepysINewton/PepyNewtonPDFs/435P2N_122193.pdf) ### The Problem * A. Throwing 6 dice, and getting at least one 6. * B. Throwing 12 dice, and getting at least two 6. * C. Throwing 18 dice, and getting at least three 6. In the case of the combinatorial approach one computes the number of favorable outcomes and divide by the number of possible outcomes to get the probability of each scenario. Each dice has 6 faces so if you throw $N$ dice we have a total of $6^N$ possible outcomes. Instead of computing the probability of each event X we compute the complement probability. Throwing dice are independent events and one can write: $\bar{A}$ is the event of not getting a 6 when we throw the 6 dice. $D_i$ the event of not throwing a six for a given throw i. For the first scenario we want: \begin{eqnarray} \bar{A} = D_1 \cap D_2 ... \cap D_6 \end{eqnarray} These events are independent and thus one can write that the probability of not getting a six for i throws is: \begin{eqnarray*} P(\bar{A}) = \Pi \, P(D_i) \end{eqnarray*} The probability of *scenario A* then becomes: \begin{eqnarray} P(A) &=& 1 - P(\bar{A}) = 1 - \left(\frac{5}{6}\right)^6 \nonumber \\ &=& 0.665 \label{res1} \end{eqnarray} For scenarion B we do the same considerations as for scenario A. Let's denote $B_1$ the event of not getting any 6 and $B_2$ the event of getting exactly one 6. \begin{eqnarray*} \bar{B} = B_1 \cup B_2 \end{eqnarray*} These are disjoinct events and we can then write that: \begin{eqnarray} P(B) = 1 - P(\bar{B}) = 1 - ( P(B_1) + P(B_2) ) \end{eqnarray} We can compute the probability of $B_1$ from equation 1: \begin{eqnarray*} P(B_1) = \Pi \, P(D_i) = \left(\frac{5}{6}\right)^{12} \end{eqnarray*} For $B_2$ we need to compute the probability of getting one 6 and 11 non-6. We have 12 dice thus the probability of throwing one 6 is $12 1/6$. We have 11 dice left that cannot roll to a 6 and thus: \begin{eqnarray*} P(B_2) = 12 \left(\frac{1}{6}\right) \left(\frac{5}{6}\right)^{11} \end{eqnarray*} The probability of event B can be written: \begin{eqnarray} P(B) &=& 1 - P(\bar{B}) = 1 - \left( \left(\frac{5}{6}\right)^{12} + 2 \left(\frac{5}{6}\right)^{11} \right) \nonumber \\ &=& 1 - \frac{17}{6} \left(\frac{5}{6}\right)^{11} \\ &=& 0.619 \label{res2} \end{eqnarray} For scenarion C we consider the outcomes of B and we need to consider one additional type of outcome. Let's denote $C_3$ the event of getting exactly two 6. \begin{eqnarray*} \bar{C} = C_1 \cup C_2 \cup C_3 \end{eqnarray*} These are disjoinct events and we can then write that: \begin{eqnarray} P(C) = 1 - P(\bar{C}) = 1 - ( P(C_1) + P(C_2) + P(C_3) ) \end{eqnarray} For $C_1$ (no 6): \begin{eqnarray*} P(C_1) = \left(\frac{5}{6}\right)^{18} \end{eqnarray*} For $C_2$ (exactly one 6): \begin{eqnarray*} P(C_2) = 18 \left(\frac{1}{6}\right) \left(\frac{5}{6}\right)^{17} \end{eqnarray*} For $C_3$ we need to consider the events where we get exactly two 6. With 18 dice we have $\frac{18*17}{2}$ ways of getting exactly two sixes and we do not want to roll any 6 on the remaining 16 dice. \begin{eqnarray*} P(C_3) = \frac{18*17}{2} \left(\frac{1}{6}\right)^{2} \left(\frac{5}{6}\right)^{16} \end{eqnarray*} One can now write the total probability for scenarion C as: \begin{eqnarray} P(C) &=& 1 - \left( \left(\frac{5}{6}\right)^{18} + 3 \left(\frac{5}{6}\right)^{17} + 17/4 \left(\frac{5}{6}\right)^{16} \right) \nonumber \\ &=& 1 - \left( \left(\frac{5}{6}\right)^{18} + 3 \left(\frac{5}{6}\right)^{17} + 17/4 \left(\frac{5}{6}\right)^{16} \right) \nonumber \\ &=& 0.597 \label{res3} \end{eqnarray} From \ref{res1}, \ref{res2} and \ref{res3} we could easily determine that: \begin{eqnarray*} P(A) > P(B) > P(C) \end{eqnarray*}