The error-aware likelihood

Folding the surrogate's own uncertainty into a parameter-estimation likelihood, so that parameters are measured less confidently wherever the model is less sure. Implemented in nrhjsurrogate/driver/fisher.py.

The short version: the effective noise becomes S_eff = S_n + S_w(theta), the detector noise plus a model-error spectrum built from the surrogate's stored variances. That addition is not a heuristic inflation. It is the exact result of marginalising the likelihood over the surrogate's own posterior, and the derivation below says precisely which premises make it exact and which of them this implementation actually satisfies.

Why variances add

Adding two independent zero-mean Gaussians adds their covariances. A power spectral density is a variance density, so power spectral densities add linearly. Amplitude spectral densities do NOT add linearly; they add in quadrature, which is the same statement written differently. S_n + S_w is the linear addition of variances and is the correct form.

The derivation

Setup. The data is one true signal plus noise,

d = h_true + n ,        n ~ Normal(0, C_n) .

The surrogate predicts h_model(theta). It is a Gaussian process regression, so what it actually holds is not a single set of coefficients but a POSTERIOR over them:

C | training ~ Normal(C_hat(theta), Sigma(theta)) .

The prediction h_model is the posterior MEAN, built from C_hat. Writing a draw from the posterior as C = C_hat + epsilon and letting J be the map from coefficients to waveform samples,

h(theta) = h_model(theta) + J epsilon ,     epsilon ~ Normal(0, Sigma) .

Marginalise. Conditional on a particular draw the likelihood is ordinary:

p(d | theta, epsilon) = Normal( d ; h_model(theta) + J epsilon , C_n ) .

We do not know epsilon, we know its distribution, so integrate it out. The integral is a Gaussian convolution and is analytic:

p(d | theta) = integral Normal(d ; h_model + J epsilon, C_n)
                        Normal(epsilon ; 0, Sigma) d epsilon
             = Normal( d ; h_model(theta) , C_n + J Sigma J^T ) .

So with C_w = J Sigma J^T,

C_eff(theta) = C_n + C_w(theta) ,

ln L(theta) = -1/2 (d - h_model)^T C_eff^-1 (d - h_model)
              -1/2 ln det C_eff  +  const .

This is exact given the premises. There is no approximation in the step from "the surrogate has a Gaussian posterior" to "add the covariances". That is worth being clear about, because an inflated noise budget usually IS a heuristic, and here it is not.

The frequency-domain form, as implemented

InnerProduct (fisher.py:634) is band-limited and frequency-weighted, so C_eff is represented by a spectrum. Under the stationarity assumption discussed in the limits, the two terms become

S_eff(f) = S_n(f) + S_w(f; theta) ,        weight = 4 df / S_eff .

Because C_eff depends on theta, the normalisation is NOT a constant and cannot be dropped. It is carried (fisher.py:891):

logdet = -0.5 * np.sum(np.log(s_eff / self.ip.snf[None, :]), axis=1)
val    = -0.5 * self.ip(r, r, weight=weight) + logdet

What that term does, and why there is an alternative. It penalises regions of parameter space where the model is less sure, so the posterior is pushed away from them. Under the marginalisation above that is formally correct. It is also easy to misread: it is not the same as being uncertain about those regions, and it can present as a preference for parameters where the surrogate happens to be confident. frozen_error_psd (fisher.py:780) evaluates S_w once at a reference point and holds it fixed, which removes the penalty at the cost of no longer varying the error budget with theta. Its own docstring names this penalty as the reason it exists. Which of the two is appropriate is a choice about the claim being made, not a bug in either.

Limits

Ranked by how much they should worry a reader.

1. RETRACTED 2026-09-02, and the correction is the interesting part.

What this section said, and what was wrong with it, is kept here rather than edited away, because the mistake is easy to make again.

IT SAID: the GPR knows only interpolation uncertainty, so the merger-ringdown arm's recorded 100 to 200 percent error is a systematic with no representation in S_w, and marginalising cannot help.

THAT IS BACKWARDS. The merger-ringdown arm's error was established to be parameter-space REGRESSION, essentially in full. Comparing the GPR's predicted element coefficients against each simulation's OWN directly fitted coefficients, with the same pinned breaks, the same p_max, the same remnant map and QNM detrend, so that only the coefficient source changes: pooled e_gpr 3.44e-03 against e_own 2.55e-05, a ratio of 115x, cross term at or below 3e-4 so the split is exact. Waveform level, own fit against oracle, 1.000 median. All ten elements regression dominated at 0.9996 to 1.0000. Not the representation, not the data.

A Gaussian process posterior variance IS an estimate of its own regression uncertainty. So the merger-ringdown error is precisely the KIND of error this construction describes, not a systematic it cannot see.

The 100 to 200 percent number was also being misused. The same log explains it: the absolute regression error is nearly FLAT across the segment, 5.9e-04 rising to 1.4e-03, while local truth power falls 23x, so only the FRACTIONAL error rises. It describes late ringdown where the signal has collapsed, not the merger in general.

SO THE REAL LIMIT IS CALIBRATION, NOT BLINDNESS. The question is whether S_w is correctly SIZED, and one measurement bears on it: the in-sample full-IMR mismatch of 7.90e-07 on (2,2) implies about 1.26e-03 relative error while the envelope reports 4.77e-04, so the envelope is a factor 2.6 OPTIMISTIC. That is in sample. Out of sample it will be worse by an amount nobody has measured.

Cutler and Vallisneri still applies to whatever error the GPR genuinely cannot see, and a marginalisation broadens rather than de-biases. The point of this retraction is only that the merger-ringdown arm is not the example of that.

1b. The original concern, in the form that survives. The posterior variance describes interpolation uncertainty between training simulations. It does not describe a systematic deficiency of the fit itself. The merger-ringdown arm is recorded in this repository as a 100 to 200 percent problem at parameter-estimation masses, regression dominated rather than representation dominated, and the GPR does not know that. That component is a BIAS with no representation in S_w, and marginalising cannot help with a bias.

This is the standard result: a waveform-model error primarily SHIFTS the posterior, by of order one sigma, and only secondarily broadens it (Cutler and Vallisneri 2007, cited at scripts/pe_error_broadening.py:39). The construction on this page captures the secondary effect and is silent on the primary one. It makes parameter estimation honest about interpolation uncertainty. It does not make a biased model unbiased.

2. Dynamics uncertainty is reported but still not propagated, and it is the LARGER term. The surrogate stores variances for 'E', 'J', 'k', 'tau' and 'xpc' alongside the mode content. They used to be dropped: mode_errors kept only tuple-keyed entries and so discarded all five.

REACHABLE since 2026-09-02 through NRHJSurAA.get_dynamics_error_estimate, a separate method so that a caller who does not want them pays nothing (MEASURED: 4.9 ms, against 109 ms for a waveform and 76 ms for the dynamics values themselves, one point, numpy backend, OMP_NUM_THREADS=1). On NRHJSur3dq8_AA_v3 at q = 4.0, chi1z = 0.3, chi2z = -0.2 the medians are tau 2.3796e-01 M, x 1.6956e-05, J 8.5713e-05, E 1.5347e-06, k 7.1095e-03.

STILL EXCLUDED FROM THE BUDGET, and this is the part that matters. tau is the timing content, a timing error is a phase error, and phase is the quantity a gravitational-wave likelihood is most sensitive to. A clock error at fixed orbital phase is an orbital-phase error at fixed time, delta Phi = -omega delta tau to first order, and mode (l, m) carries m times it. MEASURED at the same point, against the per-mode content uncertainty the budget does report:

The excluded channel is the larger one on every mode, and on the dominant mode it is 28 times larger. So S_w as built today is a small fraction of the model's own phase uncertainty.

Three things stand between that table and folding the term in, and none of them is effort. First it is a LINEARISATION: the time-domain waveform needs Phi(t), which is the INVERSE of tau(Phi), and delta Phi = -omega delta tau is only the first-order term of that inversion; the marginalisation derived above is exact only for a LINEAR map from coefficients to samples. With sigma_tau up to 0.43 M the second-order term has not been measured. Second, the clock and the mode content are fitted from the SAME simulations by the same GPR, so they are not independent and the derivation adds variances only for independent contributions; this package's precedent for correlated fits is the conservative linear sum, not quadrature. Third, folding it would move four published per-mode numbers by more than an order of magnitude. It is a measurement task before it is a coding task.

3. Only the diagonal of Sigma is propagated. The derivation needs the full posterior covariance J Sigma J^T. error_envelope sums T_j(x)^2 Var[C_j] over coefficients, which is J diag(Sigma) J^T, discarding coefficient correlations. Whether that overstates or understates the result depends on the sign of the discarded terms and has not been measured.

4. C_w is not stationary, and the implementation makes it so. The model error is time localised, small early and larger near merger. Representing it as a spectrum forces stationarity. The waveform_error_psd docstring states this plainly: it ignores the true time correlation and non-stationarity of C_w, keeping the leading effect that there is more effective noise where the surrogate is less sure.

5. Two spectral conventions that are measured, not cosmetic. Both exist because a spectrum is being built from a slowly varying amplitude.

• MODULATION. sigma_lm(t) is a non-negative amplitude, while the real error oscillates at the mode frequency. Taking the periodogram of sigma directly piles power near DC, largely below f_min where the band mask discards it, and the size of that deficit MOVES with total mass. Each mode's envelope is therefore ridden on that mode's own phase, sigma_lm exp(i arg h_lm), before projection, which is also the convention the envelope's calibration was measured in.
• INCOHERENT MODE SUM. Per-mode error series are projected separately and their periodograms added. Measured across inclination and total mass, the incoherent rule tracks the true combined error to within about 1.5x everywhere, while the coherent rule over-inflates by 2 to 3x at moderate to edge-on inclinations (scripts/swsph_likelihood_envelope.py).

6. The merger-ringdown envelope is a different quantity and must not be taken from error_envelope. Under the SWSpH representation the regressed coefficients are the detrended spheroidal residual, so error_envelope returns the uncertainty of that residual, before re-winding, mixing and the remnant-map term. That is wrong rather than merely uncalibrated. spherical_plunge_spheroidal_ringdown_error_envelope is the correct source. The released arm IS SWSpH.

Since 2026-09-02 the dispatch is made in one place, waveform_error_estimate.merger_ringdown_error_envelope, which both waveform_error_psd and NRHJSurAA.get_waveform_error_estimate call, and it GUARDS rather than guesses: an arm that carries the spheroidal ringdown and plunge groups but is not the recognised class raises instead of falling through to the plain coefficient envelope.

Its remnant_sigma defaults, (1.2e-4, 5.3e-5), are SOURCED: they are loo_Mf_rms and loo_chif_rms in paper_data/swsph_remnant_map.npz, the exact-leave-one-out errors of the shipped remnant map written by scripts/swsph_remnant_map.py. The final-mass number excludes SXS:BBH:1124 as a catalogue metadata outlier; keeping it gives 2.7e-4 instead. Both are RMS, and the maxima are about 5x and 8x larger.

7. inflate is a stress-test knob only. inflate = 1 is the calibrated envelope. Any other value no longer represents the model's actual uncertainty and must not be used to produce a quoted result.

In the time domain

Nothing in the derivation is specific to a Fourier basis. Marginalising a Gaussian coefficient posterior gives C_n + C_w in any basis, and the normalisation becomes ln det C_eff.

What changes is that limits 4 and 5 stop being necessary. The modulation convention exists only because a periodogram of a slowly varying amplitude misplaces power in frequency; with the error series in hand there is no periodogram. The incoherent-sum rule was calibrated as a rule for combining periodograms and would have to be revisited rather than carried across.

The obstacle is that a time-domain whitening built on a Toeplitz solver assumes stationary noise. tdanalysis builds its Gohberg-Semencul form from an autocorrelation (tdanalysis/likelihoods/single/pe.py, via compute_invLij_from_acorr), so a non-stationary C_w breaks the fast path.

NOT YET TESTED, and recorded here as a proposal rather than a result: C_w is naturally low rank. It is J Sigma J^T with Sigma over a modest number of Chebyshev coefficients, 14 elements by 25 coefficients per mode on the released artifact, which is small against the number of samples in a segment. That is the low-rank-update case, so the Woodbury identity gives the inverse and the matrix determinant lemma gives the determinant, both cheaply, while Gohberg-Semencul continues to handle the stationary detector noise. If that holds, the exact non-stationary correlated model-error covariance becomes affordable and limits 4 and 5 go away together. The cost has not been measured, and it should be measured on a batch, since parameter estimation is batched.

Where the code is

What the per-mode envelope covers

Restated here because limit 1 and limit 2 are both about what S_w does NOT contain, and a reader needs the positive statement too. Since 2026-09-02 the envelope that feeds S_w covers ALL SEVEN served modes over the WHOLE waveform span, not four modes over the inspiral. The four-mode reach was an artifact of the AA arm spreading one mode over three content groups while the documented route pointed at the one group holding only the odd-m modes; the AdA arm, which is what waveform_error_psd actually reads, keeps all seven in one group and was never limited. The two arms' envelopes over the join disagree by factors 0.43 to 6.94 depending on mode and parameter point, worst and most systematically on (2,2); the ramp spreads that over 15 M rather than concentrating it, and the joined envelope is continuous (verified by grid refinement, the estimator converging as dt^3).