8c9a374Unbiased Monte-Carlo entropy estimation for quantized filtered Gaussian seriesJeremy Magland 1"""Unbiased Monte-Carlo entropy estimation for quantized filtered Gaussian
2time series: x iid N(0, sigma^2) -> y = h * x -> z = round(y)."""
4import math
5from dataclasses import dataclass
7import numpy as np
9from .estimator import plugin_entropy, unbiased_entropy, level_corrections
10from .model import ConditionalChain
11from . import kernels
13__all__ = [
14 'estimate_conditional_entropy', 'Estimate', 'ConditionalChain',
15 'plugin_entropy', 'unbiased_entropy', 'level_corrections', 'kernels',
16]
20class 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
27def 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.
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.
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)