mle.tools: An R Package for Maximum Likelihood Bias Correction

Abstract:

Recently, uploaded the package mle.tools to CRAN. It can be used for bias corrections of maximum likelihood estimates through the methodology proposed by . The main function of the package, coxsnell.bc(), computes the bias corrected maximum likelihood estimates. Although in general, the bias corrected estimators may be expected to have better sampling properties than the uncorrected estimators, analytical expressions from the formula proposed by are either tedious or impossible to obtain. The purpose of this paper is twofolded: to introduce the mle.tools package, especially the coxsnell.bc() function; secondly, to compare, for thirty one continuous distributions, the bias estimates from the coxsnell.bc() function and the bias estimates from analytical expressions available in the literature. We also compare, for five distributions, the observed and expected Fisher information. Our numerical experiments show that the functions are efficient to estimate the biases by the Cox-Snell formula and for calculating the observed and expected Fisher information.

Cite PDF Tweet

Published

Nov. 1, 2017

Received

Apr 4, 2017

Citation

Mazucheli, et al., 2017

Volume

Pages

9/2

268 - 290


1 Introduction

Since it was proposed by Fisher in a series of papers from 1912 to 1934, the maximum likelihood method for parameter estimation has been employed to several issues in statistical inference, because of its many appealing properties. For instance, the maximum likelihood estimators, hereafter referred to as MLEs, are asymptotically unbiased, efficient, consistent, invariant under parameter transformation and asymptotically normally distributed . Most properties that make the MLEs attractive depend on the sample size, hence such properties as unbiasedness, may not be valid for small samples or even moderate samples . Indeed, the maximum likelihood method produces biased estimators, i.e., expected values of MLEs differ from the real true parameter values providing systematic errors. In particular, these estimators typically have biases of order O(n1), thus these errors reduce as sample size increases .

Applying the corrective Cox-Snell methodology, many researchers have developed nearly unbiased estimators for the parameters of several probability distributions. Interested readers can refer to , , , , , , , , , , , , , , , , , , and references cited therein.

In general, the Cox-Snell methodology is efficient for bias corrections. However, obtaining analytical expressions for some probability distributions, mainly for those indexed by more than two parameters, can be notoriously cumbersome or impossible. presented Maple and Mathematica scripts that may be used to calculate closed form analytic expressions for bias corrections using the Cox-Snell formula. They tested the scripts for 20 two-parameter continuous probability distributions, and the results were compared with those published in earlier works. In the same direction, researchers from the University of Illinois, at Urbana-Champaign, have developed a Mathematica program, entitled “CSCK MLE Bias Calculation” that enables the user to calculate the analytic Cox-Snell MLE bias vectors for various probability distributions with up to four unknown parameters. It is important to mention that both, Maple and Mathematica , are commercial softwares.

In this paper, our objective is to introduce a new contributed R package, namely mle.tools that computes the expected/observed Fisher information and the bias corrected estimates by the methodology proposed by . The theoretical background of the methodology is presented in Section 2. Details about the mle.tools package are described in Section 3. Closed form solutions of bias corrections are collected from the literature for a large number of distributions and compared to the output from the coxsnell.bc() function, see Section 4. In Section 5, we compare various estimates of Fisher’s information, considering a real application from the literature. Finally, Section 6 contains some concluding remarks and directions for future research.

2 Overview of the Cox-Snell methodology

Let X1,,Xn be n be independent random variables with probability density function f(xiθ) depending on a p-dimensional parameter vector θ=(θ1,,θp). Without loss of generality, let l=l(θx) be the log-likelihood function for the unknown p-dimensional parameter vector θ given a sample of n observations. We shall assume some regularity conditions on the behavior of l(θx) .

The joint cumulants of the derivatives of l are given by: κij=E[2lθiθj],κijl=E[3lθiθjθl],κij,l=E[(2lθiθj)(lθl)],κij(l)=κijθl for i,j,l=1,,p.

The bias expression of the sth element of θ^, the MLEs of θ, when the sample data are independent, but not necessarily identically distributed, was proposed by : (1)B(θ^s)=i=1pj=1pl=1pκsiκjl[0.5κijl+κij,l]+O(n2), where s=1,,p and κij is the (i,j)th element of the inverse of the negative of the expected Fisher information.

Thereafter, noticed that equation (1) holds even if the data are non-independent, and it can be re-expressed as: (2)B(θ^s)=i=1pκsij=1pl=1p[κij(l)0.5κijl]κjl+O(n2).

Defining aij(l)=κij(l)0.5κijl, A(l)={aij(l)} and K=[κij], the expected Fisher information matrix for i,j,l=1,,n, the bias expression for θ^ in matrix notation is: B(θ^)=K1Avec(K1)+O(n2), where vec(K1) is the vector obtained by stacking the columns of K1 and A={A1Ap}.

Finally, the bias corrected MLE for θs can be obtained as: (3)θ~s=θ^sB^(θ^s). Alternatively, using matrix notation the bias corrected MLEs can be expressed as : (4)θ~=θ^K^1A^vec(K^1), where K^=K|θ=θ^ and A^=A|θ=θ^.

3 The mle.tools package details

The current version of the mle.tools package, uploaded to CRAN in February, 2017, has implemented three functions — observed.varcov(), expected.varcov() and coxsnell.bc() — which are of great interest in data analysis based on MLEs. These functions calculate, respectively, the observed Fisher information, the expected Fisher information and the bias corrected MLEs using the bias formula in (1). The above mentioned functions can be applied to any probability density function whose terms are available in the derivatives table of the D() function (see “deriv.c” source code for further details). Integrals, when required, are computed numerically via the integrate() function. Below are some mathematical details of how the returned values from the three functions are calculated.

Let X1,,Xn be independent and identical random variables with probability density function f(xiθ) depending on a p-dimensional parameter vector θ=(θ1,,θp). The (j,k)th element of the observed, Hjk, and expected, Ijk, Fisher information are calculated, respectively, as Hjk=i=1n2θjθklogf(xiθ)|θ=θ^ and Ijk=n×E(2θjθklogf(xθ))=n×X2θjθklogf(xθ)×f(xθ)dx|θ=θ^, where j,k=1,,p, θ^ is the MLE of θ and X denotes the support of the random variable X.

The observed.varcov() function is as follows:

function (logdensity, X, parms, mle)

where logdensity is an R expression of the log of the probability density function, X is a numeric vector containing the observations, parms is a character vector of the parameter name(s) specified in the logdensity expression and mle is a numeric vector of the parameter estimate(s). This function returns a list with two components (i) mle: the inputed MLEs and (ii) varcov: the observed variance-covariance evaluated at the inputed MLE argument. The elements of the Hessian matrix are calculated analytically.

The functions expected.varcov() and coxsnell.bc() have the same arguments and are as follows:

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

where density and logdensity are R expressions of the probability density function and its logarithm, respectively, n is a numeric scalar of the sample size, parms is a character vector of the parameter names(s) specified in the density and log-density expressions, mle is a numeric vector of the parameter estimates, lower is the lower integration limit (-Inf is the default), upper is the upper integration limit (Inf is the default) and ... are additional arguments passed to the integrate() function. The expected.varcov() function returns a list with two components:

$mle

the inputed MLEs and

$varcov

the expected covariance evaluated at the inputed MLEs.

The coxsnell.bc() function returns a list with five components:

$mle

the inputed MLEs,

$varcov

the expected variance-covariance evaluated at the inputed MLEs,

$mle.bc

the bias corrected MLEs,

$varcov.bc

the expected variance-covariance evaluated at the bias corrected MLEs

$bias

the bias estimate(s).

Furthermore, the bias corrected MLE of θs, s=1,,p denoted by θs~ is calculated as θs~=θ^sB^(θ^s), where θ^s is the MLE of θs and B^(θ^s)=j=1pk=1pl=1pκsjκkl[0.5κjkl+κjk,l]|θ=θ^, where κjk is the (j,k)th element of the inverse of the negative of the expected Fisher information, κjkl=nX3θjθkθllogf(xθ)f(xθ)dx|θ=θ^,

κjk,l=nX2θjθklogf(xθ)θllogf(xθ)f(xθ)dx|θ=θ^ and X denotes the support of the random variable X.

It is important to emphasize that first, second and third-order partial log-density derivatives are analytically calculated via the D() function, while integrals are computed numerically, using the integrate() function. Furthermore, if numerical integration fails and/or the expected/observed information is singular, an error message is returned.

4 Comparative study

In order to evaluate the robustness of the coxsnell.bc() function, we compare, through real applications, the estimated biases obtained from the package and from the analytical expressions for a total of thirty one continuous probability distributions. The analytical expressions for each distribution, named as distname.bc(), can be found in the supplementary file “analyticalBC.R”. For example, the entry lindley.bc(n, mle) evaluates the bias estimates locally at n and mle values.

In the sequel, the probability density function, the analytical Cox-Snell expressions and the bias estimates are provided for: Lindley, inverse Lindley, inverse Exponential, Shanker, inverse Shanker, Topp-Leone, Lévy, Rayleigh, inverse Rayleigh, Half-Logistic, Half-Cauchy, Half-Normal, Normal, inverse Gaussian, Log-Normal, Log-Logistic, Gamma, inverse Gamma, Lomax, weighted Lindley, generalized Rayleigh, Weibull, inverse Weibull, generalized Half-Normal, inverse generalized Half-Normal, Marshall-Olkin extended Exponential, Beta, Kumaraswamy, inverse Beta, Birnbaum-Saunders and generalized Pareto distributions.

It is noteworthy that analytical bias corrected expressions are not reported in the literature for the Lindley, Shanker, inverse Shanker, Lévy, inverse Rayleigh, half-Cauchy, inverse Weibull, inverse generalized half-normal and Marshall-Olkin extended exponential distributions.

According to all the results presented below, we observe concordance between the bias estimates given by the coxsnell.bc() function and the analytical expression(s) for 28 out the 31 distributions. The distributions which did not agree with the coxsnell.bc() function were the beta, Kumaraswamy and inverse beta distributions. Perhaps there are typos either in our typing or in the analytical expressions reported by , and . Having this view, we recalculated the analytical expressions for the biases. For the beta and inverse beta distributions, our recalculated analytical expressions agree with the results returned by the coxsnell.bc() function, so there are actually typos in the expression of and . For the Kumaraswamy, we could not evaluate the analytical expression given by the author but we compare the results from coxsnell.bc() function with a numerical evaluation in Maple and the results are exactly equals.

  1. One-parameter Lindley distribution with scale parameter θ f(xθ)=θ21+θ(1+x)exp(θx),x>0.

    Bias expression (not previously reported in the literature): (5)B(θ^)=(θ3+6θ2+6θ+2)(θ+1)θn(θ2+4θ+2)2.

    Using the data set from we have n=100, θ^=0.1866 and se^(θ^)=0.0133. Evaluating the analytical expression (5) and the coxsnell.bc() function, we have, respectively,

    lindley.bc(n = 100, mle = 0.1866)
    ##     theta
    ## 0.0009546
    pdf <- quote(theta^2 / (theta + 1) * (1 + x) * exp(-theta * x))
    lpdf <- quote(2 * log(theta) - log(1 + theta) - theta * x)
    coxsnell.bc(density = pdf, logdensity = lpdf, n = 100,
                parms = c("theta"), mle = 0.1866, lower = 0)$bias
    ##     theta
    ## 0.0009546
  2. Inverse Lindley distribution with scale parameter θ f(xθ)=θ21+θ(1+xx3)exp(θx),x>0.

    Bias expression : (6)B(θ^)=(θ+1)θ(θ3+6θ2+6θ+2)n(θ2+4θ+2)2.

    Using the data set from we have n=58, θ^=60.0016 and se^(θ^)=7.7535. Evaluating the analytical expression (6) and the coxsnell.bc() function, we have, respectively,

    invlindley.bc(n = 58, mle =  60.0016)
    ## theta
    ## 1.017
    pdf <- quote(theta^2 / (theta + 1) * ((1 + x) / x^3) *
                 exp(-theta / x))
    lpdf <- quote(2 * log(theta) - log(1 + theta) - theta / x)
    coxsnell.bc(density = pdf, logdensity = lpdf, n = 58,
                parms = c("theta"), mle =  60.0016, lower = 0)$bias
    ## theta
    ## 1.017
  3. Inverse exponential distribution with rate parameter θ f(xθ)=θx2exp(θx),x>0.

    Bias expression : (7)B(θ^)=θn.

    Using the data set from , we have n=30, θ^=11.1786 and se^(θ^)=2.0409. Evaluating the analytical expression (7) and the coxsnell.bc() function, we have, respectively,

    invexp.bc(n = 30, mle = 11.1786)
    ##  theta
    ## 0.3726
    pdf <- quote(theta / x^2 * exp(- theta / x))
    lpdf <- quote(log(theta) - theta / x)
    coxsnell.bc(density = pdf, logdensity = lpdf, n = 30,
                parms = c("theta"), mle = 11.1786, lower = 0)$bias
    ##  theta
    ## 0.3726
  4. Shanker distribution with scale parameter θ f(xθ)=θ2θ2+1(θ+x)exp(θx),x>0.

    For bias expression (not previously reported in the literature, see the “analyticalBC.R” file.

    Using the data set from , we have n=31, θ^=0.0647 and se^(θ^)=0.0082. Evaluating the analytical expression and the coxsnell.bc() function, we have, respectively,

    shanker.bc(n = 31, mle = 0.0647)
    ##    theta
    ## 0.001035
    pdf <- quote(theta^2 / (theta^2 + 1) *  (theta + x) *
                 exp(-theta * x))
    lpdf <- quote(2*log(theta) - log(theta^2 + 1) + log(theta + x) -
                  theta * x)
    coxsnell.bc(density = pdf, logdensity = lpdf, n = 31,
                parms = c("theta"), mle = 0.0647, lower = 0)$bias
    ##    theta
    ## 0.001035
  5. Inverse Shanker distribution with scale parameter θ f(xθ)=θ21+θ2(1+θxx3)exp(θx),x>0.

    Bias expression (not previously reported in the literature): (8)B(θ^)=θ3+2θn(θ2+1).

    Using the data set from , we have n=58, θ^=59.1412 and se^(θ^)=7.7612. Evaluating the analytical expression (8) and the coxsnell.bc() function, we have, respectively,

    invshanker.bc(n = 58, mle = 59.1412)
    ## theta
    ##  1.02
    pdf <- quote(theta^2 / (theta^2 + 1) * (theta * x + 1) /
                 x^3 * exp(-theta / x))
    lpdf <- quote(log(theta) - 2 * log(x) - theta / x)
    coxsnell.bc(density = pdf, logdensity = lpdf, n = 58,
                parms = c("theta"), mle = 59.1412, lower = 0)$bias
    ## theta
    ##  1.02
  6. Topp-Leone distribution with shape parameter ν f(xν)=2ν(1x)xν1(2x)ν1,0<x<1.

    Bias expression : (9)B(ν^)=νn.

    Using the data set from , we have n=107, ν^=2.0802 and se^(ν^)=0.2011. Evaluating the analytical expression (9) and the coxsnell.bc() function, we have, respectively,

    toppleone.bc(n = 107, mle = 2.0802)
    ##      nu
    ## 0.01944
    pdf <- quote(2 * nu * x^(nu - 1) * (1 - x) * (2 - x)^(nu - 1))
    lpdf <- quote(log(nu) + nu * log(x) + log(1 - x) + (nu - 1) *
                  log(2 - x))
    coxsnell.bc(density = pdf, logdensity = lpdf, n = 107,
                parms = c("nu"), mle = 2.0802, lower = 0, upper = 1)$bias
    ##      nu
    ## 0.01944
  7. One-parameter Lévy distribution with scale parameter σ f(xσ)=σ2πx32exp(σ2x),x>0.

    Bias expression (not previously reported in the literature): (10)B(σ^)=2σn.

    Using the data set from , we have n=361, σ^=4.4461 and se^(σ^)=0.3309. Evaluating the analytical expression (10) and the coxsnell.bc() function, we have, respectively,

    levy.bc(n = 361, mle = 4.4460)
    ##   sigma
    ## 0.02463
    pdf <- quote(sqrt(sigma / (2 * pi)) * exp(-0.5 * sigma / x) /
                 x^(3 / 2))
    lpdf <- quote(0.5 * log(sigma) - 0.5 * sigma / x - (3 / 2) * log(x))
    coxsnell.bc(density = pdf, logdensity = lpdf, n = 361,
                parms = c("sigma"), mle = 4.4460, lower = 0)$bias
    ##   sigma
    ## 0.02463
  8. Rayleigh distribution with scale parameter σ f(xσ)=xσ2exp(x22σ2),x>0.

    Bias expression : (11)B(σ^)=σ8n.

    Using the data set from , we have n=69, σ^=1.2523 and se^(σ^)=0.0754. Evaluating the analytical expression (11) and the coxsnell.bc() function, we have, respectively,

    rayleigh.bc(n = 69, mle = 1.2522)
    ##     sigma
    ## -0.002268
    pdf <- quote(x / sigma^2 * exp(- 0.5 * (x / sigma)^2))
    lpdf <- quote(- 2 * log(sigma) - 0.5 * x^2 / sigma^2)
    coxsnell.bc(density = pdf, logdensity = lpdf, n = 69,
                parms = c("sigma"), mle = 1.2522, lower = 0)$bias
    ##     sigma
    ## -0.002268
  9. Inverse Rayleigh distribution with scale parameter σ f(xσ)=2σ2x3exp(σx2),x>0.

    Bias expression (not previously reported in the literature): (12)B(σ^)=3σ8n.

    Using the data set from , we have n=63, σ^=2.8876 and se^(σ^)=0.1819. Evaluating the analytical expression (12) and the coxsnell.bc() function, we have, respectively,

    invrayleigh.bc(n = 63, mle = 2.8876)
    ##   sigma
    ## 0.01719
    pdf <- quote(2 * sigma^2 / x^3 * exp(-sigma^2 / x^2))
    lpdf <- quote(2 * log(sigma) - sigma^2 / x^2)
    coxsnell.bc(density = pdf, logdensity = lpdf, n = 63,
                parms = c("sigma"), mle = 2.8876, lower = 0)$bias
    ##   sigma
    ## 0.01719
  10. Half-logistic distribution with scale parameter σ f(xσ)=2exp(xσ)σ[1+exp(xσ)]2,x>0.

    Bias expressions : (13)B(σ^)=0.05256766607σn.

    Using the data set from , we have n=34, σ^=1.3926 and se^(σ^)=0.2056. Evaluating the analytical expression (12) and the coxsnell.bc() function, we have, respectively,

    halflogistic.bc(n = 34, mle = 1.3925)
    ##     sigma
    ## -0.002153
    pdf <- quote((2/sigma) * exp(-x / sigma) / (1 + exp(-x / sigma))^2)
    lpdf <- quote(-log(sigma) - x / sigma - 2 * log(1 + exp(-x / sigma)))
    coxsnell.bc(density = pdf, logdensity = lpdf, n = 34,
                parms = c("sigma"), mle = 1.3925, lower = 0)$bias
    ##     sigma
    ## -0.002153
  11. Half-Cauchy distribution with scale parameter σ f(xσ)=2πσσ2+x2,x>0.

    Bias expression (not previously reported in the literature): (14)B(σ^)=σn.

    Using the data set from , we have n=64, σ^=28.3345 and se^(σ^)=4.4978. Evaluating the analytical expression (14) and the coxsnell.bc() function, we have, respectively,

    halfcauchy.bc(n = 64, mle = 28.3345)
    ##  sigma
    ## 0.4427
    pdf <- quote( 2 / pi * sigma / (x^2 + sigma^2))
    lpdf <- quote(log(sigma) - log(x^2 + sigma^2))
    coxsnell.bc(density = pdf, logdensity = lpdf, n = 64,
                parms = c("sigma"), mle = 28.3345, lower = 0)$bias
    ##  sigma
    ## 0.4456
  12. Half-normal distribution with scale parameter σ f(xσ)=2π1σexp(x22σ2),x>0.

    Bias expressions : (15)B(σ^)=σ4n.

    Using the data set from , we have n=69, σ^=1.5323 and se^(σ^)=0.1304. Evaluating the analytical expression (15) and the coxsnell.bc() function, we have, respectively,

    halfnormal.bc(n = 69, mle = 1.5323)
    ##     sigma
    ## -0.005552
    pdf <- quote(sqrt(2) / (sqrt(pi) * sigma) * exp(-x^2 / (2 * sigma^2)))
    lpdf <- quote(-log(sigma) - x^2 / sigma^2 / 2 )
    coxsnell.bc(density = pdf, logdensity = lpdf, n = 69,
                parms = c("sigma"), mle = 1.5323, lower = 0)$bias
    ##     sigma
    ## -0.005552
  13. Normal distribution with mean μ and standard deviation σ f(xμ,σ)=12πσexp[(xμ)22σ2],x(,).

    Bias expressions : (16)B(μ^)=0 and B(σ^)=3σ4n.

    Using the data set from , we have n=23, μ^=4.1506, σ^=0.5215, se^(μ^)=0.1087 and se^(σ^)=0.0769. Evaluating the analytical expressions (16) and the coxsnell.bc() function, we have, respectively,

    normal.bc(n = 23, mle = c(4.1506, 0.5215))
    ##       mu    sigma
    ##  0.00000 -0.01701
    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)
    coxsnell.bc(density = pdf, logdensity = lpdf, n = 23,
                parms = c("mu", "sigma"), mle = c(4.1506, 0.5215))$bias
    ##         mu      sigma
    ## -4.071e-13 -1.701e-02
  14. Inverse Gaussian distribution with mean μ and shape λ f(xμ,λ)=λ2πx3exp[λ(xμ)22xμ2],x>0.

    Bias expressions : (17)B(μ^)=0 and B(λ^)=3λn.

    Using the data set from , we have n=46, μ^=3.6067, λ^=1.6584, se^(μ^)=0.7843 and se^(λ^)=0.3458. Evaluating the analytical expressions (17) and the coxsnell.bc() function, we have, respectively,

    invgaussian.bc(n = 46, mle =  c(3.6065, 1.6589))
    ##     mu lambda
    ## 0.0000 0.1082
    pdf <- quote(sqrt(lambda / (2 * pi * x^3)) *
                 exp(-lambda * (x - mu)^2 / (2 * mu^2 * x)))
    lpdf <- quote(0.5 * log(lambda) - lambda * (x - mu)^2 /
                  (2 * mu^2 * x))
    coxsnell.bc(density = pdf, logdensity = lpdf, n = 46,
                parms = c("mu", "lambda"), mle = c(3.6065, 1.6589), 
                lower = 0)$bias
    ##        mu    lambda
    ## 3.483e-07 1.082e-01
  15. Log-normal distribution with location μ and scale σ f(xμ,σ)=12πxσexp[(logxμ)2σ2],x>0.

    Bias expressions : (18)B(μ^)=0 and B(σ^)=3σ4n.

    Using the data set from , we have n=30, μ^=2.164, σ^=1.1765, se^(μ^)=0.2148 and se^(σ^)=0.1519. Evaluating the analytical expressions (18) and the coxsnell.bc() function, we have, respectively,

    lognormal.bc(n = 30, mle = c(2.1643, 1.1765))
    ##       mu    sigma
    ##  0.00000 -0.02941
    pdf <- quote(1 / (sqrt(2 * pi) * x * sigma) *
                 exp(-0.5 * (log(x) - mu)^2 / sigma^2))
    lpdf <- quote(-log(sigma) - 0.5 * (log(x) - mu)^2 / sigma^2)
    coxsnell.bc(density = pdf, logdensity = lpdf, n = 30,
                parms = c("mu", "sigma"), mle = c(2.1643, 1.1765), 
                lower = 0)$bias
    ##         mu      sigma
    ## -5.952e-09 -2.941e-02
  16. Log-logistic distribution with shape β and scale α f(xα,β)=(β/α)(x/α)β1[1+(x/α)β]2,x>0.

    For bias expressions, see .

    From we have n=19, α^=6.2542, β^=1.1732, se^(α^)=2.1352, se^(β^)= 0.2239, B^(α^)=0.3585 and B^(β^)=0.0789. Evaluating the coxsnell.bc() function, we have:

    pdf <- quote((beta / alpha) * (x / alpha)^(beta - 1) /
                 (1 + (x / alpha)^beta)^2)
    lpdf <- quote(log(beta) - log(alpha) + (beta - 1) * log(x / alpha) -
                  2 * log(1 + (x / alpha)^beta))
    coxsnell.bc(density = pdf, logdensity = lpdf, n = 19,
                parms =  c("alpha", "beta"), mle = c(6.2537, 1.1734),
                lower = 0)$bias
    ##   alpha    beta
    ## 0.35854 0.07883
  17. Gamma distribution with shape α and rate λ f(xα,λ)=λαΓ(α)xα1exp(λx),x>0.

    Bias expressions : (19)B(α^)=α[Ψ(α)αΨ(α)]22n[αΨ(α)1]2 and (20)B(λ^)=λ[2α(Ψ(α))23Ψ(α)αΨ(α)]2n[αΨ(α)1]2.

    Using the data set from , we have n=254, α^=4.0083, λ^=0.0544, se^(α^)=0.3413 and se^(λ^)=0.0049. Evaluating the analytical expressions (19), (20) and the coxsnell.bc() function, we have, respectively,

    gamma.bc(n = 254, mle = c(4.0082, 0.0544))
    ##     alpha    lambda
    ## 0.0448278 0.0006618
    pdf <- quote((lambda^alpha) / gamma(alpha) * x^(alpha - 1) *
                 exp(-lambda *x))
    lpdf <- quote(alpha * log(lambda) - lgamma(alpha) + alpha * log(x) -
                  lambda * x)
    coxsnell.bc(density = pdf, logdensity = lpdf, n = 254,
                parms = c("alpha", "lambda"), mle = c(4.0082, 0.0544),
                lower = 0)$bias
    ##     alpha    lambda
    ## 0.0448278 0.0006618
  18. Inverse gamma distribution with shape α and scale β f(xα,β)=1Γ(α)βαxα1exp(xβ),x>0.

    Bias expressions : (21)B(α^)=0.5α2Ψ(α)+0.5Ψ(α)α1nα(Ψ(α)1)2 and (22)B(β^)=β(0.5αΨ(α)1.5Ψ(α)+(Ψ(α))2α)n(Ψ(α)α1.0)2.

    Using the data set from , we have n=31, α^=1.0479, β^=5.491, se^(α^)=0.2353 and se^(β^)=1.5648. Evaluating the analytical expressions (21), (22) and the coxsnell.bc() function, we have, respectively,

    invgamma.bc(n = 31, mle = c(5.4901, 1.0479))
    ##    beta   alpha
    ## 0.60849 0.08388
    pdf <- quote(beta^alpha / gamma(alpha) * x^(-alpha - 1) *
                 exp(-beta / x))
    lpdf <- quote(alpha * log(beta) - lgamma(alpha) -
                  alpha * log(x) - beta / x)
    coxsnell.bc(density = pdf, logdensity = lpdf, n = 31,
                parms = c("beta", "alpha"), mle = c(5.4901, 1.0479), 
                lower = 0)$bias
    ##    beta   alpha
    ## 0.60847 0.08388
  19. Lomax distribution with shape α and scale β f(xα,β)=αβ(1+βx)(α+1),x>0.

    Bias expressions : (23)B(α^)=2α(α+1)(α2+α2)(α+3)n and (24)B(β^)=2β(α+1.6485)(α+0.3934)(α1.5419)nα(α+3).

    Using the data set from , we have n=179, α^=4.9103, β^=0.0028, se^(α^)=0.6208 and se^(β^)=3.4803×104. Evaluating the analytical expressions (23), (24) and the coxsnell.bc() function, we have, respectively,

    lomax.bc(n = 179, mle = c(4.9103, 0.0028))
    ##      alpha       beta
    ##  1.281e+00 -9.438e-05
    pdf <- quote(alpha * beta / (1 + beta * x)^(alpha + 1))
    lpdf <- quote(log(alpha) + log(beta) - (alpha + 1) *
                  log(1 + beta * x))
    coxsnell.bc(density = pdf, logdensity = lpdf, n = 179,
                parms = c("alpha", "beta"), mle = c(4.9103, 0.0028),
                lower = 0)$bias
    ##      alpha       beta
    ##  1.281e+00 -9.439e-05
  20. Weighted Lindley distribution with shape α and scale θ f(xα,θ)=θα+1(θ+α)Γ(α)xα1(1+x)exp(θx),x>0.

    For bias expressions, see :

    Using the data set from , we have n=69, α^=22.8889, θ^=9.6246, se^(α^)=3.9507 and se^(θ^)=1.6295. Evaluating the analytical expressions and the coxsnell.bc function, we have, respectively,

    wlindley.bc(n = 69, mle = c(22.8889, 9.6246))
    ##  alpha  theta
    ## 1.0070 0.4167
    pdf <- quote(theta^(alpha + 1) / ((theta + alpha) * gamma(alpha)) *
                 x^(alpha - 1) * (1 + x) * exp(-theta * x))
    lpdf <- quote((alpha + 1) * log(theta) + alpha * log(x) -
                  log(theta + alpha) - lgamma(alpha) - theta * x)
    coxsnell.bc(density = pdf, logdensity = lpdf, n = 69,
                parms = c("alpha", "theta"), mle = c(22.8889, 9.6246),
                lower = 0)$bias
    ##  alpha  theta
    ## 1.0068 0.4166
  21. Generalized Rayleigh with shape α and scale θ f(xβ,μ)=2θα+1Γ(α+1)x2α+1exp(θx2),x>0.

    For bias expressions, see :

    Using the data set from , we have n=384, θ^=0.5195, α^=0.0104, se^(θ^)=0.2184 and se^(α^)=0.0014. Evaluating the analytical expressions and the coxsnell.bc() function, we have, respectively,

    generalizedrayleigh.bc(n =  384, mle = c(0.5195, 0.0104))
    ##     alpha     theta
    ## 1.035e-02 8.865e-05
    pdf <- quote(2 * theta^(alpha + 1) / gamma(alpha + 1) *
                 x^(2 * alpha + 1) * exp(-theta * x^2 ))
    lpdf <- quote((alpha + 1) * log(theta) - lgamma(alpha + 1) +
                  2 * alpha * log(x) - theta * x^2)
    coxsnell.bc(density = pdf, logdensity = lpdf,  n = 384,
                parms = c("alpha", "theta"), mle = c(0.5195, 0.0104),
                lower = 0)$bias
    ##     alpha     theta
    ## 1.035e-02 8.865e-05
  22. Weibull distribution with shape β and scale μ f(xβ,μ)=βμβxβ1exp(xμ)β,x>0.

    Bias expressions (the expressions below differs from ): (25)B(μ^)=μ(0.55433244950.3698145397β)nβ2 and (26)B(β^)=1.379530692βn.

    From , we have n=50, μ^=2.5752, β^=38.0866, se^(μ^)=0.2299 and se^(β^)=2.2299. Evaluating the analytical expression (25), (26) and the coxsnell.bc() function, we have, respectively,

    weibull.bc(n = 50, mle = c(38.0866, 2.5751))
    ##       mu     beta
    ## -0.04572  0.07105
    pdf <- quote(beta / mu^beta * x^(beta - 1) *
                 exp(-(x / mu)^beta))
    lpdf <- quote(log(beta) - beta * log(mu) + beta * log(x) -
                  (x / mu)^beta)
    coxsnell.bc(density = pdf, logdensity = lpdf, n = 50,
                parms = c("mu", "beta"), mle = c(38.0866, 2.5751), 
                lower = 0)$bias
    ##       mu     beta
    ## -0.04572  0.07105
  23. Inverse Weibull distribution with shape β and scale μ f(xβ,α)=βμβx(β+1)exp[(μx)β],x>0.

    Bias expressions (not previously reported in the literature): (27)B(β^)=1.379530690βn and (28)B(μ^)=μ(0.3698145391β+0.5543324494)nβ2.

    Using the data set from , we have n=100, β^=1.769, μ^=1.8917, se^(β^)=0.1119 and se^(μ^)=0.1138. Evaluating the analytical expressions (27), (28) and the coxsnell.bc() function, we have, respectively,

    inverseweibull.bc(n = 100, mle = c(1.7690, 1.8916))
    ##     beta       mu
    ## 0.024404 0.007305
    pdf <- quote(beta * mu^beta * x^(-beta - 1) *
                 exp(-(mu / x)^beta))
    lpdf <- quote(log(beta) + beta * log(mu) - beta * log(x) -
                  (mu / x)^beta)
    coxsnell.bc(density = pdf, logdensity = lpdf, n = 100,
                parms = c("beta", "mu"), mle = c(1.7690, 1.8916), 
                lower = 0)$bias
    ##     beta       mu
    ## 0.024404 0.007305
  24. Generalized half-normal distribution with shape α and scale θ f(xα,θ)=2παθαxα1exp[12(xθ)2α].

    Bias expressions : (29)B(α^)=1.483794456αn and (30)B(θ^)=(0.29534976610.3665611957α)θnα2.

    Using the data set from , we have n=119, α^=3.8096, θ^=4.9053, se^(α^)=0.2758 and se^(θ^)=0.0913. Evaluating the analytical expressions (29), (30) and the coxsnell.bc() function, we have, respectively,

    genhalfnormal.bc(n = 119, mle = c(3.8095, 4.9053))
    ##     alpha     theta
    ##  0.047500 -0.003127
    pdf <- quote(sqrt(2 / pi) * alpha / theta^alpha * x^(alpha - 1)*
                 exp(- 0.5 * (x / theta)^(2 * alpha) ))
    lpdf <- quote(log(alpha) - alpha * log(theta) + alpha * log(x) -
                  0.5 * (x / theta)^(2 * alpha))
    coxsnell.bc(density = pdf, logdensity = lpdf, n = 119,
                parms = c("alpha", "theta"), mle = c(3.8095, 4.9053),
                lower = 0)$bias
    ##     alpha     theta
    ##  0.047500 -0.003127
  25. Inverse generalized half-normal distribution with shape α and scale θ f(xα,θ)=2π(αx)(1θx)αexp[12(1θx)2α],x>0.

    For bias expressions (not previously reported in the literature, see the “analyticalBC.R” file.

    Using the data set from , we have n=20, α^=3.0869, θ^=0.6731, se^(α^)=0.5534 and se^(θ^)=0.0379. Evaluating the analytical expressions and the coxsnell.bc() function, we have, respectively,

    invgenhalfnormal.bc(n = 20, mle = c(3.0869, 0.6731))
    ##     alpha     theta
    ##  0.229016 -0.002953
    pdf <- quote(sqrt(2) * pi^(-0.5) * alpha * x^(-alpha - 1) *
                 exp(-0.5 * x^(-2 * alpha) * (1 / theta)^(2 * alpha)) * 
                 theta^(-alpha))
    lpdf <- quote(log(alpha) - alpha  * log(x) - 0.5e0 / (x^alpha)^2*
                  theta^(-2 * alpha) - alpha * log(theta))
    coxsnell.bc(density = pdf, logdensity = lpdf, n = 20,
                parms = c("alpha", "theta"), mle = c(3.0869, 0.6731),
                lower = 0)$bias
    ##     alpha     theta
    ##  0.229016 -0.002953
  26. Marshall-Olkin extended exponential distribution with shape α and rate λ f(xα,λ)=λαexp(λx)[1(1α)exp(λx)]2,x>0.

    For bias expressions (not previously reported in the literature, see the “analyticalBC.R” file.

    Using the data set from , we have n=20, α^=0.2782, λ^=0.0078, se^(α^)=0.2321 and se^(λ^)=0.0049. Evaluating the analytical expressions and the coxsnell.bc() function, we have, respectively,

    moeexp.bc(n = 20, mle = c(0.2781, 0.0078))
    ##    alpha   lambda
    ## 0.210919 0.003741
    pdf <- quote(alpha * lambda * exp(-x * lambda) /
                 ((1- (1 - alpha) * exp(- x * lambda)))^2)
    lpdf <- quote(log(alpha) + log(lambda) - x * lambda -
                  2 * log((1 - (1-alpha) * exp(- x * lambda))))
    coxsnell.bc(density = pdf, logdensity = lpdf, n = 20,
                parms = c("alpha", "lambda"), mle = c(0.2781, 0.0078),
                lower = 0)$bias
    ##   alpha  lambda
    ## 0.21086 0.00374
  27. Beta distribution with shapes α and β f(xα,β)=Γ(α+β)Γ(α)Γ(β)xα1(1x)β1,0<x<1.

    For bias expressions, see .

    Using the data set from , we have n=48, α^=5.941, β^=21.2024, se^(α^)=1.1812 and se^(β^)=4.3462. Evaluating the analytical expressions in , our analytical expressions and the coxsnell.bc() function, we have, respectively,

    beta.gauss.bc(n = 48, mle = c(5.941, 21.2024))
    ##  alpha   beta
    ## -4.784 -4.125
    beta.bc(n = 48, mle = c(5.941, 21.2024))
    ##  alpha   beta
    ## 0.3582 1.3315
    pdf <- quote(gamma(alpha + beta) / (gamma(alpha) * gamma(beta)) *
                 x^(alpha - 1) * (1 - x)^(beta - 1))
    lpdf <- quote(lgamma(alpha + beta) - lgamma(alpha) -
                  lgamma(beta) + alpha * log(x) + beta * log(1 - x))
    coxsnell.bc(density = pdf, logdensity = lpdf, n = 48,
                parms = c("alpha", "beta"), mle = c(5.941, 21.2024),
                lower = 0, upper =  1)$bias
    ##  alpha   beta
    ## 0.3582 1.3315
  28. Kumaraswamy distribution with shapes α and β f(xα,β)=αβxα1(1xα)β1,0<x<1.

    For bias expressions, see .

    Using the data set from , we have n=20, α^=6.3478, β^=4.4898, se^(α^)=1.5576 and se^(β^)=2.0414. Evaluating the analytical expressions and the coxsnell.bc() function, we have, respectively,

    kum.bc(n = 20, mle = c(6.3478, 4.4898))
    ##   alpha    beta
    ##  -6.573 -13.323
    pdf <- quote(alpha * beta * x^(alpha - 1) *
                 (1 - x^alpha)^(beta - 1))
    lpdf <- quote(log(alpha) + log(beta) + alpha * log(x) + (beta - 1) *
                  log(1 - x^alpha))
    coxsnell.bc(density = pdf, logdensity = lpdf, n = 20,
                parms = c("alpha", "beta"), mle = c(6.3478, 4.4898),
                lower = 0, upper = 1)$bias
    ## alpha  beta
    ## 0.514 1.013
  29. Inverse beta distribution with shapes α and β f(xα,β)=Γ(α+β)Γ(α)Γ(β)xα1(1+x)(α+β),x>0.

    For bias expressions, see .

    Using the data set from , we have n=116, α^=28.5719, β^=1.3783, se^(α^)=4.0367 and se^(β^)=0.1637. Evaluating the analytical expressions and the coxsnell.bc() function, we have, respectively,

    invbeta.bc(n = 116, mle = c(28.5719, 1.3782))
    ##  alpha   beta
    ## 534.26  17.73
    pdf <- quote(gamma(alpha + beta) * x^(alpha - 1) *
                 (1 + x)^(- alpha - beta) / gamma(alpha)/gamma(beta))
    lpdf <- quote(lgamma(alpha + beta) + alpha * log(x) - 
                  (alpha + beta) * log(1 + x) - lgamma(alpha) - lgamma(beta))
    coxsnell.bc(density = pdf, logdensity = lpdf, n = 116,
                parms = c("alpha", "beta"), mle = c(28.5719, 1.3782),
                lower = 0)$bias
    ##  alpha   beta
    ## 0.8025 0.0306
  30. Birnbaum-Saunders distribution with shape α and scale β f(xα,β)=12αβ2π[(βx)1/2+(βx)3/2]exp[12,α2(xβ+βx2)],x>0.

    Bias expressions : (31)B(α^)=α4n(1+2+α2α(2π)1/2h(α)+1) and (32)B(β^)=β2α22n[α(2π)1/2h(α)+1], where h(α)=απ2πe2/α2[1Φ(2α)].

    Using the data set from , we have n=20, α^=0.3149, β^=1.8105, se^(α^)=0.0498 and se^(β^)=0.1259. Evaluating the analytical expressions (31), (32) and the coxsnell.bc() function, we have, respectively,

    birnbaumsaunders.bc(n = 20, mle = c(0.3148, 1.8104))
    ##     alpha      beta
    ## -0.011991  0.004374
    pdf <- quote(1 / (2 * alpha * beta * sqrt(2 * pi)) *
                 ((beta / x)^0.5 + (beta / x)^1.5) *
                 exp(- 1/(2 * alpha^2) * (x / beta + beta/ x - 2)))
    lpdf <- quote(-log(alpha) - log(beta) - 1 / (2 * alpha^2) *
                  (x / beta + beta/ x - 2) + log((beta / x)^0.5 + 
                  (beta / x)^1.5))
    coxsnell.bc(density = pdf, logdensity = lpdf, n = 20,
                parms = c("alpha", "beta"), mle = c(0.3148, 1.8104), 
                lower = 0)$bias
    ##     alpha      beta
    ## -0.011991  0.004374
  31. Generalized Pareto distribution with shape ξ and scale σ f(xξ,σ)=1σ(1+ξxσ)(1/ξ+1),x>0, ξ0.

    Bias expressions : (33)B(ξ^)=(1+ξ)(3+ξ)n(1+3ξ) and (34)B(σ^)=σ(3+5ξ+4ξ2)n(1+3ξ).

    Using the data set from , we have n=58, ξ^=0.736, σ^=1.709, se^(ξ^)=0.223 and se^(σ^)=0.41. Evaluating the analytical expressions (33), (34) and the coxsnell.bc() function, we have, respectively,

    genpareto.bc(n = 58, mle = c(0.736, 1.709))
    ##       xi    sigma
    ## -0.03486  0.08126
    pdf <- quote(1 / sigma * (1 + xi * x / sigma )^(-(1 + 1 / xi)))
    Rlpdf <- quote(-log(sigma) - (1 + 1 / xi) * log(1 + xi * x / sigma))
    coxsnell.bc(density = pdf, logdensity = lpdf, n = 58,
                parms = c("xi", "sigma"), mle = c(0.736, 1.709), 
                lower = 0)$bias
    ##       xi    sigma
    ## -0.03486  0.08126

5 Additional Applications

In this section, we present additional numerical results returned by cosnell.bc(), observed.varc() and expected.varcov(). For the data describing the times between successive electric pulses on the surface of isolated muscle fiber , we fitted the exponentiated Weibull, Marshall-Olkin extended Weibull, Weibull, Marshall-Olkin extended exponential and exponential distributions. These distributions were also fitted by . There are 799 observations and for each distribution we report the MLEs, the bias corrected MLEs, the observed variance-covariance obtained from the numerical Hessian H11(θ^), the observed variance-covariance obtained from the analytical Hessian H21(θ^), the expected variance-covariance I1(θ^) and the expected variance-covariance evaluated at the bias corrected MLEs I1(θ~). The MLEs and the H11(θ^) matrix were obtained by the fitdistrplus package . The R codes used to obtain the numerical results are available in the supplementary material.

It is important to emphasize that for the Marshall-Olkin extended Weibull and exponentiated Weibull distributions, it is not possible to obtain analytical expressions for bias corrections. The exponentiated-Weibull family was proposed by . Its probability density function is: f(xλ,β,α)=αβλxβ1eλxβ(1eλxβ)α1, where λ>0 is the scale parameter and β>0 and α>0 are the shape parameters. The Marshall-Olkin extended Weibull distribution was introduced by . Its probability density function is: f(xλ,β,α)=αβλxβ1eλxβ(1αeλxβ)2, where λ>0 is the scale parameter, β>0 is the shape parameter, α>0 is an additional shape parameter and α=1α.

The fitted parameter estimates and their bias corrected estimates are shown in Table 1. We see that the bias corrected MLEs for α and λ of the MOE-Weibull and exp-Weibull distributions are quite different from the original MLEs.

Table 1: MLEs and bias corrected MLEs.
Distribution α^ β^ λ^ α~ β~ λ~
MOE-Weibull 0.3460 1.3247 0.0203 0.3283 1.3240 0.0188
exp-Weibull 1.9396 0.7677 0.2527 1.8973 0.7625 0.2461
Weibull 1.0829 0.0723 1.0811 0.0723
MOE-exponential 1.1966 0.0998 1.1820 0.0994
exponential 0.0913 0.0912

It is important to assess the accuracy of MLEs. The two common ways for this are through the inverse observed Fisher information and the inverse expected Fisher information matrices. The results below show large differences between the observed H1 and expected I1 information matrices. As demonstrated by , the I1 outperforms the H1 under a mean squared error criterion, hence with mle.tools the researchers may choose one of them and not use the easier. Furthermore, in general, we observe that the bias corrected MLEs decrease the variance of estimates.

Exponentiated Weibull distribution: H11(θ^)=[0.007260.007170.035640.007170.007180.034930.035640.034930.18045],H21(θ^)=[0.007290.007200.035790.007200.007210.035090.035790.035090.18120],I1(θ^)=[0.005320.005240.026090.005240.005270.025450.026090.025450.13333],I1(θ~)=[0.005100.005100.024820.005100.005190.024540.024820.024540.12590].

Marshall-Olkin extended Weibull distribution: H11(θ^)=[0.000040.000360.000520.000360.003610.004300.000520.004300.00748],H21(θ^)=[0.000050.000470.000680.000470.004680.005820.000680.005820.00967],I1(θ^)=[0.000060.000560.000820.000560.005420.006990.000820.006990.01146],I1(θ~)=[0.000050.000510.000720.000510.005260.006510.000720.006510.01030].

Weibull distribution: H11(θ^)=[0.000040.000180.000180.00086],H21(θ^)=[0.000040.000180.000180.00087],I1(θ^)=[0.000040.000180.000180.00089],I1(θ~)=[0.000040.000180.000180.00089].

Marshall-Olkin extended exponential distribution: H11(θ^)=[0.000040.000810.000810.02022],H21(θ^)=[0.000040.000810.000810.02023],I1(θ^)=[0.000040.000830.000830.02094],I1(θ~)=[0.000040.000820.000820.02047].

Exponential distribution: H11(θ^)=0.000010433,H21(θ^)=0.000010436,I1(θ^)=0.000010436,I1(θ~)=0.000010410.

6 Concluding Remarks

As pointed out by several works in the literature, the Cox-Snell methodology, in general, is efficient for reducing the bias of the MLEs. However, the analytical expressions are either notoriously cumbersome or even impossible to deduce. To the best of our knowledge, there are only two alternatives to obtain the analytical expressions automatically, those presented in and . They use the commercial softwares Maple and Mathematica .

In order to calculate the bias corrected estimates in a simple way, developed an R package, uploaded to CRAN on 2 February, 2017. Its main function, coxsnell.bc(), evaluates the bias corrected estimates. The usefulness of this function has been tested for thirty one continuous probability distributions. Bias expressions, for most of them, are available in the literature.

It is well known that the Fisher information can be computed using the first or second order derivatives of the log-likelihood function. In our implementation, the functions expected.varcov() and coxsnell.bc() are using the second order derivatives, analytically returned by the D() function. In a future work, we intend to check if there is any gain in calculating the Fisher information from the first order derivatives of the log-hazard rate function or from the first order derivatives of the log-reversed-hazard rate function. showed that the Fisher information can be computed using the hazard rate function. computed the Fisher information from the first order derivatives of the log-reversed-hazard rate function. In general, expressions of the first order derivatives of the log-hazard rate function (log-reversed-hazard rate function) are simpler than second order derivatives of the log-likelihood function. In this sense, the integrate() function can work better. It is important to point out that the hazard rate function and the reversed hazard rate function are given, respectively, by h(xθ)=ddxlog[S(xθ)] and h(xθ)=ddxlog[F(xθ)], where S(xθ) and F(xθ) are, respectively, the survival function and the cumulative distribution function.

In the next version of mle.tools, we will include, using analytical first and second-order partial derivatives, the following:

where L is the value of the likelihood function evaluated at the MLEs, n is the number of observations, and p is the number of estimated parameters.

Also, the next version of the package will incorporate analytical expressions for the distributions studied in Section 4 implemented in the supplementary file “analyticalBC.R”.

CRAN packages used

mle.tools, fitdistrplus

CRAN Task Views implied by cited packages

ActuarialScience, Distributions, Survival

Note

This article is converted from a Legacy LaTeX article using the texor package. The pdf version is the official version. To report a problem with the html, refer to CONTRIBUTE on the R Journal homepage.

Footnotes

    References

    J. A. Achcar, S. R. Lopes, J. Mazucheli and R. R. Linhares. A Bayesian approach for stable distributions: Some computational aspects. Open Journal of Statistics, 3(04): 268, 2013.
    A. Alzaatreh, M. Mansoor, M. Tahir, M. Zubair and S. Ali. The gamma half-Cauchy distribution: Properties and applications. HACETTEPE Journal of Mathematics and Statistics, 45(4): 1143–1159, 2016.
    M. Bader and A. Priest. Statistical aspects of fibre and bundle strength in hybrid composites. Progress in Science and Engineering of Composites, 1129–1136, 1982.
    D. K. Bhaumik, K. Kapur and R. D. Gibbons. Testing parameters of a gamma distribution for small samples. Technometrics, 51(3): 326–334, 2009.
    X. Cao. Relative performance of expected and observed Fisher information in covariance estimation for maximum likelihood estimates. 2013.
    R. Chhikara and J. Folks. The inverse Gaussian distribution as a lifetime model. Technometrics, 19(4): 461–468, 1977.
    G. M. Cordeiro and F. Cribari-Neto. An introduction to Bartlett correction and bias reduction. New York: Springer-Verlag, 2014.
    G. M. Cordeiro, E. C. Da Rocha, J. G. C. Da Rocha and F. Cribari-Neto. Bias-corrected maximum likelihood estimation for the beta distribution. Journal of Statistical Computation and Simulation, 58(1): 21–35, 1997.
    G. M. Cordeiro and R. dos Santos Brito. The Beta power distribution. Brazilian Journal of Probability and Statistics, 26(1): 88–112, 2012.
    G. M. Cordeiro and R. Klein. Bias correction in ARMA models. Statistics & Probability Letters, 19(3): 169–176, 1994.
    G. M. Cordeiro and A. J. Lemonte. On the Marshall-Olkin extended Weibull distribution. Statistical Papers, 54(2): 333–353, 2013.
    D. R. Cox and D. V. Hinkley. Theoretical statistics. CRC Press, 1979.
    D. R. Cox and P. A. W. Lewis. The statistical analysis of series of events. Netherlands: Springer-Verlag, 1966.
    D. R. Cox and E. J. Snell. A general definition of residuals. Journal of the Royal Statistical Society B, 30(2): 248–275, 1968.
    F. Cribari-Neto and K. L. P. Vasconcellos. Nearly unbiased maximum likelihood estimation for the beta distribution. Journal of Statistical Computation and Simulation, 72(2): 107–118, 2002.
    D. Datta and D. Datta. Comparison of Weibull distribution and exponentiated Weibull distribution based estimation of mean and variance of wind data. International Journal of Energy, Information and Communications, 4(4): 1–11, 2013.
    M. L. Delignette-Muller, C. Dutang and A. Siberchicot. : Help to fit of a parametric distribution to non-censored or censored data. 2017. R package version 1.0-8.
    M. Delignette-Muller, M. Cornu, A. S. S. Group, et al. Quantitative risk assessment for Escherichia coli O157: H7 in frozen ground beef patties consumed by young children in french households. International Journal of Food Microbiology, 128(1): 158–164, 2008.
    A. W. F. Edwards. Likelihood (expanded edition). Baltimore: Johns Hopkins University Press, 1992.
    B. Efron and I. M. Johnstone. Fisher’s information in terms of the hazard rate. The Annals of Statistics, 18(1): 38–62, 1990.
    M. E. Ghitany, D. K. Al-Mutairi, N. Balakrishnan and L. J. Al-Enezi. Power Lindley distribution and associated inference. Computational Statistics & Data Analysis, 64: 20–33, 2013.
    M. E. Ghitany, B. Atieh and S. Nadarajah. Lindley distribution and its application. Mathematics and Computers in Simulation, 78(4): 493–506, 2008.
    D. E. Giles. A note on improved estimation for the Topp-Leone distribution. Department of Economics, University of Victoria. 2012a.
    D. E. Giles. Bias reduction for the maximum likelihood estimators of the parameters in the half-logistic distribution. Communication in Statistics - Theory and Methods, 41(2): 212–222, 2012b.
    D. E. Giles and H. Feng. Bias of the maximum likelihood estimators of the two-parameter gamma distribution revisited. Department of Economics, University of Victoria. 2009.
    D. E. Giles, H. Feng and R. T. Godwin. Bias-corrected maximum likelihood estimation of the parameters of the generalized Pareto distribution. Communications in Statistics – Theory and Methods, 45(8): 2465–2483, 2016.
    D. E. Giles, H. Feng and R. T. Godwin. On the bias of the maximum likelihood estimator for the two-parameter Lomax distribution. Communications in Statistics - Theory and Methods, 42(11): 1934–1950, 2013.
    A. E. Gomes, C. Q. da-Silva, G. M. Cordeiro and E. M. M. Ortega. A new lifetime model: The Kumaraswamy generalized Rayleigh distribution. Journal of Statistical Computation and Simulation, 84(2): 290–309, 2014.
    A. J. Gross and V. A. Clark. Survival distributions: Reliability applications in the biometrical sciences. John Wiley & Sons, 1976.
    R. D. Gupta, R. C. Gupta and P. G. Sankaran. Some characterization results based on factorization of the (reversed) hazard rate function. Communication in Statistics - Theory and Methods, 33: 2004.
    Z. Javanshiri, A. Habibi Rad and N. R. Arghami. Exp-Kumaraswamy distributions: Some properties and applications. Journal of Sciences, Islamic Republic of Iran, 26(1): 57–69, 2015.
    P. H. Johnson, Y. Qi and Y. C. Chueh. Bias-corrected maximum likelihood estimation in actuarial science. Urbana-Champaign, Champaign, IL: University of Illinois. 2012a.
    P. H. Johnson, Y. Qi and Y. C. Chueh. CSCK MLE bias calculation. Urbana-Champaign, Champaign, IL: University of Illinois. 2012b.
    B. Jørgensen. Statistical properties of the generalized inverse gaussian distribution. Springer-Verlag, 1982.
    S. Kay. Asymptotic maximum likelihood estimator performance for chaotic signals in noise. IEEE Transactions on Signal Processing, 43(4): 1009–1012, 1995.
    S. Kumagai and I. Matsunaga. Changes in the distribution of short-term exposure concentration with different averaging times. American Industrial Hygiene Association, 56(1): 24–31, 1995.
    S. Kumagai, I. Matsunaga, K. Sugimoto, Y. Kusaka and T. Shirakawa. Assessment of occupational exposures to industrial hazardous substances. III. On the frequency distribution of daily exposure averages (8-h TWA). Japanese Journal of Industrial Health, 31(4): 216–226, 1989.
    D. Kundu. Discriminating between normal and Laplace distributions. In Advances in ranking and selection, multiple comparisons, and reliability, Eds N. Balakrishnan, N. Kannan and H. N. Nagaraja pages. 65–79 2005. Basel: Birkhäuser.
    B. Lagos-Álvarez, M. D. Jiménez-Gamero and V. Alba-Fernández. Bias correction in the type I generalized logistic distribution. Communications in Statistics - Simulation and Computation, 40(4): 511–531, 2011.
    J. F. Lawless. Statistical Models and Methods for Lifetime Data, second edition. John Wiley & Sons, 2011.
    E. L. Lehmann. Elements of large-sample theory. Springer-Verlag, 1999.
    A. J. Lemonte. Improved point estimation for the Kumaraswamy distribution. Journal of Statistical Computation and Simulation, 81(12): 1971–1982, 2011.
    A. J. Lemonte, F. Cribari-Neto and K. L. P. Vasconcellos. Improved statistical inference for the two-parameter Birnbaum-Saunders distribution. Computational Statistics & Data Analysis, 51(9): 4656–4681, 2007.
    H. Linhart and W. Zucchini. Model selection. John Wiley & Sons, 1986.
    Maple. Maplesoft, a division of waterloo maple inc. Ontario: Waterloo, 2017.
    A. W. Marshall and I. Olkin. A new method for adding a parameter to a family of distributions with application to the exponential and Weibull families. Biometrika, 84(3): 641–652, 1997.
    J. Mazucheli. Mle.tools: Expected/Observed Fisher information and bias-corrected maximum likelihood estimate(s). 2017. R package version 1.0.0.
    J. Mazucheli and S. Dey. Bias-corrected maximum likelihood estimation of the parameters of the generalized half-normal distribution. submitted to Journal of Statistical Computation and Simulation, 2017.
    G. S. Mudholkar and D. K. Srivastava. Exponentiated Weibull family for analyzing bathtub failure-rate data. IEEE Transactions on Reliability, 42(2): 299–302, 1993.
    S. Nadarajah. A truncated inverted beta distribution with application to air pollution data. Stochastic Environmental Research and Risk Assessment, 22(2): 285–289, 2008a.
    S. Nadarajah. The model for fracture toughness. Journal of Mechanical Science and Technology, 22: 1255–1258, 2008b.
    S. Nadarajah, H. S. Bakouch and R. Tahmasbi. A generalized Lindley distribution. Sankhya B, 73(2): 331–359, 2011.
    M. D. Nichols and W. Padgett. A bootstrap control chart for Weibull percentiles. Quality and reliability engineering international, 22(2): 141–151, 2006.
    R Core Team. R: A language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing, 2016. ISBN 3-900051-07-0.
    M. Z. Raqab, M. T. Madi and D. Kundu. Estimation of P(Y < X) for the three-parameter generalized exponential distribution. Communications in Statistics - Theory and Methods, 37(18): 2854–2864, 2008.
    J. Reath. Improved parameter estimation of the log-logistic distribution with applications. 2016.
    T. Ross and N. Lott. A climatology of 1980 – 2003 extreme weather and climate events. 01. National Climatic Data Center. 2003.
    K. Saha and S. Paul. Bias-corrected maximum likelihood estimator of the negative binomial dispersion parameter. Biometrics, 61(1): 179–185, 2005.
    J. Schwartz and D. E. Giles. Bias-reduced maximum likelihood estimation of the zero-inflated Poisson distribution. Communications in Statistics - Theory and Methods, 45(2): 465–478, 2016.
    J. Schwartz, R. T. Godwin and D. E. Giles. Improved maximum-likelihood estimation of the shape parameter in the Nakagami distribution. Journal of Statistical Computation and Simulation, 83(3): 434–445, 2013.
    R. Shanker. Shanker distribution and its applications. International Journal of Statistics and Applications, 5(6): 338–348, 2015.
    V. K. Sharma, S. K. Singh, U. Singh and V. Agiwal. The inverse Lindley distribution: A stress-strength reliability model with application to head and neck cancer data. Journal of Industrial and Production Engineering, 32(3): 162–173, 2015.
    B. D. Stočsić and G. M. Cordeiro. Using Maple and Mathematica to derive bias corrections for two parameter distributions. Journal of Statistical Computation and Simulation, 79(6): 751–767, 2009.
    M. Tahir, M. A. Hussain, G. M. Cordeiro, G. Hamedani, M. Mansoor and M. Zubair. The Gumbel-Lomax distribution: Properties and applications. Journal of Statistical Theory and Applications, 15(1): 61–79, 2016.
    M. Teimouri and S. Nadarajah. Bias corrected MLEs for the Weibull distribution based on records. Statistical Methodology, 13: 12–24, 2013.
    M. Teimouri and S. Nadarajah. Bias corrected MLEs under progressive type-II censoring scheme. Journal of Statistical Computation and Simulation, 86(14): 2714–2726, 2016.
    B. X. Wang, X. K. Wang and K. Yu. Inference on the Kumaraswamy distribution. Communications in Statistics - Theory and Methods, 46(5): 2079–2090, 2017.
    M. Wang and W. Wang. Bias-corrected maximum likelihood estimation of the parameters of the weighted Lindley distribution. Communications in Statistics - Theory and Methods, 46(1): 530–545, 2017.
    Wentao. Wang. Bias-corrected maximum likelihood estimation of the parameters of the weighted Lindley distribution. 2015.
    Wolfram Research, Inc. Mathematica, version 8.0. Vienna, Champaign, Illinois: Wolfram Research, Inc., 2010.
    L. Xiao and D. E. Giles. Bias reduction for the maximum likelihood estimator of the parameters of the generalized Rayleigh family of distributions. Communications in Statistics - Theory and Methods, 43(8): 1778–1792, 2014.
    G. Zhang and R. Liu. Bias-corrected estimators of scalar skew normal. Communications in Statistics - Simulation and Computation, 46(2): 831–839, 2015.

    Reuse

    Text and figures are licensed under Creative Commons Attribution CC BY 4.0. The figures that have been reused from other sources don't fall under this license and can be recognized by a note in their caption: "Figure from ...".

    Citation

    For attribution, please cite this work as

    Mazucheli, et al., "mle.tools: An R Package for Maximum Likelihood Bias Correction", The R Journal, 2017

    BibTeX citation

    @article{RJ-2017-055,
      author = {Mazucheli, Josmar and Menezes, André Felipe Berdusco and Nadarajah, Saralees},
      title = {mle.tools: An R Package for Maximum Likelihood Bias Correction},
      journal = {The R Journal},
      year = {2017},
      note = {https://rjournal.github.io/},
      volume = {9},
      issue = {2},
      issn = {2073-4859},
      pages = {268-290}
    }