Monte-Carlo ground-truth script for R, with the command printed in the UI
7 changed files+366−803
README.mdmodified+5−5View file
@@ -32,11 +32,7 @@ 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. 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.
35+math section is a stub for the full derivation.
4036
4137 ## Run it
4238
@@ -56,6 +52,10 @@ src/compress/ lossless codecs run in the browser: zlib (fflate), zstd (wasm),
5652 integer LPC; borrowed from entropy-quantized-linear-transform
5753 src/worker/ the codecs run off the main thread on a debounced parameter set
5854 src/components/ controls, filter plots, signal canvas, compression chart
55+scripts/ true_rate.py — Monte-Carlo ground truth for R (Gibbs over the
56+ latent Gaussian given the observed integers, Rao-Blackwellized
57+ next-sample pmf); the app prints the exact command to run.
58+ Requires numpy only.
5959 ```
6060
6161 Every reported size round-trips through the decoder and includes whatever the
docs/mc-true-rate.mddeleted+0−397View file
@@ -1,403 +0,0 @@
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-**Implementation status.** A reference implementation exists at
18-[`scripts/true_rate.py`](../scripts/true_rate.py) (`--selftest` runs §6's
19-cheap tests; the UI shows the ready-to-run command for its current
20-parameters). It validates on `none` (matches the closed form),
21-`first-difference`, `moving-average`, and windowed-sinc kernels up to ~21
22-taps. For longer windowed-sinc kernels the innovation std falls to a few
23-hundredths of a step and the plain sequential filter hits the §5.3
24-genealogical collapse — the script detects this and aborts rather than
25-reporting garbage. Extending to that regime needs lookahead/twisted
26-proposals (Guarniero, Johansen & Lee, "The iterated auxiliary particle
27-filter", JASA 2017) or exact-HMC rejuvenation for truncated Gaussians
28-(Pakman & Paninski 2014); single-site Gibbs rejuvenation is *not* a viable
29-substitute here because the near-deterministic dynamics make its
30-conditionals microscopically narrow. Early healthy-regime results: the
31-closed-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
33-analysis.
34-
35-
36-## 1. The system under study
37-
38-### 1.1 Generating model
39-
40-$$x_n \sim \mathcal{N}(0,\sigma^2)\ \text{i.i.d.}, \qquad
41-y = h * x, \qquad
42-z_n = \operatorname{round}(y_n + d_n)$$
43-
44-- $h_0,\dots,h_{L-1}$ is a known FIR kernel (designs in §1.2).
45-- $d_n \sim \mathcal{U}[-\tfrac12,\tfrac12)$ i.i.d. when **dither** is on;
46- $d_n = 0$ otherwise. Dither is *non-subtractive*: the stored integers
47- include it and the decoder never learns $d$.
48-- Quantization step = 1 (σ is measured in steps). Rounding ties are a
49- measure-zero event, so the tie convention does not affect any law.
50-- $y$ is a stationary Gaussian process with autocovariance
51-
52- $$r_y(\ell) = \sigma^2 \sum_{j} h_j\, h_{j+\ell}
53- \qquad (r_y(\ell)=0 \text{ for } |\ell| \ge L)$$
54-
55- and spectral density $S(f) = \sigma^2 |H(f)|^2$,
56- $H(f) = \sum_n h_n e^{-2\pi i f n}$, $f \in [0,\tfrac12]$ cycles/sample.
57-
58-The object of interest is the **entropy rate of the integer process** $z$:
59-
60-$$R \;=\; \lim_{k\to\infty} H\!\left(z_0 \mid z_{-1},\dots,z_{-k}\right)
61-\qquad \text{[bits/sample]}$$
62-
63-This is the true lossless compression limit for $z$, the number the app's
64-LPC + ANS bars try to approach from above.
65-
66-### 1.2 Kernel designs (must match to compare numbers)
67-
68-All comparisons in §7 use the **bandpass** kernel below; implement at least
69-`none` and `bandpass`.
70-
71-- `none`: $h = [1]$ (then $z$ is i.i.d. and everything has closed form —
72- the primary validation case).
73-- Hamming-windowed sinc lowpass, cutoff $f_c$ cycles/sample, odd length $n$,
74- midpoint $m = (n-1)/2$:
75-
76- $$h_i = w_i \cdot \begin{cases} 2f_c & i = m \\
77- \sin(2\pi f_c (i-m)) / (\pi (i-m)) & i \ne m \end{cases},
78- \qquad w_i = 0.54 - 0.46\cos\!\left(\tfrac{2\pi i}{n-1}\right)$$
79-
80- then normalized to unit DC gain: $h \leftarrow h / \sum_i h_i$.
81-- `bandpass(f_1, f_2)` = lowpass($f_2$) − lowpass($f_1$), same length, no
82- further normalization.
83-- The **default ephys case**: bandpass 300–6000 Hz at 30 kHz sample rate
84- ($f_1 = 0.01$, $f_2 = 0.2$), 101 taps.
85- Checkable constant: $\lVert h \rVert_2 = 0.60216$, so the filtered std is
86- $\sigma_y = 3.011$ at $\sigma = 5$.
87-
88-### 1.3 The approximate formula being tested (for comparison only)
89-
90-The app's closed-form estimate, to be reimplemented for the comparison
91-plots:
92-
93-$$S_z(f) = S(f) + \sigma_q^2, \qquad
94-\sigma_e = 2^{\int_0^{1/2} \log_2 S_z(f)\,df}, \qquad
95-R_{\text{approx}} = H_\Delta(\sigma_e)$$
96-
97-with $\sigma_q^2 = 1/12$ (no dither) or $1/6$ (dither), the integral by
98-midpoint rule (8192 points suffices), and $H_\Delta(s)$ the exact entropy of
99-$\operatorname{round}(\mathcal N(0,s^2))$ on the unit lattice:
100-$H_\Delta(s) = -\sum_z p_z \log_2 p_z$,
101-$p_z = \Phi(\tfrac{z+1/2}{s}) - \Phi(\tfrac{z-1/2}{s})$.
102-
103-Known inexactness (details in the analysis that motivated this work):
104-
105-- **It is not a bound in either direction.**
106-- In the i.i.d. case (`none`) the truth is exactly $H_\Delta(\sigma)$, while
107- the formula returns $H_\Delta(\sqrt{\sigma^2 + 1/12})$ — an overestimate
108- (a double-count of the quantization variance): +0.0024 bits at σ = 5,
109- +0.171 at σ = 0.5, +0.576 at σ = 0.1 (where the truth is ≈ 0).
110-- In the stopband-dominated regime (default bandpass, σ = 5) it lands
111- remarkably close to achievable rates, plausibly partly by cancellation of
112- errors; the MC estimator exists to resolve this.
113-
114-
115-## 2. Why Monte Carlo can be exact here
116-
117-Three facts combine:
118-
119-1. **The law of $z$ is fully known**, so conditional probabilities can be
120- *computed*, not estimated from frequencies. For a window
121- $z_{1:m}$, the event probability (no dither) is a Gaussian rectangle
122- probability:
123-
124- $$P(z_{1:m}) = \Pr\!\left[\, y_t \in [z_t - \tfrac12,\, z_t + \tfrac12)
125- \ \forall t \le m \,\right]$$
126-
127- — an $m$-dimensional box integral of $\mathcal N(0, \Sigma_y)$ with the
128- banded covariance from $r_y$. (With dither the box edges are shifted by
129- the $d_t$, which are handled by sampling; §4.4.)
130-
131-2. **Predictive decomposition.**
132- $-\log_2 P(z_{1:m}) = \sum_t -\log_2 P(z_t \mid z_{1:t-1})$, and by
133- Shannon–McMillan–Breiman, $-\tfrac1m \log_2 P(z_{1:m}) \to R$ almost
134- surely along a single simulated realization. Equivalently, averaging the
135- conditional surprisal $-\log_2 P(z_0 \mid z_{-k:-1})$ over independent
136- draws estimates $h_k = H(z_0 \mid z_{-k:-1})$, and $h_k \downarrow R$
137- monotonically.
138-
139-3. **Monotone error structure.** Every systematic error is upward:
140- - finite window: $h_k \ge R$ for every $k$ (conditioning reduces
141- entropy);
142- - inner estimation: an unbiased $\hat P$ gives
143- $\mathbb E[-\log_2 \hat P] \ge -\log_2 P$ by Jensen, with bias
144- $\approx \operatorname{Var}(\hat P) / (2 P^2 \ln 2)$, vanishing as the
145- particle count grows.
146-
147- So the estimate converges to R **from above**, and the only downward
148- fluctuation is CLT noise, which gets an error bar.
149-
150-For a two-sided certificate, use the hidden-Markov sandwich
151-(Cover & Thomas, *Elements of Information Theory*, §4.5): $z$ is a function
152-of the Markov state $s_t = (x_{t-L+1},\dots,x_t)$ (plus $d_t$), so
153-
154-$$H\!\left(z_0 \mid z_{-k:-1},\, s_{-k-1}\right) \;\le\; R \;\le\;
155-H\!\left(z_0 \mid z_{-k:-1}\right)$$
156-
157-with the lower bound *increasing* in $k$ and the gap
158-$I(s_{-k-1}; z_0 \mid z_{-k:-1}) \to 0$. In simulation the true state at
159-$-k-1$ is known (we generated it), so the lower bound is estimable by the
160-same machinery (§4.5).
161-
162-
163-## 3. The estimator in one paragraph
164-
165-Simulate one long realization $z_{1:T}$ from the pipeline. Run a sequential
166-Monte Carlo over the *latent Gaussian path constrained to the observed
167-boxes*: each particle is a sampled path $y^i$ consistent with
168-$z_{1:t}$ so far; at each step the exact one-step predictive probability of
169-the observed box, given the particle's path, is a Gaussian
170-$\Phi$-difference (because $y_t \mid y_{1:t-1}$ is Gaussian with
171-linear-prediction mean and fixed innovation variance); the mean of those
172-$\Phi$-differences across particles is an unbiased estimate of
173-$P(z_t \mid z_{1:t-1})$; the particle is extended by sampling $y_t$ from
174-the corresponding truncated Gaussian (this "optimal one-step proposal" is
175-exactly Genz's separation-of-variables construction for box probabilities,
176-run sequentially with resampling). Average $-\log_2 \hat P(z_t \mid \cdot)$
177-over $t$ after a burn-in; repeat over independent replicates for error
178-bars. That average estimates R.
179-
180-
181-## 4. Algorithm specification
182-
183-### 4.1 Precomputation
184-
185-1. Build the kernel $h$ (§1.2) and the autocovariance
186- $r_y(0),\dots,r_y(k_{\max})$ (zero beyond lag $L-1$).
187-2. Run **Levinson–Durbin** on $r_y$ up to order $k_{\max}$ (suggest
188- $k_{\max} = 4L$, capped at ~1024). Store, for each order $p \le k_{\max}$
189- — or just for the terminal order — the prediction coefficients
190- $a^{(p)}_1,\dots,a^{(p)}_p$ and innovation variance $v_p$. After
191- $p \gtrsim L$ these are effectively converged ($v_p \to$ the
192- Szegő value for $S$ *without* any noise floor; e.g. for the default
193- bandpass at σ = 5, $v_\infty \approx 0.0019$, i.e. innovation std
194- ≈ 0.044). Guard against underflow of $v_p$ for extreme stopbands (floor
195- at ~1e-30 and warn).
196-
197- For $t \le k_{\max}$ use the order-$(t-1)$ coefficients; past that, the
198- terminal ones (steady state).
199-
200-### 4.2 Per-step recursion (undithered)
201-
202-State: particle paths $y^i_{1:t-1}$, $i = 1..N$ (only the last $k_{\max}$
203-values are needed), all with equal weight after resampling.
204-
205-For $t = 1, 2, \dots, T$, with observed integer $z_t$:
206-
207-1. Prediction per particle:
208- $\mu^i = \sum_{j=1}^{p} a^{(p)}_j\, y^i_{t-j}$, innovation std
209- $s = \sqrt{v_p}$, where $p = \min(t-1, k_{\max})$.
210-2. Standardized box edges:
211- $\alpha^i = (z_t - \tfrac12 - \mu^i)/s$,
212- $\beta^i = (z_t + \tfrac12 - \mu^i)/s$.
213-3. Incremental weight $w^i = \Phi(\beta^i) - \Phi(\alpha^i)$
214- (log-space; see §5.1).
215-4. **Predictive estimate**
216- $\hat p_t = \tfrac1N \sum_i w^i$. Record $-\log_2 \hat p_t$. If
217- $\hat p_t = 0$ (all particles incompatible), the run has degenerated —
218- restart the replicate with larger N; do not clamp.
219-5. Resample particles with probabilities $\propto w^i$ (systematic
220- resampling; adaptive — only when $\mathrm{ESS} < N/2$ — is fine, but then
221- $\hat p_t$ must use the standard normalized-weight form
222- $\hat p_t = \sum_i W^i_{t-1} w^i_t$ with carried weights $W$).
223-6. Extend each surviving particle:
224- $y^i_t = \mu^i + s\,\Phi^{-1}\!\big(\Phi(\alpha^i) + U^i\,(\Phi(\beta^i)-\Phi(\alpha^i))\big)$,
225- $U^i \sim \mathcal U(0,1)$ — the truncated-Gaussian draw (§5.2).
226-
227-### 4.3 Assembling the estimate
228-
229-- $\hat R = -\dfrac{1}{T - B}\sum_{t=B+1}^{T} \log_2 \hat p_t$, with burn-in
230- $B \ge 5L$ (the early steps estimate $h_{t-1}$ for small $t$, which is
231- above R).
232-- Run $J$ independent replicates (fresh $z$, fresh particles). Report
233- mean ± SE across replicates. Within-run block averaging is acceptable but
234- replicates are simpler and honest about autocorrelation.
235-- **Bias control (essential):** repeat at $N$, $2N$, $4N$ particles. The
236- estimate must decrease and plateau (Jensen bias shrinks like $1/N$).
237- Treat the plateau as the answer; optionally Richardson-extrapolate in
238- $1/N$. A per-step delta correction
239- $-\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)$
240- is a useful diagnostic but is not exact under resampling — the
241- particle-doubling plateau is authoritative.
242-
243-### 4.4 Dither
244-
245-With dither, the box for $y_t$ is shifted by the (unknown to the decoder)
246-$d_t$: $y_t \in [z_t - \tfrac12 - d_t,\, z_t + \tfrac12 - d_t)$. Extend each
247-particle with its own $d^i_t \sim \mathcal U[-\tfrac12,\tfrac12)$ drawn
248-*before* step 2, and use the shifted edges. Everything else is unchanged;
249-$\hat p_t$ remains unbiased (the dither prior is part of the proposal).
250-Optionally, integrate $d_t$ analytically in the weight — the integral of a
251-$\Phi$-difference over a unit shift is a difference of
252-$G(t) = t\Phi(t) + \varphi(t)$ terms — and then sample $(d_t, y_t)$ jointly;
253-lower variance, more code. Start with the sampled version.
254-
255-Note: when generating the *observed* $z_{1:T}$ for a dithered run, the
256-generator also draws $d_t$; those true dither values are **not** given to
257-the estimator (non-subtractive dither).
258-
259-### 4.5 Sandwich lower bound (optional but recommended)
260-
261-Estimate $H(z_0 \mid z_{-k:-1}, s_{-k-1})$: for each of many independent
262-windows, (a) simulate truth and keep the exact latent inputs up to time
263-$-k-1$; (b) condition the Gaussian law of $y_{-k:0}$ on that known state —
264-its conditional mean is the deterministic tail response
265-$\sum_{j > t+k} h_j x_{t-j}$ and its conditional covariance is that of the
266-*truncated* kernel (only the taps applied to inputs after $-k-1$); (c) run
267-the same SOV/SMC over the $k{+}1$ constrained steps and record the last
268-step's $-\log_2 \hat p$. Average over windows. Plot lower and upper curves
269-against $k$; the closing gap certifies memory-length convergence
270-independently of any modeling argument.
271-
272-
273-## 5. Numerical hazards (all known, all manageable)
274-
275-### 5.1 Φ-differences in the tails
276-
277-When both edges are far in one tail, $\Phi(\beta)-\Phi(\alpha)$ underflows
278-in the naive form. Compute in log space:
279-use `scipy.special.log_ndtr`; for $\alpha, \beta > 0$ switch to the upper
280-tail, $\log(\Phi(-\alpha) - \Phi(-\beta))$, and combine as
281-`log_ndtr(hi) + log1p(-exp(log_ndtr(lo) - log_ndtr(hi)))`. These
282-low-probability steps are precisely the ones that dominate the surprisal,
283-so they must not be clamped to zero.
284-
285-### 5.2 Truncated-normal sampling
286-
287-The inverse-CDF trick in §4.2(6) loses precision when the interval sits
288-beyond ~6σ. Use `scipy.stats.truncnorm` (which handles tails), or Robert's
289-exponential-rejection sampler for one-sided extreme tails. A particle
290-sitting exactly on an edge after sampling is harmless.
291-
292-### 5.3 Weight degeneracy and path collapse
293-
294-The innovation std $s$ can be much smaller than the box width (deep
295-stopbands, no dither). Then most particles get $w^i \approx 1$ and the rest
296-$\approx 0$: weights are Bernoulli-like with success probability
297-$\approx 2^{-R}$ per step. This is workable — the truncated draw in
298-step 6 spans the box, i.e. many innovation-σ's, so resampled duplicates
299-re-diversify quickly — but it costs particles. If ESS collapses persistently:
300-
301-- raise N (first resort; the estimator is embarrassingly parallel across
302- particles and replicates);
303-- use a lookahead / auxiliary proposal (weight by the next few
304- observations before resampling);
305-- as a fallback formulation, run the SMC in **x-space** (particles over the
306- i.i.d. inputs, observation $z_t$ constraining
307- $\sigma \sum_j h_j x_{t-j}$). Beware: after causal reindexing, windowed-
308- sinc kernels have tiny leading taps, which makes the naive one-step
309- x-space proposal much worse than the y-space one specified here. x-space
310- only becomes attractive with block/lookahead proposals.
311-- dithered runs are *easier* (the dither smooths the likelihood); debug
312- there first.
313-
314-### 5.4 Sanity invariants to assert in code
315-
316-- $\hat p_t \le 1$ always; running $\hat R$ finite.
317-- Undithered, `none` kernel, any σ: the SMC reduces to i.i.d. draws and
318- $\hat p_t$ must equal $\Phi$-difference of the marginal for every
319- particle identically (zero variance across particles).
320-- Generated $z$ statistics: sample std of $z \approx \sqrt{\sigma^2\lVert h\rVert_2^2 + \sigma_q^2}$
321- for $\sigma_y \gtrsim 1$.
322-
323-
324-## 6. Test plan
325-
326-Ordered from cheap-and-exact to expensive-and-comparative.
327-
328-1. **i.i.d. exactness (primary).** `none` kernel. True
329- $R = H_\Delta(\sigma)$ in closed form. Require agreement within 3 SE at
330- σ ∈ {0.5, 1, 5}, and that the reimplemented $R_{\text{approx}}$
331- reproduces its known overestimates (§7 table).
332-2. **Marginal check.** $h_0$ (no conditioning) equals the marginal quantized
333- entropy for any kernel: compare a short-window run against the
334- 1-D formula with $s = \sqrt{r_y(0)}$ (undithered).
335-3. **Brute force, small windows.** For MA(1)-like kernels
336- ($h = [1, a]$, σ ≤ 2), compute $h_k$ for $k \le 3$ exactly:
337- enumerate integer tuples within ±8σ_y and evaluate box probabilities
338- with `scipy.stats.multivariate_normal.cdf` (Genz), then
339- $h_k = H_{k+1} - H_k$. The SMC at the same $k$ (fresh independent
340- windows, fixed $k$) must match within error bars.
341-4. **Monotonicity.** $\hat h_k$ non-increasing in $k$ (within noise);
342- estimate decreasing and plateauing in particle count N; sandwich gap
343- (if implemented) shrinking in $k$.
344-5. **Dead zone.** `none`, σ = 0.1: $\hat R \approx 0$ (truth ~1e-5), vs
345- $R_{\text{approx}} = 0.576$. With dither, σ → 0: truth → 0 as well
346- (uniform ±½ dither alone never flips the integer), vs
347- $R_{\text{approx}} \approx 0.56$.
348-6. **Achievability cross-check.** For the default bandpass, σ = 5:
349- $\hat R$ must be ≤ the achievable rates in §7 (LPC residual entropy
350- ≈ 1.947) plus noise, and plausibly close to $R_{\text{approx}} = 1.940$.
351- This is the headline number this whole exercise exists to pin down.
352-7. **Sweep deliverable.** σ ∈ {0.5, 1, 2, 5, 10, 20, 50, 100} × {dither
353- on, off} for the default bandpass: plot $\hat R$ with error bars against
354- $R_{\text{approx}}$, plus the i.i.d. sweep where truth is closed-form.
355- Deliver the numbers as CSV alongside the plot.
356-
357-Suggested starting parameters: N = 4096 particles, T = 4000 steps,
358-burn-in B = 1000, J = 16 replicates, $k_{\max} = 512$ for the 101-tap
359-kernel. Rough cost: O(N·k_max) per step ⇒ ~10⁹–10¹⁰ flops per replicate —
360-seconds to a minute in vectorized numpy (vectorize across particles; the
361-prediction is a matrix–vector product against the shared coefficient
362-vector).
363-
364-## 7. Reference values (default bandpass = 300–6000 Hz @ 30 kHz, 101 taps)
365-
366-Closed-form / app-measured values for validating reimplementations. "Exact"
367-rows are analytic; others were measured with the app's codecs on a
368-100k-sample block (seed-dependent in the third decimal).
369-
370-| Quantity | Value | Status |
371-|---|---|---|
372-| $H_\Delta(s)$ at s = 0.1 / 0.3 / 1 / 5 / 50 | 0.00001 / 0.55042 / 2.10483 / 4.37142 / 7.69098 | exact |
373-| True R, `none`, σ = 5 / 0.5 / 0.1 | 4.37142 / 1.24174 / ≈1e-5 | exact ($=H_\Delta(\sigma)$) |
374-| $R_{\text{approx}}$, `none`, σ = 0.1 / 1 / 5 / 100 | 0.57611 / 2.15829 / 4.37382 / 8.69096 | exact given §1.3 |
375-| $R_{\text{approx}}$, bandpass, σ = 0.5 / 5 / 20 / 100 | 0.8620 / 1.9397 / 2.7282 / 3.7166 | exact given §1.3 |
376-| $R_{\text{approx}}$, bandpass + dither, σ = 5 | 2.2111 | exact given §1.3 |
377-| $\lVert h \rVert_2$, default bandpass | 0.60216 | exact given §1.2 |
378-| Order-0 entropy of LPC(32) / LPC(128) residual, bandpass σ = 5 | 1.965 / 1.947 | measured (achievable ⇒ upper bounds on R) |
379-| LPC(32)+ANS achieved, bandpass σ = 5 | 2.002–2.011 | measured |
380-| Pure high-res rate (no floor), bandpass σ = 5 | −2.4604 | exact; demonstrates why the floor exists |
381-
382-## 8. References
383-
384-- A. Genz, "Numerical computation of multivariate normal probabilities,"
385- *J. Comput. Graph. Statist.*, 1992 — the separation-of-variables
386- construction that §4.2 runs sequentially.
387-- A. Genz, F. Bretz, *Computation of Multivariate Normal and t
388- Probabilities*, Springer, 2009.
389-- Z. I. Botev, "The normal law under linear restrictions: simulation and
390- estimation via minimax tilting," *JRSS-B*, 2017 — variance reduction if
391- §5.3 becomes limiting.
392-- T. M. Cover, J. A. Thomas, *Elements of Information Theory*, 2nd ed.:
393- §4.5 (HMM entropy-rate sandwich), §16.8 / AEP (Shannon–McMillan–Breiman).
394-- Companion theory note in the sibling repository
395- `entropy-quantized-linear-gaussian` (`entropy_quantized_linear_gaussian.md`):
396- the block-entropy version of the same box-probability machinery, including
397- the inner-bias caveat this spec's §4.3 addresses.
scripts/true_rate.pymodified+289−306View file
@@ -1,335 +1,318 @@
11 #!/usr/bin/env python3
2-"""Monte Carlo estimate of the true entropy rate of the quantized filtered
3-Gaussian pipeline, alongside the closed-form approximation the UI plots.
4-
5-Implements docs/mc-true-rate.md: sequential Monte Carlo (Genz separation of
6-variables with resampling) over the latent Gaussian path constrained to the
7-observed integer boxes. All systematic errors are upward, so the estimate
8-converges to the true R from above; double --particles and compare to check
9-convergence.
10-
11-The UI shows the exact command for the current parameter set. Requires
12-numpy and scipy.
13-
14-Examples:
15- python scripts/true_rate.py --sigma 5 --filter bandpass --low-hz 300 \
16- --high-hz 6000 --taps 101 --sample-rate 30000
17- python scripts/true_rate.py --sigma 0.5 --filter none
18- python scripts/true_rate.py --selftest
2+"""Monte-Carlo ground truth for the theoretical rate R shown in the app.
3+
4+Model (matching the app): x iid N(0, sigma^2) -> y = h * x -> optional
5+additive uniform dither on [-1/2, 1/2) -> z = round(.). The entropy rate
6+H = E[-log2 P(z_next | past)] is the true lossless limit in bits/sample.
7+(The app applies the kernel zero-phase; a time shift does not change the
8+law of the process, so causal convolution is used here.)
9+
10+Method
11+ 1. Draw a past: sample x (and dither) from the prior and push it through
12+ the pipeline to get z_1..z_M.
13+ 2. The generating x is itself an exact draw from p(x | z_1..z_M), so a
14+ Gibbs chain started there is already in stationarity - no burn-in
15+ bias, only autocorrelation.
16+ 3. Gibbs-sample x | z: this posterior is a box-truncated multivariate
17+ normal, and each conditional x_i | rest is N(0, sigma^2) truncated to
18+ an interval read off the <= L constraint boxes x_i appears in. With
19+ dither on, the dither values are extra latents with uniform
20+ conditionals. Coordinates a multiple of L apart share no constraint,
21+ so each of the L "colors" is updated as one vectorized block.
22+ 4. Rao-Blackwellization: given a chain state, the next sample is
23+ z_next = round(c + h_0 * x_free (+ d)) with x_free ~ N(0, sigma^2)
24+ still unconstrained, so P(z_next = j | state) has a closed form.
25+ Averaging these pmfs over the chain gives P(z_next | past) exactly in
26+ the limit; its entropy is the conditional entropy for that past.
27+ 5. Average over independent pasts; report mean +/- standard error.
28+
29+The estimate targets H(z_{M+1} | z_1..z_M), which is an upper bound on the
30+rate and decreases toward it as --past grows beyond the memory of the
31+process. More --sweeps reduces the (downward) plug-in bias from noise in
32+the averaged pmf; narrowband filters and large sigma mix more slowly and
33+deserve more sweeps.
34+
35+Requires numpy only.
1936 """
2037
2138 import argparse
2239 import math
23-import sys
2440
2541 import numpy as np
26-from scipy.special import log_ndtr, logsumexp, ndtr, ndtri
42+
43+SQRT2PI = math.sqrt(2 * math.pi)
44+
2745
2846 # ---------------------------------------------------------------- kernels
47+# Ported from src/model/filters.ts; must stay in step with it.
2948
30-def windowed_sinc_lowpass(fc: float, taps: int) -> np.ndarray:
31- n = taps if taps % 2 == 1 else taps + 1
49+def windowed_sinc_lowpass(fc, taps):
50+ n = taps | 1
51+ mid = (n - 1) / 2
3252 i = np.arange(n)
33- t = i - (n - 1) / 2
34- with np.errstate(invalid="ignore", divide="ignore"):
35- sinc = np.where(t == 0, 2 * fc, np.sin(2 * np.pi * fc * t) / (np.pi * t))
53+ t = i - mid
54+ sinc = np.where(t == 0, 2 * fc, np.sin(2 * np.pi * fc * t) / (np.pi * np.where(t == 0, 1, t)))
3655 w = 0.54 - 0.46 * np.cos(2 * np.pi * i / (n - 1))
3756 h = sinc * w
3857 return h / h.sum()
3958
4059
41-def design_kernel(args) -> np.ndarray:
42- f = args.filter
43- if f == "none":
60+def design_kernel(args):
61+ if args.filter == 'none':
4462 return np.array([1.0])
45- if f == "first-difference":
46- return np.array([1.0, -1.0])
47- if f == "moving-average":
63+ if args.filter == 'moving-average':
4864 return np.full(args.width, 1.0 / args.width)
49- if f == "lowpass":
50- return windowed_sinc_lowpass(args.cutoff_hz / args.sample_rate, args.taps)
51- if f == "bandpass":
52- lo = windowed_sinc_lowpass(args.low_hz / args.sample_rate, args.taps)
53- hi = windowed_sinc_lowpass(args.high_hz / args.sample_rate, args.taps)
65+ if args.filter == 'lowpass':
66+ return windowed_sinc_lowpass(args.high / args.rate, args.taps)
67+ if args.filter == 'bandpass':
68+ lo = windowed_sinc_lowpass(args.low / args.rate, args.taps)
69+ hi = windowed_sinc_lowpass(args.high / args.rate, args.taps)
5470 return hi - lo
55- raise ValueError(f)
71+ if args.filter == 'first-difference':
72+ return np.array([1.0, -1.0])
73+ raise ValueError(args.filter)
5674
57-# ------------------------------------------------- closed-form quantities
5875
59-def h_delta(s: float) -> float:
60- """Exact entropy (bits) of round(N(0, s^2)) on the unit lattice."""
61- if s <= 0:
76+# ------------------------------------------------- the app's formula for R
77+# Ported from src/model/theory.ts. Phi via math.erf (machine precision).
78+
79+def quantized_gaussian_entropy(s):
80+ if s <= 0.02:
6281 return 0.0
63- zmax = int(np.ceil(8 * s + 4))
64- z = np.arange(-zmax, zmax + 1)
65- p = ndtr((z + 0.5) / s) - ndtr((z - 0.5) / s)
66- p = p[p > 1e-300]
67- return float(-(p @ np.log2(p)))
68-
69-
70-def magnitude_sq(h: np.ndarray, f: np.ndarray) -> np.ndarray:
71- n = np.arange(len(h))
72- e = np.exp(-2j * np.pi * np.outer(f, n))
73- H = e @ h
74- return np.abs(H) ** 2
75-
76-
77-def r_approx(h: np.ndarray, sigma: float, dither: bool, points: int = 8192) -> float:
78- """The UI's formula: noise-floored Szego prediction error through H_delta."""
79- floor = 1 / 6 if dither else 1 / 12
80- f = (np.arange(points) + 0.5) * 0.5 / points
81- s_z = sigma**2 * magnitude_sq(h, f) + floor
82- integral = float(np.log2(s_z).mean() * 0.5)
83- return h_delta(2.0**integral)
84-
85-# ----------------------------------------------------------- MC estimator
86-
87-def autocovariance(h: np.ndarray, sigma: float, kmax: int) -> np.ndarray:
88- full = sigma**2 * np.correlate(h, h, "full")
89- r = np.zeros(kmax + 1)
90- take = min(kmax + 1, len(h))
91- r[:take] = full[len(h) - 1 : len(h) - 1 + take]
92- return r
93-
94-
95-def levinson_all(r: np.ndarray, kmax: int):
96- """Prediction coefficients for every order 0..kmax and innovation variances.
97-
98- A[p] holds a_1..a_p with yhat_t = sum_j a_j y_{t-j}; v[p] is the order-p
99- prediction error variance.
100- """
101- A = [np.zeros(0)]
102- v = np.zeros(kmax + 1)
103- v[0] = r[0]
104- a = np.zeros(kmax)
105- for p in range(1, kmax + 1):
106- acc = r[p] - (a[: p - 1] @ r[p - 1 : 0 : -1] if p > 1 else 0.0)
107- k = acc / v[p - 1]
108- if p > 1:
109- a[: p - 1] = a[: p - 1] - k * a[p - 2 :: -1]
110- a[p - 1] = k
111- v[p] = max(v[p - 1] * (1 - k * k), 1e-30)
112- A.append(a[:p].copy())
113- return A, v
114-
115-
116-def log_phi_diff(alpha: np.ndarray, beta: np.ndarray) -> np.ndarray:
117- """log(Phi(beta) - Phi(alpha)) elementwise, safe in both tails."""
118- out = np.empty_like(alpha)
119- hi_tail = alpha >= 0 # reflect to the lower tail
120- a = np.where(hi_tail, -beta, alpha)
121- b = np.where(hi_tail, -alpha, beta)
122- straddle = b >= 0 # a < 0 <= b: safe in linear space
123- with np.errstate(divide="ignore"):
124- out[straddle] = np.log(ndtr(b[straddle]) - ndtr(a[straddle]))
125- lo = ~straddle # both below 0: log-space difference
126- la, lb = log_ndtr(a[lo]), log_ndtr(b[lo])
127- out[lo] = lb + np.log1p(-np.exp(np.minimum(la - lb, -1e-12)))
128- return out
129-
130-
131-def sample_truncated(mu, s, lo, hi, rng):
132- """Sample N(mu, s^2) truncated to [lo, hi], vectorized, tail-safe by
133- reflection. Callers guarantee the interval has nonzero mass."""
134- alpha = (lo - mu) / s
135- beta = (hi - mu) / s
136- flip = alpha > 0
137- a = np.where(flip, -beta, alpha)
138- b = np.where(flip, -alpha, beta)
139- pa, pb = ndtr(a), ndtr(b)
140- u = pa + rng.random(len(mu)) * (pb - pa)
141- x = ndtri(np.clip(u, 1e-320, 1 - 1e-16))
82+ zmax = int(math.ceil(8 * s + 4))
83+ H = 0.0
84+ prev = 0.5 * (1 + math.erf((-zmax - 0.5) / (s * math.sqrt(2))))
85+ for z in range(-zmax, zmax + 1):
86+ cur = 0.5 * (1 + math.erf((z + 0.5) / (s * math.sqrt(2))))
87+ p = cur - prev
88+ prev = cur
89+ if p > 0:
90+ H -= p * math.log2(p)
91+ return H
92+
93+
94+def formula_rate(kernel, sigma, dither, points=8192):
95+ noise_var = 1 / 6 if dither else 1 / 12
96+ f = (0.5 * (np.arange(points) + 0.5)) / points
97+ w = -2j * np.pi * np.outer(f, np.arange(len(kernel)))
98+ S = sigma * sigma * np.abs(np.exp(w) @ kernel) ** 2
99+ integral = float(np.log2(S + noise_var).mean()) * 0.5
100+ return quantized_gaussian_entropy(2.0 ** integral)
101+
102+
103+# ------------------------------------------------ vectorized normal helpers
104+
105+def norm_pdf(t):
106+ return np.exp(-0.5 * t * t) / SQRT2PI
107+
108+
109+def norm_cdf(t):
110+ """Abramowitz-Stegun 26.2.17; |error| < 7.5e-8, plenty for sampling and
111+ for pmf weights whose entropy is wanted to ~1e-4 bits."""
112+ t = np.asarray(t, dtype=float)
113+ z = np.abs(t)
114+ k = 1.0 / (1.0 + 0.2316419 * z)
115+ poly = k * (0.319381530 + k * (-0.356563782 + k * (1.781477937 + k * (-1.821255978 + k * 1.330274429))))
116+ tail = norm_pdf(z) * poly
117+ return np.where(t >= 0, 1.0 - tail, tail)
118+
119+
120+def norm_ppf(p):
121+ """Acklam's rational approximation to the standard normal quantile."""
122+ a = (-3.969683028665376e+01, 2.209460984245205e+02, -2.759285104469687e+02,
123+ 1.383577518672690e+02, -3.066479806614716e+01, 2.506628277459239e+00)
124+ b = (-5.447609879822406e+01, 1.615858368580409e+02, -1.556989798598866e+02,
125+ 6.680131188771972e+01, -1.328068155288572e+01)
126+ c = (-7.784894002430293e-03, -3.223964580411365e-01, -2.400758277161838e+00,
127+ -2.549732539343734e+00, 4.374664141464968e+00, 2.938163982698783e+00)
128+ d = (7.784695709041462e-03, 3.224671290700398e-01, 2.445134137142996e+00,
129+ 3.754408661907416e+00)
130+ p = np.asarray(p, dtype=float)
131+ x = np.empty_like(p)
132+ plow, phigh = 0.02425, 1 - 0.02425
133+
134+ lo = p < plow
135+ hi = p > phigh
136+ mid = ~(lo | hi)
137+
138+ if mid.any():
139+ q = p[mid] - 0.5
140+ r = q * q
141+ x[mid] = (((((a[0] * r + a[1]) * r + a[2]) * r + a[3]) * r + a[4]) * r + a[5]) * q / \
142+ (((((b[0] * r + b[1]) * r + b[2]) * r + b[3]) * r + b[4]) * r + 1)
143+ if lo.any():
144+ q = np.sqrt(-2 * np.log(p[lo]))
145+ x[lo] = (((((c[0] * q + c[1]) * q + c[2]) * q + c[3]) * q + c[4]) * q + c[5]) / \
146+ ((((d[0] * q + d[1]) * q + d[2]) * q + d[3]) * q + 1)
147+ if hi.any():
148+ q = np.sqrt(-2 * np.log(1 - p[hi]))
149+ x[hi] = -(((((c[0] * q + c[1]) * q + c[2]) * q + c[3]) * q + c[4]) * q + c[5]) / \
150+ ((((d[0] * q + d[1]) * q + d[2]) * q + d[3]) * q + 1)
151+ return x
152+
153+
154+def trunc_std_normal(lo, hi, rng):
155+ """Standard normal truncated to [lo, hi], by inverse CDF. Mirrored into
156+ the lower tail so the CDF differences keep precision."""
157+ flip = (lo + hi) > 0
158+ a = np.where(flip, -hi, lo)
159+ b = np.where(flip, -lo, hi)
160+ Fa = norm_cdf(a)
161+ Fb = norm_cdf(b)
162+ u = Fa + (Fb - Fa) * rng.random(a.shape)
163+ x = norm_ppf(np.clip(u, 1e-300, 1 - 1e-16))
142164 x = np.where(flip, -x, x)
143- return mu + s * np.clip(x, alpha, beta)
165+ return np.clip(x, lo, hi)
144166
145167
146-def generate_z(h, sigma, dither, T, rng):
147- L = len(h)
148- x = rng.standard_normal(T + L - 1) * sigma
149- y = np.convolve(x, h, "valid")
168+# --------------------------------------------------------- the RB next-pmf
169+
170+def big_g(t):
171+ """G(t) = t Phi(t) + phi(t), the antiderivative of Phi."""
172+ return t * norm_cdf(t) + norm_pdf(t)
173+
174+
175+def next_pmf(c, s0, dither):
176+ """P(z_next = j | chain state): round(c + N(0, s0^2) (+ U(-1/2,1/2)))."""
177+ s = max(s0, 1e-12)
178+ half = 8 * s + (1.0 if dither else 0.0) + 1.0
179+ js = np.arange(math.floor(c - half), math.ceil(c + half) + 1)
150180 if dither:
151- y = y + rng.uniform(-0.5, 0.5, T)
152- return np.round(y).astype(np.int64)
153-
154-
155-def smc_replicate(h, sigma, dither, T, burn, N, kmax, rng):
156- """One replicate: -(1/(T-burn)) sum log2 phat(z_t | z_<t) after burn-in."""
157- A, v = levinson_all(autocovariance(h, sigma, kmax), kmax)
158- s_by_order = np.sqrt(v)
159- z = generate_z(h, sigma, dither, T, rng)
160-
161- W = np.zeros((N, kmax), dtype=np.float32) # each particle's last kmax values
162- log_p = np.zeros(T)
163- bad_steps = 0
164- for t in range(T):
165- p = min(t, kmax)
166- a = A[p]
167- s = s_by_order[p]
168- mu = (W[:, kmax - p :] @ a[::-1].astype(np.float32)).astype(np.float64) if p else np.zeros(N)
169- d = rng.uniform(-0.5, 0.5, N) if dither else 0.0
170- lo = z[t] - 0.5 - d
171- hi = z[t] + 0.5 - d
172- logw = log_phi_diff((lo - mu) / s, (hi - mu) / s)
173- lse = logsumexp(logw)
174- if not np.isfinite(lse):
175- raise RuntimeError(
176- f"particle collapse at step {t}: no particle is consistent with "
177- f"the observation — rerun with more --particles"
178- )
179- log_p[t] = lse - math.log(N)
180- # A per-step surprisal beyond ~40 bits means the particle cloud has
181- # drifted away from every path consistent with the data — the
182- # genealogical-collapse failure mode of docs/mc-true-rate.md §5.3,
183- # not a property of the data. Fail loudly rather than average it in.
184- if -log_p[t] / math.log(2) > 40:
185- bad_steps += 1
186- if bad_steps > 25:
187- raise RuntimeError(
188- "the sequential filter degenerated (near-deterministic dynamics; "
189- "see docs/mc-true-rate.md §5.3) — this kernel needs the "
190- "lookahead/twisted-proposal extension. Reduce --taps to study "
191- "the trend with a shallower stopband."
192- )
193- # Systematic resampling, then extend the chosen ancestors.
194- probs = np.exp(logw - logw.max())
195- cdf = np.cumsum(probs)
196- cdf /= cdf[-1]
197- ancestors = np.searchsorted(cdf, (rng.random() + np.arange(N)) / N)
198- mu_a = mu[ancestors]
199- lo_a = lo[ancestors] if dither else np.full(N, lo)
200- hi_a = hi[ancestors] if dither else np.full(N, hi)
201- y_new = sample_truncated(mu_a, s, lo_a, hi_a, rng)
202- W = W[ancestors]
203- W[:, :-1] = W[:, 1:]
204- W[:, -1] = y_new.astype(np.float32)
205- return float(-log_p[burn:].mean() / math.log(2))
206-
207-
208-def estimate_true_rate(h, sigma, dither, args, seed):
209- rates = []
210- for j in range(args.replicates):
211- rng = np.random.default_rng(seed + j)
212- r = smc_replicate(h, sigma, dither, args.steps, args.burn, args.particles, args.kmax, rng)
213- rates.append(r)
214- print(f" replicate {j + 1}/{args.replicates}: {r:.4f}", flush=True)
215- rates = np.array(rates)
216- se = rates.std(ddof=1) / math.sqrt(len(rates)) if len(rates) > 1 else float("nan")
217- return rates.mean(), se
218-
219-# ------------------------------------------------------------- self-test
220-
221-def selftest() -> int:
222- failures = 0
223-
224- def check(name, got, want, tol):
225- nonlocal failures
226- ok = abs(got - want) <= tol
227- failures += 0 if ok else 1
228- print(f" {'PASS' if ok else 'FAIL'} {name}: got {got:.5f}, want {want:.5f} (tol {tol})")
229-
230- print("H_delta against reference values:")
231- for s, want in [(0.1, 0.00001), (0.3, 0.55042), (1, 2.10483), (5, 4.37142), (50, 7.69098)]:
232- check(f"H_delta({s})", h_delta(s), want, 2e-4)
233-
234- print("Formula against reference values (iid):")
235- none = np.array([1.0])
236- for s, want in [(0.1, 0.57611), (1, 2.15829), (5, 4.37382), (100, 8.69096)]:
237- check(f"R_approx none sigma={s}", r_approx(none, s, False), want, 1e-3)
238-
239- print("Formula against reference values (default bandpass, end-to-end kernel check):")
240- lo = windowed_sinc_lowpass(300 / 30000, 101)
241- hi = windowed_sinc_lowpass(6000 / 30000, 101)
242- bp = hi - lo
243- check("||h||_2", float(np.sqrt((bp**2).sum())), 0.60216, 1e-4)
244- for s, want in [(0.5, 0.8620), (5, 1.9397), (20, 2.7282), (100, 3.7166)]:
245- check(f"R_approx bandpass sigma={s}", r_approx(bp, s, False), want, 1e-3)
246- check("R_approx bandpass sigma=5 dither", r_approx(bp, 5, True), 2.2111, 1e-3)
247-
248- print("Quick MC on the exact iid case (truth = H_delta(sigma)):")
249- rng = np.random.default_rng(1)
250- est = smc_replicate(none, 5.0, False, T=1500, burn=300, N=512, kmax=8, rng=rng)
251- check("MC none sigma=5", est, h_delta(5.0), 0.05)
252- est = smc_replicate(none, 0.5, False, T=1500, burn=300, N=512, kmax=8, rng=rng)
253- check("MC none sigma=0.5", est, h_delta(0.5), 0.05)
254-
255- print("FAILED" if failures else "all tests passed")
256- return 1 if failures else 0
257-
258-# ------------------------------------------------------------------ main
259-
260-def main() -> int:
261- ap = argparse.ArgumentParser(
262- description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter
263- )
264- ap.add_argument("--sigma", type=float, default=5.0, help="input std in quantization steps")
265- ap.add_argument(
266- "--filter",
267- choices=["none", "moving-average", "lowpass", "bandpass", "first-difference"],
268- default="none",
269- )
270- ap.add_argument("--low-hz", type=float, default=300, help="bandpass low edge")
271- ap.add_argument("--high-hz", type=float, default=6000, help="bandpass high edge")
272- ap.add_argument("--cutoff-hz", type=float, default=6000, help="lowpass cutoff")
273- ap.add_argument("--taps", type=int, default=101, help="windowed-sinc kernel length")
274- ap.add_argument("--width", type=int, default=8, help="moving-average width")
275- ap.add_argument("--sample-rate", type=float, default=30000)
276- ap.add_argument("--dither", action="store_true")
277- ap.add_argument("--particles", type=int, default=2048)
278- ap.add_argument("--steps", type=int, default=3000)
279- ap.add_argument("--burn", type=int, default=800)
280- ap.add_argument("--replicates", type=int, default=8)
281- ap.add_argument("--kmax", type=int, default=0, help="predictor memory; 0 = auto (4L, capped 512)")
282- ap.add_argument("--seed", type=int, default=0)
283- ap.add_argument("--selftest", action="store_true")
284- args = ap.parse_args()
181+ # Integrating the Gaussian bin probability over the dither gives a
182+ # second difference of G; as s -> 0 it degrades gracefully to the
183+ # uniform-overlap width.
184+ p = s * (big_g((js + 1 - c) / s) - 2 * big_g((js - c) / s) + big_g((js - 1 - c) / s))
185+ else:
186+ edges = norm_cdf((np.append(js, js[-1] + 1) - 0.5 - c) / s)
187+ p = np.diff(edges)
188+ return js, np.maximum(p, 0.0)
189+
285190
286- if args.selftest:
287- return selftest()
191+# ---------------------------------------------------------------- one past
288192
289- h = design_kernel(args)
193+def conditional_entropy_of_one_past(kernel, sigma, dither, M, sweeps, rng):
194+ h = np.asarray(kernel, dtype=float)
290195 L = len(h)
291- if args.kmax == 0:
292- args.kmax = min(max(4 * L, 64), 512)
293- args.burn = max(args.burn, 2 * args.kmax)
294- if args.steps <= args.burn + 500:
295- args.steps = args.burn + 2000
296-
297- norm = float(np.sqrt((h**2).sum()))
298- formula = r_approx(h, args.sigma, args.dither)
299- _, v = levinson_all(autocovariance(h, args.sigma, args.kmax), args.kmax)
300- s_inn = math.sqrt(v[args.kmax])
301- print(f"kernel: {args.filter}, L={L}, ||h||2={norm:.5f}, sigma_y={args.sigma * norm:.3f}, "
302- f"innovation std={s_inn:.4f}")
303- if s_inn < 0.12:
304- print("warning: innovation std << quantization step — near-deterministic dynamics. "
305- "The plain sequential filter (docs/mc-true-rate.md §5.3) will likely degenerate "
306- "here; reduce --taps for a shallower stopband, or implement the "
307- "lookahead/twisted-proposal extension.")
308- print(f"formula (as in the UI): R = {formula:.4f} bits/sample"
309- + (f", ideal ratio {16 / formula:.3f}x vs int16" if formula > 0 else ""))
310- if args.filter == "none":
311- exact = h_delta(args.sigma) if not args.dither else None
312- if exact is not None:
313- print(f"exact truth (iid closed form): R = {exact:.4f} bits/sample"
314- + (f", ideal ratio {16 / exact:.3f}x vs int16" if exact > 0 else ""))
315-
316- print(f"MC (N={args.particles} particles, T={args.steps} steps, burn={args.burn}, "
317- f"kmax={args.kmax}, {args.replicates} replicates):")
318- try:
319- mean, se = estimate_true_rate(h, args.sigma, args.dither, args, args.seed)
320- except RuntimeError as e:
321- print(f"aborted: {e}")
322- return 2
323- line = f"MC true rate: R = {mean:.4f} +/- {se:.4f} bits/sample"
324- if mean > 0:
325- # First-order error propagation: d(16/R) = 16 dR / R^2.
326- line += f", ideal ratio {16 / mean:.3f}x +/- {16 * se / mean**2:.3f} vs int16"
327- print(line)
328- print(f"difference (formula - MC): {formula - mean:+.4f} bits/sample")
329- print("note: finite window and inner MC both bias the estimate upward; "
330- "double --particles and compare to confirm convergence.")
331- return 0
332-
333-
334-if __name__ == "__main__":
335- sys.exit(main())
196+ hr = h[::-1]
197+ N = M + L - 1 # latents covering the windows of z_1..z_M
198+
199+ # The past, with its true latents as the (stationary) chain start.
200+ x = sigma * rng.standard_normal(N)
201+ y = np.convolve(x, h, mode='valid')
202+ d = (rng.random(M) - 0.5) if dither else None
203+ z = np.floor(y + (d if dither else 0.0) + 0.5)
204+
205+ # Boxes and y live in padded arrays so that every coordinate x_i sees
206+ # exactly L constraint rows (rows outside the data are unconstrained).
207+ P = L - 1
208+ ypad = np.zeros(M + 2 * P)
209+ lo = np.full(M + 2 * P, -np.inf)
210+ hi = np.full(M + 2 * P, np.inf)
211+
212+ def set_boxes():
213+ dd = d if dither else 0.0
214+ lo[P:P + M] = z - 0.5 - dd
215+ hi[P:P + M] = z + 0.5 - dd
216+
217+ set_boxes()
218+
219+ # Color classes: coordinates L apart share no constraint row, so a class
220+ # updates as one vectorized block. Row i+j (padded) carries coefficient
221+ # h[j] for coordinate i.
222+ classes = [np.arange(c0, N, L) for c0 in range(L)]
223+ rowmats = [idx[:, None] + np.arange(L)[None, :] for idx in classes]
224+ nonzero = h != 0
225+
226+ s0 = sigma * abs(h[0])
227+ pmf = {}
228+
229+ for _ in range(sweeps):
230+ ypad[P:P + M] = np.convolve(x, h, mode='valid') # kill fp drift
231+ if dither:
232+ ycur = ypad[P:P + M]
233+ dlo = np.maximum(-0.5, z - 0.5 - ycur)
234+ dhi = np.minimum(0.5, z + 0.5 - ycur)
235+ d = dlo + np.maximum(dhi - dlo, 0.0) * rng.random(M)
236+ set_boxes()
237+ for idx, rows in zip(classes, rowmats):
238+ r = ypad[rows] - np.outer(x[idx], h)
239+ with np.errstate(divide='ignore', invalid='ignore'):
240+ b1 = (lo[rows] - r) / h[None, :]
241+ b2 = (hi[rows] - r) / h[None, :]
242+ xlo = np.where(h[None, :] > 0, b1, b2)
243+ xhi = np.where(h[None, :] > 0, b2, b1)
244+ xlo[:, ~nonzero] = -np.inf
245+ xhi[:, ~nonzero] = np.inf
246+ xlo = xlo.max(axis=1)
247+ xhi = xhi.min(axis=1)
248+ xnew = trunc_std_normal(xlo / sigma, xhi / sigma, rng) * sigma
249+ delta = xnew - x[idx]
250+ x[idx] = xnew
251+ ypad[rows] += delta[:, None] * h[None, :]
252+
253+ c = float(hr[:-1] @ x[M:]) if L > 1 else 0.0
254+ js, p = next_pmf(c, s0, dither)
255+ for j, pj in zip(js, p):
256+ if pj > 1e-15:
257+ pmf[int(j)] = pmf.get(int(j), 0.0) + pj
258+
259+ total = sum(pmf.values())
260+ return -sum((p / total) * math.log2(p / total) for p in pmf.values() if p > 0)
261+
262+
263+# --------------------------------------------------------------------- main
264+
265+def main():
266+ ap = argparse.ArgumentParser(
267+ description='Monte-Carlo estimate of the true entropy rate, for '
268+ 'checking the analytic rate R shown in the app.')
269+ ap.add_argument('--sigma', type=float, required=True, help='input std, in quantization steps')
270+ ap.add_argument('--filter', required=True,
271+ choices=['none', 'moving-average', 'lowpass', 'bandpass', 'first-difference'])
272+ ap.add_argument('--low', type=float, help='bandpass low edge, Hz')
273+ ap.add_argument('--high', type=float, help='lowpass cutoff / bandpass high edge, Hz')
274+ ap.add_argument('--taps', type=int, default=101, help='windowed-sinc kernel length')
275+ ap.add_argument('--width', type=int, default=8, help='moving-average width')
276+ ap.add_argument('--rate', type=float, default=30000, help='sample rate, Hz')
277+ ap.add_argument('--dither', action='store_true')
278+ ap.add_argument('--past', type=int, help='conditioning window M (default max(512, 4·taps))')
279+ ap.add_argument('--pasts', type=int, default=24, help='independent pasts to average')
280+ ap.add_argument('--sweeps', type=int, default=600, help='Gibbs sweeps per past')
281+ ap.add_argument('--seed', type=int, default=0)
282+ args = ap.parse_args()
283+
284+ if args.filter == 'bandpass' and (args.low is None or args.high is None):
285+ ap.error('bandpass needs --low and --high')
286+ if args.filter == 'lowpass' and args.high is None:
287+ ap.error('lowpass needs --high')
288+
289+ kernel = design_kernel(args)
290+ L = len(kernel)
291+ M = args.past if args.past is not None else max(512, 4 * L)
292+
293+ R = formula_rate(kernel, args.sigma, args.dither)
294+ print(f'model: sigma={args.sigma} filter={args.filter} L={L} dither={args.dither}')
295+ print(f'formula R (as shown in the app): {R:.4f} bits/sample'
296+ f' (ratio vs int16: {16 / R:.3f}x)' if R > 0 else f'formula R: {R:.4f} bits/sample')
297+ print(f'MC: {args.pasts} pasts x {args.sweeps} sweeps, conditioning on M={M} samples')
298+
299+ rng = np.random.default_rng(args.seed)
300+ Hs = []
301+ for i in range(args.pasts):
302+ Hs.append(conditional_entropy_of_one_past(
303+ kernel, args.sigma, args.dither, M, args.sweeps, rng))
304+ mean = float(np.mean(Hs))
305+ se = float(np.std(Hs, ddof=1) / math.sqrt(len(Hs))) if len(Hs) > 1 else float('nan')
306+ print(f' past {i + 1:3d}/{args.pasts}: H = {Hs[-1]:.4f} running mean {mean:.4f} +/- {se:.4f}')
307+
308+ mean = float(np.mean(Hs))
309+ se = float(np.std(Hs, ddof=1) / math.sqrt(len(Hs)))
310+ print(f'\nMC entropy rate: {mean:.4f} +/- {se:.4f} bits/sample'
311+ f' (ratio vs int16: {16 / mean:.3f}x)')
312+ print(f'formula R: {R:.4f} bits/sample (formula - MC = {R - mean:+.4f})')
313+ print('note: the MC value estimates H(z_next | M past samples), an upper bound '
314+ 'on the rate that tightens as --past grows.')
315+
316+
317+if __name__ == '__main__':
318+ main()
src/App.tsxmodified+35−2View file
@@ -1,10 +1,10 @@
11 import { useEffect, useMemo, useRef, useState } from 'react'
2+import CopyableCommand from './components/CopyableCommand'
23 import Controls from './components/Controls'
34 import FilterViz from './components/FilterViz'
45 import ScrollingView from './components/ScrollingView'
56 import CompressionChart from './components/CompressionChart'
67 import MathSection from './components/MathSection'
7-import TrueRateCommand from './components/TrueRateCommand'
88 import { DEFAULT_SPEC, clampSpec, designKernel, kernelNorm } from './model/filters'
99 import { theoreticalRateBits } from './model/theory'
1010 import { LATENT_SEED } from './model/latent'
@@ -85,6 +85,36 @@ function useCompression(
8585 return state
8686 }
8787
88+/** The terminal command for scripts/true_rate.py at the current settings. */
89+function mcCommand(sigma: number, spec: ReturnType<typeof clampSpec>, rate: number, dither: boolean): string {
90+ const parts = ['python scripts/true_rate.py', `--sigma ${sigma}`]
91+ switch (spec.family) {
92+ case 'none':
93+ parts.push('--filter none')
94+ break
95+ case 'movingAverage':
96+ parts.push('--filter moving-average', `--width ${spec.width}`)
97+ break
98+ case 'lowpass':
99+ parts.push('--filter lowpass', `--high ${spec.highHz}`, `--taps ${spec.taps}`, `--rate ${rate}`)
100+ break
101+ case 'bandpass':
102+ parts.push(
103+ '--filter bandpass',
104+ `--low ${spec.lowHz}`,
105+ `--high ${spec.highHz}`,
106+ `--taps ${spec.taps}`,
107+ `--rate ${rate}`,
108+ )
109+ break
110+ case 'firstDifference':
111+ parts.push('--filter first-difference')
112+ break
113+ }
114+ if (dither) parts.push('--dither')
115+ return parts.join(' ')
116+}
117+
88118 export default function App() {
89119 const [sigma, setSigma] = useState(5)
90120 const [sampleRateHz, setSampleRateHz] = useState(30000)
@@ -203,7 +233,10 @@ export default function App() {
203233 formula in the math section — approximate where quantization dominates the spectrum
204234 (see the S(f) = 1 threshold on the response plot).
205235 </p>
206- <TrueRateCommand sigma={sigma} spec={spec} sampleRateHz={sampleRateHz} dither={dither} />
236+ <CopyableCommand
237+ label="check R against a Monte-Carlo ground truth (runs locally, ~2 min):"
238+ command={mcCommand(sigma, spec, sampleRateHz, dither)}
239+ />
207240 </section>
208241
209242 <section className="card">
src/app.cssmodified+13−19View file
@@ -463,46 +463,40 @@ body {
463463 font-weight: 500;
464464 }
465465
466-/* ---- MC cross-check command ---- */
467-
468-.cli-row {
466+.command-row {
469467 display: flex;
468+ flex-wrap: wrap;
470469 align-items: center;
471- gap: 8px;
470+ gap: 6px 10px;
472471 margin-top: 10px;
473- min-width: 0;
474-}
475-
476-.cli-label {
477472 font-size: 12px;
478473 color: var(--muted);
479- white-space: nowrap;
480474 }
481475
482-.cli-row code {
483- font-family: ui-monospace, 'SF Mono', Menlo, Consolas, monospace;
476+.command-row code {
477+ background: var(--page);
478+ border: 1px solid var(--border);
479+ border-radius: 6px;
480+ padding: 3px 8px;
484481 font-size: 11.5px;
485482 color: var(--ink-2);
486- background: color-mix(in srgb, var(--grid) 45%, transparent);
487- border-radius: 6px;
488- padding: 4px 8px;
489483 overflow-x: auto;
490484 white-space: nowrap;
485+ max-width: 100%;
491486 }
492487
493-.copy-btn {
488+.command-row button {
494489 background: var(--surface);
495490 color: var(--ink-2);
496491 border: 1px solid var(--baseline);
497492 border-radius: 6px;
498- padding: 3px 10px;
493+ padding: 2px 10px;
499494 font: inherit;
500- font-size: 12px;
495+ font-size: 11.5px;
501496 cursor: pointer;
502- white-space: nowrap;
503497 }
504498
505-.copy-btn:hover {
499+.command-row button:hover {
506500 color: var(--ink);
507501 }
508502
src/components/CopyableCommand.tsxadded+24−0View file
@@ -0,0 +1,24 @@
1+import { useEffect, useState } from 'react'
2+
3+/** A one-line terminal command with a copy button. */
4+export default function CopyableCommand({ label, command }: { label: string; command: string }) {
5+ const [copied, setCopied] = useState(false)
6+ useEffect(() => {
7+ if (!copied) return
8+ const t = setTimeout(() => setCopied(false), 1500)
9+ return () => clearTimeout(t)
10+ }, [copied])
11+ return (
12+ <div className="command-row">
13+ <span className="command-label">{label}</span>
14+ <code>{command}</code>
15+ <button
16+ onClick={() => {
17+ navigator.clipboard.writeText(command).then(() => setCopied(true))
18+ }}
19+ >
20+ {copied ? 'copied' : 'copy'}
21+ </button>
22+ </div>
23+ )
24+}
src/components/TrueRateCommand.tsxdeleted+0−74View file
@@ -1,74 +0,0 @@
1-import { useState } from 'react'
2-import type { FilterSpec } from '../model/filters'
3-
4-/** The scripts/true_rate.py invocation matching the current parameters. */
5-export function trueRateCommand(
6- sigma: number,
7- spec: FilterSpec,
8- sampleRateHz: number,
9- dither: boolean,
10-): string {
11- const parts = ['python scripts/true_rate.py', `--sigma ${sigma}`]
12- switch (spec.family) {
13- case 'none':
14- parts.push('--filter none')
15- break
16- case 'firstDifference':
17- parts.push('--filter first-difference')
18- break
19- case 'movingAverage':
20- parts.push('--filter moving-average', `--width ${spec.width}`)
21- break
22- case 'lowpass':
23- parts.push(
24- '--filter lowpass',
25- `--cutoff-hz ${spec.highHz}`,
26- `--taps ${spec.taps}`,
27- `--sample-rate ${sampleRateHz}`,
28- )
29- break
30- case 'bandpass':
31- parts.push(
32- '--filter bandpass',
33- `--low-hz ${spec.lowHz}`,
34- `--high-hz ${spec.highHz}`,
35- `--taps ${spec.taps}`,
36- `--sample-rate ${sampleRateHz}`,
37- )
38- break
39- }
40- if (dither) parts.push('--dither')
41- return parts.join(' ')
42-}
43-
44-/**
45- * The Monte Carlo cross-check, as a copy-pasteable command. The script
46- * estimates the true entropy rate (docs/mc-true-rate.md) so the dashed R can
47- * be compared against ground truth for the current parameters.
48- */
49-export default function TrueRateCommand(props: {
50- sigma: number
51- spec: FilterSpec
52- sampleRateHz: number
53- dither: boolean
54-}) {
55- const [copied, setCopied] = useState(false)
56- const command = trueRateCommand(props.sigma, props.spec, props.sampleRateHz, props.dither)
57- return (
58- <div className="cli-row">
59- <span className="cli-label">check R by Monte Carlo:</span>
60- <code>{command}</code>
61- <button
62- className="copy-btn"
63- onClick={() => {
64- navigator.clipboard.writeText(command).then(() => {
65- setCopied(true)
66- setTimeout(() => setCopied(false), 1500)
67- })
68- }}
69- >
70- {copied ? 'copied ✓' : 'copy'}
71- </button>
72- </div>
73- )
74-}