What's the probability we roll the string "66" if we roll a dice $n$ times? I decided to try and solve this with a recurrence. Let $P_E(k)$ represent the probability we roll a double 6 somewhere in those $k$ rolls starting from a $1,...,5$ and $P_6(k)$ represent the probability we roll a double 6 starting from 6. Then $P_E(k) = 1/36 + P_E(k-2)$. This recurrence gives me an answer that doesn't agree with the smaller cases of $k=2,3,$ etc. How should this problem be approached?
Generally how can we answer the probability that we roll any string of faces somewhere in $k$ rolls?
 A: I try to define $f(k)$ as the probability that double 6 is achieved at k rolls.
$g(k)$ as the probability the double 6 is not yet achieved and last roll is a 6.
Then,
$$ f(1)=0, g(1)=1/6\\
f(k+1)=f(k)+g(k)/6\\
g(k+1)=[1-f(k)-g(k)]/6\\
where \quad k\geqslant1
$$
A: Let the probability of having at least one $66$ by throw $n$ be $P(H_{n})$ and the probability of having no $66$ by throw $n$ be $P(M_{n})=1-P(H_n)$, where $H$ and $M$ stand for "hit" and "miss".
Label the probability of dice $n$ being $6$ be $P(D_n)$, then the recurrence relation is
$$
P(H_{n+1})=P(H_n)+P(D_{n+1})P(M_n)P(D_n|M_n).
$$
with $P(M_1)=1$, $P(H_1)=0$, $P(M_2)=35/36$, $P(H_2)=1/36$.
Now $P(D_n)=\frac16$ and is independent of $n$.
$P(D_n|M_n)$ is the probability that dice $n$, the final dice of the $n$-dice combination which did not contain two sixes, was a 6. This is slightly tricky.
If the number of non-double-six sequences ending with $6$ is called $P_n$ and the number of non-double sequences not ending with $6$ is $Q_n$ then $P_n=Q_{n-1}$, since we can only add a six on to the end of a sequence which does not end in a six.
If the total number of sequences is $R_n$, then $$R_n=P_n+Q_n.$$
Now $$P(D_n|M_n)={P_n\over R_n}$$
and $R_n=P_n+Q_n$. $Q_n$ has the recurrence relation $Q_n=5Q_{n-1}+5Q_{n-2}$ because we can either add a non-six to an $n-1$ sequence not ending in six, or we can add a six followed by a non-six to a sequence $n-2$ in length.
Graphically
$$
\begin{align}
Q_n=
&|...\bar⚅|&\times&|⚀⚁⚂⚃⚄|&+
&|...\bar⚅|&\times&|⚅|&\times&|⚀⚁⚂⚃⚄|\\
&Q_{n-1}&\times &5&+
&Q_{n-2}&\times &1&\times &5\\
\end{align}
$$
The first few values are $Q_1=5$ and $Q_2=30$ by counting, so $Q_3=175$ and $R_3=P_3+Q_3=Q_2+Q_3=205$ so $P(D_3|M_3)={30\over205}$.
The recurrence relation is then
$$
\begin{align}
P(H_{n+1})
&=P(H_n)+\frac16P(M_n)\times P(D_n|M_n)\\
\end{align}
$$

Here is a demonstration in computer code:
package main

import (
    "fmt"
    "math"
    "math/rand"
)

func main() {
    max := 20
    prob := makeprobs(max)
    ntries := 1000000
    for n := 2; n <= max; n++ {
        ok := 0
        for try := 0; try < ntries; try++ {
            hit := 0
            roll := make([]int, n)
            for i := 0; i < n; i++ {
                roll[i] = rand.Int()%6 + 1
                if i > 0 && roll[i] == 6 && roll[i-1] == 6 {
                    hit++
                }
            }
            if hit > 0 {
                ok++
            }
        }
        p := prob[n]
        pc := float64(ok) / float64(ntries)
        fmt.Printf("%d %g %g %.4f %d %d\n", n, pc, p, pc/p, ok, ntries)
    }
}

const hitp = float64(1) / float64(36)
const missp = float64(35) / float64(36)

func makeprobs(max int) (prob []float64) {
    if max < 2 {
        panic("Max too low")
    }
    // From 1 to max is used here.
    hit := make([]float64, max+1)
    miss := make([]float64, max+1)
    q := make([]int64, max+1)
    cond := make([]float64, max+1)
    // Ensure these are not used accidentally
    miss[0] = math.Inf(1)
    hit[0] = math.Inf(1)
    miss[1] = 1
    miss[2] = missp
    hit[1] = 0
    hit[2] = hitp
    q[0] = 1
    q[1] = 5
    q[2] = 30
    for i := 3; i <= max; i++ {
        q[i] = 5 * (q[i-1] + q[i-2])
        cond[i] = float64(q[i-1]) / float64(q[i]+q[i-1])
        fmt.Printf("%g\n", cond[i])
    }
    for i := 3; i <= max; i++ {
        // The chance of a new hit is the chance of a final six being
        // on the end of the previous miss, 1/6, times the chance of
        // another six.
        g := miss[i-1] / 6.0 * cond[i]
        hit[i] = g + hit[i-1]
        miss[i] = 1 - hit[i]
    }
    return hit
}

Run it online here.
