#' Kumaraswamy distribution maximum likelihood estimation
#'
#' Uses Newton-Raphson to estimate the parameters of the Kumaraswamy
#' distribution.
#'
#' For the density function of the Kumaraswamy distribution see
#' [Kumaraswamy][extraDistr::Kumaraswamy].
#'
#' @param x a (non-empty) numeric vector of data values.
#' @param na.rm logical. Should missing values be removed?
#' @param ... `a0` is an optional starting value for the `a` parameter.
#' `rel.tol` is the relative accuracy requested, defaults
#' to `.Machine$double.eps^0.25`. `iterlim` is a positive integer
#' specifying the maximum number of iterations to be performed before the
#' program is terminated (defaults to `100`).
#'
#' @return `mlkumar` returns an object of [class][base::class] `univariateML`.
#' This is a named numeric vector with maximum likelihood estimates for `a`
#' and `b` and the following attributes:
#' \item{`model`}{The name of the model.}
#' \item{`density`}{The density associated with the estimates.}
#' \item{`logLik`}{The loglikelihood at the maximum.}
#' \item{`support`}{The support of the density.}
#' \item{`n`}{The number of observations.}
#' \item{`call`}{The call as captured my `match.call`}
#' @seealso [Kumaraswamy][extraDistr::Kumaraswamy] for the Kumaraswamy density.
#' @examples
#' AIC(mlkumar(USArrests$Rape / 100))
#' @references Jones, M. C. "Kumaraswamy's distribution: A beta-type
#' distribution with some tractability advantages." Statistical Methodology
#' 6.1 (2009): 70-81.
#'
#' Kumaraswamy, Ponnambalam. "A generalized probability density function
#' for double-bounded random processes." Journal of Hydrology 46.1-2 (1980):
#' 79-88.
#' @export
mlkumar <- function(x, na.rm = FALSE, ...) {
if (na.rm) x <- x[!is.na(x)] else assertthat::assert_that(!anyNA(x))
ml_input_checker(x)
assertthat::assert_that(min(x) > 0)
assertthat::assert_that(max(x) < 1)
dots <- list(...)
a0 <- if (!is.null(dots$a0)) dots$a0 else 1
rel.tol <- if (!is.null(dots$rel.tol)) {
dots$rel.tol
} else {
.Machine$double.eps^0.25
}
iterlim <- if (!is.null(dots$iterlim)) dots$iterlim else 100
logs <- log(x)
for (i in 1:iterlim) {
xa <- x^a0
T1 <- a0 * mean(logs / (1 - xa))
T2 <- a0 * mean(logs * (xa / (1 - xa)))
T3 <- mean(log(1 - xa))
f <- 1 / a0 * (1 + T1 + T2 / T3)
C <- mean(logs^2 * (xa / (1 - xa)^2))
D <- mean(logs^2 * (xa / (1 - xa)))
T1diff <- 1 / a0 * T1 + a0 * C
T2diff <- 1 / a0 * T2 + 1 / a0 * T2^2 + a0 * D
fdiff <- -1 / a0^2 * f + 1 / a0 *
(T1diff + T2diff / T3 + 1 / a0 * (T2 / T3)^2)
a <- a0 - f / fdiff
if (abs((a0 - a) / a0) < rel.tol) break
a0 <- a
}
if (i == iterlim) {
warning(paste0(
"The iteration limit (iterlim = ", iterlim, ") was reached",
" before the relative tolerance requirement (rel.tol = ",
rel.tol, ")."
))
}
## Given the shape, the scale is easy to compute.
b <- -1 / mean(log(1 - x^a))
object <- c(a = a, b = b)
class(object) <- "univariateML"
attr(object, "model") <- "Kumaraswamy"
attr(object, "density") <- "extraDistr::dkumar"
attr(object, "logLik") <-
length(x) * (log(a) + log(b) + (a - 1) * mean(log(x)) + -1 + 1 / b)
attr(object, "support") <- c(0, 1)
attr(object, "n") <- length(x)
attr(object, "call") <- match.call()
object
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.