v1.31.0

This release makes the Markov-chain Monte Carlo module ran.mc a first-class part of the public API: it adds ten samplers and coordinators — RWM, AdaptiveMetropolis, Slice, HMC, MALA, NUTS, Gibbs, ARS, ParallelTempering, and the runChains multi-chain driver — alongside the gelmanRubin R-hat diagnostic, per-sampler seed()/ess()/exact-resume state(), and tree-shakeable ranjs/mc/* subpath imports. It also lands a broad round of MCMC correctness and robustness fixes and clears the outstanding npm audit devDependency advisories.

Added

• ran.mc.MCMC.state() now round-trips the sampler's complete PRNG-and-adaptation state, making a resumed sampler's subsequent draws bit-for-bit identical to an uninterrupted run instead of merely statistically equivalent (#1033, ADR-0035). state() gains a top-level prng key (the Xoshiro128+ stream position, restored by the constructor via Xoshiro128p.save()/.load(), mirroring ran.dist.Distribution.save()/.load()'s existing prngState precedent), and every subclass's _internal() now also serializes its own proposal/momentum generator (ran.mc.RWM, ran.mc.AdaptiveMetropolis, ran.mc.HMC, ran.mc.NUTS, ran.mc.MALA) and adaptation-batch accumulators — Robbins-Monro counters (RWM, MALA, ran.mc.Slice), the running covariance accumulator (AdaptiveMetropolis), and the dual-averaging and mass-matrix accumulators (HMC, NUTS) — superseding the prior "serialize effective state only" convention (ADR-0020 §2, ADR-0021, ADR-0029) for these specific fields. ran.mc.Gibbs needed no changes: restoring the base class's PRNG alone is sufficient since Gibbs has no subclass-owned generator or adaptation state. Old snapshots (missing the new fields) continue to construct valid instances with the pre-#1033 behavior — additive, non-breaking. Two scope boundaries are documented rather than silently unmet: samplingRate (thinning) is not guaranteed to reproduce exactly across a mid-warm-up resume for any subclass, since it depends on the base-class autocorrelation accumulator, which stays intentionally unserialized (ADR-0023); and RWM's per-dimension proposal scale (_base) is not guaranteed bit-for-bit reproducible across a resume that lands strictly mid-batch-window, because _refreshBase() depends on the same never-serialized base-class Welford accumulator — every other subclass's adaptation state is unaffected by this gap.
• ran.mc.NUTS now reports sampler-health diagnostics, matching the per-iteration divergent/maxTreeDepthReached signals Stan/PyMC/NumPyro expose. Every iterate() result carries a divergent boolean (a leapfrog leaf whose Hamiltonian drifted past the energy-divergence threshold — step size too large or target geometry too extreme) and a maxDepthHit boolean (the doubling tree saturated MAX_TREE_DEPTH without a U-turn — step size too small), and two aggregate accessors, divergenceCount() and maxDepthCount(), report the per-sampling-phase totals. The counts ride the same accumulator lifecycle as ar() (reset at construction and at each sample() start, so a read afterwards reflects the sampling phase only); a well-behaved run reports both as 0. Diagnostic-only — no change to sampling behavior (#1037, ADR-0035).
• ran.mc.NUTS now supports Euclidean metric (mass matrix) adaptation via config.metric, matching ran.mc.HMC: 'diag' (default) adapts a per-dimension variance and 'dense' adapts the full covariance matrix (factored via Matrix.ldl()) during warm-up. Momentum is resampled from N(0, M) instead of a standard Normal, the leapfrog integrator and kinetic energy apply the metric, and the no-U-turn criterion is evaluated on the velocity M⁻¹r; the adapted metric round-trips through state()/_internal(). This removes the previous capability regression where poorly-scaled or correlated targets mixed better under HMC than NUTS (#1035, ADR-0034).
• All 11 ran.mc samplers (AdaptiveMetropolis, ARS, gelmanRubin, Gibbs, HMC, MALA, NUTS, ParallelTempering, runChains, RWM, Slice) are now available as tree-shakeable subpath imports under a dedicated mc namespace (import RWM from 'ranjs/mc/rwm', import gelmanRubin from 'ranjs/mc/gelman-rubin'), matching the per-distribution (ranjs/dist/<name>) and per-process (ranjs/process/<name>) subpath export patterns. Each resolves to a self-contained ESM bundle at dist/mc/<name>.esm.js (#1036).
• ran.mc.AdaptiveMetropolis(logDensity, config, initialState): full-covariance adaptive Metropolis sampler (Haario, Saksman & Tamminen, 2001). Adapts the joint proposal covariance Sigma_proposal = (2.38^2 / dim) * Cov(x) + epsilon * I from the chain's own history during warm-up via an online covariance accumulator and Matrix.ldl(), then freezes the covariance for the sampling phase. Mixes substantially better than RWM's diagonal-only adaptation for correlated multi-dimensional targets (#823).
• ran.mc namespace (RWM, gelmanRubin) is now exported from the library's entry point, wiring it up to ran.mc after it was inadvertently left unexported during PR #615's cleanup (#617).
• seed(value) method on ran.mc.MCMC (and ran.mc.RWM, which additionally reseeds its internal proposal distribution) for deterministic, reproducible sampling. Internally, both classes now use a per-instance Xoshiro128p PRNG instead of the shared module-level generator, so seeding a sampler no longer affects unrelated code sharing that singleton. If the initial position was not explicitly supplied, seed() also redraws it from the newly seeded generator so that .seed(s).sample(n) is fully reproducible (#912).
• ran.mc.runChains(logDensity, config, options): runs multiple independently-seeded RWM chains and computes the gelmanRubin() diagnostic across them in one call — the recommended workflow (ADR-0024) for gating MCMC convergence, since no signal computable from a single chain can distinguish "converged" from "stuck". Defaults to 2 chains seeded [1, 2]; options.chains, options.warmUpBatches, options.sampleSize, options.seeds, and options.maxLength are all configurable. Returns { samples, rhat } (#935).
• ran.mc.Gibbs(conditionals, config, initialState): component-wise (systematic-scan) Gibbs sampler, implemented as an MCMC subclass. Cycles through each dimension in order, replacing it with a draw from the caller-supplied full conditional given the current state. Every draw comes directly from the exact conditional, so there is no accept/reject step and ar() is always 1.0 (#821).
• ran.mc.HMC(logDensity, gradLogDensity, config, initialState): Hamiltonian Monte Carlo sampler, implemented as an MCMC subclass. Uses a leapfrog integrator over config.pathLength steps of size config.stepSize to propose distant moves along Hamiltonian trajectories, with momenta resampled from N(0, I) each iteration and Metropolis accept/reject on the augmented (position, momentum) system. Step size is adapted during warm-up via Robbins-Monro dual averaging (Hoffman & Gelman 2014) toward a target acceptance probability, and jittered multiplicatively (ε ~ Uniform(0.9ε, 1.1ε)) each iteration to avoid periodicity artifacts (#824). Now also supports Euclidean metric (mass matrix) adaptation via config.metric: 'diag' (default) estimates a per-dimension variance online during warm-up so the sampler mixes efficiently on targets whose parameters span very different scales; 'dense' (opt-in) estimates the full covariance matrix via an online accumulator, regularized and factored through Matrix.ldl(), so the sampler also compensates for correlated parameters. The adapted metric round-trips through state()/_internal() alongside stepSize/pathLength (#826).
• ran.mc.MALA({ logDensity, gradLogDensity, config, initialState }): Metropolis-Adjusted Langevin Algorithm sampler, implemented as an MCMC subclass. Proposes a single gradient-informed Langevin step per iteration (x' = x + (stepSize² / 2) * ∇log p(x) + stepSize * z, z ~ N(0, I)) and applies a Metropolis-Hastings correction for the proposal's asymmetry, making the chain exact. Step size is adapted during warm-up via batch Robbins-Monro (Roberts & Rosenthal 2009), the same scheme RWM uses, toward the MALA-optimal 0.574 acceptance rate (Roberts & Rosenthal 1998) (#828). Ships with an options-object-only constructor from its first release — there is no positional form and no deprecation warning, unlike RWM/Slice/AdaptiveMetropolis (#970).
• ran.mc.NUTS({ logDensity, gradLogDensity, config, initialState }): No-U-Turn Sampler, implemented as an MCMC subclass using the identity-mass leapfrog integrator extracted to src/mc/_leapfrog.js, combined with Hoffman & Gelman's (2014) doubling-tree algorithm. Automatically tunes the trajectory length each iteration by recursively extending a leapfrog trajectory forward or backward in a random direction until the trajectory's outer endpoints start turning back toward each other (the U-turn criterion) or a maximum tree depth is reached, selecting the transition via slice sampling over the tree's valid states — eliminating the need to hand-tune pathLength. Step size is adapted during warm-up via the same Robbins-Monro dual averaging as HMC, driven by the tree-averaged acceptance statistic, toward a target acceptance probability of 0.8 (#825). Ships with the options-object constructor form from its first release — no positional form, no deprecation cycle (#972).
• ran.mc.ARS(logDensity, support, derivative): Gilks-Wild (1992) Adaptive Rejection Sampling for univariate log-concave densities on a finite support bracket. Builds a piecewise-exponential upper envelope (and a secant lower "squeeze" hull) from tangent lines to the log-density, adaptively tightening on every rejection so acceptance probability increases monotonically; throws Error for non-log-concave targets. Unlike the rest of ran.mc, it is not an MCMC subclass — it produces exact i.i.d. draws directly, with no warm-up or accept/reject Markov-chain machinery (#820).
• ran.mc.Slice(logDensity, config, initialState): coordinate-wise slice sampler (Neal 2003), implemented as an MCMC subclass. Requires only logDensity — no proposal tuning, no gradient. Each dimension is updated via stepping-out and shrinkage; the interval width w (default 1.0) is the only tunable parameter and is adapted per dimension during warm-up. Every sweep produces an accepted draw, so ar() is always 1.0. A prior, non-functional slice.js (100% commented out, never wired into the base class) was removed as dead code in PR #615; this is a fresh implementation (#822).
• ran.mc.ParallelTempering(logDensity, options): Parallel Tempering / Replica Exchange MCMC (Geyer 1991) for multimodal targets. Runs N independent replica samplers (default RWM, or a caller-supplied sampler factory) at inverse temperatures beta_1 = 1 > beta_2 > ... > beta_n — an explicit options.temperatures array, or an auto-generated geometric ladder from options.nReplicas/options.tempMax. warmUp() tunes every replica independently; sample() runs all replicas in lockstep and, after each thinned step, proposes a swap between one alternating-parity set of adjacent replica pairs, accepted with probability min(1, exp((beta_i - beta_j)(log p(x_j) - log p(x_i)))) per detailed balance on the joint replica distribution, returning the cold (beta = 1) replica's samples; swapRate() reports the accepted/attempted fraction per adjacent pair. Unlike the rest of ran.mc, it is not an MCMC subclass — it coordinates an array of replicas rather than driving a single chain, and does not support state()/resumption (ADR-0028, #830).
• ess() method on ran.mc.MCMC: computes the Effective Sample Size per dimension using Geyer's initial positive monotone sequence estimator (IPSM), N / (-1 + 2 * sum_m Gamma_m), where Gamma_m = rho[2m] + rho[2m+1] pairs consecutive lags starting at lag 0 (from the existing ac() accumulators, so the first pair always includes rho[0] = 1) and is clamped to be non-increasing, summed until the first pair whose clamped value is not positive (falling back to ess = N if even the first pair is non-positive). A fully stuck (zero-variance) chain, where ac() returns NaN at every lag, reports ess = 1 rather than saturating to N. Reads directly from the online accumulators already backing ac() and statistics() — no new accumulator state (#827, #975).

Changed

• ran.mc.runChains() is generalized to drive any ran.mc.MCMC subclass instead of hardcoding RWM: the new signature is runChains(Sampler, samplerOptions, runOptions), where samplerOptions is forwarded verbatim to new Sampler(samplerOptions) for every chain — the same options-object shape that sampler's own constructor accepts ({logDensity, config, initialState} for RWM/AdaptiveMetropolis/Slice, {logDensity, gradLogDensity, config, initialState} for HMC/MALA/NUTS, {conditionals, config, initialState} for Gibbs). runOptions keeps the previous {chains, warmUpBatches, sampleSize, seeds, maxLength} shape. ran.mc.gelmanRubin() is unaffected, since it only ever consumed the returned per-chain sample arrays, never the sampler that produced them. See ADR-0033 (#967).
• ran.mc.RWM now uses a consistent joint diagonal adaptive-Metropolis proposal in both warm-up and sampling instead of tuning per-component (Metropolis-within-Gibbs, 0.44 target) during warm-up and switching to joint proposals for sampling. Warm-up adapts a single global step scale via batch Robbins-Monro toward the optimal acceptance rate (0.44 for dim = 1, 0.234 for dim > 1) and tracks per-component scales from the running marginal standard deviations, so the proposal that is tuned is the proposal that samples. Behavior for dim = 1 is unchanged (the two schemes coincide); multi-dimensional targets are now correctly tuned. See ADR-0022.
• ran.mc.MCMC.ar() now reports the acceptance rate over a sliding window of the most recent config.arWindow iterations (default 1000) instead of the cumulative rate since the last reset, so mid-warmUp() reads aren't dragged down by early untuned batches. During the partial-fill phase (fewer than arWindow iterations since reset) the value is unchanged from before. See ADR-0021 (#920, #926).
• Code Health improved across three files by extracting shared/named helpers: test/special.js (8.28 → 9.09, shared check/checkBesselIdentity/checkF11Recurrence helpers), src/special/marcum-q.js (8.67 → 10.0, _expansionSum/_transitionBand/_initPhi helpers eliminating three Complex Method smells), and test/dist.js (8.76 → 9.09, assertFitSpec/assertParamRecovery helpers eliminating a Complex Method and Excess Arguments smell).
• ran.mc.HMC's class-level documentation and pathLength parameter docs now disclose that a fixed pathLength can still produce genuine resonance-driven negative lag-1 autocorrelation at certain target correlations, even with the existing per-iteration stepSize jitter — confirmed empirically (an investigation swept both target correlation and pathLength, finding resonance bands as narrow as 2-3 integer pathLength steps that a ±10%-scale jitter cannot reliably escape) — and point affected users to ran.mc.NUTS, which adapts trajectory length automatically. No behavior change (#1005).

Deprecated

• ran.mc.ParallelTempering's positional constructor form new ParallelTempering(logDensity, options) is deprecated in favor of the options-object form new ParallelTempering({ logDensity, ...options }), bringing it in line with every other ran.mc sampler and coordinator (RWM, AdaptiveMetropolis, Slice, HMC, MALA, NUTS, Gibbs per ADR-0030; ARS per ADR-0031) and removing the last positional-constructor wart in ran.mc. The positional form still constructs and samples correctly but emits a one-time console.warn on first use; it will be removed in v1.32.0 (#1034).

Removed

• ran.mc.RWM's, ran.mc.AdaptiveMetropolis's, ran.mc.Slice's, ran.mc.HMC's, and ran.mc.Gibbs's deprecated positional constructor arguments (e.g. new RWM(logDensity, config, initialState)) are removed; only the options-object form (new RWM({ logDensity, config, initialState })) remains. ran.mc.runChains()'s deprecated legacy call form, runChains(logDensity, config, options), is likewise removed; only the generalized runChains(Sampler, samplerOptions, runOptions) form remains. This closes out the deprecation cycle introduced in #962–#967 (ADR-0030, ADR-0031, ADR-0033). Note on the deprecation cycle: CLAUDE.md's normal deprecation-cycle rule requires a released minor version containing the warning to ship and hold for a full release before the removal lands; here the #962–#967 ### Deprecated entry never left the [Unreleased] section of this changelog (v1.30.0 predates it), so the removal is landing without that hold, at explicit maintainer request overriding the standard process (#968). No published version of ranjs ever carried the positional forms as deprecated-but-working, so no downstream user is exposed to a behavior change without warning — the positional forms and their console.warn deprecation notices are simply gone, as if they had never been introduced.

Fixed

• ran.mc.RWM, ran.mc.AdaptiveMetropolis, ran.mc.Slice, ran.mc.HMC, and ran.mc.Gibbs constructors now throw a clear, class-specific Error (e.g. "RWM: constructor requires an options object: new RWM({ logDensity, config, initialState })") when called with null, any other non-plain-object argument, or no argument at all, instead of either a generic, engine-dependent TypeError: Cannot destructure property '...' of 'null' as it is not an object. or, for a zero-argument call, silently constructing an unusable instance that only fails later inside _iter(). These five constructors became options-object-only in #968 but never received the guard ran.mc.MALA/ran.mc.NUTS already had (#970, #972), so they regressed to the confusing destructuring error MALA and NUTS were already fixed to avoid (#1029).
• ran.mc.MCMC (and all subclasses, e.g. ran.mc.RWM) now reject config.dim above 10000, config.maxLag above 10000, config.arWindow above 10000, and any individually-valid dim/maxLag combination whose combined accumulator footprint (dim * maxLag * 16 bytes) exceeds 100MB, throwing a clear Error instead of allocating oversized arrays and crashing the process with an out-of-memory error (#916, #922, #928). ran.mc.HMC now rejects config.pathLength above 1024 (2^10, matching the NUTS sampler's own literature-derived MAX_TREE_DEPTH ceiling — the Stan/PyMC/NumPyro default), throwing instead of letting warmUp()/sample() hang indefinitely on the per-iteration leapfrog cost of an unreasonably large path length (#947, #989).
• ran.mc.MCMC warm-up thinning no longer inverts for slow-mixing chains: when a dimension's autocorrelation never decays to ≤ 0.05 within maxLag, _thinningLag() now falls back to the largest measured lag instead of reporting 0. Previously a chain that mixed slower than maxLag could resolve was treated as already-decorrelated, driving samplingRate down toward 1 and under-thinning sample() — the opposite of the intended "slowest-mixing dimension wins" rule (ADR-0020 §3).
• ran.mc.MCMC.warmUp(progress, maxBatches) now runs exactly maxBatches batches (was maxBatches + 1 due to a batch <= maxBatches loop bound) and reports 100 at completion instead of firing a redundant 0% callback at the start.
• ran.mc.MCMC.sample(progress, size) now reports each integer percentage of progress exactly once. Previously the i % (iMax/100) check used a fractional modulus whenever the total iteration count was not a multiple of 100, silently skipping most progress callbacks.
• ran.mc.Gibbs's conditionals now receive the sampler's own PRNG as a second argument (conditionals[d](x, rng)), so seed() can produce reproducible chains for conditionals that draw their randomness from rng.next() instead of an independently-seeded generator. Previously Gibbs._iter() never read this.r, so gibbs.seed(42).sample(null, N) silently failed to reproduce, violating the contract documented on MCMC.seed() (ADR-0026, #938).
• ran.mc.RWM, ran.mc.AdaptiveMetropolis, and ran.mc.Gibbs constructors now have a dedicated JSDoc @param block directly on constructor(), so tsc's generated .d.ts resolves the true parameter types (Function/Function[]) instead of any; Gibbs previously had no constructor signature in the generated declaration at all (#944).
• ran.mc.ARS's hull segment-mass (_build), envelope inverse-CDF (_sampleEnvelope), and tangent-intersection breakpoint (_tangentIntersection) formulas now treat a tangent slope, or a difference between two tangent slopes, as zero once it falls below a Math.cbrt(EPS)-scaled tolerance (matching the noise floor already used elsewhere in the same file for finite-difference-derived slopes), instead of only below raw Number.EPSILON. A small-but-nonzero slope (or near-indistinguishable slope pair) previously fell through to a general-case formula — differencing two nearly-equal exponentials in _build/_sampleEnvelope, or dividing by a near-cancelled denominator in _tangentIntersection — a catastrophic-cancellation pattern that could place a hull breakpoint far outside its valid bracket (#941, #957).
• The docs build's assembleLinks() (docs/src/desc-parser.js) no longer skips every second {@link} construct in a JSDoc paragraph: an off-by-one advanced the loop index by 3 instead of 2 after each converted (text, link) pair, silently dropping the next pair instead of converting it to a hyperlink. src/mc/adaptive-metropolis.js's previously-masked bare {@link ran.mc.RWM} reference is now written in the codebase's standard bracketed [RWM]{@link ran.mc.RWM} form so it renders as a working link instead of tripping the #980 bare-link guard (#997).
• ran.mc.RWM, ran.mc.AdaptiveMetropolis, ran.mc.HMC, ran.mc.MALA, and ran.mc.NUTS no longer alias their proposal/momentum generator with their accept/reject generator after seed(). MCMC._reseedCachedLogDensity() seeded the subclass-owned _q generator with the same raw value passed to this.r; since Xoshiro128p.seed() is a pure function of its argument with no per-instance salt, both generators produced byte-identical streams, so the Metropolis acceptance uniform was a value already consumed to build an earlier proposal — violating the independence the MH ratio assumes. _q is now seeded from a derived value (`${value}-q`), mirroring ParallelTempering's per-replica seeding. Reproducibility is preserved (deterministic derivation).
• ran.mc.HMC and ran.mc.NUTS no longer permanently freeze when the caller's gradient returns NaN at a visited state (e.g. a hand-written gradient that yields NaN near a support boundary instead of a rigorous -Infinity). A non-finite acceptance statistic previously flowed unchecked through the Robbins-Monro dual-averaging recursion into _daHbar → _daLogEpsBar → stepSize, becoming a sticky NaN that silently stopped the sampler from ever moving again. _adjust() now treats a non-finite acceptance statistic as a fully-rejected (divergent) step, driving the step size down so warm-up recovers, matching Stan's divergent-proposal handling.
• ran.mc.Slice now throws for a w (in initialState.internal.w) that is neither a number nor an array (e.g. a string, boolean, object, or null). Such values were silently coerced to the 1.0 default before validation ran, so a documented-parameter type error passed unchecked instead of failing fast per the library's return-value conventions.
• ran.mc.Slice now throws (rather than hanging indefinitely) when logDensity returns NaN at the current point: logY then becomes NaN, lnp(candidate) > logY is always false, and the shrinkage loop narrows forever without accepting. _shrink() is now bounded by a MAX_SHRINK cap — the shrink analogue of the existing w: Infinity stepping-out guard.
• ran.mc.AdaptiveMetropolis's proposal-covariance regularization now scales the epsilon term by s_d, matching Haario, Saksman & Tamminen (2001)'s C_n = s_d * (Cov(x) + epsilon * I). Previously the fixed epsilon = 1e-6 sat outside the s_d = 2.38^2/dim factor (s_d * Cov(x) + epsilon * I), so the regularization floor grew relatively more influential as dimension increased — the opposite of the reference formula's proportional shrinkage.
• ran.mc.gelmanRubin() now throws when the supplied chains do not all have the same length. _gri used chains[0].length as the sample-variance divisor for every chain, so an unequal-length chain (reachable via direct calls, though never via runChains) silently read past its end (undefined → NaN) and mismatched its divisor, producing a wrong or NaN R-hat instead of an error.

Security

• Remediated npm audit findings (#960): @babel/core patched to a version above the arbitrary-file-read range (GHSA-4x5r-pxfx-6jf8) via npm audit fix; nyc bumped ^15.1.0 → ^18.0.0, which pulls a fixed istanbul-lib-processinfo that no longer depends on the vulnerable uuid (GHSA-w5hq-g745-h8pq) — verified against the full test suite, including its coverage-threshold gate; serialize-javascript pinned to ^7.0.7 via a new overrides entry to close mocha's transitive RCE/DoS vulnerabilities (GHSA-5c6j-r48x-rmvq, GHSA-qj8w-gfj5-8c6v), since mocha's own package.json range (^6.0.2) predates the fix even on its latest release. All three changes are devDependency-only; none affect src/ or the published package. Accepted risk, documented and left unresolved because no upstream fix exists: documentation@14.0.3 (latest release) bundles vue-template-compiler@2.7.16 (latest ever published, Vue 2 tooling is EOL) which has an XSS advisory (GHSA-g3ch-rx76-35fx) exploitable only via untrusted template compilation — not applicable here, since npm run docs only compiles maintainer-authored templates locally; the mathjax-node-page toolchain (mathjax, mathjax-node, jsdom, request, request-promise-core/-native, form-data, qs, tough-cookie, nested uuid, yargs/yargs-parser) is abandoned upstream (last publish 2022, itself depending on the long-deprecated request library), so npm audit fix --force's suggested resolution is an older mathjax-node-page release that carries the identical vulnerable subtree — it doesn't fix anything. Both chains are used exclusively by docs/index.js for local, maintainer-invoked API doc generation (npm run docs); docs/ is excluded from the package's files field, and neither dependency runs during npm test, npm run build, or at library runtime. Since nyc@18 declares engines.node: "20 || >=22", the CI test-job matrix (.github/workflows/ci.yml) drops Node 18 (now [20, 22]); this only affects the project's own contributor/CI tooling — the published dist/ bundle carries no Node version requirement.