R/hermite_estimator_univar.R

Defines functions kendall.hermite_estimator_univar spearmans.hermite_estimator_univar summary.hermite_estimator_univar print.hermite_estimator_univar quant.hermite_estimator_univar quantile_helper_interpolate quantile_helper_bisection cum_prob.hermite_estimator_univar dens.hermite_estimator_univar calculate_running_std.hermite_estimator_univar update_batch.hermite_estimator_univar update_sequential.hermite_estimator_univar merge_hermite_univar merge_pair.hermite_estimator_univar merge_standardized_helper_univar merge_moments_and_count_univar hermite_estimator_univar

Documented in cum_prob.hermite_estimator_univar dens.hermite_estimator_univar hermite_estimator_univar merge_hermite_univar merge_moments_and_count_univar merge_pair.hermite_estimator_univar merge_standardized_helper_univar print.hermite_estimator_univar quant.hermite_estimator_univar summary.hermite_estimator_univar update_batch.hermite_estimator_univar update_sequential.hermite_estimator_univar

#' A class to sequentially estimate univariate pdfs, cdfs and quantile functions
#'
#' This method constructs an S3 object with associated methods for univariate
#' nonparametric estimation of pdfs, cdfs and quantiles.
#'
#' The hermite_estimator_univar class allows the sequential or one-pass batch
#' estimation of the full probability density function, cumulative distribution
#' function and quantile function. It is well suited to streaming data (both
#' stationary and non-stationary) and to efficient estimation in the context of
#' massive or distributed data sets. Indeed, estimators constructed on different
#' subsets of a distributed data set can be consistently merged.
#'
#' @author Michael Stephanou <michael.stephanou@gmail.com>
#'
#' @param N An integer between 0 and 75. The Hermite series based estimator
#' is truncated at N+1 terms.
#' @param standardize A boolean value. Determines whether the observations are
#' standardized, a transformation which often improves performance.
#' @param exp_weight_lambda A numerical value between 0 and 1. This parameter
#' controls the exponential weighting of the Hermite series based estimator.
#' If this parameter is NA, no exponential weighting is applied.
#' @return An S3 object of class hermite_estimator_univar, with methods for
#' density function, distribution function and quantile function estimation.
#' @export
#' @examples
#' hermite_est <- hermite_estimator_univar(N = 30, standardize = TRUE)
hermite_estimator_univar <-
  function(N = 30,
           standardize = TRUE,
           exp_weight_lambda = NA) {
    if (!is.numeric(N)) {
      stop("N must be numeric.")
    }
    if (N < 0 | N > 75) {
      stop("N must be >= 0 and N <= 75.")
    }
    if (!(standardize == TRUE | standardize == FALSE)) {
      stop("standardize can only take on values TRUE or FALSE.")
    }
    if (!is.na(exp_weight_lambda)) {
      if (!is.numeric(exp_weight_lambda)) {
        stop("exp_weight_lambda must be numeric.")
      }
      if (exp_weight_lambda <= 0 | exp_weight_lambda > 1) {
        stop("exp_weight_lambda must be a real number > 0 and <= 1.")
      }
    }
    this <-
      list(
        N_param = N,
        coeff_vec = rep(0, N + 1),
        num_obs = 0,
        standardize_obs = standardize,
        running_mean = 0,
        running_variance = 0,
        exp_weight = exp_weight_lambda,
        normalization_hermite_vec = c(),
        h_int_lower_serialized = c(),
        h_int_upper_serialized = c()
      )
    this$normalization_hermite_vec <- h_norm_serialized[1:(this$N_param+1)]
    this$h_int_lower_serialized <- h_int_lower_serialized[1:(this$N_param+1),
                                                          ,drop = FALSE]
    this$h_int_upper_serialized <- h_int_upper_serialized[1:(this$N_param+1),
                                                          ,drop = FALSE]
    class(this) <- c("hermite_estimator_univar", "list")
    return(this)
  }

#' Internal method to consistently merge the number of observations, means and 
#' variances of two Hermite estimators
#' 
#' The algorithm to merge the variances consistently comes from
#' Schubert, Erich, and Michael Gertz. "Numerically stable parallel computation 
#' of (co-) variance." Proceedings of the 30th International Conference on 
#' Scientific and Statistical Database Management. 2018.
#'
#' @param hermite_estimator1 A hermite_estimator_univar object.
#' @param hermite_estimator2 A hermite_estimator_univar object.
#' @return An object of class hermite_estimator_univar.
merge_moments_and_count_univar <- function(hermite_estimator1, 
                                             hermite_estimator2){
  num_obs_1 <- hermite_estimator1$num_obs
  num_obs_2 <- hermite_estimator2$num_obs
  hermite_merged <- hermite_estimator_univar(hermite_estimator1$N_param, 
                                        hermite_estimator1$standardize_obs)
  hermite_merged$num_obs <- num_obs_1 + num_obs_2
  hermite_merged$running_mean <- (num_obs_1*hermite_estimator1$running_mean +
                num_obs_2*hermite_estimator2$running_mean)/(num_obs_1+num_obs_2)
  hermite_merged$running_variance <- (hermite_estimator1$running_variance +
      hermite_estimator2$running_variance) + ((num_obs_1 * num_obs_2) / 
     (num_obs_1 + num_obs_2)) *(hermite_estimator1$running_mean 
                                - hermite_estimator2$running_mean)^2
  return(hermite_merged)
}

#' Internal method to merge a list of standardized Hermite estimators
#'
#'
#' @param hermite_estimators A list of hermite_estimator_univar objects.
#' @return An object of class hermite_estimator_univar.
merge_standardized_helper_univar <- function(hermite_estimators) {
  all_N <- lapply(hermite_estimators, FUN =
                          function(x){return(x$N_param)})
  if (length(unique(all_N)) >1) {
    stop("List must contain Hermite estimators with a consistent N")
  }
  N <- hermite_estimators[[1]]$N_param
  hermite_estimator_merged <- base::Reduce(f=merge_moments_and_count_univar, 
                                             x = hermite_estimators)
  hermite_estimator_merged$coeff_vec <-
    vapply(1:(N+1),FUN=function(k){sum(vapply(hermite_estimators, 
          FUN=function(x){(x$num_obs / hermite_estimator_merged$num_obs) *
      gauss_hermite_quad_100(function(t){integrand_coeff_univar(t, x, 
                           hermite_estimator_merged, k)})},
      FUN.VALUE=numeric(1)))}, FUN.VALUE=numeric(1))
  return(hermite_estimator_merged)
}

#' Merges two Hermite estimators
#'
#' This method allows a pair of Hermite based estimators of class
#' hermite_estimator_univar to be consistently merged.
#'
#' Note that the N and standardize arguments must be the same for the two
#' estimators in order to merge them. In addition, note that exponentially
#' weighted estimators cannot be merged. If the Hermite estimators are not
#' standardized, the merged estimator will be exactly equivalent to
#' constructing a single estimator on the data set formed by combining the
#' data sets used to update the respective hermite_estimator_univar inputs.
#' If the input Hermite estimators are standardized however, then the
#' equivalence will be approximate but still accurate in most cases.
#'
#' @param this A hermite_estimator_univar object. The first Hermite series based
#' estimator.
#' @param hermite_estimator_other A hermite_estimator_univar object. The 
#' second Hermite series based estimator.
#' @return An object of class hermite_estimator_univar.
#' @export
#' @examples
#' hermite_est_1 <- hermite_estimator_univar(N = 10, standardize = FALSE)
#' hermite_est_1 <- update_batch(hermite_est_1, rnorm(30))
#' hermite_est_2 <- hermite_estimator_univar(N = 10, standardize = FALSE)
#' hermite_est_2 <- update_batch(hermite_est_2, rnorm(30))
#' hermite_merged <- merge_pair(hermite_est_1, hermite_est_2)
merge_pair.hermite_estimator_univar <-
  function(this, hermite_estimator_other) {
    if (!is(hermite_estimator_other, "hermite_estimator_univar")) {
      stop("merge_pair.hermite_estimator_univar can only be applied to 
           hermite_estimator_univar objects.")
    }
    if (this$N_param != hermite_estimator_other$N_param) {
      stop("N must be equal to merge estimators.")
    }
    if (this$standardize_obs != hermite_estimator_other$standardize_obs) {
      stop("Standardization setting must be the same to merge estimators.")
    }
    if (!is.na(this$exp_weight) |
        !is.na(hermite_estimator_other$exp_weight)) {
      stop("Cannot merge exponentially weighted estimators.")
    }
    if (this$standardize_obs == FALSE) {
      hermite_estimator_merged <- merge_moments_and_count_univar(this, 
                                                      hermite_estimator_other)
      hermite_estimator_merged$coeff_vec <-
        (
          this$coeff_vec * this$num_obs + hermite_estimator_other$coeff_vec 
          * hermite_estimator_other$num_obs
        ) / (this$num_obs + hermite_estimator_other$num_obs)
    } else {
      hermite_estimator_merged <- 
        merge_standardized_helper_univar(list(this,hermite_estimator_other))
    }
    return(hermite_estimator_merged)
  }

#' Merges a list of Hermite estimators
#'
#' This method allows a list of Hermite based estimators of class
#' hermite_estimator_univar to be consistently merged.
#'
#' Note that the N and standardize arguments must be the same for all estimators
#' in order to merge them. In addition, note that exponentially weighted
#' estimators cannot be merged. If the Hermite estimators are not
#' standardized, the merged estimator will be exactly equivalent to
#' constructing a single estimator on the data set formed by combining the
#' data sets used to update the respective hermite_estimator_univar inputs.
#' If the input Hermite estimators are standardized however, then the
#' equivalence will be approximate but still accurate in most cases.
#'
#' @param hermite_estimators A list of hermite_estimator_univar objects.
#' @return An object of class hermite_estimator_univar.
#' @export
#' @examples
#' hermite_est_1 <- hermite_estimator_univar(N = 10, standardize = FALSE)
#' hermite_est_1 <- update_batch(hermite_est_1, rnorm(30))
#' hermite_est_2 <- hermite_estimator_univar(N = 10, standardize = FALSE)
#' hermite_est_2 <- update_batch(hermite_est_2, rnorm(30))
#' hermite_merged <- merge_hermite(list(hermite_est_1, hermite_est_2))
merge_hermite_univar <- function(hermite_estimators) {
  if (length(hermite_estimators) == 1) {
    return(hermite_estimators[[1]])
  }
  if (hermite_estimators[[1]]$standardize_obs==FALSE){
   hermite_estimator_merged <- base::Reduce(merge_pair,hermite_estimators)
  } else {
   hermite_estimator_merged <- 
     merge_standardized_helper_univar(hermite_estimators)
  }
  return(hermite_estimator_merged)
}

#' Updates the Hermite series based estimator sequentially
#'
#' This method can be applied in sequential estimation settings.
#' 
#'
#' @param this A hermite_estimator_univar object.
#' @param x A numeric value. An observation to be incorporated into the
#' estimator.
#' @return An object of class hermite_estimator_univar.
#' @export
#' @examples
#' hermite_est <- hermite_estimator_univar(N = 10, standardize = TRUE)
#' hermite_est <- update_sequential(hermite_est, x = 2)
update_sequential.hermite_estimator_univar <- function(this, x) {
  if (!is.numeric(x)) {
    stop("x must be numeric.")
  }
  if (length(x) != 1) {
    stop("The sequential method is only 
         applicable to one observation at a time.")
  }
  if (is.na(x) | !is.finite(x)){
    stop("The sequential method is only 
         applicable to finite, non NaN, non NA values.")
  }
  this$num_obs <- this$num_obs + 1
  if (this$standardize_obs == TRUE) {
    if (is.na(this$exp_weight)) {
      prev_running_mean <- this$running_mean
      this$running_mean <-  (this$running_mean * (this$num_obs - 1) + x) /
        this$num_obs
      if (this$num_obs < 2) {
        return(this)
      }
      this$running_variance <- this$running_variance + (x - prev_running_mean) *
        (x - this$running_mean)
      x <- (x - this$running_mean) /  sqrt(this$running_variance /
                                            (this$num_obs - 1))
    } else {
      if (this$num_obs < 2){
        this$running_mean <- x
        this$running_variance <- 1
        return(this)
      }
      this$running_mean <-  (1 - this$exp_weight) * this$running_mean +
        this$exp_weight * x
      this$running_variance <- (1 - this$exp_weight) * this$running_variance +
        this$exp_weight * (x - this$running_mean)^2
      x <- (x - this$running_mean) / sqrt(this$running_variance)
    }
  }
  h_k <-
    as.vector(hermite_function_N(this$N_param, x, 
                                 this$normalization_hermite_vec))
  if (is.na(this$exp_weight)) {
    this$coeff_vec <-
      (this$coeff_vec * (this$num_obs - 1) + h_k) / this$num_obs
  } else {
    this$coeff_vec <-
      this$coeff_vec * (1 - this$exp_weight) + h_k * this$exp_weight
  }
  return(this)
}

#' Updates the Hermite series based estimator with a batch of data
#'
#' This method can be applied in one-pass batch estimation settings. This
#' method cannot be used with an exponentially weighted estimator.
#'
#' @param this A hermite_estimator_univar object.
#' @param x A numeric vector. A vector of observations to be incorporated
#' into the estimator.
#' @return An object of class hermite_estimator_univar.
#' @export
#' @examples
#' hermite_est <- hermite_estimator_univar(N = 10, standardize = TRUE)
#' hermite_est <- update_batch(hermite_est, x = c(1, 2))
update_batch.hermite_estimator_univar <- function(this, x) {
  if (!is.numeric(x)) {
    stop("x must be numeric.")
  }
  if (!is.na(this$exp_weight)) {
    stop("The Hermite estimator cannot be exponentially weighted.")
  }
  if (length(x) < 1) {
    stop("x must contain at least one value.")
  }
  if (any(is.na(x)) | any(!is.finite(x))){
    stop("The batch update method is only 
         applicable to finite, non NaN, non NA values.")
  }
  this$num_obs <- length(x)
  if (this$standardize_obs == TRUE) {
    this$running_mean <- mean(x)
    this$running_variance <- stats::var(x) * (length(x) - 1)
    x <-
      (x - this$running_mean) / sqrt(this$running_variance / (this$num_obs - 1))
  }
  this$coeff_vec <- hermite_function_sum_N(this$N_param, x,
                                 this$normalization_hermite_vec) / this$num_obs
  return(this)
}

calculate_running_std.hermite_estimator_univar<- function(this){
  if (is.na(this$exp_weight)) {
    running_std <- sqrt(this$running_variance / (this$num_obs - 1))
  } else {
    running_std <- sqrt(this$running_variance)
  }
  return(running_std)
}

#' Estimates the probability density for a vector of x values
#'
#' This method calculates the probability density values at a vector of
#' x values using the hermite_estimator_univar object (this).
#'
#' The object must be updated with observations prior to the use of the method.
#'
#' @param this A hermite_estimator_univar object.
#' @param x A numeric vector. Values at which to estimate the probability
#' density.
#' @param clipped A boolean value. This value determines whether
#' probability densities are clipped to be bigger than zero.
#' @param accelerate_series A boolean value. This value determines whether
#' Hermite series acceleration is applied.
#' @return A numeric vector of probability density values.
#' @export
#' @examples
#' hermite_est <- hermite_estimator_univar(N = 10, standardize = TRUE)
#' hermite_est <- update_batch(hermite_est, rnorm(30))
#' pdf_est <- dens(hermite_est, c(0, 0.5, 1))
dens.hermite_estimator_univar <- function(this, x, clipped = FALSE, 
                                          accelerate_series = TRUE) {
  if (!is.numeric(x)) {
    stop("x must be numeric.")
  }
  if (length(x) < 1) {
    stop("x must contain at least one value.")
  }
  if (this$num_obs < 2) {
    return(rep(NA, length(x)))
  }
  factor <- 1
  if (this$standardize_obs == TRUE) {
    running_std <- calculate_running_std(this)
    x <- (x - this$running_mean) / running_std
    factor <- 1 / running_std
  }
  h_k <-
    hermite_function_N(this$N_param, x, this$normalization_hermite_vec)
  # pdf_val <- crossprod(h_k, this$coeff_vec) * factor
  pdf_val <- series_calculate(h_k, this$coeff_vec, accelerate_series = 
                                accelerate_series) * factor
  if (clipped == TRUE) {
    pdf_val <- pmax(pdf_val, 1e-08)
  }
  return(as.vector(pdf_val))
}

#' Estimates the cumulative probability for a vector of x values
#'
#' This method calculates the cumulative probability values at a vector of
#' x values using the hermite_estimator_univar object (this).
#'
#' The object must be updated with observations prior to the use of this method.
#'
#' @param this A hermite_estimator_univar object.
#' @param x A numeric vector. Values at which to estimate the cumulative
#' probability
#' @param clipped A boolean value. This value determines whether cumulative
#' probabilities are clipped to lie within the range [0,1].
#' @param accelerate_series A boolean value. This value determines whether
#' Hermite series acceleration is applied.
#' @return A numeric vector of cumulative probability values.
#' @export
#' @examples
#' hermite_est <- hermite_estimator_univar(N = 10, standardize = TRUE)
#' hermite_est <- update_batch(hermite_est, rnorm(30))
#' cdf_est <- cum_prob(hermite_est, c(0, 0.5, 1))
cum_prob.hermite_estimator_univar <- function(this, x, clipped = FALSE, 
                                              accelerate_series = TRUE) {
  if (!is.numeric(x)) {
    stop("x must be numeric.")
  }
  if (length(x) < 1) {
    stop("x must contain at least one value.")
  }
  if (this$num_obs < 2) {
    return(rep(NA, length(x)))
  }
  if (this$standardize_obs == TRUE) {
    if (this$running_variance == 0) {
      return(ifelse(this$running_mean <= x,1,0))
    }
    running_std <- calculate_running_std(this)
    x <- (x - this$running_mean) / running_std
  }
  h_k <-
    hermite_function_N(this$N_param, x, this$normalization_hermite_vec)
  integrals_hermite <- hermite_int_lower(this$N_param, x, 
                                         hermite_function_matrix = h_k)
  cdf_val <- series_calculate(integrals_hermite, this$coeff_vec, 
                              accelerate_series = accelerate_series)
  if (clipped == TRUE) {
    cdf_val <- pmin(pmax(cdf_val, 1e-08), 1)
  }
  return(as.vector(cdf_val))
}

# Estimates the quantile at a vector of probability values using a vectorized
# implementation of the bisection search root finding algorithm.
#
# This helper method is intended for internal use by the 
# hermite_estimator_univar class.
quantile_helper_bisection <- function(this, p_vec, accelerate_series) {
  f_est <- function(x,p) {
    lower_idx <- which(x < x_lower)
    upper_idx <- which(x > x_upper)
    ambig_idx <- which(x >= x_lower & x <= x_upper)
    res <- rep(NA,length(x))
    if (length(lower_idx)>0){
      res[lower_idx] <- series_calculate(hermite_int_lower(this$N_param,
         x[lower_idx], normalization_hermite = this$normalization_hermite_vec), 
         this$coeff_vec, accelerate_series) - p[lower_idx]
    }
    if (length(upper_idx)>0){
      res[upper_idx] <- 1 - series_calculate(hermite_int_upper(this$N_param,
        x[upper_idx], normalization_hermite = this$normalization_hermite_vec), 
         this$coeff_vec, accelerate_series) - p[upper_idx]
    }
    if (length(ambig_idx)>0){
      if (p_upper < p_lower){
        res[ambig_idx] <- (p_upper + (x[ambig_idx]-x_lower) * 
            as.numeric((p_lower-p_upper)/(x_upper-x_lower))) - p[ambig_idx]
      }
      else if (p_upper > p_lower){
        res[ambig_idx] <- (p_lower + (x[ambig_idx]-x_lower) * 
                        (p_upper-p_lower)/(x_upper-x_lower)) - p[ambig_idx]
      } else if (p_upper == p_lower){
        res[ambig_idx] <- p_upper - p[ambig_idx]
      }
    }
    return(res)
  }
  h_int_lower_zero_serialized_mat <- 
    h_int_lower_zero_serialized[1:(this$N_param+1)]
  dim(h_int_lower_zero_serialized_mat) <- c(this$N_param+1,1)
  h_int_upper_zero_serialized_mat <- 
    h_int_upper_zero_serialized[1:(this$N_param+1)]
  dim(h_int_upper_zero_serialized_mat) <- c(this$N_param+1,1)
  p_lower <- as.numeric(series_calculate(h_int_lower_zero_serialized_mat, 
                                         this$coeff_vec, accelerate_series))
  p_upper <- 1-as.numeric(series_calculate(h_int_upper_zero_serialized_mat, 
                                           this$coeff_vec, accelerate_series))
  if (is.na(p_lower) | is.na(p_upper)){
    return(rep(NA, length(p_vec)))
  }
  if (p_upper < p_lower){
    x_lower <- tryCatch({
      stats::uniroot(
        f = function(x) {
          series_calculate(hermite_int_lower(this$N_param,x, 
               normalization_hermite = this$normalization_hermite_vec), 
               this$coeff_vec, accelerate_series) - p_upper
        },
        interval = c(-100, 100)
      )$root
    },
    error = function(e) {NA})
    x_upper <- tryCatch({
      stats::uniroot(
        f = function(x) {
          1- series_calculate(hermite_int_upper(this$N_param,x,
             normalization_hermite = this$normalization_hermite_vec), 
             this$coeff_vec, accelerate_series) - p_lower
        },
        interval = c(-100, 100)
      )$root
    },
    error = function(e) {NA})
  } else if (p_upper > p_lower) {
    x_lower <- -1e-6
    x_upper <- 1e-6
  } else if (p_upper == p_lower){
    x_lower <- 0
    x_upper <- 0
  }
  if (is.na(x_lower) | is.na(x_upper)){
    return(rep(NA, length(p_vec)))
  }
  max_steps <- 25
  eps_quant <- 2e-4
  x_0 <- rep(-50,length(p_vec))
  x_1 <- rep(50,length(p_vec))
  f_0 <- rep(0,length(p_vec)) - p_vec
  f_1 <- rep(1,length(p_vec)) - p_vec
  for (step in seq_len(max_steps)) {
    est  <- (x_0 + x_1)/2
    f_mid <- f_est(est,p_vec)
    mask_0 <- sign(f_mid) == sign(f_0)
    mask_1 <-  sign(f_mid) == sign(f_1)
    x_0 <- ifelse( mask_0, est, x_0)
    x_1 <- ifelse( mask_1, est, x_1)
    f_0 <- ifelse( mask_0, f_mid, f_0 )
    f_1 <- ifelse( mask_1, f_mid, f_1 )
    error_max <- max(abs(x_1 - x_0))
    if (error_max <= eps_quant) {break}
  }
  if (is.na(this$exp_weight)) {
    est <-
    est * sqrt(this$running_variance / (this$num_obs - 1)) + this$running_mean
  } else {
    est <- est * sqrt(this$running_variance) + this$running_mean
  }
  return(est)
}

# Estimates the quantile at a vector of probability values using an 
# interpolation approximation.
#
# This helper method is intended for internal use by the 
# hermite_estimator_univar class.
quantile_helper_interpolate <- function(this, p_vec, accelerate_series = TRUE){
  result <- rep(NA,length(p_vec))
  p_lower_vals <- series_calculate(this$h_int_lower_serialized,
                                   this$coeff_vec, accelerate_series)
  p_upper_vals <- 1 - series_calculate(this$h_int_upper_serialized,
                                       this$coeff_vec, accelerate_series)
  p_all_vals <- cummax(c(p_lower_vals,p_upper_vals))
  result <- stats::approx(p_all_vals,x_full_domain_serialized,
                   xout = p_vec,method="linear",ties = "ordered")$y
  if (is.na(this$exp_weight)) {
    result <-
      result * sqrt(this$running_variance / (this$num_obs - 1)) +
      this$running_mean
  } else {
    result <- result * sqrt(this$running_variance) + this$running_mean
  }
  return(result)
}

#' Estimates the quantiles at a vector of probability values
#' 
#' This method utilizes the estimator (13) in paper Stephanou, Michael, 
#' Varughese, Melvin and Iain Macdonald. "Sequential quantiles via Hermite 
#' series density estimation." Electronic Journal of Statistics 11.1 (2017): 
#' 570-607 <doi:10.1214/17-EJS1245>, with some modifications to improve the 
#' stability of numerical root finding. 
#'
#' @param this A hermite_estimator_univar object.
#' @param p A numeric vector. A vector of probability values.
#' @param algorithm A string. Two possible values 'interpolate' which is faster
#' but may be less accurate or 'bisection' which is slower but potentially more
#' accurate.
#' @param accelerate_series A boolean value. If set to TRUE, the series 
#' acceleration methods described in:
#'
#' Boyd, John P., and Dennis W. Moore. "Summability methods for 
#' Hermite functions." Dynamics of atmospheres and oceans 10.1 (1986): 51-62. 
#'  
#' are applied. If set to FALSE, then standard summation is applied.
#' @return A numeric vector. The vector of quantile values associated with the
#' probabilities p.
#' @export
#' @examples
#' hermite_est <- hermite_estimator_univar(N = 10, standardize = TRUE)
#' hermite_est <- update_batch(hermite_est, rnorm(30))
#' quant_est <- quant(hermite_est, c(0.25, 0.5, 0.75))
quant.hermite_estimator_univar <- function(this, p, algorithm="interpolate", 
                                           accelerate_series = TRUE) {
  if (!is.numeric(p)) {
    stop("p must be numeric.")
  }
  if (length(p) < 1) {
    stop("p must contain at least one value.")
  }
  if (any(p>1) | any(p<0)) {
    stop("p must contain probabilities i.e. p>=0 and p<=1.")
  }
  if (this$standardize_obs != TRUE) {
    stop("Quantile estimation requires standardization to be true.")
  }
  if (!(algorithm=="interpolate" | algorithm=="bisection") ) {
    stop("Algorithm must be either 'interpolate' or 'bisection'.")
  }
  result <- rep(NA, length(p))
  if (this$num_obs < 2) {
    return(result)
  }
  if (this$running_variance == 0) {
    return(rep(this$running_mean, length(p)))
  }
  if (algorithm=="interpolate"){
    result <- quantile_helper_interpolate(this,p,accelerate_series)
  }
  if (algorithm=="bisection"){
    result <- quantile_helper_bisection(this,p,accelerate_series)
  }
  return(result)
}

#' Prints univariate hermite_estimator object.
#' 
#'
#' @param x A hermite_estimator_univar object.
#' @param ... Other arguments passed on to methods used in printing.
#' @export
#' @examples
#' hermite_est <- hermite_estimator_univar(N = 10, standardize = TRUE)
#' print(hermite_est)
print.hermite_estimator_univar <- function(x, ...) {
  describe_estimator(x,"univariate")
}

#' Summarizes univariate hermite_estimator object.
#' 
#' Outputs key parameters of a univariate hermite_estimator object along with
#' estimates of the mean, standard deviation and deciles of the data that
#' the object has been updated with.
#'
#' @param object A hermite_estimator_univar object.
#' @param digits A numeric value. Number of digits to round to.
#' @param ... Other arguments passed on to methods used in summary.
#' @export
#' @examples
#' hermite_est <- hermite_estimator_univar(N = 10, standardize = TRUE)
#' hermite_est <- update_batch(hermite_est, rnorm(30))
#' summary(hermite_est)
summary.hermite_estimator_univar <- function(object, 
                              digits = max(3, getOption("digits") - 3),...) {
  describe_estimator(object,"univariate")
  if (object$num_obs > 2){
    cat("\n")
    cat(paste0("Mean = ",round(object$running_mean,digits), "\n"))
    cat(paste0("Standard Deviation = ",
               round(calculate_running_std(object),digits), "\n"))
    cat("Estimated Quantiles:\n")
    cumulative_probs <- seq(10,90,10)
    cum_probs_size <- length(cumulative_probs)
    quantile_values <- matrix(round(quant(object,p=cumulative_probs/100),
                                    digits),
                              nrow = 1, ncol = cum_probs_size, byrow = TRUE)
    colnames(quantile_values) <- paste0(cumulative_probs ,"%")
    rownames(quantile_values) <- ""
    print(quantile_values, quote = FALSE)
  }
}

spearmans.hermite_estimator_univar <- function(this, clipped) {
  stop("Spearman's Rho correlation estimation is not defined for the univariate 
       Hermite estimator")
}

kendall.hermite_estimator_univar <- function(this, clipped) {
  stop("Kendall Tau correlation estimation is not defined for the univariate 
       Hermite estimator")
}

Try the hermiter package in your browser

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

hermiter documentation built on Nov. 17, 2021, 1:07 a.m.