Skills › Scientific Agent Skills › uncertainty-and-units › uncertainty_budget.py
uncertainty_budget.py
373 lines · skills/uncertainty-and-units/scripts/uncertainty_budget.py at 5f2bc03
1#!/usr/bin/env python3
2"""Turn a list of uncertainty components into a GUM uncertainty budget.
4Each component arrives the way it is actually stated on a certificate, a data
5sheet, or a repeatability worksheet. The divisor that converts it to a standard
6uncertainty depends on which of those it is, and getting that divisor wrong is
7the most common defect in a real budget. This CLI applies JCGM 100:2008 4.3,
8combines the components, and reports the Welch-Satterthwaite effective degrees
9of freedom and the resulting coverage factor.
10"""
12from __future__ import annotations
14import argparse
15import math
16import sys
17from typing import Any
19import _common
20from _common import CliError
23TEMPLATE: dict[str, Any] = {
24 "measurand": "example flow rate",
25 "unit": "L/min",
26 "value": 12.345,
27 "coverage_probability": 0.95,
28 "components": [
29 {
30 "label": "repeatability of ten readings",
31 "type": "A",
32 "distribution": "normal",
33 "value": 0.052,
34 "dof": 9,
35 "note": "experimental standard deviation of the mean",
36 },
37 {
38 "label": "calibration certificate",
39 "type": "B",
40 "distribution": "expanded",
41 "value": 0.10,
42 "coverage_factor": 2.0,
43 },
44 {
45 "label": "display resolution",
46 "type": "B",
47 "distribution": "rectangular",
48 "value": 0.005,
49 "note": "half-width equals half of the last displayed digit",
50 },
51 {
52 "label": "long-term drift since calibration",
53 "type": "B",
54 "distribution": "rectangular",
55 "value": 0.04,
56 },
57 {
58 "label": "temperature correction",
59 "type": "B",
60 "distribution": "triangular",
61 "value": 0.03,
62 "sensitivity": 0.8,
63 },
64 ],
65}
68def normalize_component(raw: Any, index: int, measurand: float | None) -> dict[str, Any]:
69 """Validate one budget component and reduce it to a standard uncertainty."""
71 if not isinstance(raw, dict):
72 raise CliError(f"component {index + 1} must be a JSON object")
73 label = raw.get("label") or f"component {index + 1}"
74 if not isinstance(label, str) or len(label) > 200:
75 raise CliError(f"component {index + 1}: label must be a string under 200 chars")
77 kind = str(raw.get("type", "B")).strip().upper()
78 if kind not in {"A", "B"}:
79 raise CliError(f"{label}: evaluation type must be 'A' or 'B'")
81 shape = raw.get("distribution")
82 shape = "expanded" if str(shape).strip().lower() == "expanded" else shape
83 if shape == "expanded":
84 distribution = "expanded"
85 else:
86 distribution = _common.normalize_distribution(
87 shape, label=f"{label}.distribution"
88 )
90 value = _common.as_finite(raw.get("value"), label=f"{label}.value")
91 if value < 0:
92 raise CliError(f"{label}: value must not be negative")
94 relative = raw.get("relative", False)
95 if not isinstance(relative, bool):
96 raise CliError(f"{label}: relative must be a JSON boolean")
97 if relative:
98 if measurand is None:
99 raise CliError(
100 f"{label}: a relative component needs a top-level measurand 'value'"
101 )
102 value = value * abs(measurand)
104 if distribution == "expanded":
105 factor = _common.as_finite(
106 raw.get("coverage_factor"), label=f"{label}.coverage_factor"
107 )
108 if factor <= 0:
109 raise CliError(f"{label}: coverage_factor must be greater than zero")
110 divisor = factor
111 else:
112 divisor = _common.DISTRIBUTION_DIVISORS[distribution]
114 if "divisor" in raw:
115 divisor = _common.as_finite(raw["divisor"], label=f"{label}.divisor")
116 if divisor <= 0:
117 raise CliError(f"{label}: divisor must be greater than zero")
119 sensitivity = _common.as_finite(
120 raw.get("sensitivity", 1.0), label=f"{label}.sensitivity"
121 )
122 dof = _common.as_degrees_of_freedom(raw.get("dof"), label=f"{label}.dof")
123 standard = value / divisor
124 if distribution == "exact" and value != 0:
125 raise CliError(f"{label}: an exact component must have zero uncertainty")
126 note = raw.get("note")
127 if note is not None and (not isinstance(note, str) or len(note) > 500):
128 raise CliError(f"{label}: note must be a string under 500 chars")
130 return {
131 "label": label,
132 "type": kind,
133 "distribution": distribution,
134 "stated_value": value,
135 "divisor": divisor,
136 "standard_uncertainty": standard,
137 "sensitivity": sensitivity,
138 "contribution": sensitivity * standard,
139 "dof": dof,
140 "note": note,
141 }
144def build_budget(spec: dict[str, Any], coverage_override: float | None) -> dict[str, Any]:
145 """Combine components into u_c, effective degrees of freedom, and U."""
147 if spec.get("correlations") or spec.get("covariance") or spec.get("covariance_matrix"):
148 raise CliError("this budget combines independent components only; use a covariance-aware model")
149 raw_components = spec.get("components")
150 if not isinstance(raw_components, list) or not raw_components:
151 raise CliError("the spec must contain a non-empty 'components' list")
152 if len(raw_components) > _common.MAX_COMPONENTS:
153 raise CliError(f"at most {_common.MAX_COMPONENTS} components are supported")
155 measurand_value = spec.get("value")
156 if measurand_value is not None:
157 measurand_value = _common.as_finite(measurand_value, label="value")
159 components = [
160 normalize_component(raw, index, measurand_value)
161 for index, raw in enumerate(raw_components)
162 ]
164 variance = sum(item["contribution"] ** 2 for item in components)
165 combined = math.sqrt(variance)
166 for item in components:
167 item["variance_fraction"] = (
168 (item["contribution"] ** 2) / variance if variance > 0 else 0.0
169 )
171 coverage = coverage_override
172 if coverage is None:
173 coverage = _common.as_finite(
174 spec.get("coverage_probability", 0.95), label="coverage_probability"
175 )
176 if not 0 < coverage < 1:
177 raise CliError("coverage probability must be strictly between 0 and 1")
179 effective_dof = _common.welch_satterthwaite(
180 combined, [(abs(item["contribution"]), item["dof"]) for item in components]
181 )
182 factor = _common.coverage_factor(effective_dof, coverage)
183 expanded = factor * combined
185 result: dict[str, Any] = {
186 "measurand": spec.get("measurand"),
187 "unit": spec.get("unit"),
188 "value": measurand_value,
189 "components": components,
190 "combined_standard_uncertainty": combined,
191 "effective_degrees_of_freedom": (
192 None if not math.isfinite(effective_dof) else effective_dof
193 ),
194 "coverage_probability": coverage,
195 "coverage_factor": factor,
196 "expanded_uncertainty": expanded,
197 }
198 if measurand_value is not None:
199 result["coverage_interval"] = [
200 measurand_value - expanded,
201 measurand_value + expanded,
202 ]
203 if measurand_value != 0:
204 result["relative_standard_uncertainty"] = combined / abs(measurand_value)
205 result["warnings"] = collect_warnings(result)
206 for item in components:
207 item["dof"] = _common.json_dof(item["dof"])
208 return result
211def collect_warnings(budget: dict[str, Any]) -> list[str]:
212 """Flag budget defects that change how the result should be reported."""
214 warnings: list[str] = []
215 warnings.append(
216 "the normal/Student-t coverage approximation requires a suitable output "
217 "distribution; this budget does not verify coverage or measurement-model validity"
218 )
219 components = budget["components"]
220 for item in components:
221 if item["type"] == "A" and not math.isfinite(item["dof"]):
222 warnings.append(
223 f"{item['label']}: a Type A component evaluated from n readings has "
224 "n-1 degrees of freedom; leaving it infinite understates the "
225 "coverage factor"
226 )
227 if (
228 item["type"] == "B"
229 and item["distribution"] == "normal"
230 and item["divisor"] == 1.0
231 ):
232 warnings.append(
233 f"{item['label']}: a Type B component declared normal is divided by "
234 "1, so the stated value is taken as a standard uncertainty; "
235 "certificates normally quote an expanded U, which needs "
236 "distribution 'expanded' with its coverage_factor"
237 )
238 dof = budget["effective_degrees_of_freedom"]
239 if dof is not None and dof < 6:
240 warnings.append(
241 f"effective degrees of freedom is {dof:.1f}; k = "
242 f"{budget['coverage_factor']:.2f} rather than the customary 2, and the "
243 "interval is sensitive to one small sample"
244 )
245 if components:
246 dominant = max(components, key=lambda item: item["variance_fraction"])
247 if dominant["variance_fraction"] > 0.9:
248 warnings.append(
249 f"{dominant['label']} contributes "
250 f"{dominant['variance_fraction'] * 100:.1f}% of the variance; "
251 "improving other components has limited effect on the result"
252 )
253 largest = max(abs(item["contribution"]) for item in components)
254 negligible = [
255 item["label"]
256 for item in components
257 if largest > 0 and abs(item["contribution"]) < largest / 3.0
258 ]
259 if negligible:
260 warnings.append(
261 "each of these components is below one third of the largest; "
262 f"their aggregate effect may still matter: {', '.join(negligible)}"
263 )
264 return warnings
267def render_markdown(budget: dict[str, Any]) -> str:
268 """Render the budget as a Markdown table plus the combined result."""
270 unit = f" {budget['unit']}" if budget.get("unit") else ""
271 lines = [
272 f"# Uncertainty budget: {budget.get('measurand') or 'measurand'}",
273 "",
274 "| Component | Type | Distribution | Stated | Divisor | u(x) | c | "
275 "c*u(x) | % variance | dof |",
276 "| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |",
277 ]
278 for item in budget["components"]:
279 dof = _common.format_dof(item["dof"])
280 lines.append(
281 f"| {item['label']} | {item['type']} | {item['distribution']} | "
282 f"{item['stated_value']:.6g} | {item['divisor']:.4g} | "
283 f"{item['standard_uncertainty']:.4g} | {item['sensitivity']:.4g} | "
284 f"{item['contribution']:.4g} | {item['variance_fraction'] * 100:.1f} | "
285 f"{dof} |"
286 )
287 effective = budget["effective_degrees_of_freedom"]
288 lines += [
289 "",
290 "## Combined result",
291 "",
292 f"- Combined standard uncertainty u_c: "
293 f"{budget['combined_standard_uncertainty']:.6g}{unit}",
294 f"- Effective degrees of freedom: "
295 f"{'infinite' if effective is None else f'{effective:.1f}'}",
296 f"- Coverage factor k: {budget['coverage_factor']:.4g} at p = "
297 f"{budget['coverage_probability']:.3g}",
298 f"- Expanded uncertainty U: {budget['expanded_uncertainty']:.6g}{unit}",
299 ]
300 if budget.get("value") is not None:
301 low, high = budget["coverage_interval"]
302 lines.append(
303 f"- Result: {budget['value']:.6g} +/- "
304 f"{budget['expanded_uncertainty']:.6g}{unit} "
305 f"[{low:.6g}, {high:.6g}]"
306 )
307 if budget.get("warnings"):
308 lines += ["", "## Warnings", ""]
309 lines += [f"- {message}" for message in budget["warnings"]]
310 lines += [
311 "",
312 "Components are combined in quadrature, which assumes they are "
313 "uncorrelated (JCGM 100:2008 equation 10). Use "
314 "`propagate_uncertainty.py --correlation` when they are not.",
315 ]
316 return "\n".join(lines) + "\n"
319def build_parser() -> argparse.ArgumentParser:
320 parser = argparse.ArgumentParser(
321 description=(
322 "Combine stated uncertainty components into a GUM budget with "
323 "distribution-aware divisors and a Welch-Satterthwaite coverage factor."
324 )
325 )
326 parser.add_argument("--spec", help="JSON file describing the budget")
327 parser.add_argument(
328 "--template",
329 action="store_true",
330 help="print a worked example spec and exit",
331 )
332 parser.add_argument(
333 "--coverage",
334 type=_common.probability,
335 help="override the coverage probability in the spec",
336 )
337 parser.add_argument("--format", choices=("json", "markdown"), default="json")
338 parser.add_argument("--output", help="write the report to this file")
339 parser.add_argument(
340 "--force", action="store_true", help="overwrite an existing output file"
341 )
342 return parser
345def main(argv: list[str] | None = None) -> int:
346 arguments = build_parser().parse_args(argv)
347 try:
348 if arguments.template:
349 _common.emit_json(
350 TEMPLATE, output=arguments.output, force=arguments.force
351 )
352 return 0
353 spec = TEMPLATE
354 if arguments.spec:
355 loaded = _common.load_json(arguments.spec)
356 if not isinstance(loaded, dict):
357 raise CliError("the spec file must hold a JSON object")
358 spec = loaded
359 budget = build_budget(spec, arguments.coverage)
360 if arguments.format == "markdown":
361 _common.emit_text(
362 render_markdown(budget), output=arguments.output, force=arguments.force
363 )
364 else:
365 _common.emit_json(budget, output=arguments.output, force=arguments.force)
366 except CliError as error:
367 print(f"error: {error}", file=sys.stderr)
368 return 2
369 return 0
372if __name__ == "__main__":
373 raise SystemExit(main())