1# timeseries-entropy
3Unbiased Monte-Carlo estimation of the entropy of a quantized filtered
4Gaussian time series:
6 x iid N(0, sigma^2) -> y = h * x -> z = round(y)
8The estimand is the conditional entropy H(z_{M+1} | z_1..z_M) in bits, which
9decreases to the entropy rate of z — the true lossless compression limit in
10bits/sample — as the past window M grows beyond the memory of the process.
12Companion to [timeseries-compressibility](https://github.com/concept-collection/timeseries-compressibility),
13which carries a hand-synced TypeScript port of this estimator in `src/entropy/`
14(run in the browser from a worker) — when changing the algorithm here, change
15it there too.
17## Method
191. **Stationary conditional sampling.** Draw x from the prior and push it
20 through the pipeline to get a past z_1..z_M. The generating x is an exact
21 draw from p(x | z), so a Gibbs chain started there is already in
22 stationarity — no burn-in bias. Each Gibbs conditional is a box-truncated
23 normal; after each sweep the free tail latent is drawn fresh, emitting one
24 exact sample of z_{M+1}. The samples form a stationary, autocorrelated
25 discrete chain.
262. **Unbiased entropy of the chain's marginal.** Plug-in entropies of blocks
27 whose sizes double per level are combined by Rhee–Glynn randomized
28 telescoping with antithetic half-block corrections
29 Delta_m = h(B_m) - [h(B_m^1) + h(B_m^2)]/2, truncated at a random level N
30 with P(N >= m) = 2^(-r m) and reweighted. The expectation is exactly
31 H(z_{M+1} | that past) despite the plug-in bias at every finite block size
32 and despite the autocorrelation, which affects only the variance.
333. **Average over independent pasts** to get H(z_{M+1} | z_1..z_M) with a
34 valid standard error. The only remaining approximation to the entropy rate
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.
40## Install
42 pip install -e .
44Requires numpy and scipy.
46## Usage
48```python
49from timeseries_entropy import estimate_conditional_entropy, kernels
51est = estimate_conditional_entropy(kernels.moving_average(8), sigma=4.0)
52print(est.mean, est.se) # bits/sample, over independent pasts
53```
55Lower level: `ConditionalChain(kernel, sigma, past, rng).draw(k)` yields the
56stationary chain of z_{M+1} samples, and `unbiased_entropy(draw, n0, r, rng)`
57is one randomized-telescoping realization for any stationary discrete chain.
59## CLI
61 timeseries-entropy --sigma 4 --filter moving-average --width 8
62 timeseries-entropy --sigma 8 --filter lowpass --high 3000 --rate 30000
63 timeseries-entropy --sigma 2 --filter none --pasts 8
65Filters match the web app: `none`, `moving-average`, `lowpass`, `bandpass`,
66`first-difference`.
68## Tuning
70- `r` (default 1.5) sets the truncation tail P(N >= m) = 2^(-r m). Finite
71 expected work needs r > 1; finite variance needs E[Delta_m^2] to decay
72 faster than 2^(-r m). Run `--pilot 6` to see the RMS Delta_m decay before
73 trusting a value.
74- `n0` (default 128) is the base block size; `thin` inserts extra Gibbs
75 sweeps per emitted sample to cut autocorrelation for slowly mixing
76 (narrowband, large-sigma) settings.
77- `--past` sets M; increase it until the estimate stops moving to approach
78 the rate.
79- `--workers` caps the process pool (default: all cores).
81## Analytic prediction
83`timeseries_entropy.theory` predicts the rate from the Fourier modes
84H(f) = sum_j h_j e^(-2 pi i f j) of the kernel; the CLI prints it before
85each run. Szego's theorem gives the Gaussian one-step prediction error
86sigma_inf^2 = sigma^2 exp(int_0^1 ln|H(f)|^2 df) and the high-resolution
87rate (1/2) log2(2 pi e sigma_inf^2), which fails when H(f) has near-zero
88modes (the integral diverges negative while the true rate stays >= 0). Two
89quantization corrections fix it: the observed past is quantized, so the
90uniform roundoff power 1/12 is added to the spectrum before the geometric
91mean — which also keeps the integral finite in stopbands — and the next
92sample is quantized, so the Gaussian-uniform convolution entropy
93G(s) = h(N(0, s^2) + U(-1/2, 1/2)) replaces the Gaussian log term:
95 H_rate ~ G(s*), s*^2 = exp( int_0^1 ln(sigma^2 |H(f)|^2 + 1/12) df ) - 1/12
97See [THEORY.md](THEORY.md) for the full derivation, its assumptions, and
98Monte-Carlo validation across filters and sigmas.