Probability of exactly $R$ runs of $M$ consecutive success I know how to compute the probability that $0$ runs (or one or more) of $M$
consecutive successes will occur in a sequence of $N$ trials, each with  probability $P$ of success.
Is there an explicit formula (mass function) or a method of calculation (for example in Maple) to find the probability that exactly $R$ runs of $M$ consecutive successes will occur in a sequence of $N$ trials, each trial with probability $P$ of success?
I have the same question about the cumulative distribution function?
 A: Let's first count the probability that exactly $R$ runs of $M$ successes happen in $N$ trials, and the last trial fails. Denote by $p$ the probability of success, and by $q$ the probability of failure.
Each set of runs can be "tokenized" into parts of the form $S^kF$ (this is because the last trial fails). If $k < M$ then we count no runs, and otherwise we count $k - M + 1$ runs. We're going to form a bivariate generating function, with $x$ being the total number of trials, and $y$ being the number of runs. Each part contributes the following:
$$P(x,y) = xq(1 + px + \cdots + (px)^{M-1} + (px)^My + (px)^{M+1}y^2 + \cdots).$$
We can get an explicit formula for $P(x,y)$ using geometric series identities:
$$P(x,y) = xq\left(\frac{(px)^{M-1}-1}{px-1} + \frac{(px)^{M-1}}{1-pxy}\right).$$
The required probability is the coefficient of $x^N y^R$ in $1/(1-P(x,y))$.
Now let's do the same when the last trial succeeds. The last piece is then $P(x,y)/xq$, and so we want the coefficient of $x^N y^R$ in $(P(x,y)/xq)/(1-P(x,y))$. In total, the required probability is the coefficient of $x^N y^R$ in $$ \frac{1+P(x,y)/xq}{1-P(x,y)}. $$
