R/binom_param.R

In twolodzko/twextras: twextras

#' Estimation of binomial parameters when both n, p are unknown
#'
#' @param x            a numeric vector.
#' @param na.rm        a logical value indicating whether \code{NA} values
#'                     should be stripped before the computation proceeds.
#' @param alpha,lambda parameters for methods described by DasGupta and Rubin (2005).
#'
#' @details
#'
#' Methods of estimating size of binomial population \eqn{n} are based on
#' paper by DasGupta and Rubin (2005) who describe the "new moment estimator"
#' \code{n1} and the second estimator \code{n2}, additionally, the method of
#' moments estimator \code{mmes} as described by Grevstad (2006) is available.
#'
#' For estimating probability of success \eqn{p}, the "common moments estimator"
#' \code{p.tilde} and and new estimator \code{p.hat} proposed by DasGupta and
#' Rubin (2005) are available.
#'
#' @references
#'
#' Grevstad, N. (2006).
#' Binomial Distribution: Sample Size Estimation.
#' Encyclopedia of Statistical Sciences.
#'
#' @references
#'
#' DasGupta, A., and Rubin, H. (2005).
#' Estimation of binomial parameters when both n, p are unknown.
#' Journal of Statistical Planning and Inference, 130(1), 391-404.
#'
#' @export

binom_n <- function(x, alpha = 1, na.rm = FALSE) {
  if (na.rm)
    x <- na.omit(x)

  k <- length(x)
  s2 <- var(x)
  xK <- max(x)
  xbar <- mean(x)

  n1 <- (xK^(alpha+1) * s2^alpha) / (xbar^alpha * (xK - xbar)^alpha)
  n2 <- xK + sum( qbeta(1/k, 1:(n1+1), n1 - (0:n1)) )

  ntilde <- (xbar^2)/(xbar - s2)
  w <- max(1 + sqrt(2), (xK-xbar)/s2)
  mmes <- if (xbar/s2 < 1 + 1/sqrt(2)) (s2 * w^2)/(w-1) else ntilde

  list( n1 = n1, n2 = n2, mmes = max(xK, mmes) )
}

#' @rdname binom_n
#' @export

binom_p <- function(x, alpha = 1, lambda = 1, na.rm = FALSE) {
  if (na.rm)
    x <- na.omit(x)

  k <- length(x)
  s2 <- var(x)
  xK <- max(x)
  xbar <- mean(x)

  phat <- ( xbar^((alpha+2)*lambda-1) * (xbar - s2)^(1-lambda) * (xK-xbar)^(alpha*lambda) ) /
    ( xK^((alpha+1)*lambda) * s2^(alpha*lambda) )

  list( p.tilde = max(0, 1 - s2/xbar), p.hat = phat )
}