Contents
Introducing the FDA Toolbox
Disclaimer -- the FDA toolbox is not an official Matlab toolbox; it is a collection of privately-written functions and comes with no garuantees.
To get the toolbox
1. Visit functionaldata.org
2. Click on software and follow the link at the bottom of the page
3. Alternatively, go straight to
<ftp://ego.psych.mcgill.ca/pub/ramsay/FDAfuns/Matlab>
4. Dowload Matlabfunctions.zip
5. Follow the instruction in INSTALL.MATLAB.WIN.txt
Linux users; the installation is analagous
I keep my fda package in
F://work/External/matlab_packages/fdaM
Now in order to make use of these functions, you need to add them to your path:
addpath('C://work/External/matlab_packages/fdaM') % Alternatively, you can add the directory to your path permanently using % the menus: % % File -> Set Path ... % % and using the "Add Folder" option -- you do not need to select "add % subfolders". % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Now let's load some data load 'handwriting.mat' % Now we would like to start smoothing the data
Warning: Name is nonexistent or not a directory: C:\work\External\matlab_packages\fdaM.
Defining a Basis
help basis % All basis objects require a range to be specified: basis_range = [min(fdatime) max(fdatime)]; % Starting with a B-spline basis: breaks = linspace(basis_range(1),basis_range(2),11); figure(1) plot(breaks,'.') % I'll choose cubic B-splines norder = 4; % And that gives nbasis = length(breaks) + norder - 2; % This is enough to define the basis bspline_obj = create_bspline_basis(basis_range,nbasis,norder,breaks); % Now that I've got a basis, I can evaluate it: basisvals = eval_basis(fdatime,bspline_obj); size(basisvals) figure(2) plot(fdatime,basisvals) % There is also a shortcut figure(3) plot(bspline_obj) % Lets look at just one figure(4) plot(fdatime,basisvals(:,4),'linewidth',2) hold on plot([breaks; breaks],[-ones(1,length(breaks)); ones(1,length(breaks))],'r--') % Now I can evaluate some derivatives; the third argument in eval_basis % gives an order of derivative to evaluate. bbasisvals = eval_basis(fdatime,bspline_obj,1); b2basisvals = eval_basis(fdatime,bspline_obj,2); plot(fdatime,bbasisvals(:,7)/3,'m','linewidth',2) plot(fdatime,b2basisvals(:,10)/40,'k','linewidth',2) hold off % Here I have had to scale the derivatives in order to make them comparable % in size. This is important numerically -- using a very fine basis can % mean that taking derivatives leads to numerical overflow or underflow. % To give a comparison, let's try a linear basis norder = 2; nbasis = length(breaks); bspline_obj = create_bspline_basis(basis_range,nbasis,norder,breaks); basisvals = eval_basis(fdatime,bspline_obj); figure(5) plot(fdatime,basisvals) axis([0 2.3 0 2]) % There is probably not enough resolving power in this basis, so I'll try % more knots norder = 4; breaks = linspace(basis_range(1),basis_range(2),101); nbasis = length(breaks) + norder - 2; bspline_obj = create_bspline_basis(basis_range,nbasis,norder,breaks); basisvals = eval_basis(fdatime,bspline_obj); figure(6) plot(fdatime,basisvals) axis([0 2.3 0 1]) % Lets try a fourier basis: nbasis_f = 10; period = basis_range(2)-basis_range(1); fourier_obj = create_fourier_basis(basis_range,nbasis_f,period); fbasisvals = eval_basis(fdatime,fourier_obj); figure(7) plot(fdatime,fbasisvals(:,1:5)); % And maybe a polynomial basis on a centered range nbasis_p = 10; exponents = 0:9; % Note that I do not need to use integer exponents here monomial_obj = create_monomial_basis(basis_range,nbasis_p,exponents); mbasisvals = eval_basis(fdatime,monomial_obj); figure(8) plot(fdatime,mbasisvals) % Accessing aspects of a basis object getbasisrange(bspline_obj) % range of the basis getnbasis(bspline_obj) % number of basis functions bvals = getbasismatrix(fdatime,fourier_obj); % same as eval_basis size(bvals)
BASIS Creates a functional data basis.
CREATE_BASIS An alternative call to CREATE_BASIS_fd.
DISPLAY Display a functional data basis object.
EQ assesses whether two bases are equivalent.
GETBASISMATRIX Computes the basis matrix evaluated at arguments in
GETBASISPAR Extracts the basis parameters from basis object BASISOBJ.
GETBASISRANGE Extracts the range from basis object BASISOBJ.
GETBASISTYPE Extracts the type of basis from basis object BASISOBJ.
GETDROPIND Extracts the indices of basis functions to be dropped
GETNBASIS Extracts the type of basis from basis object BASISOBJ.
GETQUADVALS Extracts the quadrature points and weights
Return values of the derivative of order NDERIV of basis
ISA_BASIS checks a struct object for fields for basis objects
Plot a basis object.
PUQUADVALS Enters drop indices for
PUQUADVALS Enters the quadrature points and weights
PUTVALUES Enters a cell array of values of basis functions
TIMES for two basis objects sets up a basis suitable for
basis is both a directory and a function.
BASIS Creates a functional data basis.
Arguments:
BASISTYPE ... a string indicating the type of basis. This may be one of
'Fourier', 'fourier', 'Fou', 'fou',
'Bspline', 'bspline', 'Bsp', 'bsp',
'pol', 'poly', 'polynomial',
'mon', 'monom', 'monomial',
'con', 'const', 'constant'
'exp', 'exponen', 'exponential'
'polyg' 'polygon', 'polygonal'
RANGEVAL ... an array of length 2 containing the lower and upper
boundaries for the rangeval of argument values
NBASIS ... the number of basis functions
PARAMS ... If the basis is 'fourier', this is a single number indicating
the period. That is, the basis functions are periodic on
the interval (0,PARAMS) or any translation of it.
If the basis is 'bspline', the values are interior points at
which the piecewise polynomials join.
Note that the number of basis functions NBASIS is equal
to the order of the Bspline functions plus the number of
interior knots, that is the length of PARAMS.
This means that NBASIS must be at least 1 larger than the
length of PARAMS.
DROPIND ... A set of indices in 1:NBASIS of basis functions to drop
when basis objects are arguments. Default is []; Note
that argument NBASIS is reduced by the number of indices,
and the derivative matrices in VALUES are also clipped.
QUADVALS .. A NQUAD by 2 matrix. The first column contains quadrature
points to be used in a fixed point quadrature. The second
contains quadrature weights. For example, for Simpson's
rule for NQUAD = 7, the points are equally spaced and the
weights are delta.*[1, 4, 2, 4, 2, 4, 1]/3. DELTA is the
spacing between quadrature points. The default is [].
VALUES ... A cell array, with entries containing the values of
the basis function derivatives starting with 0 and
going up to the highest derivative needed. The values
correspond to quadrature points in QUADVALS and it is
up to the user to decide whether or not to multiply
the derivative values by the square roots of the
quadrature weights so as to make numerical integration
a simple matrix multiplication.
Values are checked against QUADVALS to ensure the correct
number of rows, and against NBASIS to ensure the correct
number of columns.
The default is VALUES{1} = [];
Returns
BASISOBJ ... a basis object with slots
type
rangeval
nbasis
params
dropind
quadvals
values
Slot VALUES contains values of basis functions and derivatives at
quadrature points weighted by square root of quadrature weights.
These values are only generated as required, and only if slot
quadvals is not empty.
An alternative name for this function is CREATE_BASIS, but PARAMS argument
must be supplied.
Specific types of bases may be set up more conveniently using functions
CREATE_BSPLINE_BASIS ... creates a b-spline basis
CREATE_FOURIER_BASIS ... creates a fourier basis
CREATE_POLYGON_BASIS ... creates a polygonal basis
CREATE_MONOM_BASIS ... creates a monomial basis
CREATE_POLYNOMIAL_BASIS ... creates a polynomial basis
CREATE_CONSTANT_BASIS ... creates a constant basis
ans =
1401 13
ans =
0 2.3000
ans =
103
ans =
1401 11
Creating Functional Data Objects
fd objects combine a basis with coefficients
coefs = randn(getnbasis(fourier_obj),1); fd_obj = fd(coefs,fourier_obj); figure(9) plot(fd_obj) % Alternatively fvals = eval_fd(fdatime,fd_obj); figure(10) plot(fdatime,fvals); % We can also define collections of fd objects coefs = randn(getnbasis(fourier_obj),10); fdnames = {'fdatime','reps','height'}; fd_obj = fd(coefs,fourier_obj,fdnames); figure(11) plot(fd_obj) % I can evaluate derivatives by fvals = eval_fd(fdatime,fd_obj,3); size(fvals) figure(12) plot(fdatime,fvals);
ans =
1401 10
Smoothing Data
% We saw in class that smoothing with a basis expansion just requires % solving the normal equations. There is a function to do this % automatically: fdax_b = data2fd(fdarray(:,:,1),fdatime,bspline_obj); figure(13) plot(fdax_b) figure(12) plotfit_fd(fdarray(:,:,1),fdatime,fdax_b) % Now if I want to plot a derivative: figure(14) plot(fdax_b,1) % compare this with figure(15) plot(fdatime(1:1400),diff(fdarray(:,:,1))) % Let's look at a fourrier basis: fdax_f = data2fd(fdarray(:,:,1),fdatime,fourier_obj); plotfit_fd(fdarray(:,:,1),fdatime,fdax_f); % Or a monomial basis fdax_m = data2fd(fdarray(:,:,1),fdatime,monomial_obj); plotfit_fd(fdarray(:,:,1),fdatime,fdax_m); % I can recover various coef_f = getcoef(fdax_f); basis_f = getbasis(fdax_f); basis_f
Basis: Type: fourier Range: 0 to 2.3 Number of basis functions: 11 Period: 2.3
Numerical Quadrature
% For approximating quadrature; usually. % I can set quadrature points and weights. r_b = getbasisrange(basis_f); m = 1000; quadpts = linspace(r_b(1),r_b(2),2*m+1); % Define a fine grid over the range quadwts = ones(2*m+1,1); % Simpson's rule quadrature weights quadwts(2*(1:(m-1))) = 2; quadwts(2*(1:m)+1) = 4; quadwts = quadwts*(r_b(2)-r_b(1))/(6*m); basis_f = putquadvals(basis_f,[quadpts' quadwts]); quads = getquadvals(basis_f); size(quads) % There are times when it may be useful to have quick access to derivatives % of the basis values at the quadrature points. k = 2; % number of derivatives I'm interested in values = cell(3,1); for i = 0:k, values{i+1} = diag(sqrt(quadwts))*eval_basis(quadpts,basis_f,i); end basis_f = putvalues(basis_f,values); vals = getvalues(basis_f,1); size(vals) figure(16) plot(vals) % Note that I multiplied by the square-root of quadwts. This is because I % will usually be interested in estimating the integral of one basis % against another. % % If I have B1 and B2 evaluated at quadrature points, then their inner % product is % % B1'*diag(quadwts)*B2 % % multiplying by the square root of quadwts avoids using the middle term.
ans =
2001 2
ans =
2001 11
Multiple Knots
% In order to demonstrate the use of B-splines for known discontinuities fda = [fdarrayx; fdarrayy]; fdatime2 = [fdatime; max(fdatime)+ fdatime(2) + fdatime]; % Now set up a order-4 B-splines with three knots at 2.3. r_f = [min(fdatime2) max(fdatime2)]; breaks = linspace(r_f(1),r_f(2),201); knots = sort([breaks 2.3*ones(1,2)]); norder = 4; nbasis = length(knots)+norder-2; basis_b = create_bspline_basis(r_f,nbasis,norder,knots); fd_obj = data2fd(fda(:,1),fdatime2,basis_b); fda1 = eval_fd(fdatime2,fd_obj); figure(17) plot(fdatime2,fda1)
