R/cov.R

In linulysses/mcfda: Mean and Covariance Estimation for Functional Data Analysis

#' Estimate Covariance Function
#' @param t a list of vectors (for irregular design) or a vector (for regular design) containing time points of observations for each individual; each vector should be in ascending order.
#' @param y a list of vectors (for irregular design) or a matrix (for regular design) containing the observed values at \code{t}; if it is a matrix, the columns correspond to the time points in the vector \code{t}.
#' @param newt  a list of vectors or a vector containing time points of observations to be evaluated; if NULL, then \code{newt} is treated as \code{t}.
#' @param mu the known or estimated mean function object; it must be a scalar (viewed as a constant function), a function handle, or an object obtained by calling \code{meanfunc}
#' @param method estimation method, 'PACE' or 'FOURIER' or 'SP' (for semiparametric method)
#' @param weig a vector of \code{length(t)} of weight for each subject, or 'OBS' or 'SUBJ' for weighting scheme
#' @param tuning tuning method to select possible tuning parameters
#' @param ... other parameters required depending on the \code{method} and \code{tuning}; see details
#' @details
#'     \itemize{
#'         \item{When \code{method='PACE'}, additional parameters are}
#'         \describe{
#'             \item{\code{kernel}}{kernel type; supported are 'epanechnikov', "rectangular", "triangular", "quartic", "triweight", "tricube", "cosine", "gauss", "logistic", "sigmoid" and "silverman"; see https://en.wikipedia.org/wiki/Kernel_(statistics) for more details.}
#'             \item{\code{deg}}{degree of the local polynomial regression; currently only \code{deg=1} is supported.}
#'             \item{\code{bw}}{bandwidth}
#'         }
#'         \item{When \code{method='FOURIER'}, additional parameters are}
#'         \describe{
#'             \item{\code{q}}{number of basis functions; if \code{NULL} then selected by \code{tuning} method}
#'             \item{\code{rho}}{roughness penalty parameter; if \code{NULL} then selected by \code{tuning} method}
#'             \item{\code{ext}}{extension margin of Fourier extension; if \code{NULL} then selected by \code{tuning} method}
#'             \item{\code{domain}}{time domain; if \code{NULL} then estimated by \code{(min(t),max(t))}}
#'         }
#'         \item{When \code{method='SP'}, additional parameters are} 
#'         \describe{
#'             \item{\code{domain}}{time domain; if \code{NULL} then estimated by \code{(min(t),max(t))}}
#'             \item{\code{corf}}{function of the form \code{function(theta,x,y)} that specifies the correlation structure with parameter \code{theta}; If NULL then Matern correlation is used.}
#'             \item{\code{sig2e}}{variance of measurement error; if \code{NULL} the automatically calculated by the calling \code{sigma2}}
#'             \item{\code{sig2x}}{variance function; a function or an object generated by \code{varfunc}; if \code{NULL} then automatically estimated by calling \code{varfunc}}
#'             \item{\code{pfunc}}{penalty function on the estimation; Not used yet.}
#'             \item{\code{theta0}}{Initial value for the parameter \code{theta} to be estimated; \code{NULL} by default}
#'             \item{\code{lb}}{vector of lower bound of \code{theta}; if \code{NULL} then set to \code{-Inf} for all coordinate}
#'             \item{\code{ub}}{vector of upper bound of \code{theta}; if \code{NULL} then set to \code{Inf} for all coordinate}
#'             \item{\code{D}}{dimension of \code{theta}}
#'         }
#'     }
#'
#' @return an object of the class 'covfunc' containing necessary information to predict/evaluate the estimated covariance function and the following output:
#'     \itemize{
#'         \item{When \code{method='PACE'}, additional parameters are}
#'         \describe{
#'             \item{\code{fitted}}{fitted value at the grid spanned by \code{newt}}
#'             \item{\code{delta}}{the largest span among all subjects; note that it is not normalized by the span of the whole study.}
#'             \item{\code{bw}}{selected bandwidth by \code{tuning} method if \code{NULL} is the input for \code{bw}.}
#'             \item{\code{mu}}{estimated mean function if \code{NULL} is the input.}
#'         }
#'         \item{When \code{method='FOURIER'}, additional parameters are}
#'         \describe{
#'             \item{\code{fitted}}{fitted value at the grid spanned by \code{newt}}
#'             \item{\code{q}}{selected \code{q} if \code{NULL} is the input}
#'             \item{\code{rho}}{selected \code{rho} if \code{NULL} is the input}
#'             \item{\code{ext}}{selected \code{ext} if \code{NULL} is the input}
#'             \item{\code{C}}{estimated coefficients}
#'             \item{\code{mu}}{estimated mean function if \code{NULL} is the input.}
#'         }
#'         \item{When \code{method='SP'}, additional parameters are} 
#'         \describe{
#'             \item{\code{fitted}}{fitted value at the grid spanned by \code{newt}}
#'             \item{\code{domain}}{time domain; if \code{NULL} then estimated by \code{(min(t),max(t))}.}
#'             \item{\code{rho}}{estimated function of the form \code{function(x,y)} of the correlation structure.}
#'             \item{\code{sig2e}}{estiamted variance of measurement error if \code{NULL} is the input.}
#'             \item{\code{sig2x}}{estiamted variance function if \code{NULL} is the input.}
#'             \item{\code{theta}}{estimated parameters for the correlation structure.}
#'             \item{\code{mu}}{estimated mean function if \code{NULL} is the input.}
#'         }
#'     }
#' 
#' @importFrom Rdpack reprompt
#' @references 
#' \insertRef{Lin2020b}{mcfda}
#' 
#' \insertRef{Lin2020}{mcfda}
#' 
#' \insertRef{Yao2005}{mcfda}
#' 
#' @examples
#' mu <- function(s) sin(2*pi*s)
#' D <- synfd::sparse.fd(mu=mu, X=synfd::gaussian.process(), n=100, m=5)
#' mu.obj <- meanfunc(D$t,D$y,newt=NULL,method='PACE',
#'                    tuning='cv',weig=NULL,kernel='gauss',deg=1)
#' cov.obj <- covfunc(D$t,D$y,newt=NULL,mu=mu.obj,method='FOURIER',
#'                    tuning='cv',weig=NULL,domain=c(0,1))
#' cov.hat <- predict(cov.obj,regular.grid())
#' @export covfunc
#' @useDynLib mcfda
covfunc <- function( t,y,newt=NULL,mu=NULL,weig=NULL,
                      method=c('FOURIER','PACE','SP'),...)
{
    method <- match.arg(method)
    R <- NULL
    if(is.list(t)) # irregular design
    {
        if(!is.list(y)) stop('y must be list when t is list')
        if(length(y) != length(t))
            stop('the length of y must match the length of t')
        if(method=='PACE')
            return(cov.pace(t,y,mu=mu,newt=newt,...))
        else if(method=='SP')
            return(cov.sp(t,y,mu=mu,newt=newt,...))
        else
            return(cov.basis(t,y,mu=mu,newt=newt,...))
    }
    else if(is.vector(t)) # regular design
    {
        if(method=='PACE')
            stop('method=PACE for regular design is not supported yet')
        else if(method=='SP')
            stop('method=SP for regular design is not supported yet')
        else
            return(cov.basis(t,y,mu=mu,newt=newt,...))
    }
    else stop('t must be a list of vectors (for irregular design) or a vector (for regular design')

    class(R) <- 'covfunc'
    return(R)
}



#' predict cov functions at new locations
#' @param covobj the object obtained by calling \code{covfunc}
#' @param newt a vector or a list of vectors of real numbers
#' @return the estimated mean function evaluated at \code{newt}. It has the same format of \code{newt}
#' @export
predict.covfunc <- function(covobj,newt)
{
    pred <- NULL
    if(toupper(covobj$method) == 'PACE')
    {

        pred <- function(newx)
        {
            R <- Lwls2D(bw=covobj$bw,
                                     kern=covobj$kernel,
                                     xin=covobj$x,
                                     yin=covobj$y,
                                     xout1=newx,
                                     xout2=newx,
                                     win=covobj$weig,
                            delta=covobj$delta)

        }

    }
    else if(toupper(covobj$method) == 'FOURIER')
    {
        pred <- function(newt)
        {
            stopifnot(is.vector(newt))
            domain <- covobj$domain
            if(domain[1]!=0 || domain[2]!=1)
            {
                newt <- (newt-domain[1])/(domain[2]-domain[1])
            }

            p <- dim(covobj$C)[1]
            if(covobj$ext > 0) newt <- shrink(newt,covobj$ext)
            B <- evaluate.basis(p,grid=newt)
            return(B %*% covobj$C %*% t(B))
        }

    }
    else if(covobj$method == 'SP')
    {
        pred <- function(newt) # newt is a vector
        {
            stopifnot(is.vector(newt))

            corr.est <- covobj$rho(newt,newt)
            var.hat <- predict(covobj$sig2x,newt)
            sig.hat <- sqrt(var.hat)
            cov.fitted <- corr.est * (sig.hat %*% t(sig.hat))
            return(cov.fitted)
        }
    }
    else stop(paste0('Method ',covobj$method, ' is not recognized.'))


    if(is.list(newt))
    {
        return(lapply(newt, function(newx) pred(newx)))
    }
    else if(is.vector(newt))
    {
        return(pred(newt))
    }
    else stop('newt must be a vector or a list of vectors of real numbers')
}



#' estimate the window width of snippets
#' @export
estimate.delta <- function(Lt)
{

    if(is.list(Lt))
    {
        tmp <- lapply(Lt, function(v) max(v)-min(v))
        return(max(unlist(tmp)))
    }
    else if(is.vector(Lt)) return(max(Lt)-min(Lt))
    else stop('unsupported type of Lt')
}