Version: v3.8.1 (2026-05) · Scope: T0 fitting → Mass Ratio → Air Ratio → J → Age → Datum Publication
v3.8.0 update: §10 (DiagramPlot) rewritten to reflect corrected inverse-isochron F = −b/a, inverse-variance-weighted WMA, and WMA-referenced MSWD per Vermeesch (2024, 2018) and Schaen et al. (2021). Previous §10 description matched the buggy v3.7.x code (§10.x note retains the bug history for traceability).
v3.8.1 update (2026-05): Documentation-only changes — no math code modified. (1) §1 T0 sigma now properly flagged as math bug (was previously described as "NTNU lab convention" — too soft); cross-references the correct intercept-SE formula and AutoPipeline
_fit_onev3.8.2 patch. (2) §6 J-Volume now flags the σ_J bracket bug atUtilities.py:2864(operator precedence drops a critical pair). (3) §10.7 flags that the plateau "avg_age" still uses arithmeticnp.mean, not the v3.8-corrected WMA in §10.5. (4) §10.8 adds a warning about isochron-intercept SE computed by running OLS on the error bars (curve_fit(linear, x_std, y_std)) atUtilities.py:389-390, 965-966— mathematically meaningless. See §⚠ for the consolidated outstanding-bug list synced with the v5 math-audit report.
This document is a one-stop reference for every numeric formula pyADR evaluates, written as math (LaTeX) plus the corresponding .py line. Error propagation is derived from first-order partials of each expression as actually coded — pyADR mixes proper quadrature (sqrt(Σ(∂F/∂x · σ_x)²)) with linear-sum approximations (F·Σ|σ_x/x|); both are documented faithfully. Where the code's choice diverges from the McDougall & Harrison (1999) / Koppers (2002) convention, a Note flags it.
Symbol convention:
mᵢ= measured intensity of ⁴ⁱAr (already T0-extrapolated, blank-corrected, decay-corrected)aᵢ/kᵢ/caᵢ/clᵢ/rᵢ= atmospheric / K-derived / Ca-derived / Cl-derived / radiogenic componentσ_x= 1σ uncertainty onxλ_total= total ⁴⁰K decay constant (parameterconstants[14], /yr)
Utilities.py:124–128
Ratio (proper quadrature):
∂(y/x)/∂y = 1/x, ∂(y/x)/∂x = −y/x². The implementation expects callers to pass the ratio value as the 5th arg.
Difference / sum (proper quadrature):
∂(y±x)/∂y = ±1, ∂(y±x)/∂x = ±1.
Note: Several places in pyADR (e.g.
getJVolumeStatistics,calcAge,calculateSlatCa) use the linear-sum approximation|F|·(σ_a/|a| + σ_b/|b|)instead of quadrature. This is conservative (slightly overestimates σ when terms are independent) and is preserved here as-is for backward compatibility with NTNU lab outputs.
Utilities.py:142–305
For each isotope i ∈ {36,37,38,39,40} the raw voltage trace (t, V) is fit with one of:
Linear (fit_func_list[0], Utilities.py:132):
Average (fit_func_list[1], Utilities.py:135):
T0 sigma (as coded — Utilities.py:187, 216, 275, 294):
where n_out is the number of cycles masked out by the outlier rule
(at most 4 cycles removed; up to 1 retry).
Goodness: R² from sklearn.metrics.r2_score.
⚠ Math bug — σ_T0 is not the intercept SE. The coded formula above is the standard error of the mean (
std(residuals)/√N), not the SE of the y-intercept of a linear fit. For typical step-heating cycle timest ∈ [320, 600] switht̄ ≈ 460 s≫ 0, the intercept extrapolated tot=0carries a large lever-arm term that is completely missing here.Correct intercept SE (Li et al. 2019 Eq. 1):
$$ \sigma_{T_0,i}^{(\text{correct})} ;=; \sigma_r \sqrt{\frac{1}{N} + \frac{\bar t^{,2}}{\sum_j (t_j - \bar t)^2}} ;=; \sqrt{\mathrm{pcov}[-1,-1]}\quad\text{from}\quad\mathrm{curve_fit}(f, t, V) $$ where
σ_r = std(residuals, ddof=N−p)is the unbiased residual SD andpcovis the regression covariance matrix.Empirical impact (NTNU 0621-01C 1100 °C step, 10 cycles, t ∈ [320, 601] s):
Isotope T0 σ_coded ( std/√N)σ_correct ( pcov)ratio ³⁶Ar 3.58 × 10⁻⁴ 4.6 × 10⁻⁶ 4.6 × 10⁻⁵ 10.1× ³⁷Ar 1.62 × 10⁻⁴ 5.6 × 10⁻⁶ 5.5 × 10⁻⁵ 9.8× ³⁸Ar 6.44 × 10⁻⁴ 6.9 × 10⁻⁶ 6.9 × 10⁻⁵ 10.0× ³⁹Ar 1.91 × 10⁻² 9.6 × 10⁻⁶ 9.8 × 10⁻⁵ 10.2× ⁴⁰Ar 2.04 × 10⁻² 4.8 × 10⁻⁶ 5.8 × 10⁻⁵ 11.9× The coded form underestimates σ_T0 by ~10× for every isotope. This propagates through
getJVolumeStatistics,calcAge, plane-fit weighting, and downstream WMA / MSWD computations — every reported σ in pyADR currently inherits this underestimate.Status: Fixed in AutoPipeline
_fit_one(v3.8.2 patch — usessqrt(pcov[-1,-1])with Li 2019 closed-form fallback). LegacyUtilities.calculateT0/REcalculateT0(this section) still uses the buggy form pending advisor review; the upstreamAndrewLiu0725/pyADRGitHub repo uses a different but also-incorrect formulastd(|residuals|)(no/√Ndivisor). Reference: Li, X., Naeher, U. and Pross, J. (2019). Mass spectrometric data processing in stable isotope analysis: regression-based standard errors of the y-intercept. J. Mass Spectrom. 54, 145–152, Eq. 1.
Utilities.py:2529–2625
For each isotope i, given N saved T0 files:
Auto-mask rule (getT0Statistics:2550):
After masking, recompute mean/std over the surviving subset → restatistics.
Caveat: the threshold
σ/2 + meanmixes a half-σ scale with the mean itself; this is what the code does. Standard practice would be|x − μ| > kσ. Keep in mind when interpreting outlier rejection.
Utilities.py:2629–2687
Inputs: mass file T0 (5 isotopes, with σ), preline/background T0 (5 isotopes, with σ), days T = SPD − OGD (sample-prep date minus original-gas date).
Step 1 — blank subtraction (:2670–2671):
Step 2 — short-lived isotope decay correction (:2673–2678):
σ scales by the same factor: ∂(m·eᵏᵀ)/∂m = eᵏᵀ.
Step 3 — five reported ratios (:2682–2685, pair_indices at :2627):
| idx | y / x | meaning |
|---|---|---|
| 0 | 39 / 40 | Ar(39ₘ)/Ar(40ₘ) |
| 1 | 36 / 40 | Ar(36ₘ)/Ar(40ₘ) |
| 2 | 39 / 36 | Ar(39ₘ)/Ar(36ₘ) |
| 3 | 40 / 36 | Ar(40ₘ)/Ar(36ₘ) — for atmospheric check |
| 4 | 38 / 36 | Ar(38ₘ)/Ar(36ₘ) — for atmospheric check |
(ratio-quadrature; abs() added in v3.7.x bug-fix.)
Utilities.py:2066–2109
For N air-shot files, extracts the (40/36)ₐ and (38/36)ₐ ratios (rows 4 and 5, col 9), then trims any |40/36| > 313 outlier:
for p ∈ {40/36, 38/36}. No weighted mean here — pure population mean / std.
Utilities.py:2111–2210, 2435+
Two estimators are computed in parallel and both reported:
Arithmetic mean & SE-of-the-mean (used for outlier mask):
Mask points where J_j falls outside [J̄ − σ_SE, J̄ + σ_SE]. After masking, the second pass uses std(J) (no /√N), i.e. switches to population std.
Inverse-variance weighted mean (:2137–2146):
Returned tuple [avg, σ_pop, μ, σ_μ]. σ_J from the J-volume calc is propagated via 1/σ² weights.
Note: pyADR uses
avg(arithmetic) for the outlier filter but reports the inverse-variance pair(μ, σ_μ)separately. Decide explicitly which to forward to AgeCalc.
Utilities.py:2747–2809
Constants used (irradiation parameter array, indexed by constants[i]):
| idx | symbol | description |
|---|---|---|
| 0 | (³⁹Ar/³⁷Ar)_Ca | Ca-derived ³⁹Ar production ratio |
| 1 | σ of [0] | |
| 2 | (³⁶Ar/³⁷Ar)_Ca | |
| 3 | σ of [2] | |
| 4 | (⁴⁰Ar/³⁹Ar)_K | K-derived ⁴⁰Ar production |
| 12 | (⁴⁰Ar/³⁶Ar)_air | atmospheric, default 295.5 or 298.56 |
| 14 | λ_total | total ⁴⁰K decay (/yr), used in T = ln(1+JF)/λ |
| 16 | λ_total /1e6 | same, scaled to Ma |
Ar component partition (:2764–2785):
Standard ratios (:2789–2796):
each with linear-sum σ: σ_G = G·(σ_{40m}/m_{40} + σ_{39m}/m_{39}) etc.
F (radiogenic ⁴⁰Ar per ³⁹Arₖ) propagation (:2796):*
pyADR codes:
This matches the partials of
i.e. the K-corrected F = (⁴⁰* − air − K-interference) / ³⁹ₖ form, with air- and Ca-correction terms.
J value & uncertainty (:2799–2804):
For an irradiation standard of known age t (years):
Hard-coded constants in this function: λ = 5.531e-10 /yr, σ_λ = 0.0135e-10 /yr (Steiger & Jäger 1977 ⁴⁰K → ⁴⁰Ar branch). For Min et al. (2000) values, change l, l_std at Utilities.py:2748–2749.
⚠ Bug — σ_J v₃ term missing a pair of parentheses (
Utilities.py:2864).The code reads:
v3 = F_std**2 * ((np.exp(l*t)) - 1 / Ar_39_K_40_r_ratio**2) ** 2Python operator precedence binds
1 / F²before(e^λt) − …, so this evaluates to
$$ v_3^{(\text{coded})} = F_{\text{std}}^{2} \cdot \bigl(e^{\lambda t} - \tfrac{1}{F_{r/k}^{2}}\bigr)^{2} $$ instead of the intended
$$ v_3^{(\text{correct})} = F_{\text{std}}^{2} \cdot \left(\frac{e^{\lambda t}-1}{F_{r/k}^{2}}\right)^{2}. $$ Fix: one extra pair of parentheses
v3 = F_std**2 * ((np.exp(l*t) - 1) / Ar_39_K_40_r_ratio**2) ** 2Impact: When
λt ≪ 1(typical young samples),e^{λt} ≈ 1so the correct form givesv₃ ≈ (λt/F²)² · σ_F²(small), while the coded form givesv₃ ≈ (1 − 1/F²)² · σ_F²(dominated by the spurious1). σ_J is massively over-estimated for young samples — sign and rough magnitude reversed from the correct error budget. 5-minute fix; deferred only because Utilities.py is on the advisor-owned upstream.
Ca/K (:2807–2808):
The factor 0.52 is the lab calibration (production cross-section for ³⁷Ar from ⁴⁰Ca / ³⁹Ar from ³⁹K). v3.0.1 fixed an inverted Ca/K and a wrong constant (was pr_ratio=0.000377).
Utilities.py:2811–2920 — restored from V3.4.1 archive in v3.7.4 (was truncated mid-function in v3.7.0–v3.7.3 release HEAD).
Returns 59-element list. Key indices: [18,19] Ar_39_K ± σ; [24,25] Ar_40_r ± σ; [36,37] F ± σ; [46,47] T ± σ (years); [48,49] J_int, T_int.
Constants used (extends §6 table):
| idx | symbol | description |
|---|---|---|
| 0,1 | (39/37)_Ca, σ | Ca-derived 39 production |
| 2,3 | (36/37)_Ca, σ | |
| 4,5 | (40/39)_K, σ | K-derived 40 production |
| 6,7 | (38/39)_K, σ | K-derived 38 production |
| 12,13 | (40/36)_air, σ | atmospheric, default 298.56 |
| 14 | λ_total (/yr) | used in degas plot per-step T = ln(1+JF)/λ |
| 16 | λ_total (/yr) | used in calcAge AND getStackPlot total age |
⚠️ Code inconsistency: per-step T ingetDegasPlot:444usesconstants[14], butcalcAge:2904andgetStackPlot:1972useconstants[16]. If PS sets [14] ≠ [16] (e.g. Steiger 5.531e-10 vs Min 5.463e-10), per-step ages and total age use different λ — check PS values match.
Step 1 — Ar component partition (:2858–2890):
³⁸Ar partition (:2876–2880):
⚠️ Misleading variable name: the field labelledAr_38_Airin the AgeCalc table actually contains ³⁸Ar(air) + ³⁸Ar(Cl). The proper split (³⁸Ar_air = ³⁶Ar_air · (³⁸/³⁶)_a, ³⁸Ar_Cl = ³⁸m − ³⁸_air − ³⁸_K) is only done at the Datum Publication stage intoDP:4648–4664. Don't quote the AgeCalc page's "Ar_38_Air" as pure atmospheric ³⁸Ar.
⁴⁰Ar partition (:2884–2891):
Step 2 — F = ⁴⁰Ar/³⁹Arₖ (:2904):*
F uncertainty — code uses the explicit-partial form (NOT the simple ratio quadrature on F):
where, using the calcAge naming (C₁,C₂,C₃,C₄) = (c₁₂, c₂, c₄, c₀):
with linear-sum σ on each (e.g. σ_G = G·(σ_{40m}/m_{40} + σ_{39m}/m_{39})). Same form as §6's σ_F. Note: code defines C3=c₄ but C3 is not used in the σ_F expression above (only C1, C2, C4 appear).
Step 3 — Age equation (:2906):
T uncertainty as coded (:2907):
i.e. partials ∂T/∂J = F/[λ(1+JF)], ∂T/∂F = J/[λ(1+JF)] in quadrature.
T uncertainty per McDougall & Harrison 1999 eq. 4.7 (full):
⚠️ Code drops the σ_λ term: pyADR'sT_stdformula at:2907omits the third partial(∂T/∂λ)²σ_λ². For Steiger & Jäger λ with σ_λ ≈ 0.24% relative, this under-estimates σ_T by ≈ 0.5% relative. For a 50 Ma age with σ_T ≈ 0.5 Ma, the missing contribution is ≈ 0.12 Ma — small but systematic. Flag for advisor discussion before publication.
Internal age uncertainty (:2908):
Same form as σ_T but with σ_J replaced by σ_J^(int) (the v₃-only J uncertainty from §6). σ_λ also dropped here.
Auxiliary ratios returned (:2895–2900):
each with linear-sum σ. Used by inverse-isochron / atmospheric-correction plots.
Utilities.py:2689–2745
CaF₂ salt → (³⁶/³⁷)_Ca and (³⁹/³⁷)_Ca (:2708–2711):
The 298.56 is the assumed atmospheric ⁴⁰/³⁶ used to back-correct the air component; this is lab-fixed (Lee et al. 2006). To use 295.5 (Steiger & Jäger 1977) edit :2708–2738.
K salt → (⁴⁰/³⁹)_K, (³⁸/³⁹)_K, (³⁹/³⁷)_K (:2737–2742):
σ analogous to the CaF case (linear-sum on the numerator's air subtraction).
NTNU_DataReduction.py:4507–4834
For each AgeCalc CSV row, the publication table re-derives:
Step heating fractions (:4636–4637):
(Note: σ not propagated for these — pure denominators.)
Per-step atmospheric ⁴⁰/³⁶ ratio (:4692):
³⁸Ar component breakdown (:4648–4664):
(quadrature, Utilities → toDP:4664)
³⁶Ar(Cl) from production ratio (:4744–4748):
Ca/K (rederived, lab calibration R = 0.52, :4710–4711):
v3.0.1 fix-record: prior version inverted Ca/K and used pr_ratio = 0.000377 instead of 0.52.
⁴⁰Ar(r+a):
Age in Ma (:4686–4687):
Utilities.py:307–490 (getDFStatistics_ls), 462–1535 (getDFStatistics_sh), 928–1050 (DFN group fit), 1340–1395 (DFI group fit), 1500–1535 (WMA & MSWD)
All formulas in this section use the York convention Y = a + b·X, where a = y-intercept and b = slope (per Vermeesch 2024 Eq. p.398, Vermeesch 2018 IsoplotR §11). pyADR's curve_fit(linear, x, y) returns popt = [b, a] (slope first, intercept second).
Linear fit gives Y = a_N + b_N·X where:
So slope = F, intercept = trapped 40/36 (direct read-off).
Linear fit gives Y = a_I + b_I·X where (Vermeesch 2024, eq. p.398):
Therefore (Utilities.py:1479):
Equivalently, the X-intercept of the inverse isochron equals 1/F (Vermeesch 2018 Geosci Frontiers, IsoplotR p.8).
Error propagation for F = −b/a including slope-intercept covariance from pcov_inv:
Implemented at Utilities.py:455-465 (LS) and Utilities.py:1479-1502 (SH).
v3.8.0 fix: Previously the code used
F = 1/inv_slope(i.e.F = 1/b). This is dimensionally and geometrically wrong — for typical samplesb < 0(X = 39/40 increases as Y = 36/40 decreases), so1/bis a small negative number, producinglog(1 + J·F) ≈ 0and absurdly young ages (e.g. SYL31 Sylhet Trap basalt, 115.4 ± 3.9 Ma per NTU thesis R94224113, returned T = 0.903 Ma in v3.7.4).
with the standard three-source error propagation (Utilities.py:466, 1503):
(σ_λ omitted; see §10.x note 3.)
where F_j = ⁴⁰Ar*_j / ³⁹Ar_K,j is computed per-step in calcAge (§7).
For N step ages {T_i ± σ_i}:
This is the maximum-likelihood estimator for a Normal model with two variance components (Vermeesch 2018, IsoplotR Eq. 5) and the standard form assumed by Schaen et al. (2021) GSA Bull. p.470 MSWD definition.
Implemented at Utilities.py:468-477 (LS) and Utilities.py:1514-1523 (SH).
v3.8.0 fix: The previous code wrote:
wma += (1/σ²·T) / (1/σ²)The numerator and denominator cancel inside the loop, so each iteration adds
T_iand the finalwma = Σ T_iis the sum, not the weighted mean. SYL31 with 40 spots reported WMA = 4481.9 (= 40 × ~112 Ma ≈ Σ T) instead of ~115 Ma. Output is internally inconsistent with the MSWD (which assumes WMA is the reference point).
For plateau / weighted-mean MSWD (Schaen et al. 2021 p.470):
For isochron-regression MSWD (group fits at Utilities.py:1029-1037, 1372-1380):
with ŷ_i = a + b·x_i from the linear fit. The N−2 reflects two fitted parameters (slope, intercept); the plateau form uses N−1 (one fitted parameter, the WMA itself).
v3.8.0 fix: Previously the plateau-MSWD reference was the arithmetic mean
T_sum/Ninstead of WMA. With a non-uniformσ_i, this gives a different — and inconsistent with the reported WMA — measure of dispersion. Schaen 2021 specifies WMA as the reference, matching the maximum-likelihood form.
Utilities.py:2033-2034. Computed before any plateau/isochron filtering.
⚠ Plateau "avg_age" still uses arithmetic mean (
Utilities.py:2031). Despite the §10.5 WMA being fixed in v3.8, the plateau display computesavg_age = np.mean([y_age[i] for i in range(n) if mask[i] == 1])i.e. unweighted arithmetic mean of unmasked step ages. The "weighted plateau age" shown in the plateau-step plot is therefore not consistent with the WMA reported in the regression summary (§10.5). For an external
±σtable, use the WMA from §10.5; treat the plot's avg-age annotation as visual only. Same pattern as the v3.7.x WMA bug — should be unified.
pyADR currently uses scipy.optimize.curve_fit(linear, x, y), which is ordinary least squares (OLS) with no weights. This:
- Ignores σ_y (no inverse-variance weighting of the fit)
- Ignores σ_x (assumes all uncertainty is in y)
- Ignores correlations between σ_x and σ_y (which are large for inverse isochrons because ⁴⁰Ar is the denominator of both axes)
The standard for ⁴⁰Ar/³⁹Ar isochron fitting is York regression (York et al. 2004; Vermeesch 2018, IsoplotR), which accounts for σ_x, σ_y, and ρ(x,y). Replacing OLS with York is a planned v3.9+ task (see CHANGELOG.md "仍待處理" list).
The MSWD formulas in §10.6 are unaffected by the OLS vs York choice as long as σ_y dominates and ρ(x,y) is small.
⚠ Intercept SE bug — running OLS on the error bars (
Utilities.py:389-390, 965-966). Two normal-isochron paths compute the intercept uncertainty by re-fitting:popt_std, _ = curve_fit(linear, x_std, y_std) # ← fit to the σ values themselves n_std = linear(0, *popt_std) # ← "intercept" of that fit, called n_stdThis treats the σ_x, σ_y arrays as if they were a new (x, y) dataset and reads its intercept. Mathematically meaningless — the resulting
n_stdis not the SE of the original-fit intercept; it is a slope-of-error-bars number with no statistical interpretation. The inverse-isochron paths (:441, :1286) correctly readpcov[1,1]from the main fit's covariance matrix — that's the form that should be used everywhere. Fix is one-line: replace eachcurve_fit(linear, x_std, y_std) → linear(0, *popt_std)pair withnp.sqrt(pcov[1, 1])from the main fit at the same call site.
For context: the v3.7.x F, WMA, and MSWD formulas in getDFStatistics_sh / getDFStatistics_ls were all buggy in distinct ways:
| Item | v3.7.x code | v3.8 fix | Reference |
|---|---|---|---|
| F (inverse isochron) | F = 1/inv_slope |
F = −b/a |
Vermeesch 2024 p.398 |
| F (LS isochron) | T = log(1 + J·iv) using Y-intercept |
F = −b/a |
Vermeesch 2024 p.398 |
| WMA | (1/σ²·T)/(1/σ²) inside loop → Σ T |
Σ(T/σ²)/Σ(1/σ²) |
Vermeesch 2018 Eq. 5 |
| MSWD reference | arithmetic mean | WMA | Schaen 2021 p.470 |
Validation: SYL31 LS (NTU thesis R94224113, Sylhet Trap basalt 115.4 ± 3.9 Ma).
PlaneFit3D.py:1–791
Reference implementation of Kent et al. (1990) maximum-likelihood plane regression for ⁴⁰Ar/³⁹Ar, following Wu (2007) NTU master thesis (R94224113, advisor: Ching-Hua Lo). Provides an alternative to the 2D isochron projection that retains all three isotope axes directly, avoiding error accumulation from ratio computation.
For a well-behaved sample with two end-members (trapped + radiogenic), the three Ar isotopes lie on a plane in (³⁶Ar, ³⁹Ar, ⁴⁰Ar) space:
Physical meaning:
-
$\alpha = (^{40}\text{Ar}/^{36}\text{Ar})_0$ = trapped (initial / atmospheric) composition (≈ 298.56 for pure air) -
$\beta = ^{40}\text{Ar}^*/^{39}\text{Ar}_K$ = radiogenic-to-K ratio (= F in the standard age equation)
Age from β: $$ T ;=; \frac{1}{\lambda} \ln(1 + \beta \cdot J) $$
Following Wu (2007) eq 3-7 and Kent et al. (1990), with each data point
Define the plane normal vector
After Lagrange-multiplier elimination of
Maximizing
PlaneFit3D.py:_Lp L126-129.
PlaneFit3D.py:build_cov L56-63.
where $k_0 = $ PR(⁴⁰Ar/³⁹Ar)_K production ratio (default 0.025004).
The off-diagonal
PlaneFit3D.py:_newton_raphson L174-198. PlaneFit3D.py:_grad, _hess L132-158.
Solve _ols_initial, L161-171):
where $$ g = -\frac{\partial L_p}{\partial \delta}\Big|{\delta_k}, \qquad H = -\frac{\partial^2 L_p}{\partial \delta \partial \delta^T}\Big|{\delta_k} $$
At the MLE,
Backtracking line search (v3.4 addition): halve the step size up to 20× until
Sign-convention note: pyADR's
_grad()returns$+\partial L_p/\partial \delta$ (opposite sign to Wu 2007 eq 3-10). Combined with$H = -\partial^2 L_p$ , the update$\delta_{k+1} = \delta_k + H^{-1} g$ moves uphill on$L_p$ (signs cancel). Documented inline in_grad's docstring.
PlaneFit3D.py:_compute_mswd L201-208.
(Wu 2007 eq 3-24; Mahon 1996.) Two free parameters (α, β) → df = n−2.
95% confidence interval for MSWD computed exactly via scipy.stats.chi2), not the normal approximation of Wendt & Carl (1991):
$$
\text{CI}{95%}(\text{MSWD}) ;=; \left[\frac{\chi^2{df}(0.025)}{df}, \frac{\chi^2_{df}(0.975)}{df}\right]
$$
If MSWD exceeds the upper bound, pyADR applies a Wendt-Carl style σ-expansion: $$ \tau^2 ;=; \frac{S^2}{df} ;=; \text{MSWD} \qquad \Longrightarrow \qquad \sigma_{\delta} \to \sqrt{\tau^2},\sigma_{\delta} $$
PlaneFit3D.py:_param_cov L211-226.
where
Wu (2007) eq 3-27 sign-error correction: The thesis writes
$\text{cov}(\hat{\delta}) = \tau^2 \cdot (\partial^2 L_p/\partial \delta \partial \delta^T)^{-1}$ . At the MLE,$\partial^2 L_p$ is negative-definite, so its inverse gives negative variances — clearly wrong. The correct ML asymptotic covariance is$\tau^2 \cdot (-\partial^2 L_p/\partial \delta \partial \delta^T)^{-1} = \tau^2 \cdot H^{-1}$ . pyADR uses the corrected form and documents this inline.
Marginal 1σ uncertainties: $$ \sigma_\alpha ;=; \sqrt{\text{cov}(\hat\delta){11}}, \qquad \sigma\beta ;=; \sqrt{\text{cov}(\hat\delta)_{22}} $$
PlaneFit3D.py:age_from_beta, age_error_1sigma L229-243.
Renne (1998) / Min et al. (2000) error propagation:
with partial derivatives: $$ \frac{\partial T}{\partial \beta} = \frac{J}{\lambda(1 + \beta J)}, \qquad \frac{\partial T}{\partial J} = \frac{\beta}{\lambda(1 + \beta J)}, \qquad \frac{\partial T}{\partial \lambda} = -\frac{\ln(1 + \beta J)}{\lambda^2} $$
PlaneFit3D.py:_apply_sigma_cap L66-71 (v3.4.3 addition).
For background-dominated steps where
with typical sigma_cap_rel parameter; None disables (classical Kent), per-axis values supported.
PlaneFit3D.py:find_subplanes L325+.
Sliding-window search over consecutive heating steps to identify the longest run with MSWD inside the 95% χ² CI. Used to detect sub-plateaus when not all steps share a single isotope system. Returns ranked candidates by window length.
Location: NTNU_DataReduction.py:3901-3902 (single fit) and :4489+ (group fit).
x40 = (df3["40Ar(r)"] + df3["40Ar(a)"]).values[m3]
s40 = np.hypot(df3["40Ar(r)_std"].values, df3["40Ar(a)_std"].values)[m3]Math problem:
$x_{40} = ^{40}\text{Ar}(r) + ^{40}\text{Ar}(a)$ - $^{40}\text{Ar}(r) = ^{40}\text{Ar}m - ^{40}\text{Ar}(a) - ^{40}\text{Ar}(K)$, so $^{40}\text{Ar}(r)$ and $^{40}\text{Ar}(a)$ are anti-correlated: $\text{cov}(^{40}\text{Ar}(r), ^{40}\text{Ar}(a)) = -\sigma{40a}^2$
- Correct variance:
$\sigma_{x_{40}}^2 = \sigma_{40r}^2 + \sigma_{40a}^2 + 2,\text{cov} = \sigma_{40r}^2 - \sigma_{40a}^2$ - pyADR's
np.hypotgives$\sigma_{40r}^2 + \sigma_{40a}^2$ — over by$2\sigma_{40a}^2$
Same pattern as v3.8.1 fix for σ_36(m)/σ_39(m). Effect: σ_40 inflated, MSWD systematically low, α/β/T potentially biased toward 36Ar/39Ar.
Proposed fix: s40 = np.sqrt(max(σ_40r² − σ_40a², 0)) (clip-to-zero in the rare case σ_40a > σ_40r). Validation sample: SYL31 (Wu 2007, 115.4 ± 3.9 Ma).
| Topic | Reference |
|---|---|
| 3D plane ML algorithm | Kent J.T. et al. (1990) Maximum likelihood estimation of a plane in three dimensions. Statistics 21: 411–426 |
| Implementation (NTU thesis) | Wu C.-Y. (2007) 3-D Plane-fitting Program in 40Ar/39Ar Dating. MSc thesis, NTU Geosciences (R94224113), advisor Ching-Hua Lo. Math derivations in Chapter 3 |
| Newton-Raphson + Lagrange | Titterington D.M., Halliday A.N. (1979) On the fitting of parallel isochrons and the method of maximum likelihood. Chem. Geol. 26: 183–195 |
| MSWD / residual analysis | Mahon K.I. (1996) The new "York" regression: application of an improved statistical method to geochemistry. Int. Geol. Rev. 38: 293–303 |
| σ-cap modified weighting | Mahon (1996); pyADR v3.4.3 addition |
| Age + σ_T propagation | Renne P.R. et al. (1998) Intercalibration of standards, absolute ages and uncertainties in 40Ar/39Ar dating. Chem. Geol. 145: 117–152. Min K. et al. (2000) A test for systematic errors in 40Ar/39Ar geochronology. GCA 64: 73–98 |
| Koppers age error tables | Koppers A.A.P. (2002) ArArCALC — software for 40Ar/39Ar age calculations. Comput. Geosci. 28: 605–619 |
| Validation sample | SYL31 (Sylhet Trap basalt, India), Rajmahal-Sylhet eruption 119–116 Ma, Kerguelen Plume early product. Wu 2007 chapter 5 |
-
calcAgetruncation — RESOLVED in v3.7.4. Function was truncated mid-statement in v3.7.0–v3.7.3; v3.7.4 restored the full implementation from the V3.4.1 archive. -
DiagramPlot F / WMA / MSWD bugs — RESOLVED in v3.8.0. Three distinct bugs in
getDFStatistics_sh/_ls: (a) inverse-isochronF = 1/slope(should be−b/a); (b) WMA loop with(1/σ²·T)/(1/σ²)cancels to Σ T; (c) MSWD reference was arithmetic mean instead of WMA. All three corrected. Validated against SYL31 LS (115.4 ± 3.9 Ma).
-
σ_T0 underestimated ~10× (§1) —
Utilities.calculateT0/REcalculateT0lines 187, 216, 275, 294. Usesstd(|residuals|)/√Ninstead of intercept SEsqrt(pcov[-1,-1]). AutoPipeline_fit_onewas patched in v3.8.2; legacy functions still pending. -
σ_J v₃ bracket bug (§6) —
Utilities.py:2864. Operator precedence parses((np.exp(l*t)) - 1 / F²)²as(e^λt − 1/F²)²instead of((e^λt − 1)/F²)². σ_J massively over-estimated for young samples. 5-min fix: add one pair of parentheses. -
Isochron intercept SE — OLS on the error bars (§10.8) —
Utilities.py:389-390(LS),965-966(SH normal).curve_fit(linear, x_std, y_std) → linear(0, *popt_std)is mathematically meaningless. Inverse-isochron paths (441,1286) correctly usepcov[1,1]; propagate that pattern.
-
Plateau "avg_age" uses arithmetic mean (§10.7) —
Utilities.py:2031.np.mean(y_age[mask==1])instead of the v3.8-corrected WMA from §10.5. Plot annotation inconsistent with the WMA in the regression summary. -
Linear-sum σ propagation in
getJVolumeStatistics/calcAge—Utilities.py:2839, 2926, 2942etc. UsesF·Σ|σ_x/x|instead of√Σ(∂F/∂x · σ_x)². Inconsistent with AutoPipeline_propagatequadrature. -
J/Salt outlier mask uses ±1·σ_SE —
getJStatistics:2209-2210,getSaltStatistics:2250. Too tight; masks ~32% of normal points. Use 2σ or Chauvenet. -
Isochron regression — OLS not York (§10.8) —
Utilities.py:360, 905, 1234. Ignores σ_x and ρ(x,y). York regression planned for v3.9+.
-
Ar_38_Airmis-named in calcAge (Utilities.py:2879). Variable contains ³⁸Ar(air) + ³⁸Ar(Cl), but is labelled "Ar_38_Air". Proper split only intoDP:4648–4664. -
T_stdin calcAge omits σ_λ (:2907). McDougall & Harrison 1999 eq. 4.7 has three partials; pyADR drops(∂T/∂λ)²σ_λ². ~0.5% relative under-estimate. -
λ source split:
constants[14]vsconstants[16]—getDegasPlot:444uses [14]; calcAge:2904,getStackPlot:1972use [16]. Audit PS to confirm [14]==[16] or unify. -
J calc uses hardcoded λ, Age calc reads
constants[16](getJVolumeStatistics:2748-2749). λ=5.531e-10 written in source. PS change creates systematic offset. -
Atmospheric ⁴⁰/³⁶ hardcoded 298.56 in salt functions (
:2708, 2737-2738). Rest of program readsconstants[12]. Code style: replace literal withconstants[12]. -
§2 T0 outlier mask threshold uses
|T_{0,j} − μ| > σ/2 + μ(getT0Statistics:2550). The+ μterm is dimensionally suspect. Verify with advisor. -
Salt σ on numerator subtraction (
calculateSlatCa:2709,calculateSlatK:2738) uses linear sum onAr36 − air/298.56. Should be quadrature. -
Datum CSV does not write cov[X,Y] (
NTNU_DR.py:3493). Without ρ(x,y), downstream York / IsoplotR refitting from the CSV is impossible. Add 4 columns per Vermeesch (2018) Eq. 2.
Primary ⁴⁰Ar/³⁹Ar data-reduction:
- McDougall, I., Harrison, T.M. (1999). Geochronology and Thermochronology by the ⁴⁰Ar/³⁹Ar Method, 2nd ed. Oxford UP.
- Koppers, A.A.P. (2002). ArArCALC — software for ⁴⁰Ar/³⁹Ar age calculations. Computers & Geosciences 28, 605–619.
- Schaen, A.J. et al. (2021). Interpreting and reporting ⁴⁰Ar/³⁹Ar geochronologic data. GSA Bulletin 133(3/4), 461–487.
Isochron regression & WMA statistics (§10):
- Vermeesch, P. (2024). Errorchrons and anchored isochrons in IsoplotR. Geochronology 6, 397–407.
- Vermeesch, P. (2018). IsoplotR: free open toolbox for geochronology. Geoscience Frontiers 9, 1479–1493.
- Vermeesch, P. (2015). Revised error propagation of ⁴⁰Ar/³⁹Ar data, including covariances. GCA 171, 325–337.
- Li, Y., Vermeesch, P. (2021). Inverse isochron regression for Re–Os, K–Ca, and other chronometers. Geochronology 3, 415–420.
- Powell, R., Green, E.C.R., Marillo Sialer, E., Woodhead, J. (2020). Robust Isochron Calculation. Geochronology.
- York, D., Evensen, N.M., Lopez-Martinez, M., De Basabe Delgado, J. (2004). Unified equations for the slope, intercept, and standard errors of the best straight line. Am. J. Phys. 72(3), 367–375.
- Kuiper, K.F. (2002). The interpretation of inverse isochron diagrams in ⁴⁰Ar/³⁹Ar geochronology. EPSL 203, 499–506.
Intercept SE for linear regression (§1):
- Li, X., Naeher, U., Pross, J. (2019). Mass spectrometric data processing in stable isotope analysis: regression-based standard errors of the y-intercept. J. Mass Spectrom. 54, 145–152.
Modified weighting & 3D plane fit (PlaneFit3D):
- Mahon, K.I. (1996). The New "York" Regression. Int. Geol. Rev. 38(4), 293–303.
- Kent, J.T., Watson, G.S., Onstott, T.C. (1990). Fitting straight lines and planes with an application to radiometric dating. EPSL 97, 1–17.
- Wu, M.-W. (2007). 3-D Plane-fitting Program in 40Ar/39Ar Dating. NTU MSc thesis R94224113.