End-to-End Bayesian Particle Analysis of a Nonlinear PK SDE

From a PK SDE to Bayesian summaries

A particle analysis links several distinct calculations: a continuous-time pharmacokinetic model, a discrete-time approximation used by the filter, an estimate of the data likelihood, and a posterior distribution for unknown parameters. Keeping those objects separate matters. A particle filter produces a random likelihood estimate, not the exact likelihood; a Bayesian parameter sampler must account for that randomness; and a latent-state estimate must use the particle weights appropriate to the observation being summarized.

The construction below uses a one-compartment amount model with known infusion input. It illustrates an end-to-end analysis without implying that one model or numerical scheme is suitable for every drug or sampling schedule. The particle-filter likelihood estimate is unbiased for the likelihood of the specified discretized state-space model under the assumptions stated below. That is not, by itself, a claim of unbiasedness for the exact continuous-time SDE likelihood.

Model specification and the numerical transition

Let \(A_t>0\) be the amount of drug in a compartment, \(R(t)\) a known infusion rate, \(k>0\) an elimination-rate parameter, and \(\sigma>0\) a diffusion scale. Between dose impulses, consider

The multiplicative diffusion makes the amount equation stochastic while allowing a positive state to be represented by \(X_t=\log A_t\). Itô's formula gives

The term \(-\sigma^2/2\) is required by the transformation; omitting it changes the model. At a known bolus dose \(D>0\), amount jumps from \(A_{t^-}\) to \(A_{t^-}+D\), so the corresponding log-state update is \(X_{t^+}=\log(e^{X_{t^-}}+D)\). The filter applies this known transformation at the event time rather than treating the dose as an unknown random disturbance.

Suppose observations are positive measured concentrations \(C_j^{\mathrm{obs}}\) at times \(t_j\), and define \(Z_j=\log C_j^{\mathrm{obs}}\). With known volume \(V\) and log-scale assay standard deviation \(\tau>0\), use the observation model

with observation errors conditionally independent given the states. Thus the observation density is \(g_j(z\mid x)=\{ \tau\sqrt{2\pi}\}^{-1}\exp[-(z-x+\log V)^2/(2\tau^2)]\). A lognormal assay model for concentration is being used; it is not the same as additive Gaussian error on the concentration scale.

For inference, let \(\theta=(k,\sigma,V,\tau)\), or a subset of these quantities if some are known. Specify a proper prior density \(p(\theta)\) on its stated support, an initial-state distribution, the observation model, known inputs and event times, and a numerical time grid. Assume the SDE has a nonexplosive solution over the study period and the coefficients are sufficiently regular between input events for the stated model to be meaningful. The Euler filter below defines an approximate state-space model with Gaussian transitions; posterior results based on it are conditional on that discretization.

Particle-filter likelihood estimation

Write the discrete-time transition density, including any Euler substeps and known event updates, as \(q_j(u\mid x;\theta)\). At observation time \(t_j\), propagate particles through this transition and assign the incremental observation weight \(g_j(Z_j\mid U_j^{(i)};\theta)\). Resampling uses the normalized weights. The particle estimate of the likelihood is the product of the mean incremental weights:

Here \(N\) is the particle count and \(U_j^{(i)}\) are the propagated particles immediately before weighting at time \(t_j\). The initial particle population is sampled from the specified initial-state distribution. The formula assumes the bootstrap construction; other proposal distributions require the corresponding importance correction in each weight. Resampling changes the particle population used at later steps, but does not replace the mean-weight factor already accumulated for the likelihood estimate.

To prove the result, let \(\widehat Z_j\) be the product of mean weights through time \(t_j\), with \(\widehat Z_0=1\). Let \(\widehat\pi_{j-1}\) be the normalized weighted empirical distribution after observing \(Z_{j-1}\). Define

Condition on the particle system after time \(j-1\). For any nonnegative measurable \(f\), the product of the new mean weight and the new post-observation weighted average of \(f\) is exactly \(N^{-1}\sum_i g_j(Z_j\mid U_j^{(i)})f(U_j^{(i)})\). The conditional expectation of this sum is \(\widehat\pi_{j-1}(Q_jf)\): each resampled ancestor has the current normalized weights in expectation, and propagation then follows \(q_j\). Hence

Taking expectations and applying this identity successively gives the ordinary iterated-integral expression for the unnormalized filtering distribution at time \(j\), integrated against \(f\). At \(j=0\), the identity is the initial-state expectation by construction. Thus induction establishes the expression at every \(j\). Setting \(f=1\) yields \(E[\widehat Z_J]=L(\theta)\), as claimed. Nonnegativity follows because all factors are nonnegative. \(\square\)

This argument uses the post-observation weighted average. It does not assert that the mean weight multiplied by an unweighted average of \(f\) over the same propagated particles has the required expectation. Those two random averages are dependent, and that product generally has an extra covariance contribution.

Sampling parameters with a random likelihood estimate

A Metropolis-Hastings algorithm that substitutes \(\widehat L_N(\theta)\) for \(L(\theta)\) can target the exact posterior of the discretized model, provided the nonnegative estimate is unbiased and the algorithm retains its random numbers or particles as part of its state. This is the pseudo-marginal method. Let \(U\) denote all random variables used by the particle filter, with density \(m_\theta(u)\), and write the estimate as \(\widehat L(\theta,u)\). The extended target is

Integrating over \(u\) gives \(p(\theta)L(\theta)\), because \(E_{m_\theta}[\widehat L(\theta,U)]=L(\theta)\). Thus its marginal in \(\theta\) is the desired parameter posterior. For a proposal \(\theta'\sim r(\theta'\mid\theta)\), draw fresh \(u'\sim m_{\theta'}\). The Metropolis-Hastings acceptance probability is

The \(m_\theta\) terms cancel because the proposed random variables are drawn from their parameter-dependent filter distribution. Crucially, when a proposal is rejected, retain the current \(u\) and likelihood estimate. Recomputing the current estimate at every iteration would not implement this extended-target Markov chain. In practice, compute log likelihood estimates by summing log mean weights to avoid numerical underflow; this does not remove the randomness of the estimate.

Latent-state summaries and diagnostics

At an observation time, a filtering expectation of a function \(h\) is estimated by \(\sum_i w_j^{(i)}h(U_j^{(i)})\), where \(w_j^{(i)}\) are the normalized post-observation weights. For the state mean, use \(h(x)=x\); for a variance, use the weighted second central moment. Weighted quantiles can summarize asymmetry. If summaries are needed for amount or concentration rather than log amount, transform each particle first: for example, the estimated mean amount is \(\sum_iw_j^{(i)}e^{U_j^{(i)}}\), not the exponential of the mean log amount.

Several diagnostics help distinguish Monte Carlo trouble from a model or discretization problem. Track effective sample size, \(1/\sum_i(w_j^{(i)})^2\), alongside the number of distinct particles after resampling. A sudden collapse indicates that a few particles dominate; repeated resampling can also leave few distinct ancestors and make past-state summaries unreliable. Increase \(N\), assess repeated runs with independent random seeds, and inspect whether likelihood estimates are stable enough for parameter sampling. In pseudo-marginal sampling, very noisy likelihood estimates can cause sticky chains even when the acceptance rule is implemented correctly.

A particle filter is not the only option. Linear Gaussian state-space models admit Kalman filtering; nonlinear Gaussian approximations such as the extended Kalman filter can be cheaper but approximate the filtering distribution. Transition-density approximations can support likelihood-based inference without particle resampling, but their accuracy depends on the approximation regime. Particle Gibbs targets latent paths by a different MCMC construction and can be useful for path inference, but it is not the same operation as marginalizing the path with a particle likelihood estimate in pseudo-marginal parameter sampling.

The particle approach also has limits. A correct algorithm for the Euler model does not eliminate Euler discretization error; grid refinement and sensitivity checks are needed when that error may matter. Finite particle populations can yield high-variance likelihood estimates, while poor observation models or incompatible initial-state assumptions can produce weight collapse. Finally, a narrow posterior from one run is not proof that the model is adequate: compare runs, inspect posterior predictive behavior, and state whether reported latent summaries are filtering or smoothing summaries.