# R/fdata2basis.R In fda.usc: Functional Data Analysis and Utilities for Statistical Computing

#### Documented in fdata2basissummary.basis.fdata

#' @title Compute fucntional coefficients from functional data represented in a base of functions
#'
#' @aliases fdata2basis summary.basis.fdata
#' @description Compute fucntional coefficients from functional data (\code{\link{fdata}} class object)
#' represented in a basis (fixed of data-driven basis).
#'
#' @param fdataobj \code{\link{fdata}} class object.
#' @param basis  	a functional basis object defining the basis
#' @param method character string, if it is "grid" the fdata object is evaluated in the grid (\code{argvals} of fdata),
#'  if it is "inprod" the basis representation of functional data is computed by inner product
#' @param object \code{basis.fdata} class object calculated by: \code{\link{fdata2basis}}
#' @param draw logical, original curves  and their basis representation are plotted
#' @param index vector, by default (if NULL) the first n curves are plotted, where n = min(4, length(fdataobj)).
#' Otherwise, index vector indicates taht curvesare plotted.
#' @param \dots Further arguments passed to or from other methods.
#' @return fdata2basis \code{fdata2bais} function return:
#' \itemize{
#' \item {coef}{a matrix or two-dimensional array of coefficients.}
#' \item {basis}{basis of \code{\link{fdata}} class evaluated in the same grid of \code{fdataobj}. }
#' }
#'
#' summary function return:
#' \itemize{
#' \item {R}{a matrix with a measure similar to R-sq for each curve aproximation (by row) and number of basis elements (by column).}
#' }
#' @author Manuel Febrero-Bande, Manuel Oviedo de la Fuente  \email{manuel.oviedo@@usc.es}
#' @seealso  Inverse function: \code{\link{gridfdata}}.
# @references
#' @keywords multivariate
#' @examples
#' \dontrun{
#' T <- 71
#' S <- 51
#' tj <- round(seq(0,1,len=T),3)
#' si <- round(seq(0,1,len=S),3)
#' beta1 <- outer(si,tj,function(si,tj){exp(-5*abs((tj-si)/5))})
#' nbasis.s =7
#' nbasis.t=11
#' base.s <- create.fourier.basis(c(0,1),nbasis=nbasis.s)
#' base.t <- create.fourier.basis(c(0,1),nbasis=nbasis.t)
#' y1 <- fdata(rbind(log(1+tj),1-5*(tj-0.5)^2),argvals=tj,rangeval=c(0,1))
#' aa <- fdata2basis(y1,base.t,method="inprod")
#' summary(aa)
#' plot(gridfdata(aa$coefs,aa$basis))
#' lines(y1,lwd=2,col=c(3,4),lty=2)
#' }
#' @export
fdata2basis <- function(fdataobj, basis, method=c("grid","inprod")){
xmean <- NULL
if (is.basis(basis)){
# print(1)
bb <- fdata(t(eval.basis(fdataobj$argvals,basis)), argvals=fdataobj$argvals,rangeval=fdataobj$rangeval) fdobj <- Data2fd(fdataobj$argvals,t(fdataobj$data),basis,lambda=0) xmean <- fdata(fd(apply(fdobj$coefs,1,mean),basis),argvals=fdataobj$argvals) fdataobj <- fdataobj -xmean } else{ if (class(basis) %in% c("fdata")){ bb <- basis Xcen <- fdata.cen(fdataobj) fdataobj <- Xcen$Xcen
xmean <- Xcen$meanX #xmean <- func.mean(fdataobj) } else { bb = basis$basis
xmean <- basis$mean fdataobj$data <- sweep(fdataobj$data,2,xmean$data,"-")
}
}
#print(2)
if (method[1]=="grid"){
#print(3)
A=t(fdataobj$data) B=t(bb$data)
coefs=t(Minverse(t(B)%*%B)%*%t(B)%*%A)
} else {
# print(4)
A=t(inprod.fdata(fdataobj,bb))
B=inprod.fdata(bb)
coefs=t(solve(B)%*%A)
}
#print(5)
if (!is.null(basis$type)) type <- basis$type
#if (!is.null(bb$names$main)) type <- bb$names$main
out <- list(coefs=coefs,basis=bb,
fdataobj=fdataobj,
type=type,
mean = xmean
)
class(out) <- "basis.fdata"
return(out)
}

#' @rdname fdata2basis
#' @export
summary.basis.fdata=function(object, draw=TRUE, index=NULL,...) {
#object <- bb
cat("\n     - SUMMARY -\n")
le <- length(object$basis) R <- numeric(le) n <- nrow(object$fdataobj)
R2 <- matrix(NA,n,le)
colnames(R2)<-paste0(substr(object$type,1,3),"(1:",1:le,")") rownames(R2)<-rownames(object$fdataobj)
#  Xcen <- fdata.cen(object$fdataobj)$Xcen
type=FALSE
cat("Type of basis:",object$type,"\nNum. of basis:",le,"\nRangeval:",object$basis$rangeval,"\n") for (l in 1:le){ # if (object$type =="pc" | object$type == "pls"){ # type=TRUE #xmean <- func.mean(object$fdataobj)
#      xmean <- object$mean fdata.est <- gridfdata(object$coefs[,1:l,drop=F],
object$basis[1:l],mu=object$mean)
#    } else{ fdata.est <- gridfdata(object$coefs[,1:l,drop=F], # object$basis[1:l],mu=object$mean) #R2[l] <- 1 - sum(norm.fdata(fdata.est-object$fdataobj)^2)/
#   sum(norm.fdata(fdata.cen(object$fdataobj)$Xcen)^2)
#    }
# fdataobj from fdata2basis comes recentred
R[l] <- 1-sum(norm.fdata(fdata.est-object$fdataob-object$mean)^2)/sum(norm.fdata(object$fdataobj)^2) R2[,l] <- 1 - norm.fdata(fdata.est-object$fdataobj-object$mean)^2/norm.fdata(object$fdataobj)^2
}
if (is.null(index)) {
ymin <- min(4,n);     index<- 1:ymin
}
txt <- paste(length(index),"curves (in grey) and their basis representation (in red) are plotted\n")
if (draw) {
# yl <- c(min(object$fdataobj,fdata.est),max(object$fdataobj,fdata.est))
yl <- c(min(object$fdataobj[index]+object$mean,fdata.est[index]),max(object$fdataobj[index]+object$mean,fdata.est[index]))
# mn <- expression( 1 - frac(paste( "||X - ",hat(X),"||"),paste( "||X - ",bar(X),"||")))
mn <-  bquote( paste("X(t) vs ",hat(X),"(t)-(",.(min(index)),":",.(max(index)),")"))
plot(object$fdataobj[index]+object$mean,main=mn,col="grey",ylim=yl)
lines(fdata.est[index],lty=2,col=2)
}

cat("\nMeasure of fit: 1 - Sum||X - hat(X)||^2  / Sum||X - bar(X)||^2:\n")
names(R) <- paste0(object$type,"(1:",1:le,")",sep="") # if (!type) R <- R[le] #print(expression( 1 - paste( "||X - ",hat(X),"|| / ||X - ",bar(X),"||"))) print(R) # if (draw) cat("\n",txt) return(invisible(R2)) } #wrapper of gridfdata basis2fdata<-function (coef, fdataobj, mu) { return(gridfdata(coef, fdataobj, mu) ) } # Devolver la matriz entera R2 para todas las curvas en un summary # parece excesivamente prolijo. Como texto habría que devolver un R2 # conjunto con el tamaño global de la base: # 1-sum(norm.fdata(.,.)^2)/sum(norm.fdata(Xcen)^2). # La función debe devolver la matriz R2 como invisible por si # se quiere usar posteriormente pero devolver el mismo objeto # que entra no tiene mucha utilidad. Respecto al texto que se # incluye al principio (solo la llamada de la función) no sé # si sería mejor poner un texto con el object$type, la longitud
# y alguna info sobre los argvals y el rangeval.

# no ayuda por el momento
fdata2basis2d=function(fdata2d,basis.s,basis.t){
n=dim(fdata2d$data)[1] if (is.basis(basis.s)) { basis.s=fdata(t(eval.basis(fdata2d$argvals[[1]],basis.s)),
argvals=fdata2d$argvals[[1]],rangeval=fdata2d$rangeval[[1]])
}
if (is.basis(basis.t)) {
basis.t=fdata(t(eval.basis(fdata2d$argvals[[2]],basis.t)), argvals=fdata2d$argvals[[2]],rangeval=fdata2d$rangeval[[2]]) } n.s=nrow(basis.s) n.t=nrow(basis.t) eta=basis.s$data;theta=basis.t$data Jeta=tcrossprod(eta);Jtheta=tcrossprod(theta) AA=Minverse(Jeta)%*%eta BB=t(theta)%*% Minverse(Jtheta) B=array(NA,dim=c(n,n.s,n.t)) #B=solve(Jeta)%*%eta%*%fdata2d%*%t(theta)%*%solve(Jtheta) for (i in 1:n){B[i,,]=AA%*%fdata2d$data[i,,]%*%BB}
return(list(coefs=B,b.s=basis.s,b.t=basis.t))
}

# no ayuda por el momento
basis2d2fdata=function(coefs,b.s,b.t){
if (!is.array(coefs))
stop("coefs must be an array")
if (dim(coefs)[2]!=length(b.s) | dim(coefs)[3]!=length(b.t))
stop("Dimensions of coefs and basis.s or basis.t do not match")
n=dim(coefs)[1]
n.s=length(b.s$argvals) n.t=length(b.t$argvals)
Beta=array(NA,dim=c(n,n.s,n.t))
for (i in 1:n){
Beta[i,,]=t(b.s$data)%*%coefs[i,,]%*%b.t$data
}
Beta=fdata(Beta,argvals=list(b.s$argvals,b.t$argvals),
rangeval=list(b.s$rangeval,b.t$rangeval),
names=list(main="B(s,t)",xlab=list("s","t")),
fdata2d=TRUE)
}

# no ayuda
intXbeta=function(X,Beta,arg.t=1:ncol(Beta),
rangeval=range(arg.t),
equi=argvals.equi(X$argvals), method=NULL){ if (is.null(method)) { par.fda.usc=eval(parse(text="fda.usc:::par.fda.usc"), envir=.GlobalEnv) method=par.fda.usc$int.method
}
if (inherits(X,"fdata")) {
nX=NROW(X)
arg.s=X$argvals equi=argvals.equi(arg.s) n.s=length(X$argvals)
if (!is.matrix(Beta)) { stop("Only implemented for Beta being a matrix")} else {
if (nrow(Beta)!=n.s) stop ("Dimensions of X and Beta do not coincide")
if (ncol(Beta)!=length(arg.t)) stop("Beta and argvals.t do not have the same length")
n.t=ncol(Beta);arg.t=arg.t
}
if (equi) {
res=X$data%*%Beta } else{ res=matrix(NA,nrow=nX,ncol=n.t) for (i in 1:nX){ for (j in 1:n.t){ res[i,j]=int.simpson2(arg.s,X$data[i,]*Beta[,j],equi=equi,method=method)
}}
}
nms <- list(main="X(s)Beta(s,t)",xlab="t",ylab="Y(t)")
res=fdata(res,argvals=arg.t,rangeval=rangeval,names=nms)
}
}

# gridfdata2d=function(coefs,fdataobj,mu){
# if (!is.fdata(fdataobj) | length(fdataobj$argvals)!=2) stop("fdataobj is not an fdata2d object") # if (is.matrix(coefs)) { # n=1 # if (nrow(coefs)!=length(fdataobj$argvals[[1]]) | ncol(coefs)!=length(fdataobj$argvals[[2]])) stop("Dimensions of coefs and fdataobj do not coincide") # n.s=nrow(coefs);n.t=ncol(coefs) # coefs=array(coefs,dim=c(1,nrow=n.s,ncol=n.t)) # } else if (is.array(coefs)){ # n=dim(coefs)[1] # if (dim(coefs)[2]!=length(fdataobj$argvals[[1]]) | dim(coefs)[3]!=length(fdataobj$argvals[[2]])) stop("Dimensions of coefs and fdataobj do not coincide") # n.s=dim(coefs)[2] # n.t=dim(coefs)[3] # } else stop("coefs must be a matrix or an array") # res=array(NA,dim=c( # if (!missing(mu)) {mu=fdata(array(rep(0,n.s*n.t),dim=c(1,n.s,n.t)),argvals=fdataobj$argvals,
# rangeval=fdataobj$rangeval, names=fdataobj$names,fdata2d=TRUE)}

# }

#' Minverse <- fda.usc:::Minverse
#' T <- 71
#' S <- 51
#' tj <- round(seq(0,1,len=T),3)
#' si <- round(seq(0,1,len=S),3)
#' beta1 <- outer(si,tj,function(si,tj){exp(-5*abs((tj-si)/5))})
#' nbasis.s =31
#' nbasis.t=51
#' base.s <- create.fourier.basis(c(0,1),nbasis=nbasis.s)
#' base.t <- create.fourier.basis(c(0,1),nbasis=nbasis.t)
#----------------------------------------------
#' y1 <- fdata(rbind(log(1+tj),1-5*(tj-0.5)^2),argvals=tj,rangeval=c(0,1))
#' aa <- fdata2basis(y1,base.t,method="inprod")
#' plot(gridfdata(aa$coefs,aa$basis))
#' lines(y1,lwd=2,col=c(3,4),lty=2)
#----------------------------------------------
#' z <- array(NA,dim=c(2,S,T))
#' for (i in 1:2) {z[i,,] <- outer(si, tj, function(u, v) {(i*u)*(v-0.5)^i})}
#' z.fdata <- fdata(z,argvals=list(si,tj),rangeval=list(c(0,1),c(0,1)),fdata2d=TRUE)
#' aa <- fdata2basis2d(z.fdata,base.s,base.t)
#' z2.fdata <- basis2d2fdata(aa$coefs,aa$b.s,aa$b.t) #' plot(z.fdata) # Original #' plot(z2.fdata) #Reconstruida #' par(mfrow=c(1,2)) #' plot(z.fdata,type="persp") #' plot(z2.fdata,type="persp") #------------------------------------------- #------------------------------------------- #' xx=fdata(rbind(-1/((si-0.5)^2+0.25),si,si^2),argvals=si,rangeval=c(0,1)) #' miBeta=outer(si,tj,function(u,v){-(v-0.25)^2}) #' miBeta=fdata(miBeta,argvals=list(si,tj),rangeval=list(c(0,1),c(0,1)),fdata2d=TRUE) #' yy=intXbeta(xx,miBeta$data[1,,],tj,rangeval=c(0,1))
#' par(mfrow=c(1,2))
#' plot(xx);plot(yy)

# eta=t(eval.basis(si,base.s))
# theta=t(eval.basis(tj,base.t))

# Jeta=tcrossprod(eta)
# Jtheta=tcrossprod(theta)

# B=solve(Jeta)%*%eta%*%beta1%*%t(theta)%*%solve(Jtheta)

# betatilde=t(eta)%*%B%*%theta
# par(mfrow=c(1,2))

# image(si,tj,beta1,main="beta")
# image(si,tj,betatilde,main="betatilde")
#--------------------------
# npc.x=5
# npc.y=5

# Sigma.X=outer(si,si,function(u,v){exp(-abs(u-v))})
# Sigma.Y=outer(tj,tj,function(u,v){exp(-abs(u-v))})
# res.X=eigen(Sigma.X)
# bpc.X=fdata(t(res.X$vector[,1:npc.x]),argvals=si) # bpc.X[["data"]]=sweep(bpc.X[["data"]],1,norm.fdata(bpc.X),"/") # res.Y=eigen(Sigma.Y) # bpc.Y=fdata(t(res.Y$vector[,1:npc.y]),argvals=tj)
# bpc.Y[["data"]]=sweep(bpc.Y[["data"]],1,norm.fdata(bpc.Y),"/")

# X=rproc2fdata(500,t=si,sigma="vexponential")
# Y=rproc2fdata(500,t=tj,sigma="vexponential")
# bpc.X=create.pc.basis(X,1:npc.x,lambda=10)$basis # bpc.Y=create.pc.basis(Y,1:npc.y,lambda=10)$basis
# eta.pc=bpc.X$data # theta.pc=bpc.Y$data
# Jeta.pc=tcrossprod(eta.pc)
# Jtheta.pc=tcrossprod(theta.pc)

# B2=solve(Jeta.pc)%*%eta.pc%*%beta1%*%t(theta.pc)%*%solve(Jtheta.pc)
# betatilde2=t(eta.pc)%*%B2%*%theta.pc
# par(mfrow=c(1,2))
# image(si,tj,beta1,main="beta")
# image(si,tj,betatilde2,main="betatilde2")


## Try the fda.usc package in your browser

Any scripts or data that you put into this service are public.

fda.usc documentation built on Oct. 17, 2022, 9:06 a.m.