create.exponential.basis <- function (rangeval=c(0,1), nbasis=NULL, ratevec=NULL, dropind=NULL, quadvals=NULL, values=NULL, basisvalues=NULL, names='exp', axes=NULL) { # This function creates an exponential functional data basis # Arguments # RANGEVAL ... An array of length 2 containing the lower and upper # boundaries for the rangeval of argument values # NBASIS ... The number of basis functions. If this conflicts with # the length of RATEVEC, the latter is used. # RATEVEC ... The rate parameters defining exp(ratevec[i]*x) # DROPIND ... A vector of integers specifying the basis functions to # be dropped, if any. # QUADVALS .. A NQUAD by 2 matrix. The firs t 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 # matrix("numeric",0,0). # VALUES ... A list, 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 value of is VALUES is vector("list",0). # 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 matrix("numeric",0,0). # BASISVALUES ... A vector of lists, allocated by code such as # vector("list",1). # This field is designed to avoid evaluation of a # basis system repeatedly at a set of argument values. # Each list within the vector corresponds to a specific set # of argument values, and must have at least two components, # which may be tagged as you wish. # The first component in an element of the list vector contains the # argument values. # The second component in an element of the list vector # contains a matrix of values of the basis functions evaluated # at the arguments in the first component. # The third and subsequent components, if present, contain # matrices of values their derivatives up to a maximum # derivative order. # Whenever function getbasismatrix is called, it checks # the first list in each row to see, first, if the number of # argument values corresponds to the size of the first dimension, # and if this test succeeds, checks that all of the argument # values 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. # For example, you might set up a vector of argument values # called "evalargs" along with a matrix of basis function # values for these argument values called "basismat". # You might want too use tags like "args" and "values", # respectively for these. You would then assign them # to BASISVALUES with code such as # basisobj$basisvalues <- vector("list",1) # basisobj$basisvalues[[1]] <- # list(args=evalargs, values=basismat) # Returns # BASISOBJ ... a functional data basis object of type "expon" # Last modified 9 November 2008 by Spencer Graves # Last modified 6 January 2008 by Jim Ramsay # Default basis for missing arguments ## ## 1. Check RANGEVAL ## if(!is.numeric(rangeval)) stop('rangaval must be numeric; class(rangeval) = ', class(rangeval) ) if(length(rangeval)<1) stop('rangeval must be a numeric vector of length 2; ', 'length(rangeval) = 0.') if (length(rangeval) == 1) { if( rangeval <= 0) stop("rangeval a single value that is not positive: ", rangeval) rangeval <- c(0,rangeval) } if(length(rangeval)>2) stop('rangeval must be a vector of length 2; ', 'length(rangeval) = ', length(rangeval)) if(diff(rangeval)<=0) stop('rangeval must cover a positive range; diff(rangeval) = ', diff(rangeval) ) ## ## 2. check nbasis and ratevec ## { if(is.null(nbasis)){ if(is.null(ratevec)){ nbasis <- 2 ratevec <- 0:1 } else { nbasis <- length(ratevec) if(nbasis<1) stop('ratevec must have positive length; length(ratevec) = 0') if(!is.numeric(ratevec)) stop('ratevec must be numeric; class(ratevec) = ', class(ratevec) ) if(length(unique(ratevec)) != nbasis) stop('ratevec contains duplicates; not allowed.') } } else { if(is.null(ratevec)) ratevec <- 0:(nbasis-1) else{ if(length(ratevec) != nbasis) stop('length(ratevec) must equal nbasis; length(ratevec) = ', length(ratevec), ' != ', 'nbasis = ', nbasis) if(length(unique(ratevec)) != nbasis) stop('ratevec contains duplicates; not allowed.') } } } ## ## 3. check DROPIND ## if (length(dropind) > 0){ if(!is.numeric(dropind)) stop('dropind must be numeric; is ', class(dropind)) doops <- which((dropind%%1)>0) if(length(doops)>0) stop('dropind must be integer; element ', doops[1], " = ", dropind[doops[1]], '; fractional part = ', dropind[doops[1]] %%1) # doops0 <- which(dropind<=0) if(length(doops0)>0) stop('dropind must be positive integers; element ', doops0[1], ' = ', dropind[doops0[1]], ' is not.') doops2 <- which(dropind>nbasis) if(length(doops2)>0) stop("dropind must not exceed nbasis = ", nbasis, '; dropind[', doops2[1], '] = ', dropind[doops2[1]]) # dropind <- sort(dropind) if(length(dropind) > 1) { if(min(diff(dropind)) == 0) stop("Multiple index values in DROPIND.") } } ## ## 4. set up the basis object ## type <- "expon" params <- as.vector(ratevec) basisobj <- basisfd(type=type, rangeval=rangeval, nbasis=nbasis, params=params, dropind=dropind, quadvals=quadvals, values=values, basisvalues=basisvalues) ## ## 5. names ## { if(length(names) == nbasis) basisobj$names <- names else { if(length(names)>1) stop('length(names) = ', length(names), '; must be either ', '1 or nbasis = ', nbasis) basisobj$names <- paste(names, 0:(nbasis-1), sep="") } } ## ## 6. Done ## if(!is.null(axes))basisobj$axes <- axes basisobj }