1import math
3import numpy as np
5from timeseries_entropy import ConditionalChain, estimate_conditional_entropy
6from timeseries_entropy.model import truncated_std_normal
9def 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
23def 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)
32def 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)
42def test_auto_thin_deterministic_across_workers():
43 kwargs = dict(sigma=4.0, past=32, pasts=3, reps=4, n0=32, thin='auto',
44 probe=64, seed=9)
45 kernel = np.full(4, 0.25)
46 a = estimate_conditional_entropy(kernel, workers=1, **kwargs)
47 b = estimate_conditional_entropy(kernel, workers=3, **kwargs)
48 assert np.array_equal(a.per_past, b.per_past)
49 assert np.array_equal(a.thin, b.thin)
50 assert np.array_equal(a.reps, b.reps)
51 assert np.all(a.thin >= 1) and np.all(a.reps >= 1)
54def test_no_filter_matches_exact_entropy():
55 # kernel [1]: z is iid round(N(0, sigma^2)), so the conditional entropy
56 # equals the exact marginal entropy for any past length.
57 sigma = 1.5
58 exact = quantized_gaussian_entropy_bits(sigma)
59 est = estimate_conditional_entropy(
60 [1.0], sigma, past=8, pasts=12, reps=24, n0=64, seed=5)
61 assert abs(est.mean - exact) < max(5 * est.se, 0.02)