R/depth.RPD.R

In ftsa: Functional Time Series Analysis

depth.RPD <- function (data, nproj = 50, deriv = c(0, 1), trim = 0.25, 
    dfunc2 = depth.mode, x = NULL, spline = TRUE, graph = FALSE, 
    ...) 
{
    functions = t(data$y)
    if (is.data.frame(functions)) 
        functions = as.matrix(functions)
    nr = nrow(functions)
    nc = ncol(functions)
    modulo = function(z) {
        sqrt(sum(z^2))
    }
    if (is.null(nr) || is.null(nc)) 
        stop("Input must be a matrix")
    if (is.null(x)) 
        x = 1:nc
    newfunc = array(NA, dim = c(nr, nc, length(deriv)))
    for (k in 1:length(deriv)) {
        if (deriv[k] == 0) {
            newfunc[, , k] = functions
        }
        else {
            for (i in 1:nr) {
                if (spline) {
                  fn = splinefun(x, functions[i, ])
                  newfunc[i, , k] = fn(x, deriv = deriv[k])
                }
                else {
                  newfunc[i, , k] = c(rep(0, deriv[k]), diff(functions[i, 
                    ], differences = deriv[k]))
                }
            }
        }
    }
    prof = rep(0, nr)
    vproject = matrix(0, nrow = nr, ncol = length(deriv))
    z = rnorm(nc * nproj)
    z = matrix(z, nrow = nproj, ncol = nc)
    modu = apply(z, 1, modulo)
    z = z/modu
    for (j in 1:nproj) {
        for (k in 1:length(deriv)) {
            matrx = newfunc[, , k]
            vproject[, k] = matrx %*% z[j, ]
        }
        vproject_mat = t(vproject)
        colnames(vproject_mat) = colnames(data$y)
        resul = dfunc2(fts(1:dim(vproject)[2], vproject_mat), ...)
        prof = prof + resul$prof
    }
    prof = prof/nproj
    k = which.max(prof)
    med = functions[k, ]
    lista = which(prof >= quantile(prof, probs = trim))
    mtrim = apply(functions[lista, ], 2, mean)
    if (graph) {
        dev.new()
        cgray = 1 - (prof - min(prof))/(max(prof) - min(prof))
        if (nc == 2) {
            plot(range(functions[, 1]), range(functions[, 2]), 
                type = "n", xlab = colnames(functions)[1], ylab = colnames(functions)[2])
            points(functions[, 1], functions[, 2], col = gray(cgray))
            points(rbind(mtrim), pch = 19, col = gray(2 * trim), 
                cex = 2)
            points(rbind(med), col = 3, pch = 20, cex = 2)
        }
        else {
            plot(range(x), range(functions), type = "n", xlab = "t", 
                ylab = "X(t)", main = "RPD Depth")
            for (i in 1:nr) {
                lines(x, functions[i, ], col = gray(cgray[i]))
            }
            lines(x, mtrim, lwd = 2, col = "yellow")
            lines(x, med, col = "red", lwd = 2)
            legend("bottomleft", legend = c(paste("Trim", trim * 
                100, "%", sep = ""), "Median"), lwd = 2, col = c("yellow", 
                "red"), cex = 0.8)
        }
    }
    return(list(median = med, lmed = k, mtrim = mtrim, ltrim = lista, 
        prof = prof, proj = z))
}