Fast Stable Parameter Estimation for Linear Dynamical Systems

Dynamical systems describe the changes in processes that arise naturally from their underlying physical principles, such as the laws of motion or the conservation of mass, energy or momentum. These models facilitate a causal explanation for the drivers and impediments of the processes. But do they describe the behaviour of the observed data? And how can we quantify the models' parameters that cannot be measured directly? This paper addresses these two questions by providing a methodology for estimating the solution; and the parameters of linear dynamical systems from incomplete and noisy observations of the processes. The proposed procedure builds on the parameter cascading approach, where a linear combination of basis functions approximates the implicitly defined solution of the dynamical system. The systems' parameters are then estimated so that this approximating solution adheres to the data. By taking advantage of the linearity of the system, we have simplified the parameter cascading estimation procedure, and by developing a new iterative scheme, we achieve fast and stable computation. We illustrate our approach by obtaining a linear differential equation that represents real data from biomechanics. Comparing our approach with popular methods for estimating the parameters of linear dynamical systems, namely, the non-linear least-squares approach, simulated annealing, parameter cascading and smooth functional tempering reveals a considerable reduction in computation and an improved bias and sampling variance.

The proposed procedure builds on the parameter cascading approach, where a linear combination of basis functions approximates the implicitly defined solution of the dynamical system.The systems' parameters are then estimated so that this approximating solution adheres to the data.By taking advantage of the linearity of the system, we have simplified the parameter cascading estimation procedure, and by developing a new iterative scheme, we achieve fast and stable computation.
We illustrate our approach by obtaining a linear differential equation that represents real data from biomechanics.Comparing our approach with popular methods for estimating the parameters of linear dynamical systems, namely, the non-linear least-squares approach, simulated annealing, parameter cascading

Introduction
Dynamical systems typically translate the natural phenomena into a set of equations based on the motion or equilibrium of the system as determined by its mechanics, chemistry, biology, etc.These models explain the underlying mechanisms that drive or hinder a processes behaviour.A set of linear differential equations denotes a linear dynamical system.Let the p th derivative of the function x at time t be D p x(t).A p th order differential equation specifies how the behaviour of the p th derivative depends on the lower order derivatives, D 0 x(t), . . ., D p−1 x(t), and other external variables, u 1 (t), . . ., u Q (t), that is, where t ∈ [t 1 , t N ], the coefficient functions β r (t|θ) and α q (t|θ) are functions of t that are dependent on a vector of parameters θ and u q (t) is the q th function at time t representing the q th external variable.The differential equation is linear if the functions β r (t|θ), α q (t|θ) and u q (t) do not depend on the values of x. 1 This formulation encompasses a broad range of phenomena, including those observed in climate science, biology and ecology.See for example, [1,2,3] and the references therein.The main challenge is determining the values of the parameters θ, defining β r (t|θ) and α q (t|θ) in (1), that ensure the approximating solution of (1) evaluated at the observed times, adheres to the observed behaviour of the process.We illustrate this problem by presenting an example of a linear differential equation for modelling head acceleration.Figure (1) depicts 133 observations of head acceleration (in cm/msec 2 ) measured 14 milliseconds before and 42.6 milliseconds after a blow to the cranium.The dashed line represents the unit pulse function which denotes the strike to the cranium that lasted one millisecond.The experiment, a simulated motor-cycle crash, is described in detail in [5].Mechanical principles imply that the acceleration x(t) convey the period of the oscillation, the change in its amplitude, as t → ∞ the oscillations decay exponentially to zero, and the size of the impact from the unit pulse respectively.Our objective is to estimate the acceleration x and the parameters θ = [β 0 , β 1 , α] in (2) so that the approximated solution of (2), x(t|θ), evaluated at the observed times, t = [t 1 , . . ., t N ], adheres to the data in Figure (1).
Non-linear least squares (NLS) is the most common approach [6,7,8] for estimating the solution x(t) and the parameters θ for the differential equation in (1).Given a set of initial conditions, D r x(0) for r = 0, . . ., p− 1, and a period of the domain over which the solution is sought, [t 1 , t N ], the solution of (1) can be approximated by a numerical iterative method (e.g.Runge-Kutta methods).
NLS then estimates the parameters θ by minimising the difference between the approximated numerical solution of (1) and the observed data values.Typically this minimisation problem has many local minima.As shown in [9] simulated annealing (SA), introduced by [10], can be used to overcome the topological difficulties in the minimisation problem.The NLS and SA approaches are both computationally intensive as a numerical approximation to the solution of the differential equation in ( 1) is required for each update of θ.Additionally, the initial values D r x(0) for r = 0, . . ., p − 1 are not usually available in exact form.
Therefore, we often need to minimise the objective function with respect to θ and D r x(0) for r = 0, . . ., p−1.This adds a great deal of extra computation and complexity to the optimisation.Parameter cascading (PC) attributable to [11] alleviates the computational cost associated with repeatedly numerically solving the differential equation and does not require the initial values D r x(0) for r = 0, . . ., p − 1 to be available in exact form.PC uses a linear combination of basis functions to approximate the solution of (1) and the estimated parameters θ are obtained by ensuring that the approximating basis function expansion adheres to the data.Similar to NLS, PC has topological difficulties in minimising the data misfit.Smooth functional tempering (SFT) proposed by [12] implements a Bayesian version of PC and borrows insights from parallel tempering [13,14] to overcome the topological difficulties.SFT produces accurate estimates of θ, but it is very computationally expensive.Quick and easy procedures have been proposed by [15,16,17,18] and [19].These methods do not account for the hierarchical structure of the parameters.The parameters that approximate the solution of (1) are dependent on the parameters θ, which determine the shape of the solution.As a consequence, each method reports a considerable increase in the bias of θ when compared to the estimates produced by the SFT, PC, SA or NLS approaches.
Motivated by the drawbacks of the existing methods discussed above, we introduce a version of the PC approach called data to linear dynamics "Data2LD".
Data2LD is a fast and stable version of the PC approach for estimating the parameters of linear dynamical systems.First, we reduce the complexity of the PC estimation procedure, which has the advantages of speed and ease of use.
Then analogous to SA and SFT, we propose an iterative scheme to overcome the topological difficulties in minimising the data misfit.One of the primary benefits of this algorithm is that it facilitates an accurate and stable estimation of the solution, x(t), and the parameters, θ defining β r (t|θ) and α q (t|θ) in (1).In comparison to other techniques, namely NLS, SA, PC and SFT our proposed method benefits from estimates of θ and x(t), with an improved bias and sampling variance obtained at a fraction of the computational cost.
Section (2) briefly reviews the existing approaches for estimating the solution x(t) and the parameters θ from data.Section (3) describes our approach detailing the dynamic model-fitting criteria, proposing an iterative scheme for estimating θ and providing formulae for approximating the sampling variance of θ and x(t).Section (4) illustrates the estimation of the solution and the parameters θ from noisy incomplete data by obtaining a linear differential equation for modelling head acceleration.Section (5) presents a simulated data example and its performance.

Background
For a detailed account of modern methods for estimating parameters in linear and non-linear differential equations see [4].Sections (2.1) to (2.4) briefly outlines four popular approaches for estimating the solution x(t) and the parameters θ of the differential equation in (1).

Non-linear least squares (NLS)
Given an initial estimate θ 0 of θ, and a set of p initial values, D r x(0), r = 0, . . ., p − 1, a numerical approximation to the solution of the differential equation in (1), x(t, θ 0 , D 0 x(0), . . ., D p−1 x(0)), evaluated at the observed times t, is computed using a method for initial value problems such as a Runge-Kutta method.The estimated parameters are then obtained by minimising, with the initial condition θ = θ 0 , using a gradient-based optimisation method (e.g. the trust-region-reflective algorithm).Typically, one obtains the gradient and hessian of (3) evaluated at the current estimate of θ using numerical differentiation (e.g.finite difference approximations).Often, the topology of the objective function in (3) is undesirable with local minima, ridges, ripples and large flat segments.

Simulated annealing (SA)
The Metropolis-Hastings (MH) algorithm draws samples of θ from the conditional distribution π(θ|T, y, t, D 0 x(0), . . ., D p−1 x(0 for a fixed temperature T .As T → ∞, π tends to a uniform probability density function and as T → 0, π tends to a delta function located at the global maximum of (4), which is equivalent to the global minimum of (3).A temperature ladder T i for i = 1, . . ., M , is constructed and the MH algorithm progressively samples from the conditional distribution in (4) as the temperature is adjusted from high to low values.At the u th iteration if the value of π u is greater than , the new θu is accepted.Otherwise, the new θu is accepted at random with a probability 1/(1+exp((π u−1 −π u )/T i ).A smaller temperature or a larger distance between π u−1 and π u will lead to a smaller acceptance probability.SA provides a means to escape local optima by accepting steps which decrease π in hopes of finding a global optimum.

Parameter cascading (PC)
PC approximates the solution of (1) by a linear combination of basis func- where Φ is the N × K matrix containing the basis function φ k (t) evaluated at the locations t and c is a vector of length K containing the corresponding coefficients.PC defines the coefficient c k as a smooth function of the parameters θ and λ, where λ is a regularity parameter that determines the trade-off between x ′ s fit to the data and adherence to the differential equation in (1).For fixed λ, the estimated parameters ĉ(θ, λ), are produced by minimising with respect to c each time the parameter vector θ is updated.The penalty term, L(θ, Φc), in ( 5) is the square L 2 norm of the differential equation in (1) with x replaced by Φc.The parameters, θ, for fixed λ are then estimated so the resulting approximating solution, x, adheres to the data, which is achieved by minimising with respect to θ.The regularity parameter λ is typically chosen by minimising generalized cross validation.Similar to NLS, PC has topological difficulties in the minimisation of (6).

Smooth Functional Tempering (SFT)
SFT implements a Bayesian version of PC for fixed λ, where π(θ) and π(σ 2 ) are prior distributions defined for θ and σ

Data2LD
Here we present a version of PC, called Data2LD, designed for the estimation of the parameters of linear differential equations as in (1).

The dynamic model-fitting criterion
Approximate the solution of the differential equation in (1) by a basis function expansion Assume the coefficients c k in (8) are smooth functions of the parameters θ and ρ, that is, c k (θ, ρ), where ρ controls x's approximation adherence to the differential equation in (1).The basis functions φ k (t) are chosen to reflect the characteristics of the data.For example, if the data exhibit cyclical behaviour then Fourier basis may be desirable.We recommend B-spline basis functions due to their flexibility and computational efficiency.The number of basis functions must be large enough to guarantee that the regularisation is controlled by the choice of the regulating parameter ρ and to ensure a satisfactory approximation of the highest order derivative D p x(t).An exact representation or interpolation of the data is achieved when K = N .As advised in [20], we typically set where O is the order of the B-spline basis functions.We recommend setting the order O > p + 3 to ensure that the highest order derivative D p x(t) is at least approximated with piece-wise cubic functions.For further details on possible basis functions and choices of K see [20].
Let R(θ) be the K × K matrix with entries and S(θ) be a K × 1 vector with entries Then the square L 2 norm of the differential equation in (1) with x(t) replaced by the basis function expansion in (8), can be written as The coefficients ĉ(θ, ρ) for fixed θ and ρ are obtained by minimising the penalised least squares criterion where y is a vector of length N containing the measured observations, Φ is an and c is a vector of length K containing the coefficients of the basis functions.
Equation (10) combines two sources of information about x(t), its fidelity to the data, as measured by the residual sum of squares in the first term in (10), and its adherence to the linear differential equation in (1), as quantified by the L 2 norm of (1) in the second and third term in (10).To facilitate a comparable scale the fist term in (10) is divided by N to obtain the average of the squared residuals and the second and third term in (10) is divided by (t N − t 1 ) to obtain the average of the adherence to the linear differential equation.The regulating parameter ρ can now be defined within the domain [0, 1).If ρ = 0 then the corresponding estimated function, x = Φĉ, is the least squares approximation of the data and hence does not depend on the differential equation.However, as ρ → 1, minimising ( 10) is equivalent to minimising (9).If ( 9) is approximately zero, x is an approximation of the solution of (1).The coefficient values that minimise (10) with respect to c for fixed θ and ρ are given analytically by See the supplementary material for the full derivation of (11).
The estimated parameters of the differential equation θ for fixed ρ are obtained by minimising a dynamic model-fitting criterion where ĉ(θ, ρ) is given in (11).Numerical optimisation methods such as Gauss-Newton can be used to minimise (12) with respect to θ for fixed ρ.The gradient of H(θ|ρ) is required for the Gauss-Newton algorithm and is given analytically by See the supplementary material for the formulae for evaluating dĉ(θ,ρ) dθ .

The iterative scheme to acquire an optimal estimate of θ
The regulating parameter ρ controls the complexity of the surface in (12).
For low values of ρ, x is a least-squares approximation of the data and H(θ|ρ) is convex.Thus, the minimum of H(θ|ρ) with respect to θ, is easy for the Gauss-Newton algorithm to locate.Low values of ρ do not require L( θ) to be small.Consequently x is not a satisfactory approximation of the solution of the differential equation in (1) and θ is an inaccurate estimate of θ.For high values of ρ, H(θ|ρ) is a non-convex surface with flat plains, ripples and a long narrow ridge around the global minimum.In this instance, unless the initial estimate for θ is within the narrow basin of attraction of the global minimum the Gauss-Newton algorithm will converge to a local minima.High values of ρ require L( θ) to be small and therefore x is a better approximation to the solution of the differential equation in (1) and θ is a more precise estimate of θ.
To illustrate this we examine the estimates of θ obtained by minimising  present in the initial optimisation of H(θ|ρ 0 ).The difference between consecutive values of ρ can be larger for ρ < 0.9 for which H(θ|ρ) is relatively convex.
For ρ > 0.9 the difference between consecutive values of ρ must be small as the surface of H(θ|ρ) can change substantially from one value of ρ to the next.
As a consequence, H(θ|ρ) is minimised with logistic values of ρ chosen so that the minimum of each H(θ|ρ) is "close" to the previous one.This will help to preclude difficulties in finding the global minimum of H(θ|ρ) from one iteration to the next.Analogous to the many choices for the temperature ladder in simulated annealing, see [21] for details, one could envisage many possible methods for reducing ρ from one iteration to the next.Our approach proposed herein is a rather conservative approach and we acknowledge that an optimal reduction of ρ is an area for future research.The estimation procedure stops when the estimated parameters converge.The details of the Data to linear dynamics iterative scheme are below: Step 1 Specify the initial values: θ0 (e.g.θ(0) = [0.01, . . ., 0.01]) and γ 0 = −4.
Let u max be the value of u when the iterative scheme stopped.The iterative scheme produces θ = θumax the estimated parameters of the differential equation that best approximate the data.Substituting θ and ρ = ρumax into (11) produces an estimate of the coefficients of the basis function expansion ĉ( θ, ρ).
The approximated solution of the differential equation is x = Φĉ( θ, ρ) and its degrees of freedom are where the . Here degrees of freedom refers to the effective dimensionality of x, as ρ → 1 it tends to the dimensionality of the solution space for the differential equation.

Approximating the sampling variation for θ and x
Assuming that y is normally distributed with variance σ 2 y .The conditional sampling variance of the estimated parameters of the differential equation can be approximated using the delta method [22]: .The formula for evaluating d θ dy is given in the supplementary material.
The point-wise conditional sampling variance of the estimated solution of the differential equation evaluated at the data points x is also approximated using the delta method: The formula for evaluating dĉ( θ, ρ) dy is given in the supplementary material.

A differential equation for modelling head acceleration
We used three order one B-splines over the knots [0, 14, 15, 56] with coefficient vector [0,1,0] to represent the unit pulse function u(t) in ( 2).For the basis expansion of x(t) we used order five B-spline functions, which by their nature have discontinuous third derivatives if all knots are singletons.The unit pulse function u(t) is discontinuous, which implies a discontinuity in D 2 x(t).
To achieve curvature discontinuity at the impact point and at that point plus one, we placed three knots at these locations.We put no knots between the first observation and the impact point, where the data indicate a flat trajectory and eleven equally spaced knots between the impact point plus one and the final observation time.We estimated the coefficients c and the parameters θ = [β 0 , β 1 , α] in (2) using Data2LD.plying that the acceleration is an under-damped process; after the blow to the cranium, the acceleration will oscillate with a decreasing amplitude that will  quickly decay to zero.The parameters of the differential equation suggest that it will take approximately 66 milliseconds after impact for the average acceleration to return to zero.The estimated function x has an effective degrees of freedom that is equal to 2.45, and a root mean squared error that is 0.05.
Figure (4) shows the accelerometer readings of the brain tissue, the fitted curve produced by Data2LD (solid line), the approximated 95% point-wise confidence interval for the fitted curve (dashed line) and the approximated 95% point-wise prediction interval for the fitted curve (grey band).
The differential equation captures the trend in the acceleration of the brain tissue.It conveys that the acceleration peaks at approximately 6.2 milliseconds after impact and troughs at around 19.8 milliseconds after impact.

Simulation Study: head acceleration model
Consider the differential equation  (solid line), the approximated 95% point-wise confidence interval for the curve (dashed line) and the approximated 95% point-wise prediction interval for the curve (grey region).
where u is one for 14 ≤ t ≤ 15 and zero otherwise and the initial conditions are x(0) = 0 and Dx(0) = 0. Let x be the analytic solution of ( 14), evaluated at N equally-spaced observations over the domain [0, 60].The data y are generated by y = x + ε where ε is a vector of N independent normally distributed random values with mean 0 and standard deviation σ × range(x).One thousand simulated samples are generated for σ = 0.01, 0.05, 0.10, and for sample sizes N = 21, 51, 101.The basis functions for approximating x and u are set up as described in Section 4.
To implement SFT we let the priors for β 0 , β 1 and α be normal distributions with mean 0 and variance 1.The prior for σ 2 was chosen to be 1 σ 2 .Four parallel chains were used with four different temperatures {10, 100, 1000, 10000}.We ran fifty thousand parallel MCMC chains and each chain was initialised with the same values.As suggested in [10] for SA we set the initial temperature to T = 100 and reduced it using a Boltzmann schedule.
Table (1) provides the root mean squared error of the estimated parameters θ with respect to the true parameters θ true , RMSE( θ) for the estimates obtained by Data2LD, simulated annealing (SA), smooth functional tempering (SFT), non-linear least squares (NLS) and parameter cascading (PC).The minimum RMSE( θ) for each of the nine simulated data configurations involving three sample sizes and three levels of error is highlighted in grey.Data2LD had the best performance for all three parameters, and all three are well determined by the data with the RMSE's being less than .5% of the parameter magnitudes.
SA and SFT are expected to produce lower RMSEs relative to NLS and PC, respectively, as these approaches are designed to deal with the complex topology of the parameter space.While SFT yielded lower RMSEs relative to PC for β 0 and β 1 , it increased the RMSE for α across all configurations.SFT and PC showed considerably higher RMSE(α) across all sample sizes relative to the other approaches.This indicates that neither SFT nor PC adequately estimated the sharp change in the process due to the impact of u(t).SA showed a substantial increase in RMSE relative to NLS for N = 21.This suggests that SA does not perform well when the sample size of the data set is small.
Table (2) provides the 95% coverage probability of the estimated confidence intervals CP( θ) for the estimates obtained by Data2LD, SA, SFT, NLS and PC.
The coverage probability measures the proportion of the estimated confidence intervals for θ that contained the true parameters θ true over 1000 simulations.
The closest to 95% for each of the nine simulated data configurations involving three sample sizes and three levels of error is highlighted in grey.In most  SA showed an increase in RMSE relative to NLS for N = 21 and N = 51.Indicating that SA provides an improvement in the estimate of the solution of the differential equation only when N is large.SFT showed an increase in RMSE relative to PC for σ = 0.1.Indicating that SFT does not provide an improve-

Discussion and Conclusions
Dynamical systems can provide a conceptual understanding of how processes evolve, which can help guide their management and prediction or can simply provide a tractable, flexible and parsimonious model of the processes.The parameters of a dynamical system determine the interrelationships between the processes which describe how these objects be it physical, engineering or demographic behave.These parameters are often unknown and must be estimated from the observed data.The most popular approaches for parameter estimation for dynamical systems are smooth functional tempering (SFT), parameter cascading (PC), simulated annealing (SA) and non-linear least squares (NLS).
The NLS and PC approach involves obtaining the minimum of ( 3) and (6) with respect to the parameters' of the differential equation.These parameter spaces can exhibit complex topology including multi-modality, ripples and nar-row ridges and as such, can be difficult to navigate.SA and SFT are popular approaches for finding the global minimum of ( 3) and (6).As shown in [9,12] and herein, SA and SFT often provide improved estimates of the parameters and the solution of the differential equation relative to NLS and PC.However, both SA and SFT are very computationally expensive and do not provide an adequate estimate of the uncertainty associated with the estimated parameters of the differential equation.
We propose Data2LD a version of the PC approach that has been tailored for linear systems.First, we reduce the complexity of the PC estimation procedure, which has the advantages of speed and ease of use.Then analogous to SA, we propose an iterative scheme to overcome the topological difficulties in minimising the data misfit.One of the primary benefits of this algorithm is that it facilitates accurate and stable estimation of the solution, x(t), and the parameters, θ defining β r (t|θ) and α q (t|θ) in (1).
We compared Data2LD with the popular existing approaches, namely SFT, PC, SA and NLS.In terms of statistical measures of performance such as rootmean-squared error for parameters estimates and the estimates of the solution of the differential equation, our simulations suggest an advantage for Data2LD.
For large sample sizes, Data2LD and SA have a similar estimation accuracy with Data2LD having a computational advantage of about five orders of magnitude.
SA does not perform well when the sample size of the data set is small.PC and SFT has difficulty estimating the sharp change in the solution at the impact point resulting in poor estimates of the parameters' of the differential equation.
We are extending Data2LD to linear dynamic systems along with data observed over space and time where processes can be denoted by a set of linear partial differential equations such as reaction-diffusion-transport family.
A Matlab package with source code and datasets for the examples presented in this article is available at https://github.com/mcareyucd/Data2LD-Matlab.
An R-package "Data2LD" that contains functions for using differential equations as modelling objects can be obtained from CRAN at (https://cran.r-project.org/web/packages/Data

Figure 1 :
Figure 1: The circles illustrate the accelerometer readings of the head acceleration before and after a blow to the cranium of a cadaver.The dashed line represents a unit pulse function which denotes the blow to the skull.This function initiates at 14 milliseconds and lasts for 1 millisecond.

2
respectively.Similar to PC, the sampling problem becomes difficult due to the multi-modality of the posterior surface π(θ|y).Gaps between modes can be traversed at lower values of λ, while individual modes can be efficiently explored at higher values of λ.In contrast to SA, SFT replaces the unidirectional reduction of the temperature ladder by a set of concurrent simulations at M different temperatures {λ m |m = 1, . . ., M }.At the u th iteration, each of the M chains independently performs a MH step to update θ u 1 , . . ., θ u M .Let U be a uniform random variable on [0, 1].If U is less than a threshold value (e.g.0.5) then a randomly selected neighbouring pair of chains, refereed to as m and m + 1, exchange states, that is, θ u m → θ u m+1 and θ u m+1 → θ u m .This exchange is accepted with probability min 1, πm(θ u m+1 |y)πm+1(θ u m |y) πm(θ u m |y)πm+1(θ u m+1 |y) .SFT increases the efficiency of the sampling of π(θ|y) and thus can improve the convergence to a global minimum.

H
(θ|ρ) with respect θ for various values of ρ for a simulated data set.The data are obtained by evaluating the analytic solution of (2) with θ = [−0.05,−0.15, 0.39] at 101 equally spaced points within the domain [0, 60] and adding a vector of 101 independently normally distributed random variables with mean 0 and σ = 0.05. Figure (2) shows the surfaces and contours of H(θ|ρ) for θ = [β 0 , β 1 , 0.39], with β 0 ranging from −0.55 to 0.45 and β 1 ranging from −0.65 to 0.35 with ρ = 0.99, 0.95, 0.71 and 0.50.The circle is the true value for (β 0 , β 1 ) = (−0.05,−0.15) and the asterisk represents the minimum of H(θ|ρ) for the respective ρ.For ρ = 0.50 the surface is quadratic, but the minima θ = [−0.05,−0.10, 0.39], is not an accurate estimate of θ.For ρ = 0.99 the surface has flat plains, ripples and a long narrow ridge which is difficult for gradient decent methods to navigate, but the minima θ = [−0.05,−0.15, 0.39], is an accurate estimate of θ. Figure (2) also shows the changes in the topology of H(θ|ρ) with respect ρ.As ρ reduces the non-convex surface, which is a delta function located at the global minimum, illustrated by H(θ|0.99), is smoothed to a convex surface with a flat area around the global minimum as illustrated by H(θ|0.5).Thus, the regulating parameter ρ controls the prominence of local and global minima and as such has the same role as the temperature parameter T in simulated annealing described in Section (2.2).Data to linear dynamics, navigates the complicated topology of H(θ|ρ) by starting with a relatively small value of ρ (e.g.ρ 0 = 0.04) and an infeasible (exterior) point θ(0) (e.g.θ(0) = [.01, . . ., .01])so that no steep valleys are

Figure 2 :
Figure 2: The optimisation surface H(θ|ρ) for the simulation of the head acceleration analysis discussed in Section (??) with ρ = 0.5, 0.71, 0.95 and ρ = 0.99 with β 0 ranging from −0.55, 0.45 and β 1 ranging from −0.65, 0.35.In each figure the global minimum is identified by the large dot and the local minima is identified by the asterisk.

Figure ( 3 )
shows how the three parameters of the differential equation in (2) and their approximated confidence intervals vary as ρ increased to 0.99.As shown in Figure (3) when the influence of the differential equation increases to the point where it is the primary determinant of the parameters, the parameter values stabilise, and the approximated confidence intervals reduce.The final parameter estimates with 95% confidence intervals are, β0 = −0.057± 0.005 for the stiffness, β1 = −0.15± 0.03 for the damping and α = 0.40 ± 0.06 for the force from the unit pulse function.Im-

Figure 3 :
Figure 3: The values of the three parameters and their approximated 95% confidence intervals for the head impact data over values of ρ converging to 0.99.

Figure 4 :
Figure 4: The accelerometer readings of the brain tissue before and after a blow to the cranium are indicated by the circles.The fitted curve produced by Data2LD with ρ = 0.99 cases, Data2LD obtained the most accurate approximation of the uncertainty associated with the estimate θ.PC substantially underestimated the coverage probability in all cases except for σ = 0.01 and N = 101.SFT overestimated the coverage probability in all cases.Table (3) assess the accuracy of the solution of the differential equation by reporting RMSE(x), the root mean squared error of the estimated solution of the differential equation x with respect to the true solution.The minimum RMSE for each of the nine configurations is highlighted in grey.Data2LD has

Table 2 :
The 95% coverage probability of the estimated confidence intervals averaged over

Table 3 :
The root mean squared errors (RMSE) times 100 averaged over 1000 simulations for the estimated solution of the differential equation for the head acceleration analysis estimated