#!/usr/bin/env python3
"""Operating characteristic of a single sampling plan, or design a plan from producer and consumer risk.
Subcommands:
oc probability of acceptance Pa(p) for plan (n, c); the p at Pa = 0.95 (producer's quality level),
at Pa = 0.10 (LTPD / consumer's quality level), and AOQ / AOQL under rectifying inspection.
Binomial model, or hypergeometric with --lot (small lots, sampling without replacement).
design smallest (n, c) with Pa(AQL) >= 1 - alpha and Pa(LTPD) <= beta (binomial).
c0 zero-acceptance plan: smallest n with Pa(LTPD) <= beta at c = 0 (n = ln(beta) / ln(1 - LTPD)).
Usage:
python oc_curve.py oc --n 80 --c 2
python oc_curve.py oc --n 80 --c 2 --lot 1000 --svg > oc.svg
python oc_curve.py design --aql 0.01 --ltpd 0.05 [--alpha 0.05 --beta 0.10]
python oc_curve.py c0 --ltpd 0.05 --beta 0.10
Exit codes: 0 computed, 1 no plan found within limits, 2 usage error.
"""
import argparse
import math
import sys
def binom_cdf(c, n, p):
if p <= 0:
return 1.0
if p >= 1:
return 0.0 if c < n else 1.0
lp, lq = math.log(p), math.log1p(-p)
return min(1.0, sum(math.exp(math.lgamma(n + 1) - math.lgamma(i + 1) - math.lgamma(n - i + 1) + i * lp + (n - i) * lq)
for i in range(0, min(c, n) + 1)))
def hyper_cdf(c, n, N, p):
D = round(p * N)
tot = math.comb(N, n)
return sum(math.comb(D, i) * math.comb(N - D, n - i) for i in range(0, min(c, D, n) + 1)) / tot
def pa(c, n, p, N=None):
return hyper_cdf(c, n, N, p) if N else binom_cdf(c, n, p)
def solve(c, n, target, N=None):
lo, hi = 0.0, 1.0
for _ in range(80):
mid = (lo + hi) / 2
if pa(c, n, mid, N) > target:
lo = mid
else:
hi = mid
return (lo + hi) / 2
def svg(points, n, c):
W, H, M = 520, 320, 50
pmax = points[-1][0]
X = lambda p: M + (W - 2 * M) * p / pmax
Y = lambda v: H - M - (H - 2 * M) * v
path = " ".join(f"{X(p):.1f},{Y(v):.1f}" for p, v in points)
ticks = "".join(f'{100 * pmax * i / 5:.1f}%' for i in range(6))
yt = "".join(f'{v:.1f}' for v in (0, 0.25, 0.5, 0.75, 1.0))
return (f'')
def main():
ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
sub = ap.add_subparsers(dest="cmd", required=True)
o = sub.add_parser("oc")
o.add_argument("--n", type=int, required=True)
o.add_argument("--c", type=int, required=True)
o.add_argument("--lot", type=int, help="lot size for the hypergeometric model and AOQ")
o.add_argument("--svg", action="store_true", help="print the OC curve as SVG instead of a table")
d = sub.add_parser("design")
d.add_argument("--aql", type=float, required=True, help="producer's quality level as a fraction, e.g. 0.01")
d.add_argument("--ltpd", type=float, required=True, help="consumer's quality level as a fraction, e.g. 0.05")
d.add_argument("--alpha", type=float, default=0.05)
d.add_argument("--beta", type=float, default=0.10)
z = sub.add_parser("c0")
z.add_argument("--ltpd", type=float, required=True)
z.add_argument("--beta", type=float, default=0.10)
a = ap.parse_args()
if a.cmd == "oc":
if a.n < 1 or a.c < 0 or a.c >= a.n or (a.lot and a.lot < a.n):
print("error: need n >= 1, 0 <= c < n and lot >= n", file=sys.stderr)
return 2
p95, p10, p50 = solve(a.c, a.n, 0.95, a.lot), solve(a.c, a.n, 0.10, a.lot), solve(a.c, a.n, 0.50, a.lot)
grid = [p10 * 1.6 * i / 40 for i in range(41)]
if a.svg:
print(svg([(p, pa(a.c, a.n, p, a.lot)) for p in grid], a.n, a.c))
return 0
print(f"Plan n = {a.n}, c = {a.c}" + (f", lot {a.lot} (hypergeometric)" if a.lot else " (binomial)"))
print(f" Pa = 0.95 at p = {p95:.4%} (producer's quality level)")
print(f" Pa = 0.50 at p = {p50:.4%} (indifference quality)")
print(f" Pa = 0.10 at p = {p10:.4%} (LTPD / consumer's quality level)\n")
print(f"{'p':>9}{'Pa':>9}" + (f"{'AOQ':>10}" if a.lot else ""))
aoqs = []
for p in grid[::4]:
v = pa(a.c, a.n, p, a.lot)
aoq = p * v * (a.lot - a.n) / a.lot if a.lot else None
aoqs.append(aoq or 0)
print(f"{p:>9.3%}{v:>9.3f}" + (f"{aoq:>10.3%}" if a.lot else ""))
if a.lot:
fine = [p * pa(a.c, a.n, p, a.lot) * (a.lot - a.n) / a.lot for p in grid]
print(f"\nAOQL (worst average outgoing quality with rejected lots 100% sorted): {max(fine):.3%}")
return 0
if a.cmd == "design":
if not 0 < a.aql < a.ltpd < 1:
print("error: need 0 < aql < ltpd < 1 (fractions)", file=sys.stderr)
return 2
for c in range(0, 60):
# smallest n that meets the consumer risk at this c
n = c + 1
while binom_cdf(c, n, a.ltpd) > a.beta:
n += 1 if n < 200 else max(1, n // 50)
if n > 50000:
break
if n > 50000:
break
if binom_cdf(c, n, a.aql) >= 1 - a.alpha:
print(f"Plan: n = {n}, c = {c}")
print(f" Pa at AQL {a.aql:.3%} = {binom_cdf(c, n, a.aql):.3f} (needs >= {1 - a.alpha:.2f})")
print(f" Pa at LTPD {a.ltpd:.3%} = {binom_cdf(c, n, a.ltpd):.3f} (needs <= {a.beta:.2f})")
return 0
print("No plan with n <= 50000 meets both risks; widen the gap between AQL and LTPD.", file=sys.stderr)
return 1
if not 0 < a.ltpd < 1 or not 0 < a.beta < 1:
print("error: --ltpd and --beta must be fractions between 0 and 1", file=sys.stderr)
return 2
n = math.ceil(math.log(a.beta) / math.log(1 - a.ltpd))
print(f"c = 0 plan: n = {n}; accept the lot only if no nonconforming unit is found.")
print(f" Lots at {a.ltpd:.2%} nonconforming are accepted with probability {binom_cdf(0, n, a.ltpd):.3f}")
print(f" A lot at {a.ltpd / 10:.2%} is still rejected {1 - binom_cdf(0, n, a.ltpd / 10):.1%} of the time: c = 0 plans are strict on producers.")
return 0
if __name__ == "__main__":
sys.exit(main())