import numpy as np, pandas as pd, matplotlib.pyplot as plt
from scipy import stats
AMBER="#d97706"; TEAL="#0d9488"; INK="#1a2138"; GRID="#e6e9f2"; PINK="#db2777"
plt.rcParams.update({"figure.facecolor":"white","axes.facecolor":"white","figure.dpi":110,"font.size":11,
   "axes.edgecolor":GRID,"axes.grid":True,"grid.color":GRID,"axes.axisbelow":True,"axes.spines.top":False,
   "axes.spines.right":False,"axes.titlesize":12,"axes.titleweight":"bold","legend.frameon":False})
BASE="https://raw.githubusercontent.com/johnfisher-ai/Statistics-Data-Science-AI-Visual-Book/main/data/"
try:    df = pd.read_csv("../../data/software_reliability_testing.csv")
except FileNotFoundError: df = pd.read_csv(BASE+"software_reliability_testing.csv")
print("loaded:", df.shape)
df.head()

loaded: (1000, 5)

r = int(df.target_successes_r.iloc[0])
fails = df["failures_encountered"]; total = df["total_tests_run"]
# total tests = r successes + failures; per-test pass prob p = r / mean(total)
p = r / total.mean()
print(f"target successes r   : {r}")
print(f"mean total tests run  : {total.mean():.3f}")
print(f"per-test pass prob p = r/mean(total) = {p:.4f}")
print(f"mean failures observed = {fails.mean():.3f}")
print(f"NegBinom mean failures = r(1-p)/p = {r*(1-p)/p:.3f}")

target successes r   : 10
mean total tests run  : 12.640
per-test pass prob p = r/mean(total) = 0.7911
mean failures observed = 2.640
NegBinom mean failures = r(1-p)/p = 2.640

fig,ax=plt.subplots(figsize=(7,3.4))
kmax=int(fails.max())+1; k=np.arange(0,kmax+1)
obs=fails.value_counts(normalize=True).sort_index()
ax.bar(obs.index,obs.values,color=AMBER,alpha=0.7,width=0.9,label="observed failures")
ax.plot(k, stats.nbinom.pmf(k,r,p),"o-",color=TEAL,ms=4,label=f"NegBinom(r={r}, p={p:.2f})")
ax.set_xlabel("failures before the 10th pass"); ax.set_ylabel("probability"); ax.set_title("Failures before 10 passes: observed vs negative binomial"); ax.legend()
plt.tight_layout(); plt.show()

# Formal goodness-of-fit: chi-square on observed vs expected counts.
# Rule of thumb: merge adjacent bins so every EXPECTED count is at least 5, then
# chi-square = sum((obs-exp)**2 / exp), with df = (bins - 1 - number of fitted parameters).
def chisq_gof(counts, pmf, n_params):
    N = len(counts); kmax = int(counts.max())
    obs = np.array([(counts == k).sum() for k in range(kmax + 1)], float)
    exp = np.array([N * pmf(k) for k in range(kmax + 1)], float)
    obs = np.append(obs, 0.0); exp = np.append(exp, N - exp.sum())   # tail bin for k > kmax
    mo, me, co, ce = [], [], 0.0, 0.0
    for o, e in zip(obs, exp):
        co += o; ce += e
        if ce >= 5:
            mo.append(co); me.append(ce); co = ce = 0.0
    if ce > 0 and mo:
        mo[-1] += co; me[-1] += ce
    mo, me = np.array(mo), np.array(me)
    chi2 = ((mo - me) ** 2 / me).sum(); dof = len(mo) - 1 - n_params
    return chi2, dof, stats.chi2.sf(chi2, dof)

chi2, dof, pval = chisq_gof(fails.values, lambda k: stats.nbinom.pmf(k, r, p), 1)
print(f"Chi-square goodness-of-fit: chi2 = {chi2:.2f}, df = {dof}, p = {pval:.3f}  -> "
      + ("cannot reject the fit (NegativeBinomial is a good model)" if pval > 0.05 else "reject: not NegativeBinomial"))

Chi-square goodness-of-fit: chi2 = 1.77, df = 8, p = 0.987  -> cannot reject the fit (NegativeBinomial is a good model)

print(f"mean total tests to 10 passes = r/p = {r/p:.2f}")
for budget in [12,14,16,18]:
    extra = budget - r            # failures allowed within the budget
    print(f"P(suite clears within {budget} tests) = P(failures <= {extra}) = {stats.nbinom.cdf(extra,r,p):.4f}")

mean total tests to 10 passes = r/p = 12.64
P(suite clears within 12 tests) = P(failures <= 2) = 0.5271
P(suite clears within 14 tests) = P(failures <= 4) = 0.8504
P(suite clears within 16 tests) = P(failures <= 6) = 0.9667
P(suite clears within 18 tests) = P(failures <= 8) = 0.9942

rng = np.random.default_rng(53)
# sum of r geometric waits (trials), minus r, = failures
geos = rng.geometric(p, size=(200_000, r)).sum(axis=1)   # total trials
sim_fail = geos - r
print(f"sum-of-geometrics: mean failures {sim_fail.mean():.3f}  (theory {r*(1-p)/p:.3f})")
print(f"scipy nbinom mean {stats.nbinom.mean(r,p):.3f}, var {stats.nbinom.var(r,p):.3f}")
print(f"a tail run: P(> 8 failures) = {stats.nbinom.sf(8,r,p):.4f}")

sum-of-geometrics: mean failures 2.639  (theory 2.640)
scipy nbinom mean 2.640, var 3.337
a tail run: P(> 8 failures) = 0.0058

	test_suite_id	environment	target_successes_r	failures_encountered	total_tests_run
0	TS_8000	Production_Mirror	10	4	14
1	TS_8001	UAT	10	2	12
2	TS_8002	UAT	10	1	11
3	TS_8003	Staging	10	4	14
4	TS_8004	Production_Mirror	10	1	11

⚙️ Setup & data¶