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);  %  time in milliseconds
motoy = motorcycledata(:,3);  %  deformation in 0.1 millimeters

% Adjust the data for baseline and plot

impact  = 14;  %  impact time
baseind = motot < impact;
basey   = mean(motoy(baseind));

% Remove baseline, change time, and convert to centimeters
motoy   = (basey - motoy)/100.0;
n       = length(motoy);

% Create a cell containing the data
yCell = cell(1);
yCell{1}(:,1)  = motot;
yCell{1}(:,2)  = motoy;

% Plot the data
figure(1)
phdl=plot(yCell{1}(:,1), yCell{1}(:,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);

%Create a cell for the basis for each variable
XbasisCell = cell(1);
XbasisCell{1} = Xbasisobj;

%Plot the basis
figure()
plot(Xbasisobj)

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

%The known external function is a unit impulse
Ubasis = create_bspline_basis(motorng, 3, 1, [motorng(1),14,15,motorng(2)]);
Ufd = fd([0;1;0],Ubasis);

%Plot the external function
figure()
plot(Ufd)

%Create a cell for the external functions
EFCell = cell(1);
EFCell{1} = Ufd;

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

%Intial parameter estimates
intial_theta = {0.1,0.1,0.1};

%Which parameters would we like to estimate. In this case all parameters.
estimate = [1,1,1];

%Run Data2LD
[Res,Results_cell] = D2LD_ODE(yCell,XbasisCell,intial_theta,...
    estimate,EFCell);


Elapsed time is 0.159864 seconds.
Elapsed time is 0.158650 seconds.
Elapsed time is 0.013674 seconds.
Rho = 0.5

Iter.        Criterion   Grad Length
  0           0.045955     0.001943
  1           0.045726     0.000386
  2           0.045676     0.000061
  3           0.045673     0.000008
  4           0.045673     0.000001
  5           0.045673     0.000000
Convergence reached
Rho = 0.73106

Iter.        Criterion   Grad Length
  0           0.045718     0.000116
  1           0.045716     0.000054
  2           0.045715     0.000002
  3           0.045715     0.000000
Convergence reached
Rho = 0.8808

Iter.        Criterion   Grad Length
  0           0.045854     0.000433
  1           0.045851     0.000265
  2           0.045841     0.000085
  3           0.045841     0.000017
  4           0.045841     0.000006
  5           0.045841     0.000001
Convergence reached
Rho = 0.95257

Iter.        Criterion   Grad Length
  0           0.046300     0.002776
  1           0.046286     0.000403
  2           0.046281     0.000120
  3           0.046281     0.000033
  4           0.046281     0.000011
  5           0.046281     0.000003
Convergence reached
Rho = 0.97069

Iter.        Criterion   Grad Length
  0           0.046799     0.003304
  1           0.046789     0.000611
  2           0.046780     0.000160
  3           0.046780     0.000041
  4           0.046779     0.000010
  5           0.046779     0.000002
Convergence reached
Rho = 0.98201

Iter.        Criterion   Grad Length
  0           0.047602     0.006326
  1           0.047583     0.001372
  2           0.047556     0.000473
  3           0.047553     0.000150
  4           0.047553     0.000041
  5           0.047553     0.000012
  6           0.047553     0.000003
Convergence reached
Rho = 0.98901

Iter.        Criterion   Grad Length
  0           0.048725     0.010017
  1           0.048698     0.002738
  2           0.048626     0.000630
  3           0.048623     0.000004
  4           0.048623     0.000000
Convergence reached
Rho = 0.99331

Iter.        Criterion   Grad Length
  0           0.050119     0.013939
  1           0.050085     0.005582
  2           0.049909     0.000710
  3           0.049906     0.000099
  4           0.049906     0.000002
  5           0.049906     0.000000
Convergence reached
Rho = 0.99593

Iter.        Criterion   Grad Length
  0           0.051423     0.017094
  1           0.051388     0.010112
  2           0.051152     0.002012
  3           0.051150     0.000102
  4           0.051150     0.000010
  5           0.051150     0.000004
Convergence reached
Rho = 0.99753

Iter.        Criterion   Grad Length
  0           0.052353     0.016602
  1           0.052329     0.013570
  2           0.052175     0.002229
  3           0.052175     0.000187
  4           0.052175     0.000086
  5           0.052175     0.000007
Convergence reached
Rho = 0.9985

Iter.        Criterion   Grad Length
  0           0.053037     0.013911
  1           0.053024     0.015415
  2           0.052946     0.001457
  3           0.052946     0.000280
  4           0.052946     0.000045
Convergence reached
Rho = 0.99909

Iter.        Criterion   Grad Length
  0           0.053532     0.010652
  1           0.053525     0.014624
  2           0.053490     0.000924
  3           0.053490     0.000239
  4           0.053490     0.000062
Convergence reached
Rho = 0.99945

Iter.        Criterion   Grad Length
  0           0.053873     0.007592
  1           0.053870     0.010930
  2           0.053855     0.000854
  3           0.053855     0.000173
  4           0.053855     0.000068
Convergence reached
Rho = 0.99966

Iter.        Criterion   Grad Length
  0           0.054099     0.005132
  1           0.054097     0.006063
  2           0.054091     0.001206
  3           0.054091     0.000146
  4           0.054091     0.000067
Convergence reached
Rho = 0.9998

Iter.        Criterion   Grad Length
  0           0.054243     0.003350
  1           0.054243     0.003442
  2           0.054241     0.001068
  3           0.054241     0.000109
  4           0.054241     0.000053
Convergence reached
Rho = 0.99988

Iter.        Criterion   Grad Length
  0           0.054335     0.002153
  1           0.054335     0.001837
  2           0.054334     0.000916
  3           0.054334     0.000086
  4           0.054334     0.000041
Convergence reached
Rho = 0.99993

Iter.        Criterion   Grad Length
  0           0.054392     0.001454
  1           0.054392     8.390675
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$ increase and the overlap of the final estimates of the ODE parameters, which illustrates that the parameter values have stabilised 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;

% Our qqplot looks pretty good and indicates a normal distribution, as it’s generally a straight line.

Check the fitted vrs residuals plot

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

% Should see random scatter in a constant band around zero that is no patterns, and no funnels.