Inference for nonlinear dynamical systems

Maximum Likelihood via Iterated Filtering

A state space model consists of an unobserved Markov process, x_t, called the state process and an observation process y_t. Here, x_t takes values in the state space ℝ^dx, and y_t in the observation space ℝ^dy. The processes depend on an (unknown) vector of parameters, θ, in ℝ^dθ. Observations take place at discrete times, t = 1, …, T; we write the vector of concatenated observations as y_1:T = (y₁, …, y_T); y_1:0 is defined to be the empty vector. A model is completely specified by the conditional transition density f(x_t|x_t−1, θ), the conditional distribution of the observation process f(y_t|y_1:t−1, x_1:t, θ) = f(y_t|x_t, θ), and the initial density f(x₀|θ). Throughout, we adopt the convention that f(·|·) is a generic density specified by its arguments, and we assume that all densities exist. The likelihood is given by the identity f(y_1:T|θ) = ∏_t=1^T f(y_t|y_1:t−1, θ). The state process, x_t, may be defined in continuous or discrete time, but only its distribution at the discrete times t = 1, …, T directly affects the likelihood. The challenge is to find the maximum of the likelihood as a function of θ.

The basic idea of our method is to replace the original model with a closely related model, in which the time-constant parameter θ is replaced by a time-varying process θ_t. The densities f(x_t|x_t−1, θ), f(y_t|x_t, θ), and f(x₀|θ) of the time-constant model are replaced by f(x_t|x_t−1, θ_t−1), f(y_t|x_t, θ_t), and f(x₀|θ₀). The process θ_t is taken to be a random walk in ℝ^dθ. Our main algorithm (Procedure 1 below) and its justification (Theorem 1 below) depend only on the mean and variance of the random walk, which are defined to be

In practice, we use the normal distributions specified by Eq. 1. Here, σ and c are scalar quantities, and the new model in Eq. 1 is identical to the fixed-parameter model when σ = 0. The objective is to obtain an estimate of θ by taking the limit as σ → 0. Σ is typically a diagonal matrix giving the respective scales of each component of θ; more generally, it can be taken to be an arbitrary positive–definite symmetric matrix. Procedure 1 below is standard to implement, as the computationally challenging step 2(ii) requires using only well studied filtering techniques (1, 13) to calculate

for t = 1, …, T. We call this procedure maximum likelihood via iterated filtering (MIF).

Procedure 1. (MIF)

1. Select starting values θ̂⁽¹⁾, a discount factor 0 < α < 1, an initial variance multiplier c², and the number of iterations N.

2. For n in 1, …, N

   1. Set σ_n = αⁿ⁻¹. For t = 1, …, T, evaluate θ̂_t⁽ⁿ⁾ = θ̂_t(θ̂⁽ⁿ⁾, σ_n) and V_t,n = V_t(θ̂⁽ⁿ⁾, σ_n).
   2. Set θ̂⁽ⁿ⁺¹⁾ = θ̂⁽ⁿ⁾ + V_1,n Σ_t=1^T V_t,n⁻¹(θ̂_t⁽ⁿ⁾ − θ̂_t−1⁽ⁿ⁾), where θ̂₀⁽ⁿ⁾ = θ̂⁽ⁿ⁾.

3. Take θ̂^(N+1) to be a maximum likelihood estimate of the parameter θ for the fixed parameter model.

The quantities θ̂_t⁽ⁿ⁾ can be considered local estimates of θ, in the sense that they depend most heavily on the observations around time t. The updated estimate is a weighted average of the values θ̂_t⁽ⁿ⁾, as explained below and in Supporting Text, which is published as supporting information on the PNAS web site. A weighted average of local estimates is a heuristically reasonable estimate for the fixed “global” parameter θ. In addition, taking a weighted average and iterating to find a fixed point obviates the need for a separate optimization algorithm. Theorem 1 asserts that (under suitable conditions) the weights in Procedure 1 result in a maximum likelihood estimate in the limit as σ → 0. Taking a weighted average is not so desirable when the information about a parameter is concentrated in a few observations: this occurs for initial value parameters, and modifications to Procedure 1 are appropriate for these parameters (Supporting Text).

Procedure 1, with step 2(i) implemented using a sequential Monte Carlo method (see ref. 13 and Supporting Text), permits flexible modeling in a wide variety of situations. The methodology requires only that Monte Carlo samples can be drawn from f(x_t| x_t−1), even if only at considerable computational expense, and that f(y_t| x_t, θ) can be numerically evaluated. We demonstrate this below with an analysis of cholera data, using a mechanistic continuous-time model. Sequential Monte Carlo is also known as “particle filtering” because each Monte Carlo realization can be viewed as a particle's trajectory through the state space. Each particle filtering step prunes particles in a way analogous to Darwinian selection. Particle filtering for fixed parameters, like natural selection without mutation, is rather ineffective. This explains heuristically why Procedure 1 is necessary to permit inference for fixed parameters via particle filtering. However, Procedure 1 and the theory given below apply more generally, and could be implemented using any suitable filter.

Example: A Compartment Model for Cholera

In a standard epidemiological approach (29, 30), the population is divided into disease status classes. Here, we consider classes labeled susceptible (S), infected and infectious (I), and recovered (R¹, …, R^k). The k recovery classes allow flexibility in the distribution of immune periods, a critical component of cholera modeling (20). Three additional classes B, C, and D allow for birth, cholera mortality, and death from other causes, respectively. S_t denotes the number of individuals in S at time t, with similar notation for other classes. We write N_t^SI for the integer-valued process (or its real-valued approximation) counting transitions from S to I, with corresponding definitions of N_t^BS, N_t^SD, etc. The model is shown diagrammatically in Fig. 1. To interpret the diagram in Fig. 1 as a set of coupled stochastic equations, we write

The population size P_t is presumed known, interpolated from census data. Transmission is stochastic, driven by Gaussian white noise

In Eq. 3, we ignore stochastic effects at a demographic scale (infinitesimal variance proportional to S_t). We model the remaining transitions deterministically

Time is measured in months. Seasonality of transmission is modeled by log(β_t) = Σ_k=0⁵ b_ks_k(t), where {s_k(t)} is a periodic cubic B-spline basis (31) defined so that s_k(t) has a maximum at t = 2k and normalized so that Σ_k=0⁵ s_k(t) = 1; ε is an environmental stochasticity parameter (resulting in infinitesimal variance proportional to S_t²); ω corresponds to a non-human reservoir of disease; β_tI_t/P_t is human-to-human transmission; 1/γ gives mean time to recovery; 1/r and 1/(kr²) are respectively the mean and variance of the immune period; 1/m is the life expectancy excluding cholera mortality, and m_c is the mortality rate for infected individuals. The equation for dN_t^BS in Eq. 4 is based on cholera mortality being a negligible proportion of total mortality. The stochastic system was solved numerically using the Euler–Maruyama method (32) with time increments of 1/20 month. The data on observed mortality were modeled as y_t ∼ 𝒩[C_t − C_t−1, τ²(C_t − C_t−1)²], where C_t = N_t^IC. In the terminology given above, the state process x_t is a vector representing counts in each compartment.

Results

Discussion

Procedure 1 depends on the viability of solving the filtering problem, i.e., calculating θ̂_t and V_t in Eq. 2. This is a strength of the methodology, in that the filtering problem has been extensively studied. Filtering does not require stationarity of the stochastic dynamical system, enabling covariates (such as P_t) to be included in a mechanistically plausible way. Missing observations and data collected at irregular time intervals also pose no obstacle for filtering methods. Filtering can be challenging, particularly in nonlinear systems with a high-dimensional state space (d_x large). One example is data assimilation for atmospheric and oceanographic science, where observations (satellites, weather stations, etc.) are used to inform large spatio-temporal simulation models: approximate filtering methods developed for such situations (4) could be used to apply the methods of this paper.

The goal of maximum likelihood estimation for partially observed data is reminiscent of the expectation–maximization (EM) algorithm (37), and indeed Monte Carlo EM methods have been applied to nonlinear state space models (24). The Monte Carlo EM algorithm, and other standard Monte Carlo Markov Chain methods, cannot be used for inference on the environmental noise parameter ε for the model given above, because these methods rely upon different sample paths of the unobserved process x_t having densities with respect to a common measure (38). Diffusion processes, such as the solution to the system of stochastic differential equations above, are mutually singular for different values of the infinitesimal variance. Modeling using diffusion processes (as in above) is by no means necessary for the application of Procedure 1, but continuous-time models for large discrete populations are well approximated by diffusion processes, so a method that can handle diffusion processes may be expected to be more reliable for large discrete populations.

Procedure 1 is well suited for maximizing numerically estimated likelihoods for complex models largely because it requires neither analytic derivatives, which may not be available, nor numerical derivatives, which may be unstable. The iterated filtering effectively produces estimates of the derivatives smoothed at each iteration over the scale at which the likelihood is currently being investigated. Although general stochastic optimization techniques do exist for maximizing functions measured with error (39), these methods are inefficient in terms of the number of function evaluations required (40). General stochastic optimization techniques have not to our knowledge been successfully applied to examples comparable to that presented here.

Each iteration of MIF requires similar computational effort to one evaluation of the likelihood function. The results in Fig. 3 demonstrate the ability of Procedure 1 to optimize a function of 13 variables using 50 function evaluations, with Monte Carlo measurement error and without knowledge of derivatives. This feat is only possible because Procedure 1 takes advantage of the state-space structure of the model; however, this structure is general enough to cover relevant dynamical models across a broad range of disciplines. The EM algorithm is similarly “only” an optimization trick, but in practice it has led to the consideration of models that would be otherwise intractable. The computational efficiency of Procedure 1 is essential for the model given above, where Monte Carlo function evaluations each take ≈15 min on a desktop computer.

Implementation of Procedure 1 using particle filtering conveniently requires little more than being able to simulate paths from the unobserved dynamical system. The new methodology is therefore readily adaptable to modifications of the model, allowing relatively rapid cycles of model development, model fitting, diagnostic analysis, and model improvement.

Theoretical Basis for MIF

Recall the notation above, and specifically the definitions in Eqs. 1 and 2.

Theorem 1.

Assuming conditions (R1–R3) below,

where ∇g is defined by [∇g]_i = ∂g/∂θ_i and θ̂₀ = θ. Furthermore, for a sequence σ_n → 0, define θ̂⁽ⁿ⁾ recursively by

where θ̂_t⁽ⁿ⁾ = θ̂_t(θ̂⁽ⁿ⁾, σ_n) and V_t,n = V_t(θ̂⁽ⁿ⁾, σ_n). If there is a θ̂ with |θ̂⁽ⁿ⁾ − θ̂|/σ_n² → 0 then∇ log f(y_1:T|θ = θ̂, σ = 0) = 0.

Theorem 1 asserts that (for sufficiently small σ_n), Procedure 1 iteratively updates the parameter estimate in the direction of increasing likelihood, with a fixed point at a local maximum of the likelihood. Step 2(ii) of Procedure 1 can be rewritten as θ̂⁽ⁿ⁺¹⁾ = V_1,n{Σ_t=1^T−1 (V _t,n⁻¹ − V _t+1,n⁻¹)θ̂_t⁽ⁿ⁾ + (V _T,n⁻¹)θ̂ _T⁽ⁿ⁾}. This makes θ̂⁽ⁿ⁺¹⁾ a weighted average, in the sense that V₁{Σ_t=1^T−1 (V _t⁻¹ − V _t+1⁻¹) + V _T⁻¹)} = I_dθ, where I_dθ is the _θ × d_θ identity matrix. The weights are necessarily positive for sufficiently small σ_n (Supporting Text).

The exponentially decaying σ_n in step 2(i) of Procedure 1 is justified by empirical demonstration, provided by convergence plots (Fig. 3). Slower decay, σ_n² = n^−β with 0 < β < 1, can give sufficient conditions for a Monte Carlo implementation of Procedure 1 to converge successfully (Supporting Text). In our experience, exponential decay yields equivalent results, considerably more rapidly. Analogously, simulated annealing provides an example of a widely used stochastic search algorithm where a geometric “cooling schedule” is often more effective than slower, theoretically motivated schedules (41).

In the proof of Theorem 1, we define f_t(ψ) = f(y_t|y_1:t−1, θ_t = ψ). The dependence on σ may be made explicit by writing f_t(ψ) = f_t(ψ, σ). We assume that y_1:T, c and Σ are fixed: for example, the constant B in R1 may depend on y_1:t. We use the Euclidean norm for vectors and the corresponding norm for matrices, i.e., |M| = sup_|u|≤1 |u′Mu|, where u′ denotes the transpose of u. We assume the following regularity conditions.

• R1. The Hessian matrix is bounded, i.e., there are constants B and σ₀ such that, for all σ < σ₀ and all θ_t ∈ ℝ^dθ, |∇²f_t(θ_t, σ)| < B.

• R2. E[|θ_t − θ̂_t−1|²| y_1:t−1] = O(σ²).

• R3. E[|θ_t−θ̂_t−1|³| y_1:t−1] = o(σ²).

R1 is a global bound over θ_t ∈ ℝ^dθ, comparable to global bounds used to show the consistency and asymptotic normality of the MLE (42, 43). It can break down, for example, when the likelihood is unbounded. This problem can be avoided by reparameterizing to keep the model away from such singularities, as is common practice in mixture modeling (44). R2–R3 require that a new observation cannot often have a large amount of new information about θ. For example, they are satisfied if θ_0:t, x_1:t, and y_1:t are jointly Gaussian. We conjecture that they are satisfied whenever the state space model is smoothly parametrized and the random walk θ_t does not have long tails.

Proof of Theorem 1.

Suppose inductively that |V_t| = O(σ²) and |θ̂_t−1 − θ| = O(σ²). This holds for t = 1 by construction. Bayes' formula gives

The numerator in Eq. 8 comes from a Taylor series expansion of f_t(θ̂_t), and R1 implies |R_t| ≤ B|θ_t − θ̂_t−1|²/2. The denominator then follows from applying this expansion to the integral in Eq. 7, invoking R2, and observing that Eq. 1 implies E[θ_t| y_1:t−1] = θ̂_t−1. We now calculate

Eq. 11 follows from Eq. 10 using Eq. 9 and R3. Eq. 12 follows from Eq. 11 using the induction assumptions on θ̂_t−1 and V_t; Eq. 12 then justifies this assumption for θ̂_t. A similar argument gives

where Eq. 14 follows from Eq. 13 via Eqs. 9 and 12 and the induction hypothesis on V_t. Eq. 14 in turn justifies this hypothesis. Summing Eq. 12 over t produces

which leads to Eq. 5. To see the second part of the theorem, note that Eq. 6 and the requirement that |θ̂⁽ⁿ⁾ − θ̂|/σ_n² → 0 imply that

Continuity then gives

which, together with Eq. 5, yields the required result.

Heuristics, Diagnostics, and Confidence Intervals

Our main MIF diagnostic is to plot parameter estimates as a function of MIF iteration; we call this a convergence plot. Convergence is indicated when the estimates reach a single stable limit from various starting points. Convergence plots were also used for simulations with a known true parameter, to validate the methodology. The investigation of quantitative convergence measures might lead to more refined implementations of Procedure 1.

Heuristically, α can be thought of as a “cooling” parameter, analogous to that used in simulated annealing (39). If α is too small, the convergence will be “quenched” and fail to locate a maximum. If α is too large, the algorithm will fail to converge in a reasonable time interval. A value of α = 0.95 was used above.

Supposing that θ_i has a plausible range [θ_i^lo, θ_i^hi] based on prior knowledge, then each particle is capable of exploring this range in early iterations of MIF (unconditional on the data) provided $\sqrt{{\Sigma }_{ii}T}$ is on the same scale as θ_i^hi − θ_i^lo. We use Σ_ii^1/2 = (θ_i^hi − θ_i^lo)/2$\sqrt{T}$ with Σ_ij = 0 for i ≠ j.

Although the asymptotic arguments do not depend on the particular value of the dimensionless constant c, looking at convergence plots led us to take c² = 20 above. Large values c² ≈ 40 resulted in increased algorithmic instability, as occasional large decreases in the prediction variance V_t resulted in large weights in Procedure 1 step 2(ii). Small values c² ≈ 10 were diagnosed to result in appreciably slower convergence. We found it useful, in choosing c, to check that [V_t]_ii plotted against t was fairly stable. In principle, a different value of c could be used for each dimension of θ; for our example, a single choice of c was found to be adequate.

If the dimension of θ is even moderately large (say, d_θ ≈ 10), it can be challenging to investigate the likelihood surface, to check that a good local maximum has been found, and to get an idea of the standard deviations and covariance of the estimators. A useful diagnostic, the “sliced likelihood” (Fig. 3B), plots ℓ(θ̂ + hδ_i) against θ̂_i + h, where δ_i is a vector of zeros with a one in the ith position. If θ̂ is located at a local maximum of each sliced likelihood, then θ̂ is a local maximum of ℓ(θ), supposing ℓ(θ) is continuously differentiable. Computing sliced likelihoods requires moderate computational effort, linear in the dimension of θ. A local quadratic fit is made to the sliced log likelihood (as suggested by ref. 35), because ℓ(θ̂ + hδ_i) is calculated with a Monte Carlo error. Calculating the sliced likelihood involves evaluating log f(y_t| y_1:t−1, θ̂ + hδ_i), which can then be regressed against h to estimate (∂/∂θ_i) log f(y_t| y_1:t−1, θ̂). These partial derivatives may then be used to estimate the Fisher information (ref. 34 and Supporting Text) and corresponding SEs. Profile likelihoods (34) can be calculated by using MIF, but at considerably more computational expense than sliced likelihoods. SEs and profile likelihood confidence intervals, based on asymptotic properties of MLEs, are particularly useful when alternate ways to find standard errors, such as bootstrap simulation from the fitted model, are prohibitively expensive to compute. Our experience, consistent with previous advice (45), is that SEs based on estimating Fisher information provide a computationally frugal method to get a reasonable idea of the scale of uncertainty, but profile likelihoods and associated likelihood based confidence intervals are more appropriate for drawing careful inferences.

As in regression, residual analysis is a key diagnostic tool for state space models. The standardized prediction residuals are {u_t(θ̂)} where θ̂ is the MLE and u_t(θ) = [Var(y_t| y_1:t−1, θ)]^−1/2 (y_t − E[y_t| y_1:t−1, θ]). Other residuals may be defined for state space models (8), such as E|∫ _t−1^t dW_s| y_1:T, θ̂] for the model in Eqs. 3 and 4. Prediction residuals have the property that, if the model is correctly specified with true parameter vector θ*, {u_t(θ*)} is an uncorrelated sequence. This has two useful consequences: it gives a direct diagnostic check of the model, i.e., {u_t(θ̂)} should be approximately uncorrelated; it means that prediction residuals are an (approximately) prewhitened version of the observation process, which makes them particularly suitable for using correlation techniques to look for relationships with other variables (7), as demonstrated above. In addition, the prediction residuals are relatively easy to calculate using particle-filter techniques (Supporting Text).

Supplementary Material