#!/usr/bin/env python3 """Process capability: Cp, Cpk, Pp, Ppk, Cpm, expected and observed ppm, a confidence interval for Cpk, and an Anderson-Darling normality check, from a column of measurements. Within-subgroup sigma (Cp, Cpk): R-bar/d2 with --subgroup or --size (n = 2..10), pooled S for n > 10, or MR-bar/1.128 for individuals. Overall sigma (Pp, Ppk): sample standard deviation. One-sided specifications are supported (give only --lsl or --usl). --method percentile uses the 0.135% and 99.865% empirical percentiles (ISO 22514-2 style) for non-normal data; it needs several hundred values. Usage: python capability.py bore.csv --value diameter --lsl 12.000 --usl 12.018 --subgroup hour python capability.py flatness.csv --value flatness --usl 0.05 --target-index 1.67 python capability.py data.csv --value x --lsl 10 --usl 20 --target 15 --method percentile --json Exit codes: 0 Cpk (or Ppk) meets --target-index, 1 below target or not normal, 2 usage or data error. """ import argparse import csv import json import math import sys from statistics import NormalDist, mean, stdev D2 = {2: 1.128, 3: 1.693, 4: 2.059, 5: 2.326, 6: 2.534, 7: 2.704, 8: 2.847, 9: 2.970, 10: 3.078} N = NormalDist() def c4(n): return math.sqrt(2 / (n - 1)) * math.exp(math.lgamma(n / 2) - math.lgamma((n - 1) / 2)) def within_sigma(values, groups): if groups and len(groups[0]) > 1: n = len(groups[0]) if n <= 10: return mean(max(g) - min(g) for g in groups) / D2[n], f"R-bar/d2 (n={n})" s2 = sum(stdev(g) ** 2 * (len(g) - 1) for g in groups) / sum(len(g) - 1 for g in groups) return math.sqrt(s2) / c4(sum(len(g) - 1 for g in groups) + 1), f"pooled S (n={n})" mr = [abs(b - a) for a, b in zip(values, values[1:])] return mean(mr) / D2[2], "MR-bar/1.128 (individuals)" def anderson_darling(values): """A^2* (small-sample adjusted) and approximate p-value (D'Agostino & Stephens 1986).""" x = sorted(values) n = len(x) m, s = mean(x), stdev(x) if s == 0: return float("inf"), 0.0 cdf = [min(max(N.cdf((v - m) / s), 1e-12), 1 - 1e-12) for v in x] a2 = -n - sum((2 * i + 1) * (math.log(cdf[i]) + math.log(1 - cdf[n - 1 - i])) for i in range(n)) / n a = a2 * (1 + 0.75 / n + 2.25 / n ** 2) if a >= 0.6: p = math.exp(1.2937 - 5.709 * a + 0.0186 * a * a) elif a >= 0.34: p = math.exp(0.9177 - 4.279 * a - 1.38 * a * a) elif a >= 0.2: p = 1 - math.exp(-8.318 + 42.796 * a - 59.938 * a * a) else: p = 1 - math.exp(-13.436 + 101.14 * a - 223.73 * a * a) return a, min(max(p, 0.0), 1.0) def percentile(sorted_x, q): pos = q * (len(sorted_x) - 1) lo = math.floor(pos) hi = min(lo + 1, len(sorted_x) - 1) return sorted_x[lo] + (sorted_x[hi] - sorted_x[lo]) * (pos - lo) def indices(mu, sigma, lsl, usl): cp = (usl - lsl) / (6 * sigma) if lsl is not None and usl is not None else None cpl = (mu - lsl) / (3 * sigma) if lsl is not None else None cpu = (usl - mu) / (3 * sigma) if usl is not None else None k = min(v for v in (cpl, cpu) if v is not None) ppm = 1e6 * ((N.cdf((lsl - mu) / sigma) if lsl is not None else 0) + (1 - N.cdf((usl - mu) / sigma) if usl is not None else 0)) return cp, k, ppm def main(): ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) ap.add_argument("csv") ap.add_argument("--value", required=True) ap.add_argument("--lsl", type=float) ap.add_argument("--usl", type=float) ap.add_argument("--target", type=float, help="nominal for Cpm") ap.add_argument("--subgroup", help="subgroup column for within sigma") ap.add_argument("--size", type=int, help="consecutive subgroup size for within sigma") ap.add_argument("--target-index", type=float, default=1.33, help="required Cpk/Ppk (default 1.33; 1.67 for new processes)") ap.add_argument("--confidence", type=float, default=0.95) ap.add_argument("--method", choices=["normal", "percentile"], default="normal") ap.add_argument("--json", action="store_true") a = ap.parse_args() if a.lsl is None and a.usl is None: print("error: give --lsl and/or --usl", file=sys.stderr) return 2 if a.lsl is not None and a.usl is not None and a.lsl >= a.usl: print("error: --lsl must be below --usl", file=sys.stderr) return 2 try: with open(a.csv, newline="", encoding="utf-8-sig") as fh: rows = list(csv.DictReader(fh)) vals = [float(r[a.value]) for r in rows if (r.get(a.value) or "").strip()] except (OSError, KeyError, ValueError) as e: print(f"error: {e}", file=sys.stderr) return 2 if len(vals) < 25: print(f"error: {len(vals)} values; capability needs at least 25 (100+ recommended)", file=sys.stderr) return 2 groups = None if a.subgroup: g = {} for r in rows: if (r.get(a.value) or "").strip(): g.setdefault(r[a.subgroup], []).append(float(r[a.value])) groups = [v for v in g.values() if len(v) == len(next(iter(g.values())))] elif a.size and a.size > 1: groups = [vals[i:i + a.size] for i in range(0, len(vals) - a.size + 1, a.size)] n, mu, s_all = len(vals), mean(vals), stdev(vals) s_within, how = within_sigma(vals, groups) cp, cpk, ppm_w = indices(mu, s_within, a.lsl, a.usl) pp, ppk, ppm_o = indices(mu, s_all, a.lsl, a.usl) obs = sum((a.lsl is not None and v < a.lsl) or (a.usl is not None and v > a.usl) for v in vals) z = N.inv_cdf(1 - (1 - a.confidence) / 2) half = z * math.sqrt(1 / (9 * n) + cpk ** 2 / (2 * (n - 1))) # Bissell (1990) ad, p_ad = anderson_darling(vals) res = {"n": n, "mean": mu, "sigma_within": s_within, "sigma_within_method": how, "sigma_overall": s_all, "Cp": cp, "Cpk": cpk, "Pp": pp, "Ppk": ppk, "ppm_within": ppm_w, "ppm_overall": ppm_o, "ppm_observed": 1e6 * obs / n, "Cpk_ci": [cpk - half, cpk + half], "AD": ad, "AD_p": p_ad} if a.target is not None and cp is not None: res["Cpm"] = (a.usl - a.lsl) / (6 * math.sqrt(s_all ** 2 + (mu - a.target) ** 2)) if a.method == "percentile": xs = sorted(vals) p_lo, p50, p_hi = percentile(xs, 0.00135), percentile(xs, 0.5), percentile(xs, 0.99865) parts = [] if a.lsl is not None: parts.append((p50 - a.lsl) / (p50 - p_lo)) if a.usl is not None: parts.append((a.usl - p50) / (p_hi - p50)) res["Ppk_percentile"] = min(parts) if cp is not None: res["Pp_percentile"] = (a.usl - a.lsl) / (p_hi - p_lo) key = "Ppk_percentile" if a.method == "percentile" else "Cpk" ok = res[key] >= a.target_index and (a.method == "percentile" or p_ad >= 0.05) if a.json: print(json.dumps(res, indent=2)) return 0 if ok else 1 f = lambda v: "-" if v is None else f"{v:.2f}" print(f"n = {n} mean = {mu:.6g} sigma within = {s_within:.5g} [{how}] sigma overall = {s_all:.5g}") print(f"Spec: LSL {a.lsl if a.lsl is not None else '-'} USL {a.usl if a.usl is not None else '-'}\n") print(f"Cp {f(cp):>6} Cpk {f(cpk):>6} ({a.confidence:.0%} CI {cpk - half:.2f} to {cpk + half:.2f})") print(f"Pp {f(pp):>6} Ppk {f(ppk):>6}" + (f" Cpm {res['Cpm']:.2f}" if "Cpm" in res else "")) if a.method == "percentile": print(f"Percentile method: Pp {f(res.get('Pp_percentile'))} Ppk {res['Ppk_percentile']:.2f}") print(f"Expected ppm: within {ppm_w:,.0f} overall {ppm_o:,.0f} observed {res['ppm_observed']:,.0f}") print(f"Normality (Anderson-Darling): A2* = {ad:.3f}, p = {p_ad:.3f}" + (" -> NOT normal" if p_ad < 0.05 else "")) if p_ad < 0.05 and a.method == "normal": print(" Normal-based indices and ppm are unreliable: find the cause (mixture, drift, skew), or use " "--method percentile / a Box-Cox transform.") if cp is not None and cpk < 0.75 * cp: print(f" Process is off-centre: centring at {(a.lsl + a.usl) / 2:.6g} would raise Cpk to about {cp:.2f}.") if ppk < 0.85 * cpk: print(" Ppk well below Cpk: between-subgroup variation (shifts, drift) - check the control chart.") if obs and res["ppm_observed"] > 3 * max(ppm_o, 1): print(f" Observed {obs} out-of-spec value(s) - far more than the normal model predicts: special causes or heavy tails.") print(" Indices are valid only for a process in statistical control: check the control chart first.") label = "Ppk (percentile)" if a.method == "percentile" else "Cpk" print(f"\n{label} {res[key]:.2f} vs required {a.target_index:g}: " + ("CAPABLE" if ok else "NOT CAPABLE / NOT VALID")) return 0 if ok else 1 if __name__ == "__main__": sys.exit(main())