handwrit.m

手写签名数据

%  -----------------------------------------------------------------------
%                 Registered Handwriting Data
%  -----------------------------------------------------------------------

% These data are the X-Y coordinates of 20 replications of writing
% the script "fda".  The subject was Jim Ramsay.  Each replication
% is represented by 1401 coordinate values.  The scripts have been 
% extensively pre-processed.  They have been adjusted to a common
% length that corresponds to 2.3 seconds or 2300 milliseconds, and
% they have already been registered so that important features in
% each script are aligned.
% 
% This analysis is designed to illustrate techniques for working
% with functional data having rather high frequency variation and
% represented by thousands of data points per record.  Comments
% along the way explain the choices of analysis that were made.
% 
% The final result of the analysis is a third order linear 
% differential equation for each coordinate forced by a 
% constant and by time.  The equations are able to reconstruct
% the scripts to a fairly high level of accuracy, and are also
% able to accommodate a substantial amount of the variation in
% the observed scripts across replications.  by contrast, a 
% second order equation was found to be completely inadequate.
% 
% An interesting suprise in the results is the role placed by
% a 120 millisecond cycle such that sharp features such as cusps
% correspond closely to this period.  This 110-120 msec cycle
% seems is usually seen in human movement data involving rapid
% movements, such as speech, juggling and so on.

%  Last modified 17 June 2008

% Add paths to functional data functions and to the data

addpath ('c:\Program Files\matlab\fdaM')
addpath ('c:\Program Files\matlab\fdaM\examples\handwrit')

%  Input the data.  These 20 records have already been
%  normalized to a common time interval of 2300 milliseconds
%  and have beeen also registered so that prominent features
%  occur at the same times across replications.
%  Time will be measured in milliseconds and space in meters.
%  The data will require a small amount of smoothing, since
%  an error of 0.5 mm is characteristic of the OPTOTRAK 3D
%  measurement system used to collect the data.

fid      = fopen('fdareg.dat','rt');
fdaarray = reshape(fscanf(fid,'%f'), [20,2,1401]);
fdaarray = permute(fdaarray,[3,1,2]);
fdaarray = fdaarray/1000;   %  convert spatial unit to meters

%  Set up time values and range.
%  It is best to choose milliseconds as a time scale
%  in order to make the ratio of the time
%  unit to the inter-knot interval not too
%  far from one.  Otherwise, smoothing parameter values
%  may be extremely small or extremely large.

fdatime  = linspace(0, 2300, 1401)';
fdarange = [0, 2300];

%  The basis functions will be B-splines, with a spline 
%  placed at each knot.  One may question whether so many
%  basis functions are required, but this decision is found to 
%  be  essential for stable derivative estimation up to the 
%  third order at and near the boundaries.

%  Order 7 was used to get a smooth third derivative, which
%  requires penalizing the size of the 5th derivative, which
%  in turn requires an order of at least 7.
%  This implies norder + no. of interior knots = 1399 + 7 = 1406 
%  basis functions.  

nbasis   = 1406;
norder   =    7;
fdabasis = create_bspline_basis(fdarange, nbasis, norder);

%  The smoothing parameter value 1e8 was chosen to obtain a
%  fitting error of about 0.5 mm, the known error level in
%  the OPTOTRACK equipment.

fdafd  = fd(zeros(nbasis,20,2), fdabasis);
lambda = 1e9;
fdaPar = fdPar(fdafd, 5, lambda);

%  set up the functional data structure 

[fdafd, df, gcv] = smooth_basis(fdatime, fdaarray, fdaPar);

%  Add suitable names for the dimensions of the data.

fdafd_fdnames{1} = 'Milliseconds';
fdafd_fdnames{2} = 'Replications';
fdafd_fdnames{3} = 'Metres';
fdafd = putnames(fdafd, fdafd_fdnames);

%  display degrees of freedom and total GCV criterion

disp(['degrees of freedom = ',num2str(df)])  %  about 92
totalgcv = sum(sum(gcv));
disp(['total GCV = ',num2str(totalgcv)])  
RMSgcv = sqrt(totalgcv)*1000; % about 0.3 mm
disp(['RMS GCV = ',num2str(RMSgcv)])  

%  plot the fit to the data.  It's good to within plotting
%  accuracy, but the title displays RMSE values

plotfit_fd(fdaarray, fdatime, fdafd);

%  plot all curves

plot(fdafd)

%  compute values of curves and the values of the curve

fdameanfd  = mean(fdafd);
fdamat     = eval_fd(fdatime, fdafd);
fdameanmat = squeeze(eval_fd(fdatime, fdameanfd));

%  Evaluate the three derivatives and their means

D1fdamat = eval_fd(fdatime, fdafd, 1);
D2fdamat = eval_fd(fdatime, fdafd, 2);
D3fdamat = eval_fd(fdatime, fdafd, 3);

D1fdameanmat = eval_fd(fdatime, fdameanfd, 1);
D2fdameanmat = eval_fd(fdatime, fdameanfd, 2);
D3fdameanmat = eval_fd(fdatime, fdameanfd, 3);

%  -----------------------------------------------------
%            Principal Components Analysis
%  -----------------------------------------------------

%  Set up a more compact basis using the df identified in
%  the spline smoothing above.  Use order 6 splines

nbasis92   = 92;
norder92   =  6;
fdabasis92 = create_bspline_basis(fdarange, nbasis92, norder92);

%  set up the functional data structure 

fdafd92 = smooth_basis(fdatime, fdaarray, fdabasis92);
%  Add suitable names for the dimensions of the data.
fdafd92_fdnames{1} = 'Milliseconds';
fdafd92_fdnames{2} = 'Replications';
fdafd92_fdnames{3} = 'Metres';
fdafd92 = putnames(fdafd92, fdafd92_fdnames);

%  principal components analysis

lambda    = 0;
fdPar92   = fdPar(fdabasis92, 4, lambda);
nharm     = 3;
fdapcastr = pca_fd(fdafd92, nharm, fdPar92);

%  plot unrotated harmonics

subplot(1,1,1)
plot_pca(fdapcastr)

%  Varimax rotation

fdarotpcastr = varmx_pca(fdapcastr);

%  plot rotated harmonics

subplot(1,1,1)
plot_pca(fdarotpcastr)

fdameanfd  = mean(fdafd92);
fdameanmat = squeeze(eval_fd(fdatime, fdameanfd));

harmfd  = fdarotpcastr.harmfd;
harmmat = eval_fd(fdatime, harmfd);

nfine = 201;

fdapointtime = linspace(0,2300,nfine)';
fdameanpoint = squeeze(eval_fd(fdapointtime, fdameanfd));
harmpointmat = squeeze(eval_fd(fdapointtime, harmfd));

fac = 0.1;
harmplusmat = zeros(nfine,3,2);
harmminsmat = zeros(nfine,3,2);
for j=1:3
    harmplusmat(:,j,:) = fdameanpoint + fac.*squeeze(harmpointmat(:,j,:));
    harmminsmat(:,j,:) = fdameanpoint - fac.*squeeze(harmpointmat(:,j,:));
end

for j=1:2
%     figure(j)
    subplot(1,2,j)
    phdl = plot(fdameanmat(:,1), fdameanmat(:,2), 'b-');
    set(phdl, 'LineWidth', 2)
%     text(harmplusmat(:,j,1), harmplusmat(:,j,2), '+');
%     text(harmminsmat(:,j,1), harmminsmat(:,j,2), '-');
    hold on
    plot(harmplusmat(:,j,1), harmplusmat(:,j,2), 'b--');
    plot(harmminsmat(:,j,1), harmminsmat(:,j,2), 'b--');
    hold off
    axis('square')
end

%  plot log eigenvalues

fdaharmeigval = fdapcastr.values;
x = ones(16,2);
x(:,2) = reshape((4:19),[16,1]);
y = log10(fdaharmeigval(4:19));
c = x\y;
subplot(1,1,1)
plot(1:19,log10(fdaharmeigval(1:19)),'-o', ...
     1:19, c(1)+ c(2).*(1:19), ':')
xlabel('Eigenvalue Number')
ylabel('Log10 Eigenvalue')

%  plot factor scores

harmscr = fdarotpcastr.harmscr;

plot(harmscr(:,1), harmscr(:,2), 'o')
xlabel('Harmonic I')
ylabel('Harmonic II')

%  -----------------------------------------------------
%  Display the neural clock cycle of 119 milliseconds
%  -----------------------------------------------------

%  Set up motor control clock cycle times at every 
%  119 milliseconds. 

cyclelength = 117;
cycle  = 0:cyclelength:2300;
ncycle = length(cycle);

%  evaluate curves at cycle times

fdamatcycle     = eval_fd(cycle, fdafd);
fdameanmatcycle = squeeze(eval_fd(cycle, fdameanfd));

%  Indices of cycle times corresponding to important features:
%  -- the cusp in "f", 
%  -- the the cusp in "d", 
%  -- the first cusp in "a", 
%  -- the rest after the first cusp in "a", and 
%  -- the second cusp in "a".
%  It is remarkable that these features correspond so closely
%  with clock cycle times!

featureindex   = [3, 5, 7, 10, 13, 16, 19];
fdafeature     = fdamatcycle(featureindex,:,:);
fdameanfeature = fdameanmatcycle(featureindex,:);

%  Plot mean, including both sampling points and fit
%  Points at cycle times are plotted as blue circles, and
%  points at feature times are plotted as red circles.

subplot(1,1,1);
lhdl = plot(fdameanmat(:,1),      fdameanmat(:,2),      'b-', ...
            fdameanmatcycle(:,1), fdameanmatcycle(:,2), 'bo', ...
            fdameanfeature(:,1),  fdameanfeature(:,2),  'ro');
set(lhdl, 'LineWidth', 2);
xlabel('\fontsize{16} Metres')
ylabel('\fontsize{16} Metres')
axis([-.040, .040,-.040, .040]);
title('\fontsize{16} Mean script');

%  Plot individual curves, including both sampling points and fit
%  also plot the mean curve in the background.
%  Note how closely individual curve features are tied to the
%  feature cycle times.

subplot(1,1,1);
for i = 1:20
  lhdl = plot(fdamat(:,i,1),        fdamat(:,i,2),        'b-', ...
              fdamatcycle(:,i,1),   fdamatcycle(:,i,2),   'bo', ...
              fdafeature(:,i,1),    fdafeature(:,i,2),    'go', ...
              fdameanmat(:,1),      fdameanmat(:,2),      'r--', ...
              fdameanmatcycle(:,1), fdameanmatcycle(:,2), 'ro', ...
              fdameanfeature(:,1),  fdameanfeature(:,2),  'mo');
  set(lhdl, 'LineWidth', 2);
  xlabel('\fontsize{16} Metres')
  ylabel('\fontsize{16} Metres')
  axis([-.040, .040,-.040, .040]);
  title(['Record ', num2str(i)]);
  pause;
end

%  Plot the individual acceleration records.
%  In these plots, acceleration is displayed as
%  metres per second per second.

%  Cycle and feature times are plotted as vertical
%  dashed lines, un-featured cycle times as red 
%  dotted lines, and cycle times of features as 
%  heavier magenta solid lines.

subplot(1,1,1);
for i=1:20
    lhdl = plot(fdatime, 1e6*D2fdamat(:,i,1), '-', ...
        fdatime, 1e6*D2fdamat(:,i,2), '-', ...
        [0,2300],[0,0], 'r:');
    set(lhdl, 'Linewidth', 2)
    hold on
    plotrange = [-12,12];
    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} Milliseconds')
    ylabel('\fontsize{16} Meters/msec/msec')
    title(['\fontsize{16} Curve ',num2str(i)])
    legend('X', 'Y')
    axis([0, 2300, -12, 12]);
    pause
end

%  Compute and plot the acceleration magnitudes,  
%  also called the tangential accelerations.  

D2mag  = sqrt(D2fdamat(:,:,1).^2 + D2fdamat(:,:,2).^2);
D2magmean = mean(D2mag,2);

subplot(1,1,1);
plot(fdatime, 1e6*D2mag)
hold on
plotrange = [0,12];
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} Milliseconds')
ylabel('\fontsize{16} Metres/sec/sec')
title('\fontsize{16} Acceleration Magnitude');
axis([0,2300,0,12]);

%  Plot the mean acceleration magnitude as well as
%  those for each curve.
%  Note the two rest cycles, one in "d" and one in "a"

subplot(1,1,1);
lhdl = plot(fdatime, 1e6*D2magmean);
set(lhdl, 'LineWidth', 2')
hold on
plotrange = [0,8];
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} Milliseconds')
ylabel('\fontsize{16} Metres/sec/sec')
title('\fontsize{16} Mean acceleration Magnitude');
axis([0,2300,0,8]);

%  Plot each individual acceleration magnitude, along
%  with the mean magnitude as a green dashed line

subplot(1,1,1)
for i=1:20
    lhdl = plot(fdatime, 1e6*D2mag(:,i), '-', ...
        fdatime, 1e6*D2magmean, '--');
    set(lhdl, 'LineWidth', 2)
    hold on
    plotrange = [0,12];
    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} Milliseconds')
    ylabel('\fontsize{16} Metres/sec/sec')
    title(['\fontsize{16} Script ',num2str(i)])
    axis([0,2300,0,12])
    pause
end

%  ------------------------------------------------------------
%                 Principal Differential Analysis 
%  A third order equation forced by a constant and time is
%          estimated for X and Y coordinates separately.
%  Forcing with constant and time is required to allow for 
%  an arbitrary origin and the left-right motion in X.
%  ------------------------------------------------------------

difeorder  = 3;  %  order of equation

%  set up the two forcing functions

%  constant forcing
constbasis = create_constant_basis(fdarange);
constfd    = fd(ones(1,20), constbasis);
ufdcell{1} = constfd;               
% time forcing
linbasis   = create_monomial_basis(fdarange, 2);
lincoef    = zeros(2,20);
lincoef(2,:) = 1;
ufdcell{2} = fd(lincoef, linbasis); 

%  set up the corresponding weight functions

constfd    = fd(1, constbasis);
constfdPar = fdPar(constfd);
awtcell{1} = constfdPar;
awtcell{2} = constfdPar;

%  Define two basis systems for the derivative weight
%  functions.  One for a background analysis used as a
%  baseline and using constant weight functions, and
%  another using 125 basis functions that will be used
%  to generate the equaations.

%  Set the number of basis functions to match the 
%  119 msec clock time estimated above, so that the
%  gap between knots is 11.9 msec.

delta   = 11.9;
breaks  = [0:delta:2300, 2300];
nbreaks = length(breaks);
nwbasis = nbreaks + 2;
wbasis  = create_bspline_basis(fdarange, nwbasis);

%  ------------------------------------------------------------
%                   Analysis for coordinate X
%  ------------------------------------------------------------

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

xbasis = create_bspline_basis(fdarange, 245);
xfdcell{1} = data2fd(squeeze(fdamat(:,:,1)), fdatime, xbasis);

%  Set up the derivative weight functions