Deriving Fisher Information for ODE Observation Models

Why Derive Information From an ODE?

An ODE model predicts how a biological state changes over time, but a fitted parameter is informed only through the observations made of that state. If a small parameter change produces a substantial, measurable change in predicted observations, those observations can carry information about that parameter. If different parameter changes produce nearly indistinguishable predictions, the observations may constrain some parameter combinations much more strongly than others.

The Fisher information matrix makes this local relationship precise. For independent Gaussian observations with specified variances, it is built from derivatives of the predicted outputs with respect to the parameters, weighted by the observation variances. This tutorial derives that result, shows how the derivatives can be obtained from the ODE, and explains the assumptions required to interpret the matrix as a measure of local precision.

Model and Assumptions

Let the model state \(x(t;\theta)\) satisfy an ODE

\frac{dx}{dt}=f(x(t;\theta),\theta,t), \qquad x(0;\theta)=x_0(\theta),

where \(\theta=(\theta_1,\ldots,\theta_p)^\mathsf{T}\) is the parameter vector. At observation occasion \(j\), the model predicts the scalar output

\mu_j(\theta)=h(x(t_j;\theta),\theta,t_j),

and the observed value is \(Y_j=\mu_j(\theta)+\varepsilon_j\). Assume the errors are independent and \(Y_j\sim N(\mu_j(\theta),\sigma_j^2)\), where each variance \(\sigma_j^2\) is specified, positive, and does not depend on \(\theta\).

Definition: Under this observation model, the Fisher information matrix at parameter value \(\theta\) is the expected outer product of the likelihood score, or equivalently the expected negative Hessian of the log likelihood. It summarizes the local curvature of the expected log likelihood for the specified observation design and error variances.

For the derivation below, assume the ODE solution and observation function are differentiable with respect to the parameters near the parameter value of interest, the solution is well-defined at each observation time, and the stated independent Gaussian error model is appropriate. If the initial condition depends on a parameter, its parameter derivative must be included. Multiple independent scalar outputs can be handled by indexing each observed component separately, with its own specified variance.

From the Likelihood to the Information Matrix

Write \(g_j(\theta)=\nabla_\theta\mu_j(\theta)\), the column vector of derivatives of the \(j\)th predicted observation. Apart from a constant that does not depend on \(\theta\), the log likelihood is

\ell(\theta) =-\frac{1}{2}\sum_j \frac{\{Y_j-\mu_j(\theta)\}^2}{\sigma_j^2}.

Differentiating gives the score, the gradient of the log likelihood:

\nabla_\theta\ell(\theta) =\sum_j \frac{Y_j-\mu_j(\theta)}{\sigma_j^2}g_j(\theta).

The residual \(Y_j-\mu_j(\theta)\) has mean zero and variance \(\sigma_j^2\) when the model is evaluated at the true parameter. Independence makes the expected cross-products between distinct residuals zero. Therefore the expected outer product of the score is

\mathcal I(\theta) =E\!\left[ \nabla_\theta\ell(\theta)\nabla_\theta\ell(\theta)^\mathsf{T} \right] =\sum_j\frac{1}{\sigma_j^2}g_j(\theta)g_j(\theta)^\mathsf{T}.

This derivation also explains why each derivative is weighted by inverse variance: an observation with larger specified noise contributes less information for the same change in its predicted mean. In matrix terms, each observation contributes a positive semidefinite outer product, so the total information matrix is positive semidefinite.

The expected negative Hessian gives the same result. Differentiating the score produces a term \(g_jg_j^\mathsf{T}/\sigma_j^2\) and a term proportional to the residual times the Hessian of \(\mu_j\). The latter has expectation zero because the residual has mean zero. This equality applies under the stated model; the observed negative Hessian for a particular data set need not equal the Fisher information exactly.

Obtaining Output Derivatives From the ODE

The output derivatives in the information formula depend on how the entire ODE trajectory changes with each parameter. For parameter \(\theta_l\), define the state sensitivity \(S_l(t)=\partial x(t;\theta)/\partial\theta_l\). Differentiating the ODE with respect to \(\theta_l\) gives the sensitivity equation

\frac{dS_l}{dt} =f_x(x,\theta,t)S_l+f_{\theta_l}(x,\theta,t), \qquad S_l(0)=\frac{\partial x_0(\theta)}{\partial\theta_l}.

Here \(f_x\) is the Jacobian of \(f\) with respect to the state, and \(f_{\theta_l}\) is its partial derivative with respect to the parameter, holding the state fixed. By the chain rule, the corresponding derivative of the predicted output is

\frac{\partial \mu_j}{\partial\theta_l} =h_x(x(t_j),\theta,t_j)S_l(t_j) +h_{\theta_l}(x(t_j),\theta,t_j).

Thus the calculation proceeds from ODE sensitivities to output derivatives and then to the information matrix. The direct term \(h_{\theta_l}\) is needed when the observation function itself depends on a parameter, even if that parameter has no direct effect on the state equation.

Worked Example: Information for a First-Order Decay Model

Consider the hypothetical ODE \(dx/dt=-kx\), with \(x(0)=A\), and observations \(Y_j=x(t_j)+\varepsilon_j\) at \(t=0,1,2\) time units. Suppose the independent errors have specified variance \(\sigma^2=1\) concentration-unit squared. Calculate the information at \(A=10\) concentration units and \(k=1\) per time unit. The predicted output is \(\mu(t)=Ae^{-kt}\), so its derivative vector in the parameter order \((A,k)\) is

g(t)= \begin{pmatrix} e^{-kt}\\ -At e^{-kt} \end{pmatrix}.

At \(t=0\), this vector is \((1,0)^\mathsf{T}\); at \(t=1\), it is \((e^{-1},-10e^{-1})^\mathsf{T}\); and at \(t=2\), it is \((e^{-2},-20e^{-2})^\mathsf{T}\). Since the variance is one, summing the three outer products gives

\begin{aligned} \mathcal I(A,k) &= \begin{pmatrix} 1+e^{-2}+e^{-4} & -10e^{-2}-20e^{-4}\\ -10e^{-2}-20e^{-4} & 100e^{-2}+400e^{-4} \end{pmatrix}\\ &\approx \begin{pmatrix} 1.1537 & -1.7197\\ -1.7197 & 20.8598 \end{pmatrix}. \end{aligned}

For example, the off-diagonal entry is \(-10(0.1353353)-20(0.0183156)\approx-1.7197\), and the lower-right entry is \(100(0.1353353)+400(0.0183156)\approx20.8598\). The determinant is approximately \(21.1077\), so the matrix is nonsingular at this parameter value. Its inverse is approximately

\mathcal I(A,k)^{-1} \approx \begin{pmatrix} 0.9883 & 0.0815\\ 0.0815 & 0.05465 \end{pmatrix}.

Under the usual regularity conditions and a sufficiently informative design, the inverse information is a local approximation to the covariance matrix of an efficient parameter estimate. In this approximation, the marginal variance for \(A\) is \(0.9883\) concentration units squared, and the marginal variance for \(k\) is \(0.05465\) per time unit squared. The off-diagonal entry is their approximate covariance, in concentration units per time unit. These are local approximations, not guaranteed finite-sample variances or confidence intervals.

Worked Examples With Other Observation Designs

Worked Example: A Baseline and a Constant Production Rate

Suppose a hypothetical state follows \(dx/dt=v\), with \(x(0)=b\), and is observed at times \(0,1,2\) time units. Then \(\mu(t)=b+vt\), with \(b\) in concentration units and \(v\) in concentration units per time unit. Assume independent errors with known variance \(4\) concentration-unit squared. At the parameter point \(b=2\), \(v=1\), suppose the observed values are \(2.1,3.0,4.2\) concentration units. The predicted-output derivative vector is \(g(t)=(1,t)^\mathsf{T}\). The observations determine the residuals, but for this specified-variance Gaussian model the information depends on the derivatives and design, not on the realized residuals.

The three derivative vectors are \((1,0)^\mathsf{T}\), \((1,1)^\mathsf{T}\), and \((1,2)^\mathsf{T}\). Weighting each outer product by \(1/4\) gives

\mathcal I(b,v) =\frac{1}{4} \left[ \begin{pmatrix}1&0\\0&0\end{pmatrix} +\begin{pmatrix}1&1\\1&1\end{pmatrix} +\begin{pmatrix}1&2\\2&4\end{pmatrix} \right] = \begin{pmatrix}0.75&0.75\\0.75&1.25\end{pmatrix}.

Its determinant is \(0.75(1.25)-0.75^2=0.375\), and its inverse is \(\begin{pmatrix}3.333&-2\\-2&2\end{pmatrix}\). The negative covariance in this local approximation reflects the fact that a higher baseline can be offset, in predictions over time, by a lower production rate. The entries have the units implied by the parameter pair \((b,v)\); for example, the approximate variance of \(v\) is \(2\) concentration-unit squared per time-unit squared.

Worked Example: Sensitivity Equations and One Observation

Return to \(dx/dt=-kx\), \(x(0)=A\), and observe the state once at \(t=1\) time unit. The sensitivity equations, in the parameter order \((A,k)\), are

\frac{dS_A}{dt}=-kS_A,\quad S_A(0)=1; \qquad \frac{dS_k}{dt}=-x-kS_k,\quad S_k(0)=0.

The first equation gives \(S_A(t)=e^{-kt}\). Substituting \(x(t)=Ae^{-kt}\) into the second gives \(S_k(t)=-At e^{-kt}\), which can also be verified by differentiating that expression and checking its initial value. At \(A=10\), \(k=1\), and \(t=1\), the output derivative vector is approximately \((0.367879,-3.67879)^\mathsf{T}\). If the specified observation variance is \(0.25\) concentration-unit squared, the information is

\mathcal I(A,k) =\frac{1}{0.25} \begin{pmatrix} 0.367879^2 & 0.367879(-3.67879)\\ 0.367879(-3.67879) & (-3.67879)^2 \end{pmatrix} \approx \begin{pmatrix} 0.5413&-5.413\\ -5.413&54.13 \end{pmatrix}.

This matrix has determinant zero because it is a single weighted outer product of one derivative vector. One observation contributes information in only one parameter-space direction; by itself it cannot provide independent local information about both \(A\) and \(k\). This is a property of the design and model derivatives, not evidence that the observed value was unusually noisy.

Interpretation, Diagnostics, and Limits

The information matrix describes local curvature at a particular parameter value. Its interpretation depends on the assumptions: observations must have the stated independent Gaussian errors, the variances must be correctly specified and fixed with respect to the parameters, and the parameter-to-output map must be sufficiently smooth. If observations are correlated, the independent-error sum is not the appropriate formula. If the variance depends on the mean or on parameters, differentiating the full likelihood introduces additional terms; the formula derived here should not be reused unchanged.

Numerical calculations also require care. Sensitivities should correspond to the same ODE solution and observation function used to compute the predictions. Solver error, inaccurate derivatives, inconsistent units, or poorly chosen parameter scales can distort the entries and the apparent conditioning of the matrix. A useful basic check is to compare sensitivity-based output derivatives with small, scale-aware perturbations of the parameters and verify that the ODE solver tolerances are sufficiently tight. This comparison checks a calculation; it does not replace the assumed statistical model.

When the information matrix is nonsingular, its inverse can provide a local covariance approximation under regularity conditions. That approximation is most persuasive when the likelihood is well described by a quadratic near its maximum and the parameter value is away from problematic boundaries. It may be misleading with strong nonlinearity, weakly informed directions, parameter boundaries, or small samples. A singular matrix has no inverse, so an inverse-information covariance approximation is unavailable without additional information or a different parameterization.

The main alternatives answer related but distinct questions. The observed Hessian uses the curvature of the realized log likelihood, whereas Fisher information averages curvature under the model; in finite data, they can differ. Direct likelihood-based uncertainty methods assess the likelihood beyond its local quadratic approximation, while simulation-based approaches can examine repeated-sample behavior under a specified data-generating model. Neither changes the assumptions required for the Fisher information formula itself. The information matrix is valuable precisely as a compact local summary, not as a universal substitute for checking the likelihood or the observation model.

Common Mistakes

Omitting the variance weight. The contribution is \(g_jg_j^\mathsf{T}/\sigma_j^2\), not just the outer product. Larger observation variance means less information.
Using state derivatives in place of output derivatives. If the measurement function transforms the state or depends directly on parameters, apply the chain rule to obtain derivatives of the predicted observation.
Forgetting the initial-condition derivative. A parameter-dependent initial state contributes through \(S_l(0)=\partial x_0/\partial\theta_l\).
Treating an inverse as an exact finite-sample covariance. The inverse information is generally a local approximation under regularity conditions, not a guarantee about a particular fitted data set.
Using the fixed-variance formula when the error variance depends on parameters. In that case the likelihood’s variance terms also change with the parameters, so the derivation must be modified.

Key takeaway: For independent Gaussian observations with specified, parameter-independent variances, each predicted output contributes its parameter-derivative outer product divided by its variance. The ODE sensitivity equations supply those derivatives; the resulting information is a local summary whose precision interpretation depends on the model assumptions and regularity conditions.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions.

Why is the contribution of an observation with variance \(\sigma_j^2\) weighted by \(1/\sigma_j^2\)?
For the decay model with parameters \((A,k)\), what information does a single observation at \(t=0\) provide about \(A\) and \(k\), if its specified variance is \(\sigma^2\)?
In the baseline-and-rate model, add an observation at \(t=3\) to the design with times \(0,1,2\). What is the new information matrix if each observation variance remains \(4\)?
Why is it not generally valid to use the fixed-variance formula unchanged when the observation variance depends on a model parameter?

Solutions. (1) The score contribution is the residual times the output derivative divided by the variance; taking its expected outer product uses \(E[(Y_j-\mu_j)^2]=\sigma_j^2\), leaving the derivative outer product divided by \(\sigma_j^2\). (2) At \(t=0\), \(g=(1,0)^\mathsf{T}\), so the contribution is \(\begin{pmatrix}1/\sigma^2&0\\0&0\end{pmatrix}\); the observation locally informs \(A\), but not \(k\). (3) The four derivative vectors have sums \(\sum 1=4\), \(\sum t=6\), and \(\sum t^2=14\). Dividing their assembled outer-product sum by \(4\) gives \(\begin{pmatrix}1&1.5\\1.5&3.5\end{pmatrix}\). (4) A parameter-dependent variance changes the likelihood through both its residual term and its variance term, so differentiating the full likelihood produces additional contributions not included in the fixed-variance derivation.