Sampling

Prelude.Sampling provides the moves of the common Monte Carlo samplers. Metropolis–Hastings, MALA, Hamiltonian Monte Carlo, slice sampling, Gibbs sampling, parallel tempering and simulated annealing are all the same loop with a different proposal: carry the state, read a proposal draw and an acceptance draw, and accept by a 0/1 mask rather than a branch. Each definition here is one step of that loop. It also covers rejection and importance sampling with the effective sample size, antithetic and control variates, the pieces of a bootstrap particle filter, and the stick-breaking construction behind Dirichlet-process and related nonparametric priors.

Every random number is an argument. A step takes, for example, a standard normal draw z and a uniform draw u in [0, 1), and returns the next state; you generate the draws in C# and run the chain by calling the step repeatedly (or wrap it in a recursive definition of your own). Every definition works elementwise, so passing vectors of states and draws advances many independent chains in one call.

The important limitation: the target is a Gaussian. This language has no function-valued parameters, so the samplers cannot take an arbitrary density. Each one targets the normal distribution \(N(\mu, \sigma^2)\) whose mean mu and standard deviation sigma are arguments (Gibbs targets a standard bivariate normal with correlation rho). That is enough to learn, test and calibrate the samplers against a known answer; it is not a way to sample an arbitrary posterior.

Notation used in the statements below. Unless a definition sums over a vector, it acts entrywise, so each statement is written for real scalars and applies to every entry of a vector argument. Many divisors pass through the guard nonzero of Prelude.Chains,

\[\operatorname{nonzero}(v) = v + 10^{-300} \, [\, \lvert v \rvert < 10^{-300} \,],\]

which leaves every divisor of magnitude at least \(10^{-300}\) unchanged. For the Gaussian target we write \(\ell(x) = \operatorname{logTarget}(x, \mu, \sigma)\) and \(g(x) = \operatorname{slopeOf}(x, \mu, \sigma)\).

open Prelude.Sampling

The target and its gradient

The target’s log density, up to an additive constant, and its derivative:

\[\log \pi(x) = -\frac{(x - \mu)^2}{2\sigma^2}, \qquad \frac{d}{dx}\log \pi(x) = \frac{\mu - x}{\sigma^2}\]

A standard deviation of zero is nudged to a tiny positive value rather than dividing by zero.

logTarget

logTarget(x, mu, sigma) · x state scalar or vector · mu mean scalar · sigma standard deviation scalar → same shape as x

Let \(x \in \mathbb{R}\) be a state, \(\mu \in \mathbb{R}\) the mean and \(\sigma > 0\) the standard deviation. Then

\[\operatorname{logTarget}(x, \mu, \sigma) = \ell(x) = -\frac{(x - \mu)^2}{\operatorname{nonzero}(2\sigma^2)} .\]

The unnormalised log density \(\log \pi(x)\) of \(N(\mu, \sigma^2)\). Its maximum is 0, at \(x = \mu\).

Intuition MCMC never needs a density’s normalising constant, only how much more likely one point is than another. Working in logs turns those ratios into differences and avoids underflow far in the tails; dropping the constant makes the peak exactly 0, so values read as “log-likelihood relative to the best point”.

Use it when you compute an acceptance ratio by hand, a tempered density with tempered, or a Gaussian likelihood (it is exactly what particleWeights exponentiates). It is the only target the samplers here know, so it is most useful for testing a sampler against a known answer rather than for real posteriors.

density-evaluationsampler-validation

def logTarget(x, mu, sigma)  ≔ 0 − (x − mu) · (x − mu) / nonzero(2 · sigma · sigma)

slopeOf

slopeOf(x, mu, sigma) · x state scalar or vector · mu mean scalar · sigma standard deviation scalar → same shape as x

Let \(x \in \mathbb{R}\), \(\mu \in \mathbb{R}\) and \(\sigma > 0\). Then

\[\operatorname{slopeOf}(x, \mu, \sigma) = g(x) = \frac{\mu - x}{\operatorname{nonzero}(\sigma^2)} .\]

The gradient of the log density, \((\mu - x)/\sigma^2\), which MALA and Hamiltonian Monte Carlo follow uphill.

Intuition An arrow pointing toward higher probability, longer the further you are from the mode and the narrower the target. For a Gaussian it is a spring pulling back to \(\mu\) with stiffness \(1/\sigma^2\).

Use it when you want proposals that head where the probability is instead of wandering blindly — it is the force in malaNext and hmcNext — or when you check a gradient-based sampler’s arithmetic against a closed form. Random-walk hastingsNext and sliceNext do not need it.

gradient-based-samplingsampler-validation

def slopeOf(x, mu, sigma)    ≔ (mu − x) / nonzero(sigma · sigma)

keeps

keeps(logRatio, u) · logRatio log acceptance ratio scalar or vector · u uniform draw in [0, 1) → mask

Let \(\lambda \in \mathbb{R}\) (the argument logRatio) be a log acceptance ratio and let \(u \in [0, 1)\) be a uniform draw. Then

\[\operatorname{keeps}(\lambda, u) = \big[\, u < e^{\min(\lambda, 0)} \,\big] .\]

The Metropolis acceptance test: 1 when \(u < \min(1, e^{\text{logRatio}})\), 0 otherwise. Every sampler below accepts through it.

Intuition Always take a move that is at least as good; take a worse one with probability equal to how much worse it is. That balance — every step into a region matched by steps out in proportion to its probability — is what makes the chain spend time exactly according to the target. Returning a 0/1 mask lets a vector of chains accept or reject independently without a branch.

Use it when you write any Metropolis-style step, including ones not on this page: pass the log ratio you computed (with any Hastings correction) and a fresh uniform. It also accepts parallel-tempering swaps (swapRatio). For a cost to minimise at a temperature use annealAccepts.

bayesian-posterior-samplingaccept-reject-decisions

def keeps(logRatio, u)     ≔ u < exp(min(logRatio, 0))

Metropolis–Hastings

Random-walk Metropolis proposes \(x' = x + \text{step}\cdot z\) and accepts with probability \(\min(1, \pi(x')/\pi(x))\). The proposal is symmetric, so the Hastings correction cancels and the ratio is the target alone.

hastingsNext

hastingsNext(x, step, mu, sigma, z, u) · x current state · step proposal scale scalar · mu mean · sigma standard deviation · z standard normal draw · u uniform draw → next state

Let \(x \in \mathbb{R}\) be the current state, \(h > 0\) (the argument step) the proposal scale, \(\mu \in \mathbb{R}\) and \(\sigma > 0\), let \(z \in \mathbb{R}\) be a standard normal draw and \(u \in [0, 1)\) a uniform draw, and let \(x' = x + h z\). Then

\[\operatorname{hastingsNext}(x, h, \mu, \sigma, z, u) = x + (x' - x) \, \big[\, u < e^{\min(\ell(x') - \ell(x),\, 0)} \,\big] .\]

One random-walk Metropolis step: the proposal if it is accepted, otherwise x unchanged. Shapes follow x; pass vectors of x, z and u to step many chains at once.

Intuition A blind hiker who steps in a random direction and keeps the step if it goes uphill, or sometimes if it goes downhill. Only the target’s ratio is used, so it works where nothing but the (unnormalised) density can be evaluated.

Use it when you want the simplest, most robust MCMC baseline, or a reference to compare other samplers against. Tune step so that roughly a quarter to a half of proposals are accepted; in one dimension about \(2.4\sigma\) is near optimal. Too small and the chain crawls, too large and it rarely moves. It mixes slowly in many dimensions, where malaNext or hmcNext use the gradient to do better.

bayesian-posterior-samplingmonte-carlo-estimationstep-size-tuning

def hastingsNext(x, step, mu, sigma, z, u) ≔ x + (walkProposal(x, step, z) − x) · keeps(hastingsRatio(x, walkProposal(x, step, z), mu, sigma), u)

Building blocks

walkProposal(x, step, z) — Let \(x \in \mathbb{R}\), \(h > 0\) (the argument step) and \(z \in \mathbb{R}\) a standard normal draw. Then \(\operatorname{walkProposal}(x, h, z) = x + h z\). The random-walk proposal, x plus step times a normal draw.
hastingsRatio(x, y, mu, sigma) — Let \(x, y \in \mathbb{R}\), \(\mu \in \mathbb{R}\) and \(\sigma > 0\). Then \(\operatorname{hastingsRatio}(x, y, \mu, \sigma) = \ell(y) - \ell(x)\). The log acceptance ratio of moving from x to y.

MALA

The Metropolis-adjusted Langevin algorithm drifts the proposal uphill along the gradient, so more proposals are worth keeping:

\[x' = x + \frac{h}{2} \nabla \log \pi(x) + \sqrt{h}\, z\]

with step size \(h\). This proposal is not symmetric, so the Hastings correction \(q(x \mid x') / q(x' \mid x)\) must be carried. Leaving it out is the classic MALA mistake: the chain still moves and looks plausible, but converges to the wrong distribution.

malaNext

malaNext(x, step, mu, sigma, z, u) · x current state · step step size h scalar · mu mean · sigma standard deviation · z standard normal draw · u uniform draw → next state

Let \(x \in \mathbb{R}\) be the current state, \(h > 0\) (the argument step) the step size, \(\mu \in \mathbb{R}\) and \(\sigma > 0\), let \(z \in \mathbb{R}\) be a standard normal draw and \(u \in [0, 1)\) a uniform draw. Let \(m(a) = a + \tfrac{h}{2} g(a)\), \(\log q(b \mid a) = -\dfrac{(b - m(a))^2}{\operatorname{nonzero}(2h)}\), and \(x' = m(x) + \sqrt{h}\, z\). Then

\[\operatorname{malaNext}(x, h, \mu, \sigma, z, u) = x + (x' - x) \, \big[\, u < e^{\min(\rho,\, 0)} \,\big], \qquad \rho = \ell(x') - \ell(x) + \log q(x \mid x') - \log q(x' \mid x) .\]

One MALA step with the full Hastings correction: the drifted proposal if accepted, otherwise x.

Intuition A random walk with a tilt: first slide downhill in energy (uphill in probability) along the gradient, then add noise. That is a discretised Langevin diffusion, whose stationary distribution is the target; the Metropolis test with the Hastings correction repairs the error the discretisation introduces.

Use it when the gradient is available and a random walk mixes too slowly, especially as dimension grows. Aim for an acceptance rate around 50–60% and shrink step if it falls well below. The step is a variance, so it scales like \(\sigma^2\), not \(\sigma\). For longer moves per step use hmcNext; if no gradient exists use hastingsNext or sliceNext.

bayesian-posterior-samplinggradient-based-samplingstep-size-tuning

def malaNext(x, step, mu, sigma, z, u) ≔ x + (malaProposal(x, step, mu, sigma, z) − x) · keeps(malaRatio(x, malaProposal(x, step, mu, sigma, z), step, mu, sigma), u)

Building blocks

malaMean(x, step, mu, sigma) — Let \(x \in \mathbb{R}\), \(h > 0\) (the argument step), \(\mu \in \mathbb{R}\) and \(\sigma > 0\). Then \(\operatorname{malaMean}(x, h, \mu, \sigma) = m(x) = x + \tfrac{h}{2} g(x)\). The centre of the proposal, x moved half a step along the gradient.
malaProposal(x, step, mu, sigma, z) — Let \(x\), \(h\), \(\mu\), \(\sigma\) be as for malaMean and \(z \in \mathbb{R}\) a standard normal draw. Then \(\operatorname{malaProposal}(x, h, \mu, \sigma, z) = m(x) + \sqrt{h}\, z\). The drifted proposal, the centre plus sqrt(step) times a normal draw.
malaLogQ(from, to, step, mu, sigma) — Let \(a, b \in \mathbb{R}\) (the arguments from, to), \(h > 0\), \(\mu \in \mathbb{R}\) and \(\sigma > 0\). Then \(\operatorname{malaLogQ}(a, b, h, \mu, \sigma) = \log q(b \mid a) = -(b - m(a))^2 / \operatorname{nonzero}(2h)\). The unnormalised log density of proposing to from from.
malaRatio(x, y, step, mu, sigma) — Let \(x, y \in \mathbb{R}\), \(h > 0\), \(\mu \in \mathbb{R}\) and \(\sigma > 0\). Then \(\operatorname{malaRatio}(x, y, h, \mu, \sigma) = \ell(y) - \ell(x) + \log q(x \mid y) - \log q(y \mid x)\). The log acceptance ratio including the Hastings correction.

Hamiltonian Monte Carlo

Hamiltonian Monte Carlo gives the state a momentum \(p\) and simulates the dynamics of the energy

\[H(x, p) = -\log \pi(x) + \tfrac{1}{2} p^2\]

with the leapfrog integrator (half kick, drift, half kick), then accepts with probability \(\min(1, e^{H(x,p) - H(x',p')})\). Each call takes one leapfrog step from a momentum draw you supply (a fresh standard normal per call). Leapfrog is symplectic, so for a sensible step size the energy wanders but does not drift; a step that is too large makes it diverge, which a good acceptance rate alone can hide. Use diverged to detect that.

In this section let \(x \in \mathbb{R}\) be the position, \(p \in \mathbb{R}\) the momentum, \(h > 0\) (the argument step) the leapfrog step size, \(\mu \in \mathbb{R}\) and \(\sigma > 0\), and write one leapfrog step as

\[x' = x + h \big( p + \tfrac{h}{2} g(x) \big), \qquad p' = p + \tfrac{h}{2} g(x) + \tfrac{h}{2} g(x') .\]

energyOf

energyOf(x, p, mu, sigma) · x position · p momentum · mu mean · sigma standard deviation → energy

Let \(x, p \in \mathbb{R}\), \(\mu \in \mathbb{R}\) and \(\sigma > 0\). Then

\[\operatorname{energyOf}(x, p, \mu, \sigma) = H(x, p) = -\ell(x) + \tfrac{1}{2} p^2 .\]

The total energy \(H(x, p)\): potential \(-\log \pi(x)\) plus kinetic \(p^2/2\).

Intuition Picture a frictionless puck on a bowl whose height is \(-\log \pi\): low where probability is high. Its speed trades height for motion, and the total energy is conserved by exact dynamics. Any change in \(H\) after a leapfrog step is purely integration error, which is why the acceptance test is built from it.

Use it when you diagnose Hamiltonian Monte Carlo: track the energy change across steps (small and unbiased means the step size is fine), or compute the acceptance ratio for a longer trajectory you integrate yourself. diverged wraps the check for one step.

sampler-validationgradient-based-sampling

def energyOf(x, p, mu, sigma)           ≔ 0 − logTarget(x, mu, sigma) + p · p / 2

hmcNext

hmcNext(x, p, step, mu, sigma, u) · x current position · p momentum, a standard normal draw · step leapfrog step size scalar · mu mean · sigma standard deviation · u uniform draw → next position

Let \(x\), \(p\), \(h\), \(\mu\), \(\sigma\) and \(x', p'\) be as above, with \(p\) a standard normal draw, and let \(u \in [0, 1)\) be a uniform draw. Then

\[\operatorname{hmcNext}(x, p, h, \mu, \sigma, u) = x + (x' - x) \, \big[\, u < e^{\min(H(x, p) - H(x', p'),\, 0)} \,\big] .\]

One Hamiltonian Monte Carlo step with a single leapfrog step: the new position if the energy test accepts it, otherwise x.

Intuition Give the state a random kick of momentum and let it coast along the bowl of \(-\log \pi\) for a moment. Physics carries it toward and across regions of high probability instead of diffusing, and the energy test corrects the integrator’s error.

Use it when you study or test HMC mechanics against a known Gaussian. Draw a fresh standard normal p every call. Note that with a single leapfrog step HMC is equivalent to MALA (with MALA step \(h^2\)), so the long, far-reaching moves that make HMC powerful need many leapfrog steps chained by your own code. Watch diverged when tuning step.

bayesian-posterior-samplinggradient-based-sampling

def hmcNext(x, p, step, mu, sigma, u)   ≔ x + (leapX(x, p, step, mu, sigma) − x) · keeps(hmcRatio(x, p, step, mu, sigma), u)

diverged

diverged(x, p, step, mu, sigma) · x position · p momentum · step step size · mu mean · sigma standard deviation → mask

Let \(x\), \(p\), \(h\), \(\mu\), \(\sigma\) and \(x', p'\) be as above. Then

\[\operatorname{diverged}(x, p, h, \mu, \sigma) = \big[\, \lvert H(x, p) - H(x', p') \rvert > 10 \,\big] .\]

1 when the leapfrog step changes the energy by more than 10, the sign of a divergent step size. Count the ones across chains or steps to diagnose a step that is too large.

Intuition Exact dynamics conserve energy, so a large energy jump means the integrator has stopped following the physics. An error of 10 corresponds to an acceptance probability of about \(e^{-10}\) when energy rises, and a step producing such errors is unstable, not merely inefficient.

Use it when you tune the HMC step size or check a run’s trustworthiness: a nonzero share of divergences means the step is too large for some region of the target, and the samples there are unreliable. Shrink step until the count is zero. The threshold 10 is fixed; compute energyOf differences yourself for a different cutoff.

sampler-validationstep-size-tuning

def diverged(x, p, step, mu, sigma)     ≔ abs(hmcRatio(x, p, step, mu, sigma)) > 10

Building blocks

kickMomentum(p, x, step, mu, sigma) — Let \(p, x \in \mathbb{R}\), \(h > 0\) (the argument step), \(\mu \in \mathbb{R}\) and \(\sigma > 0\). Then \(\operatorname{kickMomentum}(p, x, h, \mu, \sigma) = p + \tfrac{h}{2} g(x)\). A half-step momentum update along the gradient.
driftPosition(x, p, step) — Let \(x, p \in \mathbb{R}\) and \(h > 0\). Then \(\operatorname{driftPosition}(x, p, h) = x + h p\). A full-step position update along the momentum.
leapX(x, p, step, mu, sigma) — Let \(x\), \(p\), \(h\), \(\mu\), \(\sigma\) be as above. Then \(\operatorname{leapX}(x, p, h, \mu, \sigma) = x' = x + h\big(p + \tfrac{h}{2} g(x)\big)\). The position after one leapfrog step.
leapP(x, p, step, mu, sigma) — Let \(x\), \(p\), \(h\), \(\mu\), \(\sigma\) be as above. Then \(\operatorname{leapP}(x, p, h, \mu, \sigma) = p' = p + \tfrac{h}{2} g(x) + \tfrac{h}{2} g(x')\). The momentum after one leapfrog step.
hmcRatio(x, p, step, mu, sigma) — Let \(x\), \(p\), \(h\), \(\mu\), \(\sigma\) be as above. Then \(\operatorname{hmcRatio}(x, p, h, \mu, \sigma) = H(x, p) - H(x', p')\). The energy before the leapfrog step minus the energy after it: the log acceptance ratio.

Slice sampling

Slice sampling has no acceptance rate and no step size to tune. Draw a height uniformly under the density at the current point, then draw uniformly from the set of points where the density is above that height. For a Gaussian that set is an interval with a closed form,

\[\mu \pm \sigma \sqrt{-2 \log y}\]

where \(\log y = \log \pi(x) + \log u\), so no stepping-out procedure is needed.

sliceNext

sliceNext(x, mu, sigma, u, v) · x current state · mu mean · sigma standard deviation · u uniform draw for the height · v uniform draw along the slice → next state

Let \(x \in \mathbb{R}\) be the current state, \(\mu \in \mathbb{R}\) and \(\sigma > 0\), and let \(u, v \in [0, 1)\) be independent uniform draws. Let \(y = \ell(x) + \ln \operatorname{nonzero}(u)\) and \(w = \sigma \sqrt{\max(-2y, 0)}\). Then

\[\operatorname{sliceNext}(x, \mu, \sigma, u, v) = \mu - w + 2 w v .\]

One slice-sampling step. Every step moves; there is nothing to reject.

Intuition Sampling from a density is the same as sampling a point uniformly from the area under its curve. Alternate: pick a random height below the curve at the current point, then pick uniformly among all points where the curve is higher than that. The slice adapts its own width to the target, so there is no step size.

Use it when you want a sampler with nothing to tune and no wasted proposals. Here the slice is exact because the target is Gaussian; for general densities the slice must be found by stepping out and shrinking, which this definition does not do. Successive states are still correlated through the height draw, so thin or run several chains as for any MCMC.

bayesian-posterior-samplingsampler-validation

def sliceNext(x, mu, sigma, u, v) ≔ mu − sliceHalfWidth(sliceHeight(x, mu, sigma, u), sigma) + 2 · sliceHalfWidth(sliceHeight(x, mu, sigma, u), sigma) · v

Building blocks

sliceHeight(x, mu, sigma, u) — Let \(x \in \mathbb{R}\), \(\mu \in \mathbb{R}\), \(\sigma > 0\) and \(u \in [0, 1)\) a uniform draw. Then \(\operatorname{sliceHeight}(x, \mu, \sigma, u) = \ell(x) + \ln \operatorname{nonzero}(u)\). The log height of the slice, drawn under the density at x.
sliceHalfWidth(h, sigma) — Let \(y \in \mathbb{R}\) (the argument h) be a log height and \(\sigma > 0\). Then \(\operatorname{sliceHalfWidth}(y, \sigma) = \sigma \sqrt{\max(-2y, 0)}\). Half the width of the slice at log height h.

Gibbs sampling

A Gibbs move samples one coordinate exactly from its conditional given the others, so every draw is accepted. For a standard bivariate normal with correlation \(\rho\),

\[x \mid y \sim N(\rho y, \; 1 - \rho^2)\]

A sweep draws the first coordinate given the second, then the second given the new first.

gibbsDraw

gibbsDraw(other, rho, z) · other current value of the other coordinate · rho correlation in (−1, 1) scalar · z standard normal draw → new value of this coordinate

Let \(y \in \mathbb{R}\) (the argument other) be the current value of the other coordinate, let \(-1 < \rho < 1\) and let \(z \in \mathbb{R}\) be a standard normal draw. Then

\[\operatorname{gibbsDraw}(y, \rho, z) = \rho y + \sqrt{\max(1 - \rho^2, 0)} \; z .\]

One Gibbs update of a standard bivariate normal: \(\rho \cdot \text{other} + \sqrt{1 - \rho^2}\, z\).

Intuition A hard joint distribution broken into easy one-dimensional pieces: freeze one coordinate, draw the other exactly from its conditional, then swap roles. The chain moves only parallel to the axes, so when the variables are strongly correlated it zig-zags in tiny steps along a narrow diagonal ridge.

Use it when each conditional has a closed form — conjugate hierarchical models, latent-variable models — or to demonstrate how correlation slows mixing (compare \(\rho = 0.5\) with \(\rho = 0.99\)). This definition is only for zero means and unit variances; for other bivariate normals standardise first. When correlation is high, reparametrise or use a joint move such as hmcNext.

bayesian-posterior-samplingmultivariate-sampling

def gibbsDraw(other, rho, z) ≔ gibbsMean(other, rho) + gibbsSpread(rho) · z

Building blocks

gibbsMean(other, rho) — Let \(y \in \mathbb{R}\) (the argument other) and \(-1 < \rho < 1\). Then \(\operatorname{gibbsMean}(y, \rho) = \rho y\). The conditional mean, rho times the other coordinate.
gibbsSpread(rho) — Let \(-1 < \rho < 1\). Then \(\operatorname{gibbsSpread}(\rho) = \sqrt{\max(1 - \rho^2, 0)}\). The conditional standard deviation, the square root of 1 − rho squared.

Tempering and annealing

Raising the temperature \(T\) flattens a target to \(\pi(x)^{1/T}\), so a hot chain wanders freely between modes that a cold chain cannot cross. Parallel tempering runs chains at several temperatures and proposes swapping their states; the swap has its own Metropolis ratio, which keeps every chain correct:

\[\log \alpha = \Big(\frac{1}{T_c} - \frac{1}{T_h}\Big)\big(\log \pi(x_h) - \log \pi(x_c)\big)\]

Simulated annealing is tempering with a falling temperature and no swap: always accept a move to a lower cost, accept a move to a higher cost sometimes, and less often as it cools. A schedule that cools too fast freezes in the first basin it finds.

tempered

tempered(logp, temperature) · logp log density · temperature T → log density

Let \(\lambda \in \mathbb{R}\) (the argument logp) be a log density value and \(T > 0\) the temperature. Then

\[\operatorname{tempered}(\lambda, T) = \frac{\lambda}{\operatorname{nonzero}(T)} .\]

The tempered log density \(\log \pi / T\).

Intuition Heat melts the landscape: dividing the log density by \(T > 1\) shrinks every height difference, so peaks get lower and valleys shallower and a chain can walk between modes. \(T < 1\) sharpens the peaks instead, and \(T \to 0\) concentrates everything on the maximum.

Use it when you build a ladder of chains for parallel tempering, anneal toward a mode, or flatten a posterior during warm-up. For a Gaussian, temperature \(T\) is the same as scaling \(\sigma\) by \(\sqrt{T}\). The swap test between two tempered chains is swapRatio.

multimodal-samplingglobal-optimisation

def tempered(logp, temperature) ≔ logp / nonzero(temperature)

swapRatio

swapRatio(xCold, xHot, mu, sigma, cold, hot) · xCold state of the cold chain · xHot state of the hot chain · mu mean · sigma standard deviation · cold cold temperature · hot hot temperature → log acceptance ratio

Let \(x_c, x_h \in \mathbb{R}\) (the arguments xCold, xHot) be the states of the two chains, \(\mu \in \mathbb{R}\) and \(\sigma > 0\), and let \(T_c, T_h > 0\) (the arguments cold, hot) be their temperatures. Then

\[\operatorname{swapRatio}(x_c, x_h, \mu, \sigma, T_c, T_h) = \Big(\frac{\ell(x_h)}{T_c} + \frac{\ell(x_c)}{T_h}\Big) - \Big(\frac{\ell(x_c)}{T_c} + \frac{\ell(x_h)}{T_h}\Big) = \Big(\frac{1}{T_c} - \frac{1}{T_h}\Big)\big(\ell(x_h) - \ell(x_c)\big),\]

with each temperature passed through \(\operatorname{nonzero}\).

The log acceptance ratio for exchanging the two chains’ states. Accept the swap with keeps(swapRatio(…), u).

Intuition The hot chain scouts, the cold chain exploits. A swap hands a good state found by the hot chain down to the cold one, and it is accepted readily when the hot chain’s state is more probable than the cold chain’s. Because the swap is itself a Metropolis move on the joint system, every chain keeps its own correct distribution.

Use it when the target has several separated modes and a single chain gets stuck in one. Space temperatures so neighbouring chains accept swaps a reasonable share of the time (roughly 20–40%); if swaps almost never succeed, add intermediate temperatures. Only samples of the chain at \(T = 1\) belong to the target.

multimodal-samplingbayesian-posterior-sampling

def swapRatio(xCold, xHot, mu, sigma, cold, hot) ≔ (tempered(logTarget(xHot, mu, sigma), cold) + tempered(logTarget(xCold, mu, sigma), hot)) − (tempered(logTarget(xCold, mu, sigma), cold) + tempered(logTarget(xHot, mu, sigma), hot))

coolingAt

coolingAt(start, rate, k) · start initial temperature · rate cooling rate · k step number → temperature

Let \(T_0 > 0\) (the argument start) be the initial temperature, \(r \ge 0\) (the argument rate) the cooling rate and \(k \ge 0\) the step number. Then

\[\operatorname{coolingAt}(T_0, r, k) = T_k = T_0 \, e^{-r k} .\]

The exponential cooling schedule \(T_k = T_0 \, e^{-\text{rate}\cdot k}\).

Intuition Start hot enough to roam freely, then cool by the same fraction every step, so the search spends as many steps going from \(T_0\) to \(T_0/10\) as from \(T_0/10\) to \(T_0/100\). It is the same shape as decayed in Policies, written with the rate in the exponent.

Use it when you schedule simulated annealing (feed it to annealAccepts). Choose \(T_0\) comparable to typical cost differences so most uphill moves are accepted at first, and a rate so the temperature is small relative to cost differences by the last step, about \(\ln(T_0/T_{\text{end}})/k_{\max}\). Cooling too fast freezes in the first basin found; the logarithmic schedules with convergence guarantees are far slower.

global-optimisationhyperparameter-choice

def coolingAt(start, rate, k) ≔ start · exp(0 − rate · k)

annealAccepts

annealAccepts(here, there, temperature, u) · here cost of the current state · there cost of the proposed state · temperature current temperature · u uniform draw → mask

Let \(c, c' \in \mathbb{R}\) (the arguments here, there) be the costs of the current and proposed states, \(T > 0\) the temperature and \(u \in [0, 1)\) a uniform draw. Then

\[\operatorname{annealAccepts}(c, c', T, u) = \Big[\, u < \exp\Big(\min\Big(\frac{c - c'}{\operatorname{nonzero}(T)}, 0\Big)\Big) \,\Big] .\]

The annealing acceptance test for a cost to be minimised: 1 when the proposal is no worse, otherwise 1 with probability \(e^{-(\text{there} - \text{here})/T}\).

Intuition Greedy descent that is allowed to climb out of a hole. At high temperature a step uphill in cost is nearly free, so the search crosses ridges; as the temperature falls it becomes greedy and settles into a deep basin. Temperature is measured in the units of the cost.

Use it when you minimise a cost with many local minima and no useful gradient — scheduling, routing, layout, combinatorial assignment. Propose a neighbour in C#, test it with this mask at coolingAt(start, rate, k), and keep the best state ever seen as well as the current one. It does not sample a distribution; for that use keeps with a log-density ratio.

global-optimisationcombinatorial-optimisation

def annealAccepts(here, there, temperature, u) ≔ u < exp(min((here − there) / nonzero(temperature), 0))

Rejection and importance sampling

Rejection sampling is exact and wasteful: it keeps a draw when a uniform falls under the ratio of the target to an envelope, so the survivors come exactly from the target and the price is everything thrown away. Importance sampling is approximate and keeps everything: each draw is weighted by \(\pi(x)/q(x)\) instead. An importance estimate is only as good as its weights are even, and the effective sample size

\[\mathrm{ESS} = \frac{\left(\sum_i w_i\right)^2}{\sum_i w_i^2}\]

says how many equally weighted draws it is worth: \(n\) when the weights are equal, and 1 when a single draw carries everything.

underEnvelope

underEnvelope(x, mu, sigma, bound, u) · x proposed points vector · mu mean · sigma standard deviation · bound envelope height scalar · u uniform draws vector → mask vector

Let \(x \in \mathbb{R}^n\) be proposed points, \(\mu \in \mathbb{R}\), \(\sigma > 0\), let \(B \ge 1\) (the argument bound) be the envelope height and \(u \in [0, 1)^n\) independent uniform draws. Then

\[\operatorname{underEnvelope}(x, \mu, \sigma, B, u)_i = \big[\, u_i B < e^{\ell(x_i)} \,\big] .\]

Rejection sampling under a flat envelope of height bound: 1 for each proposal to keep, when \(u \cdot \text{bound} < \pi(x)\). The target is the unnormalised Gaussian whose peak is 1, so the envelope needs bound of at least 1; draw x uniformly over the region of interest.

Intuition Throw darts uniformly at a box that contains the curve and keep those that land under it. The survivors’ positions are distributed exactly like the density, and the empty space in the box is the price.

Use it when you need independent, exact draws in low dimension, or a ground truth to test an MCMC sampler against. Keep bound as close to 1 as possible and the region only as wide as the target’s support you care about: mass outside the region is never proposed, so a region of \(\mu \pm 5\sigma\) loses almost nothing. In more than a few dimensions nearly every dart misses; use MCMC or importanceWeight there. acceptedShare measures the waste.

random-variate-generationsampler-validation

def underEnvelope(x, mu, sigma, bound, u) ≔ u · bound < exp(logTarget(x, mu, sigma))

acceptedShare

acceptedShare(x, mu, sigma, bound, u) · arguments as for underEnvelope → scalar

Let \(x\), \(\mu\), \(\sigma\), \(B\) and \(u \in [0, 1)^n\) be as for underEnvelope. Then

\[\operatorname{acceptedShare}(x, \mu, \sigma, B, u) = \frac{1}{n} \sum_{i=1}^{n} \big[\, u_i B < e^{\ell(x_i)} \,\big] .\]

The fraction of proposals accepted, a measure of how tight the envelope is.

Intuition The share of darts under the curve equals the area under the curve divided by the area of the box. So it is both an efficiency score and a Monte Carlo estimate of an integral.

Use it when you tune the envelope (the share falls as bound or the region grows, so keep both as small as correctness allows), budget how many proposals are needed for a given number of samples (divide by the share), or estimate an unnormalised density’s integral over a region of length \(L\) as share × bound × \(L\) (hit-or-miss integration).

sampler-validationmonte-carlo-estimationnormalising-constants

def acceptedShare(x, mu, sigma, bound, u) ≔ einsum("i->", underEnvelope(x, mu, sigma, bound, u)) / count(u)

importanceWeight

importanceWeight(x, mu, sigma, proposalSpread) · x draws from the proposal · mu target mean · sigma target standard deviation · proposalSpread proposal standard deviation → weights

Let \(x \in \mathbb{R}\) be a draw from the proposal, \(\mu \in \mathbb{R}\) and \(\sigma > 0\) the target’s mean and standard deviation, and \(\tau > 0\) (the argument proposalSpread) the proposal’s standard deviation. Then

\[\operatorname{importanceWeight}(x, \mu, \sigma, \tau) = \exp\big( \operatorname{logTarget}(x, \mu, \sigma) - \operatorname{logTarget}(x, 0, \tau) \big) = \exp\Big( -\frac{(x - \mu)^2}{2\sigma^2} + \frac{x^2}{2\tau^2} \Big) .\]

Importance weights for draws from the proposal \(N(0, \text{proposalSpread}^2)\) (centred at zero) towards the target \(N(\mu, \sigma^2)\). Both densities are unnormalised, so the weights are correct up to one common factor; that factor cancels in weightedMeanOf and effectiveSize.

Intuition You sampled from the wrong distribution, so correct for it: each draw counts in proportion to how much more (or less) often the target would have produced it than the proposal did. Draws the proposal over-produced are down-weighted, under-produced ones up-weighted.

Use it when sampling the target is hard but a simple proposal is easy — estimating expectations, reweighting simulation results to new parameters, or off-policy evaluation. The proposal here is fixed at mean 0; make it wider than the target: with \(\tau\) below \(\sigma/\sqrt{2}\) the weights have infinite variance and a few huge ones dominate. Always check effectiveSize. When you need unweighted samples, use underEnvelope or MCMC.

importance-weightingmonte-carlo-estimation

def importanceWeight(x, mu, sigma, proposalSpread) ≔ exp(logTarget(x, mu, sigma) − logTarget(x, 0, proposalSpread))

weightedMeanOf

weightedMeanOf(x, w) · x values vector · w weights vector → scalar

Let \(x \in \mathbb{R}^n\) be values and \(w \in [0, \infty)^n\) weights. Then

\[\operatorname{weightedMeanOf}(x, w) = \frac{\sum_{i=1}^{n} x_i w_i}{\operatorname{nonzero}\big(\sum_{i=1}^{n} w_i\big)} .\]

The self-normalised weighted mean \(\sum_i x_i w_i / \sum_i w_i\).

Intuition A centre of mass: each value pulls on the average with its own weight. Dividing by the total weight means the weights need only be right up to a common factor, which is all unnormalised densities can give.

Use it when you turn importance weights into an estimate of an expectation, average particle states in a filter, or compute any weighted average. Self-normalising adds a small bias that vanishes as the sample grows, in exchange for not needing normalising constants. Pass \(f(x)\) rather than \(x\) to estimate \(\mathbb{E}[f]\). Its reliability is only as good as effectiveSize.

importance-weightingmonte-carlo-estimation

def weightedMeanOf(x, w) ≔ einsum("i,i->", x, w) / nonzero(einsum("i->", w))

effectiveSize

effectiveSize(w) · w weights vector → scalar

Let \(w \in [0, \infty)^n\) be weights. Then

\[\operatorname{effectiveSize}(w) = \frac{\big(\sum_{i=1}^{n} w_i\big)^2}{\operatorname{nonzero}\big(\sum_{i=1}^{n} w_i^2\big)} .\]

The effective sample size of a set of weights. It does not depend on their scale.

Intuition How many equally weighted draws your weighted sample is really worth. A thousand draws where three carry nearly all the weight behave like about three: the estimate has the noise of three samples, whatever the count says.

Use it when you judge an importance-sampling estimate, choose a proposal (larger ESS is better), or decide when a particle filter must resample (degenerate is a fixed-threshold version). Report it next to any weighted estimate. A high ESS does not prove the proposal covers the target: regions the proposal never reaches leave no trace in the weights.

sampler-validationimportance-weighting

def effectiveSize(w) ≔ einsum("i->", w) · einsum("i->", w) / nonzero(einsum("i,i->", w, w))

degenerate

degenerate(w) · w weights vector → 1 or 0

Let \(w \in [0, \infty)^n\) be weights. Then

\[\operatorname{degenerate}(w) = \Big[\, \operatorname{effectiveSize}(w) < \frac{n}{10} \,\Big] .\]

1 when the effective sample size has fallen below a tenth of the number of draws: a weighted estimate that is effectively a handful of samples and should not be trusted (or a particle cloud that needs resampling).

Intuition An alarm on weight collapse: it fires when fewer than one in ten draws is effectively contributing. Past that point the estimate is driven by a handful of lucky samples.

Use it when you want a yes/no trigger — resample a particle cloud only when it fires, reject a proposal distribution, or flag an importance estimate as unreliable. The threshold \(n/10\) is fixed; many particle filters resample earlier, at \(\mathrm{ESS} < n/2\), which you can test directly with effectiveSize.

sampler-validationstate-estimation

def degenerate(w)    ≔ effectiveSize(w) < count(w) / 10

Monte Carlo estimators

Two classical ways to lower the variance of a Monte Carlo average. Antithetic variates pair each uniform draw with its mirror \(1 - u\), so the errors of a monotone integrand cancel instead of accumulating. Control variates subtract a correlated quantity \(c\) whose mean is known, \(y - \beta (c - \mathbb{E}[c])\). reduction reports the variance ratio actually achieved: a value at or above one means the pairing or control was badly chosen.

antithetic

antithetic(u) · u uniform draws → mirrored draws

Let \(u \in [0, 1)\) be a uniform draw. Then

\[\operatorname{antithetic}(u) = 1 - u .\]

The antithetic partner \(1 - u\) of each uniform draw.

Intuition Balance every draw with its mirror image. If a draw near 0 gives a result that is too low, its partner near 1 gives one that is too high, and their average lands close to the truth, so pairs disagree less than independent draws would.

Use it when you estimate an expectation of a monotone function of a uniform — option prices, inverse-CDF simulations, anything that increases with the draw. Evaluate \(f(u)\) and \(f(1-u)\), average each pair, and check with reduction against an equal number of independent evaluations. For a symmetric function such as \((u - 0.5)^2\) the pair gives identical values and the trick is worse than useless.

variance-reductionmonte-carlo-estimation

def antithetic(u)         ≔ 1 − u

controlled

controlled(y, c, cmean, beta) · y samples of the quantity · c samples of the control · cmean known mean of the control · beta coefficient → adjusted samples

Let \(y, c \in \mathbb{R}\) be paired samples of the quantity and the control, \(\bar{c} = \mathbb{E}[c]\) (the argument cmean) the control’s known mean and \(\beta \in \mathbb{R}\) the coefficient. Then

\[\operatorname{controlled}(y, c, \bar{c}, \beta) = y - \beta (c - \bar{c}) .\]

The control-variate adjusted samples \(y - \beta (c - \mathbb{E}[c])\). Their mean estimates the mean of y.

Intuition If a quantity you can simulate moves together with one whose mean you already know, the known one tells you how lucky this particular sample was. Subtract that luck, scaled by \(\beta\), and the noise shared by the two cancels while the mean of y is untouched.

Use it when you have a cheap, correlated quantity with a known expectation — the analytic price of a simpler option, the input’s own mean, a linearised model. The best \(\beta\) is \(\operatorname{Cov}(y, c)/\operatorname{Var}(c)\), which cuts the variance by a factor \(1 - \rho^2\); estimating it from the same samples adds a small bias. A weakly correlated control helps little, so confirm with reduction.

variance-reductionmonte-carlo-estimation

def controlled(y, c, cmean, beta) ≔ y − beta · (c − cmean)

varianceOf

varianceOf(v) · v samples vector → scalar

Let \(v \in \mathbb{R}^n\) be samples and \(\bar{v} = \frac{1}{n} \sum_i v_i\). Then

\[\operatorname{varianceOf}(v) = \frac{\sum_{i=1}^{n} (v_i - \bar{v})^2}{\operatorname{nonzero}(n - 1)} .\]

The sample variance, with the \(n - 1\) divisor.

Intuition The average squared distance from the sample mean, corrected for the fact that the mean was fitted to these same points and so sits closer to them than the true mean does. Dividing by \(n - 1\) instead of \(n\) makes the estimate unbiased.

Use it when you need the spread of Monte Carlo samples: the standard error of their mean is \(\sqrt{\text{varianceOf}(v)/n}\), which tells you how many more samples a target accuracy needs. It assumes independent samples; for correlated MCMC output it understates the error of the mean. To compare two estimators use reduction.

uncertainty-quantificationmonte-carlo-estimation

def varianceOf(v)         ≔ einsum("i,i->", v − mean(v), v − mean(v)) / nonzero(count(v) − 1)

reduction

reduction(plain, improved) · plain samples from the plain estimator vector · improved samples from the variance-reduced estimator vector → scalar

Let \(a \in \mathbb{R}^n\) (the argument plain) and \(b \in \mathbb{R}^{n'}\) (the argument improved) be samples from the two estimators. Then

\[\operatorname{reduction}(a, b) = \frac{\operatorname{varianceOf}(b)}{\operatorname{nonzero}\big(\operatorname{varianceOf}(a)\big)} .\]

The ratio of sample variances, improved over plain. Below one is an improvement.

Intuition How many times fewer samples the improved estimator needs for the same accuracy: a ratio of 0.1 means one tenth. It measures what a variance-reduction trick actually achieved, not what it promised.

Use it when you validate antithetic pairs, controlled samples or a new proposal against the plain estimator. Compare at equal cost: an antithetic pair uses two function evaluations, so set its variance against plain samples that are averages of two independent draws. With few samples the ratio is itself noisy.

variance-reductionestimator-comparison

def reduction(plain, improved) ≔ varianceOf(improved) / nonzero(varianceOf(plain))

The bootstrap particle filter

A particle filter tracks a state through a model that need not be linear or Gaussian, where the Kalman filter cannot go (on a linear-Gaussian model the two should agree). Each step weights the particles by how well they explain the observation, accumulates the evidence from the average weight, and resamples. Resampling here is systematic: one uniform per step lays a comb of \(n\) evenly spaced teeth \((k + u)/n\) over the cumulative weights, which has lower variance than drawing \(n\) independent indices.

particleWeights

particleWeights(cloud, observed, noise) · cloud particle states vector · observed observation scalar · noise observation noise standard deviation scalar → weights vector

Let \(x \in \mathbb{R}^n\) (the argument cloud) be the particle states, \(y \in \mathbb{R}\) (the argument observed) the observation and \(\sigma > 0\) (the argument noise) the observation noise. Then

\[\operatorname{particleWeights}(x, y, \sigma)_i = \exp\Big( -\frac{(x_i - y)^2}{\operatorname{nonzero}(2\sigma^2)} \Big) .\]

The unnormalised Gaussian likelihood \(\exp(-(x_i - y)^2 / 2\sigma^2)\) of the observation under each particle.

Intuition Each particle is a guess of the hidden state; its weight is how well that guess explains what the sensor just reported. Guesses near the observation score close to 1, guesses many noise-widths away almost 0.

Use it when you track a hidden quantity from noisy measurements — position, level, a drifting parameter — with a particle filter. The observation model is fixed: the sensor reads the state directly plus Gaussian noise of standard deviation noise. Propagate the particles through your dynamics in C#, weight them here, then feed the weights to logEvidence and resample. Too small a noise collapses the weights onto one particle; check effectiveSize. For linear-Gaussian models a Kalman filter is exact and cheaper.

state-estimationnonlinear-filtering

def particleWeights(cloud, observed, noise) ≔ exp(logTarget(cloud, observed, noise))

drawnAgain

drawnAgain(cloud, w, u) · cloud particle states vector n · w weights vector n · u one uniform draw scalar → resampled states vector n

Let \(x \in \mathbb{R}^n\) (the argument cloud) be the particle states, \(w \in [0, \infty)^n\) the weights and \(u \in [0, 1)\) one uniform draw; with positions \(j = 0, 1, \dots, n-1\) let \(\bar{w} = \operatorname{normalised}(w)\), \(W_j = \sum_{j' \le j} \bar{w}_{j'}\) and \(c_j = (j + u)/n\). Then, as the code computes it, for every \(i\)

\[\operatorname{drawnAgain}(x, w, u)_i = \sum_{j=0}^{n-1} \big[\, W_j \ge c_j \,\big] \, \big[\, W_j - \bar{w}_j < c_j \,\big] \, x_j ,\]

which does not depend on \(i\). Systematic resampling would instead set entry \(k\) to \(x_j\) for the one \(j\) with \(W_j - \bar{w}_j < c_k \le W_j\): the code compares tooth \(j\) of the comb with particle \(j\) rather than every tooth with every particle.

Known issue: in this release the result is not a resampled cloud. Every entry comes back equal to the same sum of particles (for example, particles [10, 20, 30, 40] with equal weights and u = 0.5 return [100, 100, 100, 100], where [10, 20, 30, 40] is expected). Do not rely on drawnAgain until this is corrected; resample in C# from particleWeights instead.

Intuition Resampling kills particles that explain the data poorly and copies those that explain it well, so effort is not wasted tracking hopeless guesses. Systematic resampling lays \(n\) evenly spaced teeth over the cumulative weights with one random offset: a particle with weight \(\bar w\) receives almost exactly \(n \bar w\) copies, with less randomness than \(n\) independent draws.

Use it when — not in this release: as the known issue above states, it does not return a resampled cloud. Implement systematic resampling in C# from normalised particleWeights, and resample only when degenerate (or your own ESS threshold) fires, since resampling itself adds noise.

state-estimationparticle-degeneracy

def drawnAgain(cloud, w, u) ≔ einsum("ij,j->i", (runningTotal(normalised(w)) ≥ laid(comb(w, u))) · (runningTotal(normalised(w)) − normalised(w) < laid(comb(w, u))), cloud)

logEvidence

logEvidence(w) · w unnormalised weights of one step vector → scalar

Let \(w \in [0, \infty)^n\) be the unnormalised weights of one step. Then

\[\operatorname{logEvidence}(w) = \ln \operatorname{nonzero}\Big( \frac{1}{n} \sum_{i=1}^{n} w_i \Big) .\]

The log of the average weight, this step’s contribution to the log marginal likelihood of the observations; sum it over steps. With weights from particleWeights, which omit the Gaussian’s normalising constant, add \(-\log(\sigma\sqrt{2\pi})\) per observation for the true value.

Intuition The average particle weight estimates how probable this observation was given everything seen before. Multiplying those across steps — adding their logs — gives how probable the whole observation sequence is under the model: a score for the model itself, not for any one state.

Use it when you compare models or parameter settings (noise levels, dynamics) on the same data — higher is better — or drive parameter estimation such as particle MCMC. The plain average assumes the particles carried equal weight before this step, that is, you resampled; otherwise use the weighted average. Add the normalising constant when comparing different noise values, since it changes with \(\sigma\).

model-comparisonstate-estimation

def logEvidence(w) ≔ ln(nonzero(einsum("i->", w) / count(w)))

Building blocks

normalised(w) — Let \(w \in [0, \infty)^n\). Then \(\operatorname{normalised}(w)_i = w_i / \operatorname{nonzero}\big(\sum_{j} w_j\big)\). The weights divided by their sum.
comb(w, u) — Let \(w \in \mathbb{R}^n\) and \(u \in [0, 1)\), with positions \(k = 0, 1, \dots, n-1\). Then \(\operatorname{comb}(w, u)_k = (k + u)/n\). The systematic resampling comb, positions (k + u)/n for k = 0 … n−1.

Stick breaking and nonparametric priors

A Dirichlet process is a stick broken forever: break off a Beta(1, α) fraction of what remains, again and again. The pieces are the mixture weights,

\[\pi_k = v_k \prod_{j<k} (1 - v_j), \qquad v_k \sim \mathrm{Beta}(1, \alpha)\]

and a small concentration α breaks off large pieces early, so fewer clusters are effectively in play. Beta(1, α) has a closed-form inverse CDF, \(v = 1 - (1 - u)^{1/\alpha}\), so each break comes exactly from a uniform draw, with no rejection.

The expected number of clusters grows like \(\alpha \log n\) under the Chinese restaurant process and like \(n^d\) under the Pitman–Yor process with discount \(d\) (the power law that real vocabularies follow). The formulas below are the standard large-\(n\) approximations, not exact expectations.

In this section \(\alpha > 0\) is the argument concentration, \(n \ge 0\) the argument customers, and positions of a vector of \(K\) sticks are counted \(k = 0, 1, \dots, K-1\).

breakAt

breakAt(u, concentration) · u uniform draws · concentration α > 0 → break fractions

Let \(u \in [0, 1)\) be a uniform draw and \(\alpha > 0\). Then

\[\operatorname{breakAt}(u, \alpha) = 1 - (1 - u)^{1/\operatorname{nonzero}(\alpha)} .\]

Converts uniform draws into Beta(1, α) draws by the inverse CDF.

Intuition How much of the remaining stick to snap off. The average fraction is \(1/(1+\alpha)\): with \(\alpha\) small most of the stick goes at once (a few big clusters), with \(\alpha\) large each break is a sliver (many small clusters).

Use it when you need Beta(1, α) variates exactly from uniforms — the building block of stickWeights, or any model with such a prior. Because it is a monotone transform of \(u\), seeded draws give reproducible breaks and pair well with antithetic. For a vector of weights call stickWeights rather than chaining breaks by hand.

random-variate-generationmixture-modelling

def breakAt(u, concentration) ≔ 1 − (1 − u) ^ (1 / nonzero(concentration))

stickWeights

stickWeights(u, concentration) · u uniform draws, one per stick vector K · concentration α → weights vector K

Let \(u \in [0, 1)^K\) be independent uniform draws, one per stick, and \(\alpha > 0\), and let \(v_k = \operatorname{breakAt}(u_k, \alpha)\). Then

\[\operatorname{stickWeights}(u, \alpha)_k = v_k \prod_{j<k} (1 - v_j), \qquad k = 0, 1, \dots, K-1,\]

where the product is computed as \(\exp\big(\sum_{j \le k} \ln \operatorname{nonzero}(1 - v_j) - \ln \operatorname{nonzero}(1 - v_k)\big)\).

The first \(K\) stick-breaking weights of a Dirichlet process, from \(K\) uniform draws.

Intuition Mixture weights for an unbounded number of clusters: take a fraction of a unit stick for cluster 0, a fraction of what is left for cluster 1, and so on. Early clusters tend to be large and later ones ever smaller, so a finite number carry almost all the mass without fixing that number in advance.

Use it when you simulate from or initialise a Dirichlet-process mixture, a truncated variational approximation, or any model where the number of components should be learned. Choose \(K\) large enough that the leftover mass is negligible (check with stickMass). The weights come out stochastically decreasing, not sorted, and \(K\) weights sum to slightly less than one.

mixture-modellingunknown-cluster-countrandom-variate-generation

def stickWeights(u, concentration) ≔ breakAt(u, concentration) · unbroken(breakAt(u, concentration))

stickMass

stickMass(u, concentration) · u uniform draws vector K · concentration α → scalar

Let \(u \in [0, 1)^K\) be independent uniform draws and \(\alpha > 0\). Then

\[\operatorname{stickMass}(u, \alpha) = \sum_{k=0}^{K-1} \operatorname{stickWeights}(u, \alpha)_k .\]

The total of the first \(K\) weights. One minus it is the mass still left on the stick, a check on whether \(K\) is large enough.

Intuition What the truncation has captured. Whatever remains, \(1 - \text{stickMass}\), is the probability of a cluster beyond the first \(K\) — the part of the infinite model you have thrown away.

Use it when you pick the truncation level of a stick-breaking prior: raise \(K\) until the leftover mass is below your tolerance for typical draws. On average the leftover is \((\alpha/(1+\alpha))^K\), so larger \(\alpha\) needs a larger \(K\). You can also put the leftover on a final “everything else” component so the weights sum to one.

truncation-errormixture-modelling

def stickMass(u, concentration)    ≔ einsum("i->", stickWeights(u, concentration))

crpTables

crpTables(concentration, customers) · concentration α · customers n → scalar

Let \(\alpha > 0\) and \(n \ge 0\). Then

\[\operatorname{crpTables}(\alpha, n) = \alpha \ln\Big(1 + \frac{n}{\operatorname{nonzero}(\alpha)}\Big) .\]

The approximate expected number of occupied tables in a Chinese restaurant process after \(n\) customers, \(\alpha \log(1 + n/\alpha)\).

Intuition Each new customer joins an occupied table in proportion to its size or opens a new one with weight \(\alpha\). New tables keep appearing, but ever more rarely, so their number grows only logarithmically: doubling the data adds about \(\alpha \ln 2\) clusters.

Use it when you choose the concentration \(\alpha\) of a Dirichlet-process prior: solve for the \(\alpha\) that gives the number of clusters you expect for your data size, or check whether a fitted model’s cluster count is plausible under its prior. It is an approximation, not the exact expectation. If clusters in your data keep appearing like a power law (words, species), use pitmanTables.

prior-elicitationunknown-cluster-count

def crpTables(concentration, customers) ≔ concentration · ln(1 + customers / nonzero(concentration))

pitmanTables

pitmanTables(concentration, discount, customers) · concentration α · discount d in (0, 1) · customers n → scalar

Let \(\alpha > 0\), \(0 < d < 1\) (the argument discount) and \(n \ge 0\). Then

\[\operatorname{pitmanTables}(\alpha, d, n) = \frac{\alpha \, n^{d}}{\operatorname{nonzero}(d)} .\]

The power-law growth of the number of tables in a Pitman–Yor process, \(\alpha n^d / d\).

Intuition The discount \(d\) takes a little from every occupied table and hands it to new ones, so new clusters keep arriving much faster than in the Chinese restaurant process: their number grows like \(n^d\) rather than \(\log n\). That heavy tail of rare clusters is what word frequencies and species counts look like.

Use it when you set the discount of a Pitman–Yor prior for data with a power-law vocabulary (language models, topic models, ecology): fit \(d\) to how the number of distinct items grows with sample size. Treat it as the growth rate; the constant in front is a rough scale. \(d\) close to 0 behaves like crpTables.

prior-elicitationpower-law-modellingunknown-cluster-count

def pitmanTables(concentration, discount, customers) ≔ concentration · customers ^ discount / nonzero(discount)

buffetFeatures

buffetFeatures(concentration, customers) · concentration α · customers n → scalar

Let \(\alpha > 0\) and \(n \ge 0\). Then

\[\operatorname{buffetFeatures}(\alpha, n) = \alpha \ln(1 + n) .\]

The approximate expected total number of features (dishes) in an Indian buffet process after \(n\) customers, \(\alpha \log(1 + n)\).

Intuition In the Indian buffet each object can have many features at once. Customer \(i\) tries popular dishes and about \(\alpha/i\) new ones, so the total number of distinct features grows like a harmonic sum — logarithmically in the number of objects.

Use it when you choose \(\alpha\) for a latent-feature model (for example overlapping topics, hidden causes, binary factor analysis) so the prior’s feature count fits your data size, or sanity-check a fitted model’s number of features. Together with buffetOwned it fixes both the total and the per-object number of features. For exclusive cluster membership use crpTables instead.

latent-feature-modellingprior-elicitation

def buffetFeatures(concentration, customers) ≔ concentration · ln(1 + customers)

buffetOwned

buffetOwned(concentration) · concentration α → scalar

Let \(\alpha > 0\). Then

\[\operatorname{buffetOwned}(\alpha) = \alpha .\]

The expected number of features each customer takes in an Indian buffet process, which is α.

Intuition The same \(\alpha\) that controls how many features exist overall is also how many each object has on average, however many objects there are. It is the density of each row of the binary feature matrix.

Use it when you set \(\alpha\) from domain knowledge about how many features one object typically has (“each document touches about three topics”), then read the implied total from buffetFeatures. It is a single parameter doing two jobs; if the two counts you expect disagree, the one-parameter buffet prior cannot fit both.

latent-feature-modellingprior-elicitation

def buffetOwned(concentration) ≔ concentration

Building blocks

unbroken(v) — Let \(v \in [0, 1)^K\), with positions \(k = 0, 1, \dots, K-1\). Then \(\operatorname{unbroken}(v)_k = \exp\big(\sum_{j \le k} \ln \operatorname{nonzero}(1 - v_j) - \ln \operatorname{nonzero}(1 - v_k)\big) = \prod_{j<k} (1 - v_j)\) (the last equality whenever every \(1 - v_j \ge 10^{-300}\)). The length of stick left before each break, the running product of 1 − v over earlier breaks.

Examples

All examples use one file, Chainsteps.hb, so the C# class is Chainsteps. Vectors are passed as Tensor so the general overload is called.

Advance two Metropolis chains by one step

Both chains start at 0 targeting \(N(0, 1)\) with step 1 and the same normal draw 0.5, so both propose 0.5, whose log ratio is \(-0.125\) (acceptance probability 0.8825). The first chain’s uniform draw 0.9 rejects; the second’s 0.5 accepts.

open Prelude.Sampling

fn metropolisStep(x, step, mu, sigma, z, u) ↦ hastingsNext(x, step, mu, sigma, z, u)

using MLambda.Hilbert.Runtime;
using Examples;

var x = (Tensor)new[] { 0d, 0d };
var z = (Tensor)new[] { 0.5, 0.5 };
var u = (Tensor)new[] { 0.9, 0.5 };

var next = Chainsteps.MetropolisStep(x, 1.0, 0.0, 1.0, z, u);

Output (verified):

hastingsNext = [0, 0.5]

Take gradient-guided steps with MALA and HMC

Target \(N(0, 1)\). MALA from 2 with step 0.5 and draw 0 proposes \(2 + 0.25 \cdot (-2) = 1.5\), and the corrected ratio accepts it. HMC from 1 with momentum 0.5 and step 0.1 kicks the momentum to 0.45 and drifts to 1.045, which is accepted.

open Prelude.Sampling

fn malaStep(x, step, mu, sigma, z, u) ↦ malaNext(x, step, mu, sigma, z, u)
fn hamiltonStep(x, p, step, mu, sigma, u) ↦ hmcNext(x, p, step, mu, sigma, u)

using Examples;

var mala = Chainsteps.MalaStep(2.0, 0.5, 0.0, 1.0, 0.0, 0.5);
var hmc  = Chainsteps.HamiltonStep(1.0, 0.5, 0.1, 0.0, 1.0, 0.5);

Output (verified):

malaNext = 1.5
hmcNext  = 1.045

Draw with a slice sampler and a Gibbs sweep

Slice: from 0 on \(N(0, 1)\) with height draw 0.5 the slice is \(\pm\sqrt{2 \ln 2} = \pm 1.17741\), and the draw 0.75 lands three quarters of the way across it. Gibbs: with correlation 0.6, the other coordinate at 1 and a normal draw of 1, the update is \(0.6 + 0.8 = 1.4\).

open Prelude.Sampling

fn sliceStep(x, mu, sigma, u, v) ↦ sliceNext(x, mu, sigma, u, v)
fn gibbsStep(other, rho, z) ↦ gibbsDraw(other, rho, z)

using Examples;

var slice = Chainsteps.SliceStep(0.0, 0.0, 1.0, 0.5, 0.75);
var gibbs = Chainsteps.GibbsStep(1.0, 0.6, 1.0);

Output (verified, rounded to 6 significant digits):

sliceNext = 0.588705
gibbsDraw = 1.4

Judge importance weights

Equal weights are worth all four draws; a single dominant weight is worth one. A weighted mean of 1, 2 and 3 with weights 1, 1 and 2 is \(9/4\).

open Prelude.Sampling

fn essOf(w) ↦ effectiveSize(w)
fn weightedAverage(x, w) ↦ weightedMeanOf(x, w)

using MLambda.Hilbert.Runtime;
using Examples;

var even = Chainsteps.EssOf((Tensor)new[] { 1d, 1d, 1d, 1d });
var one  = Chainsteps.EssOf((Tensor)new[] { 1d, 0d, 0d, 0d });
var avg  = Chainsteps.WeightedAverage((Tensor)new[] { 1d, 2d, 3d }, (Tensor)new[] { 1d, 1d, 2d });

Output (verified):

effectiveSize (even) = 4
effectiveSize (one)  = 1
weightedMeanOf       = 2.25

Break a stick for a Dirichlet process

With α = 1 a uniform draw of 0.5 breaks off exactly half of what remains. The expected number of tables for 100 customers at α = 1 is \(\ln 101\).

open Prelude.Sampling

fn sticks(u, concentration) ↦ stickWeights(u, concentration)
fn stickTotal(u, concentration) ↦ stickMass(u, concentration)
fn tablesExpected(concentration, customers) ↦ crpTables(concentration, customers)

using MLambda.Hilbert.Runtime;
using Examples;

var u = (Tensor)new[] { 0.5, 0.5, 0.5 };

var weights = Chainsteps.Sticks(u, 1.0);
var mass    = Chainsteps.StickTotal(u, 1.0);
var tables  = Chainsteps.TablesExpected(1.0, 100.0);

Output (verified, rounded to 6 significant digits):

stickWeights = [0.5, 0.25, 0.125]
stickMass    = 0.875
crpTables    = 4.61512

Compare the variance of two estimators

A plain estimator’s samples [1, 3, 5, 7] have variance 20/3; an improved estimator’s [3, 4, 4, 5] have variance 2/3, a tenfold reduction. The same file also gives a cooling schedule: starting at 10 with rate 0.5, the temperature after 2 steps is \(10/e\).

open Prelude.Sampling

fn varianceRatio(plain, improved) ↦ reduction(plain, improved)
fn temperatureAt(start, rate, k) ↦ coolingAt(start, rate, k)

using MLambda.Hilbert.Runtime;
using Examples;

var ratio = Chainsteps.VarianceRatio((Tensor)new[] { 1d, 3d, 5d, 7d }, (Tensor)new[] { 3d, 4d, 4d, 5d });
var temp  = Chainsteps.TemperatureAt(10.0, 0.5, 2.0);

Output (verified, rounded to 6 significant digits):

reduction = 0.1
coolingAt = 3.67879