""" 04-05 Binomial & Poisson Calculator Chapters 06-01, 06-02, 12-01 Run from the project folder: python python/analysis.py """ import numpy as np import pandas as pd from scipy import stats import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt defects = pd.read_csv("data/batch_defects.csv") calls = pd.read_csv("data/calls_per_hour.csv") # ── 1. REUSABLE CALCULATORS ───────────────────────────────────────────────── def binom_calc(n, p, x, a=None, b=None): out = { "P(X = x)": stats.binom.pmf(x, n, p), "P(X <= x)": stats.binom.cdf(x, n, p), "P(X >= x)": stats.binom.sf(x - 1, n, p), "P(X < x)": stats.binom.cdf(x - 1, n, p), "P(X > x)": stats.binom.sf(x, n, p), "mean": n * p, "sd": np.sqrt(n * p * (1 - p)), } if a is not None and b is not None: out["P(a <= X <= b)"] = stats.binom.pmf(np.arange(a, b + 1), n, p).sum() return pd.Series(out) def pois_calc(lam, x, interval=1.0): L = lam * interval return pd.Series({ "lambda'": L, "P(X = x)": stats.poisson.pmf(x, L), "P(X <= x)": stats.poisson.cdf(x, L), "P(X >= x)": stats.poisson.sf(x - 1, L), "mean": L, "variance": L, "sd": np.sqrt(L), }) print("\n===== 1. BINOMIAL CALCULATOR (n = 20, p = 0.06) =====") print(binom_calc(20, 0.06, x=3, a=1, b=4).round(6).to_string()) print("\n===== 2. POISSON CALCULATOR (lambda = 3.4 per hour) =====") print("Full hour:") print(pois_calc(3.4, 4).round(6).to_string()) print("\n20-minute window (lambda rescaled by 20/60):") print(pois_calc(3.4, 1, interval=20 / 60).round(6).to_string()) # ── 3. FIT THE BINOMIAL TO THE INSPECTION DATA ────────────────────────────── print("\n===== 3. BINOMIAL FIT — batch_defects.csv =====") n_items = int(defects["items_inspected"].iloc[0]) p_hat = defects["defective"].sum() / defects["items_inspected"].sum() print(f"batches inspected {len(defects)}") print(f"items per batch {n_items}") print(f"p-hat (pooled) {p_hat:.5f}") print(f"observed mean {defects['defective'].mean():.4f} " f"binomial np = {n_items * p_hat:.4f}") print(f"observed variance {defects['defective'].var(ddof=1):.4f} " f"binomial npq = {n_items * p_hat * (1 - p_hat):.4f}") obs_b = defects["defective"].value_counts().reindex(range(n_items + 1), fill_value=0) exp_b = stats.binom.pmf(np.arange(n_items + 1), n_items, p_hat) * len(defects) fit_b = pd.DataFrame({"defective": range(n_items + 1), "observed": obs_b.to_numpy(), "expected": exp_b.round(2)}) print(fit_b[(fit_b["observed"] > 0) | (fit_b["expected"] > 0.05)].to_string(index=False)) # ── 4. FIT THE POISSON TO THE CALL DATA ───────────────────────────────────── print("\n===== 4. POISSON FIT — calls_per_hour.csv =====") lam_hat = calls["calls"].mean() print(f"hours observed {len(calls)}") print(f"lambda-hat (mean) {lam_hat:.4f}") print(f"observed variance {calls['calls'].var(ddof=1):.4f}") print(f"dispersion index {calls['calls'].var(ddof=1)/lam_hat:.4f} " "(near 1 supports Poisson)") # ── 5. GOODNESS-OF-FIT TEST FOR THE POISSON ───────────────────────────────── print("\n===== 5. CHI-SQUARE GOODNESS OF FIT (Poisson) =====") N = len(calls) kmax = int(calls["calls"].max()) obs = calls["calls"].value_counts().reindex(range(kmax + 1), fill_value=0).to_numpy().astype(float) p_k = stats.poisson.pmf(np.arange(kmax + 1), lam_hat) p_k[-1] += stats.poisson.sf(kmax, lam_hat) # lump the upper tail exp_k = p_k * N lab = [str(i) for i in range(kmax + 1)] # combine sparse classes until every EXPECTED count is >= 5 while len(exp_k) > 2 and exp_k[-1] < 5: exp_k[-2] += exp_k[-1]; obs[-2] += obs[-1]; lab[-2] = lab[-2] + "+" exp_k, obs, lab = exp_k[:-1], obs[:-1], lab[:-1] while len(exp_k) > 2 and exp_k[0] < 5: exp_k[1] += exp_k[0]; obs[1] += obs[0]; lab[1] = "<=" + lab[1] exp_k, obs, lab = exp_k[1:], obs[1:], lab[1:] contrib = (obs - exp_k) ** 2 / exp_k chi2 = contrib.sum() df = len(obs) - 1 - 1 # -1 for the ESTIMATED lambda pval = stats.chi2.sf(chi2, df) print(pd.DataFrame({"calls": lab, "observed": obs.astype(int), "expected": exp_k.round(2), "contribution": contrib.round(4)}).to_string(index=False)) print(f"\nchi-square = {chi2:.4f} df = {df} p-value = {pval:.4f}") print(f"critical chi-square(0.05, {df}) = {stats.chi2.ppf(0.95, df):.4f}") print("=> FAIL TO REJECT. The Poisson model FITS -- a large p-value is the GOOD outcome here." if pval > 0.05 else "=> REJECT. The Poisson model does not fit these counts.") # ── 6. POISSON APPROXIMATES THE BINOMIAL ──────────────────────────────────── print("\n===== 6. POISSON AS A BINOMIAL APPROXIMATION =====") print("n = 2000, p = 0.0015 -> lambda = np = 3\n") xs = np.arange(9) cmp = pd.DataFrame({ "x": xs, "binomial": stats.binom.pmf(xs, 2000, 0.0015).round(6), "poisson": stats.poisson.pmf(xs, 3).round(6), }) cmp["difference"] = (cmp["binomial"] - cmp["poisson"]).round(6) print(cmp.to_string(index=False)) print(f"\nLargest absolute difference: {cmp['difference'].abs().max():.6f} -- negligible.") # ── 7. PLOTS ──────────────────────────────────────────────────────────────── fig, axes = plt.subplots(2, 2, figsize=(11, 8.5)) w = 0.4 k = np.arange(9) axes[0, 0].bar(k, stats.binom.pmf(k, 20, 0.06), color="#5B2A86") axes[0, 0].set(title="Binomial(20, 0.06)", xlabel="Defective", ylabel="P(x)") axes[0, 1].bar(k - w/2, obs_b.to_numpy()[:9], w, label="observed", color="#5B2A86") axes[0, 1].bar(k + w/2, exp_b[:9], w, label="expected", color="#0FA3A3") axes[0, 1].set(title="Defects: observed vs binomial", xlabel="Defective", ylabel="Batches") axes[0, 1].legend() kk = np.arange(13) axes[1, 0].bar(kk, stats.poisson.pmf(kk, lam_hat), color="#0B7A7A") axes[1, 0].set(title=f"Poisson({lam_hat:.2f})", xlabel="Calls", ylabel="P(x)") kf = np.arange(kmax + 1) obs_full = calls["calls"].value_counts().reindex(kf, fill_value=0).to_numpy() axes[1, 1].bar(kf - w/2, obs_full, w, label="observed", color="#5B2A86") axes[1, 1].bar(kf + w/2, stats.poisson.pmf(kf, lam_hat) * N, w, label="expected", color="#0FA3A3") axes[1, 1].set(title="Calls: observed vs Poisson", xlabel="Calls per hour", ylabel="Hours") axes[1, 1].legend() plt.tight_layout() plt.savefig("distribution_plots.png", dpi=110) print("\nWrote distribution_plots.png")