timeseries-entropy#

Unbiased Monte-Carlo estimation of the entropy of a quantized filtered Gaussian time series:

x iid N(0, sigma^2)  ->  y = h * x  ->  z = round(y)

The estimand is the conditional entropy H(z_{M+1} | z_1..z_M) in bits, which decreases to the entropy rate of z — the true lossless compression limit in bits/sample — as the past window M grows beyond the memory of the process.

Companion to timeseries-compressibility, which carries a hand-synced TypeScript port of this estimator in src/entropy/ (run in the browser from a worker) — when changing the algorithm here, change it there too.

Method#

Stationary conditional sampling. Draw x from the prior and push it through the pipeline to get a past z_1..z_M. The generating x is an exact draw from p(x | z), so a Gibbs chain started there is already in stationarity — no burn-in bias. Each Gibbs conditional is a box-truncated normal; after each sweep the free tail latent is drawn fresh, emitting one exact sample of z_{M+1}. The samples form a stationary, autocorrelated discrete chain. By default each chain is auto-thinned: a short probe measures the draws' integrated autocorrelation time tau and the chain then takes ceil(tau) sweeps per emitted sample (capped), so downstream levels see approximately independent draws.
Unbiased entropy of the chain's marginal. Plug-in entropies of blocks whose sizes double per level are combined by Rhee–Glynn randomized telescoping with antithetic half-block corrections Delta_m = h(B_m) - [h(B_m^1) + h(B_m^2)]/2, truncated at a random level N with P(N >= m) = 2^(-r m) and reweighted. The expectation is exactly H(z_{M+1} | that past) despite the plug-in bias at every finite block size and despite the autocorrelation, which affects only the variance.
Average over independent pasts to get H(z_{M+1} | z_1..z_M) with a valid standard error. The only remaining approximation to the entropy rate is the finite window M. Pasts are independent, so they run in parallel across processes (workers=, default all cores); each past has its own spawned RNG stream, making results deterministic per seed for any worker count.

Install#

pip install -e .

Requires numpy and scipy.

Usage#

from timeseries_entropy import estimate_conditional_entropy, kernels

est = estimate_conditional_entropy(kernels.moving_average(8), sigma=4.0)
print(est.mean, est.se)   # bits/sample, over independent pasts

Lower level: ConditionalChain(kernel, sigma, past, rng).draw(k) yields the stationary chain of z_{M+1} samples, and unbiased_entropy(draw, n0, r, rng) is one randomized-telescoping realization for any stationary discrete chain.

CLI#

timeseries-entropy --sigma 4 --filter moving-average --width 8
timeseries-entropy --sigma 8 --filter lowpass --high 3000 --rate 30000
timeseries-entropy --sigma 2 --filter none --pasts 8

Filters match the web app: none, moving-average, lowpass, bandpass, first-difference.

Tuning#

r (default 1.5) sets the truncation tail P(N >= m) = 2^(-r m). Finite expected work needs r > 1; finite variance needs E[Delta_m^2] to decay faster than 2^(-r m). Run --pilot 6 to see the RMS Delta_m decay before trusting a value.
thin (default 'auto') inserts extra Gibbs sweeps per emitted sample to cut autocorrelation. This is not optional for slowly mixing (narrowband, large-sigma) settings: without thinning their Delta_m decay too slowly for r = 1.5 and the estimator's variance is infinite — still unbiased, but a rare deep truncation level then returns a value of hundreds of bits (one such draw was traced producing 305 bits on the bandpass sigma=32 grid cell). 'auto' probes each chain's integrated autocorrelation time and thins by ceil(tau), capped at thin_cap (default 64); reps then acts as a per-past budget (realizations = max(1, reps / thin)), so per-past cost stays roughly flat and accuracy accumulates over pasts instead.
n0 (default 128) is the base block size.
--past sets M; increase it until the estimate stops moving to approach the rate.
--workers caps the process pool (default: all cores).

Cached estimates#

The estimates branch caches a grid of estimates so the common settings need not be recomputed: four filters (none, moving-average 8, lowpass 3000 Hz, bandpass 300-6000 Hz at 30 kHz) crossed with sigma in {1, 2, 4, 8, 16, 32}.

https://raw.githubusercontent.com/concept-collection/timeseries-entropy/estimates/estimates.json

The estimates workflow fills it every two hours, and can also be dispatched by hand. Each run draws 8 fresh, independent pasts per cell, appends one record per cell to runs.jsonl (including the per-past values and resolved thinning, so any outlier can be traced to its seed and replayed), and rebuilds estimates.json by pooling every past ever drawn — so the means keep tightening the more often it runs. The cache was reset on 2026-08-01: draws made before auto-thinning landed had heavy-tailed outliers on the narrowband large-sigma cells (see Tuning) and were discarded. The grid lives in scripts/grid.py; adding a cell does not invalidate the cache, since each record stores its own parameters. To fill the cache locally instead:

python scripts/run_sweep.py --data-dir data

Analytic prediction#

timeseries_entropy.theory predicts the rate from the Fourier modes H(f) = sum_j h_j e^(-2 pi i f j) of the kernel; the CLI prints it before each run. Szego's theorem gives the Gaussian one-step prediction error sigma_inf^2 = sigma^2 exp(int_0^1 ln|H(f)|^2 df) and the high-resolution rate (1/2) log2(2 pi e sigma_inf^2), which fails when H(f) has near-zero modes (the integral diverges negative while the true rate stays >= 0). Two quantization corrections fix it: the observed past is quantized, so the uniform roundoff power 1/12 is added to the spectrum before the geometric mean — which also keeps the integral finite in stopbands — and the next sample is quantized, so the Gaussian-uniform convolution entropy G(s) = h(N(0, s^2) + U(-1/2, 1/2)) replaces the Gaussian log term:

H_rate ~ G(s*),   s*^2 = exp( int_0^1 ln(sigma^2 |H(f)|^2 + 1/12) df ) - 1/12

See THEORY.md for the full derivation, its assumptions, and Monte-Carlo validation across filters and sigmas.