utilities.R 5.3 KB
Newer Older
davidkep's avatar
davidkep committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160
#' Compute the Tau-Scale of Centered Values
#'
#' Compute the tau-scale without centering the values.
#'
#' @param x numeric values.
#' @return the tau scale.
#' @export
tau_size <- function (x) {
  # No checks for NA values!
  .Call(C_tau_size, as.numeric(x))
}

#' Compute the M-Scale of Centered Values
#'
#' Compute the M-scale without centering the values.
#'
#' @param x numeric values.
#' @param bdp desired breakdown point (between 0 and 0.5).
#' @param cc cutoff value for the bisquare rho function. By default, chosen
#'           for a consistent estimate under the Normal model.
#' @param opts a list of options for the M-scale equation, see [mest_options]
#'             for details.
#' @return the m-scale estimate.
#' @export
mscale <- function (x, bdp = 0.25, cc = consistency_const(bdp, 'bisquare'),
                    opts = mest_options()) {
  # No checks for NA values!
  opts$delta <- as.numeric(bdp[[1L]])
  opts$cc <- as.numeric(cc[[1L]])
  .Call(C_mscale, as.numeric(x), opts)
}

#' Compute the M-estimate of Location
#'
#' @param x numeric values.
#' @param scale scale of the `x` values.
#' @param rho the rho function to use.
#' @param cc value of the tuning constant for the chosen rho function.
#'           By default, chosen to achieve 95% efficiency under the Normal
#'           model.
#' @param opts a list of options for the M-estimating equation, see
#'             [mest_options] for details.
#' @return the m-scale estimate.
#' @importFrom stats mad
#' @export
mloc <- function (x, scale = mad(x), rho, cc, opts = mest_options()) {
  # No checks for NA values!
  opts$rho <- rho_function(rho)
  if (!missing(cc)) {
    opts$cc <- as.numeric(cc[[1L]])
  }
  .Call(C_mloc, as.numeric(x), scale, opts)
}

#' Compute the M-Location and M-Scale
#'
#' Simultaneously estiamte the M-Location and the M-Scale.
#'
#' @param x numeric values.
#' @param bdp desired breakdown point (between 0 and 0.5).
#' @param cc cutoff value for the bisquare rho function. By default, chosen
#'           for a consistent estimate under the Normal model.
#' @param opts a list of options for the M-estimating equation,
#'             see [mest_options] for details.
#' @return a vector with two elements, the M-location and the M-scale estimate.
#' @export
mlocscale <- function (x, bdp = 0.25, location_rho = c('bisquare', 'huber'),
                       scale_cc = consistency_const(bdp, 'bisquare'),
                       location_cc, opts = mest_options()) {
  # No checks for NA values!
  opts$delta <- as.numeric(bdp[[1L]])
  opts$cc <- as.numeric(scale_cc[[1L]])
  loc_opts <- list(rho = rho_function(location_rho))
  if (!missing(location_cc)) {
    loc_opts$cc <- as.numeric(location_cc[[1L]])
  }
  .Call(C_mlocscale, as.numeric(x), opts, loc_opts)
}

#' Get the Constant for Consistency for the M-Scale Using the Bisquare Rho Function
#' @param delta desired breakdown point (between 0 and 0.5)
#'
#' @return consistency constant
.bisquare_consistency_const <- function (delta) {
  ##
  ## Pre-computed values for some delta values
  ##
  eps <- sqrt(.Machine$double.eps)
  if (!isTRUE(delta < 0.5 + eps && delta > -eps)) {
    stop("`delta` is outside valid bounds")
  }

  if (abs(delta - 0.5) < eps) {
    return(1.5476450)
  } else if (abs(delta - 0.25) < eps) {
    return(2.937015)
  } else if (abs(delta - 0.1) < eps) {
    return(5.182361)
  } else if (delta < 0.005) {
    return(50) # ~.1% bdp for bisquare
  }
  integral_interval <- if (delta > 0.1) {
    c(1.5, 5.5)
  } else {
    c(5, 25)
  }

  # For bisquare we have the closed form solution to the expectation
  expectation <- function(cc, delta) {
    pnorm.mcc <- 2 * pnorm(-cc)
    1/cc^6 * exp(-(cc^2/2)) * (
      -cc * (15 - 4 * cc^2 + cc^4) * sqrt(2 / pi) +
          3 * (5 - 3 * cc^2 + cc^4) * exp(cc^2/2) * (1 - pnorm.mcc) +
          cc^6 * exp(cc^2/2) * pnorm.mcc
    ) - delta
  }
  uniroot(expectation, interval = integral_interval, delta)$root
}

#' Get the Constant for Consistency for the M-Scale
#'
#' @param delta desired breakdown point (between 0 and 0.5)
#' @param rho the name of the chosen rho function.
#'
#' @return consistency constant
#' @export
consistency_const <- function (delta, rho) {
  return(switch(rho_function(rho),
                bisquare = .bisquare_consistency_const(delta),
                huber = stop("Huber's rho function not supported for scale ",
                             "estimation!")))
}

#' List or check available rho functions.
#'
#' @param rho the name of the rho function to check.
#' @return if `rho` is missing returns a vector of rho function names, otherwise
#'         the integer representation of the rho function.
#' @export
rho_function <- function (rho) {
  available <- c('bisquare', 'huber')
  if (missing(rho)) {
    return(available)
  }
  return(match(match.arg(rho, available), available))
}

#' Approximate Value Matching
#'
#' @param x,table see [base::match] for details.
#' @param eps numerical tolerance for matching.
#' @return a vector the same lenght of `x` with integers giving the position in
#'         `table` of the first match if there is a match, or `NA_integer_`
#'         otherwise.
.approx_match <- function(x, table,
                          eps = min(sqrt(.Machine$double.eps),
                                    0.5 * min(x, table))) {
  .Call('C_approx_match', as.numeric(x), as.numeric(table),
        as.numeric(eps[[1L]]))
}