Analytic entropy-rate prediction; parallelize pasts across processes
THEORY.md derives H ~ G(s*) from Szego's theorem with a 1/12 quantization
floor on the spectrum, validated against the estimator across filters and
sigmas; theory.py implements it and the CLI prints it before each run.
Independent pasts (and the pilot) now run on a process pool with spawned
per-past RNG streams, deterministic per seed for any worker count.
6 changed files+524−18
README.mdmodified+24−1View file
@@ -32,7 +32,10 @@ it there too.
3232 and despite the autocorrelation, which affects only the variance.
3333 3. **Average over independent pasts** to get H(z_{M+1} | z_1..z_M) with a
3434 valid standard error. The only remaining approximation to the entropy rate
35- is the finite window M.
35+ is the finite window M. Pasts are independent, so they run in parallel
36+ across processes (`workers=`, default all cores); each past has its own
37+ spawned RNG stream, making results deterministic per seed for any worker
38+ count.
3639
3740 ## Install
3841
@@ -73,3 +76,23 @@ Filters match the web app: `none`, `moving-average`, `lowpass`, `bandpass`,
7376 (narrowband, large-sigma) settings.
7477 - `--past` sets M; increase it until the estimate stops moving to approach
7578 the rate.
79+- `--workers` caps the process pool (default: all cores).
80+
81+## Analytic prediction
82+
83+`timeseries_entropy.theory` predicts the rate from the Fourier modes
84+H(f) = sum_j h_j e^(-2 pi i f j) of the kernel; the CLI prints it before
85+each run. Szego's theorem gives the Gaussian one-step prediction error
86+sigma_inf^2 = sigma^2 exp(int_0^1 ln|H(f)|^2 df) and the high-resolution
87+rate (1/2) log2(2 pi e sigma_inf^2), which fails when H(f) has near-zero
88+modes (the integral diverges negative while the true rate stays >= 0). Two
89+quantization corrections fix it: the observed past is quantized, so the
90+uniform roundoff power 1/12 is added to the spectrum before the geometric
91+mean — which also keeps the integral finite in stopbands — and the next
92+sample is quantized, so the Gaussian-uniform convolution entropy
93+G(s) = h(N(0, s^2) + U(-1/2, 1/2)) replaces the Gaussian log term:
94+
95+ H_rate ~ G(s*), s*^2 = exp( int_0^1 ln(sigma^2 |H(f)|^2 + 1/12) df ) - 1/12
96+
97+See [THEORY.md](THEORY.md) for the full derivation, its assumptions, and
98+Monte-Carlo validation across filters and sigmas.
THEORY.mdadded+232−0View file
@@ -0,0 +1,232 @@
1+# Predicting the entropy rate of a quantized filtered Gaussian series
2+
3+The model is
4+
5+$$
6+x_t \overset{\text{iid}}{\sim} \mathcal N(0,\sigma^2)
7+\;\longrightarrow\;
8+y_t = \sum_j h_j\, x_{t-j}
9+\;\longrightarrow\;
10+z_t = \mathrm{round}(y_t),
11+$$
12+
13+with unit quantization step. The quantity of interest is the entropy rate
14+
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+$$
19+
20+the true lossless compression limit of $z$, which the Monte-Carlo estimator
21+in this package approaches from above as the window $M$ grows. This note
22+derives the analytic prediction implemented in
23+[`theory.py`](src/timeseries_entropy/theory.py),
24+
25+$$
26+\boxed{\;\bar H \;\approx\; G(s_*),
27+\qquad
28+s_*^2 \;=\; \exp\!\left(\int_0^1 \ln\!\big(\sigma^2\lvert H(f)\rvert^2 + \tfrac1{12}\big)\, df\right) \;-\; \frac1{12},\;}
29+$$
30+
31+where $H(f) = \sum_j h_j e^{-2\pi i f j}$ are the Fourier modes of the
32+kernel and $G(s) = h\big(\mathcal N(0,s^2) + \mathcal U(-\tfrac12,\tfrac12)\big)$
33+is the differential entropy (in bits) of a Gaussian convolved with a unit
34+uniform,
35+
36+$$
37+G(s) \;=\; -\int_{-\infty}^{\infty} g_s(v)\,\log_2 g_s(v)\; dv,
38+\qquad
39+g_s(v) \;=\; \Phi\!\left(\frac{v + \tfrac12}{s}\right) - \Phi\!\left(\frac{v - \tfrac12}{s}\right),
40+$$
41+
42+with $\Phi$ the standard normal CDF. The formula is built in three steps,
43+each repairing a failure of the previous one.
44+
45+## Step 1 — Szegő–Kolmogorov: prediction from the exact past
46+
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
49+Szegő's theorem says the one-step prediction error variance from the
50+infinite (exact) past is the *geometric mean* of the spectrum:
51+
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+$$
57+
58+Since a Gaussian process's entropy rate is the entropy of its innovation,
59+
60+$$
61+\bar h(y) = \tfrac12\log_2\!\big(2\pi e\, \sigma_\infty^2\big),
62+$$
63+
64+and in the high-resolution regime ($\sigma_\infty \gg$ 1 bin) the usual
65+approximation $H(\mathrm{round}(Y)) \approx h(Y) - \log_2\Delta$ with
66+$\Delta = 1$ gives the naive prediction
67+
68+$$
69+\bar H \;\approx\; \tfrac12\log_2\!\big(2\pi e\,\sigma_\infty^2\big).
70+$$
71+
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
74+geometric mean of $\lvert w - b\rvert$ over the unit circle as
75+$\max(1, \lvert b\rvert)$, so
76+
77+$$
78+\sigma_\infty = \sigma\, \lvert c\rvert \prod_k \max\big(1, \lvert b_k\rvert\big).
79+$$
80+
81+Examples: the moving average of width $W$ has all zeros on the unit circle
82+and leading coefficient $1/W$, so $\sigma_\infty = \sigma/W$; the first
83+difference $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$.
85+
86+**Two failures.** (i) Wherever $\lvert H(f)\rvert \approx 0$ — the stopband
87+of a lowpass or bandpass filter — the log integral dives toward $-\infty$
88+and 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
90+predictor only sees the *quantized* past, which carries strictly less
91+information than the exact past.
92+
93+## Step 2 — Quantized past: the $1/12$ noise floor
94+
95+Model roundoff as additive dither: $z_t = y_t + u_t$ with
96+$u_t \overset{\text{iid}}{\sim} \mathcal U(-\tfrac12,\tfrac12)$,
97+independent of $y$ (Bennett's approximation; exact under subtractive
98+dither). The observed process $w = y + u$ then has spectrum
99+
100+$$
101+S_w(f) = S(f) + \tfrac1{12}.
102+$$
103+
104+Kolmogorov's theorem is a statement about *linear* prediction and needs no
105+Gaussianity, so the one-step linear prediction error of $w$ from its past
106+is $\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:
108+
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
113+s_*^2 = \exp\!\left(\int_0^1 \ln\!\big(S(f) + \tfrac1{12}\big) df\right) - \frac1{12}.
114+$$
115+
116+This $s_*$ is the effective uncertainty of the next sample given the
117+quantized past. Three properties worth noting:
118+
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*).
130+
131+## Step 3 — Quantized next sample: Gaussian ⊛ uniform entropy
132+
133+Given the quantized past, $y_{t+1} \approx \mathcal N(m, s_*^2)$ with a
134+conditional mean $m$ that varies from past to past. When the marginal
135+spread of $y$ covers many bins, $m \bmod 1$ equidistributes, so
136+
137+$$
138+\bar H \;\approx\; \mathbb E_{c\sim\mathcal U(0,1)}\,
139+H\!\big(\mathrm{round}(c + \mathcal N(0, s_*^2))\big).
140+$$
141+
142+This average has a closed form — the standard dithered-quantization
143+identity. 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
146+the density of $X + U$, $U \sim \mathcal U(-\tfrac12, \tfrac12)$. The
147+intervals $\{k - c : c \in (0,1)\}$ tile the line, so
148+
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+$$
155+
156+For Gaussian $X$ define
157+
158+$$
159+G(s) = h\big(\mathcal N(0,s^2) + U\big)
160+= -\int_{-\infty}^{\infty} g_s\log_2 g_s\,dv,
161+\qquad
162+g_s(v) = \Phi\!\left(\frac{v + \tfrac12}{s}\right) - \Phi\!\left(\frac{v - \tfrac12}{s}\right).
163+$$
164+
165+Its limits are exactly the right ones:
166+
167+$$
168+G(s) \;\to\; \tfrac12\log_2\!\big(2\pi e\,(s^2 + \tfrac1{12})\big)
169+\quad (s \gg 1),
170+\qquad
171+G(s) \;\sim\; C\,s \;\to\; 0
172+\quad (s \to 0),
173+$$
174+
175+with $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
177+when quantization is fine and saturates to $0$ — instead of diverging to
178+$-\infty$ — when the conditional distribution concentrates inside one bin.
179+
180+## Validity and failure modes
181+
182+Monte-Carlo validation with this package's estimator (24+ independent
183+pasts; `--thin 4` for the slowly mixing narrowband cases):
184+
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 |
200+
201+The approximations, and where they bite:
202+
203+1. **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.
209+2. **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.
214+3. **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.
217+4. **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.
221+
222+## Numerical notes
223+
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.
src/timeseries_entropy/__init__.pymodified+38−10View file
@@ -2,6 +2,8 @@
22 time series: x iid N(0, sigma^2) -> y = h * x -> z = round(y)."""
33
44 import math
5+import os
6+from concurrent.futures import ProcessPoolExecutor, as_completed
57 from dataclasses import dataclass
68
79 import numpy as np
@@ -24,9 +26,16 @@ class Estimate:
2426 per_past: np.ndarray
2527
2628
29+def _one_past(kernel, sigma, M, thin, n0, r, reps, seed_seq):
30+ rng = np.random.default_rng(seed_seq)
31+ chain = ConditionalChain(kernel, sigma, M, rng, thin)
32+ return float(np.mean(
33+ [unbiased_entropy(chain.draw, n0, r, rng) for _ in range(reps)]))
34+
35+
2736 def estimate_conditional_entropy(kernel, sigma, past=None, pasts=24, reps=8,
2837 n0=128, r=1.5, thin=1, seed=None,
29- progress=None):
38+ progress=None, workers=None):
3039 """Unbiased estimate of H(z_{M+1} | z_1..z_M) in bits.
3140
3241 For each of `pasts` independent pasts, a stationary Gibbs chain of
@@ -37,19 +46,38 @@ def estimate_conditional_entropy(kernel, sigma, past=None, pasts=24, reps=8,
3746
3847 past defaults to max(512, 4 * len(kernel)); the estimand decreases toward
3948 the entropy rate as it grows. progress, if given, is called as
40- progress(i, values) after each past.
49+ progress(i, values) after each past finishes (completion order when
50+ parallel).
51+
52+ Pasts run in parallel across `workers` processes (default: all cores).
53+ Each past gets its own spawned RNG stream, so a given seed yields the
54+ same result for any worker count.
4155 """
4256 kernel = np.asarray(kernel, dtype=float)
4357 M = int(past) if past is not None else max(512, 4 * kernel.size)
44- rng = np.random.default_rng(seed)
58+ seeds = np.random.SeedSequence(seed).spawn(pasts)
59+ if workers is None:
60+ workers = min(pasts, os.cpu_count() or 1)
61+ per_past = np.empty(pasts)
4562 values = []
46- for i in range(pasts):
47- chain = ConditionalChain(kernel, sigma, M, rng, thin)
48- values.append(float(np.mean(
49- [unbiased_entropy(chain.draw, n0, r, rng) for _ in range(reps)])))
50- if progress is not None:
51- progress(i, values)
52- per_past = np.array(values)
63+ if workers <= 1:
64+ for i in range(pasts):
65+ per_past[i] = _one_past(kernel, sigma, M, thin, n0, r, reps,
66+ seeds[i])
67+ values.append(per_past[i])
68+ if progress is not None:
69+ progress(i, values)
70+ else:
71+ with ProcessPoolExecutor(max_workers=workers) as pool:
72+ futures = {
73+ pool.submit(_one_past, kernel, sigma, M, thin, n0, r, reps,
74+ seeds[i]): i
75+ for i in range(pasts)}
76+ for done, fut in enumerate(as_completed(futures)):
77+ per_past[futures[fut]] = fut.result()
78+ values.append(per_past[futures[fut]])
79+ if progress is not None:
80+ progress(done, values)
5381 se = (float(per_past.std(ddof=1) / math.sqrt(len(per_past)))
5482 if len(per_past) > 1 else float('nan'))
5583 return Estimate(float(per_past.mean()), se, per_past)
src/timeseries_entropy/cli.pymodified+24−7View file
@@ -1,11 +1,14 @@
11 """Command-line interface: timeseries-entropy [options]."""
22
33 import argparse
4+import os
5+from concurrent.futures import ProcessPoolExecutor
46
57 import numpy as np
68
79 from . import estimate_conditional_entropy, level_corrections, kernels
810 from .model import ConditionalChain
11+from .theory import predict_entropy_rate
912
1013
1114 def design_kernel(args):
@@ -56,6 +59,8 @@ def main():
5659 help='truncation exponent: P(N >= m) = 2^(-r m)')
5760 ap.add_argument('--thin', type=int, default=1,
5861 help='Gibbs sweeps per emitted sample')
62+ ap.add_argument('--workers', type=int,
63+ help='parallel processes over pasts (default: all cores)')
5964 ap.add_argument('--seed', type=int, default=0)
6065 ap.add_argument('--pilot', type=int, metavar='LEVELS',
6166 help='instead of estimating, print RMS Delta_m over the '
@@ -66,6 +71,11 @@ def main():
6671 L = len(kernel)
6772 M = args.past if args.past is not None else max(512, 4 * L)
6873 print(f'model: sigma={args.sigma} filter={args.filter} L={L} M={M}')
74+ pred = predict_entropy_rate(kernel, args.sigma)
75+ print(f'predicted rate: {pred["corrected"]:.4f} bits '
76+ f'(quantization-corrected; s*={pred["s_star"]:.4g}) '
77+ f'high-res Szego: {pred["highres"]:.4f} '
78+ f'(sigma_inf={pred["sigma_inf"]:.4g})')
6979
7080 if args.pilot is not None:
7181 run_pilot(kernel, args, M)
@@ -84,7 +94,7 @@ def main():
8494 est = estimate_conditional_entropy(
8595 kernel, args.sigma, past=M, pasts=args.pasts, reps=args.reps,
8696 n0=args.n0, r=args.r, thin=args.thin, seed=args.seed,
87- progress=progress)
97+ progress=progress, workers=args.workers)
8898 ratio = f' (ratio vs int16: {16 / est.mean:.3f}x)' if est.mean > 0 else ''
8999 print(f'\nH(z_next | {M} past samples) = {est.mean:.4f} +/- {est.se:.4f} '
90100 f'bits/sample{ratio}')
@@ -92,14 +102,21 @@ def main():
92102 'grows.')
93103
94104
105+def _one_pilot(kernel, sigma, M, thin, n0, levels, seed_seq):
106+ rng = np.random.default_rng(seed_seq)
107+ return level_corrections(
108+ ConditionalChain(kernel, sigma, M, rng, thin).draw, n0, levels)
109+
110+
95111 def run_pilot(kernel, args, M):
96- rng = np.random.default_rng(args.seed)
112+ seeds = np.random.SeedSequence(args.seed).spawn(args.pasts)
113+ workers = args.workers or min(args.pasts, os.cpu_count() or 1)
97114 print(f'pilot: {args.pasts} pasts, levels 1..{args.pilot}, n0={args.n0}')
98- deltas = np.array([
99- level_corrections(
100- ConditionalChain(kernel, args.sigma, M, rng, args.thin).draw,
101- args.n0, args.pilot)
102- for _ in range(args.pasts)])
115+ with ProcessPoolExecutor(max_workers=workers) as pool:
116+ deltas = np.array(list(pool.map(
117+ _one_pilot,
118+ *zip(*[(kernel, args.sigma, M, args.thin, args.n0, args.pilot, s)
119+ for s in seeds]))))
103120 rms = np.sqrt((deltas ** 2).mean(axis=0))
104121 for m in range(args.pilot):
105122 note = ''
src/timeseries_entropy/theory.pyadded+137−0View file
@@ -0,0 +1,137 @@
1+"""Analytic prediction of the entropy rate of z = round(h * x), x iid
2+N(0, sigma^2), from the Fourier modes H(f) = sum_j h_j e^{-2 pi i f j}.
3+
4+Derivation, in three steps:
5+
6+1. Szego / Kolmogorov. y = h * x is stationary Gaussian with power spectrum
7+ S(f) = sigma^2 |H(f)|^2. Its one-step prediction error variance from the
8+ exact past is the geometric mean of the spectrum,
9+
10+ sigma_inf^2 = exp( int_0^1 ln S(f) df ),
11+
12+ so in the high-resolution regime (sigma_inf >> 1 quantization step)
13+
14+ H_rate ~ 1/2 log2(2 pi e sigma_inf^2).
15+
16+ This fails in two ways tied to quantization: it ignores that only the
17+ *quantized* past is observed, and it diverges to -inf wherever H(f) has
18+ zeros or deep stopbands (the log integral blows up while the true entropy
19+ stays >= 0).
20+
21+2. Quantization noise floor. Model roundoff as additive dither: z_t ~ y_t +
22+ u_t with u_t iid U(-1/2, 1/2), variance 1/12, independent of y. The best
23+ linear one-step prediction of w = y + u has error variance
24+ exp(int ln(S + 1/12)) (Szego again, on S_w = S + 1/12), and u_{t+1} being
25+ unpredictable splits off exactly, leaving for y_{t+1} itself
26+
27+ s*^2 = exp( int_0^1 ln(S(f) + 1/12) df ) - 1/12.
28+
29+ The 1/12 inside the log is what keeps the integral finite at spectral
30+ zeros; s*^2 >= sigma_inf^2 always (geometric means are superadditive),
31+ and s*^2 -> 0 as sigma -> 0.
32+
33+3. Coarse-quantization saturation. Given the quantized past, y_{t+1} ~
34+ N(m, s*^2) with a conditional mean m that equidistributes modulo the
35+ quantization grid (valid when the marginal spread of y is >> 1 bin).
36+ The average discrete entropy of round(m + N(0, s^2)) over a uniform grid
37+ offset equals *exactly* the differential entropy of the Gaussian-uniform
38+ convolution (the standard dithered-quantization identity):
39+
40+ G(s) = h( N(0, s^2) + U(-1/2, 1/2) ) [bits],
41+
42+ which tends to 1/2 log2(2 pi e s^2) for s >> 1 and to 0 for s -> 0
43+ instead of diverging negative.
44+
45+Final formula:
46+
47+ H_rate ~ G( sqrt( exp( int_0^1 ln(sigma^2 |H(f)|^2 + 1/12) df ) - 1/12 ) )
48+
49+Known approximations: the dither term is really deterministic roundoff, not
50+independent noise; w is not Gaussian (linear prediction is not optimal); and
51+for kernels whose conditional mean does not equidistribute mod 1 (e.g. the
52+identity kernel, where it is constant) G slightly overestimates at small s.
53+"""
54+
55+import math
56+
57+import numpy as np
58+from scipy.special import ndtr
59+
60+TWO_PI_E = 2.0 * math.pi * math.e
61+
62+
63+def spectrum(kernel, n=1 << 18):
64+ """|H(f)|^2 on the rfft grid f = k/n, k = 0..n/2."""
65+ return np.abs(np.fft.rfft(np.asarray(kernel, dtype=float), n)) ** 2
66+
67+
68+def log_spectrum_mean(kernel, sigma, floor=1.0 / 12.0, n=1 << 18):
69+ """Mean over f in [0, 1) of ln(sigma^2 |H(f)|^2 + floor).
70+
71+ Trapezoid on [0, 1/2] (the spectrum is symmetric). floor > 0 keeps the
72+ integrand finite at zeros of H.
73+ """
74+ logS = np.log(sigma * sigma * spectrum(kernel, n) + floor)
75+ return float((logS.sum() - 0.5 * (logS[0] + logS[-1])) / (logS.size - 1))
76+
77+
78+def sigma_infinity(kernel, sigma):
79+ """Szego one-step prediction error std of y = h * x from its exact past:
80+ the root geometric mean of the spectrum, sigma * |c_lead| * prod max(1, |b_k|)
81+ over the zeros b_k of the tap polynomial (exact, robust to zeros of H on
82+ the unit circle, where the log integral is still convergent)."""
83+ c = np.trim_zeros(np.asarray(kernel, dtype=float)[::-1], 'f')
84+ b = np.roots(c) if c.size > 1 else np.array([])
85+ return sigma * abs(c[0]) * float(np.prod(np.maximum(1.0, np.abs(b))))
86+
87+
88+def gauss_uniform_entropy(s):
89+ """G(s) = h(N(0, s^2) + U(-1/2, 1/2)) in bits — the exact average entropy
90+ of round(c + N(0, s^2)) over a uniform grid offset c."""
91+ s = float(s)
92+ if s <= 0:
93+ return 0.0
94+ if s < 1e-3:
95+ return _edge_constant() * s
96+ dv = min(s / 8.0, 0.01)
97+ v = np.arange(0.0, 0.5 + 8.0 * s + 1.0, dv)
98+ g = ndtr((v + 0.5) / s) - ndtr((v - 0.5) / s)
99+ term = np.where(g > 0, -g * np.log2(np.clip(g, 1e-300, None)), 0.0)
100+ return float(2.0 * dv * (term.sum() - 0.5 * term[0]))
101+
102+
103+_EDGE_C = None
104+
105+
106+def _edge_constant():
107+ """int_-inf^inf h2(Phi(t)) dt: small-s slope of G(s)."""
108+ global _EDGE_C
109+ if _EDGE_C is None:
110+ t = np.linspace(-12.0, 12.0, 20001)
111+ p = np.clip(ndtr(t), 1e-300, 1 - 1e-16)
112+ h2 = -(p * np.log2(p) + (1 - p) * np.log2(1 - p))
113+ _EDGE_C = float(np.trapezoid(h2, t))
114+ return _EDGE_C
115+
116+
117+def predict_entropy_rate(kernel, sigma):
118+ """Both predictions of the entropy rate of z = round(h * x), in bits.
119+
120+ Returns a dict:
121+ sigma_inf Szego prediction error std from the exact past
122+ highres 1/2 log2(2 pi e sigma_inf^2) — valid when sigma_inf >~ 1
123+ s_star prediction error std of y_next from the quantized past
124+ corrected G(s_star) — stays valid at spectral zeros and coarse
125+ quantization
126+ """
127+ s_inf = sigma_infinity(kernel, sigma)
128+ highres = (0.5 * math.log2(TWO_PI_E * s_inf * s_inf)
129+ if s_inf > 0 else -math.inf)
130+ gm_w = math.exp(log_spectrum_mean(kernel, sigma))
131+ s_star = math.sqrt(max(gm_w - 1.0 / 12.0, 0.0))
132+ return {
133+ 'sigma_inf': s_inf,
134+ 'highres': highres,
135+ 's_star': s_star,
136+ 'corrected': gauss_uniform_entropy(s_star),
137+ }
tests/test_theory.pyadded+69−0View file
@@ -0,0 +1,69 @@
1+import math
2+
3+import numpy as np
4+import pytest
5+
6+from timeseries_entropy import kernels
7+from timeseries_entropy.theory import (
8+ gauss_uniform_entropy, log_spectrum_mean, predict_entropy_rate,
9+ sigma_infinity)
10+
11+
12+def test_sigma_infinity_closed_forms():
13+ # MA(W): all tap-polynomial zeros on the unit circle -> sigma / W.
14+ assert sigma_infinity(kernels.moving_average(8), 4.0) == pytest.approx(0.5)
15+ # first difference: 1 - w, zero at w = 1 -> sigma.
16+ assert sigma_infinity(kernels.first_difference(), 8.0) == pytest.approx(8.0)
17+ assert sigma_infinity(kernels.identity(), 2.0) == pytest.approx(2.0)
18+
19+
20+def test_log_spectrum_mean_matches_szego_when_floor_negligible():
21+ # With a floor far below the spectrum minimum, the FFT integral must
22+ # agree with the exact roots-based geometric mean.
23+ h = np.array([1.0, -0.5]) # min |H|^2 = 0.25, floor 1e-12 negligible
24+ gm = math.exp(log_spectrum_mean(h, 1.0, floor=1e-12))
25+ assert gm == pytest.approx(sigma_infinity(h, 1.0) ** 2, rel=1e-6)
26+
27+
28+def test_gauss_uniform_entropy_limits():
29+ # Large s: h(N + U) -> 1/2 log2(2 pi e (s^2 + 1/12)) (Gaussian limit).
30+ for s in [4.0, 16.0]:
31+ expect = 0.5 * math.log2(2 * math.pi * math.e * (s * s + 1 / 12))
32+ assert gauss_uniform_entropy(s) == pytest.approx(expect, abs=2e-3)
33+ # Small s: -> 0 linearly, never negative.
34+ assert gauss_uniform_entropy(0.0) == 0.0
35+ assert 0 < gauss_uniform_entropy(0.01) < 0.05
36+
37+
38+def test_gauss_uniform_entropy_is_offset_averaged_quantized_entropy():
39+ # G(s) = E_c H(round(c + N(0, s^2))), c ~ U(0, 1) (dither identity).
40+ s = 0.7
41+ cs = (np.arange(200) + 0.5) / 200
42+ zs = np.arange(-12, 13)
43+ edges = (zs[None, :] - 0.5 - cs[:, None]) / s
44+ from scipy.special import ndtr
45+ p = np.diff(ndtr(np.concatenate([edges, edges[:, -1:] + 1 / s], axis=1)))
46+ p = np.clip(p, 1e-300, None)
47+ mean_h = float((-p * np.log2(p)).sum(axis=1).mean())
48+ assert gauss_uniform_entropy(s) == pytest.approx(mean_h, abs=1e-3)
49+
50+
51+def test_predict_identity_kernel_high_resolution():
52+ # For h = [1], z is iid round(N(0, sigma^2)); at sigma >> 1 both
53+ # predictions must approach the exact marginal entropy.
54+ sigma = 16.0
55+ zmax = int(8 * sigma + 4)
56+ z = np.arange(-zmax, zmax + 1)
57+ from scipy.special import ndtr
58+ p = ndtr((z + 0.5) / sigma) - ndtr((z - 0.5) / sigma)
59+ p = p[p > 0]
60+ exact = float(-(p * np.log2(p)).sum())
61+ pred = predict_entropy_rate(kernels.identity(), sigma)
62+ assert pred['corrected'] == pytest.approx(exact, abs=2e-3)
63+ assert pred['highres'] == pytest.approx(exact, abs=2e-3)
64+
65+
66+def test_corrected_prediction_stays_finite_at_spectral_zeros():
67+ pred = predict_entropy_rate(kernels.windowed_sinc_lowpass(0.1, 101), 8.0)
68+ assert pred['highres'] < 0 # Szego diverges toward -inf
69+ assert 0 < pred['corrected'] < 16 # corrected saturates sensibly