Hybrid modeling and prediction of dynamical systems
Franz Hamilton Alun Lloyd Kevin Flores







Hybrid Modeling and Prediction of Dynamical Systems 





Franz Hamilton1,2*,
Alun Lloyd 1,2,3,
Kevin Flores 1,2,4


1 Department of Mathematics, North Carolina State University, Raleigh, NC, USA


2 Center for Quantitative Sciences in Biomedicine, North Carolina State University, Raleigh, NC, USA


3 Biomathematics Graduate Program, North Carolina State University, Raleigh, NC, USA


4 Center for Research in Scientific Computation, North Carolina State University, Raleigh, NC, USA








* fwhamilt@ncsu.edu






§ ABSTRACT

Scientific analysis often relies on the ability to make accurate predictions of a system's dynamics. Mechanistic models, parameterized by a number of unknown parameters, are often used for this purpose. Accurate estimation of the model state and parameters prior to prediction is necessary, but may be complicated by issues such as noisy data and uncertainty in parameters and initial conditions. At the other end of the spectrum exist nonparametric methods, which rely solely on data to build their predictions. While these nonparametric methods do not require a model of the system, their performance is strongly influenced by the amount and noisiness of the data. In this article, we consider a hybrid approach to modeling and prediction which merges recent advancements in nonparametric analysis with standard parametric methods. The general idea is to replace a subset of a mechanistic model's equations with their corresponding nonparametric representations, resulting in a hybrid modeling and prediction scheme. Overall, we find that this hybrid approach allows for more robust parameter estimation and improved short-term prediction in situations where there is a large uncertainty in model parameters. We demonstrate these advantages in the classical Lorenz-63 chaotic system and in networks of Hindmarsh-Rose neurons before application to experimentally collected structured population data. 







§ AUTHOR SUMMARY

The question of how best to predict the evolution of a dynamical system has received substantial interest in the scientific community. While traditional mechanistic modeling approaches have dominated, data-driven approaches which rely on data to build predictive models have gained increasing popularity. The reality is, both approaches have their drawbacks and limitations. In this article we ask the question of whether or not a hybrid approach to prediction, which combines characteristics of both mechanistic modeling and data-driven modeling, can offer improvements over the standalone methodologies. We analyze the performance of these methods in two model systems and then evaluate them on experimentally collected population data.






§ INTRODUCTION

Parametric modeling involves defining an underlying set of mechanistic equations which describe a system's dynamics. These mechanistic models often contain a number of unknown parameters as well as an uncertain state, both of which need to be quantified prior to use of the model for prediction. The success of parametric prediction is tied closely to the ability to construct accurate estimates of the model parameters and state. This can be particularly challenging in high dimensional estimation problems as well as in chaotic systems <cit.>. Additionally, there is often a degree of model error, or a discrepancy between the structure of the model and that of the system, further complicating the estimation process and hindering prediction accuracy.

Despite these potential issues, mechanistic models are frequently utilized in data analysis. The question we aim to address is when is it advantageous to use them? Under suitable conditions where model error is relatively small and parameters can be reliably estimated, parametric predictions can provide a great deal of accuracy. However, as we will see in the subsequent examples, a large uncertainty in the initial parameter values often leads to inaccurate estimates resulting in poor model-based predictions.

An alternative approach to modeling and prediction abandons the use of any mechanistic equations, instead relying on predictive models built from data. These nonparametric methods have received considerable attention, in particular those methods based on Takens' delay-coordinate method for attractor reconstruction <cit.>.  The success of nonparametric methods is strongly influenced by the amount of data available as well as the dimension of the dynamical system. If only a sparse amount of training data is available, the result is often inaccurate predictions due to the lack of suitable nearby neighbors in delay-coordinate space. Furthermore, as the dimension and complexity of the dynamical system increases, nonparametric prediction becomes significantly more difficult due to the necessary data requirements <cit.>.

Several recent works have investigated the situation where only a portion of a mechanistic model is known <cit.>. Our motivation here though is to explore how best to use a full mechanistic model when it is available. We consider a hybrid methodology to modeling and prediction that combines the complementary features of both parametric and nonparametric methods. In our proposed hybrid method, a subset of a mechanistic model's equations are replaced by nonparametric evolution. These nonparametrically advanced variables are then incorporated into the remaining mechanistic equations during the data fitting and prediction process. The result of this approach is a more robust estimation of model parameters as well as an improvement in short-term prediction accuracy when initial parameter uncertainty is large.

The utility of this method is demonstrated in several example systems. The assumption throughout is that noisy training data from a system are available as well as a mechanistic model that describes the underlying dynamics. However, several of the model parameters are unknown and the model state is uncertain due to the noisy measurements. The goal is to make accurate predictions of the system state up to some forecast horizon beyond the end of the training data. We compare the prediction accuracy of the standard parametric and nonparametric methodologies with the novel hybrid method presented here.

We begin our analysis by examining prediction in the classical Lorenz-63 system <cit.>, which exhibits chaotic dynamics. Motivated by the success of the hybrid method in the Lorenz-63 system, we consider a more sophisticated example of predicting the spiking dynamics of a neuron in a network of Hindmarsh-Rose <cit.> cells. Finally, we examine the prediction problem in a well-known experimental dataset from beetle population dynamics <cit.>.



§ MATERIALS AND METHODS

The assumption throughout is that a set of noisy data is available over the time interval [t(0),t(T)]. This is referred to as the training data of the system. Using these training data, the question is how best to predict the system dynamics over the interval [t(T+1),t(T+T_F)], known as the prediction interval. Standard parametric and nonparametric methods are presented before our discussion of the novel hybrid method which blends the two approaches. 



 §.§ Parametric Modeling and Prediction

When a full set of mechanistic equations is used for modeling and prediction, we refer to this as the parametric approach. Assume a general nonlinear system of the form

    𝐱(k+1)    =   𝐟(t(k),𝐱(k),𝐩)+𝐰(k)
    𝐲(k)    =   𝐡(t(k),𝐱(k),𝐩)+𝐯(k)

where x= [x_1,x_2,,x_n]^ T is an n-dimensional vector of model state variables and p = [p_1,p_2,,p_l]^ T is an l-dimensional vector of model parameters which may be known from first principles, partially known or completely unknown. f represents our system dynamics which describe the evolution of the state 𝐱 over time and h is an observation function which maps x to an m-dimensional vector of model observations, y = [y_1,y_2,,y_m]^ T. To simplify the description of our analysis, we assume that the training data maps directly to some subset of 𝐱. w(k) and v(k) are assumed to be mean 0 Gaussian noise terms with covariances 𝐐 and 𝐑 respectively. While discrete notation is used in Eq. <ref> for notational convenience,  the evolution of x is often described by continuous-time systems. In this situation numerical solvers, such as Runge-Kutta or Adams-Moulton methods, are used to obtain solutions to the continuous-time system at discrete time points.

When the state of a system is uncertain due to noisy or incomplete observations, nonlinear Kalman filtering can be used for state estimation <cit.>. Here we choose the unscented Kalman filter (UKF), which approximates the propagation of the mean and covariance of a random variable through a nonlinear function using a deterministic ensemble selected through the unscented transformation <cit.>. We initialize the filter with state vector 𝐱^+(0) and covariance matrix 𝐏^+(0). At the kth step of the filter there is an estimate of the state 𝐱^+(k-1) and the covariance matrix 𝐏^+(k-1). In the UKF, the singular value decomposition is used to find the square root of the matrix 𝐏^+(k-1), which is used to form an ensemble of 2n+1 state vectors.

The model 𝐟 is applied to the ensemble, advancing it forward one time step, and then observed with 𝐡. The weighted average of the resulting state ensemble gives the prior state estimate 𝐱^-(k) and the weighted average of the observed ensemble is the model-predicted observation 𝐲^-(k).  Covariance matrices 𝐏^-(k) and 𝐏^𝐲(k) of the resulting state and observed ensemble, and the cross-covariance matrix 𝐏^𝐱𝐲(k) between the state and observed ensembles,  are formed and the equations

    𝐊(k)    =   𝐏^𝐱𝐲(k)(𝐏^𝐲(k))^-1
    𝐏^+(k)    =   𝐏^-(k)-𝐏^xy(k)(𝐏^y(k))^-1𝐏^yx(k)
    𝐱^+(k)    =   𝐱^-(k)+𝐊(k)(𝐲(k)-𝐲^-(k) ).

are used to update the state and covariance estimates with the observation 𝐲(k).

The UKF algorithm described above can be extended to include the joint estimation problem allowing for parameter estimation. In this framework, the parameters 𝐩 are considered as auxiliary state variables with trivial dynamics, namely 𝐩_k+1 = 𝐩_k. An augmented n+l dimensional state vector can then be formed consisting of the original n state variables and l model parameters allowing for simultaneous state and parameter estimation <cit.>.

To implement parametric prediction, the UKF is used to process the training data and obtain an estimate of 𝐩, as well as the state at the end of the training set, 𝐱(T). The parameter values are fixed and Eq. <ref> is forward solved from t(T) to generate predictions of the system dynamics over the prediction interval [t(T+1),t(T+T_F)]. Namely, predictions x(T+1),x(T+2),,x(T+T_F) are calculated.



 §.§ Takens' Method for Nonparametric Prediction

Instead of using the mechanistic model described by Eq. <ref>, the system can be represented nonparametrically. Without loss of generality consider the observed variable x_j. Using Takens' theorem <cit.>, the d+1 dimensional delay coordinate vector x_j^d(T) = [x_j(T),x_j(T-τ),x_j(T-2τ), x_j(T-dτ)] is formed which represents the state of the system at time t(T). Here d is the number of delays and τ is the time-delay. 

The goal of nonparametric prediction is to utilize the training data in the interval [t(0),t(T) ] to build local models for predicting the dynamics over the interval [t(T+1),t(T+T_F)]. Here, the method of direct prediction is chosen. Prior to implementation of the direct prediction, a library of delay vectors is formed from the training data of x_j. 

Direct prediction begins by finding the κ nearest neighbors, as a function of Euclidean distance, to the current delay-coordinate vector x_j^d(T) within the library of delay vectors. Neighboring delay vectors

    x_j^d(T')    =   [x_j(T'),x_j(T'-τ),x_j(T'-2τ), x_j(T'-dτ)]
    
    x_j^d(T”)    =   [x_j(T”),x_j(T”-τ),x_j(T”-2τ), x_j(T”-dτ)]
       ⋮   
    
    x_j^d(T^κ)    =   [x_j(T^κ),x_j(T^κ-τ),x_j(T^κ-2τ), x_j(T^κ-dτ)]

are found within the training data and the known x_j(T'+i), x_j(T”+i), …, x_j(T^κ+i) points are used in a local model to predict the unknown value x_j(T+i) where i = 1,2,,T_F. In this article, a locally constant model is chosen

    x_j(T+i) ≈ w_j'x_j(T'+i) + w_j”x_j(T”+i) +  + w_j^κx_j(T^κ+i)

where w_j',w_j”,,w_j^κ are the weights for the j^th state that determine the contribution of each neighbor in building the prediction. In its simplest form, Eq. <ref> is an average of the nearest neighbors where w_j' = w_j” =  = w_j^κ = 1/κ. More sophisticated weighting schemes can be chosen, for example assigning the weights based on the Euclidean distance from each neighbor to the current delay vector <cit.>. Selection of values for d, τ and κ is necessary for implementation of the direct prediction algorithm. These values were optimized, within each example, to give the lowest prediction error (results not shown).

The accuracy of the predicted x_j(T+i) is subject to several factors.  The presence of noise in the training data plays a substantial role in decreasing prediction accuracy. However, recent advancements in nonparametric analysis have addressed the problem of filtering time series without use of a mechanistic model. In <cit.>, a nonparametric filter was developed which merged Kalman filtering theory and Takens' method. The resulting Kalman-Takens filter was demonstrated to be able to reduce significant amounts of noise in data. Application of the method was extended in <cit.> to the case of filtering stochastic variables without a model. In the results presented below, the training data used for nonparametric prediction are filtered first using the method of <cit.>.



 §.§ Hybrid Modeling and Prediction: Merging Parametric and Nonparametric Methods


As an alternative to the parametric and nonparametric methods described above, we propose a hybrid approach which blends the two methods together. In this framework, we assume that a full mechanistic model as described by Eq. <ref> is available. However, rather than using the full model, a subset of the mechanistic equations are used and the remainder of the variables are represented nonparametrically using delay-coordinates.

In formulating this method it is convenient to first think of Eq. <ref> without vector notation

    x_1(k+1)    =    f_1(t(k),x_1(k),x_2(k),,x_n(k),p_1,p_2,,p_l)
    
    x_2(k+1)    =    f_2(t(k),x_1(k),x_2(k),,x_n(k),p_1,p_2,,p_l)
       ⋮   
    
    x_n(k+1)    =    f_n(t(k),x_1(k),x_2(k),,x_n(k),p_1,p_2,,p_l)

Now assume only the first n-1 equations of Eq. <ref> are used to model state variables x_1,x_2,…,x_n-1, while x_n is described nonparametrically

    x_1(k+1)    =    f_1(t(k),x_1(k),x_2(k),,x_n-1(k),x_n(k),p_1,p_2,,p_l)
    
    x_2(k+1)    =    f_2(t(k),x_1(k),x_2(k),,x_n-1(k),x_n(k),p_1,p_2,,p_l)
       ⋮   
    
    x_n-1(k+1)    =    f_n-1(t(k),x_1(k),x_2(k),,x_n-1(k),x_n(k),p_1,p_2,,p_l) 
    
    x_n(k+1)    ≈    w_n'x̃_n(T'+k+1) + w_n”x̃_n(T”+k+1) +  + w_n^κx̃_n(T^κ+k+1)

We refer to Eq. <ref> as the hybrid model. Note, in Eq. <ref> only x_n is assumed to be advanced nonparametrically. This is done purely for ease of presentation and the hybrid model can instead contain several variables whose equations are replaced by nonparametric advancement.

The hybrid model has several distinguishing features. Notice, in this framework nonparametrically advanced dynamics are incorporated into mechanistic equations, essentially merging the two lines of mathematical thought. Furthermore, equations for state variables within Eq. <ref> can be replaced only if there are observations which map directly to them, otherwise their dynamics can not be nonparametrically advanced. Finally, the process of replacing equations in the hybrid method will generally result in a reduction in the number of unknown model parameters to be estimated.

In this hybrid scheme, obtaining an estimate of the unknown parameters in the n-1 mechanistic equations and an estimate of x(T) requires a combination of the nonparametric analysis developed in <cit.> and traditional parametric methodology. The state variable x_n, which is not defined by a mechanistic equation in Eq. <ref>, is represented by delay coordinates within the UKF. Therefore at step k we have the hybrid state

    𝐱^ H(k) = [x_1(k),x_2(k) ,…, x_n-1(k),x_n(k), x_n(k-τ) , x_n(k-2τ),… , x_n(k-dτ)]^ T

The UKF as described above is implemented with this hybrid state 𝐱^ H(k) and the model described by Eq. <ref>. Notice that in the case of the hybrid model when we have to advance the state dynamics and form the prior estimate in the UKF, the advancement is done parametrically for the first n-1 states and nonparametrically for the n^th state. Similarly to before, we can augment 𝐱^ H with the unknown parameters in the n-1 mechanistic equations allowing for simultaneous parameter estimation.

Once the training data are processed and an estimate of 𝐱^ H(T) and the parameters are obtained, the hybrid model in Eq. <ref> is implemented to generate predictions 𝐱^ H(T+1), 𝐱^ H(T+2),…, 𝐱^ H(T+T_F).








§ RESULTS

We demonstrate the utility of the hybrid methodology, with comparison to standard parametric and nonparametric modeling and prediction, in the following example systems. When conducting this analysis, two types of error are considered. The first, error in the observations, manifests itself as noise in the training data which all three methods will have to confront. The second type, error in the parameters, takes the form of an uncertainty in the initial parameter values used by the UKF for parameter estimation. Only the parametric and hybrid methods will have to deal with this parameter error. Throughout, we will refer to a percentage uncertainty which corresponds to the standard deviation of the distribution from which the initial parameter value is drawn relative to the mean. For example, if the true value for a parameter p_1 is 12 and we have 50% uncertainty in this value, then the initial parameter value used for estimating p_1 will be drawn from the distribution N(12,(0.5*12)^2).

To quantify prediction accuracy, the normalized root-mean-square-error, or SRMSE, is calculated for each prediction method as a function of forecast horizon. Normalization is done with respect to the standard deviation of the variable as calculated from the training data. In using the SRMSE metric, the goal is to be more accurate than if the prediction was simply the mean of the training data (corresponding to SRMSE = 1). Thus a prediction is better than a naive prediction when SRMSE < 1, though for chaotic systems prediction accuracy will eventually converge to this error level since only short-term prediction is possible.





 §.§ Prediction in the Lorenz-63 System


As a demonstrative example, consider the Lorenz-63 system <cit.>

    ẋ   =   σ(y-x)
    ẏ   =    x(ρ-z)-y
    ż   =    xy-β z

where σ = 10, ρ = 28, β = 8/3. Data are generated from this system using a fourth-order Adams-Moulton method with sample rate h = 0.05. We assume that 500 training data points of the x, y and z variables are available, or 25 units of time. The Lorenz-63 system oscillates approximately once every unit of time, meaning the training set consists of about 25 oscillations. The goal is to accurately predict the dynamics of x, y and z one time unit after the end of the training set. However,  the observations of each variable are corrupted by Gaussian observational noise with mean zero and variance equal to 4. Additionally the true value of parameters σ, ρ and β are unknown. Fig. <ref> shows an example simulation of this system. 


The parametric method utilizes Eq. <ref> to estimate the model state and parameters, and to predict the x, y and z dynamics. For the nonparametric method, delay coordinates of the variables are formed with d = 9 and τ = 1. The local constant model for prediction is built using κ = 20 nearest neighbors. For the hybrid method, the mechanistic equation governing the dynamics of y are replaced nonparametrically resulting in the reduced Lorenz-63 model

    ẋ   =   σ(y-x)
    ż   =    xy-β z

Note, the hybrid model does not require estimation of the ρ parameter since the mechanistic equation for y is removed.




Fig. <ref> shows a comparison of parametric (black), nonparametric (blue) and hybrid (red) prediction error as a function of forecast horizon. SRMSE results averaged over 500 system realizations. Various parameter uncertainty levels are shown: 80% uncertainty (solid lines), 50% uncertainty (dashed-dotted lines) and 20% uncertainty (dashed line). The hybrid method with 80% uncertainty offers improved short-term prediction of the Lorenz-63 x (Fig. <ref>a) and z (Fig. <ref>c) variables over standalone nonparametric prediction as well as parametric prediction with 80% uncertainty. Hybrid and nonparametric prediction of y (Fig. <ref>b) are comparable, which is to be expected since the hybrid approach is using nonparametric advancement of y in its formulation. Note that parametric prediction at this uncertainty level does very poorly and in the cases of y and z its result is not shown due to the scale of the error. As the uncertainty decreases for parametric prediction, its performance improves. However, hybrid prediction with 80% uncertainty still outperforms parametric prediction with 50% uncertainty in the short-term. At a small uncertainty level, parametric prediction outperforms both hybrid and nonparametric methods which is to be expected since it has access to the true model equations and starts out with close to optimal parameter values.

The success of the hybrid method at higher uncertainty levels can be traced to more accurate estimates of the model parameters in the mechanistic equations that it uses. Table <ref> shows the resulting hybrid and parametric estimation of the Lorenz-63 parameters. The hybrid method with 80% uncertainty is able to construct accurate estimates of both σ and β, with a mean close to the true value and a small standard deviation of the estimates. The parametric method with 80% and 50% uncertainty is unable to obtain reliable estimates, exemplified by the large standard deviation of the estimates. Only when the parametric method has a relatively small uncertainty of 20% is it able to accurately estimate the system parameters.






 §.§ Predicting Neuronal Network Dynamics

We now consider the difficult high dimensional estimation and prediction problem posed by neuronal network studies. If we are only interested in predicting a portion of the network, then we can use the proposed hybrid method to refine our estimation and prediction while simultaneously reducing estimation complexity. As an example
of this potential network application we consider the prediction of spiking dynamics in a network of M Hindmarsh-Rose neurons <cit.>

    ẋ_i    =    y_i-a_ix_i^3+b_ix_i^2-z_i+1.2+∑_i≠ m^M β_im/1+9e^-10x_mx_m 
    ẏ_i    =    1-c_ix_i^2 
    ż_i    =    5× 10^-5[4(x_i-(-8/5))-z_i ]
 
where i = 1,2,,M. x_i corresponds to the spiking potential while y_i and z_i describe the fast and slow-scale dynamics, respectively, of neuron i. Each individual neuron in the network has parameters a_i =1, b_i = 3 and c_i = 5 which are assumed to be unknown. β_im represents the connectivity coefficient from neuron i to neuron m. For a network of size M, we have M^2-M possible connection parameters since neuron self connections are not allowed (i.e. β_ii = 0). These connection parameters are also assumed to be unknown.






For this example we examine networks of size M = 3 with 5 random connections. Data from these networks are generated using a fourth-order Adams-Moulton method with sample rate h = 0.08 ms.  We assume that the training data consists of 3000 observations, or 240 ms, of the x_1, x_2,x_3 variables each of which are corrupted by Gaussian noise with mean 0 and variance of 0.2. Under the stated parameter regime, the neurons in the network spike approximately every 6 ms, meaning our training set has on average around 40 spikes per neuron. In this example, we restrict our focus to predicting 8 ms of the x_3 variable (though a similar analysis follows for the prediction of x_1 and x_2). Fig. <ref>a shows a representative realization of this problem. Given our interest in x_3, the hybrid method only assumes a mechanistic equation for neuron 3

    ẋ_3    =    y_3-a_3x_3^3+b_3x_3^2-z_3+1.2+∑_3 ≠ m^M β_3m/1+9e^-10x_mx_m 
    ẏ_3    =    1-c_3x_3^2
    ż_3    =    5× 10^-5[4(x_3-(-8/5))-z_3 ]

and nonparametrically represents neuron 1 and neuron 2.

Fig. <ref>b shows the resulting accuracy in predicting x_3 when using parametric (black), nonparametric (blue) and hybrid (red) methods with 80% (solid line) and 50% (dashed-dotted line) uncertainty in parameter values. The parametric approach uses the full mechanistic model described by Eq. <ref> for modeling and prediction, requiring estimation of the x,y and z state variables and parameters a,b and c for each neuron, as well as the full connectivity matrix. Notice that once again with 80% uncertainty, the scale of error for the parametric method is much larger compared to the other methods. Only with 50% uncertainty is the parametric method able to provide reliable predictions of x_3. Note that unlike in the Lorenz-63 example, we do not consider the parametric method with 20% uncertainty since reasonable parameter estimates and predictions are obtained with 50% uncertainty. The nonparametric method (τ = 1, d = 9) uses κ = 10 neighbors for building the local model for prediction. Again we observe that the hybrid method, even with a large parameter uncertainty of 80%, provides accurate predictions of x_3 compared to the other methods. Table <ref> shows the robustness of the hybrid method in estimating the individual parameters for neuron 3.










 §.§ Predicting Flour Beetle Population Dynamics


We now investigate the prediction problem in a well-known data set from an ecological study involving the cannibalistic red flour beetle Tribolium castaneum. In <cit.>, the authors present experimentally collected data and a mechanistic model describing the life cycle dynamics of T. castaneum. Their discrete time model describing the progression of the beetle through the larvae, pupae, and adult stages is given by 

    L(t+1)    =    bA(t) e^-c_elL(t) - c_eaA(t)
    
    P(t+1)    =    L(t)(1-μ_l)
    
    A(t+1)    =    P(t) e^-c_paA(t)+A(t)(1-μ_a)

where L, P and A correspond to larvae, pupae and adult populations, respectively. The essential interactions described by this model are (i) flour beetles become reproductive only in the adult stage, (ii) adults produce new larvae, (iii) adults and larvae can both cannibalize larvae, and (iv) adults cannibalize pupae. We note that since Eq. <ref> only approximates the life cycle dynamics of the beetle, there is a degree of model error in the proposed system, unlike the previous examples.




The authors of <cit.> experimentally set the adult mortality rate (μ_a) to 0.96 and the recruitment rate (c_pa) from pupae to adult to seven different values (0, 0.05, 0.10, 0.25, 0.35, 0.50, 1.0). Experiments at each recruitment rate value were replicated three times resulting in 21 different datasets. Each dataset consists of total numbers of larvae, pupae, and adults measured bi-weekly over 82 weeks resulting in 41 measurements for each life stage. These data were fit to Eq. <ref> in <cit.> and parameter estimates b = 6.598, c_el = 1.209 × 10^-2, c_ea = 1.155 × 10^-2 and μ_l = 0.2055 were obtained. We treat these parameter values as ground truth when considering the different parameter uncertainty levels for fitting the data to the model.

In our analysis of this system, we treat the first 37 measurements (or 74 weeks) within an experiment as training data and use the remaining 4 time points (or 8 weeks) for forecast evaluation. Fig. <ref> shows an example of this setup for a representative dataset. Fig. <ref> shows the results of predicting the larvae (Fig. <ref>a), pupae (Fig. <ref>b) and adult (Fig. <ref>c) populations using parametric (black), nonparametric (blue) and hybrid prediction methods with 80% (solid line) and 50% (dashed-dotted line) parameter uncertainty levels. Error bars correspond to the standard error over the 21 datasets. The parametric method uses the full mechanistic model described in Eq. <ref> to estimate the population state and parameters b, c_el, c_ea and μ_l before prediction. We note in Fig. <ref> that the parametric method with 80% uncertainty is not shown due to the scale of the error, and is significantly outperformed by the nonparametric prediction (τ = 1, d = 2, κ = 5). For the hybrid method, we only consider the mechanistic equations for pupae and adult population dynamics

    P(t+1)    =    L(t)(1-μ_1)
    
    A(t+1)    =    P(t) e^-c_paA(t)+A(t)(1-μ_a)

and nonparametrically represent larvae. Hybrid prediction with 80% uncertainty outperforms both nonparametric and parametric with 80% uncertainty for pupae and adult population levels, and is comparable to parametric with 50% uncertainty.








§ CONCLUSION


By blending characteristics of parametric and nonparametric methodologies, the proposed hybrid method for modeling and prediction offers several advantages over standalone methods. From the perspective of model fitting and the required parameter estimation that arises in this process, we have shown that the hybrid approach allows for a more robust estimation of model parameters. Particularly for situations where there is a large uncertainty in the true parameter values, the hybrid method is able to construct accurate estimates of model parameters when the standard parametric model fitting fails to do so. At first this may seem counter-intuitive, but in fact it is not that surprising. The replacement of mechanistic equations with their nonparametric representations in effect reduces the dimension of the parameter space that we have to optimize in, resulting in better parameter estimates. As we have demonstrated in the above examples, this refinement in the parameter estimates leads to an improvement in short-term prediction accuracy.

The limitations of the hybrid method are similar to those of parametric and nonparametric methods in that if not enough training data are available then accurate estimation and prediction becomes difficult. However, the demonstrated robustness of the hybrid method to large parameter uncertainty is encouraging, particularly when considering experimental situations where we may not have a good prior estimate of the model parameters. One could consider implementing the hybrid method in an iterative fashion, estimating the parameters of each equation separately, then piecing the model back together for prediction. We can think of this as an iterative hybrid method, and is the subject of future work.

We view this work as complementary to recent publications on forecasting <cit.>. The authors of <cit.> advocate nonparametric methods over parametric methods in general, while a letter <cit.> addressing the work of <cit.> showed that a more sophisticated method for model fitting results in better parameter estimates and therefore model-based predictions which outperform model-free methods. Our results support the view that no one method is uniformly better than the other. As we showed in the above examples, in situations where the model error and uncertainty in initial parameters are relatively small, the parametric approach outperforms other prediction methods. Often in experimental studies though, we are not operating in this ideal situation and instead are working with a model that has substantial error with a large uncertainty in parameters which can lead to inaccurate system inference. In situations such as these, nonparametric methods are particularly useful.

The main appeal of the hybrid method is that we can confront these situations without having to completely abandon the use of the mechanistic equations. This is important since mechanistic models often provide valuable information about the underlying processes governing the system dynamics. While we explored in detail the robustness of the hybrid method to large levels of parameter uncertainty, its usefulness stretches well beyond that. In some instances, we may only have a model for some of the states or portions of the model may have higher error than others. By supplementing these parts with their nonparametric representation, the hybrid method would allow us to only use the parts of the model we are confident in and thus improve our analysis.








10

voss
Voss H, Timmer J, Kurths J.
Nonlinear dynamical system identification from uncertain and indirect
  measurements.
Int J Bif Chaos. 2002;14:1905–1924.

baake
Baake E, Baake M, Bock H, Briggs K.
Fitting ordinary differential equations to chaotic data.
Physical Review A. 1992;45:5524–5529.

farmer
Farmer J, Sidorowich J.
Predicting chaotic time series.
Phys Rev Lett. 1987;59:845–848.

casdagli1989nonlinear
Casdagli M.
Nonlinear prediction of chaotic time series.
Physica D: Nonlinear Phenomena. 1989;35(3):335–356.

Sugihara:1990aa
Sugihara G, May RM.
Nonlinear forecasting as a way of distinguishing chaos from
  measurement error in time series.
Nature. 1990;344(6268):734–741.

smith1992identification
Smith LA.
Identification and prediction of low dimensional dynamics.
Physica D: Nonlinear Phenomena. 1992;58(1):50–76.

jimenez1992forecasting
Jimenez J, Moreno J, Ruggeri G.
Forecasting on chaotic time series: A local optimal
  linear-reconstruction method.
Phys Rev A. 1992;45(6):3553.

sauer94
Sauer T.
Time series prediction by using delay coordinate embedding.
In: Time Series Prediction: Forecasting the Future and Understanding
  the Past. Addison Wesley; 1994. p. 175–193.

sugihara1994nonlinear
Sugihara G.
Nonlinear forecasting for the classification of natural time series.
Philosophical Transactions of the Royal Society of London Series A:
  Physical and Engineering Sciences. 1994;348(1688):477–495.

schroer1998predicting
Schroer CG, Sauer T, Ott E, Yorke JA.
Predicting chaos most of the time from embeddings with
  self-intersections.
Phys Rev Lett. 1998;80(7):1410.

kugiumtzis1998regularized
Kugiumtzis D, Lingjærde O, Christophersen N.
Regularized local linear prediction of chaotic time series.
Physica D: Nonlinear Phenomena. 1998;112(3):344–360.

yuan
Yuan G, Lozier M, Pratt L, Jones C, Helfrich K.
Estimating the predicability of an oceanic time series using linear
  and nonlinear methods.
J Geophys Res. 2004;109:C08002.

hsieh2005distinguishing
Hsieh CH, Glaser SM, Lucas AJ, Sugihara G.
Distinguishing random environmental fluctuations from ecological
  catastrophes for the North Pacific Ocean.
Nature. 2005;435(7040):336–340.

strelioff2006medium
Strelioff CC, Hübler AW.
Medium-term prediction of chaos.
Phys Rev Lett. 2006;96(4):044101.

regonda
Regonda S, Rajagopalan B, Lall U, Clark M, Moon YI.
Local polynomial method for ensemble forecast of time series.
Nonlin Proc in Geophys. 2005;12:397–406.

schelter2006handbook
Schelter B, Winterhalder M, Timmer J.
Handbook of time series analysis: recent theoretical developments and
  applications.
John Wiley and Sons; 2006.

hamilton2016
Hamilton F, Berry T, Sauer T.
Ensemble Kalman filtering without a model.
Physcial Review X. 2016;6:011021.

hamilton2
Hamilton F, Berry T, Sauer T.
Predicting chaotic time series with a partial model.
Physical Review E. 2015;92(1):010902.

berry2016
Berry T, Harlim J.
Semiparametric forecasting and filtering: correcting low-dimensional
  model error in parametric models.
Journal of Computational Physics. 2016;308:305–321.

lorenz63
Lorenz E.
Deterministic nonperiodic flow.
J Atmos Sci. 1963;20:130–141.

hindmarsh
Hindmarsh J, Rose R.
A model of neuronal bursting using three coupled first order
  differential equations.
Proc Roy Soc. 1984;221:87–102.

constantino
Constantino RF, Desharnais RA, Cushing JM, Dennis B.
Chaotic dynamics in an insect population.
Science. 1997;276:1881–1882.

enkf7
Kalnay E.
Atmospheric modeling, data assimilation, and predictability.
Cambridge Univ. Press; 2003.

evensen
Evensen G.
Data assimilation: The Ensemble Kalman Filter.
Springer: Heidelberg; 2009.

rabier
Rabier F.
Overview of global data assimilation developments in numerical
  weather-prediction centres.
Quarterly Journal of the Royal Meteorological Society.
  2005;131(613):3215–3233.

cummings
Cummings JA.
Operational multivariate ocean data assimilation.
Quarterly Journal of the Royal Meteorological Society.
  2005;131(613):3583–3604.
doi:10.1256/qj.05.105.

yoshida
Yoshida K, Yamaguchi J, Kaneda Y.
Regeneration of Small Eddies by Data Assimilation in Turbulence.
Phys Rev Lett. 2005;94:014501.

stuart
Law K, Stuart A.
Evaluating data stimulation algorithms.
Mon Wea Rev. 2012;140:3757–3782.

schiffbook
Schiff SJ.
Neural control engineering.
MIT Press; 2012.

berry2
Berry T, Sauer T.
Adaptive ensemble Kalman filtering of nonlinear systems.
Tellus A. 2013;65:20331.

hamiltonEPL
Hamilton F, Cressman J, Peixoto N, Sauer T.
Reconstructing neural dynamics using data assimilation with multiple
  models.
Europhysics Letters. 2014;107:68005.

hamiltonPRE
Hamilton F, Berry T, Peixoto N, Sauer T.
Real-time tracking of neuronal network structure using data
  assimilation.
Physical Review E. 2013;88:052715.

ghanim
Ullah G, Schiff S.
Tracking and control of neuronal Hodgkin-Huxley dynamics.
Phys Rev E. 2009;79:040901.

ghanim2
Ullah G, Schiff S.
Assimilating seizure dynamics.
PLoS Computational Biology. 2010;6:e1000776.

sitz2002
Sitz A, Schwarz U, Kurths J, Voss H.
Estimation of parameters and unobserved components for nonlinear
  systems from noisy time series.
Physical Review E. 2002;66:16210.

simon
Simon D.
Optimal State Estimation: Kalman, H_∞, and Nonlinear
  Approaches.
John Wiley and Sons; 2006.

julier1
Julier S, Uhlmann J, Durrant-Whyte H.
A new method for the nonlinear transformation of means and
  covariances in filters and estimators.
IEEE Trans Automat Control. 2000;45:477–482.

julier2
Julier S, Uhlmann J, Durrant-Whyte H.
Unscented filtering and nonlinear estimation.
Proc IEEE. 2004;92:401–422.

takens
Takens F.
Detecting strange attractors in turbulence.
Lecture Notes in Math Springer-Verlag: Berlin. 1981;898:366–381.

SYC
Sauer T, Yorke JA, Casdagli M.
Embedology.
J Stat Phys. 1991;65:579–616.

perretti
Perretti C, Munch S, Sugihara G.
Model-free forecasting outperforms the correct mechanistic model for
  simulated and experimental data.
Proceedings of the National Academy of Sciences. 2013;110:5253–5257.

perretti2
Perretti C, Sugihara G, Munch S.
Nonparametric forecasting outperforms parametric methods for a
  simulated multispecies system.
Ecology. 2013;94:794–800.

hamiltonEPJ
Hamilton F, Berry T, Sauer T.
Kalman-Takens filtering in the presence of dynamical noise.
Submitted to European Physical Journal;.

hartig
Hartig F, Dormann C.
Does model-free forecasting really outperform the true model?
Proceedings of the National Academy of Sciences. 2013;110:E3975.






