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?

  • 1
    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

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.

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
  • 1
    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$

Your Answer

 
discard

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.