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.
17---
19## 1. The system under study
21### 1.1 Generating model
23$$x_n \sim \mathcal{N}(0,\sigma^2)\ \text{i.i.d.}, \qquad
24y = h * x, \qquad
25z_n = \operatorname{round}(y_n + d_n)$$
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
35 $$r_y(\ell) = \sigma^2 \sum_{j} h_j\, h_{j+\ell}
36 \qquad (r_y(\ell)=0 \text{ for } |\ell| \ge L)$$
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.
41The object of interest is the **entropy rate of the integer process** $z$:
43$$R \;=\; \lim_{k\to\infty} H\!\left(z_0 \mid z_{-1},\dots,z_{-k}\right)
44\qquad \text{[bits/sample]}$$
46This is the true lossless compression limit for $z$, the number the app's
47LPC + ANS bars try to approach from above.
49### 1.2 Kernel designs (must match to compare numbers)
51All comparisons in §7 use the **bandpass** kernel below; implement at least
52`none` and `bandpass`.
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$:
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)$$
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$.
71### 1.3 The approximate formula being tested (for comparison only)
73The app's closed-form estimate, to be reimplemented for the comparison
74plots:
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
78R_{\text{approx}} = H_\Delta(\sigma_e)$$
80with $\sigma_q^2 = 1/12$ (no dither) or $1/6$ (dither), the integral by
81midpoint 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})$.
86Known inexactness (details in the analysis that motivated this work):
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.
97---
99## 2. Why Monte Carlo can be exact here
101Three facts combine:
1031. **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:
108 $$P(z_{1:m}) = \Pr\!\left[\, y_t \in [z_t - \tfrac12,\, z_t + \tfrac12)
109 \ \forall t \le m \,\right]$$
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.)
1152. **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.
1233. **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.
131 So the estimate converges to R **from above**, and the only downward
132 fluctuation is CLT noise, which gets an error bar.
134For a two-sided certificate, use the hidden-Markov sandwich
135(Cover & Thomas, *Elements of Information Theory*, §4.5): $z$ is a function
136of the Markov state $s_t = (x_{t-L+1},\dots,x_t)$ (plus $d_t$), so
138$$H\!\left(z_0 \mid z_{-k:-1},\, s_{-k-1}\right) \;\le\; R \;\le\;
139H\!\left(z_0 \mid z_{-k:-1}\right)$$
141with 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
144same machinery (§4.5).
146---
148## 3. The estimator in one paragraph
150Simulate one long realization $z_{1:T}$ from the pipeline. Run a sequential
151Monte Carlo over the *latent Gaussian path constrained to the observed
152boxes*: 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
154the observed box, given the particle's path, is a Gaussian
155$\Phi$-difference (because $y_t \mid y_{1:t-1}$ is Gaussian with
156linear-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
159the corresponding truncated Gaussian (this "optimal one-step proposal" is
160exactly Genz's separation-of-variables construction for box probabilities,
161run sequentially with resampling). Average $-\log_2 \hat P(z_t \mid \cdot)$
162over $t$ after a burn-in; repeat over independent replicates for error
163bars. That average estimates R.
165---
167## 4. Algorithm specification
169### 4.1 Precomputation
1711. Build the kernel $h$ (§1.2) and the autocovariance
172 $r_y(0),\dots,r_y(k_{\max})$ (zero beyond lag $L-1$).
1732. 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).
183 For $t \le k_{\max}$ use the order-$(t-1)$ coefficients; past that, the
184 terminal ones (steady state).
186### 4.2 Per-step recursion (undithered)
188State: particle paths $y^i_{1:t-1}$, $i = 1..N$ (only the last $k_{\max}$
189values are needed), all with equal weight after resampling.
191For $t = 1, 2, \dots, T$, with observed integer $z_t$:
1931. 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})$.
1962. Standardized box edges:
197 $\alpha^i = (z_t - \tfrac12 - \mu^i)/s$,
198 $\beta^i = (z_t + \tfrac12 - \mu^i)/s$.
1993. Incremental weight $w^i = \Phi(\beta^i) - \Phi(\alpha^i)$
200 (log-space; see §5.1).
2014. **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.
2055. 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$).
2096. 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).
213### 4.3 Assembling the estimate
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.
229### 4.4 Dither
231With 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
233particle 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).
236Optionally, 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;
239lower variance, more code. Start with the sampled version.
241Note: when generating the *observed* $z_{1:T}$ for a dithered run, the
242generator also draws $d_t$; those true dither values are **not** given to
243the estimator (non-subtractive dither).
245### 4.5 Sandwich lower bound (optional but recommended)
247Estimate $H(z_0 \mid z_{-k:-1}, s_{-k-1})$: for each of many independent
248windows, (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 —
250its 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
253the same SOV/SMC over the $k{+}1$ constrained steps and record the last
254step's $-\log_2 \hat p$. Average over windows. Plot lower and upper curves
255against $k$; the closing gap certifies memory-length convergence
256independently of any modeling argument.
258---
260## 5. Numerical hazards (all known, all manageable)
262### 5.1 Φ-differences in the tails
264When both edges are far in one tail, $\Phi(\beta)-\Phi(\alpha)$ underflows
265in the naive form. Compute in log space:
266use `scipy.special.log_ndtr`; for $\alpha, \beta > 0$ switch to the upper
267tail, $\log(\Phi(-\alpha) - \Phi(-\beta))$, and combine as
268`log_ndtr(hi) + log1p(-exp(log_ndtr(lo) - log_ndtr(hi)))`. These
269low-probability steps are precisely the ones that dominate the surprisal,
270so they must not be clamped to zero.
272### 5.2 Truncated-normal sampling
274The inverse-CDF trick in §4.2(6) loses precision when the interval sits
275beyond ~6σ. Use `scipy.stats.truncnorm` (which handles tails), or Robert's
276exponential-rejection sampler for one-sided extreme tails. A particle
277sitting exactly on an edge after sampling is harmless.
279### 5.3 Weight degeneracy and path collapse
281The innovation std $s$ can be much smaller than the box width (deep
282stopbands, 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
285step 6 spans the box, i.e. many innovation-σ's, so resampled duplicates
286re-diversify quickly — but it costs particles. If ESS collapses persistently:
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.
301### 5.4 Sanity invariants to assert in code
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$.
310---
312## 6. Test plan
314Ordered from cheap-and-exact to expensive-and-comparative.
3161. **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).
3202. **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).
3233. **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.
3294. **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$.
3325. **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$.
3366. **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.
3407. **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.
345Suggested starting parameters: N = 4096 particles, T = 4000 steps,
346burn-in B = 1000, J = 16 replicates, $k_{\max} = 512$ for the 101-tap
347kernel. Rough cost: O(N·k_max) per step ⇒ ~10⁹–10¹⁰ flops per replicate —
348seconds to a minute in vectorized numpy (vectorize across particles; the
349prediction is a matrix–vector product against the shared coefficient
350vector).
352## 7. Reference values (default bandpass = 300–6000 Hz @ 30 kHz, 101 taps)
354Closed-form / app-measured values for validating reimplementations. "Exact"
355rows are analytic; others were measured with the app's codecs on a
356100k-sample block (seed-dependent in the third decimal).
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 |
370## 8. References
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.