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
5import os
6from concurrent.futures import ProcessPoolExecutor, as_completed
7from dataclasses import dataclass
9import numpy as np
11from .estimator import plugin_entropy, unbiased_entropy, level_corrections
12from .model import ConditionalChain
13from . import kernels
15__all__ = [
16 'estimate_conditional_entropy', 'Estimate', 'ConditionalChain',
17 'plugin_entropy', 'unbiased_entropy', 'level_corrections', 'kernels',
18]
22class Estimate:
23 """mean +/- se (over independent pasts) of H(z_{M+1} | z_1..z_M), bits."""
24 mean: float
25 se: float
26 per_past: np.ndarray
29def _one_past(kernel, sigma, M, thin, n0, r, reps, seed_seq):
30 rng = np.random.default_rng(seed_seq)
31 chain = ConditionalChain(kernel, sigma, M, rng, thin)
32 return float(np.mean(
33 [unbiased_entropy(chain.draw, n0, r, rng) for _ in range(reps)]))
36def estimate_conditional_entropy(kernel, sigma, past=None, pasts=24, reps=8,
37 n0=128, r=1.5, thin=1, seed=None,
38 progress=None, workers=None):
39 """Unbiased estimate of H(z_{M+1} | z_1..z_M) in bits.
41 For each of `pasts` independent pasts, a stationary Gibbs chain of
42 z_{M+1} draws feeds `reps` randomized-telescoping realizations (on
43 consecutive segments of the chain); their average is one unbiased value
44 per past. Returns the mean and standard error over pasts — valid because
45 pasts are independent.
47 past defaults to max(512, 4 * len(kernel)); the estimand decreases toward
48 the entropy rate as it grows. progress, if given, is called as
49 progress(i, values) after each past finishes (completion order when
50 parallel).
52 Pasts run in parallel across `workers` processes (default: all cores).
53 Each past gets its own spawned RNG stream, so a given seed yields the
54 same result for any worker count.
55 """
56 kernel = np.asarray(kernel, dtype=float)
57 M = int(past) if past is not None else max(512, 4 * kernel.size)
58 seeds = np.random.SeedSequence(seed).spawn(pasts)
59 if workers is None:
60 workers = min(pasts, os.cpu_count() or 1)
61 per_past = np.empty(pasts)
62 values = []
63 if workers <= 1:
64 for i in range(pasts):
65 per_past[i] = _one_past(kernel, sigma, M, thin, n0, r, reps,
66 seeds[i])
67 values.append(per_past[i])
68 if progress is not None:
69 progress(i, values)
70 else:
71 with ProcessPoolExecutor(max_workers=workers) as pool:
72 futures = {
73 pool.submit(_one_past, kernel, sigma, M, thin, n0, r, reps,
74 seeds[i]): i
75 for i in range(pasts)}
76 for done, fut in enumerate(as_completed(futures)):
77 per_past[futures[fut]] = fut.result()
78 values.append(per_past[futures[fut]])
79 if progress is not None:
80 progress(done, values)
81 se = (float(per_past.std(ddof=1) / math.sqrt(len(per_past)))
82 if len(per_past) > 1 else float('nan'))
83 return Estimate(float(per_past.mean()), se, per_past)