coxsnell.bc {mle.tools}

Bias-Corrected Maximum Likelihood Estimate(s)

Description

coxsnell.bc calculates the bias-corrected maximum likelihood estimate(s) using the bias formula introduced by Cox and Snell (1968).

Usage

coxsnell.bc(density, logdensity, n, parms, mle, lower = "-Inf",
  upper = "Inf", ...)

Arguments

density

An expression with the probability density function.

logdensity

An expression with the logarithm of the probability density function.

n

A numeric scalar with the sample size.

parms

A character vector with the parameter name(s) specified in the density and logdensity expressions.

mle

A numeric vector with the parameter estimate(s).

lower

The lower integration limit (lower = “-Inf” is the default).

upper

The upper integration limit (upper = “Inf” is the default).

...

Additional arguments passed to integrate function.

Details

The first, second and third-order partial log-density derivatives are analytically calculated via D function. The expected values of the partial log-density derivatives are calculated via integrate function.

Value

coxsnell.bc returns a list with five components (i) mle: the inputted maximum likelihood estimate(s), (ii) varcov: the expected variance-covariance evaluated at the inputted mle argument, (iii) mle.bc: the bias-corrected maximum likelihood estimate(s), (iv) varcov.bc: the expected variance-covariance evaluated at the bias-corrected maximum likelihood estimate(s) and (v) bias: the bias estimate(s).

If the numerical integration fails and/or the expected information is singular an error message is returned.

Author(s)

Josmar Mazucheli [e-mail]

See Also

deriv, D, expected.varcov, integrate, observed.varcov.

Examples

{library(mle.tools); library(fitdistrplus); set.seed(1)};

## Normal distribution
pdf <- quote(1 / (sqrt(2 * pi) * sigma) * exp(-0.5 / sigma ^ 2 * (x - mu) ^ 2))
lpdf <- quote(- log(sigma) - 0.5 / sigma ^ 2 * (x - mu) ^ 2)

x <- rnorm(n = 100, mean = 0.0, sd = 1.0)
{mu.hat <- mean(x); sigma.hat = sqrt((length(x) - 1) * var(x) / length(x))}

coxsnell.bc(density = pdf, logdensity = lpdf, n = length(x), parms = c("mu", "sigma"),
 mle = c(mu.hat, sigma.hat), lower = '-Inf', upper = 'Inf')

################################################################################

## Weibull distribution
pdf <- quote(shape / scale ^ shape * x ^ (shape - 1) * exp(-(x / scale) ^ shape))
lpdf <- quote(log(shape) - shape * log(scale) + shape * log(x) -
 (x / scale) ^ shape)

x <- rweibull(n = 100, shape = 1.5, scale = 2.0)

fit <- fitdist(data = x, distr = 'weibull')
fit$vcov

coxsnell.bc(density = pdf, logdensity = lpdf, n = length(x), parms = c("shape", "scale"),
 mle = fit$estimate, lower = 0)

################################################################################

## Exponentiated Weibull distribution
pdf <- quote(alpha * shape / scale ^ shape * x ^ (shape - 1) * exp(-(x / scale) ^ shape) *
 (1 - exp(-(x / scale) ^ shape)) ^ (alpha - 1))
lpdf <- quote(log(alpha) + log(shape) - shape * log(scale) + shape * log(x) -
 (x / scale) ^ shape + (alpha - 1) * log((1 - exp(-(x / scale) ^ shape))))

coxsnell.bc(density = pdf, logdensity = lpdf, n = 100, parms = c("shape", "scale", "alpha"),
 mle = c(1.5, 2.0, 1.0), lower = 0)

################################################################################

## Exponetial distribution
pdf <- quote(rate * exp(-rate * x))
lpdf <- quote(log(rate) - rate * x)

x <- rexp(n = 100, rate = 0.5)

fit <- fitdist(data = x, distr = 'exp')
fit$vcov

coxsnell.bc(density = pdf, logdensity = lpdf, n = length(x), parms = c("rate"),
 mle = fit$estimate, lower = 0)

################################################################################

## Gamma distribution
pdf <- quote(1 /(scale ^ shape * gamma(shape)) * x ^ (shape - 1) * exp(-x / scale))
lpdf <- quote(-shape * log(scale) - lgamma(shape) + shape * log(x) -
 x / scale)

x <- rgamma(n = 100, shape = 1.5, scale = 2.0)

fit <- fitdist(data = x, distr = 'gamma', start = list(shape = 1.5, scale =  2.0))
fit$vcov

coxsnell.bc(density = pdf, logdensity = lpdf, n = length(x), parms = c("shape", "scale"),
 mle = fit$estimate, lower = 0)

################################################################################

## Beta distribution
pdf <- quote(gamma(shape1 + shape2) / (gamma(shape1) * gamma(shape2)) * x ^ (shape1 - 1) *
 (1 - x) ^ (shape2 - 1))
lpdf <- quote(lgamma(shape1 + shape2) - lgamma(shape1) - lgamma(shape2) +
 shape1 * log(x) + shape2 * log(1 - x))

x <- rbeta(n = 100, shape1 = 2.0, shape2 = 2.0)

fit <- fitdist(data = x, distr = 'beta', start = list(shape1 = 2.0, shape2 =  2.0))
fit$vcov

coxsnell.bc(density = pdf, logdensity = lpdf, n = length(x), parms = c("shape1", "shape2"),
mle = fit$estimate, lower = 0, upper = 1)