Data 2 Linear Dynamics

What is Data 2 Linear Dynamics (Data2LD)?

Data 2 Linear Dynamics (Data2LD) estimates the solution and the parameters of linear dynamical systems from incomplete and noisy observations of the underlying processes.

It uses a linear combination of spline basis functions to approximate the implicitly defined solution of the dynamical system. While also estimating the systems' parameters by requiring that this approximating solution adheres to the data.

Papers:

Carey, M., Ramsay, J.O. (2020) 'Fast Stable Parameter Estimation for Linear Dynamical Systems'. Computational Statistics & Data Analysis

Carey, M., Gath, E., Hayes, K. (2016) 'A generalized smoother for linear ordinary differential equations'. Journal of Computational and Graphical Statistics.

Carey, M., Gath, E., Hayes, K. (2014) 'Frontiers in financial dynamics'. Research in International Business and Finance

Let's consider a study on traumatic brain injury (TBI), which contributes to just under a third (30.5%) of all injury-related deaths in the US and is caused by a blow to the head. Figure (1) shows the 133 accelerometer readings taken over 55.2 milliseconds. The dashed line represents the impulse function which denotes the blow to the head.

Figure (1)

The laws of motion tell us that the acceleration x(t) can be modelled by a second-order linear ordinary differential equation (ODE) with input a unit pulse u(t) representing the blow to the head and shown in the dashed lines in Figure (1).

This ODE
\begin{equation}
\frac{\textrm{d}^2x(t)}{\textrm{d}t^2} + \beta_{0} x(t) + \beta_{1} \frac{\textrm{d}x(t)}{\textrm{d}t} + \alpha u(t)
\end{equation}

contains three parameters $\beta_{0},\beta_{1}$ and $\alpha,$ and these convey the rate of the restoring force (as $t \rightarrow \infty,$ the acceleration will tend to revert back zero), the rate of the friction force (as $t \rightarrow \infty,$ the oscillations in the acceleration reduce to zero) and the rate of force from the unit pulse.

While there are several methods for estimating ODE parameters with partially observed data, they are invariably subject to several problems including high computational cost, sensitivity to initial values or large sampling variability.

Ā Data to linear dynamics (Data2LD) overcomes these issues and produces estimates of the ODE parameters that have less bias, a smaller sampling variance and a ten-fold improvement in computation in comparsion to existing approaches.

Demonstration of Data2LD: Head Impact Analysis

This page provides a detailed description of the MATLAB calculations necessary to run Data2LD code for estimating the parameters and the solution of linear differential equations.

Contents

Set-up

clear
clear functions

Load the data

load motorcycledata.txt

motot = motorcycledata(:,2);  'ko', [0,60], [0,0], 'k:');
set(phdl, 'LineWidth', 2)
axis([min(yCell{1}(:,1))-0.01,max(yCell{1}(:,1))+0.01,-1.0,1.5])
xlabel('fontsize{16} Time (milliseconds)')
ylabel('fontsize{16} Acceleration (cm/msec^2)')

Figure (1) illustrates $133$ observations of head acceleration measured 14 milliseconds before and 42.6 milliseconds after a blow to the cranium of a cadaver. The dashed line represents the impulse function which denotes the blow to the cranium that lasted 1 millisecond. The experiment, a simulated motor-cycle crash, is described in detail by Schmidt et al (1981).

The second order linear differential equation

The second order linear differential equation, as used in this example is

\begin{equation}
\frac{\textrm{d}^2x(t)}{\textrm{d}t^2} + \beta_{0} x(t) + \beta_{1} \frac{\textrm{d}x(t)}{\textrm{d}t} + \alpha u(t)
\end{equation}

The three parameters $beta_{0}$, $beta_{1} $ and $alpha$ in the linear ODE convey the period of the oscillation, the change in its amplitude, as $t rightarrow infty$ the oscillations decay exponentially to zero, and the size of the impact from the unit impulse respectively.

Our aim is to estimate the acceleration $x(t)$ and the parameters $theta=[beta_{0},beta_{1},alpha]$ from the data in Figure (1).

Set up basis functions for X

motorng   = [yCell{1}(1,1),yCell{1}(end,1)];
knots     = [motorng(1),14,14,14,15,15,linspace(15,motorng(end),11)];
norder    = 5;
nbasis    = length(knots)+ (norder - 2);
Xbasisobj = create_bspline_basis(motorng,nbasis,norder,knots);


For the expansion of $x(t)$ we used order five B-spline functions, which by their nature have discontinuous third derivatives if all knots are singletons. 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 time of observation.

Create a functional data object for the external/forcing function


Three order one B-splines over the knots [0, 14, 15, 56] with coefficient vector [0,1,0] represented the unit pulse function $u(t)$

Run Data2LD

...
    estimate,EFCell);

Rho = 0.99995

Iter.        Criterion   Grad Length
  0           0.054430     0.002710
  1           0.054430    12.076390
Elapsed time is 2.639805 seconds.

Optimal values:
Rho       = 1      
Theta1    = -0.057    
Theta2    = -0.155    
Theta3    = 0.402    
sdTheta1  = 0.003 
sdTheta2  = 0.019 
sdTheta3  = 0.032 
MSE       = 0.054
ISE       = 0
Df        = 2.234
GCV       = 0.056    

Plot the fitted curve, the confidence and prediction intervals

figure()
hp = patch([yCell{1}(:,1); yCell{1}(end:-1:1,1);yCell{1}(1)],...
    [Results_cell{4}(:,4); Results_cell{4}(end:-1:1,5); Results_cell{4}(1)],[0.85 0.85 0.85]);
hold on;
plot(yCell{1}(:,1),[Results_cell{4}(:,2),Results_cell{4}(:,3)],'k--')
plot(yCell{1}(:,1),Results_cell{4}(:,1),'k-')
plot(yCell{1}(:,1), yCell{1}(:,2), 'ko');
xlim([min(yCell{1}(:,1))-0.1,max(yCell{1}(:,1))+0.1])
ylim([min(Results_cell{4}(:,4))-0.1,max(Results_cell{4}(:,5))+0.1])
xlabel('fontsize{16} Time (milliseconds)')
ylabel('fontsize{16} Acceleration (cm/msec^2)')
legend('Prediction Interval','Confidence Interval','Confidence Interval','Fitted Curve','Data')
set(findall(gca, 'Type','text'), 'FontSize', 20)

This Figure illustrates 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 $hat{rho} = 0.99$ (solid line), the approximated 95% pointwise confidence interval for the curve (dashed line) and the approximated 95% pointwise prediction interval for the curve (grey region)

Plot the estimated parameters of the ODE and their approximated CI with respect to rho

Res_Mat  = table2array(Res);
for i=1:3
subplot(3,1,i)
indexOfInterest = (Res_Mat(:,1) > 0.95);
plot(Res_Mat(indexOfInterest,1),Res_Mat(indexOfInterest,i+1),'-*k')
hold on;
errorbar(Res_Mat(indexOfInterest,1),Res_Mat(indexOfInterest,i+1),...
    t_v.*Res_Mat(indexOfInterest,i+4),'ko')
xlabel('$rho

The values of the three parameters and their approximated 95% confidence intervals for the head impact data over values of $rho$ converging to 1.0. Note the reduction in the width of the approximated 95% confidence intervals as $rho$ increases and the overlap of the final estimates of the ODE parameters, which illustrates that the parameter values have stabilized for high values of $rho$.

Check if the residuals are approx. normal

figure()
h = normplot(yCell{1}(:,2)-Results_cell{4}(:,1));
h(1).LineWidth = 2;
h(2).LineWidth = 2;
h(3).LineWidth = 2;


Check the fitted vrs residuals plot

figure()
plot(Results_cell{4}(:,1),yCell{1}(:,2)-Results_cell{4}(:,1),'ko')