Second and Third Order Bias Reduction for One–Parameter Family Models

Silvia L.P. Ferrari

Dept. de Estatística, Univ. de São Paulo, Caixa Postal 66281, Ag. Cidade de São Paulo, São Paulo SP, 05389-970, Brazil

Denise A. Botter

Dept. de Estatística, Univ. de São Paulo, Caixa Postal 66281, Ag. Cidade de São Paulo, São Paulo SP, 05389-970, Brazil

Gauss M. Cordeiro

Dept. de Estatística, Universidade Federal de Pernambuco, Cidade Universitária, Recife PE, 50740-540, Brazil

Francisco Cribari–Neto

Department of Economics, Southern Illinois University, Carbondale IL, 62901-4515, USA

Abstract. In this paper we derive second and third order bias-corrected maximum likelihood estimates in general uniparametric models. We compare the corrected estimates and the usual maximum likelihood estimate in terms of their mean squared errors. We also obtain closed-form expressions for bias-corrected estimates in one-parameter exponential family models. Our results cover many important and commonly used distributions.

Some key words: Asymptotic expansion; bias correction; exponential family; maximum likelihood estimate.

1. Introduction

Consider the probability or density function

π(y; θ) = exp{t(y; θ)},  (1)

where θ is a scalar parameter, and let θ̂ denote the maximum likelihood estimate of θ. Expansions for the bias and variance of θ̂ to order n^{−2}, where n is the sample size, can be written as

B(θ) = \frac{B_{1}(θ)}{n} + \frac{B_{2}(θ)}{n^{2}} and V (θ) = \frac{V_{1}(θ)}{n} + \frac{V_{2}(θ)}{n^{2}},

respectively. We define three bias-corrected estimates, namely:

θ̂_{1} = θ̂ − \frac{B_{1}(θ̂)}{n},  (2)

θ̂_{2} = θ̂ − \frac{B_{1}(θ̂)}{n} − \frac{B_{2}(θ̂)}{n^{2}}  (3)

and

θ̃_{2} = θ̂ − \frac{B_{1}(θ̂)}{n} − \frac{B∗_{2}(θ̂)}{n^{2}},  (4)

where

B∗_{2}(θ) = B_{2}(θ) − B_{1}(θ)B′_{1}(θ) − \frac{1}{2}B′′_{1}(θ)V_{1}(θ),

primes denoting derivatives with respect to θ. These three modified estimates are bias-free to order n^{−1}, but only θ̃_{2} has no bias to order n^{−2}. Although θ̂_{2} seems to be a third order bias-corrected estimate, its third order bias is in general nonnull, as we shall prove later. It will also be proved that the modified estimates in (2)-(4) have the same variance and mean squared error to order n^{−2}. Closed-form formulas for B_{1}(θ), B_{2}(θ), B∗_{2}(θ), V_{1}(θ) and V_{2}(θ) for general uniparametric models are given in the next section together with a comparison of the variances and mean squared errors of the maximum likelihood estimate and its bias-corrected versions. Section 3 presents closed-form expressions for the bias and variance coefficients for one-parameter exponential family models. It also gives bias-corrected estimates for natural exponential family models.

2. Derivation of the second and third order biases

Let y_{1}, …, y_{n} be a set of n independent and identically distributed random variables with regular density (or probability) function given by (1). The asymptotic expansions for the bias and variance of the maximum likelihood estimate of θ up to order n^{−2} were written in the previous section as B(θ) = n^{−1}B_{1}(θ) + n^{−2}B_{2}(θ) and V (θ) = n^{−1}V_{1}(θ) + n^{−2}V_{2}(θ), respectively. In order to obtain closed-form expressions for the coefficients in these expansions, we need to introduce some notation for cumulants of log-likelihood derivatives with respect to θ for a single observation. Let

κ_{θθ} = E{t′′(y; θ)}, κ_{θ,θ} = −κ_{θθ} = E[{t′(y; θ)}^{2}], κ_{θθθ} = E{t′′′(y; θ)}, κ_{θθθθ} = E{t^{iv}(y; θ)}, κ_{θθθθθ} = E{t^{v}(y; θ)}, κ_{θθ,θ} = E{t′′(y; θ)t′(y; θ)},

κ_{θθ,θθ} = E[{t′′(y; θ)}^{2}] − κθθ^{2}, κ_{θ,θθθ} = E{t′(y; θ)t′′′(y; θ)},

κ_{θ,θθθθ} = E{t′(y; θ)t^{iv}(y; θ)}, κ_{θ,θθ,θθ} = E[t′(y; θ){t′′(y; θ)}^{2}] − 2κ_{θ,θθ}κ_{θθ},

κ_{θθ,θθθ} = E\{t′′(y; θ)t′′′(y; θ)\} − κ_{θθ}κ_{θθθ}, κθθ^{(θ)} = \frac{dκ_{θθ}}{dθ} , κ_{θθ}^{(θθ)} = \frac{d^{2}κ_{θθ}}{dθ^{2}} , etc.

Using the results in Shenton and Bowman (1977, pp.44-47) together with the notation for cumulants introduced above and the regularity equations (Lawley, 1956) κ_{θ,θ} = −κ_{θθ}, κ_{θ,θ,θ} =

2κ_{θθθ} − 3κθθ^{(θ)}, κ_{θ,θθ} = κθθ^{(θ)} − κ_{θθθ}, κ_{θ,θ,θθ} = κ_{θθθθ} − 2κθθθ^{(θ)} + κ_{θθ}^{(θθ)} − κ_{θθ,θθ}, we get, after some algebra,

B_{1}(θ) = −\frac{1}{2κθθ^{2}} (κ_{θθθ} − 2κθθ^{(θ)}),
B_{2}(θ) = −\frac{1}{8κθθ^{3}} (κ_{θθθθθ} + 12κ_{θθ,θθθ} + 4κ_{θ,θθθθ} + 8κ_{θ,θθ,θθ} + 4κ_{θ,θ,θθθ}) + \frac{1}{12κθθ^{4}} (13κ_{θθθθ}κ_{θθθ}

− 18κ_{θθ,θθ}κ_{θθθ} − 18κ_{θθθθ}κθθ^{(θ)} + 36κ_{θθ,θθ}κθθ^{(θ)} − 36κ_{θθθ}κθθθ^{(θ)} + 36κθθ^{(θ)}κθθθ^{(θ)} + 18κ_{θθθ}κ_{θθ}^{(θθ)})

− \frac{1}{8κθθ^{5}} \{11κθθθ^{3} − 48κθθθ^{2}κθθ^{(θ)} + 48κ_{θθθ}(κθθ^{(θ)})^{2}\},

V_{1}(θ) = −\frac{1}{κ_{θθ}} ,
V_{2}(θ) = −\frac{1}{κθθ^{3}} (κ_{θθ,θθ} − κθθθ^{(θ)} + 2κ_{θθ}^{(θθ)}) + \frac{1}{2κθθ^{4}} (10κθθ^{(θ)^{2}} − κθθθ^{2} − 4κθθ^{(θ)}κ_{θθθ}).

Bias-corrected estimates can be obtained from (2)-(4) and the equations above.

In what follows, we shall use some expansions to obtain the bias, variance and mean squared error of the estimates θ̂_{1}, θ̂_{2} and θ̃_{2}. Let H(·) be a function defined in the parameter space that can be expanded in a Taylor series. Expanding H(θ̂) around θ we obtain

H(θ̂) = H(θ) + H′(θ)(θ̂ − θ) + \frac{1}{2}H′′(θ)(θ̂ − θ)^{2} + …  (5)

Taking H(θ) = θ − B_{1}(θ)/n and using (5) we get, after some algebra, that E(θ̂_{1} − θ) = B∗_{2}(θ)/n^{2} + O(n^{−3}). Similarly, taking H(θ) = θ − B_{1}(θ)/n − B_{2}(θ)/n^{2}, we obtain E(θ̂_{2} − θ) = {B∗_{2}(θ) − B_{2}(θ)}/n^{2} + O(n^{−3}). Thus, θ̂_{2} is, like θ̂_{1}, an unbiased estimate to order n^{−1}. Next, letting H(θ) = θ − B_{1}(θ)/n − B∗_{2}(θ)/n^{2}, we obtain E(θ̃_{2} − θ) = O(n^{−3}). It then follows that only θ̃_{2} is unbiased to order n^{−2}. Although one would expect θ̂_{2} to be a third order bias-corrected estimate of θ, it was shown above that the correction term B_{1}(θ̂)/n + B_{2}(θ̂)/n^{2} does not make the n^{−2} term in the bias expansion vanish. This leads us to consider the correction term B_{1}(θ̂)/n + B∗_{2}(θ̂)/n^{2}, which, as shown, delivers a third order bias-corrected estimate.

By making use of similar developments, it is possible to show that, to order n^{−2}, the mean squared errors (MSE) of the three corrected estimates are the same and coincide with their respective variances. They are given by

MSE(θ̂_{1}) = MSE(θ̂_{2}) = MSE(θ̃_{2}) = \frac{V_{1}(θ)}{n} + \frac{1}{n^{2}} \{V_{2}(θ) − 2V_{1}(θ)B′_{1}(θ)\} + O(n^{−3}).

Therefore, from the three bias-corrected estimates, the estimate θ̃_{2} is to be preferred since they all have the same variance to order n^{−2} and θ̃_{2} is the only one that is bias-free to this order.

Finally, it remains to compare θ̂ and θ̃_{2}. The mean squared error of θ̂ is given by

MSE(θ̂) = \frac{V_{1}(θ)}{n} + \frac{1}{n^{2}} \{V_{2}(θ) + B_{1}(θ)^{2}\} + O(n^{−3}),

and hence

MSE(θ̂) − MSE(θ̃_{2}) = \frac{1}{n^{2}} \{B_{1}(θ)^{2} + 2V_{1}(θ)B′_{1}(θ)\} + O(n^{−3}).

The choice between θ̂ and θ̃_{2} is to be based on the sign of Δ = Δ(θ) = B_{1}(θ)^{2} + 2V_{1}(θ)B′_{1}(θ). For example, when B_{1}(θ) = c, a constant, we have that, to order n^{−2}, MSE(θ̂) ≥ MSE(θ̃_{2}).

Oftentimes one is interested in the estimation of a function of θ, say τ = G(θ), which does not depend on n. One can reparameterize the model in terms of τ and then use the results above to obtain τ̂_{1}, τ̂_{2} and τ̃_{2}, that is, bias-corrected maximum likelihood estimates for τ.

3. One-parameter exponential family

Let y_{1}, …, y_{n} be a set of n independent and identically distributed random variables, and assume an exponential family form for π(·; ·): π(y; θ) is as defined in (1) with t(y; θ) = − log ζ(θ) − α(θ)d(y) + v(y), that is,

π(y; θ) = \frac{1}{ζ(θ)} exp\{−α(θ) d(y) + v(y)\},  (6)

where θ is a scalar parameter and ζ = ζ(·), α = α(·), d(·) and v(·) are known functions, and it is assumed that the support of π(y; θ) does not depend upon θ. It is also assumed that α and ζ have continuous first five derivatives with respect to θ, and that ζ is positive valued, dα(θ)/dθ and dβ(θ)/dθ being different from zero for all values of θ in the parameter space, where β = β(θ) is defined as β = (dζ(θ)/dθ)(ζ(θ)dα(θ)/dθ)^{−1}. Many commonly used distributions in applied research are special cases of the family of distributions in (6). Also, this family of distributions enjoys important mathematical properties; see Bickel and Doksum (1977). We then have that t′(y; θ) = −α′{β + d(y)}, and since E{t′(y; θ)} = 0 it follows that μ = E{d(y)} = −β. The maximum likelihood estimate θ̂ comes from n^{−1} ∑ d(y_{i}) = β(θ̂). This equation may not be easy to solve. It might even require the use of numerical methods.

It can also be shown that t′′(y; θ) = −α′′{β + d(y)} − α′β′, t′′′(y; θ) = −α′′′{β + d(y)} − 2α′′β′ − α′β′′, t^{iv}(y; θ) = −α^{iv}{β + d(y)} − 3(α′′′β′ + α′′β′′) − α′β′′′, and t^{v}(y; θ) = −α^{v}{β + d(y)} − 4α^{iv}β′ − α′′β^{iv} − α′′′β′′′ − 3(α′′β′′′ + α′′′β′′). Using the results above, we can show, after some algebra, that

B_{1}(θ) = −\frac{β′′}{2α′β^{′2}} ,  (7)
B_{2}(θ) = \frac{−12β′α′′β^{′′2} − 33α′β^{′′3} + 4β^{′2}α′′β′′′ + 26α′β′β′′β′′′ − 3α′β^{′2}β^{iv}}{24α^{′3}β^{′5}} ,  (8)
B∗_{2}(θ) = \frac{1}{24α^{′4}β^{′5}} (12β^{′2}α^{′′2}β′′ + 18α′β′α′′β^{′′2} + 15α^{′2}β^{′′3} − 6α′β^{′2}β′′α′′′ − 8α′β^{′2}α′′β′′′
− 16α^{′2}β′β′′β′′′ + 3α^{′2}β^{′2}β^{iv}),  (9)
V_{1}(θ) = \frac{1}{α′β′} ,  (10)
V_{2}(θ) = \frac{2β′α′′β′′ + 5α′β′′ − 2α′β′β′′′}{2α^{′3}β^{′4}} .  (11)

It is noteworthy that the expressions for B_{1}(θ), B_{2}(θ), B∗_{2}(θ), V_{1}(θ) and V_{2}(θ) given in (7)-(11) only require knowledge of α, β and their first few derivatives with respect to θ. Unlike the general expressions given in the previous section, these expressions do not require the evaluation of moments or cumulants. Formulas (7)-(9) can be implemented in a computer algebra system such as Maple (Abell and Braselton, 1994) or Mathematica (Wolfram, 1991) to obtain bias-corrected estimates with minimal effort. They are capable of generating both simple and complex expressions for different special cases. Also, it is possible to use the expressions above to examine under what conditions the maximum likelihood estimate of θ is unbiased to a certain order of magnitude. For example, equations (7) and (8) imply that the second and third order biases vanish when β′′ = 0.

It is possible to check equations (7)-(9) from first principles for some special distributions. The simplest special case is the normal distribution with mean μ and variance θ for which B_{1}(θ) and B_{2}(θ) vanish. Here, θ̂ = ∑(y_{i} − μ)^{2}/n ∼ θχn^{2}/n and μ̂ = ∑ y_{i}/n are clearly unbiased estimates. It is also easy to verify that the same happens for the following distributions: binomial, Poisson, exponential, Laplace, truncated extreme value and inverse Gaussian with known scale parameter θ. For the inverse Gaussian distribution with known mean μ > 0 and scale parameter θ > 0, for which α(θ) = θ, ζ(θ) = θ^{−1/2}, d(y) = (y − μ)^{2}/(2μ^{2}y) and v(y) = −{log(2πy^{3})}/2, θ̂ = nμ^{2}{∑(y_{i} − μ)^{2}/y_{i}}^{−1} ∼ nθ/χn^{2}. By Taylor series expansion we can prove that E(θ̂ − θ) = 2θ/n + 4θ/n^{2} + O(n^{−3}), in agreement with the corresponding expressions for B_{1}(θ) and B_{2}(θ) that are obtained using (7) and (8). For the gamma distribution with known index k > 0 and scale parameter θ > 0, for which α(θ) = θ, ζ(θ) = θ^{−k}, d(y) = y and v(y) = (k − 1) log(y) − log{Γ(k)}, Γ(·) being the gamma function, the distribution of θ̂ is 2nkθ/χ2nk^{2}. Thus, by direct expansion E(θ̂ − θ) = θ/(nk) + θ/(nk)^{2} + O(n^{−3}), as can be obtained from (7) and (8).

For the Pareto distribution with k > 0 known (y > k) and scale parameter θ > 0, where α(θ) = θ + 1, ζ(θ) = (θk^{θ})^{−1}, d(y) = log(y) and v(y) = 0, 2nθ/θ̂ ∼ χ2n^{2}, we have that E(θ̂ − θ) = θ/n + θ/n^{2} + O(n^{−3}), leading to the corresponding expressions for B_{1}(θ) and B_{2}(θ). The maximum likelihood estimate θ̂ of θ ∈ IR for the extreme value distribution, ϕ > 0 known, for which α(θ) = exp{θ/ϕ}, ζ(θ) = ϕ exp{−θ/ϕ}, d(y) = exp{−y/ϕ} and v(y) = −y/ϕ, reduces to θ̂ = −ϕ log(n^{−1} ∑ e^{−y_{i}/ϕ}) and we can also show that E(θ̂ − θ) = ϕ/(2n) + ϕ/(12n^{2}) + O(n^{−3}). This result again agrees with the expressions of B_{1}(θ) and B_{2}(θ) obtained from (7) and (8). Finally, it is of interest to note that for the lognormal distribution with known location parameter μ > 0 and scale parameter θ > 0, where α(θ) = θ^{−2}, ζ(θ) = θ, d(y) = {log(y) − μ}^{2}/2 and v(y) = − log(y) + {log(2π)}/2, θ̂^{2} ∼ θ^{2}χn^{2}/n. The resulting expansion E(θ̂ − θ) = −θ/(4n) + θ/(32n^{2}) to order n^{−2} follows after some algebra, which agrees with the results from (7) and (8).

Next, we obtain the bias-corrected maximum likelihood estimates of the canonical parameter α = α(θ) in uniparametric natural exponential family models. The one-parameter exponential family (6) parameterized in terms of the canonical parameter α = α(θ) has the form

π(y; α) = \frac{1}{δ(α)} exp\{−αd(y) + v(y)\},  (12)

where −d(y) is the canonical statistic and δ(α) is its cumulant generating function. Similar conditions to the ones stated for the probability or density function in (6) are assumed to hold.

The rth cumulant of −d(y) is κ_{r} = d^{r}δ(α)/dα^{r}, for r = 1, 2, … Letting β(α) = δ′(α)/δ(α) we write the (r + 1)th cumulant of −d(y) as κ_{r+1} = d^{r}β/dα^{r}, for r = 1, 2, … From now on primes will denote derivatives with respect to α. The mean and variance of −d(y) are κ_{1} = β and κ_{2} = dβ/dα = β′. The variance function β′ uniquely characterizes the distribution in (12) and does not depend upon the particular parameterization used.

Using the notation introduced in Section 1, let α̂_{1}, α̂_{2} and α̃_{2} be the bias corrected maximum likelihood estimates defined by equations (2), (3) and (4), respectively. The n^{−1} and n^{−2} biases of the maximum likelihood estimate α̂ of α are obtained from equations (7) and (8) as B_{1}(α) = −β′′/(2β^{′2}) and B_{2}(α) = (−33β^{′′3}+26β′β′′β′′′−3β^{′2}β^{iv})/(24β^{′5}). Also, formulas (10) and (11) yield, to order n^{−2}, var(α̂) = V_{1}(α)/n + V_{2}(α)/n^{2} with V_{1}(α) = 1/β′ and V_{2}(α) = (5β^{′′2} − 2β′β′′′)/(2β^{′4}). Let γ_{1}^{2} = β^{′′2}/β^{′3} and γ_{2} = β′′′/β^{′2} be the third and fourth standardized cumulants of −d(y). From the equations for the mean squared errors given in Section 2 we find, after some algebra, that to order n^{−2}

MSE(α̂_{1}) = MSE(α̂_{2}) = MSE(α̃_{2}) = \frac{1}{nβ′} + \frac{γ_{1}^{2}}{2n^{2}β′}

and

MSE(α̂) − MSE(α̃_{2}) = \frac{9γ_{1}^{2} − 4γ_{2}}{4n^{2}β′} .

Since β′ > 0, it is clear that MSE(α̂) < MSE(α̃_{2}) when 9γ_{1}^{2} − 4γ_{2} < 0 and MSE(α̂) > MSE(α̃_{2}) when 9γ_{1}^{2} − 4γ_{2} > 0. The values of γ_{1} and γ_{2} must satisfy the inequalities γ_{1}^{2} + γ_{2} + 2 ≥ 0 and γ_{2} − γ_{1}^{2} ≥ 1. Then, the sign of 9γ_{1}^{2} − 4γ_{2} determines the best estimate of α.

We shall now derive simple expressions for the biases B_{1}(α), B_{2}(α) and B∗_{2}(α) for a special family of variance functions: the power variance function defined by β′ = β^{p}/c, where c > 0 and we assume that β > 0 for p ≠ 0. Here, p = 0, 1, 2, 3 for the normal, Poisson, gamma and inverse Gaussian distributions, respectively. Other values of p define a number of distributions which have been classified by Jørgensen (1987). Cases with power 0 < p < 1 do not correspond to distributions in (12).

For the power variance function we get B_{1}(α) = −p/(2β), B_{2}(α) = p(p − 6)(p + 1)β^{p−3}/(24c) and B∗_{2}(α) = p(p − 3)(p − 2)β^{p−3}/(24c). Note that B_{1}(α) > 0 for p < 0 and B_{1}(α) < 0 for p ≥ 1. Also, B_{2}(α) is positive for −1 < p < 0 or p > 6 and negative for p < −1 or 1 < p < 6. Clearly, p = 0 yields B_{1}(α) = B_{2}(α) = 0, as expected, and p = −1 and p = 6 lead to B_{2}(α) = 0. B∗_{2}(α) vanishes when p = 0, 2, 3, being positive for 1 < p < 2 or p > 3 and negative for p < 0 or 2 < p < 3. A check of the expressions for B_{1}(α) and B_{2}(α) is provided for the Poisson and gamma distributions. For the Poisson distribution (p = 1) α = log(θ), c = 1 and d(y) = −y (since β = θ > 0). By Taylor series expansion we can find E(α̂ − α) = −1/(2nθ) − 5/(12n^{2}θ^{2}) + O(n^{−3}), which agrees with the formulas of B_{1}(α) and B_{2}(α) for p = 1. For the gamma distribution with known scale parameter k (p = 2) α = −θ, c = k and d(y) = −y (since β = k/θ > 0). Direct expansion yields E(α̂ − α) = α/(nk) + α/(n^{2}k^{2}) + O(n^{−2}), in agreement with the values of B_{1}(α) and B_{2}(α) for p = 2.

For the power variance function, the n^{−2} term in the variance of α̂ reduces to V_{2}(α) = p(p + 2)/(2β^{2}), which vanishes for p = 0 and p = −2. The mean squared errors can be expressed in the forms

MSE(α̂_{1}) = MSE(α̂_{2}) = MSE(α̃_{2}) = \frac{c}{nβ^{p}} + \frac{p^{2}}{2n^{2}β^{2}} + O(n^{−3})

and

MSE(α̂) = \frac{c}{nβ^{p}} + \frac{p(3p + 4)}{4n^{2}β^{2}} + O(n^{−3}).

Therefore,

MSE(α̂) − MSE(α̃_{2}) = \frac{p(p + 4)}{4n^{2}β^{2}} + O(n^{−3}).

This equation shows that, to order n^{−2}, MSE(α̂) > MSE(α̃_{2}) for p < −4 or p > 1, and MSE(α̂) < MSE(α̃_{2}) for −4 < p < 0. For p = −4 or p = 0, MSE(α̂) = MSE(α̃_{2}).

Bias-corrected maximum likelihood estimates for other families of variance functions can also be obtained using our results. For example, one can use the results in this paper to obtain the second and third order biases of the maximum likelihood estimate for quadratic (Morris, 1982), cubic (Letac and Mora, 1990) and Babel variance functions.

Acknowledgements

The partial financial support of CNPq/Brazil is gratefully acknowledged.

References

Abell, M.L. & Braselton, J.P. (1994). The Maple V Handbook. New York: AP Professional. Bickel, P.J. & Doksum, K.A. (1977). Mathematical Statistics: Basic Ideas and Selected Topics.

Oakland: Holden-Day. Jørgensen, B. (1987). Exponential dispersion models. Journal of the Royal Statistical Society B

49, 127-162. Lawley, D.N. (1956). A general method for approximating to the distribution of the likelihood

ratio criteria. Biometrika 71, 233-244. Letac, G. & Mora, M. (1990). Natural real exponential families with cubic variance function.

Annals of Statistics 18, 1-37. Morris, C.N. (1982). Natural exponential families with quadratic variance functions. Annals of

Statistics 10, 65-80. Shenton, L.R. & Bowman, K. (1977). Maximum Likelihood Estimation in Small Samples. New

York: MacMillan. Wolfram, S. (1991). Mathematica: A System for Doing Mathematics by Computer, 2nd ed. New

York: Addison-Wesley.