R/unbGen.R

In innager/Umoments: Unbiased Central Moment Estimates

#' Unbiased central moment estimates
#'
#' Calculate unbiased estimates of central moments and their powers and products
#' up to specified order.
#'
#' Unbiased estimates up to the 6th order can be calculated. Second and third
#' orders contain estimates of the variance and third central moment, fourth
#' order includes estimates of fourth moment and squared variance
#' (\eqn{\mu_2^2}{\mu[2]^2}), fifth order - of fifth moment and a product of
#' second and third moments (\eqn{\mu_2 \mu_3}{\mu[2] \mu[3]}), sixth order - of
#' sixth moment, a product of second and fourth moments (\eqn{\mu_2
#' \mu_4}{\mu[2] \mu[4]}), squared third moment (\eqn{\mu_3^2}{\mu[3]^2}), and
#' cubed variance (\eqn{\mu_2^3}{\mu[2]^3}).
#'
#' @param smp sample.
#' @param order highest order of the estimates to calclulate. Estimates of lower
#'   orders will be included.
#'
#' @return A named vector of estimates of central moments and their powers and
#'   products up to \code{order}. The highest order available is 6th. The names
#'   of the elements are \code{"M2", "M3", "M4", "M5", "M6"} for corresponding
#'   central moments, \code{"M2M3", "M2M4"} for products of the moments (second
#'   and third, second and fourth), and \code{"M2pow2", "M2pow3", "M3pow2"} for
#'   powers of the moments - corresponding to estimates of squared variance,
#'   cubed variance, and squared third moment.
#' @seealso \code{\link{uMpool}} for two-sample pooled estimates.
#' @examples
#' smp <- rgamma(10, shape = 3)
#' uM(smp, 6)
#' @export
#'
#' @importFrom stats var
#' @importFrom utils combn
#'
uM <- function(smp, order) {
  n <- length(smp)
  order <- trunc(order)
  if (order < 0) {
    return()
  }
  if (order == 0) {
    return(1)
  }
  if (order == 1) {
    return(0)
  }
  res <- c(M2 = var(smp))
  if (order == 2) {
    return(res)
  }
  m1 <- mean(smp)
  m3 <- mean((smp - m1)^3)
  M3 <- uM3(m3, n)
  res <- c(res, M3 = M3)
  if (order == 3) {
    return(res)
  }
  m2 <- mean((smp - m1)^2)
  m4 <- mean((smp - m1)^4)
  M2pow2 <- uM2pow2(m2, m4, n)
  M4     <- uM4(    m2, m4, n)
  res <- c(res, M2pow2 = M2pow2, M4 = M4)
  if (order == 4) {
    return(res)
  }
  m5 <- mean((smp - m1)^5)
  M2M3 <- uM2M3(m2, m3, m5, n)
  M5   <- uM5(  m2, m3, m5, n)
  res <- c(res, M2M3 = M2M3, M5 = M5)
  if (order == 5) {
    return(res)
  }
  m6 <- mean((smp - m1)^6)
  M2pow3 <- uM2pow3(m2, m3, m4, m6, n)
  M3pow2 <- uM3pow2(m2, m3, m4, m6, n)
  M2M4   <- uM2M4(  m2, m3, m4, m6, n)
  M6     <- uM6(    m2, m3, m4, m6, n)
  res <- c(res, M2pow3 = M2pow3, M3pow2 = M3pow2, M2M4 = M2M4, M6 = M6)
  if (order == 6) {
    return(res)
  }
  if (order > 6) {
    warning("orders higher than 6th are not available")
    return(res)
  }
}

#' Pooled central moment estimates - two-sample
#'
#' Calculate unbiased pooled estimates of central moments and their powers and
#' products up to specified order.
#'
#' Pooled estimates up to the 6th order can be calculated. Second and third
#' orders contain estimates of the variance and third central moment, fourth
#' order includes estimates of fourth moment and squared variance
#' (\eqn{\mu_2^2}{\mu[2]^2}), fifth order - of fifth moment and a product of
#' second and third moments (\eqn{\mu_2 \mu_3}{\mu[2] \mu[3]}), sixth order - of
#' sixth moment, a product of second and fourth moments (\eqn{\mu_2
#' \mu_4}{\mu[2] \mu[4]}), squared third moment (\eqn{\mu_3^2}{\mu[3]^2}), and
#' cubed variance (\eqn{\mu_2^3}{\mu[2]^3}).
#'
#' @inherit uM params return
#' @param a vector of the same length as \code{smp} specifying categories of
#'   observations (should contain two unique values).
#' @seealso \code{\link{uM}} for one-sample unbiased estimates.
#' @examples
#' nsmp <- 23
#' smp <- rgamma(nsmp, shape = 3)
#' treatment <- sample(0:1, size = nsmp, replace = TRUE)
#' uMpool(smp, treatment, 6)
#' @export
uMpool <- function(smp, a, order) {
  if (length(unique(a)) != 2 | length(a) != length(smp)) {
    stop("design does not match sample")
  }
  smpx <- smp[a == max(a)]
  smpy <- smp[a == min(a)]
  nx <- length(smpx)
  ny <- length(smpy)
  order <- trunc(order)
  if (order < 0) {
    return()
  }
  if (order == 0) {
    return(1)
  }
  if (order == 1) {
    return(0)
  }
  mx1 <- mean(smpx)
  my1 <- mean(smpy)
  m2 <- mean(c((smpx - mx1)^2, (smpy - my1)^2))
  M2 <- uM2pool(m2, nx, ny)
  res <- c(M2 = M2)
  if (order == 2) {
    return(res)
  }
  m3 <- mean(c((smpx - mx1)^3, (smpy - my1)^3))
  M3 <- uM3pool(m3, nx, ny)
  res <- c(res, M3 = M3)
  if (order == 3) {
    return(res)
  }
  m4 <-  mean(c((smpx - mx1)^4, (smpy - my1)^4))
  M2pow2 <- uM2pow2pool(m2, m4, nx, ny)
  M4     <- uM4pool(    m2, m4, nx, ny)
  res <- c(res, M2pow2 = M2pow2, M4 = M4)
  if (order == 4) {
    return(res)
  }
  m5 <- mean(c((smpx - mx1)^5, (smpy - my1)^5))
  M2M3 <- uM2M3pool(m2, m3, m5, nx, ny)
  M5   <- uM5pool(  m2, m3, m5, nx, ny)
  res <- c(res, M2M3 = M2M3, M5 = M5)
  if (order == 5) {
    return(res)
  }
  m6 <- mean(c((smpx - mx1)^6, (smpy - my1)^6))
  M2pow3 <- uM2pow3pool(m2, m3, m4, m6, nx, ny)
  M3pow2 <- uM3pow2pool(m2, m3, m4, m6, nx, ny)
  M2M4   <- uM2M4pool(  m2, m3, m4, m6, nx, ny)
  M6     <- uM6pool(    m2, m3, m4, m6, nx, ny)
  res <- c(res, M2pow3 = M2pow3, M3pow2 = M3pow2, M2M4 = M2M4, M6 = M6)
  if (order == 6) {
    return(res)
  }
  if (order > 6) {
    warning("orders higher than 6th are not available")
    return(res)
  }
}