Spec the Monte Carlo estimator of the true entropy rate
2 changed files+390−1
README.mdmodified+5−1View file
@@ -32,7 +32,11 @@ Where the spectrum sits well above one step² this reduces to the classical
3232 Gaussian entropy rate ½log₂(2πe) + ∫log₂S df; the noise floor keeps it finite
3333 and positive where a deep stopband would send that integral to −∞. LPC + ANS
3434 should approach R; probing where the approximation holds is the point. The
35-math section is a stub for the full derivation.
35+math section is a stub for the full derivation. The formula is *not* exact —
36+[docs/mc-true-rate.md](docs/mc-true-rate.md) analyzes its error sources and
37+specifies a sequential Monte Carlo estimator of the true entropy rate
38+(implementable standalone, e.g. in Python) with a two-sided convergence
39+certificate.
3640
3741 ## Run it
3842
docs/mc-true-rate.mdadded+385−0View file
@@ -0,0 +1,385 @@
1+# Monte Carlo estimation of the true entropy rate of a quantized filtered Gaussian process
2+
3+**Purpose.** This document specifies, in enough detail for independent
4+implementation and testing (Python + numpy/scipy is the intended target), a
5+Monte Carlo method for estimating the *exact* entropy rate R of the integer
6+process produced by the pipeline below — the quantity that the closed-form
7+formula in the companion app ([timeseries-compressibility](../README.md))
8+only approximates. The estimator is consistent (converges to the true R as
9+particle count, window length, and sample count grow), all of its systematic
10+errors push in a known direction (upward), and it admits a two-sided
11+sandwich so convergence can be certified rather than assumed.
12+
13+No code from the app needs to be reused; everything required is specified
14+here. Reference values measured with the app are given in §7 for
15+cross-checking.
16+
17+---
18+
19+## 1. The system under study
20+
21+### 1.1 Generating model
22+
23+$$x_n \sim \mathcal{N}(0,\sigma^2)\ \text{i.i.d.}, \qquad
24+y = h * x, \qquad
25+z_n = \operatorname{round}(y_n + d_n)$$
26+
27+- $h_0,\dots,h_{L-1}$ is a known FIR kernel (designs in §1.2).
28+- $d_n \sim \mathcal{U}[-\tfrac12,\tfrac12)$ i.i.d. when **dither** is on;
29+ $d_n = 0$ otherwise. Dither is *non-subtractive*: the stored integers
30+ include it and the decoder never learns $d$.
31+- Quantization step = 1 (σ is measured in steps). Rounding ties are a
32+ measure-zero event, so the tie convention does not affect any law.
33+- $y$ is a stationary Gaussian process with autocovariance
34+
35+ $$r_y(\ell) = \sigma^2 \sum_{j} h_j\, h_{j+\ell}
36+ \qquad (r_y(\ell)=0 \text{ for } |\ell| \ge L)$$
37+
38+ and spectral density $S(f) = \sigma^2 |H(f)|^2$,
39+ $H(f) = \sum_n h_n e^{-2\pi i f n}$, $f \in [0,\tfrac12]$ cycles/sample.
40+
41+The object of interest is the **entropy rate of the integer process** $z$:
42+
43+$$R \;=\; \lim_{k\to\infty} H\!\left(z_0 \mid z_{-1},\dots,z_{-k}\right)
44+\qquad \text{[bits/sample]}$$
45+
46+This is the true lossless compression limit for $z$, the number the app's
47+LPC + ANS bars try to approach from above.
48+
49+### 1.2 Kernel designs (must match to compare numbers)
50+
51+All comparisons in §7 use the **bandpass** kernel below; implement at least
52+`none` and `bandpass`.
53+
54+- `none`: $h = [1]$ (then $z$ is i.i.d. and everything has closed form —
55+ the primary validation case).
56+- Hamming-windowed sinc lowpass, cutoff $f_c$ cycles/sample, odd length $n$,
57+ midpoint $m = (n-1)/2$:
58+
59+ $$h_i = w_i \cdot \begin{cases} 2f_c & i = m \\
60+ \sin(2\pi f_c (i-m)) / (\pi (i-m)) & i \ne m \end{cases},
61+ \qquad w_i = 0.54 - 0.46\cos\!\left(\tfrac{2\pi i}{n-1}\right)$$
62+
63+ then normalized to unit DC gain: $h \leftarrow h / \sum_i h_i$.
64+- `bandpass(f_1, f_2)` = lowpass($f_2$) − lowpass($f_1$), same length, no
65+ further normalization.
66+- The **default ephys case**: bandpass 300–6000 Hz at 30 kHz sample rate
67+ ($f_1 = 0.01$, $f_2 = 0.2$), 101 taps.
68+ Checkable constant: $\lVert h \rVert_2 = 0.60216$, so the filtered std is
69+ $\sigma_y = 3.011$ at $\sigma = 5$.
70+
71+### 1.3 The approximate formula being tested (for comparison only)
72+
73+The app's closed-form estimate, to be reimplemented for the comparison
74+plots:
75+
76+$$S_z(f) = S(f) + \sigma_q^2, \qquad
77+\sigma_e = 2^{\int_0^{1/2} \log_2 S_z(f)\,df}, \qquad
78+R_{\text{approx}} = H_\Delta(\sigma_e)$$
79+
80+with $\sigma_q^2 = 1/12$ (no dither) or $1/6$ (dither), the integral by
81+midpoint rule (8192 points suffices), and $H_\Delta(s)$ the exact entropy of
82+$\operatorname{round}(\mathcal N(0,s^2))$ on the unit lattice:
83+$H_\Delta(s) = -\sum_z p_z \log_2 p_z$,
84+$p_z = \Phi(\tfrac{z+1/2}{s}) - \Phi(\tfrac{z-1/2}{s})$.
85+
86+Known inexactness (details in the analysis that motivated this work):
87+
88+- **It is not a bound in either direction.**
89+- In the i.i.d. case (`none`) the truth is exactly $H_\Delta(\sigma)$, while
90+ the formula returns $H_\Delta(\sqrt{\sigma^2 + 1/12})$ — an overestimate
91+ (a double-count of the quantization variance): +0.0024 bits at σ = 5,
92+ +0.171 at σ = 0.5, +0.576 at σ = 0.1 (where the truth is ≈ 0).
93+- In the stopband-dominated regime (default bandpass, σ = 5) it lands
94+ remarkably close to achievable rates, plausibly partly by cancellation of
95+ errors; the MC estimator exists to resolve this.
96+
97+---
98+
99+## 2. Why Monte Carlo can be exact here
100+
101+Three facts combine:
102+
103+1. **The law of $z$ is fully known**, so conditional probabilities can be
104+ *computed*, not estimated from frequencies. For a window
105+ $z_{1:m}$, the event probability (no dither) is a Gaussian rectangle
106+ probability:
107+
108+ $$P(z_{1:m}) = \Pr\!\left[\, y_t \in [z_t - \tfrac12,\, z_t + \tfrac12)
109+ \ \forall t \le m \,\right]$$
110+
111+ — an $m$-dimensional box integral of $\mathcal N(0, \Sigma_y)$ with the
112+ banded covariance from $r_y$. (With dither the box edges are shifted by
113+ the $d_t$, which are handled by sampling; §4.4.)
114+
115+2. **Predictive decomposition.**
116+ $-\log_2 P(z_{1:m}) = \sum_t -\log_2 P(z_t \mid z_{1:t-1})$, and by
117+ Shannon–McMillan–Breiman, $-\tfrac1m \log_2 P(z_{1:m}) \to R$ almost
118+ surely along a single simulated realization. Equivalently, averaging the
119+ conditional surprisal $-\log_2 P(z_0 \mid z_{-k:-1})$ over independent
120+ draws estimates $h_k = H(z_0 \mid z_{-k:-1})$, and $h_k \downarrow R$
121+ monotonically.
122+
123+3. **Monotone error structure.** Every systematic error is upward:
124+ - finite window: $h_k \ge R$ for every $k$ (conditioning reduces
125+ entropy);
126+ - inner estimation: an unbiased $\hat P$ gives
127+ $\mathbb E[-\log_2 \hat P] \ge -\log_2 P$ by Jensen, with bias
128+ $\approx \operatorname{Var}(\hat P) / (2 P^2 \ln 2)$, vanishing as the
129+ particle count grows.
130+
131+ So the estimate converges to R **from above**, and the only downward
132+ fluctuation is CLT noise, which gets an error bar.
133+
134+For a two-sided certificate, use the hidden-Markov sandwich
135+(Cover & Thomas, *Elements of Information Theory*, §4.5): $z$ is a function
136+of the Markov state $s_t = (x_{t-L+1},\dots,x_t)$ (plus $d_t$), so
137+
138+$$H\!\left(z_0 \mid z_{-k:-1},\, s_{-k-1}\right) \;\le\; R \;\le\;
139+H\!\left(z_0 \mid z_{-k:-1}\right)$$
140+
141+with the lower bound *increasing* in $k$ and the gap
142+$I(s_{-k-1}; z_0 \mid z_{-k:-1}) \to 0$. In simulation the true state at
143+$-k-1$ is known (we generated it), so the lower bound is estimable by the
144+same machinery (§4.5).
145+
146+---
147+
148+## 3. The estimator in one paragraph
149+
150+Simulate one long realization $z_{1:T}$ from the pipeline. Run a sequential
151+Monte Carlo over the *latent Gaussian path constrained to the observed
152+boxes*: each particle is a sampled path $y^i$ consistent with
153+$z_{1:t}$ so far; at each step the exact one-step predictive probability of
154+the observed box, given the particle's path, is a Gaussian
155+$\Phi$-difference (because $y_t \mid y_{1:t-1}$ is Gaussian with
156+linear-prediction mean and fixed innovation variance); the mean of those
157+$\Phi$-differences across particles is an unbiased estimate of
158+$P(z_t \mid z_{1:t-1})$; the particle is extended by sampling $y_t$ from
159+the corresponding truncated Gaussian (this "optimal one-step proposal" is
160+exactly Genz's separation-of-variables construction for box probabilities,
161+run sequentially with resampling). Average $-\log_2 \hat P(z_t \mid \cdot)$
162+over $t$ after a burn-in; repeat over independent replicates for error
163+bars. That average estimates R.
164+
165+---
166+
167+## 4. Algorithm specification
168+
169+### 4.1 Precomputation
170+
171+1. Build the kernel $h$ (§1.2) and the autocovariance
172+ $r_y(0),\dots,r_y(k_{\max})$ (zero beyond lag $L-1$).
173+2. Run **Levinson–Durbin** on $r_y$ up to order $k_{\max}$ (suggest
174+ $k_{\max} = 4L$, capped at ~1024). Store, for each order $p \le k_{\max}$
175+ — or just for the terminal order — the prediction coefficients
176+ $a^{(p)}_1,\dots,a^{(p)}_p$ and innovation variance $v_p$. After
177+ $p \gtrsim L$ these are effectively converged ($v_p \to$ the
178+ Szegő value for $S$ *without* any noise floor; e.g. for the default
179+ bandpass at σ = 5, $v_\infty \approx 0.0019$, i.e. innovation std
180+ ≈ 0.044). Guard against underflow of $v_p$ for extreme stopbands (floor
181+ at ~1e-30 and warn).
182+
183+ For $t \le k_{\max}$ use the order-$(t-1)$ coefficients; past that, the
184+ terminal ones (steady state).
185+
186+### 4.2 Per-step recursion (undithered)
187+
188+State: particle paths $y^i_{1:t-1}$, $i = 1..N$ (only the last $k_{\max}$
189+values are needed), all with equal weight after resampling.
190+
191+For $t = 1, 2, \dots, T$, with observed integer $z_t$:
192+
193+1. Prediction per particle:
194+ $\mu^i = \sum_{j=1}^{p} a^{(p)}_j\, y^i_{t-j}$, innovation std
195+ $s = \sqrt{v_p}$, where $p = \min(t-1, k_{\max})$.
196+2. Standardized box edges:
197+ $\alpha^i = (z_t - \tfrac12 - \mu^i)/s$,
198+ $\beta^i = (z_t + \tfrac12 - \mu^i)/s$.
199+3. Incremental weight $w^i = \Phi(\beta^i) - \Phi(\alpha^i)$
200+ (log-space; see §5.1).
201+4. **Predictive estimate**
202+ $\hat p_t = \tfrac1N \sum_i w^i$. Record $-\log_2 \hat p_t$. If
203+ $\hat p_t = 0$ (all particles incompatible), the run has degenerated —
204+ restart the replicate with larger N; do not clamp.
205+5. Resample particles with probabilities $\propto w^i$ (systematic
206+ resampling; adaptive — only when $\mathrm{ESS} < N/2$ — is fine, but then
207+ $\hat p_t$ must use the standard normalized-weight form
208+ $\hat p_t = \sum_i W^i_{t-1} w^i_t$ with carried weights $W$).
209+6. Extend each surviving particle:
210+ $y^i_t = \mu^i + s\,\Phi^{-1}\!\big(\Phi(\alpha^i) + U^i\,(\Phi(\beta^i)-\Phi(\alpha^i))\big)$,
211+ $U^i \sim \mathcal U(0,1)$ — the truncated-Gaussian draw (§5.2).
212+
213+### 4.3 Assembling the estimate
214+
215+- $\hat R = -\dfrac{1}{T - B}\sum_{t=B+1}^{T} \log_2 \hat p_t$, with burn-in
216+ $B \ge 5L$ (the early steps estimate $h_{t-1}$ for small $t$, which is
217+ above R).
218+- Run $J$ independent replicates (fresh $z$, fresh particles). Report
219+ mean ± SE across replicates. Within-run block averaging is acceptable but
220+ replicates are simpler and honest about autocorrelation.
221+- **Bias control (essential):** repeat at $N$, $2N$, $4N$ particles. The
222+ estimate must decrease and plateau (Jensen bias shrinks like $1/N$).
223+ Treat the plateau as the answer; optionally Richardson-extrapolate in
224+ $1/N$. A per-step delta correction
225+ $-\log_2\hat p_t \to -\log_2 \hat p_t - \widehat{\operatorname{Var}}_i(w^i) / (2 N \hat p_t^2 \ln 2)$
226+ is a useful diagnostic but is not exact under resampling — the
227+ particle-doubling plateau is authoritative.
228+
229+### 4.4 Dither
230+
231+With dither, the box for $y_t$ is shifted by the (unknown to the decoder)
232+$d_t$: $y_t \in [z_t - \tfrac12 - d_t,\, z_t + \tfrac12 - d_t)$. Extend each
233+particle with its own $d^i_t \sim \mathcal U[-\tfrac12,\tfrac12)$ drawn
234+*before* step 2, and use the shifted edges. Everything else is unchanged;
235+$\hat p_t$ remains unbiased (the dither prior is part of the proposal).
236+Optionally, integrate $d_t$ analytically in the weight — the integral of a
237+$\Phi$-difference over a unit shift is a difference of
238+$G(t) = t\Phi(t) + \varphi(t)$ terms — and then sample $(d_t, y_t)$ jointly;
239+lower variance, more code. Start with the sampled version.
240+
241+Note: when generating the *observed* $z_{1:T}$ for a dithered run, the
242+generator also draws $d_t$; those true dither values are **not** given to
243+the estimator (non-subtractive dither).
244+
245+### 4.5 Sandwich lower bound (optional but recommended)
246+
247+Estimate $H(z_0 \mid z_{-k:-1}, s_{-k-1})$: for each of many independent
248+windows, (a) simulate truth and keep the exact latent inputs up to time
249+$-k-1$; (b) condition the Gaussian law of $y_{-k:0}$ on that known state —
250+its conditional mean is the deterministic tail response
251+$\sum_{j > t+k} h_j x_{t-j}$ and its conditional covariance is that of the
252+*truncated* kernel (only the taps applied to inputs after $-k-1$); (c) run
253+the same SOV/SMC over the $k{+}1$ constrained steps and record the last
254+step's $-\log_2 \hat p$. Average over windows. Plot lower and upper curves
255+against $k$; the closing gap certifies memory-length convergence
256+independently of any modeling argument.
257+
258+---
259+
260+## 5. Numerical hazards (all known, all manageable)
261+
262+### 5.1 Φ-differences in the tails
263+
264+When both edges are far in one tail, $\Phi(\beta)-\Phi(\alpha)$ underflows
265+in the naive form. Compute in log space:
266+use `scipy.special.log_ndtr`; for $\alpha, \beta > 0$ switch to the upper
267+tail, $\log(\Phi(-\alpha) - \Phi(-\beta))$, and combine as
268+`log_ndtr(hi) + log1p(-exp(log_ndtr(lo) - log_ndtr(hi)))`. These
269+low-probability steps are precisely the ones that dominate the surprisal,
270+so they must not be clamped to zero.
271+
272+### 5.2 Truncated-normal sampling
273+
274+The inverse-CDF trick in §4.2(6) loses precision when the interval sits
275+beyond ~6σ. Use `scipy.stats.truncnorm` (which handles tails), or Robert's
276+exponential-rejection sampler for one-sided extreme tails. A particle
277+sitting exactly on an edge after sampling is harmless.
278+
279+### 5.3 Weight degeneracy and path collapse
280+
281+The innovation std $s$ can be much smaller than the box width (deep
282+stopbands, no dither). Then most particles get $w^i \approx 1$ and the rest
283+$\approx 0$: weights are Bernoulli-like with success probability
284+$\approx 2^{-R}$ per step. This is workable — the truncated draw in
285+step 6 spans the box, i.e. many innovation-σ's, so resampled duplicates
286+re-diversify quickly — but it costs particles. If ESS collapses persistently:
287+
288+- raise N (first resort; the estimator is embarrassingly parallel across
289+ particles and replicates);
290+- use a lookahead / auxiliary proposal (weight by the next few
291+ observations before resampling);
292+- as a fallback formulation, run the SMC in **x-space** (particles over the
293+ i.i.d. inputs, observation $z_t$ constraining
294+ $\sigma \sum_j h_j x_{t-j}$). Beware: after causal reindexing, windowed-
295+ sinc kernels have tiny leading taps, which makes the naive one-step
296+ x-space proposal much worse than the y-space one specified here. x-space
297+ only becomes attractive with block/lookahead proposals.
298+- dithered runs are *easier* (the dither smooths the likelihood); debug
299+ there first.
300+
301+### 5.4 Sanity invariants to assert in code
302+
303+- $\hat p_t \le 1$ always; running $\hat R$ finite.
304+- Undithered, `none` kernel, any σ: the SMC reduces to i.i.d. draws and
305+ $\hat p_t$ must equal $\Phi$-difference of the marginal for every
306+ particle identically (zero variance across particles).
307+- Generated $z$ statistics: sample std of $z \approx \sqrt{\sigma^2\lVert h\rVert_2^2 + \sigma_q^2}$
308+ for $\sigma_y \gtrsim 1$.
309+
310+---
311+
312+## 6. Test plan
313+
314+Ordered from cheap-and-exact to expensive-and-comparative.
315+
316+1. **i.i.d. exactness (primary).** `none` kernel. True
317+ $R = H_\Delta(\sigma)$ in closed form. Require agreement within 3 SE at
318+ σ ∈ {0.5, 1, 5}, and that the reimplemented $R_{\text{approx}}$
319+ reproduces its known overestimates (§7 table).
320+2. **Marginal check.** $h_0$ (no conditioning) equals the marginal quantized
321+ entropy for any kernel: compare a short-window run against the
322+ 1-D formula with $s = \sqrt{r_y(0)}$ (undithered).
323+3. **Brute force, small windows.** For MA(1)-like kernels
324+ ($h = [1, a]$, σ ≤ 2), compute $h_k$ for $k \le 3$ exactly:
325+ enumerate integer tuples within ±8σ_y and evaluate box probabilities
326+ with `scipy.stats.multivariate_normal.cdf` (Genz), then
327+ $h_k = H_{k+1} - H_k$. The SMC at the same $k$ (fresh independent
328+ windows, fixed $k$) must match within error bars.
329+4. **Monotonicity.** $\hat h_k$ non-increasing in $k$ (within noise);
330+ estimate decreasing and plateauing in particle count N; sandwich gap
331+ (if implemented) shrinking in $k$.
332+5. **Dead zone.** `none`, σ = 0.1: $\hat R \approx 0$ (truth ~1e-5), vs
333+ $R_{\text{approx}} = 0.576$. With dither, σ → 0: truth → 0 as well
334+ (uniform ±½ dither alone never flips the integer), vs
335+ $R_{\text{approx}} \approx 0.56$.
336+6. **Achievability cross-check.** For the default bandpass, σ = 5:
337+ $\hat R$ must be ≤ the achievable rates in §7 (LPC residual entropy
338+ ≈ 1.947) plus noise, and plausibly close to $R_{\text{approx}} = 1.940$.
339+ This is the headline number this whole exercise exists to pin down.
340+7. **Sweep deliverable.** σ ∈ {0.5, 1, 2, 5, 10, 20, 50, 100} × {dither
341+ on, off} for the default bandpass: plot $\hat R$ with error bars against
342+ $R_{\text{approx}}$, plus the i.i.d. sweep where truth is closed-form.
343+ Deliver the numbers as CSV alongside the plot.
344+
345+Suggested starting parameters: N = 4096 particles, T = 4000 steps,
346+burn-in B = 1000, J = 16 replicates, $k_{\max} = 512$ for the 101-tap
347+kernel. Rough cost: O(N·k_max) per step ⇒ ~10⁹–10¹⁰ flops per replicate —
348+seconds to a minute in vectorized numpy (vectorize across particles; the
349+prediction is a matrix–vector product against the shared coefficient
350+vector).
351+
352+## 7. Reference values (default bandpass = 300–6000 Hz @ 30 kHz, 101 taps)
353+
354+Closed-form / app-measured values for validating reimplementations. "Exact"
355+rows are analytic; others were measured with the app's codecs on a
356+100k-sample block (seed-dependent in the third decimal).
357+
358+| Quantity | Value | Status |
359+|---|---|---|
360+| $H_\Delta(s)$ at s = 0.1 / 0.3 / 1 / 5 / 50 | 0.00001 / 0.55042 / 2.10483 / 4.37142 / 7.69098 | exact |
361+| True R, `none`, σ = 5 / 0.5 / 0.1 | 4.37142 / 1.24174 / ≈1e-5 | exact ($=H_\Delta(\sigma)$) |
362+| $R_{\text{approx}}$, `none`, σ = 0.1 / 1 / 5 / 100 | 0.57611 / 2.15829 / 4.37382 / 8.69096 | exact given §1.3 |
363+| $R_{\text{approx}}$, bandpass, σ = 0.5 / 5 / 20 / 100 | 0.8620 / 1.9397 / 2.7282 / 3.7166 | exact given §1.3 |
364+| $R_{\text{approx}}$, bandpass + dither, σ = 5 | 2.2111 | exact given §1.3 |
365+| $\lVert h \rVert_2$, default bandpass | 0.60216 | exact given §1.2 |
366+| Order-0 entropy of LPC(32) / LPC(128) residual, bandpass σ = 5 | 1.965 / 1.947 | measured (achievable ⇒ upper bounds on R) |
367+| LPC(32)+ANS achieved, bandpass σ = 5 | 2.002–2.011 | measured |
368+| Pure high-res rate (no floor), bandpass σ = 5 | −2.4604 | exact; demonstrates why the floor exists |
369+
370+## 8. References
371+
372+- A. Genz, "Numerical computation of multivariate normal probabilities,"
373+ *J. Comput. Graph. Statist.*, 1992 — the separation-of-variables
374+ construction that §4.2 runs sequentially.
375+- A. Genz, F. Bretz, *Computation of Multivariate Normal and t
376+ Probabilities*, Springer, 2009.
377+- Z. I. Botev, "The normal law under linear restrictions: simulation and
378+ estimation via minimax tilting," *JRSS-B*, 2017 — variance reduction if
379+ §5.3 becomes limiting.
380+- T. M. Cover, J. A. Thomas, *Elements of Information Theory*, 2nd ed.:
381+ §4.5 (HMM entropy-rate sandwich), §16.8 / AEP (Shannon–McMillan–Breiman).
382+- Companion theory note in the sibling repository
383+ `entropy-quantized-linear-gaussian` (`entropy_quantized_linear_gaussian.md`):
384+ the block-entropy version of the same box-probability machinery, including
385+ the inner-bias caveat this spec's §4.3 addresses.