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

Recently, (Mazucheli 2017) uploaded the package mle.tools to CRAN. It can be used for bias corrections of maximum likelihood estimates through the methodology proposed by (Cox and Snell 1968). 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 (Cox and Snell 1968) 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.

Josmar Mazucheli (Department of Statistics, Universidade Estadual de Maringá) , André Felipe Berdusco Menezes (Department of Statistics, Universidade Estadual de Maringá) , Saralees Nadarajah (School of Mathematics, University of Manchester)
2017-11-01

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 (Edwards 1992; Lehmann 1999). 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 (Kay 1995). 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 \(\mathcal{O}\left(n^{-1}\right)\), thus these errors reduce as sample size increases (Cordeiro and Cribari-Neto 2014).

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 (Cordeiro et al. 1997), (Cribari-Neto and Vasconcellos 2002), (Saha and Paul 2005), (Lemonte et al. 2007), (Giles and Feng 2009) (Lagos-Álvarez et al. 2011), (Lemonte 2011), (Giles 2012b), (Giles 2012a), (Schwartz et al. 2013), (Giles et al. 2013), (Teimouri and Nadarajah 2013), (Xiao and Giles 2014), (Zhang and Liu 2015), (Teimouri and Nadarajah 2016), (Reath 2016), (Giles et al. 2016), (Schwartz and Giles 2016), (Wang and Wang 2017), (Mazucheli and Dey 2017) 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. (Stočsić and Cordeiro 2009) 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” (Johnson et al. 2012b) 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 (Maple 2017) and Mathematica (Wolfram Research, Inc. 2010), are commercial softwares.

In this paper, our objective is to introduce a new contributed R (R Core Team 2016) package, namely mle.tools that computes the expected/observed Fisher information and the bias corrected estimates by the methodology proposed by (Cox and Snell 1968). 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 \(X_1, \ldots, X_n\) be \(n\) be independent random variables with probability density function \(f\left( x_i \mid \mathbf{\theta} \right)\) depending on a \(p\)-dimensional parameter vector \(\mathbf{\theta} = \left(\theta_1, \ldots, \theta_p\right)\). Without loss of generality, let \(l =l\left(\mathbf{\theta} \mid \mathbf{x}\right)\) be the log-likelihood function for the unknown \(p\)-dimensional parameter vector \(\mathbf{\theta}\) given a sample of \(n\) observations. We shall assume some regularity conditions on the behavior of \(l\left(\mathbf{\theta} \mid \mathbf{x}\right)\) (Cox and Hinkley 1979).

The joint cumulants of the derivatives of \(l\) are given by: \[\begin{aligned} \kappa_{ij} &=& \mathbb{E}\left[\dfrac {\partial^2\, l}{\partial\,\mathbf{\theta}_i\, \partial\,\mathbf{\theta}_j} \right], \\ \nonumber \\ \kappa_{ijl} &=& \mathbb{E}\left[\dfrac {\partial^3\, l}{\partial\,\mathbf{\theta}_i\, \partial\,\mathbf{\theta}_j\, \partial\,\mathbf{\theta}_l} \right], \\ \nonumber \\ \kappa_{ij,l} &=& \mathbb{E}\left[\left(\dfrac {\partial^2\, l}{\partial\,\mathbf{\theta}_i\, \partial\,\mathbf{\theta}_j}\right)\,\left(\dfrac {\partial\, l}{\partial\,\mathbf{\theta}_l}\right)\right], \\ \nonumber \\ \kappa_{ij}^{(l)} &=& \dfrac {\partial\,\kappa_{ij}}{\partial\,\mathbf{\theta}_l} \end{aligned}\] for \(i, j, l = 1, \ldots, p\).

The bias expression of the \(s\)th element of \(\widehat{\mathbf{\theta}}\), the MLEs of \(\mathbf{\theta}\), when the sample data are independent, but not necessarily identically distributed, was proposed by (Cox and Snell 1968): \[\begin{aligned} \label{eq:coxsnell} \mathcal{B}\left(\widehat{\theta}_s\right) = \sum_{i=1}^{p}\,\sum_{j=1}^{p}\,\sum_{l=1}^{p}\,\kappa^{si}\,\kappa^{jl}\,\left[0.5 \kappa_{ijl} + \kappa_{ij,l}\right] + \mathcal{O} \left(n^{-2}\right), \end{aligned} \tag{1}\] where \(s = 1, \ldots, p\) and \(\kappa^{ij}\) is the \((i, j)\)th element of the inverse of the negative of the expected Fisher information.

Thereafter, (Cordeiro and Klein 1994) noticed that equation (1) holds even if the data are non-independent, and it can be re-expressed as: \[\begin{aligned} \label{eq:cordeiro-klein} \mathcal{B} \left(\widehat{\theta}_s\right) = \sum_{i=1}^{p}\,\kappa^{si}\sum_{j=1}^{p}\,\sum_{l=1}^{p}\left[\kappa_{ij}^{(l)} - 0.5 \kappa_{ijl}\right]\,\kappa^{jl} + \mathcal{O} \left(n^{-2}\right). \end{aligned} \tag{2}\]

Defining \(a_{ij}^{(l)} = \kappa_{ij}^{(l)} - 0.5 \kappa_{ijl}\), \(A^{(l)} = \left\{a_{ij}^{(l)}\right\}\) and \(K = \left[ -\kappa_{ij}\right]\), the expected Fisher information matrix for \(i, j, l = 1, \ldots, n\), the bias expression for \(\widehat{\mathbf{\theta}}\) in matrix notation is: \[\begin{aligned} \mathcal{B}\left(\widehat{\mathbf{\theta}}\right) = K^{-1} A \text{vec}\left(K^{-1}\right) + \mathcal{O}\left(n^{-2}\right), \end{aligned}\] where \(\text{vec} \left(K^{-1}\right)\) is the vector obtained by stacking the columns of \(K^{-1}\) and \(A = \left\{A^{1} \mid \cdots \mid A^{p}\right\}\).

Finally, the bias corrected MLE for \(\theta_s\) can be obtained as: \[\begin{aligned} \label{eq:mle-coxsnell} \widetilde{\theta}_s = \widehat{\theta}_s - \widehat{\mathcal{B}}\left(\widehat{\theta}_s\right). \end{aligned} \tag{3}\] Alternatively, using matrix notation the bias corrected MLEs can be expressed as (Cordeiro and Klein 1994): \[\begin{aligned} \label{eq:mle-cordeiro-klein} \widetilde{\mathbf{\theta}} = \widehat{\mathbf{\theta}} - \widehat{K}^{-1} \widehat{A} \text{vec} \left(\widehat{K}^{-1}\right), \end{aligned} \tag{4}\] where \(\widehat{K} = K\big|_{\mathbf{\theta}=\widehat{\mathbf{\theta}}}\) and \(\widehat{A} = A\big|_{\mathbf{\theta}=\widehat{\mathbf{\theta}}}\).

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 \(X_{1}, \ldots, X_{n}\) be independent and identical random variables with probability density function \(f\left(x_{i}\mid \mathbf{\theta}\right)\) depending on a \(p\)-dimensional parameter vector \(\mathbf{\theta} = \left( \theta_1,\ldots,\theta_p\right)\). The \((j,k)\)th element of the observed, \(H_{jk}\), and expected, \(I_{jk}\), Fisher information are calculated, respectively, as \[\begin{aligned} H_{jk} =\left. {-\sum\limits_{i=1}^{n}\frac {\partial^{2}}{\partial \theta_{j}\partial \theta_{k}}\log f\left(x_{i}\mid {\mathbf{\theta} }\right) } \right\vert_{\mathbf{\theta }=\widehat{\mathbf{\theta}}} \end{aligned}\] and \[\begin{aligned} I_{jk}=-n\times E\left( \frac {\partial^{2}}{\partial \theta_{j}\partial \theta_{k}}\log f\left( x\mid \mathbf{\theta }\right) \right) = \left.-n \times\int\limits_{\mathcal{X}} \frac {\partial^2}{\partial\theta_j\partial \theta_k}\log f\left(x\mid \mathbf{\theta}\right) \times f\left(x\mid \mathbf{\theta }\right) \textrm{d}x\right\vert_{\mathbf{\theta }=\widehat{\mathbf{\theta}}}, \end{aligned}\] where \(j,k = 1,\ldots,p\), \(\mathbf{\widehat{\theta}}\) is the MLE of \(\mathbf{\theta}\) and \(\mathcal{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 \(\theta_s\), \(s=1,\ldots, p\) denoted by \(\widetilde{\theta_s}\) is calculated as \(\widetilde{\theta_s} = \widehat{\theta}_s - \widehat{\mathcal{B}} \left( \widehat{\theta}_s \right)\), where \(\widehat{\theta}_s\) is the MLE of \({\theta}_s\) and \[\begin{aligned} {\widehat{\mathcal{B}}\left({\widehat{\theta }}_{s}\right) = } \left. {\sum\limits_{j=1}^{p}\sum\limits_{k=1}^{p} \sum\limits_{l=1}^{p}\kappa^{sj}\kappa^{kl}\left[ 0.5\kappa_{{jkl}}+\kappa_{{jk,l}}\right]} \right\vert_{\mathbf{\theta }=\widehat{\mathbf{\theta }}}, \end{aligned}\] where \(\kappa^{jk}\) is the (\(j,k\))th element of the inverse of the negative of the expected Fisher information, \[\begin{aligned} {\kappa_{jkl}=} \left.n \int\limits_{\mathcal{X}} \frac {\partial^3}{\partial \theta_j \partial \theta_k \partial \theta_l} \log f\left(x\mid \mathbf{\theta}\right) f\left(x\mid \mathbf{\theta }\right) \textrm{d}x\right\vert_{\mathbf{\theta }=\widehat{\mathbf{\theta}}}, \end{aligned}\]

\[\begin{aligned} \kappa_{jk,l}= \left.n \int\limits_{\mathcal{X}} \frac {\partial^2}{\partial\theta_j\partial \theta_k} \log f\left(x\mid \mathbf{\theta}\right) \frac {\partial }{{\theta }_{l}}\log f\left( x\mid\mathbf{\theta }\right) f\left(x\mid \mathbf{\theta }\right) \textrm{d}x\right\vert_{\mathbf{\theta }=\widehat{\mathbf{\theta}}} \end{aligned}\] and \(\mathcal{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 (Cordeiro et al. 1997), (Lemonte 2011) and (Stočsić and Cordeiro 2009). 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 (Cordeiro et al. 1997) and (Stočsić and Cordeiro 2009). 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 (Maple 2017) and the results are exactly equals.

  1. One-parameter Lindley distribution with scale parameter \(\theta\) \[\begin{aligned} f(x\mid\theta) = \frac {\theta^2}{1 + \theta} (1 + x) \exp(-\theta x), \quad x > 0. \end{aligned}\]

    \(\bullet\) Bias expression (not previously reported in the literature): \[\begin{aligned} \label{bc-lindley} \mathcal{B}\left(\widehat{\theta}\right) = {\frac { \left( {{\theta}}^{3}+ 6\,{{\theta}}^{2}+ 6\,{\theta}+ 2 \right) \left( {\theta}+ 1\right) {\theta}}{n \left( {{\theta}}^{2}+ 4\,{\theta}+ 2 \right)^{2}}}. \end{aligned} \tag{5}\]

    Using the data set from (Ghitany et al. 2008) we have \(n = 100\), \(\widehat{\theta} = 0.1866\) and \(\widehat{se}\left(\widehat{\theta}\right)= 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 \(\theta\) \[\begin{aligned} f(x\mid\theta) = \frac {\theta^2}{1 + \theta} \left(\frac {1 + x}{x^3}\right)\exp\left(-\frac {\theta}{x}\right), \quad x > 0. \end{aligned}\]

    \(\bullet\) Bias expression (Wang 2015): \[\begin{aligned} \label{bc-invlindley} \mathcal{B}\left(\widehat{\theta}\right) ={\frac { \left( {\theta}+ 1 \right) {\theta}\, \left( {{\theta}}^{3}+ 6\,{{\theta}}^{2}+ 6\,{\theta}+ 2 \right) }{n\left( {{\theta}}^{2}+ 4\,{\theta}+ 2\right)^{2}}}. \end{aligned} \tag{6}\]

    Using the data set from (Sharma et al. 2015) we have \(n = 58\), \(\widehat{\theta} = 60.0016\) and \(\widehat{se} \left(\widehat{\theta}\right) = 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 \(\theta\) \[\begin{aligned} f(x \mid \theta) = \dfrac {\theta}{x^2}\,\exp\left(-\dfrac {\theta}{x}\right), \quad x > 0. \end{aligned}\]

    \(\bullet\) Bias expression (Johnson et al. 2012b): \[\begin{aligned} \label{bc-invexp} \mathcal{B}\left(\widehat{\theta}\right) = \dfrac {\theta}{n}. \end{aligned} \tag{7}\]

    Using the data set from (Lawless 2011), we have \(n = 30\), \(\widehat{\theta} = 11.1786\) and \(\widehat{se}\left(\widehat{\theta}\right) = 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 \(\theta\) \[\begin{aligned} f(x \mid \theta) = \dfrac {\theta^2}{\theta^2 + 1}\,(\theta + x)\,\exp(-\theta\,x), \quad x > 0. \end{aligned}\]

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

    Using the data set from (Shanker 2015), we have \(n = 31\), \(\widehat{\theta} = 0.0647\) and \(\widehat{se} \left(\widehat{\theta}\right) = 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 \(\theta\) \[\begin{aligned} f(x \mid \theta) = \frac {\theta^2}{1 + \theta^2} \left(\frac {1 + \theta\,x}{x^3}\right)\exp\left(-\frac {\theta}{x}\right), \quad x > 0. \end{aligned}\]

    \(\bullet\) Bias expression (not previously reported in the literature): \[\begin{aligned} \label{bc-invshanker} \mathcal{B} \left(\widehat{\theta}\right) = \dfrac {\theta^3 + 2\,\theta}{n\,\left(\theta^2 + 1\right)}. \end{aligned} \tag{8}\]

    Using the data set from (Sharma et al. 2015), we have \(n = 58\), \(\widehat{\theta} = 59.1412\) and \(\widehat{se}\left(\widehat{\theta}\right) = 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 \(\nu\) \[\begin{aligned} f(x \mid \nu) = 2\,\nu\,(1 - x)\,x^{\nu-1}\,(2 - x)^{\nu - 1}, \quad 0 < x < 1. \end{aligned}\]

    \(\bullet\) Bias expression (Giles 2012a): \[\begin{aligned} \label{bc-toppleone} \mathcal{B}\left(\widehat{\nu}\right) = \dfrac {\nu}{n}. \end{aligned} \tag{9}\]

    Using the data set from (Cordeiro and dos Santos Brito 2012), we have \(n = 107\), \(\widehat{\nu} = 2.0802\) and \(\widehat{se}\left(\widehat{\nu}\right) = 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 \(\sigma\) \[\begin{aligned} f(x \mid \sigma) = \sqrt{\dfrac {\sigma}{2\,\pi}}\,x^{-\frac {3}{2}}\,\exp\left(-\dfrac {\sigma}{2\,x}\right), \quad x > 0. \end{aligned}\]

    \(\bullet\) Bias expression (not previously reported in the literature): \[\begin{aligned} \label{bc-levy} \mathcal{B} \left(\widehat{\sigma}\right) = \dfrac {2\,\sigma}{n}. \end{aligned} \tag{10}\]

    Using the data set from (Achcar et al. 2013), we have \(n = 361\), \(\widehat{\sigma} = 4.4461\) and \(\widehat{se} \left(\widehat{\sigma}\right) = 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 \(\sigma\) \[\begin{aligned} f(x \mid \sigma) = \frac {x}{\sigma^2}\,\exp\left(-\dfrac {x^2}{2\,\sigma^2}\right), \quad x > 0. \end{aligned}\]

    \(\bullet\) Bias expression (Xiao and Giles 2014): \[\begin{aligned} \label{bc-rayleigh} \mathcal{B} \left(\widehat{\sigma}\right) = -\dfrac {\sigma}{8\,n}. \end{aligned} \tag{11}\]

    Using the data set from (Bader and Priest 1982), we have \(n = 69\), \(\widehat{\sigma} = 1.2523\) and \(\widehat{se} \left(\widehat{\sigma}\right) = 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 \(\sigma\) \[\begin{aligned} f(x \mid \sigma) = \frac {2\,\sigma^2}{x^3}\,\exp\left(-\dfrac {\sigma}{x^2}\right), \quad x > 0. \end{aligned}\]

    \(\bullet\) Bias expression (not previously reported in the literature): \[\begin{aligned} \label{bc-irayleigh} \mathcal{B} \left(\widehat{\sigma}\right) = \dfrac {3\sigma}{8\,n}. \end{aligned} \tag{12}\]

    Using the data set from (Bader and Priest 1982), we have \(n = 63\), \(\widehat{\sigma} = 2.8876\) and \(\widehat{se} \left(\widehat{\sigma}\right) = 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 \(\sigma\) \[\begin{aligned} f(x \mid \sigma) = \dfrac {2\,\exp\left(-\dfrac {x}{\sigma}\right)}{\sigma\,\left[ 1 + \exp\left(-\dfrac {x}{\sigma}\right) \right]^2}, \quad x > 0. \end{aligned}\]

    \(\bullet\) Bias expressions (Giles 2012b): \[\begin{aligned} \label{bc-half-logistic} \mathcal{B}\left(\widehat{\sigma}\right) = -\dfrac {0.05256766607\,\sigma}{n}. \end{aligned} \tag{13}\]

    Using the data set from (Bhaumik et al. 2009), we have \(n = 34\), \(\widehat{\sigma} = 1.3926\) and \(\widehat{se}\left(\widehat{\sigma}\right) = 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 \(\sigma\) \[\begin{aligned} f(x \mid \sigma) = \dfrac {2}{\pi}\,\dfrac {\sigma}{\sigma^2 + x^2}, \quad x > 0. \end{aligned}\]

    \(\bullet\) Bias expression (not previously reported in the literature): \[\begin{aligned} \label{bc-half-cauchy} \mathcal{B}\left(\widehat{\sigma}\right) = -\dfrac {\sigma}{n}. \end{aligned} \tag{14}\]

    Using the data set from (Alzaatreh et al. 2016), we have \(n = 64\), \(\widehat{\sigma} = 28.3345\) and \(\widehat{se} \left(\widehat{\sigma}\right) = 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 \(\sigma\) \[\begin{aligned} f(x \mid \sigma) = \sqrt{\dfrac {2}{\pi}}\,\dfrac {1}{\sigma}\,\exp\left(-\dfrac {x^2}{2\,\sigma^2}\right), \quad x > 0. \end{aligned}\]

    \(\bullet\) Bias expressions (Xiao and Giles 2014): \[\begin{aligned} \label{bc-half-normal} \mathcal{B} \left(\widehat{\sigma}\right) = -\dfrac {\sigma}{4\,n}. \end{aligned} \tag{15}\]

    Using the data set from (Raqab et al. 2008), we have \(n = 69\), \(\widehat{\sigma} = 1.5323\) and \(\widehat{se} \left(\widehat{\sigma}\right) = 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 \(\mu\) and standard deviation \(\sigma\) \[\begin{aligned} f(x \mid \mu, \sigma) = \dfrac {1}{\sqrt{2\,\pi}\,\sigma}\,\exp\left[-\dfrac {(x - \mu)^2}{2\,\sigma^2}\right], \quad x \in (-\infty, \infty). \end{aligned}\]

    \(\bullet\) Bias expressions (Stočsić and Cordeiro 2009): \[\begin{aligned} \label{bc-normal} \mathcal{B} \left(\widehat{\mu}\right) = 0 \mbox{ and } \mathcal{B}\left(\widehat{\sigma}\right) = -\dfrac {3\,\sigma}{4\,n}. \end{aligned} \tag{16}\]

    Using the data set from (Kundu 2005), we have \(n = 23\), \(\widehat{\mu} = 4.1506\), \(\widehat{\sigma} = 0.5215\), \(\widehat{se} \left(\widehat{\mu}\right) = 0.1087\) and \(\widehat{se}\left(\widehat{\sigma}\right) = 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 \(\mu\) and shape \(\lambda\) \[\begin{aligned} f(x \mid \mu, \lambda) = \sqrt{\dfrac {\lambda}{2\,\pi\,x^3}}\,\exp\left[-\dfrac {\lambda\,(x - \mu)^2}{2\,x\,\mu^2}\right], \quad x > 0. \end{aligned}\]

    \(\bullet\) Bias expressions (Stočsić and Cordeiro 2009): \[\begin{aligned} \label{bc-invgauss} \mathcal{B}\left(\widehat{\mu}\right) = 0 \text{ and } \mathcal{B}\left(\widehat{\lambda}\right) = \frac {3\lambda}{n}. \end{aligned} \tag{17}\]

    Using the data set from (Chhikara and Folks 1977), we have \(n = 46\), \(\widehat{\mu} = 3.6067\), \(\widehat{\lambda} = 1.6584\), \(\widehat{se} \left(\widehat{\mu}\right) = 0.7843\) and \(\widehat{se}\left(\widehat{\lambda}\right) = 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 \(\mu\) and scale \(\sigma\) \[\begin{aligned} f(x \mid \mu, \sigma) = \dfrac {1}{\sqrt{2\,\pi}\,x\,\sigma}\,\exp\left[-\dfrac {(\log x - \mu)^2}{\sigma^2}\right], \quad x > 0. \end{aligned}\]

    \(\bullet\) Bias expressions (Stočsić and Cordeiro 2009): \[\begin{aligned} \label{bc-lognormal} \mathcal{B} \left(\widehat{\mu}\right) = 0 \text{ and } \mathcal{B} \left(\widehat{\sigma}\right) = -\dfrac {3\,\sigma}{4\,n}. \end{aligned} \tag{18}\]

    Using the data set from (Kumagai et al. 1989), we have \(n = 30\), \(\widehat{\mu} = 2.164\), \(\widehat{\sigma} = 1.1765\), \(\widehat{se} \left(\widehat{\mu}\right) = 0.2148\) and \(\widehat{se} \left(\widehat{\sigma}\right) = 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 \(\beta\) and scale \(\alpha\) \[\begin{aligned} f(x\mid \alpha, \beta) = \dfrac {(\beta / \alpha) \, (x / \alpha)^{\beta - 1}}{\left[1 + (x / \alpha)^\beta\right]^2}, \quad x > 0. \end{aligned}\]

    \(\bullet\) For bias expressions, see (Reath 2016).

    From (Reath 2016) we have \(n = 19\), \(\widehat{\alpha} = 6.2542\), \(\widehat{\beta} = 1.1732\), \(\widehat{se} \left(\widehat{\alpha}\right) = 2.1352\), \(\widehat{se} \left(\widehat{\beta}\right) =\) 0.2239, \(\mathcal{\widehat{B}} \left(\widehat{\alpha} \right)=0.3585\) and \(\mathcal{\widehat{B}} \left(\widehat{\beta} \right) = 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 \(\alpha\) and rate \(\lambda\) \[\begin{aligned} f(x \mid \alpha, \lambda) = \dfrac {\lambda^\alpha}{\Gamma(\alpha)}\,x^{\alpha - 1}\,\exp(-\lambda\,x), \quad x > 0. \end{aligned}\]

    \(\bullet\) Bias expressions (Giles and Feng 2009): \[\begin{aligned} \label{bc-gamma-alpha} \mathcal{B} \left(\widehat{\alpha}\right) = \dfrac {\alpha\,\left[\Psi^{\prime}(\alpha) - \alpha\Psi^{\prime\prime}(\alpha)\right] - 2}{2\,n\left[\alpha\Psi^{\prime}(\alpha) - 1\right]^2} \end{aligned} \tag{19}\] and \[\begin{aligned} \label{bc-gamma-lambda} \mathcal{B} \left(\widehat{\lambda}\right) = \dfrac {\lambda\,\left[2\,\alpha\,\left(\Psi^{\prime}(\alpha)\right)^2 - 3\,\Psi^{\prime}(\alpha) - \alpha\,\Psi^{\prime\prime}(\alpha)\right]}{2\,n\left[\alpha\Psi^{\prime}(\alpha) - 1\right]^2}. \end{aligned} \tag{20}\]

    Using the data set from (Delignette-Muller et al. 2008), we have \(n = 254\), \(\widehat{\alpha} = 4.0083\), \(\widehat{\lambda} = 0.0544\), \(\widehat{se} \left(\widehat{\alpha}\right) = 0.3413\) and \(\widehat{se} \left(\widehat{\lambda}\right) = 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 \(\alpha\) and scale \(\beta\) \[\begin{aligned} f(x\mid \alpha, \beta) = \dfrac {1}{\Gamma(\alpha)\,\beta^\alpha}x^{\alpha - 1}\,\exp\left(-\dfrac {x}{\beta}\right), \quad x > 0. \end{aligned}\]

    \(\bullet\) Bias expressions (Stočsić and Cordeiro 2009): \[\begin{aligned} \label{bc-invgamma-alpha} \mathcal{B} \left(\widehat{\alpha}\right) = {\frac {- 0.5\,{{\alpha}}^{2}\,\Psi^{\prime\prime}\left(\alpha\right) + 0.5\,\Psi^{\prime}\left( \alpha \right) {\alpha}- 1}{n\, \alpha\,\left(\Psi^{\prime} \left(\alpha\right)- 1 \right)^{2}}} \end{aligned} \tag{21}\] and \[\begin{aligned} \label{bc-invgamma-beta} \mathcal{B}\left(\widehat{\beta}\right) ={\frac { \beta\,\left( - 0.5\,\alpha \,\Psi^{\prime\prime} \left( \alpha \right) - 1.5\,\Psi^{\prime} \left(\alpha \right) + \left( \Psi^{\prime} \left( \alpha \right) \right)^{2}\alpha \right) }{n \left( \Psi^{\prime} \left( \alpha \right) \alpha - 1.0 \right)^{2}}}. \end{aligned} \tag{22}\]

    Using the data set from (Kumagai and Matsunaga 1995), we have \(n = 31\), \(\widehat{\alpha} = 1.0479\), \(\widehat{\beta} = 5.491\), \(\widehat{se} \left(\widehat{\alpha}\right) = 0.2353\) and \(\widehat{se} \left(\widehat{\beta}\right) = 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 \(\alpha\) and scale \(\beta\) \[\begin{aligned} f(x \mid \alpha, \beta) = \alpha\,\beta\,(1 + \beta\,x)^{-(\alpha + 1)}, \quad x > 0. \end{aligned}\]

    \(\bullet\) Bias expressions (Giles et al. 2013): \[\begin{aligned} \label{bc-lomax-alpha} \mathcal{B}\left(\widehat{\alpha}\right) = {\frac {2\,\alpha\, \left( \alpha+ 1 \right) \left( {\alpha}^{2} +\alpha - 2 \right) }{ \left( \alpha+ 3 \right) n}} \end{aligned} \tag{23}\] and \[\begin{aligned} \label{bc-lomax-beta} \mathcal{B}\left(\widehat{\beta}\right) = - {\frac {2\,\beta\, \left( \alpha+ 1.6485\right) \left( \alpha+ 0.3934\right) \left( \alpha- 1.5419\right) }{n\,\alpha\, \left( \alpha+ 3 \right) }}. \end{aligned} \tag{24}\]

    Using the data set from Tahir et al. (2016), we have \(n = 179\), \(\widehat{\alpha} = 4.9103\), \(\widehat{\beta} = 0.0028\), \(\widehat{se} \left(\widehat{\alpha}\right) = 0.6208\) and \(\widehat{se} \left(\widehat{\beta}\right) = {3.4803\times 10^{-4}}\). 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 \(\alpha\) and scale \(\theta\) \[\begin{aligned} f(x \mid \alpha, \theta) = \dfrac {\theta^{\alpha + 1}}{(\theta + \alpha)\,\Gamma(\alpha)}\,x^{\alpha - 1}\,(1 + x)\,\exp(-\theta x), \quad x > 0. \end{aligned}\]

    \(\bullet\) For bias expressions, see (Wang and Wang 2017):

    Using the data set from (Ghitany et al. 2013), we have \(n = 69\), \(\widehat{\alpha} = 22.8889\), \(\widehat{\theta} = 9.6246\), \(\widehat{se} \left(\widehat{\alpha}\right) = 3.9507\) and \(\widehat{se}\left(\widehat{\theta}\right) = 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 \(\alpha\) and scale \(\theta\) \[\begin{aligned} f(x \mid \beta, \mu) = \dfrac {2\,\theta^{\alpha + 1}}{\Gamma(\alpha + 1)}\,x^{2\,\alpha + 1}\,\exp\left(-\theta\,x^2 \right), \quad x > 0. \end{aligned}\]

    \(\bullet\) For bias expressions, see (Xiao and Giles 2014):

    Using the data set from (Gomes et al. 2014), we have \(n = 384\), \(\widehat{\theta} = 0.5195\), \(\widehat{\alpha} = 0.0104\), \(\widehat{se} \left(\widehat{\theta}\right)= 0.2184\) and \(\widehat{se} \left(\widehat{\alpha}\right)= 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 \(\beta\) and scale \(\mu\) \[\begin{aligned} f(x \mid \beta, \mu) = \dfrac {\beta}{\mu^\beta}\,x^{\beta - 1}\,\exp\left( -\dfrac {x}{\mu}\right)^\beta, \quad x > 0. \end{aligned}\]

    \(\bullet\) Bias expressions (the expressions below differs from (Stočsić and Cordeiro 2009)): \[\begin{aligned} \label{bc-weibull-mu} \mathcal{B} \left(\widehat{\mu}\right) = \dfrac {\mu\,(0.5543324495-0.3698145397\,\beta )}{n\,\beta^2} \end{aligned} \tag{25}\] and \[\begin{aligned} \label{bc-weibull-beta} \mathcal{B} \left(\widehat{\beta}\right) =\dfrac {1.379530692\,\beta}{n}. \end{aligned} \tag{26}\]

    From (Datta and Datta 2013), we have \(n = 50\), \(\widehat{\mu} = 2.5752\), \(\widehat{\beta} = 38.0866\), \(\widehat{se} \left(\widehat{\mu}\right)= 0.2299\) and \(\widehat{se} \left(\widehat{\beta}\right)= 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 \(\beta\) and scale \(\mu\) \[\begin{aligned} f(x \mid \beta, \alpha) = \beta\,\mu^{\beta}\,x^{-(\beta + 1)}\,\exp\left[- \left(\dfrac {\mu}{x}\right)^\beta \right], \quad x > 0. \end{aligned}\]

    \(\bullet\) Bias expressions (not previously reported in the literature): \[\begin{aligned} \label{bc-invweibull-beta} \mathcal{B} \left(\widehat{\beta}\right) = \frac {1.379530690\,\beta}{n} \end{aligned} \tag{27}\] and \[\begin{aligned} \label{bc-invweibull-mu} \mathcal{B} \left(\widehat{\mu}\right) = {\frac {\mu\, \left( 0.3698145391\,\beta+ 0.5543324494 \right) }{n{\beta}^{2}}}. \end{aligned} \tag{28}\]

    Using the data set from (Nichols and Padgett 2006), we have \(n = 100\), \(\widehat{\beta} = 1.769\), \(\widehat{\mu} = 1.8917\), \(\widehat{se} \left(\widehat{\beta}\right)= 0.1119\) and \(\widehat{se} \left(\widehat{\mu}\right)= 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 \(\alpha\) and scale \(\theta\) \[\begin{aligned} f(x \mid \alpha, \theta) = \sqrt{\dfrac {2}{\pi}}\,\dfrac {\alpha}{\theta^\alpha}\,x^{\alpha - 1}\,\exp\left[-\dfrac {1}{2} \left(\dfrac {x}{\theta} \right)^{2\,\alpha}\right]. \end{aligned}\]

    \(\bullet\) Bias expressions (Mazucheli and Dey 2017): \[\begin{aligned} \label{bc-genhalfnormal-alpha} \mathcal{B} \left(\widehat{\alpha}\right) = 1.483794456 \, \dfrac {\alpha}{n} \end{aligned} \tag{29}\] and \[\begin{aligned} \label{bc-genhalfnormal-theta} \mathcal{B} \left(\widehat{\theta}\right) = (0.2953497661 - 0.3665611957\,\alpha) \, \dfrac {\theta}{n\,\alpha^2}. \end{aligned} \tag{30}\]

    Using the data set from (Nadarajah 2008b), we have \(n = 119\), \(\widehat{\alpha} = 3.8096\), \(\widehat{\theta} = 4.9053\), \(\widehat{se} \left(\widehat{\alpha}\right) = 0.2758\) and \(\widehat{se} \left(\widehat{\theta}\right) = 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 \(\alpha\) and scale \(\theta\) \[\begin{aligned} f(x \mid \alpha, \theta) = \sqrt{\dfrac {2}{\pi}}\,\left(\dfrac {\alpha}{x}\right)\,\left(\dfrac {1}{\theta\,x}\right)^\alpha\,\exp\left[-\dfrac {1}{2} \left(\dfrac {1}{\theta\,x} \right)^{2\,\alpha}\right], \quad x > 0. \end{aligned}\]

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

    Using the data set from (Nadarajah et al. 2011), we have \(n = 20\), \(\widehat{\alpha} = 3.0869\), \(\widehat{\theta} = 0.6731\), \(\widehat{se} \left(\widehat{\alpha}\right) = 0.5534\) and \(\widehat{se} \left(\widehat{\theta}\right) = 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 \(\alpha\) and rate \(\lambda\) \[\begin{aligned} f(x \mid \alpha, \lambda) = \dfrac {\lambda\,\alpha\,\exp\left(-\lambda\,x\right)}{\left[1 - (1 - \alpha)\,\exp\left(-\lambda\,x\right)\right]^2}, \quad x > 0. \end{aligned}\]

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

    Using the data set from (Linhart and Zucchini 1986), we have \(n = 20\), \(\widehat{\alpha} = 0.2782\), \(\widehat{\lambda} = 0.0078\), \(\widehat{se} \left(\widehat{\alpha}\right) = 0.2321\) and \(\widehat{se} \left(\widehat{\lambda}\right) = 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 \(\alpha\) and \(\beta\) \[\begin{aligned} f(x \mid \alpha, \beta) = \dfrac {\Gamma(\alpha + \beta)}{\Gamma(\alpha)\,\Gamma(\beta)}\,x^{\alpha - 1}\,(1 - x)^{\beta - 1}, \quad 0 < x < 1. \end{aligned}\]

    \(\bullet\) For bias expressions, see (Cordeiro et al. 1997).

    Using the data set from (Javanshiri et al. 2015), we have \(n = 48\), \(\widehat{\alpha} = 5.941\), \(\widehat{\beta} = 21.2024\), \(\widehat{se} \left(\widehat{\alpha}\right) = 1.1812\) and \(\widehat{se} \left(\widehat{\beta}\right) = 4.3462\). Evaluating the analytical expressions in (Cordeiro et al. 1997), 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 \(\alpha\) and \(\beta\) \[\begin{aligned} f(x \mid \alpha, \beta) = \alpha\,\beta\,x^{\alpha - 1}\,(1 - x^\alpha)^{\beta - 1}, \quad 0 < x < 1. \end{aligned}\]

    \(\bullet\) For bias expressions, see (Lemonte 2011).

    Using the data set from (Wang et al. 2017), we have \(n = 20\), \(\widehat{\alpha} = 6.3478\), \(\widehat{\beta} = 4.4898\), \(\widehat{se} \left(\widehat{\alpha}\right) = 1.5576\) and \(\widehat{se} \left(\widehat{\beta}\right) = 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 \(\alpha\) and \(\beta\) \[\begin{aligned} f(x \mid \alpha, \beta) = \dfrac {\Gamma(\alpha + \beta)}{\Gamma(\alpha)\,\Gamma(\beta)}\,x^{\alpha - 1}\,(1 + x)^{-(\alpha + \beta)}, \quad x > 0. \end{aligned}\]

    \(\bullet\) For bias expressions, see (Stočsić and Cordeiro 2009).

    Using the data set from (Nadarajah 2008a), we have \(n = 116\), \(\widehat{\alpha} = 28.5719\), \(\widehat{\beta} = 1.3783\), \(\widehat{se} \left(\widehat{\alpha}\right) = 4.0367\) and \(\widehat{se} \left(\widehat{\beta}\right) = 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 \(\alpha\) and scale \(\beta\) \[\begin{aligned} f(x \mid \alpha, \beta) = \dfrac {1}{2\,\alpha\,\beta\,\sqrt{2\,\pi}}\,\left[\left(\dfrac {\beta}{x}\right)^{1/2} + \left(\dfrac {\beta}{x}\right)^{3/2}\right]\, \exp\left[-\dfrac {1}{2,\alpha^2}\,\left(\dfrac {x}{\beta} + \dfrac {\beta}{x} - 2\right)\right], \quad x > 0. \end{aligned}\]

    \(\bullet\) Bias expressions (Lemonte et al. 2007): \[\begin{aligned} \label{bc-bs-alpha} \mathcal{B} \left(\widehat{\alpha}\right) = -\dfrac {\alpha}{4\,n}\,\left(1 + \dfrac {2 + \alpha^2}{\alpha\,(2\,\pi)^{-1/2}\,h(\alpha) + 1}\right) \end{aligned} \tag{31}\] and \[\begin{aligned} \label{bc-bs-beta} \mathcal{B} \left(\widehat{\beta}\right) = \dfrac {\beta^2\,\alpha^2}{2\,n\,\left[\alpha\,(2\,\pi)^{-1/2}\,h(\alpha) + 1\right]}, \end{aligned} \tag{32}\] where \[\begin{aligned} h(\alpha) = \alpha\,\sqrt{\dfrac {\pi}{2}} - \pi\,e^{2/\alpha^2}\,\left[1 - \Phi\left(\dfrac {2}{\alpha}\right)\right]. \end{aligned}\]

    Using the data set from (Gross and Clark 1976), we have \(n = 20\), \(\widehat{\alpha} = 0.3149\), \(\widehat{\beta} = 1.8105\), \(\widehat{se} \left(\widehat{\alpha}\right) = 0.0498\) and \(\widehat{se} \left(\widehat{\beta}\right) = 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 \(\xi\) and scale \(\sigma\) \[\begin{aligned} f(x \mid \xi, \sigma) = \dfrac {1}{\sigma}\,\left(1 + \dfrac {\xi\,x}{\sigma} \right)^{-(1/\xi + 1)}, \quad x > 0, \ \xi \neq 0. \end{aligned}\]

    \(\bullet\) Bias expressions (Giles et al. 2016): \[\begin{aligned} \label{bc-gp-xi} \mathcal{B} \left(\widehat{\xi}\right) = - \dfrac {(1 + \xi)\,(3 + \xi)}{n\,(1 + 3\,\xi)} \end{aligned} \tag{33}\] and \[\begin{aligned} \label{bc-gp-sigma} \mathcal{B} \left(\widehat{\sigma} \right) = - \dfrac {\sigma\, \left(3 + 5\,\xi + 4\,\xi^2\right)}{n\,(1 + 3\,\xi)}. \end{aligned} \tag{34}\]

    Using the data set from (Ross and Lott 2003), we have \(n = 58\), \(\widehat{\xi} = 0.736\), \(\widehat{\sigma} = 1.709\), \(\widehat{se} \left(\widehat{\xi}\right) = 0.223\) and \(\widehat{se} \left(\widehat{\sigma}\right) = 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 (Cox and Lewis 1966; Jørgensen 1982), we fitted the exponentiated Weibull, Marshall-Olkin extended Weibull, Weibull, Marshall-Olkin extended exponential and exponential distributions. These distributions were also fitted by (Cordeiro and Lemonte 2013). 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 \(\mathbf{H}_1^{-1} \left(\mathbf{\widehat{\theta}}\right)\), the observed variance-covariance obtained from the analytical Hessian \(\mathbf{H}_2^{-1} \left(\mathbf{\widehat{\theta}} \right)\), the expected variance-covariance \(\mathbf{I}^{-1} \left(\mathbf{\widehat{\theta}}\right)\) and the expected variance-covariance evaluated at the bias corrected MLEs \(\mathbf{I}^{-1} \left(\mathbf{\widetilde{\theta}} \right)\). The MLEs and the \(\mathbf{H}_1^{-1} \left(\mathbf{\widehat{\theta}} \right)\) matrix were obtained by the fitdistrplus package (Delignette-Muller et al. 2017). 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 (Mudholkar and Srivastava 1993). Its probability density function is: \[\begin{aligned} f(x \mid \lambda, \beta, \alpha) = \alpha\,\beta\,\lambda\,x^{\beta - 1}\,\textrm{e}^{-\lambda\,x^\beta}\,\left(1 - \textrm{e}^{-\lambda\,x^\beta} \right)^{\alpha - 1}, \end{aligned}\] where \(\lambda > 0\) is the scale parameter and \(\beta > 0\) and \(\alpha > 0\) are the shape parameters. The Marshall-Olkin extended Weibull distribution was introduced by (Marshall and Olkin 1997). Its probability density function is: \[\begin{aligned} f(x \mid \lambda, \beta, \alpha) = \dfrac {\alpha\,\beta\,\lambda\,x^{\beta - 1}\,\textrm{e}^{-\lambda\,x^\beta}}{\left(1 - \overline{\alpha}\,\textrm{e}^{-\lambda\,x^\beta}\right)^2}, \end{aligned}\] where \(\lambda > 0\) is the scale parameter, \(\beta > 0\) is the shape parameter, \(\alpha > 0\) is an additional shape parameter and \(\overline{\alpha} = 1 - \alpha\).

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

Table 1: MLEs and bias corrected MLEs.
Distribution \(\widehat{\alpha}\) \(\widehat{\beta}\) \(\widehat{\lambda}\) \(\widetilde{\alpha}\) \(\widetilde{\beta}\) \(\widetilde{\lambda}\)
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 \(\mathbf{H}^{-1}\) and expected \(\mathbf{I}^{-1}\) information matrices. As demonstrated by (Cao 2013), the \(\mathbf{I}^{-1}\) outperforms the \(\mathbf{H}^{-1}\) 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.

\(\bullet\) Exponentiated Weibull distribution: \[\begin{aligned} \mathbf{H}_1^{-1} \left(\mathbf{\widehat{\theta}}\right) &= \left[\begin{array}{rrr} 0.00726 & -0.00717 & 0.03564 \\ -0.00717 & 0.00718 & -0.03493 \\ 0.03564 & -0.03493 & 0.18045 \\ \end{array}\right], & \mathbf{H}_2^{-1} \left(\mathbf{\widehat{\theta}}\right) &= \left[\begin{array}{rrr} 0.00729 & -0.00720 & 0.03579 \\ -0.00720 & 0.00721 & -0.03509 \\ 0.03579 & -0.03509 & 0.18120 \\ \end{array}\right], \\ \\ \mathbf{I}^{-1}\left(\mathbf{\widehat{\theta}}\right) &= \left[\begin{array}{rrr} 0.00532 & -0.00524 & 0.02609 \\ -0.00524 & 0.00527 & -0.02545 \\ 0.02609 & -0.02545 & 0.13333 \\ \end{array}\right], & \mathbf{I}^{-1}\left(\mathbf{\widetilde{\theta}}\right) &= \left[\begin{array}{rrr} 0.00510 & -0.00510 & 0.02482 \\ -0.00510 & 0.00519 & -0.02454 \\ 0.02482 & -0.02454 & 0.12590 \\ \end{array}\right]. \end{aligned}\]

\(\bullet\) Marshall-Olkin extended Weibull distribution: \[\begin{aligned} \mathbf{H}_1^{-1}\left(\mathbf{\widehat{\theta}}\right) &= \left[\begin{array}{rrr} 0.00004 & -0.00036 & 0.00052 \\ -0.00036 & 0.00361 & -0.00430 \\ 0.00052 & -0.00430 & 0.00748 \\ \end{array}\right], & \mathbf{H}_2^{-1}\left(\mathbf{\widehat{\theta}}\right) &= \left[\begin{array}{rrr} 0.00005 & -0.00047 & 0.00068 \\ -0.00047 & 0.00468 & -0.00582 \\ 0.00068 & -0.00582 & 0.00967 \\ \end{array}\right], \\ \\ \mathbf{I}^{-1}\left(\mathbf{\widehat{\theta}}\right) &= \left[\begin{array}{rrr} 0.00006 & -0.00056 & 0.00082 \\ -0.00056 & 0.00542 & -0.00699 \\ 0.00082 & -0.00699 & 0.01146 \\ \end{array}\right], & \mathbf{I}^{-1}\left(\mathbf{\widetilde{\theta}}\right) &= \left[\begin{array}{rrr} 0.00005 & -0.00051 & 0.00072 \\ -0.00051 & 0.00526 & -0.00651 \\ 0.00072 & -0.00651 & 0.01030 \\ \end{array}\right]. \end{aligned}\]

\(\bullet\) Weibull distribution: \[\begin{aligned} \mathbf{H}_1^{-1}\left(\mathbf{\widehat{\theta}}\right) &= \left[\begin{array}{rr} 0.00004 & -0.00018 \\ -0.00018 & 0.00086 \\ \end{array}\right], & \mathbf{H}_2^{-1}\left(\mathbf{\widehat{\theta}}\right) &= \left[\begin{array}{rr} 0.00004 & -0.00018 \\ -0.00018 & 0.00087 \\ \end{array}\right], \\ \\ \mathbf{I}^{-1}\left(\mathbf{\widehat{\theta}}\right) &= \left[\begin{array}{rr} 0.00004 & -0.00018 \\ -0.00018 & 0.00089 \\ \end{array}\right], & \mathbf{I}^{-1}\left(\mathbf{\widetilde{\theta}}\right) &= \left[\begin{array}{rr} 0.00004 & -0.00018 \\ -0.00018 & 0.00089 \\ \end{array}\right]. \end{aligned}\]

\(\bullet\) Marshall-Olkin extended exponential distribution: \[\begin{aligned} \mathbf{H}_1^{-1}\left(\mathbf{\widehat{\theta}}\right) &= \left[\begin{array}{rr} 0.00004 & 0.00081 \\ 0.00081 & 0.02022 \\ \end{array}\right], & \mathbf{H}_2^{-1}\left(\mathbf{\widehat{\theta}}\right) &= \left[\begin{array}{rr} 0.00004 & 0.00081 \\ 0.00081 & 0.02023 \\ \end{array}\right], \\ \\ \mathbf{I}^{-1}\left(\mathbf{\widehat{\theta}}\right) &= \left[\begin{array}{rr} 0.00004 & 0.00083 \\ 0.00083 & 0.02094 \\ \end{array}\right], & \mathbf{I}^{-1}\left(\mathbf{\widetilde{\theta}}\right) &= \left[\begin{array}{rr} 0.00004 & 0.00082 \\ 0.00082 & 0.02047 \\ \end{array}\right]. \end{aligned}\]

\(\bullet\) Exponential distribution: \[\begin{aligned} \mathbf{H}_1^{-1}\left(\mathbf{\widehat{\theta}}\right) &= 0.000010433, & \mathbf{H}_2^{-1}\left(\mathbf{\widehat{\theta}}\right) &= 0.000010436, \\ \\ \mathbf{I}^{-1}\left(\mathbf{\widehat{\theta}}\right) &= 0.000010436, & \mathbf{I}^{-1}\left(\mathbf{\widetilde{\theta}}\right) &= 0.000010410. \end{aligned}\]

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 (Stočsić and Cordeiro 2009) and (Johnson et al. 2012a). They use the commercial softwares Maple (Maple 2017) and Mathematica (Wolfram Research, Inc. 2010).

In order to calculate the bias corrected estimates in a simple way, (Mazucheli 2017) developed an R (R Core Team 2016) 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. (Efron and Johnstone 1990) showed that the Fisher information can be computed using the hazard rate function. (Gupta et al. 2004) 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 \left(x\mid\mathbf{\theta}\right) = -\frac {d}{dx} \log \left[S(x\mid\mathbf{\theta})\right]\) and \(\overline{h} \left(x\mid\mathbf{\theta} \right) = \frac {d}{dx} \log \left[F(x\mid\mathbf{\theta})\right]\), where \(S\left(x\mid\mathbf{\theta}\right)\) and \(F\left(x\mid\mathbf{\theta}\right)\) 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”.

6.1 Supplementary materials

Supplementary materials are available in addition to this article. It can be downloaded at RJ-2017-055.zip

6.2 CRAN packages used

mle.tools, fitdistrplus

6.3 CRAN Task Views implied by cited packages

ActuarialScience, Distributions, Survival

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.

References

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://doi.org/10.32614/RJ-2017-055},
  doi = {10.32614/RJ-2017-055},
  volume = {9},
  issue = {2},
  issn = {2073-4859},
  pages = {268-290}
}

© The R Foundation, web page contact.