'distribution of success runs of length k or more
Parameters
x : float
count of runs of length n
k : int
length of runs
n : int
total number of observations or trials
p : float
probability of success in each Bernoulli trial
Returns
pdf : float
probability that x runs of length of k are observed
Notes
not yet vectorized
References
Muselli 1996, theorem 3'

def pdf(self, x, k, n, p):

    q = (1 - p)
    m = np.arange(x, (((n + 1) // (k + 1)) + 1))[:, None]
    terms = ((((((-1) ** (m - x)) * comb(m, x)) * (p ** (m * k))) * (q ** (m - 1))) * (comb((n - (m * k)), (m - 1)) + (q * comb((n - (m * k)), m))))
    return terms.sum(0)
