# Probability of a run of k successes in a sequence of n Bernoulli trials

I'm trying to find the probability of getting 8 trials in a row correct in a block of 25 trials, you have 8 total blocks (of 25 trials) to get 8 trials correct in a row. The probability of getting any trial correct based on guessing is 1/3, after getting 8 in a row correct the blocks will end (so getting more than 8 in a row correct is technically not possible). How would I go about finding the probability of this occurring? I've been thinking along the lines of using (1/3)^8 as the probability of getting 8 in a row correct, there are 17 possible chances to get 8 in a row in a block of 25 trials, if I multiply 17 possibilities * 8 blocks I get 136, would 1-(1-(1/3)^8)^136 give me the likelihood of getting 8 in a row correct in this situation or am I missing something fundamental here?

• I believe the problem with the argument given is that the events considered are not independent. For example, consider a single block. If I tell you that (a) there is no run of eight that starts at position 6, (b) there is a run starting at position 7 and (c) there is no run starting at position 8, what does that tell you about the probability of a run starting at positions, say, 9 through 15? Commented Jan 28, 2012 at 15:05
• I believe that the fastest solution in the form of the analytical formula can be found here: math.stackexchange.com/a/59749/51275 Commented Apr 11, 2023 at 12:55

By keeping track of things you can get an exact formula.

Let $p=1/3$ be the probability of success and $k=8$ be the number of successes in a row you want to count. These are fixed for the problem. Variable values are $m$, the number of trials left in the block; and $j$, the number of successive successes already observed. Let the chance of eventually achieving $k$ successes in a row before $m$ trials are exhausted be written $f_{p,k}(j,m)$. We seek $f_{1/3,8}(0,25)$.

Suppose we have just seen our $j^\text{th}$ success in a row with $m\gt0$ trials to go. The next trial is either a success, with probability $p$--in which case $j$ is increased to $j+1$--; or else it is a failure, with probability $1-p$--in which case $j$ is reset to $0$. In either case, $m$ decreases by $1$. Whence

$$f_{p,k}(j,m) = p f_{p,k}(j+1,m-1) + (1-p)f_{p,k}(0,m-1).$$

As starting conditions we have the obvious results $f_{p,k}(k,m)=1$ for $m \ge 0$ (i.e., we have already seen $k$ in a row) and $f_{p,k}(j,m)=0$ for $k-j \gt m$ (i.e., there aren't enough trials left to get $k$ in a row). It is now fast and straightforward (using dynamic programming or, because this problem's parameters are so small, recursion) to compute

$$f_{p,8}(0,25) = 18p^8 - 17p^9 - 45p^{16} + 81p^{17}-36p^{18}.$$

When $p=1/3$ this yields $80897 / 43046721 \approx 0.0018793$.

Relatively fast R code to simulate this is

hits8 <- function() {
x <- rbinom(26, 1, 1/3)                # 25 Binomial trials
x[1] <- 0                              # ... and a 0 to get started with diff
if(sum(x) >= 8) {                      # Are there at least 8 successes?
max(diff(cumsum(x), lag=8)) >= 8   # Are there 8 successes in a row anywhere?
} else {
FALSE                              # Not enough successes for 8 in a row
}
}
set.seed(17)
mean(replicate(10^5, hits8()))


After 3 seconds of calculation, the output is $0.00213$. Although this looks high, it's only 1.7 standard errors off. I ran another $10^6$ iterations, yielding $0.001867$: only $0.3$ standard errors less than expected. (As a double-check, because an earlier version of this code had a subtle bug, I also ran 400,000 iterations in Mathematica, obtaining an estimate of $0.0018475$.)

This result is less than one-tenth the estimate of $1-(1-(1/3)^8)^{136} \approx 0.0205$ in the question. But perhaps I have not fully understood it: another interpretation of "you have 8 total blocks ... to get 8 trials correct in a row" is that the answer being sought equals $1 - (1 - f_{1/3,8}(0,25))^8) = 0.0149358...$.

While @whuber's excellent dynamic programming solution is well worth a read, its runtime is $\mathcal O(k^2m)$ with respect to total number of trials $m$ and the desired trial length $k$ whereas the matrix exponentiation method is $\mathcal O(k^3\log(m))$. If $m$ is much larger than $k$, the following method is faster.

Both solutions consider the problem as a Markov chain with states representing the number of correct trials at the end of the string so far, and a state for achieving the desired correct trials in a row. The transition matrix is such that seeing a failure with probability $p$ sends you back to state 0, and otherwise with probability $1-p$ advances you to the next state (the final state is an absorbing state). By raising this matrix to the $n$th power, the value in the first row, and last column is the probability of seeing $k=8$ heads in a row. In Python:

import numpy as np

a = np.zeros((want + 1, want + 1))
for i in range(want):
a[i, 0] = 1 - p
a[i, i + 1] = p
a[want, want] = 1.0
return np.linalg.matrix_power(a, flips)[0, want]



yields 0.00187928367413 as desired.

According to this answer, I will explain the Markov-Chain approach by @Neil G a bit more and provide a general solution to such problems in R. Let's denote the desired number of correct trials in a row by $$k$$, the number of trials as $$n$$ and a correct trial by $$W$$ (win) and an incorrect trial by $$F$$ (fail). In the process of keeping track of the trials, you want to know whether you already had a streak of 8 correct trials and the number of correct trials at the end of your current sequence. There are 9 states ($$k+1$$):

$$A$$: We have not had $$8$$ correct trials in a row yet, and the last trial was $$F$$.

$$B$$: We have not had $$8$$ correct trials in a row yet, and the last two trials were $$FW$$.

$$C$$: We have not had $$8$$ correct trials in a row yet, and the last three trials were $$FWW$$.

$$\ldots$$

$$H$$: We have not had $$8$$ correct trials in a row yet, and the last eight trials were $$FWWWWWWW$$.

$$I$$: We've had $$8$$ correct trials in a row!

The probability of moving to state $$B$$ from state $$A$$ is $$p=1/3$$ and with probability $$1-p=2/3$$ we stay in state $$A$$. From state $$B$$, the probability of moving to state $$C$$ is $$1/3$$ and with probability $$2/3$$ we move back to $$A$$. And so on. If we are in state $$I$$, we stay there.

From this, we can construct a $$9\times9$$ transition matrix $$M$$ (as each column of $$M$$ sums to $$1$$ and all entries are positive, $$M$$ is called a left stochastic matrix):

$$M= \begin{pmatrix} 2/3 & 2/3 & 2/3 & 2/3 & 2/3 & 2/3 & 2/3 & 2/3 & 0 \\ 1/3 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1/3 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1/3 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1/3 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1/3 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1/3 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1/3 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1/3 & 1 \end{pmatrix}$$

Each column and row corresponds to one state. After $$n$$ trials, the entries of $$M^{n}$$ give the probability of getting from state $$j$$ (column) to state $$i$$ (row) in $$n$$ trials. The rightmost column corresponds to the state $$I$$ and the only entry is $$1$$ in the right lower corner. This means that once we are in state $$I$$, the probability to stay in $$I$$ is $$1$$. We are interested in the probability of getting to state $$I$$ from state $$A$$ in $$n=25$$ steps which corresponds to the lower left entry of $$M^{25}$$ (i.e. $$M^{25}_{91}$$). All we have to do now is calculating $$M^{25}$$. We can do that in R with the matrix power function from the expm package:

    library(expm)

k <- 8   # desired number of correct trials in a row
p <- 1/3 # probability of getting a correct trial
n <- 25  # Total number of trials

# Set up the transition matrix M

M <- matrix(0, k+1, k+1)

M[ 1, 1:k ] <- (1-p)

M[ k+1, k+1 ] <- 1

for( i in 2:(k+1) ) {

M[i, i-1] <- p

}

# Name the columns and rows according to the states (A-I)

colnames(M) <- rownames(M) <- LETTERS[ 1:(k+1) ]

round(M,2)

A    B    C    D    E    F    G    H I
A 0.67 0.67 0.67 0.67 0.67 0.67 0.67 0.67 0
B 0.33 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0
C 0.00 0.33 0.00 0.00 0.00 0.00 0.00 0.00 0
D 0.00 0.00 0.33 0.00 0.00 0.00 0.00 0.00 0
E 0.00 0.00 0.00 0.33 0.00 0.00 0.00 0.00 0
F 0.00 0.00 0.00 0.00 0.33 0.00 0.00 0.00 0
G 0.00 0.00 0.00 0.00 0.00 0.33 0.00 0.00 0
H 0.00 0.00 0.00 0.00 0.00 0.00 0.33 0.00 0
I 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.33 1

# Calculate M^25

Mn <- M%^%n
Mn[ (k+1), 1 ]
[1] 0.001879284


The probability of getting from state $$A$$ to state $$I$$ in 25 steps is $$0.001879284$$, as established by the other answers.

Here is some R code that I wrote to simulate this:

tmpfun <- function() {
x <- rbinom(25, 1, 1/3)
rx <- rle(x)
any( rx$lengths[ rx$values==1 ] >= 8 )
}

tmpfun2 <- function() {
any( replicate(8, tmpfun()) )
}

mean(replicate(100000, tmpfun2()))


I am getting values a little smaller than your formula, so one of us may have made a mistake somewhere.

• Does your function include trials where it is impossible to get 8 in a row right, e.g. where the "run" started on trial 20? Commented Jan 27, 2012 at 21:03
• Most likely me, my R simulation is giving me smaller values as well. I'm just curious if there is an algebraic solution to solve this as a simple probability issue in case someone disputes a simulation. Commented Jan 28, 2012 at 0:48
• I think this answer would be improved by providing the output you obtained so that it can be compared. Of course, including something like a histogram in addition would be even better! The code looks right to me at first glance. Cheers. :) Commented Jan 28, 2012 at 15:06

Here is a Mathematica simulation for the Markov chain approach, note that Mathematica indexes by $1$ not $0$:

M = Table[e[i, j] /. {
e[9, 1] :> 0,
e[9, 9] :> 1,
e[_, 1] :> (1 - p),
e[_, _] /; j == i + 1 :> p,
e[_, _] :> 0
}, {i, 1, 9}, {j, 1, 9}];

x = MatrixPower[M, 25][[1, 9]] // Expand


This would yield the analytical answer: $$18 p^8 - 17 p^9 - 45 p^{16} + 81 p^{17} - 36 p^{18}$$

Evaluating at $p=\frac{1.0}{3.0}$

x /. p -> 1/3 // N


Will return $0.00187928$

This can also be evaluated directly using builtin Probability and DiscreteMarkovProcess Mathematica functions:

Probability[k[25] == 9, Distributed[k, DiscreteMarkovProcess[1, M /. p -> 1/3]]] // N


Which will get us the same answer: $0.00187928$

As stated in my comment, there is an analytical formula for this problem in https://math.stackexchange.com/a/59749/51275.

The accepted answer, translated into Python looks like this:

import numpy as np
from scipy.special import binom

def theoretical_distr2(N, p, run_length)->float:
# Mathematica code: Sum[(-1)^(j + 1) * (p + (n - j*m + 1)/j*(1 - p))*Binomial[n - j m, j - 1] p^(j*m)*(1 - p)^(j - 1), {j, 1, Floor[n/m]}]
# where n = N, m = run_length, p = p

ans = 0
for j in range(1, N // run_length + 1):
ans += (-1) ** (j + 1) * (p + (N - j * run_length + 1) / j * (1 - p)) * binom(N - j * run_length, j - 1) * p ** (j * run_length) * (1 - p) ** (j - 1)

return ans


I don't deserve a credit. I just found the other solution and translated it into a code. There's an obvious issue with the code: for large N binom explodes to infinity, yielding a NaN answer. The @Neil's code is more stable.

For large N we can use the approximate formula (source: https://www.maa.org/sites/default/files/pdf/upload%5flibrary/22/Polya/07468342.di020742.02p0021g.pdf)

def theoretical_distr3(N, p, run_length)->float:
# Mathematica code: -E^-p^(0.5 + m - Log[n - n p]/Log[1/p]) p^(0.5 + m - Log[n - n p]/Log[1/p]) Log[p]
# -Exp[-p^(0.5 + i - Log[n (1 - p)]/Log[1/p])] p^(0.5 + i - Log[n (1 - p)]/Log[1/p]) Log[p]
ans = 1
for i in range(1, run_length):
ans -= -np.exp(-p ** (0.5 + i - np.log(N * (1 - p)) / np.log(1 / p))) * p ** (0.5 + i - np.log(N * (1-p)) / np.log(1 / p)) * np.log(p)
return ans