2 * The process and its conditional Gibbs sampler: x iid N(0, σ²) → y = h*x
3 * (causal FIR) → z = round(y), and a stationary chain of exact draws of
4 * z_{M+1} given a fixed past z_1..z_M.
5 *
6 * TypeScript port of model.py from the sibling timeseries-entropy package —
7 * keep the two in step. The numpy version updates same-color coordinate
8 * blocks vectorized; here a sweep is a plain sequential scan over the
9 * coordinates, an equally valid systematic-scan Gibbs sweep. This CPU sweep
10 * is the piece a WebGPU backend would replace.
11 */
12import { ndtr, ndtri } from './normal'
13import type { Rng } from './rng'
15/** Standard normal truncated to [lo, hi], by inverse CDF. Mirrored into the
16 * lower tail so the CDF differences keep precision. */
17export function truncatedStdNormal(lo: number, hi: number, rng: Rng): number {
18 const flip = lo > -hi // midpoint above 0 (robust to (-inf, inf) intervals)
19 const a = flip ? -hi : lo
20 const b = flip ? -lo : hi
21 const fa = ndtr(a)
22 const fb = ndtr(b)
23 const u = Math.min(Math.max(fa + (fb - fa) * rng.uniform(), 1e-300), 1 - 1e-16)
24 let x = ndtri(u)
25 if (flip) x = -x
26 return Math.min(Math.max(x, lo), hi)
27}
29/**
30 * Constructing the chain draws the past from the prior; the generating
31 * latents are themselves an exact draw from p(x | z), so the Gibbs chain
32 * starts in stationarity — no burn-in bias, only autocorrelation. Each Gibbs
33 * conditional x_i | rest is N(0, σ²) truncated to the interval read off the
34 * ≤ L constraint boxes x_i appears in. draw(k) advances the chain k steps
35 * (thin sweeps each) and returns k sampled z_{M+1} values, each marginally
36 * distributed exactly as z_{M+1} | z_1..z_M.
37 */
38export class ConditionalChain {
39 private readonly h: Float64Array
40 private readonly sigma: number
41 private readonly thin: number
42 private readonly rng: Rng
43 private readonly L: number
44 private readonly M: number
45 /** Latents x_0..x_{M+L-2}; boxes live in padded rows so that coordinate i
46 * sees exactly L constraint rows i..i+L-1 (rows outside the data are
47 * unconstrained), with coefficient h[j] in row i+j. */
48 private readonly x: Float64Array
49 private readonly ypad: Float64Array
50 private readonly lo: Float64Array
51 private readonly hi: Float64Array
53 constructor(kernel: Float64Array, sigma: number, past: number, rng: Rng, thin = 1) {
54 if (kernel.length === 0) throw new Error('kernel must be nonempty')
55 if (!(sigma > 0)) throw new Error('sigma must be positive')
56 if (past < 1) throw new Error('past must be >= 1')
57 this.h = kernel
58 this.sigma = sigma
59 this.thin = thin
60 this.rng = rng
61 const L = (this.L = kernel.length)
62 const M = (this.M = past)
64 const x = (this.x = new Float64Array(M + L - 1))
65 for (let i = 0; i < x.length; i++) x[i] = sigma * rng.normal()
67 const P = L - 1
68 this.ypad = new Float64Array(M + 2 * P)
69 this.lo = new Float64Array(M + 2 * P).fill(-Infinity)
70 this.hi = new Float64Array(M + 2 * P).fill(Infinity)
71 this.refreshY()
72 for (let m = 0; m < M; m++) {
73 const z = Math.floor(this.ypad[P + m] + 0.5)
74 this.lo[P + m] = z - 0.5
75 this.hi[P + m] = z + 0.5
76 }
77 }
79 /** ypad[P + m] = y_m = Σ_j h[j] x_{m+L-1-j}, recomputed to kill fp drift. */
80 private refreshY(): void {
81 const { h, x, ypad, L, M } = this
82 const P = L - 1
83 for (let m = 0; m < M; m++) {
84 let y = 0
85 for (let j = 0; j < L; j++) y += h[j] * x[m + L - 1 - j]
86 ypad[P + m] = y
87 }
88 }
90 private sweep(): void {
91 const { h, x, ypad, lo, hi, sigma, rng, L } = this
92 this.refreshY()
93 for (let i = 0; i < x.length; i++) {
94 const xi = x[i]
95 let xlo = -Infinity
96 let xhi = Infinity
97 for (let j = 0; j < L; j++) {
98 const hj = h[j]
99 if (hj === 0) continue
100 const row = i + j
101 const res = ypad[row] - hj * xi
102 const b1 = (lo[row] - res) / hj
103 const b2 = (hi[row] - res) / hj
104 if (hj > 0) {
105 if (b1 > xlo) xlo = b1
106 if (b2 < xhi) xhi = b2
107 } else {
108 if (b2 > xlo) xlo = b2
109 if (b1 < xhi) xhi = b1
110 }
111 }
112 const xn = truncatedStdNormal(xlo / sigma, xhi / sigma, rng) * sigma
113 const d = xn - xi
114 if (d !== 0) {
115 for (let j = 0; j < L; j++) ypad[i + j] += d * h[j]
116 x[i] = xn
117 }
118 }
119 }
121 /** The next k samples of z_{M+1}, continuing the chain. */
122 draw = (k: number): Int32Array => {
123 const { h, x, sigma, rng, L, M } = this
124 const out = new Int32Array(k)
125 for (let s = 0; s < k; s++) {
126 for (let t = 0; t < this.thin; t++) this.sweep()
127 // z_{M+1} = round(c + h[0] · x_free) with x_free ~ N(0, σ²) fresh.
128 let c = 0
129 for (let i = 0; i < L - 1; i++) c += h[L - 1 - i] * x[M + i]
130 out[s] = Math.floor(c + sigma * h[0] * rng.normal() + 0.5)
131 }
132 return out
133 }
134}