Commit beb99e56 authored by KaKiLa's avatar KaKiLa
Browse files

Initial commit

To cite GNU Octave in publications use:
Ramsay, J., Hooker, G., & Graves, S. (2009).
Functional Data Analysis with R and MATLAB.
New York, NY: Springer New York.
DOI: 10.1007/978-0-387-98185-7
A BibTeX entry for LaTeX users is:
address = {New York, NY},
author = {Ramsay, James and Hooker, Giles and Graves, Spencer},
doi = {10.1007/978-0-387-98185-7},
isbn = {978-0-387-98184-0},
publisher = {Springer New York},
title = {{Functional Data Analysis with R and MATLAB}},
url = {},
year = {2009}
This diff is collapsed.
Name: fda
Version: 1.0.0
Date: 23-02-2018
Author: J. O. Ramsay <ramsay at>, Hadley Wickham, Spencer Graves, Giles Hooker
Maintainer: Juan Pablo Carbajal <>
Title: Functional Data Analysis
Description: These functions were developed to support functional data analysis as described in Ramsay, J. O. and Silverman, B. W. (2005) Functional Data Analysis. New York: Springer.
Depends: octave (>= 4.2.0)
Autoload: no
License: GPLv3+,
## Makefile to simplify Octave Forge package maintenance tasks
PACKAGE = $(shell $(SED) -n -e 's/^Name: *\(\w\+\)/\1/p' DESCRIPTION | $(TOLOWER))
VERSION = $(shell $(SED) -n -e 's/^Version: *\(\w\+\)/\1/p' DESCRIPTION | $(TOLOWER))
DEPENDS = $(shell $(SED) -n -e 's/^Depends[^,]*, \(.*\)/\1/p' DESCRIPTION | $(SED) 's/ *([^()]*),*/ /g')
HTML_TARBALL = $(PACKAGE)-html.tar.gz
MD5SUM ?= md5sum
SED ?= sed
TAR ?= tar
# Follow jwe suggestion on not inheriting these vars from
# the enviroment, so they can be set as command line arguemnts
MKOCTFILE := mkoctfile
OCTAVE := octave --no-gui
TOLOWER = $(SED) -e 'y/ABCDEFGHIJKLMNOPQRSTUVWXYZ/abcdefghijklmnopqrstuvwxyz/'
FIX_PERMISSIONS ?= chmod -R a+rX,u+w,go-w,ug-s
.PHONY: help dist html release install all check run doc clean maintainer-clean get-upstream
@echo "Targets:"
@echo " dist - Create $(RELEASE_TARBALL) for release"
@echo " html - Create $(HTML_TARBALL) for release"
@echo " release - Create both of the above and show md5sums"
@echo " install - Install the package in GNU Octave"
@echo " all - Build all oct files"
@echo " check - Execute package tests (w/o install)"
@echo " run - Run Octave with development in PATH (no install)"
@echo " doc - Build Texinfo package manual"
@echo " clean - Remove releases, html documentation, and oct files"
@echo " maintainer-clean - Additionally remove all generated files"
$(RELEASE_DIR): .hg/dirstate
@echo "Creating package version $(VERSION) release ..."
-rm -rf $@
hg archive --exclude ".hg*" --exclude Makefile --type files $@
$(TAR) cf - --posix $< | gzip -9n > $@
-rm -rf $<
$(HTML_DIR): install
@echo "Generating HTML documentation. This may take a while ..."
-rm -rf $@
$(OCTAVE) --silent \
--eval 'graphics_toolkit ("gnuplot");' \
--eval 'pkg load generate_html $(PACKAGE);' \
--eval 'generate_package_html ("$(PACKAGE)", "$@", "octave-forge");'
$(TAR) cf - --posix $< | gzip -9n > $@
-rm -rf $<
release: dist html
@echo "Upload @"
@echo " and note the changeset the release corresponds to"
@echo "Installing package locally ..."
$(OCTAVE) --silent --eval 'pkg install $(RELEASE_TARBALL);'
check: all
$(OCTAVE) --silent \
--eval 'if(!isempty("$(DEPENDS)")); pkg load $(DEPENDS); endif;' \
--eval 'addpath (fullfile ([pwd filesep "inst"]));' \
--eval 'runtests ("inst");'
run: all
$(OCTAVE) --silent --persist \
--eval 'if(!isempty("$(DEPENDS)")); pkg load $(DEPENDS); endif;' \
--eval 'addpath (fullfile ([pwd filesep "inst"]));'
maintainer-clean: clean
get-upstream: .hg/distate
mkdir -p upstream
unzip -d upstream
rm -f
Summary of important user-visible changes for releases of the fda package
fda-1.0.0 Release Date: xx-xx-xxxx Release Manager: Juan Pablo Carbajal
** First official release.
function Lfdobj = Lfd(nderiv, bwtcell)
% LFD creates a linear differential operator object of the form
% Lx(t) = w_0(t) x(t) + ... + w_{m-1}(t) D^{m-1}x(t) +
% \exp[w_m(t) D^m x(t) where nderiv = nderiv.
% Function x(t) is operated on by this operator L, and the operator
% computes a linear combination of the function and its first nderiv
% derivatives. The function x(t) must be scalar.
% The linear combination of derivatives is defined by the weight
% or coefficient functions w_j(t), and these are assumed to vary
% over t, although of course they may also be constant as a
% special case.
% The weight coefficient for D^m is special in that it must
% be positive to properly identify the operator. This is why
% it is exponentiated. In most situations, it will be 0,
% implying a weight of one, and this is the default.
% The inner products of the linear differential operator L
% applied to basis functions is evaluated in the functions
% called in function EVAL_PENALTY().
% Some important functions also have the capability of allowing
% the argument that is an LFD object be an integer. They convert
% the integer internally to an LFD object by INT2LFD(). These are:
% Arguments:
% NDERIV ... the order of the operator
% the highest order of derivative.
% BWTCELL ... A cell vector object with either NDERIV or
% NDERIV+1 cells.
% If there are NDERIV cells, then the coefficient of
% D^NDERIV is set to 1; otherwise, cell NDERIV+1
% contains a function that is exponentiated to define
% the actual coefficient.
% Simple cases:
% All this generality may not be needed, and, for example,
% often the linear differential operator will be
% simply L = D^m, defining Lx(t) = D^mx(t). Or the weights and
% forcing functions may all have the same bases, in which case
% it is simpler to use a functional data objects to define them.
% These situations cannot be accommodated within Lfd(), but
% there is function int2Lfd(m) that converts a nonnegative
% integer nderiv into an Lfd object equivalent to D^m.
% There is also fd2cell(fdobj) and that converts a functional
% data object into cell object, which can then be used as
% an argument of Lfd().
% Returns:
% LFDOBJ ... a functional data object
% last modified 3 January 2008
% check nderiv
if ~isnumeric(nderiv)
error('Order of operator is not numeric.');
if nderiv ~= round(nderiv)
error('Order of operator is not an integer.');
if nderiv < 0
error('Order of operator is negative.');
% check that BWTCELL is a cell object
if ~iscell(bwtcell)
error('BWTCELL not a cell object.');
if isempty(bwtcell)
if nderiv > 0
error(['Positive derivative order accompanied by ', ...
'empty weight cell.']);
bwtsize = size(bwtcell);
bfdPar = bwtcell{1};
bfd = getfd(bfdPar);
brange = getbasisrange(getbasis(bfd));
% BWTCELL two-dimensional.
% Only possibilities are (1) N > 1, nderiv == 1, and
% (2) N = 1, nderiv >= 1;
if length(bwtsize) == 2
if bwtsize(1) > 1 && bwtsize(2) > 1
error('BWTCELL is not a vector.');
if nderiv > 0
if bwtsize(1) ~= nderiv && bwtsize(2) ~= nderiv
error(['Dimension of BWTCELL not ', ...
'compatible with NDERIV.']);
% BWTCELL has more than two dimensions.
if length(bwtsize) > 2
error('BWTCELL has more than two dimensions.');
% Check the ranges for compatibility
for j=2:nderiv
brangej = getbasisrange(getbasis(getfd(bwtcell{j})));
if any(brangej ~= brange)
error('Incompatible ranges in weight functions.');
% set up the Lfd object.
Lfdobj.nderiv = nderiv;
Lfdobj.bwtcell = bwtcell;
Lfdobj = class(Lfdobj, 'Lfd');
function display(Lfd)
% Last modified 20 July 2006
nderiv = Lfd.nderiv;
fprintf(['NDERIV = ', num2str(nderiv),'\n']);
if nderiv > 0
for ideriv=1:nderiv
fprintf(['\n\nWFD(',num2str(ideriv-1),') fdPar object:\n'])
function bwtcell = getbwtcell(Lfdobj)
% GETWFDCELL Extracts the weight function cell object from LFDOBJ.
% last modified 18 June 2013
if ~isa_Lfd(Lfdobj)
error('Argument is not a linear differential operator object');
bwtcell = Lfdobj.bwtcell;
function nderiv = getnderiv(Lfdobj)
% GETNDERIV Extracts the order of the operator from LFDOBJ.
% last modified 3 January 2008
if ~isa_Lfd(Lfdobj)
error('Argument is not a linear differential operator object');
nderiv = Lfdobj.nderiv;
function wfdcell = getwfdcell(Lfdobj)
% GETWFDCELL Extracts the weight function cell object from LFDOBJ.
% last modified 18 November 2003
if ~isa_Lfd(Lfdobj)
error('Argument is not a linear differential operator object');
wfdcell = Lfdobj.bwtcell;
function intwrd = isinteger(Lfdobj)
% ISINTEGER returns 1 of WFD and AFD are both zero functions.
% Last modified 3 January 2008
% check WFDCELL for emptyness or all zero
wfdcell = getwfdcell(Lfdobj);
wintwrd = 1;
if ~isempty(wfdcell)
nderiv = Lfdobj.nderiv;
for j=1:nderiv
wfdParj = wfdcell{j};
wfdj = getfd(wfdParj);
if any(getcoef(wfdj) ~= 0.0)
wintwrd = 0;
intwrd = wintwrd;
\ No newline at end of file
function plot(Lfdobj)
% Plots a linear differential operator object.
% The coefficients \beta_j defining the homogeneous
% part of the operator are plotted first, followed
% by the coefficients for the forcing functions.
m = Lfdobj.nderiv;
wfdcell = Lfdobj.bwtcell;
nplot = m;
for j=1:m
function basisobj = basis(basistype, rangeval, nbasis, params, ...
dropind, quadvals, values, basisvalues)
% BASIS Creates a functional data basis.
% Arguments or slots:
% BASISTYPE ... a string indicating the type of basis. This may be one of
% 'Bspline', 'bspline', 'Bsp', 'bsp',
% 'con', 'const', 'constant'
% 'exp', 'exponen', 'exponential'
% 'Fourier', 'fourier', 'Fou', 'fou',
% 'mon', 'monom', 'monomial',
% 'Nspline', 'nspline', 'Nsp', 'nsp',
% 'polyg' 'polygon', 'polygonal'
% 'pol', 'poly', 'polynomial'
% 'power', 'pow'
% 'QW'
% 'QWM'
% 'QS'
% 'slide'
% 'fd'
% 'FEM'
% 'TP'
% 'fdVariance', 'fdV', fdVar'
% 'IRT3PL'
% 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 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 cells 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. It is
% up to the user 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} = [];
% BASISVALUES ... A cell array. The cell array must be 2-dimensional,
% with a variable number of rows and two or more columns.
% This field is designed to avoid evaluation of a
% basis system repeatedly at a set of argument values.
% Each row corresponds to a specific set of argument values.
% The first cell in that row contains the argument values.
% The second cell in that row contains a matrix of values of
% the basis functions.
% The third and subsequent cells contain matrices of values
% their derivatives up to a maximum derivative order.
% Whenever function getbasismatrix is called, it checks
% the first cell in each row to see, first, if the number of
% argument values corresponds to the size of the first
% imension, and if this test succeeds, checks that all of
% the arguments match. This takes time, of course, but is
% much faster than re-evaluation of the basis system. Even
% this time can be avoided by direct retrieval of the
% desired array.
% Returns
% BASISOBJ ... a basis object with slots
% type
% rangeval
% nbasis
% params
% dropind
% quadvals
% values
% basisvalues
% 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_CONSTANT_BASIS ... creates a constant basis
% CREATE_FOURIER_BASIS ... creates a fourier basis
% CREATE_MONOM_BASIS ... creates a monomial basis
% CREATE_NSPLINE_BASIS ... creates a n-spline basis
% CREATE_POLYGON_BASIS ... creates a polygonal basis
% CREATE_POLYNOMIAL_BASIS ... creates a polynomial basis
% CREATE_POWER_BASIS ... creates a polygonal basis
% CREATE_QW_BASIS ... creates a Weibull W basis
% CREATE_QWM_BASIS ... creates a modified Weibull W basis
% CREATE_SLIDE_BASIS ... creates a slide basis
% CREATE_FD_BASIS ... creates a functional data object basis
% CREATE_FEM_BASIS ... creates a finite element basis
% CREATE_TP_BASIS ... creates a tensor product basis
% CREATE_fdVariance_BASIS ... creates a fdVariance basis
% CREATE_IRT3PL_BASIS ... creates a IRT3PL basis
% Steps for setting up new bases:
% 1. Code a "create_..._basis" function, which specifies the basis type,
% parameter structure and so forth. create_bspline_basis is a good
% example of such a function.
% 2. Set up an evaluation function, like bsplineM, for evaluating a
% a basis or one of its derivatives at a set of argument
% value.
% 3. Set up a function for evaluating the penalty matrix by computing
% the inner products of basis functions. If these are not available
% analytically, function inprod can be used if step 2 has been done.
% This step is only required if a roughness penalty is ever needed,
% since very simple or low-dimensional basis probably won't such
% things.
% 4. Add the basis to function getbasismatrix in folder @basis
% 5. Add the basis to function eval_penalty if a roughness penalty is
% ever needed.
% 6. Add the basis to function use_proper_basis inside function basis
% 7 Add the basis to function getbasistype in folder @basis as well
% Last modified 7 July 2014
% Set up default basis if there are no arguments
if nargin==0
basisobj.type = 'bspline';
basisobj.rangeval = [0,1];
basisobj.nbasis = 2;
basisobj.params = [];
basisobj.dropind = [];
basisobj.quadvals = [];
basisobj.values = {};
basisobj.basisvalues = {};
basisobj = class(basisobj, 'basis');
% If arguments are supplied, at least the first four must be supplied.
if nargin < 4
error('Less than four arguments found.');
% if first argument is a basis object, return
if isa(basistype, 'basis')
basisobj = basistype;
% check basistype
basistype = use_proper_basis(basistype);
if strcmp(basistype,'unknown')
error ('TYPE unrecognizable.');
% check if QUADVALS is present, and set to default if not
if nargin < 6
quadvals = [];
if ~isempty(quadvals)
[nquad,ncol] = size(quadvals);
if nquad == 2 && ncol > 2
quadvals = quadvals';
[nquad,ncol] = size(quadvals);
if nquad < 2
error('Less than two quadrature points are supplied.');
if ncol ~= 2
error('QUADVALS does not have two columns.');
% check VALUES if present, and set to a single empty cell
% if not.
if nargin < 7
values = {};
if ~iscell(values)
error('VALUES argument is not a cell array.');
if ~isempty(values)
if ~isempty(values{1})
nvalues = length(values);
for ivalue=1:nvalues
[n,k] = size(full(values{ivalue}));
if n ~= nquad
error(['Number of rows in VALUES not equal to ', ...
'number of quadrature points.']);
if k ~= nbasis
error(['Number of columns in VALUES not equal to ', ...
'number of basis functions.']);
% check BASISVALUES are present, and if not set to empty cell array
if nargin < 8
basisvalues = {};
if ~isempty(basisvalues)
if ~iscell(basisvalues)
error('BASISVALUES is not a cell object.')
sizevec = size(basisvalues);
if length(sizevec) > 2
error('BASISVALUES is not 2-dimensional.')
for i=1:sizevec(1)
if length(basisvalues{i,1}) ~= size(basisvalues{i,2},1)
error(['Number of argument values not equal number ', ...
'of values.']);
% check if DROPIND is present, and set to default if not
if nargin < 5
dropind = [];
if ~isempty(dropind)
% check DROPIND
ndrop = length(dropind);
if ndrop >= nbasis
error('Too many index values in DROPIND.');
dropind = sort(dropind);
if ndrop > 1
if any(diff(dropind)) == 0
error('Multiple index values in DROPIND.');