Skip to content

Commit fcffd4a

Browse files
JamesRamsay5JamesRamsay5
committed
fRegress issues
1 parent 543d17d commit fcffd4a

14 files changed

+1056
-2568
lines changed

‎.Rhistory

+129-129
Original file line numberDiff line numberDiff line change
@@ -1,132 +1,3 @@
1-
Knots <- growth$age
2-
norder <- 6
3-
nbasis <- length(Knots) + norder - 2
4-
hgtbasis <- create.bspline.basis(range(Knots), nbasis, norder, Knots)
5-
Lfdobj <- 4
6-
lambda <- 1e-2
7-
growfdPar <- fdPar(hgtbasis, Lfdobj, lambda)
8-
hgtfd <- smooth.basis(growth$age,
9-
cbind(growth$hgtm,growth$hgtf),growfdPar)$fd
10-
cbasis = create.constant.basis(range(Knots))
11-
maleind = c(rep(1,ncol(growth$hgtm)),rep(0,ncol(growth$hgtf)))
12-
constfd = fd( matrix(1,1,length(maleind)),cbasis)
13-
maleindfd = fd( matrix(maleind,1,length(maleind)),cbasis)
14-
xfdlist = list(constfd,maleindfd)
15-
betalist = list(fdPar(hgtbasis,2,1e-6),fdPar(hgtbasis,2,1e-6))
16-
Fres = Fperm.fd(hgtfd,xfdlist,betalist,nperm=10)
17-
F.res2 = Fperm.fd(annualprec, xfdlist, betalist, nperm=10)
18-
annualprec <- log10(apply(
19-
CanadianWeather$dailyAv[,,"Precipitation.mm"], 2,sum))
20-
# set up a smaller basis using only 40 Fourier basis functions
21-
# to save some computation time
22-
smallnbasis <- 40
23-
smallbasis <- create.fourier.basis(c(0, 365), smallnbasis)
24-
tempfd <- smooth.basis(day.5, CanadianWeather$dailyAv[,,"Temperature.C"],
25-
smallbasis)$fd
26-
constantfd <- fd(matrix(1,1,35), create.constant.basis(c(0, 365)))
27-
xfdlist <- vector("list",2)
28-
xfdlist[[1]] <- constantfd
29-
xfdlist[[2]] <- tempfd[1:35]
30-
betalist <- vector("list",2)
31-
# set up the first regression function as a constant
32-
betabasis1 <- create.constant.basis(c(0, 365))
33-
betafd1 <- fd(0, betabasis1)
34-
betafdPar1 <- fdPar(betafd1)
35-
betalist[[1]] <- betafdPar1
36-
nbetabasis <- 35
37-
betabasis2 <- create.fourier.basis(c(0, 365), nbetabasis)
38-
betafd2 <- fd(matrix(0,nbetabasis,1), betabasis2)
39-
lambda <- 10^12.5
40-
harmaccelLfd365 <- vec2Lfd(c(0,(2*pi/365)^2,0), c(0, 365))
41-
betafdPar2 <- fdPar(betafd2, harmaccelLfd365, lambda)
42-
betalist[[2]] <- betafdPar2
43-
# Should use the default nperm = 200
44-
# but use 10 to save test time for illustration
45-
F.res2 = Fperm.fd(annualprec, xfdlist, betalist, nperm=10)
46-
q()
47-
q()
48-
source('~/Documents/R/fda_revision/fda/R/smooth.basis.R')
49-
library(fda)
50-
b1.1 <- create.bspline.basis(nbasis=1, norder=1)
51-
y12 <- 1:2
52-
fd1.1 <- Data2fd(y12, basisobj=b1.1)
53-
source('~/Documents/R/fda_revision/fda/R/Data2fd.R')
54-
fd1.1 <- Data2fd(y12, basisobj=b1.1)
55-
plot(fd1.1)
56-
fd1.1.5 <- Data2fd(y12, basisobj=b1.1, lambda=0.5)
57-
eval.fd(seq(0, 1, .2), fd1.1.5)
58-
b1.2 <- create.bspline.basis(nbasis=2, norder=1)
59-
# fit the data without smoothing
60-
fd1.2 <- Data2fd(1:2, basisobj=b1.2)
61-
op <- par(mfrow=c(2,1))
62-
plot(b1.2, main='bases')
63-
plot(fd1.2, main='fit')
64-
par(op)
65-
b1.1 <- create.bspline.basis(nbasis=1, norder=1)
66-
# data values: 1 and 2, with a mean of 1.5
67-
y12 <- 1:2
68-
# smooth data, giving a constant function with value 1.5
69-
fd1.1 <- Data2fd(y12, basisobj=b1.1)
70-
plot(fd1.1)
71-
fd1.1.5 <- Data2fd(y12, basisobj=b1.1, lambda=0.5)
72-
eval.fd(seq(0, 1, .2), fd1.1.5)
73-
source('~/Documents/R/fda_revision/fda/R/Data2fd.R')
74-
source('~/Documents/R/fda_revision/fda/R/Data2fd.R')
75-
source('~/Documents/R/fda_revision/fda/R/smooth.basisPar.R')
76-
source('~/Documents/R/fda_revision/fda/R/Data2fd.R')
77-
help("with")
78-
source('~/Documents/R/fda_revision/fda/R/Data2fd.R')
79-
source('~/Documents/R/fda_revision/fda/R/Data2fd.R')
80-
source('~/Documents/R/fda_revision/fda/R/Data2fd.R')
81-
source('~/Documents/R/fda_revision/fda/R/Data2fd.R')
82-
source('~/Documents/R/fda_revision/fda/R/Data2fd.R')
83-
source('~/Documents/R/fda_revision/fda/R/smooth.basisPar.R')
84-
source('~/Documents/R/fda_revision/fda/R/smooth.basisPar.R')
85-
source('~/Documents/R/fda_revision/fda/R/smooth.basis.R')
86-
source('~/Documents/R/fda_revision/fda/R/smooth.basis.R')
87-
source('~/Documents/R/fda_revision/fda/R/smooth.basis.R')
88-
source('~/Documents/R/fda_revision/fda/R/smooth.basisPar.R')
89-
source('~/Documents/R/fda_revision/fda/R/Data2fd.R')
90-
source('~/Documents/R/fda_revision/fda/R/Data2fd.R')
91-
source('~/Documents/R/fda_revision/fda/R/Data2fd.R')
92-
source('~/Documents/R/fda_revision/fda/R/Data2fd.R')
93-
source('~/Documents/R/fda_revision/fda/R/Data2fd.R')
94-
library(fda)
95-
invasion1 <- as.Date('1775-09-04')
96-
invasion2 <- as.Date('1812-07-12')
97-
earlyUS.Canada <- c(invasion1, invasion2)
98-
BspInvasion <- create.bspline.basis(earlyUS.Canada)
99-
earlyYears <- seq(invasion1, invasion2, length.out=7)
100-
(earlyQuad <- (as.numeric(earlyYears-invasion1)/365.24)^2)
101-
fitQuad <- Data2fd(earlyYears, earlyQuad, BspInvasion)
102-
source('~/Documents/R/fda_revision/fda/R/smooth.basis.R')
103-
source('~/Documents/R/fda_revision/fda/R/smooth.basis.R')
104-
source('~/Documents/R/fda_revision/fda/R/monomialpen.R')
105-
source('~/Documents/R/fda_revision/fda/R/monomialpen.R')
106-
source('~/Documents/R/fda_revision/fda/R/monomialpen.R')
107-
source('~/Documents/R/fda_revision/fda/R/basisfd.R')
108-
source('~/Documents/R/fda_revision/fda/R/basisfd.R')
109-
source('~/Documents/R/fda_revision/fda/R/monomialpen.R')
110-
source('~/Documents/R/fda_revision/fda/R/smooth.basis.R')
111-
source('~/Documents/R/fda_revision/fda/R/monomialpen.R')
112-
is.integer(1.2)
113-
source('~/Documents/R/fda_revision/fda/R/basisfd.R')
114-
source('~/Documents/R/fda_revision/fda/R/basisfd.R')
115-
source('~/Documents/R/fda_revision/fda/R/basisfd.R')
116-
source('~/Documents/R/fda_revision/fda/R/basisfd.R')
117-
round(1.2)
118-
source('~/Documents/R/fda_revision/fda/R/basisfd.R')
119-
source('~/Documents/R/fda_revision/fda/R/smooth.basis.R')
120-
source('~/Documents/R/fda_revision/fda/R/smooth.basis.R')
121-
source('~/Documents/R/fda_revision/fda/R/smooth.basis.R')
122-
class(earlyYears)
123-
help("Date")
124-
earlyYears
125-
as.vector(earlyYears)
126-
source('~/Documents/R/fda_revision/fda/R/basisfd.R')
127-
AmRev.ct <- as.POSIXct1970(c('1776-07-04', '1789-04-30'))
128-
BspRev.ct <- create.bspline.basis(AmRev.ct)
129-
AmRevYrs.ct <- seq(AmRev.ct[1], AmRev.ct[2], length.out=14)
1301
(AmRevLin.ct <- as.numeric(AmRevYrs.ct-AmRev.ct[1]))
1312
AmRev.ct <- as.POSIXct1970(c('1776-07-04', '1789-04-30'))
1323
BspRev.ct <- create.bspline.basis(AmRev.ct)
@@ -510,3 +381,132 @@ q()
510381
q()
511382
getwd()
512383
q()
384+
library(fds)
385+
data("Adelaide")
386+
data()
387+
help("SAelectdemand")
388+
Rcpp::sourceCpp('~/Documents/R/Data2LD/src/prodIntfd.cpp')
389+
Rcpp::sourceCpp('~/Documents/R/Data2LD/src/prodIntfd.cpp')
390+
Rcpp::sourceCpp('~/Documents/R/Data2LD/src/prodIntfd.cpp')
391+
Rcpp::sourceCpp('~/Documents/R/Data2LD/src/prodIntfd.cpp')
392+
Rcpp::sourceCpp('~/Documents/R/Data2LD/src/prodIntfd.cpp')
393+
Rcpp::sourceCpp('~/Documents/R/Data2LD/src/prodIntfd.cpp')
394+
Rcpp::sourceCpp('~/Documents/R/Data2LD/src/prodIntfd.cpp')
395+
Rcpp::sourceCpp('~/Documents/R/Data2LD/src/prodIntfd.cpp')
396+
Rcpp::sourceCpp('~/Documents/R/Data2LD/src/prodIntfd.cpp')
397+
Rcpp::sourceCpp('~/Documents/R/Data2LD/src/prodIntfd.cpp')
398+
Rcpp::sourceCpp('~/Documents/R/Data2LD/src/prodIntfd.cpp')
399+
Rcpp::sourceCpp('~/Documents/R/Data2LD/src/prodIntfd.cpp')
400+
Rcpp::sourceCpp('~/Documents/R/Data2LD/src/prodIntfd.cpp')
401+
Rcpp::sourceCpp('~/Documents/R/Data2LD/src/prodIntfd.cpp')
402+
Rcpp::sourceCpp('~/Documents/R/Data2LD/src/prodIntfd.cpp')
403+
Rcpp::sourceCpp('~/Documents/R/Data2LD/src/prodIntfd.cpp')
404+
Rcpp::sourceCpp('~/Documents/R/Data2LD/src/prodIntfd.cpp')
405+
Rcpp::sourceCpp('~/Documents/R/Data2LD/src/prodIntfdCpp.cpp')
406+
Rcpp::sourceCpp('~/Documents/R/Data2LD/src/prodIntfdCpp.cpp')
407+
Rcpp::sourceCpp('~/Documents/R/Data2LD/src/prodIntfdCpp.cpp')
408+
Rcpp::sourceCpp('~/Documents/R/Data2LD/src/prodIntfdCpp.cpp')
409+
q()
410+
q()
411+
annualprec <- log10(apply(CanadianWeather$dailyAv[,,"Precipitation.mm"],
412+
2,sum))
413+
# The simplest 'fRegress' call is singular with more bases
414+
# than observations, so we use a small basis for this example
415+
smallbasis <- create.fourier.basis(c(0, 365), 25)
416+
# There are other ways to handle this,
417+
# but we will not discuss them here
418+
tempfd <- smooth.basis(day.5,
419+
CanadianWeather$dailyAv[,,"Temperature.C"], smallbasis)$fd
420+
getwd()
421+
library(fda)
422+
annualprec <- log10(apply(CanadianWeather$dailyAv[,,"Precipitation.mm"],
423+
2,sum))
424+
tempfd <- smooth.basis(day.5,
425+
CanadianWeather$dailyAv[,,"Temperature.C"], smallbasis)$fd
426+
smallbasis <- create.fourier.basis(c(0, 365), 25)
427+
# There are other ways to handle this,
428+
# but we will not discuss them here
429+
tempfd <- smooth.basis(day.5,
430+
CanadianWeather$dailyAv[,,"Temperature.C"], smallbasis)$fd
431+
precip.Temp.f <- fRegress(annualprec ~ tempfd)
432+
precip.Temp.mdl <- fRegress(annualprec ~ tempfd, method="model")
433+
precip.Temp.m <- do.call('fRegress', precip.Temp.mdl)
434+
nbetabasis <- 21
435+
betabasis2. <- create.fourier.basis(c(0, 365), nbetabasis)
436+
betafd2. <- fd(rep(0, nbetabasis), betabasis2.)
437+
# add smoothing
438+
betafdPar2. <- fdPar(betafd2., lambda=10)
439+
precip.Temp.mdl2 <- precip.Temp.mdl
440+
precip.Temp.mdl2[['betalist']][['tempfd']] <- betafdPar2.
441+
precip.Temp.m2 <- do.call('fRegress', precip.Temp.mdl2)
442+
precip.Temp.f[['df']] # 26
443+
precip.Temp.m2[['df']]# 22 = saved 4 degrees of freedom
444+
(var.e.f <- mean(with(precip.Temp.f, (yhatfdobj-yfdPar)^2)))
445+
(var.e.m2 <- mean(with(precip.Temp.m2, (yhatfdobj-yfdPar)^2)))
446+
xfdlist <- list(const=rep(1, 35), tempfd=tempfd)
447+
betabasis1 <- create.constant.basis(c(0, 365))
448+
betafd1 <- fd(0, betabasis1)
449+
betafdPar1 <- fdPar(betafd1)
450+
betafd2 <- with(tempfd, fd(basisobj=basis, fdnames=fdnames))
451+
# convert to an fdPar object
452+
betafdPar2 <- fdPar(betafd2)
453+
betalist <- list(const=betafdPar1, tempfd=betafdPar2)
454+
precip.Temp <- fRegress(annualprec, xfdlist, betalist)
455+
TempRgn.mdl <- fRegress(tempfd ~ region, CanadianWeather, method='m')
456+
%names(TempRgn.mdl)
457+
# make desired modifications here
458+
# then run
459+
TempRgn.m <- do.call('fRegress', TempRgn.mdl)
460+
TempRgn.mdl <- fRegress(tempfd ~ region, CanadianWeather, method='m')
461+
names(TempRgn.mdl)
462+
TempRgn.m <- do.call('fRegress', TempRgn.mdl)
463+
class(TempRgn.mdl)
464+
class(TempRgn.mdl.y)
465+
class(TempRgn.mdl$y)
466+
source('~/Documents/R/fda_work/fda/R/fRegress.R')
467+
TempRgn.m <- do.call('fRegress', TempRgn.mdl)
468+
source('~/Documents/R/fda_work/fda/R/fRegress.R')
469+
TempRgn.m <- do.call('fRegress', TempRgn.mdl)
470+
source('~/Documents/R/fda_work/fda/R/fRegress.R')
471+
TempRgn.m <- do.call('fRegress', TempRgn.mdl)
472+
rm("fRegress.fdPar.R")
473+
source('~/Documents/R/fda_work/fda/R/fRegress.R')
474+
source('~/Documents/R/fda_work/fda/R/fRegress.R')
475+
TempRgn.m <- do.call('fRegress', TempRgn.mdl)
476+
library(fda)
477+
q()
478+
source('~/Documents/R/fda_work/fda/R/fRegress.R')
479+
source('~/Documents/R/fda_work/fda/R/fRegress.R')
480+
q()
481+
q()
482+
library(fda)
483+
print(names(CanadianWeather))
484+
# set up log10 of annual precipitation for 35 weather stations
485+
annualprec <- log10(apply(CanadianWeather$dailyAv[,,"Precipitation.mm"], 2,sum))
486+
# The simplest 'fRegress' call is singular with more bases
487+
# than observations, so we use only 25 basis functions, for this example
488+
smallbasis <- create.fourier.basis(c(0, 365), 25)
489+
# The covariate is the temperature curve for each station.
490+
tempfd <- smooth.basis(day.5,
491+
CanadianWeather$dailyAv[,,"Temperature.C"], smallbasis)$fd
492+
TempRgn.mdl <- fRegress(tempfd ~ region, CanadianWeather, method='f')
493+
source('~/Documents/R/fda_work/fda/R/fRegress.R')
494+
TempRgn.mdl <- fRegress(tempfd ~ region, CanadianWeather, method='f')
495+
TempRgn.mdl <- fRegress(tempfd ~ region, CanadianWeather, method='f')
496+
TempRgn.mdl <- fRegress(tempfd ~ region, CanadianWeather, method='m')
497+
TempRgn.mdl <- fRegress(tempfd ~ region, CanadianWeather, method='f')
498+
library(fda)
499+
source('~/Documents/R/fda_work/fda/R/fRegress.R')
500+
# data are in Canadian Weather object
501+
# print the names of the data
502+
print(names(CanadianWeather))
503+
# set up log10 of annual precipitation for 35 weather stations
504+
annualprec <- log10(apply(CanadianWeather$dailyAv[,,"Precipitation.mm"], 2,sum))
505+
# The simplest 'fRegress' call is singular with more bases
506+
# than observations, so we use only 25 basis functions, for this example
507+
smallbasis <- create.fourier.basis(c(0, 365), 25)
508+
# The covariate is the temperature curve for each station.
509+
tempfd <- smooth.basis(day.5,
510+
CanadianWeather$dailyAv[,,"Temperature.C"], smallbasis)$fd
511+
TempRgn.mdl <- fRegress(tempfd ~ region, CanadianWeather, method='f')
512+
q()

‎DESCRIPTION

+3-3
Original file line numberDiff line numberDiff line change
@@ -1,12 +1,12 @@
11
Package: fda
2-
Version: 5.1.6
3-
Date: 2020-08-7
2+
Version: 5.1.7
3+
Date: 2020-09-04
44
Title: Functional Data Analysis
55
Author: J. O. Ramsay <ramsay@psych.mcgill.ca> [aut,cre],
66
Spencer Graves <spencer.graves@effectivedefense.org> [ctb],
77
Giles Hooker <gjh27@cornell.edu> [ctb]
88
Maintainer: J. O. Ramsay <ramsay@psych.mcgill.ca>
9-
Depends: R (>= 3.5), Matrix
9+
Depends: R (>= 3.5), Matrix, fds
1010
Suggests: deSolve, lattice
1111
Description: These functions were developed to support functional data
1212
analysis as described in Ramsay, J. O. and Silverman, B. W.

‎NAMESPACE

+12-5
Original file line numberDiff line numberDiff line change
@@ -40,6 +40,7 @@ export(AmpPhaseDecomp,
4040
deriv.fd,
4141
df2lambda,
4242
dirs,
43+
eigchk,
4344
Eigen,
4445
eigen.pda,
4546
eval.basis,
@@ -60,8 +61,13 @@ export(AmpPhaseDecomp,
6061
fourierpen,
6162
Fperm.fd,
6263
fRegress,
64+
fRegress.fd,
65+
fRegress.fdPar,
6366
fRegress.CV,
6467
fRegress.stderr,
68+
fRegressArgCheck,
69+
fRegress.formula,
70+
fRegress.character,
6571
Fstat.fd,
6672
geigen,
6773
getbasismatrix,
@@ -156,6 +162,7 @@ import("graphics")
156162
import("stats")
157163
import("splines")
158164
import("Matrix")
165+
import("fds")
159166

160167
importFrom("grDevices", "dev.new", "heat.colors")
161168
importFrom("utils", "data")
@@ -184,14 +191,12 @@ S3method(predict, fd)
184191
S3method(sqrt, fd)
185192
S3method(sum, fd)
186193
S3method(summary, fd)
187-
S3method(fRegress, fd)
188194
S3method(boxplot, fdPar)
189195
S3method(coef, fdPar)
190196
S3method(coefficients, fdPar)
191197
S3method(plot, fdPar)
192198
S3method(predict, fdPar)
193199
S3method(summary, fdPar)
194-
S3method(fRegress, fdPar)
195200
S3method(as.fd, fdSmooth)
196201
S3method(coef, fdSmooth)
197202
S3method(boxplot, fdSmooth)
@@ -209,12 +214,14 @@ S3method(residuals, monfd)
209214
S3method(predict, monfd)
210215
S3method(as.fd, 'function')
211216
S3method(as.fd, smooth.spline)
217+
S3method(fRegress, fd)
218+
S3method(fRegress, fdPar)
212219
S3method(fRegress, numeric)
220+
S3method(fRegress, double)
221+
S3method(fRegress, formula)
222+
S3method(fRegress, character)
213223
S3method(fRegress, stderr)
214224
S3method(fRegress, CV)
215-
S3method(fRegress, character)
216-
S3method(fRegress, formula)
217-
S3method(matplot, numeric)
218225
S3method(fitted, posfd)
219226
S3method(predict, posfd)
220227
S3method(residuals, posfd)

‎R/eigchk.R

+14
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,18 @@
11
eigchk <- function(Cmat) {
2+
3+
# Last modified 25 August 2020 by Jim Ramsay
4+
5+
# Cmat for NA's
6+
7+
if (any(is.na(Cmat))) stop("Cmat has NA values.")
8+
9+
# check Cmat for Cmatmetry
10+
11+
if (max(abs(Cmat-t(Cmat)))/max(abs(Cmat)) > 1e-10) {
12+
stop('CMAT is not symmetric.')
13+
} else {
14+
Cmat <- (Cmat + t(Cmat))/2
15+
}
216

317
# check Cmat for singularity
418

0 commit comments

Comments
 (0)