8c9a374Unbiased Monte-Carlo entropy estimation for quantized filtered Gaussian seriesJeremy Magland 1import math
3import numpy as np
4import pytest
efed18aAuto-thin each chain to its measured autocorrelation timeJeremy Magland 6from timeseries_entropy import (plugin_entropy, unbiased_entropy,
7 level_corrections, integrated_autocorr_time)
8c9a374Unbiased Monte-Carlo entropy estimation for quantized filtered Gaussian seriesJeremy Magland 8
10def test_plugin_entropy_exact():
11 assert plugin_entropy([0, 1]) == pytest.approx(1.0)
12 assert plugin_entropy([3, 3, 3]) == pytest.approx(0.0)
13 assert plugin_entropy([0, 0, 0, 1]) == pytest.approx(
14 -(0.75 * math.log2(0.75) + 0.25 * math.log2(0.25)))
17def test_unbiased_on_iid_categorical():
18 p = np.array([0.5, 0.3, 0.2])
19 true_h = -np.sum(p * np.log2(p))
20 rng = np.random.default_rng(1)
21 draw = lambda k: rng.choice(3, size=k, p=p)
22 vals = np.array([unbiased_entropy(draw, n0=32, r=1.5, rng=rng)
23 for _ in range(4000)])
24 se = vals.std(ddof=1) / math.sqrt(len(vals))
25 assert abs(vals.mean() - true_h) < 5 * se
26 assert se < 0.01
29def test_unbiased_on_autocorrelated_chain():
30 # Stationary two-state Markov chain with sticky transitions; the marginal
31 # is uniform, so H = 1 bit despite strong autocorrelation.
32 rng = np.random.default_rng(2)
33 state = [int(rng.random() < 0.5)]
35 def draw(k):
36 out = np.empty(k, dtype=int)
37 for i in range(k):
38 if rng.random() < 0.1:
39 state[0] = 1 - state[0]
40 out[i] = state[0]
41 return out
43 vals = np.array([unbiased_entropy(draw, n0=64, r=1.5, rng=rng)
44 for _ in range(3000)])
45 se = vals.std(ddof=1) / math.sqrt(len(vals))
46 assert abs(vals.mean() - 1.0) < 5 * se
47 assert se < 0.02
efed18aAuto-thin each chain to its measured autocorrelation timeJeremy Magland 50def test_autocorr_time_iid():
51 rng = np.random.default_rng(6)
52 tau = integrated_autocorr_time(rng.integers(0, 4, size=4096))
53 assert 0.8 < tau < 1.5
56def test_autocorr_time_sticky_markov():
57 # Two-state chain flipping w.p. eps: rho_k = (1 - 2 eps)^k, so
58 # tau = 1 + 2 rho / (1 - rho) = 19 at eps = 0.05.
59 rng = np.random.default_rng(7)
60 flips = rng.random(200000) < 0.05
61 x = np.cumsum(flips) % 2
62 tau = integrated_autocorr_time(x)
63 assert 12 < tau < 28
66def test_autocorr_time_degenerate():
67 assert integrated_autocorr_time([3]) == 1.0
68 assert integrated_autocorr_time([2, 2, 2, 2]) == 1.0
8c9a374Unbiased Monte-Carlo entropy estimation for quantized filtered Gaussian seriesJeremy Magland 71def test_level_corrections_decay():
72 rng = np.random.default_rng(3)
73 p = np.array([0.5, 0.3, 0.2])
74 deltas = np.array([
75 level_corrections(lambda k: rng.choice(3, size=k, p=p), n0=64, levels=5)
76 for _ in range(200)])
77 rms = np.sqrt((deltas ** 2).mean(axis=0))
78 # E[Delta_m^2] should decay roughly like 2^(-2m); demand clear decay.
79 assert rms[-1] < rms[0] / 4