On Selection Criteria for the Tuning Parameter in Robust Divergence

Although robust divergence, such as density power divergence and γ-divergence, is helpful for robust statistical inference in the presence of outliers, the tuning parameter that controls the degree of robustness is chosen in a rule-of-thumb, which may lead to an inefficient inference. We here propose a selection criterion based on an asymptotic approximation of the Hyvarinen score applied to an unnormalized model defined by robust divergence. The proposed selection criterion only requires first and second-order partial derivatives of an assumed density function with respect to observations, which can be easily computed regardless of the number of parameters. We demonstrate the usefulness of the proposed method via numerical studies using normal distributions and regularized linear regression.


Introduction
Data with outliers naturally arise in diverse areas. In the analysis of data containing outliers, statistical models with robust divergence are known to be powerful and have been used regularly. In particular, the density power divergence [1] and γ-divergence [2] have been routinely used in this context due to their robustness properties while there now exist others. In these studies, the theoretical properties of the robustness of robust divergence against outliers are also clarified through the analysis of influence functions. For its interesting applications, see for example [3,4] and references therein. Robust divergence, in general, holds a tuning parameter that controls robustness under model misspecification or contamination. Ref. [1] noted that there is a trade-off between estimation efficiency and strength of robustness; thereby, a suitable choice of the tuning parameter seems crucial in practice. However, a well-known selection strategy such as cross-validation is not straightforward under contamination, so that we need to rely on a trial-and-error way to find a reasonable value of the tuning parameter.
To select a turning parameter, we here propose a simple but novel selection criterion for the tuning parameter by using the asymptotic approximation of Hyvarinen score [5,6] with unnormalized models based on robust divergence. Typical existing methods [7,8] choose a tuning parameter based on the asymptotic approximation of the mean square error but have the drawback of requiring some pilot estimators and an analytical expression of the asymptotic variance. In addition, their works are essentially limited to the simple normal distribution and simple linear regression. Our proposed method has the following advantages over the existing studies.

1.
Our method does not require an explicit representation of the asymptotic variance. Therefore, our method can be applied to rather complex statistical models, such as multivariate models, which seems difficult to be handled by the previous methods; 2.
In the existing studies, it is necessary to determine a certain value as a pilot estimate to optimize a tuning parameter. Thus, the estimates may strongly depend on the pilot estimate. On the other hand, our method does not require a pilot estimate and is stable and statistically efficient; 3.
Although our proposed method is based on a simple asymptotic expansion, it is more statistically meaningful and easier to interpret the results statistically than existing methods because it is based on the theory of parameter estimation for unnormalized statistical models.
Through numerical studies under simple settings, we show that the existing methods can be sensitive to a pilot estimate and tends to select an unnecessarily larger value of a tuning parameter, leading to loss of efficiency compared with the proposed method. Moreover, we still apply the proposed selection method, an estimation procedure in which the asymptotic variance is difficult to compute. As an illustrative example of such a case, we consider robust linear regression with γ-divergence and 1 -regularization, where the existing approach is infeasible to apply.
As related works, there are two information criteria using the Hyvarinen score. [9] proposed AIC-type information criteria for unnormalized models by deriving an asymptotic unbiased estimator of the Hyvarinen score, but it does not allow unnormalized models whose normalizing constants do not exist. Hence, the criterion cannot be applied to the current situation. On the other hand, [10] proposed an information criterion via Laplace approximation of the marginal likelihood in which the potential function is constructed by the Hyvarinen score. Although [10] covers unnormalized models with possibly diverging normalizing constants, the estimator used in the criterion is entirely different from one defined as the maximizer of robust divergence; thereby, the criterion does not apply to the tuning parameter selection of robust divergence either. Moreover, ref. [11] developed an robust sequential Monte Carlo sampler based on robust divergence in which γ is adaptively selected. However, it does not provide selection of γ in a frequentist framework.
The rest of the paper is organized as follows. Section 2 introduces a new selection criterion based on the Hyvarinen score. We then provide concrete expressions of the proposed criterion under density power divergence and γ-divergence in Section 3. We numerically illustrate the proposed method in two situations in Section 4. Concluding remarks are given in Section 5.

Tuning Parameter Selection of Robust Divergence
Suppose we observe y 1 , . . . , y n as realizations from a true distribution or data generating process G, and we want to fit a statistical model { f θ : θ ∈ Θ} where Θ ⊆ R d for some d ≥ 1. Furthermore, assume that the density of G is expressed as (1 − ω) f θ * + ωδ, where δ is a contaminated distribution that produces outliers in observations. Our goal is to make statistical inference on θ * by successfully eliminating information of outliers. To this end, robust divergence such as density power divergence [1] and γ-divergence [2] is typically used for robust inference on θ * . Let y = (y 1 , . . . , y n ) be a vector of observations and D γ (y; θ) be a (negative) robust divergence with a tuning parameter γ. We assume that the robust divergence has a additive form, namely, D γ (y; θ) = ∑ n i=1 D γ (y i ; θ). This constraint is necessary when using the H-score, but it is not a strong constraint in the context of this study because the well-known robust divergence, as presented in Section 3, satisfies this property.
For selecting the tuning parameter γ, our main idea is to regard L γ (y i ; θ) ≡ exp{D γ (y i ; θ)} as an unnormalized statistical model whose normalizing constant may not exist. Recently, ref. [10] pointed out that the role of such unnormalized models can be recognized in terms of relative probability. For such model, we employ the Hyvarinen score (H-score) in terms of Bayesian model selection [5,6], defined as where L (m) γ (y) is the marginal likelihood given by with some prior distribution π(θ). We consider an asymptotic approximation of the H-score (1) under large sample sizes. Under regularity conditions [12], the Laplace approximation of (2) is where θ γ is the M-estimator given by and H( θ γ ) is the Hessian matrix at θ = θ γ . Then, we have the following approximation, where the proof is deferred to Appendix A.

Proposition 1. Under some regularity conditions, it follows that
The above results give the following approximation of the original H-score: which satisfies H n (γ) = H * n (γ) + o p (1) under n → ∞. We then define the optimal γ as γ opt = argmin γ H n (γ).
Let θ * γ be the quantity that θ γ converges to. Then, from the argument given in [5,10], the criterion (4) would converge to the Fisher divergence between the unnormalized model exp{D γ (y; θ * γ )} and the true data generating process. Although the unnormalized model is not a valid statistical model in the sense that it does not have a finite normalizing constant, the Fisher divergence can be recognized as the distance in terms of relative probabilities [10].
Existing selection strategies for γ mostly use the asymptotic variance of θ γ . For example, under the density power divergence, refs. [7,8] suggested using asymptotic approximation of the mean squared errors of θ γ . However, computation of the asymptotic variance is not straightforward, especially when an additional penalty function is incorporated into the objective function or the dimension of θ is large. On the other hand, the proposed criterion (4) does not require the computation of asymptotic variance but only needs the derivatives of robust divergence concerning y i . Furthermore, it should be noted that the proposed criterion (4) can be applied to a variety of robust divergence.

Possible Robust Divergences to Consider
We here provide detailed expressions for the proposed criterion (4) under some robust divergences. For simplicity, we focus on two robust divergences which can be empirically estimated from the data. Still, the proposed method could be applied to other divergences such as Hellinger divergence [13] or αβ-divergence [14]. In what follows, we shall use the notations, f (y i ; θ) = ∂ f (y i ; θ)/∂y i and f (y i ; θ) = ∂ 2 f (y i ; θ)/∂y 2 i .

Density Power Divergence
The density power divergence [1] for a statistical model f (y i ; θ) is It can be seen that D γ (y i ; θ) + 1 − 1/γ → log f (y i ; θ) as γ → 0, so the above function can be regarded as an extension of the standard log-likelihood. Then, a straightforward calculation leads to the expression of (4), given by

γ-Divergence
The original form of γ-divergence [2] for a statistical model f (y i ; θ) is given by , which is not an additive form. However, the maximization of the above function with respect to θ is equivalent to the maximization of the transformed version of γ-divergence, . Then, we have

Normal Distribution with Density Power Divergence
We first consider a simple example of robust estimation of the normal population mean under unknown variance. Let y 1 , . . . , y n be sampled observations and we fit N(µ, σ 2 ) to the data. The density power divergence [1] of the model is given by where φ(y i ; µ, σ 2 ) is the density function of N(µ, σ 2 ). In this case, the criterion (4) is expressed as where µ γ and σ γ are the estimator based on the density power divergence. We first demonstrate the proposed selection strategy through simulation studies. We simulated y 1 , . . . , y n from the normal distribution with true parameters, µ = 2, and σ = 1, and then replace the first nω observations by y i + 7. We adopted four settings for ω ∈ {0, 0.05, 0.1, 0.15}. Using the simulated dataset, the optimal γ is selected among {0, 0.01, . . . , 0.69, 0.70} through the criterion H n (γ), and we obtain the adaptive estimator µ γ opt . For comparison, we also employed two selection methods, OWJ [7] and IWJ [8], in which the optimal value of γ is selected via asymptotic approximation of mean squared errors of the estimator. We set γ = 0.5 to compute a pilot estimator that must be specified in the two methods. Furthermore, we also computed µ γ with γ = 0.1, 0.3, and 0.5. Using an estimator of the asymptotic variance of µ γ given in [8], we also computed the Wald-type 95% confidence interval of µ. Based on 5000 simulated datasets, we obtained the squared root of mean squared error (RMSE) of the point estimator, as well as coverage probability (CP) and average length (AL) of the interval estimation. The results are reported in Table 1. It is observed that the use of small γ (such as γ = 0.1) may lead to unsatisfactory results when the contamination is heavy. It can also be seen that with the use of relatively large γ, the estimation results can be inefficient. On the other hand, the proposed method can adaptively select a suitable value of γ since the averaged value of γ opt increases with the contamination ratio ω, as confirmed in Table 2. This would be the main reason that the proposed method provides reasonable performance in all the scenarios. We next apply the proposed method to Simon Newcomb's measurements of the speed of light data, motivated by applications in Basu et al. [1], Basak et al. [8], Stigler [15]. We searched the optimal γ among {0, 0.01, . . . , 0.69, 0.70} and the H-sores are shown in left panel in Figure 1. The obtained optimal value is γ opt = 0.09, which is substantially smaller thanγ = 0.23 selected by the existing methods as reported in [8]. Since the method proposed in [8] requires a pilot estimate and the estimation results depend significantly on it, we believe that our estimation results are more reasonable. In fact, it is unlikely that we will have to use a value of γ = 0.23 for a dataset that contains only two outliers. As shown in the right panel in Figure 1, the estimated density functions are almost the same when γ = 0.09 and when γ = 0.23. However, it would be preferable to adopt the smaller value of γ = 0.09 if the estimates are almost identical in terms of statistical efficiency.

Gamma Distribution with Density Power Divergence
We next consider robust estimation of the gamma distribution. Let y 1 , . . . , y n be sampled observations and we fit Ga(α, β) to the data, where α is a shape parameter and β is a rate parameter, so that the expectation of the gamma distribution is α/β. The density power divergence of the model is given by where f Ga (y i ; α, β) is the density function of Ga(α, β). In this case, the criterion of γ is one given in Section 3.1 in which the first and second order derivatives of the density are given as We demonstrate the proposed selection criterion through simulation studies. We simulated y 1 , . . . , y n from the gamma distribution with true parameters, α = 2, and β = 4, and then replace the first nω observations by y i + 5. We adopted three settings for ω ∈ {0, 0.05, 0.1} and two scenarios for n ∈ {100, 200}. Using the simulated dataset, the optimal γ is selected among {0, 0.01, . . . , 0.49, 0.50} through the HS criterion H n (γ), and we obtain the adaptive estimators, α γ opt and β γ opt . For comparison, we applied the standard maximum likelihood (ML) estimator, as well as the robust estimator with fixed γ ∈ {0.1, 0.5}. In this study, we do not include OWJ or IWJ methods since the asymptotic variance formula in this case has not been investigated before and the derivation would require tedious algebraic calculation.
Based on 5000 simulated datasets, we obtained the squared root of mean squared error (RMSE) of the point estimator, where the results are shown in Table 3. We also provide the average values of the selected γ in Table 4. It is observed that the (non-robust) ML and the robust method using the small γ (γ = 0.1) leads to unsatisfactory results when the data are contaminated. It can also be confirmed that γ = 0.5 does not hold reasonable accuracy when the contamination ratio is small or 0, which indicates that a suitable selection step is substantially related to the estimation result. The proposed method, however, has some adaptive property. When there is not contamination, the estimation performance is almost identical to the the ML estimator which is the most efficient in this scenario since a small value (sometimes exactly zero) of γ is selected in this scenario, as shown in Table 4. On the other hand, the estimation performance is still better than the other methods when the data are contaminated, by successfully finding a suitable value (increasing according to ω) of γ.

Regularized Linear Regression with γ-Divergence
Note that the proposed criterion can be used when some regularized terms are introduced in the objective function, while the existing method requiring an asymptotic variance of the estimator is not simply applicable. We demonstrate the advantage of the proposed method through regularized linear regression with γ-divergence [16]. Let y i and x i be a response variable and a p-dimensional vector of covariates, respectively, for i = 1, . . . , n. The model is y i ∼ N(x i β, σ 2 ). Then, the transformed γ-divergence is , and the H-score is expressed as Here, β γ and σ 2 γ are estimated as the minimizer of the following regularized γ-divergence: where λ is an additional tuning parameter that can be optimized via 10-fold cross-validation. We use the R package gamreg [16] to estimate β and σ 2 under given γ.
We apply the aforementioned method to the well-known Boston housing dataset [17]. In this analysis, we included the original 13 covariates and 12 quadratic terms of the covariates except for one binary covariate, resulting in 25 covariates in total. We searched the optimal γ among {0, 0.02, . . . , 0.68, 0.70}, and the estimated H-scores are shown in the left panel in Figure 2, where the optimal value of γ was 0.16. For comparison, we estimated the regression coefficients with γ = 0 and γ = 0.5. Note that γ = 0 reduces to the (non-robust) standard regularized linear regression. The scatter plots of the estimated standardized coefficients under γ = 0.16 against ones under the two choices of γ are shown in the right panel of Figure 2. It is confirmed that the estimates with γ = 0.16 and γ = 0.5 are comparable while there are substantial differences between estimates with γ = 0.16 and γ = 0, indicating that a certain amount of robustness is required for the dataset.

Concluding Remarks
We proposed a new criterion for selecting the optimal tuning parameter in robust divergence, using the Hyvarinen score for unnormalized models with robust divergence. The proposed criterion does not require the asymptotic variance formula of the estimator that is needed in the existing selection methods. Although we simply focused on the univariate and continuous situation, the proposed criterion can also be applied to multivariate or discrete distribution, where finite differences under discrete distributions should replace derivatives. Applications of the proposed score to such cases would also be helpful, and we left it to future work. where we defined S γ (y i ; θ) = ∂S γ (y i ; θ)/∂y i . Note that ∂ θ γ /∂y i = O p (n −1 ) under large n. From (3), the first order partial derivative of the marginal log-likelihood can be approximated as Under the regularity conditions for π(θ), it follows that under large n. From the same argument, we can also show that ∂ log |H( θ γ )|/∂y i = o p (1). Regarding the first term in (A1), we have since S γ (y j ; θ) is a score function and ∑ n j=1 S γ (y j ; θ γ ) = O p (n 1/2 ). Using the expression of the first order derivative (A2), it holds that S γ (y j ; θ γ ) H( θ γ ) −1 ∂ ∂y i S γ (y i ; θ γ ).