# 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? – cardinal Jan 28 '12 at 15:05

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? – Michelle Jan 27 '12 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. – AcidNynex Jan 28 '12 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. :) – cardinal Jan 28 '12 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$