R/statistics.R

In DistributionalDDA: Distributional Data Analysis in Bayes Spaces

### statistics.R



#' A distributional median
#'
#' The distributional median is the distribution (from a dataset),
#' that minimizes the sum over the Bayes Space distances in the
#' dataset. The sample median selects the distribution from the
#' dataset that minimizes this distance (or maximizes the spatial
#' depth).
#'
#' @param x a BData object
#' @param na.rm if true missing distances will be set to the maximum
#'     distance
#' @param distType "meandistance" for mean distance to datapoints,
#'       "spatialdepth" for the spatial depth.
#' @param ... further parameters to dist
#'
#' @returns An object of the same type as x, but with only the median,
#'     defined as the value with the smallest distance to all others
#'     (distType "mean"), or the largest spatial depth (distType
#'     "spatialdepth").
#'
#' @family caseStatistics
#' @export
#'
#' @examples
#' data(DTD)
#' Xpdf    = BDataPDFs(DTD$shapePDF,basis=B2Unif01Basis5)
#' median(Xpdf)
#'
median.BData = function(x,na.rm=TRUE,...,distType=c("meandistance","spatialdepth")) {
  i <- switch(match.arg(distType),
              "meandistance" =  {
                D <- as.matrix(dist(x,...))
                if( na.rm ) {
                  MDist <- max(D,na.rm=TRUE)
                  D[!is.finite(D)] <- MDist
                }
                which.min(rowMeans(D))
              },
              "spatialdepth"= which.max(spatialDepth(x,...))
  )
  x[i]
}

## envelopes

#' Distance to Center
#'
#' Computes the distance of a datapoint to the center of a dataset, for the definition of
#' central regions.
#'
#'
#' @param x the dataset
#' @param na.rm boolean whether to remove na values in dist
#' @param ... further parameters to dist
#' @param distType "meandistance" for mean distance to datapoints, "median"
#'     for distance to median, "spatialdepth" for one minus spatial depth.
#' @returns A numeric vector of distances of each observation to the centre
#'     of the dataset.
#' @family caseStatistics
#' @export
#'
#' @examples
#' x = rnorm(100)
#'
#' plot(x,distanceToCenter(x,distType="median"),col="green")
#' points(x,distanceToCenter(x,distType="meandist"))
#' points(x,distanceToCenter(x,distType="spatialdepth"),col="red")
#'
#' data(DTD)
#' Xpdf    = BDataPDFs(DTD$shapePDF,basis=B2Unif01Basis5)
#' distanceToCenter(Xpdf)
distanceToCenter = function(x,na.rm=TRUE,...) UseMethod("distanceToCenter",x)
#' @describeIn distanceToCenter default method
#' @export
distanceToCenter.default =function(x,na.rm=TRUE,...,distType=c("meandistance","median","spatialdepth")) {
  distType = match.arg(distType)
  if( distType =="spatialdepth" )
    return(1-spatialDepth(x))
  D = as.matrix(dist(x,...))
  if( na.rm ) {
    MDist = max(D,na.rm=TRUE)
    D[!is.finite(D)]=MDist
  }
  switch(distType,
         median=D[,which.min(rowMeans(D))],
         meandistance=rowMeans(D),
         default=stop("Unimplemented distType for distanceToCenter")
  )
}

#' Spatial Depth in dataset
#'
#' The depth in a dataset as defined by 1 the length of the mean of
#' the unit vectors pointing to all (other) datapoints. For a point
#' with depth 0 the direction to all other points is parallel.  For a
#' point with depth 1, the point corresponds exactly to the median of
#' the dataset.
#' @param x the values for which the depth shall be computed
#' @param data the dataset against which the spatial depth of `x` is computed
#'     (the reference cloud of points). Defaults to `x` itself.
#' @param ... further arguments to norm determining computation or type of norm
#' @param norm the norm to be used
#'
#' @family caseStatistics
#' @export
#' @returns the spatial depth of the values in x
#'
#' @examples
#' data(DTD)
#' Xpdf    = BDataPDFs(DTD$shapePDF,basis=B2Unif01Basis5)
#' \donttest{ spatialDepth(Xpdf) }
#'
spatialDepth = function(x,data=x,...,norm) UseMethod("spatialDepth",x)
#' @describeIn spatialDepth method for BData objects
#' @export
#'
spatialDepth.BData = function(x,data=x,...,norm=DistributionalDDA::norm) {
  nX = numberOfCases(x)
  nD = numberOfCases(data)
  iD = rep(1:nD,nX)
  iX = rep(1:nX,each=nD)
  Delta = x[iX]-data[iD]
  normDelta = norm(Delta,...)
  # Norm 0 vectors get 0:
  unitvectors = Delta/ifelse(normDelta==0,1,normDelta)
  sapply(1:nX,\(i) 1-norm(mean(unitvectors[nD*(i-1)+(1:nD)]),...))
}
#' @describeIn spatialDepth method for numerical data
#' @export
spatialDepth.numeric = function(x,data=x,...,norm=abs ) {
  nX = numberOfCases(x)
  nD = numberOfCases(data)
  iD = rep(1:nD,nX)
  iX = rep(1:nX,each=nD)
  Delta = x[iX]-data[iD]
  normDelta = norm(Delta,...)
  # Norm 0 vectors get 0:
  unitvectors = Delta/ifelse(normDelta==0,1,normDelta)
  sapply(1:nX,\(i) 1-norm(mean(unitvectors[nD*(i-1)+(1:nD)]),...))
}
#' @describeIn spatialDepth method for multivariate data
#' @export
spatialDepth.matrix = function(x,data=x,...,norm=\(x) sqrt(rowMeans(x^2)) ) {
  nX = numberOfCases(x)
  nD = numberOfCases(data)
  iD = rep(1:nD,nX)
  iX = rep(1:nX,each=nD)
  Delta = x[iX,]-data[iD,]
  normDelta = norm(Delta,...)
  # Norm 0 vectors get 0:
  unitvectors = Delta/ifelse(normDelta==0,1,normDelta)
  sapply(1:nX,\(i) 1-norm(mean(unitvectors[nD*(i-1)+(1:nD),]),...))
}



#' centralityQuantile
#'
#' Computes quantile of distanceToCenter of data in a dataset, 0 is the median, 1 is the farthest point
#'
#' @param x the dataset
#' @param ... further parameters to [distanceToCenter]
#' @param d the distances to median
#' @param distType "meandistance" for mean distance to datapoints, "median"
#'     for distance to median, "spatialdepth" for the spatial depth.
#' @param ties.method see [rank]: the method for breaking ties in the ranks
#' @returns A numeric vector in [0,1], 0 for the most central point and 1 for
#'     the most peripheral.
#' @family caseStatistics
#' @export
#' @examples
#' x = rnorm(100)
#' plot(x,centralityQuantile(x))
#' plot(x,centralityRank(x))
#'
centralityQuantile = function(x,...) UseMethod("centralityQuantile",x)
#' @describeIn centralityQuantile default method
#' @export
centralityQuantile.default = function(x,...,
                                      d=distanceToCenter(x,...,distType=distType),
                                      distType=c("meandistance","median","spatialdepth"),
                                      ties.method="first")
  (rank(d,ties.method=ties.method)-1)/(length(d)-1)
#' @describeIn centralityQuantile returns the rank rather than the quantile:
#'     an integer vector of ranks, 1 for the most central point up to the
#'     number of observations for the most peripheral.
#' @export
centralityRank = function(x,...) UseMethod("centralityRank",x)
#' @describeIn centralityQuantile default method returning ranks
#' @export
centralityRank.default = function(x,...,d=distanceToCenter(x,...,distType=distType),
                                  distType=c("meandistance","median","spatialdepth"),
                                  ties.method="first")
  rank(d,ties.method=ties.method)



#' Compute Maximum and Minimum Functions
#'
#' @param ... arbitrary list of vectorized functions or vectorized functions, behaving returning numerical values of the same size as a result to the same arguments
#' @param na.rm boolean, should na results be discared?
#' @param bound a limiting function. All values exceeding (or smaller than for minFunc) are ignored in the maximum.
#' @param Funcs the list of all functions that should be used in the maximum or minimum
#'
#' @returns a function taking like the given functions, and returning the maximum (or minium) of the values not exceeding the boundary.
#' @export
#' @examples
#' F = maxFunc(\(x) dnorm(x),\(x) dnorm(x,mean=2),bound=\(x) rep(0.3,length(x)))
#' x=seq(-2,5,by=0.01)
#' plot(x,F(x))
#' F = minFunc(\(x) dnorm(x),\(x) dnorm(x,mean=2),bound=\(x) rep(0.1,length(x)))
#' plot(x,F(x))
maxFunc <- function(...,na.rm=FALSE, bound=NULL,Funcs=c(...)) {
  Funcs
  na.rm
  .=NULL
  if( is.null(bound) ) {
    function(...) do.call("pmax",c(Funcs %.% .(...),na.rm=na.rm))
  } else {
    function(...) {
      b = bound(...)
      suppress = function(x) ifelse(x<=b,x,-Inf)
      do.call("pmax",c(Funcs %.% suppress(.(...)),na.rm=na.rm))
    }
  }
}
#' @describeIn maxFunc computes the minimum
#' @export
minFunc <- function(...,na.rm=FALSE, bound=NULL,Funcs=c(...)) {
  Funcs
  na.rm
  .=NULL
  if( is.null(bound) ) {
    function(...) do.call("pmin",c(Funcs %.% .(...),na.rm=na.rm))
  } else {
    function(...) {
      b = bound(...)
      suppress = function(x) ifelse(x>=b,x,Inf)
      do.call("pmin",c(Funcs %.% suppress(.(...)),na.rm=na.rm))
    }
  }
}



## caseQuantiles
#' Quantiles of distributions in a data set of distributions
#'
#' Ideally this should always be called on the original samples.
#'
#' @param X the dataset, typically a BData object
#' @param ... further arguments to called functions (e.g. quantile)
#'
#' @returns a vector or a matrix (with different probabilities in colums) of quantiles of the distributions in the dataset.
#' @family caseStatistics
#' @export
#' @examples
#' data(DTD)
#' X = BDataSamples(DTD$data %.% shape,basis=B2Unif01Basis5)
#' caseQuantiles(X,probs=c(0.25,0.75))
#' caseQuantiles(X,probs=0.5)
#'
caseQuantiles <- function(X,...) UseMethod("caseQuantiles",X)
#' @describeIn caseQuantiles method for BDataSamples, simply calles [stats::quantile]
#' @export
#'
caseQuantiles.BDataSamples <- function(X,...) {
  erg = sapply(X,quantile,...)
  if( length(dim(erg)) == 2 ) {
    t(erg)
  } else if( numberOfCases(X)==1 && length(erg)>1 ) {
    t(c(erg))
  } else {
    c(erg)
  }
}

## TBD: Implement caseQuantiles for BData
#' @describeIn caseQuantiles method for BData not yet implemented
#' @export
#'
caseQuantiles.BData <- function(X,...) stop("caseQuantiles not yet implemented for class",class(X),". Try applying it directly to the original sample.")

## caseMeans
#' Expected values (means) of real-valued distributions
#'
#' Ideally this should always be called on the original samples.
#'
#'
#' @param x the dataset of distributions,
#' @param ... further arguments to called functions (e.g. caseE)
#'
#' @returns a vector or a matrix (with different cases in columns) of the means
#'     of the distributions in the dataset. For multivariate (vector-valued)
#'     means, the last dimension iterates over the elements of x.
#'
#' @family caseStatistics
#' @export
#'
#' @examples
#' data(DTD)
#' X = BDataSamples(DTD$data %.% shape,basis=B2Unif01Basis5)
#' caseMeans(X)
#' X = BDataSamples(DTD$data,basis=B2Unif01Basis5)
#' caseMeans(X)
#'
caseMeans <- function(x,...) caseE(x,...,X=function(x) x)

## caseVars
#' Variances of real-valued distributions
#'
#' Ideally this should always be called on the original samples.
#'
#'
#' @param x the dataset of distributions,
#' @param ... further arguments to called functions (e.g. caseE)
#' @param X a function representing the real or vectorial random variable on the distributions in x
#'
#' @returns A vector or tensor of variances. In case of multiple
#'     dimensional means, the last dimension iterates over the
#'     elements of x.
#'
#' @family caseStatistics
#' @export
#'
#' @examples
#' data(DTD)
#' X = BDataSamples(DTD$data %.% shape,basis=B2Unif01Basis5)
#' caseVars(X)
#' X = BDataSamples(DTD$data,basis=B2Unif01Basis5)
#' caseVars(X)
#'
caseVars <- function(x,...,X=\(x) x) UseMethod("caseVars",x)
#' @describeIn caseVars method for samples
#' @export
#'
caseVars.BDataSamples <- function(x,...,X=\(x) x) {
  vars = sapply(x,var,simplify="array")
  if(length(dim(vars))==3) {
    reportIfNot(dim(vars)[3]==numberOfCases(x),dimVars=dim(vars),lengthVars=length(vars),numberOfCases=numberOfCases(x))
    c(vars)
    vars
  } else if(length(dim(vars))<=1) {
    reportIfNot(length(vars)==numberOfCases(x),dimVars=dim(vars),lengthVars=length(vars),numberOfCases=numberOfCases(x))
    c(vars)
  } else {
    reportIfNot(FALSE,list(dimVars=dim(vars)))
  }
}
#' @describeIn caseVars method for distributions
#' @export
#'
caseVars.BData <- function(x,...,X=\(x) x) {
  x = as.BDataDensities(x)
  ref = getReferenceMeasure(getBayesSpaceBasis(x))
  E = sapply(x,simplify="array",\(f) {
    Ai = expectation(ref,f=\(x) f(x),...)
    M  = expectation(ref,f=\(x) X(x)*f(x),...)/Ai
    if( length(M)== 1 ) {
      X2 = \(x) (X(x)-M)^2
      V = expectation(ref,f=\(x) X2(x)*f(x),...)/Ai
    } else {
      NN=length(M)*length(M)
      Xc = \(x) t(t(X(x))-M)
      XXtc = \(x) apply(Xc(x),length(M),function(X) outer(X,X))
      V  = expectation(ref,XXtc(x)*rep(f(x),each=NN),...)
    }
    V
  }
  )
  D = dim(E)
  if( length(D)>1 ) {
    if( numberOfCases(x)==1 )
      D=c(D,1)
    dim(E)=D
    reportIfNot(D[length(D)]==numberOfCases(x),D=D,n=numberOfCases(x))
    E
  } else {
    if( numberOfCases(x)>1 ) {
      reportIfNot(length(E)==numberOfCases(x),lengthE=length(E),n=numberOfCases(x))
      c(E)
    } else {
      if( length(E)>1 ) {
        dim(E)=c(length(E),1)
        E
      } else {
        c(E)
      }
    }
  }
}



## caseSds
#' Standard deviations of real-valued distributions
#'
#' Ideally this should always be called on the original samples.
#'
#'
#' @param x the dataset of distributions,
#' @param ... further arguments to called functions (e.g. caseE)
#'
#' @returns a vector the standard deviations of the distributions in
#'     the dataset.
#'
#' @family caseStatistics
#' @export
#'
#' @examples
#' data(DTD)
#' X = BDataSamples(DTD$data %.% shape,basis=B2Unif01Basis5)
#' caseSds(X)
#'
caseSds <- function(x,...) sapply(caseVars(x,...),sqrt)



## caseE
#' Expected values of functions of distributions
#'
#' Ideally this should always be called on the original samples.
#'
#'
#' @param x the dataset of distributions,
#' @param X the real random variable, i.e. a real function define on whats random in the distributions
#' @param ... further arguments to called functions (e.g. quantile)
#'
#' @returns a vector, matrix or with the Expected values for each
#'     element of x.  if the expected value is multivariate the last
#'     dimension iterates over the cases in x.
#' @family caseStatistics
#' @export
#'
#' @examples
#' data(DTD)
#' X = BDataSamples(DTD$data %.% shape,basis=B2Unif01Basis5)
#' caseE(X,X=function(x) x,probs=c(0.25,0.75))
#' caseE(X,X=function(x) x^2,probs=0.5)
#'
caseE <- function(x,X,...) UseMethod("caseE",x)
#' @describeIn caseE method for BDataSamples simply uses means
#' @export
#'
caseE.BDataSamples <- function(x,X,...) {
  E = sapply(x,simplify="array",function(x) {
    Y = X(x)
    if(numberOfCases(x)<2)
      return(Y)
    if( length(dim(Y))==2 ) {
      return(colMeans(Y))
    } else {
      return(mean(Y))
    }
  })
  D = dim(E)
  if( length(D)>1 ) {
    if( numberOfCases(x)==1 )
      D=c(D,1)
    dim(E)=D
    reportIfNot(D[length(D)]==numberOfCases(x),D=D,n=numberOfCases(x))
    E
  } else {
    if( numberOfCases(x)>1 ) {
      reportIfNot(length(E)==numberOfCases(x),lengthE=length(E),n=numberOfCases(x))
      c(E)
    } else {
      if( length(E)>1 ) {
        dim(E)=c(length(E),1)
        E
      } else {
        c(E)
      }
    }
  }
}

#' @describeIn caseE method for distributions
#' @export
#'
caseE.BData <- function(x,X,...) {
  x = as.BDataDensities(x)
  ref = getReferenceMeasure(getBayesSpaceBasis(x))
  E = sapply(x,simplify="array",\(f) expectation(ref,f=\(x) X(x)*f(x),...)/
               expectation(ref,f=\(x) f(x),...))
  D = dim(E)
  if( length(D)>1 ) {
    if( numberOfCases(x)==1 )
      D=c(D,1)
    dim(E)=D
    reportIfNot(D[length(D)]==numberOfCases(x),D=D,n=numberOfCases(x))
    E
  } else {
    if( numberOfCases(x)>1 ) {
      reportIfNot(length(E)==numberOfCases(x),lengthE=length(E),n=numberOfCases(x))
      c(E)
    } else {
      if( length(E)>1 ) {
        dim(E)=c(length(E),1)
        E
      } else {
        c(E)
      }
    }
  }
}



## caseSim


## caseLogRatios
#' Log-density ratios for each distribution (case)
#'
#' Density log ratios independent of the reference measure and any other values.
#'
#' @param x a BData dataset of distributions
#' @param numerator numerator, i.e. the same number of possible values of the random varianble
#' @param denominator denominator values, i.e. possible values of the random variable
#' @param ... further unused arguments for future generic use
#' @returns a dataset of the same number of cases as x of the corresponding log rations
#'
#' @family caseStatistics
#' @export
#'
#' @examples
#' data(DTD)
#' X  = BDataPDFs(DTD$shapePDF,basis=B2Unif01Basis5)
#' caseLogRatios(X,0.8,0.9)
#'
caseLogRatios <- function(x,numerator,denominator,...) UseMethod("caseLogRatios",x)
#' @export
#'
caseLogRatios.BDataDensities <- function(x,numerator,denominator,...) {
  stopifnot(length(numerator)==length(denominator))
  . = NULL # avoid warning
  E = log(x %..% .(numerator) / x %..% .(denominator))
  D = dim(E)
  if( length(D)>1 ) {
    if( numberOfCases(x)==1 )
      D=c(D,1)
    dim(E)=D
    reportIfNot(D[length(D)]==numberOfCases(x),D=D,n=numberOfCases(x))
    E
  } else {
    if( numberOfCases(x)>1 ) {
      reportIfNot(length(E)==numberOfCases(x),lengthE=length(E),n=numberOfCases(x))
      c(E)
    } else {
      if( length(E)>1 ) {
        dim(E)=c(length(E),1)
        E
      } else {
        c(E)
      }
    }
  }
}
#' @export
#'
caseLogRatios.BDataPDFs = caseLogRatios.BDataDensities
#' @export
#'
caseLogRatios.BDataPDFs = caseLogRatios.BDataDensities
#' @export
caseLogRatios.BDataLogDensities <- function(x,numerator,denominator,...) {
  stopifnot(length(numerator)==length(denominator))
  . = NULL # avoid warning
  E = x %..% .(numerator) - x %..% .(denominator)
  D = dim(E)
  if( length(D)>1 ) {
    if( numberOfCases(x)==1 )
      D=c(D,1)
    dim(E)=D
    reportIfNot(D[length(D)]==numberOfCases(x),D=D,n=numberOfCases(x))
    E
  } else {
    if( numberOfCases(x)>1 ) {
      reportIfNot(length(E)==numberOfCases(x),lengthE=length(E),n=numberOfCases(x))
      c(E)
    } else {
      if( length(E)>1 ) {
        dim(E)=c(length(E),1)
        E
      } else {
        c(E)
      }
    }
  }
}
#' @export
caseLogRatios.BDataLogPDFs <- caseLogRatios.BDataLogDensities
#' @export
caseLogRatios.BDataSplines <- function(x,numerator,denominator,...) {
  caseLogRatios(as.BDataLogDensities(x),numerator=numerator,denominator=denominator,...)
}