# 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](https://github.com/concept-collection/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

1. **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.
2. **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.
3. **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

```python
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](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](THEORY.md) for the full derivation, its assumptions, and
Monte-Carlo validation across filters and sigmas.