Pavlina Jordanova

Zdeněk Fabián

Philipp Hermann

Lubos Strelec

Andrés Rivera

Stéphane Girard

Sebastián Torres

Milan Stehlík
email: mlnstehlik@gmail.com

Luboš Střelec

Weak properties and robustness of t-Hill estimators

Keywords: Point estimation, Asymptotic properties of estimators, t-Hill estimator, t-lgHill estimator Mathematics Subject Classification (2000) 62F10, 62F12

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The aim of this paper is to introduce a novel method of heavy tails estimation based on transformed score. Especially for small or/and contaminated samples such a method can be better than standard maximum likelihood, since it relates to the method of moments for transformed scores (see [START_REF] Fabián | Score function of distribution and revival of the moment method[END_REF]; [START_REF] Hosking | Parameter and quantile estimation for the generalized pareto distribution[END_REF]). Let us denote RV a the class of regularly varying functions at infinity, with an index of regular variation equal to a ∈ R, i.e. positive measurable functions g(•) such that for all x > 0, g(tx)/g(t) → x a , as t → ∞.

Here we suppose that X 1 , X 2 , . . . , X n are possibly dependent copies of X with distribution function (d.f.) F, upper order statistics

X (1,n) ≤ X (2,n) ≤ • • • ≤ X (n,n) and 1 -F ∈ RV -α , α > 0. (1)
We define α to be the tail parameter and γ := 1/α to be the extreme value index. [START_REF] Hill | A simple general approach to inference about the tail of a distribution[END_REF] derived a procedure of Pareto tail estimation by MLE, obtaining the following Hill estimator

H k n ,n = 1 k n k n ∑ i=1 log { X (n-i+1,n) X (n-k n ,n) } , k n = 1, 2, . . . , n -1. ( 2
)
Later on, many authors tried to robustify the Hill estimator, but they still relied on maximum likelihood method. [START_REF] Alves | A location invariant hill-type estimator[END_REF] has introduced a new lower bound and [START_REF] Gomes | Censoring estimators of a positive tail index[END_REF] respectively [START_REF] Li | Asymptotic normality of location invariant heavy tail index estimator[END_REF] have introduced powers of original statistics. However, the influence function of Hill estimator is slowly increasing but unbounded, thus, the Hill procedure is not robust. Further approaches of robustifying the original Hill estimator were given in [START_REF] Beran | On robust tail index estimation[END_REF]; Vandewalle et al. (2007). In [START_REF] Fabián | Induced cores and their use in robust parametric estimation[END_REF] a new score method of score moment estimators has been proposed. It appeared that these score moment estimators are robust for very heavy tailed distributions, see Stehlík et al. (2010a). [START_REF] Jordanova | Weak asymptotic results for t-hill estimator[END_REF] consider an independent identically distributed (i.i.d.) sample and find the limit distribution of the t-Hill estimator for a fixed number k of the threshold order statistics. They prove that a sample of Pareto distributed observations does not have to be large in order to receive the corresponding limit distribution for fixed k. In that case under suitable normalizations and a large sample the t-Hill estimator is asymptotically normal for k(n) → ∞. Under the more general conditions, the t-Hill estimator is asymptotically normal for k(n) = o(n). The Hill estimator procedure with the score moment estimator has been investigated in Stehlík et al. (2012) for optimal testing for normality against Pareto tail. Recently, nice generalizations of t-Hill have been published, see [START_REF] Billingsley | A simple generalization of the hill estimator[END_REF]; [START_REF] Paulauskas | On the improvement of hill and some other estimators[END_REF]; [START_REF] Beran | The harmonic moment tail index estimator: asymptotic distribution and robustness[END_REF]. Resnick and Starica (1993) generalize the Hill estimator for more general settings with possibly dependent data. In this paper we continue these investigations, since dependencies are expected in real data applications. We obtain weak consistency of the t-Hill estimator for a special class of dependent data, the infinite moving average model. Moreover, we provide some examples showing that in contrast to the i.i.d. case the t-Hill and the Hill estimators applied to the moving average model are not robust with respect to large observations. Under the concept of the Hill estimator we understand the successive averaging of ordered values up to given k. In this paper we understand "The Hill estimator" as a specific procedure for studying the tail of Pareto like distributions. Instead of implementing "The Hill estimator" procedure, we implement the t-score moment procedure, which can give different estimators for different families. We illustrate t-Hill and we also quantify the robustness and compare efficiency with other competitors.

The paper is organized as follows: in the next section we recall the theory of score function of the distribution introduced firstly by [START_REF] Fabián | Induced cores and their use in robust parametric estimation[END_REF] and studied later in series of works. In section 3 we provide asymptotical results (asymptotic consistency and asymptotic normality) for t-Hill estimators for Pareto and log-gamma distributions, under iid and moving average settings. In section 4 we compare maximum likelihood, t-Hill, transformed score moment and Hill estimators.We also introduce t-lgHill estimator, which is a transformed score moment estimator for case of loggamma distribution. Contamination of underlying data is controlled by the extreme value index, shifting parameter or the means of transformed score variance of contaminating Pareto distribution. Comparisons show that in cases when the contamination has heavier tail than the original distribution, t-Hill estimator outperforms Hill estimator and several robust tail estimators (Integrated Squared Error Estimator (ISE), Partial Density Component (PDC) estimator (Vandewalle et al. (2007)), Least Squares estimator (LS), moment estimator (ME, see Dekkers et al. (1989)), QQ estimator [START_REF] Kratz | The qq-estimator and heavy tails[END_REF]) or Weighted MLE estimator (WMLE, see [START_REF] Dupuis | Robust weighted likelihood estimators with an application to bivariate extreme value problems[END_REF]; [START_REF] Dupuis | A robust prediction error criterion for pareto modelling of upper tails[END_REF]). We investigate robustness of t-Hill estimator for a particular case of moving average sequence. In section 5 we provide an application of t-Hill in comparison to Hill and Zipf's estimator for a small sample data set from Guanaco Glacier. Summary concludes the paper. For reasons of readability we provide all proofs in the Appendix.

2 Transformed Score Function and Transformed Score Moment Estimators [START_REF] Pearson | of frequency constants and on the influence of random se[END_REF], Edgeworth (1908a), Edgeworth (1908b), and[START_REF] Fisher | Statistical Methods For Research Workers. Cosmo study guides[END_REF] developed the basic statistical inference function, the so called score function. Their score is the gradient with respect to parameter θ of the logarithm of the likelihood function S(θ , X) = ∂ ∂ θ log L(θ ; X), indicating the sensitivity of L (its derivative normalized by its value). The variance of the score is the Fisher information

I (θ ) = E θ (S 2 ).
Fabián (2001,2007,2010,2015) has introduced a more general scalar-valued inference function, in this paper called transformed score function (or t-score) of distribution, reflecting main features of continuous probability distributions and enabling an introduction of their new relevant characteristics.

The t-score of distribution G with support R and unimodal differentiable density g is a function

S G (y) = - g ′ (y) g(y) (3)
expressing the relative rate of the change of g. For a distribution G µ with location parameter µ ∈ R, that is, with density in the form g(x -µ), the generalized t-score function equals the Fisher score for µ as

S G (y -µ) = ∂ ∂ µ log g(y -µ).
We consider that t-score in equation ( 3) is a significant function of parametric distribution G θ , θ ∈ Θ ⊆ R, and the mode y * : S G (y; θ ) = 0 its 'central' point.

Example 1 The t-score of N (µ, σ ) is S N (y) = y-µ σ 2 and y * = µ.

In general, the t-score of distribution F with arbitrary interval support X ⊆ R is defined by the following construction.

Definition 1 Let η : X → R be a strictly increasing smooth mapping. The t-score of distribution F with interval support X ⊆ R and (almost surely) two-times differentiable density f (x) is

T F (x) = - 1 f (x) d dx [ 1 η ′ (x) f (x) ] . ( 4
)
Supposing that the solution x * of equation

T F (x) = 0 (5)
is unique, it is called the t-score mean and function

S F (x) = η ′ (x * )T F (x) (6)
the generalized t-score function of distribution F.

Definition 1 can be explained as follows:

F is a transformed distribution F(x) = G(η(x)), where G(y) is a 'prototype' with support R. F has density f (x) = g(η(x))η ′ (x) (7)
where η ′ (x) = dη(x)/dx is the Jacobian of the transformation. It was shown in [START_REF] Fabián | Induced cores and their use in robust parametric estimation[END_REF] that the t-score is actually the transformed generalized t-score function of the prototype, T F (x) = S G (η(x)). The Jacobian of the transformation does not carry any information about the distribution. The term in brackets of equation ( 4) is the density without the Jacobian. The justification of expression ( 6) is the following: consider a location prototype G µ with density g(y -µ). The density of the transformed distribution

F τ = G µ • η(x) is f (x; τ) = g(η(x) -η(τ))η ′ (x),
where τ = η -1 (µ) is a so called transformed location parameter. It was proven in [START_REF] Fabián | Estimation of simple characteristics of samples from skewed and heavy-tailed distribution[END_REF] that the generalized t-score function with transformed location parameter is identical with the Fisher score for this parameter,

S F (x; τ) = d dx log f (x; τ).
Generalized t-score function ( 6) and its parametric form S F (x; θ ) generalize the concept of the Fisher score for distributions with arbitrary support without a 'central' parameter or without parameters at all. If the prototype G has mode y * , the score mean of the transformed distribution is the transformed mode of the prototype G, x * = η -1 (y * ). The generalized t-score function is actually the (generalized) Fisher score for the score mean, which may not be a parameter of the distribution and is unique if the prototype is unimodal.

The mapping η(x) occurring in Definition 1 is often the inner part of f (x), so that η(x) and/or η ′ (x) in equation ( 7) are clearly identifiable. This η is called an innate mapping; [START_REF] Fabián | Score function of distribution and revival of the moment method[END_REF] discusses how to choose a most suitable innate mapping if it is not apparent from the density formula.

Let us consider two examples, related to log-gamma and Pareto distributions studied later in this paper. For other examples, e.g. beta-prime distribution see [START_REF] Fabián | Score function of distribution and revival of the moment method[END_REF].

Example 2 The log-gamma distribution with support

X = (1, ∞) has density f (x) = α c Γ (c) (log x) c-1 1 x α+1 = α c Γ (c) (log x) c x -α 1 x log x (8) so that η ′ (x) = 1
x log x and the innate mapping is η(x) = log log x. By (4) one obtains

T F (x) = - 1 f (x) d dx [(log x) α x -α ] = α log x -c, with score mean x * = e c/α . Example 3 The density f (x) = α x α+1 (9)
of the Pareto distribution with X = (1, ∞) does not contain any 'visible' Jacobian of any transformation. By using η ′ (x) from the foregoing example we would receive

T F (x) = - 1 f (x) d dx [x log x f (x)] = α log x -1,
which is the Fisher score for α. By using η(x) = log(x -1), η ′ (x) = 1/(x -1), the t-score ( 4) is

T F (x) = - 1 f (x) d dx [(x -1) f (x)] = α -(α + 1)/x = α(1 -x * /x), ( 10
)
with score mean x * = (α +1)/α. Since the 'shifted' Pareto with X = R + and density f (x) = α (x+1) α+1 is the particular case of the beta type II distribution with bounded generalized t-score, the latter one is to be taken as the t-score of the Pareto distribution.

The score moments (SM) were introduced for any k ∈ N instead of the ordinary moments by

M k (θ ) = ES k F (θ ) = ∫ X S F (x; θ ) k f (x; θ ) dx, ( 11
)
existing only if f satisfies the usual regularity requirements. It appeared that the score moments are often expressed by elementary functions of parameters. It is easy to see that M 1 = 0. The value M 2 = ES 2 F of distributions with (transformed) location parameter is the Fisher information for (transformed) location parameter. Accordingly, ES 2 F is the Fisher information for the t-mean. The reciprocal value

ω 2 = 1 ES 2 F , ( 12
)
the t-score variance, appeared to be a measure of the variability (dispersion) of the distribution even in cases where the usual variance does not exist (see [START_REF] Fabián | Estimation of simple characteristics of samples from skewed and heavy-tailed distribution[END_REF]). Notice, that in case of the Pareto distribution

ET 2 F = α/(α + 2), ES 2 F = (x * - 1) -2 ET 2 F so that ω 2 F = α + 2 α 3 . ( 13
)
Let X 1 , . . . , X n be i.i.d. random variables according to some F. Assuming F as a member of the model family [START_REF] Fabián | Induced cores and their use in robust parametric estimation[END_REF][START_REF] Fabián | Score moment estimators[END_REF] introduced the t-score moment estimate θSM as the solution of equations θSM :

{F θ , θ ∈ Θ }, Θ ∈ R m ,
1

n n ∑ i=1 S k F (x i ; θ ) = ES F k (θ ), k = 1, . . . , m, (14)
the statistical counterpart of (11

∑ i=1 (α log x i -c) = 0 (15) 1 n n ∑ i=1 (α log x i -c) 2 = c ( 16
)
with solution α = s 1 s 2 -s 2 1 , c = s 2 1 s 2 -s 2 1 where s 1 = 1 n ∑ n i=1 log x i and s 2 = 1 n ∑ n i=1 (log x i ) 2 .
Example 5 Pareto distribution: From the first moment equation

∑ n i=1 (1 -x * /x i ) = 0 we obtain x * = xH , where xH = n/ n ∑ i=1 1/x i
is the harmonic mean, and

α = 1/( xH -1). ( 17
)
We use the same notations of order statistics as in section 1, see (1). Hence, due to (17), the t-Hill estimator of γ = α -1 was suggested in the form

H * k n ,n = { 1 k n k n ∑ i=1 X (n-k n ,n) X (n-i+1,n) } -1 -1, k n = 1, 2, ..., n -1. ( 18
)
This estimator was firstly published in technical report of [START_REF] Fabián | On robust and distribution sensitive hill like method[END_REF].

Since it is based on harmonic mean, which is a generalized t-Hill estimator, it is expected to be resistant to large observations to a certain extent so that it could yield more realistic values than the ordinary Hill estimator. Further generalizations have been made by [START_REF] Billingsley | A simple generalization of the hill estimator[END_REF]; [START_REF] Beran | The harmonic moment tail index estimator: asymptotic distribution and robustness[END_REF]. In particular, the tradeoff between efficiency and robustness has been studied in [START_REF] Beran | The harmonic moment tail index estimator: asymptotic distribution and robustness[END_REF], a mean of order p ≥ 0 (MOP) generalization is given by [START_REF] Billingsley | A simple generalization of the hill estimator[END_REF].

For the definition of the t-lgHill estimator we use the same notations of order statistics as in section 1, see (1), and

M ( j) n = 1 k n k n ∑ i=1 ( log X (n-i+1,n) X (n-k n ,n) ) j , j = 1, 2, ...
Hence the t-lgHill estimator of γ = α -1 has the form

H L k n ,n = M (2) n -(M (1) n ) 2 M (1) n . ( 19
)
If we understand the t-Hill estimator as an algorithm depending on the assumed distribution, the t-Hill for log-gamma distribution is, according to the solution of ( 15) and ( 16)

H L k n ,n = ∑ k n i=1 log 2 { X (n-i+1,n) X (n-kn,n) } ∑ k n i=1 log { X (n-i+1,n) X (n-kn,n) } - 1 k n k n ∑ i=1 log { X (n-i+1,n) X (n-k n ,n) } k n = 1, 2, ..., n -1. ( 20
)
Let us denote H L k n ,n sequence by t-lgHill. Referring to Section 2, the t-lgHill estimate of the tail index is given by ( 16) and has a closed-form expression as

H L k n ,n = s 2 s 1 -s 1 .
3 Asymptotical results

t-Hill Estimator in Case of Moving Average Sequence

Suppose at least one of the real numbers c j , j = 0, 1, ... is positive and there exists δ ∈ (0, 1), δ < α such that

∞ ∑ j=0 |c j | δ < ∞. ( 21
)
Consider the moving average sequence

X n = ∞ ∑ j=0 c j Z n-j , -∞ < n < ∞, (22)
where Z i , -∞ < i < ∞, are non-negative i.i.d. innovations with d.f. G, such that G ∈ RV -α , α > 0. We will use the t-Hill estimator as introduced before and consider the point measure

µ X,k n ,n (•) := 1 k n n ∑ i=1 ε { X i b( n k n ) ∈ (•) } ,
as a random element in the space E + of positive Radon measures on (0, ∞] endowed with the vague topology. Here

b(t) := F ← (1 - 1 t ) = ( 1 F ) ← (t) (23)
is the tail function of F. Let

k n → ∞, k n n → 0. ( 24
)
By Proposition 3.3 of Resnick and Starica (1993)

1 k n n ∑ i=1 I [ X i b( n kn ) >0 ] ε { X i b( n k n ) ∈ (•) } ⇒µ, n → ∞, in E + , where µ : σ ((0, ∞]) → [0, ∞) and µ(x, ∞] = x -α , x > 0. However, for all n ∈ N, k n = 1, 2, ..., n and A ∈ E + 1 k n n ∑ i=1 I [ X i b( n kn ) >0 ] ε { X i b( n k n ) ∈ A } =µ X,k n ,n (A), (25) in distribution, therefore µ X,k n ,n ⇒µ, n → ∞, (26)
in E + . By Proposition 2.1. in Resnick and Starica (1993), having an intermediate sequence,

X (n-k n ,n) b ( n k n ) P →1 ( 27
)
and by Proposition 2.2. in Resnick and Starica (1993)

µ X,k n ,n (•) := 1 k n n ∑ i=1 ε { X i X (n-k n ,n) ∈ (•) } P →µ(•), n → ∞, (28)
in E + . We will use the following Potter's inequality for distribution functions with regularly varying tails. If a function H is regularly varying with exponent γ ∈ R, then for ε > 0 there will exist t 0 (ε), such that for t > t 0 (ε) and x ≥ 1,

(1 -ε)x γ-ε ≤ H(tx) H(t) ≤ (1 + ε)x γ+ε . ( 29
)
Its proof can be found in de [START_REF] De Haan | On regular variation and its application to the weak convergence of sample extremes[END_REF].

Lemma 1 Let X 1 , X 2 , ... be random variables with d.f. F such that

F(x) ∈ RV -α . ( 30
) Assume that k n → ∞, k n /n → 0 and µ X,k n ,n ⇒ µ as n → ∞. Define Φ n = 1 k n k n ∑ i=1 φ ( X (n-i+1,n) X (n-k n ,n) ) , where φ is a differentiable function on (1, +∞) such that x -α φ(x) → ∞ as x → ∞ and ∫ ∞ t x -α+δ 1 |φ ′ (x)| dx < +∞, for all t > 1, 0 < δ 1 < α. Then, Φ n is a weakly consistent estimator of φ(1) + ∫ +∞ 1 x -α φ ′ (x)dx.
As a consequence of Lemma 1, considering φ(x) = 1/x yields that H * k n ,n , defined in ( 18) is a weakly consistent estimator of 1/α. Now we are ready to prove the weak consistency of the t-Hill estimator in case of infinite moving average (MA) sequence.

Proposition 1 Suppose at least one of the real numbers c j , j = 0, 1, ... is positive and there exists δ ∈ (0, 1), δ < α such that condition ( 21) is satisfied. Consider the moving average sequence ( 22) where Z 24), then H * k,n is a weakly consistent estimator for α -1 .

i , -∞ < i < ∞, are non-negative i.i.d. innova- tions with d.f. G, such that G ∈ RV -α , α > 0. If (

t-lgHill Estimator in Case of Moving Average Sequence

Suppose at least one of the real numbers c j , j = 0, 1, ... is positive and there exists δ ∈ (0, 1), δ < α such that ( 21) is satisfied.

Consider the moving average sequence ( 22), where

Z i , -∞ < i < ∞, are non- negative i.i.d. innovations with d.f. G, such that G ∈ RV -α , c = α > 0.
Assume (24) holds. By weak consistency of Hill estimator we obtain that

M (1) n P → 1 α . Now we need to prove that M (2) n -(M (1) n ) 2 P → 1 α 2 , i.e. M (2) n P → 2 α 2 (31)
and then to use Slutsky arguments. Consider the point measure

µ X,k n ,n (•) := 1 k n n ∑ i=1 ε { X i b( n k n ) ∈ (•) } ,
as a random element in the space E + of positive Radon measures on (0, ∞] endowed with the vague topology, where b is defined in (23).

In section 3.1 we proved that under these conditions

µ X,k n ,n ⇒µ, n → ∞, (32)
in E + and ( 26) and ( 27) are true. As a consequence of Lemma 1, letting φ(x) = ln(x) entails that M

(2)

n is a weakly consistent estimator for 2/α 2 .

Proposition 2 Suppose at least one of the real numbers c j , j = 0, 1, ... is positive and there exists δ ∈ (0, 1), δ < α such that condition ( 21) is satisfied. Consider the moving average sequence ( 22

) where Z i , -∞ < i < ∞, are non-negative i.i.d. innovations with d.f. G, such that G ∈ RV -α , α > 0. If (24), then i.) M (2)
n is a weakly consistent estimator for 2α -2 ; ii.) t-lgHill estimator H L k,n is a weakly consistent estimator for α -1 .

t-lgHill Estimator is Asymptotically Normal in iid case

Definition 2 If the tail function of a non-negative random variable X is F := 1 -F and F : R → [0; 1] satisfies that F ∈ RV -α with α > 0. Then F is said to be of secondorder regular variation with parameter ρ ≤ 0, if there exists a function Q(t) that ultimately has a constant sign with lim t→∞ Q(t) = 0 and a constant c ̸ = 0 such that

lim t→∞ F(tx) F(t) -x -α Q(t) = H α,ρ (x) = cx -α ∫ x 1 u ρ-1 du, x > 0 ( 33
)
Then it is written as F ∈ 2RV -α,ρ and Q(t) is referred to as the auxiliary function of F.

It is known from de [START_REF] De Haan | Generalized regular variation of second order[END_REF] or a more relevant form in Geluk et al. [START_REF] Geluk | Second-order regular variation, convolution and the central limit theorem[END_REF]

that if H α,ρ (x) is not a multiple of x -α then ρ < 0 implies that there exists a c ̸ = 0 such that H α,ρ (x) = cx -α x ρ -1 ρ
and |A| ∈ RV ρ and no other choices of ρ are consistent with Q(t) → 0. There are many distributions which satisfy the second order RV condition. These are (see e.g. [START_REF] Drees | How to make a hill plot[END_REF]):

Cauchy γ = 1, ρ = -2, Fréchet(1) γ = 1, ρ = -1, Student t(4) γ = 1/4, ρ = -1/2, t(10) γ = 1/10, ρ = -1/5 or loggamma γ = 1/3, ρ = 0.
The most seminal results on this topic could be found in de Haan and Ferreira (2006).

Theorem 1 Let X 1 , X 2 , . . . be iid random variables. Suppose that the second order condition 33 holds with α > 0, and auxiliary function

Q. If √ k n Q(n/k n ) → 0 as n → ∞ then, √ k n ( αH L k n ,n -1 ) d - → N(0, 8).
4 Empirical study of the effect of contamination on Pareto tail index.

The maximum likelihood estimator αML,n := 1

H kn,n
of Pareto tail index is very sensitive to deviations from the theoretical distributions, namely in the heavy-tailed class of distributions, see Stehlík et al. (2010b). It is unbiased, asymptotically consistent and has smallest variance in the set of all unbiased estimators for α. Its variance is equal to α 2 /n. However, it is not robust with respect to large contamination because its finite sample upper breakdown point is 0. Therefore, we are looking for an estimator αn with asymptotic relative efficiency

Var αML,n

Var αn = α 2 nVar αn ≤ 1

and close to 1. The latter means that the estimator's variance should be close to α 2 /n, although it will be larger than this value. It seems difficult to theoretically compare mean values and variances of Hill and t-Hill estimators for all distribution functions with regularly varying tails. In this section we empirically compare their properties.

A good discussion of this topic can be found in [START_REF] Finkelstein | Pareto tail index estimation revisited[END_REF]. First they make the Pareto probability integral transform of the initial data and then they use the properties of the Uniform distribution in order to obtain their estimator of α and to examine its properties. The biggest advantage of their approach is that in this way "even infinite contamination has a bounded effect on the transformed data". Robust methods considering contaminated distributions work much better, but sweep out the differences among distributions so that they do not use the complete prior information. The main concepts about the robustness are very well presented e.g. in [START_REF] Huber | Robust statistics, series in probability and mathematical statistics[END_REF]; Rousseeuw and Stahel (1986).

In order to compare the quality of the estimators we calculated the percentage relative bias of Hill, t-Hill and t-lgHill estimators (RB). It is defined as

RB = 1 m ∑ m i=1 ( γi -γ) γ 100,
where γi is correspondingly estimator of γ = 1/α for the ith simulated sample, i = 1, 2, ..., m.

Another characteristic that we use in order to determine the quality of these estimators is their relative root-mean-square error (RRMSE) [START_REF] Brzezinski | Robust estimation of the pareto tail index: a monte carlo analysis[END_REF].

RRMSE = √ 1 m ∑ m i=1 ( γi -γ) 2 γ 100, see e.g.
Analogously to the Hill plot we consider the set of points with coordinates

( k n , H * k n ,n ) , k n ∈ {1, 2, ..., n}.
Further on we call this plot 't-Hill plot'. Similarly for t-lgHill plot. The straight line represents the true value of γ = 1/α in all following figures in this section.

4.1 Independent observations.

t-Hill versus other robust estimators

In our Monte Carlo simulation we focus on parameter α of the classical (or Type I) Pareto distribution P(α, δ ), which is defined using cumulative distribution function as follows:

F(x) = 1 - ( δ x ) α , x ≥ δ > 0, ( 34
)
where δ is a scale parameter and α > 0 is a shape parameter, generally known as the Pareto tail index.

The existing literature offers various estimators for the Pareto tail index α. In our simulation study we focus primarily on the following estimators: Hill estimator (maximum likelihood estimator introduced by [START_REF] Hill | A simple general approach to inference about the tail of a distribution[END_REF], which we use as a nonrobust benchmark for the Pareto tail index), Integrated Squared Error Estimator (ISE, estimator introduced by Vandewalle et al. (2007), which is based on the relative excesses of observations above a certain treshold), Partial Density Component estimator (PDC, estimator introduced by Vandewalle et al. (2007), which is also based on the relative excesses of observations above a certain treshold), Least Squares estimator (LS), moment estimator (ME, see Dekkers et al. (1989)), QQ estimator (QQ, estimator based on the Quantile-Quantile approach, see [START_REF] Kratz | The qq-estimator and heavy tails[END_REF]), and finally Weighted MLE (WMLE, see [START_REF] Dupuis | Robust weighted likelihood estimators with an application to bivariate extreme value problems[END_REF]; [START_REF] Dupuis | A robust prediction error criterion for pareto modelling of upper tails[END_REF]). For the purpose of comparison we use the percentage relative bias (RB) and the percentage relative root-mean-square error (RRMSE) as defined in section 4.

Therefore, the data sets are simulated from Pareto distribution P(α, 1), for various values of α, as well as from contaminated Pareto distribution with function

F = (1 -ε)P(α 1 , 1) + εP(α 2 , 1) (35)
where ε = 0.05, 0.10 and α 1 = 1, 1.7, α 2 = 0.3, 0.5, 1. For all cases we assume the appropriate sample sizes n. The next investigations refine the conclusions made in [START_REF] Fabián | On robust and distribution sensitive hill like method[END_REF]. All estimators are calculated by means of their corresponding values in m = 100 samples. First we observe samples with size n = 40. The means of the Hill, t-lgHill and t-Hill estimators, separately for any fixed k n = 1, 2, ..., 39, together with their averages are given in Figure 1. They show that both estimators have similar properties for small samples, although the t-Hill oscillates stronger. The values of RB and RRMSE are given in Table 2. We also made analogous observations for samples with size n = 500. However, this figure is not provided here. The same conclusions could be made also from the values of RB and RRMSE, given in Table 3.

t-Hill versus Hill and t-lgHill in Pareto sample

t-Hill versus Hill and t-lgHill for distribution function with regularly varying tail at a logarithmic rate

Hill, t-lgHill and t-Hill estimators may perform very poorly if the slowly varying function in the tail is faraway from a constant. In (4.16) of [START_REF] Embrechts | Modeling extremal events for insurance and finance[END_REF] the authors consider

F ← (p) = (1 -p) -1/α (-ln(1 -p)), p ∈ (0, 1), ( 36
)
with respect to the Hill estimator. We simulated samples of n = 10000 observations of random variables with quantile functions (36), separately for α = 0.3, 1, 1.7 and plotted the Hill, t-lgHill and the t-Hill plots for k n = 5, 6, ..., 499. The rate of convergence of these estimators can be observed in Figure 2. Notice that the t-lgHill outperforms both Hill and the t-Hill estimators. Although Hill estimator seems to be slightly closer to the estimated value, it is clear that three estimators are not good for such c.d.fs, see

t-Hill versus t-lgHill and Hill for log-gamma sample

Figure 3 shows the t-lgHill for different α equal to 0.3, 1, 1.7 and log-gamma distribution (8), c = 1. Green line is for t-Hill, blue is Hill, red is lgHill.

t-Hill versus Hill for contaminated data from Pareto sample

Practical examples, see [START_REF] Fabián | On robust and distribution sensitive hill like method[END_REF] among others, have shown that data are expected to be contaminated in specific situations. Frequently there are outliers in the right tail of the distribution. It is known that the Hill estimator is not robust. This is due to the fact that data are entered by their logarithm in its calculation. But what about the t-Hill estimator? In (18) data are involved by their reciprocal value. Therefore, it is not difficult to deduce that these estimators are robust with respect to large values. They are sensitive to the center of the distribution.

Contamination, governed by α 2 . t-Hill estimator provides the superior results with respect to the relative bias (RB) and the relative root-mean-square error (RRMSE) for the contaminated Pareto distribution with distribution function ( 35) among all employed estimators ISE, PDC, LS, ME, QQ and WMLE.

Here we simulated small n = 40 and bigger samples of n = 200 data from contaminated Pareto distribution as given in (35) separately for α 1 = 0.3, 1 or 1.7 and α 2 = 0.5 or 1. The relative part of the contamination is governed by the parameter ε. We worked with two values 0.1 and 0.05. When considering a sample size with 10% contamination, this would yield ε = 0.1 . For sample sizes of 40, see Figure 4.

To obtain proper estimations for 1/α 1 we need k n to be large and k n < n. When α 1 < 1 and α 2 ≥ 1 the Hill estimators are better than the t-Hill estimators. When α 1 < 1 and α 2 < 1, or α 1 ≥ 1 and α 2 > 1 the t-Hill estimators are comparable to the corresponding Hill estimators. Finally if we observe the means of the estimators we can conclude that t-Hill estimators are better than Hill estimators for α 1 ≥ 1 and α 2 ≤ 1.

The same conclusions could be made on RB and RRMSE basis of both Hill and t-Hill estimators. They are given in Table 5 for k = 10 or 20 and for k = 100 or 200, and different values of α 1 , α 2 and ε.

Comparison of maximum likelihood with t-score moment estimators Let us compare maximum likelihood (ML) estimators with t-score moment (SM) estimators for two different cases, on the one hand based on a sample of one distribution and on the other hand in presence of contamination, i.e. distribution of the form

F ε = (1 -ε)F(ω) + εF(ω 1 ), (37)
where ω 1 > ω are score variances. Note: The score variance of the Pareto distribution is ω 2 = (α + 2)/α 3 . If ω = 1, we have 1/α = 0.657. Pareto distribution with score variance ω 2 will be denoted by Pa(ω).

Random samples were generated from heavy-tailed distributions, with increasing level ω 1 of contamination and α = 1.

Ratios of r(k) = H k n ,n /γ(ω * ) -1 and r * (k) = H * k n ,n /γ(ω * ) -1 for k n = 100, 250 and 900 and n = 1000 are given in Table 6.

The variances of SM estimates are somewhat higher in accordance with the values of theoretical asymptotic efficiencies, but SM estimates are less biased and thus, do not seem to be much worse than the ML estimates even for not contaminated samples.

However, when data contamination comes into consideration, i.e. ε = 0.1 of F( 10), the estimation with t-scores gives better results than standard ML-method. In addition to a smaller bias for very small sample sizes (n ≤ 25), variances are smaller for SM-method compared with ML-method the smaller the number of observations. Recall that for the two experiments new simulations have been performed such that a comparison of the goodness of the results between contaminated and noncontaminated data is obsolete.

t-Hill versus Hill for data from the log-gamma distribution

It is apparent that both, t-Hill and Hill, show systematic decrease with expecting heavier Pareto tail. Figure 5 shows paths of averages for H k n ,n , H * k n ,n and H L k n ,n from 50 experiments with contaminated log-gamma data from (37) with ω = 5.546 and ε = 0.05. The values of ω 1 in plots in Figure 5 are computed from log-gamma distri- bution.

Figure 5 shows paths of averages for H k n ,n , H * k n ,n and H L k n ,n from 50 experiments with contaminated log-gamma data from (37) with ω = 5.546 and ε = 0.05. The values of ω 1 in plots in Figure 5 It is apparent from Figure 5 that for small contamination t-lgHill (red) estimates γ from log-gamma data perform better than both Hill (blue) and t-Hill (green). However, for larger contamination it is definitely not robust, perhaps due to the squared

Infinite moving average process

In this subsection we consider the infinite moving average process with the following autoregressive form

X i = 1.3X i-1 -0.7X i-2 + Z i , i = 1, 2, ..., n. ( 38
)
It is discussed in Resnick and Starica (1993) with respect to the Hill estimator. Here we compare Hill estimator with the corresponding t-Hill estimator.

t-Hill versus Hill and lgHill for Pareto noise

The t-Hill and Hill plots of the averages of the corresponding estimators, calculated on samples of n = 2000 observations on moving average sequence (38) with noise Z i that is P(α, 1) distributed and k n = 10, 11, ..., 499 are given on Figure 6. The corresponding values of α are 0.3 in the first column, 1 in the second column and 1.7 in the third column. At first glance, we see that both estimators have a similar behavior for a fixed number of upper order statistics. However, the values of RRMSE and RB in Table 8 show that in this case the t-Hill estimators are better.

t-Hill versus Hill for contaminated Pareto noise

Contamination, governed by ε.

Here we contaminated the distribution of the noise component, calculated the means of the corresponding Hill and t-Hill estimators for k n = 10, 11, ..., 9999 and plotted them in Figure 7 and Figure 8 for relative parts on contamination ε = 0.1 and ε = 0.05. The observations in each sample are n = 10000. More precisely the distribution of the noise component Z i in ( 38) is ( 35). The corresponding values of RRMSE and RB are given in Table 9. Although both estimators give almost the same results, in most of these cases the t-Hill estimator is slightly better than the Hill estimator.

Conclusions of the Simulations

Having in mind the observations in section 4 we can make the following conclusions:

1. when the contamination has heavier tail than the original distribution, t-Hill estimator outperforms Hill estimator and several robust tail estimators (ISE, PDC, LS, ME, QQ and WMLE). In contaminated i.i.d. case, where the severity of the contamination is governed by the extremal index, the Hill estimator has more narrow confidence intervals than t-Hill estimator. When α 1 < 1 and α 2 ≥ 1 the Hill estimators are better than the t-Hill estimators. When α 1 < 1 and α 2 < 1, or α 1 ≥ 1 and α 2 > 1 the t-Hill estimators are comparable to the corresponding Hill estimators. Finally we can conclude that t-Hill estimators are better than Hill estimators for α 1 ≥ 1 and α 2 ≤ 1. This effect increases for small samples.

2. Both Hill and t-Hill estimators are asymptotically consistent, however the t-Hill estimator has bigger variance. 3. The reason that the Hill estimator is not robust seems to be the fact that in its formula the data are entered by their logarithm. But what about the t-Hill estimator? In (18) data are involved by their reciprocal value. Therefore, it is not difficult to deduce that these estimators are robust with respect to large values. They are sensitive to the center of the distribution. 4. When the i.i.d. data contamination comes into consideration, i.e. ε = 0.1 of F( 10), the estimation with t-scores gives better results than standard ML-method. In addition to a smaller bias for very small sample sizes (n ≤ 25), variances are smaller for SM-method compared with ML-method, especially for small samples. 5. In moving average case without contamination and with Pareto noise and α 1 ≤ 1 the t-Hill estimator is better with respect to the unbiasedness when n → ∞ and k n /n → 0, but it has bigger variance than the Hill estimator. 6. In moving average cases with contaminated Pareto noise the t-Hill estimators are better than Hill only in case α 1 ≥ 1 and α 2 ≤ 1 and small relative part of contamination. The smaller the relative part of the contamination, the faster the rate of convergence. Moreover in any cases the estimators converge to 1/min(α 1 , α 2 ) for n → ∞, k n → ∞ and k n /n → 0. 7. In case when the observed random variable have log-gamma distribution, the t-lgHill estimator is more robust that t-Hill and Hill estimators. This is due to the fact that the t-lgHill estimator is distribution sensitive while the Hill and t-Hill estimators by construction in this case are useful only due to the information about the regular variation of the distribution. 8. From all of the charts it visible that the t-lgHill estimator has relatively small variance. It seems to be almost the same as of the Hill estimator. 9. In Moving Average case with Pareto noise the t-Hill outperforms both Hill and t-lgHill estimators. When alpha increases, the variance of all three estimator decreases.

The larger the number of samples m, e.g. m = 1000, the more visible are the conclusions about Hill, t-lgHill and t-Hill estimators, however computational time is very large. We shall recall a general consensus that we shall have a multi-criterial approach in order to estimate the parameters of heavy tail. Consistency is an important property, but we shall also consider importance of approximation of finite sample distribution and robustness.

5 Application to a small sample data from Guanaco glacier A glacier is a solid ice mass which is fed by solid water (snow, hail or hoarfrost), transforms this solid water in ice and restores it via steam (evaporation/sublimation) or via liquid (water drained by the runoff stream-flow), see [START_REF] Francou | Métodos de observación de glaciares en los Andes tropicales, mediciones de campo y procesamiento de datos. Versión 1. Great-Ice, IRD Gascoin S[END_REF]. Such gain and loss of mass can be analyzed as a balance (or budget). Guanaco glacier is located in the III region of Chile (latitude 29S) in the semi-arid Andes, it has a surface area of 1.86[km 2 ] and a maximum thickness of 120[m], see [START_REF] Francou | Métodos de observación de glaciares en los Andes tropicales, mediciones de campo y procesamiento de datos. Versión 1. Great-Ice, IRD Gascoin S[END_REF], the area contributing water to the Pacific Ocean (included in the above calculation of 23%), is only 1.26[km 2 ], in other words, less than 35% of the total ice area in this zone. The rest of the Guanaco glacier contributes water to the Atlantic ocean. In this sense, the relative importance of glacier meltwater to the river basin as a complete unit is very low compared to the snow and rain contributions. Its geographical location is shown in Figure 9, this glacier (jointly with other minor glaciers and glacierets in the zone) is highly relevant for studies, because it is an important water reservoir in a semi-arid zone of Chile, contributing with up to the 23% of the streamflow in the corresponding valley, see Gascoin et al. (2011). All glaciers in the Pascua Lama region (3.64 km2) can contribute between 3 and 23 % to the very high basin area of the Huasco river. The glacier meltwater contribution goes from 100% at the foot of the glacier, down to near 0 at the sea outlet of this river (Huasco). Given that this zone is one of the world's most richly endowed territories of copper, there is high interest for the occupation of these water resources, resulting in a uneven mining-agricultural competition for water rights, see [START_REF] Oyarzún | Sustainable development threats, inter-sector conflicts and environmental policy requirements in the arid, mining rich, northern chile territory[END_REF]. The mass balance observations over the glacier consists of annual records for the period 2002-2014, the dataset corresponds to the negated mass balance because we are interested in extremely large losses of volume in the glacier (see Figure 10).

In order to study a threshold model for this data, estimations for the tail index γ are computed using Hill, t-Hill (as previously defined) and Zipf's estimator (see [START_REF] Kratz | The qq-estimator and heavy tails[END_REF]) which is a form of least square estimator (Schultze and Steinebach (1996)) of the tail index. It follows the fact that if we suspect that X 1,n ≤ ... ≤ X n,n are the order statistics from a Pareto family, then the plot of should be roughly linear with intercept 0 and slope γ. The Zipf's estimator is computed for a given k n , as

{( -log ( 1 - i n + 1 ) , log X i,n ) , 1 ≤ i ≤ n } ,
γzipf (k n ) = slope ({ -log ( 1 - i k n + 1 ) , log X n-k n +i,n , 1 ≤ i ≤ k n }) .
We used σ = X n-k n +1,n γ as estimator of the scale parameter. The asymptotic standard deviation is computed as

H kn,n √ k n for Hill's estimator, as γzipf (k n ) √ k n /2
for Zipf's, and for t-Hill as

H * k n ,n √ k n   1 + H * k n ,n √ 1 + 2H * k n ,n   .
Zipf's estimator plot (see Figure 11) shows linear stability up to k n = 7, which makes this value of k n a good option for further computations. Hill and t-Hill estimators were computed along with Zipf's estimator, taking into account k n = 6 and k n = 7. The stability of the parameter estimation in a small sample context was tested using a parametric bootstrap approach, where 1000 samples of size 10 were generated from GP distribution with shape parameter computed with the three previously discussed estimators. The results are displayed in form of histograms in Figure 12.

Table 10 summarizes the information about the EVI obtained from the simulation: 10 The initial value is the obtained using the estimator over the original sample, the mean and standard deviation accounts for the simulated samples.

Given the results for each estimator, it is clear that k n = 7 is a better choice than k n = 6 because estimations are less biased and all estimators are more close to each other, also, Zipf's estimator should be treated with more care in this context because it shows an important bias in the parametric bootstrap without an improvement in the variance. The main advantage of usage of t-Hill estimator is its robustness and good small sample properties.

In Figure 13 we plot Hill, t-Hill, Zipf and t-lgHill for 40+ data sample of the Pascualama zone.

In Figure 14 we provide a QQ-plot with assumed theoretical distribution of GPD with parameters: µ = 0, σ = 1, ξ ∈ {-0.5, -0.25, 0.25, 0.5}, with cdf defined for ξ

̸ = 0 as G(z) = 1 -[1 + ξ (z -µ)/σ ] -1/ξ .

Summary

In this paper we introduce and study properties of t-score estimator under different model settings. In particular, we derived asymptotic normality for t-lgHill estimator for iid case. We also derived consistency for t-Hill and t-lgHill under i.i.d. and moving average setup. We can conclude that:

i) The t-score is an unbounded function in cases of light-tailed distributions and a bounded function when distributions are heavy-tailed. Therefore, it seems to be extraordinary suitable for the estimation of parameters by the generalized moment method. It secures robustness 'when it is needed by the assumed distribution'. In this sense, the suggested estimates are 'naturally' robust without need of using any further function of robust statistics. In particular robustness of t-Hill can be well applied in cases when existence of first or second moment is questionable. Such an application for LATAM returns, i.e. daily log-returns (in percent) of an equity fund investing in Latin America (LATAM), can be found in Stehlík and Hermann (2015).

ii) It appeared that in cases of very small contamination, the maximum likelihood and t-moment method give similar results. However, the larger the contamination, the better are the t-moment estimates.

iii) It is apparent that an introduction of the t-mean and t-variance for distributions without ordinary mean and variance enables to describe the behavior of the estimates in a unique way and thus, compares results of the estimation for various assumed families of distributions having different parameters. iv) We propose t-Hill estimator for possibly dependent data and investigate its weak consistency and asymptotic behavior. Particularly we consider the infinite moving average model. Besides the i.i.d. case there are cases where the t-Hill estimator is more robust than the Hill estimator. We applied t-Hill for small sample size real data for Glacier and received satisfactory results despite a small sample size. Proof The result is proved analogously to Resnick and Starica (1993). An integration by parts entails The positive measures µ X,kn,n are Radon and the measure µ X,kn,n converges to µ and α > 0. Therefore, in view of (28) and the definition of µ, the Hospital rule yields Then, for all t ≥ 1,

Φ n = ∫ ∞ 1 φ(x) { 1 k n n ∑ i=1 ε { X i X (n-kn,n) ∈ dx }} = ∫ ∞ 1 φ(x) µ X,
lim n→∞ Φ n = lim n→∞ { ∫ t 1 µ X,kn,n (x, ∞] φ ′ (x) dx + ∫ ∞ t µ X,kn,n (x, ∞] φ ′ (x) dx } .
From the definition of b, n/k n → ∞ implies b(n/k n ) → ∞ and

n k n F(b(n/k n )) → 1
as n → ∞. Thus, in view of (30) and Potter's inequality for regularly varying functions, for any 0 < δ 1 < α there exist x 0 and n δ 1 large enough such that, for all x > max(1, x 0 ) and n > n δ 1 , Eµ X,kn,n ((1 -δ 0 )x, ∞] ≤ (1 + δ 1 )((1 -δ 0 )x) -α+δ 1 , and for all t > 1,

∫ ∞ t Eµ X,kn,n ((1 -δ 0 )x, ∞] |φ ′ (x) dx ≤ (1 + δ 1 )(1 -δ 0 ) -α+δ 1 ∫ ∞ t x -α+δ 1 |φ ′ (x)| dx.
We apply Fubini's theorem and obtain that for all t > 1 and n > n δ ,

P { ∫ ∞ t µ X,kn,n ((1 -δ 0 )x, ∞] |φ ′ (x)| dx > ε } ≤ 1 ε ∫ ∞ t Eµ X,kn,n ((1 -δ 0 )x, ∞] |φ ′ (x) dx ≤ 1 ε (1 + δ 1 )(1 -δ 0 ) -α+δ 1 ∫ ∞ t x -α+δ 1 |φ ′ (x)| dx.
As a consequence,

lim t→∞ lim n→∞ P { ∫ ∞ t µ X,kn,n ( X (n-kn,n) b(n/k n ) x, ∞ ] |φ ′ (x)| dx > ε, X (n-kn,n) b(n/k n ) -1 < δ } = 0
and (39) is satisfied.

A.2 Proof of Proposition 1

Proof The random variables X n , -∞ < n < ∞, are identically distributed. [START_REF] Cline | Infinite series of random variables with regularly varying tails[END_REF] proves that under these settings

F(x) = P ( ∞ ∑ j=0 c j Z j > x ) ∼ ∞ ∑ j:c j >0 c α j G(x) ∈ RV -α .
We use (26), apply the above Lemma 1 and complete the proof.

A.3 Proof of Proposition 2

Proof i.) The random variables X n , -∞ < n < ∞, are identically distributed. [START_REF] Cline | Infinite series of random variables with regularly varying tails[END_REF] proves that under these settings

F(x) = P ( ∞ ∑ j=0 c j Z j > x ) ∼ ∞ ∑ j:c j >0 c α j G(x) ∈ RV -α .
We use (32), apply Lemma 1 and complete the proof.

ii.) We use Slutsky's arguments and obtain weak consistency of the t-lgHill estimator in case of infinite moving average (MA) sequence.

  Fig. 1 t-Hill plots (green --), Hill plots (blue) and t-lgHill (red -.-) for α = 0.3 (left), α = 1(middle) and α = 1.7 (right), n = 40.

  Fig. 2 t-Hill plots (green --), Hill plots (blue) and t-lgHill (red -.-) for α = 0.3 (left), α = 1(middle) and α = 1.7 (right), n = 10000.

  Fig.3t-Hill plots (green --), Hill plots (blue) and t-lgHill (red -.-) for α = 0.3 (left), α = 1(middle) and α = 1.7 (right), n = 400.

Fig. 4

 4 Fig.4Rate of convergence of the t-Hill (first line) and Hill (second line) estimators with α 1 = 1, α 2 = 0.5 (left), α 1 = 1.7, α 2 = 0.5 (middle) and α 1 = 1.7, α 2 = 1 (right), ε = 0.1 for n = 40.

  α (see 8) in addition to their variances with different sample sizes averaged over 5000 experiments without (rows 1-4) and with contaminated Pareto distribution (rows 5-8), i.e. ε = 0.1.

Fig. 5

 5 Fig. 5 Graphs of averages for H kn,n (blue), H * kn,n (green) and H L kn,n (red) from 50 experiments with contaminated log-gamma data and increasing ω 1 , n = 1000.

  logarithmic term. Notice that t-Hill is more robust than the classical Hill even for non-Pareto data.

  Fig. 6 t-Hill plots (green --), Hill plots (blue) and t-lgHill (red -.-) for α = 0.3 (left), α = 1(middle) and α = 1.7 (right), n = 2000.

Fig. 7

 7 Fig.7Rate of convergence of the t-Hill (first line) and Hill (second line) estimators with α 1 = 1, α 2 = 0.5 (left), α 1 = 1.7, α 2 = 0.5 (middle) and α 1 = 1.7, α 2 = 1 (right), ε = 0.1, n = 10000.

Fig. 8

 8 Fig.8Rate of convergence of the t-Hill (first line) and Hill (second line) estimators with α 1 = 1, α 2 = 0.5 (left), α 1 = 1.7, α 2 = 0.5 (middle) and α 1 = 1.7, α 2 = 1 (right), ε = 0.05, n = 10000.

Fig. 9

 9 Fig. 9 Geographical location of Guanaco glacier and stakes where mass balance was measured in the period 2011-2014 along with aerial changes of the glacier computed via geodetic methods. Left image by Abermann et al. (2014), right image courtesy of the Center for Scientific Studies (CECs).

Fig. 10

 10 Fig. 10 Histogram for the dataset of mass balance in Guanaco glacier, n = 12.

Fig. 11

 11 Fig. 11 Left: Zipf's estimator plot for several values of k n . Right: Extreme value index computed over the dataset using 3 estimators and different values for k n .

Fig. 12

 12 Fig.12Histograms showing stability of the estimated tail index using different estimators, bin-width is chosen using the Freedman-Diaconis rule. Estimation over the original sample is shown as a vertical red line. Length of x-axis is fixed.

Fig. 13

 13 Fig. 13 Hill, t-Hill, Zipf and t-lgHill estimators (right) for 40+ data sample of the Pascualama zone.

  port of ANR project Desire FWF I 833-N18. The authors are very grateful to the Editor, Associated Editor and Referees for their valuable comments. lection on variation and correlation. Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 191:229-311, DOI 10.1098/rsta.1898.0007 Rabatel A, Castebrunet H, Favier V, Nicholson L, Kinnard C (2011) Glacier changes in the pascua-lama region, chilean andes (29 s): recent mass balance and 50 yr surface area variations. The Cryosphere 5(4):1029-1041 Resnick S, Starica C (1993) Consistency of hill's estimator for dependent data. Tech. rep., Technical report No. 1077, School of Operational Research and Industrial Engeneering Cornell University Itaca NY 14853 Rousseeuw FHERP, Stahel W (1986) Robust statistics: the approach based on influence functions. J Wiley, New York Schultze J, Steinebach J (1996) On least squares estimates of an exponential tail coefficient. Statistics & Risk Modeling 14(4):353-372 Stehlík M, Hermann P (2015) Letter to the editor. Annals of Applied Statistics 9:2051 Stehlík M, Potocký R, Fabián Z (2010a) On the favourable estimation of fitting heavy tailed data. Computational Statistics 25:485-503 Stehlík M, Potocký R, Waldl H, Fabián Z (2010b) On the favourable estimation of fitting heavy tailed data. Tech. rep., IFAS res. report Nr. 32 JKU Linz Stehlík M, Fabián Z, Střelec L (2012) Small sample robust testing for normality against pareto tails. Communications in Statistics -Simulation and Computation 41:1167-1194 Vandewalle B, Beirlant J, Christmann A, Hubert M (2007) A robust estimator for the tail index of pareto-type distributions. Computational Statistics & Data Analysis 51

  kn,n (dx) = -∫ ∞ 1 φ(x) d µ X,kn,n (x, ∞] = lim x→1 φ(x) µ X,kn,n (x, ∞]lim x→∞ φ(x) µ X,kn,n (x, ∞] + ∫ ∞ 1 µ X,kn,n (x, ∞] d φ(x).

  µ X,kn,n (x, ∞] = lim x→1 lim n→∞ φ(x) µ X,kn,n (x, ∞] = φ(1).

Table 1

 1

				n = 1000
		α = 0.5		α = 1		α = 3
		RB	RRMSE	RB	RRMSE	RB	RRMSE
	Hill	0.14	3.18	0.10	3.20	0.08	3.20
	tHill	0.14	4.27	0.08	3.67	0.06	3.30
	t-lgHill	0.43	3.48	0.78	7.03	2.69	21.83
	ISE	0.23	5.52	0.16	5.07	0.09	4.72
	LS	0.83	4.59	0.79	4.59	0.78	4.54
	ME	0.31	3.59	0.48	4.60	1.93	10.73
	PDC	1.02	14.07	0.53	9.40	0.25	6.68
	QQ	-0.90	4.49	-0.94	4.51	-0.95	4.46
	WML	3.48	8.30	3.43	8.20	3.35	8.08

presents the percentage relative bias (RB) and the percentage relative rootmean-square error (RRMSE) for the uncontamined Pareto distribution P(α, 1). If we suppose large sample size n = 1000, the results of RB and RRMSE of the t-Hill, the Hill estimator, ISE estimators are comparable. The largest relative bias as well as relative root-mean-square error have the ME and WML estimators.

Table 1

 1 Results of the percentage relative bias (RB) and the percentage relative root-mean-square error (RRMSE) for the various estimators of the Pareto tail index α with data drawn from uncontaminated Pareto distribution P(α, 1).

Table 2

 2 RRMSE and RB of Hill end t-Hill estimators for n = 40 i.i.d. data with Pareto distribution, for different α and k n .

	α	k n	estimator RRMSE	RB
	0.3	100	Hill tHill t-lgHill	9.92 16.27 21.65	-0.07 1.78 -1.94
		200	Hill tHill	7.04 11.35	-0.04 0.96
			t-lgHill	15.59	-1.02
	1	100	Hill tHill t-lgHill	9.90 11.65 21.68	0.05 0.73 -1.97
		200	Hill tHill	7.01 8.14	-0.01 0.31
			t-lgHill	15.48	-0.97
	1.7	100	Hill tHill t-lgHill	10.01 10.83 21.76	0.10 0.53 -1.74
		200	Hill tHill	7.13 7.66	0.00 0.17
			t-lgHill	15.62	-0.78

Table 3

 3

RRMSE and RB of Hill end t-Hill estimators for n = 500 i.i.d. data with Pareto distribution, for different α and k n .

Table 4

 4 , however the t-lgHill has the best performance.

	α	k n	estimator	RRMSE	RB
	0.3	100	Hill tHill t-lgHill	11.93 19.32 22.54	5.41 8.35 1.99
		200	Hill tHill	9.62 14.86	6.19 8.31
			t-lgHill	16.39	3.31
	1	100	Hill tHill t-lgHill	21.72 25.42 26.00	18.34 20.79 11.40
		200	Hill tHill	22.61 25.78	21.01 23.55
			t-lgHill	21.64	13.69
	1.7	100	Hill tHill t-lgHill	33.49 37.00 32.02	31.00 33.89 20.07
		200	Hill tHill	36.78 40.28	35.62 38.82
			t-lgHill	29.63	23.59

Table 4

 4

RRMSE and RB of Hill end t-Hill estimators for n = 10000 i.i.d. data with Pareto distribution, for different α and k n , m = 10000.

Table 5

 5 are computed from log-gamma distribution. RRMSE and RB of Hill end t-Hill estimators for n = 40 and 200 i.i.d. data with contaminated Pareto distribution, for different α 1 , α 2 , k n and ε.

	n	α 1	α 2	k n	ε	estimator	RRMSE	RB
				10	0.1	Hill tHill	45.39 56.03	18.42 22.61
					0.05	Hill	37.89	9.33
		1	0.5	20	0.1	tHill Hill tHill	47.62 31.02 34.57	14.63 13.800 14.14
					0.05	Hill	25.96	6.63
						tHill	30.74	8.35
				10	0.1	Hill tHill	41.01 37.50	-6.05 -14.73
					0.05	Hill	39.70	-22.78
	40		0.5	20	0.1	tHill Hill tHill	38.87 30.76 32.82	-26.63 -18.69 -26.24
					0.05	Hill	35.59	-29.53
		1.7		10	0.1	tHill Hill tHill	37.19 40.65 40.47	-32.92 -34.28 -32.34
					0.05	Hill	42.81	-37.59
			1	20	0.1	tHill Hill tHill	42.26 38.84 38.80	-35.43 -35.80 -35.26
					0.05	Hill	40.97	-38.44
						tHill	40.64	-37.59
				100	0.1	Hill tHill	18.49 16.87	13.80 10.51
					0.05	Hill	13.34	7.02
		1	0.5	200	0.1	tHill Hill tHill	13.61 13.00 11.57	5.60 9.99 7.32
					0.05	Hill	9.27	5.00
						tHill	9.35	3.74
				100	0.1	Hill tHill	21.55 29.11	-18.51 -27.88
					0.05	Hill	31.16	-29.86
	201		0.5	200	0.1	tHill Hill tHill	35.42 27.79 33.56	-34.63 -26.98 -33.14
					0.05	Hill	34.60	-34.16
		1.7		100	0.1	tHill Hill tHill	37.54 36.46 37.19	-37.22 -35.84 -36.53
					0.05	Hill	38.94	-38.43
			1	200	0.1	tHill Hill tHill	39.20 37.36 38.01	-38.62 -37.07 -37.70
					0.05	Hill	39.37	-39.12
						tHill	39.66	-39.38

Table 6

 6 Values of r(k) = H kn,n /γ(ω * ) -1 and r * (k) = H * kn,n /γ(ω * ) -1 are computed for contaminated distributions (P cont ) with three different ω * and ε = 0.05 and 0.1, respectively.

	ε	ω * γ(ω * ) r(100) r * (100) r(250) r * (250) r(900) r * (900)
		1.5	0.879	.0363	.0376	.0235	.0187	.0184	.0170
	.05 3	1.500	.2027	.1517	.1333	.0982	.0679	.0473
		5	2.165	.4982	.3049	.2792	.1573	.1203	.0656
		1.5	0.879	.0683	.0707	.0502	1.0458	.0339	.0296
	.10 3	1.500	.3842	.2828	.2512	.1799	.1340	.0970
		5	2.165	.9298	.6250	.5536	.3352	.2383	.1310
		contamination estimate		n=12	n=25	n=50	n=75
				1/H kn,n		1.086	1.041 1.023 1.014
		ε = 0		1/H L kn,n Var (1/H kn,n ) 0.342 1.052	1.028 1.010 1.009 0.215 0.144 0.116
				Var (1/H L kn,n ) 0.371 0.241	0.166	0.135
				1/H kn,n		1.298	1.188 1.129 1.127
		ε = 0.1	1/H L kn,n Var (1/H kn,n ) 0.203 1.011	1.090 1.198 1.144 0.073 0.030 0.020
				Var (1/H L kn,n ) 0.105 0.061	0.037	0.022

Table 7

 7 Comparison of ML and SM estimators of

Table 8

 8 RRMSE and RB of Hill end t-Hill estimators for n = 10000 moving average data with Pareto distribution, for different α and k n .

	α	k n	estimator	RRMSE	RB
	0.3	100	Hill tHill t-lgHill	33.11 33.87 45.93	-10.26 -0.69 -35.38
		200	Hill tHill	24.45 23.68	-4.09 3.38
			t-lgHill	39.80	-23.62
	1	100	Hill tHill t-lgHill	20.52 18.77 36.65	-2.19 0.01 -14.47
		200	Hill tHill	16.03 15.77	-1.11 -0.05
			t-lgHill	30.49	-7.75
	1.7	100	Hill tHill t-lgHill	20.58 18.02 45.27	-3.69 -4.98 -1.84
		200	Hill tHill	15.25 14.27	-6.56 -8.15
			t-lgHill	37.13	-0.03

Table 9

 9 RRMSE and RB of Hill end t-Hill estimators for n = 10000 moving average data with contaminated Pareto distribution, for different α 1 , α 2 , k n and ε.

	α 1	α 2	k n	ε	estimator	RRMSE	RB
			2000	0.1	Hill tHill	55.92 47.05	54.32 46.10
				0.05	Hill	32.79	30.88
	1	0.5	3000	0.1	tHill Hill tHill	25.38 52.85 47.37	24.30 51.80 46.77
				0.05	Hill	31.15	29.97
					tHill	27.50	26.85
			2000	0.1	Hill tHill	197.34 166.95	195.51 165.59
				0.05	Hill	120.22	117.81
		0.5	3000	0.1	tHill Hill tHill	79.98 169.79 133.34	78.63 168.45 132.47
				0.05	Hill	95.57	93.88
	1.7		2000	0.1	tHill Hill tHill	61.26 11.52 8.10	60.42 9.67 6.41
				0.05	Hill	5.64	-0.84
		1	3000	0.1	tHill Hill tHill	5.31 11.13 9.41	-2.98 9.96 8.50
				0.05	Hill	4.47	0.44
					tHill	3.67	-0.06

  Estimation over the original sample is shown as a vertical red line. Length of x-axis is fixed.

		Estimator Initial Mean Std.Dev. Asymp.sd.
		Hill	0.436	0.429	0.171	0.178
	k n = 6	t-Hill	0.515	0.541	0.241	0.224
		Zipf	0.315	0.401	0.210	0.257
		Hill	0.412	0.413	0.153	0.156
	k n = 7	t-Hill	0.466	0.485	0.198	0.186
		Zipf	0.359	0.441	0.208	0.271
	Table

Acknowledgements Zdeněk Fabián thanks to institutional support RVO:67985807. Milan Stehlík acknowledges the support of Fondecyt Proyecto Regular No. 1151441. Philipp Hermann thanks to the sup-

Let us first remark that

Second, in order to apply Theorem 4.2. of [START_REF] Billingsley | A simple generalization of the hill estimator[END_REF], we have to check that lim t→∞

Our aim is to replace µ X,kn,n with µ X,kn,n . To this aim, we choose δ > 0 and consider

Clearly,

, converges to zero when n → ∞. In view of ( 27), for any 0 < δ 0 < 1 + α there exists n δ 0 large enough such that for all t > 1 and n > n δ 0 ,

Chebishev's inequality guarantees that

A.4 Proof of Theorem 1

Proof Recall that the t-lgHill estimator H L kn,n of γ = α -1 can be written as

while the Moment estimator is defined as

see (3.5.9), p.102, in de Haan and Ferreira (2006), also (1.7) in page 1834 of Dekkers and de Haan (1989)).

The difference between the two estimators appears clearly when comparing ( 40) and ( 41). In the following, we shall also make use of γ -:= γ M -M

(1) n , see Remark 3.5.7 in de Haan and Ferreira ( 2006)), which is an estimator of γ -:= min(0, γ). It is thus clear that the t-lgHill estimator (40) c an be rewritten as

Corollary 3.5.6 of de Haan and Ferreira (2006) states that

where Q, P are Gaussian distributions given in Lemma 3.5.5 of de Haan and Ferreira (2006). Using the previous notations and remarking that γ -= 0, this result can be simplified as

n -γ + ) where γ + := max(0, γ) = γ = 1/α. It is easily seen that ξ ′ n can be expanded as ξ

The limit distribution of the first term is established in (3.5.24), page 109 of de Haan and Ferreira (2006):

The limit distribution of the second term has already been established to be 1 2 Q -2P and therefore ξ ′ n converges in distribution to γP. Plugging ξ n and ξ ′ n in (42) yields

). According to Lemma 3.5.5 of de Haan and Ferreira (2006)


Benoît Chappet De Vangel

Cesar Torres-Huitzil

Cinvestav Tamaulipas

Bernard Girau

Stochastic and Asynchronous Spiking Dynamic Neural Fields

Stochastic computing was extensively studied for artificial neural networks (ANN) implementation with a good time/area trade-off on up-to-date FPGAs. We propose here to use the same paradigm for the hardware implementation of dynamic neural fields (DNF) on FPGAs. The all-to-all connectivity of these neural population models make straight-forward hardware mappings impossible for high density fields. It is necessary to adapt the architecture to fit the cellular nature of computing substrates such as FPGAs. Following the previous work on randomly spiking dynamic neural fields, we propose here a new implementation inspired by stochastic ANNs. We introduce here the Cellular Array of Stochastic Asynchronous Spiking DNF model, or CASAS-DNF. While keeping the fully decentralized cellular characteristics, this new approach is much more competitive in terms of speed and area. We also show that the basic behaviors of DNFs are preserved. The low hardware cost and the cellular design of this model make it easily scalable.

I. INTRODUCTION

Artificial neural networks (ANNs) are well known as a tool inspired from biology and applied successfully in several fields such as pattern recognition, language processing, classification, etc. Numerous hardware options have been developed to rise the processing power of these ANNs. Different implementations include analogue or digital hardware substrate [START_REF] Misra | Artificial neural networks in hardware: A survey of two decades of progress[END_REF]. Among digital hardware FPGAs are often retained as the best trade-off between logic density and design time.

The usual limiting factor of hardware implementation of ANN is the synaptic density. Usual ANN have 10 to 1000 time more synapses than neurons and their hardware implementation is costly as a synapse requires multiplying a neural activity by a specific synaptic weight. Several options were proposed to address this issue: time multiplexing, serial arithmetic or stochastic arithmetic. Stochastic arithmetic (reviewed in [START_REF] Alaghi | Survey of stochastic computing[END_REF]) describes a real value with a bit-stream where a bit has a probability p to be high and 1 -p to be low. The advantages are numerous as the operators for complex operations are very light. The multiplication can be computed with a simple AND gate while the sum can be approximated by an OR gate (to some extent) or by a multiplexer. These principles have been successfully applied to FPGA-based ANN implementations (see [START_REF] Van Daalen | A stochastic neural architecture that exploits dynamically reconfigurable fpgas[END_REF], [START_REF] Bade | Fpga-based stochastic neural networksimplementation[END_REF]).

We focus here on using the same principles to implement a widely used tool for bio-inspired cognition: the dynamic neural fields (DNF) introduced by Beurle [START_REF] Beurle | Properties of a mass of cells capable of regenerating pulses[END_REF] and Wilson [START_REF] Wilson | Excitatory and inhibitory interactions in localized populations of model neurons[END_REF]. The DNF is a mathematical tool built on the continuous neural field theory (CNFT), which models the layers of cortical columns in the brain as a continuous field. Many applications exist, especially for perceptual tasks and autonomous robotics [START_REF] Erlhagen | The dynamic neural field approach to cognitive robotics[END_REF].

In [START_REF] Chappet De | Randomly spiking dynamic neural fields[END_REF], we introduced the randomly spiking dynamic neural fields (RSDNF). The idea was to adapt DNF to hardware in a cellular way. We used randomly propagated spikes in a regular mesh of von Neumann connected neurons so as the spike propagation and randomization would induce a virtually fully connected neural network with a Gaussian shaped lateral weights function. If functionally the random behavior was not altering the model abilities, the hardware performance of this decentralized model was not competitive compared to a centralized spike routing approach [START_REF]Spiking dynamic neural fields architectures on fpga[END_REF]. Besides this attempt, no dedicated hardware implementation has been proposed for dynamic neural fields. In this paper we propose a new approach based on stochastic computation. It maintains the idea of using a regular mesh but instead of propagating spikes randomly, stochastic bit streams are propagated taking advantage of the simplicity of their operators. This approach has resulted in the cellular array-based stochastic and asynchronous spiking dynamic neural field (CASAS-DNF) model that we describe in this paper.

II. RELATED WORK

A. Hardware implementations of stochastic artificial neural networks

In the 90's, two works explained the potential of stochastic neurons (bit-stream neurons or pulse mode neuron) for ANN FPGA implementation, see [START_REF] Van Daalen | A stochastic neural architecture that exploits dynamically reconfigurable fpgas[END_REF] and [START_REF] Bade | Fpga-based stochastic neural networksimplementation[END_REF]. But the area consumed by pseudo random number generators was limiting the applications at that time. With the rising power of FPGAs, stochastic ANNs were studied more recently by [START_REF] Hikawa | A digital hardware pulse-mode neuron with piecewise linear activation function[END_REF], [START_REF] Nedjah | Fpga-based hardware architecture for neural networks: binary radix vs. stochastic[END_REF] and [START_REF] Rossello | Hardware implementation of stochastic-based neural networks[END_REF].

The idea of stochastic computing is to represent numbers with the frequency of '1's in a bit-stream. For instance we can represent the number v = 1/3 with a bit stream of size s = 9 with 3 '1's and 6 '0's. More generally the bit stream of arbitrary size s and value v will be generated with a probability p = v. The definition can then be extended to encode values v > 1 (unipolar representation) and even negative values (bipolar representation). The main advantage of stochastic bit streams is that multiplication can be implemented with a single AND gate. Neuron activation functions can also be implemented with simple operators. The sum is less straight forward than the multiplication: the OR gate with two bit streams A and B gives a bit stream A + B -AB (see fig 1). Another solution is to use a multiplexer which will return a bit stream encoding (A + B)/2. To encode numbers 0,0,1,1,0,1,1,0 (4/8) A 1,1,1,0,1,0,1,0 (5/8) B 0,0,1,0,0,0,1,0 (2/8) AB into stochastic bit streams, stochastic number generator (SNG) requires a pseudo random number generator (PRNG). The pseudo random numbers are then compared to a probability p and the result will generate a bit stream of probability p.

(a) Stochastic Multiplication 0,0,1,1,0,1,1,0 (4/8) A 1,1,1,0,1,0,1,0 (5/8) B 1,1,1,1,1,1,1,0 (7/8) A + B -AB (
The quality of the PRNG and the precision of the comparison with p are two parameters which will influence the precision of stochastic computation. The number of PRNG has to be scaled with the number of SNG as the random numbers correlation leads to critical drop of stochastic bit stream quality. The historical approach to produce pseudo random numbers used linear-feedback shift registers [START_REF] Golomb | Shift register sequences[END_REF]. Their quality can be very good but their area cost has motivated research for new biomorphic cellular hardware (see [START_REF] Suri | Bio-inspired stochastic computing using binary cbram synapses[END_REF] and [START_REF] Kavehei | Highly scalable neuromorphic hardware with 1-bit stochastic nano-synapses[END_REF]). To decode a stochastic number, the stochastic number decoder (SND) has to count the number of ones in a bit stream with a simple counter.

In summary the operators to manipulate stochastic bitstreams are very competitive in term of area except for the SNG which is heavier and critical for computation quality. The quality of the stochastic representation depends on the size of the bit-streams, the quality of the stochastic number generators, the quality of the operators and the bit stream correlation. This limits will have to be tuned carefully to optimize the speedarea trade-off.

B. Cellular computing and spiking neural networks

Despite their intrinsic parallelism and the simplicity of their computation, neural networks are not easy to map onto parallel hardware devices. With the rising computational power needed for precise simulations of neuro-biological models, a whole field of science, namely neuromorphic engineering [START_REF] Indiveri | Frontiers in neuromorphic engineering[END_REF], is trying to find new hardware paradigms for largescale spiking neural networks simulation. The main problem for vast neuron dynamic simulation is the non-locality of the computation which is highly incompatible with state of the art computing devices. The spiking neural networks (SNN) address partly this issue by reducing the non-local informations size to 1 bit. Specific networks on chip (NoC) for SNN were developed to improve the speed of spike propagation based on the Address Event Representation paradigm (AER) [START_REF] Mahowald | Vlsi analogs of neuronal visual processing: A synthesis of form and function[END_REF]. The idea is to sequentially address every spike to its destinations using a dedicated high speed bus or a binary tree [START_REF] Vainbrand | Scalable network-on-chip architecture for configurable neural networks[END_REF]. These solutions can be an inspiration for hardware artificial network implementations, and the speed-area trade-off is generally very competitive. In [START_REF]Spiking dynamic neural fields architectures on fpga[END_REF], we derived an AER-based hardware implementation of spiking DNFs with good performance.

However, this sequential processing of spikes is not satisfying as it lacks the main characteristic of the brain: its decentralization. We are more interested in the advantages of decentralized cellular hardware for their numerous advantages [START_REF] Sipper | The emergence of cellular computing[END_REF]. It perfectly fits our current hardware paradigm with the local computation and good parallelism. It is fault tolerant as the loss of one element does not alter the general behavior. And finally, the computation is an emergent property of the system exactly like in the brain. This emergence confers a very competitive robustness to noise and is thus compatible with the introduction of stochastic computation.

Attempts to mix cellular computing and neural networks already exist: locality in hardware implementation of ANNs was achieved with FPNAs (in [START_REF] Girau | Neural networks on FPGAs : a survey[END_REF] and in [START_REF] Bohrn | Field programmable neural array for feed-forward neural networks[END_REF] with up-to date boards), while cellular neural networks are already well-known [START_REF] Chua | Cellular neural networks: Applications[END_REF] but they handle purely local communications between neurons.

III. FROM DNFS TO HARDWARE FRIENDLY DNFS

A. Dynamic neural fields

Dynamic neural fields (DNFs) were first introduced to model neural activity at a mesoscopic level of a cortical column. It is an approximation initially based on the very simple assumption of the continuum neural field theory (CNFT) introduced by Wilson [6] and developed by Amari [START_REF] Amari | Dynamics of pattern formation in lateral-inhibition type neural fields[END_REF]. In this theory, the neural population of a cortical column is approximated by a locally cooperative and globally competitive neural field. The neurons are rate-coded and their potential (u x,t for neuron x at time t) depends on the whole potential field M , following the equation (already discretized with an Euler method):

u x,t+dt = u x,t + dt τ (-u x,t + M ω(||x -x ||)σ(u x ,t )dx +I x,t + h), (1)
where dt is the discretizing step, I x,t the input feeding, h the resting potential, τ the time constant and ω is the lateral weights function. σ is the activation function which can be a sigmoid and ω models inhibitory and excitatory synaptic weights according distance d between neurons. This lateral weights function ω is a difference of Gaussians (DoG):

ω(d) = Ae d 2 a 2 -Be d 2 b 2 , A, a, B, b ∈ R * + . (2)
The dynamic of the DNF (reviewed in [START_REF] Coombes | Waves, bumps, and patterns in neural field theories[END_REF]) is very rich and interesting. The characteristics we are interested in for visual processing applications are the selection and tracking abilities [START_REF] Rougier | Emergence of attention within a neural population[END_REF]. It relies on a trade-off between the afferent feeding (the inputs I) and the lateral feeding term M ω(||xx ||)σ(u x ,t )dx . The global competition within the field is the negative part of the lateral weights function (representing inhibitory connections), while the local cooperation is the positive part of ω (representing local excitatory connections).

Taken as it is, a straightforward hardware implementation of a DNF would be problematic as it would need a very costly convolution operator to compute the lateral feeding.

B. Spiking dynamic neural fields

The first step towards an easier hardware implementation was to introduce spiking neurons which would exchange only binary informations. Spiking dynamic neural fields have been introduced by [START_REF] Chevallier | Visual focus with spiking neurons[END_REF] and developed by [START_REF] Vazquez | Visual attention using spiking neural maps[END_REF]. This model is based on Leaky Integrate-and-Fire (LIF) neurons which have the property to fire a spike and reset their potential when the potential reaches a threshold. The discretized version of the SDNF equation follows:

u x,t+dt = u x,t + dt τ (-u x,t + I x,t + h) + I syn x,t . (3)
The lateral influence is computed by applying the synaptic weight (ω M ) of eq. 2 to each received spike(S):

I syn x,t = x ω(||x -x ||)σ(u x ,t ). ( 4
)
Finally the neurons emit a spike when their potential reaches a threshold θ.

σ(x) = 1 if x ≥ θ 0 if x < θ. ( 5
)
If a spike is fired, the potential is reset:

u x,t = h.
The dynamics being harder to describe with analytical tools, more work need to be done to fully understand their behavior. However, basic abilities like selection, tracking and robustness were demonstrated to surpass DNF abilities [START_REF] Vazquez | Visual attention using spiking neural maps[END_REF].

C. Randomly spiking dynamic neural fields

Despite their reduced 1-bit bandwidth, the spiking dynamic neural fields are fully connected neural networks. Therefore their hardware implementation remains problematic as we need a personalized all-to-all connectivity in the network.

The first attempt to reduce neuron connectivity to 4 neighbors was introduced with the randomly spiking dynamic neural field model (RSDNF) [START_REF] Chappet De | Randomly spiking dynamic neural fields[END_REF]. To simulate global connectivity with local 4-neighbors connectivity, we need to be able to simulate as precisely as possible the lateral weights function ω. As the connections are local, information between neurons (in our case spikes) can only be transferred from cell to cell. This induces a propagation delay to wait for every information to reach its destination.

In RSDNFs, the idea is to simulate the DoG's shape by using probabilities. For that, every spike activation information is propagated from the activated neuron to the borders of the map with a probabilistic network on chip (NoC). A spike is subdivided in N packets of one bit, which are then sent through the network where every hop depends on a probability p. The NoC is designed such that there can only be one path between two neurons a and b and such that its distance (in number of hop) is the Manhattan (or XY) distance between the neurons:

dist XY (a, b) = |a x -b x | + |a y -b y |.
There is one NoC for the excitatory weights (the probability of each hop is p a ) and another one for the inhibitory weights (the probability of each hop is p b ). Consequently the decentralized lateral weights function is, on average,

ω(a, b) = ω(d = dist XY (a, b)) = Ap d e -Bp d i . ( 6
)
The shape of the corresponding curve is different from equation 1 (see Fig. 2), but new optimal parameters where derived from optimization methods (see [START_REF] Quinton | Exploring and Optimizing Dynamic Neural Fields Parameters Using Genetic Algorithms[END_REF] and [START_REF] Fix | Template based black-box optimization of dynamic neural fields[END_REF]). Behavioral simulation results confirmed the validity of the approach and the scalability of the cellular architecture of this model. However, the RSDNF speed×area product is not very competitive, and asynchronous neuron updates are not possible as a neuron has to wait for the end of spikes propagation to compute its next potential. This can be a problem to extend the RSDNF to synaptic delayed DNF [START_REF] Hutt | Activity spread and breathers induced by finite transmission speeds in two-dimensional neural fields[END_REF], or to implement asynchronously computing DNFs [START_REF] Taouali | No clock to rule them all[END_REF].

IV. THE CASAS-DNF MODEL

As in the RSDNF paradigm, the hardware mesh is regular and composed of relatively simple identical cells. Each cell contains one neuron, two dendrites, eight synapses and eight routers (one inhibitory and one excitatory for each neighbor). The cells are disposed on a regular mesh connecting each other to their north, south, east and west neighbors.

A. Stochastic decentralized lateral weights

The challenge of cellular simulation for a globally connected neural networks is to have a good simulation of the synaptic weights between every neuron: in our case it is the lateral weights function ω. To assert the quality of the model we introduce in this work two criteria: the curve comparison (we compare our theoretic lateral weights function with the DNF one), and behavioral comparison (we compare our model dynamic with the standard model dynamic on a set of simulation scenarios). In this paragraph we present the analytical results for the shape of the CASAS-DNF lateral weights.

As we saw in eq. 2, the DNF lateral weights function ω only depends on the distance between neurons. Thus, following the idea introduced in RSDNF, we propose to simulate the lateral weights using probabilities that decrease with distance, since it was proved to have a good curve similarity (Fig. 2) and a similar behavior despite the introduced randomness (see [START_REF] Chappet De | Randomly spiking dynamic neural fields[END_REF]).

We are using stochastic computation (or bit-stream arithmetic) to reproduce the lateral weights decreasing with the distance. The representation of a real number a is made by a bit-stream a bs generated with a probability a p = a, such that P (a bs = 1) = a p = a. To simplify the demonstration, we consider that the value of a bit stream is its probability.

Every activated neuron generates a bit-stream δ bs of size s representing a spike of real value δ with a probability δ p = δ. The spike will then be propagated throughout the network, being multiplied by synaptic weights at each hop. The synaptic weights do not depend on the neuron or synapse: the constant weight W e is used for the excitatory layer, and W i for the inhibitory layer.

To compute the lateral weights function for one spike we need to model the excitatory bit-stream received by the excitatory dendrites (E) and by inhibitory dendrites (I), assuming these dendrites belong to a neuron at a distance d from the activated neuron:

P (E bs = 1) = E = δW d e P (I bs = 1) = I = δW d i (7)
Hence, the probabilistic lateral feeding equation integrated by the neuron is:

ω M (d) = I e δW d e -I i δW d i , (8)
where I e and I i are respectively the intensity of excitatory spikes and inhibitory spikes, defined by:

I e = A δ and
I i = B δ (9)
As this lateral weights equation is very similar to (eq. 6), we will use similar parameters. However, the bad quality of stochastic addition performed by the OR gate will bias these lateral weights.

The cellular nature of the network and the connectivity graph of the routers result in horizontal routers (east and west) having 4 predecessors, while vertical routers (north and south) have 2 predecessors (see fig. 4). Every router generates its outputs by summing the predecessors bit-streams, and as we saw in the stochastic computation paragraph, there is no perfect way to perform stochastic bit-stream addition. The reason is that as we sum numbers in the interval [0, 1], the sum can be greater than the output bit-stream interval which is also [0, 1].

Let A be an activation patch of n neurons (a set of n contiguous neurons that spike simultaneously). Let b represent one neuron at distance d from the patch (to simplify we assume that the neuron b is at distance d from every neuron of the patch). If we adapt the bit-stream integration in eq. 7 for the integration of n spike bit-streams instead of one, we have:

E = nδW d e I = nδW d i (10)
It is important to note that the bit-stream additions are performed in several places: in the routers when several bitstreams are going through the same router and in the dendrite module as it receives the 4 neighbors bit-streams. Equation 10is valid for a standard addition. However, as we are using OR gate for addition, we need to adapt it.

The generalized stochastic sum on n stochastic bit-streams of probability p using OR gate is

n p bs = 1 -(1 -p) n . (11)
Consequently eq. 10 becomes

E = 1 -(1 -δW d e ) n I = 1 -(1 -δW d i ) n (12)
Modifying eq. 8 to take into account the integration of n bit-streams using OR gates, the lateral influence of the patch of neurons A on neuron b becomes:

ω M (A, b) = I e (1 -(1 -δW d e ) n ) -I i (1 -(1 -δW d i ) n ). ( 13
)
This is an approximation as we assumed that all the neurons of the patch A were at the same distance from neuron b and we did not take into account the precision of the bit-streams (their length s) which influences the quality of encoding and decoding. But it is enough to predict the influence of the different parameters on the lateral weights quality of the model. On Fig. 3 we see that the error rises with δ and n. When δ is greater than 0.01 and n is greater than 5 the error with the RSDNF lateral weights function is too important to be neglected. Consequently a new set of parameter has to be found. In this work we used particle swarm optimization techniques to find optimal parameters for the set of task we are presenting in the results section.

B. Spiking neuron

The dendrite module is a stochastic number decoder (SND). This modules receives the 4 neighbors bit-streams, sum them with an OR gate and add a quantum of excitation (q e ) or a quantum of inhibition (q i ) to a register when the output of the OR gate is high (q e = I e /s and q i = I i /s) . When the neuron computation is activated, this lateral influence is integrated in the neuron potential, and the dendrite register is reset. The spiking neuron computation was made as close to the spiking neuron definition as possible (see Fig. IV-B). The only modification was made to transform the leaking factor division by a right bit shift to avoid a division. Details are available in [START_REF] Chappet De | Randomly spiking dynamic neural fields[END_REF].

C. Pseudo random numbers generation

The stochastic number generator module (SNG) uses pseudo random numbers (PRNs). We have shown in [START_REF] Chappet De | Randomly spiking dynamic neural fields[END_REF] that an 8-bit precision is enough for synaptic weight probabilities. However, the efficient generation of pseudo random numbers is often a complex question. Standard implementations use linear feedback shift registers but we are using massively replicated cellular automata which are a scalable and distributed way to generate multiple high quality PRNs [START_REF] Girau | Tiled cellular automata for area-efficient distributed random number generators[END_REF]. To be fair with 3a show the error with a rising value for spike probability δ and a fixed number of activated neurons n = 9 and Fig. 3b show the error with a rising number of activated neuron n and with δ = 0.1.

RSDNFs, we never use twice a random bit, but it can be noted that a high random numbers re-utilization is compatible with our simple applications. For this architecture we need to generate four random numbers per synapse and per propagation cycle and one random number per neuron for the spike bit-stream generation. For a network of 35x35 neurons with excitatory and inhibitory synapses, it results in 11025 random numbers with 8 bits minimum. These 88200 random bits are simultaneously generated by Cellular Automata-based Pseudo-Random Number Generators (CAPRNG) that are very well adapted to FPGA devices and to our cellular computation requirement: non-homogeneous cellular automata with high random quality emerging from simple logic rules:

σ t+1 (i,j) = X ⊕ C.σ t (i,j) ⊕ N.σ t (i-1,j) ⊕ W σ t (i,j-1) ⊕S.σ t (i+1,j) ⊕ Eσ t (i,j+1) (14)
where σ t c is the state of the cell c at time t, and where Design of the routers and their connectivity. The CASAS cells (white circles) are organized on a 2D-grid. Every router's direction (square) contains one OR gate performing a sum approximation of the inputs stochastic bit-streams. The vertical OR get input from the opposite direction neighbor while horizontal OR get input from opposite direction and vertical directions neighbors. Every direction receive as well an input from the spike's stochastic bit stream generated by the CASAS Cell's neural activation (gray circle). The CASAS map contains two layer of this routers: one for the excitatory bit streams and one for the inhibitory bit streams. XCN W SE ∈ [0, 1] 6 are the control variables determining the rule of each cell in the cellular automaton. In [START_REF]Evolution of 2-dimensional cellular automata as pseudo-random number generators[END_REF], [START_REF] Vlassopoulos | A Metric for Evolving 2-D Cellular Automata As Pseudo-Random Number Generators[END_REF], genetic algorithm were used to find the best rules for 8x8 CA grid. This PRNG shows very good behavior in standard tests while offering an outstanding compacity for FPGA implementation.

V. RESULTS

The CASAS-DNF model was implemented in software and hardware to study more precisely its dynamics and performance. Software simulation was programmed with a Python/C++ framework simulating the cellular computation of CASAS-DNF. The simulation parameters are shown in Table I.

A. Behavioral results

Behavioral software simulations of the two models SDNF and RSDNF have already been reported. In [START_REF] Vazquez | Visual attention using spiking neural maps[END_REF] a behavioral study shows the properties of spiking DNF for basic visual attention tasks. In [START_REF] Chappet De | Randomly spiking dynamic neural fields[END_REF], the random behavior of the RSDNF is compared to the behavior of spiking DNF to assert the validity of the cellular approach in basic visual tasks. A short summary of the main results is presented here.

We consider three different scenarios to test the error during basic visual selective attention and tracking tasks. The error is measured with the euclidean distance between the barycenter The leak is applied to the old potential with a division by 8 using a shift right logical (SRL) of 3 bits. Once the leak has been removed from the potential, the new potential is compared with the activation threshold δ. If the neuron is activated, the potential will be reset. Otherwise, the new potential is output and will update the potentials register.

of the field's activated neuron coordinates and the stimulus the field is focusing at. For instance we can see on Figure 7 the tracking scenario and the measured error. After nine computation iterations the neural field is spiking close to one stimulus: this is the selective attention or focus. When the stimulus circles, the activation of the field follows the stimulus. Figure 8a shows 50 iterations of tracking scenario for different sizes of the bit-stream s for CASAS-DNF and RSDNF. We can see that for CASAS-DNF s is optimal around 100. For Fig. 8b, we chose s = 100 and we repeated 50 experiments on the three different scenarios. The error bars show the bootstrap 95% confidence intervals. The model seems to be less tolerant to noise and distracters than RSDNF. It can be explained by the OR summation sensibility to the number of activated neurons n (see previous section discussions on lateral weights function quality).

B. Hardware utilization

The synthesis and placement was done on a FPGA Xilinx xc6vlx760t-3ff1156, the results on Fig. 9 show that the speed×area product is much better and scalable for the CASAS-DNF model. The speed is not only computed with the maximum clock frequency of the design but also taking into account the computation speed of the model. All these models need time to propagate the spikes information throughout the neural map. This time is expressed in terms of propagation iterations. For CASAS-DNF, the propagation time lasts s + 2res as a bit-stream of size s needs s iterations to be fully generated, and needs (at worse) 2res iteration to go through the whole map (2res is the maximum Manhattan distance for our XY propagation with a resolution res). The update rate is computed with the division of the maximal frequency by the simulated propagation time and the SpeedxArea product is then computed with the division of the number of slice by the update rate.

We are comparing our architecture to RSDNF and to a completely centralized architecture based on the address event representation (AER) that we introduced in [START_REF]Spiking dynamic neural fields architectures on fpga[END_REF].

We can see on Table II that while CASAS-DNF has a slightly slower propagation time than RSDNF, it takes much less hardware resources. However the centralized AER-based implementation is still more competitive in terms of speed and area, though future CASAS-DNF optimizations (by partially sharing PRNs for example) might reverse this conclusion.

VI. CONCLUSION

We presented in this work a new DNF model with the main aim of facilitating its FPGA implementation. Using principles from stochastic computing, the CASAS-DNF model is a cellular array substratum which is able to simulate an all-to-all connectivity between thousands of neurons with a hardware friendly 4-neighbors connectivity.

Neuronal computations are by nature very robust. Stochastic bit stream operators are thus perfectly adapted as their intrinsic noise will not disturb the emergent behavior of the dynamic neural field. Therefore it is not surprising that they represent a very competitive speed-area trade-off for cellular computation of a globally connected neural network. It is also expected that precision and quality of the random numbers are not critical to the behavior. To be fair with the previous implementations, we used very good random numbers quality (one 8bit pseudo random number generation by propagation cycle and by cell routing layer direction: 88200 bits per cycle). But first simulations show that this number can be dramatically reduced. It seems that 5 bits precision is enough and that using the same 5 bits for every SNG of each CASAS cell is enough to maintain an identical behavior. Thus only 6125 pseudo random bits by cycle of bit propagation would be enough reducing the amount of CA-PRNG from a 37x38 grid to a 10x10 grid of 64 cells cellular automaton. Thus if the CA-PRNG are taking now more than 60% of the area a decreasing of the random number quality would reduce the overall area by 50%.

Beyond the improvement for hardware resources utilization compared to the RSDNF model, this architecture is adapted to much more DNF extensions than RSDNF. In this architecture the lateral-weights bit propagation is symmetric in every direction. It means that it will be possible to update neuron potentials on every bit reception instead of waiting the whole lateral-weights propagation time before updating the neurons. With this architecture the pipelined simulation of DNF with synaptic delays will be possible, inducing a highly satisfactory time×area product. Future work will also focus on stronger predictions for lateral-weights quality and on the implementation of learning abilities.

Fig. 2 .

 2 Fig. 2. Difference between the DNF lateral weights function which is a difference of Gaussian (DoG) and the RSDNF lateral weights function which is a difference of exponentials (DoE).

Fig. 3 .

 3 Fig.3. The error between the RSDNF lateral weights function (ωRSDN F ) and the CASAS DNF lateral weights function (ωCASAS) as a function of the neural distance from the activation patch (in number of neurons). Fig.3ashow the error with a rising value for spike probability δ and a fixed number of activated neurons n = 9 and Fig.3bshow the error with a rising number of activated neuron n and with δ = 0.1.

  Fig. 4.Design of the routers and their connectivity. The CASAS cells (white circles) are organized on a 2D-grid. Every router's direction (square) contains one OR gate performing a sum approximation of the inputs stochastic bit-streams. The vertical OR get input from the opposite direction neighbor while horizontal OR get input from opposite direction and vertical directions neighbors. Every direction receive as well an input from the spike's stochastic bit stream generated by the CASAS Cell's neural activation (gray circle). The CASAS map contains two layer of this routers: one for the excitatory bit streams and one for the inhibitory bit streams.

Fig. 5 .Fig. 6 .

 56 Fig.5. Design of the CASAS cell. The 4 neighbors' bit stream and the cell's spike bit stream are routed by the bit stream router (Bs Router) before being multiplied by the synaptic weight generated using a stochastic number generator (SNG). They will then be transmitted to their destination neighbors. The input bit streams will also feed a so called dendrite module (Dent. Exc and Dent Inh.) which will simply sum (with an OR gate) decode and convert the bitstreams and send it to the neuron. The spiking neuron module compute the potential of the neuron. When the potential reaches a threshold it activates the activation module which will stay activated during s cycles. A AND gate will mutliply this activation signal with a stochastic number of probability δ. Note that every SNG modules are fed with pseudo random numbers (PRNIn) generated with a cellular automata.

Fig. 7 .

 7 Fig. 7. Example of the competition scenario. The Inputs maps consist of two Gaussian bubbles rotating clock-wise. The Focus map which represents the neural field activations, tracks one bubble. The error distance between the barycenter of the field activation and the tracked bubble center is displayed on the bottom.

Fig. 8 .

 8 Fig. 8. 50 repetitions of simulation for RSDNF and CASAS-DNF. On 8b we compare RSDNF and CASAS-DNF behaviors for 3 different scenarios. On 8a we show the distance error for rising qualities of stochastic bit-streams.

Fig. 9 .

 9 Fig. 9. Time×area performance for the different models and for different neural map sizes.

TABLE I .

 I SIMULATION PARAMETERS FOR BEHAVIORAL COMPARISON BETWEEN RSDNF AND CASAS-DNF MODELS

	Model	A	B	pe	pi	s
	CASAS-DNF	1.16	0.87	0.035	0.99	100
	RSDNF	1.25	0.70	0.004	0.90	100

TABLE II .

 II RESOURCE UTILIZATION AND PERFORMANCE OF THE ARCHITECTURE ON A XILINX FPGA DEVICE, XC6VLX760T-3FF1156, FOR35X35 NEURONS NEURAL MAP.

	Model	Registers	LUTs	Slices	Freq (Mhz)
	CASAS-DNF	90497	112207	31851	127
	RSDNF	233155	261011	86106	95.7
	AER	74870	121266	46543	54.4
		Prop. it. theory and simu.	Update Rate (Mhz) SpeedxArea
	CASAS-DNF	O(s + 2res)	170	0.36	42635
	RSDNF	O(nN + 2res)	120	0.7975	107970
	AER	n	10	5.44	8556