Revision fd92910f232786c2531c1c74e02ece6b49a2f76f authored by J. O. Ramsay on 03 November 2009, 00:00:00 UTC, committed by Gabor Csardi on 03 November 2009, 00:00:00 UTC
1 parent 8ce7a35
fdarm-ch06.m
%
% Ramsay, Hooker & Graves (2009)
% Functional Data Analysis with R and Matlab (Springer)
% Remarks and disclaimers
% These R commands are either those in this book, or designed to
% otherwise illustrate how R can be used in the analysis of functional
% data.
% We do not claim to reproduce the results in the book exactly by these
% commands for various reasons, including:
% -- the analyses used to produce the book may not have been
% entirely correct, possibly due to coding and accuracy issues
% in the functions themselves
% -- we may have changed our minds about how these analyses should be
% done since, and we want to suggest better ways
% -- the R language changes with each release of the base system, and
% certainly the functional data analysis functions change as well
% -- we might choose to offer new analyses from time to time by
% augmenting those in the book
% -- many illustrations in the book were produced using Matlab, which
% inevitably can imply slightly different results and graphical
% displays
% -- we may have changed our minds about variable names. For example,
% we now prefer "yearRng" to "dayrange" for the weather data.
% -- three of us wrote the book, and the person preparing these scripts
% might not be the person who wrote the text
% Moreover, we expect to augment and modify these command scripts from time
% to time as we get new data illustrating new things, add functionality
% to the package, or just for fun.
%
% ch. 6. Descriptions of Functional Data
%
% Set up some strings for constructing paths to folders.
% These strings should be modified so as to provided access
% to the specified folders on your computer.
% Path to the folder containing the Matlab functional data analysis
% software
fdaMPath = 'c:/Program Files/MATLAB/R2009a/fdaM';
addpath(fdaMPath)
% Path to the folder containing the examples
examplesPath = [fdaMPath,'/examples'];
addpath(examplesPath)
%
% Section 6.1 Some Functional Descriptive Statistics
%
% set up a functional data object for log precipitation
% ---------- Statistics for the log precipitation data ---------------
% path to the daily weather data example folder
weatherPath = [examplesPath,'/weather'];
addpath(weatherPath)
load daily
tempav = daily.tempav;
precav = daily.precav;
yearRng = daily.rng;
time = daily.time;
place = daily.place;
dayperiod = daily.period;
% define the harmonic acceleration operator
Lcoef = [0,(2*pi/dayperiod)^2,0];
harmaccelLfd = vec2Lfd(Lcoef, yearRng);
% organize data to have winter in the center of the plot
dayOfYearShifted = [182:365, 1:181];
% change 0's to 0.05 mm in precipitation data
prectmp = precav;
for j=1:35
index = prectmp(:,j)==0;
prectmp(index,j) = 0.05;
end
% set up functional data object for log precipitation
logprecav = log10(prectmp(dayOfYearShifted,:));
% set up a saturated basis: as many basis functions as observations
daybasis = create_fourier_basis(yearRng, 365);
% smooth data with lambda that minimizes GCV
lambda = 1e6;
fdParobj = fdPar(daybasis, harmaccelLfd, lambda);
logprec_fd = smooth_basis(time, logprecav, fdParobj);
fdnames{1} = 'Day (July 1 to June 30)';
fdnames{2,1} = 'Weather Station';
fdnames{2,2} = place;
fdnames{3} = 'Log 10 Precipitation (mm)';
logprec_fd = putnames(logprec_fd, fdnames);
% elementary pointwise mean and standard deviation
meanlogprec = mean(logprec_fd);
stddevlogprec = std_fd(logprec_fd);
% Section 6.1.1 The Bivariate Covariance Function v(s; t)
logprecvar_bifd = var_fd(logprec_fd);
weektime = linspace(0,365,53);
logprecvar_mat = eval_bifd(weektime, weektime, logprecvar_bifd);
% Figure 6.1
surf(weektime, weektime, logprecvar_mat)
xlabel('\fontsize{13} Day (July 1 to June 30)')
ylabel('\fontsize{13} Day (July 1 to June 30)')
zlabel('\fontsize{13} variance(log10 precip)')
contour(weektime, weektime, logprecvar_mat);
xlabel('\fontsize{13} Day (July 1 to June 30)')
ylabel('\fontsize{13} Day (July 1 to June 30)')
% Figure 6.2
% Matlab doesn't indicate contour level values in the plot,
% but instead uses an optional colorbar to indicate height.
contour(weektime, weektime, logprecvar_mat);
xlabel('\fontsize{13} Day (July 1 to June 30)')
ylabel('\fontsize{13} Day (July 1 to June 30)')
colorbar
%
% Section 6.2 The Residual Variance-Covariance Matrix Se
%
% (no computations in this section)
%
% Section 6.3 Functional Probes rho[xi]
%
% see section 6.5 below
%
% Section 6.4 Phase-plane Plots of Periodic Effects
%
% ----------- Phase-plane plots for the nondurable goods data ---------
% path to the nondurable goods index example folder
goodsindexPath = [examplesPath,'/goodsindex'];
addpath(goodsindexPath)
% load the goodsindex data from the .mat file goodsindex.mat
load goodsindex
ndur = goodsindex.ndur;
yeartime = goodsindex.yeartime;
nondurables = goodsindex.nondurables;
lognondur = goodsindex.lognondur;
% set up a basis for smoothing log nondurable goods index
goodsRng = [1919,2000];
nbasis = 979;
norder = 8;
goodsbasis = create_bspline_basis(goodsRng, nbasis, norder);
% set up fdPar object
LfdobjNonDur= int2Lfd(4);
lambda = 1e-11;
goodsLfdPar = fdPar(goodsbasis, LfdobjNonDur, lambda);
logNondur_fd = smooth_basis(yeartime, lognondur, goodsLfdPar);
% Fig. 6.3 The log nondurable goods index for 1964 to 1967
t64_67 = linspace(1964, 1967, 601);
logNondur_vec = eval_fd(t64_67, logNondur_fd);
sel64_67 = find(1964 <= yeartime & yeartime <= 1967);
ymin = min(lognondur(sel64_67));
ymax = max(lognondur(sel64_67));
phdl=plot(yeartime(sel64_67), lognondur(sel64_67), 'bo', ...
t64_67, logNondur_vec, 'b-', ...
[1965,1965], [ymin,ymax], 'b--', [1966,1966], [ymin,ymax], 'b--');
set(phdl, 'LineWidth', 2)
xlabel('\fontsize{13} Year')
ylabel('\fontsize{13} Log_{10} Nondurable Goods Index')
axis([1964,1967,1.615,1.725])
% Section 6.4.1. Phase-plane Plots Show Energy Transfer
% Figure 6.4. Phase-plane plot for a simple harmonic function
tvec = linspace(0,1,101);
sinvec = sin(2*pi*tvec);
cosvec = cos(2*pi*tvec);
Dsinvec = 2*pi*cosvec;
D2sinvec = -(2*pi).^2*sinvec;
phdl=plot(Dsinvec, D2sinvec, 'b-', ...
[min(Dsinvec),max(Dsinvec)], [0,0], 'b:', ...
[0,0], [min(D2sinvec),max(D2sinvec)], 'b:');
set(phdl, 'LineWidth', 2)
xlabel('\fontsize{13} Velocity')
ylabel('\fontsize{13} Acceleration')
text(-4, 45, '\fontsize{12} Max. potential energy')
text(-4, -45, '\fontsize{12} Max. potential energy')
text(-9.5, 5, '\fontsize{12} Max.')
text(-9.5, 0, '\fontsize{12} kinetic')
text(-9.5, -5, '\fontsize{12} energy')
text( 6.5, 5, '\fontsize{12} Max.')
text( 6.5, 0, '\fontsize{12} kinetic')
text( 6.5, -5, '\fontsize{12} energy')
axis([-10,10,-50,50])
% Section 6.4.2 The Nondurable Goods Cycles
% ----------- Phase-plane plots for the nondurable goods data ---------
% Figure 6.5
startyr = 64; % starting year
yearindex = startyr + 1900 - 1919; % indices of year
% Set up the derivatives for phase-plane plotting.
index = (1:(13)) + yearindex*12;
durrange = [min(yeartime(index)), max(yeartime(index))];
durfine = linspace(durrange(1),durrange(2),401)';
durcrse = linspace(durrange(1),durrange(2),13)';
D1f = eval_fd(durfine, logNondur_fd, 1);
D2f = eval_fd(durfine, logNondur_fd, 2);
D1c = eval_fd(durcrse, logNondur_fd, 1);
D2c = eval_fd(durcrse, logNondur_fd, 2);
% The phase-plane plotting.
monthlabs = ['jan';'F ';'mar';'A ';'M ';'Jun'; ...
'Jly';'Aug';'S ';'Oct';'N ';'D '];
plot([0,0], [-12,12], 'b:', [-.75,.75], [0,0], 'b:')
x0 = durrange(1);
hold on
indj = find(durfine >= x0 & durfine <= x0+1);
phdl = plot(D1f(indj), D2f(indj), 'b-');
set(phdl, 'LineWidth', 2);
indexc = 1:12;
for k=1:12
plot(D1c(indexc(k)), D2c(indexc(k)), 'o')
text(D1c(indexc(k))+0.05, D2c(indexc(k))+0.5, monthlabs(k,:))
end
hold off
axis([-.5, .5, -12, 12])
xlabel('\fontsize{13} Velocity')
ylabel('\fontsize{13} Acceleration')
% sec. 6.4.3. Phase-Plane Plotting the Growth of Girls
% ------------- Phase-plane diagrams for the growth data ---------------
% path to the Berkeley growth data example folder
growthPath = [examplesPath,'/growth'];
addpath(growthPath)
% load the growth data struct object from the .mat file growth.mat
load growth
% the dimensions of the data
ncasem = growth.ncasem;
ncasef = growth.ncasef;
nage = growth.nage;
% the heights of the 54 girls
heightmat = growth.hgtfmat;
% the 31 ages of measurement
age = growth.age;
% define the range of the ages and set up a fine mesh of ages
ageRng = [1,18];
norder = 6;
breaks = age;
nbasis = length(breaks) + norder - 2;
gr_basis = create_bspline_basis(ageRng, nbasis, norder, breaks);
children = 1:10;
coef = zeros(nbasis,10);
Wfd0 = fd(coef, gr_basis);
gr_Lfd = 3;
gr_lambda = 10^(-1.5);
gr_fdPar1_5 = fdPar(Wfd0, gr_Lfd, gr_lambda);
[Wfd, beta] = smooth_monotone(age, heightmat(:,children), gr_fdPar1_5);
agefine = linspace(1,18,101);
i11_7 = find(abs(agefine-11.7) == min(abs(agefine-11.7)));
betamat = repmat(beta(2,:),101,1);
velffine = betamat.*eval_monfd(Wfd, agefine, 1);
accffine = betamat.*eval_monfd(Wfd, agefine, 2);
phdl=plot(velffine(:,1), accffine(:,1), 'b-', ...
velffine(i11_7, 1), accffine(i11_7, 1), 'bo', ...
[0,12], [0,0], 'b:');
set(phdl, 'LineWidth', 2);
hold on
for i = 2:10
phdl=plot(velffine(:,i), accffine(:,i), 'b-', ...
velffine(i11_7, i), accffine(i11_7, i), 'bo');
set(phdl, 'LineWidth', 2);
end
hold off
xlabel('\fontsize{13} Velocity (cm/yr)')
ylabel('\fontsize{13} Acceleration (cm/yr^2)')
axis([0,12,-6,3])
%
% Section 6.5 Confidence Intervals for Curves and their Derivatives
%
% sec. 6.5.1. Two Linear mappings Defining a Probe Value
% ----------------- Probe values for Canadian weather data ------------
% define a probe to emphasize mid-winter
dayvec = linspace(0,365,101);
xivec = exp(20*cos(2*pi*(dayvec-197)/365));
xibasis = create_bspline_basis(yearRng,13);
xifd = smooth_basis(dayvec, xivec, xibasis);
plot(xifd)
% define bases for temperature and precipitation
tempbasis = create_fourier_basis(yearRng, 65);
precbasis = create_fourier_basis(yearRng,365);
% probe values for basis functions with respect to xifd
tempLmat = inprod(tempbasis, xifd);
precLmat = inprod(precbasis, xifd);
% sec. 6.5.3. Confidence Limits for Prince Rupert's Log Precipitation
% smooth data with lambda that minimizes GCV getting
% all of the output up to matrix y2cMap
lambda = 1e6;
fdParobj = fdPar(daybasis, harmaccelLfd, lambda);
[logprec_fd, df, gcv, SSE, penmat, y2cMap] = ...
smooth_basis(time, logprecav, fdParobj);
fdnames{1} = 'Day (July 1 to June 30)';
fdnames{2,1} = 'Weather Station';
fdnames{2,2} = place;
fdnames{3} = 'Log 10 Precipitation (mm)';
logprec_fd = putnames(logprec_fd, fdnames);
% compute the residual matrix and variance vector
logprecmat = eval_fd(time, logprec_fd);
logprecres = logprecav - logprecmat;
logprecvar = sum(logprecres.^2, 2)/(35-1);
% smooth log variance vector
lambda = 1e8;
resfdParobj = fdPar(daybasis, harmaccelLfd, lambda);
logvar_fd = smooth_basis(time, log(logprecvar), resfdParobj);
% evaluate the exponentiated log variance vector and
% set up diagonal error variance matrix SigmaE
varvec = exp(eval_fd(time, logvar_fd));
SigmaE = diag(varvec);
% compute variance covariance matrix for fit
c2rMap = eval_basis(time, daybasis);
Sigmayhat = c2rMap * y2cMap * SigmaE * y2cMap' * c2rMap';
% extract standard error function for yhat
logprec_stderr= sqrt(diag(Sigmayhat));
% plot Figure 6.6
logprecvec29 = eval_fd(time,logprec_fd(29));
phdl = plot(time, logprecvec29, 'b-', ...
time, logprec29 + 2*logprec_stderr, 'b--', ...
time, logprec29 - 2*logprec_stderr, 'b--', ...
time, logprecav(:,29), 'bo');
set(phdl, 'LineWidth', 2);
xlabel('\fontsize{13} Day (July 1 to June 30)')
ylabel('\fontsize{13} Log 10 Precipitation (mm)')
title('\fontsize{13} Prince Rupert')
axis([0,365,0.2,1.3])
%
% Section 6.6 Some Things to Try
%
% (exercises for the reader)
![swh spinner](/static/img/swh-spinner.gif)
Computing file changes ...