handwrit.m

手写签名数据

bfd     = fd(zeros(1,1), constbasis);
bfdPar  = fdPar(bfd, 1, 0);
bwtcell{1,1} = bfdPar;
bwtcell{1,2} = bfdPar;
bwtcell{1,3} = bfdPar;

%  carry out principal differential analysis

[bestwtcell, aestwtcell, resfdcell] = ...
    pda_fd(xfdcell, bwtcell, awtcell, ufdcell, difeorder);

%  evaluate forcing functions

resfd  = resfdcell{1};
resmat = eval_fd(fdatime, resfd);

MSY = mean(mean(resmat.^2));

%  Set the number of basis functions to 125,

bfd     = fd(zeros(nwbasis,1), wbasis);
bfdPar  = fdPar(bfd, 1, 0);
bwtcell{1,1} = bfdPar;
bwtcell{1,2} = bfdPar;
bwtcell{1,3} = bfdPar;

%  carry out principal differential analysis

[bestwtcell, aestwtcell, resfdcell] = ...
    pda_fd(xfdcell, bwtcell, awtcell, ufdcell, difeorder);

%  evaluate forcing functions

resfd  = resfdcell{1};
resmat = eval_fd(fdatime, resfd);

%  compute a squared multiple correlation measure of fit
%  MSY = mean(mean(resmat.^2)); %  Used only with constant basis
%  Uncomment this line when the constant basis is used, and
%  comment it out otherwise.

MSE = mean(mean(resmat.^2));
RSQ = (MSY-MSE)/MSY;

disp(['R-squared = ',num2str(RSQ)])

%  Plot the weight functions

for j=1:3
    subplot(3,1,j)
    plot(getfd(bestwtcell{1,j}));
    ylabel(['Weight function ',num2str(j-1)]);
end

%  Plot the first derivative weight function.

b2fdX = getfd(bestwtcell{1,2});
b2vecX = eval_fd(fdatime, b2fdX);
b2meanX = mean(b2vecX);

subplot(1,1,1)
lhdl = plot(fdatime, b2vecX, '-', ...
            [0, 2300], [b2meanX, b2meanX], 'g--');
set(lhdl, 'LineWidth', 2);
hold on
plotrange = [0,6e-3];
for k=1:length(cycle)
    lhdl = plot([cycle(k),cycle(k)],plotrange,'r:');
    set(lhdl, 'LineWidth', 1')
end
for j=1:length(featureindex)
    k = featureindex(j);
    lhdl = plot([cycle(k),cycle(k)],plotrange,'m-');
    set(lhdl, 'LineWidth', 2')
end
hold off
xlabel('\fontsize{16} Time (milliseconds)')
ylabel('\fontsize{16} \beta_1(t)')
axis([0, 2300, 0, 6e-3])

%  plot the weight coefficient for the second derivative

b3fdX = getfd(bestwtcell{1,3});
b3vecX = eval_fd(fdatime, b3fdX);
b3meanX = mean(b3vecX);

subplot(1,1,1)
lhdl = plot(fdatime, b3vecX, '-', ...
            [0, 2300], [b3meanX, b3meanX], 'g--');
set(lhdl, 'LineWidth', 2);
hold on
plotrange = [-.05,.05];
for k=1:length(cycle)
    lhdl = plot([cycle(k),cycle(k)],plotrange,'r:');
    set(lhdl, 'LineWidth', 1')
end
for j=1:length(featureindex)
    k = featureindex(j);
    lhdl = plot([cycle(k),cycle(k)],plotrange,'m-');
    set(lhdl, 'LineWidth', 2')
end
hold off
xlabel('\fontsize{16} Time (milliseconds)')
ylabel('\fontsize{16} \beta_2(t)')
axis([0, 2300, plotrange])

%  display coefficients for forcing weight functions

disp(getcoef(getfd(aestwtcell{1})))
disp(getcoef(getfd(aestwtcell{2})))
 
%  plot all forcing functions

subplot(1,1,1)
plot(fdatime, 1e9*resmat, '-', fdatime, 1e9*D3fdameanmat(:,1), 'r:')
xlabel('\fontsize{16} Milliseconds')
ylabel('\fontsize{16} Meters/sec/sec/sec')
axis([0,2300,-200,200])

%  plot the mean forcing function along with third deriv.

resmeanfd  = mean(resfd);
resmeanvec = eval_fd(fdatime, resmeanfd);

subplot(1,1,1)
plot(fdatime, 1e9*resmeanvec, 'b-', ...
     fdatime, 1e9*D3fdameanmat(:,1), 'r:')
xlabel('\fontsize{16} Milliseconds')
ylabel('\fontsize{16} Meters/sec/sec/sec')
axis([0,2300,-200,200])

% Define a functional data object for the 
%  three derivative weight functions

wcoef1 = getcoef(getfd(bestwtcell{1}));
wcoef2 = getcoef(getfd(bestwtcell{2}));
wcoef3 = getcoef(getfd(bestwtcell{3}));

wcoef  = [wcoef1, wcoef2, wcoef3];
wfd    = fd(wcoef,wbasis);

%  solve the equation

xstart = eye(3);
xstart(1,1) =   fdameanmat(1,1);
xstart(2,2) = D1fdameanmat(1,1);
xstart(3,3) = D2fdameanmat(1,1);

odeoptions = odeset('RelTol', 1e-6);

[tp1, xp1] = ode45(@derivs, fdatime, xstart(:,1), odeoptions, wfd);
[tp2, xp2] = ode45(@derivs, fdatime, xstart(:,2), odeoptions, wfd);
[tp3, xp3] = ode45(@derivs, fdatime, xstart(:,3), odeoptions, wfd);

%  plot the three solutions

umatx   = [xp1(:,1),xp2(:,1),xp3(:,1)];
Dumatx  = [xp1(:,2),xp2(:,2),xp3(:,2)];
D2umatx = [xp1(:,3),xp2(:,3),xp3(:,3)];

subplot(3,1,1)
plot(fdatime, umatx, '-', [0,2300], [0,0], 'r:');
title('Function')

subplot(3,1,2)
plot(fdatime, Dumatx, [0,2300], [0,0], 'r:')
title('First Derivative');

subplot(3,1,3)
plot(fdatime, D2umatx, [0,2300], [0,0], 'r:')
title('Second Derivative');

%  plot fit to each curve ... hit any key after each plot

subplot(1,1,1)
index  = 1:20;
fdamat = eval_fd(fdatime, fdafd);
zmat   = [ones(1401,1),fdatime-1150,umatx];
for i = index
   xhat = zmat * (zmat\fdamat(:,i,1));
   lhdl = plot(fdatime, xhat, '-', ...
               fdatime, fdamat(:,i,1), '--');
   set(lhdl, 'LineWidth', 2)
   title(['X curve ',num2str(i)])
   axis([0, 2300, -0.04, 0.04])
   legend('\fontsize{13} Fit', '\fontsize{13} Observed', ...
          'Location', 'NorthWest')
   pause;
end

%  ------------------------------------------------------------
%  analysis for Y: a third order equation forced by a 
%  constant and time
%  ------------------------------------------------------------

%  Define the variable

yfdcell{1} = fdafd(:,2);

%  Define the variable with a reduced basis to save
%  computation time

ybasis = xbasis;
yfdcell{1} = data2fd(squeeze(fdamat(:,:,2)), fdatime, ybasis);

%  A constant basis is used for a background level of error

bfd     = fd(zeros(1,1), constbasis);
bfdPar  = fdPar(bfd, 1, 0);
bwtcell{1,1} = bfdPar;
bwtcell{1,2} = bfdPar;
bwtcell{1,3} = bfdPar;

%  carry out principal differential analysis

[bestwtcell, aestwtcell, resfdcell] = ...
    pda_fd(yfdcell, bwtcell, awtcell, ufdcell, difeorder);

%  evaluate forcing functions

resfd  = resfdcell{1};
resmat = eval_fd(fdatime, resfd);

MSY = mean(mean(resmat.^2));

%  Use 125 basis functions.

bfd     = fd(zeros(nwbasis,1), wbasis);
bfdPar  = fdPar(bfd, 1, 0);
bwtcell{1,1} = bfdPar;
bwtcell{1,2} = bfdPar;
bwtcell{1,3} = bfdPar;

%  carry out principal differential analysis

[bestwtcell, aestwtcell, resfdcell] = ...
    pda_fd(yfdcell, bwtcell, awtcell, ufdcell, difeorder);

%  evaluate forcing functions

resfd  = resfdcell{1};
resmat = eval_fd(fdatime, resfd);

%  compute a squared multiple correlation measure of fit
%  MSY = mean(mean(resmat.^2)); %  Used only with constant basis
%  Uncomment this line when the constant basis is used, and
%  comment it out otherwise.

MSE = mean(mean(resmat.^2));
RSQ = (MSY-MSE)/MSY;

disp(['R-squared = ',num2str(RSQ)])

%  Plot the weight functions

for j=1:3
    subplot(3,1,j)
    plot(getfd(bestwtcell{1,j}));
    ylabel(['Weight function ',num2str(j-1)]);
end

%  Plot the second derivative weight, defining the period
%  of a harmonic oscillator.

b2fdY = getfd(bestwtcell{1,2});
b2vecY = eval_fd(fdatime, b2fdY);
b2meanY = mean(b2vecY);

subplot(1,1,1)
lhdl = plot(fdatime, b2vecY, '-', ...
            [0, 2300], [b2meanY, b2meanY], 'g--');
set(lhdl, 'LineWidth', 2);
hold on
plotrange = [0,6e-3];
for k=1:length(cycle)
    lhdl = plot([cycle(k),cycle(k)],plotrange,'r:');
    set(lhdl, 'LineWidth', 1')
end
for j=1:length(featureindex)
    k = featureindex(j);
    lhdl = plot([cycle(k),cycle(k)],plotrange,'m-');
    set(lhdl, 'LineWidth', 2')
end
hold off
axis([0, 2300, 0, 6e-3])

%  Plot the third derivative weight

b3fdY = getfd(bestwtcell{1,3});
b3vecY = eval_fd(fdatime, b3fdY);
b3meanY = mean(b3vecY);

subplot(1,1,1)
lhdl = plot(fdatime, b3vecY, '-', ...
            [0, 2300], [b3meanY, b3meanY], 'g--');
set(lhdl, 'LineWidth', 2);
hold on
plotrange = [-0.07,0.07];
for k=1:length(cycle)
    lhdl = plot([cycle(k),cycle(k)],plotrange,'r:');
    set(lhdl, 'LineWidth', 1')
end
for j=1:length(featureindex)
    k = featureindex(j);
    lhdl = plot([cycle(k),cycle(k)],plotrange,'m-');
    set(lhdl, 'LineWidth', 2')
end
hold off
axis([0, 2300, -0.07, 0.07])

%  display coefficients for forcing weight functions

disp(getcoef(getfd(aestwtcell{1})))
disp(getcoef(getfd(aestwtcell{2})))
 
%  plot forcing functions

subplot(1,1,1)
plot(fdatime, 1e9*resmat, '-', ...
     fdatime, 1e9*D3fdameanmat(:,2), 'r:')
xlabel('\fontsize{16} Milliseconds')
ylabel('\fontsize{16} Meters/sec/sec/sec')
axis([0,2300,-200,200])

%  plot the mean forcing function along with second deriv.

resmeanfd  = mean(resfd);
resmeanvec = eval_fd(fdatime, resmeanfd);

plot(fdatime, 1e9*resmeanvec, 'b-', ...
     fdatime, 1e9*D3fdameanmat(:,2), 'r:')
xlabel('\fontsize{16} Milliseconds')
ylabel('\fontsize{16} Meters/sec/sec/sec')
axis([0,2300,-200,200])

% Define a functional data object for the 
%  three derivative weight functions

wcoef1 = getcoef(getfd(bestwtcell{1}));
wcoef2 = getcoef(getfd(bestwtcell{2}));
wcoef3 = getcoef(getfd(bestwtcell{3}));

wcoef  = [wcoef1, wcoef2, wcoef3];
wfd    = fd(wcoef,wbasis);

%  Set up a linear differential operator.
%  This isn't used in these analyses.

fdaLfd = Lfd(difeorder, fd2cell(wfd));

%  solve equation

ystart = eye(3);
ystart(1,1) =   fdameanmat(1,2);
ystart(2,2) = D1fdameanmat(1,2);
ystart(3,3) = D2fdameanmat(1,2);

odeoptions = odeset('RelTol', 1e-6);

[tp1, yp1] = ode45(@derivs, fdatime, ystart(:,1), odeoptions, wfd);
[tp2, yp2] = ode45(@derivs, fdatime, ystart(:,2), odeoptions, wfd);
[tp3, yp3] = ode45(@derivs, fdatime, ystart(:,3), odeoptions, wfd);

%  plot the three solutions

umaty   = [yp1(:,1),yp2(:,1),yp3(:,1)];
Dumaty  = [yp1(:,2),yp2(:,2),yp3(:,2)];
D2umaty = [yp1(:,3),yp2(:,3),yp3(:,3)];

subplot(3,1,1)
plot(fdatime, umaty, '-', [0,2300], [0,0], 'r:');
title('Function')

subplot(3,1,2)
plot(fdatime, Dumaty, [0,2300], [0,0], 'r:')
title('First Derivative');

subplot(3,1,3)
plot(fdatime, D2umaty, [0,2300], [0,0], 'r:')
title('Second Derivative');

%  plot fit to each curve ... hit any key after each plot

subplot(1,1,1)
index  = 1:20;
fdamat = eval_fd(fdatime, fdafd);
zmat   = [ones(1401,1),fdatime-1150,umaty];
for i = index
   yhat = zmat * (zmat\fdamat(:,i,2));
   lhdl = plot(fdatime, yhat, '-', ...
               fdatime, fdamat(:,i,2), '--');
   set(lhdl, 'LineWidth', 2)
   title(['X curve ',num2str(i)])
   axis([0, 2300, -0.04, 0.04])
   legend('\fontsize{13} Fit', '\fontsize{13} Observed', ...
          'Location', 'SouthEast')
   pause;
end

%  plot fit to each script ... hit any key after each plot

subplot(1,1,1)
index  = 1:20;
fdamat = eval_fd(fdatime, fdafd);
zmatx   = [ones(1401,1),fdatime-1150,umatx];
zmaty   = [ones(1401,1),fdatime-1150,umaty];
for i = index
   xhat = zmatx * (zmatx\fdamat(:,i,1));
   yhat = zmaty * (zmaty\fdamat(:,i,2));
   lhdl = plot(xhat, yhat, '-', ...
               fdamat(:,i,1), fdamat(:,i,2), '--');
   set(lhdl, 'LineWidth', 2)
   xlabel('\fontsize{16} X')
   ylabel('\fontsize{16} Y')
   title(['\fontsize{16} Script ',num2str(i)])
   axis([-0.05,0.05,-0.05,0.05])
   legend('\fontsize{13} Fit', '\fontsize{13} Observed', ...
          'Location', 'SouthEast')
   pause;
end

%  plot the two weight functions for the second derivative

subplot(1,1,1)
lhdl = plot(fdatime, b2vecX, 'b-', fdatime, b2vecY, 'g-', ...
            [0, 2300], [b2meanX, b2meanX], 'b--', ...
            [0, 2300], [b2meanY, b2meanY], 'g--');
set(lhdl, 'LineWidth', 2);
hold on
plotrange = [0,6e-3];
for k=1:length(cycle)
    lhdl = plot([cycle(k),cycle(k)],plotrange,'r:');
    set(lhdl, 'LineWidth', 1')
end
for j=1:length(featureindex)
    k = featureindex(j);
    lhdl = plot([cycle(k),cycle(k)],plotrange,'m-');
    set(lhdl, 'LineWidth', 2')
end
hold off
axis([0, 2300, 0, 6e-3])

%  Conclusions:

%  For each coordinate, a single differential equation was
%  estimated to describe all 20 coordinate records 
%  simulataneously.
%  The differential equation solutions are readable.
%  Moreover, they accommodate a good deal of the curve-to-curve
%  variation in the scripts.  
%  The cusp first in the "a" is not handled as well as the other cusps.
%  This may be due to poor registration at this point, or to the fact
%  that this cusp has a lot of variation from one replication to another.