From 16d80d1a87b26458f5a2760b6442526efc7319f8 Mon Sep 17 00:00:00 2001 From: Ramzi Abujamra <31161729+rabujamra@users.noreply.github.com> Date: Thu, 10 Sep 2026 11:50:48 -0500 Subject: [PATCH] add short-series sampling and band scripts MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Adds the two reproducibility scripts from #114: the original short-series sampling study and the n=3–30 band extension. --- validation/short_series_bands.py | 88 ++++++++++++++++++++++++ validation/short_series_sampling.py | 100 ++++++++++++++++++++++++++++ 2 files changed, 188 insertions(+) create mode 100644 validation/short_series_bands.py create mode 100644 validation/short_series_sampling.py diff --git a/validation/short_series_bands.py b/validation/short_series_bands.py new file mode 100644 index 0000000..cdb00b9 --- /dev/null +++ b/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.") diff --git a/validation/short_series_sampling.py b/validation/short_series_sampling.py new file mode 100644 index 0000000..9ee209e --- /dev/null +++ b/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()