Hitting a prime in 2.43 dice rolls (on average)
What is the number of rolls of fair 6-sided dice until the first time the total sum of all rolls is a prime? We compute the expectation and the variance of this random variable up to an additive error of less than 10^{-4}. This is a solution to a puzzle suggested by DasGupta (2017) in the Bulletin of the Institute of Mathematical Statistics, where the published solution is incomplete. The proof is simple, combining a basic dynamic programming algorithm with a quick Matlab computation and basic facts about the distribution of primes.
\section{Description of the problem}

The following puzzle appears in the Bulletin of the Institute of
Mathematical Statistics \citep{D2017}:
Let $X_1,X_2,\ldots$ be independent uniform random variables
on the integers $1,2,\ldots,6$,
and define $S_n=X_1+\ldots+X_n$ for $n=1,2,\ldots$.
Denote by $\tau$ the discrete time in which $S_n$ first hits
the set of prime numbers $P$:
\begin{equation*}
\label{eq:aaa}
\tau=\min\left\{n\geq 1: S_n \in P\right\}.
\end{equation*}

The contributing Editor \citep{D2017} provides a lower bound of
$2.3479$ for the expectation $E(\tau)$ and mentions the following heuristic
approximation for it: $E\left(\tau\right)\approx 7.6.$
He also adds that it is not known whether or not $\tau$ has finite variance.

In this note we compute the value of $E(\tau)$ up to an additive error
of less than $10^{-7}$, showing  that it is much closer to the
lower bound mentioned above than to $7.6$. We also show that the variance
is finite and compute its value up to an additive error of less than
$10^{-4}$. It will be clear from the discussion that it is not
difficult to get better approximation for both
quantities by increasing the amount of computation performed.
The proof is very simple, and the note is written mainly
in order to illustrate the power of the combination of a simple
dynamic programming algorithm combined with a quick computer-aided
computation and basic facts about the distribution of primes in the study
of problems of this type.

Before describing the rigorous argument, we present, in Table 1 below,
the outcomes of Monte-Carlo simulations of the process.

\begin{table}[H]
\label{eq:t}
\caption{Monte-Carlo simulations}
\begin{center}
\begin{tabular}{lllllll}
  number of repetitions  & $mean(\tau)$ & $variance(\tau)$ & $\max(\tau)$\\
  \hline
  $10^6$     &  2.4316    & 6.2735      & 49\\
  $2*10^6$   &  2.4274    & 6.2572      &67\\
  $3*10^6$   &  2.4305    & 6.2372      &70\\
  $5*10^6$   &  2.4287    & 6.2418      &64\\
  $10^7$     &  2.4286     &6.2463       &65\\
\end{tabular}
\end{center}
\end{table}

We proceed with a rigorous computation of $E(\tau)$ and $Var(\tau)$
up to an additive error smaller than $1/10,000$. Not surprisingly,
this computation shows that
the simulations supply accurate values.

Note first that
\begin{equation}
\label{eq:E}
E\left(\tau\right)=\sum_{k\geq 1}P(\tau\geq k),\,\,\,\,\,\,\,
E\left(\tau^2\right)=\sum_{k\geq 1}(2k-1)P(\tau\geq k).
\end{equation}

We next apply
dynamic programming to compute $p(k)\equiv P(\tau\geq k)$ exactly
for every $k$ up to $1000$, and then provide an upper bound for the
sums of the remaining terms.
\smallskip

\noindent
{\bf Dynamic-Programming Algorithm}
\smallskip

For each integer $k\geq 1$ and each non-prime $n$ satisfying
$k\leq n\leq 6k$, let $p(k,n)$ denote the probability that
$X_1+\ldots+X_k=n$ and that for every $i<k$, $X_1+\ldots+X_i$ is non prime.
Fix a parameter $K$ (in our computation we later take
$K=1000$.)
By the definition of $p(k,n)$ and the rule of total probability
we have the following dynamic programming algorithm for computing
$p(k,n)$ precisely for all $1 \leq k \leq K$ and $k \leq n \leq 6k$.
\begin{enumerate}
\item[1.]  $p(1,1)=p(1,4)=p(1,6)=1/6.$
\item[2.]
For $k=2,\ldots,K$ and any non-prime $n$ between $k$ and $6k$
\begin{equation*}
p(k,n)=\frac{1}{6}\sum_{i}p(k-1,n-i),
\end{equation*}
where the sum ranges over all $i$ between $1$ and $6$ so that
$n-i$ is non-prime.
\end{enumerate}

Denote by $E_{K}$ and $Var_{K}$ the estimators (lower bounds)
of $E(\tau)$ and $Var(\tau)$
based on the values of $p(k)$ for the first $K$ values,
which we obtain as follows:
\begin{align*}
&
E\left(\tau\right)=E_{K}+R_{K},\,\,\,\,
E\left(\tau^2\right)=E^{(2)}_{K}+R^{(2)}_{K},
\end{align*}
where $E_K=\sum_{k=1}^{K}p(k)$, $R_K=\sum_{k\geq K+1}p(k)$,
$E^{(2)}_K=\sum_{k=1}^{K}(2k-1)p(k)$,
$R^{(2)}_K=\sum_{k\geq K+1}(2k-1)p(k)$.
We also have
\begin{align*}
&
Var{(\tau)}=Var_K+RV_{K},
\end{align*}
where $Var_K=E^{(2)}_{K}-(E_K)^2$, $RV_{K}=R^{(2)}_{K}-2E_K R_K-(R_K)^2$.

Applying the dynamic-programming algorithm in Matlab, with
an execution time of less than $5$ seconds, we obtain
$$E_{1000}=2.428497913693504,\,\,\,\, Var_{1000}= 6.242778668279075.$$

\noindent
{\bf Bounding the remainders}
\smallskip

It remains to show that each of the sums of the remaining terms is
bounded by $10^{-4}$. To this end, we prove the following simple
result by induction on $k$.

\begin{pro}
For every $k$ and every non-prime $n$,
\begin{equation}
\label{eq:cl1}
p(k,n)<\frac{1}{3}\left(\frac{5}{6}\right)^{\pi(n)},
\end{equation}
where $\pi(n)$ is the number of primes smaller than $n$.
\end{pro}

\begin{proof}
Note first that \eqref{eq:cl1} holds for $k=1$, as $1/6=p(1,6)<(1/3) (5/6)^3$,
$1/6=p(1,4)<(1/3)(5/6)^2$  and $1/6=p(1,1)<(1/3)(5/6)^{0}$, with
room to spare.
Assuming the inequality holds for $k-1$ (and every relevant $n$)
we prove it for $k$. Suppose there are $q$ primes in the set
$\left\{n-6,\ldots,n-1\right\}$, then $\pi(n-i)\geq \pi(n)-q$
for all non-prime $n-i$ in this set. Thus, by the induction hypothesis
\begin{align*}
&
p(k,n) \leq \frac{1}{6}(6-q)\frac{1}{3}
\left(\frac{5}{6}\right)^{\pi(n)-q}
\leq \left(\frac{5}{6}\right)^q \frac{1}{3}
\left(\frac{5}{6}\right)^{\pi(n)-q}
=\frac{1}{3} \left(\frac{5}{6}\right)^{\pi(n)}.
\end{align*}
\end{proof}

By the prime number theorem (cf.,  e.g., \cite{HW2008}),
for every $n>1000$
$\pi(n) > 0.9 \frac{n}{\ln n}$ (again, with room to spare).
Therefore, by the above estimate we get:

\begin{cor}
\label{eq:cor1}
For every $k>1000$ and every non-prime $n (n\geq k)$,
\begin{equation}
\label{eq:cl2}
p(k,n)<\frac{1}{3}\left(\frac{5}{6}\right)^{0.9 \frac{n}{\ln n}}.
\end{equation}
\end{cor}

Using the above, we can now bound the sums of the remaining terms.

\begin{align*}
&
R_{1000}\equiv\sum_{k>1000}P(\tau\geq k)
=\sum_{k>1000} \sum_{\left\{n:\,\, k \leq n \leq 6k\right\}}
p(k,n) \nonumber\\
&
<
 \sum_{k >1000} \sum_{\left\{n:\,\, k \leq n \leq 6k\right\}}
\frac{1}{3} \left(\frac{5}{6}\right)^{0.9 n/ \ln n}
=\sum_{n\geq 1001}\sum_{k=\max(1001,n/6)}^{n}\frac{1}{3}
\left(\frac{5}{6}\right)^{0.9 n/ \ln n}\nonumber \\
 &
 <
 \sum_{n\geq 1000}(n-1000)\frac{1}{3}
\left(\frac{5}{6}\right)^{0.9 n/ \ln n},
\end{align*}
where the first inequality is obtained from Corollary \ref{eq:cor1}.

Define $f(n)=(n-1000)\frac{1}{3} \left(\frac{5}{6}\right)^{0.9 n/ \ln n},$
where $n$ is an integer $\geq 1000$. It is easy to check that for
$n \geq 1000$ the function
$f(n)$ has a unique maximum at $n=1050$.
(To see it, it suffices to compute $f(n)$ precisely for all
$1000 \leq n \leq 1100$ and observe that for $n>1100$
$f(n)$ is far smaller than $f(1050)$.)
It is also easy to check that for any $n\geq1050$ $f(n+13\ln n)/f(n)<1/2$.
Therefore,

\begin{align}
\label{eq:RE}
&
R_{1000}<\sum_{n\geq 1000}(n-1000)\frac{1}{3}
\left(\frac{5}{6}\right)^{0.9 n/ \ln n}=\sum_{n\geq 1000}f(n) \nonumber\\
&<
50f(1050)+2(13\ln 1050) f(1050)<7\cdot10^{-8}.
\end{align}

Similarly
$$
R^{(2)}_{1000}\equiv\sum_{k>1000}(2k-1)P(\tau\geq k)
=\sum_{k>1000} \sum_{\left\{n:\,\, k \leq n \leq 6k\right\}} (2k-1)p(k,n)
$$
$$
<\sum_{k>1000} \sum_{\left\{n:\,\, k \leq n \leq 6k\right\}}
(2k-1)\left(\frac{5}{6}\right)^{0.9 n/ \ln n}
\leq
\sum_{n\geq 1001}\frac{1}{3} \left(\frac{5}{6}\right)^{0.9 n/ \ln n}
\sum_{k=\max(1001,n/6)}^{n}(2k-1)
$$
$$
\leq
\sum_{n\geq 1000}\frac{1}{3} \left(\frac{5}{6}\right)^{0.9 n/ \ln n}
(n^2-1000^2),
$$
where the first inequality is obtained from Corollary \ref{eq:cor1}.

Denote by $g(n)=(n^2-1000^2)\frac{1}{3}
\left(\frac{5}{6}\right)^{0.9 n/ \ln n},$
where $n$ is an integer $\geq 1000$. For $n \geq 1000$ the function
$g(n)$ has an unique maximum at $n=1051$
and also for any $n\geq1051$ $g(n+13\ln n)/g(n)<1/2$.
Therefore,

\begin{align}
\nonumber
&
R^{(2)}_{1000}<\sum_{n\geq 1000}(n^2-1000^2)\frac{1}{3}
\left(\frac{5}{6}\right)^{0.9 n/ \ln n}=\sum_{n\geq 1000}g(n)\\
&
<
51g(1051)+2(13\ln 1051)g(1051)< 3.1\cdot 10^{-5}.
\end{align}

Therefore, the error $RV_{1000}$ in the variance estimation based on
the first $1000$ values
is  below $1/10,000$.

\iffalse
\section*{Acknowledgements}
The research of YM was supported by grant no. 2020063 from the
United States--Israel Binational Science Foundation (BSF),
Jerusalem, Israel.
\fi


2022-09-19T02:07:05.000
2209
2209.07698
en
https://arxiv.org/abs/2209.07698
What is the number of rolls of fair 6-sided dice until the first time the total sum of all rolls is a prime? We compute the expectation and the variance of this random variable up to an additive error of less than 10^{-4}. This is a solution to a puzzle suggested by DasGupta (2017) in the Bulletin of the Institute of Mathematical Statistics, where the published solution is incomplete. The proof is simple, combining a basic dynamic programming algorithm with a quick Matlab computation and basic facts about the distribution of primes.
Probability (math.PR); Number Theory (math.NT); Other Statistics (stat.OT)
Hitting a prime in 2.43 dice rolls (on average)
Qwen/Qwen-72B
1. YES
2. YES