# R/POIFD.R In fdaPOIFD: Partially Observed Integrated Functional Depth

#### Documented in POIFD

```#' Integrated Depth for Partially Observed Functional Data
#'
#' Compute the depth measures of a partially observed functional data set evaluated in a common grid.
#'
#' @param data matrix p by n, being n the number of functions and p the number of grid points.
#' Rownames are the dense grid x and colnames the identifier of each functional data.
#' @param type chosen depth measure. Fraiman and Muniz depth (\code{"FMD"}), Modified band depth (\code{"MBD"})
#' or Modified Half Region Depth and Modified Epigraph/Hipograph Index \code{"MHRD"})
#' @param phi phi function of weights for the POIFD. The default value is as in the paper, i.e. the proportion of observed functions
#' at each time point.
#'
#' @return Ordered vector of depths from the deepest to outward. The names are the functions names (if provided)
#' or the column position.
#'
#' @examples
#' data(exampleData)
#' data <- exampleData\$PoFDintervals
#' poifd <- POIFD(data,  type = c("MBD"))
#'
#' @export
POIFD <- function(data, type = c("MBD", "FMD", "MHRD"), phi){

# if(!missing(xgrid)){
#   data <- getBinnedData_PACE
# }

n <- dim(data)[2]
if(is.null(colnames(data))){colnames(data) <- c(1:n)}

depth <- switch(type, MBD = mbd(data, phi), FMD = fmd(data, phi), MHRD = mhrd(data, phi))

return(depth)
}

mbd <- function(data, phi){
n <- dim(data)[2]
p <- dim(data)[1]

if(missing("phi") ){
phi_num <- (n - rowSums(is.na(data)))
}else{
phi_num <- phi*n
}

MBD_it <- t(mapply(aux.mbd, as.data.frame(t(data)), phi_num, USE.NAMES = FALSE))

phi_den <- apply(data, 2, function(x) sum(phi_num[!is.na(x)]) )

depth <- colSums(MBD_it, na.rm = TRUE)/phi_den
names(depth) <- colnames(data)

MBD <- sort(depth, decreasing=TRUE, index.return=TRUE); names(MBD\$x) <- colnames(data)[MBD\$ix]
return(MBD\$x)
}

aux.mbd <- function(data.t, phi_num){
if(phi_num != 0){
empiricalDistribution <- ecdf(data.t)(data.t)

2*empiricalDistribution*(1 - empiricalDistribution)*phi_num
}else{rep(NA, length(data.t))}
}

fmd <- function(data, phi){
n <- dim(data)[2]
p <- dim(data)[1]

if( missing("phi") ){
phi_num <- (n - rowSums(is.na(data)))
}else{
phi_num <- phi*n
}

FM_it <- t(mapply(aux.fmd, as.data.frame(t(data)), phi_num, USE.NAMES = FALSE))

phi_den <- apply(data, 2, function(x) sum(phi_num[!is.na(x)]) )

depth <- colSums(FM_it, na.rm = TRUE)/phi_den
names(depth) <- colnames(data)

FMd <- sort(depth, decreasing=TRUE, index.return=TRUE); names(FMd\$x) <- colnames(data)[FMd\$ix]
return(FMd\$x)
}

aux.fmd <- function(data.t, phi_num){
if(phi_num != 0){(1 - abs(1/2 - ecdf(data.t)(data.t)))*phi_num}else{rep(NA, length(data.t))}
}

mhrd <- function(data, phi){
n <- dim(data)[2]
p <- dim(data)[1]

if( missing("phi") ){
phi_num <- (n - rowSums(is.na(data))) #observed functions at each t
}else{
phi_num <- phi*n
}

HRD_it <- mapply(aux.hrd, as.data.frame(t(data)), phi_num, USE.NAMES = FALSE)

phi_den <- apply(data, 2, function(x) sum(phi_num[!is.na(x)]) )

mhipo <- colSums(matrix(unlist(HRD_it["hipo",]), nrow = p, ncol = n, byrow = TRUE) , na.rm = TRUE)/phi_den
mepi <- colSums(matrix(unlist(HRD_it["epi",]), nrow = p, ncol = n, byrow = TRUE) , na.rm = TRUE)/phi_den
hrd <- apply(cbind(mepi, mhipo), 1, min)

mepi_pofd <- sort(mepi, decreasing=TRUE, index.return=TRUE); names(mepi_pofd\$x) <- colnames(data)[mepi_pofd\$ix]
mhipo_pofd <- sort(mhipo, decreasing=TRUE, index.return=TRUE); names(mhipo_pofd\$x) <- colnames(data)[mhipo_pofd\$ix]
hrd_pofd <- sort(hrd, decreasing=TRUE, index.return=TRUE); names(hrd_pofd\$x) <- colnames(data)[hrd_pofd\$ix]

return(list(mepi = mepi_pofd\$x, mhipo = mhipo_pofd\$x, hrd = hrd_pofd\$x))
}

aux.hrd <- function(data.t, phi_num.t){
if(phi_num.t != 0){
rmat <- rank(data.t, na.last = "keep" , ties.method = "max")  #in ties, both curves get higher rank

up <- phi_num.t - rmat
down <- rmat - 1
mepi <- up/phi_num.t
mhipo <- down/phi_num.t

output <- list(epi = mepi*phi_num.t, hipo = mhipo*phi_num.t)
}else{
output <- list(epi = rep(NA, length(data.t)), hipo = rep(NA, length(data.t)))
}
return(output)
}
```

## Try the fdaPOIFD package in your browser

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

fdaPOIFD documentation built on May 16, 2022, 5:05 p.m.