Asymptotics of probability of Newton-Pepys problem

Question

In Newton-Pepys problem one is interested in probability $p_n$ of getting at least $n$ sixes in $6 n$ independent throws of regular 6-sided die.

The number of sixes $S_m$ obtained in $m$ throws follows a Binomial distribution $\operatorname{Bin}(m,p)$, where $p=\frac{1}{6}$ is the probability of getting a six in a single throw. Thus $$ p_n = \Pr\left(S_{6n} \geqslant n\right) $$ Notice that $\mathsf{E}(S_{6n}) = 6 n p = n$, and variance $\mathbb{Var}(S_{6n}) = \sqrt{5 n/6}$, hence in the large $n$ limit $$ \lim_{n \to \infty} p_n = \lim_{n \to \infty} \Pr\left(\frac{S_{6n}-n}{\sqrt{5 n/6} }\geqslant 0\right) \stackrel{\text{CLT}}{=} \Pr(Z \geqslant 0) = \frac{1}{2} $$ In fact $p_n$ is monotonically decreasing sequence:

Q.: How can one find the large $n$ asymptotics of $p_n$?

Answer 1

By summing the series for $p_n$: $$ p_n = \Pr(S_{6n} \geqslant n) = \sum_{m=n}^{6n} \binom{6n}{m} p^m (1-p)^{6n-m} = \frac{1}{\operatorname{B}\left(n,5n+1\right)} \int_0^p t^{n-1} (1-t)^{5n} \mathrm{d}t $$ The maximum of the integrand is achieved at $t = \frac{n-1}{6n-1}$ which is $\frac{5}{6n-1}$ away from the upper limit of integration. Making a change of variables: $$ t = \frac{n-1}{6n-1}\left(1+ u \sqrt{\frac{5n}{(6n-1)(n-1)}} \right) $$ The integral will become: $$ p_n = \mathcal{N}_n \int_{-\infty}^{u_\max(n)} \mathrm{e}^{-u^2/2} \left(1 + c_3(n) u^3 + \cdots\right) \mathrm{d}u $$ for some quite unwieldy expressions of $\mathcal{N}_n$, $u_\max(n)$ and coefficients $c_k(n)$.

Doing all of this analysis with the help of CAS (Mathematica in my case) I got $$ p_n = \frac{1}{2} + \frac{7}{6\sqrt{15\pi}}\frac{1}{\sqrt{n}} \left(1 - \frac{113}{7!} \frac{1}{n} - \frac{15097}{7! \cdot 360} \frac{1}{n^2} - \frac{1321931}{2^{10} \cdot 3^7 \cdot 5^4} \frac{1}{n^3} + \mathcal{o}\left(n^{-3}\right)\right) $$ Here is the visualization of the error term:

Interestingly, the expression for $p_n$ satisfies a second order recurrence equation:

which can be used to easily obtain higher order expansion terms.

Answer 2

The probability $p_n$ has a representation in terms of regularized beta-function, which incidentally is the cumulative distribution function of beta distribution. $$ p_n = \Pr(S_{6n} \geqslant n) = I_{p}\left(n, 5n+1\right) = \Pr(X < p ) = \Pr\left(X < \frac{1}{6} \right)) $$ where $X$ follows $\mathcal{Be}\left(n,5n+1\right)$. It is well-known that $X \stackrel{d} = \frac{Y}{Y+Z}$ where $Y$ and $Z$ are independent random variables following gamma-distirbutions $\Gamma(n,1)$ and $\Gamma(5n+1,1)$ respectively. $$ \begin{eqnarray} p_n &=& \Pr\left(\frac{Y}{Y+Z} < \frac{1}{6}\right) = \Pr(5 Y < Z) = \Pr\left(5 (Y-n) < Z -5 n\right) \\ &=& \Pr\left(5 \frac{Y-n}{\sqrt{n}} < \frac{Z -5 n}{\sqrt{n}}\right) = \Pr\left(5 \frac{Y-n}{\sqrt{n}} < \frac{Z -5 n-1}{\sqrt{n}} + \frac{1}{\sqrt{n}}\right) \\ &=& \Pr\left(5 \frac{Y-n}{\sqrt{n}} < \frac{Z -5 n-1}{\sqrt{5n+1}} \sqrt{5 + \frac{1}{n}} + \frac{1}{\sqrt{n}}\right) \end{eqnarray} $$ Since mean and standard deviation of $Y$ equal $n$ and $\sqrt{n}$ and mean and standard deviation of $Z$ equal $5n+1$ and $\sqrt{5n+1}$.

Approximating $\frac{Y-n}{\sqrt{n}}$ and $\frac{Z -5 n-1}{\sqrt{5n+1}}$ with iid standard normal random variables: $$ p_n \approx \Pr\left(5 Z_1 - Z_2 \sqrt{5+\frac{1}{n}} < \frac{1}{\sqrt{n}}\right) = \Phi\left(\frac{1}{\sqrt{30 n + 1}}\right) $$ Unfortunately, the above approximation is not very accurate.