Bases funcionales
#install.packages("fda")
library(fda)
## Warning: package 'fda' was built under R version 4.4.3
## Cargando paquete requerido: splines
## Cargando paquete requerido: fds
## Cargando paquete requerido: rainbow
## Cargando paquete requerido: MASS
## Warning: package 'MASS' was built under R version 4.4.3
## Cargando paquete requerido: pcaPP
## Cargando paquete requerido: RCurl
## Cargando paquete requerido: deSolve
## Warning: package 'deSolve' was built under R version 4.4.3
##
## Adjuntando el paquete: 'fda'
## The following object is masked from 'package:graphics':
##
## matplot
tobs = seq(0,1,0.01)
nobs = length(tobs)
knots = c(seq(0,1,0.1));
nknots = length(knots);
norder = 4;
nbasis = length(knots) + norder - 2;
basis = create.bspline.basis(c(min(tobs),max(tobs)),nbasis,norder,knots);
basis
## $call
## basisfd(type = type, rangeval = rangeval, nbasis = nbasis, params = params,
## dropind = dropind, quadvals = quadvals, values = values,
## basisvalues = basisvalues)
##
## $type
## [1] "bspline"
##
## $rangeval
## [1] 0 1
##
## $nbasis
## [1] 13
##
## $params
## [1] 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
##
## $dropind
## NULL
##
## $quadvals
## NULL
##
## $values
## list()
##
## $basisvalues
## list()
##
## $names
## [1] "bspl4.1" "bspl4.2" "bspl4.3" "bspl4.4" "bspl4.5" "bspl4.6"
## [7] "bspl4.7" "bspl4.8" "bspl4.9" "bspl4.10" "bspl4.11" "bspl4.12"
## [13] "bspl4.13"
##
## attr(,"class")
## [1] "basisfd"
# basis values at samplincurv points
basismat = eval.basis(tobs, basis);
dim(basismat)
## [1] 101 13
##quartz()
plot(tobs,basismat[,1],type = "l",col=1,lwd=3)
lines(tobs,basismat[,2],type = "l",col=2,lwd=3)
##quartz()
matplot(tobs,basismat,type='l',lwd=2,lty=1, xlab='day',ylab='basis',cex.lab=1.5,cex.axis=1.5)
for (i in 1:nknots)
{
abline(v=knots[i],type="l", lty=2, lwd=3)
}
## Warning in int_abline(a = a, b = b, h = h, v = v, untf = untf, ...): el
## parámetro gráfico "type" es obsoleto
## Warning in int_abline(a = a, b = b, h = h, v = v, untf = untf, ...): el
## parámetro gráfico "type" es obsoleto
## Warning in int_abline(a = a, b = b, h = h, v = v, untf = untf, ...): el
## parámetro gráfico "type" es obsoleto
## Warning in int_abline(a = a, b = b, h = h, v = v, untf = untf, ...): el
## parámetro gráfico "type" es obsoleto
## Warning in int_abline(a = a, b = b, h = h, v = v, untf = untf, ...): el
## parámetro gráfico "type" es obsoleto
## Warning in int_abline(a = a, b = b, h = h, v = v, untf = untf, ...): el
## parámetro gráfico "type" es obsoleto
## Warning in int_abline(a = a, b = b, h = h, v = v, untf = untf, ...): el
## parámetro gráfico "type" es obsoleto
## Warning in int_abline(a = a, b = b, h = h, v = v, untf = untf, ...): el
## parámetro gráfico "type" es obsoleto
## Warning in int_abline(a = a, b = b, h = h, v = v, untf = untf, ...): el
## parámetro gráfico "type" es obsoleto
## Warning in int_abline(a = a, b = b, h = h, v = v, untf = untf, ...): el
## parámetro gráfico "type" es obsoleto
## Warning in int_abline(a = a, b = b, h = h, v = v, untf = untf, ...): el
## parámetro gráfico "type" es obsoleto
# Comments: If x(0) = 1, set the coefficient to the first basis function = 1;
# evaluate the first derivative of the basis functions
Dbasismat = eval.basis(tobs, basis,1);
#quartz()
matplot(tobs,Dbasismat,type='l',lwd=2,lty=1, xlab='day',ylab='basis',cex.lab=1.5,cex.axis=1.5)
for (i in 1:nknots)
{
abline(v=knots[i],type="l", lty=2, lwd=3)
}
## Warning in int_abline(a = a, b = b, h = h, v = v, untf = untf, ...): el
## parámetro gráfico "type" es obsoleto
## Warning in int_abline(a = a, b = b, h = h, v = v, untf = untf, ...): el
## parámetro gráfico "type" es obsoleto
## Warning in int_abline(a = a, b = b, h = h, v = v, untf = untf, ...): el
## parámetro gráfico "type" es obsoleto
## Warning in int_abline(a = a, b = b, h = h, v = v, untf = untf, ...): el
## parámetro gráfico "type" es obsoleto
## Warning in int_abline(a = a, b = b, h = h, v = v, untf = untf, ...): el
## parámetro gráfico "type" es obsoleto
## Warning in int_abline(a = a, b = b, h = h, v = v, untf = untf, ...): el
## parámetro gráfico "type" es obsoleto
## Warning in int_abline(a = a, b = b, h = h, v = v, untf = untf, ...): el
## parámetro gráfico "type" es obsoleto
## Warning in int_abline(a = a, b = b, h = h, v = v, untf = untf, ...): el
## parámetro gráfico "type" es obsoleto
## Warning in int_abline(a = a, b = b, h = h, v = v, untf = untf, ...): el
## parámetro gráfico "type" es obsoleto
## Warning in int_abline(a = a, b = b, h = h, v = v, untf = untf, ...): el
## parámetro gráfico "type" es obsoleto
## Warning in int_abline(a = a, b = b, h = h, v = v, untf = untf, ...): el
## parámetro gráfico "type" es obsoleto
# true curve
ytru = (tobs-0.3)^2
plot(tobs,ytru,type = "l")
# put noise to the true curve and generate noisy data
nobs = length(tobs)
noise = 0.03*rnorm(nobs)
yobs = ytru + noise
points(tobs,yobs)
# estimate basis coefficient
Mmat = ginv(t(basismat)%*%basismat)%*%t(basismat)
chat = Mmat%*%yobs
# fitted curve
yhat = basismat%*%chat;
lines(tobs,yhat,type = "l",col="red")
# estimate the variance of noise
SSE = t(yhat-yobs)%*%(yhat-yobs)
sigma2 = SSE/(nobs-nbasis)
sigma2
## [,1]
## [1,] 0.0007502305
sqrt(sigma2)
## [,1]
## [1,] 0.02739034
# estimate the variance of the fitted curve
Smat = basismat%*%Mmat
varYhat = diag(Smat%*%Smat*matrix(sigma2,nobs,nobs))
# 95% confidence interval
yhat025 = yhat-1.96*sqrt(varYhat)
yhat975 = yhat+1.96*sqrt(varYhat)
lines(tobs,yhat025,type="l", lty=2, lwd=3,col="blue")
lines(tobs,yhat975,type="l", lty=2, lwd=3,col="blue")
###########################
# Smoothing Splines
#########################
# Use quadrature to get integral - Composite Simpson's Rule
delta <- 0.02
quadpts <- seq(0,1,delta)
nquadpts <- length(quadpts)
quadwts <- as.vector(c(1,rep(c(4,2),(nquadpts-2)/2),4,1),mode="any")
quadwts <- c(1,rep(c(4,2),(nquadpts-1)/2))
quadwts[nquadpts] <- 1
quadwts <- quadwts*delta/3
# Second derivative of basis functions at quadrature points
Q2basismat = eval.basis(quadpts, basis,2);
# estimates for basis coefficients
Rmat = t(Q2basismat)%*%(Q2basismat*(quadwts%*%t(rep(1,nbasis))))
dim(Rmat)
## [1] 13 13
basismat2 = t(basismat)%*%basismat;
lambda = 0.05 # smoothing parameter
Bmat = basismat2 + lambda*Rmat;
chat = ginv(Bmat)%*%t(basismat)%*%yobs;
# fitted value
yhat = basismat%*%chat;
yhat2 = basismat%*%ginv(t(basismat)%*%basismat)%*%t(basismat)%*%yobs;
#quartz()
plot(tobs,ytru,type = "l")
points(tobs,yobs)
lines(tobs,yhat,type = "l",col="red")
lines(tobs,yhat2,type = "l",col="blue")
# degrees of freedom
Mmat = ginv(Bmat)%*%t(basismat)
Smat = basismat%*%Mmat
df = sum(diag(Smat))
c(df,nbasis)
## [1] 3.37606 13.00000
# estimate the variance of noise
SSE = t(yhat-yobs)%*%(yhat-yobs)
sigma2 = SSE/(nobs-df)
sigma2
## [,1]
## [1,] 0.0009588998
sqrt(sigma2)
## [,1]
## [1,] 0.03096611
# estimate the variance of the fitted curve
varYhat = diag(Smat%*%Smat*matrix(sigma2,nobs,nobs))
# 95% confidence interval
yhat025 = yhat-1.96*sqrt(varYhat)
yhat975 = yhat+1.96*sqrt(varYhat)
lines(tobs,yhat025,type="l", lty=2, lwd=3,col="blue")
lines(tobs,yhat975,type="l", lty=2, lwd=3,col="blue")
