2 * Analytic prediction of the entropy rate of z = round(h * x), x iid
3 * N(0, σ²) — the quantization-corrected formula from theory.py of the
4 * sibling timeseries-entropy package (see its docstring for the
5 * derivation); keep the two in step:
6 *
7 * R ≈ G( √( exp( ∫₀¹ ln(σ² |H(f)|² + 1/12) df ) − 1/12 ) ),
8 *
9 * where G(s) is the differential entropy of N(0, s²) + U(−½, ½). Szegő's
10 * one-step prediction error with roundoff as a 1/12 dither floor, fed
11 * through the dithered-quantization entropy — finite at spectral zeros,
12 * saturating to 0 at coarse quantization. Only this corrected prediction
13 * is ported; the high-resolution Szegő form needs polynomial roots and is
14 * not shown in the app.
15 */
16import { ndtr } from './normal'
18const FLOOR = 1 / 12
20/**
21 * Mean over f in [0, 1) of ln(σ² |H(f)|² + 1/12): trapezoid on [0, 1/2]
22 * (the spectrum is symmetric), with |H(f)|² = r₀ + 2 Σ_k r_k cos(2πfk)
23 * from the kernel autocorrelation r, the cosines by Chebyshev recurrence.
24 * n = 2¹⁴ intervals matches the Python 2¹⁸-point FFT integral to ~1e-8
25 * over the app's whole kernel/σ range.
26 */
27export function logSpectrumMean(kernel: Float64Array, sigma: number, n = 1 << 14): number {
28 const L = kernel.length
29 const r = new Float64Array(L)
30 for (let k = 0; k < L; k++) {
31 let s = 0
32 for (let j = 0; j + k < L; j++) s += kernel[j] * kernel[j + k]
33 r[k] = s
34 }
35 const S = new Float64Array(n + 1).fill(r[0])
36 for (let k = 1; k < L; k++) {
37 const w = 2 * r[k]
38 if (w === 0) continue
39 const t = 2 * Math.cos((Math.PI * k) / n) // c_{i+1} = t·c_i − c_{i−1}
40 let cPrev = 1
41 let c = t / 2
42 S[0] += w
43 for (let i = 1; i <= n; i++) {
44 S[i] += w * c
45 const cNext = t * c - cPrev
46 cPrev = c
47 c = cNext
48 }
49 }
50 const s2 = sigma * sigma
51 let sum = 0
52 for (let i = 0; i <= n; i++) {
53 const v = Math.log(s2 * Math.max(S[i], 0) + FLOOR)
54 sum += i === 0 || i === n ? v / 2 : v
55 }
56 return sum / n
57}
59/** G(s) = h(N(0, s²) + U(−½, ½)) in bits — the exact average entropy of
60 * round(c + N(0, s²)) over a uniform grid offset c. */
61export function gaussUniformEntropy(s: number): number {
62 if (s <= 0) return 0
63 if (s < 1e-3) return edgeConstant() * s
64 const dv = Math.min(s / 8, 0.01)
65 const vMax = 0.5 + 8 * s + 1
66 let sum = 0
67 for (let i = 0; i * dv < vMax; i++) {
68 const v = i * dv
69 const g = ndtr((v + 0.5) / s) - ndtr((v - 0.5) / s)
70 const term = g > 0 ? -g * Math.log2(g) : 0
71 sum += i === 0 ? term / 2 : term
72 }
73 return 2 * dv * sum
74}
76let EDGE_C: number | null = null
78/** ∫ h₂(Φ(t)) dt: the small-s slope of G(s). */
79function edgeConstant(): number {
80 if (EDGE_C === null) {
81 const n = 20001
82 let sum = 0
83 for (let i = 0; i < n; i++) {
84 const t = -12 + (24 * i) / (n - 1)
85 const p = Math.min(Math.max(ndtr(t), 1e-300), 1 - 1e-16)
86 const h2 = -(p * Math.log2(p) + (1 - p) * Math.log2(1 - p))
87 sum += i === 0 || i === n - 1 ? h2 / 2 : h2
88 }
89 EDGE_C = (sum * 24) / (n - 1)
90 }
91 return EDGE_C
92}
94/** The predicted entropy rate G(s*) of z = round(h * x), in bits/sample. */
95export function predictEntropyRate(kernel: Float64Array, sigma: number): number {
96 const gmW = Math.exp(logSpectrumMean(kernel, sigma))
97 return gaussUniformEntropy(Math.sqrt(Math.max(gmW - FLOOR, 0)))
98}