# What is the expected length of the largest run of heads if we make 1,000 flips?

Is there a way to calculate on average, the maximum amount of times we can expect a coin to land heads during 1,000 flips?

So the answer (and formula if one exists) I am looking for would be something like: during 1,000 flips we can expect a maximum run of 12 heads in a row.

• I think that the precise question you're asking is difficult to give an analytical solution for. There's a duplicate floating around somewhere on math.SE; if I can find it, I'll add a link. Some notion of the general picture, though, can be obtained by observing that one expects to reach a streak of $2$ about once every other flip; a streak of $3$ about once every fourth flip; a streak of $4$ about once every eighth flip; a streak of $5$ about once every sixteenth flip; and so on. – Brian Tung Aug 25 '15 at 17:17
• I should point out that in such streaks, the coin is not really deviating from its $50$-$50$ probability; it is, in fact, adhering to it. That's just the nature of random behavior. – Brian Tung Aug 25 '15 at 17:18
• Makes sense. And now that I think about it, your right, the more times the coin is flipped the closer it keeps getting to an actual 50-50 number. I am going to remove the last sentence from my questions as the coin is not deviating as you have pointed out. – Juan Velez Aug 25 '15 at 17:22
• One way to phrase your question is "What is the minimal $k$ so that the probability of getting a maximum of $k$ or less heads in a row out of 1000 flips is at least $1/2$?" Or at least $0.95$? Etc. Another way to phrase your question is "What is the expected length of the largest run of heads if we make 1000 flips?" Which version/phrasing of the question are you interested in? – user2566092 Aug 25 '15 at 17:36
• I went with your 2nd phrasing. thanks. – Juan Velez Aug 25 '15 at 17:55

Some intuition about what you can expect can be found here: Longest Run of Heads.

Let the random variable $L_n$ be the largest contiguous heads sequence in $n$ coin tosses (suppose the coin is biased with heads probability $p$).

In the paper you can find the following intuitive approximation to the expectation of $L_n$: denote by $N_k$ the expected number of heads sequences having length$\ge k$. Since each tails outcome is a possible beginning of a heads sequence(ignoring edges), the expected number of heads sequences with length$\ge 1$ is $N_1\approx n(1-p)p$. Similarly, for length$\ge 2$ sequences the expectation is $N_2\approx n(1-p)p^2$ and generally $N_k\approx p(1-p)p^k$. Now you can approximate the expectation of $L_n$ by the solution to $N_k=1$ and this yields: $L_n$

$\langle L_n\rangle\approx -\log_pn(1-p)=\log_\frac{1}{p}n(1-p)$.

Although this appears extremely loose, it gives you an idea about the asymptotic behaviour $\langle L_n\rangle$ (logarithmic growth).

More accurately, you have $\frac{L_n}{\log_\frac{1}{p}n} \rightarrow 1$ in probability, i.e. $\forall \epsilon>0 \lim\limits_{n\to\infty}\mathbb{P} \left(\left|\frac{L_n}{\log_{1/p}n}-1\right|>\epsilon \right)=0$.

You may want to look at the following plot I got a while back. For $1\le n\le 1000$ $\langle L_n\rangle$ is calculated using $1000$ trails of $n$ unbiased coin tosses:

• Nice graphic. It actually helped me realize that as you flip a coin more and more, the max number of head coin flips in a row you witness will increase gradually over time. It's also interesting to see that it will also start to flat out as the amount of flips increases. I'm assuming the plot will eventually start to look like a flat horizontal line with very little movement after a while. I'm guessing we could probably also calculate how many times on avg. you would have to flip a coin to get a certain number of heads in a row. mind blown – Juan Velez Aug 25 '15 at 20:02
• "I'm assuming the plot will eventually start to look like a flat horizontal line with very little movement after a while" - well, it grows logarithmically, so clearly the length of the longest sequence tends to infinity (although slower than you might have expected). As to your last comment, suppose you managed to calculate the probability $p_k$ of having a largest sequence of size $k$ (in $n$ tosses), then the number of $n$ tosses trials you need in order to witness a $k$ longest sequence is geometrically distributed and it's expectation is $\frac{1}{p_k}$. – Ariel Aug 25 '15 at 20:24
• @Ariel: yes, and yet in another sense it does tend to a flat line because the slope approaches zero! Such is the way of limits. – wchargin Aug 26 '15 at 1:12
• I like this answer because it gives a justification about the growth rate, but it is based on approximation twice. Perhaps the paper you cited explains why this approach should be a particularly good approximation of the expected value of the maximum length run of heads. I would note, that even without the justification, at the very least, SOME kind of logarithmic growth should be expected because the probability of getting $k$ heads in a row for a given $k$ flips is $1/p^k$. Anyway, kudos on finding a paper that talks about it and you gave a nice accessible explanation. +1 – user2566092 Aug 26 '15 at 17:01
• I only meant to provide some intuition on $\langle L_n\rangle$, I think it's useful, if only to make sure that the values you're calculating are making sense. In any case, I voted you up because you gave a nice efficient way of computing the values of $P_k$. – Ariel Aug 26 '15 at 18:19

If you can calculate the probability of getting no more than $k$ heads in a row, call it $P_k$, then the probability that the maximum number of heads in a row is exactly $k$ is $P_k - P_{k-1}$. Then your expectation is

$$\sum_{k=1}^{1000} k(P_k - P_{k-1}) = 1000P_{1000} - \sum_{k=0}^{999} P_k$$

So it suffices to calculate all the $P_k$. Fortunately there is an easy way to do this with a recurrence on the total number of flips $n$. Let $Q_{k,m}(n)$ be the number of ways to have $n$ flips so that there are no more than $k$ heads in a row, and the flips end with a sequence of $m$ heads in a row. Then

$$Q_{k,0}(n+1) = \sum_{m=0}^k Q_{k,m}(n)$$

and

$$Q_{k,m}(n+1) = Q_{k,m-1}(n)$$ for all $m$ between $1$ and $k$, and your desired $P_k$ (as a function of $n = 1000$) is $P_k = (1/2^{1000}) \sum_{m=0}^k Q_{k,m}(1000)$.

You can calculate the $Q$ values in order of increasing $n$ up to 1000, fairly quickly, and then $P_k$ is as described. With a computer program code this would compute all $P_k$, and hence your expectation, in a fraction of a second.

Note: Here I'd like to present some information of an instructive and detailed treatment of this problem provided in Analytic Combinatorics by P. Flajolet and R. Sedgewick.

At first the authors derive a generating function very similar to the answer of @CountIblis. But the exciting thing (for me) is, that they also provide a rather detailed in depth analysis to find the asymptotic behaviour. It's great to see the interplay of many different facts and techniques in order to derive a precise asymptotic estimation.

Generating function:

We consider words (or strings) over a binary alphabet $\mathcal{A}=\{0,1\}$. A sequence $\mathop{SEQ}(0)=\{\varepsilon, 0,00,000,\ldots\}$ contains the set of all finite words containing $0$'s together with the empty word $\varepsilon$. The corresponding generating function is $\frac{1}{1-z}=1+z+z^2+\ldots$, the power of $z^n$ indicating the word with $n$ zeros and the coefficient $1$ giving the number of occurrences of this word.

All binary words can be characterized by all words starting with zero or more $0$'s followed by zero or more blocks each starting with a $1$ and followed by zero or more $0$'s. In symbols and as generating function we derive \begin{align*} \mathcal{W}\cong \mathop{SEQ}(0)\mathop{SEQ}(1\mathop{SEQ}(0))\qquad\Longrightarrow\qquad W(z)=\frac{1}{1-z}\frac{1}{1-z\frac{1}{1-z}}\tag{1} \end{align*} Observe that $W(z)$ reduces to $W(z)=\frac{1}{1-2z}=1+2z+4z^2+\ldots$ reflecting the fact that there are $2^n$ binary words of length $n$. The representation (1) is convenient for further development.



Considering longest runs of zeros, we denote with ${\mathop{SEQ}}^{<k}(0)$ the set of binary words with length less than $k$. The corresponding generating function is $\frac{1-z^{k}}{1-z}$. With respect to (1) all binary words having longest run of zeros less than $k$ are \begin{align*} \mathcal{W}^{\langle k\rangle}\cong {\mathop{SEQ}}^{<k}(0)\mathop{SEQ}(1{\mathop{SEQ}}^{<k}(0))\qquad\Longrightarrow \end{align*} \begin{align*}\qquad W(z)^{\langle k\rangle}&=\frac{1-z^k}{1-z}\frac{1}{1-z\frac{1}{1-z^k}}\tag{2}\\ &=\frac{1-z^k}{1-2z+z^{k+1}}\tag{3} \end{align*}

In the following asymptotic analysis we will use both representations (2) and (3).

The derivation of the generating function was the purely formal aspect without any needs for analytical means. It can be found in section I.4.1 of the book. But now we carry on on a somewhat deeper level. The treatment of this problem is part of chapter V: Applications of rational and meromorphic asymptotics.

Expectation value of longest run: We claim according to Proposition V.1

The longest run parameter $L$ taken over the set of binary words of length $n$ (endowed with the uniform distribution) satisfies the uniform estimate

\begin{align*} \mathbb{E}_n(L)=\lg(n)+\frac{\gamma}{\log(2)}-\frac{3}{2}+P(\lg n)+O\left(\frac{\log^2n}{\sqrt{n}}\right)\tag{4} \end{align*} where $P$ is a continuous periodic function whose Fourier expansion is given by \begin{align*} P(\omega)=-\frac{1}{\log 2}\sum_{k\in \mathbb{Z}\backslash\{0\}}\Gamma\left(\frac{2ik\pi}{\log 2}\right)e^{-2ik\pi\omega}\tag{5} \end{align*}

The oscillating function $P(\omega)$ is found to have tiny fluctuations, of the order of $10^{-6}$.

The symbol $\lg(n)\equiv\log_2(n)$ is the binary logarithm and $\gamma$ is the Euler-Mascheroni constant

Note: The expression (4) is a good approximation. We obtain for $n=1000$ the expectation value \begin{align*} \mathbb{E}_{1000}(L)=\lg(1000)+\frac{\gamma}{\log(2)}-\frac{3}{2}\doteq 9.29853\tag{6} \end{align*} and the first Fourier coefficient has amplitude: \begin{align*}\frac{\left|\Gamma\left(\frac{2i\pi}{\log 2}\right)\right|}{\log 2}\doteq 7.86\cdot 10^{-7}\end{align*} and since the fluctuations of the Fourier expansion is of order $10^{-6}$, the value stated at (6) is rather accurate.

Note, that the value $9.29853$ at $n=1000$ nicely matches with the $graph$ presented in the answer from @Ariel and the first order approximation $\lg(n)$ is the one provided in the answer by @leonbloy.

How to find the asymptotics:

The authors provide a detailed proof subdivded in four steps:

• locate the dominant pole

• estimate the corresponding contribution

• separate the dominant pole from the other poles in order to derive constructive error terms

• finally approximate the main quantitites of interest



Step 1: Locate the dominant pole

The dominant pole is the one nearest to the origin, responsible for the main contribution of the asymptotic approximation. By the first form (2) of the OGF $\mathcal{W}^{\langle k\rangle}$ the dominant pole $\rho_k$ is a root of the equation \begin{align*} 1=s(\rho_k)\qquad\text{where}\qquad s(z)=\frac{z(1-z^k)}{1-z} \end{align*} We consider $k>2$. Since $s(z)$ is an increasing polynomial for $z\geq 0$, we conclude from $s(0)=0, s(\frac{1}{2})<1, s(1)=k$, the root $\rho_k$ must lie in the open interval $\left(\frac{1}{2},1\right)$. In fact it's easy to verify that $s(0.6)>1$, hence the first estimate \begin{align*} \frac{1}{2}<\rho_k<\frac{3}{5}\qquad\qquad(k>2)\tag{7} \end{align*} We can now derive precise estimates of $\rho_k$ by bootstrapping, an iteration technique for approaching a fixed point. We write $1=s(z)$ in the form \begin{align*} z=\frac{1}{2}\left(1+z^{k+1}\right) \end{align*} and by using (7) we derive \begin{align*} \frac{1}{2}\left(1+{\left(\frac{1}{2}\right)}^{k+1}\right)<\rho_k <\frac{1}{2}\left(1+{\left(\frac{3}{5}\right)}^{k+1}\right) \end{align*} Thus, $\rho_k$ is exponentially close to $\frac{1}{2}$ and further iteration yields \begin{align*} \rho_k=\frac{1}{2}+\frac{1}{2^{k+2}}+O\left(\frac{k}{2^{2k}}\right)\tag{8} \end{align*}

Note, that $\rho_k$ corresponds to the calculation of the pole $p$ with $n=k+1$ from the answer of @CountIblis.

Step 2: Contribution from the dominant pole

The (expected) main contribution of the asymptotic estimation is the contribution from the dominant pole $\rho_k$. We calculate the residue of $W(z)^{\langle k\rangle}(z)=\frac{g(z)}{h(z)}$ and obtain \begin{align*} R_{n,k}&:=-\mathop{Res}\left[\frac{W(z)^{\langle k\rangle}(z)}{z^{n+1}};z=\rho_k\right]=-\frac{g(\rho_k)}{h^{\prime}(\rho_k)}\rho_k^{-n-1}\\ &=\frac{1-\rho_k^k}{2-(k+1)\rho_k^k}\rho_k^{-n-1} \end{align*}

This quantitiy is expected to provide the main approximation to the coefficients of $W(z)^{\langle k\rangle}(z)$ as $n\rightarrow \infty$.

Step 3: Separation from the subdominant poles

Next we want to show that $\rho_k$ is the only zero of the denominator of $W(z)^{\langle k\rangle}(z)$ within the disc $|z|\leq \frac{3}{4}$. We apply Rouché's theorem and regard the denominator of $W(z)^{\langle k\rangle}(z)$ \begin{align*} 1-2z+z^{k+1} \end{align*} as sum of two functions $f(x)=1-2z$ and $g(z)=z^{k+1}$. The term $f(z)$ has on the circle $|z|=\frac{3}{4}$ a modulus that varies between $\frac{1}{2}$ and $\frac{5}{2}$, while the term $g(z)\leq \frac{27}{64}$ for any $k\geq 2$.

Thus on the circle $|z|=\frac{3}{4}$, one has $|g(z)|<|f(z)|$, so that $f(z)$ and $f(z)+g(z)$ have the same number of zeros inside the circle. Since $f(z)$ admits $z=\frac{1}{2}$ as only zero there, the denominator must also have a unique root in $|z|\leq \frac{3}{4}$, and that root must coincide with $\rho_k$.

Similar arguments are further applied to derive a representation $q_{n,k}$ of the probability that the longest run in a random word of length $n$ is less than $k$. The following main estimate is obtained for $k>2$

\begin{align*} q_{n,k}:=\mathbb{P}_n(L<k)=\frac{1-\rho_k^k}{1-(k+1)\frac{\rho_k^k}{2}}\left(\frac{1}{2\rho_k}\right)^{n+1} +O\left(\left(\frac{2}{3}\right)^n\right)\tag{9} \end{align*} which holds uniformly with respect to $k$.

Step 4: Final approximations

Note: I present this last step without going into details (which are partly beyond my skills). Remarkable is (at least for me) the kind of art and mastery to find an appropriate partitioning of the domain in order to identify tail regions and central regions, the latter of them relevant for the main approximation. I would appreciate if someone could provide fruitful hints regarding this theme.

The first part of this step is to estimate the tail inequalities comprising the regions which do not contribute to the main approximation.

They are derived by applying the representation (8) of $\rho_k$ in the main approximation (9). Thus, for asymptotic purposes, only a relative small region around $\ln n$ needs to be considered.

Regarding the central regime, for $k=\lg n+x$ and $x\in[-\frac{1}{4}\lg n,\lg n]$, again with the approximation (8) of $\rho_{k}$ and related quantities, one finds \begin{align*} (2\rho_k)^{-n}=\exp\left(-\frac{n}{2^{k+1}}+O(kn2^{-2k})\right) =e^{-\frac{n}{2^{k+1}}}\left(1+O\left(\frac{\log n}{\sqrt{n}}\right)\right) \end{align*} and it follows (uniformly in $k$) \begin{align*} q_{n,k}=e^{-\frac{n}{2^{k+1}}}\left(1+O\left(\frac{\log n}{\sqrt{n}}\right)\right)\tag{11} \end{align*}

The following mean estimate is derived from the fact that the distribution quickly decays at values away from $\lg n$ by (10) while it satisfies equation (11) in the central region. \begin{align*} \mathbb{E}_n(L):=\sum_{h\geq 1}[1-\mathbb{P}_n(L<h)]=\Phi\left(\frac{n}{2}\right)-1+O\left(\frac{\log ^2 n}{n}\right) \end{align*} with \begin{align*} \Phi(x):=\sum_{h\geq 0}\left[1-e^{-\frac{x}{2^h}}\right] \end{align*} We consider the three cases $h<h_0,h\in[h_0,h_1]$, and $h>h_1$, with $h_0=\lg x-\lg \lg x$ and $h_1=\lg x+\lg \lg x$, where the general term is (respectively) close to $1$, between $0$ and $1$, and close to $0$. By summing, one finds elementarily \begin{align*} \Phi(x)=\lg x+O(\lg \lg x) \qquad \qquad x\rightarrow\infty \end{align*}

According to the authors the method of choice for precise asymptotics is to treat $\Phi(x)$ as a harmonic sum and apply Mellin transform techniques. The Mellin transform of $\Phi(x)$ is \begin{align*} \Phi^{\star}(s):=\int_0^{\infty}\Phi(x)x^{s-1}dx=\frac{\Gamma(s)}{1-2^s}\qquad\qquad \mathcal{R}(s)\in (-1,0) \end{align*}

They finally conclude:

The double pole of $\Phi^{\star}$ at $0$ and the simple poles at $s=\frac{2ik\pi}{\log 2}$ are reflected by an asymptotic expansion that involves a Fourier series \begin{align*} \Phi(x)=\lg x+\frac{\gamma}{\log 2}+\frac{1}{2}+P(\lg x)+O(x^{-1}) \end{align*} with $P(\omega)$ as stated in (5). \begin{align*} \qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\Box \end{align*}

You can also attack this problem using generating functions. If we give a sequence of n coin throws a weight of $x$, then the generating function for a single string of at most $n$ heads is:

$$G_n(x) = \frac{1-x^{n+1}}{1-x}$$

A general sequence of coin throws such that there are at most $n$ heads in any sequence of heads, can consist of any arbitrary number of such sequences, the generating function $f(x)$ is obtained by summing over all possible ways to join together such sequences with a singe tail inbetween:

$$f(x) = G_n(x) + G_n(x) x G_n(x) + G_n(x) x G_n(x) x G_n(x)+\cdots$$

We thus have:

$$f(x) = G_n(x)\sum_{k=0}^{\infty}x^kG_n^k(x) = \frac{G_n(x)}{1-x G_n(x)}=\frac{1-x^{n+1}}{1-2x+x^{n+2}}$$

The coefficient of $x^r$ can be written as:

$$c_r = \frac{1}{2\pi i}\oint_C\frac{1}{z^{r+1}}\frac{1-z^{n+1}}{1-2z+z^{n+2}}dz$$

where the contour $C$ is a small counterclockwise circle around the origin that avoids the poles of $f(z)$. Since the contour integral over a circle of radius $R$ with center the origin will tend to zero in the limit of $R\to\infty$, the sum of all the residues of the integrand equals zero, therefore $c_r$ is equal to minus the sum of all the resides at all the poles other than the one at the origin. The dominant contribution to $c_r$ for large $r$ comes from the pole near $z=\frac{1}{2}$. To evaluate it, one can write the zero of the denominator as $\frac{1}{2} + t$ expand to first order in $t$ and solve for $t$. The pole $p$ is then found to be approximately located at:

$$p \approx \frac{1}{2} + \frac{1}{2^{n+3} - 2 n -4}$$

This then yields:

$$c_r\approx 2^r \left(1+\frac{1}{2^{n+2} -n-2}\right)^{-(r+1)}\frac{1-p^{n+1}}{1-\frac{(n+2)}{2} p^{(n+1)}}$$

The probability for there being a sequence of at most $n$ heads is thus approximately:

$$P_n \approx \exp\left(-\frac{r+1}{2^{n+1}-n-2}\right)\frac{2^{n+1}-1}{2^{n+1}-\frac{n+2}{2}}$$

A quick approximation (works asymptotically) would be to model that the succesive runs as iid geometric variables, say $r_n \sim Geom(1/2)$, so $E(r_n)=2$, so we'd have approximately $N/2$ runs, $N/4$ being head runs. Then the problem would be equivalent to finding the maximum of $N/4$ iid geometric rv... which is not an easy problem. Using the (again, approximate) results of the linked answer, we'd get

$$E[\max_{N/4}\{ r_n\}] \approx \frac{1}{2} + \frac{1}{\log 2} H_{N/4}\approx \frac{1}{2} +\frac{\log(N/4) +\gamma}{\log 2} \approx \log_2(N)$$