diff --git a/CHANGELOG.md b/CHANGELOG.md index 99749dc..807369b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -35,6 +35,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 the residuals guide. ### Changed +- **Manual citations follow the 11-7 edition.** Dr. Bishop's manual of 2026-11-07 has 20 chapters + (chapters 7 and up are one lower than before) and dot-numbered captions; every equation and + chapter citation in the code, docs and tests now uses its numbering, in the hyphen form the + manual repository checks ("11-7 manual Eq 13-8"). Seven equations it removed are cited by + their replacements. No behaviour changes. - **mR charts check Rule 1 only by default.** Wheeler reads the moving range chart for points above its upper limit. The zone and run rules assume symmetric, independent points, and successive moving ranges share a reading, so with all eight rules on, 95% of 100-point diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index 5e06a20..0229394 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -108,7 +108,7 @@ compares outputs, not prose. ## Citing the methodology Code and docs cite Bishop's VAS documentation manual by edition and equation number: -`(10-12 manual Eq 14-8)`, `Eqs 14-4..14-13`, `Eqs 11-16/17`. Name the edition once per file and use +`(11-7 manual Eq 13-8)`, `Eqs 13-4..13-13`, `Eqs 10-16/17`. Name the edition once per file and use the chapter-number form; `tests/test_manual_citations.py` checks both. The citations are read against the manual itself, so a new edition reports every place in the library whose equation changed or moved. diff --git a/docs/appendix/wheeler-terminology.md b/docs/appendix/wheeler-terminology.md index 7799a42..5ddbc38 100644 --- a/docs/appendix/wheeler-terminology.md +++ b/docs/appendix/wheeler-terminology.md @@ -232,11 +232,11 @@ Wheeler uses the standard SPC constants from Shewhart's work: | Constant | Purpose | Formula Source | |----------|---------|----------------| -| c₄ | Unbiasing s | Related to gamma function; Bishop's αN (10-12 manual Eqs 10-11, 10-12) | +| c₄ | Unbiasing s | Related to gamma function; Bishop's αN (11-7 manual Eqs 9-8, 9-9) | | A₃ | Xbar limits from s | 3 / (c₄√n) | -| B₃, B₄ | S chart limits | Functions of c₄ and n (Eqs 10-14, 11-8, 11-9) | -| d₂ | Unbiasing range | Tabulated; 1.128 for moving ranges of two (Eqs 12-2, 12-4) | -| D₃, D₄ | mR chart limits | Functions of d₂ and d₃; D₄ = 1 + 3(0.8525)/1.128 for moving ranges of two (Eq 12-5) | +| B₃, B₄ | S chart limits | Functions of c₄ and n (Eqs 9-11, 10-8, 10-9) | +| d₂ | Unbiasing range | Tabulated; 1.128 for moving ranges of two (Eqs 11-2, 11-4) | +| D₃, D₄ | mR chart limits | Functions of d₂ and d₃; D₄ = 1 + 3(0.8525)/1.128 for moving ranges of two (Eq 11-5) | --- diff --git a/docs/reference/api.md b/docs/reference/api.md index ba9e3c8..6ea2f16 100644 --- a/docs/reference/api.md +++ b/docs/reference/api.md @@ -420,7 +420,7 @@ From `study.capability()`. Fields include `pp`, `ppk` (with `ppk_lower`/`ppk_upp `stability_warning`, and the capability `window` when one was requested. - `sigma_hat_r2` — SR2 / c4(NR2), the potential σ̂ from R2 - (10-12 manual Eq 16-10). + (11-7 manual Eq 15-10). - `potential_center` — the mean of the potential values ȳ + R2, i.e. ȳ + mean(R2). Potential CPL/CPU are measured from it and the potential chart's center line sits there; it equals `y_bar` in ADS 1. `None` when potential capability is unavailable. @@ -445,10 +445,10 @@ From `study.loss_function()`. Taguchi loss decomposition: `centering`, counterparts, `pdc_by_factor`, and the `target` used (`target_is_default` says whether it was derived). -`unexplained` is the average cell variance in ADS 1 (10-12 manual Eq 15-16) and the -sample variance of R2 in ADS 2/3 (Eq 15-17). `interaction` is the average squared +`unexplained` is the average cell variance in ADS 1 (11-7 manual Eq 14-14) and the +sample variance of R2 in ADS 2/3 (Eq 14-15). `interaction` is the average squared interaction effect in ADS 1; in ADS 2/3 it is what the other parts leave of -(Ȳ − T)² + S², set to 0 if negative (Eqs 15-18..15-20), so there `total` equals +(Ȳ − T)² + S², set to 0 if negative (Eqs 14-16..14-18), so there `total` equals (Ȳ − T)² + S² unless the interaction was set to 0. `loss_function()` raises `ValidationError` when R2 is unavailable for the layout. @@ -467,8 +467,8 @@ loss.plot(*, structured=False, orientation='vertical', theme=None, From `study.maximum_information()`. Fields: `n`, `r2_mean`, `r2_mR`, `sigma_hat`, `upl`, `lpl`, `n_signals`. An X chart of R2: limits -`r2_mean ± (3/1.128)·r2_mR`, `sigma_hat = r2_mR / 1.128` (10-12 manual Eqs 12-4, 12-10, -12-11). Raises +`r2_mean ± (3/1.128)·r2_mR`, `sigma_hat = r2_mR / 1.128` (11-7 manual Eqs 11-4, 11-8, +11-9). Raises `ValidationError` when R2 is unavailable for the layout. ```python @@ -614,9 +614,9 @@ Secondary surface. The control-chart constants are functions of subgroup size: ```python from processbehavior.spc_constants import c4, b3, b4, VALID_BASE_CHARTS -c4(n) # bias-correction factor for s (10-12 manual Eqs 10-11, 10-12); exact for every n -b3(n, sigma_multiplier=3) # S-chart lower-limit factor (Eqs 10-14, 11-8) -b4(n, sigma_multiplier=3) # S-chart upper-limit factor (Eqs 10-14, 11-9) +c4(n) # bias-correction factor for s (11-7 manual Eqs 9-8, 9-9); exact for every n +b3(n, sigma_multiplier=3) # S-chart lower-limit factor (Eqs 9-11, 10-8) +b4(n, sigma_multiplier=3) # S-chart upper-limit factor (Eqs 9-11, 10-9) VALID_BASE_CHARTS # {'Histogram', 'Xbar', 'S', 'X', 'mR'} ``` diff --git a/docs/reference/sds_definitions.md b/docs/reference/sds_definitions.md index 3c80bbd..4a0df9f 100644 --- a/docs/reference/sds_definitions.md +++ b/docs/reference/sds_definitions.md @@ -41,7 +41,7 @@ For background on the lineage itself (Planned → Observed → Analytical), see - **Grid Status**: Semi-Complete - **Description**: Every (factor × time) cell has exactly 1 observation - **Implications**: - - No within-cell variance: R2 is estimated from the condition-and-period-adjusted series, differenced along the condition-then-time sequence and scaled by the R2 scale factor (Eqs 14-4..14-8) + - No within-cell variance: R2 is estimated from the condition-and-period-adjusted series, differenced along the condition-then-time sequence and scaled by the R2 scale factor (Eqs 13-4..13-8) - Interaction effects confounded with pure error - VAS residuals available when the R2 scale factor is defined: at least 2 process design conditions and M ≥ K + 2 observations (see [R2 calculation](#r2-calculation-methods-by-ds)) - Common in unreplicated factorial designs @@ -52,7 +52,7 @@ For background on the lineage itself (Planned → Observed → Analytical), see - **Grid Status**: Semi-Complete - **Description**: Mix of cells - some have n=1, others have n≥2 - **Implications**: - - R2 by the same scaled difference as DS 2, over every observation (Eqs 14-9..14-13); any one-observation cell selects `ma2` for every observation + - R2 by the same scaled difference as DS 2, over every observation (Eqs 13-9..13-13); any one-observation cell selects `ma2` for every observation - Xbar/S is the default chart - **Most common in real-world data** - Partial interaction effect estimation @@ -123,12 +123,12 @@ Equation numbers refer to Bishop's VAS documentation manual dated 10-12-2026. | DS | R2 Method | Description | |-----|-----------|-------------| | 0 | N/A | No VAS decomposition | -| 1 | Exact (within-cell) | `R2 = Y - Ȳ_kt` (Eq 14-3) | -| 2 | Scaled difference (`ma2`) | Condition-and-period-adjusted series, differenced and divided by 2·c(K, M), M = KT (Eqs 14-4..14-8) | -| 3 | Scaled difference (`ma2`) | Any one-observation cell: the same over every observation (Eqs 14-9..14-13) | +| 1 | Exact (within-cell) | `R2 = Y - Ȳ_kt` (Eq 13-3) | +| 2 | Scaled difference (`ma2`) | Condition-and-period-adjusted series, differenced and divided by 2·c(K, M), M = KT (Eqs 13-4..13-8) | +| 3 | Scaled difference (`ma2`) | Any one-observation cell: the same over every observation (Eqs 13-9..13-13) | | 4 | Exact (present cells) | Collapses to ADS 1 after cleaning; within-cell for available data | -| 5 | Scaled difference (`ma2`) | Collapses to ADS 2 after cleaning; all present cells are unreplicated (Eqs 14-9..14-13) | -| 6 | Scaled difference (`ma2`) | Collapses to ADS 3 after cleaning (Eqs 14-9..14-13) | +| 5 | Scaled difference (`ma2`) | Collapses to ADS 2 after cleaning; all present cells are unreplicated (Eqs 13-9..13-13) | +| 6 | Scaled difference (`ma2`) | Collapses to ADS 3 after cleaning (Eqs 13-9..13-13) | The scaled difference first removes the process mean and the condition and period effects, Z = Y − Ȳ_k − Ȳ_t + Ȳ (unweighted means of cell means), then differences Z along the diff --git a/docs/reference/validation.md b/docs/reference/validation.md index 5d15859..515f666 100644 --- a/docs/reference/validation.md +++ b/docs/reference/validation.md @@ -70,7 +70,7 @@ effects). It is a report, not a gate. CI runs it on every push and keeps its out - The limits of the two factor S charts (pages 25 and 27) in the known-effects run, whose subgroups have 1000 and 2000 values. Its VAS run of 6 October 2026 held αN - (10-12 manual Eq 10-11) at 0.9995 from N = 343 up, where Minitab's GAMMA function overflows, + (11-7 manual Eq 9-8) at 0.9995 from N = 343 up, where Minitab's GAMMA function overflows, so its limits follow N = 500 instead of the subgroup's own N. Dr. Bishop's reruns of the PM runs (13 October 2026) compute αN for every N, and every R6 S chart in them is checked and matches; the known-effects run waits for the same rerun. diff --git a/docs/tutorials/loss-function.ipynb b/docs/tutorials/loss-function.ipynb index a638f9c..389d088 100644 --- a/docs/tutorials/loss-function.ipynb +++ b/docs/tutorials/loss-function.ipynb @@ -7,7 +7,7 @@ "source": [ "# Taguchi Loss Function Analysis\n", "\n", - "Process behavior charts answer *\"Is this process stable?\"* The Taguchi Loss Function (Bishop Chapter 15) answers a different question: *\"Where is the variation coming from, and how much does each source contribute to expected loss?\"*\n", + "Process behavior charts answer *\"Is this process stable?\"* The Taguchi Loss Function (Bishop Chapter 14) answers a different question: *\"Where is the variation coming from, and how much does each source contribute to expected loss?\"*\n", "\n", "It decomposes the expected squared loss into **5 components**:\n", "\n", @@ -2304,7 +2304,7 @@ "\n", "## DS 2: ADS 2/3 Unexplained and Interaction\n", "\n", - "DS 2 has singleton cells, so there is no within-cell spread to average. The unexplained component is the sample variance of R2 (10-12 manual Eq 15-17), and the interaction is what the other four components leave of the total loss $(\\bar{Y}_{..} - T)^2 + S^2$, set to 0 if negative (Eqs 15-18 to 15-20). With one observation per subgroup, interaction and noise cannot be separated cell by cell, so the interaction is backed out of the total rather than estimated directly:" + "DS 2 has singleton cells, so there is no within-cell spread to average. The unexplained component is the sample variance of R2 (11-7 manual Eq 14-15), and the interaction is what the other four components leave of the total loss $(\\bar{Y}_{..} - T)^2 + S^2$, set to 0 if negative (Eqs 14-16 to 14-18). With one observation per subgroup, interaction and noise cannot be separated cell by cell, so the interaction is backed out of the total rather than estimated directly:" ] }, { @@ -2503,9 +2503,9 @@ "\n", "## The Equations\n", "\n", - "For reference, here are the Bishop Chapter 15 equations implemented by `loss_function()`:\n", + "For reference, here are the Bishop Chapter 14 equations implemented by `loss_function()`:\n", "\n", - "**Expected Loss Decomposition** (Eq. 15-13):\n", + "**Expected Loss Decomposition** (Eq. 14-13):\n", "\n", "$$EL = (\\bar{Y}_{..} - T)^2 + \\hat{\\sigma}^2_w + \\frac{1}{K}\\sum_k \\rho_k^2 + \\frac{1}{T}\\sum_t \\tau_t^2 + \\frac{1}{KT}\\sum_{k,t} (\\rho\\tau)_{kt}^2$$\n", "\n", @@ -2514,18 +2514,18 @@ "| Term | Component | Equation |\n", "|------|-----------|----------|\n", "| $(\\bar{Y}_{..} - T)^2$ | Centering | Squared distance from grand mean to target |\n", - "| $\\hat{\\sigma}^2_w$ | Unexplained | Noise variance (Eq. 15-16 / 15-17) |\n", - "| $\\frac{1}{K}\\sum_k \\rho_k^2$ | PDC | Mean squared PDC effect (Eq. 15-22) |\n", + "| $\\hat{\\sigma}^2_w$ | Unexplained | Noise variance (Eq. 14-14 / 14-15) |\n", + "| $\\frac{1}{K}\\sum_k \\rho_k^2$ | PDC | Mean squared PDC effect (Eq. 14-20) |\n", "| $\\frac{1}{T}\\sum_t \\tau_t^2$ | Time | Mean squared time effect |\n", "| $\\frac{1}{KT}\\sum_{k,t} (\\rho\\tau)_{kt}^2$ | Interaction | Mean squared PDC x time interaction (DS 1/4); remainder (DS 2/3/5/6) |\n", "\n", "Equation numbers are those of Bishop's VAS manual dated 10-12-2026.\n", "\n", - "**Replicated cells** (ADS 1, Eq. 15-16) — the average cell variance:\n", + "**Replicated cells** (ADS 1, Eq. 14-14) — the average cell variance:\n", "\n", "$$\\hat{\\sigma}^2_w = \\frac{1}{KT}\\sum_{k,t} S_{kt}^2$$\n", "\n", - "**Singleton cells** (ADS 2/3) — the variance of R2 (Eq. 15-17), and the interaction as the remainder of the total loss (Eqs. 15-18 to 15-20), floored at 0:\n", + "**Singleton cells** (ADS 2/3) — the variance of R2 (Eq. 14-15), and the interaction as the remainder of the total loss (Eqs. 14-16 to 14-18), floored at 0:\n", "\n", "$$\\hat{\\sigma}^2_w = S_{R2}^2 \\qquad EL_{PDC\\times PT} = \\max\\left(S^2 - \\frac{1}{K}\\sum_k \\rho_k^2 - \\frac{1}{T}\\sum_t \\tau_t^2 - S_{R2}^2,\\ 0\\right)$$\n", "\n", diff --git a/docs/tutorials/process-capability.ipynb b/docs/tutorials/process-capability.ipynb index 608cd6f..cbc07ef 100644 --- a/docs/tutorials/process-capability.ipynb +++ b/docs/tutorials/process-capability.ipynb @@ -9,7 +9,7 @@ "\n", "Process capability answers a different question from process behavior charts. Charts ask *\"Is this process stable?\"* Capability asks *\"Can this process consistently meet specifications?\"*\n", "\n", - "Wheeler & Bishop (Chapter 16) define two families of capability indices:\n", + "Wheeler & Bishop (Chapter 15) define two families of capability indices:\n", "\n", "| Index | Measures | Based on |\n", "|-------|----------|----------|\n", @@ -1703,7 +1703,7 @@ "source": [ "## The Equations\n", "\n", - "For reference, here are the Bishop Chapter 16 equations implemented by `capability()`:\n", + "For reference, here are the Bishop Chapter 15 equations implemented by `capability()`:\n", "\n", "**Current Capability** (overall sigma):\n", "\n", diff --git a/docs/tutorials/sds-validation.ipynb b/docs/tutorials/sds-validation.ipynb index 216e3a5..d90b34e 100644 --- a/docs/tutorials/sds-validation.ipynb +++ b/docs/tutorials/sds-validation.ipynb @@ -73055,7 +73055,7 @@ "---\n", "## PM SDS 3: Full Validation Against Bishop's Reference\n", "\n", - "This section replicates Tom Bishop's VAS analysis of PM SDS 3 — an SDS 3 (Partial Replication) dataset with a mix of n=1 and n≥2 cells — against his VAS run of 3 October 2026 (for the 10-1-2026 manual). In that manual Xbar/S is VAS's default chart for SDS 3, so the original data and the interaction are charted as Xbar/S. We validate each chart against Bishop's reference output.\n", + "This section replicates Tom Bishop's VAS analysis of PM SDS 3 — an SDS 3 (Partial Replication) dataset with a mix of n=1 and n≥2 cells — against his VAS run of 3 October 2026 (for the 11-7 manual). In that manual Xbar/S is VAS's default chart for SDS 3, so the original data and the interaction are charted as Xbar/S. We validate each chart against Bishop's reference output.\n", "\n", "Spec limits: **LSL = 232, USL = 242, Target = 237**" ] diff --git a/docs/tutorials/staged-charts.ipynb b/docs/tutorials/staged-charts.ipynb index c1c51bd..28f05c0 100644 --- a/docs/tutorials/staged-charts.ipynb +++ b/docs/tutorials/staged-charts.ipynb @@ -7291,7 +7291,7 @@ "source": [ "Each stage is computed as a separate analysis. Its center is the unweighted mean of its subgroup\n", "means, its S̄ is the mean of its subgroup standard deviations, and its σ̂ and limits follow from those\n", - "(10-12 manual Eqs 11-5, 11-8, 11-9, 11-16, 11-17). Nothing is pooled across stages, so each stage is\n", + "(11-7 manual Eqs 10-5, 10-8, 10-9, 10-16, 10-17). Nothing is pooled across stages, so each stage is\n", "judged against its own data." ] }, diff --git a/docs/user-guide/chart-types.md b/docs/user-guide/chart-types.md index 3cdf92e..031e1f9 100644 --- a/docs/user-guide/chart-types.md +++ b/docs/user-guide/chart-types.md @@ -266,8 +266,8 @@ staged.phase_statistics('Xbar') # one row per stage: center, S-bar, sigma-hat, Each stage is computed as a separate analysis, exactly as `by=[time]` computes one chart per stage: its center is the unweighted mean of its subgroup means, its S̄ is the mean of its -subgroup standard deviations, and its σ̂ and limits follow from those (10-12 manual Eqs 11-5, -11-8, 11-9, 11-16, 11-17). Nothing is pooled across stages. On Wheeler's four-cavity +subgroup standard deviations, and its σ̂ and limits follow from those (11-7 manual Eqs 10-5, +10-8, 10-9, 10-16, 10-17). Nothing is pooled across stages. On Wheeler's four-cavity thickness data, with the cavity as the stage, the per-cavity centers and limits match Minitab's staged Xbar-S chart to the digit; the per-cavity charts of `by=['Hour']` give the same numbers, one chart at a time. @@ -314,8 +314,8 @@ factor pooled or charted on its own. Plots each individual observation. -- **Centerline**: Average of all observations (X̅) (10-12 manual Eqs 12-7, 12-9) -- **Control Limits**: X̅ ± (3/1.128) × R̅, about X̅ ± 2.66 × R̅ (R̅ is the average of the T − 1 moving ranges; Eqs 12-1, 12-10, 12-11) +- **Centerline**: Average of all observations (X̅) (11-7 manual Eq 11-7) +- **Control Limits**: X̅ ± (3/1.128) × R̅, about X̅ ± 2.66 × R̅ (R̅ is the average of the T − 1 moving ranges; Eqs 11-1, 11-8, 11-9) - **Interpretation**: Points beyond limits indicate special causes ### The mR (Moving Range) Chart @@ -323,7 +323,7 @@ Plots each individual observation. Plots the absolute difference between consecutive observations. - **Centerline**: Average moving range (R̅) -- **UCL**: (1 + 3 × 0.8525/1.128) × R̅, about 3.267 × R̅ (10-12 manual Eq 12-5) +- **UCL**: (1 + 3 × 0.8525/1.128) × R̅, about 3.267 × R̅ (11-7 manual Eq 11-5) - **LCL**: 0 (range cannot be negative) - **Interpretation**: Large ranges indicate sudden changes @@ -443,7 +443,7 @@ result = study.execute(chart='X', by=['lane'], value='R4') result = study.execute(chart='X', by=['lane'], value='R4', recentered=True) ``` -Re-centering adds back the grand mean, RCR = R + Y̅ (Bishop's 10-12 manual, Eq 14-26, for R1 +Re-centering adds back the grand mean, RCR = R + Y̅ (Bishop's 11-7 manual, Eq 13-26, for R1 and R3–R6), so a re-centered chart plots the effect plus R2 on the measurement scale. The chart shifts by Y̅; the limit width is unchanged. RCR2 is the exception: Y̅_kt + R2. Chart period effects from R4; an Xbar of RCR3 by period is flat, because the interaction sums to diff --git a/docs/user-guide/formulation.md b/docs/user-guide/formulation.md index 8b8cb13..7a83b17 100644 --- a/docs/user-guide/formulation.md +++ b/docs/user-guide/formulation.md @@ -417,7 +417,7 @@ This tells you: the plan called for full replication across 480 cells. Three cel ### `study.capability()` -Assesses process capability against specification limits (Bishop Ch. 16). Returns both **current capability** (Pp/Ppk, based on overall variation) and **potential capability** (Cp/Cpk, based on the R2 noise only). +Assesses process capability against specification limits (Bishop Ch. 15). Returns both **current capability** (Pp/Ppk, based on overall variation) and **potential capability** (Cp/Cpk, based on the R2 noise only). ```python from processbehavior import SpecLimits @@ -438,7 +438,7 @@ cap.plot(values=study.dataset[study.response].dropna().values) | Metric | Based On | Meaning | |--------|----------|---------| | **Pp/Ppk** | Overall σ̂ (all variation) | Current capability — process as-is | -| **Cp/Cpk** | R2 σ̂ = SR2 / c4(NR2) (noise only; 10-12 manual Eq 16-10) | Potential capability — achievable if all assignable causes are eliminated | +| **Cp/Cpk** | R2 σ̂ = SR2 / c4(NR2) (noise only; 11-7 manual Eq 15-10) | Potential capability — achievable if all assignable causes are eliminated | The potential values are ȳ + R2. Potential CPL/CPU are measured from their mean, `cap.potential_center` (ȳ + mean of R2), and the potential chart's center line sits there. @@ -449,7 +449,7 @@ Potential capability (Cp/Cpk) requires VAS residuals (factors + time). If unavai ### `study.loss_function()` -Decomposes expected loss into five components using the Taguchi Loss Function (Bishop Ch. 15; equation numbers below are from Bishop's VAS documentation manual dated 10-12-2026). Identifies the largest sources of variation as a Pareto analysis. It needs R2, so it raises `ValidationError` when R2 is unavailable for the layout. +Decomposes expected loss into five components using the Taguchi Loss Function (Bishop Ch. 14; equation numbers below are from Bishop's VAS documentation manual dated 11-7-2026). Identifies the largest sources of variation as a Pareto analysis. It needs R2, so it raises `ValidationError` when R2 is unavailable for the layout. ```python # Use grand mean as target (centering = 0) @@ -478,10 +478,10 @@ loss.plot(structured=True) # Expands PDC into per-factor components | Component | Formula | What It Captures | |-----------|---------|-----------------| | **Centering** | (Ȳ - Target)² | Loss from being off-target | -| **Unexplained** | ADS 1: average of the cell variances, (1/KT)·ΣS²kt (Eq 15-16). ADS 2/3: S²R2, the sample variance of R2 (Eq 15-17) | Irreducible noise (R2) | +| **Unexplained** | ADS 1: average of the cell variances, (1/KT)·ΣS²kt (Eq 14-14). ADS 2/3: S²R2, the sample variance of R2 (Eq 14-15) | Irreducible noise (R2) | | **PDC** | Between-condition variance | Process design condition main effects (R5) | | **Time** | Between-time variance | Time period effects | -| **Interaction** | ADS 1: average squared interaction effect (Eq 15-16). ADS 2/3: what the other parts leave of the total, S² − PDC − Time − Unexplained, set to 0 if negative (Eqs 15-18..15-20) | How factor effects change over time | +| **Interaction** | ADS 1: average squared interaction effect (Eq 14-14). ADS 2/3: what the other parts leave of the total, S² − PDC − Time − Unexplained, set to 0 if negative (Eqs 14-16..14-18) | How factor effects change over time | S² is the sample variance of all observations. The total is the sum of the five components; in ADS 2 and 3 it equals (Ȳ - Target)² + S², unless the interaction remainder was negative diff --git a/docs/user-guide/residuals.md b/docs/user-guide/residuals.md index 9019438..271a2b9 100644 --- a/docs/user-guide/residuals.md +++ b/docs/user-guide/residuals.md @@ -15,10 +15,10 @@ Dr. Thomas A. Bishop's **Variation Analysis System (VAS)** decomposes total vari Where: - Y = individual observation -- Y̅ = grand mean: the mean of the cell means, each cell counted once (10-12 manual Eq 10-9) +- Y̅ = grand mean: the mean of the cell means, each cell counted once (11-7 manual Figure 13.3, the estimates of the VAS model parameters) - Y̅k = mean for factor level k: the mean of that level's cell means - Y̅t = mean at time t: the mean of that period's cell means -- Y̅kt = mean for factor k at time t (cell mean, Eq 10-7) +- Y̅kt = mean for factor k at time t (cell mean, Eq 9-6) ## Accessing Residuals @@ -44,7 +44,7 @@ print(result.residuals.head()) **Purpose**: View the original measurements as ± about zero, so the total range of variation is read directly rather than against an arbitrary process mean. -**Formula**: R1 = Y - Y̅ (Eq 14-2) +**Formula**: R1 = Y - Y̅ (Eq 13-2) R1 is a pure location shift — it subtracts the overall mean and changes nothing else. It is also the building block the other residuals are derived from, which makes it the @@ -83,10 +83,10 @@ exactly in DS 1, approximately in DS 2 and 3, where differencing the adjusted se interaction that changes smoothly over time). **Formula by DS**: -- **DS 1 (Full Replication)**: R2 = Y - Y̅kt, the deviation from the cell mean (Eq 14-3) +- **DS 1 (Full Replication)**: R2 = Y - Y̅kt, the deviation from the cell mean (Eq 13-3) - **DS 2 (No Replication)** and **DS 3 (Partial)**: with any one-observation subgroup there is no within-cell deviation to use, so R2 is built from the whole series in three steps - (Eqs 14-4..14-8 for DS 2, 14-9..14-13 for DS 3, 5 and 6): + (Eqs 13-4..13-8 for DS 2, 14-9..14-13 for DS 3, 5 and 6): 1. Remove the process mean and the condition and period effects: Z = Y - Y̅k - Y̅t + Y̅ (Y̅k, Y̅t and Y̅ are unweighted means of cell means). @@ -121,7 +121,7 @@ fig = result.plot(show_zones=True, title='Unexplained Variation (R2)') **Purpose**: Detect factor × time interactions. -**Formula**: R3 = (Y̅kt - Y̅k - Y̅t + Y̅) + R2 (Eq 14-14) +**Formula**: R3 = (Y̅kt - Y̅k - Y̅t + Y̅) + R2 (Eq 13-14) This removes both main effects, leaving the interaction plus the unexplained noise. In DS 1, where R2 = Y - Y̅kt, it reduces to Y - Y̅k - Y̅t + Y̅. @@ -287,9 +287,9 @@ result = study.execute(chart='X', by=['lane'], value='R4', recentered=True) ``` Re-centering formulas: -- RCR = R + Y̅ for R1, R3, R4, R5 and R6 (Eq 14-26) — each adds back the grand mean only, so a +- RCR = R + Y̅ for R1, R3, R4, R5 and R6 (Eq 13-26) — each adds back the grand mean only, so a re-centered chart is the zero-centered chart shifted onto the measurement scale -- RCR2 = Y̅kt + R2 — not part of Eq 14-26; in DS 1 this reconstructs Y +- RCR2 = Y̅kt + R2 — not part of Eq 13-26; in DS 1 this reconstructs Y So an effect chart of RCR5 by condition, RCR4 by period, or RCR3 by subgroup plots that effect plus R2 on the measurement scale. Chart period effects from R4: an Xbar of RCR3 by @@ -314,10 +314,10 @@ R6 requires factors (it is computed from R5 and R2 at `execute()` time). **Note on R2 calculation**: R2 is the one residual whose formula depends on structure, and the choice is made on the analytical design state (ADS), after cleaning: -- **ADS 1** (every cell n ≥ 2): within-cell deviation, `R2 = Y - Ȳ_kt` (Eq 14-3) +- **ADS 1** (every cell n ≥ 2): within-cell deviation, `R2 = Y - Ȳ_kt` (Eq 13-3) - **ADS 2 and 3** (any one-observation cell): the condition-and-period-adjusted series Z, differenced along the full canonical sequence and divided by twice the R2 scale factor, - for every observation (Eqs 14-4..14-13); no per-cell mixing. See [R2](#r2-unexplained-residuals). + for every observation (Eqs 13-4..13-13); no per-cell mixing. See [R2](#r2-unexplained-residuals). DS 4, 5 and 6 are sampling states; they collapse to ADS 1, 2 and 3 once empty cells are dropped. @@ -429,7 +429,7 @@ has the same number of observations. R2 = Y - Ȳ_kt ``` -Where Y̅kt is the mean of observations in cell (k, t) (Eq 14-3). +Where Y̅kt is the mean of observations in cell (k, t) (Eq 13-3). ### R2: Scaled Difference (DS 2 and 3, any one-observation cell) @@ -440,7 +440,7 @@ c(K, M) = sqrt( (K - 1) / (2K) · (1 - K / (M - 1)) ) ``` K is the number of process design conditions present and M the number of observations -(M = KT in DS 2); c(K, M) is the R2 scale factor (Eqs 14-4..14-8 for DS 2, 14-9..14-13 for +(M = KT in DS 2); c(K, M) is the R2 scale factor (Eqs 13-4..13-8 for DS 2, 14-9..14-13 for DS 3, 5 and 6). The first observation in the sequence has no value (NaN). The scale factor is undefined, and R2 unavailable, when K < 2 or M < K + 2. @@ -451,7 +451,7 @@ R3 = (Ȳ_kt - Ȳ_k - Ȳ_t + Ȳ) + R2 = InteractionEffect + R2 ``` -(Eq 14-14.) In DS 1 this equals Y - Ȳ_k - Ȳ_t + Ȳ, since Y = Ȳ_kt + R2 there. +(Eq 13-14.) In DS 1 this equals Y - Ȳ_k - Ȳ_t + Ȳ, since Y = Ȳ_kt + R2 there. ### R4: Time Main Effect + Unexplained diff --git a/docs/user-guide/sds-detection.md b/docs/user-guide/sds-detection.md index e77e65c..f3e026c 100644 --- a/docs/user-guide/sds-detection.md +++ b/docs/user-guide/sds-detection.md @@ -265,7 +265,7 @@ Printed, the report adds a line for each of these when it applies: ``` **Capabilities**: -- ⚠️ No within-cell variance: R2 is the condition-and-period-adjusted series, differenced along the condition-then-time sequence and divided by twice the R2 scale factor (Eqs 14-4..14-8) +- ⚠️ No within-cell variance: R2 is the condition-and-period-adjusted series, differenced along the condition-then-time sequence and divided by twice the R2 scale factor (Eqs 13-4..13-8) - ⚠️ R2 needs at least 2 process design conditions (K) and M ≥ K + 2 observations; otherwise it is unavailable - ✅ Main effects analysis - ✅ Xbar-S of the effect residuals (R4, R5), with limits from R2 @@ -293,7 +293,7 @@ Printed, the report adds a line for each of these when it applies: **Valid Charts**: Histogram, Xbar, S, X, mR :::{note} Why DS 3 does not mix methods -For R2 in DS 3 the library does not combine an exact estimate from the replicated cells with a separate estimate for the one-observation cells. Bishop's R2 for any layout with a one-observation cell (Eqs 14-9..14-13) is defined over the whole ordered sequence, with no grouping by cell: remove the process mean and the condition and period effects, difference what is left along the canonical condition-then-time sort, and divide by twice the R2 scale factor c(K, M). Only the first observation in the sequence has no value. The replicated cells' within-cell deviations are not used for R2. This is the same calculation as DS 2, and the ADS 3 reference assertions in the [validation page](../reference/validation.md) hold against it. The recommended chart is Xbar (with S), VAS's default for DS 3 and 6. Every subgroup mean, one-observation subgroups included, counts once toward the Xbar center line (per stratum too). One-observation subgroups have no within-subgroup spread, so they take no part in σ̂ (Eq 11-5); they are still charted, with limits center ± 3σ̂ (Eqs 11-16/17 with N = 1), as VAS draws them. The companion S chart leaves them out. +For R2 in DS 3 the library does not combine an exact estimate from the replicated cells with a separate estimate for the one-observation cells. Bishop's R2 for any layout with a one-observation cell (Eqs 13-9..13-13) is defined over the whole ordered sequence, with no grouping by cell: remove the process mean and the condition and period effects, difference what is left along the canonical condition-then-time sort, and divide by twice the R2 scale factor c(K, M). Only the first observation in the sequence has no value. The replicated cells' within-cell deviations are not used for R2. This is the same calculation as DS 2, and the ADS 3 reference assertions in the [validation page](../reference/validation.md) hold against it. The recommended chart is Xbar (with S), VAS's default for DS 3 and 6. Every subgroup mean, one-observation subgroups included, counts once toward the Xbar center line (per stratum too). One-observation subgroups have no within-subgroup spread, so they take no part in σ̂ (Eq 10-5); they are still charted, with limits center ± 3σ̂ (Eqs 10-16/17 with N = 1), as VAS draws them. The companion S chart leaves them out. ::: ## DS 4: Incomplete, No Singletons @@ -386,9 +386,9 @@ The R2 method is determined by the tidy data structure (ADS), not the raw DS. Eq | ADS | R2 Method | Description | |-----|-----------|-------------| -| 1 | exact | Within-cell deviation (`R2 = Y - Ȳ_kt`, Eq 14-3) | -| 2 | ma2 | `Z = Y - Ȳ_k - Ȳ_t + Ȳ`, differenced along the condition-then-time sequence and divided by twice the R2 scale factor c(K, M) (Eqs 14-4..14-8) | -| 3 | ma2 | Any one-observation cell: the same calculation over every observation (Eqs 14-9..14-13) | +| 1 | exact | Within-cell deviation (`R2 = Y - Ȳ_kt`, Eq 13-3) | +| 2 | ma2 | `Z = Y - Ȳ_k - Ȳ_t + Ȳ`, differenced along the condition-then-time sequence and divided by twice the R2 scale factor c(K, M) (Eqs 13-4..13-8) | +| 3 | ma2 | Any one-observation cell: the same calculation over every observation (Eqs 13-9..13-13) | The token `ma2` names Bishop's size-2 moving average of Z; R2 is Z minus that average, divided by the R2 scale factor. The factor is computed from the layout, with K the number of diff --git a/docs/user-guide/series-length.md b/docs/user-guide/series-length.md index 90def48..55f82ec 100644 --- a/docs/user-guide/series-length.md +++ b/docs/user-guide/series-length.md @@ -57,7 +57,7 @@ The same numbers appear in code as `study.series_length`, a frozen record with ` ## How precise a moving-range sigma is -For a stable normal process, the moving-range estimate of sigma (MR̄/d₂, 10-12 manual Eq 12-4) scatters around +For a stable normal process, the moving-range estimate of sigma (MR̄/d₂, 11-7 manual Eq 11-4) scatters around the truth by an amount that depends only on how many points there are. The table gives the 10th, 50th and 90th percentiles of the estimate relative to the true sigma, and the ratio of the 90th to the 10th, for every series length from 3 to 30. It was generated by diff --git a/processbehavior/analysis.py b/processbehavior/analysis.py index cb68c48..5051cf3 100644 --- a/processbehavior/analysis.py +++ b/processbehavior/analysis.py @@ -254,7 +254,7 @@ def _mr_family_limits(mr_spec: _MRChartSpec, *, mean, mR) -> pd.DataFrame: """Limits for one MR-family chart per row of a per-stratum or per-phase table. The array form of the scalar ``calculate_limits`` these sites used to call through - ``DataFrame.apply(axis=1)``: the same formulae, 10-12 manual Eqs 12-10/12-11 (X) and + ``DataFrame.apply(axis=1)``: the same formulae, 11-7 manual Eqs 11-8/11-9 (X) and 12-5 (mR). The scalar call built a ``pd.Series`` per stratum, about 0.1 ms each, so it cost 30s per chart at 300,000 strata. ``mean`` is ignored for the mR chart. """ @@ -921,7 +921,7 @@ def _residual_grain(self, value_col: str) -> list[str]: """Return the cell-grid columns at which to compute Bishop's grand mean for an Xbar chart on this column. - Each re-centred residual is R + Ybar (10-12 manual Eq 14-26; see + Each re-centred residual is R + Ybar (11-7 manual Eq 13-26; see `AnalysisDataSet._calculate_centered_residuals` and `residual_calculator.calculate_r6_residuals`). The canonical grand mean of an RCRk chart averages over the grain its effect lives at — equal weight per @@ -1216,10 +1216,10 @@ def _build_output_columns(self, df: pd.DataFrame, value_cols: list[str]) -> pd.D def _xbar_limits(self, out: pd.DataFrame, s_bar: float, n_to_use: str) -> pd.DataFrame: """Xbar natural process limits for the subgroups in ``out`` (columns center, s, n, N). - Equal subgroup sizes (or n_mode="average", Eq 11-20): S-bar / c4(N) for every row. - Unequal sizes: sigma_hat is the average of each subgroup's S_r / c4(N_r) (10-12 manual - Eq 11-5) and each subgroup's limits are center ± multiplier · sigma_hat / sqrt(N_r) - (Eqs 11-16, 11-17), so they step with N_r, as VAS draws them. A one-reading subgroup + Equal subgroup sizes (or n_mode="average", Eq 10-20): S-bar / c4(N) for every row. + Unequal sizes: sigma_hat is the average of each subgroup's S_r / c4(N_r) (11-7 manual + Eq 10-5) and each subgroup's limits are center ± multiplier · sigma_hat / sqrt(N_r) + (Eqs 10-16, 10-17), so they step with N_r, as VAS draws them. A one-reading subgroup has no S_r: it takes no part in sigma_hat and is charted at center ± multiplier · sigma_hat. """ if n_to_use == 'n': @@ -1283,12 +1283,12 @@ def _fit_xbar_group(self, out: pd.DataFrame) -> _GroupFit | None: Every subgroup with a value is charted. A one-reading subgroup has no S, so sigma-hat comes from the subgroups with n >= 2, and its limits are center ± 3·sigma-hat (10-12 - manual Eqs 11-16, 11-17 with N_r = 1), as VAS draws them; the companion S chart keeps + manual Eqs 10-16, 10-17 with N_r = 1), as VAS draws them; the companion S chart keeps only n >= 2. n counts the limits column, so a reading without one (R2's first value) is not charted, as in VAS. Returns None when no subgroup has n >= 2; the caller decides whether that raises (one chart) or is recorded (a stratum, a phase). Under ``n_mode='average'`` the ``N`` column becomes the mean fitted n and both charts take - the constant-N path (Eq 11-20). + the constant-N path (Eq 10-20). """ out = out[out['n'].ge(1)].copy() # n=0: every charted value missing fit = out['n'].ge(2) @@ -2007,9 +2007,9 @@ def _calculate_xbar_s_phased(self, value_col: str | None = None) -> dict[str, Ch A phase is one combination of the grouping factors: a contiguous run of subgroups in time order on the stage-by-stage layout that ``by=None`` already draws. Each phase is computed as the per-condition chart of ``by=[time]`` computes one stratum — the - unweighted mean of its subgroup means as centre (10-12 manual Eq 10-9 within the - phase), its own S-bar and sigma-hat (Eq 11-5), its own limits (Eqs 11-8, 11-9, 11-16, - 11-17) — and the phases are concatenated in display order with the stepped columns. + unweighted mean of its subgroup means as centre (11-7 manual Figure 13.3 within the + phase), its own S-bar and sigma-hat (Eq 10-5), its own limits (Eqs 10-8, 10-9, 10-16, + 10-17) — and the phases are concatenated in display order with the stepped columns. Nothing is pooled across phases: Dr. Bishop (2026-10-18), "treat each stage separately as if it were a single analysis". @@ -2177,10 +2177,10 @@ def _calculate_mr_chart( The X and mR charts share >85% of their calculation logic. The behavioral differences are encoded in ``mr_spec`` — no boolean flags needed. The - statistics are the 10-12 manual's Chapter 12: mR-bar over the T - 1 - consecutive ranges (Eq 12-1), X centre Y-bar (Eqs 12-7, 12-9), X limits - Y-bar ± 3·mR-bar/d2 (Eqs 12-10, 12-11) and mR limits 0 … D4·mR-bar - (Eq 12-5), within each stratum or phase when the chart is split. + statistics are the 11-7 manual's Chapter 11: mR-bar over the T - 1 + consecutive ranges (Eq 11-1), X centre Y-bar (Eq 11-7), X limits + Y-bar ± 3·mR-bar/d2 (Eqs 11-8, 11-9) and mR limits 0 … D4·mR-bar + (Eq 11-5), within each stratum or phase when the chart is split. Parameters ---------- @@ -2895,11 +2895,11 @@ def _calculate_mr_chart_ungrouped( # noqa: C901 return result # === Global limits path (unchanged) === - # Moving range: mR-bar over the T - 1 consecutive ranges (10-12 manual Eq 12-1). A re-centred + # Moving range: mR-bar over the T - 1 consecutive ranges (11-7 manual Eq 11-1). A re-centred # R2 chart takes them over R2 (_resolve_mr_source_column), not the plotted RCR2. out['mr'] = out[mr_source_col].diff().abs() mR = out['mr'].mean() - mean_ = out[value_col].mean() # Y-bar, the X chart's centre (Eqs 12-7, 12-9) + mean_ = out[value_col].mean() # Y-bar, the X chart's centre (Eq 11-7) # R chart: drop first observation if mr_spec.drops_first_mr: diff --git a/processbehavior/analysis_dataset.py b/processbehavior/analysis_dataset.py index 6e984b2..0e212b0 100644 --- a/processbehavior/analysis_dataset.py +++ b/processbehavior/analysis_dataset.py @@ -172,7 +172,7 @@ def r2_unavailable_reason(self) -> str | None: def _check_r2_scale_factor(self, r2_method: str) -> None: """Record (and warn once) when the layout leaves the R2 scale factor undefined. - The scaled-difference R2 (10-12 manual Eqs 14-4..14-13) needs at least two process + The scaled-difference R2 (11-7 manual Eqs 13-4..13-13) needs at least two process design conditions and more than K + 1 observations; outside that it has no value. """ if r2_method != 'ma2': @@ -244,7 +244,7 @@ def _compute_structure_stats(self) -> None: # Compute n_per_cell and cell means once self._n_per_cell = df.groupby('cell_key', observed=True)[y].transform('size') - self._ybar_kt = df.groupby('cell_key', observed=True)[y].transform('mean') # Ȳ_kt (10-12 manual Eq 10-7) + self._ybar_kt = df.groupby('cell_key', observed=True)[y].transform('mean') # Ȳ_kt (11-7 manual Eq 9-6) # Handle edge case: empty data or all NaN if len(self._n_per_cell) == 0 or self._n_per_cell.isna().all(): @@ -356,12 +356,12 @@ def _calculate_centered_residuals(self): Calculate the re-centred residuals (RCR) and the R1 cell/marginal means (Rbar). Re-centring puts a residual back on the measurement scale without changing what it - shows: RCR = R + Ȳ.. (Bishop's 10-12 manual Eq 14-26, for R1, R3, R4, R5 and R6; RCR6 + shows: RCR = R + Ȳ.. (Bishop's 11-7 manual Eq 13-26, for R1, R3, R4, R5 and R6; RCR6 is built per request in residual_calculator.calculate_r6_residuals). An effect chart of RCR5 by condition, RCR4 by period or RCR3 by subgroup therefore plots that effect plus R2, centred on the grand mean. - RCR2 is not part of Eq 14-26; it is kept as Ȳ_kt + R2 (which reconstructs Y in ADS 1). + RCR2 is not part of Eq 13-26; it is kept as Ȳ_kt + R2 (which reconstructs Y in ADS 1). Note: This method intentionally mutates self.analysis_dataset in-place (adding columns directly to the DataFrame). This differs from @@ -382,9 +382,9 @@ def _calculate_centered_residuals(self): df['Rbar_k'] = df.groupby([self.spec.rsg_var_name], observed=True)['R1'].transform('mean') df['Rbar_t'] = df.groupby([self.spec.time_var], observed=True)['R1'].transform('mean') - # Re-centred residuals, Eq 14-26: RCR = R + Ybar + # Re-centred residuals, Eq 13-26: RCR = R + Ybar df['RCR1'] = df['Ybar'] + df['R1'] - df['RCR2'] = df['Ybar_kt'] + df['R2'] # not in Eq 14-26; reconstructs Y in ADS 1 + df['RCR2'] = df['Ybar_kt'] + df['R2'] # not in Eq 13-26; reconstructs Y in ADS 1 df['RCR3'] = df['Ybar'] + df['R3'] df['RCR4'] = df['Ybar'] + df['R4'] df['RCR5'] = df['Ybar'] + df['R5'] diff --git a/processbehavior/capability.py b/processbehavior/capability.py index a3b13ce..53017ca 100644 --- a/processbehavior/capability.py +++ b/processbehavior/capability.py @@ -1,5 +1,5 @@ """ -Process Capability Analysis — Bishop Chapter 16. +Process Capability Analysis — Bishop Chapter 15 (11-7 manual). Computes capability indices given specification limits. Core logic lives as pure module-level functions (composition pattern matching ResidualCalculator / @@ -19,7 +19,7 @@ References ---------- -Bishop, T. *Variation Analysis System* — Chapter 16. +Bishop, T. *Variation Analysis System* — Chapter 15 (11-7 manual). """ from __future__ import annotations @@ -89,7 +89,7 @@ def is_two_sided(self) -> bool: @dataclass(frozen=True) class CapabilityResult: """ - Result of a process capability assessment (Bishop Ch. 16). + Result of a process capability assessment (Bishop Ch. 15). Raw ``float`` values are stored unrounded. Use :meth:`as_dict` for a rounded dictionary or ``repr()`` for a formatted summary. @@ -111,7 +111,7 @@ class CapabilityResult: pp, ppk_lower, ppk_upper, ppk : float | None Current capability indices (from overall sigma). sigma_hat_r2 : float | None - Unbiased sigma from R2 residuals: S_R2 / c4(N_R2) (10-12 manual Eq 16-10). + Unbiased sigma from R2 residuals: S_R2 / c4(N_R2) (11-7 manual Eq 15-10). cp, cpk_lower, cpk_upper, cpk : float | None Potential capability indices (from R2 sigma), with CPL/CPU measured from ``potential_center``. @@ -403,9 +403,9 @@ def compute_sigma_hat(values: np.ndarray) -> tuple[float, float]: Notes ----- - Bishop Ch. 16 — ddof=1 is mandatory. The c4 correction + Bishop Ch. 15 — ddof=1 is mandatory. The c4 correction removes the small-sample bias inherent in the sample standard - deviation: sigma_hat = S / alpha_N (10-12 manual Eq 10-12). + deviation: sigma_hat = S / alpha_N (11-7 manual Eq 9-9). """ n = len(values) s = float(np.std(values, ddof=1)) @@ -560,8 +560,8 @@ def _process_mean(frame: pd.DataFrame, response_var: str) -> float: Each experimental condition counts once, whatever its cell size, as for the Xbar centre line and the loss function's centring term. VAS prints this as the capability charts' PROCESS MEAN (PM SDS 6: 237.86, where the plain average of - the readings is 237.834), and the 10-12 manual's one-sided indices (Eqs 16-5, 16-6, - 16-12, 16-13) centre on it; on balanced data the two are equal. + the readings is 237.834), and the 11-7 manual's one-sided indices (Eqs 15-5, 15-6, + 15-12, 15-13) centre on it; on balanced data the two are equal. A study without cells (no factors or time) has one cell: the plain average. """ valid = frame[frame[response_var].notna()] @@ -714,9 +714,9 @@ def assess_capability( n_r2 = len(r2_values) if n_r2 >= 2: - _, sigma_hat_r2 = compute_sigma_hat(r2_values) # Eq 16-10 + _, sigma_hat_r2 = compute_sigma_hat(r2_values) # Eq 15-10 potential_values = y_bar + r2_values - # Centre the potential indices on the potential values themselves (Eqs 16-11..16-13 as + # Centre the potential indices on the potential values themselves (Eqs 15-11..15-13 as # VAS computes them): y_bar + mean(R2). Tom's Medicare run of 10/3/2026 prints # "PROCESS MEAN 10831.3" on the potential slide (y_bar 10832.4, mean R2 -1.15). potential_center = float(np.mean(potential_values)) diff --git a/processbehavior/loss_function.py b/processbehavior/loss_function.py index 8835cff..3e45340 100644 --- a/processbehavior/loss_function.py +++ b/processbehavior/loss_function.py @@ -1,5 +1,5 @@ """ -Taguchi Loss Function Analysis — Bishop Chapter 15. +Taguchi Loss Function Analysis — Bishop Chapter 14 (11-7 manual). Decomposes expected loss into 5 components that identify where variation comes from: centering, unexplained (within-cell), process design conditions @@ -9,7 +9,7 @@ References ---------- -Bishop, T. *Variation Analysis System* — Chapter 15. +Bishop, T. *Variation Analysis System* — Chapter 14 (11-7 manual). """ from __future__ import annotations @@ -38,7 +38,7 @@ @dataclass(frozen=True) class LossResult: """ - Result of a Taguchi Loss Function analysis (Bishop Ch. 15). + Result of a Taguchi Loss Function analysis (Bishop Ch. 14). Decomposes expected loss into 5 components: centering, unexplained, PDC, time, and PDC×time interaction. When multiple factors exist, @@ -62,8 +62,8 @@ class LossResult: Analytical SDS. centering, unexplained, pdc, time, interaction : float The 5 loss components (squared-loss units). In ADS 1 each is estimated - directly (10-12 manual Eq 15-16). In ADS 2/3 the interaction is what the - other four leave of the total loss (Ȳ.. - T)² + S² (Eqs 15-18 to 15-20), + directly (11-7 manual Eq 14-14). In ADS 2/3 the interaction is what the + other four leave of the total loss (Ȳ.. - T)² + S² (Eqs 14-16 to 14-18), set to 0 when that remainder is negative. total : float Sum of 5 components. In ADS 2/3 this equals (Ȳ.. - T)² + S² unless the @@ -84,7 +84,7 @@ class LossResult: independent readings with a common variance, estimated by the unexplained part, and is exact for any pattern of occupied cells and cell sizes; with every cell holding n readings it is σ²(K-1)/(KTn), σ²(T-1)/(KTn) and σ²(K-1)(T-1)/(KTn). In ADS 2/3 the interaction is - what the other parts leave (10-12 manual Eq 15-20), so its reference is None. Not valid + what the other parts leave (11-7 manual Eq 14-18), so its reference is None. Not valid when readings are autocorrelated. """ @@ -335,13 +335,13 @@ def __repr__(self) -> str: def _compute_centering(y_bar: float, target: float) -> float: - """10-12 manual Eq 15-16 / 15-21, first term: (Ȳ.. - T)².""" + """11-7 manual Eq 14-14 / 14-19, first term: (Ȳ.. - T)².""" return (y_bar - target) ** 2 def _compute_unexplained_replicated(df, response_var: str) -> float: """ - 10-12 manual Eq 15-14/15-16: within-cell variance for fully replicated designs (ADS 1). + 11-7 manual Eq 14-14: within-cell variance for fully replicated designs (ADS 1). unexplained = (1/KT) Σ S_kt² @@ -354,7 +354,7 @@ def _compute_unexplained_replicated(df, response_var: str) -> float: def _compute_unexplained_pooled(df) -> float: """ - 10-12 manual Eq 15-17: unexplained loss for ADS 2/3 is the sample variance of R2. + 11-7 manual Eq 14-15: unexplained loss for ADS 2/3 is the sample variance of R2. unexplained = S²_R2 (n - 1 divisor). R2 is already on the noise scale (the R2 scale factor puts it there), so no further correction is applied. @@ -364,7 +364,7 @@ def _compute_unexplained_pooled(df) -> float: def _compute_pdc(df, rsg_var_name: str) -> float: - """10-12 manual Eq 15-22: (1/K) Σ (Ȳ_k. - Ȳ..)².""" + """11-7 manual Eq 14-20: (1/K) Σ (Ȳ_k. - Ȳ..)².""" ybar = df['Ybar'].iloc[0] ybar_k = df.groupby(rsg_var_name, observed=True)['Ybar_k'].first() rho_k = ybar_k - ybar @@ -372,7 +372,7 @@ def _compute_pdc(df, rsg_var_name: str) -> float: def _compute_time(df, time_var: str) -> float: - """10-12 manual Eq 15-16 / 15-21, time term: (1/T) Σ (Ȳ_.t - Ȳ..)².""" + """11-7 manual Eq 14-14 / 14-19, time term: (1/T) Σ (Ȳ_.t - Ȳ..)².""" ybar = df['Ybar'].iloc[0] ybar_t = df.groupby(time_var, observed=True)['Ybar_t'].first() tau_t = ybar_t - ybar @@ -380,7 +380,7 @@ def _compute_time(df, time_var: str) -> float: def _compute_interaction(df, time_var: str) -> float: - """10-12 manual Eq 15-16, last term (ADS 1): (1/KT) Σ (Ȳ_kt - Ȳ_k. - Ȳ_.t + Ȳ..)².""" + """11-7 manual Eq 14-14, last term (ADS 1): (1/KT) Σ (Ȳ_kt - Ȳ_k. - Ȳ_.t + Ȳ..)².""" ybar = df['Ybar'].iloc[0] cell = df.groupby('cell_key', observed=True)[['Ybar_kt', 'Ybar_k', 'Ybar_t']].first() rho_tau = cell['Ybar_kt'] - cell['Ybar_k'] - cell['Ybar_t'] + ybar @@ -389,18 +389,18 @@ def _compute_interaction(df, time_var: str) -> float: def _compute_interaction_remainder(df, response_var: str, pdc: float, time_loss: float, unexplained: float) -> float: """ - 10-12 manual Eqs 15-18 to 15-20 (ADS 2/3): interaction is what the other parts leave of the total. + 11-7 manual Eqs 14-16 to 14-18 (ADS 2/3): interaction is what the other parts leave of the total. The total loss is (Ȳ.. - T)² + S², with S² the sample variance of all - observations (Eq 15-18). Subtracting the mean, condition, time and unexplained - parts leaves the PDC×PT interaction (Eq 15-20), set to 0 when negative. With + observations (Eq 14-16). Subtracting the mean, condition, time and unexplained + parts leaves the PDC×PT interaction (Eq 14-18), set to 0 when negative. With one observation per subgroup, interaction and noise cannot be told apart cell by cell, so the interaction is backed out rather than estimated directly. """ s2 = float(np.var(df[response_var].to_numpy(dtype=float), ddof=1)) remainder = s2 - pdc - time_loss - unexplained if remainder < 0.0: - logger.debug('Eq 15-20 interaction remainder %.6g is negative; set to 0', remainder) + logger.debug('Eq 14-18 interaction remainder %.6g is negative; set to 0', remainder) return 0.0 return remainder @@ -413,7 +413,7 @@ def _no_effect_reference( Each cell mean is μ + e_kt with Var(e_kt) = σ²/n_kt, independent across cells. Ȳ.., Ȳ_k. and Ȳ_.t are unweighted means of the occupied cells' means, so the expected square of each - deviation in Eq 15-16 follows from these variances and covariances. + deviation in Eq 14-14 follows from these variances and covariances. """ cells = df.groupby('cell_key', observed=True).agg( k=(rsg_var_name, 'first'), t=(time_var, 'first'), n=(rsg_var_name, 'size') @@ -496,7 +496,7 @@ def assess_loss( round_to: int = 3, ) -> LossResult: """ - Assess Taguchi Loss Function decomposition (Bishop Ch. 15). + Assess Taguchi Loss Function decomposition (Bishop Ch. 14). Parameters ---------- @@ -555,14 +555,14 @@ def assess_loss( pdc = _compute_pdc(df, rsg_var_name) time_loss = _compute_time(df, time_var) - # Unexplained and interaction: branch on cell replication (10-12 manual §15.4) + # Unexplained and interaction: branch on cell replication (11-7 manual §15.4) replicated = ads._structure_stats.n_cell_min >= 2 if replicated: - unexplained = _compute_unexplained_replicated(df, response_var) # Eq 15-16 - interaction = _compute_interaction(df, time_var) # Eq 15-16 + unexplained = _compute_unexplained_replicated(df, response_var) # Eq 14-14 + interaction = _compute_interaction(df, time_var) # Eq 14-14 else: - unexplained = _compute_unexplained_pooled(df) # Eq 15-17 - interaction = _compute_interaction_remainder(df, response_var, pdc, time_loss, unexplained) # Eq 15-20 + unexplained = _compute_unexplained_pooled(df) # Eq 14-15 + interaction = _compute_interaction_remainder(df, response_var, pdc, time_loss, unexplained) # Eq 14-18 total = centering + unexplained + pdc + time_loss + interaction diff --git a/processbehavior/maximum_information.py b/processbehavior/maximum_information.py index d32db0f..1094ea1 100644 --- a/processbehavior/maximum_information.py +++ b/processbehavior/maximum_information.py @@ -200,13 +200,13 @@ def assess_maximum_information( if n < 2: raise ValidationError(f'Maximum information analysis requires at least 2 R2 values, got {n}.') - # --- XmR statistics on R2 (10-12 manual Chapter 12) --- + # --- XmR statistics on R2 (11-7 manual Chapter 11) --- r2_mean = float(np.mean(r2_values)) mr_values = np.abs(np.diff(r2_values)) - r2_mR = float(np.mean(mr_values)) # Eq 12-1 - sigma_hat = r2_mR / D2_N2 # Eq 12-4 + r2_mR = float(np.mean(mr_values)) # Eq 11-1 + sigma_hat = r2_mR / D2_N2 # Eq 11-4 - # Limits: mean ± E2 * mR (Eqs 12-10, 12-11) + # Limits: mean ± E2 * mR (Eqs 11-8, 11-9) upl = r2_mean + XMR_LIMIT_MULTIPLIER * r2_mR lpl = r2_mean - XMR_LIMIT_MULTIPLIER * r2_mR @@ -222,7 +222,7 @@ def assess_maximum_information( mr_series[0] = np.nan mr_series[1:] = mr_values - # R chart limits (Eq 12-5) + # R chart limits (Eq 11-5) r_upl = R_UPPER_LIMIT_MULTIPLIER * r2_mR xmr_data = pd.DataFrame( diff --git a/processbehavior/plotting/renderers.py b/processbehavior/plotting/renderers.py index 12bf9e1..80b852b 100644 --- a/processbehavior/plotting/renderers.py +++ b/processbehavior/plotting/renderers.py @@ -693,7 +693,7 @@ def _add_histogram_stats( vline_kw['annotation_font_size'] = theme.annotation_font_size fig.add_vline(**vline_kw, **subplot_kw) - # Std deviation lines (±1σ, ±2σ, ±3σ): Y-bar ± k·S; k = 3 is the manual's 3-sigma band (10-12 manual Eq 10-15) + # Std deviation lines (±1σ, ±2σ, ±3σ): Y-bar ± k·S; k = 3 is the manual's 3-sigma band (11-7 manual Eq 9-12) if std is None or n < 2 or not np.isfinite(std) or std <= 0: return diff --git a/processbehavior/residual_calculator.py b/processbehavior/residual_calculator.py index d2c9a35..e1935ea 100644 --- a/processbehavior/residual_calculator.py +++ b/processbehavior/residual_calculator.py @@ -4,8 +4,8 @@ This module calculates the Bishop VAS residuals (R1-R5) that decompose total variation into interpretable components. They estimate the terms of the re-parameterized VAS model Y_ktn = mu + rho_k + tau_t + rhotau_kt + e'_ktn -(10-12 manual Eq 10-3), with the means estimated as unweighted means of the cell -means (Eq 10-9; on an incomplete grid, SDS 4-6, over the cells present): +(11-7 manual Eq 9-2), with the means estimated as unweighted means of the cell +means (Figure 13.3, the estimates of the model parameters; on an incomplete grid, SDS 4-6, over the cells present): - R1: Total deviation from grand mean (pure algebra) - R2: Within-cell (unexplained) variation (structure-dependent) @@ -14,13 +14,13 @@ - R5: Factor effects + unexplained (pure algebra given R2) R2 is the ONLY residual whose calculation varies by structure (equation numbers -are Bishop's 10-12 manual): -- exact (ADS 1): R2 = Y - Ȳ_kt (Eq 14-3, the error estimate of Eq 10-10), when every cell has n >= 2 +are Bishop's 11-7 manual): +- exact (ADS 1): R2 = Y - Ȳ_kt (Eq 13-3, the estimate of the model's error term in Eq 9-1), when every cell has n >= 2 - ma2 (ADS 2 & 3): when any cell has n = 1. The data are first stripped of the process mean and the condition and period effects, Z = Y - Ȳ_k - Ȳ_t + Ȳ - (Eq 14-4 / 14-9). Z is then differenced along the condition-then-time stream + (Eq 13-4 / 13-9). Z is then differenced along the condition-then-time stream and divided by twice the R2 scale factor c(K, M), which puts it back on the - noise scale: R2_j = (Z_j - Z_{j-1}) / (2c) (Eqs 14-5..14-8 / 14-10..14-13). + noise scale: R2_j = (Z_j - Z_{j-1}) / (2c) (Eqs 13-5..13-8 / 13-10..13-13). "ma2" names Bishop's size-2 moving average of Z; R2_j = Z_j minus that average. All other residuals (R1, R3, R4, R5) are pure algebraic transformations @@ -162,7 +162,7 @@ def calculate_time_means(df: pd.DataFrame, response_var: str, time_var: str) -> def calculate_cell_means(df: pd.DataFrame, response_var: str, rsg_var_name: str, time_var: str) -> pd.Series: """ - Calculate cell means (Ȳ_kt) - average for each (factor × time) cell (10-12 manual Eq 10-7). + Calculate cell means (Ȳ_kt) - average for each (factor × time) cell (11-7 manual Eq 9-6). Broadcasts the mean for each cell to all rows in that cell. @@ -208,7 +208,7 @@ def calculate_r1_residual(df: pd.DataFrame, response_var: str, grand_mean: float """ Calculate R1 residual: total deviation from grand mean. - R1 = Y - Ȳ (10-12 manual Eq 14-2) + R1 = Y - Ȳ (11-7 manual Eq 13-2) R1 represents the total variation of each observation around the overall average. It's the foundation for all other residuals. @@ -248,7 +248,7 @@ def calculate_r3_residual( """ Calculate R3 residual: interaction effects (factor × time). - R3 = Y - Ȳ_k - Ȳ_t + Ȳ (10-12 manual Eq 14-14, interaction part) + R3 = Y - Ȳ_k - Ȳ_t + Ȳ (11-7 manual Eq 13-14, interaction part) R3 captures the interaction between factors and time. It represents variation that can't be explained by factor effects or time effects alone. @@ -285,7 +285,7 @@ def calculate_r4_residual(time_means: pd.Series, grand_mean: float, r2: pd.Serie """ Calculate R4 residual: time effects + unexplained. - R4 = Ȳ_t - Ȳ + R2 (10-12 manual Eq 14-16) + R4 = Ȳ_t - Ȳ + R2 (11-7 manual Eq 13-16) R4 represents time effects plus within-cell variation. Used to assess if time contributes meaningful variation. @@ -311,7 +311,7 @@ def calculate_r5_residual(factor_means: pd.Series, grand_mean: float, r2: pd.Ser """ Calculate R5 residual: factor effects + unexplained. - R5 = Ȳ_k - Ȳ + R2 (10-12 manual Eq 14-19) + R5 = Ȳ_k - Ȳ + R2 (11-7 manual Eq 13-19) R5 represents factor effects plus within-cell variation. Used to assess if factors contribute meaningful variation and to calculate @@ -343,7 +343,7 @@ def r2_scale_factor(n_conditions: int, n_observations: int) -> float: """ R2 scale factor c(K, M) for the singleton design states. - Bishop's 10-12 manual Eq 14-8 (SDS 2, where M = KT) and Eq 14-13 (SDS 3, 5, 6):: + Bishop's 11-7 manual Eq 13-8 (SDS 2, where M = KT) and Eq 13-13 (SDS 3, 5, 6):: c(K, M) = sqrt( (K - 1) / (2K) * (1 - K / (M - 1)) ) @@ -382,15 +382,15 @@ def calculate_r2( R2 is the ONLY structure-dependent residual. ``r2_method`` comes from ``SDSRegistry.get_r2_method()`` (exact iff every cell has n >= 2) and drives - the branch; equation numbers are Bishop's 10-12 manual. + the branch; equation numbers are Bishop's 11-7 manual. - - ``'exact'`` (ADS 1): R2 = Y - Ȳ_kt (Eq 14-3). - - ``'ma2'`` (ADS 2 and 3): Z = Y - Ȳ_k - Ȳ_t + Ȳ (Eq 14-4 / 14-9), using the + - ``'exact'`` (ADS 1): R2 = Y - Ȳ_kt (Eq 13-3). + - ``'ma2'`` (ADS 2 and 3): Z = Y - Ȳ_k - Ȳ_t + Ȳ (Eq 13-4 / 13-9), using the unweighted cell-mean averages already on the frame, so the process mean and the condition and period effects are removed before anything is differenced. In the canonical condition-then-time order (``sort_key``), - R2_j = (Z_j - Z_{j-1}) / (2 c(K, M)) for j >= 2 (Eqs 14-5..14-8 / - 14-10..14-13). The stream runs across condition boundaries with no grouping; + R2_j = (Z_j - Z_{j-1}) / (2 c(K, M)) for j >= 2 (Eqs 13-5..13-8 / + 13-10..13-13). The stream runs across condition boundaries with no grouping; only the first observation has no value (NaN). Because Z carries no condition or period effect, the step from one condition's last period to the next condition's first period no longer carries a level difference. @@ -515,7 +515,7 @@ def calculate_vas_residuals( # Step 2: Derive marginal means from cell means (unweighted means analysis) # Bishop VAS uses mean of cell means, giving each experimental - # condition equal weight regardless of sample size within cells (10-12 manual Eq 10-9). + # condition equal weight regardless of sample size within cells (11-7 manual Figure 13.3). cell_means_unique = out.groupby([spec.rsg_var_name, spec.time_var], observed=True)['Ybar_kt'].first() grand_mean = cell_means_unique.mean() @@ -540,7 +540,7 @@ def calculate_vas_residuals( raise RuntimeError("R2 method 'exact' requires every cell to have n >= 2") out['R2'] = calculate_r2(out, y, r2_method, n_conditions=int(out[spec.rsg_var_name].nunique())) - # Step 5: Calculate R3 (interaction effects), Eq 14-14 + # Step 5: Calculate R3 (interaction effects), Eq 13-14 # Unified formula: Ybar_kt - Ybar_k - Ybar_t + Ybar + R2 # For exact (ADS 1): R2 = Y - Ybar_kt, so this simplifies to # Y - Ybar_k - Ybar_t + Ybar. @@ -611,9 +611,9 @@ def calculate_r6_residuals(df: pd.DataFrame, groupby_key: str | list[str], recen """Return a new frame carrying R6 (and RCR6 when recentered). Never mutates ``df``. R6 = α_i + R2, where α_i is the factor-level effect: the condition effects ρ̂_k = R5 - R2 - averaged over the level's rows (10-12 manual Eqs 14-19, 14-24, 14-25). + averaged over the level's rows (11-7 manual Eqs 13-19, 13-24, 13-25). - The manual writes α_i as the level's mean of R5 (Eq 14-24, "≈ α_i"). R5 = ρ̂_k + R2, so + The manual writes α_i as the level's mean of R5 (Eq 13-24, "≈ α_i"). R5 = ρ̂_k + R2, so that mean also carries the level's average R2, which VAS leaves out: Tom's Medicare run (10/3/2026, slide 52) centres the organisation-effects chart at 10831.1, the effect without it, where the mean of R5 gives 10829.8. On ADS 1, R2 averages to zero inside diff --git a/processbehavior/sds_detector.py b/processbehavior/sds_detector.py index b0e9008..9a281d6 100644 --- a/processbehavior/sds_detector.py +++ b/processbehavior/sds_detector.py @@ -227,7 +227,7 @@ def residual_charts(self) -> list[tuple[str, str]]: if not self.vas_residuals_supported: return [] - # R1 (10-12 manual Eq 14-2): the response centred at 0 — a location shift, so only the + # R1 (11-7 manual Eq 13-2): the response centred at 0 — a location shift, so only the # location charts. X is unconditional: charting the centred raw data is the whole # diagnostic. No S/mR — R1 differs from the response by a constant, so its # dispersion charts would duplicate the response's exactly. @@ -379,11 +379,11 @@ class SDSRegistry: **SDS 2**: Semi-Complete - No replication - Min N_kt = 1 and Max N_kt = 1 (all cells are singletons) - - R2 by the scaled difference of the effect-adjusted series (10-12 manual Eqs 14-4..14-8) + - R2 by the scaled difference of the effect-adjusted series (11-7 manual Eqs 13-4..13-8) **SDS 3**: Semi-Complete - Partial replication - Min N_kt = 1 and Max N_kt ≥ 2 (mixed) - - R2 by the same scaled difference over every observation (Eqs 14-9..14-13); + - R2 by the same scaled difference over every observation (Eqs 13-9..13-13); Xbar/S is the default chart **Incomplete** (has empty cells, requires sampling plan to detect): @@ -697,7 +697,7 @@ def get_sds_characteristics(self, sds: int) -> dict: 3: { 'description': 'Partial replication (mixed n=1 and n≥2)', 'replication_type': 'partial', - 'r2_method': 'ma2', # Any singleton: scaled difference over every observation (Eqs 14-9..14-13) + 'r2_method': 'ma2', # Any singleton: scaled difference over every observation (Eqs 13-9..13-13) 'capabilities': ['all_residuals', 'partial_interactions', 'main_effects'], 'interaction_analysis': 'partial', 'variance_decomposition': True, @@ -713,7 +713,7 @@ def get_sds_characteristics(self, sds: int) -> dict: 5: { 'description': 'Incomplete grid without replication (has 0s, max=1)', 'replication_type': 'none', - 'r2_method': 'ma2', # No replication: scaled difference (Eqs 14-9..14-13) + 'r2_method': 'ma2', # No replication: scaled difference (Eqs 13-9..13-13) 'capabilities': ['partial_vas', 'stratification'], 'interaction_analysis': False, 'variance_decomposition': True, # VAS supported via the scaled-difference R2 @@ -758,9 +758,9 @@ def get_r2_method(self, stats: StructureStats) -> R2Method: Notes ----- - Deterministic rule based on Bishop's 10-12 manual: - - Eq 14-3 (exact): R2 = Y - Ȳ_kt, requires replication in every cell - - Eqs 14-4..14-13 (ma2): Z = Y - Ȳ_k - Ȳ_t + Ȳ, differenced over the whole + Deterministic rule based on Bishop's 11-7 manual: + - Eq 13-3 (exact): R2 = Y - Ȳ_kt, requires replication in every cell + - Eqs 13-4..13-13 (ma2): Z = Y - Ȳ_k - Ȳ_t + Ȳ, differenced over the whole condition-then-time sequence (j = 2..M, no grouping by cell) and divided by twice the R2 scale factor c(K, M). Once any cell is a singleton this applies to every observation; the replicated cells' within-cell @@ -1321,7 +1321,7 @@ def get_analysis_plan(sds: int, min_cell_size: int = 0) -> SDSAnalysisPlan: has_time=True, has_replication='partial', valid_charts=['Histogram', 'Xbar', 'S', 'X', 'mR'], - recommended_chart='Xbar', # 10-12 manual: Xbar/S is VAS's default for ADS 3 + recommended_chart='Xbar', # 11-7 manual: Xbar/S is VAS's default for ADS 3 invalid_charts=[], vas_residuals_supported=True, residuals_available=['R1', 'R2', 'R3', 'R4', 'R5', 'R6'], @@ -1336,7 +1336,7 @@ def get_analysis_plan(sds: int, min_cell_size: int = 0) -> SDSAnalysisPlan: 'Pilot studies with targeted replication', ], limitations=[ - 'R2 uses the scaled difference over the full sequence (10-12 manual Eqs 14-9..14-13); ' + 'R2 uses the scaled difference over the full sequence (11-7 manual Eqs 13-9..13-13); ' 'the replicated cells contribute no separate within-cell estimate', 'Variance estimates less precise than SDS 1', 'May have unequal subgroup sizes', diff --git a/processbehavior/spc_constants.py b/processbehavior/spc_constants.py index 3eb50ec..4a9b2fd 100644 --- a/processbehavior/spc_constants.py +++ b/processbehavior/spc_constants.py @@ -30,9 +30,9 @@ SIGMA_MULTIPLIER = 3 # Moving-range constants for n = 2 (consecutive pairs), as in Bishop's VAS manual -# (10-12 manual Eq 12-2: E(R) = d2·sigma, sigma_R = d3·sigma): sigma = mR / d2 (Eq 12-4), -# X limits = X ± 3·mR/d2 (Eqs 12-10, 12-11), and the moving-range upper limit -# = (1 + 3·d3/d2)·mR (Eq 12-5). The manual prints the +# (11-7 manual Eq 11-2: E(R) = d2·sigma, sigma_R = d3·sigma): sigma = mR / d2 (Eq 11-4), +# X limits = X ± 3·mR/d2 (Eqs 11-8, 11-9), and the moving-range upper limit +# = (1 + 3·d3/d2)·mR (Eq 11-5). The manual prints the # rounded values 2.66 and 3.268; Bishop's VAS software computes them from d2 and # d3 without rounding, and so does this library. D2_N2 = 1.128 @@ -55,13 +55,13 @@ def c4(n: int) -> float: Calculate c4 bias constant for Xbar and S charts. The c4 constant corrects for bias in the standard deviation estimate - when using subgroups: E(S) = c4(n)·sigma, the manual's alpha_N (10-12 manual - Eq 10-11), so S / c4(n) is unbiased for sigma (Eq 10-12). It approaches 1.0 + when using subgroups: E(S) = c4(n)·sigma, the manual's alpha_N (11-7 manual + Eq 9-8), so S / c4(n) is unbiased for sigma (Eq 9-9). It approaches 1.0 as n increases. c4 is exact for every n, as the manual's alpha_N is: it is computed with log-gamma, so it has no overflow limit. Since sigma_S uses sqrt(1 - alpha_N²) - (Eq 10-14), an S chart's limits close in on S-bar gradually as n grows: + (Eq 9-11), an S chart's limits close in on S-bar gradually as n grows: S-bar ± 3% at n = 5000, ± 1.5% at 20,000. Parameters @@ -132,8 +132,8 @@ def b3(n: int, sigma_multiplier: float = 3) -> float: Notes ----- Formula: b3(n) = 1 - sigma_multiplier/c4(n) * sqrt(1 - c4(n)²), i.e. S-bar - less sigma_multiplier estimated standard deviations of S (10-12 manual - Eq 10-14), the S chart's lower limit (Eq 11-8). + less sigma_multiplier estimated standard deviations of S (11-7 manual + Eq 9-11), the S chart's lower limit (Eq 10-8). For small subgroup sizes (n < 6), the raw b3 formula yields negative values. Since standard deviations cannot be negative, these values are @@ -183,8 +183,8 @@ def b4(n: int, sigma_multiplier: float = 3) -> float: Notes ----- Formula: b4(n) = 1 + sigma_multiplier/(c4(n)) * sqrt(1 - c4(n)^2), i.e. S-bar - plus sigma_multiplier estimated standard deviations of S (10-12 manual - Eq 10-14), the S chart's upper limit (Eq 11-9). + plus sigma_multiplier estimated standard deviations of S (11-7 manual + Eq 9-11), the S chart's upper limit (Eq 10-9). The b4 constant decreases as subgroup size increases, approaching the value of 1 + 3*sqrt(1-1) = 1 for very large subgroups. @@ -315,7 +315,7 @@ def calculate_limits( f'The limits calculation for {limits_type} requires (mean, and mR). Got: mean={mean}, mR={mR}' ) - # LPL = X̄ - (E2 * mR), UPL = X̄ + (E2 * mR) (10-12 manual Eqs 12-10, 12-11) + # LPL = X̄ - (E2 * mR), UPL = X̄ + (E2 * mR) (11-7 manual Eqs 11-8, 11-9) lpl = mean - (XMR_LIMIT_MULTIPLIER * mR) upl = mean + (XMR_LIMIT_MULTIPLIER * mR) @@ -324,7 +324,7 @@ def calculate_limits( raise ValueError(f'The limits calculation for {limits_type} requires (mR). Got: mR={mR}') # LPL = 0 (ranges cannot be negative) - # UPL = mR * D4 (10-12 manual Eq 12-5) + # UPL = mR * D4 (11-7 manual Eq 11-5) lpl = 0 upl = mR * R_UPPER_LIMIT_MULTIPLIER @@ -335,10 +335,10 @@ def calculate_limits( def subgroup_sigma_hat(s, n) -> float: - """Bishop's estimate of sigma from R rational subgroups (10-12 manual Eq 11-5). + """Bishop's estimate of sigma from R rational subgroups (11-7 manual Eq 10-5). The average of each subgroup's own unbiased estimate: sigma_hat = (1/R) sum S_r / c4(N_r). - With equal subgroup sizes this is S-bar / c4(N) (Eq 11-6); with unequal sizes it is not, + With equal subgroup sizes this is S-bar / c4(N) (Eq 10-6); with unequal sizes it is not, because each S_r is unbiased-corrected with its own c4 before averaging. VAS computes it this way (Tom's Medicare run: one subgroup of 3 among 23 of 4). """ @@ -361,8 +361,8 @@ def calculate_limits_vectorized( """Array form of :func:`calculate_limits` — same formulae, whole columns at once. ``sigma`` (Xbar only) gives sigma_hat directly, for subgroups of unequal size: the limits - are then mean ± multiplier · sigma / sqrt(N_r) for each subgroup (10-12 manual Eqs 11-16, - 11-17), with sigma from :func:`subgroup_sigma_hat`. Without it, ``sd`` is S-bar and each + are then mean ± multiplier · sigma / sqrt(N_r) for each subgroup (11-7 manual Eqs 10-16, + 10-17), with sigma from :func:`subgroup_sigma_hat`. Without it, ``sd`` is S-bar and each row's limits use S-bar / c4(N), which is the same thing when every N is equal. :func:`calculate_limits` stays the scalar reference and is unchanged: it is what the @@ -413,7 +413,7 @@ def _per_n(func, sizes, *args): ) sizes = np.asarray(N) if sigma is not None: - # Eqs 11-16/11-17: half-width = multiplier * sigma_hat / sqrt(N_r) + # Eqs 10-16/10-17: half-width = multiplier * sigma_hat / sqrt(N_r) half = (sigma_multiplier * np.asarray(sigma, dtype=float)) / np.sqrt(sizes) else: # Wd = S / c4(n), then half-width = (multiplier * Wd) / sqrt(n) @@ -435,13 +435,13 @@ def _per_n(func, sizes, *args): f'The limits calculation for {limits_type} requires (mean, and mR). ' f'Got: mean={mean}, mR={mR}' ) - half = XMR_LIMIT_MULTIPLIER * np.asarray(mR) # 10-12 manual Eqs 12-10, 12-11 + half = XMR_LIMIT_MULTIPLIER * np.asarray(mR) # 11-7 manual Eqs 11-8, 11-9 lpl, upl = np.asarray(mean) - half, np.asarray(mean) + half elif limits_type == 'R': if mR is None: raise ValueError(f'The limits calculation for {limits_type} requires (mR). Got: mR={mR}') - upl = np.asarray(mR) * R_UPPER_LIMIT_MULTIPLIER # 10-12 manual Eq 12-5 + upl = np.asarray(mR) * R_UPPER_LIMIT_MULTIPLIER # 11-7 manual Eq 11-5 lpl = np.zeros_like(upl) else: @@ -523,19 +523,19 @@ def calibrated_limits( lims = pd.Series({'lpl': mean - half, 'upl': mean + half}, index=['lpl', 'upl']) elif limits_type == 'S': assert N is not None - center = c4(N) * sigma # E(S) = alpha_N·sigma (10-12 manual Eq 10-11) + center = c4(N) * sigma # E(S) = alpha_N·sigma (11-7 manual Eq 9-8) lims = calculate_limits( limits_type='S', mean=0, sd=c4(N) * sigma, N=N, round_to=round_to, sigma_multiplier=n_sigma, ) elif limits_type == 'XmR': - center = mean # mean ± 3·sigma (10-12 manual Eq 12-6) + center = mean # mean ± 3·sigma (11-7 manual Eq 11-6) lims = calculate_limits( limits_type='XmR', mean=mean, sd=0, N=0, mR=D2_N2 * sigma, round_to=round_to, ) elif limits_type == 'R': - center = D2_N2 * sigma # E(R) = d2·sigma; limits 0 … D4·d2·sigma (10-12 manual Eqs 12-2, 12-3) + center = D2_N2 * sigma # E(R) = d2·sigma; limits 0 … D4·d2·sigma (11-7 manual Eqs 11-2, 11-3) lims = calculate_limits( limits_type='R', mean=0, sd=0, N=0, mR=D2_N2 * sigma, round_to=round_to, @@ -705,7 +705,7 @@ def suggest_chart_name(name: str) -> str: # labelled "Noise / Unexplained variation", which is R2's meaning, in a table that # also mapped `within_cell` to R2 and so contradicted itself. RESIDUAL_LABELS: dict[str, str] = { - # 10-12 manual Eq 14-2 (centring the original PM data at 0): R1 = Y_ktn - Ybar.., the + # 11-7 manual Eq 13-2 (centring the original PM data at 0): R1 = Y_ktn - Ybar.., the # response re-expressed as +/- about zero. Not a within-subgroup quantity. 'R1': 'Response Centered at 0', 'R2': 'Within-Cell', diff --git a/processbehavior/study.py b/processbehavior/study.py index 4647e94..2e47754 100644 --- a/processbehavior/study.py +++ b/processbehavior/study.py @@ -1091,7 +1091,7 @@ def dataset(self) -> pd.DataFrame: - rsg: Rational subgroup identifier - Ybar, Ybar_k, Ybar_t, Ybar_kt: Hierarchical means - R1-R5: VAS residuals (where applicable for the SDS) - - RCR1-RCR5: Re-centred residuals, R + Ybar (10-12 manual Eq 14-26; RCR2 = Ybar_kt + R2) + - RCR1-RCR5: Re-centred residuals, R + Ybar (11-7 manual Eq 13-26; RCR2 = Ybar_kt + R2) Returns a copy to preserve immutability. The frame itself is immutable after formulate() — execute() never adds columns to it; request @@ -1372,12 +1372,12 @@ def available_analysis_methods(self) -> pd.DataFrame: { 'method': 'Capability', 'available': self.analytical_design_state.sds > 0, - 'description': 'Process capability against specification limits (Ch. 16)', + 'description': 'Process capability against specification limits (Ch. 15)', }, { 'method': 'Loss Function', 'available': vas_available, - 'description': 'Taguchi loss decomposition into 5 components (Ch. 15)', + 'description': 'Taguchi loss decomposition into 5 components (Ch. 14)', }, { 'method': 'Maximum Information', @@ -1796,7 +1796,7 @@ def capability( window: tuple | None = None, ) -> CapabilityResult: """ - Assess process capability against specification limits (Ch. 16). + Assess process capability against specification limits (Ch. 15). Parameters ---------- @@ -1847,7 +1847,7 @@ def capability( def loss_function(self, target: float | None = None) -> LossResult: """ - Assess Taguchi Loss Function decomposition (Ch. 15). + Assess Taguchi Loss Function decomposition (Ch. 14). Decomposes expected loss into 5 components: centering, unexplained, PDC, time, and PDC×time interaction. diff --git a/tests/test_ads3_charts.py b/tests/test_ads3_charts.py index f54793e..4c9da88 100644 --- a/tests/test_ads3_charts.py +++ b/tests/test_ads3_charts.py @@ -1,4 +1,4 @@ -"""ADS 3 (partial replication) charts follow the 10-12 manual: Xbar/S is the default chart. +"""ADS 3 (partial replication) charts follow the 11-7 manual: Xbar/S is the default chart. Tom Bishop's VAS run of 10/3/2026 draws PM SDS 3 as Xbar/S throughout (original data and the PDC×PT interaction); the per-condition centre lines on that deck are the unweighted mean of every diff --git a/tests/test_capability.py b/tests/test_capability.py index 8dd6782..13fc290 100644 --- a/tests/test_capability.py +++ b/tests/test_capability.py @@ -1,5 +1,5 @@ """ -Tests for process capability analysis (Bishop Chapter 16). +Tests for process capability analysis (Bishop Chapter 15, 11-7 manual). Covers: - SpecLimits validation @@ -964,7 +964,7 @@ def test_custom_round_to(self): # ============================================================================ -# Potential centre (y_bar + mean R2) — 10-12 manual potential capability +# Potential centre (y_bar + mean R2) — 11-7 manual potential capability # ============================================================================ diff --git a/tests/test_loss_function.py b/tests/test_loss_function.py index 539903e..be75b53 100644 --- a/tests/test_loss_function.py +++ b/tests/test_loss_function.py @@ -1,5 +1,5 @@ """ -Tests for Taguchi Loss Function Analysis (Bishop Ch. 15). +Tests for Taguchi Loss Function Analysis (Bishop Ch. 14, 11-7 manual). Uses validation dataset (PBTESTDATABASE_T100.csv) as ground truth. """ @@ -75,7 +75,7 @@ def study_no_time(pb): class TestValidationPMSDS1: - """Match Tom's VAS run of 10/3/2026 (10-1 manual), PM SDS 1 deck slides 30-31.""" + """Match Tom's VAS run of 10/3/2026 (11-7 manual), PM SDS 1 deck slides 30-31.""" def test_unstructured_percentages(self, study_sds1): """5-component Pareto matches Tom's reference.""" @@ -172,7 +172,7 @@ def test_decomposition_default_target(self, study_sds1): class TestUnexplained: def test_sds1_uses_per_cell(self, study_sds1): - """SDS 1 (replicated): unexplained is the average cell variance (Eq 15-16), no c4.""" + """SDS 1 (replicated): unexplained is the average cell variance (Eq 14-14), no c4.""" result = study_sds1.loss_function() assert result.sds == 1 df = study_sds1.dataset @@ -180,7 +180,7 @@ def test_sds1_uses_per_cell(self, study_sds1): assert result.unexplained == pytest.approx(expected, rel=1e-12) def test_sds2_percentages(self, study_sds2): - """SDS 2 Pareto (Eqs 15-17, 15-20) matches Tom's VAS run of 10/3/2026 (PM SDS 2 deck, slide 30).""" + """SDS 2 Pareto (Eqs 14-15, 14-18) matches Tom's VAS run of 10/3/2026 (PM SDS 2 deck, slide 30).""" result = study_sds2.loss_function(target=237.0) assert result.sds == 2 @@ -191,7 +191,7 @@ def test_sds2_percentages(self, study_sds2): assert round(result.pct_time, 1) == 2.9 def test_sds3_percentages(self, study_sds3): - """SDS 3 Pareto (Eqs 15-17, 15-20) matches Tom's VAS run of 10/3/2026 (PM SDS 3 deck, slide 30).""" + """SDS 3 Pareto (Eqs 14-15, 14-18) matches Tom's VAS run of 10/3/2026 (PM SDS 3 deck, slide 30).""" result = study_sds3.loss_function(target=237.0) assert result.sds == 3 @@ -202,14 +202,14 @@ def test_sds3_percentages(self, study_sds3): assert round(result.pct_time, 1) == 2.8 def test_ads23_unexplained_is_r2_variance(self, study_sds2): - """Eq 15-17: ADS 2/3 unexplained is the sample variance of R2, nothing more.""" + """Eq 14-15: ADS 2/3 unexplained is the sample variance of R2, nothing more.""" result = study_sds2.loss_function(target=237.0) r2 = study_sds2.dataset['R2'].dropna() assert result.unexplained == pytest.approx(r2.var(ddof=1), rel=1e-12) @pytest.mark.parametrize('response', ['PM SDS 2', 'PM SDS 3']) def test_ads23_parts_add_to_total_loss(self, pb, response): - """Eqs 15-18 to 15-20: with no floor, the parts add to (Ybar - T)^2 + S^2 exactly.""" + """Eqs 14-16 to 14-18: with no floor, the parts add to (Ybar - T)^2 + S^2 exactly.""" study = pb.formulate(response=response, factors=['FACTOR 1', 'FACTOR 2'], time='PRODUCTION TIME') result = study.loss_function(target=237.0) assert result.interaction > 0.0 @@ -217,7 +217,7 @@ def test_ads23_parts_add_to_total_loss(self, pb, response): assert result.total == pytest.approx(result.centering + s2, rel=1e-12) def test_ads23_interaction_floored_at_zero_on_pure_noise(self): - """PM INERT (pure noise) under the PM SDS 5 pattern: Eq 15-20 remainder is negative, set to 0. + """PM INERT (pure noise) under the PM SDS 5 pattern: Eq 14-18 remainder is negative, set to 0. Tom's VAS run of 10/3/2026 (PM5 INERT deck, slide 30) reports 84.0 / 14.0 / 1.5 / 0.5 / 0.0. """ diff --git a/tests/test_maximum_information.py b/tests/test_maximum_information.py index 0d11878..7f3aaf1 100644 --- a/tests/test_maximum_information.py +++ b/tests/test_maximum_information.py @@ -201,7 +201,7 @@ def test_plot_with_bins(self): # --------------------------------------------------------------------------- -# 10-1 manual R2: limits match Tom's VAS run of 10/3/2026 (slide 33 of each deck) +# 11-7 manual R2: limits match Tom's VAS run of 10/3/2026 (slide 33 of each deck) # --------------------------------------------------------------------------- diff --git a/tests/test_phased_xbar_s.py b/tests/test_phased_xbar_s.py index 4eb51c7..2a1b37a 100644 --- a/tests/test_phased_xbar_s.py +++ b/tests/test_phased_xbar_s.py @@ -241,7 +241,7 @@ def test_unequal_n_uses_eq_11_5_within_the_phase(self): a = stages['A'] stds = np.array([np.std(g, ddof=1) for g in a]) sizes = np.array([len(g) for g in a]) - sigma_hat = np.mean(stds / np.array([c4(n) for n in sizes])) # Eq 11-5 within the phase + sigma_hat = np.mean(stds / np.array([c4(n) for n in sizes])) # Eq 10-5 within the phase rows = xbar[xbar['phase'] == 'A'] assert (rows['upl'] - rows['center']).to_numpy() == pytest.approx(3 * sigma_hat / np.sqrt(sizes), abs=1e-9) assert phases.loc['A', 'sigma_hat'] == pytest.approx(sigma_hat, abs=1e-9) @@ -267,7 +267,7 @@ def test_average_n_mode_is_applied_per_phase(self): for stage, subgroups in stages.items(): n_avg = np.mean([len(g) for g in subgroups]) s_bar = np.mean([np.std(g, ddof=1) for g in subgroups]) - half = 3 * (s_bar / c4(int(n_avg))) / np.sqrt(n_avg) # Eq 11-20, as the unphased chart applies it + half = 3 * (s_bar / c4(int(n_avg))) / np.sqrt(n_avg) # Eq 10-20, as the unphased chart applies it rows = xbar[xbar['phase'] == stage] assert phases.loc[stage, 'N'] == pytest.approx(n_avg) assert rows['upl'].nunique() == 1 @@ -298,7 +298,7 @@ def test_one_reading_subgroups_sit_on_xbar_and_not_on_s(self, reference_studies) assert len(s) == len(xbar) - len(singles) phases = result.phase_statistics('Xbar').set_index('phase') for _, row in singles.iterrows(): - # N_r = 1: center ± n_sigma · sigma_hat (10-12 manual Eqs 11-16, 11-17) + # N_r = 1: center ± n_sigma · sigma_hat (11-7 manual Eqs 10-16, 10-17) assert row['upl'] - row['center'] == pytest.approx(3 * phases.loc[row['phase'], 'sigma_hat'], abs=3e-3) def test_single_subgroup_phase_has_finite_limits(self): diff --git a/tests/test_r2_unavailable.py b/tests/test_r2_unavailable.py index b300546..5bedacc 100644 --- a/tests/test_r2_unavailable.py +++ b/tests/test_r2_unavailable.py @@ -1,7 +1,7 @@ """R2 is unavailable, with a stated reason, when the layout leaves the R2 scale factor undefined. With a one-observation subgroup, R2 is the condition-and-period-adjusted series differenced and divided -by the R2 scale factor c(K, M) (10-12 manual Eqs 14-4..14-13). c needs at least two process design +by the R2 scale factor c(K, M) (11-7 manual Eqs 13-4..13-13). c needs at least two process design conditions; with one condition the adjusted series is identically zero and R2 carries no information. """ diff --git a/tests/test_r6_effect.py b/tests/test_r6_effect.py index d415fb2..fc4fbed 100644 --- a/tests/test_r6_effect.py +++ b/tests/test_r6_effect.py @@ -1,7 +1,7 @@ """R6's factor-level effect is the condition effect averaged over the level, without R2's level mean. -R6 = α_i + R2 (10-12 manual Eq 14-25). The manual takes α_i as the level's mean of R5 (Eq 14-24, "≈ α_i"); -R5 = ρ̂_k + R2 (Eq 14-19), so that mean also carries the level's average R2. VAS leaves it out: Tom's +R6 = α_i + R2 (11-7 manual Eq 13-25). The manual takes α_i as the level's mean of R5 (Eq 13-24, "≈ α_i"); +R5 = ρ̂_k + R2 (Eq 13-19), so that mean also carries the level's average R2. VAS leaves it out: Tom's Medicare run (10/3/2026, slide 52) centres the organisation-effects Xbar chart at 10831.1, which is the effect without it (the mean of R5 gives 10829.8). On ADS 1 the two agree, because R2 averages to zero in every cell. diff --git a/tests/test_recentering.py b/tests/test_recentering.py index 91d6b40..665a438 100644 --- a/tests/test_recentering.py +++ b/tests/test_recentering.py @@ -1,4 +1,4 @@ -"""Re-centred residuals follow Bishop's 10-12 manual Eq 14-26: RCR = R + Ybar. +"""Re-centred residuals follow Bishop's 11-7 manual Eq 13-26: RCR = R + Ybar. Re-centring moves a residual onto the measurement scale without changing what it shows, so an effect chart of a re-centred residual plots that effect plus R2, centred on the grand mean. diff --git a/tests/test_residual_calculator.py b/tests/test_residual_calculator.py index 150f96f..dd46160 100644 --- a/tests/test_residual_calculator.py +++ b/tests/test_residual_calculator.py @@ -148,7 +148,7 @@ def test_calculate_r1_residual_sum_is_zero(): ], ) def test_r2_scale_factor_matches_manual_layouts(k, m, expected): - """c(K, M) = sqrt((K-1)/(2K) * (1 - K/(M-1))), Eq 14-8 / 14-13.""" + """c(K, M) = sqrt((K-1)/(2K) * (1 - K/(M-1))), Eq 13-8 / 13-13.""" assert r2_scale_factor(k, m) == pytest.approx(expected, abs=5e-6) @@ -177,7 +177,7 @@ def test_calculate_r2_exact(): pd.testing.assert_series_equal(result, expected, check_names=False) -# R2 for the singleton design states (10-12 manual Eqs 14-4..14-13). +# R2 for the singleton design states (11-7 manual Eqs 13-4..13-13). # Hand-worked 2 x 3 grid, one observation per cell: # A: 10, 12, 11 B: 14, 13, 18 # Ybar = 13; Ybar_k = 11 (A), 15 (B); Ybar_t = 12, 12.5, 14.5 @@ -557,7 +557,7 @@ def test_no_nan_residuals_all_r2_methods(r2_method, df_key, sds1_df, sds2_df, sd # ============================================================================ -# Test: R2 properties for the singleton design states (10-12 manual) +# Test: R2 properties for the singleton design states (11-7 manual) # ============================================================================ diff --git a/tests/test_sds.py b/tests/test_sds.py index 07c6c55..95a7cbe 100644 --- a/tests/test_sds.py +++ b/tests/test_sds.py @@ -269,8 +269,8 @@ def test_rcr_centered_residuals(self): ds = ads.analysis_dataset TOL = 1e-10 - # Re-centred residuals, Bishop's 10-12 manual Eq 14-26: RCR = R + Ybar (R1, R3, R4, R5); - # RCR2 is not part of Eq 14-26 and stays Ybar_kt + R2. + # Re-centred residuals, Bishop's 11-7 manual Eq 13-26: RCR = R + Ybar (R1, R3, R4, R5); + # RCR2 is not part of Eq 13-26 and stays Ybar_kt + R2. expected_rcr1 = ds['Ybar'] + ds['R1'] assert (ds['RCR1'] - expected_rcr1).abs().max() <= TOL, 'RCR1 formula incorrect' @@ -279,7 +279,7 @@ def test_rcr_centered_residuals(self): for k in (3, 4, 5): expected = ds['Ybar'] + ds[f'R{k}'] - assert (ds[f'RCR{k}'] - expected).abs().max() <= TOL, f'RCR{k} must equal R{k} + Ybar (Eq 14-26)' + assert (ds[f'RCR{k}'] - expected).abs().max() <= TOL, f'RCR{k} must equal R{k} + Ybar (Eq 13-26)' # End-to-end reconstruction identity # The fixture is SDS 1 with n_min=2, guaranteeing the exact R2 method diff --git a/tests/test_sds_detector.py b/tests/test_sds_detector.py index 41cb8a8..9d6843d 100644 --- a/tests/test_sds_detector.py +++ b/tests/test_sds_detector.py @@ -648,7 +648,7 @@ def test_realistic_scenario_designed_experiment_no_replication(detector): class TestR2ChartAvailability: """Tests for R2_S vs R2_X selection based on actual cell sizes. - Per Bishop Section 20.6.1: + Per Bishop's 11-7 manual (Chapter 12, the analytic design states, and Section 13.3.9): - R2_S is available when rational subgroups have n>=2 - R2_X is used when n=1 (no within-cell variation) """ @@ -1317,7 +1317,7 @@ class TestR2MethodSelection: """``get_r2_method`` has two answers, and structure alone decides between them. The docs once described a third, per-cell "hybrid" for mixed cell sizes. It was removed in - cb037f8: with any singleton cell, R2 is Bishop's ungrouped MA2 (10-12 manual Eqs 14-9..14-13) for every + cb037f8: with any singleton cell, R2 is Bishop's ungrouped MA2 (11-7 manual Eqs 13-9..13-13) for every observation, which the ADS 3 reference assertions validate. These pin the vocabulary to that. """ diff --git a/tests/test_spc_constants.py b/tests/test_spc_constants.py index acd933e..7569429 100644 --- a/tests/test_spc_constants.py +++ b/tests/test_spc_constants.py @@ -131,7 +131,7 @@ def test_b3_b4_raises_on_invalid_n(func): class TestMovingRangeConstants: - """The XmR and mR constants follow Bishop's VAS manual (10-12 manual Eqs 12-4, 12-5, 12-10, 12-11) + """The XmR and mR constants follow Bishop's VAS manual (11-7 manual Eqs 11-4, 11-5, 11-8, 11-9) and are computed from d2 = 1.128 and d3 = 0.8525 without rounding, as VAS does.""" def test_building_blocks_are_the_manuals(self): diff --git a/tests/test_unequal_subgroups.py b/tests/test_unequal_subgroups.py index 1d09f9b..e5e8c7c 100644 --- a/tests/test_unequal_subgroups.py +++ b/tests/test_unequal_subgroups.py @@ -1,4 +1,4 @@ -"""Xbar/S charts with subgroups of unequal size follow the 10-12 manual (Eqs 11-5, 11-7..11-9, 11-16, 11-17). +"""Xbar/S charts with subgroups of unequal size follow the 11-7 manual (Eqs 10-5, 10-7..10-9, 10-16, 10-17). sigma_hat is the average of each subgroup's own unbiased S_r / c4(N_r); each subgroup's Xbar limits are centre ± 3 sigma_hat / sqrt(N_r), and its S limits are B3(N_r) S-bar and B4(N_r) S-bar. N_r counts the @@ -20,7 +20,7 @@ def test_subgroup_sigma_hat_is_average_of_unbiased_estimates(): s, n = np.array([1.0, 2.0, 1.5]), np.array([3, 4, 5]) assert subgroup_sigma_hat(s, n) == pytest.approx(np.mean(s / np.array([c4(3), c4(4), c4(5)]))) - # equal sizes: S-bar / c4(N) (Eq 11-6) + # equal sizes: S-bar / c4(N) (Eq 10-6) assert subgroup_sigma_hat(s, np.array([4, 4, 4])) == pytest.approx(s.mean() / c4(4)) @@ -36,7 +36,7 @@ def test_subgroup_size_counts_the_charted_values(medicare): table = medicare.execute(chart='S', by=['ACO'], value='R6').chart_table('S') sizes = medicare.dataset.groupby('ACO')['R2'].count() assert sizes['ACO-001'] == 3 and (sizes.drop('ACO-001') == 4).all() - # ACO-001's S limit uses B4(3), the others B4(4) (Eqs 11-8, 11-9) + # ACO-001's S limit uses B4(3), the others B4(4) (Eqs 10-8, 10-9) s_bar = table['center'].iloc[0] assert table['upl'].iloc[0] == pytest.approx(b4(3) * s_bar, rel=1e-9) assert table['upl'].iloc[1] == pytest.approx(b4(4) * s_bar, rel=1e-9) @@ -74,8 +74,8 @@ def with_one_reading_cells(): def test_one_reading_subgroups_are_charted_at_three_sigma(with_one_reading_cells): - """A one-reading subgroup has no S_r, so it takes no part in sigma_hat (Eq 11-5), but it is charted with - limits centre ± 3 sigma_hat (Eqs 11-16/17 with N_r = 1), as VAS draws it (PM SDS 3 and 6, 10/3/2026).""" + """A one-reading subgroup has no S_r, so it takes no part in sigma_hat (Eq 10-5), but it is charted with + limits centre ± 3 sigma_hat (Eqs 10-16/17 with N_r = 1), as VAS draws it (PM SDS 3 and 6, 10/3/2026).""" study = with_one_reading_cells result = study.execute(chart='Xbar', by=[], companion=True) xbar, s = result.chart_table('Xbar'), result.chart_table('S') diff --git a/tests/test_vas_chart_constants.py b/tests/test_vas_chart_constants.py index b4a8459..2cc8f4f 100644 --- a/tests/test_vas_chart_constants.py +++ b/tests/test_vas_chart_constants.py @@ -2,8 +2,8 @@ The expected values are read from Bishop's VAS output decks of 29 September 2026 (Medicare and PM SDS 2) and, for the R2 chart, 3 October 2026, at the precision the charts print. They pin the -XmR constants E2 = 3/1.128 and D4 = 1 + 3(0.8525)/1.128 (10-12 manual Eqs 12-4, 12-5, 12-10, -12-11): with the rounded 2.66 and 3.268 the Medicare limits miss by 0.6 and 1.0, and the PM SDS 2 +XmR constants E2 = 3/1.128 and D4 = 1 + 3(0.8525)/1.128 (11-7 manual Eqs 11-4, 11-5, 11-8, +11-9): with the rounded 2.66 and 3.268 the Medicare limits miss by 0.6 and 1.0, and the PM SDS 2 moving-range limit prints 2.84. """ @@ -88,7 +88,7 @@ def test_moving_range_limits(self, pm_sds_2): class TestMedicareR2Chart: - """VAS run of 10/3/2026 (10-1 manual R2), MEDICARE PCE slide 60: R2 chart -1.15 (-1045.19 / 1042.88). + """VAS run of 10/3/2026 (11-7 manual R2), MEDICARE PCE slide 60: R2 chart -1.15 (-1045.19 / 1042.88). R2 here is the condition-and-period-adjusted series differenced along the ACO-then-year stream and divided by twice the R2 scale factor c(24, 96); the raw-difference R2 of diff --git a/validation/e2e_bishop_report.py b/validation/e2e_bishop_report.py index eb547f0..186d52c 100644 --- a/validation/e2e_bishop_report.py +++ b/validation/e2e_bishop_report.py @@ -271,7 +271,7 @@ def run_validation(run, pb, study, ref): # noqa: C901 # Pre-compute results we'll need computed = {} - # Xbar/S wherever some subgroups are replicated: ADS 1, and ADS 3, where the 10-12 manual makes + # Xbar/S wherever some subgroups are replicated: ADS 1, and ADS 3, where the 11-7 manual makes # Xbar/S VAS's default (Tom's 10/3/2026 run draws PM SDS 3 that way). ADS 2 uses X/mR. xbar_s = run.ads in (1, 3) @@ -293,7 +293,7 @@ def run_validation(run, pb, study, ref): # noqa: C901 chart='S', by=factors, value='R6' ) # PT effects (pages 22-23) — Xbar/S of R4 by time for all SDS (R4 carries the period effect, - # 10-12 manual Eq 14-16). When charted by=[time], each time subgroup has multiple factor + # 11-7 manual Eq 13-16). When charted by=[time], each time subgroup has multiple factor # levels, giving n>1 subgroups even in SDS 2, so Xbar/S is correct. computed['pt_effects_xbar'] = study.execute( chart='Xbar', by=[time], value='R4', recentered=True @@ -378,7 +378,7 @@ def stats_from(result_obj, chart_type): lpl_val = stats.get('lpl') upl_val = stats.get('upl') if stats.get('limits_vary') and lpl_val is None and upl_val is None: - # Limits step with subgroup size (10-12 manual Eqs 11-16/17). VAS/Minitab labels the + # Limits step with subgroup size (11-7 manual Eqs 10-16/17). VAS/Minitab labels the # limits of the last subgroup, so that is the value its slide carries. table = safe_chart_table(result_obj, chart_type) if table is not None and len(table): @@ -888,7 +888,7 @@ def generate_myst_summary(all_results, aux_by_run, runs): '', '- The limits of the two factor S charts (pages 25 and 27) in the known-effects run, whose', ' subgroups have 1000 and 2000 values. Its VAS run of 6 October 2026 held αN', - " (10-12 manual Eq 10-11) at 0.9995 from N = 343 up, where Minitab's GAMMA function overflows,", + " (11-7 manual Eq 9-8) at 0.9995 from N = 343 up, where Minitab's GAMMA function overflows,", " so its limits follow N = 500 instead of the subgroup's own N. Dr. Bishop's reruns of the PM", ' runs (13 October 2026) compute αN for every N, and every R6 S chart in them is', ' checked and matches; the known-effects run waits for the same rerun.', diff --git a/validation/truth_report.py b/validation/truth_report.py index a26ac29..7293d4e 100644 --- a/validation/truth_report.py +++ b/validation/truth_report.py @@ -15,7 +15,7 @@ For each design state it reports the noise estimate behind potential capability, potential and current Pp against their true values, the loss shares against the true shares of the expected loss -(Ybar - T)^2 + sigma_total^2 (10-1 manual Eq 15-18), and how far the estimated condition and time +(Ybar - T)^2 + sigma_total^2 (11-7 manual Eq 14-16), and how far the estimated condition and time effects sit from the true ones. """ from __future__ import annotations