Tutorials › Stochastic Differential Equations in Pharmacokinetics › Ancestor Sampling for Particle Gibbs

Particle Filters and Bayesian Estimation · Tutorial 72 of 100

Ancestor Sampling for Particle Gibbs

Learn how to calculate and interpret ancestor-sampling probabilities for a fixed reference trajectory, and why they can help Particle Gibbs update latent paths.

Advanced 13 min read

What You'll Learn

  • Write the latent-path smoothing target with all observation likelihood factors, including the time-zero factor.
  • Derive the ancestor-sampling distribution from particle weights and transition densities.
  • Calculate ancestor probabilities for Gaussian and discrete-state examples.
  • Explain how future-state compatibility can redirect a reference path to a different particle genealogy.
  • Diagnose degenerate ancestor weights and distinguish ancestor sampling from backward simulation.

Why a Reference Path Needs an Ancestor Choice

Particle Gibbs updates a latent trajectory while retaining a previously selected trajectory as a reference. In a conditional particle filter, the reference states are inserted into the particle population, but their connection to earlier particles need not be fixed. Ancestor sampling uses the next reference state to choose which earlier particle should be its parent. This gives a path update a chance to change the reference trajectory’s earlier history instead of keeping the same genealogy.

The key calculation is local: among particles at time \(t-1\), how plausible is each one as the parent of the fixed reference state at time \(t\)? The answer combines its current particle weight with the model’s transition density into that reference state. Observation likelihoods for the fixed future segment do not distinguish candidate parents once the reference state is fixed.

Target and Conditional Particle Filter

Let \(X_{0:T}=(X_0,\ldots,X_T)\) be a Markov latent process, with initial density \(\mu(x_0)\) and transition density \(f_t(x_t\mid x_{t-1})\). Conditional on the states, observations \(Y_{0:T}\) have likelihood contributions \(g_t(y_t\mid x_t)\). The smoothing target for a realized observation sequence \(y_{0:T}\) is

$$ \pi(x_{0:T}\mid y_{0:T}) \propto \mu(x_0)g_0(y_0\mid x_0) \prod_{t=1}^{T} f_t(x_t\mid x_{t-1})g_t(y_t\mid x_t). $$

The factor \(g_0(y_0\mid x_0)\) is included because an observation at time zero is part of the stated data. If the model has no observation at time zero, set that likelihood contribution to one. Known inputs, such as doses, may be included in the transition conditioning variables; they are held fixed in the ancestor calculation.

For the bootstrap particle filter, particles at time zero are drawn from \(\mu\), and their initial weights are proportional to \(g_0(y_0\mid X_0^i)\). At each later time, particles are propagated using \(f_t\) and weighted by \(g_t(y_t\mid X_t^i)\), with resampling as needed. Write \(W_{t-1}^i\) for the normalized weights of the particles \(X_{t-1}^i\), so that \(\sum_i W_{t-1}^i=1\).

Definition: In ancestor sampling, a fixed reference state \(x_t'\) at time \(t\) is assigned an ancestor index \(A_{t-1}\) among the particles at time \(t-1\). For the bootstrap particle filter, its conditional probabilities are $$ P(A_{t-1}=i\mid x_t',X_{t-1}^{1:N},W_{t-1}^{1:N}) = \frac{W_{t-1}^i f_t(x_t'\mid X_{t-1}^i)} {\sum_{j=1}^{N} W_{t-1}^j f_t(x_t'\mid X_{t-1}^j)}, $$ provided the denominator is positive and finite.

Deriving the Ancestor Probabilities

Fix a time \(t\), a reference suffix \(x_{t:T}'\), and the particle population and weights at \(t-1\). Consider candidate parent \(i\). Its contribution to the probability of this connection is its current weight \(W_{t-1}^i\), multiplied by the probability density of reaching the fixed next state, \(f_t(x_t'\mid X_{t-1}^i)\). Once \(x_t'\) is fixed, the observation likelihood \(g_t(y_t\mid x_t')\) and every later transition and observation factor in the suffix are the same for every candidate parent.

1
Start with the candidate’s filtering mass.
Particle \(i\) represents mass \(W_{t-1}^i\) at time \(t-1\).
2
Account for reaching the reference state.
The transition contributes \(f_t(x_t'\mid X_{t-1}^i)\). A particle poorly compatible with \(x_t'\) receives little ancestor mass even if its current weight is large.
3
Normalize across candidates.
Divide each product by their sum. The fixed suffix likelihood terms cancel because they do not depend on \(i\).

This derivation also shows why the observation likelihood at time \(t\) is not an extra candidate-specific factor in the ancestor formula: the reference state \(x_t'\) is the same for all candidate ancestors. The observation likelihood is important in the particle weights used at time \(t\), and the complete smoothing target still contains every observation factor.

Worked Example: Gaussian Transition Compatibility

Suppose the particles at time \(t-1\) are \(0,1,2\), with normalized weights \(0.2,0.5,0.3\). The transition is \(X_t\mid X_{t-1}=x\sim N(0.8x,0.25)\), and the fixed reference state is \(x_t'=1.2\). Calculate the ancestor probabilities.

The transition density has a common normalizing constant across the three candidates, so it suffices to use \(\exp\{-(1.2-0.8x)^2/(2(0.25))\}\). For \(x=0\), the residual is \(1.2\), giving \(\exp(-2.88)=0.0561\). For \(x=1\), the residual is \(0.4\), giving \(\exp(-0.32)=0.7261\). For \(x=2\), the residual is \(-0.4\), also giving \(0.7261\).

Multiplying by the particle weights gives unnormalized ancestor masses \(0.01122695\), \(0.36307452\), and \(0.21784471\). Their sum is \(0.59214618\). Therefore the probabilities are approximately \(0.01896,\ 0.61315,\ 0.36789\), respectively. The middle particle is most likely because it has both a substantial weight and a transition mean close to the reference state.

Two Consequences of the Formula

The first result gives the conditional distribution directly. The second makes the role of future-state compatibility especially clear by comparing two candidate ancestors.

Result 1 (conditional ancestor distribution): Suppose the particle weights are nonnegative and normalized, each transition density is defined at the reference state, and $$ 0<\sum_{j=1}^{N} W_{t-1}^j f_t(x_t'\mid X_{t-1}^j)<\infty. $$ Then the conditional probability of ancestor \(i\) is proportional to \(W_{t-1}^i f_t(x_t'\mid X_{t-1}^i)\).

Proof. For candidate \(i\), the joint contribution from the weighted particle and the transition into the fixed state is \(W_{t-1}^i f_t(x_t'\mid X_{t-1}^i)\). Conditional on that state, the remaining suffix contributes the common factor \(g_t(y_t\mid x_t')\prod_{s=t+1}^{T}f_s(x_s'\mid x_{s-1}')g_s(y_s\mid x_s')\). The factor is independent of \(i\), so it cancels when the candidate contributions are normalized. The assumed positive finite denominator makes this normalization well-defined. \(\square\)

Result 2 (ancestor odds): If candidates \(i\) and \(j\) both have positive ancestor mass, then $$ \frac{P(A_{t-1}=i\mid \cdots)}{P(A_{t-1}=j\mid \cdots)} = \frac{W_{t-1}^i}{W_{t-1}^j} \frac{f_t(x_t'\mid X_{t-1}^i)} {f_t(x_t'\mid X_{t-1}^j)}. $$ Thus the odds combine the candidates’ current weight ratio with their relative transition compatibility.

Derivation. Substitute the normalized ancestor formula for both probabilities. Their common normalizing denominator cancels, leaving the displayed ratio. If the transition density for candidate \(i\) is zero and the denominator is positive, its ancestor probability is zero: that particle cannot reach the fixed reference state in one model transition. \(\square\)

Worked Example: A High-Weight Particle Can Lose

Two particles have states \(0\) and \(2\), each with weight \(0.5\). The transition is \(N(x,0.25)\), and the reference state is \(0.2\). Compare the ancestor probabilities with the probabilities obtained by resampling on weights alone.

The residuals are \(0.2\) for the first particle and \(-1.8\) for the second. The ratio of their transition densities is \(\exp\{-[1.8^2-0.2^2]/(2(0.25))\}=\exp(-6.4)\approx0.001662\). Since the weights are equal, the ratio of ancestor probabilities is also \(0.001662\). The unnormalized masses, after dropping a common factor, are \(0.5\) and \(0.5(0.001662)=0.000831\). Normalizing gives probabilities about \(0.99834\) and \(0.00166\).

Resampling on weights alone would assign probability \(0.5\) to each particle. Ancestor sampling instead strongly favors the particle whose state is transition-compatible with the fixed reference state.

Worked Example: A Discrete Transition with an Impossible Parent

Suppose the state space is \(\{A,B\}\). At time \(t-1\), the particles are in states \(A\) and \(B\), with weights \(0.4\) and \(0.6\). The transition probabilities into the fixed reference state \(A\) are \(P(A\mid A)=0.9\) and \(P(A\mid B)=0.1\). Find the ancestor probabilities.

The unnormalized masses are \(0.4(0.9)=0.36\) for parent \(A\) and \(0.6(0.1)=0.06\) for parent \(B\). Their sum is \(0.42\), so the probabilities are \(0.36/0.42=6/7\approx0.8571\) and \(0.06/0.42=1/7\approx0.1429\). Although the \(B\) particle has the larger filtering weight, its transition to the reference state is less plausible, so it receives the smaller ancestor probability.

Using Ancestor Sampling in Particle Gibbs

In a conditional bootstrap particle filter, one particle is constrained to equal the reference state \(x_t'\) at each time. For the other particles, the filter performs its usual resampling and propagation. At each \(t\geq1\), ancestor sampling draws the parent index of the constrained particle from the formula above. At the final time, an index is drawn according to the final particle weights, and its stored parent links determine the new trajectory.

The key change is that the fixed states need not inherit the old reference path’s parent links. Their states remain fixed during the conditional filter, but their ancestors can be selected anew. In Particle Gibbs with ancestor sampling, the conditional particle construction and these ancestor updates form a valid particle-Gibbs update: the selected trajectory has the smoothing distribution as its invariant target. The extended-target justification rests on the same conditional calculation: in the augmented particle system, the parent choice conditional on the fixed state and particle population has mass proportional to the particle weight times the transition density. A standard PGAS update uses this conditional rather than choosing the parent by weight alone.

This is not a guarantee that every update changes the path substantially. The selected ancestor may coincide with the old one, and resampling can still leave few distinct lineages. The benefit is that the reference path can reconnect to a particle whose state is compatible with its next fixed state, rather than being forced to follow a parent chosen without that information.

Diagnostics, Alternatives, and Limitations

A useful diagnostic is the effective number of ancestor candidates, \(1/\sum_i (W_{t-1}^i f_t(x_t'\mid X_{t-1}^i)/C_t)^2\), where \(C_t=\sum_i W_{t-1}^i f_t(x_t'\mid X_{t-1}^i)\). This value lies between 1 and \(N\): a value near 1 means one candidate dominates, while a larger value indicates that several particles can plausibly connect to the reference state. Also check whether \(C_t\) is positive and finite, and whether a few candidates repeatedly receive nearly all the mass. If all computed masses are zero, the algorithm cannot perform the stated draw; investigate support mismatch, numerical underflow, or an inconsistent transition model rather than silently normalizing zeros.

ApproachAncestor or path choiceMain distinction
Particle Gibbs without ancestor samplingThe constrained particle’s parent is fixed to the reference particle at the preceding time.Does not redraw that parent using transition compatibility.
Particle Gibbs with ancestor samplingUses \(W_{t-1}^i f_t(x_t'\mid X_{t-1}^i)\) for the fixed particle’s parent.Uses next-state compatibility during the conditional particle update.
Backward simulationSelects a path backward from a forward particle population using transition compatibility.Reconstructs a path after the forward pass rather than assigning the fixed particle’s parent online.

Ancestor sampling requires evaluating the transition density at the reference state for each candidate. For a PK SDE, this may require a tractable transition density or a chosen transition approximation. The resulting update targets the model represented by that density; an approximation is not automatically an exact transition for the underlying SDE. In high-dimensional states or over a long interval, all candidate transitions may be very unlikely, so ancestor weights can still collapse. Ancestor sampling improves the information used for the parent choice; it does not remove the need for adequate particle diversity or a suitable transition model.

Common Pitfall

Wrong reasoning: “The reference trajectory is fixed, so its ancestor should be chosen by the particle weights alone; the transition density only matters when propagating new particles.” This ignores the fact that ancestor sampling is precisely choosing which earlier state connects to a fixed later state. Weight-only resampling does not assess that connection.

Correction: Multiply each particle’s weight by its transition density into the fixed next state, then normalize. Conversely, do not multiply by a different observation likelihood for each candidate at time \(t\): the reference state is fixed, so that likelihood is common and cancels. The smoothing target itself must still include all observation factors, including \(g_0(y_0\mid x_0)\) when a time-zero observation is present.

Key takeaway: Ancestor sampling chooses a parent for a fixed reference state by combining filtering weight with transition compatibility. It can reconnect a reference suffix to a more plausible earlier particle, while its practical effectiveness depends on nondegenerate candidate weights and a valid transition model.

Check Your Understanding

Answer each question, then compare with the solutions below.

  1. Why does the smoothing target include \(g_0(y_0\mid x_0)\) when an observation is recorded at time zero?
  2. Two candidate ancestors have equal weights, but one has twice the transition density into the reference state. What is the ratio of their ancestor probabilities?
  3. What does an ancestor effective number close to one indicate?
  4. In one sentence, distinguish ancestor sampling from backward simulation.

Solutions. (1) Because the time-zero observation contributes its likelihood conditional on the initial state; omitting it changes the stated smoothing target unless there is no such observation or it has already been incorporated into the initial density. (2) The first candidate’s ancestor probability is twice the second’s, by the odds formula. (3) Nearly all normalized ancestor mass is concentrated on one candidate. (4) Ancestor sampling chooses the fixed particle’s parent online within conditional particle Gibbs, whereas backward simulation reconstructs a path backward from a forward particle population.

Clinical Biostats