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

Question

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?

Answer 1

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...$.

Answer 2

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

def heads_in_a_row(flips, p, want):
    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]

print(heads_in_a_row(flips=25, p=1.0 / 3.0, want=8))

yields 0.00187928367413 as desired.

Answer 3

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.

Answer 4

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.

Answer 5

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$

Answer 6

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