Ya N Pa

Fellow, IEEE Derek Abbott

Noise Enhancement in Robust Estimation of Location

Keywords: Noise enhancement, location, M-estimator, asymptotic efficiency, maximum bias

In this paper, we investigate the noise benefits to maximum likelihood type estimators (M-estimator) for the robust estimation of a location parameter. Two distinct noise benefits are shown to be accessible under these conditions. With symmetric heavy-tailed noise distributions, the asymptotic efficiency of the estimation can be enhanced by injecting extra noise into the Mestimators. With an asymmetric contaminated noise model having a convex cumulative distribution function, we demonstrate that addition of noise can reduce the maximum bias of the median estimator. These findings extend the analysis of stochastic resonance effects for noise-enhanced signal and information processing.

I. INTRODUCTION

A N OPTIMAL noise level, obtained by appropriately adding extra noise to a given signal processor or by tuning the existing noise level, can sometimes improve information processing [START_REF] Mitaim | Adaptive stochastic resonance[END_REF]- [START_REF] Chen | Noise-enhanced information systems[END_REF]. This effect was initially called stochastic resonance [START_REF] Benzi | The mechanism of stochastic resonance[END_REF] but is now widely referred to as noise enhancement [START_REF] Zozor | Noise-aided processing: Revisting dithering in a Σ∆ quantizer[END_REF]- [START_REF] Rousseau | Noise-improved Bayesian estimation with arrays of one-bit quantizers[END_REF], resulting in noise benefit-a term coined by Kosko [START_REF] Patel | Optimal noise benefits in Neyman-Pearson and inequality-constrained statistical signal detection[END_REF]. In recent years, noise-enhanced signal processing theory has received significant attention [START_REF] Kay | Can detectability be improved by adding noise?[END_REF]- [START_REF] Chen | Noise-enhanced information systems[END_REF], [START_REF] Zozor | Noise-aided processing: Revisting dithering in a Σ∆ quantizer[END_REF]- [START_REF] Harmer | A review of stochastic resonance: Circuits and measurement[END_REF]. In the signal estimation field, classical parameter estimation problems [START_REF] Chen | Noise-enhanced information systems[END_REF], [START_REF] Chen | Noise enhanced parameter estimation[END_REF], [START_REF] Balkan | CRLB based optimal noise enhanced parameter estimation using quantized observations[END_REF], [START_REF] Papadopoulos | Sequential signal encoding from noisy measurements using quantizers with dynamic bias control[END_REF]- [START_REF] Dabeer | Signal parameter estimation using 1-bit dithered quantization[END_REF] and the Bayesian method [START_REF] Chen | Noise enhanced parameter estimation[END_REF]- [START_REF] Rousseau | Noise-improved Bayesian estimation with arrays of one-bit quantizers[END_REF], [START_REF] Uhlich | Bayes risk reduction of estimators using artificial observation noise[END_REF]- [START_REF] Zeitler | Bayesian parameter estimation using single-bit dithered quantization[END_REF] are of considerable current research interest. The related results of [START_REF] Chen | Noise-enhanced information systems[END_REF], [START_REF] Chen | Noise enhanced parameter estimation[END_REF]- [START_REF] Rousseau | Noise-improved Bayesian estimation with arrays of one-bit quantizers[END_REF], [START_REF] Uhlich | Bayes risk reduction of estimators using artificial observation noise[END_REF]- [START_REF] Zhu | Parameter estimation from quantized observations in multiplicative noise environments[END_REF] confirm that noise enhancement can be exploited to optimize the estimator performance.

Statistical signal processing often relies on strong assumptions, e.g., a Gaussian distribution or a particular signal model, under which optimal estimators or detectors are derived [START_REF] Huber | Robust Statistics[END_REF]- [START_REF] Zoubir | Robust estimation in signal processing: A tutorial-style treatment of fundamental concepts[END_REF]. Nevertheless, many areas of engineering, for instance, Y. Pan and F. Duan are with the Institute of Complexity Science, Qingdao University, Qingdao 266071, China (e-mail: yanpn.pan@gmail.com; fabing.duan@gmail.com).

F. Chapeau-Blondeau is with the Laboratoire Angevin de Recherche en Ingénierie des Systèmes, Université d'Angers, Angers 49000, France (e-mail: chapeau@univ-angers.fr).

D. Abbott is with the Centre for Biomedical Engineering and School of Electrical and Electronic Engineering, University of Adelaide, Adelaide, SA 5005, Australia (e-mail: derek.abbott@adelaide.edu.au). mobile communication channels, radar systems and neural networks, result in outliers or atypical observations that do not obey these assumptions [START_REF] Zoubir | Robust estimation in signal processing: A tutorial-style treatment of fundamental concepts[END_REF]. This fact can lead to significantly degraded estimator performance, and also highlights the need for robust estimation techniques that are insensitive to outliers [START_REF] Huber | Robust Statistics[END_REF]- [START_REF] Zoubir | Robust estimation in signal processing: A tutorial-style treatment of fundamental concepts[END_REF].

In this paper, we mainly focus on the noise-enhanced effect in maximum likelihood type estimators (M-estimator) [START_REF] Huber | Robust Statistics[END_REF]- [START_REF] Zoubir | Robust estimation in signal processing: A tutorial-style treatment of fundamental concepts[END_REF] for estimating a location parameter. Under this scenario, we derive two key results that demonstrate two distinct noise benefits:

(i) For the symmetric heavy-tailed noise distributions, we derive the expression of the asymptotic efficiency for an array of identical M-estimators, and show the possibility of the improved asymptotic efficiency of M-estimators by injecting mutually independent noise components. For a single M-estimator, a discriminant function, elicited from the background noise probability density and the given M-estimator function, provides the condition of existence of the optimal additive noise. It is proven that the optimal additive noise, when it exists, is the symmetrical dichotomous noise given by Chen et al. [START_REF] Chen | Theory of the stochastic resonance effect in signal detection: Part I-Fixed detectors[END_REF]- [START_REF] Chen | Noise enhanced parameter estimation[END_REF]. When the optimal noise does not exist in a single M-estimator, it is observed that the asymptotic efficiency can still be enhanced to be very close to the upper bound of unity by an array of M-estimators with added noise selected from a parametric class of noise. This fact motivates us to focus on the asymptotic efficiency of an array of identical M-estimators with more than one estimator element. We further demonstrate that, for a given M-estimator function and fixed noise levels, the asymptotic efficiency is a monotonically increasing function of the M-estimator number. Thus, the asymptotic efficiency reaches its maximum value for an array with an infinite number of identical M-estimators. Aiming to maximizing the asymptotic efficiency of an infinite number of M-estimators, the optimal probability density of the added noise is proven to be the weighted minimum L 2 -norm solution of Eq. [START_REF] Perlman | Jensen's inequality for a convex vector-valued function on an infinite-dimension[END_REF]. Specially, by the Cauchy-Schwarz inequality, the upper bound of the asymptotic efficiency of an infinite number of M-estimators is proven to be unity, and the corresponding optimal noise density is the deconvolution of the maximum likelihood estimator and the given M-estimator function. However, this optimal noise density is frequently unattainable, due to absence of a solution to the deconvolution and the fact that the infinite-size array of M-estimators can only be approached in practice. Therefore, as an alternative to such a theoretical optimal distribution, we can select an aprioriparametric class of noise distributions, of interest, and optimize the parameters and variance of the added noise. The results show that, compared to the method without added noise, the addition of extra noise that obeys the optimized parametric noise distribution does provide an improved asymptotic efficiency. It is also noted that the observed noise-enhanced effect for the robust estimation of a location parameter can also be related to suprathreshold stochastic resonance (SSR) that has been demonstrated in a parallel array of identical elements with the injection of independent noise components [START_REF] Patel | Optimal mean-square noise benefits in quantizerarray linear estimation[END_REF]- [START_REF] Chapeau-Blondeau | Noise-enhanced fisher information in parallel arrays of sensors with saturation[END_REF], [START_REF] Stocks | Suprathreshold stochastic resonance in multilevel threshold systems[END_REF], [START_REF] Mcdonnell | Stochastic Resonance: From Suprathreshold Stochastic Resonance To Stochastic Signal Quantization[END_REF]. In the context of SSR, such nonlinear elements were chosen as threshold comparators [START_REF] Rousseau | Suprathreshold stochastic resonance and noise-enhanced fisher information in arrays of threshold devices[END_REF], quantizers [START_REF] Patel | Optimal mean-square noise benefits in quantizerarray linear estimation[END_REF]- [START_REF] Chen | Nonparametric one-bit quantizers for distributed estimation[END_REF], [START_REF] Xu | Adaptive recursive algorithm for optimal weighted suprathreshold stochastic resonance[END_REF], [START_REF] Mcdonnell | Stochastic Resonance: From Suprathreshold Stochastic Resonance To Stochastic Signal Quantization[END_REF] or saturating devices [START_REF] Chapeau-Blondeau | Noise-enhanced fisher information in parallel arrays of sensors with saturation[END_REF] and have been evaluated by various performance measures such as mean square error (MSE) [START_REF] Patel | Optimal mean-square noise benefits in quantizerarray linear estimation[END_REF], [START_REF] Chen | Nonparametric one-bit quantizers for distributed estimation[END_REF], [START_REF] Xu | Adaptive recursive algorithm for optimal weighted suprathreshold stochastic resonance[END_REF], [START_REF] Mcdonnell | Stochastic Resonance: From Suprathreshold Stochastic Resonance To Stochastic Signal Quantization[END_REF], Fisher information [START_REF] Chapeau-Blondeau | Fisher information and noise-aided power estimation from one-bit quantizers[END_REF], [START_REF] Rousseau | Suprathreshold stochastic resonance and noise-enhanced fisher information in arrays of threshold devices[END_REF], [START_REF] Chapeau-Blondeau | Noise-enhanced fisher information in parallel arrays of sensors with saturation[END_REF], and Cramer-Rao bound [START_REF] Chen | Performance limit for distributed estimation systems with identical one-bit quantizers[END_REF], [START_REF] Chen | Nonparametric one-bit quantizers for distributed estimation[END_REF], [START_REF] Mcdonnell | Stochastic Resonance: From Suprathreshold Stochastic Resonance To Stochastic Signal Quantization[END_REF] in the signal estimation field. By contrast here, we mainly focus on the noise-enhanced effect in a parallel array of robust M-estimators evaluated by the measure of the asymptotic efficiency.

(ii) In an asymmetric contaminated noise model, the bias dominates the variance component of M-estimators for a sufficiently large number of observations. We theoretically prove that, under the condition of a convex cumulative distribution function, the maximum bias of M-estimators can be reduced by adding extra noise. Moreover, aiming to minimize the maximum bias of the median estimator, the dichotomous noise is proved to be optimal. Theoretical and numerical results also show that the maximum bias can be distinctly diminished at an optimal dichotomous noise level.

II. M-ESTIMATOR FOR A LOCATION PARAMETER

Consider a location model of observations [START_REF] Huber | Robust Statistics[END_REF], [START_REF] Maronna | Robust Statistics: Theory and Methods[END_REF]

x i = θ + w i ,i=1, 2, ••• ,n (1)
where the location parameter θ is unknown and observation errors (noise) w i are independent and identically distributed (i.i.d.) random variables with a common probability density f w . Then, the observations x i have the distribution f w (x -θ).

Here, we assume w = σw 0 is with the scale parameter σ that represents the noise intensity, and w 0 has its probability density f w 0 and unity scale.

Assuming that the loss function ρ satisfies certain regularities [START_REF] Huber | Robust Statistics[END_REF], [START_REF] Maronna | Robust Statistics: Theory and Methods[END_REF], we define the M-estimator θ as

θ = arg min θ n i=1 ρ(x i -θ). (2)
Differentiating Eq. ( 2) with respect to θ yields

n i=1 ψ(x i -θ)=0, (3)
where ψ = dρ/dθ = ρ ′ and the root θ that corresponds to the global minimum of the loss function in Eq. ( 2) is just the M-estimator or the generalized maximum likelihood estimator [START_REF] Huber | Robust Statistics[END_REF]- [START_REF] Zoubir | Robust estimation in signal processing: A tutorial-style treatment of fundamental concepts[END_REF]. The maximum likelihood estimator is a special case of M-estimators when ψ = -f ′ w /f w ψ M and ρ = -log f w [START_REF] Huber | Robust Statistics[END_REF]- [START_REF] Zoubir | Robust estimation in signal processing: A tutorial-style treatment of fundamental concepts[END_REF].

In order to compare the approximation by θ with the true value of θ, we usually compute the MSE

MSE( θ)=E[( θ -θ) 2 ]=v ar( θ)+b 2 ( θ), (4)
with the estimator variance var( θ)=E( θ2 ) -E 2 ( θ) and the bias b( θ)=E( θ) -θ [34]- [START_REF] Zoubir | Robust estimation in signal processing: A tutorial-style treatment of fundamental concepts[END_REF]. Here, the expectation operator E(•)= •f w (x -θ)dx represents the expectation according to the shifted probability density f w (x -θ) that is the density for the noisy random signal x i in Eq. ( 1). According to the central limit theorem and for a sufficiently large observation number n,a Fisher-consistent M-estimator θ satisfies E[ψ(x -θ)] = 0 [34]- [START_REF] Zoubir | Robust estimation in signal processing: A tutorial-style treatment of fundamental concepts[END_REF], and asymptotically converges to the Gaussian distribution with mean θ and variance

var( θ) ≈ 1 n E[ψ 2 (x -θ)] E 2 [ψ ′ (x -θ)] = 1 n E w [ψ 2 (x)] E 2 w [ψ ′ (x)] , (5)
by taking the first-order Taylor expansion of ψ(x i -θ) at θ [START_REF] Huber | Robust Statistics[END_REF]- [START_REF] Zoubir | Robust estimation in signal processing: A tutorial-style treatment of fundamental concepts[END_REF] (also see Appendix B). Here, the expectation operator

E w (•)= •f w (x)
dx is according to the probability density of the random noise w i . Also in Eq. ( 5) the change of variable x -θ → x has been performed. It is noted that the minimum asymptotic variance of var( θ)=1/ nJ(θ) is achieved by the maximum likelihood estimator [START_REF] Maronna | Robust Statistics: Theory and Methods[END_REF], [START_REF] Cover | Elements of Information Theory[END_REF], and the Fisher information J(θ) with respect to the location parameter θ is defined as

J(θ)=E (df w (x -θ)/dθ) 2 f 2 w (x -θ) =E w (df w (x)/dx) 2 f 2 w (x) = J(f w ), (6)
where J(f w ) is the Fisher information of the probability density f w [START_REF] Cover | Elements of Information Theory[END_REF]. Therefore, in order to measure how near the Mestimator θ is to the optimal estimator, the asymptotic efficiency is defined as the ratio [START_REF] Maronna | Robust Statistics: Theory and Methods[END_REF] Eff

( θ)= 1 nJ(f w ) 1 var( θ) = 1 J(f w ) E 2 w [ψ ′ (x)] E w [ψ 2 (x)] . (7)
In the following section, we will use the asymptotic efficiency of Eq. ( 7) and the bias b( θ) to assess the occurrence of noise-enhanced effects in M-estimators associated with two noise distribution models.

III. NOISE-ENHANCED ASYMPTOTIC EFFICIENCY

In a number of areas of engineering, outliers in the measurement cause a heavy-tailed distribution and so a Gaussian assumption no longer holds true [START_REF] Huber | Robust Statistics[END_REF]- [START_REF] Zoubir | Robust estimation in signal processing: A tutorial-style treatment of fundamental concepts[END_REF]. In this Section, we consider the symmetric heavy-tailed distribution model of noise. Under such scenarios, the function ψ in Eq. ( 3) is assumed to be odd, so that there is no bias problem for the unbiased Mestimator θ derived from Eq. (3). Thus, we mainly discuss the noise-enhanced asymptotic efficiency effect in this situation under consideration.

When the expression for the noise distribution f w is known, the maximum likelihood estimator is the optimal option if it exists [START_REF] Huber | Robust Statistics[END_REF]- [START_REF] Zoubir | Robust estimation in signal processing: A tutorial-style treatment of fundamental concepts[END_REF]. However, in practice, the noise distribution f w may be incompletely known, or a complete closed-form of the maximum likelihood estimator is too complex to implement [START_REF] Huber | Robust Statistics[END_REF]- [START_REF] Zoubir | Robust estimation in signal processing: A tutorial-style treatment of fundamental concepts[END_REF]. Thus, other M-estimators with an easily implemented feature are frequently employed [START_REF] Huber | Robust Statistics[END_REF]- [START_REF] Zoubir | Robust estimation in signal processing: A tutorial-style treatment of fundamental concepts[END_REF]. We note that Mestimators are critically connected the ψ-function in Eq. ( 3), which can be grouped into the monotone or the redescending classes [START_REF] Maronna | Robust Statistics: Theory and Methods[END_REF]. Here, we mainly consider the monotone Huber function [START_REF] Huber | Robust Statistics[END_REF]

ψ(x)= x, |x|≤γ, γ sgn(x), |x| >γ, (8)
and the redescending bisquare function [START_REF] Huber | Robust Statistics[END_REF], [START_REF] Maronna | Robust Statistics: Theory and Methods[END_REF]

ψ(x)=x 1 -(x/γ) 2 2 , (9)
for |x|≤γ and otherwise zero. Here, γ is an adjustable parameter of the M-estimator. We further consider a scale family of the probability density f w (x)=f w 0 (x/σ)/σ, where f w 0 has unity scale and the Fisher information of Eq. ( 6) can be expressed as J(f w )=J(f w 0 )/σ 2 [START_REF] Cover | Elements of Information Theory[END_REF], [START_REF] Kassam | Signal Detection in Non-Gaussian Noise[END_REF]. Here, J(f w 0 ) is the Fisher information of the probability density f w 0 . Then, the asymptotic efficiency of Eq. ( 7) can be calculated as

Eff( θ)= 1 -2F w 0 (-γ σ ) 2 2J(f w 0 ) γ/σ 0 x 2 f w 0 (x)dx + γ 2 σ 2 F w 0 (-γ σ ) (10)
for the Huber estimator of Eq. ( 8) and

Eff( θ)= γ/σ -γ/σ 1 -σ 2 γ 2 x 2 1 -5σ 2 γ 2 x 2 f w 0 (x)dx 2 J(f w 0 ) γ/σ -γ/σ x 2 1 -σ 2 γ 2 x 2 4 f w 0 (x)dx (11)
for the bisquare estimator of Eq. ( 9). Here, F w 0 is the cumulative distribution function of noise w 0 .

For instance, we especially consider the generalized Cauchy noise with the probability density

f w (x)=B(k, v) 1+ 1 v |x| A(k) k -(v +1/k) (12)
with

A(k)=σ Γ(1/k)/Γ(3/k), B(k, v)=kv -1/k Γ(v + 1/k)/[2A(k)Γ(v)Γ(1/k)],
and parameters σ, k, v > 0 [START_REF] Kassam | Robust techniques for signal processing: A survey[END_REF], [START_REF] Kassam | Signal Detection in Non-Gaussian Noise[END_REF]. This generalized Cauchy distribution of Eq. ( 12) is capable of modelling a wide range of practical noise types. For instance, we obtain a typical heavy-tailed distribution model for v =1/2 and k =2, which is the Cauchy noise [START_REF] Kassam | Signal Detection in Non-Gaussian Noise[END_REF]. In Fig. 1(a) and (b), we plot the asymptotic efficiencies of Eq. [START_REF] Chen | Theory of the stochastic resonance effect in signal detection-Part II: Variable detectors[END_REF] and Eq. ( 11) as a function of γ/σ, respectively. It is clearly seen that, for the generalized Cauchy noise in Eq. ( 12) with v =0.5, 1, 2, 5 and k =2, the asymptotic efficiency Eff( θ) can be enhanced to its maximum as γ/σ increases. For instance, in the case of Cauchy noise with v =1/2 and k =2in Eq. ( 12), the asymptotic efficiency Eff( θ) of Huber estimator of Eq. ( 10) is maximized to the value of 0.88 at γ/σ =0.39 and for the bisquare estimator, the asymptotic efficiency Eff( θ) of Eq. [START_REF] Chen | Noise enhanced parameter estimation[END_REF] has the maximum value of 0.90 as γ/σ =3.28. For a given parameter γ, this is the noise-enhanced effect occurring in the location M-estimator, since the asymptotic efficiency is maximized at an optimal noise scale. Of course, if the noise scale σ is given, we can also optimize the estimator parameter γ to obtain the maximum asymptotic efficiency. 10) for the Huber estimator and (b) Eq. ( 11) for the bisquare estimator as a function of γ/σ for the generalized Cauchy noise in Eq. ( 12) and the Laplacian noise in Eq. ( 18) (β =1).

Here, the generalized Cauchy noise is with parameters v =1/2, 1, 2, 5, ∞ and k =2.

It is also noted that, for the Gaussian noise (v = ∞ and k =2in Eq. ( 12)), the asymptotic efficiency Eff( θ) approaches unity as γ/σ increases, as shown in Fig. 1(a) and (b). The reason is that, as the parameter γ →∞, both Huber and the bisquare functions are linear, which just corresponds to the optimal estimator for the Gaussian noise. For the Laplacian noise in Eq. ( 18) (β =1), Fig. 1(a) shows the asymptotic efficiency Eff( θ) monotonically decreases from unity as γ/σ increases. The reason is that, as γ → 0, the Huber function of Eq. ( 8) becomes ψ(x)=sgn(x) that corresponds to the maximum likelihood estimator for the Laplacian noise, i.e., the median estimator [START_REF] Huber | Robust Statistics[END_REF]- [START_REF] Zoubir | Robust estimation in signal processing: A tutorial-style treatment of fundamental concepts[END_REF].

Next, two interesting problems arise. The first problem is, when the noise scale σ can be exactly estimated or is known and the asymptotic efficiency is maximized at an optimal parameter γ opt , whether the addition of extra noise is beneficial or not?

The second problem arises because the noise scale σ is often unknown, while the optimization of the asymptotic efficiency Eff( θ) is closely tied with the ratio γ/σ, as shown in Fig. 1. Thus, the M-estimators of Eqs. ( 8) and ( 9) often estimate the location away of the maximum asymptotic efficiency obtained by the optimal parameter γ opt . Under this circumstance, can we add extra noise to the given observations and obtain an improved performance of a robust M-estimator with a fixed parameter γ? In the following discussion, the answers to these two interesting problems will be presented. 46) where addition of noise is beneficial, and the complementary region S -of the negative signs of D(ψ, f w ), versus the parameter v of the generalized Cauchy distribution f w in Eq. ( 12) (k =2) and the parameter γ of (a) the Huber and (b) the bisquare function. Here, the existing noise scale σ =1and the added uniform noise level d → 0.

We purposefully add the extra noise to the original observations

y i = x i + η i = θ + z i ,i=1, 2, ••• ,n (13)
where the composite noise components z i = w i + η i and the extra noise components η i are the i.i.d. random variables with the common probability density f η . Therefore, z i have a convolved probability density f z (z)= f w (z -η)f η (η)dη.I ti s proven by the Fisher information convolution inequality in Appendix A that the addition of extra noise to the original observations cannot improve the performance of the maximum likelihood estimator. Furthermore, inspired by the suprathreshold stochastic resonance model [START_REF] Stocks | Suprathreshold stochastic resonance in multilevel threshold systems[END_REF], [START_REF] Mcdonnell | Stochastic Resonance: From Suprathreshold Stochastic Resonance To Stochastic Signal Quantization[END_REF], we collect L statistics θl generated by identical M-estimators for l =1, 2, ••• ,L, and average them as the M-estimation of the location parameter in the context of noise benefits. It is emphasized that each M-estimator is subject to the same input x i but mutually independent noise components η li with the common probability density f η . Therefore, the input of each M-estimator is

y li = θ + w i + η li = θ + z li , (14)
where the composite noise z li = w i + η li have the probability density f z . The estimate θl satisfies n i=1 ψ(y li -θl )=0of Eq. ( 3), and the devised M-estimator is given by

θ = 1 L L l=1 θl . ( 15
)
Due to the independent characteristic of noise η li and using the Taylor expansion of θl , the asymptotic efficiency Eff a ( θ) of θ in Eq. ( 15) is theoretically derived in Appendix B as

Eff a ( θ)= 1 J(f w ) L E 2 z [ψ ′ (z)] E z [ψ 2 (z)] + (L -1)E w E 2 η [ψ(w + η)] , (16)
with the expectations

E z (•)= •f z (x)dx and E η (•)= •f η (x)dx.
For a single M-estimator (L =1), the asymptotic efficiency Eff a ( θ) of Eq. ( 16) reduces to Eq. ( 32), and the corresponding optimal noise probability density f opt η that maximizes the asymptotic efficiency Eff a ( θ) of Eq. ( 16) is discussed in Appendix C. This optimization problem is similar to the issue of finding the optimal additive noise addressed by Chen et al. [START_REF] Chen | Theory of the stochastic resonance effect in signal detection: Part I-Fixed detectors[END_REF]- [START_REF] Chen | Noise enhanced parameter estimation[END_REF] based on Carathéodory's theorem. It is proven [START_REF] Chen | Noise enhanced parameter estimation[END_REF] that, for a general estimation parameter, the optimal noise, if one exists, can be chosen as a suitable randomization of no more than two constant vectors to reduce of the estimation variance without increasing the bias of the estimator. In Appendix C, for a location parameter under the Fisher-consistent constraint E z [ψ(z)] = 0, we provide the existence condition of the optimal noise implied by the function g(η) in Eq. (34b) and prove the symmetrical dichotomous noise in Eq. ( 25) to be the optimal one on the basis of a convex real-value vector function and the Jensen's inequality [START_REF] Balkan | CRLB based optimal noise enhanced parameter estimation using quantized observations[END_REF], [START_REF] Boyd | Convex Optimization[END_REF], [START_REF] Perlman | Jensen's inequality for a convex vector-valued function on an infinite-dimension[END_REF]. The demonstration of Appendix C is in particular complementary to the approach of Ref. [START_REF] Chen | Noise enhanced parameter estimation[END_REF] for location estimation.

However, the derived optimal additive noise may not exist for a single M-estimator in some situations where the function of g(η) in Eq. (34b) has its maximum at η =0, e.g., the Cauchy background noise f w of Eq. ( 12) and the Huber M-estimator of Eq. [START_REF] Zozor | Noise-aided processing: Revisting dithering in a Σ∆ quantizer[END_REF]. It is interesting to note in Fig. 3 that, even if the added noise is not optimal, the asymptotic efficiency Eff a ( θ) can be enhanced by injecting extra noise components into the M-estimator of Eq. ( 15) with the estimator number L>1.This noise-enhanced effect motivates us to further consider the situation of L>1.

Furthermore, for a given M-estimator function ψ and at a given extra noise level, we demonstrate in Appendix D that the asymptotic efficiency Eff a (θ) in Eq. ( 16) is a monotonically increasing function of the estimator number L, as illustratively shown in Fig. 3. Here, the additive noise distribution f η is arbitrary. We recognize that, due to the nonlinear functional term E w E 2 η [ψ(w + η)] in Eq. ( 16), the theoretical determination of the optimal probability density that maximizes the asymptotic efficiency of Eq. ( 16) remains an open question for a finite Mestimator number 1 <L<∞. In the limiting case of L →∞, the asymptotic efficiency of Eq. ( 16) attains its maximum with respect to L and asymptotically approaches to

Eff ∞ ( θ) = lim L →∞ Eff a ( θ)= 1 J(f w ) E 2 z [ψ ′ (z)] E w E 2 η [ψ(w + η)] , (17)
where

E z [ψ 2 (z)] < ∞ and lim L →∞ E z [ψ 2 (z)]/L → 0.I n Appendix E
, aiming to maximize the asymptotic efficiency of Eq. ( 17), we find the optimal noise density f opt η is the solution of the weighted minimum L 2 -norm of Eq. ( 44), and the asymptotic efficiency Eff ∞ ( θ) is theoretically proven to be not more than unity by the Cauchy-Schwarz inequality. Then, for achieving the upper bound of Eff ∞ ( θ)=1, it is demonstrated that the optimal noise probability density f opt η of Eq. ( 42) is the deconvolution of ψ M (w) and ψ(x). Here, the condition of existence of this optimal probability density f opt η is that the inverse of the Fourier transform in Eq. ( 42) must exist and be nonnegative. In Appendix E, an achievable optimum as a logistic noise is worked out as an example. Unfortunately, this kind of optimal noise often does not exist. For instance, there is no solution of Eq. ( 42) for the Huber function in Eq. ( 8) (or the bisquare function in Eq. ( 9)) and the Cauchy noise distribution f w in Eq. ( 18) (v =1/2 and k =2). In addition, the infinite M-estimator number L = ∞ is also inaccessible and can only be approached in practice.

As an alternative, some prior interesting parametric classes of noise distributions can be selected and optimized to improve the asymptotic efficiency of M-estimators. For example, consider the generalized Gaussian noise η i with its probability density

f η (x)= c 1 (β) σ η exp -c 2 (β) x σ η 2/(1+β ) , (18)
where 1+β ) , the decay exponent β ≥-1 and the scale parameter σ η > 0 [START_REF] Kay | Fundamentals of Statistical Signal Processing: Detection Theory[END_REF]. This generalized Gaussian noise model of Eq. ( 18) describes many noise types encountered in real-world systems, ranging from Gaussian noise (β =0), Laplacian noise (β =1)t o uniform noise (β = -1) [START_REF] Kassam | Signal Detection in Non-Gaussian Noise[END_REF], [START_REF] Kay | Fundamentals of Statistical Signal Processing: Detection Theory[END_REF]. For β = -1 in Eq. ( 18), the probability density of uniform noise can be expressed as

c 1 (β)=Γ 1/2 (3(1 + β)/2)/(1 + β)Γ 3/2 ((1 + β)/2), c 2 (β) = [Γ(3(1 + β)/2)/Γ((1 + β)/2)] 1/(
f η (x)=1/(2d) for -d ≤ x ≤ d (d = √ 3σ η
) and zero otherwise.

The theoretical asymptotic efficiency of Eq. ( 17) is then applicable to evaluate the benefit of added noise for various background noise types. For instance, with the addition of uniform noise, the asymptotic efficiency in Eq. ( 17) can be expressed as

Eff ∞ ( θ, d)= 1 J(f w ) E 2 w d -d 1 2d ψ ′ (w + η)dη E w d -d 1 2d ψ(w + η)dη 2 = 1 J(f w ) E 2 w [ψ(w + d) -ψ(w -d)] E w [ρ(w + d) -ρ(w -d)] 2 . ( 19
)
In Appendix F, we assume the uniform noise level d → 0, and derive the condition of Eq. ( 46) to determine whether the addition of uniform noise is beneficial or not. For the generalized Cauchy noise in Eq. ( 12) (k =2) with the fixed scale σ =1,the signs of D(ψ, f w ) in Eq. ( 46) versus the parameter γ of the ψ function and the parameter v of the generalized Cauchy distribution f w in Eq. ( 12) are shown in Fig. 2(a) and (b). It is seen in Fig. 2 that, using the discriminant inequality of Eq. ( 46), the regions S + of the positive signs of D(ψ, f w ) clearly manifest that the benefit of adding uniform noise exists for substantial ranges of the parameter pair (γ, v). In addition, the benefit of added uniform noise occurs for various noise distributions characterised by the parameter v of the generalized Cauchy distribution f w in Eq. [START_REF] Chapeau-Blondeau | Noise-enhanced performance for an optimal Bayesian estimator[END_REF].

The discriminant inequality of Eq. ( 46) only theoretically demonstrates the noise benefit for a very small noise level and an infinite M-estimator number. Furthermore, for the background Cauchy noise and the estimator parameter γ indicated in Fig. 2 that elicits the benefit of adding noise, the asymptotic efficiency Eff a ( θ) as a function of the uniform noise level d is shown in Fig. 3 for different M-estimator numbers. Here, the Huber function in Fig. 3(a) has a given estimator parameter γ =0.1, and the bisquare function in Fig. 3(b) is with γ =1. The fixed scale of the Cauchy noise is σ =1. It is seen in Fig. 3 that, for a single M-estimator (L =1), the asymptotic efficiency Eff a ( θ) monotonically decreases with the increase of noise level d.While,for the estimator number L>1, the asymptotic efficiency Eff a ( θ) can be improved to a maximum as the noise level d increases. This kind of noise enhancement is distinct from the beneficial effect of noise in a single processor in that it critically relies on the collective role of the estimator number L [START_REF] Stocks | Suprathreshold stochastic resonance in multilevel threshold systems[END_REF]. Then, the noise-enhanced effect shown in Fig. 3 can be regarded as a form of SSR [START_REF] Stocks | Suprathreshold stochastic resonance in multilevel threshold systems[END_REF] for location estimation in essence.

We also numerically demonstrate the noise-enhanced asymptotic efficiency effect for the bisquare M-estimator by the iterative reweighting method [START_REF] Maronna | Robust Statistics: Theory and Methods[END_REF]. The procedure of the numerical compution is given in Appendix G-A. In Fig. 3(b), the asymptotic efficiency Eff a ( θ) is numerically evaluated for 10 4 simulation trails and the tolerance parameter ζ =10 -6 . Here, the observation size n = 7000 is large enough but finite. The numerical results of the asymptotic efficiency Eff a ( θ) are plotted in Fig. 3(b) for the M-estimator numbers L =1( * ), L =3 (▽), L =7(•), L =15( ) and L = 100 (⋆).I ti ss h o w ni n Fig. 3(b) that the numerical results of Eff a ( θ) obtained by this method are well consistent with the theoretical values of Eff a ( θ) in Eq. ( 16). Now we return to the question of whether, for an optimized M-estimator with its optimal parameter γ opt , the addition of extra noise is beneficial or not. For a fixed input Cauchy noise Fig. 5. Asymptotic efficiency Eff( θ) of Eq. ( 11) without added noise (solid line), and the maximum asymptotic efficiency Eff a ( θ) of Eq. ( 16) versus the estimator parameter γ for Gaussian (•), Laplacian ( * ) and uniform (△) types of the added noise from Eq. [START_REF] Akbay | Noise benefits in joint detection and estimation problems[END_REF]. For each point of the curves, the level σ η of added noise is optimized to achieve the maximum asymptotic efficiency Eff a ( θ). Here, the background Cauchy noise scale σ =1and the M-estimator number L = 1000.

scale σ =1 and the M-estimator number L = 1000,w ep l o t the asymptotic efficiency Eff a ( θ) of the bisquare M-estimator array as a function of both the estimator parameter γ and the added noise level d, as shown in Fig. 4(a). It is found in Fig. 4(a) that the asymptotic efficiency Eff a ( θ) of Eq. ( 16) attains its maximum value of 0.93 ( * ) at the optimized parameter pair of (d opt ,γ opt )=(1.3, 2.8). Without the help of added noise (d = 0), the maximum Eff( θ) is 0.90 ( ) with the optimal parameter γ opt =3.28, as shown in Fig. 4(a). The maximum asymptotic efficiency Eff a ( θ)=0.93 via optimizing the parameter pair (γ, d) presents a slight improvement over the maximum Eff( θ) of 0.90 without the help of added noise.

Furthermore, we can consider the possibility of adding extra noise from Eq. ( 18) to the M-estimator. For different decay exponents β of the generalized Gaussian distributions of Eq. ( 18), Fig. 4(b) shows the maximum asymptotic efficiency Eff a ( θ) in Eq. ( 16) via optimizing the parameter pair (γ, σ η ). It is interesting to note in Fig. 4(b) that, for a finite M-estimator number L = 1000, the maximum asymptotic efficiency of Eq. ( 16) can reach 0.956 (▽) for the decay exponent β =1(Laplacian noise). More interestingly, based on Eq. ( 17) and for L = ∞, Fig. 4(b) also shows the the maximum asymptotic efficiency Eff ∞ ( θ) can achieve as high as 0.999 (•) that is very close to unity. Here, the decay exponent β =1.2 and the optimal parameter pair(γ opt ,σ opt η )=(1.4, 8.1). Of course, this limit value of Eff ∞ ( θ)=0.999 can only be approached in practice. Under this circumstance, Fig. 4(c) illustrates the asymptotic efficiency Eff a ( θ) as a monotonically increasing function of the M-estimator number L with the parameters γ opt , σ opt η and β that correspond to the maximum asymptotic efficiency Eff ∞ ( θ)=0.999. It is shown in Fig. 4(c) that the asymptotic efficiency Eff a ( θ) can gradually approach the limit value of Eff ∞ ( θ)=0.999, but at the cost of the very heavy computation complexity induced by the large number L.

Finally, we answer the second question of the possibility of a noise benefit to an array of robust M-estimators with a fixed parameter γ. Fig. 5 shows the asymptotic efficiency Eff( θ) (solid line) of Eq. ( 11) without added noise. Furthermore, for each given parameter γ and the estimator number L = 1000, Fig. 5 also depicts the maximum asymptotic efficiency Eff a ( θ) of Eq. ( 16) for Gaussian (•), Laplacian ( * ) and uniform (△) types of the added noise via optimizing the noise level σ η (d = √ 3σ η for uniform noise). From Fig. 5, it is visible that the asymptotic efficiencies Eff a ( θ) obtained by optimizing the level of added noise are always superior to the efficiency with no added noise (solid line). These results show that, for a fixed parameter γ,the addition of extra noise is effectively beneficial to the estimation process.

IV. ADDING NOISE TO REDUCE THE MAXIMUM BIAS

In Section III, we discussed the noise-enhanced asymptotic efficiency in symmetric heavy-tailed noise models. In this Section, we further consider the case of ǫ-contaminated noise model h ranging in the ǫ-neighborhood [START_REF] Huber | Robust Statistics[END_REF], [START_REF] Maronna | Robust Statistics: Theory and Methods[END_REF] h ǫ (x, θ)={h| h =(1-ǫ)f w (x, θ)+ǫq(x),q ∈Q}, (20) of an assumed parametric distribution f w (x, θ), where q that represents the probability model of outliers belongs to a family of probability density functions Q and 0 ≤ ǫ<1 indicates the proportion of the observations generated by the contamination model of q [34], [START_REF] Maronna | Robust Statistics: Theory and Methods[END_REF]. Here, the contamination model of q is usually not exactly known, and depicts the statistical characteristics of the outliers produced by some unknown mechanisms [START_REF] Huber | Robust Statistics[END_REF], [START_REF] Kassam | Robust techniques for signal processing: A survey[END_REF]- [START_REF] Zoubir | Robust estimation in signal processing: A tutorial-style treatment of fundamental concepts[END_REF]. In order to analyze some extreme contaminated cases, q is chosen from the set of point mass distribution δ(x 0 ) such that the probability of x = x 0 is unity [START_REF] Huber | Robust Statistics[END_REF], [START_REF] Maronna | Robust Statistics: Theory and Methods[END_REF]. Here, δ is the Dirac function. It has been pointed out [START_REF] Huber | Robust Statistics[END_REF], [START_REF] Maronna | Robust Statistics: Theory and Methods[END_REF] that, for the location M-estimation in neighborhoods of Eq. ( 20) of a symmetric f w with a symmetric probability density q of contamination, it is easy to derive an unbiased M-estimator θ with E( θ)=θ and we still care about the asymptotic efficiency of M-estimators. Thus, in this Section, we mainly focus on the asymmetric contaminated noise model of h ǫ (x, θ). Under this condition, the bias of the estimators to the true value becomes the main factor of MSE as a measurement of estimators rather than the variance for a sufficiently large observation size n.This is because the estimation variance var( θ) of Eq. ( 4) tends to zero with the order of n -1 , while the bias does not [START_REF] Huber | Robust Statistics[END_REF]- [START_REF] Zoubir | Robust estimation in signal processing: A tutorial-style treatment of fundamental concepts[END_REF].

For the asymmetric contaminated noise model of h ǫ (x, θ) of Eq. ( 20), there are two cases: either a symmetric central probability density f w but an asymmetric point mass distribution q,or an asymmetric f w plus an asymmetric (or symmetric) q.Inthe case of the noise components w i with the asymmetric probability density f w , we can add the i.i.d. noise components w ′ i with the probability density of -w i , i.e., f w ′ (x)=f w (-x),t ot h e original noise components w i . Here, w ′ i are independent with w i . Then, the updated random variables wi = w i + w ′ i have a symmetric probability density f w (x)= f w (x -u)f w ′ (u)du = f w ′ (-x -u)f w (u)du = f w (-x). Therefore, we here only concentrate on the maximum bias induced by the contaminated noise model of h ǫ (x, θ) with a symmetric f w but an asymmetric q in the following parts.

Consider the estimate θn depending on the observation data (x 1 ,x 2 , ••• ,x n ) and for a location M-estimator with the monotonic ψ in Eq. ( 3), the asymptotic value θ∞ = θ under f w , i.e., the solution of E[ψ(x -θ)] = 0 [START_REF] Huber | Robust Statistics[END_REF], [START_REF] Maronna | Robust Statistics: Theory and Methods[END_REF]. Then, the bias for an arbitrary distribution h is b θ = θ∞ (h) -θ. For a given proportion ǫ and the assumed parametric probability density f w , the maximum bias is defined as [START_REF] Huber | Robust Statistics[END_REF], [START_REF] Maronna | Robust Statistics: Theory and Methods[END_REF]

b ǫ = max{|b θ (h, θ)| : h ∈ h ǫ }, (21)
over an arbitrary contamination model of q ∈Q. Let ψ be a nondecreasing function, instead of a redescending function, with a bound s = ψ(∞), the fraction of contamination 0 <ǫ<1/2 and θ =0(without loss of generality). It can be shown [START_REF] Huber | Robust Statistics[END_REF], [START_REF] Maronna | Robust Statistics: Theory and Methods[END_REF] that the maximum bias b ǫ of a location Mestimator is the solution of

E w [ψ(x + b)] = sǫ 1 -ǫ , ( 22
)
which is attained in the case of the point mass distribution q = δ(x 0 ) in Eq. ( 20) with x 0 →±∞. The median estimator, as a special Huber estimator in Eq. ( 8) with its parameter γ → 0, has the minimax bias property for a symmetric unimodal distribution f w of noise [START_REF] Huber | Robust Statistics[END_REF], [START_REF] Maronna | Robust Statistics: Theory and Methods[END_REF], [START_REF] Zoubir | Robust estimation in signal processing: A tutorial-style treatment of fundamental concepts[END_REF]. Since the median estimator with ψ(x)=sgn(x) is bounded by s =1, the maximum bias b ǫ in Eq. ( 22) can be calculated as [START_REF] Huber | Robust Statistics[END_REF], [START_REF] Maronna | Robust Statistics: Theory and Methods[END_REF]

b ǫ = F -1 w [1/2(1 -ǫ)], (23)
where F -1 w is the inverse function of the cumulative distribution F w . When we add noise η i to the median estimator, the observations in Eq. ( 13) are now buried in the composite noise z i = w i + η i (i =1, 2, ••• ,n). Then, the maximum bias of Eq. ( 23) is updated as

b ǫ = F -1 z [1/2(1 -ǫ)], (24)
with the inverse function F -1 z of the cumulative distribution F z . In Appendix H, we derive the condition for reducing maximum bias b ǫ of the median estimator in Eq. ( 23) with additive noise. It is demonstrated that if there exists a positive interval where the cumulative distribution function F w is strictly convex, then the maximum bias b ǫ of the median estimator in Eq. ( 23) can be reduced. Furthermore, under the condition of convexity of F w , we also demonstrate in Appendix H that the optimal noise probability density f opt η is given by

f opt η (η)= 1 2 [δ(η -a opt )+δ(η + a opt )], (25)
with the optimal noise level

a opt = arg max a 1 2 [F w (b ǫ -a)+F w (b ǫ + a)]. ( 26
)
The noise with the distribution in Eq. ( 25) is also referred to as the dichotomous noise [START_REF] Chapeau-Blondeau | Theory of stochastic resonance in signal transmission by static nonlinear systems[END_REF]. This condition applies for bimodal noise distribution models arising in diverse problems [START_REF] Baldry | Quantifying the bimodal color-magnitude distribution of galaxies[END_REF]- [START_REF] Lategahn | Texture classification by modeling joint distributions of local patterns with Gaussian mixtures[END_REF], where the cumulative distribution function F w satisfies the convexity condition given in Appendix H. For example, consider the Gaussian mixture noise with density

f w (x)= 1 2 [κ(x; -µ, τ )+κ(x; µ, τ )] , (27)
where κ(x; µ, τ

)=(1/ √ 2πτ 2 )exp[-(x -µ) 2 /(2τ 2 )
] with parameters µ, τ ≥ 0. The dichotomous noise of Eq. ( 25) with the optimal noise level a opt of Eq. ( 26) is added to the observations. For the fraction of contamination ǫ =0.1,F i g .6 ( a )p r e s e n t s that the maximum bias b ǫ (solid line) of Eq. ( 23) without added noise increases as the ratio of µ/τ grows, while the maximum bias b ǫ (•) of Eq. ( 24) obtained by optimizing the added noise level a opt maintains a lower value with the increase of µ/τ . The corresponding optimal noise level a opt is also plotted in Fig. 6(b). The results clearly show the noise-reduced bias effect in the median estimator for the asymmetric ǫ-contaminated noise model.

V. D ISCUSSION

In this paper, we study the enhancement of the asymptotic efficiency and the reduction of the maximum bias by adding noise in robust location M-estimators. For symmetric heavy-tailed noise models, we show that the asymptotic efficiencies of two commonly used M-estimators with the Huber function and the bisquare function are non-monotonic functions of the ratio of the noise scale and the estimator parameter. With a summing array of M-estimators, by injecting extra noise into the array, we show the possibility of noise enhancement of the asymptotic efficiency that is maximized by an optimal nonzero amount of extra noise.

Furthermore, in order to maximize the asymptotic efficiency of a single M-estimator with added noise, a discriminant function is proposed to establish whether an optimal noise condition exists or not. When the optimal additive noise exists, its probability density function is proven to be the same one given by Chen et al. [START_REF] Chen | Theory of the stochastic resonance effect in signal detection: Part I-Fixed detectors[END_REF]- [START_REF] Chen | Noise enhanced parameter estimation[END_REF]. Nevertheless, even when the optimal additive noise does not exist, the noise-enhanced effect occurs for an array of M-estimators with estimator number larger than one. Moreover, it is proven that the asymptotic efficiency of an array of M-estimators is a monotonically increasing function of the M-estimator number, and achieves its maximum for an infinite M-estimator number. Then, aiming to maximize the asymptotic efficiency of an infinite number of M-estimators, we find that the optimal noise distribution is the weighted minimum L 2norm solution of Eq. [START_REF] Perlman | Jensen's inequality for a convex vector-valued function on an infinite-dimension[END_REF]. It is also theoretically demonstrated that the upper bound of the asymptotic efficiency of an infinite number of M-estimators is unity, and the corresponding optimal noise is analytically derived. Unfortunately, this optimal noise distribution often does not exist. Therefore, we select a class of parametric noise distributions, and optimize the noise distribution parameters as well as the noise variance. Specially, we present a condition for a benefit by adding uniform noise, and show the noise benefit is effective for a number of background heavy-tailed noise models.

It is also shown that, compared with the method without extra noise, the addition of noise can provide an improved asymptotic efficiency with a local optimal probability density chosen from the given parametric noise distributions. In addition, upon increasing the M-estimator number to a sufficiently large value, we can gradually obtain an excellent asymptotic efficiency close to unity. The asymmetric contaminated model is also considered, which contains a proportion of the observations generated by the true distribution, while another proportion of outliers is mixed by an unknown mechanism. In this case, we consider the extreme situation of contamination with an infinite point-mass distribution, and characterize the noise type for which the maximum bias of the median estimator can be minimized by adding an optimal amount of the dichotomous noise.

Several questions of interest remain open: Although the optimal noise density f opt η that maximizes the asymptotic efficiency of an infinite number of M-estimators is proven to be the weighted minimum L 2 -norm solution of Eq. ( 44), an explicit form of f opt η remains an interesting open question for future study. We also note that, for achieving the upper bound of unity of the asymptotic efficiency, the optimal additive noise exists only if the inverse of the Fourier transform of Eq. ( 42) is nonnegative over all the probability space. Otherwise, the optimal noise does not exist. In cases where the inverse of the Fourier transform of Eq. ( 42) exists but has negative parts, one can envisage to clip the negative parts to zero and renormalize the positive parts as a suboptimal solution for the noise density. We test this approach with the following example. Consider a zero-mean Student distribution f w (x)= c v (1 + x 2 /v) -(v +1)/2 with the degree of freedom v>0 and the constant c v =Γ((v +1)/2)/( √ vπΓ(v/2)). The corresponding maximum likelihood estimator ψ M = -f ′ w (x)/f w (x)= (v +1)x/(v + x 2 ). When the median estimator ψ(x)=sgn(x) is employed, the inverse of the Fourier transform of Eq. ( 42) is given by I(x)= π/2(v + 1)(v -x 2 )/(v + x 2 ) 2 , which has negative parts as |x| > √ v and is actually not a probability density function. We can clip these negative parts of I(x) to zero and renormalize the nonnegative parts as a probability density function

f c η (x)= √ v(v -x 2 )/(v + x 2 ) 2 for |x|≤ √ v and
zero otherwise. It is very interesting to note that, with this clipped probability density function f c η and the degree of freedom v =4, the asymptotic efficiency of Eq. ( 17) is 0.995 that is very close to unity. Thus, this suboptimal solution f c η (x) provides a suitable noise density to improve the asymptotic efficiency of an array of M-estimators. Adding noise in an array as we propose is another alternative to improve the efficiency. However, even if this optimal noise distribution exists, the infinite-size array can only be approached in practice. As an alternative, it is meaningful to further explore the optimal noise type that achieves the maximum asymptotic efficiency for a finite M-estimator number, and ideally with a noise type which is not restricted to a prior class of parametric noise distributions. In addition, for the special case of the median estimator, the optimal dichotomous noise is demonstrated to minimize the maximum bias for estimating a location buried in bimodal noise models with convex cumulative distributions. It is interesting to analyze the reduction of maximum bias with addition of noise for other robust M-estimators, and prove whether there are other contaminated models exhibiting the noise-enhanced effects or not. We here consider the asymmetrical contaminated model of Eq. ( 20) for the point mass distribution of extreme cases, and derive the optimal symmetrical noise to reduce the maximum bias at a given proportion of the contamination. In this contaminated model, the contaminating probability density q is actually arbitrary and unknown. But, with sufficient knowledge of the contamination distribution in Eq. [START_REF] Soganci | Optimal signal design for multiparameter estimation problems[END_REF], it is very meaningful to further investigate the optimal asymmetrical noise obeying the distribution of h ǫ (-x). Then, the asymmetrical contaminated model of Eq. ( 20) can be simplified as an unbiased location estimation. For robust statistics of multi-parameter problems and M-estimators for regression models [START_REF] Huber | Robust Statistics[END_REF]- [START_REF] Zoubir | Robust estimation in signal processing: A tutorial-style treatment of fundamental concepts[END_REF], the benefits of noise also deserve to be further studied. These questions are of interest for robust estimation where addition of noise properly optimized may reveal, as in the present study, to be a beneficial option.

APPENDIX A NO NOISE-ENHANCED EFFECT FOR THE MAXIMUM LIKELIHOOD ESTIMATOR

In this Appendix, we show that adding extra noise to observations cannot reduce the variance of the maximum likelihood estimator. Based on Eqs. ( 5) and ( 6), the asymptotic variance of the maximum likelihood estimator is 1/(nJ(f w )) that is closely related to the original noise distribution f w . When we add extra noise η to the observations and the background noise becomes the composite noise z = w + η, the asymptotic variance of the maximum likelihood estimator that is optimal for noise z is updated as 1/(nJ(f z )). According to the Fisher information convolution inequality J(f z ) ≤ min{J(f w ),J(f η )} in Refs. [START_REF] Cover | Elements of Information Theory[END_REF] (Chapter 16.6) and [START_REF] Blachman | The convolution inequality for entropy power[END_REF], we can derive that 1/(nJ(f z )) ≥ 1/(nJ(f w )), which indicates no noise-enhanced effect for the maximum likelihood estimator.

APPENDIX B ASYMPTOTIC EFFICIENCY OF AN ARRAY OF M-ESTIMATORS

In this Appendix, we present the derivation of the asymptotic efficiency of Eq. ( 16) for a summing array of M-estimators. The location estimate θl of each M-estimator satisfies n i=1 ψ(y liθl )=0for l =1, 2, ••• ,L. We can perform a Taylor expansion of ψ as a function of θl around the true value of θ,sotheresult yields [START_REF] Huber | Robust Statistics[END_REF], [START_REF] Maronna | Robust Statistics: Theory and Methods[END_REF] 0= n i=1 ψ(y li -θl )

≈ n i=1 ψ(y li -θ) -( θl -θ) n i=1 ∂ ∂θ ψ(y li -θ) = A n -( θl -θ)B n , (28)
where A n = n i=1 ψ(y li -θ)/n and B n = n i=1 (∂ψ(y liθ)/∂θ)/n. Hence, by Slutsky's lemma and according to the central limit theorem for a sufficiently large n, B n in probability converges to B =E[∂ψ(y -θ)/∂θ]=E z [ψ ′ (z)], and A n converges to a Gaussian distribution with zero mean and variance E[ψ 2 (y -θ)]/n =E z [ψ 2 (z)]/n [START_REF] Huber | Robust Statistics[END_REF], [START_REF] Maronna | Robust Statistics: Theory and Methods[END_REF]. Therefore, each estimate θl converges to a Gaussian distribution with mean θ and variance var( θl )=n -1 E z [ψ 2 (z)]/E 2 z [ψ ′ (z)] of Eq. [START_REF] Zozor | Stochastic resonance in locally optimal detectors[END_REF]. In line with this general proof, the average estimate θ = L l=1 θl /L, also converges to a Gaussian distribution with mean θ and variance

var a ( θ)=v ar 1 L L l=1 θl = 1 L 2 L var( θl )+ L l,p=1 l = p E z [( θl -θ)( θp -θ)] = 1 L 2 L var( θl )+L(L -1)E z [( θl -θ)( θp -θ)] = 1 L var( θl )+(L -1)E w E η [( θl -θ)( θp -θ)] = 1 L var( θl )+ L -1 L E w E η [( θl -θ)( θp -θ)] . (29)
Since for a given w i , w i + η li are i.i.d. random variables for l =1, 2, ••• ,L and the Taylor expansion of ψ i nE q .( 2 8 ) ,s o we can obtain

E w E η [( θl -θ)( θp -θ)] = E w E 2 η [ψ(w + η)] n E 2 z [ψ ′ (z)] . (30)
Substituting Eqs. ( 5) and (30) into Eq. ( 29), the asymptotic variance of θ can be calculated as

var a ( θ)= 1 n E z [ψ 2 (z)] + (L -1)E w E 2 η [ψ(w + η)] L E 2 z [ψ ′ (z)]
.

Then, the asymptotic efficiency of Eq. ( 16) is obtained.

APPENDIX C OPTIMIZATION OF THE ASYMPTOTIC EFFICIENCY OF A SINGLE M-ESTIMATOR

For a single M-estimator (L =1), the asymptotic efficiency of Eff a ( θ) of Eq. ( 16) becomes

Eff a ( θ)= 1 J(f w ) E 2 z [ψ ′ (z)] E z [ψ 2 (z)] . (32)
Thus, it is interesting to find an optimal noise probability density f opt η to maximize the asymptotic efficiency of Eff a ( θ) in Eq. ( 32) under the constraint of unbiasedness. Since the Fisher information J(f w ) > 0, then this optimization problem can be expressed as

max f η E 2 z [ψ ′ (z)] E z [ψ 2 (z)] , (33)
with the constraint of E z [ψ(z)] = 0. Consider a real-value vector function f

(Z)=Z 2 1 /Z 2 with Z =[Z 1 ,Z 2 ] T and Z 2 > 0. It is seen that f (Z) is convex, be- cause its Hessian matrix ∇ 2 f (Z)=2[Z 2 , -Z 1 ] T [Z 2 , -Z 1 ]/Z 3 2
is positive semidefinite [START_REF] Boyd | Convex Optimization[END_REF]. Then, Jensen's inequality [START_REF] Cover | Elements of Information Theory[END_REF], [START_REF] Boyd | Convex Optimization[END_REF], [START_REF] Perlman | Jensen's inequality for a convex vector-valued function on an infinite-dimension[END_REF]. Therefore, defining

E 2 η (Z 1 )/E η (Z 2 2 ) ≤ E η (Z 2 1 /Z 2 ) is valid in the definition do- main {(Z 1 ,Z 2 ) ∈ R 2 |Z 2 > 0} [17],
Z 1 =E w [ψ ′ (w + η)] and Z 2 =E w [ψ 2 (w + η)],w e have E 2 z [ψ ′ (z)] E z [ψ 2 (z)] = E 2 η {E w [ψ ′ (w + η)]} E η {E w [ψ 2 (w + η)]} = E 2 η (Z 1 ) E η (Z 2 ) ≤ E η Z 2 1 /Z 2 (34a) =E η [g(η)] (34b) ≤ max η g(η), (34c)
where g(η)=Z 2 1 /Z 2 , Z 2 1 and Z 2 are even functions for the odd function ψ and the even probability density f w . It is noted that E η [g(η)] cannot be larger than the maximum value of g(η) in Eq. (34c) and the optimal added noise has a density function of one or multiple mass points (depending on the number of points of a opt = arg max η g(η)). Furthermore, since g(η) is an even function, then it must have an extremum value g(0) at the origin. If g(0) is the absolute maximum over the real space η ∈ R, then the optimal solution of Eq. ( 33) is f opt η (η)=δ(η). This indicates that there is no benefit in adding noise. Meanwhile, if g(η) has two absolute maximum values of g(a opt )=g(-a opt ) > g(0) at ±a opt (a opt =0), then the optimal dichotomous noise density exists and can be expressed by Eq. ( 25) for obtaining the maximum of g(a opt ). Substituting Eq. ( 25) with the optimal noise level a opt into Eq. (34a) and Eq. (34c), both equalities hold and the maximum value of Eq. ( 33) can be achieved. This result accords with the determination of the optimal noise proven by Chen et al. [START_REF] Chen | Theory of the stochastic resonance effect in signal detection: Part I-Fixed detectors[END_REF]- [START_REF] Chen | Noise enhanced parameter estimation[END_REF] based on Carathéodory's theorem [START_REF] Roberts | Convex Functions[END_REF].

APPENDIX D ASYMPTOTIC EFFICIENCY VERSUS THE M-ESTIMATOR NUMBER

For the convex function x 2 and by the Jensen inequality, we obtain

E η [ψ 2 (w + η)] ≥ E 2 η [ψ(w + η)], (35)
for any fixed variable w [START_REF] Cover | Elements of Information Theory[END_REF]. Then, we have

E z [ψ 2 (z)] = E w E η [ψ 2 (w + η)] ≥ E w E 2 η [ψ(w + η)] . (36)
Immediately, we find

1 L +1 E z [ψ 2 (z)] + L L +1 E w E 2 η [ψ(w + η)] ≤ 1 L E z [ψ 2 (z)] + L -1 L E w E 2 η [ψ(w + η)] . (37)
From Eq. ( 16) and Eq. ( 37), we derive

1 J(f w ) (L + 1)E 2 z [ψ ′ (z)] E z [ψ 2 (z)] + LE w E 2 η [ψ(w + η)] ≥ 1 J(f w ) LE 2 z [ψ ′ (z)] E z [ψ 2 (z)] + (L -1)E w E 2 η [ψ(w + η)] , which indicates Eff a ( θ, L +1)≥ Eff a ( θ, L). (38)
Thus, for the given M-estimator function ψ and fixed noise components w and η, the asymptotic efficiency Eff a ( θ, L) in Eq. ( 16) is a monotonically increasing function of the M-estimator number L.

APPENDIX E OPTIMIZATION OF THE ASYMPTOTIC EFFICIENCY FOR AN INFINITE NUMBER OF M-ESTIMATORS

By using the Cauchy-Schwarz inequality, we find that the asymptotic efficiency Eff ∞ ( θ) in Eq. ( 17) satisfies

Eff ∞ ( θ)= 1 J(f w ) E 2 w {dE η [ψ(w + η)]/dw} E w E 2 η [ψ(w + η)] = 1 J(f w ) E 2 w {E η [ψ(w + η)] ψ M (w)} E w E 2 η [ψ(w + η)] ≤ 1 J(f w ) E w ψ 2 M (w) =1, (39)
where the equality occurs as

E η [ψ(w + η)] = Cψ M (w), (40)
with an arbitrary constant C.

Thereby, for the given function ψ and the background noise distribution f w , the optimal additive noise probability density f opt η is the solution of Eq. [START_REF] Mcdonnell | Stochastic Resonance: From Suprathreshold Stochastic Resonance To Stochastic Signal Quantization[END_REF]. With the assumptions of symmetrical distribution function f w and the odd M-estimator function ψ, the added noise probability density function is also even, i.e., f η (-η)=f η (η), due to the Fisher-consistent M-estimator satisfying E z [ψ(z)] = 0 [START_REF] Maronna | Robust Statistics: Theory and Methods[END_REF]. Then, Eq. ( 40) can be rewritten as

E η [ψ(w + η)] = + ∞ -∞ ψ(w + η)f η (-η)dη = -∞ + ∞ ψ(w -u)f η (u)d(-u) = ψ(w) * f η (w) = Cψ M (w), (41)
which indicates that the optimal noise probability density f opt η is the deconvolution of Cψ M (x) and ψ(x). If the Fourier transformations of ψ(x) and ψ M (x) both exist [START_REF] Corduneanu | Integral Equations and Applications.N e wY o r k[END_REF], [START_REF] Oppenheim | Signals and Systems[END_REF], we have

f opt η (x)=C F -1 F [ψ M (x)] F [ψ(x)] , (42)
where F and F -1 denote the Fourier transform and its inverse transform, respectively. Here, we emphasize that, due to the probability density f η (η) ≥ 0, the optimal density f opt η exists only if the inverse of the Fourier transform in Eq. ( 42) exists and is a nonnegative function. The constant C can be calculated by the normalization condition f η (η)dη =1. In addition, both ψ M (x) and ψ(x) are odd functions. According to the convolution property [START_REF] Oppenheim | Signals and Systems[END_REF], f opt η is an even function when it exists. For example, consider a kind of heavy-tailed noise with the hyperbolic secant distribution f w (x)=sec h(πx/2) [START_REF] Kassam | Signal Detection in Non-Gaussian Noise[END_REF] and the median estimator ψ(x) = sign(x). Then, from Eq. ( 42), the optimal probability density becomes f opt η (x)=πe πx /[4(1 + e πx ) 2 ], which indicates logistic noise as a feasible optimal noise.

It is noted that the optimal noise indicated in Eq. ( 42) corresponding to the upper bound of unity is only a special case of

f opt η (x) = arg max f η Eff ∞ ( θ). (43)
An achievable optimal noise density that achieves the maximum of the asymptotic efficiency (not more than unity) of Eq. ( 17) among all feasible densities is practically significant. Define the inner product of two functions f and g as f, g = E w [f (w)g(w)], and the L 2 -norm of a function g as g = g, g 1/2 . Here, we define φ(w)=E η [ψ(w + η)] and its normalized function can be expressed as φ(w)=φ(w)/ φ(w) . Correspondingly, the maximum likelihood estimator ψ M (w) can be also normalized as ψ M (w)=ψ M (w)/ ψ M (w) with its norm ψ M (w) = J(fw). From Eqs. ( 17) and ( 41), we have

Eff ∞ ( θ)= E η [ψ(w + η)],ψ M (w) 2 E η [ψ(w + η)] 2 ψ M (w) 2 = φ(w), ψ M (w) 2 .

Note that φ(w) -ψ M (w) 2 =2-2 φ(w), ψ M (w) . Thus, the optimal noise density indicated in Eq. ( 43) is just the minimum L 2 -norm solution of φ(w) -ψ M (w) 2 , expressed as

f opt η (x) = arg min f η φ(w) -ψ M (w) 2 , (44)
with the constrains of f η (x)dx =1and f η (x) ≥ 0 for x ∈ R.

A further point is to consider the set of all possible noise densities, and then in this set there is necessarily at least one that realizes the maximum efficiency, which is equivalent to solving the constrained minimization of Eq. ( 44). How to calculate such an achievable optimal density of the additive noise remains an interesting open question for future study.

APPENDIX F BENEFIT CONDITION OF ADDING UNIFORM NOISE

Based on Eq. ( 19), the derivative of Eff ∞ ( θ, d) with respect to the uniform noise level d is given by w [ρ(w + d) -ρ(w -d)] 2 > 0, then the optimal noise level d opt is the solution of A =0. We further consider the small noise level d → 0 (d>0). In this situation, the derivative of Eff ∞ ( θ, d) in Eq. ( 45) can be approximated with the second-order Taylor expansion as

∂Eff ∞ ∂d ≈ d J(f w ) E 2 w [ψ ′ (w)] E w [ψ 2 (w)] E w [ψ ′′′ (w)] E w [ψ ′ (w)] - E w [ψ(w)ψ ′′ (w)] E w [ψ 2 (w)]
.

Since E 2 w [ψ ′ (w)]/E w [ψ 2 (w)] > 0 and the noise level d>0, then the condition of a benefit by adding uniform noise is given by

D(ψ, f w )= E w [ψ ′′′ (w)] E w [ψ ′ (w)] - E w [ψ(w)ψ ′′ (w)] E w [ψ 2 (w)] > 0, (46)
which indicates when the addition of uniform noise leads to the increase of the asymptotic efficiency Eff ∞ ( θ, d). If the second or third derivative of ψ in Eq. [START_REF] Baldry | Quantifying the bimodal color-magnitude distribution of galaxies[END_REF] does not exist, we can reduce the derivative order of ψ, and explore the derivative of f w by applying the subsection integral method. For example, E w [ψ ′′′ (w)] = -ψ ′′ (w)f ′ w (w)dw and note f w (±∞)=0.In practice, we can simply calculate the sign of D(ψ, f w ) in Eq. ( 46) to determine whether the addition of uniform noise is beneficial or not for a background noise distribution f w ,a s illustrated in Fig. 2. the optimal distribution f opt η is given by arg max

f η F z (b ǫ )=E η [F w (b ǫ -η)], (53)
with the constraint E z [ψ(z)] = E η [2F w (η) -1] = 0.Next,we will demonstrate that the solution of Eq. ( 53) is just the noise distribution f η of Eq. [START_REF] Chen | Performance limit for distributed estimation systems with identical one-bit quantizers[END_REF]. Proof: Consider the set of all pairs of (F w (b ǫ -η), 2F w (η) -1) named P = {(p 1 ,p 2 )|p 1 = F w (b ǫ -η),p 2 =2F w (η) -1}. T = {(t 1 ,t 2 )|t 1 = F z (b ǫ ),t 2 =E z [ψ(z)]} denotes the set of pairs (F z (b ǫ ), E z [ψ(z)]).L e tE be the convex hull of P and the dimension Dim(E) ≤ 2. As discussed in [START_REF] Chen | Theory of the stochastic resonance effect in signal detection: Part I-Fixed detectors[END_REF]- [START_REF] Chen | Noise enhanced parameter estimation[END_REF], [START_REF] Soganci | Optimal stochastic parameter design for estimation problems[END_REF], we have T = E. From Carathéodory's theorem [START_REF] Roberts | Convex Functions[END_REF], any point in T can be expressed as a combination of at most three elements of P . Since our aim is to maximize F z (b ǫ ) with the constraint of E z [ψ(z)] = 0, the optimum point can only occur on the boundary of T . Therefore, the optimal pair (F z (b ǫ ), E z [ψ(z)]) can be expressed as the convex combination of two points of P . Based on the constraint such that E z [ψ(z)] = 0, i.e., f η (-η)=f η (η), we derive the optimal noise distribution f η of Eq. ( 25).

Fig. 1 .

 1 Fig.1. Asymptotic efficiency Eff( θ) of (a) Eq. (10) for the Huber estimator and (b) Eq. (11) for the bisquare estimator as a function of γ/σ for the generalized Cauchy noise in Eq. (12) and the Laplacian noise in Eq. (18) (β =1).Here, the generalized Cauchy noise is with parameters v =1/2, 1, 2, 5, ∞ and k =2.

Fig. 2 .

 2 Fig.2. Region S + (gray) of the positive signs of D(ψ, f w ) in Eq.[START_REF] Baldry | Quantifying the bimodal color-magnitude distribution of galaxies[END_REF] where addition of noise is beneficial, and the complementary region S -of the negative signs of D(ψ, f w ), versus the parameter v of the generalized Cauchy distribution f w in Eq. (12) (k =2) and the parameter γ of (a) the Huber and (b) the bisquare function. Here, the existing noise scale σ =1and the added uniform noise level d → 0.

Fig. 3 .

 3 Fig. 3. Asymptotic efficiency Eff a ( θ) of Eq. (16) for an array of (a) Huber estimators with γ =0.1 and (b) bisquare estimators with γ =1as a function of uniform noise level d and the estimator number L. Here, the existing Cauchy noise is with the fixed scale σ =1. Using the iterative reweighting algorithm (see Appendix G-A), the numerical results computed by 10 4 Monte Carlo trails are plotted in (b) for the M-estimator number L =1( * ), L =3( ▽), L =7 (•), L =15( ) and L = 100 (⋆). Here, the observation size takes n = 7000, the number of simulation trails is 10 4 and the tolerance parameter ζ =10 -6 .

Fig. 4 .

 4 Fig. 4. (a) Asymptotic efficiency Eff a ( θ) of Eq. (16) for an array of bisquare estimators versus the estimator parameter γ and uniform noise level d. Here, the M-estimator number L = 1000 and the background Cauchy noise scale σ =1; (b) Maximum asymptotic efficiency Eff a ( θ) of Eq. (16) as a function of the decay exponents β of generalized Gaussian distributions via the optimized parameter pair of (σ η ,γ) for L = 1000 (solid line) and L = ∞ (dashed line). (c) Asymptotic efficiency Eff a ( θ) as a monotonically increasing function of the M-estimator number L with the parameters γ opt , σ opt η and β that correspond to the maximum asymptotic efficiency Eff ∞ ( θ)=0.999 (•) in Fig. 4(b).

Fig. 6 .

 6 Fig. 6. (a) Maximum bias b ǫ (•) of Eq. (24) obtained by optimizing the added noise level a, and maximum bias b ǫ (solid line) of Eq. (23) without the added noise as a function of µ/τ . The numerical results of maximum bias b ǫ ( )and b ǫ ( ) are produced by Monte Carlo simulation (see Appendix G-B). (b) The optimal noise level a opt corresponding to (a) versus µ/τ . Here, the fraction of contamination ǫ =0.1.

  (w + d) -ρ(w -d)] 2 , (45) with A =2E w [ψ(w + d) -ψ(w -d)]E w [ψ ′ (w + d)+ψ ′ (wd)]E w [ρ(w + d) -ρ(w -d)] 2 -2E 2 w [ψ (w + d) -ψ (wd)]E w [ρ(w + d) -ρ(w -d)][ψ(w + d)+ψ(w -d)].T h e optimal noise level d opt is solved by ∂Eff ∞ ( θ, d)/∂d =0. Since J(f w ) > 0 and E 2

APPENDIX G NUMERICAL SIMULATION

A. Numerical Computation of the Asymptotic Efficiency

We adopt the iterative reweighting method to obtain the numerical estimation θ of a location parameter [START_REF] Maronna | Robust Statistics: Theory and Methods[END_REF]. Assuming

then Eq. ( 3) can be rewritten as [START_REF] Maronna | Robust Statistics: Theory and Methods[END_REF]

From Eq. ( 48), the estimated parameter can be obtained as a weighted mean [START_REF] Maronna | Robust Statistics: Theory and Methods[END_REF]

with

For the method of adding extra noise to an array of M-estimators, we first generate L groups of updated samples with the addition of L samples of uniform noise components η li with the level of d. Then, compute L values of parallel estimators ( θ1 , ••• , θL ) by using L groups of observations respectively. For each group of observations y l =(y l1 ,y l2 , ••• ,y ln ),thelth estimator θl is obtained by the iterative reweighting method as follows:

1) Set the initial estimator θl (0) as the median of the updated observations y l .

2) For k =0, 1, 2,..., compute the location parameter θl (k + 1) at iteration k +1as

with the weight function

where ζ is a small tolerance parameter, and the numerical M-estimator θl (k +1) isassumedtobethelth expected M-estimator θl .

4) Compute the estimator θ of Eq. ( 15) by averaging

Then, we can realize the numerical M-estimator θ for M =10 4 Monte Carlo trails by the above iterative reweighting method, and the variance of θ is evaluated as var( θ)= M m =1 ( θmθ) 2 /(M -1) with the mean θ = M m =1 θm /M . Substituting the variance var( θ) into Eq. ( 7), the numerical asymptotic efficiency Eff a ( θ) can be obtained, as shown in Fig. 3(c).

B. Numerical Simulation of the Maximum Bias

We use 10 4 Monte Carlo trails to evaluate the maximum bias with optimal noise and without added noise, respectively. For each trail, the observation size n =5× 10 4 and the fraction of contamination ǫ =0.1. We generate n observations obeying the Gaussian mixture distribution, then replacing 0.1n of observations with the outliers, resulting in the contaminated data x i . According to the contamination model of Eq. ( 20), the maximum bias is theoretically evaluated with the extreme outliers obeying the point mass distribution q = δ(x 0 ) at x 0 →∞. In practice, the outliers are taken as real numbers with sufficiently large magnitude as 10 6 ≫ µ +3τ . Then, we generate n random components η i composed by two values of a opt and -a opt with equal probabilities. Then, the updated observations y i = x i + η i are obtained. Next, calculate the median of the samples by the original data x i and the updated observations y i , respectively. Finally, the maximum bias with optimal noise and without added noise are obtained for 10 4 numerical trials, respectively, as shown in Fig. 6.

APPENDIX H NOISE-REDUCED MAXIMUM BIAS OF THE MEDIAN ESTIMATOR

For the median estimator and a given proportion 0 <ǫ<1/2, we expect that the updated bias b ǫ of Eq. ( 24) is smaller than the bias b ǫ of Eq. ( 23), i.e.,

Here, the median estimator has the breakdown point of ǫ * =1/2 that yields an unbounded estimator θ →∞, thus requiring 0 < ǫ<1/2. Since the probability density f w or f z is symmetric, the cumulative distribution function is then analyzed for the nonnegative domain of x ≥ 0. It is known that the cumulative distribution function is nondecreasing, then finding a smaller updated bias of b ǫ <b ǫ is equivalent to satisfying the condition of F z (b ǫ ) >F w (b ǫ ) for b ǫ = b ǫ . The cumulative distribution function F z of the composite noise z can be calculated as

For a given positive value of b ǫ ∈ (U, V ) for 0 ≤ U<Vand a function F w (x) convex over the interval (U, V ), then, according to the Jensen inequality [START_REF] Cover | Elements of Information Theory[END_REF], we obtain

with the zero mean E η [η]=0. This indicates that if there exists a positive interval (U, V ) where the cumulative distribution function F w is strictly convex, then the maximum bias of the median estimator can be reduced by adding extra noise η.

Under the above convexity condition of F w over a certain positive interval (U, V ), a naturally emerging problem is to find an optimal noise distribution f opt η , aiming to minimize the maximum bias of Eq. ( 24) without changing the consistency of the estimator under f z , i.e., E z [ψ(z)] = 0 with ψ(x)=sgn(x). It is noted that minimizing the maximum bias is equivalent to maximizing F z (b ǫ ) for a given positive bias b ǫ ∈ (U, V ). Thus,


Fabio Fichera

Andrea Mariscotti
email: andrea.mariscotti@unige.it

Ieee Member

IEEE Senior Member Ade Ogunsola
email: ade.ogunsola@parsons.com

Evaluating Stray Current from DC Electrified Transit Systems with Lumped Parameter and Multi-Layer Soil Models 1

Keywords: Corrosion, Electromagnetic conductive interference, Guideway transportation power systems, Power distribution electromagnetic interference, Stray current

DC electrified traction systems are a potential source of stray current, that may corrode both internal and external structures and installations. The stray current intensity depends on several factors (traction current, rail insulation, concrete mat and piers resistance), that are accounted for in the proposed model. Moreover, the effect of a uniform or multi-layer soil model is considered, showing the differences and the impact on the countermeasures with respect to a uniform soil model simplifying approach.

I. INTRODUCTION

In DC electrified railway, the running rails are used for traction current return and form part of the signaling circuit for the control of train movements. However it is necessary to consider the consequences that this solution involves: as the running rails are isolated from ground, current will leak to ground, depending on the rail-to-earth conductance term for each running rail, with the ground now forming an unintended current return path.

In practice, limitation of the current leaking to ground and thus causing corrosion is implemented by decreasing the rail resistance and increasing the rail to earth resistance -achieved by adequate insulation and coating. A lower running rail resistance has the advantage of increasing the proportion of current propagating through the rails, while a greater rail to earth resistance has the advantage of limiting the amount of leakage current.

The earthing scheme implemented on a dc electrified railway largely influences the current distribution and thus the stray current [START_REF] Jamali | Effects of different earthing schemes on the stray current in rail transit systems[END_REF]- [START_REF] Lee | Assessment of grounding schemes on rail potential and stray currents in a DC transit system[END_REF]. When a metallic structure is electrically influenced by stray current, the potential of the structure shifts in the positive or in the negative direction where the current leaves or enters the metal surface. The positive potential deviation of the structure is a measurable effect which is used as electrolytic corrosion hazard criterion in standards and regulations of many countries. However such measurements are usually post construction and thus may lead to an expensive retrofit program, in the case of undesirable results.

Determining the response of the structure is not an easy task, because it is a function of three things: the location of the structure with respect to the electric flow field generated by the stray current source (protection system), the magnitude of the electric field, and the electrochemical response of the structure to the interference (protection).

In the electric circuit approach, based on the earth-return circuit theory, which is more suitable for modeling extended structures (pipelines, cables, etc.), the polarization phenomenon is not taken into account. Perhaps the first publication that provides a detailed analytical approach for the assessment of stray current is [START_REF] Mccollum | Leakage of currents from electric railways[END_REF]. The formulation demonstrates the importance of keeping the leakage factor, the square root of the value of the rail resistance per unit length divided by the track-to-earth leakage resistance per unit length, as low as possible, and provides an analytical solution for the estimation of stray current for grounded and ungrounded dc electrified railways. It however does not consider the presence of a stray current collector mat and thus is only applicable to primary current leakages from rail to ground and not to modern metro and rapid transit systems.

Due to the limitation expressed above, analytical models have been developed by various researchers [START_REF] Jamali | Effects of different earthing schemes on the stray current in rail transit systems[END_REF][5][7]- [START_REF] Charalambous | Influence of soil structures on corrosion performance of floating-DC transit systems[END_REF] for the estimation of stray current from dc-electrified railways. Explicit in these models is the assumption that the dc electrified railway is directly laid above ground; furthermore the analytical solutions gives no indication on the potential impressed on a buried metallic object nor an estimate of the corrosion current on a buried pipeline.

In [START_REF] White | Earthing the London underground[END_REF] a distributed parameter model is presented, where the resistance of the various interposed structural elements (supporting element and stray current concrete mat, for example), as well as the grounding resistance of the viaduct piers, may be included in the mesh of resistive elements, that is part of the model.

II. OVERVIEW OF THE MODEL STRUCTURE

The considered model is a generic model describing a dcelectrified railway on a viaduct consisting of two Traction Power Stations (sourcing the currents I 1 and I 2 , respectively) and one train, at a generic position x between 0 and L and absorbing a constant current I [START_REF] Ogunsola | Estimation of Stray Current from a dc Electrified Railway and Impressed Potential on a Buried Pipe[END_REF]. The distributed parameter model in [START_REF] Ogunsola | Estimation of Stray Current from a dc Electrified Railway and Impressed Potential on a Buried Pipe[END_REF] for the same system represents the reference for the validation of the lumped parameter model developed successively. A lumped parameter model, following a circuit approach, seems a less elegant and a less performing solution to the same problem. Yet, the reason for investigating it resides in the possibility of introducing easily discontinuities and modeling different and complicated situations, for example introducing local defects in the rail fasteners, the sleepers, the ballast, the soil, buried metal structures and the grounding system. The overall accuracy is not sacrificed: stray current phenomena are low frequency, or dc in the present case, phenomena and do not require a fine grain representation in the longitudinal axis, thus keeping the spatial step size at a convenient value. For the results shown in the following 10 m were used, after a validation of the chosen value by means of a sensitivity analysis.

The lumped parameter model proceeds writing down the Kirchhoff equations and the constitutive equations for the series and shunt circuit elements; then the "sparse tableau method" is applied, obtaining a circuit representation as follows:

[ ] [ ] [ ] [ ] U I Z V Y V A V I A n T = + = - = 0 0 (1)
where I and V represent the branch current and voltage vectors for all the elements of the circuit model, V n is the vector of the node voltages, U is the vector of the external variables (either currents or voltages), [Y] and [Z] are admittance and impedance matrices and [A] is an incidence matrix implementing the basic Kirchhoff equations.

This method leads to a square matrix that can be inverted, of size 2R+N-1, where R is the number of branches and N is the number of the nodes of the circuit:

[ ] [ ] [ ] [ ] [ ]                   - = 0 0 0 0 1 A Z Y A T T (2)
The solution of the circuit model is automatic by inversion of T. The model has a graphical user interface to ease the data input for simple cases and to display the main results, as shown in Fig. 2.

III. MODEL VERIFICATION AND SIMULATION RESULTS

The verification has proceeded through initial simple cases, featuring only one layer of circuit elements, simulating the case of rails directly laid on a uniform ground (that includes the effects of the track ballast and soil), with no buried structures. Then, a stray current mat and later a buried pipe are added. Last, the uniform soil model is changed to a multilayer soil model.

A. Simple case of track on uniform soil, no buried structures

This configuration was considered to perform a sensitivity analysis and a direct comparison with the expressions and values found in [1][2], where a similar configuration is assumed. Assuming the rails fastened to sleepers, sunk in the ballast and soil, with the rails electrically insulated from the sleepers themselves by an insulating pad, the following conductance terms can be identified: G rp for the rail to pad term, G ps for the pad to sleeper term, G sb for the sleeper to ballast term, and G be for the ballast to earth term. For each term a minimum and maximum value, based on sensible assumptions and practical knowledge, has been assumed as shown in Table I. The resulting min and max values for the overall rail to earth conductance G re (given by the series combination of the four terms above) are about 0.001 and 1.0 S/km. It can be noted that in the best track conditions (low conductance values for the four terms) G re is strongly influenced by G rp , whilst in the worst case of high conductance to earth G re is limited by the ballast conductance value G be .

With reference also to Fig. 1, the other parameter that is needed for this simplified configuration is the rail longitudinal resistance R r that is normally in the range of 0.03 -0.08 Ω/km, the track resistance R t for the two rails in parallel being half of it.

For a 1000 A absorbed train current and a TPS separation L = 3 km, it is possible to calculate and plot the stray current per unit length as a functions of G re and R r (see Fig. 3). The calculated value is below the limit of 2.5 mA/m [2] if it lies in the portion below the horizontal arrow, that is for the maximum G re value if R t is below 0.01 Ω/km and for the maximum R t value if G re is below 0.19 S/km.

It is interesting to observe that the rail to earth voltage U re can be higher than the limit of 5 V prescribed in [1] for a G re > 0.5 S/km, without exceeding the limit on the stray current. This means that the two formulations for rail voltage and stray current given in [1][2] are not consistent for all parameter values, but bear an unexpressed safety margin. The expression for stray current reported in the EN 50122-2 standard [2] has also been considered and found in agreement with the results shown in Fig. 3 with a maximum error of 16%.

The longitudinal train position that produces the largest stray current is that in the middle; positions closer to either TPS make a larger current flow in the rails of one of the two sections, but over a shorter length, thus not leading to the maximum stray current intensity.

B. Inclusion of the stray current concrete mat

This configuration corresponds to the one shown in Fig. 1 and considered in [START_REF] Ogunsola | Estimation of Stray Current from a dc Electrified Railway and Impressed Potential on a Buried Pipe[END_REF]. The resistance of the concrete mat R cm may be estimated by different methods, but the one proposed in [START_REF] Thapar | Ground resistance of concrete foundations in substation yards[END_REF] is followed.

      +       - ρ ρ ρ = 1 1 4 r r r R c c c cm (3)
where r c is the external radius of the concrete and r is the internal radius of the iron mass, estimated as √(π/A), A being the surface.

The concrete resistivity range can be very wide, depending on environment conditions of the trackbed, and normally it is assumed in the range 50-250 Ωm. The influence of the concrete resistivity is shown in Fig. 6 where the calculations are made for a soil resistivity of 100 Ωm, R r = 0.02 Ω/km, G r,cm = 0.5 S/km and 1000 A of absorbed train current.

C. Efficiency of stray current mat and impact on buried pipe

A buried metallic pipe is added (as shown in Fig. 7) and the efficiency of the stray current mat system is verified by simulation. Under the assumption of an ideal connection of the stray current mat to the negative terminal of the TPS by means of a diode (or a rail potential control device) that helps having cathodic protection and making the current flow unidirectional, the impact on an external buried metallic pipe is considered. The resistance to earth of the pipe has been estimated using the Dwight's expression [START_REF] Ogunsola | Estimation of Stray Current from a dc Electrified Railway and Impressed Potential on a Buried Pipe[END_REF] [START_REF] Dwight | Calculation of resistance to ground[END_REF].

For brevity the efficiency E of the stray current mat is introduced, The uniform soil model, featuring only one layer of soil with constant resistivity and infinite depth, may lead in some cases to an overestimation of the pipe impressed potential, thus requiring unnecessary countermeasures. It was found that in particular soil resistivity configurations, where different layers can be distinguished, a multi-layer soil model is needed [START_REF] Jabbehdari | Sensitivity Analysis of Multilayer Earth Models[END_REF]. For a top soil layer with lower resistivity (that is the most common case) the results show a dramatic reduction of the calculated pipe potential (as shown in Fig. 9). If, on the contrary, the top layer has higher resistivity than the bottom one, the stray current mat efficiency decreases and the pipe potential increases (as shown in Fig. 10). The two cases above (shown in Fig. 9 and 10) were using purposely extreme values of the resistivity values of the two soil layers (10 and 1000 Ωm, exchanged in the two considered cases). A comparison between the uniform soil model and the two layer soil model for the same measured apparent resistivity values (taken from one reference case), gives the following results (see Fig. 11). The difference between the values predicted by either of the two models is remarkable, has relevance for the required countermeasures and thus poses the problem of the best modeling approach when assessing the stray current phenomenon.

IV. CONCLUSION

The problem of stray current evaluation in a DC transit system by means of simulation has been considered and solved with a lumped parameter model, compared for validation with a previous distributed parameter model [START_REF] Ogunsola | Estimation of Stray Current from a dc Electrified Railway and Impressed Potential on a Buried Pipe[END_REF]. This model has the possibility of including any discontinuity or simulating local defects or deviations of some parameters from the nominal or general values. Moreover, this model has been use in conjunction with both a uniform soil model and a multi-layer soil model, demonstrating the degree of approximation of the former and the reduction of the undue margins for countermeasures with the latter.

F

  . Fichera, A. Mariscotti, and A. Ogunsola, "Evaluating Stray Current from DC Electrified Transit Systems with Lumped Parameter and Multi-Layer Soil Models", IEEE EUROCON 2013, July 1-4, 2013, Zagreb, Croatia. doi: 10.1109/EUROCON.2013.6625131

Fig. 1 .

 1 Fig.1. Track scheme including the stray current mat and the viaduct piers[START_REF] Ogunsola | Estimation of Stray Current from a dc Electrified Railway and Impressed Potential on a Buried Pipe[END_REF]

Fig. 2 .

 2 Fig. 2. Graphical user interface of the stray current calculator

  /km] Gre [S/km] max. stray current per unit length [mA/m]

Fig. 3 .

 3 Fig. 3. Stray current per unit length as a functions of G re and R r

Fig. 4 .

 4 Fig. 4. Rail to earth voltage U re for R t = 0.01 Ω/km and G re = 0.8 S/km

Fig. 5 .

 5 Fig. 5. Stray current as a function of train longitudinal position

Fig. 6 .

 6 Fig. 6. Current in the (a) concrete mat and (b) soil as a function of concrete resistivity

Fig. 7 .

 7 Fig. 7. Stray current mat and buried metallic pipe

Fig. 8 .

 8 Fig. 8. Pipe potential for the uniform soil model and two stray current efficiency values: (a) 50%, (b) 90%

Fig. 9 .

 9 Fig. 9. Pipe potential for the two layer soil model (top layer with lowest resistivity)

Fig. 10 .

 10 Fig. 10. Pipe potential for the two layer soil model (top layer with highest resistivity)

Fig. 11 .

 11 Fig. 11. Pipe potential comparing the uniform (1) and two layer soil models (2)