diff --git a/CHANGELOG.md b/CHANGELOG.md index 722f505..560e99f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,21 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +### Added +- **The design report states where sigma comes from at the observed structure and + series length.** One line under ``Structure``: ``Series length: T=4. Sigma for X/mR + rests on 3 moving ranges.`` With every cell replicated it reads ``Series length: T=1. + Sigma rests on within-cell replication.`` — a single time period is a complete study. + Exposed as ``study.series_length`` (``SeriesLengthPrecision``, frozen) and + ``study.series_length_description``. The result also carries ``within_cell_df`` and, + for the moving-range case, ``mr_interval_80`` — the p10/p90 of MRbar/d2 by n = 3..30 + from ``MR_SIGMA_INTERVAL_80`` (``spc_constants``), generated by + ``validation/short_series_bands.py`` (contributed by @rabujamra in #114 / #119) — for + callers who want the numbers; the report prints only the sentence. It is a fact beside + the design state, not a judgment on it: no threshold, label, or warning, and the design + state, recommended chart and analysis menu do not move with T (pinned by test). Closes + the series-length half of #114. + ### Changed - **The R2 method for partial replication is reported as ``ma2``, not ``hybrid``.** The docs described a per-cell hybrid for design state 3 (exact deviation where a cell has diff --git a/docs/getting-started/key-concepts.md b/docs/getting-started/key-concepts.md index 284b514..18a2c2c 100644 --- a/docs/getting-started/key-concepts.md +++ b/docs/getting-started/key-concepts.md @@ -47,6 +47,11 @@ See [DS Definitions](../reference/sds_definitions.md) for the formal classificat The DS determines: - Which chart types are valid - How within-group variance is estimated (R2 method: exact or ma2) +- Where the sigma behind the limits comes from at the observed structure and series length. + The design report states this as one line under `Structure` (`study.series_length`): with + any singleton cell, `Sigma for X/mR rests on T − 1 moving ranges`; with every cell + replicated, `Sigma rests on within-cell replication`, and T does not enter. No threshold + or warning is attached; the analyst weighs it. - Which VAS residuals can be computed - What conclusions you can draw diff --git a/docs/reference/api.md b/docs/reference/api.md index 2dcca3d..f90460c 100644 --- a/docs/reference/api.md +++ b/docs/reference/api.md @@ -163,6 +163,8 @@ study.analytical_design_state # SDSResult — ADS: what the analysis runs at study.plan_design_state # SDSResult | None — PDS, when plan= was given study.ads_reason # str — machine-readable reason ('full_replication', ...) study.ads_description # str — human-readable description +study.series_length # SeriesLengthPrecision — where sigma comes from at this structure, and how precise at this T +study.series_length_description # str — the 'Series length:' sentence from the design report (may be empty) ``` `SDSResult.sds` is the design state on Bishop's 1–6 reference scale — see the diff --git a/processbehavior/series_length.py b/processbehavior/series_length.py new file mode 100644 index 0000000..19cdf3a --- /dev/null +++ b/processbehavior/series_length.py @@ -0,0 +1,105 @@ +""" +Series-length precision of the sigma estimate. + +One fact for the design report: where the sigma behind the natural process +limits comes from at the observed structure and series length. No threshold, +label, or warning is attached; the analyst weighs it (Bishop: judgment belongs +to the analyst, not to a rule). + +Two sources of sigma: + +- Every cell replicated (ADS 1): within-cell deviation. The number of time + periods does not enter; T = 1 is a complete study. +- Any singleton cell (ADS 2, ADS 3): the 2-point moving range over the ordered + sequence, T - 1 ranges for T time points. + +The sentence names the source and the count. The within-cell degrees of freedom +and the 80% interval of the moving-range estimate at this T +(``MR_SIGMA_INTERVAL_80``, from #114) ride on the result for callers who want +them; the report does not print them. +""" + +from __future__ import annotations + +from dataclasses import dataclass + +from .spc_constants import MR_SIGMA_INTERVAL_80 + +_TABLE_MAX = max(MR_SIGMA_INTERVAL_80) + + +@dataclass(frozen=True) +class SeriesLengthPrecision: + """ + Where sigma comes from, and how precise it is, at the observed series length. + + Attributes + ---------- + T : int or None + Distinct time points on the analysis dataset; None when no time variable. + n_moving_ranges : int or None + T - 1 when the moving range carries sigma and T >= 2; else None. + within_cell_df : int + Sum over cells of (N_kt - 1); 0 when no cell is replicated. + sigma_from_time : bool + True when sigma rests on moving ranges (any singleton cell); False when + it rests on within-cell replication. + mr_interval_80 : tuple of (float, float) or None + (p10, p90) of the moving-range sigma estimate relative to the truth on a + stable process, from ``MR_SIGMA_INTERVAL_80``; only when + ``sigma_from_time`` and 3 <= T <= 30. + description : str + The sentence printed in the design report, e.g. + ``T=4. Sigma for X/mR rests on 3 moving ranges.`` Empty when there is + nothing to say (no time variable and no replication). + """ + + T: int | None + n_moving_ranges: int | None + within_cell_df: int + sigma_from_time: bool + mr_interval_80: tuple[float, float] | None + description: str + + +def assess_series_length(T: int | None, within_cell_df: int, sigma_from_time: bool) -> SeriesLengthPrecision: + """ + Build the precision statement for a study. + + The sentence names only the source of sigma at this structure and T. The + numbers behind it (``within_cell_df``, ``mr_interval_80``) ride on the + result for callers who want them; the report does not print them. + + Parameters + ---------- + T : int or None + Distinct time points on the analysis dataset (None if no time variable). + within_cell_df : int + Sum over cells of (N_kt - 1) on the analysis dataset. + sigma_from_time : bool + True when the limits rest on the moving range (any singleton cell). + + Returns + ------- + SeriesLengthPrecision + """ + if not sigma_from_time: + if T is None: + description = 'Sigma rests on within-cell replication.' + else: + description = f'T={T}. Sigma rests on within-cell replication.' + return SeriesLengthPrecision(T, None, within_cell_df, False, None, description) + + if T is None: + return SeriesLengthPrecision(None, None, within_cell_df, True, None, '') + + n_mr = T - 1 if T >= 2 else 0 + interval = MR_SIGMA_INTERVAL_80.get(T) if T <= _TABLE_MAX else None + + if T < 2: + description = f'T={T}. Sigma for X/mR rests on the moving range, which needs at least 2 time points.' + elif T == 2: + description = 'T=2. Sigma for X/mR rests on a single moving range.' + else: + description = f'T={T}. Sigma for X/mR rests on {n_mr} moving ranges.' + return SeriesLengthPrecision(T, n_mr, within_cell_df, True, interval, description) diff --git a/processbehavior/spc_constants.py b/processbehavior/spc_constants.py index b6a65cb..9e9d38f 100644 --- a/processbehavior/spc_constants.py +++ b/processbehavior/spc_constants.py @@ -690,3 +690,47 @@ def suggest_chart_name(name: str) -> str: 'design_condition_main_effects': 'R5', 'design_factor_main_effects': 'R6', } + + +# ============================================================================ +# Series-length precision of the moving-range sigma estimate +# ============================================================================ + +# 80% interval (p10, p90) of MRbar/d2 relative to the true sigma, for a stable +# normal process observed at n points, n = 3..30. Generated by +# validation/short_series_bands.py (seed 20260903, 100,000 replicates, an +# independent stream per n), contributed by @rabujamra in #114 / PR #119. Read +# as: at n=4, 80% of moving-range sigma estimates fall between 0.43x and 1.68x +# the truth. The spread narrows monotonically with n, so the n=30 row bounds +# every longer series. Stated as a fact in the design report; the library +# attaches no threshold, label, or warning to it. +MR_SIGMA_INTERVAL_80: dict[int, tuple[float, float]] = { + 3: (0.337, 1.810), + 4: (0.430, 1.676), + 5: (0.488, 1.593), + 6: (0.528, 1.526), + 7: (0.565, 1.484), + 8: (0.595, 1.453), + 9: (0.618, 1.418), + 10: (0.637, 1.396), + 11: (0.656, 1.375), + 12: (0.670, 1.360), + 13: (0.680, 1.344), + 14: (0.694, 1.330), + 15: (0.703, 1.315), + 16: (0.714, 1.310), + 17: (0.722, 1.300), + 18: (0.730, 1.288), + 19: (0.739, 1.280), + 20: (0.745, 1.274), + 21: (0.749, 1.266), + 22: (0.756, 1.258), + 23: (0.761, 1.253), + 24: (0.766, 1.249), + 25: (0.771, 1.243), + 26: (0.775, 1.238), + 27: (0.780, 1.234), + 28: (0.784, 1.229), + 29: (0.787, 1.224), + 30: (0.790, 1.220), +} diff --git a/processbehavior/study.py b/processbehavior/study.py index 1cabd2c..b2e343c 100644 --- a/processbehavior/study.py +++ b/processbehavior/study.py @@ -29,6 +29,7 @@ from .exceptions import ChartNotAvailableError, FactorNotFoundError, ValidationError from .residual_calculator import resolve_r6_groupby from .sds_detector import SDSRegistry +from .series_length import SeriesLengthPrecision, assess_series_length from .spc_constants import ( ALL_RESIDUALS, RESIDUAL_ALIASES, @@ -173,6 +174,7 @@ class DesignReport: _pds_result: SDSResult | None = None # Plan Design State _ods_result: SDSResult | None = None # Observed Design State _ads_result: SDSResult | None = None # Analytical Design State + _series_length: SeriesLengthPrecision | None = None # Where sigma comes from, and how precise @property def factors_table(self) -> pd.DataFrame: @@ -765,6 +767,12 @@ def __repr__(self) -> str: lines.append('') lines.append(f' Structure: {self.structure_summary}') + # Series length: where sigma comes from at this structure and T. A fact beside the + # design state, not a verdict on it — ADS 2 is ADS 2 whether T is 4 or 400 (#114). + # No threshold, label, or warning is attached. + if self._series_length is not None and self._series_length.description: + lines.append(f' Series length: {self._series_length.description}') + # Available analyses (derived from ADS) if self._ads_result and self._ads_result.sds > 0: lines.extend(self._repr_available_analyses()) @@ -1111,6 +1119,39 @@ def ads_reason(self) -> str | None: """ return self.analytical_design_state.reason + @property + def series_length(self) -> SeriesLengthPrecision: + """ + Where sigma comes from at this structure and series length. + + A fact stated beside the design state, not a judgment on it: no + threshold, label, or warning. With any singleton cell the limits rest on + T - 1 moving ranges; with every cell replicated they rest on within-cell + replication and T does not enter. The result also carries + ``within_cell_df`` and, for the moving-range case, ``mr_interval_80`` + (from ``MR_SIGMA_INTERVAL_80``) for callers who want the numbers; the + report prints only the sentence. + + Returns + ------- + SeriesLengthPrecision + """ + ads_df = self._ads.analysis_dataset + T = None + if self._spec.time_var and self._spec.time_var in ads_df.columns: + T = int(ads_df[self._spec.time_var].nunique()) + within_cell_df = 0 + if 'cell_key' in ads_df.columns and len(ads_df): + sizes = ads_df.groupby('cell_key', observed=True).size() + within_cell_df = int((sizes - 1).sum()) + sigma_from_time = self.analytical_design_state.sds != 1 + return assess_series_length(T, within_cell_df, sigma_from_time) + + @property + def series_length_description(self) -> str: + """The series-length sentence printed in the design report (may be empty).""" + return self.series_length.description + @property def ads_description(self) -> str: """ @@ -1667,6 +1708,7 @@ def design(self) -> DesignReport: _pds_result=self._pds_result, _ods_result=self._sds_result, _ads_result=self._ads.analytical_design_state, + _series_length=self.series_length, ) def capability( diff --git a/tests/test_series_length.py b/tests/test_series_length.py new file mode 100644 index 0000000..2325ede --- /dev/null +++ b/tests/test_series_length.py @@ -0,0 +1,153 @@ +"""Series-length precision statement (#114). + +One line under Structure in the design report: where sigma comes from at the observed +structure, and how precise it is at the observed T. A fact beside the design state, never a +verdict on it. These tests pin the sentences, the table they draw on, and the invariant that +the design state, recommendation and chart menu do not move with T. +""" + +import dataclasses + +import numpy as np +import pandas as pd +import pytest + +import processbehavior as pb +from processbehavior.series_length import SeriesLengthPrecision, assess_series_length +from processbehavior.spc_constants import D2_N2, MR_SIGMA_INTERVAL_80 + + +def _stable(T, seed=7): + rng = np.random.default_rng(seed) + return pd.DataFrame({'t': range(T), 'y': rng.normal(5.5, 1.2, T)}) + + +# --------------------------------------------------------------------------- +# The pure function +# --------------------------------------------------------------------------- + + +class TestAssessSeriesLength: + def test_unreplicated_in_table_range_states_the_interval(self): + r = assess_series_length(T=4, within_cell_df=0, sigma_from_time=True) + assert r.n_moving_ranges == 3 + assert r.mr_interval_80 == MR_SIGMA_INTERVAL_80[4] + assert r.description == 'T=4. Sigma for X/mR rests on 3 moving ranges.' + + def test_unreplicated_beyond_table_is_bounded_by_the_last_row(self): + r = assess_series_length(T=400, within_cell_df=0, sigma_from_time=True) + assert r.n_moving_ranges == 399 and r.mr_interval_80 is None + assert r.description == 'T=400. Sigma for X/mR rests on 399 moving ranges.' + + def test_two_points_is_a_single_moving_range(self): + r = assess_series_length(T=2, within_cell_df=0, sigma_from_time=True) + assert r.n_moving_ranges == 1 and r.mr_interval_80 is None + assert r.description == 'T=2. Sigma for X/mR rests on a single moving range.' + + def test_replicated_does_not_depend_on_T(self): + r1 = assess_series_length(T=1, within_cell_df=25, sigma_from_time=False) + assert r1.n_moving_ranges is None and r1.mr_interval_80 is None + assert r1.within_cell_df == 25 + assert r1.description == 'T=1. Sigma rests on within-cell replication.' + r_none = assess_series_length(T=None, within_cell_df=25, sigma_from_time=False) + assert r_none.description == 'Sigma rests on within-cell replication.' + + def test_partial_replication_names_the_moving_range_and_keeps_the_df(self): + r = assess_series_length(T=8, within_cell_df=32, sigma_from_time=True) + assert r.description == 'T=8. Sigma for X/mR rests on 7 moving ranges.' + assert r.within_cell_df == 32 and r.mr_interval_80 == MR_SIGMA_INTERVAL_80[8] + + def test_no_time_and_no_replication_says_nothing(self): + r = assess_series_length(T=None, within_cell_df=0, sigma_from_time=True) + assert r.description == '' + + def test_no_verdict_words(self): + for T in (2, 3, 4, 6, 12, 20, 30, 31, 400): + text = assess_series_length(T, 0, True).description.lower() + for word in ( + 'short', + 'adequate', + 'provisional', + 'warning', + 'insufficient', + 'too few', + 'degrees of freedom', + 'estimate', + 'stable process', + ): + assert word not in text, (T, word) + + def test_result_is_frozen(self): + r = assess_series_length(T=4, within_cell_df=0, sigma_from_time=True) + assert isinstance(r, SeriesLengthPrecision) + with pytest.raises(dataclasses.FrozenInstanceError): + r.T = 5 # type: ignore[misc] + + +# --------------------------------------------------------------------------- +# The table +# --------------------------------------------------------------------------- + + +class TestIntervalTable: + def test_covers_3_through_30_and_narrows_monotonically(self): + assert sorted(MR_SIGMA_INTERVAL_80) == list(range(3, 31)) + spreads = [hi / lo for lo, hi in (MR_SIGMA_INTERVAL_80[n] for n in range(3, 31))] + assert all(a > b for a, b in zip(spreads, spreads[1:], strict=False)) + + @pytest.mark.parametrize('n', [4, 20]) + def test_rows_regenerate_from_the_validation_script_seed(self, n): + """validation/short_series_bands.py: default_rng([20260903, n]), 100,000 reps.""" + rng = np.random.default_rng([20260903, n]) + z = rng.standard_normal((100_000, n)) + v = np.abs(np.diff(z, axis=1)).mean(axis=1) / D2_N2 + p10, p90 = np.percentile(v, [10, 90]) + assert (round(p10, 3), round(p90, 3)) == MR_SIGMA_INTERVAL_80[n] + + +# --------------------------------------------------------------------------- +# Through Study and the design report +# --------------------------------------------------------------------------- + + +class TestStudySurface: + def test_unreplicated_series_reports_moving_ranges(self): + st = pb.formulate(_stable(4), response='y', time='t') + assert st.series_length.T == 4 and st.series_length.n_moving_ranges == 3 + assert st.series_length_description == st.series_length.description + assert 'Series length: T=4. Sigma for X/mR rests on 3 moving ranges' in repr(st.design()) + + def test_design_state_and_menu_do_not_move_with_T(self): + """The invariant #114's first table demonstrates: ADS 2 is ADS 2 at T=3 and T=300.""" + short, long = (pb.formulate(_stable(T), response='y', time='t') for T in (3, 300)) + assert short.analytical_design_state.sds == long.analytical_design_state.sds == 2 + assert short.ads_reason == long.ads_reason + assert short.recommended_chart == long.recommended_chart == 'X' + assert short.valid_charts == long.valid_charts + assert short.series_length.n_moving_ranges == 2 and long.series_length.n_moving_ranges == 299 + + def test_replicated_single_period_reports_within_cell_df(self): + """Tom's T=1 file: 5 conditions x 6 replicates, one time point — a complete study.""" + df = pd.DataFrame( + { + 'TIME': [1] * 30, + 'TEMP': np.repeat([0, 25, 50, 75, 100], 6), + 'LIFE': np.arange(30, dtype=float) + 55, + } + ) + st = pb.formulate(df, response='LIFE', factors=['TEMP'], time='TIME') + assert st.analytical_design_state.sds == 1 + sl = st.series_length + assert sl.T == 1 and not sl.sigma_from_time and sl.within_cell_df == 25 + assert 'Series length: T=1. Sigma rests on within-cell replication.' in repr(st.design()) + + def test_no_time_variable_with_replication_still_states_the_source(self): + df = pd.DataFrame({'TEMP': np.repeat([0, 25, 50], 4), 'LIFE': np.arange(12, dtype=float)}) + st = pb.formulate(df, response='LIFE', factors=['TEMP']) + assert st.series_length.T is None + assert st.series_length.within_cell_df == 9 + assert 'Series length: Sigma rests on within-cell replication.' in repr(st.design()) + + def test_no_warning_is_emitted_for_a_short_series(self, recwarn): + pb.formulate(_stable(3), response='y', time='t').design() + assert not [w for w in recwarn if 'Series length' in str(w.message) or 'moving range' in str(w.message)]