Tutorials › Stochastic Differential Equations in Pharmacokinetics › Exact Densities for Time-Inhomogeneous Linear Diffusions

Likelihood for Discretely Observed Diffusions · Tutorial 57 of 100

Exact Densities for Time-Inhomogeneous Linear Diffusions

Learn to obtain exact conditional transition means and variances for time-inhomogeneous linear diffusions, then use them to evaluate a Gaussian transition likelihood.

Advanced 11 min read

What You'll Learn

  • State conditions that ensure the scalar linear SDE has an exact Gaussian transition law.
  • Derive its transition mean using the fundamental solution of the time-varying drift.
  • Derive its transition variance from the deterministic Brownian-integral kernel.
  • Include a known, time-varying input in the conditional mean.
  • Evaluate a transition density and log-likelihood from an explicit PK-style model.
  • Identify when the Gaussian density is degenerate or unsuitable for concentration data.

Why Exact Transition Densities Matter

For discretely observed diffusion data, a likelihood can be built from the conditional density of each observation given the preceding state. The earlier material on constructing a discrete-time likelihood introduces this principle; here the key task is to obtain that density exactly for a useful class of models. In a scalar linear diffusion, coefficients may change with time and a known input may vary during a dosing interval, yet the transition remains Gaussian under suitable conditions.

This exact result is valuable for likelihood calculation: it avoids replacing the continuous-time model by a time-stepping approximation. It also clarifies which parts of the density depend on the drift, the diffusion, the observation interval, and a known input.

Model and Transition Formula

Let \(W_t\) be standard Brownian motion, and consider the scalar SDE on an interval \([s,t]\):

$$ dX_r=\{a(r)X_r+b(r)+c(r)u(r)\}\,dr+\sigma(r)\,dW_r,\qquad s\leq r\leq t. $$

Here \(a,b,c,\sigma\) are known deterministic functions, and \(u\) is a known deterministic input, such as an infusion rate. The term \(c(r)u(r)\) could be combined with \(b(r)\); writing it separately highlights the input contribution. Assume \(a,b,cu\) are integrable on \([s,t]\), and \(\sigma\) is square-integrable there. These conditions ensure the integrals below are finite and the linear SDE has a well-defined solution. We condition on \(X_s=x\).

Definition: The fundamental solution of the scalar drift on \([s,t]\) is \(\Phi(t,r)=\exp\!\left(\int_r^t a(v)\,dv\right)\), for \(s\leq r\leq t\). Define the transition mean \(m(t\mid s,x)\) and variance \(q(t\mid s)\) by \(m(t\mid s,x)=\Phi(t,s)x+\int_s^t\Phi(t,r)\{b(r)+c(r)u(r)\}\,dr\) and \(q(t\mid s)=\int_s^t\Phi(t,r)^2\sigma(r)^2\,dr\).

The factor \(\Phi(t,r)\) describes how a perturbation to the state at time \(r\) is carried forward to time \(t\). It also weights both the accumulated input and the Brownian noise. In particular, noise introduced earlier in the interval is scaled by the drift over the remaining time.

Derivation and Gaussian Transition Result

Theorem: Suppose \(a,b,cu\) are deterministic and integrable on \([s,t]\), and \(\sigma\) is deterministic and square-integrable there. For a fixed starting value \(X_s=x\), the solution has conditional distribution \(X_t\mid X_s=x\sim N\!\left(m(t\mid s,x),q(t\mid s)\right)\). If \(q(t\mid s)>0\), its transition density at \(y\) is \(p(y,t\mid x,s)=\{2\pi q(t\mid s)\}^{-1/2}\exp[-\{y-m(t\mid s,x)\}^2/\{2q(t\mid s)\}]\). If \(q(t\mid s)=0\), the transition is degenerate at its deterministic mean and has no ordinary Gaussian density with respect to Lebesgue measure.

To derive the result, use the integrating factor for the linear drift. Multiplying the SDE by \(\exp(-\int_s^r a(v)\,dv)\) and applying the product rule gives a differential whose drift no longer contains \(X_r\). Integrating from \(s\) to \(t\) and rearranging yields

$$ X_t=\Phi(t,s)x+\int_s^t\Phi(t,r)\{b(r)+c(r)u(r)\}\,dr+\int_s^t\Phi(t,r)\sigma(r)\,dW_r. $$

The first two terms are deterministic given \(X_s=x\). The last term is an Itô integral with a deterministic square-integrable integrand. Such an integral is a centered Gaussian random variable with variance equal to the integral of the squared integrand. Therefore its variance is \(q(t\mid s)\), proving the stated conditional distribution. When that variance is positive, the normal density formula gives the transition density. When it is zero, the stochastic integral is zero almost surely, so \(X_t=m(t\mid s,x)\) almost surely. This proves the result. \(\square\)

The deterministic-integrand condition matters: a stochastic integral with a state-dependent integrand is not generally Gaussian. Likewise, the formula is conditional on a known starting state. If the starting state is random, this conditional transition law still applies, but the unconditional distribution requires averaging over the distribution of that state.

Worked Examples

Worked Example: Time-Varying Drift and Diffusion

On \(0\leq r\leq1\), take \(a(r)=-1/(1+r)\), \(b(r)=0\), \(c(r)u(r)=0\), and \(\sigma(r)=1/(1+r)\). Let \(X_0=2\). Find the distribution of \(X_1\).

First, \(\Phi(1,r)=\exp\{-\int_r^1(1+v)^{-1}dv\}=(1+r)/2\). Thus \(\Phi(1,0)=1/2\), and the mean is \(m(1\mid0,2)=2(1/2)=1\).

For the variance, the squared transition kernel times the squared diffusion coefficient is \(\Phi(1,r)^2\sigma(r)^2=\{(1+r)/2\}^2\{1/(1+r)\}^2=1/4\). Therefore \(q(1\mid0)=\int_0^1(1/4)\,dr=1/4\), and \(X_1\mid X_0=2\sim N(1,1/4)\). The standard deviation is \(1/2\). Both coefficients vary with time, but the transition law is still exactly Gaussian.

Worked Example: A Known Input and One-Transition Likelihood

Suppose \(dX_r=(-0.4X_r+0.3)\,dr+0.2\,dW_r\) over an interval of length \(0.5\), with \(X_s=1\). The constant \(0.3\) is a known input contribution to the drift. Evaluate the transition density at the observed endpoint \(y=1.083073\).

The fundamental solution is \(\Phi(s+0.5,r)=e^{-0.4(s+0.5-r)}\). The mean is \(m=e^{-0.2}(1)+0.3\int_0^{0.5}e^{-0.4(0.5-v)}\,dv =e^{-0.2}+0.75(1-e^{-0.2})\). Using \(e^{-0.2}\approx0.818730753\), this gives \(m\approx0.954682688\).

The variance is \(q=0.2^2\int_0^{0.5}e^{-0.8(0.5-v)}\,dv =0.04(1-e^{-0.4})/0.8\approx0.016483998\). Hence the standard deviation is approximately \(0.128390\). The residual is \(y-m\approx1.083073-0.954682688=0.128390312\), approximately one standard deviation. Substitution into the normal density gives \(p(y,s+0.5\mid1,s)=\{2\pi q\}^{-1/2}\exp[-(y-m)^2/(2q)]\approx1.8847\). The corresponding one-transition log-likelihood is \(\log p\approx0.6337\). A density can exceed one; it is probability per unit of the state’s measurement scale, not the probability of a single exact value.

Worked Example: A Piecewise Time-Varying Elimination Rate

Let \(dX_r=a(r)X_r\,dr+0.1\,dW_r\) on \([0,1]\), with \(X_0=3\), where \(a(r)=-1\) for \(0\leq r<0.5\) and \(a(r)=-2\) for \(0.5\leq r\leq1\). Compute the transition mean and variance.

The total drift integral is \(-0.5-1=-1.5\), so \(\Phi(1,0)=e^{-1.5}\) and \(m=3e^{-1.5}\approx0.669390\). For \(0\leq r<0.5\), \(\Phi(1,r)=e^{-(0.5-r)-1}=e^{-1.5+r}\); for \(0.5\leq r\leq1\), \(\Phi(1,r)=e^{-2(1-r)}\). Splitting the variance integral at the coefficient change gives

$$ q=0.01\left\{\int_0^{0.5}e^{-3+2r}\,dr+\int_{0.5}^1 e^{-4(1-r)}\,dr\right\} =0.01\left\{\frac{e^{-2}-e^{-3}}{2}+\frac{1-e^{-2}}{4}\right\} \approx0.0025894. $$

Thus \(X_1\mid X_0=3\sim N(0.669390,0.0025894)\), with standard deviation approximately \(0.05089\). Splitting the integral is essential: using either drift coefficient over the whole interval would produce the wrong attenuation and variance.

From Transition Densities to Likelihoods

Suppose the state is observed exactly at times \(t_0<t_1<\cdots<t_n\), with values \(x_0,\ldots,x_n\). For fixed model parameters, the Markov property and the transition densities give the conditional likelihood

$$ L=\prod_{i=1}^{n}p(x_i,t_i\mid x_{i-1},t_{i-1}),\qquad \ell=\sum_{i=1}^{n}\left[-\frac12\log\{2\pi q_i\} -\frac{\{x_i-m_i\}^2}{2q_i}\right], $$

where \(m_i=m(t_i\mid t_{i-1},x_{i-1})\) and \(q_i=q(t_i\mid t_{i-1})\). If parameters are estimated, these quantities are evaluated at each candidate parameter value. The initial-state contribution must also be included if \(x_0\) is not conditioned on and has a modelled distribution. This likelihood concerns exact state observations; an observation-error model would require additional terms and is not part of this transition formula.

Diagnostics, Alternatives, and Limitations

Before using the density, check that the coefficients and input used in each interval match the actual dosing and infusion schedule, that the integrals are evaluated over the full observation interval, and that every \(q_i\) is positive when a continuous transition density is required. A computed negative variance signals an error, since \(q_i\) is an integral of a nonnegative function. Very small variances can also make likelihood calculations numerically unstable; evaluate the log-density directly rather than multiplying many densities.

The main alternative in this setting is a time-discretization approximation, such as Euler–Maruyama, followed by an approximate transition likelihood. That approach can be useful when the model is nonlinear or its transition law is unavailable, but it introduces discretization error. For a scalar linear SDE with deterministic coefficients, the integral formulas above give the exact transition and should be preferred when those integrals can be computed accurately. For piecewise-defined coefficients, integrate separately on each piece. For more complicated time-varying functions, numerical quadrature may be needed; its accuracy should be checked by refining the quadrature.

Exact Gaussian transitions do not guarantee a realistic concentration model. A Gaussian distribution has support on the entire real line, so it assigns positive probability to negative states. If concentrations must remain nonnegative, this may be scientifically unacceptable even though the transition calculation is exact. The result also does not apply unchanged to nonlinear state-dependent diffusion, random inputs, or a model whose coefficients depend on the evolving state in a way that destroys the linear structure.

Common pitfall: Using \(\exp\{\int_s^t a(v)\,dv\}\) as the noise multiplier for every Brownian increment gives the wrong variance. Noise added at time \(r\) is propagated only from \(r\) to \(t\), so its multiplier is \(\Phi(t,r)\), and the variance integrates \(\Phi(t,r)^2\sigma(r)^2\) over \(r\).

Check Your Understanding

Try each question, then check the brief solutions below.

  1. For a scalar linear SDE with deterministic coefficients, what two integrals determine the conditional transition mean and variance?
  2. Why does a zero transition variance not yield the ordinary normal density formula?
  3. For the piecewise example, what is the transition mean if the starting value is changed from \(3\) to \(2\), with all coefficients unchanged?
  4. What is the log-likelihood contribution of one observed transition when its residual \(y-m\) is zero and its variance is \(q>0\)?
  5. Give one reason an exact Gaussian transition might still be inappropriate for a concentration model.

Solutions. (1) The mean is \(\Phi(t,s)x+\int_s^t\Phi(t,r)\{b(r)+c(r)u(r)\}\,dr\); the variance is \(\int_s^t\Phi(t,r)^2\sigma(r)^2\,dr\). (2) The state is deterministic at the mean, a point mass, rather than a distribution with a Lebesgue density of the stated normal form. (3) The mean is \(2e^{-1.5}\approx0.446260\). (4) It is \(-\tfrac12\log(2\pi q)\), since the squared-residual term is zero. (5) A Gaussian transition assigns positive probability to negative concentrations.

Clinical Biostats