effectiveness.py

#!/usr/bin/env python3
"""Decide with statistics whether a corrective action worked.

Subcommands:
  rate    compare the defect rate before and after the action (exact conditional binomial test for
          Poisson counts: given k = k1 + k2 defects, k2 ~ Binomial(k, n2 / (n1 + n2)) when nothing changed)
  zero    units that must run with zero defects to show the rate is below a target at a confidence
          (n = ln(1 - C) / ln(1 - p)), and the upper bound reached so far
  upper   exact (Clopper-Pearson) upper confidence bound on the defect rate after the action

Verdicts:
  EFFECTIVE               rate fell significantly (p < alpha) and by at least --min-reduction
  NOT YET DEMONSTRATED    rate fell but not significantly - collect more data (zero shows how much)
  NOT EFFECTIVE           rate did not fall

Usage:
    python effectiveness.py rate --before 12 --before-units 4000 --after 1 --after-units 3500
    python effectiveness.py zero --target 0.001 --confidence 0.95 --so-far 1200
    python effectiveness.py upper --defects 1 --units 3500 --confidence 0.95

Exit codes: 0 effective / requirement met, 1 not effective or not yet demonstrated, 2 usage error.
"""

import argparse
import math
import sys


def binom_cdf(k, n, p):
    """P(X <= k) for X ~ Binomial(n, p), summed in log space for stability."""
    if p <= 0:
        return 1.0
    if p >= 1:
        return 1.0 if k >= n else 0.0
    total, lp, lq = 0.0, math.log(p), math.log1p(-p)
    for i in range(0, k + 1):
        total += math.exp(math.lgamma(n + 1) - math.lgamma(i + 1) - math.lgamma(n - i + 1) + i * lp + (n - i) * lq)
    return min(total, 1.0)


def upper_bound(k, n, conf):
    """Clopper-Pearson one-sided upper bound: smallest p with P(X <= k | p) <= 1 - conf."""
    if k >= n:
        return 1.0
    lo, hi = k / n, 1.0
    for _ in range(100):
        mid = (lo + hi) / 2
        if binom_cdf(k, n, mid) > 1 - conf:
            lo = mid
        else:
            hi = mid
    return hi


def usage(msg):
    print(f"error: {msg}", file=sys.stderr)
    return 2


def fmt_rate(r):
    return f"{r * 1e6:,.0f} ppm ({r:.4%})"


def cmd_rate(a):
    if a.before_units <= 0 or a.after_units <= 0 or a.before < 0 or a.after < 0:
        return usage("units must be > 0 and defects >= 0")
    r1, r2 = a.before / a.before_units, a.after / a.after_units
    k = a.before + a.after
    share = a.after_units / (a.before_units + a.after_units)
    p = binom_cdf(a.after, k, share) if k else 1.0
    red = 1 - r2 / r1 if r1 else 0.0
    print(f"Before: {a.before} defects in {a.before_units} units = {fmt_rate(r1)}")
    print(f"After:  {a.after} defects in {a.after_units} units = {fmt_rate(r2)}")
    print(f"Reduction: {red:.0%}   one-sided p-value (after < before): {p:.4f}")
    ub = upper_bound(a.after, a.after_units, 1 - a.alpha)
    print(f"After-rate upper {1 - a.alpha:.0%} bound: {fmt_rate(ub)}")
    if r2 >= r1:
        print("Verdict: NOT EFFECTIVE - the defect rate did not fall; re-open the root cause analysis")
        return 1
    if p < a.alpha and red >= a.min_reduction:
        print(f"Verdict: EFFECTIVE - significant at alpha {a.alpha:g} with {red:.0%} reduction (>= {a.min_reduction:.0%})")
        return 0
    if p < a.alpha:
        print(f"Verdict: NOT YET DEMONSTRATED - significant but reduction {red:.0%} < required {a.min_reduction:.0%}")
        return 1
    # units needed after the action, at the current after-rate pattern, to reach significance with 0 further defects
    extra = max(1, a.after_units // 20)
    while extra < 100 * a.after_units:
        s = (a.after_units + extra) / (a.before_units + a.after_units + extra)
        if binom_cdf(a.after, k, s) < a.alpha:
            break
        extra = int(extra * 1.1) + 1
    print(f"Verdict: NOT YET DEMONSTRATED - p {p:.3f} >= {a.alpha:g}; about {extra} more units with no further "
          "defects would show it")
    return 1


def cmd_zero(a):
    if not 0 < a.target < 1 or not 0 < a.confidence < 1:
        return usage("--target and --confidence must be between 0 and 1")
    n = math.ceil(math.log(1 - a.confidence) / math.log(1 - a.target))
    print(f"To show the defect rate is below {fmt_rate(a.target)} with {a.confidence:.0%} confidence: "
          f"{n} consecutive units with zero defects")
    if a.so_far is not None:
        ub = 1 - (1 - a.confidence) ** (1 / a.so_far) if a.so_far else 1.0
        left = max(0, n - a.so_far)
        print(f"So far {a.so_far} defect-free units: rate below {fmt_rate(ub)} at {a.confidence:.0%}; "
              + ("requirement met" if not left else f"{left} more units needed"))
        return 0 if not left else 1
    return 0


def cmd_upper(a):
    if a.units <= 0 or not 0 <= a.defects <= a.units:
        return usage("need 0 <= defects <= units and units > 0")
    ub = upper_bound(a.defects, a.units, a.confidence)
    print(f"{a.defects} defects in {a.units} units: observed {fmt_rate(a.defects / a.units)}; "
          f"{a.confidence:.0%} upper bound {fmt_rate(ub)}")
    if a.target is not None:
        ok = ub <= a.target
        print(f"Target {fmt_rate(a.target)}: " + ("demonstrated" if ok else "not demonstrated yet"))
        return 0 if ok else 1
    return 0


def main():
    ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    sub = ap.add_subparsers(dest="cmd", required=True)
    r = sub.add_parser("rate", help="before/after defect rate comparison")
    r.add_argument("--before", type=int, required=True)
    r.add_argument("--before-units", type=int, required=True)
    r.add_argument("--after", type=int, required=True)
    r.add_argument("--after-units", type=int, required=True)
    r.add_argument("--alpha", type=float, default=0.05)
    r.add_argument("--min-reduction", type=float, default=0.5, help="required fractional reduction (default 0.5)")
    z = sub.add_parser("zero", help="zero-defect demonstration size")
    z.add_argument("--target", type=float, required=True, help="maximum acceptable defect rate, e.g. 0.001")
    z.add_argument("--confidence", type=float, default=0.95)
    z.add_argument("--so-far", type=int, help="defect-free units already run")
    u = sub.add_parser("upper", help="upper confidence bound on the defect rate")
    u.add_argument("--defects", type=int, required=True)
    u.add_argument("--units", type=int, required=True)
    u.add_argument("--confidence", type=float, default=0.95)
    u.add_argument("--target", type=float)
    a = ap.parse_args()
    return {"rate": cmd_rate, "zero": cmd_zero, "upper": cmd_upper}[a.cmd](a)


if __name__ == "__main__":
    sys.exit(main())