https://github.com/cran/fda
Tip revision: aac3e3f7207f4732067b6e1d1b735f5d7855664f authored by J. O. Ramsay on 06 July 2007, 00:00:00 UTC
version 1.2.4
version 1.2.4
Tip revision: aac3e3f
quadset.R
quadset <- function(nquad=5, basisobj=NULL, breaks){
#function [basisobj, quadpts, quadwts] = quadset(nquad, basisobj, breaks)
# last modified 2007 April 25 by Spencer Graves
#% Matlab version last modified 14 August 2006
##
## 1. Check nquad
##
{
if(nquad<5){
warning("nquad must be at least 5; increase to this minimum.")
nquad <- 5
}
else {
if((nquad%%2)!=1){
warning("nquad must be an odd integer; increased to enforce this.")
nquad <- 1+2*floor(nquad/2)
}
}
}
##
## 2. check basisobj
##
if(!is.null(basisobj) && !is.basis(basisobj))
stop('basisobj is not a basis object.')
##
## 3. check breaks
##
if(missing(breaks)) {
if(is.null(basisobj) || !is.basis(basisobj))
stop("Either 'breaks' or 'basisobj' must be provided.")
#
# type = getbasistype(basisobj);
type <- basisobj$type
if(type != 'bspline')
stop(
"'breaks' not supplied and 'basisobj' is not a spline basis.")
#
# rangeval = getbasisrange(basisobj);
rangeval = basisobj$range #?????
# params = getbasispar(basisobj);
params = basisobj$params
# knots = [rangeval(1), params, rangeval(2)];
knots <- c(rangeval[1], params, rangeval[2])
breaks = unique(knots)
}
##
## 4. quadpts and quadwts
##
nbreaks = length(breaks);
#
db <- diff(breaks)
nquad1 <- nquad-1
nbreaks1 <- nbreaks-1
# 4.1. First create quadpts as a matrix
quadpts. <- array(NA, dim=c(nbreaks1, nquad) )
quadpts.[, 1] <- breaks[-nbreaks]
db. <- db/nquad1
for(i in 2:nquad)
quadpts.[, i] <- (quadpts.[, i-1]+db.)
# 4.2. Now convert quadpts matrix to the desired vector
quadpts <- as.vector(t(quadpts.))
# 4.3. Similarly first create quadwts as a matrix
quadwts. <- outer(c(1, 4, rep(c(2, 4), (nquad1-2)/2), 1),
db/(nquad1*3) )
# 4.4. Now convert quadwts matrix to the desired vector
quadwts <- as.vector(quadwts.)
# quadpts1 = linspace(breaks(1),breaks(2),nquad)';
# quadwts1 = ones(nquad,1);
# quadwts1(2:2:nquad-1) = 4;
# quadwts1(3:2:nquad-2) = 2;
# quadwts1 = ((breaks(2)-breaks(1))/(nquad-1)).*quadwts1/3;
# quadvals = [quadpts1, quadwts1];
# for( i in 3:nbreaks )
# quadptsi = linspace(breaks(i-1),breaks(i),nquad)';
# quadwtsi = ones(nquad,1);
# quadwtsi(3:2:nquad-2) = 2;
# quadwtsi(2:2:nquad-1) = 4;
# quadwtsi = ((breaks(i)-breaks(i-1))/(nquad-1)).*quadwtsi/3;
# quadvals = [quadvals;[quadptsi, quadwtsi]];
#end
quadvals <- cbind(quadpts=quadpts, quadwts=quadwts)
#quadpts = quadvals(:,1);
#quadwts = quadvals(:,2);
#basisobj = putquadvals(basisobj, quadvals);
if(is.null(basisobj))return(quadvals)
#
basisobj$quadvals <- quadvals
values <- vector('list', 2)
for( ivalue in 1:2){
# basisvalues = eval_basis(quadpts, basisobj, ivalue-1);
# values{ivalue} = basisvalues;
values[[ivalue]] <- eval.basis(quadpts, basisobj, ivalue-1)
}
# basisobj = putvalues(basisobj, values);
basisobj$values <- values
basisobj
}