Guillaume Sicard

Camille Salaün

Serena Ivaldi

Vincent Padois

Olivier Sigaud

Learning the velocity kinematics of ICUB for model-based control: XCSF versus LWPR

The model-based control of humanoid robots requires the availability of accurate mechanical models that can be hard to obtain in practice. One approach to this problem consists in calling upon machine learning methods. In this paper, using a standard control approach based on visual servoing, we compare the accuracy of two supervised learning methods, namely LWPR and XCSF, to extract the forward velocity kinematics of the upper body of the ICUB robot. Experiments are performed in simulation, using one arm and the head for reaching tasks. We show that both methods provide accurate models of the robot, with a slight advantage to XCSF over LWPR.

I. INTRODUCTION

Given the complexity of humanoid robots and of their associated missions, using machine learning techniques for the control of these platforms is becoming mandatory (see [START_REF] Sigaud | From motor learning to interaction learning in robots[END_REF] for an overview). Among different approaches, trying to learn models of the robot at the kinematic or the dynamic level has recently received a lot of interest (e.g. [START_REF] Vijayakumar | LWPR: A scalable method for incremental online learning in high dimensions[END_REF], [START_REF] Peters | Learning to control in operational space[END_REF], [START_REF] Nguyen-Tuong | Real-time local gaussian processes model learning[END_REF]). The idea consists in extracting these models with general function approximation methods based on input/output data that are recorded while using the robot.

In [START_REF] Salaün | Control of redundant robots using learned models: an operational space control approach[END_REF], a control framework combining the Resolved Motion Rate Control (RMRC) framework and learning the Forward Velocity Kinematics (FVK) of the system with LWPR [START_REF] Vijayakumar | LWPR: A scalable method for incremental online learning in high dimensions[END_REF] is presented. The specificity of this framework is that it can combine several tasks ranked by priority, either compatible or not. However, this study is limited in several respects. First, it is based on a simple planar arm model. Here, we evaluate the potential of this approach with the ICUB humanoid robot [START_REF] Metta | The iCub humanoid robot: an open platform for research in embodied cognition[END_REF]. In that respect, our robotic set-up and goal share a lot of similarities with the framework presented in [START_REF] Natale | Learning precise 3d reaching in a humanoid robot[END_REF] that uses JAMES, an ancestor of ICUB, and RFWR, an ancestor of LWPR. Second, in [START_REF] Salaün | Control of redundant robots using learned models: an operational space control approach[END_REF] an important motor babbling stage is required to initialize the learned FVK, which would be hardly feasible on a real robot. Here, the motor babbling stage is not necessary anymore. Finally, the focus of the paper is on an empirical comparison between LWPR and XCSF [START_REF] Butz | Computational Complexity of the XCS Classifier System[END_REF], another state-of-the-art function approximation algorithm.

In order to elaborate this comparison, we first study the performance of both learning techniques on 8 sequential reaching tasks with the simulator of the ICUB robot [START_REF] Tikhanoff | The icub humanoid robot simulator[END_REF], using a FVK based on the CAD parameters of the robot as baseline.

The paper is organized as follows. In Section II, we present the task, set-up and methods used to perform our comparison. In Section III, we describe the series of evaluations performed and present the results. We discuss these results in Section IV before concluding.

II. METHODS

In this section, we first describe the robotic set-up. Then, we give a formal presentation of the visual servoing and inverse velocity kinematics problem that we address. Finally, we briefly describe LWPR and XCSF, the learning algorithms used to extract the model out of experimental data.

A. Robotics set-up and objectives

Our goal is to compare the accuracy of two supervised learning algorithms, LWPR and XCSF, when approximating the FVK of ICUB, a 53 degrees of freedom and 104cm high humanoid robot. ICUB is actuated with DC brushless motors. It was designed in the context of the ROBOTCUB project 1and is present in several research laboratories in Europe. The ICUB software architecture is based on YARP [START_REF] Fitzpatrick | Towards long-lived robot genes[END_REF] which provides a powerful abstraction layer allowing, for example, to transparently test most of the code in simulation before actually executing it on the real robot without any additional implementation effort. This is possible because the ICUB rigid body dynamics simulator also uses YARP as its software interface.

In the following experiments, we use an arm and the head of ICUB in different contexts during a reaching task. The shoulder is a 3 degrees of freedom joint with some coupling between degrees of freedom. The decoupling is done with a low level controller. The elbow is a one degree of freedom joint. The neck is composed of two consecutive joints, namely the head yaw and pitch angles. It is actuated with two direct drive DC motors whereas for the arm, transmissions based on cables and pulleys are used.

In total, as illustrated in Fig. 1, we control 4 degrees of freedom for the arm and 2 degrees of freedom for the head 2 while the torso of the robot is attached to a fixed base. Fig. 1: Parametrization of the arm and head of ICUB. Head, shoulder and elbow controlled axes are represented as well as the head and torso frames.

To learn the FVK in simulation, we simulate the vision process as it would be implemented on the real robot. Furthermore, Gaussian noise is added to the control layer to emulate a physical setup. In order to provide a baseline for the comparison, an analytical expression of this model is obtained using the Kinematic Dynamics Library (KDL) from the Orocos Project [START_REF] Smits | Orocos: A software framework for complex sensordriven robot tasks[END_REF] using the parameters of the CAD model, which are available in the ICUB wiki 3 .

B. Visual servoing problem formulation

In terms of control, the goal of this paper is to derive, at the velocity kinematics level, a closed-loop controller for the end-effector of the considered kinematic chain (ICUB hand here) based on visual feedback from cameras of the robot. Our main focus being not vision, we define the end effector of the robot as a green ball that the robot holds with a stick (cf. Fig. 1). We use the Cartesian position of the center of the ball expressed in the head frame as features ξ = (ξ x , ξ y , ξ z )

T .

This position is reconstructed based on a reprojection software module that uses a ball tracker algorithm based on the OpenCV software [START_REF] Bradski | The OpenCV Library[END_REF] as well as on the camera calibration matrices.

The visual servoing problem consists in ensuring the reaching or the tracking of the desired values of the features ξ = ξ x , ξ y , ξ z T . In that case, both the target and the end effector need to be visible at any time. An exponential decrease of the error is achieved by a proportional controller ξ = λ (ξ -ξ) where λ is a positive definite matrix. In order to design a velocity-level controller for our system, we thus need to relate ξ to the articular velocity of the robot.

3 http://eris.liralab.it/wiki/ICubForwardKinematics

C. Inverse Velocity kinematics

From the structure of the robot, we know that ξ is a function of both the head configuration q h = q hp , q hy T and the arm configuration q a = q sp , q sr , q sy , q e T , i.e. we have a surjective non-linear function f ξ , namely the forward kinematics, that relates the configuration of the robot q = q h T , q a T T

to features ξ = f ξ (q).

The considered kinematic chain being redundant with respect to the end-effector task, there is an infinite number of possible inverses for f ξ . However, there is no simple method to span the set of possible solutions at the geometric level and this type of mapping is often described at the velocity level. A a matter of fact, we can derive this expression and obtain

ξ = ∂f ξ (q) ∂q h ∂f ξ (q) ∂q a J (q) qh qa q . ( 1
)
The terms

∂f ξ (q)
∂q h and

∂f ξ (q)
∂q a define two Jacobian matrices: J ξ,h (q), which relates ξ to the head velocity, and J ξ,a (q), which relates ξ to the arm velocity. Thus (1) can be rewritten

ξ = J ξ,h (q) qh + J ξ,a (q) qa . (2)
In order to keep the target ξ visible by the cameras of the robot, we choose to place it along the line of sight and reach it with the end effector. For that purpose we need to control the head and the arm independently from one another. Since we only control the head pitch rate qhp and yaw rate qhy , we can implement the strategy of [START_REF] Natale | Learning precise 3d reaching in a humanoid robot[END_REF] to control the head in order to place the target along the line of sight

q hp = K hp ξ y + K D hp ξ y q hy = K hy ξ x + K D hy ξ x , (3)
which leads asymptotically to ξ x → 0 and ξ y → 0.

Defining q h = q hp , q hy T and ξ a = ξ -J ξ,h (q) q h , the arm FVK (2) can be written

ξ a = J ξ,a (q) q a , (4)
where q a is the velocity to apply to the arm in order to exponentially decrease the feature errors given a head velocity q h used to place the target along the line of sight.

In the non singular case, this mapping is redundant [START_REF] Salaün | Control of redundant robots using learned models: an operational space control approach[END_REF] and there is an infinite number of solutions q a satisfying (4). Among them, the one provided by the Moore-Penrose pseudo-inverse J ξ,a (q) + provides minimum norm solutions (see [START_REF] Doty | A theory of generalized inverses applied to Robotics[END_REF]). To avoid critical effects due to singular configurations, the Damped Least Square Pseudo-inverse (DLS-PINV) of J ξ,a (q) is preferred. It can be computed based on the SVD of J ξ,a (q) but an intuitive expression of the DLS-PINV is given by J ξ,a (q) + = J ξ,a (q) T J ξ,a (q) J ξ,a (q

) T + k 2 I -1 (5)
where I is the identity matrix. The parameter k 2 can be determined adaptively in different ways and we have chosen to use the method proposed in [START_REF] Chiaverini | Achieving userdefined accuracy with damped least-squares inverse kinematics[END_REF],

k 2 = 0 σ min > σ [1 -( σmin σ ) 2 ]k 2 max σ min < σ (6)
where σ min is the smallest singular value of the Jacobian matrix, k max is a regularization term that we set to 0.316 × (180/π) and σ is a threshold which determines if the plant is in singularity or not, that we set to 7.10 -4 . The obtained solution is not exact but is robust to singular configurations while ensuring minimal errors in any other case. Redundancy can also be to partially achieve a secondary task while not disturbing the main one (in our case, reaching). This decoupling requires the computation of P Jξ,a (q), the projector onto the null space of J ξ,a (q). It can be done using the SVD of J ξ,a and the column vectors v i of the resulting matrix V . As a consequence, the general form of the solution to (4) can be written

q a = J ξ,a (q) + ξ a + P Jξ,a (q) q0 , (7)
where q0 = -k null ∂H ∂q is chosen accordingly to [START_REF] Marchand | Using the task function approach to avoid robot joint limits and kinematic singularities in visual servoing[END_REF] in order to maintain the configuration far from joint limits (k null > 0).

The use of the inverse velocity kinematics scheme described by [START_REF] Natale | Learning precise 3d reaching in a humanoid robot[END_REF] requires to have access to J (q). The next subsections give a brief overview of the learning methods used in this work to incrementally approximate J (q) based on joint positions and velocities as well as visual features measurements.

D. LWPR in a nutshell

The Locally Weighted Projection Regression (LWPR) algorithm [START_REF]Locally Weighted Projection Regression[END_REF] is a recursive function approximator, which provides accurate approximation in very large spaces at low computational cost. It uses a sum of linear models weighted by normalized Gaussians. These Gaussians, also called receptive fields, define a zone of influence for each corresponding linear model. The receptive fields and the corresponding linear models are both updated incrementally to match the training data. LWPR reduces the input dimensionality using the Partial Least Squares (PLS) algorithm [START_REF] Wold | Soft modelling by latent variables: the non-linear iterative partial least squares (NIPALS) approach[END_REF], [START_REF] Elden | Partial least-squares vs. Lanczos bidiagonalization-I: analysis of a projection method for multiple regression[END_REF], [START_REF] Tenenhaus | La régression PLS: théorie et pratique[END_REF]. The global algorithm provides as output the weighted sum of all outputs of each receptive field

ŷ (x) = K k=1 w k ŷk (x) K k=1 w k (8)
where K is the number of receptive fields. We refer the reader to [START_REF] Schaal | Scalable techniques from nonparametric statistics for real time robot learning[END_REF] or [START_REF] Vijayakumar | Locally weighted projection regression: An o (n) algorithm for incremental real time learning in high dimensional space[END_REF] for a presentation of the incremental version of the algorithm.

E. XCSF in another nutshell

XCSF [START_REF] Wilson | Classifiers that Approximate Functions[END_REF] is another function approximator that shares some similarities with LWPR but comes from Learning Classifier Systems (LCSs) [START_REF] Sigaud | Learning classifier systems: A survey[END_REF]. As any LCS, XCSF manages a population of rules, called classifiers. These classifiers contain a condition part and a prediction part. In XCSF, the condition part defines the region of validity of a local model whereas the prediction part contains the local model itself.

A classifier defines a domain φ i (z) and uses a corresponding linear model β i to predict a local output vector y i relative to an input vector x i . The linear model is updated using the Recursive Least Squares (RLS) algorithm, the incremental version of the Least Squares method. The classifiers in XCSF form a population P that clusters the condition space into a set of overlapping prediction models. XCSF uses only a subset of the classifiers to generate an approximation. Indeed, at each learning iteration, XCSF generates a match set M that contains all classifiers in the population P whose condition part Z matches the input data z i.e., for which φ i (z) is above a threshold φ 0 4 .

In XCSF, the output ŷ is given for a (x, z) pair as the sum of the linear models of each matching classifier i weighted by its fitness

F i ŷ (x, z) = n M k=1 F k (z) ŷk (x) n M k=1 F k (z) (9)
where n M is the number of classifiers in the match set M . In all other respects, the mechanisms that drive the evolution of the population of classifiers are directly inherited from XCS and are described in [START_REF] Butz | Context-dependent predictions and cognitive arm control with XCSF[END_REF].

An important process in the context of this study is compaction. At the end of a learning process, the final population is composed of highly overlapping classifiers. To reduce the size of the population, XCSF uses a Closest Classifier Matching (CCM) rule to have a fixed size match set M [START_REF] Butz | Context-dependent predictions and cognitive arm control with XCSF[END_REF].

III. EXPERIMENTS AND RESULTS

We use LWPR or XCSF to learn J (q), the overall FVK expressed in the head frame of a ball handed by the robot with a stick. More specifically, the retained control scheme requires the access to J ξ,h (q) and J ξ,a (q) which are blocks of J (q). As a matter of fact, we learn the FVK under a matrix form. Learning the FVK directly with LWPR would not provide such a matrix. Thus, instead, when using LWPR we learn f ξ (q) and use the first order derivative of the learned function (in our case, J (q)). Using XCSF, the Jacobian matrix is learned directly providing q as condition parameter and q, ξ as prediction parameter. The capacity to separate those spaces leads also to the possibility of predicting the learned model with only one condition parameter. In Fig. 4, one can see the trajectory, in 2D, performed by the end effector during the first, 17th and 34th experiments, with the KDL model, LWPR, and XCSF, respectively. One can see that the trajectory does not change with the KDL model, apart from some minor variability due to control and vision noise. Moreover, the trajectory is initially worse with LWPR and XCSF, it improves faster with LWPR than with XCSF, but finally the XCSF model gets even more accurate after condensation. Those results are confirmed in Fig. 5, where one can see the evolution of the time necessary to perform a complete asterisk motion. Actually, the 34th asterisk is performed in 77 seconds with KDL, 71 with LWPR and 56 with XCSF. With LWPR, the learning output is ξ whereas XCSF uses q as condition parameter and q, ξ as prediction parameter. Fig. 4: Simulation results. First line: trajectory of the end effector using KDL. Second line: using LWPR. Third line: using XCSF. The columns represent the evolution between the first, the 17th and the 34th asterisk experiment respectively.

IV. DISCUSSION

First, independently from the model used for control, the reaching results and associated Cartesian trajectories are far from perfect. This is partly due to the fact that at the joint torque level of ICUB, control is performed in a decentralized manner: the whole-body dynamics is not accounted for and torque at joint i is computed using a PID like control structure only relying on the error in velocity at joint i.

Moreover, it is often objected to the use of machine learning asterisk duration [s] asterisk index Fig. 5: Evolution of the time necessary to perform an entire asterisk task with the three used algorithms: KDL in red, LWPR in green and XCSF in blue. methods for FVK identification that it is easy to get this model analytically from the CAD model, which is itself generally very accurate. Our experiments show that the controllers using the FVK learned from both LWPR and XCSF outperform the one relying on the KDL model. This surprising result is explained by three factors. First, the control loop relies on visual information that is slightly wrong mainly because of poor calibration of the cameras of the robot. Indeed, in the case of ICUB, the origin of the visual frame of reference is not very precisely specified. Second, a detailed analysis of the KDL model revealed that it is slightly erroneous due to the evolution of the prototype over the years 5 . Third, in the case of real robot experiments, the position of the stick in the hand of the robot may not be precisely specified and may change from an experiment to another, which may induce errors in the estimation of the ball position relatively to the hand. The ball can be considered as a generic tool which is usually difficult to model accurately in the FVK of the robot. Learning algorithms can compensate for to all those uncertainties. Finally, XCSF slightly outperforms LWPR on the tasks tested in this paper. Furthermore, as Fig. 2 shows, the final model is about ten times smaller with XCSF than with LWPR (about 100 classifiers against about 1350 receptive fields). LWPR is one of the most used machine learning methods. It has been evaluated on higher dimensional mechanical systems than XCSF, but always learning the model only along specific trajectories [START_REF] Sigaud | On-line regression algorithms for learning mechanical models of robots: a survey[END_REF], whereas here we learn over the whole joint space of the studied system. Recently, [START_REF] Stalph | A Comparative Study: Function Approximation with LWPR and XCSF[END_REF] concluded that XCSF was outperforming LWPR on simple function approximation problems. Our work shows that this is also true for learning the FVK of the upper body of ICUB under realistic simulation conditions, which is a more significant result. Nevertheless, a gap exists between the simulated and the real robot and optimal parameters obtained in simulation may not be optimal for real world experiments, thus we must now perform evaluations on the real robot.

V. CONCLUSION AND PERSPECTIVES

In this paper we have evaluated the applicability of a control framework based on visual servoing and model learning to the problem of controlling the ICUB humanoid robot. More precisely, we have compared the performance of LWPR and XCSF, showing that both were able to learn a model that is appropriate for control, with a slightly better performance and a much smaller model for XCSF.

Preliminary trials have shown that the approach can be transferred to the real robot, despite a greater variability due to additional sources of noise. In our immediate research agenda, we have to quantitatively evaluate the results on the real robot. Then we want to show that our approach to visual servoing with learning a separate model for the head and for the arm can be exploited to deal with the case where the robot is seeing the target, but not its arm. On a longer term, we will make profit of the access to the dynamics resulting from the availability the IDYN library [START_REF] Ivaldi | Computing robot internal/external wrenches by means of inertial, tactile and f/t sensors: theory and implementation on the icub[END_REF] to try to learn the dynamics model of the robot in the context of interaction with objects.

Fig. 3 :

 3 Fig. 3: On-line learning and closed-loop control. Each box is a YARP module

Fig. 2 :

 2 Fig.2: Evolution of the normalized mean square error (left blue circle) and number of receptive fields (right red plain line) during 40000 learning steps with LWPR (a) and XCSF (b). With LWPR, the learning output is ξ whereas XCSF uses q as condition parameter and q, ξ as prediction parameter.

  All the authors are with: Université Pierre et Marie Curie, Institut des Systèmes Intelligents et de Robotique -CNRS UMR 7222, Pyramide Tour 55 -Boite Courrier 173, 4 Place Jussieu, 75252 Paris CEDEX 5, France.

Contact: firstname.lastname@isir.upmc.fr This work is supported by the French ANR program (ANR 2010 BLAN 0216 01), more at http://macsi.isir.upmc.fr

http://www.robotcub.org

see http://eris.liralab.it/wiki/ICub joints for more information about iCub joints

This threshold is named θm in[START_REF] Butz | Context-dependent predictions and cognitive arm control with XCSF[END_REF]

The model has been updated shortly after the realization of the research described here.

A. Tuning the learning algorithms LWPR with XCSF come with a large set of parameters and the performance is highly dependent on the effort put in tuning these parameters. To circumvent this difficulty, our methodology consists in the systematic exploration of these parameters with an equal amount of time attributed to tuning both algorithms. All parameters are tested on random joint positions and velocities as input. The output is computed using an existing KDL model as a baseline which provides a good estimation of the behavior of the system while speeding-up the exploration process. We consider having a sufficiently accurate prediction with 3.10 The parameters used for subsequent experiments are in bold face. They are chosen so that the NMSE converges to the lowest asymptotic value while decreasing fast in the first steps. The input dimension is 6 and the output dimension is 3. The norm_in parameter is set to 180. The corresponding performance evolution is shown in Fig. 2(a).

As XCSF may use different types of prediction space (quadratic or linear) and condition spaces (ellipsoid, rectangle, sphere, rotating or not), we choose linear predictions and rotating ellipsoid conditions in order to easily compare XCSF with LWPR. The condition and prediction input dimensions are set to 6. The output dimension is 3. The maximum learning iteration parameter is set to 150000. Since the error is approximately the same for startCompaction = 0.4 and 0.6, we finally choose startCompaction = 0.5. The corresponding performance evolution is shown in Fig. 2(b).

B. Performing the asterisk experiment with the iCub simulator

In this experiment, we compare the results obtained with LWPR and XCSF in a reaching task. Results obtained with KDL are provided as a baseline. Our approach is tested on the 'star-like' asterisk task (see [START_REF] Nakanishi | Operational space control: A theoretical and empirical comparison[END_REF]) which consists in a go and return from a center point to eight target points that are visited sequentially in a clockwise manner. Unlike [START_REF] Nakanishi | Operational space control: A theoretical and empirical comparison[END_REF] where the whole trajectory is specified, only the target points are given. The center of the asterisk is placed in ξ * = (ξ * x = 0.13m, ξ * y = 0.3m, ξ * z = 0.24m) T relatively to the torso frame (cf. Fig. 1). The radius of the asterisk is 9cm and a target is considered to be reached under a threshold of 1cm.

The closed loop control and learning implementation are described in Fig. 3. The learning algorithm module provides J (q) after a learning period. The visualization module is a display using a virtual camera. The reprojection module computes the 3D position of the end effector with respect to the head. The control loop module computes the joint velocities needed to track the target with the head and the arm independently, as described in Section II-C.


Mariem Slimani
email: slimani.mariem@telecom-paristech.fr

Fernando Silveira
email: silveira@fing.edu.uy

Philippe Matherat
email: philippe.matherat@telecom-paristech.fr

Variability-Speed-Consumption Trade-off in Near Threshold Operation

Keywords: Sub-threshold logic, Near-threshold operation, Variability, Modeling

Sub-threshold operation is an efficient solution for ultra low power applications. However, it is very sensitive to process variability which can impact the robustness and effective performance of the circuit. On the other hand this sensitivity decreases as we move towards near-threshold operation. In this paper, the impact of variability on sub-threshold and near-threshold circuit performance is investigated through analytical modeling and circuit simulation in a 65 nm industrial low power CMOS process. We show that variability moves the effective minimum energy point towards the near threshold region. Thus, we demonstrate that when variability is taken into account, a complete model that includes the near threshold (moderate inversion) region is necessary in order to correctly model circuit performance around the minimum energy point. Finally, we present the resulting speed-consumption trade-off in a variability-aware analysis of sub-threshold and near-threshold operation.

Introduction

Over the last decade, sub-threshold logic has been used as an ideal option to achieve Ultra Low Energy consumption for applications with low demand in speed requirements. Here, sub-threshold term refers to the weak inversion (WI) region where the minimum energy can be achieved using a supply voltage V DD well below the threshold voltage.

As the interest for ultra low energy has increased, research related to subthreshold logic has attained considerable importance. Modeling and characterization of sub-threshold operation for standard CMOS cell designs have been investigated for energy and performance analysis [START_REF] Vittoz | Weak inversion for ultimate low-power logic: in Low-Power CMOS Circuits[END_REF] [START_REF] Wang | A 180-mV subthreshold FFT processor using a minimum energy design methodology: Solid-State Circuits[END_REF]. However, several works have shown that variability in subthreshold logic, especially intrinsic one due to Random Dopant Fluctuations (RDF), is a critical limit to achieve robust ultra low energy devices as it cannot be compensated through adaptive body biasing (ABB) [START_REF] Hanson | Energy optimality and variability in subthreshold design: in Low Power Electronics and Design[END_REF][4] [START_REF] Bol | Interests and limitations of technology scaling for subthreshold logic: Very Large Scale Integration (VLSI) Systems[END_REF]. As currents in weak inversion region exponentially depends on threshold voltage (V t ), random V t variations significantly affect the on and off currents and gate delays and can affect the output swings and result in functional failures of some gates. Moreover, the minimum energy point can be strongly affected by variability. As reported in [START_REF] Bol | Nanometer MOSFET effects on the minimum-energy point of 45nm subthreshold logic[END_REF], device variability leads to 20% of minimum energy increase.

Models considering variability have been developed in previous works specifically for use in subthreshold, considering, therefore, just the weak inversion region. In this paper, we will show through 1K-point Monte Carlo Spice simulations in a Low Power (LP) technology from an industrial foundry that variability moves the effective minimum energy point towards the near-threshold region (Moderate Inversion). Thus, restrictive weak inversion models can no longer model circuit performance around the minimum energy point. Therefore, we apply a complete and compact transistor model valid from weak to strong inversion in order to model the delay and energy performance over the weak and moderate inversion regions. This will allow us to correctly model the circuit performance in a variability aware analysis.

This paper is organized as follows. Section 2 presents the concept of subthreshold operation and investigates the impact of process variations on minimum energy point. In section 3, we propose a complete model that includes the near threshold region. We show through Monte Carlo Spice simulations that the new model is necessary to correctly model the circuit performance in a variability aware analysis.

2

Sub-threshold circuit design

In this section, the concept of sub-threshold operation is briefly explained. Variability impact on sub-threshold circuit performance through Monte Carlo Spice simulations of an inverter chain implemented in a Low Power technology from an industrial foundry is then presented.

Sub-threshold operation

Sub-threshold operation consist in reducing the power supply voltage V DD below the threshold voltage V T in order to achieve minimum energy consumption. The concept is simple: as P dyn is proportional to the square of V DD , a small reduction in supply voltage causes quadratic decrease in dynamic power consumption at the cost of a significant increase in delay. This results in an increase in leakage energy that leads to a minimum energy point achieved at an optimum supply voltage. In [START_REF] Calhoun | Characterizing and modeling minimum energy operation for subthreshold circuits: in Low Power Elec-tronics and Design[END_REF], authors provide an analytical solution for the optimum V DD that minimize the energy for a given operating frequency. However, lowering the power supply voltage will expose the circuit to the effect of process variations which can impact the functionality of the circuit and result in a functional failure of the gate. Thus, a process corners analysis which determines the minimum operation voltage and the minimum transistor sizing which ensures a good functionality of the gate is the first step to do, as described in [START_REF] Wang | A 180-mV subthreshold FFT processor using a minimum energy design methodology: Solid-State Circuits[END_REF]. We have applied this analysis to an inverter from a 65nm Low Power industrial library. Figure 1 shows the limits for the ratio of the pMOS width to nMOS width, which must be between the minimum set by the Fast 0-Slow 1 (FS) corner and the maximum set by the Slow 0-Fast 1 (SF) corner in order to assure an output swing of 10%-90% of the supply voltage. We observe that this restriction determines a minimum possible operating voltage, which for this technology occurs at 143mV by sizing the PMOS width to be 1.83 the NMOS one. We see that the inverter of the standard cell library is guaranteed to operate down to 160mV. Thus, if the minimum energy point is achieved by a supply voltage between 143mV and 160mV, a modified logic cell library should be created.

Variability effect on subthreshold circuit

Process variations is a critical limit of sub-threshold circuits. Intrinsic variations due to Random doping fluctuation (RDF) is considered to be the dominant source of variability in sub-threshold circuits [START_REF] Bol | Interests and limitations of technology scaling for subthreshold logic: Very Large Scale Integration (VLSI) Systems[END_REF].

To show the impact of variability on sub-threshold circuit performance, Monte Carlo Spice simulations with 1000 points have been performed for different values of V DD . The considered benchmark circuit is a chain of inverters where the first inverters are not considered in order to correctly calculate the static current, increased due to the degeneration of stable states as mentioned in [START_REF] Vittoz | Weak inversion for ultimate low-power logic: in Low-Power CMOS Circuits[END_REF]. Since delay variability depends on the logic depth of the circuit, an appropriate choice of the logic depth is essential. For our case study, we have choosen a logic depth of 23 like in [START_REF] Bol | Pushing ultra-low-power digital circuits into the nanometer era[END_REF] where the circuit is a standard 8-bit ripple-carry-array (RCA) multiplier.

Figure 2 shows the evolution of total energy consumption without and with variability effect, considering typical and 3σ worst case delay, respectively. We observe that variability leads to a considerable increase in leakage energy. This results in the increase of the minimum energy point by 50%, located now in the moderate inversion region.

Variability aware circuit model

In [START_REF] Vittoz | Weak inversion for ultimate low-power logic: in Low-Power CMOS Circuits[END_REF], [START_REF] Bol | Interests and limitations of technology scaling for subthreshold logic: Very Large Scale Integration (VLSI) Systems[END_REF] and [START_REF] Bol | Nanometer MOSFET effects on the minimum-energy point of 45nm subthreshold logic[END_REF], models in WI region are applied to describe sub-threshold operation. However, as we have seen in section 2, variability considerably affects the minimum energy point which moves in the moderate inversion region. In [START_REF] Vittoz | Weak inversion for ultimate low-power logic: in Low-Power CMOS Circuits[END_REF] the impact of operation in moderate inversion applying an all region transistor model is also considered, but the impact of variability is not analyzed in this case. [START_REF] Markovic | Ultralow-power design in nearthreshold region[END_REF] is another prior work where authors suggest to work in the near threshold voltage region in order to recover some of delay performance at the expense of a little energy increase. An energy-delay modeling framework that extends over all inversion regions is developed in this paper. But variability is still not considered in these models.

In this section, we present a variability aware model that extends over the weak and moderate inversion region. This model is based on the EKV model expressions [START_REF] Vittoz | Weak inversion for ultimate low-power logic: in Low-Power CMOS Circuits[END_REF]. The main contribution here is that we consider V T and β variations in our analysis. We show through Spectre and Matlab simulations that the WI model in no longer sufficient to model the performance of a system exposed to process variations.

Current and delay model under variability analysis

In weak and moderate inversion region, the drain current can be expressed as [START_REF] Vittoz | Weak inversion for ultimate low-power logic: in Low-Power CMOS Circuits[END_REF]:

I DS =I S (ln 2 [1 + exp V GS -V T 2nU T ]- ln 2 [1 + exp V GS -V T 2nU T exp -V DS 2U T ]) (1)
Where n is the subthreshold slope factor, V GS and V DS are respectively the gate to source and drain to source voltages and V T is the threshold voltage. I S is the specific current given by

I S = 2nµC ox U 2 T W/L = 2nβU 2 T ( 2
)
Where µ is the mobility, C ox is the oxide capacitance, U T is the thermal voltage and W/L denotes the channel width-length ratio of the transistor. Equation 1 tends to the classical exponential WI model when V GS -V T is negative.

The expression of delay can be derived from the current model as follows:

T d = C L V DD I on (3)
Where C L is the load capacitance, V DD is the supply voltage and I on is the saturated on-current.

Leakage current is also determined from the I DS expression when V GS = 0. The model of Equation 1 is a long channel model that does not include effects such as mobility reduction, velocity saturation and drain induced barrier lowering (DIBL). The former two are mainly of impact in the strong inversion region and that is why, for simplicity sake, were not considered in this work. The last one (DIBL), has some impact on the performance in the WI and MI regions, as analyzed in [START_REF] Bol | Interests and limitations of technology scaling for subthreshold logic: Very Large Scale Integration (VLSI) Systems[END_REF]. In our case we considered the effect of DIBL when extracting the parameters for the model of Equation 1 by considering the drain current data of the 65nm Standard Threshold Voltage, Low Power (SVTLP) NMOS transistor with drain voltage equal to V G . To extract model parameters, we have applied the method based on the gm/ID curve described in [START_REF] Jespers | The Gm/ID Methodology, a Sizing Tool for Low-voltage Analog CMOS Circuits: The Semi-empirical and Compact Model Approaches[END_REF]. The following values were obtained:

-n = 1.22; -V T = 0.38(V); -β = 4.83e -4 (A/V 2 ).
Figure 3 shows the I D versus V GS for NMOS transistor determined with the weak inversion model, the complete model and spice simulations. As expected, the weak inversion model is not sufficient to model the current in near-threshold region. We observe that the complete model is not so accurate in strong inversion region due to the lack of modeling of mobility reduction and velocity saturation [START_REF] Jespers | The Gm/ID Methodology, a Sizing Tool for Low-voltage Analog CMOS Circuits: The Semi-empirical and Compact Model Approaches[END_REF]. For variability analysis, we will consider just random V T and β variations, modeled as a normal distributions (µ VT , σ VT , µ β , σ β ), determined through Monte Carlo simulations or through analytical expressions given in [START_REF] Bol | Interests and limitations of technology scaling for subthreshold logic: Very Large Scale Integration (VLSI) Systems[END_REF]. Current model considering process variations is derived from Equation 1 by replacing β and V T by values that follow the normal distribution described in the previous paragraph.

The evolution of typical and 3σ Worst Case(WC) delay is presented in Figure 4 (left). We observe that the complete model with and without variability consideration follow perfectly Monte-carlo Spice simulations whereas the delay of the WI model deviates from 0.3V of V dd .

Figure 4 (right) plots delay variability versus the supply voltage. The simulated variability, obtained through Monte-Carlo simulations, shows that variability decreases as the supply voltage increases. We remark that delay variability, obtained with the WI model, remains constant at different supply voltage, while, the one obtained with the complete model presents the same shape as spectre simulations and has even close values. We conclude that the WI model is a restrictive model that can not be used to model process variations and that the model developed, considering V T and β variations, is a good model for variability analysis.

The delay model presented in Equation 3 is valid for a simple gate. For a circuit with a logic depth L D , the delay will be T circuit = L D .T d . If T d is a normal distribution defined by (µ delay , σ delay ), T circuit will be a normal distribution too defined by (L D .µ delay , (L D ).σ delay ). Thus, the delay variability of the circuit can be obtained as follow:

var circuit-delay = (1/ (L D )).var gate-delay (4)
Table 1 lists delay variability for different Logic depth at V DD = 0.2V . As expected, the delay variability of a circuit decreases as its logic depth increases, and the decrease follows perfectly 1/ (L D ) law.

Energy model under variability analysis

The total energy consumed by the circuit is the sum of the dynamic energy E dyn , required to charge and discharge parasitic capacitances during logic transitions, and static energy E stat due to leakage currents I leak . This can be summerized in the following expression :

E tot = αC L V 2 DD + V DD I leak T circuit ( 5
)
Where α is the switching activity of the circuit and C L is the load capacitance.

Here, we consider Justintime operation [START_REF] Seok | Analysis and optimization of sleep modes in subthreshold circuit design[END_REF], where the circuit works in its maximum frequency, i.e., the period is set to be the critical path delay of the circuit.

To consider variability in energy analysis, 3σ worst-case delay and mean I leak are considered in the static energy calculation as follows

E stat = V DD µ I leak µ delay (1 + 3 * σ delay µ delay ) (6)
Figure 5 shows the consumed energy under process variations consideration. We observe that the total energy consumed is slightly different from that determined by the models. The error of the energy on the minimum point is of 5% and 6%, when obtained by the complete and the WI model, respectively. Not what we expect, the error of the WI model is comparable to that of the complete model as shown in Figure 5 (right). This can be explained by the inverse tendency of variability and delay determined by the WI model. On the one hand, the variability of the WI model is constant whatever the value of the supply voltage. It is therefore overestimated in the moderate and strong inversion regions. On the other hand, the delay obtained by the WI model is underestimated with respect to the one obtained by Spice simulations as observed in Figure 4. As the energy contains the product of variability and the delay, there is a compensation that let the WI model remain a good model of energy consumption even under variability analysis.

Conclusion

In this paper, modeling under variability analysis is investigated. We show that the Weak Inversion model, normally used to describe sub-threshold circuit, is a restrictive model that can no longer model circuit delay around the minimum energy point of circuit exposed to process variations.

We develop a complete model that extends over the weak and moderate inversion regions. Through Monte Carlo Spice simulations of a chain of inverters in a Low Power technology, we show that the new model is necessary to correctly model the variability of the circuit.

Nevertheless, instead of what may be expected, the Weak Inversion model remains a good model of energy consumption even under variability analysis. This is due to the compensation of delay and variability errors resulting when the WI model is applied.

Fig. 1 .

 1 Fig. 1. Worst case corners analysis of the inverter for different supply voltage.

Fig. 2 .

 2 Fig. 2. Leakage and total energy evolution with and without variability.

Fig. 3 .

 3 Fig. 3. ID versus VGS curves for NMOS transistor.

Fig. 4 .

 4 Fig. 4. Evolution of typical and WC delay (left), and normalized delay variability (right) for different VDD.

Fig. 5 .

 5 Fig. 5. Consumed energy under process variations (left), and % of complete and WI model error compared to Spice simulations (right).

Table 1 .

 1 Delay variability for different logic depth at V dd =0.2 V.

	LD	Delay variability (σ delay /µ delay )
		Spice simulation Analytical model
	1	0.86	0.99
	7	0.28	0.36
	15	0.25	0.26
	23	0.205	0.209

Acknowledgments.

The authors would like to thank STIC AmSud Program for the financial support given to this work through the NanoRadio project.


Clément Dombry
email: clement.dombry@math.univ-poitiers.fr

Frédéric Eyi-Minko
email: frederic.eyi.minko@math.univ-poitiers.fr

Extremes of independent stochastic processes: a point process approach

Keywords: extreme value theory, partial maxima process, superextremal process, functional regular variations, weak convergence. AMS Subject classification. Primary: 60G70 Secondary: 60F17

For each n ≥ 1, let {X in , i 1} be independent copies of a nonnegative continuous stochastic process X n = (X n (t)) t∈T indexed by a compact metric space T . We are interested in the process of partial maxima

where the brackets [ • ] denote the integer part. Under a regular variation condition on the sequence of processes X n , we prove that the partial maxima process Mn weakly converges to a superextremal process M as n → ∞. We use a point process approach based on the convergence of empirical measures. Properties of the limit process are investigated: we characterize its finite-dimensional distributions, prove that it satisfies an homogeneous Markov property, and show in some cases that it is max-stable and self-similar. Convergence of further order statistics is also considered. We illustrate our results on the class of log-normal processes in connection with some recent results on the extremes of Gaussian processes established by Kabluchko.

Introduction

A classical problem in extreme value theory (EVT) is to determine the asymptotic behavior of the maximum of independent and identically distributed (i.i.d.) random variables (Z i ) i≥1 . What are the assumptions that ensure the weak convergence of the rescaled maximum max

1≤i≤n Z i -b n a n , a n > 0, b n ∈ R,
and what are the possible limit distributions ? These questions were at the basis of the development of EVT and found its answer in the Theorem by Fisher and Tippett [START_REF] Fisher | Limiting forms of the frequency of the largest or smallest member of a sample[END_REF] characterizing all the max-stable distributions and in the description of their domain of attraction by Gnedenko [START_REF] Gnedenko | Sur la distribution limite du terme maximum d'une série aléatoire[END_REF] and de Haan [START_REF] De Haan | On regular variation and its application to the weak convergence of sample extremes[END_REF]. Another stimulating point of view developped by Lamperti [START_REF] Lamperti | On extreme order statistics[END_REF] is to introduce a time variable and consider the asymptotic behavior of the partial maxima

max 1≤i≤[nu] Z i -b n a n , u ≥ 0.
The corresponding limit process is known as an extremal process: it is a pure jump Markov process wich is also max-stable, see Resnick [START_REF] Resnick | Weak convergence to extremal processes[END_REF].

Since then, EVT has known many developments. Among several other directions, the extension of the theory to multivariate and spatial settings is particularly important, as well as the statistical issues raised by the applications on real data sets. For excellent reviews of such developments, the reader is invited to refer to the monographies by Resnick [START_REF] Resnick | Extreme values, regular variation and point processes[END_REF], de Haan and Fereira [START_REF] De Haan | Extreme value theory[END_REF] or Beirlant et al [START_REF] Beirlant | Statistics of extremes[END_REF] and the references therein.

Our purpose here is to focus on the functional framework and investigate the asymptotic behavior of the partial maxima processes based on a doubly infinite array of independent random processes. Let T be a compact metric space and, for each n ≥ 1, let {X in , i 1} be independent copies of a sample continuous stochastic process (X n (t)) t∈T . Without loss of generality, we will always suppose that X n is non-negative (otherwise consider X n (t) = e Xn(t) ) and we denote by C + = C + (T ) be the set of non-negative continuous functions on T . We are mainly interested in the process of pointwise maxima

M n (t) = max{X in (t), 1 i n}, t ∈ T, (1)
and the process of partial maxima

Mn (u, t) = max{X in (t), 1 i [nu]}, u ≥ 0, t ∈ T. (2)
We use the convention that the maximum over an empty set is equal to 0, so Mn (0, t) ≡ 0.

Clearly, we also have Mn (1, t) = M n (t). In this framework, the parameter t ∈ T is thought as a space parameter and u ∈ [0, +∞) as a time variable.

Our approach relies on the convergence of the following empirical measures

β n = n i=1 δ X in and βn = i≥1 δ (X in ,i/n) (3)
on C + and C + ×[0, +∞) respectively. The maxima processes M n and Mn can be written as functionals of the empirical measure β n and βn respectively. We will prove the continuity of the underlying functionals and then use the continuous mapping Theorem to deduce convergence of the maxima processes from convergence of empirical measures.

Connections between EVT and point processes are well known. In the case when the underlying state space is locally compact, the following result holds (see Proposition 3.21 in [START_REF] Resnick | Extreme values, regular variation and point processes[END_REF]): if nP[X n ∈ • ] vaguely converges to some measure µ, the empirical measures β n and βn converge to Poisson random measures with intensity µ and µ ⊗ respectively, being the Lebesgue measure on [0, +∞). However, in our framework the state space C + is not locally compact and we need a suitable generalization of the above result. To this aim, we follow the approach by Davis and Mikosch [START_REF] Davis | Extreme value theory for space-time processes with heavy-tailed distributions[END_REF] based on the notion of boundedly finite measures and -weak convergence detailed by Daley and Vere-Jones in [START_REF] Daley | An introduction to the theory of point processes[END_REF].

The paper is organized as follows. In section 2, we introduce the technical material needed on boundedly finite measures, -weak convergence and convergence of empirical measures. Our main result is the convergence of the partial maxima process Mn which is stated and proved in section 3. Then, in section 4, we investigate some properties of the limit process M , known as a superextremal process. A brief extension of our results to further order statisctics is considered in section 5. The last section is devoted to an application of our results to the class of log-normal processes, based on a recent work by Kabluchko [START_REF] Kabluchko | Extremes of independent gaussian processes[END_REF] on the extremes of Gaussian processes.

Preliminaries on boundedly finite measures and point processes

In this section, we present general results on boundedly finite measures and point processes that will be useful in the sequel. The reader should refer to Appendix 2.6 in Daley and Vere-Jones [START_REF] Daley | An introduction to the theory of point processes[END_REF] or to section 2 of Davis and Mikosch [START_REF] Davis | Extreme value theory for space-time processes with heavy-tailed distributions[END_REF]. This has also closed connections with the theory of regular variations, see Hult and Lindskog [START_REF] Hult | Extremal behavior of regularly varying stochastic processes[END_REF][START_REF] Hult | Regular variation for measures on metric spaces[END_REF].

Boundedly finite measures

Let (E, d) be a complete separable metric space and E the Borel σ-algebra of E. We denote by M b (E) the set of all finite measures on (E, E). A sequence of finite measures {µ n , n 1} is said to converge weakly to µ ∈ M b (E) if and only if f dµ n → f dµ for all bounded continuous function f on E. We write µ n w → µ to denote weak convergence. It is well known that this notion of convergence is metrized by the Prokhorov metric

p(µ 1 , µ 2 ) = inf{ε > 0; ∀A ∈ E, µ 1 (A) µ 2 (A ε ) + ε and µ 2 (A) µ 1 (A ε ) + ε} where A ε = {x ∈ E; ∃a ∈ A, d(a, x) < ε} is the ε-neighborhood of A. Endowed with this metric, M b (E) is a complete separable metric space.
A measure µ on (E, E) is called boundedly finite if it assigned finite measure to bounded sets, i.e. µ(B) < ∞ for all bounded B ∈ E. Let M b (E) be the set of all boundedly finite measures µ on (E, E). A sequence of boundedly finite measures {µ n , n

1} is said to converge -weakly to µ ∈ M b (E) if and only if f dµ n → f dµ for all bounded continuous function f on E with bounded support. There exists a metric p on M b (E) that is compatible with this notion of -weak convergence and that makes M b (E) a complete and separable metric space. Such a metric can be constructed as follows. Fix an origin e 0 ∈ E and, for r > 0, let Br = {x ∈ E; d(x, e 0 ) ≤ r} be the closed ball of center e 0 and radius r. For any µ 1 , µ 2 ∈ M b (E), let µ

p (µ 1 , µ 2 ) = ∞ 0 e -r p r (µ (r) 1 , µ (r) 2 ) 1 + p r (µ (r) 1 , µ (r) 2 )
dr.

There are several equivalent characterizations of -weak convergence. Let {µ n , n 1} and µ be boundedly finite measures on E. The following statements are equivalent: i ) f dµ n → f dµ for all f bounded continuous real valued function on E with bounded support;

ii ) p (µ n , µ) → 0;

iii ) there exist a sequence r k +∞ such that, for all k 1, µ

(r k ) n w → µ (r k ) ; iv ) µ n (B) → µ(B) for all bounded B ∈ E such that µ(∂B) = 0.
We write µ n w → µ to denote -weak convergence.

A boundedly finite point measure on E is a measure µ of the form

µ = i∈I δ x i
with {x i , i ∈ I} a finite or countable family of points in E such that any bounded set B ⊆ E contains at most a finite number of the x i 's. The set of boundedly finite point measures on E is denoted by M (b,p) (E). It is a closed subset of M b (E) and hence is complete and separable when endowed with the induced metric.

The following characterization of -weak convergence of boundedly finite point measures will be useful. Let {µ n = i δ x n i , n 1} and µ = i δ x i be elements of M (b,p) (E). Then µ n w → µ if and only if there exist some sequence r k +∞ such that for all k ≥ 1, there exist N 0 and m 0 such that

∀n N, µ (r k ) n = m i=1 δ x n,k i and µ (r k ) = m i=1 δ x k i for some points {x n,k i , x k i ; 1 ≤ i ≤ m, n ≥ N } in Br k such that x n,k i → x k i .

Convergence of the empirical measures

Let (Ω, F, P) be a probability space. A boundedly finite point process on E is a measurable mapping N : Ω → M (b,p) (E). A typical example of boundedly finite point process on E is a Poisson point process Π ν with intensity ν ∈ M b (E). We will also consider the following empirical point processes β n and βn defined by Equation ( 3), where, for each n ≥ 1, {X in , i ≥ 1} are independent copies of an E-valued random variable X n . The random variable β n is a finite point process on E, while βn is a boundedly finite point process on

Ẽ = E × [0, +∞) endowed with the metric d((x 1 , u 1 ), (x 2 , u 2 )) = d(x 1 , x 2 ) + |u 2 -u 1 |, x 1 , x 2 ∈ E, u 1 , u 2 ∈ [0, +∞).
The following Proposition will play a key role in the sequel.

Proposition 2.1. The following statements are equivalent, where the limits are taken as n → ∞:

i ) nP[X n ∈ • ] w → ν;
ii ) β n ⇒ Π ν with Π ν a Poisson point process on E with intensity µ and ⇒ standing for weak convergence in M (b,p) (E).

iii ) βn ⇒ Πν with Πν Poisson point process on E × [0, +∞) with intensity µ ⊗ , being the Lebesgue measure on [0, +∞) and ⇒ standing for weak convergence in

M (b,p) (E × [0, +∞)).
The proof of this Proposition follows by an adapation of Proposition 3.21 in Resnick [START_REF] Resnick | Extreme values, regular variation and point processes[END_REF] which states a similar result in the case of a locally compact state space and in terms of vague convergence. As noticed by Davis and Mikosch [START_REF] Davis | Extreme value theory for space-time processes with heavy-tailed distributions[END_REF], the proof remains valid for a complete separable metric space E if we change vague convergence by -weak convergence. See also Theorem 4.3 in Davydov, Molchanov and Zuyev [START_REF] Davydov | Strictly stable distributions on convex cones[END_REF].

In the terminology of Hult and Lindskog [START_REF] Hult | Extremal behavior of regularly varying stochastic processes[END_REF][START_REF] Hult | Regular variation for measures on metric spaces[END_REF], the condition i) in Proposition 2.1 means that the sequence X n is regularly varying. A particularly important case is when E is endowed with a strucure of cone, i.e. a multiplication by positive scalars. If there exists a sequence (a n ) n≥1 of positive reals such that

nP[a -1 n X 1 ∈ • ] w → ν
for some nonzero ν ∈ M b (E), then the random variable X 1 is said to be regularly varying.

Under some technical assumptions on the structure of the convex cone E (e.g. continuity properties of the multiplication), the limit measure ν is proved to be homogeneous of order -α < 0, i.e. ν(λA) = λ -α ν(A), λ > 0, A ∈ E.

Furthermore, the sequence a n is regularly varying with index 1/α. In this framework, if {X i , i ≥ 1} is an i.i.d. sequence of regularly varying random variables and

X in = a -1 n X i , then condition i) is equivalent to nP[a -1 n X 1 ∈ • ] w → ν
where the limit measure ν is homogeneous of order -α < 0 and

M n (u, t) = a -1 n max 1≤i≤[nu] X i (t), u ≥ 0, t ∈ T.

Extremes of independent stochastic processes

We go back to our original problem where the X in 's are non negative continuous processes on T . We endow the set C + = C + (T ) of non negative continuous functions on T with the uniform norm x = sup t∈T |x(t)|, x ∈ C + . It turns out that, when working with maxima, the uniform metric is not adaptated, mainly because sets such as {x ∈ C + ; x ≥ ε} are not bounded for this metric. For this reason, we introduce C + 0 = (0, +∞] × S C + where S C + = {x ∈ C + : x ∞ = 1} is the unit sphere. We define the metric

d((r 1 , s 1 ), (r 2 , s 2 )) = |1/r 1 -1/r 2 | + s 1 -s 2 , (r 1 , s 1 ), (r 2 , s 2 ) ∈ C + 0 .
The metric space (C + 0 , d) is complete and separable. Define also C + 0 = C + \ {0} and consider the "polar decomposition":

T : C + 0 → (0, ∞) × S C + x → ( x , x/ x ) .
The mapping T is an homeomorphism and we identify in the sequel C + 0 and (0, +∞)×S C + . In this metric, a subset B of C + 0 is bounded if and only if it is bounded away from zero, i.e. included in a set of the form {x ∈ C + ; x ≥ ε} for some ε > 0.

Spatial maximum process

For the sake of clarity, we present first our results on the spatial maximum process M n defined by Equation (1).

Theorem 3.1. Assume that nP[X n ∈ • ] w → ν in M b (C + 0 ) and that ν(C + 0 \ C + 0 ) = 0.
Then, the process M n weakly converges in C(T ) as n → ∞ to the process M defined by

M (t) = sup i≥1 Y i (t), t ∈ T, (4)
where i≥1 δ Y i is a Poisson point process on C + 0 with intensity ν.

The proof of Theorem 3.1 relies on the following Lemma. With a slight abuse of notation, we define

M (b,p) (C + 0 ) = {µ ∈ M (b,p) (C + 0 ) : µ(C + 0 \ C + 0 ) = 0}. Note M (b,p) (C + 0 ) is an open subset of M (b,p) (C + 0 ).
Lemme 3.1. The mapping θ : M (b,p) (C + 0 ) → C + defined by θ(0) ≡ 0 and

θ i∈I δ (r i ,s i ) = t → sup{r i s i (t); i ∈ I}
is well-defined and continuous.

Proof of Lemma 3.1:

• First, we show that the mapping θ is well defined. This is not obvious since the pointwise supremum of a countable family of functions is not necessarily finite nor continuous.

For ε > 0, let µ (ε) be the restriction of µ on the set Bε = [ε, ∞] × S C . Since Bε is bounded for the metric d, the point measure µ (ε) has only a finite number of atoms and therefore θ(µ (ε) ) is the maximum of a finite number of non negative continuous functions and is hence a non negative continuous function.

Furthermore, for all t ∈ T , |θ(µ (ε) )(t) -θ(µ)(t)| ε. So θ(µ) is the uniform limit as ε → 0 of the continuous functions θ ε (µ), and hence θ(µ) ∈ C + and satisfies

θ(µ (ε) ) -θ(µ) ε. (5)
• Second, we show that the mapping θ is continuous.

Let {µ n = i∈In δ x n i , n 1} be a sequence of measures in M (b,p) (C + 0 ) converging to µ = i∈I δ x i . For all ε > 0, there exists ε < ε, N 0 and m 0 such that

∀n N, µ (ε ) n = m i=1 δ x nε i and µ (ε ) = m i=1 δ x ε i with x nε i , x ε i ∈ Bε and x nε i → x ε i as n → ∞.
Clearly, this implies that θ(µ

(ε ) n ) converges in C + to θ(µ (ε ) ) as n → ∞.
Hence there exists an integer n 0 such that θ(µ

(ε ) n ) -θ(µ (ε )
) ≤ ε for all n ≥ n 0 . Then, thanks to Equation (5), we obtain for n ≥ n 0 ). Note that M n = θ(β n ), and according to Lemma 3.1, the mapping θ is continuous. So the continuous mapping Theorem (see e.g. Theorem 5.1 in Billingsley [START_REF] Billingsley | Convergence of probability measures[END_REF]) entails θ(β n ) ⇒ θ(Π ν ) which is equivalent to M n ⇒ M .

θ(µ n ) -θ(µ) ≤ θ(µ n ) -θ(µ (ε ) n ) + θ(µ (ε ) n ) -θ(µ (ε ) ) + θ(µ (ε ) ) -θ(µ)) ≤ 3ε. This proves that θ(µ n ) → θ(µ)

Spatio-temporal maximum process

We consider now convergence of the space-time process Mn defined by Equation (2).

For fixed u ≥ 0, the space process t → Mn (u, t) is sample continuous and non negative, i.e. a random elements of C + . Furthermore, the time process u → Mn (u, •) can be seen as a C + -valued càd-làg process on [0, +∞); it is indeed constant on intervals of the form [k/n, k/n + 1/n), k ∈ N. Hence, we will consider the process Mn as a random element of the Skohorod space D([0, +∞), C + ) endowed with the J 1 -topology (see for example Ethier and Kurtz [START_REF] Ethier | Markov processes[END_REF] for the definition and properties of Skohorod space).

Theorem 3.2. Assume that nP[X n ∈ • ] w → ν in M b (C + 0 ) and that ν(C + 0 \ C + 0 ) = 0.
Then, the process Mn weakly converges in D([0, +∞), C + ) as n → ∞ to the superextremal process M defined by

M (u, t) = sup{Y i (t)1 [U i ,+∞) (u); i ≥ 1}, u > 0, t ∈ T, (6)
where i≥1 δ (Y i ,U i ) is a Poisson point process on C + 0 × [0, +∞) with intensity ν ⊗ .
For the proof, we will need the following analoguous of Lemma 3.1 in the space-time framework. Let M b,p (C + 0 × [0, +∞)) be the subset of measures µ ∈ M (b,p) (C

+ 0 × [0, +∞)) such that µ((C + 0 \ C + 0 ) × [0, +∞)) = 0 and define C = µ ∈ M (b,p) (C + 0 × [0, +∞)); µ(C + 0 × {u}) ≤ 1 for all u ≥ 0 .
In other words, a measure µ = i∈I δ (r i ,s i ,u i ) belongs to C if and only if r i < +∞ for all i ∈ I and the u i 's are pairwise distinct.

Lemme 3.2. The mapping θ :

M (b,p) (C + 0 ×[0, +∞)) → D([0, +∞), C + ) defined by θ(0) ≡ 0 and θ i∈I δ (r i ,s i ,u i ) = (u, t) → sup{r i s i (t)1 [u i ,+∞) (u); i ∈ I}
is well defined and continuous on C.

Proof of Lemma 3.2:

The proof is similar to the proof of Lemma 3.1 and we give only the main lines.

• First we show that θ is well defined. Recall that

C + 0 × [0, +∞) is endowed with the metric d((r 1 , s 1 , u 1 ), (r 2 , s 2 , u 2 )) = d((r 1 , s 1 ), (r 2 , s 2 )) + |u 2 -u 1 |. Let µ ∈ M (b,p) (C + 0 × [0, +∞)). For ε > 0 and M > 0, let µ (ε,M ) be its restriction to Bε,M = [ε, +∞] × S C + × [0, M ].
Since Bε,M is bounded, µ (ε,M ) has only a finite number of atoms and we easily check that θ(µ (ε,M ) ) belongs to D([0, +∞), C + ). Furthermore, for u < M and t ∈ T ,

| θ(µ (ε,M ) )(u, t) -θ(µ)(u, t)| ≤ ε so that θ(µ (ε,M ) ) converges uniformly on [0, M ] × T to θ(µ) as ε → 0.
The constant M being arbitrary, this implies that θ(µ) ∈ D([0, +∞), C + ) and that the application θ is well defined. It also holds that

sup (u,t)∈[0,M ]×T | θ(µ (ε,M ) )(u, t) -θ(µ)(u, t)| ≤ ε. (7)
• Next, we show that θ is continuous on C.

Let {µ n = i∈In δ x n i ; n 1} be a sequence of M (b,p) (C + 0 × [0, +∞)) converging to µ = i∈I δ x i ∈ C. For all ε > 0 and M > 0, there exist ε < ε, M > M , N ≥ 1 and some m ≥ 1, such that ∀n N, µ (ε ,M ) n = m i=1 δ x n i and µ (ε ,M ) = m i=1 δ x i with x n i = (r n i , s n i , u n i ), x i = (r i , s i , u i ) ∈ Bε ,M and x n i → x i as n → ∞.
The condition µ ∈ C ensures that the u i 's are pairwise distinct and we can suppose without loss of generality that u 1 < • • • < u m . For n large enough, we will also have

u n 1 < • • • < u n m . Define δ M the metric associated with the J 1 -topology on D([0, M ], C + ) by δ M (x, y) = inf λ sup (u,t)∈[0,M ]×T |x(λu, t) -x(u, t)|
where the infinimum is taken over the set of non-decreasing homeomorphisms λ of [0, M ]. Since for large n the u i 's and the u n i 's are in the same relative order, there exists a non-decreasing homeomorphism

λ n M of [0, M ] such that λ n M (u n i ) = u i for all 1 ≤ i ≤ m. We then have δ M ( θ(µ (ε ,M ) n ), θ(µ (ε ,M ) )) ≤ δ M ( θ(µ n )), θ(µ (ε ,M ) n )) + δ M ( θ(µ (ε ,M ) n ), θ(µ (ε ,M ) )) + δ M ( θ(µ (ε ,M ) ), θ(µ)) ≤ sup (u,t)∈[0,M ]×T | max 1≤i≤m (r n i s n i (•)1 [λ n M u n i ,+∞) (u)) -max 1≤i≤m (r i s i (•)1 [u i ,+∞) (u))| = max 1≤i≤m r n i s n i -r i s i → 0 as n → ∞.
Hence, for sufficiently large n, δ M ( θ(µ

(ε ,M ) n
), θ(µ (ε ,M ) )) ≤ ε and thanks to Equation [START_REF] De Haan | Extreme value theory[END_REF], this entails

δ M ( θ(µ n )), θ(µ)) ≤ δ M ( θ(µ n )), θ(µ (ε ,M ) n )) + δ M ( θ(µ (ε ,M ) n ), θ(µ (ε ,M ) )) + δ M ( θ(µ (ε ,M ) ), θ(µ)) ≤ 3ε.
Since M is arbitrary large and ε arbitrary small, this proves the convergence θ(µ n ) to θ(µ) in D([0, +∞), C + ).

Proof of Theorem 3.2 :

The proof is very similar to the proof of Theorem 3.1. Note that Mn = θ( βn ) where βn is defined by Equation (3). According to Proposition 2.1, the empirical measure βn weakly converges in M (b,p) (C

+ 0 × [0, +∞)) to a Poisson point process Πν with intensity ν ⊗ . The assumption ν(C + 0 \ C + 0
) and the fact that the Lebesgue measure has no atoms ensure that Πµ lies almost surely in C. From Lemma 3.2, the mapping θ is continuous on C, so that the continuous mapping Theorem implies θ( βn ) ⇒ θ( Πν ) which is equivalent to Mn ⇒ M .

Properties of the limit process

In this section, we give some properties of the superextremal process M defined by Equation (6) in Theorem 3.2. First, we characterize its finite dimensional distributions. We use vectorial notations: for l ≥ 1, t = (t 1 , . . . , t l ) ∈ T l , y = (y 1 , . . . , y l ) ∈ [0, +∞) l and u > 0, we write M (u, t) ≤ y if and only if M (u, t i ) ≤ y i for all i ∈ {1, . . . , l}.

Proposition 4.1. For k 1, 0 = u 0 < u 1 < • • • < u k , t ∈ T l and y 1 , . . . , y k ∈ [0, +∞) l , it holds P M (u 1 , t) y 1 , ..., M (u k , t) y k = k i=1 exp[-(u j -u j-1 )ν(A j )] with A j = {f ∈ C + 0 ; ∃i ∈ 1, . . . , l, f (t i ) ≥ min k≥j y k,i }.
Proof of Proposition 4.1 : Let j ∈ {1, . . . , k}. Note that M (u j , t) y j if and only if Πµ doesn't intersect the set

B j = {(f, u) ∈ C + 0 × [0, +∞); u ≤ u j and ∃i ∈ {1, • • • l} f (t i ) > y j,i } So, we have P M (u 1 , t) y 1 , ..., M (u k , t) y k = P Πν ∩ (∪ k j=1 B j ) = ∅] = exp[-(ν ⊗ )(∪ k j=1 B j )]. The B j 's are not pairwise disjoint. To compute the measure (ν ⊗ )(∪ k j=1 B j ), we observe that ∪ k j=1 B j = ∪ k j=1 ([u j-1 , u j ) × A j ) ∪ ({u k } × A k )
where the sets in the right hand side are pairwise disjoint. From this, we deduce

(ν ⊗ )(∪ k j=1 B j ) = k j=1 (u j -u j-1 )ν(A j ).
This proves the Proposition.

Next we prove that the process u → M (u, •) is a C + -valued homogeneous Markov process. Let F u be the σ-algebra generated by { M (s, t); s ∈ [0, u], t ∈ T }. The symbol ∨ stands for poinwise maximum. Proposition 4.2. Let u ≥ 0. The conditional distribution of ( M (u + h, •)) h≥0 given F u is equal to the distribution of ( M (u, •) ∨ M (h, •)) h≥0 where M is an independent copy of M .

In some sense, this Proposition states that the process M has "independent and stationary increments" with respect to the maximum: for 0 = u 0 < u 1 < • • • < u k , the distribution of ( M (u i , •)) 1≤i≤k is equal to the distribution of

(∨ i j=1 M j (u j -u j-1 , •)) 1≤i≤k ,
where M 1 , . . . , M k are i.i.d. copies of M . This property is similar to the property of independence and stationarity of increments (with respect to the addition) of Lévy processes. For a fixed point t ∈ T , the process {M (u, t), u ≥ 0} is known as an extremal process (see Proposition 4.7 of Resnick [START_REF] Resnick | Extreme values, regular variation and point processes[END_REF]).

Proof of Proposition 4.2:

Consider the decomposition Πµ = Π[0,u] ν ∪ Π(u,∞) ν where Π[0,u] ν = Πν ∩ (C + 0 × [0, u]) and Π(u,∞) ν = Πν ∩ (C + 0 × (u, ∞)).
• Similarly, the self-similarity is proven as follows:

P M (cu j , t) y j , ∀j ∈ {1, . . . , k}

= k i=1 exp[-c(u j -u j-1 )ν(A j )] = k i=1 exp[-(u j -u j-1 )ν(c -1/α A j )]
= P c 1/α M (u j , t) y j , ∀j ∈ {1, . . . , k} .

Further order statistics

The point process approach for extremes is powerful: the convergence of the empirical measure β n entails not only the convergence of the maxima M n but also of all the order statistics, and a similar result hold for the space-time version of these processes.

For 1 ≤ r ≤ n and t ∈ T , define M r n (t) as the r-th largest value among X 1n (t), . . . , X nn (t). Note M 1 n (t) = M n (t) is simply the maximum. For r > n, we use the convention M r n (t) = 0. The process M r n = (M r n (t)) t∈T is refered to as the r-th order statistic of the sample {X 1n , . . . , X 1n }. Similarly, for r ≥ 1, t ∈ T and u ≥ 0 define M r n (u, t) as the r-th largest value among X 1n (t), . . . , X [nu]n (t) with the convention M r n (u, t) = 0 if r > [nu]. An alternating definition in terms of the empirical measure βn is given by

M r n (u, t) = sup{y ≥ 0; βn (B y u,t ) ≥ r} with B y u,t = {(f, v) ∈ C + 0 × [0, +∞); f (t) ≥ y, v ≤
u} and the convention that the supremum over an empty set is equal to zero.

With these notations, we can strengthen Theorem 3.2 as follows:

Theorem 5.1. Assume that nP[X n ∈ • ] w → ν in M b (C + 0 ) and that ν(C + 0 \ C + 0 ) = 0. Then, for each r ≥ 1, the joint order statistics ( M 1 n , . . . , M r n ) weakly converges in D([0, +∞), (C + ) r ) as n → ∞ to ( M 1 , . . . , M r ) defined by M j (u, t) = sup{y ≥ 0; Πν ∩ B y u,t ≥ r}, 1 ≤ j ≤ r, u ≥ 0, t ∈ T,
with Πν a Poisson point process on C + 0 × [0, +∞) with intensity ν ⊗ .

• Proof of i): Let r > 0 and r = ln r. We have

nP[sup t∈T X n (t) > r] = nP[sup t∈T Y n (t) > r] = n √ 2πb n e b 2 n /2 R e -w-w 2 /2b 2 n P[sup t∈T Y w n (t) > r]dw ≤ M R e -w P[sup t∈T Y w n (t) > r]dw, so that sup n≥1 nP[sup t∈T X n (t) > r] ≤ M R e -w sup n≥1 P[sup t∈T Y w n (t) > r]dw. (13)
First, we have

+∞ 0 e -w sup n≥1 P[sup t∈T Y w n (t) > r]dw +∞ 0 e -w dw < ∞. (14)
By tightness of the family {Y w n -µ w n ; w ∈ R, n ≥ 1} (see Lemme 6.1), there exists c 1 > 0 such that for all n ≥ 1 and w ∈ R,

P[sup t∈T (Y w n (t) -µ w n (t)) > c 1 ] < 1 2 .
Furthermore, Equations ( 11)- [START_REF] Hult | Regular variation for measures on metric spaces[END_REF] and assumption ii) of Theorem 6.1 together imply that there are some c 2 , c 3 > 0 and n 0 ≥ 1 such that for all n ≥ n 0 and all w < 0, sup t∈T µ w n (t)

Y w n (t) > r] < 2ψ - r -w/2 -c 1 -c 2 c 3 ≤ 2e -( r-w/2-c 1 -c 2 c 3 ) 2 , w < 0,
where ψ is the tail of the standard Gaussian distribution. Consequently,

0 -∞ e -w sup n≥n 0 P[sup t∈T Y w n (t) > r]dw < 2 0 -∞ e -w e -( r-w/2-c 1 -c 2 c 3 ) 2 dw < ∞. (15)
Equations ( 14) and ( 15

In view of Equation ( 13), this proves i).

• Proof of ii): For δ > 0, we have We will apply Lebesgue dominated convergence Theorem to prove that this has limit zero as δ → 0. Lemma 6.1 implies that for every fixed w, the sequence {e Y w n ; n 1} is tight so that lim Equation ( 16) with r = ln ε provides the required domintation condition.

• Proof of iii):

We have to prove that for all set A ⊆ C + 0 of the form with r = min 1≤j≤k y j . Equation ( 16) yields the required domination condition.

A = {f ∈ C + 0 ; ∃j ∈ [[1, k]], f (t j ) > x k }, k ≥ 1,

  restriction of µ 1 and µ 2 to Br . Note that µ 1 and µ 2 are finite measures on Br and denote by p r the Prokhorov metric on M b ( Br ). Define

  as n → ∞ and that the mapping θ is continuous. Proof of Theorem 3.1: Under the assumption nP[X n ∈ • ] w → ν, we know from Proposition 2.1 that the empirical measure β n defined by Equation (3) weakly converges in M (b,p) (C + 0 ) to a Poisson point process Π ν with intensity ν. The assumption ν(C + 0 \ C + 0 ) ensures that Π ν lies almost surely in M (b,p) (C + 0

  ) > r]dw < ∞.

  e Y w n (δ) ε]dw.

δ→0 e -w sup n≥1 P

 n≥1 [ω e Y w n (δ) ε] = 0.It remains to prove a suitable domination condition: sinceδ → P[ω e Y w n (δ) ε] is a nondecreasing function, it is enough to prove that for some δ 0 > 0 R e -w sup n≥1 P[ω e Y w n (δ 0 ) ε]dw < ∞.Sinceω e Y w n (δ 0 ) ≤ exp[sup t∈T Y w n (t)],we have e -w sup n≥1 P[ω e Y w n (δ 0 ) ε] e -w sup n≥1 P[sup t∈T Y w n (t) ln ε].

  t 1 , . . . , t k ∈ T, x 1 , . . . , x k ∈ (0, +∞),Letting y j = ln x j , we have nP[X n ∈ A] = nP[∃j ∈ [[1, k]], X n (t j ) > x j ] = nP[∃j ∈ [[1, k]], Y n (t j ) > y j ]We will apply once again Lebesgue's dominated convergence Theorem. Lemma 6.2 entailsP[∃j ∈ [[1, k]], Y w n (t j ) > y j ] → P[∃j ∈ [[1, k]], w + W (t j ) -[[1, k]], we W (t j )-1 2 Γ(t j ,t 0 ) > x j ]w -2 dw = ν(A).We are left to prove a suitable domination condition: for all n ≥ 1, we have

						1 2	Γ(t j , t 0 ) > y j ],
	so that we expect
	nP[X n ∈ A] →	e -w P[∃j ∈ [[1, k]], w + W (t j ) -∞ R	1 2	Γ(t j , t 0 ) > y j ]dw
	= P[∃j ∈ e 0 -w-w 2 2b 2
	it holds
				lim
	=	√	n 2πb n e b 2 n /2	R	e -w-w 2 2b 2 n P[∃j ∈ [[1, k]], Y w n (t j ) > y j ]dw
	∼		e	-w-w 2 2b 2
		R

n→∞ nP[X n ∈ A] = ν(A). n P[∃j ∈ [[1, k]], Y w n (t j ) > y j ]dw. n P[∃j ∈ [[1, k]], Y w n (t j ) > y j ] ≤ e -w sup n≥1 P[sup t∈T Y w n (t) > r]

and the two terms in the r.h.s. are independent. It is easily seen that the first term is equal to M (u, •), and that the invariance of the Lebesgue measure implies that the second term has the same distribution as M (h, •). This proves the Proposition.

In the particular case when the measure ν is homogeneous of order -α < 0, the process M enjoys further interesting properties. Proposition 4.3. Suppose ν is homogeneous with index -α < 0. Then:

• M is self-similar with index 1 α , i.e. for all c > 0, the rescaled process ( M (cu, •)) u≥0 has the same distribution as (c 1/α M (cu, •)) u≥0

Proof of Proposition 4.3 :

We check that the the two processes have the same finite dimensional distributions. To this aim, we use Proposition 4.1 and the notations of the proposition.

• Let M 1 , . . . , M n be independent copies of M . By Proposition 4.1 and the homogeneity of µ,

M i (u j , t) y j , ∀j ∈ {1, . . . , k} = P M (u j , t) y j , ∀j ∈ {1, . . . , k}

This proves the max-stability.

The limit process M j corresponds to the r-th space-time order statistic associated to the point process Πν .

Proof of Theorem 5.1 : Using once again Proposition 2.1 and the continuous mapping Theorem, it is enough to prove the following generalization of Lemma 3.2: the mapping

defined by θ(0) ≡ 0 and

is well defined and continuous on C. The proof is very similar to the proof of Lemma 3.2 and the details are omitted for the sake of brevity.

Gaussian case

In geostatistics, Gaussian processes are often used as spatial models. In the finitedimensional setting, the asymptotic behavior of the maxima of Gaussian random vectors was first investigated by Huesler and Reiss [START_REF] Hüsler | Maxima of normal random vectors: between independence and complete dependence[END_REF]. Recently, Kabluchko, Schlather and de Haan [START_REF] Kabluchko | Stationary max-stable fields associated to negative definite functions[END_REF] and Kabluchko [START_REF] Kabluchko | Extremes of independent gaussian processes[END_REF] consider the functional setting and prove convergence of the maxima of i.i.d. centered Gaussian processes under suitable conditions on their covariance structure. The family of Brown-Resnick processes has quickly become a popular model for spatial extremes due to its simple characterization by a negative definite functions.

Our purpose here is to revisit their results and apply the present framework: we put the emphasis on regular variations and point processes. In the sequel, we follow the approach by Kabluchko, see Theorems 2 and 6 in [START_REF] Kabluchko | Extremes of independent gaussian processes[END_REF].

We suppose that (T, d) is a compact metric space satisfying the entropy condition

where N (ε) is the smallest number of balls of radius needed to cover T . Let Z n = {Z n (t); t ∈ T }, n ≥ 1, be a sequence of continuous zero-mean unit-variance Gaussian processes with covariance function

Define the scaling sequence

and the rescaled process

We consider the log-normal process

Theorem 6.1. Fix t 0 ∈ T and suppose that: i ) Uniformly in t 1 , t 2 ∈ T :

ii ) For all t 1 , t 2 ∈ T , there exists C > 0 that

Then,

and {W (t), t ∈ T } a centered Gaussian process such that W (t 0 ) = 0 and with incremental variance Γ.

Following [START_REF] Kabluchko | Extremes of independent gaussian processes[END_REF][START_REF] Kabluchko | Stationary max-stable fields associated to negative definite functions[END_REF], we define Y w n the process Y n conditioned by Y n (t 0 ) = w and we note µ w n (t) = E[Y w n (t)]. We will need the following two Lemmas.

Lemme 6.1. Under conditions i) and ii) of Theorem 6.1, the family {Y w n ; n ≥ 1} is tight in C(T ) for all fixed ω ∈ R, as well as the family {Y w n -µ w n ; w ∈ R, n ≥ 1}.

Lemma 6.1 follows from the proof of Theorem 6 in [START_REF] Kabluchko | Extremes of independent gaussian processes[END_REF] (see also the proof of Theorem 17 in [START_REF] Kabluchko | Stationary max-stable fields associated to negative definite functions[END_REF]). Details are omitted here.

Under assumption i), it holds

This implies that Var(Y w n (t 2 ) -Y w n (t 1 )) → Γ(t 1 , t 2 ) and proves the Lemma.

Proof of Theorem 6.1 :

Our proof relies on the following criterion for -convergence in C + 0 by Hult et Lindskog. According to Theorem 4.4 of [START_REF] Hult | Regular variation for measures on metric spaces[END_REF], it is enough to prove that: i) for every r > 0, sup n≥1 nP[sup t∈K X n (t) > r] < ∞;

ii) for every ε > 0, lim is bounded by some constant M > 0.


Marwa Hamza

Abdelhamid Hassairi
email: abdelhamid.hassairi@fss.rnu.tn

Conjugate Priors for Exponential Families Having Cubic Variance Functions

Keywords: Natural exponential family, variance function, cumulant function, Monge-Ampère equation, prior distribution

Conjugate Priors for Exponential Families Having Cubic Variance Functions

Maroua Hamza,

Introduction and preliminaries

For the convenance of the reader, we first introduce some notations and recall some facts concerning the natural exponential families and their variance functions, our notations are the ones used in [START_REF] Letac | Lectures on natural exponential families and their variance function[END_REF]. Let E be a linear vector space with finite dimension n, denote by E * its dual, and let E * × E → IR : (θ, x) → θ, x be the duality bracket. If µ is a positive radon measure on E, then

L µ (θ) = E exp( θ, x )µ(dx) (1.1)
denotes its Laplace transform. We also denote by M(E) the set of measures µ such that the set

Θ(µ) = interior{θ ∈ E * ; L µ (θ) < +∞} (1.2)
is non empty and µ is not concentrated on an affine hyperplane of E. The cumulant function of an element µ of M(E) is the function defined for θ in Θ(µ) by k µ (θ) = log L µ (θ).

To each µ in M(E) and θ in Θ(µ), we associate the probability distribution on E P (θ, µ)(dx) = exp ( θ, xk µ (θ)) µ(dx). The set

F = F (µ) = {P (θ, µ); θ ∈ Θ(µ)}
is called the natural exponential family (NEF) generated by µ. We also say that µ is a basis of F .

The function k µ is strictly convex and real analytic. Its first derivative k ′ µ defines a diffeomorphism between Θ(µ) and its image M F . Since k ′ µ (θ) = the probability of F with mean m, we have

F = {P (m, F ); m ∈ M F } ,
which is the parametrization of F by the mean. Now the covariance operator of P (m, F ) is denoted by V F (m) and the map

M F -→ L s (E * , E); m -→ V F (m) = k ′′ µ (ψ µ (m))
is called the variance function of the NEF F . It is easy proved that for all m ∈ M F ,

V F (m) = (ψ ′ µ (m)) -1 ,
and an important feature of V F is that it characterizes F in the following sense: If F and F ′ are two NEFs such that V F (m) and

V F ′ (m) coincide on a nonempty open set of M F ∩ M F ′ , then F = F ′ .
Now, let us examine the influence of an affine transformation and a power convolution on a NEF

F = F (µ). If ϕ(x) = a(x) + b, where a ∈ GL(E) and b ∈ E, is an affine transformation of E, then ϕ(F (µ)) = F (ϕ(µ)), M ϕ(F ) = ϕ(M F ), and
V ϕ(F ) (m) = a V F (ϕ -1 (m)) a * ,
where a * is the transpose of a. On the other hand the set

Λ(µ) = {λ > 0; ∃ µ λ ∈ M(E) such that L µ λ (θ) = (L µ (θ)) λ for all θ ∈ Θ(µ)}
is called the Jorgensen set of µ and the measure µ λ is the λ-power of convolution of µ.

For λ in Λ(µ), we have that

M F (µ λ ) = λM F , and V F λ (m) = λ V F ( m λ ).
A very interesting fact is that the most common real and multivariate probability distributions belong to the natural exponential families such that the variance function is a polynomial of degree less then or equal to three in the mean m. For instance, up to affine transformations and power of convolution (up to the type), the Gaussian, Poisson, gamma, binomial, negative binomial and hyperbolic cosine distributions form the class of all real NEF's whose variance function is a polynomial of degree less than or equal to 2 characterized by Morris [START_REF] Morris | Natrual exponential families with quadratic variance-function[END_REF]. Letac and Mora [START_REF] Letac | Natural real exponential familiies with cubic variance functions[END_REF] have added six others types of distributions, namely, the inverse Gaussian, Ressel, Abel, Tackàs, strict arcsine and large arcsine, to get the class of real cubic NEF's, that is the class of NEF's such that variance function is a polynomial of degree less than or equal to three. The classification of NEF's with polynomial variance function have been extended to the multivariate NEF's. The multivariate version of the Morris class, called the class of simple quadratic NEF's, has been completely described by Casalis [START_REF] Casalis | The 2d+4 simple natrual exponential families on IR d[END_REF], it contains 2n + 4 types. Hassairi [6] has defined and characterized the so-called class of multivariate simple cubic NEF's which is the natural extension of the class of real cubic NEF's. It is worth mentioning here that the simple quadratic NEF's are not the only families which have quadratic variance functions, the Wishart families on symmetric matrices have also quadratic variance functions. The classifications of NEF's by the form of the variance function provide an important tool in the study of distributions. In fact, in many important cases, the variance function is very simple and is easier to use than the distribution itself or the Laplace transform. Moreover, the fact that the variance function is quadratic or cubic, is not only a question of form, but the form corresponds to some very interesting analytical characteristic properties. In this respect, let us mention that for the Morris class of real quadratic NEF's, we have the Meixner characterization based on some families of orthogonal polynomials which generate exactly the Morris class (see [START_REF] Meixner | Orthogonal Polynomsysteme mit einer besonderen Gestalt der erzengenden function[END_REF]). Another characterization due to Feinsilver [START_REF] Feinsilver | Some classes of orthogonal polynomials associated with martingales[END_REF] states that a certain class of polynomials naturally associated to a NEF is orthogonal if and only if the family is in the Morris class. This characterization has been extended to the Casalis class of simple quadratic NEF's by Labeye-Voisin, and Pommeret [START_REF] Labeye-Voisin | Polynmes orthogonaux associs aux familles exponentielles de[END_REF]. Concerning the cubic NEF's, Hassairi and Zarai [START_REF] Hassairi | Characterization of the cubic exponential families by orthogonality of polynomials[END_REF] introduced a notion of 2-orthogonality for a sequence of polynomials to give an extended version of the Meixner and Feinsilver characterization which subsume the Letac-Mora class of real cubic NEF's. Hassairi and Zarai [START_REF] Hassairi | Characterization of the simple cubic multivariate exponential families[END_REF] have also introduced a notion of trans-orthogonality for a sequence of multivariate polynomials to extend their characterization result to the class of multivariate simple cubic NEF's. Besides these characterizations based essentially on different notions of orthogonality of polynomials, it is stated in Casalis [START_REF] Casalis | The 2d+4 simple natrual exponential families on IR d[END_REF] that the simple quadratic NEF satisfies a property based on two conjugates families of prior distributions related to the NEF. For a NEF F = F (µ), consider the family of prior distributions Π introduced by Diaconis and Ylvisaker [START_REF] Diaconis | Conjugate priors for exponential families[END_REF] and defined by

Π = {π t,m 0 (dθ) = C t,m 0 exp t( m 0 , θ -k µ (θ)) 1 Θ(µ) (θ)dθ, t > 0, m 0 ∈ M F } (1.3)
where C t,m 0 is a normalizing constant. Consider also the family Π * introduced by Consonni et al [START_REF] Consonni | Reference priors for exponential families with simple quadratic variance function[END_REF], see also [START_REF] Gutierrez-Pena | Reference priors for exponential families[END_REF] and defined by

Π * = {π * t,m 0 , t > 0, m 0 ∈ M F } (1.4) where π * t,m 0 (dm) = C * t,m 0 exp t( m 0 , ψ µ (m) -k µ (ψ µ (m)) ) 1 M F (m)dm,
and the constant C * t,m 0 is a normalizing constant. Then, when F (µ) is a Wishart or a simple quadratic NEF, we have that k ′ µ (Π) = Π * . It is also shown that this property is equivalent to two other properties expressed in terms of some differential equations satisfied by the cumulant function k µ . In the real case, the property characterizes the Morris class of real quadratic NEF's, that is k ′ µ (Π) = Π * if and only if the NEF F is in the Morris class. In the present paper, we extend these results to the class of multivariate simple cubic NEF's. We construct two families of prior distributions related to a multivariate NEF, and we show that these families coincide when the NEF is simple cubic. We then show that this property is equivalent to the fact that the cumulant function is a solution of some Monge-Ampère equation and also equivalent to the fact that the variance function satisfies a differential equation. As a corollary, we obtain three new characterizations of the Letac-Mora class of real cubic NEF's.

Some equivalent properties

Throughout this section, we suppose that F = F (µ) is a NEF on a linear vector space E with dimension n. Besides the family Π of prior distributions defined in (1.3), we introduce another family Π of prior distributions. Let β be in E * such that the set

Θ = {θ ∈ Θ(µ); 1 + β, k ′ µ (θ) > 0}
is nonempty, and denote M = k ′ µ ( Θ). Consider the family of prior distributions

Π = { π t,m 0 ; t ∈ IR * + , m 0 ∈ M F }, (2.5)
where

π t,m 0 (dm) = C t,m 0 (1 + β, m ) -n-2 exp t{ m 0 , ψ µ (m) -k µ (ψ µ (m))} 1 M (m)dm.
With these notations, we next state and prove our first main result.

Theorem 2.1 The three following properties are equivalent (1) There exists (a, b, c)

∈ E × IR 2 such that for all m in M F , det V F (m) = (1 + β, m ) n+2 exp{ a, ψ µ (m) + bk µ (ψ µ (m)) + c}.
(2) There exists (a, b) ∈ E × IR such that for all m in M F and any basis, (e i ) n i=1 of E, with dual basis (e * i ) n i=1 , we have

n i=1 [V ′ F (m)(e i )]e * i = n + 2 1 + β, m V F (m)(β) + a + bm. ( 2

.6)

(3) There exists an open subset

Ω of IR * + × M F such that k ′ µ (Π) = Π Ω = { π t,m 0 ; (t, m 0 ) ∈ Ω}.
Note that (1) may be stated in terms of the cumulant function as there exists (a, b, c) ∈ E × IR 2 such that for all θ in Θ(µ),

det k ′′ µ (θ) = (1 + β, k ′ µ (θ) ) n+2 exp{ a, θ + bk µ (θ) + c},
that is the cumulant function is solution of some Monge-Ampère equation (see [START_REF] Zuily | Local existence and regularity of the Dirichlet problem for the Monge-Ampère equation, Journes "équations aux Drives Partielles[END_REF]). Proof We will show that (1) ⇔ (2) and ( 1) ⇔ (3).

(1) ⇒ (2) Suppose that V F (m) satisfies ( 1), then we have

log det V F (m) = (n + 2) log(1 + β, m ) + { ψ µ (m), a + b k µ (ψ µ (m)) + c}.
Taking the derivative, we get

trace(V -1 F (m)V ′ F (m)(.)) = (n + 2) β, . 1 + β, m + ψ ′ µ (m)(.), a + b m, ψ ′ µ (m)(.) , which is equivalent to n i=1 [V ′ F (m)(.)V -1 F (m)(e i )](e * i ) = (n + 2) β, . 1 + β, m + ψ ′ µ (m)(.), a + b m, ψ ′ µ (m)(.) .
Replacing (.) by V F (m)(.), and using the condition of symmetry [START_REF] Letac | Lectures on natural exponential families and their variance function[END_REF], page 103), we obtain

V ′ (m)(V (m)(α))(β) = V ′ (m)(V (m)(β))(α) ∀ α, β ∈ E * (see
n i=1 [V ′ F (m)(e i )(.)](e * i ) = (n + 2) β, V F (m)(.) 1 + β, m + a, (.) + b m, (.) . (2.7) As V F (m) is symmetric, we get n i=1 [V ′ F (m)(e i )](e * i ) = (n + 2) 1 + β, m V F (m)(β) + a + bm.
(2) ⇒ (1) Suppose that (2) holds. Then, we easily get (2.7). Replacing, in (2.7), (.) by V -1 F (m)(.), one obtains

n i=1 V ′ F (m)(e i )V -1 F (m)(.)(e * i ) = (n + 2) β, (.) 1 + β, m + a, ψ ′ µ (m)(.) + b m, ψ ′ µ (m)(.) .
This is equivalent to

trace(V -1 F (m)V ′ (m)(.)) = (n + 2) β, (.) 1 + β, m + a, ψ ′ µ (m)(.) + b m, ψ ′ µ (m)(.) .
Integrating, we deduce that there exists c in IR such that

log det(V F (m)) = (n + 2) log(1 + β, m ) + { ψ µ (m), a + b k µ (ψ µ (m)) + c},
and the result follows.

(1) ⇒ (3) Suppose that (1) holds, and define

Ω = {(t, m 0 ) ∈ IR * + × M F ; t > b and m 0 ∈ (1 - b t )M F - a t }.
Take (t, m 0 ) in Ω and denote ν the image of π t,m 0 by ψ µ . Then it is easy to verify that

ν(dθ) = C t,m 0 e c exp{ tm 0 + a, θ -(t -b)k µ (θ)} 1 Θ (θ)dθ. Since (t, m 0 ) is in Ω, we have that t -b > 0 and tm 0 + a t -b ∈ M F .
Thus taking t 1 = tb and m 1 = tm 0 + a tb , we obtain that

ν(dθ) = C t 1 ,m 1 exp{t( m 1 , θ -k µ (θ))} 1 Θ (θ)dθ.
Hence ψ µ ( Π Ω ) ⊂ Π, and it follows that Π Ω ⊂ k ′ µ (Π). Conversely, if π t,m 0 is an element of Π, then its image σ by k ′ µ is given by

σ(dm) = C t,m 0 e -c (1 + β, m ) -n-2 exp{ tm 0 -a, ψ µ (m) -(t + b)k µ (ψ µ (m))} 1 M (m)dm.
Taking

t 1 = t + b and m 1 = tm 0 -a t + b . Then (t 1 , m 1
) is in Ω, and we have

σ(dm) = C t 1 ,m 1 (1 + β, m ) -n-2 exp t 1 { m 1 , ψ µ (m) -k µ (ψ µ (m))}1 M (m)dm ,
which is an element of Π Ω .

(3) ⇒ (1) Suppose that k ′ µ (Π) = Π Ω . Then, for an element π t,m 0 of Π, we have on the one hand,

k ′ µ (π t,m 0 )(dm) = (det V F (m)) -1 C t,m 0 exp t{ m 0 , ψ µ (m) -k µ (ψ µ (m))} 1 M (m)dm.
On the other hand, since

k ′ µ (π t,m 0 ) is in Π Ω , there exists (t 1 , m 1 ) in Ω such that k ′ µ (π t,m 0 )(dm) = C t 1 ,m 1 (1 + β, m ) -n-2 exp t 1 { m 1 , ψ µ (m) -k µ (ψ µ (m))} 1 M (m)dm. Comparing these two expressions of k ′ µ (π t,m 0 ) gives det V (m) = (1 + β, m ) n+2 exp{ a, ψ µ (m) + bk µ (ψ µ (m)) + c}, where a = tm 0 -t 1 m 1 , b = t 1 -t, and c = log( c t,m 0 c t 1 ,m 1
).

3 Characterizations of the Letac-Mora class of real cubic NEFs

In this section, we prove that a multivariate simple cubic NEF satisfies the properties in Theorem (2.1), and that the real versions of these properties characterize the real cubic NEFs. Recall that a simple cubic NEF is obtained form a simple quadratic NEF by the so-called action of the linear group GL(IR × E) on the NEFs of E. For more details, we refer the reader to [START_REF] Hassairi | La classification des familles exponentielles naturelles sur IR n par l'action du groupe linaire de IR n+1[END_REF], where a complete description of this class is given. This action is in fact an extension of the way in which the Letac-Mora class of real cubic NEFs is obtained from the Morris class of real quadratic NEF's. For our purposes here, we need only to mention that, up to affine transformations and power of convolution, a simple cubic variance function is of the form

V (m) = (1 + β, m ) (I + m ⊗ β) V 1 ( m 1 + β, m ) (I + β ⊗ m), (3.8)
where V 1 is the variance function of a simple quadratic NEF F 1 , and m is in (M F 1 ) β , where

(M F 1 ) β = {m ∈ M F 1 ; 1 + β, m > 0 and m 1 + β, m ∈ M F 1 }.
The relation (3.8) is invertible and conversely, we have

V 1 (M ) = (1 -β, M ) (I -M ⊗ β) V ( M 1 -β, M ) (I -β ⊗ M ), (3.9)
where M is in (M F ) -β . We also mention that the relation between a simple cubic NEF F (µ) and a simple quadratic NEF F (ν) may also be expressed in terms of the cumulant functions by

     k µ (λ) = k ν (θ) -k 0 λ = -βk ν (θ) + θ -λ 0 (3.10)
or equivalently by

     k ν (θ) = k µ (λ) -k 1 θ = βk µ (λ) + λ -θ 1 (3.11)
where (k 0 , λ 0 ) and (k 1 , θ 1 ) are constants in IR × E. Note that if β = 0 in (3.8) we obtain the simple quadratic class. Then for more accuracy we exclude this case and we keep only β in E * \ {0}.

We now prove that the multivariate simple cubic NEF's satisfy the properties in Theorem 2.1.

Proposition 3.1 Let F = F (µ) be a simple cubic NEF on E, then there exists (a, b, c) in E * × IR 2 such that det(V F (m)) = (1 + β, m ) n+2 exp{ ψ µ (m), a + b k µ (ψ µ (m)) + c}.
Proof Given that the family F is simple cubic, then there exist β in E * and F 1 = F (ν) a simple quadratic NEF such that

V F (m) = (1 + β, m ) (I + m ⊗ β) V F 1 ( m 1 + β, m ) (I + β ⊗ m), see (3.8). As det(I + m ⊗ β) = 1 + β, m , we obtain det(V F (m)) = (1 + β, m ) n+2 det(V F 1 ( m 1 + β, m
)).

We now use the fact for a simple quadratic NEF F 1 (see [START_REF] Casalis | The 2d+4 simple natrual exponential families on IR d[END_REF]), there exist a ′ in E * and b ′ , c ′ in IR such that, for all As the Letac-Mora class of real cubic NEFs is nothing but the simple cubic class, when the dimension n is equal to 1, this class satisfies the real version of the properties in Theorem 2.1. We will show that, in this case, these properties are characteristic. Solving this differential equation by standard methods gives

M in M F 1 , det(V F 1 (M )) = exp{ a ′ , ψ ν (M ) + b ′ k ν (ψ ν (M )) + c ′ }.
V F (m) = λ (1 + βm) 3 - b β 2 (1 + βm) 2 + b -βa 2β 2 (1 + βm),
which is a polynomial of degree less then or equal to 3. 2

E

  xP (θ, µ)(dx), M F is called the domain of the means of F . The inverse function of k ′ µ is denoted by ψ µ and setting P (m, F ) = P (ψ µ (m), µ)

  It follows that det(V F (m)) = (1 + β, m ) n+2 exp{ ψ ν ( m 1 + β, m ), a ′ + b ′ k ν (ψ ν ( m 1 + β, m)) + c ′ }.From (3.10), putting λ = ψ µ (m) and θ = ψ ν ( m 1 + β, m), we getk ν (ψ ν ( m 1 + β, m )) = k µ (ψ µ (m)) + k 0 , ψ ν ( m 1 + β, m ) = ψ µ (m) + βk µ (ψ µ (m)) + βk 0 + λ 0 . Then det V F (m) = (1 + β, m ) n+2 exp{ ψ µ (m), a ′ + (b ′ + a ′ , β )k µ (ψ µ (m)) +(b ′ k 0 + a ′ , βk 0 + λ 0 + c ′ )}.Setting a = a ′ , b = b ′ + a ′ , β and c = b ′ k 0 + a ′ , βk 0 + λ 0 + c ′ , we obtain the desired result. 2

Theorem 3 . 2

 32 Let F = F (µ) be a NEF on the real line, then F is cubic if and only ifk ′ µ (Π) = Π.Proof Suppose that k ′ µ (Π) = Π. Then according to Theorem(2.1), the variance function V F (m) satisfies the differential equation(1 + βm)V ′ F (m) -3β V F (m) = (a + bm)(1 + βm).