Auto-thin entropy chains to their measured autocorrelation time - concept-collection/timeseries-compressibility

concept-collection / timeseries-compressibility

Auto-thin entropy chains to their measured autocorrelation time

Hand-sync with timeseries-entropy: with thin=1, slow-mixing narrowband large-sigma settings violate the telescoping estimator's r=1.5 variance condition, producing rare realizations of hundreds of bits. Each past now probes the chain's integrated autocorrelation time over 512 draws, thins by ceil(tau) capped at 64, and treats REPS_PER_PAST as a budget (realizations = max(1, budget/thin)) so per-past cost stays roughly flat. Progress messages carry the resolved rep count.

Jeremy Magland <jmagland@flatironinstitute.org> committed 18 days ago commit 581f213935ba parent 6139258 Browse files

4 changed files+69−14

src/entropy/estimator.tsmodified+30−0View file

@@ -66,3 +66,33 @@ export function unbiasedEntropy(draw: Draw, n0: number, r: number, rng: Rng): nu

6666 for (let m = 1; m <= N; m++) est += deltas[m - 1] * 2 ** (r * m)

6767 return est

6868 }

69+

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

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

src/entropy/index.tsmodified+31−11View file

@@ -5,11 +5,11 @@

55 * defaults match its CLI so browser runs and `timeseries-entropy` runs

66 * target the same estimand with the same variance behavior.

77 */

8-import { unbiasedEntropy } from './estimator'

8+import { integratedAutocorrTime, unbiasedEntropy } from './estimator'

99 import { ConditionalChain } from './model'

1010 import { Rng } from './rng'

1111

12-export { unbiasedEntropy } from './estimator'

12+export { unbiasedEntropy, integratedAutocorrTime } from './estimator'

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

1414 export { Rng } from './rng'

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

@@ -17,8 +17,13 @@ export { predictEntropyRate, gaussUniformEntropy, logSpectrumMean } from './theo

1717

1818 export const N0 = 128

1919 export const R_EXPONENT = 1.5

20+/** Per-past realization budget at thin = 1; the resolved thin divides it. */

2021 export const REPS_PER_PAST = 8

21-export const THIN = 1

22+/** Draws taken at thin = 1 to measure each chain's autocorrelation time. */

23+export const PROBE = 512

24+/** Cap on the auto-resolved thin (bounds cost; also the probe cannot

25+ * resolve times much beyond PROBE / 10). */

26+export const THIN_CAP = 64

2227

2328 /** The conditioning window M at a given kernel length, as in the CLI. */

2429 export function defaultPast(kernelLength: number): number {

@@ -31,22 +36,37 @@ export function pastSeed(baseSeed: number, pastIndex: number): number {

3136 return (baseSeed + Math.imul(0x9e3779b9, pastIndex + 1)) >>> 0

3237 }

3338

34-/** One independent past's unbiased estimate: a fresh stationary chain,

35- * averaged over reps randomized-telescoping realizations run back-to-back

36- * on it. Averaging these over pasts estimates H(z_{M+1} | z_1..z_M). */

39+/**

40+ * One independent past's unbiased estimate: a fresh stationary chain,

41+ * auto-thinned to its measured mixing, averaged over randomized-telescoping

42+ * realizations run back-to-back on it. Averaging these over pasts estimates

43+ * H(z_{M+1} | z_1..z_M).

44+ *

45+ * Auto-thinning mirrors thin='auto' in the Python package: a PROBE-draw

46+ * pilot at thin = 1 estimates the chain's integrated autocorrelation time

47+ * tau, the chain then takes ceil(tau) sweeps per draw (capped at THIN_CAP),

48+ * and REPS_PER_PAST acts as a budget — realizations = max(1, budget/thin) —

49+ * so per-past cost stays roughly flat. Without this, slowly mixing chains

50+ * (narrowband kernels x large sigma) make the level corrections decay too

51+ * slowly for R_EXPONENT and the estimator's variance is infinite: still

52+ * unbiased, but rare realizations of hundreds of bits.

53+ */

3754 export function estimateOnePast(

3855 kernel: Float64Array,

3956 sigma: number,

4057 past: number,

4158 seed: number,

42- onRep?: (repsDone: number) => void,

59+ onRep?: (repsDone: number, reps: number) => void,

4360 ): number {

4461 const rng = new Rng(seed)

45- const chain = new ConditionalChain(kernel, sigma, past, rng, THIN)

62+ const chain = new ConditionalChain(kernel, sigma, past, rng, 1)

63+ const tau = integratedAutocorrTime(chain.draw(PROBE))

64+ chain.thin = Math.min(THIN_CAP, Math.max(1, Math.ceil(tau)))

65+ const reps = Math.max(1, Math.round(REPS_PER_PAST / chain.thin))

4666 let sum = 0

47- for (let rep = 0; rep < REPS_PER_PAST; rep++) {

67+ for (let rep = 0; rep < reps; rep++) {

4868 sum += unbiasedEntropy(chain.draw, N0, R_EXPONENT, rng)

49- onRep?.(rep + 1)

69+ onRep?.(rep + 1, reps)

5070 }

51- return sum / REPS_PER_PAST

71+ return sum / reps

5272 }

src/entropy/model.tsmodified+4−1View file

@@ -38,7 +38,10 @@ export function truncatedStdNormal(lo: number, hi: number, rng: Rng): number {

3838 export class ConditionalChain {

3939 private readonly h: Float64Array

4040 private readonly sigma: number

41- private readonly thin: number

41+ /** Sweeps per emitted sample. May be reassigned between draws (e.g. probe

42+ * at thin = 1, then thin by the measured autocorrelation time);

43+ * stationarity is unaffected. */

44+ thin: number

4245 private readonly rng: Rng

4346 private readonly L: number

4447 private readonly M: number

src/worker/entropyWorker.tsmodified+4−2View file

@@ -27,9 +27,11 @@ const post = self.postMessage as (message: EntropyUpdate) => void

2727 self.onmessage = (e: MessageEvent<EntropyRequest>) => {

2828 const { kernel, sigma, past, seed, startPast } = e.data

2929 for (let i = startPast; ; i++) {

30+ // The rep count is only known once the past's mixing probe resolves the

31+ // thinning; until the first onRep, report the thin=1 budget.

3032 post({ type: 'progress', pastIndex: i, repsDone: 0, reps: REPS_PER_PAST })

31- const value = estimateOnePast(kernel, sigma, past, pastSeed(seed, i), repsDone =>

32- post({ type: 'progress', pastIndex: i, repsDone, reps: REPS_PER_PAST }),

33+ const value = estimateOnePast(kernel, sigma, past, pastSeed(seed, i), (repsDone, reps) =>

34+ post({ type: 'progress', pastIndex: i, repsDone, reps }),

3335 )

3436 post({ type: 'past', pastIndex: i, value })

3537 }