Show analytic entropy-rate prediction alongside the Monte-Carlo estimate - concept-collection/timeseries-compressibility

concept-collection / timeseries-compressibility

Show analytic entropy-rate prediction alongside the Monte-Carlo estimate

Port the quantization-corrected prediction from timeseries-entropy's theory.py to src/entropy/theory.ts (matches the Python to ~4e-5 bits across the app's parameter space) and draw it as an always-on dotted reference line on the compression chart. Rework the stat row into two paired readouts — theory and Monte-Carlo ground truth — each keyed to its chart line by a sample of that line's stroke, with implied best ratio beneath. Chart labels move to a reserved band on two rows, the Monte-Carlo line gains a +/- 1 se band, and the bits/sample axis now scales to the data instead of pinning at 16.

Jeremy Magland <jmagland@flatironinstitute.org> committed 19 days ago commit fec22841ed85 parent 9d25fea Browse files

7 changed files+283−42

README.mdmodified+6−2View file

@@ -34,7 +34,11 @@ can beat — is estimated in the browser by the method of the companion

3434 [timeseries-entropy](https://github.com/concept-collection/timeseries-entropy)

3535 package: an unbiased Monte-Carlo estimator of H(z_next | a long past), by Gibbs

3636 sampling the latent Gaussian under the rounding constraints and applying

37-Rhee–Glynn randomized telescoping to the sampled chain. A button starts a web

37+Rhee–Glynn randomized telescoping to the sampled chain. The package's analytic

38+approximation of R — Szegő's one-step prediction error with roundoff as a 1/12

39+dither floor, fed through the Gaussian⊕uniform entropy — is drawn as a dotted

40+reference line at all times, so the Monte-Carlo estimate lands beside its

41+prediction. A button starts a web

3842 worker that averages one independent past at a time (live mean ± se, dashed

3943 line on the chart) until stopped; the app also shows the exact command to run

4044 the Python original at the same settings as an independent check. The

@@ -58,7 +62,7 @@ src/model/ the latent source (fixed seeded randomness indexed by sample

5862 src/entropy/ the unbiased entropy-rate estimator: hand-synced TypeScript

5963 port of the timeseries-entropy package (Gibbs conditional

6064 sampler, Rhee–Glynn telescoping, Cody erfc / Acklam ndtri,

61- xoshiro128** RNG)

65+ xoshiro128** RNG, the analytic rate prediction of theory.py)

6266 src/compress/ lossless codecs run in the browser: zlib (fflate), zstd (wasm),

6367 ans.ts (a bit-identical port of simple_ans), FLAC-style

6468 integer LPC (borrowed from entropy-quantized-linear-transform),

src/App.tsxmodified+26−20View file

@@ -6,6 +6,7 @@ import ScrollingView from './components/ScrollingView'

66 import CompressionChart from './components/CompressionChart'

77 import MethodNote from './components/MethodNote'

88 import { useEntropyRate } from './components/useEntropyRate'

9+import { predictEntropyRate } from './entropy'

910 import { DEFAULT_SPEC, clampSpec, designKernel, kernelNorm } from './model/filters'

1011 import { LATENT_SEED } from './model/latent'

1112 import { DEFAULT_LPC_ORDER, LPC_ORDERS } from './compress/codecs'

@@ -128,6 +129,7 @@ export default function App() {

128129 const sigmaY = useMemo(() => sigma * kernelNorm(kernel), [kernel, sigma])

129130 const compression = useCompression(kernel, sigma, lpcOrder, blockSize)

130131 const entropyRate = useEntropyRate(kernel, sigma)

132+ const theoryBits = useMemo(() => predictEntropyRate(kernel, sigma), [kernel, sigma])

131133

132134 return (

133135 <div className="app">

@@ -160,27 +162,36 @@ export default function App() {

160162

161163 <section className="card">

162164 <h2>Compression</h2>

165+ {/* Two readouts of the same number, each keyed to its chart line by a

166+ sample of that line's own stroke. */}

163167 <div className="stat-row">

164168 <div className="stat">

165- <span className="label">predicted std of z</span>

166- <span className="value">

167- {sigmaY.toFixed(2)} <small>steps</small>

169+ <span className="label">

170+ <span className="line-swatch" style={{ borderTop: '2px dotted var(--theory)' }} />

171+ entropy rate R — analytic theory

168172 </span>

169- </div>

170- <div className="stat">

171- <span className="label">measured std of z</span>

172173 <span className="value">

173- {compression.results.length > 0 ? compression.empiricalStd.toFixed(2) : '…'}{' '}

174- <small>steps</small>

174+ {theoryBits.toFixed(2)} <small>bits/sample</small>

175+ </span>

176+ <span className="stat-sub">

177+ {theoryBits > 0 ? `best possible ratio ${(16 / theoryBits).toFixed(2)}×` : '—'}

175178 </span>

176179 </div>

177180 <div className="stat">

178- <span className="label">entropy rate R</span>

181+ <span className="label">

182+ <span className="line-swatch" style={{ borderTop: '2px dashed var(--ink-2)' }} />

183+ entropy rate R — Monte-Carlo ground truth

184+ </span>

179185 <span className="value">

180186 {entropyRate.mean !== null ? entropyRate.mean.toFixed(2) : '—'}

181187 {entropyRate.se !== null && <small> ± {entropyRate.se.toFixed(2)}</small>}{' '}

182188 <small>bits/sample</small>

183189 </span>

190+ <span className="stat-sub">

191+ {entropyRate.mean !== null && entropyRate.mean > 0

192+ ? `best possible ratio ${(16 / entropyRate.mean).toFixed(2)}×`

193+ : 'run the estimate to check the theory'}

194+ </span>

184195 {/* The estimate lives with its readout: start, watch it refine,

185196 stop; a model change resets it. */}

186197 <span className="stat-action">

@@ -198,14 +209,6 @@ export default function App() {

198209 </span>

199210 </span>

200211 </div>

201- <div className="stat">

202- <span className="label">implied best ratio</span>

203- <span className="value">

204- {entropyRate.mean !== null && entropyRate.mean > 0

205- ? `${(16 / entropyRate.mean).toFixed(2)}×`

206- : '—'}

207- </span>

208- </div>

209212 </div>

210213 {/* Settings of the measurement, not of the model — so they live with

211214 the chart they change rather than in the model bar. */}

@@ -237,15 +240,18 @@ export default function App() {

237240 <CompressionChart

238241 results={compression.results}

239242 rateBits={entropyRate.mean}

243+ rateSe={entropyRate.se}

244+ theoryBits={theoryBits}

240245 computing={compression.computing}

241246 />

242247 )}

243248 <p className="card-note">

244249 Measured on a {blockSize.toLocaleString()}-sample block of the same latent data the

245250 signal view shows; sizes include everything a decoder needs (ANS symbol table, LPC

246- coefficients). Baseline is raw int16 (16 bits/sample). The dashed line is the entropy

247- rate R of the process itself, estimated by the button in the stat row above — the one

248- limit no lossless method whatsoever can beat (see the method section at the bottom).

251+ coefficients). Baseline is raw int16 (16 bits/sample). The two reference lines mark

252+ the entropy rate R of the process — the one limit no lossless method whatsoever can

253+ beat: dotted for the analytic theory, dashed for the Monte-Carlo ground truth, shaded

254+ by its standard error (see the method section at the bottom).

249255 What separates the methods is the model each one codes against: ANS uses the histogram

250256 of whatever stream it is given, so a better prefilter is the only way it improves,

251257 while the conditional-Gaussian coder codes each sample against a prediction and can

src/app.cssmodified+16−0View file

@@ -14,6 +14,7 @@

1414 --series-2: #eb6834;

1515 --series-3: #1baf7a;

1616 --series-4: #8a5cd6;

17+ --theory: #e87ba4;

1718 }

1819 @media (prefers-color-scheme: dark) {

1920 :root {

@@ -30,6 +31,7 @@

3031 --series-2: #d95926;

3132 --series-3: #199e70;

3233 --series-4: #9a6fe8;

34+ --theory: #d55181;

3335 }

3436 }

3537

@@ -278,6 +280,20 @@ body {

278280 color: var(--muted);

279281 }

280282

283+/* A sample of a chart reference line's stroke, keying a stat to its line. */

284+.line-swatch {

285+ display: inline-block;

286+ width: 18px;

287+ margin: 0 6px 3px 0;

288+ vertical-align: middle;

289+}

290+

291+.stat-sub {

292+ font-size: 12px;

293+ color: var(--ink-2);

294+ font-variant-numeric: tabular-nums;

295+}

296+

281297 .stat-action {

282298 display: flex;

283299 align-items: center;

src/components/CompressionChart.tsxmodified+118−19View file

@@ -10,6 +10,9 @@ const COND_VAR = 'var(--series-4)'

1010 const LABEL_W = 100

1111 const RIGHT_PAD = 64

1212 const AXIS_H = 26

13+/** Strip under the axis reserved for the two reference-line labels, so they

14+ * never sit on top of the bars or each other. */

15+const REF_BAND = 30

1316 const GROUP_H = 20

1417 const ROW_H = 24

1518 const BAR_H = 16

@@ -47,6 +50,43 @@ function barPath(x0: number, y: number, len: number, h: number): string {

4750 return `M${x0},${y} h${len - r} a${r},${r} 0 0 1 ${r},${r} v${h - 2 * r} a${r},${r} 0 0 1 ${-r},${r} h${-(len - r)} z`

4851 }

4952

53+/** A reference-line label: a short sample of the line's own stroke, then the

54+ * value in ink, flipped to end-anchored near the right edge. */

55+function RefLabel(props: {

56+ x: number

57+ y: number

58+ width: number

59+ stroke: string

60+ dash: string

61+ text: string

62+}) {

63+ const flip = props.x > props.width - 170

64+ const dir = flip ? -1 : 1

65+ const x0 = props.x + 6 * dir

66+ return (

67+ <g>

68+ <line

69+ x1={x0}

70+ x2={x0 + 16 * dir}

71+ y1={props.y - 4}

72+ y2={props.y - 4}

73+ stroke={props.stroke}

74+ strokeWidth={1.5}

75+ strokeDasharray={props.dash}

76+ />

77+ <text

78+ x={x0 + 20 * dir}

79+ y={props.y}

80+ textAnchor={flip ? 'end' : 'start'}

81+ className="bar-value"

82+ fill="var(--ink)"

83+ >

84+ {props.text}

85+ </text>

86+ </g>

87+ )

88+}

89+

5090 interface Group {

5191 label: string

5292 rows: Row[]

@@ -95,6 +135,10 @@ export default function CompressionChart(props: {

95135 results: CodecResult[]

96136 /** The browser-estimated entropy rate R, once at least one past is in. */

97137 rateBits: number | null

138+ /** Standard error of that estimate, drawn as a band around its line. */

139+ rateSe: number | null

140+ /** The analytic prediction of R, always shown as a dotted reference. */

141+ theoryBits: number

98142 computing: boolean

99143 }) {

100144 const ref = useRef<HTMLDivElement>(null)

@@ -103,7 +147,7 @@ export default function CompressionChart(props: {

103147 const [tip, setTip] = useState<Tip | null>(null)

104148 const [hovered, setHovered] = useState<string | null>(null)

105149

106- const { results, rateBits } = props

150+ const { results, rateBits, rateSe, theoryBits } = props

107151 if (results.length === 0) {

108152 return <p className="card-note">Computing compression on the first block…</p>

109153 }

@@ -116,7 +160,7 @@ export default function CompressionChart(props: {

116160 // row count, so positions accumulate rather than being indexed.

117161 const groupLabels: { label: string; y: number }[] = []

118162 const placed: { row: Row; y: number }[] = []

119- let yCursor = AXIS_H

163+ let yCursor = AXIS_H + REF_BAND

120164 for (const g of groups) {

121165 groupLabels.push({ label: g.label, y: yCursor + 15 })

122166 yCursor += GROUP_H

@@ -131,10 +175,8 @@ export default function CompressionChart(props: {

131175 metric === 'bits' ? r.bitsPerSample : r.ratio

132176 const rateValue =

133177 rateBits !== null && rateBits > 0 ? (metric === 'bits' ? rateBits : 16 / rateBits) : null

134- const xMax =

135- metric === 'bits'

136- ? Math.max(16, ...rows.map(value), rateValue ?? 0) * 1.02

137- : Math.max(...rows.map(value), rateValue ?? 0) * 1.1

178+ const theoryValue = theoryBits > 0 ? (metric === 'bits' ? theoryBits : 16 / theoryBits) : null

179+ const xMax = Math.max(...rows.map(value), rateValue ?? 0, theoryValue ?? 0) * 1.1

138180

139181 const plotW = width - LABEL_W - RIGHT_PAD

140182 const height = yCursor + 6

@@ -151,9 +193,25 @@ export default function CompressionChart(props: {

151193 const rateLabel =

152194 rateBits !== null && rateBits > 0

153195 ? metric === 'bits'

154- ? `R = ${rateBits.toFixed(2)}`

155- : `R ⇒ ${(16 / rateBits).toFixed(2)}×`

196+ ? `Monte-Carlo = ${rateBits.toFixed(2)}`

197+ : `Monte-Carlo ⇒ ${(16 / rateBits).toFixed(2)}×`

156198 : ''

199+ const theoryX = theoryValue !== null ? xOf(theoryValue) : 0

200+ const theoryLabel =

201+ theoryValue !== null

202+ ? metric === 'bits'

203+ ? `theory ≈ ${theoryBits.toFixed(2)}`

204+ : `theory ⇒ ${(16 / theoryBits).toFixed(2)}×`

205+ : ''

206+ // ± one standard error around the Monte-Carlo line, in the plotted metric.

207+ let band: { x: number; w: number } | null = null

208+ if (rateBits !== null && rateBits > 0 && rateSe !== null && rateSe > 0) {

209+ const loBits = Math.max(rateBits - rateSe, 1e-9)

210+ const hiBits = rateBits + rateSe

211+ const x1 = xOf(metric === 'bits' ? loBits : 16 / hiBits)

212+ const x2 = Math.min(xOf(metric === 'bits' ? hiBits : 16 / loBits), LABEL_W + plotW)

213+ band = { x: x1, w: Math.max(x2 - x1, 0) }

214+ }

157215

158216 return (

159217 <div>

@@ -196,6 +254,16 @@ export default function CompressionChart(props: {

196254 </g>

197255 ))}

198256 <line x1={xOf(0)} x2={xOf(0)} y1={AXIS_H - 6} y2={height - 4} stroke="var(--baseline)" strokeWidth={1} />

257+ {band && (

258+ <rect

259+ x={band.x}

260+ y={AXIS_H - 2}

261+ width={band.w}

262+ height={height - 2 - AXIS_H}

263+ fill="var(--ink-2)"

264+ opacity={0.15}

265+ />

266+ )}

199267 {groupLabels.map(g => (

200268 <text key={g.label} x={0} y={g.y} className="bar-group-label">

201269 {g.label}

@@ -230,6 +298,27 @@ export default function CompressionChart(props: {

230298 </g>

231299 )

232300 })}

301+ {theoryValue !== null && (

302+ <g>

303+ <line

304+ x1={theoryX}

305+ x2={theoryX}

306+ y1={AXIS_H - 2}

307+ y2={height - 4}

308+ stroke="var(--theory)"

309+ strokeWidth={1.5}

310+ strokeDasharray="2 3"

311+ />

312+ <RefLabel

313+ x={theoryX}

314+ y={AXIS_H + 11}

315+ width={width}

316+ stroke="var(--theory)"

317+ dash="2 3"

318+ text={theoryLabel}

319+ />

320+ </g>

321+ )}

233322 {rateValue !== null && (

234323 <g>

235324 <line

@@ -241,15 +330,14 @@ export default function CompressionChart(props: {

241330 strokeWidth={1.5}

242331 strokeDasharray="5 4"

243332 />

244- <text

245- x={rateX + (rateX > width - 150 ? -6 : 6)}

246- y={AXIS_H + 10}

247- textAnchor={rateX > width - 150 ? 'end' : 'start'}

248- className="bar-value"

249- fill="var(--ink)"

250- >

251- {rateLabel}

252- </text>

333+ <RefLabel

334+ x={rateX}

335+ y={AXIS_H + 25}

336+ width={width}

337+ stroke="var(--ink-2)"

338+ dash="5 4"

339+ text={rateLabel}

340+ />

253341 </g>

254342 )}

255343 </svg>

@@ -312,11 +400,22 @@ export default function CompressionChart(props: {

312400 <td>{r.ratio.toFixed(3)}</td>

313401 </tr>

314402 ))}

403+ {theoryBits > 0 && (

404+ <tr>

405+ <td>entropy rate R — analytic theory</td>

406+ <td>—</td>

407+ <td>{theoryBits.toFixed(3)}</td>

408+ <td>{(16 / theoryBits).toFixed(3)}</td>

409+ </tr>

410+ )}

315411 {rateBits !== null && rateBits > 0 && (

316412 <tr>

317- <td>entropy rate R (Monte-Carlo)</td>

413+ <td>entropy rate R — Monte-Carlo ground truth</td>

318414 <td>—</td>

319- <td>{rateBits.toFixed(3)}</td>

415+ <td>

416+ {rateBits.toFixed(3)}

417+ {rateSe !== null && rateSe > 0 ? ` ± ${rateSe.toFixed(3)}` : ''}

418+ </td>

320419 <td>{(16 / rateBits).toFixed(3)}</td>

321420 </tr>

322421 )}

src/components/MethodNote.tsxmodified+18−1View file

@@ -32,8 +32,25 @@ export default function MethodNote() {

3232 autocorrelation of the Gibbs draws. Averaging over independent pasts gives R with an

3333 honest standard error.

3434 </p>

35+ <p>

36+ The dotted line on the chart is the package's analytic approximation, computed

37+ instantly from the filter and σ: Szegő's formula gives the error of linearly

38+ predicting the next sample from the past, with roundoff entering as a uniform noise

39+ floor of variance 1/12,

40+ </p>

41+ <Tex

42+ display

43+ tex="R \;\approx\; G(s_*), \qquad s_*^2 \;=\; \exp\!\Big(\int_0^1 \ln\big(\sigma^2\,|H(f)|^2 + \tfrac{1}{12}\big)\,df\Big) \;-\; \tfrac{1}{12},"

44+ />

45+ <p>

46+ where G(s) is the differential entropy of N(0, s²) + U(−½, ½) — exactly the entropy

47+ of a rounded Gaussian averaged over grid offsets. The floor keeps the integral finite

48+ where H(f) has zeros, and G saturates to zero at coarse quantization instead of

49+ diverging. Treating roundoff as independent dither and prediction as linear are

50+ approximations — the Monte-Carlo estimate is the exact check on them.

51+ </p>

3552 <p className="card-note">

36- The estimate button beside the R readout runs exactly this method in a web worker — a

53+ The estimate button under the Monte-Carlo readout runs exactly this method in a web worker — a

3754 TypeScript port of the package (src/entropy, hand-synced), one independent past at a

3855 time until stopped. The command line runs the Python original at the same settings for

3956 an independent check.

src/entropy/index.tsmodified+1−0View file

@@ -13,6 +13,7 @@ export { unbiasedEntropy } from './estimator'

1313 export { ConditionalChain, truncatedStdNormal } from './model'

1414 export { Rng } from './rng'

1515 export { ndtr, ndtri, erfc } from './normal'

16+export { predictEntropyRate, gaussUniformEntropy, logSpectrumMean } from './theory'

1617

1718 export const N0 = 128

1819 export const R_EXPONENT = 1.5

src/entropy/theory.tsadded+98−0View file

@@ -0,0 +1,98 @@

1+/**

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+ */

16+import { ndtr } from './normal'

17+

18+const FLOOR = 1 / 12

19+

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+ */

27+export 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+}

58+

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. */

61+export 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+}

75+

76+let EDGE_C: number | null = null

77+

78+/** ∫ h₂(Φ(t)) dt: the small-s slope of G(s). */

79+function 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+}

93+

94+/** The predicted entropy rate G(s*) of z = round(h * x), in bits/sample. */

95+export 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+}