Statistical Inference based on Bridge Divergences

TL;DR

This paper introduces a generalized divergence family that smoothly connects the density power divergence and the logarithmic density power divergence, enhancing robust statistical inference and clarifying their relationship.

Contribution

It develops a unified divergence framework bridging DPD and LDPD, resolving longstanding questions about their hierarchy and optimal selection in data analysis.

Findings

Provides a new family of divergences linking DPD and LDPD

Clarifies the hierarchy and relation between the two divergence families

Offers a tool for improved robust inference in statistical models

Abstract

M-estimators offer simple robust alternatives to the maximum likelihood estimator. Much of the robustness literature, however, has focused on the problems of location, location-scale and regression estimation rather than on estimation of general parameters. The density power divergence (DPD) and the logarithmic density power divergence (LDPD) measures provide two classes of competitive M-estimators (obtained from divergences) in general parametric models which contain the MLE as a special case. In each of these families, the robustness of the estimator is achieved through a density power down-weighting of outlying observations. Both the families have proved to be very useful tools in the area of robust inference. However, the relation and hierarchy between the minimum distance estimators of the two families are yet to be comprehensively studied or fully established. Given a particular set of real data, how does one choose an optimal member from the union of these two classes of divergences? In this paper, we present a generalized family of divergences incorporating the above two classes; this family provides a smooth bridge between the DPD and the LDPD measures. This family helps to clarify and settle several longstanding issues in the relation between the important families of DPD and LDPD, apart from being an important tool in different areas of statistical inference in its own right.

Tables

Table 1. Table 1: Bridge Divergence Estimators of σ 𝜎 \sigma with α = 0.8 𝛼 0.8 \alpha=0.8 . (The global minimizers are presented in the parentheses).

$\lambda$	${\hat{\theta }}_{\alpha ,\lambda }$	$\lambda$	${\hat{\theta }}_{\alpha ,\lambda }$	$\lambda$	${\hat{\theta }}_{\alpha ,\lambda }$
$0.9$	1.245808	$0.6$	1.248243	$0.3$	1.252942
	(1.4411e-5)		(1.4124e-5)		(1.4085e-5)
$0.8$	1.246477	$0.5$	1.249401	$0.2$	1.255756
	(1.4218e-5)		(1.4106e-5)		(1.4078e-5)
$0.7$	1.247269	$0.4$	1.250926	$0.1$	1.259926
	(1.4155e-5)		(1.4094e-5)		(1.4073e-5)

Equations

$${\rho }_1^{(\alpha )}(g,f_{\theta })=\int \left\{f_{\theta }^{1+\alpha }-\left(1+\frac{1}{\alpha }\right)g{f}_{\theta }^{\alpha }+\frac{1}{\alpha }{g}^{1+\alpha }\right\}.$$

$${\rho }_1^{(0)}(g,f_{\theta })=\int g\log \left(\frac{g}{f_{\theta }}\right),$$

$${\rho }_0^{(\alpha )}(g,f_{\theta })=\log \int f_{\theta }^{1+\alpha }-\left(1+\frac{1}{\alpha }\right)\log \int g{f}_{\theta }^{\alpha }+\frac{1}{\alpha }\log \int g^{1+\alpha }.$$

$${\rho }_{\varphi }^{(\alpha )}(g,f_{\theta })=\frac{1}{\varphi }{\left(\int f_{\theta }^{1+\alpha }\right)}^{\varphi }-\frac{1}{\varphi }\left(1+\frac{1}{\alpha }\right){\left(\int g{f}_{\theta }^{\alpha }\right)}^{\varphi }+\frac{1}{\alpha \varphi }{\left(\int g^{1+\alpha }\right)}^{\varphi },$$

$${\hat{\theta }}_{n1}^{\alpha }:=\underset{\theta \in \Theta }{\mathrm{argmin}}\left[\int f_{\theta }^{1+\alpha }-\left(1+\frac{1}{\alpha }\right)\frac{1}{n}\sum _{i=1}^n{f}_{\theta }^{\alpha }(X_i)\right],$$

$${\hat{\theta }}_{n0}^{\alpha }:=\underset{\theta \in \Theta }{\mathrm{argmin}}\left[\log \left(\int f_{\theta }^{1+\alpha }\right)-\left(1+\frac{1}{\alpha }\right)\log \left(\frac{1}{n}\sum _{i=1}^n{f}_{\theta }^{\alpha }(X_i)\right)\right].$$

$$\int f_{\theta }^{1+\alpha }{u}_{\theta }=\frac{1}{n}\sum _{i=1}^n{f}_{\theta }^{\alpha }(X_i)u_{\theta }(X_i),$$

$$\left(\frac{1}{n}\sum _{i=1}^n{f}_{\theta }^{\alpha }(X_i)\right)\int f_{\theta }^{1+\alpha }{u}_{\theta }=\left(\frac{1}{n}\sum _{i=1}^n{f}_{\theta }^{\alpha }(X_i)u_{\theta }(X_i)\right)\int f_{\theta }^{1+\alpha }.$$

$$t_1(\theta ):=\int f_{\theta }^{1+\alpha },t_2(\theta ):=\int g{f}_{\theta }^{\alpha },t_3(\theta ):=\int g^{1+\alpha }.$$

$$\int f_{\theta }^{1+\alpha }{u}_{\theta }=\int f_{\theta }^{\alpha }g{u}_{\theta },$$

$$\int f_{\theta }^{\alpha }g\int f_{\theta }^{\alpha }{u}_{\theta }=\int f_{\theta }^{\alpha }g{u}_{\theta }\int f_{\theta }^{\alpha }.$$

$$\left[\frac{\nabla t_1(\theta )}{\alpha +1}-\frac{\nabla t_2(\theta )}{\alpha }\right]$$

$$\left[\frac{t_2(\theta )\nabla t_1(\theta )}{\alpha +1}-\frac{t_1(\theta )\nabla t_2(\theta )}{\alpha }\right]$$

$$\lambda \left[\frac{t_1^{\prime }(\theta )}{\alpha +1}-\frac{t_2^{\prime }(\theta )}{\alpha }\right]+(1-\lambda )\left[\frac{t_2(\theta )t_1^{\prime }(\theta )}{\alpha +1}-\frac{t_1(\theta )t_2^{\prime }(\theta )}{\alpha }\right]=0.$$

$$\frac{t_1^{\prime }(\theta )}{1+\alpha }[\lambda +\bar{\lambda }{t}_2(\theta )]$$

$$\frac{t_1^{\prime }(\theta )}{\lambda +\bar{\lambda }{t}_1(\theta )}$$

$$\frac{1}{\bar{\lambda }}\log (\lambda +\bar{\lambda }{t}_1(\theta ))-\frac{1}{\bar{\lambda }}\left(\frac{1+\alpha }{\alpha }\right)\log (\lambda +\bar{\lambda }{t}_2(\theta ))+C,$$

$${\rho }^{(\alpha ,\lambda )}(g,f_{\theta })=\frac{1}{\bar{\lambda }}\log (\lambda +\bar{\lambda }{t}_1(\theta ))-\frac{1}{\bar{\lambda }}\left(\frac{1+\alpha }{\alpha }\right)\log (\lambda +\bar{\lambda }{t}_2(\theta ))+\frac{1}{\bar{\lambda }}\frac{1}{\alpha }\log (\lambda +\bar{\lambda }{t}_3(\theta )).$$

$$\frac{\int f_{\theta }^{1+\alpha }{u}_{\theta }}{\lambda +\bar{\lambda }\int f_{\theta }^{1+\alpha }}=\frac{\frac{1}{n}{\sum }_{i=1}^n{f}_{\theta }^{\alpha }(X_i)u_{\theta }(X_i)}{\lambda +\bar{\lambda }\frac{1}{n}{\sum }_{i=1}^n{f}_{\theta }^{\alpha }(X_i)}.$$

$${\hat{\theta }}_n^{(\alpha ,\lambda )}:=\underset{\theta \in \Theta }{\mathrm{argmin}}\frac{1}{\bar{\lambda }}\log (\lambda +\bar{\lambda }{t}_1(\theta ))-\frac{1}{\bar{\lambda }}\left(\frac{1+\alpha }{\alpha }\right)\log \left(\lambda +\bar{\lambda }\frac{1}{n}\sum _{i=1}^n{f}_{\theta }^{\alpha }(X_i)\right).$$

$$P_{i,\gamma }=\int g^i{f}_{\theta }^{\gamma },Q_{i,\gamma }=\int g^i{f}_{\theta }^{\gamma }{u}_{\theta }{u}_{\theta }^{\top },R_{i,\gamma }=\int g^i{f}_{\theta }^{\gamma }{u}_{\theta },S_{i,\gamma }=\int g^i{f}_{\theta }^{\gamma }\nabla u_{\theta }.$$

$$\underset{n\to \infty }{\lim }{f}_{{\theta }_n}(x)=f_{\theta }(x),$$

$$\underset{\theta \in \Theta }{\sup }\left|\frac{1}{n}\sum _{i=1}^n{f}_{\theta }^{\alpha }(X_i)-\int g{f}_{\theta }^{\alpha }\right|=o_p(1),(\text{mbox}followsfromuniformSLLN)$$

$$\Rightarrow$$

$$K(\theta )$$

$$-\bar{\lambda }\left[\lambda +\bar{\lambda }{P}_{0,\alpha +1}\right]R_{1,2\alpha }{R}_{0,\alpha +1}^{\top }+{\left[\lambda +\bar{\lambda }{P}_{0,\alpha +1}\right]}^2{Q}_{1,2\alpha }$$

$$-{\bar{\lambda }}^2{P}_{1,\alpha }^2{R}_{0,\alpha +1}{R}_{0,\alpha +1}^{\top }+\bar{\lambda }\left[\lambda +\bar{\lambda }{P}_{0,\alpha +1}\right]P_{1,\alpha }{R}_{1,\alpha }{R}_{0,\alpha +1}^{\top }$$

$$+\bar{\lambda }\left[\lambda +\bar{\lambda }{P}_{0,\alpha +1}\right]P_{1,\alpha }{R}_{0,\alpha +1}{R}_{1,\alpha }^{\top }-{\left[\lambda +\bar{\lambda }{P}_{0,\alpha +1}\right]}^2{R}_{1,\alpha }{R}_{1,\alpha }^{\top }$$

$$J(\theta )$$

$$-\alpha Q_{1,\alpha }\left[\lambda +\bar{\lambda }{P}_{0,\alpha +1}\right]-S_{1,\alpha }\left[\lambda +\bar{\lambda }{P}_{0,\alpha +1}\right]$$

$$+\bar{\lambda }\alpha R_{1,\alpha }{R}_{0,\alpha +1}^{\top }-\bar{\lambda }(\alpha +1)R_{0,\alpha +1}{R}_{1,\alpha }^{\top }.$$

$$d_{\lambda ,\alpha }(g,f)=\frac{1}{\bar{\lambda }}\frac{1}{1+\alpha }\log \left(\lambda +\bar{\lambda }\int f^{1+\alpha }\right)-\frac{1}{\alpha \bar{\lambda }}\log \left(\lambda +\bar{\lambda }\int g{f}^{\alpha }\right).$$

$$D_{\lambda ,\alpha }(g,f)=-d_{\lambda ,\alpha }(g,g)+d_{\lambda ,\alpha }(g,f).$$

$$d_{\lambda ,\alpha }(g,h)=d_{\lambda ,\alpha }(f,h)-\frac{1}{\alpha \bar{\lambda }}\log (1-\epsilon )+O(T_{ϵ,\delta }),$$

$$T_{\varepsilon,\delta}:=\frac{\varepsilon}{1-\varepsilon}\left[\lambda+\bar{\lambda}\int\delta h^{\alpha}\right]\bigg{/}\alpha\bar{\lambda}\left[\lambda+\bar{\lambda}\int fh^{\alpha}\right].$$

Full text

1 Introduction

Statistical procedures based on minimization of divergences are popular in the literature. In our context, a divergence is a distance like dissimilarity measure between two distributions which does not necessarily demand metric properties. Density based minimum divergence methods are very useful as they combine high efficiency with strong robustness properties. Basu et al., (1998) proposed the class of density power divergences (DPD). These divergences, when constructed between an arbitrary density $g$ and a parametric model density $f_{\theta}$ are indexed by a non-negative robustness tuning parameter $\alpha$ and have the form

[TABLE]

The divergence at $\alpha=0$ is defined as the continuous limit of the expression in Equation (1.1) as $\alpha\to 0$. This generates the divergence

[TABLE]

where $\log$ represents natural logarithm. Observe that the only term on the right hand side of Equation (1.1) containing both $g$ and $f_{\theta}$ has $g$ in degree one. A divergence with such a property has been referred to as a decomposable divergence by some authors; see, for example, Broniatowski et al., (2012).

Another class of divergences with some similar properties was introduced by Jones et al., (2001) and was followed up by Fujisawa and Eguchi, (2008), Broniatowski et al., (2012) and Fujisawa, (2013) among others. In spite of its formal similarity with the DPD, this family, referred to herein as the logarithmic density power divergence (LDPD) family, was originally developed following the robust model fitting idea of Windham, (1995). This class also has a non-negative robustness tuning parameter $\alpha$ and is defined as

[TABLE]

The divergence $\rho_{0}^{(0)}$ is once again defined as a limiting case and some simple algebra shows that $\rho_{0}^{(0)}(g,f_{\theta})=\rho_{1}^{(0)}(g,f_{\theta})$, the common divergence being a version of the Kullback-Leibler divergence.

Jones et al., (2001) also provided a general form of divergence measures in terms of two tuning parameters $\alpha$ and $\phi$ given by

[TABLE]

for $0\leq\phi\leq 1$ and $\alpha\geq 0$. The DPD and the LDPD measures are recovered from this general form for $\phi=1$ and $\phi=0$ , the second one being obtained as the limiting case as $\phi\rightarrow 0$. Accordingly, the DPD and the LDPD families have also been referred to as type 1 and type 0 divergences in the literature. A scaled version of the LDPD has also been called the $\gamma$-divergence by Fujisawa and Eguchi, (2008). Both the DPD and the LDPD measures have proved to be useful additions to the literature of robust parameter estimation based on divergences. Both families, as well as the corresponding minimum divergence estimators, have been heavily cited in the literature and applied to many practical problems, and a thorough exploration of their relationship and hierarchy can help us develop better compromises. In either case, a simple density power down-weighting indicates the source of robustness of the resulting estimator. In both cases, the parameter $\alpha$ controls the trade-off between robustness and efficiency; smaller values of $\alpha$ lead to greater model efficiency and larger values of $\alpha$ lead to greater outlier stability. The minimum divergence estimator corresponding to $\alpha=0$ is the maximum likelihood estimator in either case.

In terms of comparison between the minimum DPD and the minimum LDPD estimators, Jones et al., (2001) expressed a (weak) preference for the former. In particular, they observed that the presence of observations very close to zero could lead to spurious global minimum in case of the LDPD under the exponential model. On the other hand, Fujisawa and Eguchi, (2008) and Fujisawa, (2013) have claimed that the minimum LDPD estimators exhibit a greater relative stability under heavy contamination leading to smaller bias and smaller mean squared error. They argued that the small bias that the minimum LDPD estimator has is due to an approximate Pythagorean relation that the LDPD satisfies. In contrast, Broniatowski et al., (2012), based on their own simulation study, do not report any particular relative advantage for the minimum LDPD estimator.

These points raise certain unsettled issues involving these two useful classes of divergences and corresponding estimators which we hope to at least partially reconcile in the present paper. The generalized family of divergences in Equation (1.4), useful as it is, does not generate minimum divergence estimators that are legitimate $M$-estimators except when $\phi=0$ or $1$. We will construct an alternative generalized class of divergences providing a bridge between these classes where several of the intermediate divergences lead to reasonable compromises between the positives of these two families. This new family will be called the family of bridge density power divergences and the generated estimators are all $M$-estimators. This will provide the user with the flexibility of choosing a suitable estimator from a larger class.

The new family of divergences will depend on two tuning parameters which are (i) the robustness parameter $\alpha$ and (ii) the bridge parameter $\lambda$. It would be of interest to choose the tuning parameters adaptively with respect to the proportion of outliers so that the estimation is optimal in an appropriate sense. The robustness tuning parameter $\alpha$ should be close to zero for pure data and should assume a moderately large positive value in case of contaminated data. In this respect, Hong and Kim, (2001) proposed the first method of choosing the tuning parameter by minimizing an estimate of the asymptotic variance and Warwick and Jones, (2005) refined this process by using the mean squared error criterion together with a pilot estimate. We provide some justification for using a modified Hong and Kim, (2001) procedure; see Section 8.

Before concluding this section, we summarize the main achievements of this paper.

We introduce a new family of divergences, the Bridge Density Power Divergences (BDPD) which produce a smooth link between the DPD and the LDPD; each intermediate divergence is decomposable and leads to an $M$-estimator. Apart from helping to understand the relations between the DPD and the LDPD, this family is important in its own right, and in specific cases the performance of the intermediate divergences turn out to be superior than both the two marginal divergences (DPD and LDPD). 2. 2.

We demonstrate that the spurious root problem observed in case of the LDPD as noticed by Jones et al., (2001) is not an isolated problem for the exponential model, and is a more general phenomenon. 3. 3.

It is shown that the results and assertions of Fujisawa and Eguchi, (2008) and Fujisawa, (2013) are substantially correct, but only when a number of residual concerns are answered. In particular, the observed superior performance of the “minimum LDPD estimator”, is, most of the time, achieved at a local minimum of the LDPD measure (rather than a global one), and has to be viewed as a weighted method of moments estimator rather than a minimum divergence estimator. 4. 4.

We demonstrate that the weighted method of moment equations corresponding to the LDPD (and indeed many other BDPD measures) can potentially throw up multiple roots; this is a real problem with practical implications. As the desired root is not necessarily the minimizer of the corresponding divergence, there is no automatic root selection strategy. As a choice of the incorrect solution can be disastrous, it is imperative that a clear answer to this question is provided in the literature, where this issue is so far unexplored. 5. 5.

We clearly define the target parameter when using the LDPD for parametric estimation. Our results in this paper show that it cannot, in general, be the global minimizer of the LDPD, nor any arbitrary root of the weighted method of methods equation. Our description of the target depends on an actual numerical algorithm for its selection (see item 6 below). This interpretation of the target parameter actually extends to the entire BDPD family, except the ordinary DPD. 6. 6.

Finally, we provide an algorithm for the selection of the suitable root of the weighted method of moment equations for any member of the BDPD family, including the LDPD. We also provide a method for the selection of “optimal” tuning parameters for analyzing the real data.

The rest of the paper is organized as follows. In Section 2, we introduce the basic parametric set up and also provide the construction leading to the family of bridge divergences. Various properties of the estimators including strong consistency, asymptotic normality and robustness are discussed in Section 3. In Section 4, we provide a heuristic explanation for the spurious minimum behaviour of the LDPD measures in the case of small outliers. Section 5 rigorously proves that the LDPD is bound to fail in some specific examples. In Section 6, we address issues regarding multiple roots of the LDPD estimating equation and well-definedness of the estimators as discussed in the existing literature. We also introduce an algorithm for getting hold of a good root of the LDPD and all the other BDPD estimating equations. In Section 7, we provide the results of a simulation study conducted to support the claims made in previous sections. In Section 8, we give some directions on how to choose the tuning parameters. Finally we conclude with a few remarks in Section 9.

2 The Bridge Density Power Divergence (BDPD) Family

The problem of parameter estimation considered in this paper uses the following set up and notation. Let $\mathcal{G}$ denote the set of all probability distributions having densities with respect to some $\sigma$-finite base measure $\mu$. We assume that the data generating distribution $G$ and the model family $\mathcal{F}_{\Theta}=\{F_{\theta}:\theta\in\Theta\subset\mathbb{R}^{p}\}$ belong to $\mathcal{G}$. Let $g$ and $f_{\theta}$ be the corresponding densities (with respect to $\mu$). Let $X_{1},X_{2},\ldots,X_{n}$ be an independent and identically distributed (iid) random sample from $G$ which is to be modeled by the family $\mathcal{F}_{\Theta}$. Our aim is to estimate the parameter $\theta$ by choosing the model density which gives the “closest fit” to the data. In this paper, we quantify the “closeness” using the density power divergences, the logarithmic density power divergences or their generalizations as mentioned in the previous section. However we will see in the latter sections that this idea of closeness will require a refinement for most of the minimum BDPD estimators. All the integrals considered in this paper including those in the previous section are with respect to the measure $\mu.$

In order to estimate $\theta$ based on an iid sample, one needs to construct an empirical estimate of the divergence. In this regard, note that the first terms of Equations (1.1) and (1.3), depending only on the model density $f_{\theta}$, need no estimation. The third terms are independent of $\theta$ and so do not figure in any minimization process over $\theta$. For the second terms, note that the integral involved can be written as $\mathbb{E}_{g}f_{\theta}^{\alpha}(X)$ and so can be estimated by the average of $f_{\theta}^{\alpha}(X_{i}),\,1\leq i\leq n$. This is a consequence of the decomposability property referred to in Section 1. In light of this discussion, the parameter estimates based on the divergences $\rho^{(\alpha)}_{1}$ and $\rho^{(\alpha)}_{0}$ are given by

[TABLE]

and

[TABLE]

Under certain regularity conditions allowing the interchange of the derivative and the integral, the estimating equations corresponding to the above divergences are given, respectively, by

[TABLE]

and

[TABLE]

Here $u_{\theta}(x)$ represents the likelihood score function given by $u_{\theta}(x)=\nabla\log f_{\theta}(x)$. Throughout this manuscript, we use the symbols $\nabla,\nabla_{2}$ to denote the first and the second derivatives with respect to $\theta$. Equations (2.3) and (2.4) demonstrate that both the minimum DPD and the minimum LDPD estimators are legitimate $M$-estimators. Later in this section and in Section 6, we will consider the weighted method of moments representation of these estimating equations.

In this paper, we will show that the DPD and the LDPD families can be combined in a larger super-family of divergences where each intermediate divergence class leads to an $M$-estimator unlike the generalization given in Equation (1.4). For notational simplicity in the following derivation, define

[TABLE]

We consider the density versions (in terms of the true density $g$) of Equations (2.3) and (2.4) rather than the versions based on the empiricals. These equations are given by

[TABLE]

and

[TABLE]

Rewriting the estimating equations (2.5) and (2.6) in terms of $t_{1}$ and $t_{2}$, we get

[TABLE]

For $\lambda\in[0,1]$, consider an estimating equation which equates a convex combination of the above two estimating functions to zero; it is given explicitly by

[TABLE]

Rearranging the terms on the left hand side leads to the differential equation,

[TABLE]

where $\bar{\lambda}=1-\lambda$. It is now easy to derive the corresponding objective function which turns out to be

[TABLE]

for some constant $C$ independent of $\theta$. Imposing the condition that the divergence has to be zero when $g=f_{\theta}$, the constant $C$ is recovered to be $\frac{1}{\bar{\lambda}}\frac{1}{\alpha}\log(\lambda+\bar{\lambda}t_{3}(\theta))$. Hence the new divergence can be written as

[TABLE]

This is our class of bridge density power divergences, defined for $\alpha\geq 0$ and $\lambda\in[0,1]$. For $\lambda=0$, this reduces to the class of logarithmic density power divergences; as $\lambda\to 1$, we recover the class of density power divergences. It is also easy to verify that as $\alpha\rightarrow 0$, the limiting divergence is the Kullback-Leibler divergence as given in Equation (1.2) irrespective of the value of $\lambda$. Due to the consideration of efficiency, we will, in the rest of the paper, restrict the robustness parameter $\alpha$ to the range $[0,1]$.

The estimating equation (2.10) for the bridge divergence with parameters $(\alpha,\lambda)$ can be written as:

[TABLE]

This may be referred to as the BDPD moment equations. It is immediately obvious that the above equation represents a very rich class of score moment equations. When $\alpha=0$, the above represents the ordinary likelihood score equation irrespective of the value of $\lambda$. More generally it represents a score moment equation where the scores are variably weighted by powers of densities depending on the parameter $\alpha$, and the weights are variably normalized depending on the parameter $\lambda$. When $\lambda=0$, the weights are perfectly normalized (as in the case of the LDPD), while we get the non-normalized equations for $\lambda=1$, as in the case of the DPD. For intermediate values of $\lambda$ we get partially normalized estimating equations, with the degree of normalization dropping off with increasing $\lambda$. We will observe the effect of this normalization throughout the rest of the paper, most notably in Section 4, where we will give some indication of its role in the occurrence of the spurious roots. Remark 2.1 The BDPD family was constructed using an appropriate combination of the DPD and the LDPD estimating equations with the same tuning parameter $\alpha$; there is no specific reason, apart from mathematical convenience, for considering the same $\alpha$ in both the estimating equations. However, if estimating equations with different values of $\alpha$ are combined it does not appear to lead to a genuine objective function for such an estimating equation. $\diamond$

Remark 2.2 Equation (2.9) cannot be obtained by starting directly with a convex combination of the DPD and the LDPD with the same value of $\alpha$. $\diamond$

3 Properties of Minimum BDPD Estimators

Based on the bridge density power divergence defined in the previous section and an iid random sample $X_{1},X_{2},\ldots,X_{n}$ from a distribution $G$, an estimator of $\theta$ is given by

[TABLE]

By definition the BDPD leads to an $M$-estimator and so its asymptotic and robustness properties can be derived from the well-established $M$-estimation theory. Therefore, we present the results for consistency and asymptotic normality of the bridge divergence estimator without proofs. For this, we use the following notation.

[TABLE]

We now present a consistency theorem with conditions in alignment with those in Wald’s consistency theorem (Ferguson, (1996, Chapter 17)) for the maximum likelihood estimator and as the proof is essentially the same, it is moved to the Supplementary Material.

Theorem 3.1** (Consistency).**

Suppose that the following assumptions hold.

1. (C1)

$\Theta$* is a compact metric space;* 2. (C2)

There exists a function $K(x)$ (independent of $\theta$) such that $|\varphi_{\alpha}(x)|\leq K(x)$ and $K(X)$ has finite expectation (with respect to $G$). Here $\varphi_{\alpha}(x)=f_{\theta}^{\alpha}(x)$ for $\alpha>0$ and $\varphi_{0}(x)=\log f_{\theta}(x)$; 3. (C3)

For each $x$ and any sequence $\theta_{n}\to\theta$,

[TABLE]

for all $x$, except possibly on a set (which might depend on $\theta$ but not on the sequence $\{\theta_{n}\}$) of $\mu$-measure zero. Also, $\theta_{g}^{(\alpha,\lambda)}$ is the unique minimizer of the bridge density power divergence $\rho^{(\alpha,\lambda)}(g,f_{\theta})$.

Then the minimum BDPD estimator $\hat{\theta}_{n}^{(\alpha,\lambda)}$ is strongly consistent for $\theta_{g}^{(\alpha,\lambda)}$ for any fixed $\alpha\in[0,1]$ and $\lambda>0$.

Remark 3.1 Note that the conditions are independent of the value of $\lambda$ and match those of Wald’s consistency theorem when $\alpha=0$. The consistency theorem only requires $\Theta$ to be a metric space and not an Euclidean space. The proof can be extended to settings other than iid samples using various generalizations of the uniform strong law of large numbers (USLLN). $\diamond$

Remark 3.2 The main ingredient in the proof of Theorem 3.1 is Theorem 5.7 of van der Vaart, (1998) that uses uniform convergence of the sample based objective functions to their population counterparts. For $\lambda>0$, the function $\xi\mapsto\log(\lambda+\bar{\lambda}\xi)$ is Lipschitz on any bounded closed interval $I$ that does not contain [math]. This fact allows one to conclude the following implication (under (C2)):

[TABLE]

At $\lambda=0$ the Lipschitz property fails and the above implication breaks down. This is the reason for restricting the range of $\lambda$ to $\lambda>0$. Thus the above proof does not automatically cover the DPD family. To extend the consistency claim to $\lambda=0$, one needs to make the additional assumption $\inf_{\theta\in\Theta}\int gf_{\theta}^{\alpha}>0.$ $\diamond$

The proof of asymptotic normality of the minimum BDPD estimator can be obtained through a quadratic approximation of the objective function. Using Cram$\acute{e}$r-Rao type regularity conditions, we also observe that there exists a sequence of roots of the estimating equation which is consistent and asymptotically normal. We do not repeat the conditions here but refer the reader to Theorem 2 of Basu et al., (1998). The proof is similar to Lehmann’s proof (Lehmann and Casella, (1998)) of consistency and asymptotic normality of the MLE, and is hence omitted.

Theorem 3.2** (Asymptotic Normality).**

Under certain regularity conditions, there exists a sequence of roots $\theta_{n}$ of the bridge divergence estimating equation which is consistent and $\sqrt{n}({\theta}_{n}-\theta_{g}^{(\alpha,\lambda)})$ has an asymptotically normal distribution with mean 0 and variance given by $J(\theta_{g})^{-1}K(\theta_{g})J(\theta_{g})^{-1}$, where $\theta_{g}:=\theta_{g}^{(\alpha,\lambda)}$,

[TABLE]

3.1 Pythagorean Relation

Fujisawa and Eguchi, (2008) have claimed that under heavy contamination, the minimum LDPD estimator achieves a smaller bias (compared to the minimum DPD estimator); it is also suggested that this phenomenon can be partially explained by an approximate Pythagorean relation which the LDPD satisfies. We present a sequence of results which shows that such a Pythagorean relation holds in general for all BDPD (except the DPD), so that the LDPD result is a special case of the more general result. For any $0\leq t\leq 1$, let $\bar{t}=1-t$. Define the cross entropy between any two densities $g$ and $f$ for $\lambda,\alpha\in[0,1]$ by,

[TABLE]

The divergence induced by this cross entropy is given by

[TABLE]

which is a scaled version of $\rho^{\lambda,\alpha}(g,f)$, satisfying $D_{\lambda,\alpha}(\cdot,\cdot)=\rho^{\lambda,\alpha}(\cdot,\cdot)/(1+\alpha)$.

Theorem 3.3**.**

Let $f$ and $\delta$ be given probability density functions and let $0\leq\varepsilon\leq 1$. If $g(\cdot)=(1-\varepsilon)f(\cdot)+\varepsilon\delta(\cdot)$ and $h$ is any positive function, then for any $\alpha\in[0,1]$ and $\lambda\in[0,1)$,

[TABLE]

where

[TABLE]

For $\lambda=1$, we have

[TABLE]

Proof.

First consider the case $0\leq\lambda<1$. By definition of $d_{\lambda,\alpha}$, we have, using $\ell(h)=\frac{1}{\bar{\lambda}(1+\alpha)}\log\left(\lambda+\bar{\lambda}\int h^{1+\alpha}\right)$,

[TABLE]

Now by a simple Taylor series expansion, we get that

[TABLE]

from which the result follows. For the case $\lambda=1$, observe that

[TABLE]

∎

Remark 3.3 This theorem can be used to prove a Pythagorean relation (for $0\leq\lambda<1$) similar to Theorem 3.2 of Fujisawa and Eguchi, (2008). The statement of the result is given below for completeness; however the proof is very similar to the one in Fujisawa and Eguchi, (2008) and is hence omitted. The Fujisawa and Eguchi, (2008) result thus becomes a particular case of the following theorem. $\diamond$

Theorem 3.4** (Pythagorean Relation).**

Let $f$ and $\delta$ be given probability density functions and let $0\leq\varepsilon\leq 1$. If $g(\cdot)=(1-\varepsilon)f(\cdot)+\varepsilon\delta(\cdot)$ and $h$ is any positive function, then for any $0\leq\lambda<1$,

[TABLE]

where

[TABLE]

Remark 3.4 Theorem 3.4 is one of the main theoretical reasons for the behaviour of the minimum BDPD estimators to be observed in the subsequent sections over various choices of tuning parameters and various contaminating distributions. Note that the function

[TABLE]

for fixed positive values of $a$ and $b$, is necessarily increasing in $\lambda\in[0,1)$ if and only if $a\leq b$. This implies that when

[TABLE]

the error term $T_{\varepsilon,\delta}$ is an increasing function of $\lambda\in[0,1)$ and so the bias of the minimum BDPD estimator is expected to be an increasing function of $\lambda\in[0,1)$; however as the Pythagorean relation does not exist for the DPD, calculations based on $\xi(\lambda)$ are not helpful in theoretically comparing the bias of the minimum DPD estimator. The condition (3.2) with $h=f_{\theta}$ and $f=f_{\theta_{0}}$ for $\theta,\theta_{0}\in\Theta$ is implied by the usual approximate singularity condition in the robustness literature; the latter condition dictates that the contaminating distribution (in this case represented by $\delta$) in the gross error model is approximately singular with the parametric family so that $\int\delta h^{\alpha}\approx 0$ and $\int fh^{\alpha}$ is relatively large. In the numerical simulations, we will see that this reasoning generally fits in well with the observed behaviour of the estimators except for very small values of $\alpha$; in the latter case the quantities (in Equation (3.2)) can be very close and the above observations may not hold.

But when the contaminating distribution $\delta$ is concentrated at (or around) the mode of the model density $h=f_{\theta_{0}}$, the integral $\int\delta h^{\alpha}$ is often at least as large as $\int fh^{\alpha}$ so that the changing behaviour of $\xi(\lambda)$ over $\lambda$ with $a=\int\delta h^{\alpha}$ and $b=\int fh^{\alpha}$, is less predictable. We will observe in the rest of the paper that the performance of the minimum LDPD estimator (and several other minimum BDPD estimators) may suffer badly in this case. $\diamond$

In the following sections, we will consider two kinds of contamination for the parametric model. The first case will correspond to the approximate singularity idea of the gross error model where the contaminating (minor) component is well-separated from the target (major) component. We will refer to this as outer contamination; in this case the contaminating values will be surprising observations in the sense of Lindsay, (1994). In the second case of contamination, we will choose the contaminating component near the mode of the major component so that these observations are no longer surprising observations but will nevertheless distort the shape of the distribution relative to the parametric model. We will refer to this case as inner contamination.

4 Spurious Behaviour under Inner Contamination

Jones et al., (2001) reported that in case of the exponential distribution, small outliers (near the origin) made the minimum LDPD estimator non-robust. These small outliers are essentially what we have described as constituting inner contamination in Section 3. Let $\eta_{\theta}(x)$ denote the density of the exponential distribution with mean $\theta$. The non-robustness of the minimum LDPD estimator under small outliers was demonstrated by Jones et al., (2001) in simulation studies and was confirmed by determining the LDPD between the densities $\eta_{\theta}$ and the $0.85\eta_{1}+0.15\delta_{x_{0}}$ mixture with $x_{0}=0.0001$ where $\delta_{x_{0}}$ represents the indicator function at $x=x_{0}$. In Figures 1 and 2, we have exhibited the BDPD objective function between the above two densities over $\theta$ for $\alpha=0.5$ and several values of $\lambda$ in $[0,1]$.

It is clear that the global minimizer of the objective function remains stuck at a value very close to zero at least up to $\lambda=0.7$. It is noteworthy that this spurious minimum may easily be missed by any gradient descent root search algorithm since the minimum is obtained as a needle sharp notch over a very limited region. It might very well be possible (and indeed it is the case in Figures 1 and 2) that a local minimizer of the LDPD may serve as a reasonable robust solution even in this case. But the asymptotic results for a sequence of roots of the estimating equation present in the literature do not prescribe the method to choose a suitable sequence of roots in case of multiple roots. In any case, it is clear that the global minimizer of the LDPD up to (at least) $\lambda=0.7$ is a nonsensical value.

The main reason for this phenomenon in the LDPD and the bridge divergences close to it is the behaviour of the $\log$ function at [math], where it diverges to $-\infty$. If one writes the LDPD between the $\eta_{\theta}$ and $0.85\eta_{1}+0.15\delta_{0.0001}$ densities, we can easily check that there is a sharp drop in the objective function around a value very close to zero because of the above logarithm effect (although at 0, the objective function is positive infinity). This effect slowly smooths out as the $\lambda$ parameter increases since the presence of the additive $\lambda$ term in each of the logarithms eventually forces the argument to be bounded away from zero. The validity of this reasoning is confirmed by examining the behaviour of the LDPD and the BDPD close to it in case of the $N(0,\sigma^{2})$ model, where it can be expected that the estimate of $\sigma$ will be driven to zero if there are some outliers near [math] in the sample. For that matter, this behaviour can also be seen in case of the $N(\mu,\sigma^{2})$ model with some outliers near the true mean. Figure 3 exhibits the LDPD at $\alpha=0.5$ as a function of $\sigma$ when computed between the densities of the $N(0,\sigma^{2})$ and the mixture $0.85N(0,1)+0.15\delta_{0.001}$. Clearly, there is a reasonable local minimum around $\sigma=0.8$; but it is beaten hands down by the useless global minimum in the neighbourhood of zero.

While the behaviour observed in Figures 1–3 represent the patterns in the divergences between the actual densities, such behaviour can also be observed under pure data (although it is relatively rare in scale models). An actual random sample of size 20 from $N(0,1)$ with seed 129 was obtained in R as

[TABLE]

The exact value of the last observation is $1.048653\times 10^{-5}$, which forces a spurious behaviour in the LDPD objective function plotted in Figure 4 for $\alpha=0.8$. To be precise, if $f_{\sigma}(x)$ denotes the probability density of $N(0,\sigma^{2})$ with respect to the Lebesgue measure, then

[TABLE]

and the LDPD objective function satisfies

[TABLE]

where $X_{20}$ is the last observation which is close to zero and $n=20$. This inequality confirms our reasoning. As $\sigma$ slides down towards zero, for a while the first term involving the logarithm dominates, pulling the objective function down. However as $\sigma$ gets really close to zero, the third term takes over and there is an extremely sharp rise in the objective function, which generates a minimum with a razor sharp notch. Figure 4 shows the spurious (global) minimum near zero (at $\sigma=0.00001309906$ to be exact), although there is a reasonable local minimum at $\sigma=1.265882$. For this sample, the spurious global minimum phenomenon is observed for all $\alpha\geq 0.500144205$ and all $\lambda\in[0,0.228827308]$. In contrast, the global minimum of the DPD objective function (also plotted in Figure 4) for $\alpha=0.8$ is obtained at $\sigma=1.20833.$ If the last observation $1.048653\times 10^{-5}$ were removed from this data set, this spurious minimum behaviour of the LDPD objective function disappears. The minimum LDPD estimator now equals $\sigma=1.321807$ at $\alpha=0.8$, an entirely reasonable value. (The corresponding minimum DPD estimator is $1.298228$). Thus one single observation can bring about an absolutely drastic change in the minimum LDPD estimator which is against the spirit of stability that robust estimators should have; so far we (or indeed anybody else), have not detected such spurious behaviour in any scenario involving the DPD.

Thus the spurious root issue in case of the LDPD is not a isolated problem limited to the case of the exponential distribution. It is, in fact, a more serious problem in the case of the location-scale model (compared to just the scale model), as we will see in the next section.

5 Unboundedness of the Bridge Divergences

In the previous section, we explained the reason for observing a spurious minimum with LDPD in the case of inner contamination and argued that it might happen even in case of real data generated from the pure model. Intuitively, this may be explained by the fact that the probability of observing a value near the mode of the majority distribution is not too small at moderate sample sizes.

Here we will formally prove that the sample version of the BDPD objective function (given on the right hand side of Equation (2.12)) is unbounded below in case the parametric model family is a location-scale family $\mathcal{G}_{f,\Theta}$, where $f$ is a probability density function, $\Theta=\mathbb{R}\times\mathbb{R}^{+}$, and

[TABLE]

This result proved in Theorem 5.1 implies that one cannot fit a location-scale family using the minimizer of a BDPD, except the DPD.

Theorem 5.1**.**

Let $X_{1},X_{2},\ldots,X_{n}$ be independent and identically distributed observations from a density $g,$ which is modeled by the densities in the location-scale family of distributions $\mathcal{G}_{f,\Theta}$ with a fixed $f$ satisfying $f(0)>0.$ Then, for any $0\leq\lambda<1$ and $\alpha>0,$

[TABLE]

Proof.

Set the objective function on the left hand side of Equation (5.1) as $M_{n}(\theta)$. Then for $\theta=(\mu,\sigma),$

[TABLE]

We will now show that for each $1\leq j\leq n$, $M_{n}(X_{j},\sigma)$ (that is, $\theta=(X_{j},\sigma)$) converges to $-\infty$ as $\sigma\downarrow 0$. Fix $1\leq j\leq n$. Since $f(x)\geq 0$ for all $x$, we get, by taking $\mu=X_{j}$,

[TABLE]

Substituting this bound in $M_{n}(\theta)=M_{n}(X_{j},\sigma)$, we obtain

[TABLE]

where

[TABLE]

Letting $\sigma$ tend to zero implies $M_{n}(X_{j},\sigma)\to-\infty$ proving Equation (5.1). ∎

Remark 5.1 From the proof, it is seen that the global minimizer of the bridge divergence in case of any sample is at the extreme point ($\sigma=0$) and the global minimum over all $(\mu,\sigma)$ combinations, is attained for atleast $n$ points namely $(X_{j},0),\,1\leq j\leq n$. This is similar in spirit to the classical example of likelihood inference for Gaussian mixture modelling as explained in Example 2.4.5 of Bickel and Doksum, (2015). Note that the theorem does not cover the case of DPD which is given below. $\diamond$

In order to illustrate the phenomenon discussed above, we plot the LDPD objective with $\alpha=0.5$ in Figure 5, based on an artificial sample of size 4, where the sample observations are $-2,-1,0.5,2.$ The underlying model family is assumed to be $\{N(\mu,\sigma^{2}):\mu\in\mathbb{R},\sigma>0\}$. We observe that the LDPD objective drops down sharply to $-\infty$ as $\sigma$ approaches [math], when $\mu$ is either of the four sample points. See Section S.2 of supplementary materials for more details.

Theorem 5.2**.**

Let $X_{1},X_{2},\ldots,X_{n}$ be independent and identically distributed observations from a density $g,$ which is modeled by the densities of the location-scale family $\mathcal{G}_{f,\Theta}$ with a fixed $f$. Suppose that $\lim_{|x|\rightarrow\infty}f(x)=0.$ Then, for any $\alpha>0,$ the following holds for all $n$ if $\mu$ does not belong to the set $\{X_{1},...,X_{n}\}$, and for all $n>\frac{(1+\alpha)f^{\alpha}(0)}{\alpha\int f^{1+\alpha}}$ if $\mu$ is equal to some $X_{i}$ in the above set:

[TABLE]

Proof.

First, note that the quantity on the left hand side of (5.2) whose limit needs to be evaluated as $\sigma\rightarrow 0$, can be written as

[TABLE]

Observe that, if $\mu$ does not belong to the set $\{X_{1},...,X_{n}\}$, then

[TABLE]

Consequently, (5.3) goes to $\infty$ as $\sigma\rightarrow 0.$ On the other hand, if $\mu$ equals $X_{j}$ for some $1\leq j\leq n,$ then (5.3) can be written as:

[TABLE]

It is now easy to see that for $n>\frac{(1+\alpha)f^{\alpha}(0)}{\alpha\int f^{1+\alpha}}$, (5.4) goes to $\infty$ as $\sigma\rightarrow 0.$ This completes the proof. ∎

Remark 5.2 Theorem 5.2 shows that the DPD objective function for the location-scale family (for $\alpha>0$) diverges to $\infty$ as $\sigma\rightarrow 0,$ for all $n,$ if $\mu$ is different from all the observations. On the other hand, if $\mu$ equals one of the observations, the situation differs in the sense that there exists $a_{f}$ such that if $n>a_{f}$, the DPD objective function diverges to $\infty$ and if $n<a_{f}$, it diverges to $-\infty$. The case of $n=a_{f}$ (when $a_{f}$ is a natural number) depends on the rate of convergence of $f(x)$ to zero as $|x|\to\infty$. $\diamond$

Remark 5.3 Theorem 5.1, in particular, covers the case of $N(\mu,\sigma^{2})$ parametric family and out of entire BDPD family the DPD alone stands out as a practically valid method of inference. Fujisawa and Eguchi, (2008) provided an iterative algorithm to obtain the minimizer of the LDPD in case the parametric family is one of the exponential families, applicable to $N(\mu,\sigma^{2})$ family. However, while their iterative algorithm is monotone decreasing, it does not in general guarantee convergence to the global minimizer. In fact, the small bias property shown for the scale parameter in their numerical study (Table 2, Fujisawa and Eguchi, (2008)) suggests, in view of our discussion in Section 4, that their algorithm mostly converges to a local minimizer, atleast in the normal family. This is also very similar to the case of likelihood inference for Gaussian mixture modelling wherein one uses the EM algorithm to get hold of the useful local minimizer of the log-likelihood. This fact is also consistent with the asymptotic normality result presented in Section 5 of Fujisawa and Eguchi, (2008). The asymptotic normality holds for a sequence of roots of the estimating equation while the global minimizer does not exist (in the interior of the parameter space). Note that in this case the global minimizer does not appear as a root of the estimating equation owing to the fact that any global minimum is attained on the boundary of the parameter space. $\diamond$

Remark 5.4 The result as stated in Theorem 5.1 does not apply to a scale family with a fixed location (say, [math]). But if one of the observations is exactly [math], then the proof as presented above can be employed to get the same conclusion. More precisely, we have, (taking $\sigma=X_{(1)}$),

[TABLE]

where $D_{f}=\log\int f^{1+\alpha}-\left(1+{\alpha}\right)\log f(1)$ under the assumption $f(1)>0$. Here $X_{(1)}$ denotes the first order statistic in the sample $X_{1},\ldots,X_{n}$ and $f_{\sigma}(x)=f(x/\sigma)/\sigma$ for a density $f$ on $\mathbb{R}^{+}$. This indicates that if we have inner contamination (that is, contamination near the mode [math]) so that $n^{1+1/\alpha}X_{(1)}$ is small enough, then there is a possibility of a spurious global minimum near zero. Observe that the true distribution is a mixture described as

[TABLE]

where $Z\sim\mbox{Bernoulli}(p)$ and $\delta_{0}$ represents a point mass at [math]. This represents a simple case inner contamination model with contamination level $p$. Then in a sample of size $n$, we have for any fixed $\varepsilon>0,$

[TABLE]

implying convergence with probability one of $n^{1+1/\alpha}X_{(1)}$ to zero as $n\to\infty$ for any fixed $p>0$. Here $F$ represents the cumulative distribution function of $X$. In a given sample though, this spurious minimum may disappear (or may not be significantly pronounced) depending on the proportion of contamination. In this case, the BDPD do not behave so badly as in the location-scale family. As already observed, the divergences close to the DPD will lead to reasonable estimators. $\diamond$

6 Towards the Desired Estimator

As we have seen, the analysis using the LDPD can be beset with different kinds of problems, both theoretical and computational. In this section we add to the discussion of the possible flaws in the current state of the analysis based on the LDPD, eventually giving a recipe for consolidating the existing knowledge and overcoming the present difficulties so that one can arrive at a best compromise. We do not exactly refute the findings of Fujisawa and Eguchi, (2008) and Fujisawa, (2013) as we find a substantial part of this research to be valuable and useful. However, there are too many loose ends that remain unaccounted for which limit the practical applications of the method without further consolidation.

Fujisawa and Eguchi, (2008) and Fujisawa, (2013) define the latent bias of an estimator as the difference between the target parameter and the limit of the estimator. To describe the basic flaw with the minimized LDPD estimator or in general the normalized estimating equations, we rework the arguments of Fujisawa, (2013) which gives the hint of a very small latent bias under heavy contamination for the minimum LDPD estimator. Let $f(x)=f_{\theta^{*}}(x)$ be the target density within our parametric family and $\delta(x)$ be the contamination density so that the data generating density $g$ is given by

[TABLE]

Let $\ell(x;\theta)=\log f_{\theta}(x)$ and $u_{\theta}(x)=\nabla\ell(x;\theta)$ be the usual log-likelihood and the score, respectively. A general normalized estimating equation is given by

[TABLE]

where $\xi(\cdot)$ is a non-negative weight function satisfying $\xi(a)\to 0$ as $a\to-\infty$ and $\mathbb{E}_{h}[\cdot]$ represents the expectation with respect to the density $h$. Assume that

[TABLE]

in a neighborhood of $\theta=\theta^{*}$, which is usually satisfied under the classical outlier model formulation (or, as we call it, outer contamination). Substituting the form of the density $g$, we get

[TABLE]

Whenever the approximations in (6.2) hold, we get the following approximate equality

[TABLE]

Thus, $\theta^{*}$ is an approximate solution of the normalized estimating equation (6.1). For concreteness, consider the case where $f(x)=f_{0}(x)$ (i.e., $\theta^{*}=0$), where $f_{\mu}$ denotes the density of the normal distribution with mean $\mu$ and variance $1$ and $\delta(x)=f_{10}(x)$. Take $\theta=10$ in the normalized estimating equation (6.1). It is easy to see that

[TABLE]

in a neighborhood of $\theta=10$. Therefore, as above, the equality

[TABLE]

holds approximately. Thus, not only $\theta^{*}=0$ but $\theta=10$ is also an approximate root of this estimating equation. The latent bias argument for normalizing estimating equations would be one sided if one were to use it only for $\theta^{*}$, and not for the root corresponding to the contaminating component. One probably can have an estimator with small latent bias in this case, but one would have to appropriately choose between the multiple roots to ensure this small bias.

To validate our argument, we generated 20 observations randomly from the mixture distribution $0.9N(10,1)+0.1N(0,1)$. We present the entire sample for the future reproducibility of our results.

[TABLE]

The parametric family is taken to be $\mathcal{F}=\{N(\mu,1):\mu\in\mathbb{R}\}.$ The LDPD objective function in (2.2) with $\alpha=0.5$ is plotted against $\theta$ in Figure 6. It is seen to have a (local) minimum near [math] and a (possibly global) minimum around $10$, which in this case is the mean of the larger component of the mixture. Note that while there are at least two obvious roots of the corresponding estimating equation, a clear root selection strategy is unavailable. While it so happens that in this case the divergence is well behaved, in many standard cases it is not (as we have seen in Sections 4 and 5).

We are now in a position to give an analysis of the minimum LDPD estimator in particular (and the minimum BDPD estimator in general) by consolidating the different bits of results that we have presented so far in this paper. In Section 3 we have proved the consistency of the minimum BDPD estimator. What it means in layman’s terms is that the global minimizer of the empirical LDPD objective function converges to the global minimizer of the theoretical objective function obtained by replacing the empirical distribution function by the true one. However, the usefulness of these consistency results are immediately challenged by our findings in Sections 4 and 5 as well as the previous part of this section. First of all, it is quite possible, even in case of pure data, that the global minimizer of the empirical objective function is a value which is entirely useless to the statistician for all practical purposes. Yet, the procedure might actually have a very reasonable local minimizer close to our target. This is what we have observed in Figure 4, for example, Fujisawa and Eguchi, (2008) or Fujisawa, (2013) – or, indeed, anybody else for that matter – do not give any root selection criteria in such cases; and in the absence of any such criteria, one would be hard put to justify the choice of a local minimizer, when a perfectly legitimate (although statistically useless) global minimizer exists. Our results in Section 5 show that the global minimizers for the LDPD would always be at the boundary of the parameter space and will be of very little practical value to us under location-scale models. However, in almost every such case there exists – as we have seen in our simulations – reasonable roots of the estimating equations, representing appropriate local minima of the objective function. The bias and mean square error reported in Figures 2-5 of Fujisawa, (2013) correspond to the root generated by this reasonable local (but not global) minimum.

Secondly, an equally serious problem is the extreme shift in the target parameter itself under inner contamination when LDPD is the divergence of choice. As we have seen in Figure 3, a small proportion of inner contamination is enough to shift the target parameter (the global minimizer of the divergence between the theoretical densities) to such an extreme degree that the new target does not characterize the original major component in any way at all. Yet, once again, the divergence has a local minimizer which reasonably characterizes the the major component. Once again, however, there is no conceivable reason for considering a local minimizer as the target in the presence of a global minimizer (which, unfortunately, does not characterize the component of interest).

It is clear, therefore, that usual inference based on the LDPD alone can lead the inference astray, if one proceeds with the global minimizer. Other reasonable roots do exist, but root selection strategies do not. Existing literature claims the existence of “a root” which has the desired properties, without any recipe for arriving at that root. Such inadequacy exists in case of many other members of the BDPD class. The only member within this class which are entirely free – as far as it is observed in our explorations – from the anomalies of spurious roots and ill-defined targets is the DPD. In every model and every case that we have looked at, the DPD has a well defined minimizer that would reasonably represent, in theory, a properly described target, and would generate a suitable estimator for the same target under randomly generated data. Acknowledging that the LDPD does have certain roots which are superior in dealing with the bias under heavy contamination, we propose to utilize the global minimizer of the DPD to arrive at the desired root of the LDPD. The next subsection describes the chain algorithm which we propose for this purpose.

6.1 The Chain Algorithm

Under repeated simulations we have observed, under many standard parametric models, that as far as the global minimizers of the BDPDs are considered over all $\lambda\in[0,1]$ given a specific value $\alpha$, the minimum DPD estimator is almost always the best in terms of latent bias (and not the minimum LDPD estimator). At the same time we acknowledge that the minimum LDPD estimating equations (or, more generally, the minimum BDPD estimating equations) may have roots which represent some local minimum of the divergences which might improve upon the performance of the minimum DPD estimator in terms of latent bias. As indicated by the heuristics, and as we have seen in our simulations, the LDPD (more generally the BDPD) has a local minimizer that is close to the true $\theta^{*}$. This we believe is at least partially due to the Pythagorean relation that we presented in Section 3. Since the members of the BDPD family form a literal bridge between the DPD and LDPD, we have an opportunity for getting at a root of the LDPD estimating equations that is consistent and asymptotically normal. A completely rigorous proof is unavailable at this time, but we provide a strong heuristic argument and extensive simulations to support this claim.

From the form of the bridge divergence, for a fixed $\alpha$ and $\theta$, we obtain the limit

[TABLE]

Heuristically this indicates that $\rho^{(\alpha,\lambda)}(g,f_{\theta})$ has at least a local minimizer close to $\hat{\theta}^{\alpha}_{n1}$ (defined in Equation (2.1)) if $\lambda$ is close to one. For example, one can take $\lambda=1-n^{-1/2}$ which is sufficiently close to $1$ for this root to be reasonably good. We shall now present a chain algorithm for getting hold of a good root of the LDPD estimating equation (indeed, all bridge estimating equations). The chain algorithm proceeds in the following steps: Fix $\alpha\in[0,1].$

First choose a sequence $\{\lambda_{i}\}$ satisfying

[TABLE]

with $r_{n}\to 0$ at some rate. Define a function $\lambda\mapsto\hat{\theta}^{\alpha}(\lambda).$ 2. 2.

Solve the DPD problem completely, that is, find the global minimizer $\hat{\theta}^{\alpha}_{n1}$ of the sample estimate of the DPD. Set $\hat{\theta}(\lambda_{K})=\hat{\theta}^{\alpha}_{n1}$. Note that $\lambda=1$ in the bridge divergence expression in Equation (2.12) corresponds to the DPD (in the limit). 3. 3.

For $i=K-1,K-2,\ldots,0$, find the local minimizer of the sample estimate of the bridge divergence with parameters $(\alpha,\lambda_{i})$ that is closest to $\hat{\theta}(\lambda_{i+1}).$ Set this closest local minimizer as $\hat{\theta}(\lambda_{i})$. 4. 4.

Return $\{\hat{\theta}(\lambda_{i}):\,0\leq i\leq K\}$.

Step 3 above can be solved by using any off-the-shelf algorithm for minimization with $\hat{\theta}(\lambda_{i+1})$ as a starting point. To draw analogue with existing algorithms of this type, we note that the LASSO that has taken over the literature in the past few years is solved by using a chain/path algorithm as above; there the issue is about getting fast algorithm not distinguishing local and global minimizers (it is a convex problem).

Conjecture 6.1**.**

In any model under the assumptions that guarantee the asymptotic normality of $\hat{\theta}_{n1}^{\alpha}$, $\hat{\theta}^{\alpha}(\lambda_{0})$ is the consistent root of the LDPD estimating equation that is closest to $\theta^{*}$. And this is the solution claimed by Fujisawa and Eguchi, (2008) to have a small latent bias.

The simulations using this chain algorithm are very encouraging and we think that a proof of the above conjecture can be obtained by an application of the Taylor series and using (6.3) with an appropriate choice of rate of convergence of $\lambda_{k}-\lambda_{0}$ to zero. In many nicely behaved parametric densities, simply choosing the root of LDPD estimating equation closest to the minimum DPD might be good enough. However, in general cases, the difference between LDPD and DPD estimating equations can be drastic for this simple trick to work (at least for a theoretical guarantee).

In the following section we will frequently use the term “minimum bridge divergence estimator”. It is to be understood that this will refer to the local minimizer obtained by using the chain algorithm which uses the global minimizer of the DPD as the starting value. In addition, we will also frequently use the term“root of the estimating equation”. In this case it will be understood that the correspondence is with a local minimizer, and not an intervening local maximizer (or saddle point), which may also produce a legitimate root of the estimating equation.

7 Simulation Study

In this section, we apply the chain algorithm discussed in Section 6 on simulated data. We consider the exponential and the normal scale families, the former to illustrate the case of outer contamination, and the latter to illustrate the case of inner contamination. The estimators are computed from 1000 replications. The bias and mean squared error over these $1000$ replications are calculated against the target component of the underlying majority distribution.

The chain algorithm is applied for each pair $(\alpha,\epsilon)$, where the robustness parameter $\alpha$ runs from [math] to $1$ in steps of $0.2$ and the contamination level $\epsilon$ is [math], $0.05$ or $0.2$. For each $(\alpha,\epsilon)$ combination, the bridge parameter $\lambda$ runs from $1$ to [math] in steps of $0.1$ as the chain algorithm proceeds. Since the parameter $\alpha$ is somewhat better understood in the literature and our introduced parameter $\lambda$ needs to be analyzed more deeply, the latter is varied through finer steps than the former.

7.1 Exponential Scale Model

Here the model is the class of exponential distributions with mean $\sigma$ over $\sigma\in(0,\infty)$. Data are simulated from the mixture $(1-\epsilon)$Exp$(1)$ + $\epsilon$Unif$(6-10^{-4},6+10^{-4})$, where the first component is exponential with rate 1, and the second is uniform over the indicated range. The contamination level $\epsilon$ is taken to be $0,0.05$ and $0.2$. For step 2 of the chain algorithm which involves solving the DPD problem completely, $25$ starting points are drawn uniformly from the interval $[0,0.1]$ and the remaining $75$ uniformly from the interval $(0.1,10]$.

7.2 Normal Scale Model

Here the model is the class of normal distributions with mean $5$ and variance $\sigma^{2}$ over $\sigma\in(0,\infty)$. Data are simulated from the mixture $(1-\epsilon)N(5,1)$ + $\epsilon$Unif$(5-10^{-5},5+10^{-5})$, where the first and the second components are normal and uniform respectively, with parameters as indicated. The contamination level $\epsilon$ is taken to be $0,0.05$ and $0.2$. For step 2 of the chain algorithm which involves solving the DPD problem completely, $25$ starting points are drawn uniformly from the interval $[0,0.1]$ and the remaining $75$ uniformly from the interval $(0.1,10]$.

For the sake of brevity, we will only include the graphs of the mean-squared errors (scaled by sample size $n=100$) of the minimum bridge divergence estimators in the exponential scale model case, for contamination levels $0.05$ and $0.2$ (in Figure 7, respectively) in the main article. Tables presenting exact values of the bias and MSE of the minimum bridge divergence estimators for both the exponential and normal scale families, for all the three contamination levels $0,0.05$ and $0.2$ are given in the online supplement.

It is observed that for the $5\%$ contamination level, the MSE of the minimum bridge divergence estimators decreases as the chain algorithm proceeds ($\lambda$ goes from 1 to 0) for $\alpha$ upto $0.6$, and increases as the chain algorithm proceeds for $\alpha=0.8$ and $1$. On the other hand, for the $20\%$ contamination level, the MSE of the minimum bridge divergence estimators decreases as the chain algorithm proceeds from $\lambda=1$ to $\lambda=0$ for all values of $\alpha$. As a function of $\alpha$ (with $\lambda$ held fixed), for $5\%$ contamination, the MSE decreases roughly upto $\alpha=0.6$ and then increases with $\alpha$, whereas for $20\%$ contamination, a completely decreasing trend of the MSE is observed as $\alpha$ increases.

Note that the values reported here are not, except for the DPD ($\lambda=1$ case) case, based on the global minimizers of the divergences; rather they are the chain algorithm solutions. What this shows is that for heavy contamination (at the level of 20%) the LDPD solution (not necessarily the global minimizer) as obtained from the chain algorithm does seem to dominate the other BDPD solutions, and the mean square errors decrease uniformly in the direction $1\rightarrow 0$ over $\lambda$ for each $\alpha$. This provides partial confirmation of the results of Fujisawa and Eguchi, (2008) and Fujisawa, (2013). However the results are now stronger and more useful in that we identify the root for which it works; it does not work for just any root, and certainly not for the global minimizer. This is so in spite of the fact that the contamination considered here is an outer contamination. Any exponential distribution has a mode at zero, and randomly throwing up values close to zero, leading to spurious global minimizers, is not entirely uncommon.

A similar situation is observed in the normal distribution case. Here the contamination is an inner contamination, which leads to spurious global minimizers more often for the BDPDs. However the chain algorithm guides the process to a sensible root (that is also a local minimizer) in each case.

For the data of size 20 presented in Section 4 from $N(0,1)$, the minimum bridge divergence estimators of $\sigma$ under the $N(0,\sigma^{2})$ model with $\alpha=0.8$ and varying $\lambda$ from [math] to $1$ are presented in Table 1. For each value of $\lambda$, the true global minimizers of the bridge divergence objective functions are presented in the parentheses. It is remarkable that even with $\lambda$ as high as 0.9, the spurious root phenomenon is observed.

8 Choice of Tuning Parameters

Most of the robust estimators derive their robustness from down-weighting the observations suspected as outliers. As can be seen from the BDPD objective functions, the observations are proportionally weighted by powers of the model density which demonstrates that the observations are down-weighted for being outliers with respect to the density $f_{\theta}$. It is this outlier stability property which makes the corresponding estimators desirable from the point of view of robustness. If there are no outliers in the data, then the tuning parameter $\alpha$ must ideally be set to zero in order to get optimal inference. But if the data do contain outliers, then the tuning parameters $\alpha$ and $\lambda$ should be chosen so as to partially or fully eliminate the effect of the outlying observations. However the choice of the tuning parameter must be an automatic procedure where a data-driven choice of the parameters $\alpha$ and $\lambda$ is generated by the method. One of the first methods of choosing tuning parameters in case of the DPD was proposed in Hong and Kim, (2001), where the tuning parameter is chosen by minimizing the trace of an estimate of the asymptotic covariance matrix. Warwick and Jones, (2005) refined this method to find an “optimal” tuning parameter by minimizing the trace of an estimate of the MSE of the estimator. Here the MSE is computed by separately estimating the bias and the variance components. From the asymptotic variance formula, one can easily estimate the variance. But bias estimation involves the use of a pilot estimate, and the estimated mean square error has a strong dependence on it. Warwick and Jones, (2005) used the minimum $L_{2}$ distance estimator as the pilot; however, under the current state of the art, there is no universally acceptable choice of the pilot, on which the process critically depends.

We acknowledge that the method of Warwick and Jones, (2005) would be perfect if we could eliminate its dependence on the pilot estimator, or find an independent estimate of the bias as a function of the tuning parameter. However, here we present a modified version of the approach of Hong and Kim, (2001) and give an (almost theoretical) justification of the same. Firstly, one should not use the trace of an estimate of the asymptotic covariance matrix (although some of the present authors have done so in the past in the absence of a better strategy) for the construction of the objective function to be minimized. This leads to adding quantities with different units. Consider the $N(\mu,\sigma^{2})$ example in which case the asymptotic variance of $(\hat{\mu},\hat{\sigma}^{2})$ or any of the variance estimates have the units of $X_{i}$ and $X_{i}^{2}$ that should not be added. One could, of course, think of estimating $(\mu,\sigma)$ where adding the variance estimates would be more sensible. But this remedy does not work in general parameter spaces.

Why is the asymptotic variance adequate? An intuitive explanation is as follows. It is well-known that the delete-$d$ jack-knife estimator is a consistent estimator of the asymptotic variance. The delete-$d$ jack-knife estimator is obtained by calculating the difference between the estimator calculated based on $n$ observations and the estimator calculated based on $n-d$ observations. By the definition of a good robust estimator, this difference should be small for a robust estimator but for a non-robust estimator like the maximum likelihood estimator this difference would be large if the $d$ observations deleted contain some outliers. This implies that robust estimators should have “smaller” asymptotic variance than non-robust estimators. Of course, if there is no contamination then it is known that the maximum likelihood estimator would have “smaller” variance for large enough sample size. Using a bootstrap asymptotic variance estimator also gives a similar conclusion. We do ignore the bias, but the question clearly goes in favor of a robust estimator. By definition they are closer to the true parameter based on the majority of the data.

A more concrete explanation is offered by using the closeness measure introduced in Zhang and He, (2016). Let $\mathcal{V}_{\Lambda}$ define a set of variance matrices $V_{\alpha}$ indexed by a (possibly vector) parameter $\alpha\in\Lambda$.

Definition 8.1**.**

A non-negative real-valued matrix function $\mathfrak{m}(\cdot)$ defined on the set $\mathcal{V}_{\Lambda}$ is called a closeness measure if and only if the following conditions hold.

Consistency. For any $V_{1},V_{2}\in\mathcal{V}_{\Lambda}$, if $V_{1}\succeq V_{2}$, then $\mathfrak{m}(V_{1})\geq\mathfrak{m}(V_{2})$, where the equality holds only if $V_{1}=V_{2}$. 2. 2.

Continuity. For any matrix sequence $\{V_{n}\}\subseteq\mathcal{V}_{\Lambda}$ with $V_{n}\succeq B$. As $n\to\infty$, $\lVert V_{n}-B\rVert_{\infty}\to 0$ if and only if $\mathfrak{m}(V_{n})\to\mathfrak{m}(B)$.

Here by $A\succeq B$, we mean $A-B$ is non-negative definite. Definition 8.1 implies that if there is an efficient variance matrix in the set $\mathcal{V}_{\Lambda}$, then minimizing $\mathfrak{m}(V_{\alpha})$ over $\alpha\in\Lambda$ leads to such a matrix. Zhang and He, (2016) additionally prove that the trace and the Frobenious norm are valid closeness measures. We use the determinant of the matrix as the measure of closeness (see Section S.5 of the supplementary material). Note that determinant is a valid quantity to consider from the point of view of units and it also has a practical interpretation as the volume of a confidence set. These arguments give a justification for minimizing a closeness measure of an estimate of the asymptotic variance.

Summing up all these arguments, we claim that minimizing the determinant of an estimate of the asymptotic variance provides an asymptotically optimal estimator whenever such an estimator exists in the family of estimators under consideration. A detailed analysis of this procedure would be taken up in a future paper. Just as a remark, for the normal data of size 20 presented in Section 4, the optimal $\lambda$ parameter with $\alpha=0.8$ using this procedure turned out to be $\lambda=1$. A plot of the asymptotic variance over all $\lambda$ is given in the supplementary material.

9 Concluding Remarks

In this paper, the competing families of divergences, DPD and LDPD, are critically examined for their strengths and deficiencies. The bridge divergence family introduced in this paper is an attempt at combining the good of both and nullify the disadvantages of either. Unlike the DPD, the LDPD estimating equation admits roots with small latent bias. However, these roots may not be the global minimizers of the LDPD objective. The phenomenon of spurious global minimizers of the LDPD is rigorously proved in specific parametric families where the DPD provably works. The bridge divergences partly suppresses this problem for certain tuning parameter ($\lambda$) values along the bridge.

However, the Bridge divergence is not a perfect solution in that it also faces the same problem as LDPD in some cases. The point made here is that one cannot expect the global minimizer of the LDPD to generate a good estimator and the DPD estimator is a safe bet in all the cases, but whenever possible some members of the minimum bridge divergence estimators can help reduce the latent bias of the DPD estimator.

Acknowledgements

The authors would like to thank Srijata Samanta of University of Florida for her contribution towards Remark 5.