Bias correction of posterior means using MCMC outputs
Yukito Iba


 

















The Institute of Statistical Mathematics,
10-3 Midori cho, Tachikawa City, Tokyo 190-8562, Japan. 
iba@ism.ac.jp
 
§ ABSTRACT 
We propose algorithms for addressing the bias of the posterior mean when used as an estimator of parameters. These algorithms build upon the recently proposed Bayesian infinitesimal jackknife approximation (<cit.>) and can be implemented using the posterior covariance and third-order combined cumulants easily calculated from MCMC outputs. Two algorithms are introduced: The first algorithm utilises the output of a single-run MCMC with the original likelihood and prior to estimate the bias. A notable feature of the algorithm is that its ability to estimate definitional bias (<cit.>), which is crucial for Bayesian estimators. The second algorithm is designed for high-dimensional and sparse data settings, where “quasi-prior” for bias correction is introduced. The quasi-prior is iteratively refined using the output of the first algorithm as a measure of the residual bias at each step. These algorithms have been successfully implemented and tested for parameter estimation in the Weibull distribution and logistic regression in moderately high-dimensional settings.


Bias correction of posterior means using MCMC outputs
    Yukito Iba
    September 9, 2024
=====================================================


 
Keywords. 
 
Bayesian statistics;
Markov Chain Monte Carlo;
bias correction;
infinitesimal jackknife;
posterior cumulant
 




§ INTRODUCTION



This paper presents algorithms for estimating and correcting the bias of the posterior mean when used as an estimator of parameters. We discuss two algorithms, both of which are based on posterior covariance and third-order combined cumulants.

The first algorithm utilizes the output of a single run of MCMC with the original likelihood and prior. A simple formula can be applied to any likelihood and prior, and the estimator is automatically computed from the posterior samples; no model-specific analytical calculation is required.

The second algorithm is designed for high-dimensional and sparse data settings, where a “quasi-prior” for bias correction is introduced. The quasi-prior is improved iteratively using the output of the first algorithm as a measure of the residual bias at each step. Typically, after several MCMC runs, the bias correction by the first algorithm yields satisfactory results; it is important to note that we do not need to entirely remove the bias through this iteration. Experiments show that the second algorithm is effective, for example, in models with 20 or 60 parameters.

These algorithms heavily rely on the recently proposed idea of the Bayesian infinitesimal jackknife approximation (Bayesian IJK, <cit.>). While frequentist covariance is represented by the posterior covariance in <cit.>, our algorithm estimates bias using posterior cumulants; we also refer to a related idea suggested in <cit.>. The proposed algorithms, however, accommodate the following novel features:

    
  * The definitional bias (<cit.>) is included in the proposed algorithms; this is essential for the Bayesian estimators. 
 
    
  * Iterative tuning of the quasi-prior is introduced in the second algorithm. This remarkably improves the performance of the algorithm in high-dimensional and sparse-data settings.


Details of (a) and (b) are discussed in the remainder of the paper. Here we provide a simple example of the definitional bias. Consider a binomial likelihood with a beta prior Beta(α,β) on the success probability q, and the posterior mean 
[q]= 
(X+α)/(n+α+β)
as an estimator of q, where n and X is the number of the trials and successes, respectively. In this context, the definitional bias b_0 is expressed as
b_0
= 
(n[X]+α)/(n+α+β)
-q_0
= 
(nq_0+α)/(n+α+β)
-q_0,
where [X] denotes the population average of X. Since the estimator [q] is linear in X, the bias caused by the non-linearity of the estimator is zero. Definitional bias b_0 is often disregarded in discussions based on the von Mises expansion or nonparametric bootstrap; in many theoretical  analyses, the “bias” of an estimator is defined as the deviation from θ_0+b_0 rather than from the population value θ_0 of the paramete. However, numerical experiments indicate that definitional bias b_0 can be significantly large for posterior mean estimators even when seemingly noninformative priors are used. This contrasts with the case of the maximum likelihood estimator, where b_0 is always 0. In this study, we introduce a simple method to estimate b_0 using the posterior covariance. The proposed estimator may fail in sparse data settings, but this issue is addressed by the use of the quasi-prior and its iterative improvement, as employed in the second algorithm.  

The rest of the paper is organized as follows: In Sec. <ref>, we introduce the proposed algorithms. In sections <ref> and <ref>, illustrative examples are presented for these algorithms. In Sec. <ref>, the proposed algorithms are derived from the jackknife approximation to the bias. Sec. <ref> provides a summary and conclusion. Appendix <ref> defines the bias and definitional bias by means of the von Mises expansion and Appendix <ref> derives the jackknife bias correction formula. 



§ PROPOSED METHOD





 §.§ Settings and notation


Let us denote a sample of size n from the population G as X^n=(X_1,…, X_n).
The posterior distribution is defined by

    p(θ| X^n ) = 
        exp{∑_i=1^nℓ(X_i; θ)}p(θ) /∫exp{∑_i=1^nℓ(X_i;θ')}p(θ') dθ',

where p(θ) is a prior density on θ=(θ_1,⋯, θ_K) and ℓ(x;θ)=log p(x |θ) is log-likelihood of the model. 

Our objective is to correct the frequentist bias in the posterior mean defined by

    [A]  = ∫ A(θ)p(θ| X^n) dθ
 
of given statistics A(θ). We also define the posterior covariance and a third-order combined posterior cumulant of the arbitrary statistics A(θ), B(θ), and C(θ) as

    [A(θ),B(θ)]  =
    [(A(θ)-[A(θ)])( B(θ)-[B(θ)])],
    
       [A(θ),B(θ), C(θ)]  =
    [(A(θ)-[A(θ)])(B(θ)-[B(θ)])(C(θ)-[C(θ)]]).
 

  



 §.§ Algorithm 1


In the algorithm 1, the bias b(A;n) is estimated using the formula

    b̂(A;n)=  - ∑_i=1^n [A(θ), ℓ(X_i; θ)]+1/2∑_i=1^n [A(θ), ℓ(X_i; θ), ℓ(X_i; θ)].

In practice, the posterior covariance and cumulants in the above formula are estimated from posterior samples obtained from a single MCMC run. The resultant algorithm is summarized as follows:

 
 
 



 §.§ Algorithm 2


In the algorithm 2, simultaneous bias correction of all parameters is essential. For simplicity, here we restrict to the case of simultaneous bias correction of the original parameters (θ,⋯,θ_K), while general cases including an arbitrary statistics A(θ) may be dealt with a change of the parameter.  

As a basis of the algorithm,  we introduce a modified posterior distribution defined by

    p_λ(θ| X^n ) = 
        exp{∑_i=1^nℓ(X_i; θ)-∑_k=1^K λ_k θ_k }p(θ) /∫exp{∑_i=1^nℓ(X_i;θ')-∑_k=1^K λ_k θ'_k}p(θ') dθ',

where λ=(λ_1,⋯,λ_K) is a vector of constants. An average [   ], covariance [   ], and third-order combined cumulant [   ] with respect to the distribution (<ref>) are defined in a similar manner that [   ], [   ], and [   ] is defined, respectively. Using them, we define

    b̂(θ_k;n;λ)=  - ∑_i=1^n [θ_k, ℓ(X_i; θ)]+1/2∑_i=1^n [θ_k, ℓ(X_i; θ), ℓ(X_i; θ)] ,
    and    C(θ_k,θ_k^';n;λ)=[θ_k, θ_k^'].

Hereafter, the vector whose components is b̂(θ_k;n;λ) is denoted as b̂(θ;n;λ), while the matrix whose (k,k^') component is C(θ_k,θ_k^';n;λ) is expressed as C(θ;n;λ). 

The proposed procedure is an iterative improvement of the constants λ=(λ_1,⋯,λ_K), which intends to reduce the bias b̂(θ;n;λ) estimated in each step. Let us denote the value of λ at step m as λ^(m) and define λ^(m) as a vector that solves a linear equation 

    C(θ;n;λ^(m)) λ^(m) = 
        b̂(θ;n;λ^(m)).
 
Using λ^(m), the value of λ^(m) is updated as

    λ^(m+1)=λ^(m)+δ λ^(m),
   
where 0< δ≤ 1 is a constant in the algorithm, the value of which is chosen to prevent oscillation of the value of λ. Since C(θ;n;λ^(m)) is a  matrix that express the posterior covariance between parameters, the above equation has a unique solution, unless parameters are not identifiable or extremely correlated. 

An essential point is that we need not entirely remove the bias by the iteration procedure. Our strategy is to reduce the magnitude of the bias towards the range where the approximation (<ref>) is sufficiently accurate, then estimate it using (<ref>).


 

The resultant algorithm is summarized as follows:
 
 
 



 §.§ Definition of b̂_0 and b̂_2


Since the first and second term in (<ref>) have a considerably different nature, contribution of each term is checked separately in the following numerical experiments. Hereafter, we denote the first and second term in (<ref>) as

    b̂_0(A;n)=     - ∑_i=1^n [A(θ), ℓ(X_i; θ)], 
    b̂_2(A;n)=    1/2∑_i=1^n [A(θ), ℓ(X_i; θ), ℓ(X_i; θ)],
  
which can be abbreviated as b̂_0(A) and b̂_2(A), or even b̂_0 and b̂_2, respectively. By an abuse of notation, the corresponding components of b̂(A;n;λ) are also expressed by the same symbols, such as b̂_0 and b̂_2.

In appendix <ref>, we define two components b_0(A;n) and b_2(A;n) of the bias based on the von Mises expansion; b̂_0(A;n) and b̂_2(A;n) are considered as an estimator of b_0(A;n) and b_2(A;n), respectively.  



§ EXAMPLE: WEIBULL FITTING



Here we test the algorithm 1 for the parameter estimation of the Weibull distribution, the probability density of which is given by

    p^wb(x |λ, γ)=γ/λ (x/λ )^γ-1exp [- ( x/λ)^γ ],     x ≥ 0.
 
This likelihood is not an exponential family in the shape parameter γ. An improper prior uniform on [0,∞) is assumed for each of the parameters λ and γ. Experiments below are performed with 5000 sets of artificial data of sample size n=30; they are generated from the Weibull distribution with λ=1.0 and γ=0.8. 

The results are presented in Fig.<ref>. The first and second panel of Fig.<ref> correspond to the scale parameter λ and shape parameter γ, respectively. In each panel, from left to right, results for b̂_2, b̂_0, and b̂_0+b̂_2 are shown. The horizontal dotted line colored red indicates the true value of the bias estimated from the average over the sets of artificial data. 

Fig.<ref> indicates that the bias is mostly explained by b̂_0 in the case of the scale parameter λ, while it is predominantly due to b̂_2 in the case of the shape parameter γ. This exhibit the importance of both b̂_0 and b̂_2 in the estimation of the bias. Fig.<ref> also presents that the proposed estimator b̂_0+b̂_2 provides reasonable estimates in both cases, while a small but systematic deviation is observed in the case of shape parameter; it may be related to a highly nonlinear nature of the estimator of γ.






§ EXAMPLE: LOGISTIC REGRESSION



Here we consider the logistic regression in a moderately high-dimensional settings, where the algorithm 2 is effectively applied. Let us consider the  logistic regression for data Y^n=(Y_i),   Y_i ∈{0,1} as:

    p(Y^n;a)=∏_i=1^n q_i^Y_i(1-q_i)^1-Y_i,     logq_i/1-q_i= ∑_s=1^N_p a_s x_is,

where a_s,   s=1,⋯,N_p is regression coefficients and x_is,   s=1,⋯,N_p are explanatory variables corresponding to Y_i. We consider the following settings: 

    
  * The model comprises N_p=3m explanatory variables. Here, we consider the cases of N_p=60 and N_p=21.
    
  * Artificial data of sample size n are generated using the same model. The values of the regression coefficients (a_i) are 10, -5, and 0 in each of the three groups of the size m, respectively.
    
  * The value of the explanatory variables are randomly selected from the normal distribution of the variance 1/n and mean 0. When we consider the effect of confounding between explanatory variables, the values of the explanatory variables are sampled from a multivariate normal distribution; off-diagonal elements of the covariance matrix are set as ρ/n.
    
  * We select δ that defines a scale of increments of λ_is as 0.2 throughout the experiments; an adaptive choice of δ is importnat, as well as an adequate stopping criterion of the algorithm, which is left for the future study.


In figures <ref>–<ref>, medians of the estimated regression coefficients a_s,   s=1,⋯, N_p in 30 trials are plotted; each of these 30 trials employs a different set of Y_i and x_is. The true values of a_s is expressed as horizontal red lines. 

Fig. <ref> presents a result of algorithm 1, where the number of parameters is N_P=60 and the sample size n is set as 10× N_p=600. From left to right, the original estimates (posterior means), estimates corrected by b̂_2, b̂_0, and b̂_0+b̂_2 are presented; the rightmost panel corresponds to the output of algorithm 1. While the original estimates [a_s] in the leftmost panel have considerable bias, the algorithm 1 provides reasonably corrected values of a_s. In addition, results in the second and third panel indicate that both of b̂_2 and b̂_0 contribute the bias correction in algorithm 1.  

Fig. <ref> presents another result, where the sample size n reduces to 5× N_p=300 with a corresponding change of the variance of explanatory variables. The first row of Fig. <ref> presents the results of the algorithm 1. In contrast to the case n=600, the rightmost panel indicates that the proposed estimator b̂_0 +b̂_2 fails to correct the bias. Results in the second and third panel indicate that adequate corrections are not achieved even when we select one of the terms b̂_0 and b̂_2 and discard the other. Numerical attempts of evaluating the accuracy of b̂_0 and b̂_2 separately (not shown here) suggest that the majority of the error in a sparse-data case comes from b̂_0.

On the other hand, the second row of Fig. <ref> presents the result of algorithm 2 for the same sets of artificial data.  As shown in the rightmost panel shows, algorithm 2 successfully corrects the bias after only three iterations of updating λ_is using (<ref>) and (<ref>). We stress that the bias is not entirely removed even after updating λ_is; it is evident in the leftmost panel in the second row, where the posterior means defined with the updated values of λ_is are presented. This means that the role of the  bias correction terms -∑_k=1^K λ_k A_k(θ) are to reduce the bias and make the bias correction by algorithm 1 effective.  





Thus far, we assumed that the values of explanatory variables x_is are independently generated for each s; this assumption might hardly hold in a real-world problem. In Fig.<ref>, we present the results for artificial data of N_p=21 and n=5× Np=105, where the values of the explanatory variables x_is are correlated in the artificial sets of the data;  off-diagonal elements of the covariance matrix are set as 0.5/N_p. The results are essentially the same as those in Fig.<ref>: algorithm 1 fails to correct the bias as shown in the rightmost panel of the first row, while algorithm 2 successfully correct the bias as shown in the rightmost panel of the second row.





§ DERIVATION





 §.§ Settings and regularity conditions


The posterior distribution p_w(θ ;X^n) with weights (w_1,w_2,…,w_n) of the observations X^n=(X_1,X_2,… X_n) is defined by: 

    p_w(θ ;X^n) = 
        exp{∑_i=1^n w_i ℓ(X_i; θ)}p(θ) /∫exp{∑_i=1^n w_i ℓ(X_i;θ')}p(θ') dθ'.

The average of A(θ) over the distribution p_w(θ ;X^n) are expressed as ^w[A(θ)]. Hereafter w=1 is used as an abbreviation of w_i=1, i=1,…,n.

We also introduce the following notation for a fourth order combined cumulant: 

    [A(θ),B(θ)]  = 
       ^w  [(A(θ)-[A(θ)])( B(θ)-^w [B(θ)])^3 ] 
    
        -3 ×^w  [(A(θ)-^w[A(θ)])(B(θ)- ^w [B(θ)])  ]
    ×^w  [(B(θ)-^w[B(θ)])^2  ].


We assume the following conditions (C1) and (C2). Hereafter [   ] denotes a population average over an observation X_1;  𝒲^(-1) is the set defined by 0 ≤ w_1 ≤ 1 and w_i≠ 1=1. 


    
  (C1) 
    [
        sup_w∈𝒲^(-1)|
        [A(θ),ℓ(X_1;θ)]
        |
        ]
         =o(1/n^2)
    
  (C2) There exists an integrable function M(θ) such that
    for all w∈𝒲^(-1), 
    
    | p(θ)
        exp(w_1ℓ(X_1;θ)+
        ∑_k≠ 1ℓ(X_k;θ) )
        A^l(θ) {ℓ(X_1;θ) }^m |  ≤ M(θ),     0 ≤ l ≤ 1,    0 ≤ m ≤ 3.






 §.§ Local case sensitivity formulae


Let us introduce local case sensitivity formulae as a basic tool for the following derivation. The first-order local case sensitivity formula () represents the first-order derivative of ^w[A(θ)] using the posterior covariance as

    . ∂/∂ w_i^w[A(θ)]  |_w=1 
    =[A(θ),  ℓ(X_i;θ)].
     In a similar manner, special cases of second- and third- order formulae required here are expressed as   . ∂^2/∂ w_i^2^w[A(θ)]  |_w=1
    = [A(θ), ℓ(X_i;θ), ℓ(X_i;θ)].
     
       ∂^3/∂ w_i^3^w[A(θ)] 
    = [A(θ), ℓ(X_i;θ)],

where the last one is defined for an arbitrary 0 ≤ w ≤ 1. The proof of these formulae is straightforward when an exchange of integration and derivation is allowed under the condition (C2); it is presented in <cit.> in a generic form; (<ref>) is also proved in <cit.>.



 §.§ Bayesian IJK for the bias


The algorithm 1 is derived as a Bayesian IJK approximation to the jackkinfe bias correction formula:

    b_jack= (n-1)  (1/n∑_i=1^n ^-i[A] -[A]  )
        = ∑_i=1^n  (^-i[A] -[A]  ) +o_p(1/n).

This formula is well known, but we give a derivation in the appendix <ref>. In the jackknife method, the number of the observations are changed from n to n-1 by the deletion of an observation.  It is essential to deal with the difinitional bias. Seemingly more sophisticated sample-reuse approaches, such as the bootstrap method, may fail in dealing with the definitional bias, when the number of the observations is kept constant before and after disturbing the data.

Setting w_i=0 and w_j=1   (j≠ i) in the weighted posterior (<ref>) and consider a Taylor expansion with respect to w_i-1 of the posterior expectation, we obtain

    ^-i[A] = [A]
        -. ∂/∂ w_i^w[A(θ)]  |_w=1
        +1/2. ∂^2/∂ w_i^2^w[A(θ)]  |_w=1
        - 1/6∑_i=1^n . ∂^3/∂ w_i^3^w[A(θ)]  |_w=w*

where 0 ≤ w* ≤ 1. On the other hand, the expression (<ref>) and condition (C1) gives

    1/6∑_i=1^n . ∂^3/∂ w_i^3^w[A(θ)]  |_w=w* =o_p(1/n).


Substituting (<ref>) in (<ref>) and use (<ref>), b_jack is expressed as

    b_jack= - ∑_i=1^n . ∂/∂ w_i^w[A(θ)]  |_w=1
        +1/2∑_i=1^n . ∂^2/∂ w_i^2^w[A(θ)]  |_w=1
        + o_p(1/n).

Now that we use (<ref>)
and (<ref>) in the expansion (<ref>), which completes the derivation as

    b_jack= - ∑_i=1^n [A(θ),  ℓ(X_i;θ)]
        +1/2∑_i=1^n [A(θ), ℓ(X_i;θ), ℓ(X_i;θ)]
        + o_p(1/n).


 


 §.§ Derivation of the algorithm 2


The algorithm 2 is considered as an iterative improvement of the residual bias based on the following representation of the derivative

    ∂/∂λ_k[θ_k^']=
         -[θ_k, θ_k^'],

where [  ] denote the covariance with the distribution (<ref>). The expression (<ref>) is obtained by a direct computation when an exchange of integration and derivation is allowed under appropriate regularity conditions. The idea behind (<ref>) is the same as that of local case sensitivity formulae; applications of similar formulae have a long history in statistics (<cit.>). 

Let us define λ^* as a value of λ that gives zero bias for all parameters θ_k; such a value might not exist, but here we introduce it as a hypothetical target of the bias reduction procedure.

    b̂_k(θ;n;λ^*)=0.
 
If we define λ^(m)=λ^*-λ^(m), where λ^(m) is the value of λ at step m.
Then, we expand (<ref>) around λ^(m) as

    b̂_k(θ;n;λ^(m))
        +∑_k^'∂/∂λ_k^'[θ_k]
        λ^(m)_k^'
        + o(λ^(m))=0.
 
Using (<ref>) and ignore the residual terms, this is expressed as

    ∑_k^'[θ_k(θ), θ_k^'(θ)]
        λ^(m)_k^'=b̂_k(θ;n;λ^(m)).
 
In a vector form, this gives the equation (<ref>).



§ SUMMARY AND FUTURE PROBLEMS


We proposed algorithms for estimating and correcting the bias of the posterior mean as an estimator of parameters. Two algorithms are introduced, both of which rely on posterior covariance and cumulants derived from MCMC outputs, eliminating the need for additional analytical calculations. The first algorithm is based on a Bayesian infinitesimal jackknife approximation and successfully estimate the bias including the difinitional bias (<cit.>) using the result of a single MCMC run. The second algorithm involves an iterative improvement of a quasi-prior based on the output of the first algorithm; it is shown to be effective in high-dimensional and sparse settings for logistic regression.   

In addition to utilizing the recently emerging Bayesian infinitesimal jackknife approximation (<cit.>), the second algorithm is characterized by a hybrid approach that combines bias estimation and correction. An significant aspect of this strategy is that the “bias reduction terms” do not need to completely eliminate the bias. Instead, their role is to reduce the bias sufficiently so that the bias correction by the first algorithm becomes effective. However, the use of the first algorithm within the second algorithm requires justification of the first algorithm in cases where quasi-prior is non-negligible. Preliminary analysis suggests that the success of this approach may depend on the use of a quasi-prior whose logarithm is linear in θ, but further analysis and experimentation are required for the understanding.

There are a number of subjects left for future studies: First, more theoretical analysis is needed along with the development of an appropriate stopping criterion for the algorithm. Additionally, the relationship between the proposed method and existing approaches, such as the Firth method for bias correction of the maximum likelihood estimator (<cit.>), should be explored. Finally, applications to large-scale real-world data and further tests using artificial data are necessary to assess the potential and limitations of the proposed method.



§ ACKNOWLEDGEMENTS


I would like to thank Keisuke Yano for the fruitful discussions and valuable advice.





§ VON MISES EXPANSION



Here we introduce the von Mises expansion and provide an interpretation of the definitional bias b_0 in this context. Let us define a formal posterior for an arbitrary distribution F as

    p(θ ; F,n ) = 
        exp{n ∫ℓ(x; θ) dF(x) }p(θ) /∫exp{n ∫ℓ(x; θ') dF(x) }p(θ') dθ',

where F substitutes for the empirical distribution Ĝ_n in (<ref>); here as usual the empirical distribution is defined by the sum dĜ_n=(1/n)∑_i=1^n δ_X_i of the point measures δ_X_i concentrated on observations X^n=(X_1,X_2, …, X_n). 

An estimator T^A(F, n) of A(θ) is defined as a posterior mean as

    T^A(F, n)=∫ A(θ)p(θ;F,n) dθ.
 
Specifically, T^A(Ĝ_n,n)=[A(θ)]. When F=G, T^A(G,n) defines an “ideal value at sample size n” for the estimator; this does not necessarily coincides with the “true value” (or projection) defined as A_0=lim _n →∞T^A(G,n) and can have some bias even for a consistent estimator. 
 

Now that we introduce von Mises expansion (<cit.>) as

    T^A(Ĝ,n)=T^A(G,n)+1/n∑_i=1^n T^A_1(X_i;G)
    +1/2n^2∑_i=1^n ∑_j=1^n T^A_2(X_i,X_j;G)+o_p  (1/n ),

where T^A_1(X_i;G) and T^A_2(X_i,X_j;G) are influence functions for the true distribution G. These influence functions are assumed to satisfy

    [T^A_1(X_i;G)]=0,    _X_i[T^A_2(X_i,X_j;G,n)]=_X_j[T^A_2(X_i,X_j;G,n)]=0,

where _X_i means frequentist expectation over X_i, keeping other components of X^n fixed.  
The dependence of T^A_1 and T^A_2 on the sample size n is omitted here; we will find that it is crucial in the first term T^A(G,n) for our purpose of estimating the bias. If we take the frequentist expectation for both sides and  use (<ref>), it gives:

    [T^A(Ĝ,n)]= [T^A(G,n)]+1/2n^2∑_i=1^n [T^A_2(X_i,X_i;G)]+o  (1/n ).


If we restrict ourselves to a consistent estimator and regard lim_n →∞ T^A(G,n) as the true value of A, the bias of the estimator T^A(Ĝ,n) is given by:

    b_0(A,n)=[T^A(G,n)]-lim_n →∞ T^A(G,n) , 
        b_2(A,n)=1/2n^2∑_i=1^n [T^A_2(X_i,X_i;G)],
     
        b(A,n)=  [T^A(Ĝ,n)]-lim_n →∞ T^A(G,n)=b_0(A,n)+b_2(A,n)+o  (1/n ),

where b_0(A,n), b_2(A,n), and b(A,n) may be abbreviated as b_0, b_2, and b, respectively.

It is natural to estimate b_2 using b̂_2 defined by (<ref>), because the formula (<ref>) indicates that the second-order influence function T^A_2(X_i,X_i;G) can be represented by the third-order posterior cumulant. To be precise, we need to consider relations in (<ref>),

    _X_i[T^A_2(X_i,X_j;G,n)]=_X_j[T^A_2(X_i,X_j;G,n)]=0,

imposed on the second order-influence function. However, the corresponding correction appears to have little effect in examples. If we consider b̂_2 as an estimator of b_2, the rest part b̂_0 of the proposed estimator should be regarded as an estimator of b_0; we can confirm this in an example of the binomial likelihood. 



§ JACKKNIFE BIAS CORRECTION



First, we assume the bias b of the estimator Â(X^n) of statistics A is asymptotically proportional to the inverse of the sample size n as:

    [Â(X^n)]=A_0+c/n+r(n),

where A_0 is a “true value” of A_0 defined as a limit lim_n→∞Â(X^n); c is a constant and the residual term r(n) is assumed to be in the order of O(1/n^1+α), α>0.

Let us consider modified data where the observation i is removed and express it as Â(X^-i). Since sample size of X^-i is n-1 and each component comes from the same population as for X^n, (<ref>) indicates

    [Â(X^-i)]=A_0+c/n-1+r(n-1).

From (<ref>) and (<ref>), we obtain the desired result as 

    [ ∑_i=1^n  (Â(X^-i) - Â(X)  )  ]=n   (c/n-1-c/n ) +o(1/n)
        = c/n+o(1/n).



 
 
apalike

