Unbiased Estimating Equation on Inverse Divergence and Its Conditions
Masahiro Kobayashi Kazuho Watanabe

 





Unbiased Estimating Equation on Inverse Divergence and Its Conditions
    
  Masahiro Kobayashi
  Information and Media Center

                    Toyohashi University of Technology

                    Toyohashi, Aichi, Japan

                    Email: kobayashi@imc.tut.ac.jp
  Kazuho Watanabe
  Department of Computer Science and Engineering
 
                    Toyohashi University of Technology

                    Toyohashi, Aichi, Japan

                    Email: wkazuho@cs.tut.ac.jp

    April 29, 2024
============================================================================================================================================================================================================================================================================================================================================================================================================================================




This paper focuses on the Bregman divergence defined by the reciprocal function, called the inverse divergence.
For the loss function defined by the monotonically increasing function f and inverse divergence, the conditions for the statistical model and function f under which the estimating equation is unbiased are clarified. 
Specifically, we characterize two types of statistical models, an inverse Gaussian type and a mixture of generalized inverse Gaussian type distributions, to show that the conditions for the function f are different for each model. 
We also define Bregman divergence as a linear sum over the dimensions of the inverse divergence and extend the results to the multi-dimensional case.


firstpage 



§ INTRODUCTION

The maximum likelihood estimation (MLE) is a standard method in parametric estimation although it is vulnerable to outliers. 
In robust statistics, the methods studied to overcome the adverse effects of outliers <cit.> include M-estimation, which is a well-known robust estimation method. 
In M-estimation, the assumed model in the MLE is changed to another with heavy tails. 
The well-known minimum divergence method changes the Kullback–Leibler divergence corresponding to the MLE to a more robust divergence for estimation <cit.>. 
The M-estimation and minimum divergence methods obtain estimators as solutions to estimating equations, which are mainly of two types: non-normalized and normalized.
Bregman divergence and its special cases <cit.> correspond to the non-normalized estimating equation, whereas γ-divergence <cit.> corresponds to the normalized estimating equation. 
Notably, the normalized estimating equation has the distinctive feature of bringing the latent bias close to zero, even in cases heavily contaminated by outliers <cit.>. 
However, the analytically intractable integrals involved in these estimating equations limit the choice of models and weight functions that can be used. 
Recent studies on the minimization of divergence have adopted the stochastic optimization framework to avoid these intractable integrals <cit.>.


The f-separable distortion measure <cit.> was proposed as an extension of the average and maximum distortions that are commonly used in information theory. 
Recently, we proposed a parameter estimation method that minimizes the f-separable distortion using Bregman divergence as a base distortion measure <cit.>.
This method is a type of M-estimation method.
The property of the estimator is determined by the shape of the monotonically increasing function f, exhibiting robustness against outliers when f is concave. 
However, extending the loss function using function f does not always guarantee an unbiased estimating equation, which is a necessary condition for the consistency of the estimator. 
To satisfy the unbiased estimating equation, a bias correction term involving an analytically intractable integral is required. 
Consequently, to avoid this problem, combinations of the statistical model, Bregman divergence, and function f that satisfy the unbiased estimating equation without a bias correction term have been investigated <cit.>.
It has been suggested that the necessity of a bias correction term depends on the type of Bregman divergence utilized in the estimation <cit.>. 
When a bias correction term is not required, the Bregman divergence describes the corresponding statistical model used for estimation, and the available functions f are characterized by the boundedness of a specific simple integral. 
In rare cases, the bias correction term disappears; however, the combinations generating this condition are unknown.

This paper focuses on the Bregman divergence defined by the reciprocal function, called the inverse divergence.
We clarify combinations of the statistical model and function f that eliminate the bias correction term when using inverse divergence for estimation.
Furthermore, we extend the result to a multi-dimensional  case by expressing the inverse divergence as a linear sum in multiple dimensions.









§ F-SEPARABLE BREGMAN DISTORTION MEASURES



 §.§ Problem setting and loss function

In estimating the parameter θ∈Θ⊆ℝ^d of a statistical model p(x|θ) based on the given data x^n={x_1,⋯,x_n},  x_i=(x_i^(1),⋯,x_i^(d))^ T∈χ⊆ℝ^d, we assume that the true distribution can be realized by p(x|θ^*). 
When the expected value exists for the statistical model, we assume that θ=𝔼[X]=∫xp(x|θ)dx.
The loss function is defined as

    L(θ) = 1/n∑_i=1^n f(d_ϕ(x_i, θ)),

using differentiable and continuous monotonically increasing function f:ℝ_+→ℝ and the Bregman divergence d_ϕ(x,θ):χ×Θ→ℝ_+, where ℝ_+ is the set of nonnegative real numbers <cit.>.
The Bregman divergence is defined as

    d_ϕ(x, θ) = ϕ(x) - ϕ(θ) - ⟨∇ϕ(θ), x-θ⟩,

where ϕ:χ→ℝ  is a differentiable strictly convex function with gradient vector ∇ϕ and ⟨·, ·⟩ denotes the inner product. 
The property of the estimator depends on the shape of the function f. 
When the function f is concave, the estimator is robust against outliers. 
Furthermore, when the function f is linear, the estimation problem reduces to the MLE for the expected value parameter in a regular exponential family:

    p(x|θ) = r_ϕ(x)exp(-d_ϕ(x, θ)),

where the strictly convex function ϕ uniquely determines r_ϕ(x) <cit.>.




 §.§ Conditions for unbiased estimating equation

The estimator θ̂ is the solution to the stationary point, obtained by differentiating the loss function (<ref>) with respect to the parameter θ and setting it to 0, in solving the estimating equation:

    1/n∑_i=1^n f'(d_ϕ(x_i, θ))(x_i - θ) = 0,

where f' is the derivative of f. 
However, this estimating equation is generally biased. 
This is because the left-hand side of the estimating equation (<ref>) does not necessarily converge asymptotically to 0. 
In the limiting value, bias should be pre-corrected to satisfy the unbiased estimating equation. 
Two main types of unbiased estimating equations are known: non-normalized (<ref>) and normalized (<ref>). 
The latter is particularly noted for its ability to reduce latent biases to 0, even in situations with a high proportion of outliers:

    1/n∑_i=1^n f'(d_ϕ(x_i, θ))(x_i - θ) = 𝔼_p(x|θ)[f'(d_ϕ(X, θ))(X-θ)], 
    ∑_i=1^n f'(d_ϕ(x_i, θ))(x_i-θ)/∑_j=1^n f'(d_ϕ(x_j, θ)) = 𝔼_p(x|θ)[f'(d_ϕ(X, θ))(X-θ)]/𝔼_p(x|θ)[f'(d_ϕ(X, θ))],

<cit.>.
The bias correction term on the right hand side of each estimating equation is generally difficult to integrate analytically. 
However, if the following condition holds, the bias correction terms on the right sides of both the non-normalized (<ref>) and normalized (<ref>) estimating equations vanish, resulting in the same problem expressed by the estimating equation (<ref>):

    𝔼_p(x|θ)[f'(d_ϕ(X, θ))(X-θ)] = 0.

In other words, if (<ref>) is satisfied, (<ref>) is the unbiased estimating equation not requiring a bias correction term.
Additionally, as it can be interpreted as the normalized estimating equation (<ref>), the latent bias can be approximately reduced to 0, even in situations with a high proportion of outliers.




 §.§ Existing results

Table <ref> shows the combinations of Bregman divergence, statistical model, and function f for which the bias correction term vanishes, and (<ref>) is satisfied. 
Previous research has clarified the conditions for the statistical model, function f corresponding to the Mahalanobis distance, and one-dimensional Bregman divergence[1] <cit.>. 
The Mahalanobis distance and one-dimensional Bregman divergence correspond to elliptical and continuous Bregman distributions, respectively.
Elliptical distributions are a family of distributions defined by the generating function and the Mahalanobis distance <cit.>. 
They include the well-known Gaussian, Laplace, and  t- distributions as special cases.
Whereas, continuous Bregman distributions are a distribution family defined by the generating function g and a strictly convex function ϕ. 
The special cases of this family are the one-dimensional elliptical distribution and the Itakura–Saito (IS) distribution, which generalizes the gamma distribution (Table <ref>).


The loss function (<ref>) is an extension of the negative log-likelihood function under the regular exponential family (<ref>), enhanced by the function f. 
Therefore, it is preferable that statistical models that satisfy the unbiased estimating equation (<ref>) correspond to the regular exponential family. 
However, apart from the Gaussian and gamma distributions, no special cases of the regular exponential family have been identified as satisfying the unbiased estimating equation (<ref>). 
The following sections clarify that for estimations using inverse divergence, the statistical model satisfying the unbiased estimating equation (<ref>) corresponds to the inverse Gaussian distribution (and its generalized distribution family), which is a special case of the regular exponential family.


§ MODELS AND CONDITIONS UNDER INVERSE DIVERGENCE

This section discusses the Bregman divergence specified by the strictly convex function ϕ(x)=λ/x,  (x>0, λ>0) which is called the inverse divergence <cit.>, where λ is a nuisance parameter. 
The inverse divergence, defined as

    d_ Inv^λ(x, θ) ≜λ(x-θ)^2/θ^2x,

is known to correspond to the inverse Gaussian distribution <cit.>, <cit.> and is also a particular case of the β-divergence <cit.>. 
In the following subsections, we discuss the two types of statistical models satisfying the unbiased estimating equation (<ref>) when using inverse divergence for estimation and show that the conditions for the function f differ for each model.



 §.§ Inverse Gaussian type (IGT) distribution


    For x∈ℝ_++, the parameters θ∈Θ=ℝ_++ and λ∈ℝ_++, and a nonnegative generating function g:ℝ_+→ℝ_+, the IGT distribution is defined as follows:
    
    p(x|θ, λ)    = 1/C_ IGT√(λ/x^3)g(d_ Inv^λ(x, θ)), 
    
            C_ IGT   = ∫_0^∞1/√(t)g(t)dt,

    if the normalization constant C_ IGT exists.

Here, ℝ_++=ℝ_+∖{0}.
Note that the normalization constant C_ IGT is the same as that for the one-dimensional elliptical distribution <cit.>.
When the generating function is g(t)=exp(- t/2), the IGT distribution (<ref>) reduces to the inverse Gaussian distribution <cit.> as follows:

    p(x|θ, λ) = √(λ/2π x^3)exp(-λ (x-θ)^2/2θ^2x).

The expected value of the IGT distribution, 𝔼[X]=θ, is independent of the generating function g <cit.>. 
This fact can also be derived from the following Corollary <ref>, which is derived below.

    There exists an IGT distribution (<ref>) corresponding to the nonnegative generating function g, i.e., C_ IGT<∞.

 
    Under Assumption <ref>, the estimating equation without a bias correction term equivalently, (<ref>) holds if and only if
    
    ∫_0^∞ g(t)f'(t)1/√(t+1)dt<∞

    holds for the combination of the function f and the statistical model (<ref>).


    Substituting the inverse divergence (<ref>) and IGT distribution (<ref>) into the left-hand side of (<ref>), we have
    
    𝔼_p(x|θ, λ)[f'(d_ Inv^λ(X, θ))(X-θ)] 
    
        =   ∫_0^∞1/C_ IGT√(λ/x^3)g(d_ Inv^λ(x, θ))f'(d_ Inv^λ(x, θ))(x-θ)dx 
    ∝   ∫_0^θ1/√(x^3)g(d_ Inv^λ(x, θ))f'(d_ Inv^λ(x, θ))(x-θ)dx 
       + ∫_θ^∞1/√(x^3)g(d_ Inv^λ(x, θ))f'(d_ Inv^λ(x, θ))(x-θ)dx 
    
        =   θ/√(λ)∫_0^∞ g(t)f'(t) [1/√(t+4λ/θ)-1/√(t+4λ/θ)]dt = 0.

    We used integration by substitution, t=d_ Inv^λ(x, θ). 
    The details of the substitution integration from (<ref>) to the next line are provided in Appendix <ref>.
    Therefore, if the following integral exists for any a>0, the unbiased estimating equation (<ref>) holds without a bias correction term:
    
    I(a) ≜∫_0^∞ g(t)f'(t)1/√(t + a)dt<∞.

    Conversely, the above discussion also shows that 𝔼_p(x|θ, λ)[|f'(d_ Inv^λ(X,θ))(X-θ)|]∝ 2θ√(λ)^-1I(4λθ^-1).
    In other words, I(a)<∞ is also a necessary condition.

    However, since this integral includes the parameter a=4λ/θ>0, rewriting it in a form that does not depend on a is desirable. 
    The integrand of (<ref>) is a strictly monotonically decreasing and continuous function with respect to a. 
    Therefore, the I(a) obtained through integration is also a strictly monotonically decreasing and continuous function with respect to a.
    We assume that I(a^†)<∞ for ∃ a^†∈ (0,∞).
    The strict monotonic decrease of I(a) ensures that ∀ε>0, I(a^† +ε)<∞ holds. 
    Furthermore, owing to the continuity of I(a), ∀ε∈(0,a^†), I(a^† - ε)<∞ holds. 
    Therefore, if I(a) is bounded at some point a^†∈(0,∞), it is bounded for any a∈(0,∞).
    Thus, we set a=1.
    Based on the above discussion, (<ref>) is a necessary and sufficient condition for the unbiased estimating equation to hold without the bias correction term.
    


The condition ∫_0^∞ g(t)√(t+1)^-1dt<∞ for the existence of the expected value of the IGT distribution is obtained from Theorem <ref> by substituting f'(t)=1.
As this condition is the lower bound of the normalization constant (<ref>) of the IGT distribution, the following corollary is obtained.

    The expected value of the IGT distribution (<ref>) always exists, independent of the generating function g, and satisfies 𝔼[X]=θ.


In the above discussion, we showed that for estimation using the inverse divergence (<ref>), the statistical model satisfying the unbiased estimating equation (<ref>) corresponds to the IGT distribution.
The inverse Gaussian distribution, a special case of the IGT distribution, is also a special case of the regular exponential family of distributions.
On the other hand, continuous Bregman distributions correspond to the one-dimensional Bregman divergence. 
Thus, the inverse divergence also corresponds to continuous Bregman distributions generated by the reciprocal function. 
In the following, we show that generalized inverse Gaussian type (GIGT) mixture distributions, which are special cases of continuous Bregman distributions, can be generated and correspond to the inverse divergence.



 §.§ Generalized IGT (GIGT) mixture distribution


    For x∈ℝ_++, the parameters θ∈Θ=ℝ_++, λ∈ℝ_++, ν∈ℝ, and a nonnegative generating function g:ℝ_+→ℝ_+, the GIGT distribution is defined as follows:
    
    q(x|θ, λ, ν)    = 1/C_ GIGT(θ, λ, ν)x^ν-1g(d_ Inv^λ (x,θ)), 
    
            C_ GIGT(θ, λ, ν)    = ∫_0^∞ t^ν-1g(d_ Inv^λ (t,θ))dt,

    if the normalization constant C_ GIGT exists.

When ν=-1/2, the GIGT distribution (<ref>) reduces to the IGT distribution (<ref>).
In continuous Bregman distributions, when the strictly convex function is set as ϕ(x)=λ/x,   (x>0, λ>0), the corresponding statistical model is defined by a two-component mixture of GIGT distribution (<ref>) as

    p(x|θ, λ) = w q(x|θ, λ, 0) + (1-w)q(x|θ, λ, -1), 
       w = C_ GIGT(θ, λ, 0)/C_ GIGT(θ, λ, 0)+θ C_ GIGT(θ, λ, -1).

Specifically, when the generating function is g(t)=exp(-t/2), the GIGT distribution (<ref>) reduces to the generalized inverse Gaussian (GIG) distribution (<ref>) <cit.>, and the GIGT mixture distribution (<ref>) becomes the GIG mixture distribution (<ref>):

    p(x|θ, λ) = w p(x|α, θ, 0) + (1-w) p(x|α, θ, -1), 
       p(x|α, η, ν) = η^-νx^ν-1/2K_ν(α)exp(-α/2(x/η+η/x)), 
       w = K_0(α)/K_0(α)+K_-1(α),  α = λ/θ,

where K_ν(·) represents the modified Bessel function of the third kind with index λ.

The combination of the one-dimensional Bregman divergence and continuous Bregman distribution has been shown to satisfy the unbiased estimating equation without a bias correction term. 
Therefore, the following corollary can be obtained from <cit.>.

    There exists a GIGT mixture distribution (<ref>) corresponding to the nonnegative generating function g, i.e., C_ GIGT(θ, λ, 0)<∞ and C_ GIGT(θ, λ, -1)<∞.


    Under Assumption <ref>, the estimating equation without a bias correction term equivalently, (<ref>) holds if and only if
    
    ∫_0^∞ g(t)f'(t)dt < ∞

    holds for the combination of the function f and the statistical model (<ref>).





 §.§ Discussion

In this section, we showed that when using inverse divergence for estimation, two types of statistical models satisfy the unbiased estimating equation: the IGT (<ref>) and GIGT mixture (<ref>) distributions, each with different conditions for the function f. 
Comparing conditions (<ref>) and (<ref>), clearly, the factor 1/√(t+1) applied to the integrand differs. 
This implies that with a fixed generating function g, the functions f available for estimation differ.
When fixing f'(t)=1, the IGT distribution always has an expected value independent of the generating function g. 
In contrast, the existence of the expected value in the GIGT mixture distribution depends on the generating function g.
More generally, the same can be concluded  when comparing the IGT distribution with the continuous Bregman distributions, which generalize the GIGT mixture distribution.
For distributions such as continuous Bregman, elliptical, and IS distributions that are generalized by the generating function g, the existence of the expected value depends on the generating function g. 
For example, the Cauchy distribution, a particular case of elliptical distributions, does not have an expected value. 
This highlights that the IGT distribution is somewhat unique within the family of distributions generalized by the generating function g.


§ EXTENSION TO MULTI-DIMENSIONAL CASE

In this section, we extend the problem to cases of multi-dimensional data points.
In the following, to simplify the notation, the j-th dimension variable is represented by a subscript j.
Let the Bregman divergence be given as a linear sum over the dimensions of the inverse divergence:

    d_ MInv^λ(x, θ) ≜∑_j=1^d  d_ Inv^λ_j(x_j, θ_j) .

The corresponding strictly convex function is given by ϕ(x)=∑_j=1^d λ_j/x_j.

    For x∈ℝ_++^d, the parameters θ∈Θ=ℝ_++^d, λ∈ℝ_++^d, and a nonnegative generating function g:ℝ_+→ℝ_+, the MIGT distribution is defined as follows if the normalization constant C_ MIGT exists:
    
    p(x|θ, λ)    = 1/C_ MIGT∏_j=1^d [√(λ_j/x_j^3)]g(d_ MInv^λ(x, θ)), 
    
            C_ MIGT   = π^d/2/Γ(d/2)∫_0^∞ g(t)t^d-2/2dt,

    where Γ(·) represents the gamma function.

The MIGT distribution has the same normalization constant C_ MIGT as elliptical distributions <cit.> and extends the IGT distribution to multiple dimensions. 
Setting the generating function as g(t)=exp(-t/2),  we obtain the simultaneous distribution of d independent inverse Gaussian distributions.
As shown in Corollary <ref>, the expected value of the MIGT distribution exists independent of the generating function g, and 𝔼[X]=θ holds.

    There exists a MIGT distribution (<ref>) corresponding to the nonnegative generating function g, i.e., C_ MIGT<∞.

 
    Let us assume d≥ 2.
    Under Assumption <ref>, the estimating equation without a bias correction term equivalently, (<ref>) holds if and only if 
    
    ∫_0^∞∫_0^∞ g(t+s)f'(t+s)t^d-3/2/√(s+1)dtds<∞

    holds for the combination of the function f and the statistical model (<ref>).

The proof of Theorem <ref> is provided in Appendix <ref>.

    The expected value of a MIGT distribution (<ref>) exists, independent of the generating function g, and satisfies 𝔼[X]=θ.


    From Theorem <ref>, by setting f'(t)=1, we immediately obtain the following relationship:
    
    ∫_ℝ_+^2 g(t+s)t^d-3/2/√(s+1)dtds<∞𝔼_p(x|θ, λ)[X]=θ<∞.

    We consider the upper bound of left condition of (<ref>):
    
    ∫_ℝ_+^2 g(t+s)t^d-3/2/√(s+1)dtds
                < ∫_ℝ_+^2 g(t+s)t^d-3/2/√(s)dtds 
       =√(π)Γ(d-1/2)/Γ(d/2)∫_0^∞ g(t)t^d-2/2dt 
       ≤π^d/2/Γ(d/2)∫_0^∞ g(t)t^d-2/2dt.

    The first inequality is derived from a simple comparison of the integrand. 
    The second line is obtained by setting  m=2, α_1=1/2, α_2=(d-1)/2, u(t)=g(t) and applying Lemma <ref> in Appendix <ref> to the first line. 
    Equation (<ref>) is the normalization constant (<ref>) of the MIGT distribution.
    Therefore, the expected value of the MIGT distribution exists independent of the generating function g, and 𝔼[X]=θ holds.
    



§ CONCLUSION

In this paper, we discussed the conditions under which the unbiased estimating equation holds without the bias correction term for loss functions composed of a monotonically increasing function f and Bregman divergence. 
In the case of inverse divergence, this scenario was satisfied by the IGT and GIGT mixture distributions (the latter is a particular case of the continuous Bregman distribution), each with different conditions of the function f.
In estimating the IGT distribution based on robust divergence, the bias correction terms, which are analytically intractable integrals, do not vanish. 
By defining the Bregman divergence as a linear sum over the dimensions of inverse divergence, we extended the discussion to multi-dimensional cases. 
The corresponding statistical model in this case (the MIGT distribution) was newly defined in this paper. 
Similar to the IGT distribution, we clarified that if a probability density function exists, an expected value, satisfying 𝔼[X]=θ  exists as well. 
Furthermore, we showed that the conditions for the function f that can be used for estimation are provided by a double integral.

When the statistical model is the regular exponential family, the Bregman divergences ensuring the unbiased estimating equation are the squared and IS distances reported in previous studies and the inverse divergence reported in this study,  which are special cases of β-divergence <cit.>. 
Although this paper did not delve into the details of robust estimation, by combining discussions from existing research, it was proven that latent biases can be reduced to approximately zero, even in the presence of outliers, and that the consistency of the estimator can be established <cit.>.



§ ACKNOWLEDGMENT

Except the abstract, this paper was translated from Japanese into English by OpenAI ChatGPT-4.
We would like to thank Editage and Enago for English language editing.
This work was supported in part by JSPS KAKENHI Grant numbers JP23K16849 and JP19K11825.







-1.2cm 
































IEEEtran






§ DETAILED PROOF OF THEOREM <REF>

Here, we set

    t = d_ Inv^λ(x, θ)

and perform substitution integration, which requires that the integration interval for x be divided.
The integration intervals for x from 0 to θ and from θ to ∞ are transformed into the integration intervals for t from ∞ to 0 and from 0 to ∞, respectively.
The factor of the integral transform is given by

    dx = λ^-1θ^2x^2/(x-θ)(x+θ)dt.

In addition, expanding (<ref>) with respect to x, we obtain

    λ x^2-θ(θ t + 2λ)x+λθ^2=0 .

From the quadratic formula, we obtain

    x = θ/2λ((θ t+2λ)±√(θ t(θ t + 4λ))).

For any θ, let us denote the range of x_<θ as (0,θ) and the range of x_≥θ as [θ, ∞).
Then, x_<θ and x_≥θ can be expressed as the inverse function of t in (<ref>) as follows:

    x_<θ(t) = θ/2λ((θ t+2λ)- √(θ t(θ t + 4λ))), 
    
        x_≥θ(t) = θ/2λ((θ t+2λ)+ √(θ t(θ t + 4λ))) .

Although these are functions of t, in the following discussion, t may occasionally be omitted for simplicity.
Let us define the following function:

    h(x) = √(x)/x+θ.

Substituting (<ref>) and (<ref>) into this function, we obtain

    h(x_<θ(t)) = h(x_≥θ(t)) = √(λ)/θ1/√(t+4λ/θ) .

Using the previous equations (<ref>), (<ref>), and (<ref>), the substitution integral of (<ref>) is calculated as

    ∫_∞^0 1/√(x_<θ^3) g(t)f'(t)(x_<θ-θ)λ^-1θ^2x_<θ^2/(x_<θ-θ)(x_<θ+θ)dt 
       +∫_0^∞1/√(x_≥θ^3) g(t)f'(t)(x_≥θ-θ)λ^-1θ^2x_≥θ^2/(x_≥θ-θ)(x_≥θ+θ)dt 
    
        =   θ^2/λ∫_0^∞ g(t)f'(t) [h(x_≥θ(t))-h(x_<θ(t)) ]dt 
    
        =   θ/√(λ)∫_0^∞ g(t)f'(t) [1/√(t+4λ/θ)-1/√(t+4λ/θ)]dt = 0.




§ PROOF OF THEOREM <REF>
 
    Substituting the multivariate inverse divergence (<ref>) and MIGT distribution (<ref>) into the left-hand side of (<ref>), we have
    
    𝔼_p(x|θ, λ)[f'(d_ MInv^λ(X, θ))(X-θ)] 
    ∝   ∫_ℝ_+^d∏_j=1^d [1/√(x_j^3)]g(d_ MInv^λ(x, θ))f'(d_ MInv^λ(x, θ))(x-θ)dx.

    In the following, we set g̅(t)=g(t)f'(t) and focus our discussion on the k-th dimension:
    
    ∫_ℝ_+^d∏_j=1^d [1/√(x_j^3)] g̅(d_ MInv^λ(x, θ))(x_k-θ_k)dx
    
            =   ∫_ℝ_+^d-1∫_0^θ_k1/√(x_k^3)g̅(d_ MInv^λ(x, θ))(x_k-θ_k)dx_k ∏_j≠ k^d 1/√(x_j^3)dx_j  
       +∫_ℝ_+^d-1∫_θ_k^∞1/√(x_k^3)g̅(d_ MInv^λ(x, θ))(x_k-θ_k)dx_k ∏_j≠ k^d 1/√(x_j^3)dx_j 
    
            =   θ_k/√(λ_k)∫_ℝ_+^dg̅(t_k+∑_j≠ k^d  d_ Inv^λ_j(x_j, θ_j)) 
       ·[1/√(t_k+4λ_k/θ_k)-1/√(t_k+4λ_k/θ_k)]dt_k ∏_j≠ k^d 1/√(x_j^3)dx_j = 0.

    The transformation from the first equality to the second involves setting t_k= d_ Inv^λ_k(x_k, θ_k) and applying substitution integration, similar to the proof in Theorem <ref>.
    For (<ref>) to be zero, the following expression must be bounded:
    
    ∫_ℝ_+^dg̅(t_k+∑_j≠ k^d d_ Inv^λ_j(x_j, θ_j))1/√(t_k+4λ_k/θ_k)dt_k ∏_j≠ k^d 1/√(x_j^3)dx_j .

    Iteratively integrating the above expression, we obtain
    
    ∏_j≠ k^d[1/√(λ_j)]∫_ℝ_+^dg̅(∑_j=1^d t_j)1/√(t_k+4λ_k/θ_k)dt_k ∏_j≠ k^d 1/√(t_j)dt_j 
    
                =∏_j≠ k^d[1/√(λ_j)]π^d-1/2/Γ(d-1/2)∫_0^∞∫_0^∞g̅(t+t_k)t^d-3/2/√(t_k+4λ_k/θ_k)dtdt_k .

    The transformation of the equation involves setting m=d-1 and u(t)=g̅(t+t_k), and using Lemma <ref> in Appendix <ref>.
    For the same reasons as in the proof of Theorem <ref>, the parameters λ_k and θ_k can be eliminated from the integrand.
    Note that although the proof focuses on the k-th dimension, the bounded condition of the integral that must to be satisfied in each dimension is the same and independent of the dimension.

    Therefore, if (<ref>) exists, the unbiased estimating equation (<ref>) holds without a bias correction term.
    Conversely, the above discussion also shows that 𝔼_p(x|θ, λ)[|f'(d_ MInv^λ(X,θ))(X_k-θ_k)|]∝ 2×(<ref>).
    In other words, (<ref>) is also a necessary condition.
    




§ LEMMAS


Given a positive integer m≥1 and α_j>0, where (j=1,⋯,m), for a nonnegative function u, the following relation <cit.> holds:

    ∫_ℝ_+^m u(∑_j=1^m t_j) ∏_j=1^m t_j^α_j-1dt_j 
       =∏_j=1^mΓ(α_j)/Γ(∑_j=1^mα_j)∫_0^∞ u(t) t^∑_j=1^mα_j-1dt.



In Lemma <ref>, by setting α_j=1/2, (j=1,⋯,m), we obtain the following lemma.

For a positive integer m≥1 and a nonnegative function u, the following relation holds:

    ∫_ℝ_+^m u(∑_j=1^m t_j) ∏_j=1^m1/√(t_j)dt_j 
       =π^m/2/Γ(m/2)∫_0^∞ u(t) t^m/2-1dt.



