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