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
4ad3394Analytic entropy-rate prediction; parallelize pasts across processesJeremy Magland 5import os
6from concurrent.futures import ProcessPoolExecutor, as_completed
8c9a374Unbiased Monte-Carlo entropy estimation for quantized filtered Gaussian seriesJeremy Magland 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
4ad3394Analytic entropy-rate prediction; parallelize pasts across processesJeremy Magland 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)]))
8c9a374Unbiased Monte-Carlo entropy estimation for quantized filtered Gaussian seriesJeremy Magland 36def estimate_conditional_entropy(kernel, sigma, past=None, pasts=24, reps=8,
37 n0=128, r=1.5, thin=1, seed=None,
4ad3394Analytic entropy-rate prediction; parallelize pasts across processesJeremy Magland 38 progress=None, workers=None):
8c9a374Unbiased Monte-Carlo entropy estimation for quantized filtered Gaussian seriesJeremy Magland 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
4ad3394Analytic entropy-rate prediction; parallelize pasts across processesJeremy Magland 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.
8c9a374Unbiased Monte-Carlo entropy estimation for quantized filtered Gaussian seriesJeremy Magland 55 """
56 kernel = np.asarray(kernel, dtype=float)
57 M = int(past) if past is not None else max(512, 4 * kernel.size)
4ad3394Analytic entropy-rate prediction; parallelize pasts across processesJeremy Magland 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)
8c9a374Unbiased Monte-Carlo entropy estimation for quantized filtered Gaussian seriesJeremy Magland 62 values = []
4ad3394Analytic entropy-rate prediction; parallelize pasts across processesJeremy Magland 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)
8c9a374Unbiased Monte-Carlo entropy estimation for quantized filtered Gaussian seriesJeremy Magland 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)