1# Predicting the entropy rate of a quantized filtered Gaussian series
3The model is
5$$
6x_t \overset{\text{iid}}{\sim} \mathcal N(0,\sigma^2)
7\;\longrightarrow\;
8y_t = \sum_j h_j\, x_{t-j}
9\;\longrightarrow\;
10z_t = \mathrm{round}(y_t),
11$$
13with unit quantization step. The quantity of interest is the entropy rate
15$$
16\bar H \;=\; \lim_{M\to\infty} H\!\left(z_{M+1}\mid z_1,\dots,z_M\right)
17\qquad\text{[bits/sample]},
18$$
20the true lossless compression limit of $z$, which the Monte-Carlo estimator
21in this package approaches from above as the window $M$ grows. This note
22derives the analytic prediction implemented in
23[`theory.py`](src/timeseries_entropy/theory.py),
25$$
26\boxed{\;\bar H \;\approx\; G(s_*),
27\qquad
28s_*^2 \;=\; \exp\!\left(\int_0^1 \ln\!\big(\sigma^2\lvert H(f)\rvert^2 + \tfrac1{12}\big)\, df\right) \;-\; \frac1{12},\;}
29$$
31where $H(f) = \sum_j h_j e^{-2\pi i f j}$ are the Fourier modes of the
32kernel and $G(s) = h\big(\mathcal N(0,s^2) + \mathcal U(-\tfrac12,\tfrac12)\big)$
33is the differential entropy (in bits) of a Gaussian convolved with a unit
34uniform,
36$$
37G(s) \;=\; -\int_{-\infty}^{\infty} g_s(v)\,\log_2 g_s(v)\; dv,
38\qquad
39g_s(v) \;=\; \Phi\!\left(\frac{v + \tfrac12}{s}\right) - \Phi\!\left(\frac{v - \tfrac12}{s}\right),
40$$
42with $\Phi$ the standard normal CDF. The formula is built in three steps,
43each repairing a failure of the previous one.
45## Step 1 — Szegő–Kolmogorov: prediction from the exact past
47$y$ is stationary Gaussian with power spectral density
48$S(f) = \sigma^2 \lvert H(f)\rvert^2$, $f\in[0,1)$. Kolmogorov's form of
49Szegő's theorem says the one-step prediction error variance from the
50infinite (exact) past is the *geometric mean* of the spectrum:
52$$
53\sigma_\infty^2
54= \exp\!\left(\int_0^1 \ln S(f)\, df\right)
55= \sigma^2 \exp\!\left(\int_0^1 \ln \lvert H(f)\rvert^2\, df\right).
56$$
58Since a Gaussian process's entropy rate is the entropy of its innovation,
60$$
61\bar h(y) = \tfrac12\log_2\!\big(2\pi e\, \sigma_\infty^2\big),
62$$
64and in the high-resolution regime ($\sigma_\infty \gg$ 1 bin) the usual
65approximation $H(\mathrm{round}(Y)) \approx h(Y) - \log_2\Delta$ with
66$\Delta = 1$ gives the naive prediction
68$$
69\bar H \;\approx\; \tfrac12\log_2\!\big(2\pi e\,\sigma_\infty^2\big).
70$$
72**Closed form for FIR kernels.** Writing the tap polynomial
73$P(w) = \sum_j h_j w^j = c \prod_k (w - b_k)$, Jensen's formula gives the
74geometric mean of $\lvert w - b\rvert$ over the unit circle as
75$\max(1, \lvert b\rvert)$, so
77$$
78\sigma_\infty = \sigma\, \lvert c\rvert \prod_k \max\big(1, \lvert b_k\rvert\big).
79$$
81Examples: the moving average of width $W$ has all zeros on the unit circle
82and leading coefficient $1/W$, so $\sigma_\infty = \sigma/W$; the first
83difference $h = (1,-1)$ has $\int_0^1 \ln(4\sin^2\pi f)\,df = 0$, so
84$\sigma_\infty = \sigma$; the identity kernel has $\sigma_\infty = \sigma$.
86**Two failures.** (i) Wherever $\lvert H(f)\rvert \approx 0$ — the stopband
87of a lowpass or bandpass filter — the log integral dives toward $-\infty$
88and the formula predicts *negative* entropy ($-4.8$ bits for the
89$f_c = 0.1$ lowpass at $\sigma = 8$), while the truth is $\ge 0$. (ii) The
90predictor only sees the *quantized* past, which carries strictly less
91information than the exact past.
93## Step 2 — Quantized past: the $1/12$ noise floor
95Model roundoff as additive dither: $z_t = y_t + u_t$ with
96$u_t \overset{\text{iid}}{\sim} \mathcal U(-\tfrac12,\tfrac12)$,
97independent of $y$ (Bennett's approximation; exact under subtractive
98dither). The observed process $w = y + u$ then has spectrum
100$$
101S_w(f) = S(f) + \tfrac1{12}.
102$$
104Kolmogorov's theorem is a statement about *linear* prediction and needs no
105Gaussianity, so the one-step linear prediction error of $w$ from its past
106is $\exp \int_0^1 \ln S_w$. Because $u_{t+1}$ is independent of both
107$y_{t+1}$ and the past of $w$, its variance splits off exactly:
109$$
110\mathrm{Var}\big(w_{t+1}\mid w_{\le t}\big)
111= \mathrm{Var}\big(y_{t+1}\mid w_{\le t}\big) + \tfrac1{12}
112\quad\Longrightarrow\quad
113s_*^2 = \exp\!\left(\int_0^1 \ln\!\big(S(f) + \tfrac1{12}\big) df\right) - \frac1{12}.
114$$
116This $s_*$ is the effective uncertainty of the next sample given the
117quantized past. Three properties worth noting:
119- **Regularization.** The $\tfrac1{12}$ inside the logarithm is exactly
120 what keeps the integral finite at zeros of $H(f)$ — the fix for failure
121 (i) falls out of modeling failure (ii).
122- **Ordering.** $s_*^2 \ge \sigma_\infty^2$ always, because geometric means
123 are superadditive: quantizing the past can only increase the entropy
124 rate. The gap is a real, measurable effect even with no stopband (for the
125 first difference at $\sigma = 2$ it is $+0.10$ bits, confirmed by Monte
126 Carlo).
127- **Limits.** $s_*^2 \to \sigma_\infty^2$ when $S \gg \tfrac1{12}$
128 everywhere, and $s_* \to 0$ as $\sigma \to 0$ (for the identity kernel,
129 $s_* = \sigma$ *exactly*).
131## Step 3 — Quantized next sample: Gaussian ⊛ uniform entropy
133Given the quantized past, $y_{t+1} \approx \mathcal N(m, s_*^2)$ with a
134conditional mean $m$ that varies from past to past. When the marginal
135spread of $y$ covers many bins, $m \bmod 1$ equidistributes, so
137$$
138\bar H \;\approx\; \mathbb E_{c\sim\mathcal U(0,1)}\,
139H\!\big(\mathrm{round}(c + \mathcal N(0, s_*^2))\big).
140$$
142This average has a closed form — the standard dithered-quantization
143identity. For any $X$ with density, $\mathrm{round}(X + c) = k$ iff
144$X \in [k - c - \tfrac12,\, k - c + \tfrac12)$, an event of probability
145$g(k - c)$ where $g(v) = F_X(v + \tfrac12) - F_X(v - \tfrac12)$ is exactly
146the density of $X + U$, $U \sim \mathcal U(-\tfrac12, \tfrac12)$. The
147intervals $\{k - c : c \in (0,1)\}$ tile the line, so
149$$
150\mathbb E_c\, H\big(\mathrm{round}(X + c)\big)
151= -\int_0^1 \sum_k g(k - c) \log_2 g(k - c)\, dc
152= -\int_{-\infty}^{\infty} g \log_2 g
153= h(X + U).
154$$
156For Gaussian $X$ define
158$$
159G(s) = h\big(\mathcal N(0,s^2) + U\big)
160= -\int_{-\infty}^{\infty} g_s\log_2 g_s\,dv,
161\qquad
162g_s(v) = \Phi\!\left(\frac{v + \tfrac12}{s}\right) - \Phi\!\left(\frac{v - \tfrac12}{s}\right).
163$$
165Its limits are exactly the right ones:
167$$
168G(s) \;\to\; \tfrac12\log_2\!\big(2\pi e\,(s^2 + \tfrac1{12})\big)
169\quad (s \gg 1),
170\qquad
171G(s) \;\sim\; C\,s \;\to\; 0
172\quad (s \to 0),
173$$
175with $C = \int_{-\infty}^{\infty} h_2(\Phi(t))\,dt \approx 2.6061$
176($h_2$ the binary entropy). So $G$ reproduces the high-resolution formula
177when quantization is fine and saturates to $0$ — instead of diverging to
178$-\infty$ — when the conditional distribution concentrates inside one bin.
180## Validity and failure modes
182Monte-Carlo validation with this package's estimator (24+ independent
183pasts; `--thin 4` for the slowly mixing narrowband cases):
185| filter | $\sigma$ | $s_*$ | $G(s_*)$ | naive Szegő | Monte Carlo $\pm$ se |
186|---|---|---|---|---|---|
187| none | 0.5 | 0.50 | 1.2544 | 1.047 | 1.2380 ± 0.0069 |
188| none | 2 | 2.00 | 3.0620 | 3.047 | 3.0491 ± 0.0124 |
189| none | 8 | 8.00 | 5.0480 | 5.047 | 5.0009 ± 0.0372 |
190| first-diff | 0.5 | 0.60 | 1.4579 | 1.047 | 1.4579 ± 0.0053 |
191| first-diff | 2 | 2.13 | 3.1511 | 3.047 | 3.1709 ± 0.0192 |
192| first-diff | 8 | 8.14 | 5.0731 | 5.047 | 5.1004 ± 0.0489 |
193| MA(8) | 1 | 0.23 | 0.6100 | −0.953 | 0.4962 ± 0.0220 |
194| MA(8) | 2 | 0.38 | 0.9825 | 0.047 | 0.9881 ± 0.0092 |
195| MA(8) | 4 | 0.65 | 1.5529 | 1.047 | 1.5641 ± 0.0191 |
196| MA(8) | 32 | 4.18 | 4.1125 | 4.047 | 4.1420 ± 0.0429 |
197| lowpass $f_c$=0.1 | 8 | 0.50 | 1.2632 | −4.756 | 1.2232 ± 0.0190 (M=512), 1.2732 ± 0.0422 (M=1024) |
198| lowpass $f_c$=0.1 | 64 | 0.89 | 1.9513 | −1.756 | 1.8258 ± 0.0796 |
199| bandpass 0.01–0.2 | 8 | 1.03 | 2.1476 | −1.783 | 2.1362 ± 0.0217 |
201The approximations, and where they bite:
2031. **Dither independence** (Step 2) requires the marginal spread
204 $\sigma\lVert h\rVert_2$ to be at least about one bin. First-difference
205 at $\sigma = 0.5$ (spread 0.71 bins) still agrees to within its se;
206 MA(8) at $\sigma = 1$ (spread 0.35 bins) is overpredicted by
207 $\approx 0.11$ bits — when the whole signal lives inside one bin,
208 roundoff is deterministic, not dither-like, and the true rate is lower.
2092. **Equidistribution of the conditional mean** (Step 3) fails for kernels
210 with no memory: the identity kernel pins $m = 0$, and the exact answer
211 is the *centered* quantized-Gaussian entropy, below $G(\sigma)$ by
212 $\approx 0.013$ bits at $\sigma = 0.5$ (and exponentially little for
213 $\sigma \gtrsim 1$). Any kernel with real memory washes this out.
2143. **Linear prediction / Gaussianity of $w$** (Step 2): $w$ is not
215 Gaussian, and linear prediction of it is not optimal, so $s_*$ errs
216 slightly high; the effect is within the Monte-Carlo error bars above.
2174. **Near-singular spectra** (lowpass/bandpass) have long memory; the
218 Monte-Carlo column is an upper bound that decreases in $M$, and the
219 Gibbs sampler mixes slowly (hence `--thin`). The $M = 1024$ lowpass run
220 agrees with the prediction to well within its error bar.
222## Numerical notes
224- $\sigma_\infty$ is computed exactly from the roots of the tap polynomial
225 (robust to zeros of $H$ *on* the unit circle, where the log integral is
226 still convergent — an integrable singularity).
227- The Step-2 integral uses the trapezoid rule on the rfft grid over
228 $[0, \tfrac12]$ with $n = 2^{18}$ points; the $\tfrac1{12}$ floor makes
229 the integrand smooth and strictly positive, so no special handling of
230 spectral zeros is needed.
231- $G(s)$ is integrated on a grid of spacing $\min(s/8, 0.01)$; below
232 $s = 10^{-3}$ the linear asymptote $C s$ is used.