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?
 A: 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.
A: 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$
A: 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...$.
A: 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.
