Coverage for src / tracekit / core / uncertainty.py: 90%

1"""Measurement uncertainty propagation and estimation. 

2 

3This module provides data structures and utilities for tracking measurement 

4uncertainty through TraceKit analysis operations, following GUM (Guide to the 

5Expression of Uncertainty in Measurement) principles. 

6 

7References: 

8 JCGM 100:2008 - Guide to the Expression of Uncertainty in Measurement (GUM) 

9 ISO/IEC 17025:2017 - General Requirements for Competence of Testing and Calibration Laboratories 

10 NIST Technical Note 1297 - Guidelines for Evaluating and Expressing Uncertainty 

11""" 

12 

13from __future__ import annotations 

14 

15from dataclasses import dataclass 

16from typing import TYPE_CHECKING, Any 

17 

18import numpy as np 

19 

20if TYPE_CHECKING: 

21 from numpy.typing import NDArray 

22 

23 

24@dataclass 

25class MeasurementWithUncertainty: 

26 """Measurement value with uncertainty estimate. 

27 

28 Stores a measurement value along with its associated uncertainty, 

29 following GUM (Guide to the Expression of Uncertainty in Measurement) 

30 conventions. 

31 

32 Attributes: 

33 value: The measured value (best estimate). 

34 uncertainty: Standard uncertainty (1-sigma, coverage factor k=1). 

35 unit: Unit of measurement (e.g., "V", "s", "Hz") (optional). 

36 coverage_factor: Coverage factor for expanded uncertainty (default 1.0). 

37 confidence_level: Confidence level for expanded uncertainty (default 68.27% for k=1). 

38 n_samples: Number of samples used in measurement (optional). 

39 degrees_of_freedom: Effective degrees of freedom (optional). 

40 

41 Properties: 

42 expanded_uncertainty: Uncertainty * coverage_factor. 

43 relative_uncertainty: uncertainty / |value| if value != 0. 

44 lower_bound: value - expanded_uncertainty. 

45 upper_bound: value + expanded_uncertainty. 

46 

47 Example: 

48 >>> result = MeasurementWithUncertainty(value=1.234, uncertainty=0.005, unit="V") 

49 >>> print(f"{result.value:.3f} ± {result.uncertainty:.3f} {result.unit}") 

50 1.234 ± 0.005 V 

51 

52 Example with expanded uncertainty (k=2, 95.45% confidence): 

53 >>> result = MeasurementWithUncertainty( 

54 ... value=10.0e6, 

55 ... uncertainty=1000.0, 

56 ... unit="Hz", 

57 ... coverage_factor=2.0, 

58 ... confidence_level=0.9545 

59 ... ) 

60 >>> print(f"{result.value/1e6:.3f} ± {result.expanded_uncertainty/1e6:.3f} MHz (95.45%)") 

61 10.000 ± 0.002 MHz (95.45%) 

62 

63 References: 

64 JCGM 100:2008 Section 2.3.5 (standard uncertainty) 

65 JCGM 100:2008 Section 2.3.6 (expanded uncertainty) 

66 """ 

67 

68 value: float 

69 uncertainty: float 

70 unit: str | None = None 

71 coverage_factor: float = 1.0 

72 confidence_level: float = 0.6827 # 68.27% for k=1 (Gaussian) 

73 n_samples: int | None = None 

74 degrees_of_freedom: float | None = None 

75 

76 def __post_init__(self) -> None: 

77 """Validate measurement result after initialization.""" 

78 if not np.isfinite(self.value): 

79 # Allow NaN/Inf but issue warning 

80 pass 

81 if self.uncertainty < 0: 

82 raise ValueError(f"uncertainty must be non-negative, got {self.uncertainty}") 

83 if self.coverage_factor <= 0: 

84 raise ValueError(f"coverage_factor must be positive, got {self.coverage_factor}") 

85 if not 0 < self.confidence_level <= 1.0: 

86 raise ValueError(f"confidence_level must be in (0, 1], got {self.confidence_level}") 

87 

88 @property 

89 def expanded_uncertainty(self) -> float: 

90 """Expanded uncertainty (U = k * u). 

91 

92 Returns: 

93 Standard uncertainty multiplied by coverage factor. 

94 

95 References: 

96 JCGM 100:2008 Section 6.2.1 

97 """ 

98 return self.coverage_factor * self.uncertainty 

99 

100 @property 

101 def relative_uncertainty(self) -> float: 

102 """Relative standard uncertainty (u_r = u / |value|). 

103 

104 Returns: 

105 Relative uncertainty, or np.inf if value is zero. 

106 

107 References: 

108 JCGM 100:2008 Section 5.1.6 

109 """ 

110 if self.value == 0: 

111 return np.inf 

112 return abs(self.uncertainty / self.value) 

113 

114 @property 

115 def lower_bound(self) -> float: 

116 """Lower bound of uncertainty interval. 

117 

118 Returns: 

119 value - expanded_uncertainty. 

120 """ 

121 return self.value - self.expanded_uncertainty 

122 

123 @property 

124 def upper_bound(self) -> float: 

125 """Upper bound of uncertainty interval. 

126 

127 Returns: 

128 value + expanded_uncertainty. 

129 """ 

130 return self.value + self.expanded_uncertainty 

131 

132 def __str__(self) -> str: 

133 """String representation of measurement with uncertainty.""" 

134 unit_str = f" {self.unit}" if self.unit else "" 

135 if self.coverage_factor == 1.0: 135 ↛ 138line 135 didn't jump to line 138 because the condition on line 135 was always true

136 return f"{self.value:.6g} ± {self.uncertainty:.6g}{unit_str}" 

137 else: 

138 return ( 

139 f"{self.value:.6g} ± {self.expanded_uncertainty:.6g}{unit_str} " 

140 f"(k={self.coverage_factor:.2f}, {self.confidence_level * 100:.2f}%)" 

141 ) 

142 

143 def __repr__(self) -> str: 

144 """Detailed representation.""" 

145 return ( 

146 f"MeasurementWithUncertainty(value={self.value:.6g}, " 

147 f"uncertainty={self.uncertainty:.6g}, " 

148 f"unit={self.unit!r})" 

149 ) 

150 

151 

152class UncertaintyEstimator: 

153 """Utilities for estimating measurement uncertainty. 

154 

155 Provides methods for calculating Type A (statistical) and Type B 

156 (systematic) uncertainties according to GUM principles. 

157 

158 References: 

159 JCGM 100:2008 Section 4 (Evaluating Standard Uncertainty) 

160 """ 

161 

162 @staticmethod 

163 def type_a_standard_deviation(data: NDArray[np.floating[Any]]) -> float: 

164 """Type A uncertainty from sample standard deviation. 

165 

166 Args: 

167 data: Array of repeated measurements. 

168 

169 Returns: 

170 Standard deviation (Type A uncertainty component). 

171 

172 References: 

173 JCGM 100:2008 Section 4.2 (Type A evaluation) 

174 """ 

175 if len(data) < 2: 

176 return np.nan 

177 return float(np.std(data, ddof=1)) # Sample std (Bessel correction) 

178 

179 @staticmethod 

180 def type_a_standard_error(data: NDArray[np.floating[Any]]) -> float: 

181 """Type A uncertainty from standard error of the mean. 

182 

183 Args: 

184 data: Array of repeated measurements. 

185 

186 Returns: 

187 Standard error (sigma / sqrt(n)), Type A uncertainty of the mean. 

188 

189 References: 

190 JCGM 100:2008 Section 4.2.3 

191 """ 

192 if len(data) < 2: 192 ↛ 193line 192 didn't jump to line 193 because the condition on line 192 was never true

193 return np.nan 

194 return float(np.std(data, ddof=1) / np.sqrt(len(data))) 

195 

196 @staticmethod 

197 def combined_uncertainty( 

198 uncertainties: list[float], correlation_matrix: NDArray[np.float64] | None = None 

199 ) -> float: 

200 """Combine multiple uncertainty sources. 

201 

202 Combines uncorrelated or correlated uncertainty components using 

203 the law of propagation of uncertainty. 

204 

205 Args: 

206 uncertainties: List of standard uncertainties to combine. 

207 correlation_matrix: Correlation matrix for correlated inputs (optional). 

208 If None, assumes all inputs are uncorrelated. 

209 

210 Returns: 

211 Combined standard uncertainty. 

212 

213 Example (uncorrelated): 

214 >>> u1, u2, u3 = 0.01, 0.02, 0.005 

215 >>> u_combined = UncertaintyEstimator.combined_uncertainty([u1, u2, u3]) 

216 >>> print(f"Combined: {u_combined:.4f}") 

217 Combined: 0.0233 

218 

219 Example (correlated): 

220 >>> import numpy as np 

221 >>> u = [0.01, 0.02] 

222 >>> R = np.array([[1.0, 0.5], [0.5, 1.0]]) # 50% correlation 

223 >>> u_combined = UncertaintyEstimator.combined_uncertainty(u, R) 

224 

225 References: 

226 JCGM 100:2008 Section 5.1.2 (uncorrelated) 

227 JCGM 100:2008 Section 5.2.2 (correlated) 

228 """ 

229 u_array = np.array(uncertainties, dtype=np.float64) 

230 

231 if correlation_matrix is None: 

232 # Uncorrelated: u_c² = Σ u_i² 

233 return float(np.sqrt(np.sum(u_array**2))) 

234 else: 

235 # Correlated: u_c² = u^T R u 

236 u_combined_sq = u_array @ correlation_matrix @ u_array 

237 return float(np.sqrt(u_combined_sq)) 

238 

239 @staticmethod 

240 def type_b_rectangular(half_width: float) -> float: 

241 """Type B uncertainty from rectangular (uniform) distribution. 

242 

243 Used when only min/max bounds are known with equal probability. 

244 

245 Args: 

246 half_width: Half-width of the interval (a). 

247 

248 Returns: 

249 Standard uncertainty u = a / √3. 

250 

251 Example: 

252 >>> # Scope resolution ±0.5 LSB 

253 >>> lsb = 1.0e-3 # 1 mV per bit 

254 >>> u = UncertaintyEstimator.type_b_rectangular(0.5 * lsb) 

255 >>> print(f"Quantization uncertainty: {u*1e6:.3f} µV") 

256 Quantization uncertainty: 288.675 µV 

257 

258 References: 

259 JCGM 100:2008 Section 4.3.7 (rectangular distribution) 

260 """ 

261 return float(half_width / np.sqrt(3)) 

262 

263 @staticmethod 

264 def type_b_triangular(half_width: float) -> float: 

265 """Type B uncertainty from triangular distribution. 

266 

267 Used when values near the center are more likely than extremes. 

268 

269 Args: 

270 half_width: Half-width of the interval (a). 

271 

272 Returns: 

273 Standard uncertainty u = a / √6. 

274 

275 References: 

276 JCGM 100:2008 Section 4.3.9 (triangular distribution) 

277 """ 

278 return float(half_width / np.sqrt(6)) 

279 

280 @staticmethod 

281 def type_b_from_tolerance(tolerance: float, confidence: float = 0.95) -> float: 

282 """Type B uncertainty from manufacturer tolerance specification. 

283 

284 Assumes Gaussian distribution unless otherwise specified. 

285 

286 Args: 

287 tolerance: Tolerance limit (e.g., ±1% of reading). 

288 confidence: Confidence level of the tolerance (default 95%). 

289 

290 Returns: 

291 Standard uncertainty. 

292 

293 Example: 

294 >>> # Scope vertical accuracy: ±2% of reading, 95% confidence 

295 >>> reading = 1.0 # 1V 

296 >>> tolerance = 0.02 * reading # ±0.02 V 

297 >>> u = UncertaintyEstimator.type_b_from_tolerance(tolerance, 0.95) 

298 >>> print(f"Uncertainty: {u:.4f} V") 

299 Uncertainty: 0.0102 V 

300 

301 References: 

302 JCGM 100:2008 Section 4.3.4 

303 """ 

304 # For Gaussian, 95% confidence → k ≈ 1.96 

305 # For 99% confidence → k ≈ 2.58 

306 if confidence == 0.95: 

307 k = 1.96 

308 elif confidence == 0.99: 308 ↛ 310line 308 didn't jump to line 310 because the condition on line 308 was always true

309 k = 2.58 

310 elif confidence == 0.6827: 

311 k = 1.0 

312 else: 

313 # General case: approximate using normal distribution 

314 from scipy import stats 

315 

316 k = stats.norm.ppf((1 + confidence) / 2) 

317 

318 return float(tolerance / k) 

319 

320 @staticmethod 

321 def time_base_uncertainty(sample_rate: float, timebase_accuracy_ppm: float = 50.0) -> float: 

322 """Calculate time base uncertainty from oscilloscope specification. 

323 

324 Args: 

325 sample_rate: Sample rate in Hz. 

326 timebase_accuracy_ppm: Timebase accuracy in parts per million (typical: 10-50 ppm). 

327 

328 Returns: 

329 Standard uncertainty in time per sample (seconds). 

330 

331 Example: 

332 >>> # 1 GSa/s scope, 25 ppm timebase (typical for OCXO) 

333 >>> u_t = UncertaintyEstimator.time_base_uncertainty(1e9, 25.0) 

334 >>> print(f"Timebase uncertainty: {u_t*1e12:.2f} ps per sample") 

335 Timebase uncertainty: 25.00 ps per sample 

336 

337 References: 

338 IEEE 181-2011 Annex B (Measurement Uncertainty) 

339 """ 

340 time_per_sample = 1.0 / sample_rate 

341 # ppm = parts per million = 1e-6 

342 relative_uncertainty = timebase_accuracy_ppm * 1e-6 

343 return time_per_sample * relative_uncertainty 

344 

345 @staticmethod 

346 def vertical_uncertainty( 

347 reading: float, 

348 vertical_accuracy_percent: float = 2.0, 

349 offset_error_volts: float = 0.0, 

350 ) -> float: 

351 """Calculate vertical (voltage) uncertainty from oscilloscope specification. 

352 

353 Typical scope spec: ±(2% of reading + 0.1% of full scale + 1 mV) 

354 

355 Args: 

356 reading: Measured voltage value. 

357 vertical_accuracy_percent: Gain accuracy in percent (e.g., 2.0 for ±2%). 

358 offset_error_volts: Fixed offset error in volts. 

359 

360 Returns: 

361 Standard uncertainty in volts. 

362 

363 Example: 

364 >>> # Tektronix TDS3000 series: ±3% of reading 

365 >>> reading = 1.5 # 1.5 V 

366 >>> u_v = UncertaintyEstimator.vertical_uncertainty(reading, 3.0, 0.001) 

367 >>> print(f"Voltage uncertainty: {u_v*1000:.2f} mV") 

368 Voltage uncertainty: 45.01 mV 

369 

370 References: 

371 IEEE 1057-2017 Section 4.4 (Amplitude Accuracy) 

372 """ 

373 # Gain error 

374 u_gain = abs(reading) * (vertical_accuracy_percent / 100.0) 

375 

376 # Offset error 

377 u_offset = offset_error_volts 

378 

379 # Combine (assume uncorrelated) 

380 return float(np.sqrt(u_gain**2 + u_offset**2)) 

381 

382 

383__all__ = ["MeasurementWithUncertainty", "UncertaintyEstimator"]