Claude Touzet
email: claude.touzet@univ-provence.fr

The Illusion of Internal Joy

Keywords: theory of neural cognition, internal joy drive, motivation, consciousness, cortical maps, unsupervised learning, associative memories

Schmidhuber proposes a "theory of fun & intrinsic motivation & creativity" that he has developed over the last two decades. This theory is precise enough to allow the programming of artificial agents exhibiting the requested behaviors. Schmidhuber's theory relies on an explicit 'internal joy drive' implemented by an 'information compression indicator'. In this paper, we show that this indicator is not necessary as soon as the 'brain' implementation involves associative memories, i.e., hierarchical cortical maps. The 'compression factor' is replaced by the 'smallest common activation pattern' in our framework, with the advantage of an immediate and plausible neural implementation. Our conclusion states that the 'internal joy' is an illusion. This remind us of the eliminative materialism position which claims that 'free-will' is also an illusion.

Introduction

J. Schmidhuber build his "theory of fun & intrinsic motivation & creativity"

 [START_REF] Schmidhuber | Formal Theory of Creativity, Fun, and Intrinsic Motivation[END_REF]on the maximization of an 'internal joy' that drives a reinforcement learning process. He proposes an operational description of it, a necessary step in order to provide an artifact with fun, motivation and creativity. The intrinsic reward is computed as the compression progress expressed as the number of saved bits [START_REF] Schmidhuber | Driven by compression progress: A simple principle explains essential aspects of subjective beauty, novelty, surprise, interestingness, attention, curiosity, creativity, art, science[END_REF]. A number of systems have been build following this recommendation, which exhibit the desired creativity behavior [START_REF] Schmidhuber | Simple Algorithmic Principles of Discovery, Subjective Beauty, Selective Attention, Curiosity & Creativity[END_REF].

We have to be aware that, as soon as something such as 'intrinsic motivation' is defined, reinforcement learning becomes de facto the unique valid candidate for the learning process. Reinforcement learning is certainly a very efficient way to acquire behaviors [START_REF] Santos | Exploration Tuned Reinforcement Function[END_REF][START_REF] Touzet | Q-learning for robots[END_REF], but it is not the only one. Supervised learning and self-organization do exist. They are not considered as valid candidates for the learning of 'creativity' behaviors because:

-in the case of supervised learning, its 'supervision' would limit the freedom that we think is necessarily involved by 'creativity', -in the case of self-organization, its (self-)'organization' only allows to represent the data, and therefore (again) lacks the ability or freedom involved in 'creativity'.

Both opinions are misplaced. There is no link between the learning process and the ability to escape the learned knowledge. 'Generalization' is the ability to process correctly new unknown inputs (following the 'rules' extracted from previously learned knowledge). The generalization quality depends only on the learning samples (not the learning process) AND the implementation.

It is our goal to show here that the 'generalization' associated to a hierarchical cortical maps implementation is able to create any behavior involving 'intrinsic motivation'. If successful, our demonstration will also state that any learning (not only reinforcement learning) may contribute to 'creativity' behaviors, and that 'intrinsic motivation', also referred sometimes as 'internal joy', is an illusion.

Cortex organization

The cortex is a hierarchy of cortical maps, each cortical map acting as a selforganizing associative memory that preserves the topology and distribution of the input data [START_REF] Kohonen | Self-Organizing Maps[END_REF]. Each cortical map acts as a novelty filter, and stops any known situation (or part of situation). Only new unknown situation (or part of it) is allowed to proceed along the hierarchy towards maps of higher level of abstraction (Fig. 1). Behaviors are automatically generated through the cooperation of pairs of cortical maps.

Fig. 1. The cortex is organized as a hierarchy of cortical maps, from the primary cortex receiving sensory inputs, to the secondary cortex allowing for inputs fusion, to the associative cortex where map encoding is sensory modality independent. The cortical maps are novelty filters: as the information progresses in the hierarchy, it is stopped as soon as it is recognized. only uncommon (extra-ordinary) situations (for a given individual) reach the 'abstract' levels that account for the 'goals'.

Behaviors are goal-directed

A behavior is a sequence of actions that can be related to the same goal. The goal is a specific situation, which will end the behavior as soon as it is reached, after what a Secondary cortex

Primary cortex

Extra-ordinary situations

Very frequent situations

Ordinary situations

Associative cortex new behavior starts (with a new goal). During a given behavior, each action is chosen in order to reduce the distance between the present situation and the goal-situation. Within our framework, a situation is not a x-y-z vector coding for a given location in the real world, but a vector in the cortical maps space. The cortex is build of several hundreds cortical maps, each one devoted to a specific category of information. The cortical maps pattern of activations at any given time is representative of the input (i.e., experienced) situation combined with the memorized lifelong experience encoded in the synapses. A goal (situation) is therefore a multidimensional vector sharing similarities with the experienced situation.

Coordination of pairs of cortical maps

The work described in [START_REF] Touzet | Modeling and Simulation of Elementary Robot Behaviors using Associative Memories[END_REF] explains how two associative memories (cortical maps) cooperate in order to produce goal-directed behaviors (Fig. 2). To resume (on map 1), a cluster of activity -labeled 'goal' -is defined by some external (or internal) inputs. On the same map, the experienced current situation is represented by the fact that it activates a particular cluster of cortical columns. If there is a difference between both activities, then both activities (experienced situation and goal) will help activate cortical columns neighbors of the experienced situation. These neighbors columns encode for a situation that is close, but nevertheless different, from the current one. The difference between the two (neighbor) activities selects cortical columns on associated cortical maps (map 2). These cortical columns have memorized the actions associated to any variation of activity (of map 1). This action is done in the real world, and put the agent in a new situation (that should be closer to the goal), and everything starts again. Fig. 2. Two cortical maps cooperate in order to generate a sequence of actions (i.e., a behavior) in response to experienced situations. Each intersection represents a cortical column. The topology of the input situation space is preserved by the first cortical map [START_REF] Schmidhuber | Formal Theory of Creativity, Fun, and Intrinsic Motivation[END_REF]. Therefore, the situation neighbor to the input situation and closer to the goal situation (a) is an intermediate situation in the process devoted to reach this 'goal'. The variation of activity between experienced and intermediate situations (b) serves as input to the second map [START_REF] Schmidhuber | Driven by compression progress: A simple principle explains essential aspects of subjective beauty, novelty, surprise, interestingness, attention, curiosity, creativity, art, science[END_REF], where it activates the muscle command (action) associated (c) to this variation of situations.

Variation of situations

Experienced situation

Goal situation (defined by external inputs)

Illusions

The definition of an illusion that is of interest here is "a misinterpretation of a true sensation". We claim that the 'internal joy drive' and 'intrinsic joy' are both misinterpretations.

The illusion of 'joy driven' behaviors

When a behavior's goal is not seen or perceived or understood by the observer (who can be the individual himself), then the behavior is said to be 'driven by joy'. Why 'joy'? Just because it is easier to believe that the individual is looking for, or responding to, something pleasant when behaving -instead of the opposite (searching for bad things). Using the information provided in section 2, we claim that any behavior is 'goal driven' and that the 'internal progress' is in fact a 'distance-to-goal progress'. As long as there is a discrepancy between the goal and the experienced situation locations on the cortical map(s), the behavior will occur.

The illusion of 'internal joy'

If we assume that 'joy driven' behaviors are just 'goal driven' behaviors whose goals are non explicit (to a human observer), then it follows that we must get better acquainted with these goals. Moreover, since they are related to 'joy', it is of tremendous importance for our well-being to know more about them.

First things first: how are the goals defined and selected? Cortical maps are associative memories (i.e., content addressable memories). Hebbian learning induces the storage of the activity patterns generated by the (lifelong) experienced situations. The number of samples required to (self-)organize a map is several times the number of cortical columns of the map. Therefore, only the most represented situations (among the samples) are going to be memorized. For each subset of situations, the memorized information is the most redundant one, i.e., the most shared cortical column activities: the 'smallest common activation pattern'.

Next thing on the list: how does a 'joy driven' behavior start? Let's imagine that the individual is in a situation where there is no explicit 'goal'. There is nevertheless a situation to experience (even if it is sensory deprivation). This experienced situation will activate a representation of it by activating (some) cortical columns on some cortical maps. This activity acts as a probe of the associative memories and will sooner or later activate an already stored activity pattern: a 'goal'. Now, the system (i.e., the cortex) is faced with two different patterns of activity: one representing the experienced situation and one the 'goal'. A behavior will emerge.

It may happened from time to time that the experienced novel situation matches exactly an already memorized -but still never experienced (otherwise it would have been stopped shortly after the primary or secondary cortex) -pattern of activity (Fig. 3). In this case, the neural pattern of activity is minimum, which allows for a much better memorization of the event (in the episodic memory). Minimum activity means that the number of cortical columns involved is minimum -but their electric activity is maximum! It follows that this particular event will become easily remembered and may often serve as a goal (i.e., an attractor situation) in the following life of the individual. Less neural activity will be required to activate this particular representation. Immediate recognition and better memorization: all the ingredients of a very meaningful experience that could be named 'joy' or 'beauty'.

Discussion & Conclusion

We must point out the similarities between Schmidhuber's compression index and the 'smallest common activation pattern'. Both consider that the size of the representation of information by the brain is the important factor. Both also acknowledge the fact that as soon as there is no more any variation of this 'factor', then the behavior stops.

Looking at the biological plausibility of each proposal, a compression index needs to be computed and therefore requires resources. On the opposite, there is no computational resource required in our proposal.

Last but not least, our proposal allows for any learning (supervised, reinforcement and seal-organization) to occur, not only reinforcement learning.

'Curiosity' and 'boredom' are not equivalent to 'internal joy' even if they can be implemented using a reinforcement signal computed by the difference between the expected situation and the experienced situation [START_REF] Schmidhuber | A possibility for implementing curiosity and boredom in model-building neural controllers[END_REF]. To us, the discrepancy between expected and experienced situations is the root cause of the attentional processes, with the involvement of more and more higher (or abstract) cortical maps as the discrepancy hold [START_REF] Touzet | Conscience, intelligence, libre-arbitre ? Les réponses de la Théorie neuronale de la Cognition[END_REF]. In our view of the attentional processes, there is no need for a reinforcement signal computed elsewhere to achieve 'attention', and therefore also 'curiosity' and 'boredom'.

The goal situations should be considered as attraction basins since they act in twisting the behavior of the individual towards specific representations (the attractive ones). A naive observer will see an individual whose behaviors tend to favor, even research, (a) (b) some specific situations. It follows that 'internal joy' is only a side-effect of the brain learning and memory processes. These explanations define 'internal joy' as an illusion, as it is also the case of 'consciousness' and 'intelligence' in the eliminative materialism paradigm. We may feel as if something (joy) was lost in the argument, but the gain is enormous. With a clear explanation of the 'internal joy' illusion, it can now be efficiently looked-for by any individual. In particular, since goal situations depend only on the individual experiences, then the illusions of 'joy' (or 'beauty') and 'joy drive' (or 'motivation') are individual dependent. There is no need to covet what makes somebody else happy: our happiness lies in the (brain-)crystallization of our life, and will evolved as time passes in order to account for any variation in our daily experience. With a careful design of our daily experiences, we can change what makes us happy, and definitely improve our odds in an happiness quest.

Regarding in silico implementation of (human) cognition, we feel that it is definitely within our reach. The burden has been moved from the mysterious 'internal joy' definition to the more accessible sample learning. The condition to respect in order to build a human-like cognition is that the learning base must only involve human-like situations.

  and associated action map[START_REF] Schmidhuber | Driven by compression progress: A simple principle explains essential aspects of subjective beauty, novelty, surprise, interestingness, attention, curiosity, creativity, art, science[END_REF]

Fig. 3 .

 3 Fig. 3. (a) Columns activity (black triangles) associated to an experienced extra-ordinary situation, that quickly resume in an activity (b) involving less columns (each one exhibiting stronger activity, i.e., bigger triangles). The pattern (a) is dependent on the life-long learning. It contains in essence the pattern (b), which was only waiting for this extra-ordinary experience in order to emerge. The memorization and future activation of pattern (b) is much easier than was pattern (a) activation. The quick transition from (a) to (b) is certainly an unique experience, that an observer may called 'joy' .


Claude F Touzet

Q-learning for Robots

published or not. The documents may come

Introduction

Robot learning is a challenging -and somewhat uniqueresearch domain. If a robot behavior is defined as a mapping between situations that occurred in the real world and actions to be accomplished, then the supervised learning of a robot behavior requires a set of representative examples (situation, desired action). In order to be able to gather such learning base, the human operator must have a deep understanding of the robot-world interaction (i.e., a model). But, there are many application domains where such models cannot be obtained, either because detailed knowledge of the robot's world is unavailable (e.g., spatial or underwater exploration, nuclear or toxic waste management), or because it would be to costly. In this context, the automatic synthesis of a representative learning base is an important issue. It can be sought using reinforcement learning techniques -in particular Q-learning which does not require a model of the robot-world interaction. Compared to supervised learning, Q-learning examples are triplets (situation, action, Q value), where the Q value is the utility of executing the action in the situation. The supervised learning base is obtained by recruiting the triplets with the highest utility.

Because it allows the synthesis of behaviors despite the absence of a robot-world interaction model, Q-learning [START_REF] Watkins | Learning from Delayed Rewards[END_REF]) has become the most used learning algorithm for autonomous robotics. Although the convergence theorem does not apply to the robotics domain (due to the limited number of situation-action pairs that can be explored during the life-time of the robot batteries), heuristically adapted Qlearning has proved successful in applications such as obstacle avoidance, wall following, go-to-the-nest, etc. This is mostly due to neural-based implementations such as multilayer perceptrons trained with backpropagation, or self-organizing maps. Such implementations provide an efficient generalization, i.e., fast learning, and designate the critic -the reinforcement function definition -as the real issue. The articles REINFORCEMENT LEARNING and REINFORCEMENT LEARNING IN MOTOR CONTROL provide background information on reinforcement learning. [START_REF] Kaelbling | Reinforcement Learning: A Survey[END_REF] and [START_REF] Sutton | Reinforcement Learning[END_REF] are two other sources of information. For more detailed treatments, the reader should consult [START_REF] Touzet | Neural Reinforcement Learning for Behaviour Synthesis[END_REF].

Q-learning

Figure 1 shows a functional decomposition of Q-learning. Three different functions are involved: evaluation, memorization and updating. Using the information stored in the robot memory, the current situation is evaluated to select the best action to accomplish (i.e., the most reward-promising action). This proposition is modified so as to allow exploration of the situation-action space. The new situation, entered as a consequence of the execution of the action, is qualified by the reinforcement function. Its qualitative criterion (reinforcement) is used by the updating algorithm to adjust the Q values in the following way:

Q(s,a) new = Q(s,a) old + ß(r + γ.Max(Q(s',a)) -Q(s,a) old ) (1)
where s is the situation, a is the action, r is the reinforcement and s' represents all situations that can be reached from s. ß and γ are positive coefficients less than 1.

Convergence

Recent developments in the theory of reinforcement learning have allowed proof of asymptotic convergence [START_REF] Dayan | TD(λ) convergences with probability 1[END_REF]). These proofs rely on several assumptions that do not apply to robots facing real-world tasks. In particular, the asymptotic convergence requires a discrete coding of the situation-action pairs (look-up table storage) and to try out every action for every situation an infinite number of times. A robot is a mechanical device than needs at least a few hundred microseconds to execute any action. Therefore, due to batterylife that typically last less than 10 hours, only a few thousand situation-action can be visited during a given experiment. This is an extremely small number when compared to the potential number of situation-action pairs (e.g., 10 26 for the Khepera miniature mobile robot, today the most common research robot). Thus, generalization between similar situation-action pairs is mandatory. This action is the one that has the best probability of reward. However, this proposition may be modified by the evaluation function to allow an extensive exploration of the situation-action space. After the execution of the action by the robot in the real world, a reinforcement function provides a reinforcement value. This value -a simple qualitative criterion (e.g., +1, -1 or 0) -is used by the updating algorithm to adjust the utility value associated with the situation-action pair.

Generalization

Improvements emphasizing generalization have been proposed by [START_REF] Mahadevan | Automatic Programming of Behavior-based Robots using Reinforcement Learning[END_REF] who use weighted Hamming distance to generalize between similar situations. This simple method is limited to syntactic situation criteria (i.e., it is dependent on the coding of the situations). A second method, proposed by the same authors, adds the action into the syntactic criteria, using clusters to generalize across similar situation-action sets. One of the problems is that the clusters have to be handpicked.

Neural Q-learning

Neural implementations offer a compact representation (i.e., limited memory requirement) and good generalization performance (as demonstrated by numerous connectionist applications). The memorization function uses the weight set of the neural network: the memory size required by the system to store the knowledge is defined, a priori, by the number of connections in the network. It is independent of the number of explored situation-action pairs. The proposed action is the processing result of the situation by the network, plus the addition of a random component for the exploration. The update function uses the weight modification algorithm to store the utility values computed by the Q-learning rule (1). The ideal neural implementation would provide, in a given situation, the best action to undertake and its associated Q value. However, training of a such network requires the definition of an error on the output layer, i.e., knowledge of the best action to undertake in every situation. Such knowledge can be inferred if there is only two different possible actions for the robot as in the cart pole balancing problem [START_REF] Barto | Neuron like elements that can solve difficult learning control problems[END_REF]. However, in the general case, the number of actions is larger. [START_REF] Lin | Reinforcement Learning for Robots Using Neural Networks[END_REF] -who proposed the first multilayer perceptron implementation of the Q-learning (Q-Con) -uses as many perceptrons as there are actions, each network output coding for the utility of accomplishing this action in the current situation. Therefore, only one Q-Con network is updated at every time step, and generalization between networks (i.e., actions) is impossible. Other multilayer perceptron implementations have been proposed [START_REF] Ackley | Interactions Between Learning and Evolution[END_REF]Littman 1991, Touzet 1997), but they do not yet solve the output error definition problem.

Q-Kohon

Unsupervised learning models -such as SELF-ORGANIZING FEATURE MAPS: KOHONEN MAPS (q.v.), RADIAL BASIS FUNCTION NETWORKS (q.v.), ADAPTIVE RESONANCE THEORY (ART) (q.v.), CMACdo not require an error definition for updating their weight values. Q-Kohon, a Kohonen map implementation of the Qlearning, is a method of state grouping involving syntactic similarity and locality [START_REF] Mccallum | Instance-based State Identification for Reinforcement Learning[END_REF]. Each neuron codes a particular triplet (situation, action, Q value); therefore the number of neurons equals the number of stored associations. The neighborhood property of the self-organizing map accounts for the generalization across similar situation-action pairs.

Q-Kohon uses the self-organizing map as an associative memory. This associative memory stored triplets. Part of a triplet is used to probe the self-organizing map in search of the corresponding information. Here, situation and Q value are used to find the action: the best action to undertake in a world situation is given by the neuron that has the minimal distance to the input situation and to a Q value of value +1. The selected neuron corresponds to a triplet (situation, action, Q ). It is this particular action that should offer the best reward in the world situation. To update the Q value, equation (1) requires the maximum Q value of the new entered situation. This is easily obtained by probing the map with the new situation and a Q value of value + 1. The selected neuron Q value will be the maximum possible value.

A nice side effect of using clustering techniques to implement the Q-learning is that the learned behavior can be interpreted by looking (see SELF-ORGANIZING FEATURE MAPS: KOHONEN MAPS) at the network weights (something extremely difficult with multilayer implementations). Also, because the neurons of the selforganizing map approximate the probability density function of the inputs, one can predict that if a correct behavior is learned, all neurons will code positive Q values. This is most useful to determine when a correct behavior has been learned. This last fact results in the optimization of the stored knowledge.

Q-learning + Hamming + clustering Q-Comp Q-Kohon

Comparisons

Experiments aimed at comparing various implementations of the Q-learning in a task of synthesizing an obstacle avoidance behavior for the miniature robot Khepera [START_REF] Touzet | Neural Reinforcement Learning for Behaviour Synthesis[END_REF] demonstrate that neural Q-learning implementations require alot less memory, less learning examples and learn faster (cf. table 1). The Q-Kohon implementation also exhibits the best behavior after learning, i.e., less negative reinforcements received than all the other implementations.

Reinforcement Function Design

The reinforcement function quality is intrinsically limited by the expert's abilities. When a reinforcement learning experiment does not converge, it is impossible to know if this is due to the fact that the experiment was too short and more examples are needed, or if the intrinsic nature of the reinforcement function forbids convergence. Today, reinforcement learning researchers use a slow-and-painful trial and error approach to define the reinforcement function. In the meantime, efforts have been devoted to find ways to automatically define such functions. [START_REF] Santos | Dynamic Update of the Reinforcement Function during Learning, Connection Science, Special issue on Adaptive Robots[END_REF] have proposed an Update Parameter Algorithm (UPA) to automatically adjust the threshold values: θ + and θ -within a particular definition of the reinforcement function:

RF(s 1 ,...,s u ) = +1 if g 1 (s 1 ,...,s u ) > θ + -1 if g 2 (s 1 ,...,s u ) < θ - 0 otherwise     
where (s 1 ,..., s u ) is the output readings of the sensors, g 1 ( ) and g 2 ( ) are any functions linking the sensor data to the rewards.

The resulting effect is to optimize the exploration part of the learning phase by achieving and maintaining pre-defined ratios of positive and negative rewards. If there is no positive reward, the evaluation function built during the learning phase will have "0" as maximum value and the policy cannot select effective actions. If there is no negative reward, the robot can remain in a dead-end situation forever. If there is no null reward, the evaluation function will be non-continuous at the frontier between positive and negative situation-action pairs.

A dynamic version of UPA [START_REF] Santos | Dynamic Update of the Reinforcement Function during Learning, Connection Science, Special issue on Adaptive Robots[END_REF] updates the threshold values during the learning phase -exploration and exploitation (so as to take into account the improvement of the robot policy). It allows behavior performance improvements without the need of some sort of external supervisor, capable of ranking situations by difficulty and of choosing tasks of increasing difficulty [START_REF] Dorigo | Robot Shaping: An Experiment in Behavior Engineering[END_REF] Another major issue is to be able to overcome the exponentially growing number of required learning examples that comes with target behaviors of greater complexity. Battery-life time seems to impose a definite limit, and researchers tend to promote knowledge incorporation as a speed-up mechanism. The goal is to bias the exploration towards "rewarding" part of the search space -at the expense of tabula rasa methods. The drawback is that newunforeseen -solutions cannot be discovered. Lazy learning (Aha 1997), also called instance-based learning, provides a way to add samples without implying bias. In a lazy learning approach, the computation of the inputs is delayed until the necessity arises. Lazy learning samples the situation-action space, storing the succession of events in memory and, when needed, probes the associative memory for the best move. The sampling process stores the successive situation-action pairs generated by a random action selection policy. The exploration phase is done only once, stored and used later by all future experiments. The probing of the memory involves complicated computations: clustering, pattern matching, and so forth.

By storing situation-action pairs, a lazy memory builds an non-explicit model of the situation transition function, that is used as a bias to leverage the model-free following learning phase (i.e., Q-learning). [START_REF] Sheppard | A Teaching Strategy for Memory-Based Control[END_REF] propose to mix lazy learning and reinforcement learning, probing the memory with the reinforcement function. Their objective is to provide a method for predicting the rewards for some state-action pairs without explicitly generating them. They call their algorithm lazy Q-learning. For the current real world situation, a situation matcher locates all the states in the memory that are within a given distance. If the situation matcher has failed to find any nearby situations, the action comparator selects an action at random. Otherwise, the action comparator examines the expected rewards associated with each of these situations and selects the action with the highest expected reward. This action is then executed, resulting in a new situation. There is a fixed probability of generating a random action regardless of the outcome of the situation matcher. New situation-action pairs are added to the memory, along with their Q values computed in the classical way. Among similar situation-action pairs in the memory, an update of the stored Q values is made. There is a limit to the genericity of this lazy memory because the Q values associated with the situation-action pairs only apply for a particular application.

Learning is not restricted to single robots. Learning in cooperative robotics promises are beguiling: a way to program a set of robots without having to explicitly model their interactions with the world -including the other team members -to achieve cooperation. To achieve this goal, mechanisms that relay the unique information associated with the team behavior (reinforcement value) to the individual robots have to be found. Results from the multi-agent research community cannot be applied since they are usually symbolic methods, where robot Q-learning requires a sub-symbolic approach.

Despite all the efforts and success around Q-learning, there are several drawbacks associated with supervised and reinforcement learning when in comes to real applications. First, the time needed to achieve the synthesis of any behavior is prohibitive. Second, the robot behavior during the learning phase is -by definition -bad, it may even be dangerous. Third, except within the lazy learning approach, a new behavior implies a new learning phase. What is needed is a learning that instantaneously synthesizes any behavior, and which performance improves by the mere repetition of this behavior (for further discussion see [START_REF] Touzet | Programming Robots with Associative Memories[END_REF]).

Figure 1 .

 1 Figure1. Q-learning method functional decomposition. In response to the present situation, an action is proposed by the robot memory. This action is the one that has the best probability of reward. However, this proposition may be modified by the evaluation function to allow an extensive exploration of the situation-action space. After the execution of the action by the robot in the real world, a reinforcement function provides a reinforcement value. This value -a simple qualitative criterion (e.g., +1, -1 or 0) -is used by the updating algorithm to adjust the utility value associated with the situation-action pair.

  . Santos et al. have been able to synthesize -for the first time -a wall following behavior by reinforcement learning (cf. figure 2), a demonstrative support for reinforcement function design techniques. Discussion Q-learning is the most used (reinforcement) learning technique for behavior-based robots (see REACTIVE ROBOTIC SYSTEMS). Neural based Q-learning implementations provide compactness and generalization. Clustering-based neural based methods, such as Q-Kohon, allow drastic reduction of the learning time and number of examples required. Their efficiency put forward the reinforcement function definition as one of the major issues.

Figure 2 .

 2 Figure 2. The trace of the miniature Khepera robot after the synthesis of a wall following behavior using a RBF implementation of the Q-learning and Dynamic-UPA in a new environment. Only about 2000 learning iterations are needed.

Table 1 .

 1 Comparison

	Time length	55 mn	25 mn	30 mn	8 mn	2 mn
	# iterations	7500	3500	4000	2000	500
	Memory size	6400	6400	1.6 10 6	56	176

of various implementations. The learning time is the time in seconds needed to synthesize an obstacle avoidance behavior. It reflects the number of real world experiments required. The number of learning iterations is the number of updates to the memory (look-up table or neural network). The memory size is the number of floats required to store the information.


Paul Rochet

Isabel Serra

The Mean/Max Statistic in Extreme Value Analysis

Keywords: Generalized Pareto model, tail index, hypothesis testing MSC (2000): 62G32, 62F03, 62G15

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

Two fundamental results marked the starting point of Extreme Value Theory (EVT): Fisher-Tippet-Gnedenko and Pickands-Balkemax-de Haan Theorems. These theorems enabled to characterize the asymptotic behavior of extreme values by a real number k, called tail index (we refer to [START_REF] Markovich | Nonparametric analysis of univariate heavy-tailed data : research and practice[END_REF], [2], [START_REF] Coles | An introduction to statistical modeling of extreme values[END_REF][START_REF] Embrechts | Modelling extremal events: for insurance and finance[END_REF][START_REF] Finkenstadt | Extreme values in finance, telecommunications, and the environment[END_REF] and [START_REF] Alexander | Quantitative risk management: concepts, techniques, and tools[END_REF] for a survey). Ever since Gnedenko's original paper, estimation of the tail index has become a major concern in the literature, for which even the most used techniques can produce unsatisfactory results in many situations. For instance, several authors have highlighted the estimation problems that arise from the second approach, see [START_REF] Castillo | Fitting the generalized pareto distribution to data[END_REF][START_REF] Jonathan | Parameter and quantile estimation for the generalized pareto distribution[END_REF][START_REF] Zhang | A new and efficient estimation method for the generalized pareto distribution[END_REF].

Pickands-Balkemax-de Haan Theorem has widened the use of the Generalized Pareto Distribution (GPD) in extreme values theory as a model for tails. According to the behavior of the probability density functions of the GPD, we can distinguish three submodels corresponding to k < 0, 0 < k < 1 and k > 1 respectively, separated by the exponential distribution (k = 0) and the uniform distribution (k = 1). Since no method is able to provide a satisfying estimation of the tail index regardless of the true model (see [START_REF] Del | Likelihood inference for generalized pareto distribution[END_REF]), the parameter inference should be made subsequently to the choice of submodel, even more so for small samples. The first recommendation is to specify the submodel of the GPD a priori, in order to estimate the parameters. Therefore, we propose to use a test procedure to distinguish between the two transitional phases of the tail index in GPD models: the exponential and uniform distributions. Due to the monotonic behavior of the test statistic with respect to the tail index of the GPD, the test can be used for classification purposes.

Tests on separate families of hypotheses have been originally studied in [START_REF] David R Cox | Tests of separate families of hypotheses[END_REF] who proposed a generalization of the Neyman-Pearson principle, based on a likelihood ratio criterion. The so-called maximumlikelihood ratio test is not necessarily optimal, although it provides a practical methodology to test a wide range of hypotheses. Examples studied in the literature include invariant tests of exponential versus normal or uniform distributions in [START_REF] Vincent | An optimum test property of two weil-known statistics[END_REF], normal versus Cauchy distribution in [START_REF] Franck | The most powerful invariant test of normal versus cauchy with applications to stable alternatives[END_REF] and so on.

In this paper, we show that Cox's maximum likelihood ratio test for uniform versus exponential distributions is the most powerful among scale-free tests. It relies on the ratio τ n between the empirical mean and sample maximum, whose behavior is highly dictated by the tail of the distribution, thus making it a main feature of the GPD framework. The mean/max statistic τ n then provides a simple yet effective way to detect the tail behavior in a Generalized Pareto model. In fact, the tail classification method based solely on the value of τ n produces extremely conclusive results that compare favorably with the more computationally expensive inference methods such as Zhang and Stephens [START_REF] Zhang | A new and efficient estimation method for the generalized pareto distribution[END_REF], Song and Song [START_REF] Song | A quantile estimation for massive data with generalized pareto distribution[END_REF] or maximum likelihood.

The problem of bad classification in GPD estimation was recently detected [START_REF] Del | Likelihood inference for generalized pareto distribution[END_REF] but no satisfactory solution has been proposed so far. Summarizing, the quality of inference methods for the tail index k in the GPD highly depends on the underlying submode, whether it is Model A (k > 1), Model B (0 ≤ k ≤ 1) or Model C (k < 0). Because standard methods tend to specialize to certain regions of the values of the parameter, it is crucial to separate the GPD model into these three submodels reflecting the different behavior of the tails.

The test on separate families of hypotheses to distinguish between uniform and exponential distributions is described in Section 2. In Section 3, we show how to extend the test procedure for tail classification and detection in Generalized Pareto models or more general distributions. In Section 4, we propose a simple general recipe to infer on the tail of the distribution of a sample and an application on global seismic activity data is presented. Theoretical results regarding the distribution of the mean/max statistic are discussed in the Appendix.

Test to distinguish between uniform and exponential tail distributions

Let X = (X 1 , . . . , X n ) be a sample of independent identically distributed variables drawn from a distribution on [0, +∞) assumed to be either uniform on some interval [0, θ] with θ > 0 or exponential with parameter λ > 0. As it is usually the case when dealing with separate families of hypotheses, a uniformly most powerful test does not exist because the likelihood ratio statistic depends on the true value of the parameter. The maximum likelihood ratio statistic was proposed in [START_REF] David R Cox | Tests of separate families of hypotheses[END_REF] as a way to generalize Neyman-Pearson's principle to families of distributions, although the optimality of the test is no longer backed up by the theory. Nevertheless, we show that in this situation, Cox's maximum likelihood ratio test is the most powerful among a large class of tests whose distribution is invariant within each hypothesis. For such tests, both the level and power are determined by the common distributions of the statistic under the null hypothesis and the alternative. The unicity of these distributions within each hypothesis enables to deduce the most powerful test from a direct application of Neyman-Pearson's lemma.

We are interested in testing the null hypothesis H 0 : "The X i 's are uniformly distributed on an interval [0, θ] with θ > 0" against the alternative H 1 : "The X i 's are exponentially distributed."

In this particular situation, one can exploit the fact that both the uniform and exponential models are closed by positive scaling: if the sample X belongs to the uniform (resp. exponential) model, then so does tX for all t > 0. We consider the class of scale-free tests, that is binary valued functions Φ = Φ(X) of the sample X satisfying Φ(X) = Φ(tX) , ∀t > 0.

If we let Φ(X) = 1{X ∈ R} for R a Borel set of R n , saying that Φ is scale-free simply means that R is a cone. Equivalently, a statistic is scale-free if it can be expressed as a function of the normalized sample X/ X , for any norm . on R n .

Because the uniform and exponential models are stable by positive scaling, the distribution of a scalefree statistic does not depend on the parameter, whether we are under the null hypothesis (uniform) or the alternative (exponential). Therefore, we know there exists a most powerful scale-free test for these hypotheses, which we can derive from the likelihood ratio statistic of the normalized sample. As we show below, the resulting test is function of the statistic

τ n := X n X (n)
where X (n) = max{X 1 , . . . , X n } and X n = 1 n n i=1 X i , whose properties are often discussed in relations with the tail distribution, see [START_REF] José | A bootstrap goodness of fit test for the generalized pareto distribution[END_REF].

Theorem 2.1. The test Φ = 1{τ n ≤ c α } with c α such that P(τ n ≤ c α |H 0 ) = α is the most powerful scale-free test of level α ∈ (0, 1) to test H 0 against H 1 .

Proof. Let Φ = Φ(X) be a scale-free test of level α ∈ (0, 1), we can write almost surely (for

X n = 0) Φ(X 1 , . . . , X n ) = Φ X 1 X n , . . . , X n-1 X n , 1 := Ψ X 1 X n , . . . , X n-1 X n , for Ψ : R n-1 → {0, 1}. Let Y = (Y 1 , . . . , Y n-1 ) = (X 1 /X n , . . . , X n-1 /X n ). Since
L 0 (Y ) = 1 n(max{1, Y (n-1) }) n and L 1 (Y ) = (n -1)! (1 + n-1 i=1 Y i ) n .
Recall that Y i = X i /X n , Neyman-Pearson's likelihood ratio is thus given by

1 n! (1 + n-1 i=1 Y i ) n (max{1, Y (n-1) }) n = (nτ n ) n n! .
Since it is an increasing function of τ n , the likelihood ratio test writes as 1{τ n ≤ c α } for some suitable c α .

In this framework, the most powerful scale-free test actually recovers Cox's maximum likelihood ratio test, introduced in [START_REF] David R Cox | Tests of separate families of hypotheses[END_REF] in a more general context. Indeed, considering the likelihood functions

L 0 (θ, X) = θ -n 1{X i ≤ θ, i = 1, . . . , n} and L 1 (λ, X) = λ n exp(-λ n i=1 X i ), the maximum likelihood ratio statistic is given by sup λ>0 L 1 (λ, X) sup θ>0 L 0 (θ, X) = X (n) e X n n = 1 e τ n n .
Thus, Cox's maximum likelihood ratio test rejects the null hypothesis H 0 for sufficiently small values of τ n . The threshold c α corresponding to a test of significance level α can be computed easily by Monte-Carlo. Nevertheless, to get an analytical expression of the threshold and the corresponding power of the test, one needs to derive the true distribution of the statistic τ n under the null hypothesis and under the alternative. This issue is discussed in the Appendix.

Tail classification in Generalized Pareto distribution

If there is no particular reason to favor the uniform model over the exponential one, the significance level can be chosen so that the probabilities of error under H 0 and H 1 are equal, thus inducing a minimal probability of error in the worst case scenario. In this purpose, one must choose the threshold c n as the unique solution of

P(τ n ≥ c n |H 0 ) = P(τ n ≤ c n |H 1 ) := A(n).
Thus, when using the threshold c n , the probability A(n) of selecting the correct model between uniform and exponential no longer depends on the true distribution. The values of the threshold c n and percentage of accurate selections A(n) are computed in Of course, the level of accuracy is only exact if the model is well-specified, which is rarely the case in practice. Nevertheless, the test procedure may provide substantial information on the tail of the distribution. To study the behavior of extreme events, one focuses generally on the tail observations, i.e. the data over a certain threshold µ, so as to obtain independent realizations X i conditionally to X i > µ. Under regularity conditions, the tail observations (translated to the left by a factor µ so as to have their support starting at zero) converge in distribution as µ → ∞ towards independent realizations of a Generalized Pareto Distribution (GPD), with density

f k,σ (x) = 1 σ 1 -k x σ 1-k k , x ≥ 0 if k ≤ 0 0 ≤ x ≤ σ k if k > 0.
The scale σ > 0 and the shape k ∈ R (also called the tail index) are the two parameters used to describe the tail. The Generalized Pareto density can be extended by continuity to the values k = 0 and k = 1 recovering respectively the exponential and uniform distributions. Generally, the relevant information on the tail comes from the scale parameter k and less importance is given to the scale σ.

When working with tail observations, it is thus customary to assume that the data are drawn from a GPD. In this model, the behavior of the mean/max statistic τ n is highly dictated by the tail index k while its distribution does not depend on the scale parameter σ (τ n is scale-free). The mean/max statistic is thus perfectly suited to investigate the tail index of the GPD without having to concern about the scale. As a matter of fact, a crucial matter in Generalized Pareto models is to determine to what submodel the tail belongs, i.e. if the tail index k is greater than 1 (model A), between 0 and 1 (model B) or negative (model C). Because the frontiers between the different models are achieved for the exponential (k = 0) and uniform (k = 1) cases, the test procedure to distinguish between uniform and exponential distributions can be naturally extended for classification purposes.

The asymptotic distribution of τ n (see Theorem A.3 in the Appendix) reveals that τ n vanishes at a rate of n -k when k ∈ (0, 1), converges to k/(k + 1) when k > 0 and is of the order 1/ log(n) in-between, for k = 0. In particular, the median of τ n can be deduced from Theorem A.3,

med(τ n ) =                      - k k + 1 log(2) n -k + o(n k ) for -1 < k < 0 1 log(n) + log log(2) log 2 (n) + o 1 log 2 (n) for k = 0 k k + 1 + o 1 √ n for k > 0. (1)
The monotonic behavior of τ n with respect to the tail index k makes it a usefool tool to learn to which submodel the distribution of the data belongs. The idea is simple: the practitioner chooses two reals numbers 0 < a n < b n < 1 and concludes to the model A if

τ n ≥ b n , the model B if a n < τ n < b n and the model C if τ n ≤ a n .
As suitable values of a n and b n , we use the theoretical median of τ n under the exponential and uniform distributions respectively. By taking these values, the probability of selecting each model in the transitional cases k = 0 and k = 1 does not exceed 1/2 so that none of the model A, B or C is favored. The actual value of med(τ n ) in the uniform case follows directly from Lemma A.1. Because we could not obtain an analytical expression of med(τ n ) in the exponential model, we use the asymptotic approximation given in Equation (1), which is in fact extremely accurate, even for small sample sizes. Thus, the bounds a n and b n used for the classification are

a n = 1 log(n) + log log(2) log 2 (n) and b n = n + 1 2n . ( 2
)
Although it is quite straight-forward to implement, the proposed procedure performs well in terms of classification compared to other estimation methods such as Zhang and Stephens (ZSE) [START_REF] Zhang | A new and efficient estimation method for the generalized pareto distribution[END_REF], Song and Song (SSE) [START_REF] Song | A quantile estimation for massive data with generalized pareto distribution[END_REF] or maximum likelihood (MLE). A comparative study is shown in Table 3. Our classification method manages to not favor any model in particular, as the probability of correct selection depends essentially on the distance between the scale parameter k and the closest transitional phase. This is well illustrated in Figure 2 where the probability of accurate selection appears to be nearly symmetric locally around 0 and 1. In comparison, classification by MLE has a clear tendency to overestimate the scale parameter while SSE tends to underestimate it. For a small sample size, these biases are particularly clear for k = 0.9 where model C is selected 88.8% of the time and for k = 0.1 where model A is selected 73.8% of the time by SSE. In fact, for the small sample case n = 15, our method is almost objectively best compared to the three other methods. For the examples considered, the situation n = 100, k ≈ 0 is the only scenario where our method is outperformed, ZSE being clearly best. This strongly suggests that the information on the tail of the distribution is not contained entirely in the mean/max statistic in this case. Nevertheless, the results of the method are overall quite satisfying given the small computational cost.

Model A Model B Model C k = -0.1 k = 0.1 k = 0.9 k = 1.1 A B C A B C A B C A B C ZSE n =

Applications

Fit a Generalized Pareto distribution

Due to the missclassification issue in Generalized Pareto models, inference methods tend to specialize in a specific region of the set of parameters, with none being uniformly best [START_REF] Del | Likelihood inference for generalized pareto distribution[END_REF]. To solve the problem, we propose a simple recipe that involves a two-step procedure. We assume that the data X = (X 1 , . . . , X n ) are independent and identically distributed from a GPD and that no prior information is available on the parameters.

Step 1. Compute the mean/max statistic τ n = X n /X (n) and compare it to the values a n and b n obtained in Eq. ( 2). Select:

• Model A if τ n < a n • Model B if a n < τ n ≤ b n • Model C if τ n ≥ b n
Step 2. Consider the estimation method depending on the selected model and sample size by this recommendation based on the analysis in [START_REF] Del | Likelihood inference for generalized pareto distribution[END_REF].

sample size Model A Model B Model C n ≤ 30 SSE ZSE MLE n > 30 ZSE ZSE MLE
In all situations, the method must be applied with the corresponding restriction on the set of parameters. Furthermore, when the sample is classified as Model B or C, a previous estimation of the endpoint can be recommended. Non-parametric approach shows satisfactory results, see [START_REF] Alves | A general estimator for the right endpoint with an application to supercentenarian women's records[END_REF].

This recipe is meant for data independently drawn from a Generalized Pareto distribution. For general problems related to tail detection and calibration, this assumes implicitly that the observations have been already extracted from the tail values of a larger dataset and suitably shifted so as to fit the GPD corresponding to the tail distribution. Nevertheless, the mean/max statistic can be used directly for tail detection. Suppose one is interested in the tail of the distribution of a sample Y 1 , ..., Y N , the idea is to calculate the mean/max statistic over the shifted tail observations

X i := Y (N -n+i) -Y (N -n) , i = 1, . . . , n
for various values of n. This way, one obtains a trajectory

τ n = X n X (n) = n i=1 (Y (N -n+i) -Y (N -n) ) n(Y (N ) -Y (N -n) ) , n = 1, 2, ...
that can be compared to the corresponding values of the bounds a n , b n . Ideally, only the values in the tail of the distribution are to be used for the classification so that n need not be too large. The optimal threshold for tail detection may appear more or less clearly (see Figure 4.1). variable with modified tail. On the left, the distribution is truncated at the value 1.5 (left), on the right, it is extended to an exponential tail starting from 1.5. The bounds an, bn from Eq. ( 2) are represented in straight lines along with their respective 95% confidence regions. The red-striped area corresponds to potential values of τn with more than 5% uncertainty to distinguish between uniform and exponential tails.

An overview of the mean/max trajectory can be helpful in order to select the correct model for the tail. In the examples displayed in Figure 4.1, the criterion is most conclusive around n = 20 even though a sample a this size may contain observations that do not obey the tail distribution. As a possible rule of thumb, the model for the tail can be decided based on the last value of τ n that leaves the red-striped area, thus aiming for a precision of over 95%.

Distribution of extreme seismic events

The Gutenberg-Richter (GR) law states that the distribution of seismic moment corresponds a to power-law distribution [START_REF] Sornette | Rank-ordering statistics of extreme events: Application to the distribution of large earthquakes[END_REF][START_REF] Corral | Scaling and universality in the dynamics of seismic occurrence and beyond[END_REF] with a density at a value M proportional to 1/M 1+β for a β evaluated to be approximately equal to 0.65. However, it is suspected that the collected data on seismic moments show a deviation of GR law for large values of M , entailing that the power-law model might have to be extended in order to add an exponential decay above a certain threshold [START_REF] Corral | Scaling and universality in the dynamics of seismic occurrence and beyond[END_REF]. In [START_REF] Yan Y Kagan | Seismic moment distribution revisited: I. statistical results[END_REF], Kagan enumerates the requirements that an extension of the GR law should fulfill. He also argues however that available seismic catalogs do not allow the reliable estimation of the threshold, except in the global case where the faster exponential decay may take place at the highest observed values of M , for which the available data are very poor [START_REF] Zöller | Convergence of the frequency-magnitude distribution of global earthquakes: Maybe in 200 years[END_REF].

The existence of a theoretical maximal magnitude earthquake entails that the true distribution of seismic moments must be right-truncated, thus displaying a uniform tail, see [START_REF] Del | Likelihood inference for generalized pareto distribution[END_REF]. Nevertheless, the question remains as to know if the largest seismic events on record are not best modeled by an exponential decay. In order to infer on the underlying physical model, we compute the max/mean plot of the seismic moment collected data from the centroid moment tensor (CMT) catalog [START_REF] Ekström | The global cmt project 2004-2010: centroid-moment tensors for 13,017 earthquakes[END_REF]. The analysis is restricted to shallow events (depth < 70 km) from 1976 to 2013, as recommended in [START_REF] Yan Y Kagan | Seismic moment distribution revisited: I. statistical results[END_REF]. The results are displayed in Figure 4 below. This analysis shows no reason to suspect a change in the tail of the distribution, and confirms that the tail must be classified as Model A with k < 0. Consequently, the existence of a cut-off close to the maximal observed seismic moments has to be discarded.

Detection of experimental limitations

We apply this methodology in the analysis of acoustic emission data in failure under compression experiments of a nanoporous quartz sample [1,17]. A cylindrical sample is placed between to plates that approach each other at a constant velocity. Compression is done with no lateral confinement and the experiment finishes when the sample has experienced a big failure and has literally disappeared. In this compressive process, two transductors are embedded in both plates in order to detect acoustic emission activity. These signals are pre-amplified in order to measure properly different magnitudes: amplitude, duration, energy, etc. The discrete measure in dB of the signal amplitude is defined as A = 20 log 10 (A v /1µV ), where A v is the highest voltage value of each acoustic emission signal and 1µV is a reference voltage. Remark that the exponential distribution for the tail, the largest values, is the most natural model from a physical point of view. Regarding this physical magnitude, one must take into account that the signal pre-amplification can lead to a saturation for the largest acoustic emission events. Therefore, this saturation must be understood as a fictitious cut-off which is inherent to the experimental set-up. Actually, this experimental limitations are present in all measuring devices since all of them have a certain measurement range. The mean/max statistic τ n enables to detect this experimental artifact in most samples, as we see Figure 5 where V.Navas and E.Vives (private communication) experiments are shown.

Compression experiments of Vycor

number of tail observations residual mean/max q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q D2H1 D1H005 D1H025 According to the experimental intuition, one expects that acoustic emission signals do not reach saturation if the contact area between the plates is small. The mean/max methodology provides a simple and automatic way to identify if the data extracted from any experimental measurement are saturated or not. Figure 5 shows the trajectories of the mean/max statistic for the 100 largest values obtained from three compression experiments of nanoporous quartz samples. The circle trajectories correspond to the greatest material, 2mm diameter. This mean/max plot shows saturated experiment, since for the 20 largest emission signals the exponential decay is rejected. In fact, the tail of the dataset is classified as a compact support, Model C. In this case, a descriptive analysis reveals a clear suspicion of saturation experiment, since the largest value is repeated too many times. However, for small materials, 1mm diameter, it is difficult to detect if the saturation occurs. The mean/max test (for the 20 largest) reveals that the triangle trajectories was affected by the saturation, as in the previous case, but the trajectory of squares is classified as Model B, therefore a non-saturation can not be rejected, since Model B includes the exponential distribution for the tail.

A Distribution of the mean/max statistic

We investigate the distribution of the mean/max statistic τ n when the sample is drawn from a GPD. The two transitional phases corresponding to uniform and exponential distributions are given a particular interest as they are directly linked to the level and power of the test discussed in Section 2.

A.1 The uniform case

The threshold c α corresponding to a test of significance level α ∈ (0, 1) in Theorem 2.1 follows from calculating the distribution of τ n under the null hypothesis of a uniformly distributed sample X 1 , . . . , X n . Noticing that

τ n = 1 n 1 + n-1 i=1 X (i) X (n) ,
it appears that n(τ n -1) is distributed as the sum of n-1 independent standard uniform random variables. This distribution is known as the Irwin-Hall distribution of parameter n -1, whose density is given by

I n-1 (t) = n 2 n-1 k=0 sgn(t -k) (-1) k (t -k) n-2 k!(n -k)! , t ∈ [0, n -1],
where sgn(y) is the sign of y, equal to 1 if y is positive, -1 if y is negative and 0 if y is zero. Geometrically, I n-1 (t) represents the volume of the intersection of the 1 -sphere of radius t in R n-1 with the hypercube [0, 1] n-1 . For sake of completeness, we derive the distribution of τ n under the null hypothesis in the next lemma. The density of a random variable Z will be denoted by f Z .

Lemma A.1. If X = (X 1 , . . . , X n ) is an iid sample from a uniform distribution on [0, θ], with θ > 0, then τ n has density

f τn (t) = nI n-1 (nt -1) , t ∈ 1 n , 1 .
In particular, E(τ n ) = med(τ n ) = n + 1 2n .

Proof. First remark that the distribution of τ n does not depend on θ so that we can assume that θ = 1 without loss of generality. Let W 1 = X (1) , W i = X (i) -X (i-1) for i = 2, . . . , n and W n+1 = 1 -X (n) . We use that W = (W 1 , . . . , W n+1 ) has a standard Dirichlet distribution Dir(n + 1), which means that

W d = Y 1 n+1 i=1 Y i , . . . , Y n+1 n+1 i=1 Y i
where Y 1 , . . . , Y n+1 is an iid sample with exponential distribution. Thus, the vector

W = X (1) X (n) , . . . , X (n-1) -X (n-2) X (n) d = Y 1 n i=1 Y i , . . . , Y n n i=1 Y i has Dir(n) distribution and Z i = X (i) /X (n) d = i j=1 Y j / n j=1 Y j , i =

A.2 The exponential case

Calculating the power of the test requires to know the distribution of τ n under H 1 . Here again, this distribution can be computed explicitly.

Lemma A.2. If X = (X 1 , . . . , X n ) is an iid sample from an exponential distribution of parameter λ > 0, then τ n has density

f τn (t) = n! n n-1 I n-1 (nt -1) t n , t ∈ 1 n , 1 .
Proof. The distribution of τ n does not depend on λ so that we can assume that λ = 1 without loss of generality. We know that the ordered sample has density on R n given by

(X (1) , . . . , X (n) ) ∼ n! exp - n i=1 x i 1{0 ≤ x 1 ≤ • • • ≤ x n }. Let Y i = X (1) /X (n) for i = 1, . . . , n -1. By the change of variable y i = x i /x n , we get (Y 1 , . . . , Y n-1 , X (n) ) ∼ n! x n-1 n exp -x n (1 + n-1 i=1 y i ) 1{0 ≤ y 1 ≤ • • • ≤ y n-1 ≤ 1, x n > 0}.
Integrating the density over x n , we find

Y = (Y 1 , . . . , Y n-1 ) ∼ n!(n -1)! (1 + n-1 i=1 y i ) n 1{0 ≤ y 1 ≤ • • • ≤ y n-1 ≤ 1}
To compute the density of S n := nτ n -1 =

A.3 Asymptotic distribution of τ n in the GPD model

Finally, we investigate the distribution of τ n when the observations are drawn from generalized Pareto distributions with shape parameter k. Because τ n is scale-free, the following results do not depend on the scale parameter σ of the GPD, which can be taken equal to one without loss of generality. The actual distribution of τ n being difficult to compute as a function of k, we only discuss the asymptotic distribution.

Theorem A.3. Let X = (X 1 , . . . , X n ) be an iid sample drawn from a Generalized Pareto distribution with scale parameter k. Then,

• If k > 0, τ n ∼ k k + 1 + N k √ n + o P 1 √ n ,
where N k has normal distribution N 0, k 2 (1 + k) 2 (1 + 2k) .

• If k = 0 (exponential case),

τ n ∼ 1 log(n) - G log 2 (n) + o P 1 log 2 (n)
,

where G has standard Gumbel distribution.

• If -1 < k < 0,

τ n ∼ W k n -k + o P (n k )
where W k has Weibull distribution with shape parameter -1/k and scale parameter -k/(k + 1).

Proof. We tackle each case separately. For k > 0, we have by the central limit theorem,

√ n X n - 1 1 + k d ----→ n→∞ N 0, 1 (1 + k) 2 (1 + 2k)
while X (n) converges a.s. towards 1/k, the upper bound of the support of the GPD. The result follows easily by Slutsky's lemma. In the exponential case k = 0, X n converges a.s. towards 1 while for t > -log(n), P X (n) -log(n) ≤ t = 1 -e -(log(n)+t) n = 1 -1

n 1+ t log(n) n ----→ n→∞ e -e -t .
Thus X (n) -log(n) converges towards a standard Gumbel distribution. The result follows in view of

1 X (n) = 1 log(n) -(X (n) -log(n)) = 1 log(n) + X (n) -log(n) log 2 (n) + o P 1 log 2 (n)
.

For the final case -1 < k < 0, we have by the strong law of large numbers

X n = 1 1 + k + o P (1)
For t > n k , write

P n k (1 -kX (n) ) ≤ t = P X (n) ≤ 1 k 1 - t n k = 1 - t 1/k n n ----→ n→∞ e -t 1/k .

Figure 2 :

 2 Figure 2: Monte-Carlo estimated probabilities of selecting the models A,B and C in function of the shape parameter k of the GPD, for sample sizes of n = 15 (left) and n = 50 (right). The straight line gives the probability of selecting the correct model.

Figure 3 :

 3 Figure 3: Examples of trajectories of the mean/max statistic τn from a sample of size 100 drawn from a standard Gaussian

Figure 4 :

 4 Figure 4: Trajectory of the mean/max statistic computed over the n largest (shifted) global seismic events, for n = 1 to n = 100.

Figure 5 :

 5 Figure 5: Trajectories of the mean/max statistic for the 100 largest values obtained from three compression experiments of nanoporous quartz samples.

1 ,

 1 . . . , n -1 has the same distribution as an ordered iid sample of n -1 uniform random variables on [0, 1]. We deduce that n-1 i=1 Z i = n(τ n -1) has Hirwin-Hall distribution with parameter n -1 and the result follows.

n- 1 i=1

 1 Y i , we now need to integrate the joint density over the level sets of the 1 -norm D s = {y ∈ R n-1 : n-1 i=1 y i = s}. We obtainf Sn (s) = Ds n!(ny i ) n 1{0 ≤ y 1 ≤ • • • ≤ y n-1 ≤ 1} dy 1 . . . dy n-1 = n! (1 + s) n Ds (n -1)! 1{0 ≤ y 1 ≤ • • • ≤ y n-1 ≤ 1} dy 1 . . . dy n-1 . Remark that (n -1)! 1{0 ≤ y 1 ≤ • • • ≤ y n-1 ≤ 1}is the density of an ordered sample (U (1) , . . . , U (n-1) ) of independent uniform random variables on [0, 1]. Thus,Ds (n -1)! 1{0 ≤ y 1 ≤ • • • ≤ y n-1 ≤ 1} dy 1 . . . dy n-1 = I n-1 (s)to which we deduce that S n = nτ n -1 has density n!I n-1 (s)/(1 + s) n over [0, n -1]. The result follows by a simple change of variable.

  the distribution of Y is constant over H 0 and over H 1 , Neyman-Pearson's Lemma tells us that the likelihood ratio test has maximal power among tests of significance level α. Both likelihood functions L 0 and L 1 of Y under H 0 and H 1 respectively can be calculated explicitly, yielding

Table 1

 1

	n	5	10	20	50	100
	A(n) 0.7178 0.8437 0.9514 0.9984	1
	c n	0.531	0.460	0.420	0.392 0.380

for different sample sizes n. It appears that the test quickly reaches a near perfect accuracy as the sample size increases. For a sample of size n = 5, the method selects the correct model more than 71% of the time, while a sample of size n = 20 achieves a precision above 95%.

Table 1 :

 1 Optimal

threshold cn and percentage of accuracy A(n) for sample sizes n = 5, 10, 20, 50, 100, computed by Monte-Carlo simulations with 5.10 4 replications.

Table 2 :

 2 Percentages

	15	65.8 34.1 0.1 44.1 55.5	0.4	0.9	76.5 22.6	0.4	57.7 41.9
	n = 100	84.2 15.8 0.0 19.9 80.1	0.0	0.0	90.7	9.3	0.0	31.0 69.0
	SSE
	n = 15	85.3 14.5 0.3 73.8 25.6	0.7	20.7 50.0 29.4 14.7 40.7 44.6
	n = 100	70.1 30.0 0.0 29.0 70.9	0.0	0.0	77.1 22.9	0.0	27.3 72.7
	MLE
	n = 15	41.4 50.2 8.4 18.8 64.6 16.6	0.1	11.2 88.8	0.0	4.1	95.9
	n = 100	75.0 25.0 0.0	9.7	90.3	0.0	0.0	49.8 50.2	0.0	4.4	95.6
	n = 15	59.7 39.9 0.4 38.6 60.3	1.1	0.2	59.6 40.2	0.0	40.7 59.3
	n = 100	74.4 25.6 0.0 20.4 79.6	0.0	0.0	79.8 20.2	0.0	22.3 77.7

of the model classifications obtained by Zhang and Stephens Estimation (ZSE), Song and Song Estimation (SSE) and Maximum Likelihood Estimation (MLE) compared to the proposed methodology (bottom lines). In bold are the percentages of correct classification.

Acknowledgement

We are grateful to V. Navas and E. Vives for their feedback. Research expenses were founded by projects MTM2012-31118 from Spanish MINECO, 2014SGR-1307 from AGAUR, and the Collaborative Mathematics Project from La Caixa Foundation.

We deduce that n k (1-kX (n) ) ∼ -kn k X (n) converges towards a Fréchet distribution with shape parameter -1/k as n → ∞. Thus, -n -k /kX (n) converges to the inverse of a Fréchet variable whose distribution is Weibull. The result follows by Slutsky's lemma.