Thinking in Waves Episode 59 / Part XIX — Back to Physics, with the Tools
From estimating systematic error to measuring it
Measuring the systematic error ──
and what helped most was not the window
Turning Episode 58’s “change the window and the value moves” around, we built a procedure that hunts for the region where it does not move.
★★★ What turned up on the way ── that overshoot appears with zero noise too. Its source is frequencies below the record’s resolution ── the part that, inside the record, looks like nothing but a trend.
★★ So drawing one straight line removes it (\(\beta=5\): \(5.512\to4.992\)) ── Episode 58 said “raise the window order,” but there was a much cheaper move.
Episode 58 produced a systematic error with a fixed sign. ★ Here we turn it into an instrument ── extract both the true value and its uncertainty from how the estimate moves.
01Sweeping \(\hat\beta\) over \(\alpha\)
Window order \(\alpha=1\) (rectangular) through 8, averaged over 32 realizations.
| true \(\beta\) | \(a{=}1\) | 2 | 3 | 4 | 5 | 6 | 7 | 8 |
|---|---|---|---|---|---|---|---|---|
| 1 | 0.990 | 0.989 | 0.991 | 0.988 | 0.988 | 0.987 | 0.988 | 0.990 |
| 3 | 2.033 | 3.160 | 2.992 | 2.990 | 2.989 | 2.988 | 2.989 | 2.990 |
| 5 | 1.984 | 4.020 | 5.512 | 5.061 | 4.990 | 4.990 | 4.990 | 4.991 |
| 6 | 1.981 | 4.021 | 6.000 | 7.623 | 5.997 | 5.992 | 5.990 | 5.992 |
| 7 | 1.981 | 4.021 | 6.004 | 8.001 | 7.633 | 7.182 | 6.992 | 6.994 |
The core of §01
★ At small \(\alpha\) the estimate is clamped by leakage (Episode 58’s ceiling \(2\alpha\)).
★★ Just below the ceiling it overshoots ── 7.623 at \(\beta=6,\alpha=4\) and 8.001 at \(\beta=7,\alpha=4\), both pushing up towards the ceiling \(2\alpha=8\).
★ At larger \(\alpha\) it settles, lining up at a constant about 0.01 below the truth for every \(\beta\) ── that is the plateau.
★★★ That is the backbone of the procedure: find the plateau automatically and report its width as the systematic error.
02★★★ Killing the noise to see the bias alone
Convolving the expected periodogram \(E[P(k)]=\int S(f)|W(k-f)|^2df\) directly gives the bias with zero statistical error.
| true \(\beta\) | \(a{=}1\) | 2 | 3 | 4 | 5 | 6 | 7 | 8 |
|---|---|---|---|---|---|---|---|---|
| 1 | 0.005 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 |
| 3 | −1.009 | +0.125 | 0.001 | 0.001 | 0.001 | 0.001 | 0.001 | 0.002 |
| 5 | −3.012 | −0.988 | +0.625 | +0.068 | 0.003 | 0.003 | 0.004 | 0.004 |
| 7 | −5.012 | −2.988 | −1.000 | +1.002 | +0.918 | +0.252 | 0.007 | 0.008 |
★★★ The overshoot appears with zero noise (\(+0.625\) at \(\beta=5,\alpha=3\)) ── matching the \(+0.61\) measured in Episode 58. It is not a statistical effect but the shape of the expected spectrum itself.
★★ A prediction failed: I expected a main-lobe smoothing bias to grow at large \(\alpha\), but it is only \(+0.002\) to \(+0.008\). What limits the large-\(\alpha\) side is not bias but variance (§04).
★ So where does the leakage come from? Move the lower limit of the integral (\(\beta=5,\alpha=3\)):
| lower limit \(f\) [bins] | 0.0625 | 0.25 | 0.5 | 1.0 |
|---|---|---|---|---|
| estimated \(\beta\) | 5.625 | 5.062 | 5.006 | 5.002 |
The core of §02 ── the culprit identified
★★★ The overshoot appears only when frequencies below the record’s 1-bin resolution (0.0625–1 bin) are included.
★★ The leakage is dominated by “the part that looks like nothing but a trend inside the record” ── waves whose period exceeds the record. Exactly the component Episodes 5 and 56 called “indistinguishable from a constant.”
★ Not merely indistinguishable ── it was contaminating the measurement.
03★★★ Then subtract the trend
If it looks like a trend inside the record, fitting and subtracting one should remove it. With a Hann window (\(\alpha=3\)), subtracting a degree-\(q\) polynomial first:
| true \(\beta\) | none | \(q{=}0\) (mean) | \(q{=}1\) (line) | \(q{=}2\) | \(q{=}3\) |
|---|---|---|---|---|---|
| 3 | 2.992 | 2.992 | 2.991 | 2.991 | 2.991 |
| 4 | 4.000 | 4.000 | 3.992 | 3.992 | 3.991 |
| 5 | 5.512 | 5.512 | 4.992 | 4.993 | 4.992 |
| 6 | 6.000 | 6.000 | 6.035 | 5.994 | 5.992 |
| 7 | 6.004 | 6.004 | 6.504 | 6.039 | 6.574 |
The core of §03 ── the most practical thing here
One straight line removes the overshoot.
\(\beta=5\): \(5.512\to4.992\)
★★★ Episode 58 said “raise the window order”; there was a cheaper move. Raising the order costs variance (§04); subtracting a trend is nearly free.
★ Subtracting the mean alone (\(q=0\)) does nothing ── it takes degree 1.
★★ But the ceiling does not move ── \(\beta=7\) measured with \(\alpha=3\) (ceiling 6) stays at 6.04–6.57 for every \(q\). Detrending removes the overshoot but not the clamping.
★ Over-subtracting hurts too (\(\beta=7\) with \(q=3\) worsens to 6.57). The handbook entry: subtract a degree 1–2 trend, then choose a window with \(2\alpha\ge\beta+2\).
04The scatter grows with \(\alpha\)
A window reduces the number of independent points. Equivalent noise bandwidth \(\mathrm{ENBW}=N\sum w^2/(\sum w)^2\).
| \(\alpha\) | ENBW [bins] | measured sd(\(\beta\)) | Episode 51’s prediction | ratio |
|---|---|---|---|---|
| 1 | 1.000 | 0.1170 | 0.0865 | (biased) |
| 2 | 1.234 | 0.1815 | 0.0961 | (biased) |
| 3 | 1.500 | 0.0839 | 0.1059 | 0.792 |
| 5 | 1.944 | 0.1011 | 0.1206 | 0.838 |
| 8 | 2.473 | 0.1193 | 0.1360 | 0.877 |
★ The prediction is Episode 51’s \(\mathrm{sd}=1.2825/\sqrt{\text{effective points}}\) with effective points \(=S_{xx}/\mathrm{ENBW}\). At \(\alpha=1,2\) the estimator is still biased at \(\beta=3\), so those sd values are meaningless.
★★ In the unbiased range \(\alpha\ge3\) the ratio is 0.79–0.88 ── the prediction runs 15–20 % high. The periodogram’s points are correlated by the window, and ENBW alone overstates that correlation. Use it as an upper bound.
05The best \(\alpha\) ── and a minimum you must not trust
| true \(\beta\) | \(a{=}1\) | 2 | 3 | 4 | 5 | 6 | 7 | 8 |
|---|---|---|---|---|---|---|---|---|
| 3 | 1.011 | 0.242 | 0.085 | 0.093 | 0.100 | 0.106 | 0.112 | 0.116 |
| 4 | 2.019 | 0.026 | 0.085 | 0.093 | 0.099 | 0.105 | 0.111 | 0.116 |
| 6 | 4.019 | 1.978 | 0.008 | 1.686 | 0.097 | 0.104 | 0.110 | 0.115 |
| 7 | 5.019 | 2.979 | 0.998 | 1.001 | 0.885 | 0.254 | 0.110 | 0.115 |
The core of §05 (the table is RMS error)
★★★ The 0.026 at \(\beta=4,\alpha=2\) and the 0.008 at \(\beta=6,\alpha=3\) are two biases cancelling by accident. Without knowing \(\beta\) you cannot aim at them ── do not ride that minimum.
★ The neighbours give it away: at \(\beta=6\), \(\alpha=3\) gives 0.008 but \(\alpha=4\) gives 1.686. A minimum that unstable is not a usable minimum.
★★ The safe rule is \(2\alpha\ge\beta+2\) (\(\alpha\ge3\) for \(\beta=4\), \(\alpha\ge4\) for \(\beta=6\)). You cannot apply it without knowing \(\beta\), which is why the §01 sweep is needed.
06★★ Making it a procedure ── and measuring how often it is right
① measure \(\hat\beta\) for \(\alpha=1\ldots8\) ② find the longest run whose spread is within \(\mathrm{tol}=0.15\) ③ the value is its mean, the systematic error its half-spread ④ the statistical error from Episode 51’s formula (with ENBW) ⑤ add in quadrature and report
| true \(\beta\) | mean estimate | total error | \(|\)err\(|\) median | 68th pct | 95th pct | actual sd | coverage |
|---|---|---|---|---|---|---|---|
| 1 | 1.015 | 0.129 | 0.057 | 0.078 | 0.178 | 0.317 | 87.0 % |
| 3 | 3.021 | 0.129 | 0.060 | 0.092 | 0.187 | 0.366 | 84.0 % |
| 5 | 5.014 | 0.132 | 0.062 | 0.093 | 0.201 | 0.224 | 84.5 % |
| 7 | 7.009 | 0.136 | 0.068 | 0.105 | 0.216 | 0.121 | 79.5 % |
The core of §06 ── the error distribution is not Gaussian
★ Against a nominal 68 %, the measured coverage is 79–87 % ── the quoted error runs about 1.4× high (the 68th percentile of \(|\)error\(|\) is 0.08–0.11 while the quoted total is 0.13).
★★★ Yet the actual sd is far larger than the quoted error (0.32 against 0.13 at \(\beta=1\)). No contradiction ── the error distribution is not Gaussian.
★★ The core is narrow and it misses badly on rare occasions (the 95th percentile is about twice the 68th, and the sd runs 0.12–0.41). Those rare misses are what inflate the sd.
★★★ So this procedure’s error must not be written as an sd. “68 % fall within this width” ── stated as coverage ── is the right form, for exactly the reason Episode 51’s \(\ln(\chi^2_2/2)\) was asymmetric.
07Where the procedure breaks
| fitting band | true \(\beta\) | estimate | total error | coverage |
|---|---|---|---|---|
| \(k=32\ldots512\) | 3 | 3.021 | 0.129 | 84.0 % |
| \(k=32\ldots512\) | 6 | 6.006 | 0.133 | 83.0 % |
| \(k=64\ldots256\) | 3 | 3.027 | 0.352 | 75.5 % |
| \(k=64\ldots256\) | 6 | 6.015 | 0.367 | 78.5 % |
| \(k=16\ldots1024\) | 3 | 3.019 | 0.074 | 86.5 % |
| \(k=16\ldots1024\) | 6 | 6.008 | 0.075 | 85.0 % |
★ Change the band and the coverage moves (75.5–86.5 %) ── even though the estimate itself barely moves (around 3.02).
★★ So this procedure is “conditional” too (Episode 55). Three assumptions: ① the fitting band ② tol ③ the window family. Building a procedure does not make the assumptions go away.
Exercises
- §02 identified the culprit as periods longer than the record. How does that connect to Episode 56?
Show answer
Episode 56 said they are invisible; here we learn they are invisible and still doing damage.
★ Episode 56: “the knee from \(\tau\) longer than the observation cannot be seen.” ★★ Here: the tail still leaks through the window into high frequency and shifts the slope by \(+0.6\).
★★★ This is what Episode 57’s exercise meant by “a smaller \(n\) lets outside information reach further in.” There I wrote that one might in principle estimate \(\tau_{\max}\) from it; what actually happens first is contamination.
★ The order matters: subtract first (§03), then measure. Inference about the outside comes after that. - §03 found that subtracting the mean does nothing. Why does it take degree 1?
Show answer
Because the window has already killed the zeroth-order leakage. A Hann window (\(\alpha=3\)) has both value and derivative vanishing at the ends, so leakage from a constant (\(f=0\)) is already down by \(f^{-3}\).
★ What remains is low but nonzero frequency ── the part that looks like a straight line inside the record. ★★ Subtracting a line removes most of it.
★★★ Generally: subtracting degree \(q\) removes “low frequencies that look like a degree-\(q\) polynomial inside the record,” which acts like raising the window order. ★ The difference is the price: raising the order costs variance (§04), while detrending costs only \(q+1\) degrees of freedom ── nearly free at 4096 points. - §06 says not to write this as an sd. Are the series’ earlier numbers safe?
Show answer
Episode 51’s \(\mathrm{sd}(\beta)=1.2825/\sqrt{M}\) is fine ── that is the standard deviation of \(\ln(\chi^2_2/2)\) itself, not the output of a procedure.
★ The risky ones are the five Episode 55 classified as “conditional on a model” ── those are closer to procedure outputs, and their distributions need not be Gaussian.
★★ To add to Episode 55’s audit: for a “conditional” number, also state the shape of its error (Gaussian, or narrow-core with heavy tails).
★★★ And measuring the coverage, as here, is the surest way ── measure the error’s shape before its size.
Summary — systematic error is measurable, but the window was not the point
Turning Episode 58’s “change the window and the value moves” around, we built a procedure that hunts for where it does not move.
★ \(\hat\beta(\alpha)\) has a plateau ── clamped at the ceiling \(2\alpha\) for small \(\alpha\), overshooting just below it, settling for large \(\alpha\).
★★★ And what turned up on the way ──
the overshoot appears with zero noise.
Its source is the part that looks like nothing but a trend inside the record.
★★ Moving the integral’s lower limit settled it ── the \(+0.625\) overshoot appears only when frequencies below the 1-bin resolution (0.0625–1 bin) are included. Episodes 5 and 56’s “a period longer than the observation is indistinguishable from a constant” becomes “not merely indistinguishable ── it was contaminating the measurement.”
★★★ So one straight line removes it (\(\beta=5\): \(5.512\to4.992\); \(\beta=6\) needs \(q=2\), giving 5.994). Episode 58 said “raise the window order,” but there was a cheaper move ── raising the order costs variance, while detrending is nearly free at 4096 points. ★ The ceiling does not move, though (\(\beta=7\) with \(\alpha=3\) stays at 6.04–6.57 for every \(q\)).
★★ One prediction failed: I expected a smoothing bias to grow at large \(\alpha\), but it is \(+0.002\) to \(+0.008\) and negligible; what limits the large-\(\alpha\) side is variance. ★ And the \(\alpha\) with minimum RMS is sometimes set by two biases cancelling by accident (0.008 at \(\beta=6,\alpha=3\), against 1.686 at \(\alpha=4\)) ── a minimum not to ride. The safe rule is \(2\alpha\ge\beta+2\).
★★★ We also measured how often the procedure is right ── coverage 79–87 % against a nominal 68 %, while the actual sd runs 1–3× the quoted error. Not a contradiction: the error distribution is not Gaussian ── a narrow core with rare large misses. ★ So this procedure’s error belongs in coverage, not in an sd, for exactly the reason Episode 51’s \(\ln(\chi^2_2/2)\) was asymmetric.
★ And the procedure itself is conditional ── on the fitting band, on tol, and on the window family; changing the band moves coverage between 75.5 % and 86.5 %. Building a procedure does not make the assumptions go away.
★ Next we look at that part itself ── where exactly does a trend end and a signal begin?
★★ Episode 5’s “indistinguishable from a constant,” Episode 29’s renormalisation, and this episode’s detrending ── all of them cut off the long end. ★ We ask where what was cut off goes.
This document is Episode 59 of the “Thinking in Waves” series. The computations are in kenshou/keito.js (own FFT, zero dependencies). On the numbers in this episode: the \(\hat\beta(\alpha)\) sweep, the bias from the noiseless expected periodogram, the dependence on the integral’s lower limit, the effect of detrending, ENBW and scatter, the RMS errors, and the procedure’s coverage and error quantiles are all computed here. Spectral leakage in the periodogram, equivalent noise bandwidth, and the practice of detrending before taking a periodogram are all well-established and are not discoveries of this episode ── what this episode did was organise them on Episode 58’s \(\alpha\) yardstick, locate the overshoot’s source through the lower-limit dependence, and measure the procedure’s coverage over 200 realizations. Every number depends on a record length of 4096, a synthesis length of 65536, a fitting band \(k=32\)–512 and \(\mathrm{tol}=0.15\); as §07 shows, changing the band moves coverage between 75.5 % and 86.5 %. The “0.79–0.88 of the prediction” in §04 is for \(\beta=3\) over 200 realizations. The “minimum not to ride” in §05 is specific to these synthesis conditions; which \(\beta\) and \(\alpha\) happen to cancel varies with them. This episode also uses synthetic data only, and has not been applied to real measurements. ── To print, use your browser’s Print → Save as PDF.