Expected number of runs of heads of size $k$? I just completed this question, and I think I have a solution but I am not entirely sure. Here's the question, with an attempt later:

A coin having a probability $p$ of landing on heads is flipped $n$ times. Compute the expected number of runs of heads of size $1$, of size $2$, and of size $k$, $1 \leq k \leq n$.

Here's my attempt at a solution which I am not entirely sure is correct:

Let $x_i$ be $1$ if the $i$th coin flip is a tail that completes a run of size $k$ and $0$ otherwise. Doing so, we find that the following is true:
          $$E[x_i] = p^k(1-p) \text{ for } i > k \text{, 0 otherwise}$$
          Therefore, the expected number of runs of size $k$ is the following, where $y$ is $1$ if the last coin completes a run of $k$ heads and $0$ otherwise:
          $$E[x_1 + x_2 + \ldots + x_n + y] = \sum_{i = k+1}^{n}p^k(1-p) + p^k$$
          Simplifying the above sum, we find that the following is the expected number of runs of size $k$:
          $$E[x_1 + x_2 + \ldots + x_n + y] = (n-k)p^k(1-p) + p^k$$

Does this work? If not, please explain?
 A: I hate to do this, because I love generating functions and I admire Marko's agility in manipulating them in interesting ways to solve difficult problems and his readiness to share it with us, but in this case using generating functions is what we'd call "shooting at sparrows with cannons" in German.
There are $n-k-1$ opportunities for a run of $k$ heads bounded by tails on either side, and $2$ opportunities for a run of $k$ heads at either end, bounded by tails on one side. By linearity of expectation, the expected number of runs is the sum of the probabilities of these runs, which is
$$
(n-k-1)p^k(1-p)^2+2p^k(1-p)=(n-k)p^k(1-p)^2+p^k(1-p^2)\;.
$$ 
A: Here  is  the generating  function  for  runs  that are  not  maximal,
i.e.  substrings  of a  run  count towards  shorter  runs.   We use  a
trivariate generating function with $z$  and $w$ marking runs of heads
and tails resp. and $v$ marking runs of heads of size $k.$
We obtain
$$G(z,w,v)
= (1+z+\cdots+vz^k+v^2z^{k+1}+\cdots)
\\ \times \left(\sum_{q\ge 0} 
(w+w^2+\cdots)^q (z+\cdots+vz^k+v^2z^{k+1}+\cdots)^q\right)
(1+w+w^2+\cdots).$$
This is
$$\left(\frac{1-z^k}{1-z}+vz^k\frac{1}{1-vz}\right)
\left(\sum_{q\ge 0} \frac{w^q}{(1-w)^q}
\left(z\frac{1-z^{k-1}}{1-z}+vz^k\frac{1}{1-vz}\right)^q\right)
\frac{1}{1-w}$$
which simplifies to
$$\frac{1 - z^k - vz + vz^{k+1} + vz^k - vz^{k+1}}{(1-z)(1-vz)}
\\ \times \frac{1}
{1-w(z - z^k - vz^2 + vz^{k+1} + vz^{k} - vz^{k+1})
/(1-z)/(1-vz)/(1-w)}
\\ \times \frac{1}{1-w}
\\ = \frac{1 - z^k - vz + vz^{k}}
{(1-z)(1-vz)(1-w)-w(z - z^k - vz^2 + vz^k)}.$$
As a sanity check when we set $w=z$ and $v=1$ this turns into
$$\frac{1-z^k-z+z^{k}}
{(1-z)^3-z(z - z^k - z^2 + z^k)}
\\ = \frac{1-z}{(1-z)^3 - z^2(1-z)}
= \frac{1}{(1-z)^2 - z^2}
= \frac{1}{1-2z}.$$
This  says we  have $2^n$  strings of  length $n$  and the  check goes
through.

Continuing with the usual procedure we compute
$$\left.\frac{\partial}{\partial v} G(z, w, v)\right|_{v=1}$$
to get for the derivative
$$\frac{z^k-z}{(1-z)(1-vz)(1-w)-w(z - z^k - vz^2 + vz^k)}
\\ - \frac{1 - z^k - vz + vz^{k}}
{((1-z)(1-vz)(1-w)-w(z - z^k - vz^2 + vz^k))^2}
\times (-z(1-z)(1-w) + w (z^2-z^k)).$$
Evaluating at $v=1$ we have
$$(1-z)(1-vz)(1-w) - w(z - z^k - vz^2 + vz^k)
\\ = (1-z)^2 (1-w) - w(z-z^2)
= (1-z)(1-w-z)$$
and obtain
$$\frac{z^k-z}{(1-z)(1-w-z)}
\\ - \frac{1}{(1-z)(1-w-z)^2}
(-z(1-z)(1-w) + w (z^2-z^k))
\\ = \frac{z^k(1-z)}{(1-z)(1-w-z)^2}
= \frac{z^k}{(1-w-z)^2}.$$
Introducing the probabilities we find
$${\large \left.\frac{\partial}{\partial v} G(z, w, v)
\right|_{v=1, z=pz, w=(1-p)z} 
= \frac{p^k z^k}{(1-(1-p)z-pz)^2}
= \frac{p^k z^k}{(1-z)^2}}.$$
Extracting coefficients we finally have the answer
$$[z^n] \frac{p^k z^k}{(1-z)^2}
= p^k [z^{n-k}] \frac{1}{(1-z)^2} 
= p^k (n-k+1).$$
The  formulae  that  were   presented  here  were  computed  by  total
enumeration as well as symbolically.
This was the Maple code.

V :=
proc(n, k)
    option remember;
    local gf, ind, d, cnt, runs, run, runq,
    cur, len, pos;

    gf := 0;

    for ind from 2^n to 2^(n+1)-1 do
        d := convert(ind, base, 2);

        cur := d[1]; len := 1; runs := [];

        cnt := [0, 0]; cnt[cur+1] := 1;

        for pos from 2 to n+1 do
            if pos = n+1 or d[pos] <> cur then
                runs := [op(runs), [cur, len]];

                cur := d[pos]; len := 1;
            else
                len := len + 1;
            fi;

            if pos < n+1 then
                cnt[cur+1] := cnt[cur+1] + 1;
            fi;
        od;

        runq := 0;

        for run in runs do
            if run[1] = 1 and run[2] >= k then
                runq := runq + 1+run[2]-k;
            fi;
        od;

        gf := gf  + z^cnt[2]*w^cnt[1]*v^runq;
    od;

    gf;
end;

S :=
proc(k)
    local gf;

    gf := (1-z^k-v*z+v*z^k)
    /((1-z)*(1-v*z)*(1-w)
      -w*(z-z^k-v*z^2+v*z^k));

    factor(subs([w=z, v=1], gf));
end;

G :=
proc(n, k)
    option remember;
    local gf, cf;

    gf := (1-z^k-v*z+v*z^k)
    /((1-z)*(1-v*z)*(1-w)
      -w*(z-z^k-v*z^2+v*z^k));

    cf :=
    add(coeftayl(coeftayl(gf, z=0, q), w=0, n-q)
        *z^q*w^(n-q), q=0..n);

    expand(cf);
end;


X := (n, k) -> subs(v=1, diff(subs([z=p, w=1-p], V(n,k)), v));
XX := (n, k) -> p^k*(n-k+1);

