Skills › Scientific Agent Skills › analytical-method-validation › _common.py
_common.py
1000 lines · skills/analytical-method-validation/scripts/_common.py at 5f2bc03
1#!/usr/bin/env python3
2"""Shared statistics and I/O for analytical method validation checks.
4Standard library only. Every distribution function here is implemented from the
5regularised incomplete beta and gamma functions so the scripts run in any
6Python 3.11+ interpreter without numpy or scipy.
8These helpers compute and report. They never decide that a procedure is
9validated, fit for purpose, or acceptable to a regulator -- that judgement
10belongs to the analyst and the quality unit.
11"""
13from __future__ import annotations
15import argparse
16import csv
17import io
18import json
19import math
20import sys
21from dataclasses import dataclass, field
22from pathlib import Path
23from typing import Any, Iterable, Sequence
25# --------------------------------------------------------------------------
26# Limits and exit codes
27# --------------------------------------------------------------------------
29MAX_INPUT_BYTES = 5_000_000
30MAX_ROWS = 20_000
32EXIT_OK = 0
33EXIT_FINDINGS = 1
34EXIT_INPUT_ERROR = 2
36TINY = 1e-300
39class InputError(Exception):
40 """Raised for malformed or out-of-bounds user input."""
43def finite_float(value: str) -> float:
44 """Argparse converter that refuses NaN and infinity."""
45 try:
46 result = float(value)
47 except ValueError as exc:
48 raise argparse.ArgumentTypeError("expected a finite number") from exc
49 if not math.isfinite(result):
50 raise argparse.ArgumentTypeError("expected a finite number")
51 return result
54# --------------------------------------------------------------------------
55# Special functions
56# --------------------------------------------------------------------------
59def _betacf(a: float, b: float, x: float, itmax: int = 400, eps: float = 3e-16) -> float:
60 """Continued fraction for the incomplete beta function (Lentz's method)."""
61 qab, qap, qam = a + b, a + 1.0, a - 1.0
62 c = 1.0
63 d = 1.0 - qab * x / qap
64 if abs(d) < TINY:
65 d = TINY
66 d = 1.0 / d
67 h = d
68 for m in range(1, itmax + 1):
69 m2 = 2 * m
70 aa = m * (b - m) * x / ((qam + m2) * (a + m2))
71 d = 1.0 + aa * d
72 if abs(d) < TINY:
73 d = TINY
74 c = 1.0 + aa / c
75 if abs(c) < TINY:
76 c = TINY
77 d = 1.0 / d
78 h *= d * c
79 aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2))
80 d = 1.0 + aa * d
81 if abs(d) < TINY:
82 d = TINY
83 c = 1.0 + aa / c
84 if abs(c) < TINY:
85 c = TINY
86 d = 1.0 / d
87 delta = d * c
88 h *= delta
89 if abs(delta - 1.0) < eps:
90 break
91 return h
94def betainc(a: float, b: float, x: float) -> float:
95 """Regularised incomplete beta function I_x(a, b)."""
96 if x <= 0.0:
97 return 0.0
98 if x >= 1.0:
99 return 1.0
100 log_front = (
101 math.lgamma(a + b)
102 - math.lgamma(a)
103 - math.lgamma(b)
104 + a * math.log(x)
105 + b * math.log1p(-x)
106 )
107 front = math.exp(log_front)
108 if x < (a + 1.0) / (a + b + 2.0):
109 return front * _betacf(a, b, x) / a
110 return 1.0 - front * _betacf(b, a, 1.0 - x) / b
113def gammainc_lower(a: float, x: float) -> float:
114 """Regularised lower incomplete gamma P(a, x)."""
115 if x <= 0.0:
116 return 0.0
117 if x < a + 1.0:
118 # Series representation.
119 term = 1.0 / a
120 total = term
121 n = a
122 for _ in range(1000):
123 n += 1.0
124 term *= x / n
125 total += term
126 if abs(term) < abs(total) * 1e-16:
127 break
128 return total * math.exp(-x + a * math.log(x) - math.lgamma(a))
129 # Continued fraction for Q(a, x), then complement.
130 b = x + 1.0 - a
131 c = 1.0 / TINY
132 d = 1.0 / b
133 h = d
134 for i in range(1, 1000):
135 an = -i * (i - a)
136 b += 2.0
137 d = an * d + b
138 if abs(d) < TINY:
139 d = TINY
140 c = b + an / c
141 if abs(c) < TINY:
142 c = TINY
143 d = 1.0 / d
144 delta = d * c
145 h *= delta
146 if abs(delta - 1.0) < 1e-16:
147 break
148 q = math.exp(-x + a * math.log(x) - math.lgamma(a)) * h
149 return 1.0 - q
152def _bisect_ppf(cdf, target: float, lo: float, hi: float, tol: float = 1e-12) -> float:
153 """Invert a monotone CDF by bisection."""
154 for _ in range(300):
155 mid = 0.5 * (lo + hi)
156 if cdf(mid) < target:
157 lo = mid
158 else:
159 hi = mid
160 if hi - lo < tol * max(1.0, abs(mid)):
161 break
162 return 0.5 * (lo + hi)
165def t_cdf(t: float, df: float) -> float:
166 """CDF of Student's t with df degrees of freedom."""
167 if df <= 0:
168 raise InputError("t distribution needs df > 0")
169 x = df / (df + t * t)
170 half = 0.5 * betainc(0.5 * df, 0.5, x)
171 return half if t <= 0 else 1.0 - half
174def t_ppf(p: float, df: float) -> float:
175 """Quantile of Student's t."""
176 if not 0.0 < p < 1.0:
177 raise InputError("t_ppf needs 0 < p < 1")
178 return _bisect_ppf(lambda t: t_cdf(t, df), p, -1e4, 1e4)
181def chi2_cdf(x: float, df: float) -> float:
182 """CDF of the chi-square distribution."""
183 if x <= 0:
184 return 0.0
185 return gammainc_lower(0.5 * df, 0.5 * x)
188def chi2_ppf(p: float, df: float) -> float:
189 """Quantile of the chi-square distribution."""
190 if not 0.0 < p < 1.0:
191 raise InputError("chi2_ppf needs 0 < p < 1")
192 return _bisect_ppf(lambda x: chi2_cdf(x, df), p, 1e-12, 1e6)
195def f_cdf(x: float, df1: float, df2: float) -> float:
196 """CDF of the F distribution."""
197 if x <= 0:
198 return 0.0
199 return betainc(0.5 * df1, 0.5 * df2, df1 * x / (df1 * x + df2))
202def f_sf(x: float, df1: float, df2: float) -> float:
203 """Upper tail of the F distribution (the p-value for an F test).
205 Computed from the complementary incomplete beta rather than as 1 - cdf,
206 which underflows to exactly 0 for large F and would print a lack-of-fit
207 p-value of 0 in a validation report.
208 """
209 if x <= 0:
210 return 1.0
211 return betainc(0.5 * df2, 0.5 * df1, df2 / (df1 * x + df2))
214def z_ppf(p: float) -> float:
215 """Standard normal quantile."""
216 from statistics import NormalDist
218 return NormalDist().inv_cdf(p)
221# --------------------------------------------------------------------------
222# Descriptive helpers
223# --------------------------------------------------------------------------
226def mean(values: Sequence[float]) -> float:
227 if not values:
228 raise InputError("mean of an empty sequence")
229 return math.fsum(values) / len(values)
232def sample_sd(values: Sequence[float]) -> float:
233 n = len(values)
234 if n < 2:
235 return float("nan")
236 m = mean(values)
237 return math.sqrt(math.fsum((v - m) ** 2 for v in values) / (n - 1))
240def rsd_percent(values: Sequence[float]) -> float:
241 """Relative standard deviation (%CV). NaN when the mean is ~0."""
242 m = mean(values)
243 if abs(m) < 1e-15:
244 return float("nan")
245 return 100.0 * sample_sd(values) / abs(m)
248def median(values: Sequence[float]) -> float:
249 if not values:
250 raise InputError("median of an empty sequence")
251 s = sorted(values)
252 n = len(s)
253 mid = n // 2
254 return s[mid] if n % 2 else 0.5 * (s[mid - 1] + s[mid])
257def sd_confidence_interval(sd: float, df: float, level: float = 0.90) -> tuple[float, float]:
258 """Chi-square confidence interval for a standard deviation."""
259 if df <= 0 or not math.isfinite(sd):
260 return (float("nan"), float("nan"))
261 alpha = 1.0 - level
262 lo_chi = chi2_ppf(1.0 - alpha / 2.0, df)
263 hi_chi = chi2_ppf(alpha / 2.0, df)
264 return (sd * math.sqrt(df / lo_chi), sd * math.sqrt(df / hi_chi))
267# --------------------------------------------------------------------------
268# Regression
269# --------------------------------------------------------------------------
272@dataclass
273class LinearFit:
274 """Weighted least-squares straight-line fit and its diagnostics."""
276 n: int
277 slope: float
278 intercept: float
279 se_slope: float
280 se_intercept: float
281 residual_sd: float
282 r_squared: float
283 r: float
284 df: int
285 residuals: list[float] = field(default_factory=list)
286 fitted: list[float] = field(default_factory=list)
287 weights: list[float] = field(default_factory=list)
289 def predict(self, x: float) -> float:
290 return self.intercept + self.slope * x
292 def slope_ci(self, level: float = 0.95) -> tuple[float, float]:
293 t = t_ppf(0.5 + level / 2.0, self.df)
294 return (self.slope - t * self.se_slope, self.slope + t * self.se_slope)
296 def intercept_ci(self, level: float = 0.95) -> tuple[float, float]:
297 t = t_ppf(0.5 + level / 2.0, self.df)
298 return (self.intercept - t * self.se_intercept, self.intercept + t * self.se_intercept)
301def fit_linear(
302 xs: Sequence[float], ys: Sequence[float], weights: Sequence[float] | None = None
303) -> LinearFit:
304 """Fit y = a + b*x by (optionally weighted) least squares."""
305 n = len(xs)
306 if n != len(ys):
307 raise InputError("x and y must be the same length")
308 if n < 3:
309 raise InputError("a regression needs at least 3 points")
310 w = [1.0] * n if weights is None else [float(v) for v in weights]
311 if len(w) != n:
312 raise InputError("weights must match the number of points")
313 if any(v < 0 for v in w):
314 raise InputError("weights must be non-negative")
316 sw = math.fsum(w)
317 swx = math.fsum(wi * xi for wi, xi in zip(w, xs))
318 swy = math.fsum(wi * yi for wi, yi in zip(w, ys))
319 swxx = math.fsum(wi * xi * xi for wi, xi in zip(w, xs))
320 swxy = math.fsum(wi * xi * yi for wi, xi, yi in zip(w, xs, ys))
321 denom = sw * swxx - swx * swx
322 if abs(denom) < 1e-300:
323 raise InputError("x values are collinear or identical; slope is undefined")
325 slope = (sw * swxy - swx * swy) / denom
326 intercept = (swy - slope * swx) / sw
327 fitted = [intercept + slope * xi for xi in xs]
328 residuals = [yi - fi for yi, fi in zip(ys, fitted)]
329 df = n - 2
330 ss_res = math.fsum(wi * ri * ri for wi, ri in zip(w, residuals))
331 residual_sd = math.sqrt(ss_res / df)
332 se_slope = residual_sd * math.sqrt(sw / denom)
333 se_intercept = residual_sd * math.sqrt(swxx / denom)
335 ybar_w = swy / sw
336 ss_tot = math.fsum(wi * (yi - ybar_w) ** 2 for wi, yi in zip(w, ys))
337 r_squared = 1.0 - ss_res / ss_tot if ss_tot > 0 else float("nan")
338 r = math.copysign(math.sqrt(max(0.0, r_squared)), slope)
340 return LinearFit(
341 n=n,
342 slope=slope,
343 intercept=intercept,
344 se_slope=se_slope,
345 se_intercept=se_intercept,
346 residual_sd=residual_sd,
347 r_squared=r_squared,
348 r=r,
349 df=df,
350 residuals=residuals,
351 fitted=fitted,
352 weights=w,
353 )
356def runs_test(residuals: Sequence[float]) -> dict[str, Any]:
357 """Wald-Wolfowitz runs test on residual signs.
359 ICH Q2(R2) 3.2.2.1 asks for the impact of any non-random pattern in the
360 residual plot to be assessed. Curvature shows up as too few runs.
361 """
362 signs = [1 if r >= 0 else -1 for r in residuals if r != 0]
363 n = len(signs)
364 n_pos = sum(1 for s in signs if s > 0)
365 n_neg = n - n_pos
366 runs = 1 + sum(1 for i in range(1, n) if signs[i] != signs[i - 1]) if n else 0
367 if n_pos < 1 or n_neg < 1 or n < 8:
368 return {
369 "runs": runs,
370 "n_pos": n_pos,
371 "n_neg": n_neg,
372 "z": float("nan"),
373 "p_value": float("nan"),
374 "note": "too few points for a meaningful runs test (need n>=8 with both signs)",
375 }
376 exp = 2.0 * n_pos * n_neg / n + 1.0
377 var = (2.0 * n_pos * n_neg * (2.0 * n_pos * n_neg - n)) / (n * n * (n - 1.0))
378 if var <= 0:
379 return {"runs": runs, "n_pos": n_pos, "n_neg": n_neg, "z": float("nan"),
380 "p_value": float("nan"), "note": "degenerate variance"}
381 z = (runs - exp) / math.sqrt(var)
382 from statistics import NormalDist
384 p = 2.0 * NormalDist().cdf(-abs(z))
385 return {"runs": runs, "n_pos": n_pos, "n_neg": n_neg, "expected_runs": exp,
386 "z": z, "p_value": p, "note": ""}
389def lack_of_fit(xs: Sequence[float], ys: Sequence[float], fit: LinearFit) -> dict[str, Any]:
390 """ANOVA lack-of-fit F test, which needs replicate y at some x levels.
392 This is the statistically meaningful test of a linear calibration model.
393 r-squared is not: it rises with range and is insensitive to curvature.
394 """
395 groups: dict[float, list[float]] = {}
396 for x, y in zip(xs, ys):
397 groups.setdefault(round(float(x), 12), []).append(float(y))
398 k = len(groups)
399 n = len(xs)
400 replicated = sum(1 for vals in groups.values() if len(vals) > 1)
401 df_pe = n - k
402 df_lof = k - 2
403 if df_pe < 1 or df_lof < 1:
404 return {
405 "applicable": False,
406 "levels": k,
407 "replicated_levels": replicated,
408 "reason": "needs replicates at >=1 level and >=3 distinct levels",
409 }
410 # Use the same variance model as the fitted line. Weights must be constant
411 # within a concentration level for this replicate-based partition.
412 level_weights: dict[float, float] = {}
413 for x, weight in zip(xs, fit.weights):
414 level = round(float(x), 12)
415 if level in level_weights and weight != level_weights[level]:
416 return {"applicable": False, "reason": "weights differ within a level"}
417 level_weights[level] = weight
418 ss_pe = math.fsum(
419 level_weights[level] * math.fsum((v - mean(vals)) ** 2 for v in vals)
420 for level, vals in groups.items()
421 )
422 ss_res = math.fsum(w * r * r for w, r in zip(fit.weights, fit.residuals))
423 ss_lof = max(0.0, ss_res - ss_pe)
424 ms_pe = ss_pe / df_pe
425 ms_lof = ss_lof / df_lof
426 if ms_pe <= 0:
427 return {"applicable": False, "levels": k, "replicated_levels": replicated,
428 "reason": "zero pure-error variance; replicates are identical"}
429 f_stat = ms_lof / ms_pe
430 return {
431 "applicable": True,
432 "levels": k,
433 "replicated_levels": replicated,
434 "df_lack_of_fit": df_lof,
435 "df_pure_error": df_pe,
436 "ms_lack_of_fit": ms_lof,
437 "ms_pure_error": ms_pe,
438 "f_statistic": f_stat,
439 "p_value": f_sf(f_stat, df_lof, df_pe),
440 }
443def heteroscedasticity(xs: Sequence[float], residuals: Sequence[float]) -> dict[str, Any]:
444 """Compare residual spread in the lowest and highest thirds of the range.
446 A large ratio is a variance-model diagnostic. Heteroscedasticity alone does
447 not bias OLS coefficients, but can reduce efficiency and low-range precision.
448 """
449 pairs = sorted(zip(xs, residuals), key=lambda p: p[0])
450 n = len(pairs)
451 if n < 6:
452 return {"applicable": False, "reason": "needs at least 6 points"}
453 cut = max(2, n // 3)
454 low = [r for _, r in pairs[:cut]]
455 high = [r for _, r in pairs[-cut:]]
456 var_low = math.fsum(r * r for r in low) / len(low)
457 var_high = math.fsum(r * r for r in high) / len(high)
458 if var_low <= 0:
459 return {"applicable": False, "reason": "zero residual variance in the low third"}
460 ratio = var_high / var_low
461 return {
462 "applicable": True,
463 "n_low": len(low),
464 "n_high": len(high),
465 "variance_ratio_high_over_low": ratio,
466 "sd_ratio": math.sqrt(ratio),
467 }
470# --------------------------------------------------------------------------
471# Variance components (precision)
472# --------------------------------------------------------------------------
475@dataclass
476class PrecisionComponents:
477 """One-way random-effects decomposition of precision."""
479 grand_mean: float
480 n_total: int
481 n_groups: int
482 ms_between: float
483 ms_within: float
484 df_between: int
485 df_within: int
486 sd_repeatability: float
487 sd_between: float
488 sd_intermediate: float
489 balanced: bool
490 n_effective: float
492 def rsd(self, sd: float) -> float:
493 if abs(self.grand_mean) < 1e-15:
494 return float("nan")
495 return 100.0 * sd / abs(self.grand_mean)
497 def satterthwaite_df(self) -> float:
498 """Effective df for the total (intermediate) SD."""
499 var_total = self.sd_intermediate ** 2
500 if var_total <= 0:
501 return float("nan")
502 n = self.n_effective
503 c1 = 1.0 / n
504 c2 = (n - 1.0) / n
505 num = var_total ** 2
506 den = 0.0
507 if self.df_between > 0:
508 den += (c1 * self.ms_between) ** 2 / self.df_between
509 if self.df_within > 0:
510 den += (c2 * self.ms_within) ** 2 / self.df_within
511 return num / den if den > 0 else float("nan")
514def one_way_components(groups: dict[str, Sequence[float]]) -> PrecisionComponents:
515 """Decompose precision into within-group and between-group components.
517 Groups are the intermediate-precision factor -- day, analyst, instrument,
518 or a combined run. Within-group scatter estimates repeatability; the total
519 estimates intermediate precision.
520 """
521 clean = {k: [float(v) for v in vals] for k, vals in groups.items() if len(vals) >= 1}
522 if len(clean) < 2:
523 raise InputError("intermediate precision needs at least 2 groups")
524 if all(len(v) < 2 for v in clean.values()):
525 raise InputError("at least one group needs >=2 replicates to estimate repeatability")
527 counts = [len(v) for v in clean.values()]
528 n_total = sum(counts)
529 k = len(clean)
530 all_values = [v for vals in clean.values() for v in vals]
531 grand = mean(all_values)
533 ss_within = math.fsum(
534 math.fsum((v - mean(vals)) ** 2 for v in vals) for vals in clean.values()
535 )
536 ss_between = math.fsum(len(vals) * (mean(vals) - grand) ** 2 for vals in clean.values())
537 df_within = n_total - k
538 df_between = k - 1
539 ms_within = ss_within / df_within if df_within > 0 else float("nan")
540 ms_between = ss_between / df_between if df_between > 0 else float("nan")
542 balanced = len(set(counts)) == 1
543 if balanced:
544 n_eff = float(counts[0])
545 else:
546 # Standard unbalanced coefficient for the expected mean square.
547 n_eff = (n_total - math.fsum(c * c for c in counts) / n_total) / (k - 1)
549 var_within = max(0.0, ms_within) if math.isfinite(ms_within) else 0.0
550 var_between = 0.0
551 if math.isfinite(ms_between) and math.isfinite(ms_within) and n_eff > 0:
552 var_between = max(0.0, (ms_between - ms_within) / n_eff)
554 return PrecisionComponents(
555 grand_mean=grand,
556 n_total=n_total,
557 n_groups=k,
558 ms_between=ms_between,
559 ms_within=ms_within,
560 df_between=df_between,
561 df_within=df_within,
562 sd_repeatability=math.sqrt(var_within),
563 sd_between=math.sqrt(var_between),
564 sd_intermediate=math.sqrt(var_within + var_between),
565 balanced=balanced,
566 n_effective=n_eff,
567 )
570# --------------------------------------------------------------------------
571# Method comparison
572# --------------------------------------------------------------------------
575def deming(
576 xs: Sequence[float], ys: Sequence[float], lambda_ratio: float = 1.0
577) -> dict[str, Any]:
578 """Deming regression: errors in both variables.
580 `lambda_ratio` is var(error in y) / var(error in x) -- the variance of the
581 random error in the TEST (y) procedure over that in the COMPARATIVE (x) one.
582 Check the direction against the limits, which are unambiguous: as
583 lambda -> infinity the fit converges on the ordinary least-squares slope of
584 y on x (x treated as error-free), and as lambda -> 0 it converges on the
585 inverse regression (y treated as error-free). lambda = 1 means equal error
586 variances and reduces to orthogonal regression.
588 In practice lambda is estimated as (SD of y replicates / SD of x replicates)
589 squared, so equal-precision procedures give 1.
591 Ordinary least squares assumes x is error-free, which is false when
592 comparing procedures; independent error in x can attenuate the slope.
593 """
594 n = len(xs)
595 if n != len(ys):
596 raise InputError("x and y must be the same length")
597 if n < 3:
598 raise InputError("Deming regression needs at least 3 points")
599 if not math.isfinite(lambda_ratio) or lambda_ratio <= 0:
600 raise InputError("lambda_ratio must be > 0")
602 def _fit(xv: Sequence[float], yv: Sequence[float]) -> tuple[float, float]:
603 xb, yb = mean(xv), mean(yv)
604 sxx = math.fsum((x - xb) ** 2 for x in xv)
605 syy = math.fsum((y - yb) ** 2 for y in yv)
606 sxy = math.fsum((x - xb) * (y - yb) for x, y in zip(xv, yv))
607 if abs(sxy) < 1e-300:
608 raise InputError("zero covariance; Deming slope is undefined")
609 term = syy - lambda_ratio * sxx
610 slope = (term + math.sqrt(term * term + 4.0 * lambda_ratio * sxy * sxy)) / (
611 2.0 * sxy
612 )
613 return slope, yb - slope * xb
615 slope, intercept = _fit(xs, ys)
617 # Jackknife standard errors.
618 slopes, intercepts = [], []
619 for i in range(n):
620 xv = list(xs[:i]) + list(xs[i + 1 :])
621 yv = list(ys[:i]) + list(ys[i + 1 :])
622 try:
623 s, a = _fit(xv, yv)
624 except InputError:
625 continue
626 slopes.append(s)
627 intercepts.append(a)
628 if len(slopes) > 2:
629 m = len(slopes)
630 se_slope = math.sqrt((m - 1) / m * math.fsum((s - mean(slopes)) ** 2 for s in slopes))
631 se_int = math.sqrt(
632 (m - 1) / m * math.fsum((a - mean(intercepts)) ** 2 for a in intercepts)
633 )
634 df = m - 2
635 else:
636 se_slope = se_int = float("nan")
637 df = 1
639 t = t_ppf(0.975, df) if df > 0 else float("nan")
640 return {
641 "n": n,
642 "lambda_ratio": lambda_ratio,
643 "slope": slope,
644 "intercept": intercept,
645 "se_slope": se_slope,
646 "se_intercept": se_int,
647 "slope_ci95": (slope - t * se_slope, slope + t * se_slope),
648 "intercept_ci95": (intercept - t * se_int, intercept + t * se_int),
649 "df": df,
650 }
653def passing_bablok(xs: Sequence[float], ys: Sequence[float]) -> dict[str, Any]:
654 """Passing-Bablok for distinct x and a positively associated linear relation.
656 Robust to outliers and does not assume a known error-variance ratio, which
657 is why CLSI EP09-style method comparison work often prefers it.
658 """
659 n = len(xs)
660 if n != len(ys):
661 raise InputError("x and y must be the same length")
662 if n < 5:
663 raise InputError("Passing-Bablok needs at least 5 points")
664 if n > 1000:
665 raise InputError("Passing-Bablok limited to 1000 pairs because pairwise slopes are quadratic")
666 if len(set(xs)) != n:
667 raise InputError("bundled Passing-Bablok does not support tied reference values")
668 xb, yb = mean(xs), mean(ys)
669 if math.fsum((x - xb) * (y - yb) for x, y in zip(xs, ys)) <= 0:
670 raise InputError("bundled Passing-Bablok requires positive association")
672 slopes: list[float] = []
673 for i in range(n):
674 for j in range(i + 1, n):
675 dx = xs[j] - xs[i]
676 dy = ys[j] - ys[i]
677 if dx == 0 and dy == 0:
678 continue
679 if dx == 0:
680 continue # vertical pair carries no finite slope
681 slope = dy / dx
682 if slope != -1.0:
683 slopes.append(slope)
684 if not slopes:
685 raise InputError("no usable pairwise slopes")
687 slopes.sort()
688 n_slopes = len(slopes)
689 shift = sum(1 for s in slopes if s < -1.0)
691 def _shifted_median(offset: int) -> float:
692 idx = n_slopes // 2 + offset
693 if n_slopes % 2:
694 return slopes[min(max(idx, 0), n_slopes - 1)]
695 lo = slopes[min(max(idx - 1, 0), n_slopes - 1)]
696 hi = slopes[min(max(idx, 0), n_slopes - 1)]
697 return 0.5 * (lo + hi)
699 slope = _shifted_median(shift)
700 intercept = median([y - slope * x for x, y in zip(xs, ys)])
702 # Rank-based 95% CI on the slope. M1 and M2 are 1-based order statistics of
703 # the shifted slope list, so both convert to 0-based with the same -1.
704 c = z_ppf(0.975) * math.sqrt(n * (n - 1.0) * (2.0 * n + 5.0) / 18.0)
705 m1 = int(round((n_slopes - c) / 2.0))
706 m2 = n_slopes - m1 + 1
707 if m1 + shift < 1 or m2 + shift > n_slopes:
708 raise InputError("Passing-Bablok confidence ranks exceed finite slope support; use a dedicated implementation")
709 lo_idx = min(max(m1 + shift - 1, 0), n_slopes - 1)
710 hi_idx = min(max(m2 + shift - 1, 0), n_slopes - 1)
711 slope_lo, slope_hi = slopes[lo_idx], slopes[hi_idx]
712 int_lo = median([y - slope_hi * x for x, y in zip(xs, ys)])
713 int_hi = median([y - slope_lo * x for x, y in zip(xs, ys)])
715 return {
716 "n": n,
717 "n_slopes": n_slopes,
718 "slope": slope,
719 "intercept": intercept,
720 "slope_ci95": (slope_lo, slope_hi),
721 "intercept_ci95": (int_lo, int_hi),
722 }
725def bland_altman(
726 xs: Sequence[float], ys: Sequence[float], relative: bool = False
727) -> dict[str, Any]:
728 """Bias and limits of agreement between paired measurements."""
729 n = len(xs)
730 if n != len(ys):
731 raise InputError("x and y must be the same length")
732 if n < 3:
733 raise InputError("Bland-Altman needs at least 3 pairs")
734 means = [0.5 * (x + y) for x, y in zip(xs, ys)]
735 if relative:
736 diffs = []
737 for x, y, m in zip(xs, ys, means):
738 if abs(m) < 1e-15:
739 raise InputError("relative differences need non-zero pair means")
740 diffs.append(100.0 * (y - x) / m)
741 else:
742 diffs = [y - x for x, y in zip(xs, ys)]
744 bias = mean(diffs)
745 sd = sample_sd(diffs)
746 t = t_ppf(0.975, n - 1)
747 se_bias = sd / math.sqrt(n)
748 loa_lo, loa_hi = bias - 1.96 * sd, bias + 1.96 * sd
749 se_loa = sd * math.sqrt(1.0 / n + (1.96 ** 2) / (2.0 * (n - 1)))
751 # Proportional-bias check: does the difference trend with the mean?
752 trend = None
753 try:
754 tf = fit_linear(means, diffs)
755 t_stat = tf.slope / tf.se_slope if tf.se_slope > 0 else float("nan")
756 trend = {
757 "slope": tf.slope,
758 "p_value": 2.0 * (1.0 - t_cdf(abs(t_stat), tf.df)) if math.isfinite(t_stat) else float("nan"),
759 }
760 except InputError:
761 trend = None
763 return {
764 "n": n,
765 "relative": relative,
766 "bias": bias,
767 "sd_differences": sd,
768 "bias_ci95": (bias - t * se_bias, bias + t * se_bias),
769 "loa_lower": loa_lo,
770 "loa_upper": loa_hi,
771 "loa_ci95_halfwidth": t * se_loa,
772 "proportional_bias": trend,
773 }
776def tost_paired(
777 diffs: Sequence[float], margin: float, alpha: float = 0.05
778) -> dict[str, Any]:
779 """Two one-sided tests for equivalence on paired differences.
781 Absence of a significant difference is not evidence of equivalence. TOST
782 tests the hypothesis that actually matters at a method transfer: that the
783 true difference lies inside +/- margin.
784 """
785 n = len(diffs)
786 if n < 2:
787 raise InputError("TOST needs at least 2 differences")
788 if not math.isfinite(margin) or margin <= 0:
789 raise InputError("margin must be > 0")
790 if not 0 < alpha < 0.5:
791 raise InputError("alpha must be between 0 and 0.5")
792 d = mean(diffs)
793 sd = sample_sd(diffs)
794 se = sd / math.sqrt(n)
795 df = n - 1
796 if se <= 0:
797 raise InputError("zero variability; TOST is undefined")
798 t_lower = (d + margin) / se
799 t_upper = (d - margin) / se
800 p_lower = 1.0 - t_cdf(t_lower, df)
801 p_upper = t_cdf(t_upper, df)
802 p = max(p_lower, p_upper)
803 t_crit = t_ppf(1.0 - alpha, df)
804 ci = (d - t_crit * se, d + t_crit * se)
805 return {
806 "n": n,
807 "mean_difference": d,
808 "sd_difference": sd,
809 "margin": margin,
810 "alpha": alpha,
811 "p_lower": p_lower,
812 "p_upper": p_upper,
813 "p_value": p,
814 "ci_1_minus_2alpha": ci,
815 "equivalent": bool(ci[0] > -margin and ci[1] < margin),
816 }
819# --------------------------------------------------------------------------
820# I/O
821# --------------------------------------------------------------------------
824def read_input(path: str | None) -> str:
825 """Read a bounded amount of text from a path or stdin."""
826 if path in (None, "-"):
827 # Reading a terminal would block forever with no indication why, so an
828 # omitted --input becomes an error rather than an apparent hang.
829 if path is None and sys.stdin.isatty():
830 raise InputError("no input given; pass --input FILE, or '-' to read stdin")
831 data = sys.stdin.read(MAX_INPUT_BYTES + 1)
832 else:
833 p = Path(path)
834 if not p.is_file():
835 raise InputError(f"not a file: {path}")
836 if p.stat().st_size > MAX_INPUT_BYTES:
837 raise InputError(f"input larger than {MAX_INPUT_BYTES} bytes")
838 data = p.read_text(encoding="utf-8", errors="replace")
839 if len(data) > MAX_INPUT_BYTES:
840 raise InputError(f"input larger than {MAX_INPUT_BYTES} bytes")
841 return data
844def parse_rows(text: str, path_hint: str | None = None) -> list[dict[str, str]]:
845 """Parse CSV, TSV, or a JSON array of objects into a list of dicts."""
846 stripped = text.lstrip()
847 if stripped.startswith("[") or stripped.startswith("{"):
848 try:
849 payload = json.loads(stripped)
850 except json.JSONDecodeError as exc:
851 raise InputError(f"invalid JSON: {exc}") from exc
852 if isinstance(payload, dict):
853 payload = payload.get("rows", payload.get("data"))
854 if not isinstance(payload, list):
855 raise InputError("JSON input must be an array of objects, or {\"rows\": [...]}")
856 # Refuse rather than truncate: silently dropping validation data would
857 # produce a clean-looking result computed on part of the study.
858 if len(payload) > MAX_ROWS:
859 raise InputError(f"more than {MAX_ROWS} rows")
860 rows = []
861 for item in payload:
862 if not isinstance(item, dict):
863 raise InputError("JSON rows must be objects")
864 # Match CSV/TSV whitespace handling while retaining null as empty.
865 rows.append({
866 str(k).strip(): "" if v is None else str(v).strip()
867 for k, v in item.items()
868 })
869 if not rows:
870 raise InputError("no data rows found")
871 return rows
873 delimiter = "\t" if (path_hint or "").endswith((".tsv", ".tab")) else None
874 if delimiter is None:
875 first = text.splitlines()[0] if text.splitlines() else ""
876 delimiter = "\t" if first.count("\t") > first.count(",") else ","
877 reader = csv.DictReader(io.StringIO(text), delimiter=delimiter)
878 rows = []
879 for i, row in enumerate(reader):
880 if i >= MAX_ROWS:
881 raise InputError(f"more than {MAX_ROWS} rows")
882 # DictReader stores surplus fields as a list under its None restkey.
883 if None in row:
884 raise InputError(
885 f"row {i + 1}: more fields than the header ({len(row[None])} extra)"
886 )
887 rows.append({(k or "").strip(): (v or "").strip() for k, v in row.items()})
888 if not rows:
889 raise InputError("no data rows found")
890 return rows
893def require_columns(rows: list[dict[str, str]], columns: Iterable[str]) -> None:
894 columns = tuple(columns)
895 for i, row in enumerate(rows):
896 missing = [c for c in columns if c not in row]
897 if missing:
898 raise InputError(f"row {i + 1}: missing required column(s): {', '.join(missing)}")
901def to_float(value: str, column: str, row_index: int) -> float:
902 try:
903 number = float(str(value).strip())
904 if not math.isfinite(number):
905 raise ValueError("non-finite number")
906 return number
907 except (TypeError, ValueError) as exc:
908 raise InputError(
909 f"row {row_index + 1}: column '{column}' is not numeric (finite required): {value!r}"
910 ) from exc
913def fmt(value: Any, digits: int = 4) -> str:
914 """Format a number for a table cell."""
915 if value is None:
916 return ""
917 if isinstance(value, bool):
918 return "yes" if value else "no"
919 if isinstance(value, float):
920 if math.isnan(value):
921 return "n/a"
922 if math.isinf(value):
923 return "inf"
924 if value != 0 and (abs(value) < 1e-4 or abs(value) >= 1e6):
925 return f"{value:.{digits}e}"
926 return f"{value:.{digits}f}"
927 return str(value)
930def emit_table(rows: list[dict[str, Any]], stream=None) -> None:
931 """Print aligned columns."""
932 stream = stream or sys.stdout
933 if not rows:
934 print("(no rows)", file=stream)
935 return
936 headers = list(rows[0].keys())
937 cells = [[fmt(r.get(h)) for h in headers] for r in rows]
938 widths = [
939 max(len(h), *(len(c[i]) for c in cells)) if cells else len(h)
940 for i, h in enumerate(headers)
941 ]
942 print(" ".join(h.ljust(w) for h, w in zip(headers, widths)).rstrip(), file=stream)
943 for c in cells:
944 print(" ".join(v.ljust(w) for v, w in zip(c, widths)).rstrip(), file=stream)
947def emit(rows: list[dict[str, Any]], fmt_name: str, stream=None) -> None:
948 """Print rows as a table, TSV, or JSON."""
949 stream = stream or sys.stdout
950 if fmt_name == "json":
951 json.dump(_json_safe(rows), stream, indent=2, allow_nan=False)
952 print(file=stream)
953 elif fmt_name == "tsv":
954 if not rows:
955 return
956 headers = list(rows[0].keys())
957 print("\t".join(headers), file=stream)
958 for r in rows:
959 print("\t".join(fmt(r.get(h)) for h in headers), file=stream)
960 else:
961 emit_table(rows, stream=stream)
964def _json_safe(obj: Any) -> Any:
965 if isinstance(obj, float) and (math.isnan(obj) or math.isinf(obj)):
966 return None
967 if isinstance(obj, (tuple, list)):
968 return [_json_safe(value) for value in obj]
969 if isinstance(obj, dict):
970 return {key: _json_safe(value) for key, value in obj.items()}
971 return obj
974def note(message: str) -> None:
975 """Write provenance and caveats to stderr so stdout stays parseable."""
976 print(f"note: {message}", file=sys.stderr)
979def finding(message: str) -> None:
980 print(f"finding: {message}", file=sys.stderr)
983def add_common_args(parser: argparse.ArgumentParser) -> None:
984 parser.add_argument(
985 "--format",
986 choices=("table", "tsv", "json"),
987 default="table",
988 help="output format (default: table)",
989 )
992def run_cli(main_func) -> None:
993 """Wrap a main() so InputError becomes a clean exit code 2."""
994 try:
995 sys.exit(main_func())
996 except InputError as exc:
997 print(f"error: {exc}", file=sys.stderr)
998 sys.exit(EXIT_INPUT_ERROR)
999 except BrokenPipeError:
1000 sys.exit(EXIT_OK)