1/**
2 * Rhee–Glynn (randomized telescoping) unbiased entropy estimation.
3 *
4 * TypeScript port of estimator.py from the sibling timeseries-entropy
5 * package — keep the two in step. See that file for the full derivation;
6 * in short: plug-in entropies of blocks whose sizes double per level form a
7 * telescoping sum via the antithetic correction
8 * Δ_m = h(B_m) − [h(B_m¹) + h(B_m²)]/2; truncating at a random level N with
9 * P(N ≥ m) = 2^(−r m) and reweighting by the survival probabilities gives
10 * an estimator whose expectation is exactly the entropy of the stationary
11 * marginal, despite the bias of every finite-block plug-in estimate and any
12 * autocorrelation of the draws. All entropies are in bits.
13 */
14import type { Rng } from './rng'
16/** draw(k) returns the next k consecutive samples of a stationary chain. */
17export type Draw = (k: number) => Int32Array
19function entropyFromCounts(counts: Map<number, number>, n: number): number {
20 let s = 0
21 for (const c of counts.values()) s += c * Math.log2(c)
22 return Math.log2(n) - s / n
23}
25function countInto(counts: Map<number, number>, seg: Int32Array): void {
26 for (let i = 0; i < seg.length; i++) {
27 counts.set(seg[i], (counts.get(seg[i]) ?? 0) + 1)
28 }
29}
31/** h(B_0) and [Δ_1 … Δ_levels] over one growing block; counts merge upward
32 * so the cost is linear in the n0 * 2**levels samples drawn. */
33function telescope(draw: Draw, n0: number, levels: number): { h0: number; deltas: number[] } {
34 const counts = new Map<number, number>()
35 countInto(counts, draw(n0))
36 let size = n0
37 const h0 = entropyFromCounts(counts, size)
38 let hPrev = h0
39 const deltas: number[] = []
40 for (let m = 0; m < levels; m++) {
41 const half = new Map<number, number>()
42 countInto(half, draw(size))
43 const h2 = entropyFromCounts(half, size)
44 for (const [v, c] of half) counts.set(v, (counts.get(v) ?? 0) + c)
45 size *= 2
46 const hFull = entropyFromCounts(counts, size)
47 deltas.push(hFull - 0.5 * (hPrev + h2))
48 hPrev = hFull
49 }
50 return { h0, deltas }
51}
53/**
54 * One randomized-telescoping realization of the marginal entropy (bits).
55 * Consumes n0 * 2**N samples with P(N ≥ m) = 2^(−r m); average many
56 * realizations (they may continue one chain back-to-back) — each has
57 * expectation exactly H. r = 1.5 suits the typical Δ second-moment decay
58 * of 2^(−2m); finite work needs r > 1, finite variance needs decay > r.
59 */
60export function unbiasedEntropy(draw: Draw, n0: number, r: number, rng: Rng): number {
61 const rho = 2 ** -r
62 let N = 0
63 while (rng.uniform() < rho) N++
64 const { h0, deltas } = telescope(draw, n0, N)
65 let est = h0
66 for (let m = 1; m <= N; m++) est += deltas[m - 1] * 2 ** (r * m)
67 return est
68}
70/**
71 * Integrated autocorrelation time of a stationary sequence, in samples.
72 * tau = 1 + 2 Σ_k ρ_k with Sokal's automatic windowing: the sum stops at
73 * the smallest lag W ≥ c·tau(W). Resolving tau needs x.length ≫ c·tau;
74 * longer times saturate near x.length / (2c), so cap the result when the
75 * chain may mix slower than the probe can see. Returns ≥ 1.
76 */
77export function integratedAutocorrTime(x: Int32Array | Float64Array, c = 5): number {
78 const n = x.length
79 if (n < 2) return 1
80 let mean = 0
81 for (let i = 0; i < n; i++) mean += x[i]
82 mean /= n
83 const d = new Float64Array(n)
84 for (let i = 0; i < n; i++) d[i] = x[i] - mean
85 let denom = 0
86 for (let i = 0; i < n; i++) denom += d[i] * d[i]
87 if (denom === 0) return 1
88 let tau = 1
89 let csum = 0
90 for (let w = 1; w <= n >> 1; w++) {
91 let acov = 0
92 for (let i = 0; i + w < n; i++) acov += d[i] * d[i + w]
93 csum += acov / denom
94 tau = 1 + 2 * csum
95 if (w >= c * tau) break
96 }
97 return Math.max(tau, 1)
98}