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 (G. M. 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 (G. M. Cordeiro et al. 1997), (Cribari-Neto and Vasconcellos 2002), (Saha and Paul 2005), (Artur J. Lemonte, Cribari-Neto, and Vasconcellos 2007), (David E. Giles and Feng 2009) (Lagos-Álvarez, Jiménez-Gamero, and Alba-Fernández 2011), (A. J. Lemonte 2011), (David E. Giles 2012), (D. E. Giles 2012), (J. Schwartz, Godwin, and Giles 2013), (D. E. Giles, Feng, and Godwin 2013), (Teimouri and Nadarajah 2013), (Xiao and Giles 2014), (Zhang and Liu 2015), (Teimouri and Nadarajah 2016), (Reath 2016), (D. E. Giles, Feng, and Godwin 2016), (Jacob Schwartz and Giles 2016), (M. 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, Qi, and Chueh 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 (D. R. 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.

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)\) (David Roxbee 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 (D. R. 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} (\#eq:coxsnell)\] 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, (G. M. Cordeiro and Klein 1994) noticed that equation @ref(eq:coxsnell) 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} (\#eq:cordeiro-klein)\]

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} (\#eq:mle-coxsnell)\] Alternatively, using matrix notation the bias corrected MLEs can be expressed as (G. M. 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} (\#eq:mle-cordeiro-klein)\] where \(\widehat{K} = K\big|_{\mathbf{\theta}=\widehat{\mathbf{\theta}}}\) and \(\widehat{A} = A\big|_{\mathbf{\theta}=\widehat{\mathbf{\theta}}}\).

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 @ref(eq:coxsnell). 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.

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 (G. M. Cordeiro et al. 1997), (A. J. 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 (G. M. 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} (\#eq:bc-lindley)\]

    Using the data set from (Ghitany, Atieh, and Nadarajah 2008) we have \(n = 100\), \(\widehat{\theta} = 0.1866\) and \(\widehat{se}\left(\widehat{\theta}\right)= 0.0133\). Evaluating the analytical expression @ref(eq:bc-lindley) 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 (Wentao. 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} (\#eq:bc-invlindley)\]

    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 @ref(eq:bc-invlindley) 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, Qi, and Chueh 2012b): \[\begin{aligned} \label{bc-invexp} \mathcal{B}\left(\widehat{\theta}\right) = \dfrac {\theta}{n}. \end{aligned} (\#eq:bc-invexp)\]

    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 @ref(eq:bc-invexp) 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} (\#eq:bc-invshanker)\]

    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 @ref(eq:bc-invshanker) 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 (D. E. Giles 2012): \[\begin{aligned} \label{bc-toppleone} \mathcal{B}\left(\widehat{\nu}\right) = \dfrac {\nu}{n}. \end{aligned} (\#eq:bc-toppleone)\]

    Using the data set from (Gauss Moutinho 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 @ref(eq:bc-toppleone) 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} (\#eq:bc-levy)\]

    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 @ref(eq:bc-levy) 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} (\#eq:bc-rayleigh)\]

    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 @ref(eq:bc-rayleigh) 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} (\#eq:bc-irayleigh)\]

    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 @ref(eq:bc-irayleigh) 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 (David E. Giles 2012): \[\begin{aligned} \label{bc-half-logistic} \mathcal{B}\left(\widehat{\sigma}\right) = -\dfrac {0.05256766607\,\sigma}{n}. \end{aligned} (\#eq:bc-half-logistic)\]

    Using the data set from (Bhaumik, Kapur, and Gibbons 2009), we have \(n = 34\), \(\widehat{\sigma} = 1.3926\) and \(\widehat{se}\left(\widehat{\sigma}\right) = 0.2056\). Evaluating the analytical expression @ref(eq:bc-irayleigh) 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} (\#eq:bc-half-cauchy)\]

    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 @ref(eq:bc-half-cauchy) 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} (\#eq:bc-half-normal)\]

    Using the data set from (Raqab, Madi, and Kundu 2008), we have \(n = 69\), \(\widehat{\sigma} = 1.5323\) and \(\widehat{se} \left(\widehat{\sigma}\right) = 0.1304\). Evaluating the analytical expression @ref(eq:bc-half-normal) 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} (\#eq:bc-normal)\]

    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 @ref(eq:bc-normal) 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} (\#eq:bc-invgauss)\]

    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 @ref(eq:bc-invgauss) 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} (\#eq:bc-lognormal)\]

    Using the data set from (S. 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 @ref(eq:bc-lognormal) 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 (David E. 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} (\#eq:bc-gamma-alpha)\] 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} (\#eq:bc-gamma-lambda)\]

    Using the data set from (M. 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 @ref(eq:bc-gamma-alpha), @ref(eq:bc-gamma-lambda) 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} (\#eq:bc-invgamma-alpha)\] 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} (\#eq:bc-invgamma-beta)\]

    Using the data set from (Shinji 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 @ref(eq:bc-invgamma-alpha), @ref(eq:bc-invgamma-beta) 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 (D. E. Giles, Feng, and Godwin 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} (\#eq:bc-lomax-alpha)\] 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} (\#eq:bc-lomax-beta)\]

    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 @ref(eq:bc-lomax-alpha), @ref(eq:bc-lomax-beta) 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 (M. 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} (\#eq:bc-weibull-mu)\] and \[\begin{aligned} \label{bc-weibull-beta} \mathcal{B} \left(\widehat{\beta}\right) =\dfrac {1.379530692\,\beta}{n}. \end{aligned} (\#eq:bc-weibull-beta)\]

    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 @ref(eq:bc-weibull-mu), @ref(eq:bc-weibull-beta) 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} (\#eq:bc-invweibull-beta)\] 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} (\#eq:bc-invweibull-mu)\]

    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 @ref(eq:bc-invweibull-beta), @ref(eq:bc-invweibull-mu) 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} (\#eq:bc-genhalfnormal-alpha)\] 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} (\#eq:bc-genhalfnormal-theta)\]

    Using the data set from (S. Nadarajah 2008), 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 @ref(eq:bc-genhalfnormal-alpha), @ref(eq:bc-genhalfnormal-theta) 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 (Saralees Nadarajah, Bakouch, and Tahmasbi 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 (G. M. Cordeiro et al. 1997).

    Using the data set from (Javanshiri, Habibi Rad, and Arghami 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 (G. M. 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 (A. J. Lemonte 2011).

    Using the data set from (B. X. Wang, Wang, and Yu 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 (Saralees Nadarajah 2008), 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 (Artur J. Lemonte, Cribari-Neto, and Vasconcellos 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} (\#eq:bc-bs-alpha)\] 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} (\#eq:bc-bs-beta)\] 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 @ref(eq:bc-bs-alpha), @ref(eq:bc-bs-beta) 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 (D. E. Giles, Feng, and Godwin 2016): \[\begin{aligned} \label{bc-gp-xi} \mathcal{B} \left(\widehat{\xi}\right) = - \dfrac {(1 + \xi)\,(3 + \xi)}{n\,(1 + 3\,\xi)} \end{aligned} (\#eq:bc-gp-xi)\] 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} (\#eq:bc-gp-sigma)\]

    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 @ref(eq:bc-gp-xi), @ref(eq:bc-gp-sigma) 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

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 (D. R. 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 (Gauss M. 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 (M. L. Delignette-Muller, Dutang, and Siberchicot 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}\]

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, Qi, and Chueh 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, Gupta, and Sankaran 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:

  • the MLEs of \(g\left({\mathbf{\theta}}\right)\) and \({\textrm{Var}} \left[g \left({\mathbf{\theta}} \right) \right]\),

  • the negative log likelihood value \(-2\log (L)\),

  • the Akaike’s information criterion \(-2\log (L)+2p\),

  • the corrected Akaike’s information criterion \(-2\log (L)+\frac {2np}{n-p-1}\),

  • the Schwarz’s Bayesian information criterion \(-2\log (L)+p\log (n)\),

  • the Hannan-Quinn information criterion \(-2\log (L)+2\log \log(n) p\),

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”.

Achcar, Jorge A, Sílvia RC Lopes, Josmar Mazucheli, and Raquel R Linhares. 2013. “A Bayesian Approach for Stable Distributions: Some Computational Aspects.” Open Journal of Statistics 3 (04): 268.
Alzaatreh, Ayman, M Mansoor, MH Tahir, M Zubair, and Shakir Ali. 2016. “The Gamma Half-Cauchy Distribution: Properties and Applications.” HACETTEPE Journal of Mathematics and Statistics 45 (4): 1143–59.
Bader, MG, and AM Priest. 1982. “Statistical Aspects of Fibre and Bundle Strength in Hybrid Composites.” Progress in Science and Engineering of Composites, 1129–36.
Bhaumik, Dulal K., Kush Kapur, and Robert D. Gibbons. 2009. “Testing Parameters of a Gamma Distribution for Small Samples.” Technometrics 51 (3): 326–34.
Cao, Xumeng. 2013. “Relative Performance of Expected and Observed Fisher Information in Covariance Estimation for Maximum Likelihood Estimates.” PhD thesis, Johns Hopkins University.
Chhikara, RS, and JL Folks. 1977. “The Inverse Gaussian Distribution as a Lifetime Model.” Technometrics 19 (4): 461–68.
Cordeiro, G. M., and F. Cribari-Neto. 2014. An Introduction to Bartlett Correction and Bias Reduction. New York: Springer-Verlag.
Cordeiro, G. M, E. C. Da Rocha, J. G. C Da Rocha, and F. Cribari-Neto. 1997. “Bias-Corrected Maximum Likelihood Estimation for the Beta Distribution.” Journal of Statistical Computation and Simulation 58 (1): 21–35.
Cordeiro, G. M., and R. Klein. 1994. “Bias Correction in ARMA Models.” Statistics & Probability Letters 19 (3): 169–76.
Cordeiro, Gauss M., and Artur J. Lemonte. 2013. “On the Marshall-Olkin Extended Weibull Distribution.” Statistical Papers 54 (2): 333–53.
Cordeiro, Gauss Moutinho, and Rejane dos Santos Brito. 2012. “The Beta Power Distribution.” Brazilian Journal of Probability and Statistics 26 (1): 88–112.
Cox, D. R., and P. A. W. Lewis. 1966. The Statistical Analysis of Series of Events. Netherlands: Springer-Verlag.
Cox, D. R., and E. J. Snell. 1968. “A General Definition of Residuals.” Journal of the Royal Statistical Society B 30 (2): 248–75.
Cox, David Roxbee, and David Victor Hinkley. 1979. Theoretical Statistics. CRC Press.
Cribari-Neto, F., and K. L. P. Vasconcellos. 2002. “Nearly Unbiased Maximum Likelihood Estimation for the Beta Distribution.” Journal of Statistical Computation and Simulation 72 (2): 107–18.
Datta, Debanshee, and D. Datta. 2013. “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.
Delignette-Muller, Marie Laure, Christophe Dutang, and Aurélie Siberchicot. 2017. : Help to Fit of a Parametric Distribution to Non-Censored or Censored Data.
Delignette-Muller, ML, M Cornu, Afssa Stec Study Group, et al. 2008. “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–64.
Edwards, A. W. F. 1992. Likelihood (Expanded Edition). Baltimore: Johns Hopkins University Press.
Efron, B., and I. M. Johnstone. 1990. “Fisher’s Information in Terms of the Hazard Rate.” The Annals of Statistics 18 (1): 38–62.
Ghitany, M. E., D. K. Al-Mutairi, N. Balakrishnan, and L. J. Al-Enezi. 2013. “Power Lindley Distribution and Associated Inference.” Computational Statistics & Data Analysis 64: 20–33.
Ghitany, M. E., B. Atieh, and S. Nadarajah. 2008. “Lindley Distribution and Its Application.” Mathematics and Computers in Simulation 78 (4): 493–506.
Giles, D. E. 2012. “A Note on Improved Estimation for the Topp-Leone Distribution.” Department of Economics, University of Victoria.
Giles, D. E., H. Feng, and R. T. Godwin. 2013. “On the Bias of the Maximum Likelihood Estimator for the Two-Parameter Lomax Distribution.” Communications in Statistics - Theory and Methods 42 (11): 1934–50.
———. 2016. “Bias-Corrected Maximum Likelihood Estimation of the Parameters of the Generalized Pareto Distribution.” Communications in Statistics – Theory and Methods 45 (8): 2465–83.
Giles, David E. 2012. “Bias Reduction for the Maximum Likelihood Estimators of the Parameters in the Half-Logistic Distribution.” Communication in Statistics - Theory and Methods 41 (2): 212–22.
Giles, David E, and Hui Feng. 2009. “Bias of the Maximum Likelihood Estimators of the Two-Parameter Gamma Distribution Revisited.” Department of Economics, University of Victoria.
Gomes, A. E., C. Q. da-Silva, G. M; Cordeiro, and E. M. M. Ortega. 2014. “A New Lifetime Model: The Kumaraswamy Generalized Rayleigh Distribution.” Journal of Statistical Computation and Simulation 84 (2): 290–309.
Gross, Alan John, and Virginia A. Clark. 1976. Survival Distributions: Reliability Applications in the Biometrical Sciences. Edited by New York. John Wiley & Sons.
Gupta, R. D., R. C. Gupta, and P. G. Sankaran. 2004. “Some Characterization Results Based on Factorization of the (Reversed) Hazard Rate Function.” Communication in Statistics - Theory and Methods 33.
Javanshiri, Zohre, Arezou Habibi Rad, and Nasser Reza Arghami. 2015. “Exp-Kumaraswamy Distributions: Some Properties and Applications.” Journal of Sciences, Islamic Republic of Iran 26 (1): 57–69.
Johnson, P. H., Y. Qi, and Y. C. Chueh. 2012a. “Bias-Corrected Maximum Likelihood Estimation in Actuarial Science.” Working Papers. Urbana-Champaign, Champaign, IL: University of Illinois.
———. 2012b. “CSCK MLE Bias Calculation.” Working Papers. Urbana-Champaign, Champaign, IL: University of Illinois.
Jørgensen, B. 1982. Statistical Properties of the Generalized Inverse Gaussian Distribution. Edited by New York. Springer-Verlag.
Kay, S. 1995. “Asymptotic Maximum Likelihood Estimator Performance for Chaotic Signals in Noise.” IEEE Transactions on Signal Processing 43 (4): 1009–12.
Kumagai, Shinji, and Ichiro Matsunaga. 1995. “Changes in the Distribution of Short-Term Exposure Concentration with Different Averaging Times.” American Industrial Hygiene Association 56 (1): 24–31.
Kumagai, S, I Matsunaga, K Sugimoto, Y Kusaka, and T Shirakawa. 1989. “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–26.
Kundu, Debasis. 2005. “Discriminating Between Normal and Laplace Distributions.” In Advances in Ranking and Selection, Multiple Comparisons, and Reliability, edited by N Balakrishnan, Nandini Kannan, and H. N. Nagaraja, 65–79. Basel: Birkhäuser.
Lagos-Álvarez, B., M. D. Jiménez-Gamero, and V. Alba-Fernández. 2011. “Bias Correction in the Type I Generalized Logistic Distribution.” Communications in Statistics - Simulation and Computation 40 (4): 511–31.
Lawless, Jerald F. 2011. Statistical Models and Methods for Lifetime Data, Second Edition. Vol. 362. John Wiley & Sons.
Lehmann, E. L. 1999. Elements of Large-Sample Theory. Springer-Verlag.
Lemonte, A. J. 2011. “Improved Point Estimation for the Kumaraswamy Distribution.” Journal of Statistical Computation and Simulation 81 (12): 1971–82.
Lemonte, Artur J., Francisco Cribari-Neto, and Klaus L. P. Vasconcellos. 2007. “Improved Statistical Inference for the Two-Parameter Birnbaum-Saunders Distribution.” Computational Statistics & Data Analysis 51 (9): 4656–81.
Linhart, H., and W. Zucchini. 1986. Model Selection. Edited by New York. John Wiley & Sons.
Maple. 2017. Maplesoft, a Division of Waterloo Maple Inc. Ontario: Waterloo.
Marshall, Albert W., and Ingram Olkin. 1997. “A New Method for Adding a Parameter to a Family of Distributions with Application to the Exponential and Weibull Families.” Biometrika 84 (3): 641–52.
Mazucheli, J. 2017. Mle.tools: Expected/Observed Fisher Information and Bias-Corrected Maximum Likelihood Estimate(s).
Mazucheli, J., and S. Dey. 2017. “Bias-Corrected Maximum Likelihood Estimation of the Parameters of the Generalized Half-Normal Distribution.” Submitted to Journal of Statistical Computation and Simulation.
Mudholkar, G. S., and D. K. Srivastava. 1993. “Exponentiated Weibull Family for Analyzing Bathtub Failure-Rate Data.” IEEE Transactions on Reliability 42 (2): 299–302.
Nadarajah, S. 2008. “The Model for Fracture Toughness.” Journal of Mechanical Science and Technology 22: 1255–58.
Nadarajah, Saralees. 2008. “A Truncated Inverted Beta Distribution with Application to Air Pollution Data.” Stochastic Environmental Research and Risk Assessment 22 (2): 285–89.
Nadarajah, Saralees, Hassan S Bakouch, and Rasool Tahmasbi. 2011. “A Generalized Lindley Distribution.” Sankhya B 73 (2): 331–59.
Nichols, Michele D, and WJ Padgett. 2006. “A Bootstrap Control Chart for Weibull Percentiles.” Quality and Reliability Engineering International 22 (2): 141–51.
R Core Team. 2016. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing.
Raqab, Mohammad Z, Mohamed T Madi, and Debasis Kundu. 2008. “Estimation of P(Y \(<\) X) for the Three-Parameter Generalized Exponential Distribution.” Communications in Statistics - Theory and Methods 37 (18): 2854–64.
Reath, Joseph. 2016. “Improved Parameter Estimation of the Log-Logistic Distribution with Applications.” PhD thesis, Michigan Technological University.
Ross, T., and N. Lott. 2003. “A Climatology of 1980 – 2003 Extreme Weather and Climate Events.” 01. National Climatic Data Center.
Saha, Krishna, and Sudhir Paul. 2005. “Bias-Corrected Maximum Likelihood Estimator of the Negative Binomial Dispersion Parameter.” Biometrics 61 (1): 179–85.
Schwartz, Jacob, and David E. Giles. 2016. “Bias-Reduced Maximum Likelihood Estimation of the Zero-Inflated Poisson Distribution.” Communications in Statistics - Theory and Methods 45 (2): 465–78.
Schwartz, J., R. T. Godwin, and D. E. Giles. 2013. “Improved Maximum-Likelihood Estimation of the Shape Parameter in the Nakagami Distribution.” Journal of Statistical Computation and Simulation 83 (3): 434–45.
Shanker, Rama. 2015. “Shanker Distribution and Its Applications.” International Journal of Statistics and Applications 5 (6): 338–48.
Sharma, V. K., S. K. Singh, U. Singh, and V. Agiwal. 2015. “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–73.
Stočsić, B. D., and G. M Cordeiro. 2009. “Using Maple and Mathematica to Derive Bias Corrections for Two Parameter Distributions.” Journal of Statistical Computation and Simulation 79 (6): 751–67.
Tahir, MH, M Adnan Hussain, Gauss M Cordeiro, GG Hamedani, M Mansoor, and M Zubair. 2016. “The Gumbel-Lomax Distribution: Properties and Applications.” Journal of Statistical Theory and Applications 15 (1): 61–79.
Teimouri, M., and S. Nadarajah. 2013. “Bias Corrected MLEs for the Weibull Distribution Based on Records.” Statistical Methodology 13: 12–24.
———. 2016. “Bias Corrected MLEs Under Progressive Type-II Censoring Scheme.” Journal of Statistical Computation and Simulation 86 (14): 2714–26.
Wang, Bing Xing, Xiu Kun Wang, and Keming Yu. 2017. “Inference on the Kumaraswamy Distribution.” Communications in Statistics - Theory and Methods 46 (5): 2079–90.
Wang, M., and W. Wang. 2017. “Bias-Corrected Maximum Likelihood Estimation of the Parameters of the Weighted Lindley Distribution.” Communications in Statistics - Theory and Methods 46 (1): 530–45.
Wang, Wentao. 2015. “Bias-Corrected Maximum Likelihood Estimation of the Parameters of the Weighted Lindley Distribution.” PhD thesis, Michigan Technological University.
Wolfram Research, Inc. 2010. Mathematica, Version 8.0. Vienna, Champaign, Illinois: Wolfram Research, Inc.
Xiao, L., and D. E. Giles. 2014. “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–92.
Zhang, Guoyi, and Rong Liu. 2015. “Bias-Corrected Estimators of Scalar Skew Normal.” Communications in Statistics - Simulation and Computation 46 (2): 831–39.