Unbiased Monte-Carlo entropy estimation for quantized filtered Gaussian series
Gibbs-sample z_next given a fixed past (stationary start from the
generating latents), then apply Rhee-Glynn randomized telescoping with
antithetic half-block corrections to the sampled chain. Library, CLI
with the app's filter presets, and tests against exact/known entropies.
10 changed files+617−0
.gitignoreadded+6−0View file
@@ -0,0 +1,6 @@
1+__pycache__/
2+*.egg-info/
3+dist/
4+build/
5+.venv/
6+.pytest_cache/
README.mdadded+72−0View file
@@ -0,0 +1,72 @@
1+# timeseries-entropy
2+
3+Unbiased Monte-Carlo estimation of the entropy of a quantized filtered
4+Gaussian time series:
5+
6+ x iid N(0, sigma^2) -> y = h * x -> z = round(y)
7+
8+The estimand is the conditional entropy H(z_{M+1} | z_1..z_M) in bits, which
9+decreases to the entropy rate of z — the true lossless compression limit in
10+bits/sample — as the past window M grows beyond the memory of the process.
11+
12+Companion to [timeseries-compressibility](https://github.com/concept-collection/timeseries-compressibility).
13+
14+## Method
15+
16+1. **Stationary conditional sampling.** Draw x from the prior and push it
17+ through the pipeline to get a past z_1..z_M. The generating x is an exact
18+ draw from p(x | z), so a Gibbs chain started there is already in
19+ stationarity — no burn-in bias. Each Gibbs conditional is a box-truncated
20+ normal; after each sweep the free tail latent is drawn fresh, emitting one
21+ exact sample of z_{M+1}. The samples form a stationary, autocorrelated
22+ discrete chain.
23+2. **Unbiased entropy of the chain's marginal.** Plug-in entropies of blocks
24+ whose sizes double per level are combined by Rhee–Glynn randomized
25+ telescoping with antithetic half-block corrections
26+ Delta_m = h(B_m) - [h(B_m^1) + h(B_m^2)]/2, truncated at a random level N
27+ with P(N >= m) = 2^(-r m) and reweighted. The expectation is exactly
28+ H(z_{M+1} | that past) despite the plug-in bias at every finite block size
29+ and despite the autocorrelation, which affects only the variance.
30+3. **Average over independent pasts** to get H(z_{M+1} | z_1..z_M) with a
31+ valid standard error. The only remaining approximation to the entropy rate
32+ is the finite window M.
33+
34+## Install
35+
36+ pip install -e .
37+
38+Requires numpy and scipy.
39+
40+## Usage
41+
42+```python
43+from timeseries_entropy import estimate_conditional_entropy, kernels
44+
45+est = estimate_conditional_entropy(kernels.moving_average(8), sigma=4.0)
46+print(est.mean, est.se) # bits/sample, over independent pasts
47+```
48+
49+Lower level: `ConditionalChain(kernel, sigma, past, rng).draw(k)` yields the
50+stationary chain of z_{M+1} samples, and `unbiased_entropy(draw, n0, r, rng)`
51+is one randomized-telescoping realization for any stationary discrete chain.
52+
53+## CLI
54+
55+ timeseries-entropy --sigma 4 --filter moving-average --width 8
56+ timeseries-entropy --sigma 8 --filter lowpass --high 3000 --rate 30000
57+ timeseries-entropy --sigma 2 --filter none --pasts 8
58+
59+Filters match the web app: `none`, `moving-average`, `lowpass`, `bandpass`,
60+`first-difference`.
61+
62+## Tuning
63+
64+- `r` (default 1.5) sets the truncation tail P(N >= m) = 2^(-r m). Finite
65+ expected work needs r > 1; finite variance needs E[Delta_m^2] to decay
66+ faster than 2^(-r m). Run `--pilot 6` to see the RMS Delta_m decay before
67+ trusting a value.
68+- `n0` (default 128) is the base block size; `thin` inserts extra Gibbs
69+ sweeps per emitted sample to cut autocorrelation for slowly mixing
70+ (narrowband, large-sigma) settings.
71+- `--past` sets M; increase it until the estimate stops moving to approach
72+ the rate.
pyproject.tomladded+20−0View file
@@ -0,0 +1,20 @@
1+[build-system]
2+requires = ["hatchling"]
3+build-backend = "hatchling.build"
4+
5+[project]
6+name = "timeseries-entropy"
7+version = "0.1.0"
8+description = "Unbiased Monte-Carlo entropy estimation for quantized filtered Gaussian time series"
9+readme = "README.md"
10+requires-python = ">=3.9"
11+dependencies = ["numpy>=1.22", "scipy>=1.8"]
12+
13+[project.scripts]
14+timeseries-entropy = "timeseries_entropy.cli:main"
15+
16+[project.optional-dependencies]
17+dev = ["pytest"]
18+
19+[tool.hatch.build.targets.wheel]
20+packages = ["src/timeseries_entropy"]
src/timeseries_entropy/__init__.pyadded+55−0View file
@@ -0,0 +1,55 @@
1+"""Unbiased Monte-Carlo entropy estimation for quantized filtered Gaussian
2+time series: x iid N(0, sigma^2) -> y = h * x -> z = round(y)."""
3+
4+import math
5+from dataclasses import dataclass
6+
7+import numpy as np
8+
9+from .estimator import plugin_entropy, unbiased_entropy, level_corrections
10+from .model import ConditionalChain
11+from . import kernels
12+
13+__all__ = [
14+ 'estimate_conditional_entropy', 'Estimate', 'ConditionalChain',
15+ 'plugin_entropy', 'unbiased_entropy', 'level_corrections', 'kernels',
16+]
17+
18+
19+@dataclass
20+class Estimate:
21+ """mean +/- se (over independent pasts) of H(z_{M+1} | z_1..z_M), bits."""
22+ mean: float
23+ se: float
24+ per_past: np.ndarray
25+
26+
27+def estimate_conditional_entropy(kernel, sigma, past=None, pasts=24, reps=8,
28+ n0=128, r=1.5, thin=1, seed=None,
29+ progress=None):
30+ """Unbiased estimate of H(z_{M+1} | z_1..z_M) in bits.
31+
32+ For each of `pasts` independent pasts, a stationary Gibbs chain of
33+ z_{M+1} draws feeds `reps` randomized-telescoping realizations (on
34+ consecutive segments of the chain); their average is one unbiased value
35+ per past. Returns the mean and standard error over pasts — valid because
36+ pasts are independent.
37+
38+ past defaults to max(512, 4 * len(kernel)); the estimand decreases toward
39+ the entropy rate as it grows. progress, if given, is called as
40+ progress(i, values) after each past.
41+ """
42+ kernel = np.asarray(kernel, dtype=float)
43+ M = int(past) if past is not None else max(512, 4 * kernel.size)
44+ rng = np.random.default_rng(seed)
45+ 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)
53+ se = (float(per_past.std(ddof=1) / math.sqrt(len(per_past)))
54+ if len(per_past) > 1 else float('nan'))
55+ return Estimate(float(per_past.mean()), se, per_past)
src/timeseries_entropy/cli.pyadded+114−0View file
@@ -0,0 +1,114 @@
1+"""Command-line interface: timeseries-entropy [options]."""
2+
3+import argparse
4+
5+import numpy as np
6+
7+from . import estimate_conditional_entropy, level_corrections, kernels
8+from .model import ConditionalChain
9+
10+
11+def design_kernel(args):
12+ if args.filter == 'none':
13+ return kernels.identity()
14+ if args.filter == 'moving-average':
15+ return kernels.moving_average(args.width)
16+ if args.filter == 'lowpass':
17+ if args.high is None:
18+ raise SystemExit('lowpass needs --high')
19+ return kernels.windowed_sinc_lowpass(args.high / args.rate, args.taps)
20+ if args.filter == 'bandpass':
21+ if args.low is None or args.high is None:
22+ raise SystemExit('bandpass needs --low and --high')
23+ return kernels.windowed_sinc_bandpass(
24+ args.low / args.rate, args.high / args.rate, args.taps)
25+ if args.filter == 'first-difference':
26+ return kernels.first_difference()
27+ raise ValueError(args.filter)
28+
29+
30+def main():
31+ ap = argparse.ArgumentParser(
32+ prog='timeseries-entropy',
33+ description='Unbiased Monte-Carlo estimate of H(z_next | M past '
34+ 'samples), in bits, for x iid N(0, sigma^2) -> h * x '
35+ '-> round.')
36+ ap.add_argument('--sigma', type=float, required=True,
37+ help='input std, in quantization steps')
38+ ap.add_argument('--filter', required=True,
39+ choices=['none', 'moving-average', 'lowpass', 'bandpass',
40+ 'first-difference'])
41+ ap.add_argument('--low', type=float, help='bandpass low edge, Hz')
42+ ap.add_argument('--high', type=float,
43+ help='lowpass cutoff / bandpass high edge, Hz')
44+ ap.add_argument('--taps', type=int, default=101,
45+ help='windowed-sinc kernel length')
46+ ap.add_argument('--width', type=int, default=8, help='moving-average width')
47+ ap.add_argument('--rate', type=float, default=30000, help='sample rate, Hz')
48+ ap.add_argument('--past', type=int,
49+ help='conditioning window M (default max(512, 4*L))')
50+ ap.add_argument('--pasts', type=int, default=24,
51+ help='independent pasts to average')
52+ ap.add_argument('--reps', type=int, default=8,
53+ help='randomized realizations per past')
54+ ap.add_argument('--n0', type=int, default=128, help='base block size')
55+ ap.add_argument('--r', type=float, default=1.5,
56+ help='truncation exponent: P(N >= m) = 2^(-r m)')
57+ ap.add_argument('--thin', type=int, default=1,
58+ help='Gibbs sweeps per emitted sample')
59+ ap.add_argument('--seed', type=int, default=0)
60+ ap.add_argument('--pilot', type=int, metavar='LEVELS',
61+ help='instead of estimating, print RMS Delta_m over the '
62+ 'pasts for m = 1..LEVELS, to help choose --r')
63+ args = ap.parse_args()
64+
65+ kernel = design_kernel(args)
66+ L = len(kernel)
67+ M = args.past if args.past is not None else max(512, 4 * L)
68+ print(f'model: sigma={args.sigma} filter={args.filter} L={L} M={M}')
69+
70+ if args.pilot is not None:
71+ run_pilot(kernel, args, M)
72+ return
73+
74+ print(f'{args.pasts} pasts x {args.reps} reps, n0={args.n0} r={args.r} '
75+ f'thin={args.thin}')
76+
77+ def progress(i, values):
78+ mean = float(np.mean(values))
79+ se = (float(np.std(values, ddof=1) / np.sqrt(len(values)))
80+ if len(values) > 1 else float('nan'))
81+ print(f' past {i + 1:3d}/{args.pasts}: H = {values[-1]:.4f} '
82+ f'running mean {mean:.4f} +/- {se:.4f}')
83+
84+ est = estimate_conditional_entropy(
85+ kernel, args.sigma, past=M, pasts=args.pasts, reps=args.reps,
86+ n0=args.n0, r=args.r, thin=args.thin, seed=args.seed,
87+ progress=progress)
88+ ratio = f' (ratio vs int16: {16 / est.mean:.3f}x)' if est.mean > 0 else ''
89+ print(f'\nH(z_next | {M} past samples) = {est.mean:.4f} +/- {est.se:.4f} '
90+ f'bits/sample{ratio}')
91+ print('note: an upper bound on the entropy rate that tightens as --past '
92+ 'grows.')
93+
94+
95+def run_pilot(kernel, args, M):
96+ rng = np.random.default_rng(args.seed)
97+ 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)])
103+ rms = np.sqrt((deltas ** 2).mean(axis=0))
104+ for m in range(args.pilot):
105+ note = ''
106+ if m > 0 and rms[m] > 0:
107+ note = f' decay exponent {np.log2(rms[m - 1] / rms[m]) * 2:.2f}'
108+ print(f' m={m + 1}: rms Delta = {rms[m]:.5f}{note}')
109+ print('choose r safely below the E[Delta^2] decay exponent (and > 1); '
110+ 'r=1.5 suits decay near 2.')
111+
112+
113+if __name__ == '__main__':
114+ main()
src/timeseries_entropy/estimator.pyadded+96−0View file
@@ -0,0 +1,96 @@
1+"""Rhee-Glynn (randomized telescoping) unbiased entropy estimation.
2+
3+Input: a stationary, ergodic sequence of discrete draws, possibly
4+autocorrelated, each with the target marginal p. Plug-in entropies of blocks
5+whose sizes double from level to level form a telescoping sum via the
6+antithetic correction
7+
8+ Delta_m = h(B_m) - [h(B_m^1) + h(B_m^2)] / 2,
9+
10+where B_m^1, B_m^2 are the two halves of B_m. Truncating the sum at a random
11+level N with P(N >= m) = 2^(-r m) and reweighting by the survival
12+probabilities gives an estimator whose expectation is exactly the entropy of
13+p, despite the bias of every finite-block plug-in estimate.
14+
15+All entropies are in bits.
16+"""
17+
18+import numpy as np
19+
20+
21+def plugin_entropy(samples):
22+ """Plug-in entropy (bits) of a block of discrete samples."""
23+ _, counts = np.unique(np.asarray(samples), return_counts=True)
24+ return _entropy(counts.astype(float), counts.sum())
25+
26+
27+def _entropy(counts, n):
28+ return float(np.log2(n) - (counts * np.log2(counts)).sum() / n)
29+
30+
31+def _dict_entropy(counts, n):
32+ return _entropy(np.fromiter(counts.values(), dtype=float), n)
33+
34+
35+def _telescope(draw, n0, levels):
36+ """h(B_0) and [Delta_1, ..., Delta_levels] over one growing block.
37+
38+ Counts are merged upward (each sample is counted once), so the cost is
39+ linear in the n0 * 2**levels samples drawn.
40+ """
41+ counts = {}
42+ _count_into(counts, draw(n0))
43+ size = n0
44+ h0 = _dict_entropy(counts, size)
45+ h_prev = h0
46+ deltas = []
47+ for _ in range(levels):
48+ half = {}
49+ _count_into(half, draw(size))
50+ h2 = _dict_entropy(half, size)
51+ for v, c in half.items():
52+ counts[v] = counts.get(v, 0) + c
53+ size *= 2
54+ h_full = _dict_entropy(counts, size)
55+ deltas.append(h_full - 0.5 * (h_prev + h2))
56+ h_prev = h_full
57+ return h0, deltas
58+
59+
60+def _count_into(counts, seg):
61+ vals, cnts = np.unique(np.asarray(seg), return_counts=True)
62+ for v, c in zip(vals.tolist(), cnts.tolist()):
63+ counts[v] = counts.get(v, 0) + c
64+
65+
66+def unbiased_entropy(draw, n0=128, r=1.5, rng=None):
67+ """One randomized-telescoping realization of the marginal entropy (bits).
68+
69+ draw(k) must return the next k consecutive samples of a stationary
70+ discrete chain; successive calls continue the chain. The realization
71+ consumes n0 * 2**N samples with P(N >= m) = 2^(-r m). Average many
72+ realizations (they may continue one chain back-to-back) to reduce
73+ variance; each has expectation exactly H(p).
74+
75+ r trades expected work against variance: levels must decay like
76+ E[Delta_m^2] = O(2^(-r'm)) with r' > r > 1 for both to be finite.
77+ r = 1.5 suits the typical second-moment decay 2^(-2m); check with
78+ level_corrections when in doubt.
79+ """
80+ rng = np.random.default_rng() if rng is None else rng
81+ rho = 2.0 ** -r
82+ N = 0
83+ while rng.random() < rho:
84+ N += 1
85+ h0, deltas = _telescope(draw, n0, N)
86+ return h0 + sum(d * 2.0 ** (r * m) for m, d in enumerate(deltas, start=1))
87+
88+
89+def level_corrections(draw, n0=128, levels=6):
90+ """Deterministic pilot: [Delta_1, ..., Delta_levels] from one block.
91+
92+ Draws n0 * 2**levels samples. Repeat over fresh chains, look at the decay
93+ of mean(Delta_m^2) with m, and pick r below the decay exponent.
94+ """
95+ _, deltas = _telescope(draw, n0, levels)
96+ return np.array(deltas)
src/timeseries_entropy/kernels.pyadded+33−0View file
@@ -0,0 +1,33 @@
1+"""Kernel constructors matching timeseries-compressibility's filters."""
2+
3+import numpy as np
4+
5+
6+def identity():
7+ return np.array([1.0])
8+
9+
10+def moving_average(width):
11+ return np.full(width, 1.0 / width)
12+
13+
14+def first_difference():
15+ return np.array([1.0, -1.0])
16+
17+
18+def windowed_sinc_lowpass(fc, taps):
19+ """Hamming-windowed sinc, cutoff fc in cycles/sample, unit DC gain."""
20+ n = taps | 1
21+ mid = (n - 1) / 2
22+ i = np.arange(n)
23+ t = i - mid
24+ sinc = np.where(t == 0, 2 * fc,
25+ np.sin(2 * np.pi * fc * t) / (np.pi * np.where(t == 0, 1, t)))
26+ w = 0.54 - 0.46 * np.cos(2 * np.pi * i / (n - 1))
27+ h = sinc * w
28+ return h / h.sum()
29+
30+
31+def windowed_sinc_bandpass(f_lo, f_hi, taps):
32+ """Difference of two windowed-sinc lowpasses; edges in cycles/sample."""
33+ return windowed_sinc_lowpass(f_hi, taps) - windowed_sinc_lowpass(f_lo, taps)
src/timeseries_entropy/model.pyadded+115−0View file
@@ -0,0 +1,115 @@
1+"""The process and its conditional sampler.
2+
3+Model: x iid N(0, sigma^2) -> y = h * x (causal FIR, kernel length L)
4+-> z = round(y). A zero-phase or otherwise shifted application of the same
5+kernel gives the same law, so causal convolution loses no generality.
6+
7+ConditionalChain targets H(z_{M+1} | z_1..z_M): it fixes an observed past
8+z_1..z_M and Gibbs-samples the latent x under the box constraints
9+y_t in [z_t - 1/2, z_t + 1/2), emitting one exact draw of z_{M+1} per step.
10+"""
11+
12+import numpy as np
13+from scipy.special import ndtr, ndtri
14+
15+
16+def truncated_std_normal(lo, hi, rng):
17+ """Standard normal truncated to [lo, hi], elementwise, by inverse CDF.
18+ Mirrored into the lower tail so the CDF differences keep precision."""
19+ flip = lo > -hi # midpoint above 0 (robust to (-inf, inf) intervals)
20+ a = np.where(flip, -hi, lo)
21+ b = np.where(flip, -lo, hi)
22+ fa = ndtr(a)
23+ fb = ndtr(b)
24+ u = fa + (fb - fa) * rng.random(a.shape)
25+ x = ndtri(np.clip(u, 1e-300, 1 - 1e-16))
26+ x = np.where(flip, -x, x)
27+ return np.clip(x, lo, hi)
28+
29+
30+class ConditionalChain:
31+ """Stationary chain of exact draws of z_{M+1} given a fixed past z_1..z_M.
32+
33+ The constructor draws the past from the prior; the generating latents are
34+ themselves an exact draw from p(x | z), so the Gibbs chain starts in
35+ stationarity — no burn-in bias, only autocorrelation. draw(k) advances the
36+ chain k steps (thin sweeps each) and returns the k sampled z_{M+1} values,
37+ each marginally distributed exactly as z_{M+1} | z_1..z_M.
38+
39+ Each Gibbs conditional x_i | rest is N(0, sigma^2) truncated to the
40+ interval read off the <= L constraint boxes x_i appears in. Coordinates a
41+ multiple of L apart share no constraint, so each of the L "colors"
42+ updates as one vectorized block.
43+ """
44+
45+ def __init__(self, kernel, sigma, past, rng=None, thin=1):
46+ h = np.asarray(kernel, dtype=float)
47+ if h.ndim != 1 or h.size == 0:
48+ raise ValueError('kernel must be a nonempty 1-D array')
49+ if sigma <= 0:
50+ raise ValueError('sigma must be positive')
51+ if past < 1:
52+ raise ValueError('past must be >= 1')
53+ self.h = h
54+ self.sigma = float(sigma)
55+ self.thin = int(thin)
56+ self.rng = np.random.default_rng() if rng is None else rng
57+
58+ L = h.size
59+ M = int(past)
60+ self.L, self.M = L, M
61+
62+ # The past, with its true latents as the (stationary) chain start.
63+ x = self.sigma * self.rng.standard_normal(M + L - 1)
64+ y = np.convolve(x, h, mode='valid')
65+ self.z = np.floor(y + 0.5)
66+ self.x = x
67+
68+ # Boxes and y live in padded arrays so that every coordinate x_i sees
69+ # exactly L constraint rows (rows outside the data are unconstrained).
70+ P = L - 1
71+ self.ypad = np.zeros(M + 2 * P)
72+ self.lo = np.full(M + 2 * P, -np.inf)
73+ self.hi = np.full(M + 2 * P, np.inf)
74+ self.lo[P:P + M] = self.z - 0.5
75+ self.hi[P:P + M] = self.z + 0.5
76+
77+ # Row i+j (padded) carries coefficient h[j] for coordinate i.
78+ self.classes = [np.arange(c0, M + L - 1, L) for c0 in range(L)]
79+ self.rowmats = [idx[:, None] + np.arange(L)[None, :]
80+ for idx in self.classes]
81+ self.nonzero = h != 0
82+ # z_{M+1} = round(hr_head @ x[M:] + h[0] * x_free), x_free fresh.
83+ self.hr_head = h[::-1][:-1]
84+
85+ def _sweep(self):
86+ h, x, sigma = self.h, self.x, self.sigma
87+ M, L = self.M, self.L
88+ P = L - 1
89+ self.ypad[P:P + M] = np.convolve(x, h, mode='valid') # kill fp drift
90+ for idx, rows in zip(self.classes, self.rowmats):
91+ r = self.ypad[rows] - np.outer(x[idx], h)
92+ with np.errstate(divide='ignore', invalid='ignore'):
93+ b1 = (self.lo[rows] - r) / h[None, :]
94+ b2 = (self.hi[rows] - r) / h[None, :]
95+ xlo = np.where(h[None, :] > 0, b1, b2)
96+ xhi = np.where(h[None, :] > 0, b2, b1)
97+ xlo[:, ~self.nonzero] = -np.inf
98+ xhi[:, ~self.nonzero] = np.inf
99+ xlo = xlo.max(axis=1)
100+ xhi = xhi.min(axis=1)
101+ xnew = truncated_std_normal(xlo / sigma, xhi / sigma, self.rng) * sigma
102+ self.ypad[rows] += (xnew - x[idx])[:, None] * h[None, :]
103+ x[idx] = xnew
104+
105+ def draw(self, k):
106+ """The next k samples of z_{M+1}, continuing the chain."""
107+ out = np.empty(k, dtype=np.int64)
108+ M = self.M
109+ for i in range(k):
110+ for _ in range(self.thin):
111+ self._sweep()
112+ c = float(self.hr_head @ self.x[M:]) if self.L > 1 else 0.0
113+ y_next = c + self.sigma * self.h[0] * self.rng.standard_normal()
114+ out[i] = int(np.floor(y_next + 0.5))
115+ return out
tests/test_estimator.pyadded+57−0View file
@@ -0,0 +1,57 @@
1+import math
2+
3+import numpy as np
4+import pytest
5+
6+from timeseries_entropy import plugin_entropy, unbiased_entropy, level_corrections
7+
8+
9+def test_plugin_entropy_exact():
10+ assert plugin_entropy([0, 1]) == pytest.approx(1.0)
11+ assert plugin_entropy([3, 3, 3]) == pytest.approx(0.0)
12+ assert plugin_entropy([0, 0, 0, 1]) == pytest.approx(
13+ -(0.75 * math.log2(0.75) + 0.25 * math.log2(0.25)))
14+
15+
16+def test_unbiased_on_iid_categorical():
17+ p = np.array([0.5, 0.3, 0.2])
18+ true_h = -np.sum(p * np.log2(p))
19+ rng = np.random.default_rng(1)
20+ draw = lambda k: rng.choice(3, size=k, p=p)
21+ vals = np.array([unbiased_entropy(draw, n0=32, r=1.5, rng=rng)
22+ for _ in range(4000)])
23+ se = vals.std(ddof=1) / math.sqrt(len(vals))
24+ assert abs(vals.mean() - true_h) < 5 * se
25+ assert se < 0.01
26+
27+
28+def test_unbiased_on_autocorrelated_chain():
29+ # Stationary two-state Markov chain with sticky transitions; the marginal
30+ # is uniform, so H = 1 bit despite strong autocorrelation.
31+ rng = np.random.default_rng(2)
32+ state = [int(rng.random() < 0.5)]
33+
34+ def draw(k):
35+ out = np.empty(k, dtype=int)
36+ for i in range(k):
37+ if rng.random() < 0.1:
38+ state[0] = 1 - state[0]
39+ out[i] = state[0]
40+ return out
41+
42+ vals = np.array([unbiased_entropy(draw, n0=64, r=1.5, rng=rng)
43+ for _ in range(3000)])
44+ se = vals.std(ddof=1) / math.sqrt(len(vals))
45+ assert abs(vals.mean() - 1.0) < 5 * se
46+ assert se < 0.02
47+
48+
49+def test_level_corrections_decay():
50+ rng = np.random.default_rng(3)
51+ p = np.array([0.5, 0.3, 0.2])
52+ deltas = np.array([
53+ level_corrections(lambda k: rng.choice(3, size=k, p=p), n0=64, levels=5)
54+ for _ in range(200)])
55+ rms = np.sqrt((deltas ** 2).mean(axis=0))
56+ # E[Delta_m^2] should decay roughly like 2^(-2m); demand clear decay.
57+ assert rms[-1] < rms[0] / 4
tests/test_model.pyadded+49−0View file
@@ -0,0 +1,49 @@
1+import math
2+
3+import numpy as np
4+
5+from timeseries_entropy import ConditionalChain, estimate_conditional_entropy
6+from timeseries_entropy.model import truncated_std_normal
7+
8+
9+def quantized_gaussian_entropy_bits(s):
10+ """Exact entropy of round(N(0, s^2)) by direct summation."""
11+ zmax = int(math.ceil(8 * s + 4))
12+ H = 0.0
13+ prev = 0.5 * (1 + math.erf((-zmax - 0.5) / (s * math.sqrt(2))))
14+ for z in range(-zmax, zmax + 1):
15+ cur = 0.5 * (1 + math.erf((z + 0.5) / (s * math.sqrt(2))))
16+ p = cur - prev
17+ prev = cur
18+ if p > 0:
19+ H -= p * math.log2(p)
20+ return H
21+
22+
23+def test_truncated_std_normal_bounds():
24+ rng = np.random.default_rng(0)
25+ lo = np.array([-1.0, 0.5, -np.inf, -8.0])
26+ hi = np.array([-0.5, 2.0, np.inf, -7.0])
27+ for _ in range(100):
28+ x = truncated_std_normal(lo, hi, rng)
29+ assert np.all(x >= lo) and np.all(x <= hi)
30+
31+
32+def test_gibbs_preserves_constraints():
33+ rng = np.random.default_rng(4)
34+ kernel = np.full(3, 1 / 3)
35+ chain = ConditionalChain(kernel, sigma=2.0, past=64, rng=rng)
36+ z0 = chain.z.copy()
37+ chain.draw(50)
38+ y = np.convolve(chain.x, chain.h, mode='valid')
39+ assert np.array_equal(np.floor(y + 0.5), z0)
40+
41+
42+def test_no_filter_matches_exact_entropy():
43+ # kernel [1]: z is iid round(N(0, sigma^2)), so the conditional entropy
44+ # equals the exact marginal entropy for any past length.
45+ sigma = 1.5
46+ exact = quantized_gaussian_entropy_bits(sigma)
47+ est = estimate_conditional_entropy(
48+ [1.0], sigma, past=8, pasts=12, reps=24, n0=64, seed=5)
49+ assert abs(est.mean - exact) < max(5 * est.se, 0.02)