# R/band_depths.R In roahd: Robust Analysis of High Dimensional Data

#### Documented in BDBD.defaultBD.fDataBD_relativeBD_relative.defaultBD_relative.fDataMBDMBD.defaultMBD.fDataMBD_relativeMBD_relative.defaultMBD_relative.fData

#' Band Depth for univariate functional data
#'
#' This function computes the Band Depth (BD) of elements of a functional
#' dataset.
#'
#' Given a univariate functional dataset, \eqn{X_1(t), X_2(t), \ldots, X_N(t)},
#' this function computes the sample BD of each element with respect to the
#' other elements of the dataset, i.e.:
#'
#' \deqn{ BD( X( t ) ) = {N \choose 2 }^{-1} \sum_{1 \leq i_1 < i_2 \leq N} I(
#' G(X) \subset B( X_{i_1}, X_{i_2} ) ),}
#'
#' where \eqn{G(X)} is the graphic of \eqn{X(t)}, \eqn{B(X_{i_1},X_{i_2})} is
#' the envelope of \eqn{X_{i_1}(t)} and \eqn{X_{i_2}(t)}, and \eqn{X \in
#' \left\{X_1, \ldots, X_N\right\}}.
#'
#' See the References section for more details.
#'
#'
#' @param Data either an object of class \code{fData} or a matrix-like dataset
#' of functional data (e.g. \code{fData$values}), #' with observations as rows and measurements over grid points as columns. #' #' @return The function returns a vector containing the values of BD for the #' given dataset. #' #' @references #' #' Lopez-Pintado, S. and Romo, J. (2009). On the Concept of Depth for Functional #' Data, \emph{Journal of the American Statistical Association}, 104, 718-734. #' #' Lopez-Pintado, S. and Romo. J. (2007). Depth-based inference for functional #' data, \emph{Computational Statistics & Data Analysis} 51, 4957-4968. #' #' @seealso \code{\link{MBD}}, \code{\link{BD_relative}}, #' \code{\link{MBD_relative}}, \code{\link{fData}} #' #' @examples #' #' grid = seq( 0, 1, length.out = 1e2 ) #' #' #' D = matrix( c( 1 + sin( 2 * pi * grid ), #' 0 + sin( 4 * pi * grid ), #' 1 - sin( pi * ( grid - 0.2 ) ), #' 0.1 + cos( 2 * pi * grid ), #' 0.5 + sin( 3 * pi + grid ), #' -2 + sin( pi * grid ) ), #' nrow = 6, ncol = length( grid ), byrow = TRUE ) #' #' fD = fData( grid, D ) #' #' BD( fD ) #' #' BD( D ) #' #' @export #' BD = function( Data ) { UseMethod( 'BD', Data ) } #' @rdname BD #' #' @export BD.fData = function( Data ) { Data = Data$values
NextMethod()
}

#' @rdname BD
#'
#' @export
BD.default = function( Data )
{
# Number of rows
N = nrow( Data )

# Actual number of function values strictly above
N_a = apply( Data, 1, function( x ) sum( apply( Data, 1, function( s ) all( x < s ) ) ) )

# Actual number of function values strictly below
N_b = apply( Data, 1, function( x ) sum( apply( Data, 1, function( s ) all( x > s ) ) ) )

Depths = ( N_a * N_b + ( N - 1 ) ) / ( N * ( N - 1 ) / 2 )

return( Depths )
}

#' Relative Band Depth of functions in a univariate functional dataset
#'
#' This function computes Band Depth (BD) of elements of a univariate functional
#' dataset with respect to another univariate functional dataset.
#'
#' Given a univariate functional dataset of elements \eqn{X_1(t), X_2(t),
#' \ldots, X_N(t)}, and another univariate functional dataset of elements
#' \eqn{Y_1(t), Y_2(t) \ldots, Y_M(t)}, this function computes the BD of
#' elements of the former with respect to elements of the latter, i.e.:
#'
#'  \deqn{ BD( X_i( t ) ) = {M \choose 2 }^{-1} \sum_{1 \leq i_1 < i_2 \leq M} I(
#' G(X_i) \subset B( Y_{i_1}, Y_{i_2} ) ), \quad \forall i = 1, \ldots, N,}
#'
#' where \eqn{G(X_i)} is the graphic of \eqn{X_i(t)} and \eqn{B(Y_{i_1},Y_{i_2})} is
#' the envelope of \eqn{Y_{i_1}(t)} and \eqn{Y_{i_2}(t)}.
#'
#' @param Data_target is the univariate functional dataset, provided either as
#' a \code{fData} object or in matrix
#' form (N observations as rows and P measurements as columns), whose BD
#' have to be computed with respect to the reference dataset.
#' @param Data_reference is the dataset, provided either as a \code{fData}
#' object or in matrix form (N observations
#' as rows and P measurements as columns), containing the reference
#' univariate functional data that must be used to compute the band depths of
#' elements in \code{Data_target}. If \code{Data_target} is \code{fData}, it
#' must be of class \code{fData}
#'
#'
#' @return The function returns a vector containing the BD of elements in
#' \code{Data_target} with respect to elements in \code{Data_reference}.
#'
#'
#' @examples
#'
#' grid = seq( 0, 1, length.out = 1e2 )
#'
#' Data_ref = matrix( c( 0  + sin( 2 * pi * grid ),
#'                       1  + sin( 2 * pi * grid ),
#'                       -1 + sin( 2 * pi * grid ) ),
#'                    nrow = 3, ncol = length( grid ), byrow = TRUE )
#'
#' Data_test_1 = matrix( c( 0.6 + sin( 2 * pi * grid ) ),
#'                       nrow = 1, ncol = length( grid ), byrow = TRUE )
#'
#' Data_test_2 = matrix( c( 0.6 + sin( 2 * pi * grid ) ),
#'                       nrow = length( grid ), ncol = 1, byrow = TRUE )
#'
#' Data_test_3 = 0.6 + sin( 2 * pi * grid )
#'
#' Data_test_4 = array( 0.6 + sin( 2 * pi * grid ), dim = length( grid ) )
#'
#' Data_test_5 = array( 0.6 + sin( 2 * pi * grid ), dim = c( 1, length( grid ) ) )
#'
#' Data_test_6 = array( 0.6 + sin( 2 * pi * grid ), dim = c( length( grid ), 1 ) )
#'
#' Data_test_7 = matrix( c( 0.5  + sin( 2 * pi * grid ),
#'                          -0.5 + sin( 2 * pi * grid ),
#'                          1.1 + sin( 2 * pi * grid ) ),
#'                       nrow = 3, ncol = length( grid ), byrow = TRUE )
#'
#' fD_ref = fData( grid, Data_ref )
#' fD_test_1 = fData( grid, Data_test_1 )
#' fD_test_2 = fData( grid, Data_test_2 )
#' fD_test_3 = fData( grid, Data_test_3 )
#' fD_test_4 = fData( grid, Data_test_4 )
#' fD_test_5 = fData( grid, Data_test_5 )
#' fD_test_6 = fData( grid, Data_test_6 )
#' fD_test_7 = fData( grid, Data_test_7 )
#'
#' BD_relative( fD_test_1, fD_ref )
#' BD_relative( Data_test_1, Data_ref )
#'
#' BD_relative( fD_test_2, fD_ref )
#' BD_relative( Data_test_2, Data_ref )
#'
#' BD_relative( fD_test_3, fD_ref )
#' BD_relative( Data_test_3, Data_ref )
#'
#' BD_relative( fD_test_4, fD_ref )
#' BD_relative( Data_test_4, Data_ref )
#'
#' BD_relative( fD_test_5, fD_ref )
#' BD_relative( Data_test_5, Data_ref )
#'
#' BD_relative( fD_test_6, fD_ref )
#' BD_relative( Data_test_6, Data_ref )
#'
#' BD_relative( fD_test_7, fD_ref )
#' BD_relative( Data_test_7, Data_ref )
#'
#' @export
#'
BD_relative = function( Data_target, Data_reference )
{
if( any(class( Data_target ) != class( Data_reference )) )
{
if( ! all( class( Data_target ) %in% c( 'numeric', 'array', 'matrix' ) &
class( Data_reference ) %in% c( 'numeric', 'array', 'matrix' ) ) )
{
stop( 'Error in BD_relative: you have to provide target and reference data
with the same class')
}
}

UseMethod( 'BD_relative', Data_target )
}

#' @rdname BD_relative
#'
#' @export
BD_relative.fData = function( Data_target, Data_reference )
{
Data_target = Data_target$values Data_reference = Data_reference$values
NextMethod()
}

#' @rdname BD_relative
#'
#' @export
BD_relative.default = function( Data_target, Data_reference )
{
# Observations
N = nrow( Data_reference )

# Number of columns
P = ncol( Data_reference )

Data_target = toRowMatrixForm( Data_target )

if( ncol( Data_target ) != P ) stop( 'Error: you provided a Data_target with
not compliant dimensions to BD_relative')

# Number of target_functions
N_target = nrow( Data_target )

Depths = rep( 0, N_target );

for( iObs in 1 : N_target )
{

# Actual number of function values strictly above
N_a = sum( apply( Data_reference,
1,
function( x ) all( Data_target[ iObs, ] < x ) ) )

# Actual number of function values strictly below
N_b = sum( apply( Data_reference,
1,
function( x ) all( Data_target[ iObs, ] > x ) ) )

Depths[ iObs ] = ( N_a * N_b ) / ( N * ( N - 1 ) / 2 )
}

return( Depths )
}

#' Modified Band Depth for univariate functional data
#'
#' This function computes the Modified Band Depth (MBD) of elements of a
#' functional dataset.
#'
#' Given a univariate functional dataset, \eqn{X_1(t), X_2(t), \ldots, X_N(t)},
#' defined over a compact interval \eqn{I= [a,b]},
#' this function computes the sample MBD of each element with respect to the
#' other elements of the dataset, i.e.:
#'
#' \deqn{ MBD( X( t ) ) = {N \choose 2 }^{-1} \sum_{1 \leq i_1 < i_2 \leq N}
#' \tilde{\lambda}\big( {t : \min( X_{i_1}(t), X_{i_2}(t) ) \leq X(t) \leq
#' \max( X_{i_1}(t), X_{i_2}(t) ) } \big), }
#'
#' where \eqn{\tilde{\lambda}(\cdot)} is the normalized Lebesgue measure over
#' \eqn{I=[a,b]}, that is \eqn{\tilde{\lambda(A)} = \lambda( A ) / ( b - a )}.
#'
#' See the References section for more details.
#'
#'
#' @param Data either a \code{fData} object or a matrix-like dataset of functional
#' data (e.g. \code{fData$values}), with observations as rows and measurements #' over grid points as columns. #' @param manage_ties a logical flag specifying whether a check for ties and #' relative treatment must be carried out or not (default is \code{FALSE}). #' #' @return The function returns a vector containing the values of MBD for the #' given dataset. #' #' @references #' #' Lopez-Pintado, S. and Romo, J. (2009). On the Concept of Depth for Functional #' Data, \emph{Journal of the American Statistical Association}, 104, 718-734. #' #' Lopez-Pintado, S. and Romo. J. (2007). Depth-based inference for functional #' data, \emph{Computational Statistics & Data Analysis} 51, 4957-4968. #' #' @seealso \code{\link{BD}}, \code{\link{MBD_relative}}, #' \code{\link{BD_relative}}, \code{\link{fData}} #' #' @examples #' grid = seq( 0, 1, length.out = 1e2 ) #' #' #' D = matrix( c( 1 + sin( 2 * pi * grid ), #' 0 + sin( 4 * pi * grid ), #' 1 - sin( pi * ( grid - 0.2 ) ), #' 0.1 + cos( 2 * pi * grid ), #' 0.5 + sin( 3 * pi + grid ), #' -2 + sin( pi * grid ) ), #' nrow = 6, ncol = length( grid ), byrow = TRUE ) #' #' fD = fData( grid, D ) #' #' MBD( fD ) #' #' MBD( D ) #' #' @export MBD = function( Data, manage_ties = FALSE ) { UseMethod( 'MBD', Data ) } #' @rdname MBD #' #' @export MBD.fData = function( Data, manage_ties = FALSE ) { Data = Data$values
NextMethod()
}

#' @rdname MBD
#'
#' @export
MBD.default = function( Data, manage_ties = FALSE )
{
# Number of rows
N = nrow( Data )

# Number of time points in the discrete grid
P = ncol( Data )

if( manage_ties )
{
# Compute ranks of matrix-like representation of data with min'
# tie-breaking rule
rk_min = apply( Data, 2, function( v )( rank( v, ties.method = 'min' ) ) )

# Compute ranks of matrix-like representation of data with max'
# tie-breaking rule
rk_max = apply( Data, 2, function( v )( rank( v, ties.method = 'max' ) ) )

# Times each function value is repeated in the dataset, for each time point
# ( matrix is N x P)
Repetitions = rk_max - rk_min + 1

# Actual number of function values strictly above
N_a = N - rk_max

# Actual number of function values strictly below
N_b = rk_min - 1

# if( any( Repetitions > 1 ) )
# {
# Now, owing to repetitions we gain more bands containing each functional
# observation.
# In standard cases (no coincident values) we have N - 1, in general we
# have N - 1 + N - 2 + ... + N - K, where K is the number of repetitions
# of the function value in the dataset, time by time.
# If you split the sum and use Gauss formula, you can get the following
# expression

added_bands = N * Repetitions - 0.5 * ( Repetitions * Repetitions +
Repetitions )

Depths = rowSums( N_a * N_b  + added_bands ) / ( P * ( N - 1 ) * N / 2 )

} else {

rk =  apply( Data, 2, function( v )( rank( v ) )  )

Depths = ( rowSums( ( N - rk ) * ( rk - 1 ) ) / P +
N - 1 ) / ( N * ( N - 1 ) / 2 )

}

return( Depths )
}

#' Relative Modified Band Depth of functions in a univariate functional dataset
#'
#' This function computes Modified Band Depth (BD) of elements of a univariate
#' functional dataset with respect to another univariate functional dataset.
#'
#' Given a univariate functional dataset of elements \eqn{X_1(t), X_2(t),
#' \ldots, X_N(t)}, and another univariate functional dataset of elements
#' \eqn{Y_1(t), Y_2(t) \ldots, Y_M(t)}, defined over the same compact interval
#' \eqn{I=[a,b]}, this function computes the MBD of
#' elements of the former with respect to elements of the latter, i.e.:
#'
#'\deqn{ MBD( X_i( t ) ) = {M \choose 2 }^{-1} \sum_{1 \leq i_1 < i_2 \leq M}
#' \tilde{\lambda}\big( {t : \min( Y_{i_1}(t), Y_{i_2}(t) ) \leq X_i(t) \leq
#' \max( Y_{i_1}(t), Y_{i_2}(t) ) } \big),}
#'
#' \eqn{\forall i = 1, \ldots, N}, where \eqn{\tilde{\lambda}(\cdot)} is the
#' normalized Lebesgue measure over \eqn{I=[a,b]}, that is
#' \eqn{\tilde{\lambda(A)} = \lambda( A ) / ( b - a )}.
#'
#' @param Data_target is the univariate functional dataset, provided either as
#' an \code{fData} object or in matrix
#' form (N observations as rows and P measurements as columns), whose MBD
#' have to be computed with respect to the reference dataset.
#' @param Data_reference is the dataset, provided either as an \code{fData}
#' object or in matrix form (N observations
#' as rows and P measurements as columns), containing the reference
#' univariate functional data that must be used to compute the MBD of
#' elements in \code{Data_target}. If \code{Data_target} is \code{fData}, it
#' must be of class \code{fData}.
#'
#'
#' @return The function returns a vector containing the MBD of elements in
#' \code{Data_target} with respect to elements in \code{Data_reference}.
#'
#'
#' @examples
#'
#' grid = seq( 0, 1, length.out = 1e2 )
#'
#' Data_ref = matrix( c( 0  + sin( 2 * pi * grid ),
#'                       1  + sin( 2 * pi * grid ),
#'                       -1 + sin( 2 * pi * grid ) ),
#'                    nrow = 3, ncol = length( grid ), byrow = TRUE )
#'
#' Data_test_1 = matrix( c( 0.6 + sin( 2 * pi * grid ) ),
#'                       nrow = 1, ncol = length( grid ), byrow = TRUE )
#'
#' Data_test_2 = matrix( c( 0.6 + sin( 2 * pi * grid ) ),
#'                       nrow = length( grid ), ncol = 1, byrow = TRUE )
#'
#' Data_test_3 = 0.6 + sin( 2 * pi * grid )
#'
#' Data_test_4 = array( 0.6 + sin( 2 * pi * grid ), dim = length( grid ) )
#'
#' Data_test_5 = array( 0.6 + sin( 2 * pi * grid ), dim = c( 1, length( grid ) ) )
#'
#' Data_test_6 = array( 0.6 + sin( 2 * pi * grid ), dim = c( length( grid ), 1 ) )
#'
#' Data_test_7 = matrix( c( 0.5  + sin( 2 * pi * grid ),
#'                          -0.5 + sin( 2 * pi * grid ),
#'                          1.1 + sin( 2 * pi * grid ) ),
#'                       nrow = 3, ncol = length( grid ), byrow = TRUE )
#'
#' fD_ref = fData( grid, Data_ref )
#' fD_test_1 = fData( grid, Data_test_1 )
#' fD_test_2 = fData( grid, Data_test_2 )
#' fD_test_3 = fData( grid, Data_test_3 )
#' fD_test_4 = fData( grid, Data_test_4 )
#' fD_test_5 = fData( grid, Data_test_5 )
#' fD_test_6 = fData( grid, Data_test_6 )
#' fD_test_7 = fData( grid, Data_test_7 )
#'
#' MBD_relative( fD_test_1, fD_ref )
#' MBD_relative( Data_test_1, Data_ref )
#'
#' MBD_relative( fD_test_2, fD_ref )
#' MBD_relative( Data_test_2, Data_ref )
#'
#' MBD_relative( fD_test_3, fD_ref )
#' MBD_relative( Data_test_3, Data_ref )
#'
#' MBD_relative( fD_test_4, fD_ref )
#' MBD_relative( Data_test_4, Data_ref )
#'
#' MBD_relative( fD_test_5, fD_ref )
#' MBD_relative( Data_test_5, Data_ref )
#'
#' MBD_relative( fD_test_6, fD_ref )
#' MBD_relative( Data_test_6, Data_ref )
#'
#' MBD_relative( fD_test_7, fD_ref )
#' MBD_relative( Data_test_7, Data_ref )
#'
#'
#' @export
#'
MBD_relative = function( Data_target, Data_reference )
{
if( any(class( Data_target ) != class( Data_reference )) )
{
if( ! all( class( Data_target ) %in% c( 'numeric', 'array', 'matrix' ) &
class( Data_reference ) %in% c( 'numeric', 'array', 'matrix' ) ) )
{
stop( 'Error in MBD_relative: you have to provide target and reference data
with the same class')
}
}

UseMethod( 'MBD_relative', Data_target )
}

#' @rdname MBD_relative
#'
#' @export
MBD_relative.fData = function( Data_target, Data_reference )
{
Data_target = Data_target$values Data_reference = Data_reference$values
NextMethod()
}

#' @rdname MBD_relative
#'
#' @export
MBD_relative.default = function( Data_target, Data_reference )
{
# Observations
N = nrow( Data_reference )

# Number of columns
P = ncol( Data_reference )

Data_target = toRowMatrixForm( Data_target )

if( ncol( Data_target ) != P ) stop( 'Error: you provided a Data_target with
not compliant dimensions to MBD_relative')

# Number of target_functions
N_target = nrow( Data_target )

Depths = rep( 0, N_target );

for( iObs in 1 : N_target )
{
rk = apply( rbind( Data_target[ iObs, ], Data_reference ),
2, rank, ties.method = 'average' );

N_a = N + 1 - rk[ 1, ]

N_b = rk[ 1, ] - 1

Depths[ iObs ] = sum( N_a * N_b ) / (  P *  N * ( N - 1 ) / 2 )
}

return( Depths )
}


## Try the roahd package in your browser

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

roahd documentation built on Nov. 4, 2021, 1:07 a.m.