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. By default each chain is auto-thinned: a short probe
26 measures the draws' integrated autocorrelation time tau and the chain
27 then takes ceil(tau) sweeps per emitted sample (capped), so downstream
28 levels see approximately independent draws.
292. **Unbiased entropy of the chain's marginal.** Plug-in entropies of blocks
30 whose sizes double per level are combined by Rhee–Glynn randomized
31 telescoping with antithetic half-block corrections
32 Delta_m = h(B_m) - [h(B_m^1) + h(B_m^2)]/2, truncated at a random level N
33 with P(N >= m) = 2^(-r m) and reweighted. The expectation is exactly
34 H(z_{M+1} | that past) despite the plug-in bias at every finite block size
35 and despite the autocorrelation, which affects only the variance.
363. **Average over independent pasts** to get H(z_{M+1} | z_1..z_M) with a
37 valid standard error. The only remaining approximation to the entropy rate
38 is the finite window M. Pasts are independent, so they run in parallel
39 across processes (`workers=`, default all cores); each past has its own
40 spawned RNG stream, making results deterministic per seed for any worker
41 count.
43## Install
45 pip install -e .
47Requires numpy and scipy.
49## Usage
51```python
52from timeseries_entropy import estimate_conditional_entropy, kernels
54est = estimate_conditional_entropy(kernels.moving_average(8), sigma=4.0)
55print(est.mean, est.se) # bits/sample, over independent pasts
56```
58Lower level: `ConditionalChain(kernel, sigma, past, rng).draw(k)` yields the
59stationary chain of z_{M+1} samples, and `unbiased_entropy(draw, n0, r, rng)`
60is one randomized-telescoping realization for any stationary discrete chain.
62## CLI
64 timeseries-entropy --sigma 4 --filter moving-average --width 8
65 timeseries-entropy --sigma 8 --filter lowpass --high 3000 --rate 30000
66 timeseries-entropy --sigma 2 --filter none --pasts 8
68Filters match the web app: `none`, `moving-average`, `lowpass`, `bandpass`,
69`first-difference`.
71## Tuning
73- `r` (default 1.5) sets the truncation tail P(N >= m) = 2^(-r m). Finite
74 expected work needs r > 1; finite variance needs E[Delta_m^2] to decay
75 faster than 2^(-r m). Run `--pilot 6` to see the RMS Delta_m decay before
76 trusting a value.
77- `thin` (default `'auto'`) inserts extra Gibbs sweeps per emitted sample to
78 cut autocorrelation. This is not optional for slowly mixing (narrowband,
79 large-sigma) settings: without thinning their Delta_m decay too slowly for
80 r = 1.5 and the estimator's variance is infinite — still unbiased, but a
81 rare deep truncation level then returns a value of hundreds of bits (one
82 such draw was traced producing 305 bits on the bandpass sigma=32 grid
83 cell). `'auto'` probes each chain's integrated autocorrelation time and
84 thins by ceil(tau), capped at `thin_cap` (default 64); `reps` then acts as
85 a per-past budget (realizations = max(1, reps / thin)), so per-past cost
86 stays roughly flat and accuracy accumulates over pasts instead.
87- `n0` (default 128) is the base block size.
88- `--past` sets M; increase it until the estimate stops moving to approach
89 the rate.
90- `--workers` caps the process pool (default: all cores).
92## Cached estimates
94The `estimates` branch caches a grid of estimates so the common settings need
95not be recomputed: four filters (`none`, moving-average 8, lowpass 3000 Hz,
96bandpass 300-6000 Hz at 30 kHz) crossed with sigma in {1, 2, 4, 8, 16, 32}.
98 https://raw.githubusercontent.com/concept-collection/timeseries-entropy/estimates/estimates.json
100The **estimates** workflow fills it every two hours, and can also be dispatched
101by hand. Each run draws 8 fresh, independent pasts per cell, appends one record
102per cell to `runs.jsonl` (including the per-past values and resolved
103thinning, so any outlier can be traced to its seed and replayed), and
104rebuilds `estimates.json` by pooling every past ever drawn — so the means keep
105tightening the more often it runs. The cache was reset on 2026-08-01: draws
106made before auto-thinning landed had heavy-tailed outliers on the
107narrowband large-sigma cells (see Tuning) and were discarded. The grid lives in
108[scripts/grid.py](scripts/grid.py); adding a cell does not invalidate the
109cache, since each record stores its own parameters. To fill the cache locally
110instead:
112 python scripts/run_sweep.py --data-dir data
114## Analytic prediction
116`timeseries_entropy.theory` predicts the rate from the Fourier modes
117H(f) = sum_j h_j e^(-2 pi i f j) of the kernel; the CLI prints it before
118each run. Szego's theorem gives the Gaussian one-step prediction error
119sigma_inf^2 = sigma^2 exp(int_0^1 ln|H(f)|^2 df) and the high-resolution
120rate (1/2) log2(2 pi e sigma_inf^2), which fails when H(f) has near-zero
121modes (the integral diverges negative while the true rate stays >= 0). Two
122quantization corrections fix it: the observed past is quantized, so the
123uniform roundoff power 1/12 is added to the spectrum before the geometric
124mean — which also keeps the integral finite in stopbands — and the next
125sample is quantized, so the Gaussian-uniform convolution entropy
126G(s) = h(N(0, s^2) + U(-1/2, 1/2)) replaces the Gaussian log term:
128 H_rate ~ G(s*), s*^2 = exp( int_0^1 ln(sigma^2 |H(f)|^2 + 1/12) df ) - 1/12
130See [THEORY.md](THEORY.md) for the full derivation, its assumptions, and
131Monte-Carlo validation across filters and sigmas.