A serious in-depth problem of probabilities for flipping coins

Lets say I'm doing 10,000 flips of a coin. I would like to know the probability of how many flips it takes to get 4 or more consecutive heads in a row.

The count would work as the following, you would count one successive round of flips being only heads (4 heads or more). When a tails hits and breaks the streak of heads the count would start again from the next flip. This would then repeat for 10,000 flips.

I'd like to know the probability of not just 4 or more heads in a row but, 6 or more, and 10 or more. To clarify if a streak of 9 heads is achieved it would be tallied as 1 streak 4 or more (and/or 6 or more), not 2 separate streaks. For example if the coin came THTHTHTHHHHHH///THAHTHT....the count would be 13 and begin once again on the next tails.

Let's say the data comes out to be highly skewed to the right; the mean being 40 flips on average it takes to achieve a streak of 4 or more, and the distribution is u = 28. Obviously skewed.

I'm doing my best to find a way to make some sense out of descriptive data, except as of now I've found nothing.

I want to find some way of getting some sensible probability out of it. Like a normal curve where +/- 1 SD is 68% etc..I've looked into log normalizing and that's only really used for a parametric testing which is not my goal.

I have been told beta distributions but every suggestion I've had has been quite confusing. I've asked this question a year ago and got some insight but unfortunately I still do not have an answer. Thank you to any of you who have ideas.

• I probably should clarify a little. 1) I'm looking for away to make sense of descriptive data of the amount of consecutive heads above 4 in 1000 flips (similar to something like a normal curve probability +/-1 SD = 68%) from a skewed data set. 2) Recommended to use a beta distribution but ANY other suggestion would be great! – Dan Apr 16 '13 at 0:07
• Dan, I just noticed that your example set of heads and tails includes an "A". – Glen_b Mar 11 '14 at 2:54
• The edit you made are a big improvement, but we need to make some more changes. Where you say "and the distribution is u = 28", what exactly do you mean? Are you talking about the median? – Glen_b Mar 11 '14 at 2:59
• @Dan the beta could only possibly factor into this problem if you were using a Bayesian approach and estimating the probability of a heads, then applying that distribution (and its associated uncertainty) into the mathematical result of the problem you stated. – AdamO Mar 20 '14 at 16:37

If I have understood correctly, then the problem is to find a probability distribution for the time at which the first run of $n$ or more heads ends.

Edit The probabilities can be determined accurately and quickly using matrix multiplication, and it is also possible to analytically compute the mean as $\mu_-=2^{n+1}-1$ and the variance as $\sigma^2=2^{n+2}\left(\mu-n-3\right) -\mu^2 + 5\mu$ where $\mu=\mu_-+1$, but there is probably not a simple closed form for the distribution itself. Above a certain number of coin flips the distribution is essentially a geometric distribution: it would make sense to use this form for larger $t$.

The evolution in time of the probability distribution in state space can be modelled using a transition matrix for $k=n+2$ states, where $n=$ the number of consecutive coin flips. The states are as follows:

• State $H_0$, no heads
• State $H_i$, $i$ heads, $1\le i\le(n-1)$
• State $H_n$, $n$ or more heads
• State $H_*$, $n$ or more heads followed by a tail

Once you get into state $H_*$ you can't get back to any of the other states.

The state transition probabilities to get into the states are as follows

• State $H_0$: probability $\frac12$ from $H_i$, $i=0,\ldots,n-1$, i.e. including itself but not state $H_n$
• State $H_i$: probability $\frac12$ from $H_{i-1}$
• State $H_n$: probability $\frac12$ from $H_{n-1},H_n$, i.e. from the state with $n-1$ heads and itself
• State $H_*$: probability $\frac12$ from $H_n$ and probability 1 from $H_*$ (itself)

So for example, for $n=4$, this gives the transition matrix

$$X = \left\{ \begin{array}{c|cccccc} & H_0 & H_1 & H_2 & H_3 & H_4 & H_* \\ \hline H_0 & \frac12 & \frac12& \frac12& \frac12 &0 & 0 \\ H_1 & \frac12 &0 & 0 &0 & 0 & 0 \\ H_2 & 0 & \frac12 &0 &0 & 0 & 0 \\ H_3 & 0 & 0 &\frac12 &0 & 0 & 0 \\ H_4 & 0 & 0 &0 &\frac12 & \frac12 & 0 \\ H_* & 0 & 0 & 0&0 &\frac12 & 1 \end{array}\right\}$$

For the case $n=4$, the initial vector of probabilities $\mathbf{p}$ is $\mathbf{p}=(1,0,0,0,0,0)$. In general the initial vector has $$\mathbf{p}_i=\begin{cases}1&i=0\\0&i>0\end{cases}$$

The vector $\mathbf{p}$ is the probability distribution in space for any given time. The required cdf is a cdf in time, and is the probability of having seen at least $n$ coin flips end by time $t$. It can be written as $(X^{t+1}\mathbf{p})_k$, noting that we reach state $H_*$ 1 timestep after the last in the run of consecutive coin flips.

The required pmf in time can be written as $(X^{t+1}\mathbf{p})_k - (X^{t}\mathbf{p})_k$. However numerically this involves taking away a very small number from a much larger number ($\approx 1$) and restricts precision. Therefore in calculations it is better to set $X_{k,k}=0$ rather than 1. Then writing $X'$ for the resulting matrix $X'=X|X_{k,k}=0$, the pmf is $(X'^{t+1}\mathbf{p})_k$. This is what is implemented in the simple R program below, which works for any $n\ge 2$,

n=4
k=n+2
X=matrix(c(rep(1,n),0,0, # first row
rep(c(1,rep(0,k)),n-2), # to half-way thru penultimate row
1,rep(0,k),1,1,rep(0,k-1),1,0), # replace 0 by 2 for cdf
byrow=T,nrow=k)/2
X

t=10000
pt=rep(0,t) # probability at time t
pv=c(1,rep(0,k-1)) # probability vector
for(i in 1:(t+1)) {
#pvk=pv[k]; # if calculating via cdf
pv = X %*% pv;
#pt[i-1]=pv[k]-pvk # if calculating via cdf
pt[i-1]=pv[k] # if calculating pmf
}

m=sum((1:t)*pt)
v=sum((1:t)^2*pt)-m^2
c(m, v)

par(mfrow=c(3,1))
plot(pt[1:100],type="l")
plot(pt[10:110],type="l")
plot(pt[1010:1110],type="l")


The upper plot shows the pmf between 0 and 100. The lower two plots show the pmf between 10 and 110 and also between 1010 and 1110, illustrating the self-similarity and the fact that as @Glen_b says, the distribution looks like it can be approximated by a geometric distribution after a settling down period.

It's possible to investigate this behaviour further using an eigenvector decomposition of $X$. Doing so shows that the for sufficiently large $t$, $p_{t+1}\approx c(n)p_t$, where $c(n)$ is the solution of the equation $2^{n+1}c^n(c-1)+1=0$. This approximation gets better with increasing $n$ and is excellent for $t$ in the range from about 30 to 50, depending on the value of $n$, as shown in the plot of log error below for computing $p_{100}$ (rainbow colours, red on the left for $n=2$). (In fact for numerical reasons, it would actually be better to use the geometric approximation for probabilities when $t$ is larger.)

I suspect(ed) there might be a closed form available for the distribution because the means and variances as I have calculated them as follows

$$\begin{array}{r|rr} n & \text{Mean} & \text{Variance} \\ \hline 2 &7& 24& \\ 3 &15& 144& \\ 4 &31& 736& \\ 5 &63& 3392& \\ 6 &127& 14720& \\ 7 &255& 61696& \\ 8 &511& 253440& \\ 9 &1023& 1029120& \\ 10&2047& 4151296& \\ \end{array}$$

(I had to bump up the number up the time horizon to t=100000 to get this but the program still ran for all $n=2,\ldots,10$ in less than about 10 seconds.) The means in particular follow a very obvious pattern; the variances less so. I have solved a simpler, 3-state transition system in the past, but so far I'm having no luck with a simple analytic solution to this one. Perhaps there's some useful theory that I'm not aware of, e.g. relating to transition matrices.

Edit: after a lot of false starts I came up with a recurrence formula. Let $p_{i,t}$ be the probability of being in state $H_i$ at time $t$. Let $q_{*,t}$ be the cumulative probability of being in state $H_*$, i.e. the final state, at time $t$. NB

• For any given $t$, $p_{i,t}, 0\le i\le n$ and $q_{*,t}$ are a probability distribution over space $i$ , and immediately below I use the fact that their probabilities add to 1.
• $p_{*,t}$ form a probability distribution over time $t$. Later, I use this fact in deriving the means and variances.

The probability of being at the first state at time $t+1$, i.e. no heads, is given by the transition probabilities from states that can return to it from time $t$ (using the theorem of total probability). \begin{align} p_{0,t+1} &= \frac12p_{0,t} + \frac12p_{1,t} + \ldots \frac12p_{n-1,t}\\ &= \frac12\sum_{i=0}^{n-1}p_{i,t}\\ &= \frac12\left( 1-p_{n,t}-q_{*,t} \right) \end{align} But to get from state $H_0$ to $H_{n-1}$ takes $n-1$ steps, hence $p_{n-1,t+n-1} = \frac1{2^{n-1}}p_{0,t}$ and $$p_{n-1,t+n} = \frac1{2^n}\left( 1-p_{n,t}-q_{*,t} \right)$$ Once again by the theorem of total probability the probability of being at state $H_n$ at time $t+1$ is \begin{align} p_{n,t+1} &= \frac12p_{n,t} + \frac12p_{n-1,t}\\ &= \frac12p_{n,t} + \frac1{2^{n+1}}\left( 1-p_{n,t-n}-q_{*,t-n} \right)\;\;\;(\dagger)\\ \end{align} and using the fact that $q_{*,t+1}-q_{*,t}=\frac12p_{n,t} \implies p_{n,t} = 2q_{*,t+1}-2q_{*,t}$, \begin{align} 2q_{*,t+2} - 2q_{*,t+1} &= q_{*,t+1}-q_{*,t}+\frac1{2^{n+1}}\left( 1-2q_{*,t-n+1}+q_{*,t-n} \right) \end{align} Hence, changing $t\to t+n$, $$2q_{*,t+n+2} - 3q_{*,t+n+1} +q_{*,t+n} + \frac1{2^n}q_{*,t+1} - \frac1{2^{n+1}}q_{*,t} - \frac1{2^{n+1}} = 0$$

This recurrence formula checks out for the cases $n=4$ and $n=6$. E.g. for $n=6$ a plot of this formula using t=1:994;v=2*q[t+8]-3*q[t+7]+q[t+6]+q[t+1]/2**6-q[t]/2**7-1/2**7 gives machine order accuracy.

Edit I can't see where to go to find a closed form from this recurrence relation. However, it is possible to get a closed form for the mean.

Starting from $(\dagger)$, and noting that $p_{*,t+1}=\frac12p_{n,t}$, \begin{align} p_{n,t+1} &= \frac12p_{n,t} + \frac1{2^{n+1}}\left( 1-p_{n,t-n}-q_{*,t-n} \right)\;\;\;(\dagger)\\ 2^{n+1}\left(2p_{*,t+n+2}-p_{*,t+n+1}\right)+2p_{*,t+1} &= 1-q_{*,t} \end{align} Taking sums from $t=0$ to $\infty$ and applying the formula for the mean $E[X] = \sum_{x=0}^\infty (1-F(x))$ and noting that $p_{*,t}$ is a probability distribution gives \begin{align} 2^{n+1}\sum_{t=0}^\infty \left(2p_{*,t+n+2}-p_{*,t+n+1}\right) + 2\sum_{t=0}^\infty p_{*,t+1} &= \sum_{t=0}^\infty(1-q_{*,t}) \\ 2^{n+1} \left(2\left(1-\frac1{2^{n+1}}\right)-\;1\;\right) + 2 &= \mu \\ 2^{n+1} &= \mu \end{align} This is the mean for reaching state $H_*$; the mean for the end of the run of heads is one less than this.

Edit A similar approach using the formula $E[X^2] = \sum_{x=0}^\infty (2x+1)(1-F(x))\;$ from this question yields the variance. \begin{align} \sum_{t=0}^\infty(2t+1)\left (2^{n+1}\left(2p_{*,t+n+2}-p_{*,t+n+1}\right) + 2 p_{*,t+1}\right) &= \sum_{t=0}^\infty(2t+1)(1-q_{*,t}) \\ 2\sum_{t=0}^\infty t\left (2^{n+1}\left(2p_{*,t+n+2}-p_{*,t+n+1}\right) + 2 p_{*,t+1}\right) + \mu &= \sigma^2 + \mu^2 \\ 2^{n+2}\left(2\left(\mu-(n+2)+\frac1{2^{n+1}}\right)-(\mu-(n+1))\right) + 4(\mu-1) &\\+ \mu &= \sigma^2 + \mu^2 \\ 2^{n+2}\left(2(\mu-(n+2))-(\mu-(n+1))\right) + 5\mu &= \sigma^2 + \mu^2 \\ 2^{n+2}\left(\mu-n-3\right) + 5\mu &= \sigma^2 + \mu^2 \\ 2^{n+2}\left(\mu-n-3\right) -\mu^2 + 5\mu &= \sigma^2 \end{align}

The means and variances can easily be generated programmatically. E.g. to check the means and variances from the table above use

n=2:10
m=c(0,2**(n+1))
v=2**(n+2)*(m[n]-n-3) + 5*m[n] - m[n]^2


Finally, I'm not sure what you wanted when you wrote

when a tails hits and breaks the streak of heads the count would start again from the next flip.

If you meant what is the probability distribution for the next time at which the first run of $n$ or more heads ends, then the crucial point is contained in this comment by @Glen_b, which is that the process begins again after one tail (c.f. the initial problem where you could get a run of $n$ or more heads immediately).

This means that, for example, the mean time to the first event is $\mu-1$, but the mean time between events is always $\mu+1$ (the variance is the same). It's also possible to use a transition matrix to investigate the long-term probabilities of being in a state after the system has "settled down". To get the appropriate transition matrix, set $X_{k,k,}=0$ and $X_{1,k}=1$ so that the system return immediately to state $H_0$ from state $H_*$. Then the scaled first eigenvector of this new matrix gives the stationary probabilities. With $n=4$ these stationary probabilities are

$$\begin{array}{c|cccccc} & \text{probability} \\ \hline H_0 & 0.48484848 \\ H_1 & 0.24242424 \\ H_2 & 0.12121212 \\ H_3 & 0.06060606 \\ H_4 & 0.06060606 \\ H_* & 0.03030303 \end{array}$$ The expected time between states is given by the reciprocal of the probability. So the expected time between visits to $H_* = 1/0.03030303 = 33 = \mu + 1$.

Appendix: Python program used to generate exact probabilities for n = number of consecutive heads over N tosses.

import itertools, pylab

def countinlist(n, N):
count = [0] * N
sub = 'h'*n+'t'
for string in itertools.imap(''.join, itertools.product('ht', repeat=N+1)):
f = string.find(sub)
if (f>=0):
f = f + n -1 # don't count t, and index in count from zero
count[f] = count[f] +1
# uncomment the following line to print all matches
# print "found at", f+1, "in", string
return count, 1/float((2**(N+1)))

n = 4
N = 24
counts, probperevent = countinlist(n,N)
probs = [count*probperevent for count in counts]

for i in range(N):
print '{0:2d} {1:.10f}'.format(i+1,probs[i])
pylab.title('Probabilities of getting {0} consecutive heads in {1} tosses'.format(n, N))
pylab.xlabel('toss')
pylab.ylabel('probability')
pylab.plot(range(1,(N+1)), probs, 'o')
pylab.show()


I am not sure the beta will be likely to be particularly suitable as a way of dealing with this problem -- "The number of plays until ..." is clearly a count. It's an integer, and there is no upper limit on values where you get positive probability.

By contrast the beta distribution is continuous, and on a bounded interval, so it would seem to be an unusual choice. If you moment match a scaled beta, the cumulative distribution functions might perhaps approximate not so badly in the central body of the distribution. However, some other choice is likely to substantially better further into either tail.

If you have either an expression for the probabilities or simulations from the distribution (which you presumably need in order to find an approximating beta), why would you not use those directly?

If your interest is in finding expressions for probabilities or the probability distribution of the number of tosses required, probably the simplest idea is to work with probability generating functions. These are useful for deriving functions from recursive relationships among probabilities, which functions (the pgf) in turn allow us to extract whatever probabilities we require.

Here's a post with a good answer taking the algebraic approach, which explains both the difficulties and makes good use of pgfs and recurrence relations. It has specific expressions for mean and variance in the "two successes in a row" case:

https://math.stackexchange.com/questions/73758/probability-of-n-successes-in-a-row-at-the-k-th-bernoulli-trial-geometric

The four successes case will be substantially more difficult of course. On the other hand, $p = \frac{1}{2}$ does simplify things somewhat.

--

If you just want numerical answers, simulation is relatively straightforward. The probability estimates could be used directly, or alternatively it would be reasonable to smooth the simulated probabilities.

If you must use an approximating distribution, you can probably choose something that does pretty well.

It's possible a mixture of negative binomials (the 'number of trials' version rather than 'the number of successes' version) might be reasonable. Two or three components should be expected to give a good approximation in all but the extreme tail.

If you want a single continuous distribution for an approximation, there may be better choices than the beta distribution; it would be something to investigate.

Okay, I've since done a little algebra, some playing with recurrence relations, some simulation and even a little thinking.

To a very good approximation, I think you can get away with simply specifying the first four nonzero probabilities (which is easy), computing the next few handfuls of values via recurrence (also easy) and then using a geometric tail once the recurrence relation has smoothed out the initially less smooth progression of probabilities.

It looks like you can use the geometric tail to very high accuracy past k=20, though if you're only worried about say 4 figure accuracy you could bring it in earlier.

This should let you compute the pdf and cdf to good accuracy.

I'm a little concerned - my calculations give that the mean number of tosses is 30.0, and the standard deviation is 27.1; if I understood what you mean by "x" and "u", you got 40 and 28 in your tossing. The 28 looks okay but the 40 seems quite a way off what I got... which makes me worry I did something wrong.

====

NOTE: Given the complexities between first time and subsequent times that we encountered, I just want to be absolutely certain now that we are counting the same thing.

Here's a short sequence, with the ends of the '4 or more H' sequences marked (pointing to the gap between flips immediately after the last H)

       \/                     \/
TTHHHHHHTTHTTTTTHHTTHTTHHTHHHHHT...
/\                     /\


Between those two marks I count 23 flips; that is as soon as the previous sequence of (6 in this case) H's ends, we start counting at the immediately following T and then we count right to the end of the sequence of 5 H's (in this case) that ends the next sequence, giving a count of 23 in this case.

Is that how you count them?

Given the above is correct, this is what the probability function of the number of tosses after one run of at least 4 heads is complete until the next run of at least 4 heads is complete looks like:

At first glance it looks like it's flat for the first few values, then has a geometric tail, but that impression is not quite accurate - it takes a while to settle down to an effectively geometric tail.

I am working on coming up with a suitable approximation you can use to answer whatever questions you'd have about probabilities associated with this process to good accuracy that is at the same time as simple as possible. I have a pretty good approximation that should work (that I have already checked against a simulation of a billion coin tosses) but there's some (small but consistent) bias in the probabilities the approximation gives in part of the range and I'd like to see if I can get an extra digit of accuracy out of it.

It may be that the best way to do it is simply to give you a table of the probability function and cdf out to a point beyond which a geometric distribution can be used.

However, it would help if you can give some idea of the range of things you need to use the approximation for.

I hope to follow through on the pgf approach, but it's possible someone else will be more proficient with them than me and can do not just the 4-case but other cases.

• To perhaps clarify things further. A distribution that adjusts or appromixate simulation that takes into account the fluxuation of 4 more successful heads would be ideal. For example, if the populatoin mean is 150 flips for 4 consecutive heads. If 4 or more heads came at the 8th flip. It's unlikely that another 4 or more heads wouldn't come in another 20 or so flips (I'm just guessing) and perhaps be closer to the mean. Something that would get me the probability for when its probable 4 consecutive heads will occur within a certain range of tosses would be AMAZING. – Dan Apr 18 '13 at 5:04
• When you've just had 4 head, if you get a 5th head, does that most recent set of 4 count as another set of 4, or does the count reset, so you start again from the first head (as soon as you see one)? – Glen_b Apr 18 '13 at 5:11
• (I assumed so far that if you generate many sequences of four then there's no overlap - once you get 4, the count of S's resets to 0.) – Glen_b Apr 18 '13 at 9:50
• Its 4 heads or more so as soon as you'd get a tail after 4 heads the streak would stop. Then the count would restart until you see 4 heads or more again consecutively. – Dan Apr 19 '13 at 18:18
• 4 heads or more - I see that's actually what it says in the question, I just hadn't understood it quite right. So 9 heads wouldn't count as two lots of 4 heads. That totally changes the calculations I was doing. The recurrence relation I was using is wrong. The basic concept - that it should have a geometric tail - that will still hold though. – Glen_b Apr 20 '13 at 1:16

You want the geometric distribution. From Wikipedia:

The probability distribution of the number $X$ of Bernoulli trials needed to get one success, supported on the set { 1, 2, 3, ...}.

Let heads H be a failure and tails T be a success. The random variable $X$ is then the number of coin flips needed to see the first tails. For example, $X=4$ will be the sequence HHHT. Here is the probability distribution for $X$:

$$P(X=x)=(1-p)^{x-1}p$$

However, we only want the number of heads. Let's instead define $Y=X-1$ as the number of heads. Here is it's distribution:

$$P(Y+1=x) = (1-p)^{x-1}p \\ P(Y=x-1) = (1-p)^{x-1}p \\ P(Y=y) = (1-p)^yp$$

for $y=0,1,2,3...$. We assume a fair coin, making $p=0.5$. Therefore: \begin{align} P(Y=y) &= (0.5)^y(0.5) \\ &= 0.5^{y+1} \end{align}

This all assumes the number of flips $n$ is sufficiently large (like 10,000). For smaller (finite) $n$, we would need to add a normalization factor to the expression. Put simply, we need to ensure that the total sum is equal to 1. We can do this by dividing by the sum of all probabilities, defined here as $\alpha$:

$$\alpha = \sum_{i=0}^{n-1}{P(Y=i)}$$

This means the corrected form of $Y$, denoted $Z$, will be:

\begin{align} P(Z=z) &= \frac{1}{\alpha}(1-p)^zp \\ &= \frac{1}{\sum_{i=0}^{n-1}{(1-p)^ip}}(1-p)^zp \end{align}

Again, with $p=0.5$, we can reduce this even further, using a geomoetric series summation:

\begin{align} P(Z=z) &= \frac{1}{\sum_{i=0}^{n-1}{0.5^{i+1}}} 0.5^{z+1} \\ &= \frac{1}{1-0.5^n} 0.5^{z+1} \\ &= \frac{0.5^{z+1}}{1-0.5^n} \end{align}

And we can see that, as $n \rightarrow \infty$, our modified version $Z$ approaches $Y$ from earlier.

• I think there are some details of the question you may have missed. Unless I've badly misunderstood the question, it's not simply geometric. – Glen_b Mar 20 '14 at 23:01
• I updated it to handle finite n. And yes, I see now that he wanted to move a window rather than exact counts. Mine only works for chains, not the time between them. – clintonmonk Mar 20 '14 at 23:24
• A good first step is to have a look at the graph in @Glen_b's post, and see if you can replicate that. I've also added the Python program I wrote to check the exact probabilities. If you are able to run this, uncomment the line that prints matches, decreasing N to somewhere between 5 and 7, and you'll get a good feel for the events that are required (note pylab is only required for plotting). – TooTone Mar 21 '14 at 0:56
• Unfortunately I'm not around a PC that I can test that with currently. I started to use a Markov process to show the stationary sol'n was geometric (and E[time to return] = 1/$\pi_i$), but I didn't have time to fully flesh it out. – clintonmonk Mar 21 '14 at 1:35
• Yes, if by stationary solution, you're talking about the ratio of consecutive probabilities in the tail converging to a constant, then the stationary solution is indeed geometric, as both previous answers have said. – Glen_b Mar 21 '14 at 8:47