1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63

"""Bug discovery curves: fitting a reliability-growth model to weekly bug counts.

The Goel-Okumoto model says the expected number of bugs found by time t is
m(t) = a (1 - exp(-b t)): a is the total number of bugs that testing of this kind will find,
b how fast they are found (fact M7). A project (illustrative, a = 80, b = 0.12 per week) finds
bugs week by week as Poisson counts. At week 8 and at week 16 the model is fitted to the counts
so far by maximum likelihood, giving an estimate of a and so of the bugs still to find.
Then 300 more simulated projects with the same true curve show how much such estimates spread.
"""
import math
import random

rng = random.Random(21)
A_TRUE, B_TRUE, WEEKS = 80, 0.12, 16

def poisson(lam):
    k, p, l = 0, 1.0, math.exp(-lam)
    while True:
        p *= rng.random()
        if p < l:
            return k
        k += 1

m = lambda a, b, t: a * (1 - math.exp(-b * t))
weekly = [poisson(m(A_TRUE, B_TRUE, w) - m(A_TRUE, B_TRUE, w - 1)) for w in range(1, WEEKS + 1)]

def fit(counts):
    """Maximum likelihood for grouped Poisson counts; a has a closed form for each b."""
    T, n = len(counts), sum(counts)
    best = None
    for i in range(1, 1000):
        b = i / 2000
        a = n / (1 - math.exp(-b * T))
        ll = 0.0
        for w, c in enumerate(counts, start=1):
            lam = m(a, b, w) - m(a, b, w - 1)
            ll += c * math.log(lam) - lam
        if best is None or ll > best[0]:
            best = (ll, a, b)
    return best[1], best[2]

print("week " + " ".join(f"{w:3d}" for w in range(1, WEEKS + 1)))
print("bugs " + " ".join(f"{c:3d}" for c in weekly))
for upto in (8, 16):
    a, b = fit(weekly[:upto])
    found = sum(weekly[:upto])
    print(f"fit at week {upto:2d}: a = {a:5.1f}, b = {b:.3f}; found {found}, estimated still to find {a - found:5.1f}")
a16, b16 = fit(weekly)
fitted = [m(a16, b16, w) - m(a16, b16, w - 1) for w in range(1, WEEKS + 1)]
print("fit  " + " ".join(f"{f:3.0f}" for f in fitted))

def spread(upto, runs=300):
    est = []
    for _ in range(runs):
        counts = [poisson(m(A_TRUE, B_TRUE, w) - m(A_TRUE, B_TRUE, w - 1)) for w in range(1, upto + 1)]
        est.append(fit(counts)[0])
    est.sort()
    return est[runs // 10], est[runs // 2], est[9 * runs // 10]
s8, s16 = spread(8), spread(16)
print(f"300 projects, estimated total at week  8: median {s8[1]:5.1f}, 10-90% range {s8[0]:5.1f} to {s8[2]:5.1f}")
print(f"300 projects, estimated total at week 16: median {s16[1]:5.1f}, 10-90% range {s16[0]:5.1f} to {s16[2]:5.1f}")
assert (s16[2] - s16[0]) < (s8[2] - s8[0]) and abs(s16[1] - A_TRUE) < 10
print(f"PASS: estimates centre near the true {A_TRUE} but spread widely; the range narrows from week 8 to week 16")