Add short-series sampling and band scripts.

validation/short_series_bands.py

@@ -0,0 +1,88 @@
+"""Sampling distribution of MRbar/d2 for every n from 3 to 30.
+
+Companion to short_series_sampling.py, for issue #114. Self-contained: numpy
+and processbehavior. No data files.
+
+    python short_series_bands.py
+
+WHY THIS IS A SEPARATE SCRIPT WITH A DIFFERENT RNG DESIGN
+---------------------------------------------------------
+`short_series_sampling.py` draws all its lengths from ONE stream, so its rows
+depend on the order of the loop -- which is why the n=5 row has to stay in it.
+That is fine for reproducing a fixed table and bad for a table someone will cut
+thresholds out of.
+
+Here each length gets its OWN independent stream, `default_rng([SEED, n])`, so
+every row is reproducible on its own and no row depends on which other lengths
+were computed. Rows therefore agree with the original table to within Monte
+Carlo error rather than exactly; the agreement is checked and printed below.
+
+Replicates raised to 100,000 so the band edges are not chosen off noise. The
+Monte Carlo standard error of each reported quantile is given, because a
+threshold placed where two adjacent n differ by less than their MC error is a
+threshold placed on nothing.
+
+NO CRITERION IS PROPOSED HERE. Quantiles and spreads only. Which functional of
+this distribution should govern an adequacy band -- p90/p10, RSD, a coverage
+probability, something else -- and where the cut falls are design decisions for
+the maintainers. This supplies the arithmetic once that choice is made.
+"""
+import numpy as np
+from processbehavior.spc_constants import D2_N2
+
+SEED = 20260903
+REPS = 100_000
+LENGTHS = range(3, 31)
+QS = [5, 10, 25, 50, 75, 90, 95]
+
+
+def sweep(n, reps=REPS, seed=SEED):
+    rng = np.random.default_rng([seed, n])
+    z = rng.standard_normal((reps, n))
+    return np.abs(np.diff(z, axis=1)).mean(axis=1) / D2_N2
+
+
+def q_se(v, q):
+    """MC standard error of the q-th percentile, via the density at that point."""
+    p = q / 100.0
+    x = np.percentile(v, q)
+    h = 1.06 * v.std(ddof=1) * len(v) ** -0.2          # Silverman bandwidth
+    dens = np.mean(np.abs(v - x) < h) / (2 * h)
+    return np.sqrt(p * (1 - p) / len(v)) / dens if dens > 0 else float("nan")
+
+
+if __name__ == "__main__":
+    import processbehavior as pb
+    print("processbehavior %s | numpy %s" % (pb.__version__, np.__version__))
+    print("MRbar/d2 under an iid standard normal process, sigma = 1.0")
+    print("%d replicates per length, independent stream per n "
+          "(default_rng([%d, n]))\n" % (REPS, SEED))
+
+    hdr = "%4s" % "n" + "".join("%8s" % ("p%d" % q) for q in QS) + \
+          "%9s%8s%10s%10s" % ("p90/p10", "RSD", "se(p10)", "se(p90)")
+    print(hdr); print("-" * len(hdr))
+    rows = {}
+    for n in LENGTHS:
+        v = sweep(n)
+        qv = np.percentile(v, QS)
+        rsd = v.std(ddof=1) / v.mean()
+        rows[n] = dict(zip(QS, qv))
+        print("%4d" % n + "".join("%8.3f" % x for x in qv) +
+              "%8.2fx%8.0f%%%10.4f%10.4f"
+              % (qv[QS.index(90)] / qv[QS.index(10)], 100 * rsd,
+                 q_se(v, 10), q_se(v, 90)))
+
+    print("\nAGREEMENT WITH THE TABLE IN THE ISSUE (single-stream, 20k reps)")
+    pub = {3: (0.889, 0.337, 1.797), 4: (0.925, 0.429, 1.689),
+           5: (0.942, 0.482, 1.585), 8: (0.969, 0.592, 1.454),
+           12: (0.976, 0.671, 1.351), 25: (0.989, 0.769, 1.245)}
+    print("%4s %22s %22s %10s" % ("n", "published (p50/p10/p90)",
+                                  "here (p50/p10/p90)", "max |diff|"))
+    for n, (m, a, b) in pub.items():
+        h = (rows[n][50], rows[n][10], rows[n][90])
+        print("%4d %22s %22s %10.4f"
+              % (n, "%.3f / %.3f / %.3f" % (m, a, b),
+                 "%.3f / %.3f / %.3f" % h,
+                 max(abs(h[0] - m), abs(h[1] - a), abs(h[2] - b))))
+    print("\nDifferences are Monte Carlo, not a discrepancy: the two scripts use\n"
+          "different streams by design.")

validation/short_series_sampling.py

@@ -0,0 +1,100 @@
+"""Short-series behaviour of MRbar/d2, and why a drift test cannot rescue it.
+
+Reproduces every simulated number in issue #114 on cnicholas/processbehavior.
+Self-contained: numpy, pandas, processbehavior. Nothing else, no data files.
+
+    python short_series_sampling.py
+
+Sections
+  1  limit width vs series length, on one stable simulated process   (seed 7)
+  2  sampling distribution of MRbar/d2 by n                (rng 20260903)
+  3  the drift test's null under T iid normal draws       (rng 20260903)
+"""
+import numpy as np, pandas as pd
+import processbehavior as pb
+from processbehavior.spc_constants import D2_N2
+
+MEAN, SD = 5.50, 1.20
+TRUE_WIDTH = 6.0 * SD                      # natural process limits span 6 sigma
+LENGTHS = (3, 4, 6, 8, 12, 20, 30, 60, 150)
+NS = (3, 4, 5, 8, 12, 25)   # n=5 was computed but omitted from the table in #114
+REPS_SAMPLING = 20_000
+REPS_NULL = 400_000
+PRECISION = 12                             # do not let rounding move the answer
+
+
+def x_limits(y, precision=PRECISION):
+    """Ask the library for the X-chart limits on this series.
+
+    precision=12 on every call. At the default of 3 the limits are rounded
+    before the width is taken, which moves section 1 at these scales.
+    """
+    df = pd.DataFrame({"t": range(len(y)), "y": np.asarray(y, float)})
+    st = pb.formulate(df, response="y", time="t", precision=precision)
+    r = st.execute(chart="X", companion=True)
+    s = r.get_statistics("X")
+    return st, s
+
+
+def section_1():
+    print("1. LIMIT WIDTH vs SERIES LENGTH")
+    print("   one stable process, mean %.2f sd %.2f, true limit span %.2f\n"
+          % (MEAN, SD, TRUE_WIDTH))
+    y = np.random.default_rng(7).normal(MEAN, SD, max(LENGTHS))
+    print("   %4s %12s %14s %6s %12s" % ("T", "limit width", "share of true", "ADS", "recommended"))
+    for T in LENGTHS:
+        st, s = x_limits(y[:T])
+        w = float(s["upl"]) - float(s["lpl"])
+        print("   %4d %12.2f %13.0f%% %6d %12s"
+              % (T, w, 100 * w / TRUE_WIDTH,
+                 int(st.analytical_design_state.sds), str(st.recommended_chart)))
+    print()
+
+
+def mrbar_over_d2(y):
+    return np.abs(np.diff(y)).mean() / D2_N2
+
+
+def section_2():
+    print("2. SAMPLING DISTRIBUTION OF MRbar/d2, sigma fixed at 1.0")
+    print("   %d replicates per length\n" % REPS_SAMPLING)
+    rng = np.random.default_rng(20260903)
+    print("   %4s %8s %8s %8s %9s %6s" % ("n", "median", "p10", "p90", "p90/p10", "RSD"))
+    for n in NS:
+        draws = rng.standard_normal((REPS_SAMPLING, n))
+        v = np.abs(np.diff(draws, axis=1)).mean(axis=1) / D2_N2
+        p10, p50, p90 = np.percentile(v, [10, 50, 90])
+        print("   %4d %8.3f %8.3f %8.3f %8.2fx %5.0f%%"
+              % (n, p50, p10, p90, p90 / p10, 100 * v.std(ddof=1) / v.mean()))
+    print()
+
+
+def drift_share(y):
+    d = np.diff(np.asarray(y, float))
+    return abs(d.mean()) / np.abs(d).mean()
+
+
+def section_3():
+    print("3. THE DRIFT TEST'S NULL")
+    print("   %d replicates of T iid standard normal draws." % REPS_NULL)
+    print("   For T iid observations from a continuous distribution, all T!")
+    print("   orderings are equally likely, so P(monotone) = 2/T!\n")
+    rng = np.random.default_rng(20260903)
+    print("   %4s %12s %14s %14s" % ("T", "P(monotone)", "theory 2/T!", "95th pct share"))
+    for T in (4, 5, 6):
+        z = rng.standard_normal((REPS_NULL, T))
+        d = np.diff(z, axis=1)
+        mono = ((d > 0).all(axis=1) | (d < 0).all(axis=1)).mean()
+        share = np.abs(d.mean(axis=1)) / np.abs(d).mean(axis=1)
+        theory = 2.0 / __import__("math").factorial(T)
+        print("   %4d %11.4f %14.4f %14.3f"
+              % (T, mono, theory, np.percentile(share, 95)))
+    print("\n   At T=4, P(monotone)=1/12=8.3% > 5%, so the 95th percentile of the")
+    print("   null is 1.000 -- the maximum attainable. No four-point series can")
+    print("   exceed it, so the test cannot fire.\n")
+
+
+if __name__ == "__main__":
+    print("processbehavior", pb.__version__, "| numpy", np.__version__,
+          "| pandas", pd.__version__, "\n")
+    section_1(); section_2(); section_3()