James L Crowley

Patrick Reignier

Rémi Barraquand

Situation Models: A Tool for Observing and Understanding Activity

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

I. INTRODUCTION

Human activity is extremely rich. Real world scenes can contain an overwhelming number of possible agents and objects to detect and observe. As are result, systems and services based on observation of activity must, either implicitly or explicitly, be able to choose where to look next and what to look for. Designers of system for observing activity are increasingly confronted with the problem of control attention.

Attention is not the only problem confronting designers of systems for observing activity. Activity in the real world often occurs in less than ideal viewing conditions. Poor lighting, background clutter, object texture, and occlusions can degrade the reliability of even the most well designed systems. Thus systems and services must be able to detect and discard uncertain and unreliable observations, and if appropriate, substitute default information. In addition, many services require real time information from perception. In such systems it may be preferable to provide an immediate response with default information and to use background processes to verify that the response was Manuscript received March 9, 2009

correct.

Current systems for observing activities tend to be constructed in an ad-hoc manner with control structures that are hard-wired into the system design. Such systems are generally restricted to detecting a very small set of activities observed within a highly controlled environment. Adapting such systems to different operating environments or modifying such systems to observe different forms of activity can involve extensive reprogramming.

In this paper we propose an approach for constructing systems for observing activity based on a model from Cognitive Science. We propose the use of situation models to organize, control, and interpret perception of activity. We will first provide some background from Cognitive Science concerning the use of situation models as a model of human mental activity. We then describe how to adopt such a model into software systems that provide services. We propose a layered, component-oriented software architecture for building situation aware services, and examine how situation models can be used to structure perceptual components and to provide default information for understanding activity. We conclude with discussion of the problems of automatically acquiring situation models through developmental learning.

II. SITUATION MODELS FOR MENTAL PROCESSES

Situation models have been proposed by Johnson-Laird [START_REF] Johnson-Laird | Mental Models: Towards a Cognitive Science of Language, Inference, and Consciousness[END_REF], as a cognitive theory for human mental models. Over the last 25 years, Theories about situation models have been adopted and developed by a large community of cognitive psychologists. Key publications include [START_REF] Radvansky | The retrieval of situation-specific information[END_REF], [START_REF] Zwaan | Situation Models in Language Comprehension and Memory[END_REF] as well as [START_REF] Johnson-Laird | Mental models[END_REF].

Situations are defined as a set of relations between entities, where an entity is anything that can be observed, and a relation is a predicate function. According to Radansky [START_REF] Radvansky | The retrieval of situation-specific information[END_REF], a situation model is a mental representation of a described or experienced situation in a real or imaginary world.

Situation models are commonly composed of four primary types of information: 1) A spatial-temporal framework (spatial locations, time frames) 2) Entities (people, objects, ideas, etc. ) 3) Properties of entities (color, emotions, goals, shape, etc. ) 4) Relational information (spatial, temporal, causal, ownership, kinship, social, etc. )

Situation Models: A Tool for Observing and Understanding Activity James L. Crowley, Member, IEEE, Patrick Reignier and Remi Barranquand Situation models can be structured along dimensions of space, time, causality, actors and objects. Extensions of situations models have been proposed to represent intentions of actors. It is commonly assumed that both general world knowledge (knowledge about concept types, e.g., scripts, schemas, categories, etc ) and referent specific knowledge (knowledge about specific entities, independent of the situation) are used in constructing situation models.

Situation models are used for representations of: 1) Information about events.

2) Information about sequences of events.

3) Information about collections of episodes

We have adapted the concept situation models to construct system and services based on monitoring and observing human activity [START_REF] Crowley | Context Driven Observation of Human Activity[END_REF], [START_REF] Crowley | Learning Situation Models for Understanding Activity In The 5th International Conference on Development and Learning[END_REF], [START_REF] Brdiczka | Learning situation models for providing context-aware services[END_REF]. Although most of our implementations have been constructed using smart environments, such services can also be designed using robotic systems. Indeed, our approach to smart environments is to see the environment as a form of "inside out" robot, observing and interacting with occupants. Thus we maintain that models for understanding activity in smart environments may also be adapted for construction of autonomous robots.

III. SITUATION MODELS FOR OBSERVING ACTIVITY

Situation models can be used to addresses the twin problems of focus of attention, and operation with unreliable, erroneous or missing data. They can also be used to decouple services from the time constraints normally imposed by near (or not-so near) real time vision systems. We present our technique in the context of a service-oriented architecture constructed using a layered, component-based, software model. For the robotics and vision communities, these concepts may require some explanation.

The term "service" is used here in its most general form. Generally, it will refer to the "services" that informatics systems, including robotic systems, provide to people. User services can be designed as software agents that observe human activity and respond to events. Over the last few years, we have constructed a variety of user services based on dynamic assembly of perceptual components. This include services for lecture recording [START_REF] Metze | The FAME Interactive Space[END_REF], meeting services [START_REF] Danninger | MyConnector: analysis of context cues to predict human availability for communication[END_REF], monitoring of the health and well-being of elderly, and availability monitoring [START_REF] Brdiczka | Learning situation models for providing context-aware services[END_REF]. As sensor and actuator technology mature, we can expect to see the emergence of an increasing variety of such systems for domestic services (cleaning, logistics, cooking), commercial services (shopping, queue management, customer assistance), health monitoring and assisted living, security monitoring, and a variety of other application domains. All of these examples require observing and understanding the actions of humans. We believe that situation models will provide an important component of such systems.

We note that the term "service oriented" also has a more technical meaning for the software engineering community. In software engineering, a "service oriented" system is one in which software components interact according to a welldefined contract. For example, a location service integrates information from a variety of sources to estimate the current location of a user. Although the two uses of the term "service" are not incompatible, they can cause some confusion. Thus we will use the explicit term "Software services" for services that are primarily designed to interact with software components. We will interchangeably use "user services" or simply "services" for systems that interact with and assist people.

Modern software systems are generally designed using a layered architecture. A layered architecture organizes the system into a hierarchy of interchangeable components, with well-defined interfaces. The design and operation at a particular layer may proceed independently of the underlying components. Components that make up a particular layer may be reused or shared by a variety of services. Components that are temporarily inoperative may be replaced with components alternative components. A common example of this approach is provided by the current generation of location aware services on mobile devices that can interchangeably use location information from GPS, cell phone repeaters, or WIFI repeater identity. Components for providing location from WIFI, GPS or cell-phone repeaters are a form of "perceptual component" that operate in parallel using competing methods to make available a key piece of information: current location. We propose a similar approach to building components for observing activity. Perceptual components can be constructed to observe a scene with competing methods to provide information that may then be shared between different services.

A situation model falls naturally at the interface between user services and perceptual components. For user services, the situation model provides a default reasoning system that can complete or repair partial or missing information from sensing. For the perceptual components, the situation model can be used to focus attention on the objects and events that are relevant to a service, allowing irrelevant objects or events to be ignored. The situation model can be used to predict possible events based, both to focus attention, and to prime a response before the event occurs.

In the following, we describe layered architecture for observation systems and present a component-oriented approach to building observation systems. We then provide a formal definition for a situation model, and describe how such a model can be used to configure and control perceptual components as well as focus attention, predict events, and provide default reasoning for observation of activity.

A. Services, Sensors, and Components

We are interested in services that provide assistance through the observation of human activity. A service determines requirements for perception and action, without specifying how these requirements are to be met. Hardwiring the interconnection between sensor signals and actuators is possible, and can provide simplistic services that are hardware dependent and have limited utility.

Separating services from their underlying hardware makes it possible to build systems that operate in a larger range of environments, for a larger variety of functions. However such separation requires that the sensor-actuator layer provide logical interfaces, or standard API's, that are function centered and device independent. Hardware independence and generality require abstractions for perception and action. A layered architecture of user services is shown in figure 1. At the lowest layer, the service's view of the world is provided by a collection of physical sensors and actuators. This corresponds to the sensor-actuator layer. This layer depends on the technology and encapsulates the diversity of sensors and actuators by which the system interacts with the world. Information at this layer is expressed in terms of sensor signals and device commands.

Service abilities for perception and action are provided by components for perception and action. Components make observations about the world, interact with users as well as to make changes to the environment. In our systems, services maintain information about users and environment in a situation model. The situation model has the form of a network of situations. Each situation has three facets: Observation, Reaction and Prediction. The observation facet specifies the entities and relations needed to define the situation. This can act as a specification that serves to activate a set of perception components capable of providing observations about the required entities and their relations. The reaction facet specifies how the service should behave in each situation, including both the desired state of the environment, and a specification communications that the service should make with the user. The Prediction facet indicates possible changes to the current situation, by pointing to adjacent situations and describing the events that can indicate the change.

Sensors are devices that make measurements, ranging simple devices that measure temperature or humidity, to devices that capture motion (infrared motion detectors), acoustic energy (microphones) and images (cameras).

Actuators impart change on the environment. Such devices can range from information displays, control of lighting and sound systems, motorized controls for doors, windows and window blinds, as well as mobile robotic devices for cleaning and entertainment.

Components for perception and action operate at a higher level of abstraction than sensors and actuators. While sensors and actuators operate on device-specific signals, perception and action operate in terms of environmental state. Perception interprets sensor signals by detecting, recognizing and observing people, things and events. Action components alter the environment to being it to a desired state. Tightly coupling perception and action can offer many advantages. Controlling action with perception allows a service to adapt action in accordance with the effect on the environment. Action can also be used to reconfigure the environment to improve perception, or even to probe the environment as part of perception.

B. Components for Perception and Action

Perception and action components are autonomous assemblies of modules executed in a cyclic manner by a component supervisor. Components communicate via synchronous data streams and asynchronous events in order to provide software services for action or perception. We propose a data-flow process architecture for software components for perception and action [START_REF]Software Process Modeling and Technology[END_REF], [START_REF] Rasure | The Khoros application development environment[END_REF], [START_REF] Crowley | Integration and Control of Reactive Visual Processes[END_REF]. Component based architectures, as described in Shaw and Garlan [START_REF] Shaw | Software Architecture: Perspectives on an Emerging Disciplines[END_REF], are composed of auto-descriptive functional components joined by connectors. Such an architecture is well adapted to interoperability of components, and thus provides a framework by which multiple partners can explore design of specific components without having to rebuild the entire system.

Components are controlled by a supervisory module. The component supervisor interprets commands and parameters, supervises the execution of the transformation, and responds to queries with a description of the current state and capabilities of the component. The auto-critical report from modules allows a component supervisor to monitor the execution time and to adapt the schedule of modules for the next cycle so as to maintain a specified quality of service, such as execution time or number of targets tracked. Such monitoring can be used, for example, to reduce the resolution of processing by selecting 1 pixel of N [START_REF] Piater | Event-based Activity Analysis in Live Video using a Generic Object Tracker[END_REF] or to selectively delete targets judged to be uninteresting or erroneous [START_REF] Hall | An automatic approach for parameter selection in self-adaptive tracking[END_REF]. In addition to recognition the supervisory component provides execution scheduling, self-monitoring, parameter regulation, and communications. The supervisor also acts as a scheduler, invoking execution of modules in a synchronous manner. For self-monitoring, a component applies a model of its own behavior to estimate both quality of service and confidence for its outputs. Monitoring allows a process to detect and adapt to degradations in performance due to changing operating conditions by reconfiguring its component modules and operating parameters. Monitoring also enables a process to provide a symbolic description of its capabilities and state.

Homeostasis or "autonomic regulation of internal state" is a fundamental property for robust operation in an uncontrolled environment. A component is auto-regulated when processing is monitored and controlled so as to maintain a certain quality of service. The process supervisor maintains homeostasis by adapting module parameters to maximize estimated quality service. For example, processing time and precision are two important state variables for a tracking process. Quality of service measures such as cycle-time, number of targets, or precision can be maintained by dropping targets based on a priority assignment or by changing resolution for processing of some targets.

During the communication phase, the supervisor may respond to requests from other components. These requests may ask for descriptions of process state, process capabilities, or may provide specification of new recognition methods. The supervisor acts as a programmable interpreter, receiving snippets of code script that determine the composition and nature of the process execution cycle and the manner in which the process reacts to events. Recognition procedures are small procedures interpreted by a lightweight language interpreter for the lisp-like language "scheme" [START_REF] Lux | The Imalab Method for Vision Systems[END_REF]. In our implementation, such procedures may be preprogrammed, or they may be downloaded to the component during configuration as snippets of code.

For most human activities, there are a potentially infinite number of entities that could be observed and an infinite number of possible relations for any set of entities. The appropriate entities and relations must be determined with respect to the service to be provided. This is the role of the situation model, as described in the previous section. The situation model allows the system to focus perceptual attention and computing resources, in order to associate the current state of the activity, with the appropriate system action.

Perceptual components communicate using Streams, Events, and Queries. Streams are synchronous communication channels for communicating synchronous data such as image frames or acoustic signals. An important role for perceptual components is to process streams in order to observe entities and their properties. Events are asynchronous messages generated by components in response to changes in entities or their properties. Events may be sent to other components or to the situation model. Queries are communication transactions in which a service, the situation model, or another component interrogate a component about its entities and their properties.

C. Assembling Components to Provide Services

We have constructed a middle-ware environment [START_REF] Emonet | O3MiSCID: an Object Oriented Opensource Middleware for Service Connection, Introspection an Discover[END_REF] that allows us to dynamically launch and connect components on different machines. This environment, called O3MiSCID, provides an XML based interface that allows components to declare input command messages, output data structures, as well as current operational state. In this environment, a user service may be created by assembling a collection of perceptual components.

Available components are discovered by interrogating an ontology server. An open research challenge is to provide an ontological system for indexing components based on function in a manner that is sufficiently general to capture future functionalities as they emerge. In addition the ontology server is used to establish compatible communications of data. The problem of aligning ontologies of data structures is manageable when components are co-designed by a single group. The problem becomes very difficult when components are developed independently with no prior effort to agree on specifications. This problem resembles to web-ontology alignment problem that currently receives attention in software engineering.

Figure 3 shows a simple example of a service provided by an assembly of perceptual components. This service integrates information from multiple cameras to provide 3-D target tracking. A set of tracked entities is provided by a Bayesian 3D tracking process that tracks targets in 3D scene coordinates. This process specifies the predicted 2-D Region of Interest (ROI) and detection method for a set of pixel-level detection components. These components use color, motion or background difference subtraction to detect and track blobs in an image stream from a camera. The O3MICID middle-ware makes it possible to dynamically add or drop cameras to the process during tracking. Fig. 3. An example of an assembly of perceptual components. The 3D Bayesian blob tracker provides a ROI and detection method for a number of 2D entity detection components. The result is used to update a list of 3D blobs.

D. Entities and Relations

Situations are defined as relations between entities. An "entity" is anything that can be observed. This solipsistic viewpoint admits that the system can only see what it knows how to see. At same time, it sidesteps existential dilemmas related to how to define notions of "object" and "class". In this view, a chair is anything that can be used as a chair, regardless of its apparent form.

Formally, entities are correlated sets of observations. Entities are grounded in the software components for observation of activity, typically through some form of tracking process that correlates observations over time. Entities can be decorated with properties that make possible the determination of relations between entities.

A relation is a predicate or binary function computed on the properties of one or more entities. Relations have an arity, that specifies the number of properties of entities that are concerned. An arity-1 relation is true when a property is observed to be within some range of values, or is otherwise signaled as true by a sensor. Examples can include (standing person) or (running person).

Relations of Arity-2 may be any of the classical relation spaces including spatial relations, temporal relations or more abstract functions such as social-behaviour. Spatial relations can be 2D or 3D and relative or absolute, depending on the requirements of the service. Examples can include absolute position of actors (at podium person), (seated-at table person), relative position (facing person1 person2), or even refer to the posture of persons (standing person). Observing human interaction can require perceptual components that detect more abstract social behaviour, such as (talking-to person1 person2) or (smiling-at person1 person2).

As mentioned above, a key problem is that the number of potential relations that might be observed is an unbounded set. The situation model for a service specifies the entities (agents and objects) that must be observed, the properties that are required, and the relations between entities that are relevant. The task of the system designer is to provide perceptual components that can detect and track the required entities, measure the required properties, and detect when the required relations are true.

Human Attention is an important relation in social situations. In our approach, we adopt the attention model developed by Maisonnasse [START_REF] Maisonnasse | Attentional Model for Perceiving Social Context in Intelligent Environments, Artificial Intelligence Apllications and Innovations[END_REF]. In this work, attention is defined as the cognitive process of selectively concentrating on one aspect of the environment while ignoring other things. We include attention of agents as one of the fundamental relations for describing the social situation.

E. Generalizing with Roles

In most situations, the exact identity of the entity is not important. Thus we have generalized situation models by the introducing of the concept of "role" [START_REF] Crowley | Context Driven Observation of Human Activity[END_REF]. A role is a form of abstract model for an entity. In applying a situation model to describe a scene, a system will select from available entities to determine which entity can "fill" each role.

Operationally, a role is an abstract generalization for a class of entities. Role classes are typically defined based on the set of actions that entities in the class can take (actors), or the set of actions that the entities can enable (props). Formally, role is a function that selects an entity from the set of observed entities.

A "role" is NOT an intrinsic property of an entity, but rather, is an interpretation assigned to an entity by the system. Entities are assigned to roles by a role assignment process. Role assignment generally occurs by applying a set of tests to available entities. The role assignment process acts as a form of "filter" [START_REF] Brdiczka | Automatic Detection of Interaction Groups[END_REF] that sorts suitability of entities based on their properties. The most suitable entity wins the role assignment.

In our experiments for automatic learning of situation models [START_REF] Crowley | Learning Situation Models for Understanding Activity In The 5th International Conference on Development and Learning[END_REF], we have discovered that roles provide generalization in learning. Reactions learned for a situation composed of one set of entities can be used to understand a different set of entities.

F. Situations

The situation model acts as a non-linear script for interpreting activity and predicting the corresponding appropriate and inappropriate actions for services. This framework organizes the observation of interaction using a hierarchy of concepts: scenario, situation, role, entity and relations. A situation is defined as a configuration of relations over a set of entities playing roles. Thus a situation is a form of state, expressed as a logical expression (a conjunction of predicates). This logical expression is composed of predicates whose arguments are roles. This concept generalizes and extends the common practice of defining situations based on the relative position of actors and objects.

Relations test the properties of entities that have been assigned to roles. As mentioned above, situations also predict possible future situations. This is captured by the connectivity of a situation network. Changes in the logical expression of relations or in the selection of entities playing roles are represented as changes in situation. Such changes can trigger system actions.

A situation is a form of state, expressed as a logical expression (a conjunction of predicates). Situations are defined as a set of relations between entities, where entities may be agents, objects or any abstract concept observed as a correlated set of properties. This logical expression is composed of predicates. This concept generalizes and extends the common practice of defining situations based on the relative position of actors and objects. As predicates, relations are Boolean functions with one or more arguments. Relations are true or false, depending on the properties of their arguments.

G. Situation Models to Understand Activity

Situations are organized into networks, with transition probabilities, so that possible next situations may be predicted from the current situation.

In our systems, the situation model drives focus of attention by specifying the entities and relations that should be attended. When a service is initiated, a list of relevant entities and relations are provided, along with the relevant configuration information. This list is used to initiate and configure the perceptual and action components needed to maintain the situation model.

Each situation contains a list of expected relations, as well as expected observed entities and their expected properties. Transitions between situations can be triggered by events, and do not require verification for the entire set of relations, entities and properties. Thus it is possible for a situation to provide default values for relations, and properties that have not been verified. When interrogated by a service, a situation model may respond with the default values, whether or not these values can be currently verified. Such a response can be provided without waiting for an actual verification to occur. However, this verification can be used as an integrity check for the situation model When a system responds with a default value, it is good practice for the system to query the relevant perceptual components to verify the default values that have been returned. In some cases, this may indicate a divergence between the situation model and the observed properties of entities. Such a divergence can be used to trigger a diagnostic process to recover from the current error, by adapting perception to changes in the environment or by developing the situation model by adding new situations or behaviours.

IV. 5. LEARNING SITUATION NETWORKS

We distinguish the concepts of adaptation from development [START_REF] Coutaz | Context is Key[END_REF]. Adaptation allows a system to maintain consistent behaviour across variations in operating environments. The environment denotes the physical world (e.g., in the street, lighting conditions), the user (identification, location, goals and activities), social settings, and computational, communicational and interactive resources. Development refers to the acquisition of abilities, in this case encoded as situation models composed of the entities, roles and relations with which situation is described and service actions are performed.

Systems for providing services based on observing activity must both adapt and develop. Adaptation is necessary to maintain consistent behaviour while accommodating changes in the operating environment, task, user population, preferences or some other factors. At the same time, human activity is too complex to be fully captured in a pre-programmed situation model. An activity model must develop through observation and interaction with users. A fundamental challenge is to provide both automatic adaptation and automatic development without disruption.

Current learning technologies, such as hidden Markov models and neural networks, require large sets of training data -something that is difficult to obtain for an extensible environment. Non-disruptive development of context models requires new ways of looking at learning, and may ultimately require a new class of minimally supervised learning algorithms. This requires that learning be studied as part of a semi-autonomous system. It requires that systems have properties of self-description, self-evaluation and auto-regulation, and may well lead to new classes of learning algorithms specifically suitable to developing and evolving context models in a non-disruptive manner.

We are currently experimenting with techniques for adapting activity models based on pre-defined stereotypical situations [START_REF] Barraquand | Learning Polite Behavior with Situation Models[END_REF]. We are exploring different approaches to learning for development of activity models starting from a predefined stereotypical model using feedback about the system actions. Because the different components of the model (entities, roles, relations, and situations) depend on each other, these cannot be developed simultaneously. Thus we have focused on the development of the situation networks and the associated system actions.

Bayesian models (in particular Hidden Markov Models [START_REF] Rabiner | A Tutorial on Hidden Markov Models and selected Applications in Speech Recognition[END_REF] as well as algorithms based on first-order logic [START_REF] Quinlan | Learning Logical Definitions from Relations[END_REF] can be used to represent and adapt the situation network. However, these approaches do not have desirable properties concerning the extension of the number of situations. Bayesian models require a large amount of example data to extend the number of states. First-order logic algorithms cannot create new predicates (problem of higher order logic), which is necessary for the extension of situations. Thus we propose an approach for changes in the structure of the situation network, as shown in figure 4.

The input to the algorithm is a predefined situation network along with feedback from prior use mediated by a supervisor. The supervisor corrects, deletes or preserves the actions executed by the system while observing a user in the environment. Each correction, deletion, or preservation generates a training example for the learning algorithm containing current situation, roles and configuration of relations, and the (correct) (re)action. The differences between the actions given in the training examples and the actions provided in the predefined situation network will drive the different steps of the algorithm.

Initially, our approach has been to directly modify system actions using the existing situation network. If action A is associated with situation S, and all training examples indicate that action B must be executed instead of A, then B is associated to S and the association between A and S is deleted. V. CONCLUSIONS Activity models for context aware services can be expressed as a network of situations concerning a set of roles, entities and relations. Roles are abstract classes for entities. Entities may be interpreted as playing a role, based on their current properties. Relations between entities playing roles define situations. This conceptual framework provides default reasoning, focus of attention, and real time response for services that require observation of human activity. This model can also provide a basis for adaptation and development of non-disruptive software services for aiding human-to-human interaction.

Socially aware observation of activity and interaction is a key requirement for development of non-disruptive services. For this to become reality, we need methods for robust observation of activity, as well as methods to automatically learn about activity without imposing disruptions. The framework and techniques described in this paper are intended as a foundation for such observation.

Fig. 1 .

 1 Fig. 1. A layered model for systems that observe human activity.

Figure 2 .

 2 Figure 2. An example of perceptual component based on visual tracking

Fig 4 :

 4 Fig 4: Overview of the algorithm for adapting system actions

  . This work was supported in part by project ANR CASPER, as well as the European IST projects FAME (IST 2000-28323), CAVIAR (IST 2001-37540) and CHIL (IST 506909) James L. Crowley is Professor at Grenoble National Polytechnique Institute (INPG) and directs the PRIMA Research group at INRIA Grenoble Centre de Recherche, 655 Ave de l'Europe, 38334 St. Ismier, France. Patrick Reignier is a Junior Professor (MdC) at the University Joseph Fourier in Grenoble, and member of the PRIMA Research group at INRIA Grenoble Centre de Recherche, 655 Ave de l'Europe, 38334 St. Ismier, France. Remi Barraquand is a doctoral student at Grenoble National Polytechnique Institute (INPG) under the direction of James L. Crowley and member of the PRIMA Research group at INRIA Grenoble Centre de Recherche.


Isabelle Charlier

Davy Paindaveine

Jérôme Saracco

QuantifQuantile : an R Package for Performing Quantile Regression through Optimal Quantization

In quantile regression, various quantiles of a response variable Y are modelled as functions of covariates (rather than its mean). An important application is the construction of reference curves/surfaces and conditional prediction intervals for Y. Recently, a nonparametric quantile regression method based on the concept of optimal quantization was proposed. This method competes very well with k-nearest neighbor, kernel, and spline methods. In this paper, we describe an R package, called QuantifQuantile, that allows to perform quantization-based quantile regression. We describe the various functions of the package and provide examples.

Introduction

In numerous applications, quantile regression is used to evaluate the impact of a d-dimensional covariate X on a (scalar) response variable Y. Quantile regression is an interesting alternative to standard regression whenever the conditional mean does not provide a satisfactory picture of the conditional distribution. Denoting by F(•|x) the conditional distribution of Y given X = x, the conditional quantile functions

x → q α (x) = inf {y ∈ R : F(y|x) ≥ α} , α ∈ (0, 1), (1)
indeed always yield a complete description of the conditional distribution. For our purposes, it is useful to recall that the conditional quantiles in (1) can be equivalently defined as

q α (x) = arg min a∈R E[ρ α (Y -a)|X = x], (2)
where ρ α (z) = αzI [z≥0] -(1α)zI [z<0] is the so-called check function.

For fixed α, the quantile functions x → q α (x) provide reference curves (when d = 1), one for each value of α. For fixed x, they provide conditional prediction intervals of the form I α = [q α (x), q 1-α (x)] (α < 1/2). Such reference curves and prediction intervals are widely used, e.g. in economics, ecology, or lifetime analysis. In medicine, they are used to provide reference growth curves for children's height and weight given their age.

Many approaches have been developed to estimate conditional quantiles. After the seminal paper of [START_REF] Koenker | Regression quantiles[END_REF] that introduced linear quantile regression, much effort has been made to consider nonparametric quantile regression. The most classical procedures in this vein are the nearest neighbor estimators [START_REF] Bhattacharya | Kernel and nearest-neighbor estimation of a conditional quantile[END_REF], the (kernel) local linear estimators [START_REF] Yu | Local linear quantile regression[END_REF] or the spline-based estimators [START_REF] Koenker | Quantile smoothing splines[END_REF][START_REF] Koenker | Penalized triograms: total variation regularization for bivariate smoothing[END_REF]. For related work, we also refer to, e.g. [START_REF] Fan | Robust nonparametric function estimation[END_REF], [START_REF] Gannoun | Reference curves based on non-parametric quantile regression[END_REF], [START_REF] Muggeo | Estimating growth charts via nonparametric quantile regression: a practical framework with application in ecology[END_REF] and [START_REF] Yu | Quantile regression: applications and current research areas[END_REF]. There also exists a wide variety of R functions/packages dedicated to the estimation of conditional quantiles. Among them, let us cite the functions rqss (only for d ≤ 2) and gcrq (only for d = 1) from the packages quantreg [START_REF] Koenker | quantreg: Quantile Regression[END_REF] and quantregGrowth [START_REF] Muggeo | quantregGrowth: Growth charts via regression quantiles[END_REF], respectively.

Recently, Charlier et al. (2015a) proposed a nonparametric quantile regression method based on the concept of optimal quantization. Optimal quantization replaces the (typically continuous) covariate X with a discretized version X N obtained by projecting X on a collection of N points (these N points, that form the quantization grid, are chosen to minimize the L p -norm of X -X N ; see below). As shown in Charlier et al. (2015b), the resulting conditional quantile estimators compete very well with their classical competitors.

The goal of this paper is to describe an R package, called QuantifQuantile (Charlier et al., 2015c), that allows to perform quantization-based quantile regression. This includes the data-driven selection of the grid size N (that plays the role of a tuning parameter), the construction of the corresponding quantization grid, the computation of the resulting sample conditional quantiles, as well as (for d ≤ 2) their graphical representation.

The paper is organized as follows. The first section briefly recalls the construction of quantizationbased quantile regression introduced in Charlier et al. (2015a,b) and explains the various steps needed to obtain the resulting estimators. The second section lists the main functions of QuantifQuantile and describes their inputs and outputs. Finally, the third section provides several examples that illustrate the use of the various functions. We conclude the paper by comparing our method with R alternatives on a real data set. An illustration of the function computing optimal quantization grids is given in the Appendix, which can be of independent interest in numerical probability, finance or numerical integration where quantization is extensively used [START_REF] Pagès | A space quantization method for numerical integration[END_REF][START_REF] Pagès | Optimal quantization methods and applications to numerical problems in finance[END_REF]).

Quantile regression through quantization

As mentioned above, the R package we describe in this paper implements the Charlier et al. (2015a,b) quantization-based methodology to perform nonparametric quantile regression. This section describes this methodology.

Approximating population conditional quantiles through quantization

Let γ N ∈ (R d ) N be a grid of size N, that is, a collection of N points in R d . For any x ∈ R d , we will denote by xγ N = Proj γ N (x) the projection of x onto this grid, that is, the point of γ N that is closest to x (absolute continuity assumption makes ties unimportant in the sequel). This allows to approximate the d-dimensional covariate X by its quantized version X γ N .

Obviously, the choice of the grid has a significant impact on the quality of this approximation. Under the assumption that X p := E[|X| p ]1/p < ∞ (throughout, | • | denotes the Euclidean norm), optimal quantization selects the grid γ N that minimizes the L p -quantization error X -X γ N p . Such an optimal grid exists under the assumption that the distribution P X of X does not charge any hyperplane, i.e. under the assumption that P X [H] = 0 for any hyperplane H [START_REF] Pagès | A space quantization method for numerical integration[END_REF]. In practice, an optimal grid is constructed using a stochastic gradient algorithm (see the following section). For more details on quantization, the reader may refer to [START_REF] Pagès | A space quantization method for numerical integration[END_REF] and [START_REF] Graf | Foundations of quantization for probability distributions[END_REF].

Based on optimal quantization of X, we can approximate the conditional quantile q α (x) in (2) by

q N α (x) := arg min a∈R E ρ α (Y -a)| X N = xN , (3)
where X N (resp., xN ) denotes the projection of X (resp., x) onto an optimal grid. It is shown in Charlier et al. (2015a) that -under mild assumptions, 1 q N α (x) converges to q α (x) as N → ∞, uniformly in x.

Obtaining an optimal N-grid

As we will see below, whenever independent copies (X 1 , Y 1 ) , . . . , (X n , Y n ) of (X , Y) are available, the first step to obtain a sample version of (3) is to compute an optimal N-grid (we assume here that N is fixed). As already mentioned, this can be done through a stochastic gradient algorithm. This algorithm, called Competitive Learning Vector Quantization (CLVQ) when p = 2, is an iterative procedure that operates as follows :

• First, an initial grid -γN,0 , say -is chosen by sampling randomly without replacement among the X i 's. • Second, n iterations are performed (one for each observation available). The grid γN,t = ( γN,t 1 , . . . , γN,t N ) at step t is obtained through

γN,t i = γN,t-1 i -δ t | γN,t-1 i -X t | p-1 γN,t-1 i -X t | γN,t-1 i -X t | if Proj γN,t-1 (X t ) = γN,t-1 i , γN,t-1 i otherwise,
where (δ t ), t ∈ N 0 , is a deterministic sequence in (0, 1) such that ∑ t δ t = ∞ and ∑ t δ 2 t < ∞. At the t th iteration, only one point of the grid moves, namely the one that is closest to X t .

The optimal grid provided by this algorithm is then γN,n .

Estimating conditional quantiles

Assume now that a sample (X 1 , Y 1 ) , . . . , (X n , Y n ) as above is indeed available. The sample analog of (3) is then defined as follows :

(S1) First, we compute the optimal grid γN,n through the stochastic gradient algorithm just described, and we write X N i = Proj γN,n (X i ), i = 1, . . . , n. (S2) Then, the conditional quantiles are estimated by

q N,n α (x) = arg min a∈R n ∑ i=1 ρ α (Y i -a)I [ X N i = xN ] , (4)
where xN = Proj γN,n (x). In practice, q N,n α (x) is simply evaluated as the sample α-quantile of the Y i 's for which

X N i = xN .
It is shown in Charlier et al. (2015a) that -under mild assumptions, q N,n α (x), for any fixed N and x, converges in probability to q N α (x) as n → ∞, provided that quantization is based on p = 2. When the sample size n is small to moderate (n ≤ 300, say), the estimated reference curves x → q N,n α (x) typically are not smooth. To improve on this, Charlier et al. (2015a,b) introduced the following bootstrapped version of the estimator in (4). For some positive integer B, generate B samples of size n from the original sample {(X i , Y i ) } i=1,...,n with replacement. From each of these bootstrap samples, the stochastic gradient algorithm provides an "optimal" grid, using these bootstrapped samples to perform the iterations. The bootstrapped estimator of conditional quantile is then

qN,n α,B (x) = 1 B B ∑ b=1 q(b) α (x), (5)
where q(b) α (x) = q(b),N,n α (x) denotes the estimator in (4) computed on the basis of the b th optimal grid.

We stress that, when computing q(b) α (x), the original sample is used in (S2); the bootstrapped samples are only used to provide the B different grids. As shown in Charlier et al. (2015a,b), the bootstrapped reference curves are much smoother than the original ones. Of course, B should be chosen large enough to make the bootstrap useful, but also small enough to keep the computational burden under control. For d = 1, we usually choose B = 50.

Selecting the grid size N

Both for the original estimators q N,n α (x) and for their bootstrapped versions qN,n α,B (x), an appropriate value of N should be identified. If N is too small, then reference curves will have a large bias, while if N is too large, then they will show much variability. This leads to the usual bias/variance trade-off that is to be achieved when selecting the value of a smoothing parameter in nonparametric statistics. Charlier et al. (2015b) proposed the following data-driven method to choose N. Let {x 1 , . . . , x J } be a set of x-values at which we want to estimate q α (x) (the x j 's are for instance chosen equispaced on the support of X) and let N be a finite collection of N-values (this represents the values of N one allows for and should typically be chosen according to the sample size n). Ideally, we would like to select the optimal value of N as

N ᾱ;opt = arg min N∈N ISE ᾱ(N), with ISE ᾱ (N) = 1 J J ∑ j=1 qN,n α,B (x j ) -q α (x j ) 2 . (6)
This, however, is infeasible, since the population quantiles q α (x j ) are unknown. This is why we draw B extra bootstrap samples (still of size n) from the original sample and consider

Nᾱ ;opt = arg min N∈N ISE ᾱ (N), with ISE ᾱ(N) = 1 J J ∑ j=1 1 B B ∑ b=1 qN,n α,B (x j ) - q( b) α (x j ) 2 , ( 7
)
where q( b) α (x j ), for b = 1 . . . , B, makes use of this bth new bootstrap sample; more precisely, the bootstrap sample is still only used to perform the iterations of the algorithm, whereas the original sample is used in both the initial grid and in (S2).

As shown in Charlier et al. (2015b) through simulations, both N → ISE ᾱ(N) and N → ISE ᾱ(N) are essentially convex in N and lead to roughly the same minimizers. This therefore provides a feasible way to select a reasonable value of N for the estimator qN,n α,B (x) in ( 5). Note that this also applies to qN,n α (x) by simply taking B = 1 in the procedure above.

If quantiles are to be estimated at various orders α, (7) will provide an optimal N-value for each α. It may then happen, in principle, that the resulting reference curves cross, which is of course undesirable.

The R Journal Vol. XX/YY, AAAA ISSN 2073-4859

One way to guarantee that no such crossings occur is to identify a common N-value for the various α's.

In such a case, N will be chosen as

Nō pt = arg min N∈N ISE ¯(N), with ISE ¯(N) = ∑ α ISE ᾱ (N), (8)
where the sum is computed over the various α-values considered. Charlier et al. (2015b) performed extensive comparisons between the quantization-based estimators in (5) -based on the efficient data-driven selection method for N just described -and some of their main competitors, namely estimators obtained from spline, k-nearest neighbor, and kernel methods. This revealed that the quantization-based estimators compete well in all cases, and often dominates their competitors (in terms of integrated square errors), particularly so for complex link functions; see Charlier et al. (2015b) for details.

Unlike the local linear and local constant estimators from [START_REF] Yu | Local linear quantile regression[END_REF], that are usually based on a global-in-x bandwidth, our quantization-based estimators are locally adaptive in the sense that, when estimating q α (x), the "working bandwidth" -that is, the size of the quantization cell containing x -depends on x. The k-nearest neighbor (kNN) estimator is closer in spirit to quantizationbased estimators but always selects k neighbors, irrespective of the x-value considered, whereas the number of X i 's in the quantization cell of x may depend on x. This subtle local-in-x behavior may explain the good empirical performances of quantization-based estimators over kernel and nearestneighbor competitors. Finally, spline methods (implemented in R with the rqss and qss functions) tend to perform poorly for complex link functions, since they always provide piecewise linear reference curves [START_REF] Koenker | Quantile smoothing splines[END_REF]. Moreover, the current implementation of the rqss function only supports dimensions 1 and 2, whereas our package allows to compute quantization-based estimators in any dimension d.

The QuantifQuantile package

This section provides a description of the various functions offered in the R package QuantifQuantile. We first detail the three functions that allow to estimate conditional quantiles through quantization. Then we describe a function computing optimal quantization grids.

Conditional quantile estimation

QuantifQuantile is composed of three main functions that each provide estimated conditional quantiles in (4)-( 5). These functions work in a similar way but address different values of d (recall that d is the dimension of the covariate vector X) :

• The function QuantifQuantile is suitable for d = 1.

• The function QuantifQuantile.d2 addresses the case d = 2.

• Finally, QuantifQuantile.d can deal with an arbitrary value of d.

Combined with the plot function, the first two functions provide reference curves and reference surfaces, respectively. No graphical outputs can be obtained from the third function if d > 2.

The three functions share the same arguments, but not necessarily the same default values. For each function, using args() displays the various arguments and corresponding default values : We now give more details on these arguments.

• X: a d × n real array (required by all three functions, a vector of length n for QuantifQuantile).

The columns of this matrix are the X i 's, i = 1, . . . , n.

• Y: an n × 1 real array (required by all three functions). This vector collects the Y i 's, i = 1, . . . , n.

• alpha: an r × 1 array with components in (0, 1) (optional for all three functions). This vector contains the orders for which q α (x) should be estimated.

• x: a d × J real array (optional for QuantifQuantile and QuantifQuantile.d2, required by QuantifQuantile.d). The columns of this matrix are the x j 's at which the quantiles q α (x j ) are to be estimated. If x is not provided when calling QuantifQuantile, then it is set to a vector of J = 100 equispaced values between the minimum and the maximum of the X i 's. If this argument is not provided when calling QuantifQuantile.d2, then the default for x is a matrix whose J = 202 = 400 column vectors are obtained as follows: 20 equispaced values are considered between the minimum and maximum values of the (X i ) 1 's and similarly for the second component. The 400 column vectors of the default x are obtained by considering all combinations of those 20 values for the first component with the 20 values for the second one 2 .

• testN: an m × 1 vector of pairwise distinct positive integers (optional for all three functions). The entries of this vector are the elements of the set N in ( 7)-( 8), hence are the N-values for which the ISE ᾱ quantity considered will be evaluated. The default is (35, 40, . . . , 55) but it is strongly recommended to adapt it according to the sample size n at hand.

• p: a real number larger than or equal to one (optional for all three functions). This is the parameter p to be used when performing optimal quantization in L p -norm.

• B: a positive integer (optional for all three functions). This is the number of bootstrap replications B to be used in (5).

• tildeB: a positive integer (optional for all three functions). This is the number of bootstrap replications B to be used when determining Nᾱ ;opt or Nō pt . • same_N: a boolean variable (optional for all three functions). If same_N=TRUE, then a common value of N (that is, Nō pt in (8)) will be selected for all α's. If same_N=FALSE, then optimal values of N will be chosen independently for the various of α (which will provide several Nᾱ ;opt , as in ( 7)).

• ncores: number of cores to use. These functions can use parallel computation to save time by increasing this parameter. Parallel computation relies on mclapply from parallel package, hence is not available on Windows.

All three functions return the following list of objects, which is of class "QuantifQuantile" :

• hatq_opt: an r × J real array (where r is the number of α-values considered). If same_N=TRUE, then the entry (i, j) of this matrix is q Nō pt ,n α i ,B (x j ). If same_N=FALSE, then it is rather q Nᾱ i ;opt ,n α i ,B (x j ). This object can also be returned using the usual fitted.values function.

• N_opt: a positive integer (if same_N=TRUE) or an r × 1 array of positive integers (if same_N=FALSE).

Depending on same_N, this provides the value of Nō pt or the vector ( Nᾱ 1 ;opt , . . . , Nᾱ r ;opt ). • hatISE_N: an r × m real array. The entry (i, j) of this matrix is ISE ᾱi (N j ). Plotting this for fixed α or plotting its average over the various α, in both cases over testN, allows to assess the global convexity of these ISEs. Hence, it can be used to indicate whether or not the choice of testN was adequate. This will be illustrated in the examples below.

• hatq_N: an r × J × m real array. The entry (i, j, ) of this matrix is qN ,n α i ,B (x j ), where N is the th entry of the argument testN. From this output, it is easy by fixing the third entry to get the matrix of the qN,n α i ,B (x j ) values for any N in testN. • The arguments X, Y, x, alpha, and testN are also reported in this response list.

Moreover, when the optimal value N_opt selected is on the boundary of testN, which means that testN most likely was not well chosen, a warning message is printed.

The "QuantifQuantile" class response can be used as argument of the functions plot (only for d ≤ 2), summary and print. The plot function draws the observations and plots the estimated conditional quantile curves (d = 1) or surfaces (d = 2) -for d = 2, the rgl package is used [START_REF] Adler | rgl: 3D visualization device system (OpenGL)[END_REF], which allows to change the perspective in a dynamic way. In order to illustrate the selection of N, the function plot also has an optional argument ise. Setting this argument to TRUE (the default is FALSE), this function, that can be used for any dimension d, provides the plot (against N) of the ISE ᾱ and ISE ¯quantities in (7) or in (8), depending on the choice same_N=FALSE or same_N=TRUE, respectively; see the examples below for details. If d ≤ 2, it also returns the fitted quantile curves or surfaces.

Computing optimal grids

Besides the functions that allow to estimate conditional quantiles and to plot the corresponding reference curves/surfaces, QuantifQuantile provides a function that computes optimal quantization grids. This function, called choice.grid, admits the following arguments :

• X: a d × n real array (required). The columns of this matrix are the X i 's, i = 1, . . . , n, for which the optimal quantization grid should be determined. Each point of X is used as a stimulus in the stochastic gradient algorithm to get an optimal grid.

• N: a positive integer (required). The size of the desired quantization grid.

• ng: a positive integer (optional). The number of desired quantization grids. The default is 1.

• p: a real number larger than or equal to one (optional). This is the parameter p used in the quantization error. The default is 2.

In some cases, it may be necessary to have several quantization grids, such as in ( 5), where B + tildeB grids are needed. The three functions computing quantization-based conditional quantiles then call the function choice.grid with ng > 1. In such case, the various grids are obtained using as stimuli a resampling version of X (the X t 's in the previous section).

The output is a list containing the following elements :

• init_grid: a d×N×ng real array. The entry (i, j, ) of this matrix is the i th component of the j th point of the th initial grid.

• opti_grid: a d×N×ng real array. The entry (i, j, ) of this matrix is the i th component of the j th point of the th optimal grid provided by the algorithm.

Illustrations

In this section, we illustrate on several examples the use of the functions described above. Examples 1-3 restrict to QuantifQuantile/QuantifQuantile.d2 and provide graphical representations in each case.

Example 4 deals with a three-dimensional covariate, without graphical representation. An illustration of the function choice.grid is given in the Appendix.

Example 1 : simulated data with one-dimensional covariate

We generate a random sample (X i , Y i ) , i = 1, . . . , n = 300, where the X i 's are uniformly distributed over the interval (-2, 2) and where the Y i 's are obtained by adding to X 2 i a standard normal error term that is independent of X i :

set.seed(258164) n <-3 X <-runif(n, -2, 2) Y <-X^2 + rnorm(n)
We test the number N of quantizers between 10 and 30 by steps of 5 and we do not change the default values of the other arguments. We then run the function QuantifQuantile with these arguments and stock the response in res. testN <-seq(1 , 3 , by = 5) res <-QuantifQuantile(X, Y, testN = testN)

No warning message is printed, which means that this choice of testN was adequate. To assess this in a graphical way, we use the function plot with ise argument set to TRUE that plots hatISEmean_N (obtained by taking the mean of hatISE_N over alpha) against the various N-values in testN. Figure 1a provides the resulting graph, which confirms that testN was well chosen since hatISEmean_N is larger for smaller and larger values of N than N_opt. We then plot the corresponding estimated conditional quantiles curves in Figure 1b. The default colors of the points and of the curves are changed by using the col.plot argument. This argument is a vector of size 1+length(alpha), whose first component fixes the color of the data points and whose remaining components determine the colors of the various reference curves.

col.plot <-c("grey", "red", "orange", "green", "orange", "red") plot(res, col.plot = col.plot, xlab = "X", ylab = "Y")

It is natural to make the grid testN finer. Of course, the more N-values we test, the longer it takes. This is why we adopted this two-stage approach, where the goal of the first stage was to get a rough approximation of the optimal N-value. In the second stage, we can then refine the grid only in the vicinity of the value N_opt obtained in the first stage.

testN <-c(seq(1 , 2 , by = 1), seq(25, 3 , by = 5)) res_step1 <-QuantifQuantile(X, Y, testN = testN) plot(res_step1, ise = TRUE, col.plot = col.plot, xlab = "X", ylab = "Y")

The resulting graphs are provided in Figures 1c-1d, respectively. We observe that the value of N_opt is made more precise, since we now get N_opt=18 instead of 15. The resulting estimated conditional quantiles curves in Figure 1b are very similar to the ones in Figure 1d. Figure 2: The same results as in Figure 1, but when selecting optimal values of N separately for each α.

So far, we used the default value same_N=TRUE, which leads to selecting an N-value that is common to all α's. For the sake of comparison, we now explore the results for same_N=FALSE. testN <-c(seq(1 , 3 , by = 5)) res2 <-QuantifQuantile(X, Y, testN = testN, same_N = FALSE) plot(res2, ise = TRUE, col.plot = col.plot, xlab = "X", ylab = "Y") testN <-c(seq(1 , 2 , by = 1), seq(25, 3 , by = 5)) res2_step1 <-QuantifQuantile(X, Y, testN = testN, same_N = FALSE) plot(res2_step1, ise = TRUE, col.plot = col.plot, xlab = "X", ylab = "Y")

The results are provided in Figure 2. Comparing the left panels of Figures 1 and2, we see that when choosing N by steps of five, we find N_opt=15 with same_N=TRUE and N_opt= 15 or 20 (depending on alpha) for same_N=FALSE. When we refine the grid testN, we find analogously N_opt=18 for same_N=TRUE and N_opt=14, 15, or 16 for same_N=FALSE. In the present setup, thus, both methods provide relatively close optimal N-values, which explains why the corresponding estimated reference curves are so similar (see the right panels of Figures 1 and2). Therefore, the grid of N-values tested in Figure 1, that may seem too coarse at first sight, actually provides fitted curves that are as satisfactory as those associated with the finer grid in Figure 2.

Example 2 : simulated data with two-dimensional covariate

The sample considered here is made of n = 1, 000 independent realizations of a random vector (X , Y) , where X = (X 1 , X 2 ) is uniformly distributed on the square (-2, 2) 2 and where Y is obtained by adding to X 2 1 + X 2 2 an independent standard normal error term.

The R Journal Vol. XX/YY, AAAA ISSN 2073-4859

set.seed(642516) n <-1 X <-matrix(runif(n*2, -2, 2), ncol = n) Y <-apply(X^2, 2, sum) + rnorm(n)
We test N between 40 and 90 by steps of 10. We change the values of B and tildeB to reduce the computational burden, which is heavier for d = 2 than for d = 1. We keep the default values of all other arguments when running the function QuantifQuantile.d2. Here, a warning message is printed informing us that testN was not well-chosen. We confirm it with the function plot with ise argument set to TRUE.

testN <-seq(4 , 9 , by = 1 ) B <-2 tildeB <-15 res <-QuantifQuantile.d2(X, Y, testN = testN, B = B, tildeB = tildeB) plot(res, ise = TRUE)

Figure 3a provides the resulting graph. The parameter testN was not well chosen since hatISEmean_N becomes smaller and smaller as N_opt increases. We then adapt the choice of testN accordingly and rerun the procedure, which identifies an optimal N-value equal to 100; see Figure 3b.

testN <-seq(8 , 13 , by = 1 ) res <-QuantifQuantile.d2(X, Y, testN = testN, B = B, tildeB = tildeB) plot(res, ise = TRUE)

We then plot the corresponding estimated conditional quantile surfaces in Figure 4. col.plot <-c("black", "red", "orange", "green", "orange", "red") plot(res, col.plot = col.plot, xlab = "X_1", ylab = "X_2", zlab = "Y")

Example 3 : real data study and comparison with some competitors

This example aims at illustrating the proposed estimated reference curves on a real data set and at comparing them with some competitors. In this example, the ncores parameter of QuantifQuantile function was set to the number of cores detected by R minus 1. The data used here, that are included in the QuantifQuantile package, involves several variables related to employment, housing and environment associated with n = 542 towns/villages in Gironde, France. For the present illustration, we restrict to the regressions R 1 and R 2 involving (X, Y) = (percentage of owners living in their primary residence, percentage of buildings area) and (X, Y) = (percentage of middle-range employees, population density), respectively. In both cases, n = 542 observations are available and we are interested in the estimation of reference curves for α = 0.05, 0.25, 0.50, 0.75 and 0.95. For both R 1 and R 2 , we tested the number N of quantizers to be used between 5 and 15 by step of 1, using the methodology described in Example 1.

The R Journal Vol. XX/YY, AAAA ISSN 2073-4859 col.plot <-c("grey", "red", "orange", "green", "orange", "red") plot(res, col.plot = col.plot, xlab = "X", ylab = "Y")

The same exercise is repeated with (X, Y) = (percentage of middle-range employees, population density). For each α-value considered, we obtained Nᾱ ;opt = 13 for R 1 and Nᾱ ;opt = 7 for R 2 . The corresponding quantization-based reference curves are plotted in Figures 5a and5c, respectively. For the sake of comparison, spline-based curves are provided in Figures 5b and5d. These were obtained from the function rqss in the package quantreg. Since the parameter λ involved, that governs the trade-off between fidelity and smoothness, is not automatically selected by rqss, we selected it through AIC (via the AIC function), separately for each order α. rank <-rank(X, ties.method = "random") X[rank] <-X Y[rank] <-Y alpha <-c ( . 5, .25, .5, .75, .95) x <-seq(min(X), max(X), length = 1 ) n <-length(X) lambda <-array( , dim = c(length(alpha), 1)) interval = c( .2, 1 ) for(i in 1:length(alpha)){ for(i in 1:length(alpha)){ lines(fitted_matrix[, i] ~X, type = "l", col = col.plot[i+1]) }

The same exercise is repeated for R 2 , but with λ tested between 0.5 and 15. Since they are piecewise linear, the resulting spline-based reference curves are less smooth than the one based on quantization. Arguably, the latter better adapt to the samples even though they are sometimes quite wiggly.

Of course, the computational burden is also an important issue. Therefore, Table 1 gathers, for each method and each regression problem, the average and standard deviation of the computing times in a collection of 50 runs (these 50 runs were considered to make results more reliable). In each case, our method is faster than its spline-based competitor.

Example 4 : real data with three-dimensional covariate

To treat an example with d > 2, we reconsider the data set in Example 3, this time with the response Y = population density and the three covariates X 1 = percentage of farmers, X 2 = percentage of unemployed workers, and X 3 = percentage of managers. In this setup, no graphical output is available. We therefore restrict to a finite collection of x-values where conditional quantiles are to be esti- mated. Denoting by M j and X j , j = 1, 2, 3, the maximal value and the average of X ij , i = 1, . . . , n = 542, respectively, we consider the following eight values of x :

x 1 =    X 1 X 2 X 3    , x 2 =    1 2 (X 1 + M 1 ) X 2 X 3    , x 3 =    X 1 1 2 (X 2 + M 2 ) X 3    , x 4 =    X 1 X 2 1 2 (X 3 + M 3 )    , x 5 =    1 2 (X 1 + M 1 ) 1 2 (X 2 + M 2 ) X 3    , x 6 =    1 2 (X 1 + M 1 ) X 2 1 2 (X 3 + M 3 )    , x 7 =    X 1 1 2 (X 2 + M 2 ) 1 2 (X 3 + M 3 )    , x 8 =    1 2 (X 1 + M 1 ) 1 2 (X 2 + M 2 ) 1 2 (X 3 + M 3 )    .
The function QuantifQuantile.d is then evaluated for the response and covariates indicated above, and with the arguments alpha=(0.25, 0.5, 0.75) , testN=(5, 6, 7, 8, 9, 10) , x being the 3 × 8 matrix whose columns are the vectors x 1 , x 2 , . . . , x 8 just defined and ncores being the number of cores detected by R minus 1.

data(gironde) set.seed(729848) X1 <-gironde[[1]]$farmers X2 <-gironde[[1]]$unemployed X3 <-gironde[[1]]$managers Y <-gironde[[2]
]$density X <-matrix(c(X1, X2, X3), nr = 3, byrow = TRUE) n <-length(X)/3 d <-3 alpha <-c( .25, .5, .75) x1 <-round(c(mean(X1), mean(X2), mean(X3))) x2 <-round(c((mean(X1) + max(X1))/2, mean(X2), mean(X3))) x3 <-round(c(mean(X1), (mean(X2) + max(X2))/2, mean(X3))) x4 <-round(c(mean(X1), mean(X2), (mean(X3) + max(X3))/2)) x5 <-round(c((mean(X1) + max(X1))/2, (mean(X2) + max(X2))/2, mean(X3))) x6 <-round(c((mean(X1) + max(X1))/2, mean(X2), (mean(X3) + max(X3))/2)) x7 <-round(c(mean(X1), (mean(X2) + max(X2))/2, (mean(X3) + max(X3))/2)) x8 <-round(c((mean(X1) + max(X1))/2, (mean(X2) + max(X2))/2, (mean(X3) + max(X3))/2))

x <-matrix(c(x1, x2, x3, x4, x5, x6, x7, x8), nr = d) res <-QuantifQuantile.d(X, Y, x , alpha = alpha, testN = seq(5, 1 , by = 1), same_N = F, ncores = detectCores() -1) round(fitted.values(res), 2)

This provided Nᾱ

;opt = 8, 7 and 7, for α = 0.25, 0.50 and 0.75, respectively. The total computation time is 6.86 seconds. The fitted.values function then allowed to return the following matrix hatq_opt of estimated conditional quantiles : [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [1,] 44.3 22.59 39.5 71.59 25. 5 22.4 76.37 24.19 [2,] 8 . 7 32.31 81.24 161.85 35. 1 31.92 145.29 38.18 [3,] 139.16 46.5 223.13 344.92 53.73 47. 1 4 2.98 73.19 This collection of (estimated) conditional quartiles allows to appreciate the impact of a marginal perturbation of the covariates on Y's conditional median (location) or interquartile range (scale). For instance, the results suggest that Y's conditional median decreases with X 1 , is stable with X 2 , and increases with X 3 , whereas its conditional interquartile range decreases with X 1 but increases much with X 2 and with X 3 . The eight x-values considered further allow to look at the joint impact of two or three covariates on Y's conditional location and scale. Of course, other shifts in the covariates (and other orders α) should further be considered to fully appreciate the dependence of Y on X.

Conclusion

In this paper, we described the package QuantifQuantile that allows to implement the quantizationbased quantile regression method introduced in Charlier et al. (2015a,b). The package is simple to use, as the function QuantifQuantile and its multivariate versions essentially only require providing the covariate and response as arguments. Since the choice of the tuning parameter N is crucial, a warning message is printed if it is not well-chosen and the function plot can also be used as guide to change adequately the value of the parameter testN in the various functions. Moreover, a graphical illustration is directly provided by the same function plot when the dimension of the covariate is smaller than or equal to 2. Finally, this package also contains a function that provides optimal quantization grids, which might be useful in other contexts, too.

Finally, we stress that quantization-based estimators, like most nonparametric smoothing procedures, are likely to perform poorly in high-dimensional situations due to the curse of dimensionality. For large d, it is therefore unclear how to assess whether a given covariate has a significant impact on the response variable. For small d, however, it is always possible, in the absence of a formal testing procedure, to resort to visual inspection. In the simplest case of a single covariate (d = 1), this would lead to looking whether or not fitted curves approximately are horizontal lines. This can be extended to the case d = 2. Figure 6: For n = 500 (left) and n = 5, 000 (right), a random sample of size n from the uniform distribution over (-2, 2) (in grey), an initial grid of size 15 obtained by sampling without replacement among these n observations (in red), and the "optimal" grid returned by choice.grid (in green). The R Journal Vol. XX/YY, AAAA ISSN 2073-4859

Finally, we turn to the bivariate case and generate two random samples of size n = 2, 000 and size n = 20, 000 from the uniform distribution over the square (-2, 2) 2 . The function choice.grid was applied to these samples with N=30 and ng=1. The resulting couple of initial and optimal grids are plotted in Figure 7. As in the univariate case, we observe an improvement when going from the initial grids to the corresponding optimal grids provided by the function choice.grid (here as well, the population optimal grid should be uniformly spread over the support of the underlying distribution). Also, it is still the case that the resulting optimal grid is better when based on a larger sample size n.

set.seed(345689) n <-2 X <-matrix(runif(n*2, -2, 2), nc = n) N <-3 ng <-1 res <-choice.grid(X, N, ng) col <-c("red", "forestgreen") plot(res$init_grid[1,,1], res$init_grid[2,,1], col = col[1], xlab = "", ylab = "") plot(res$opti_grid[1,,1], res$opti_grid[2,,1], col = col[2], xlab = "", ylab = "")

  > args(QuantifQuantile) function(X, Y, alpha = c( . 5, .25, .5, .75, .95), x = seq(min(X), max(X), length = 1 ), testN = c(35, 4 , 45, 5 , 55), p = 2, B = 5 , tildeB = 2 , same_N = TRUE, ncores = 1) > args(QuantifQuantile.d2) function(X, Y, alpha = c( . 5, .25, .5, .75, .95), x = matrix(c(rep(seq(min(X[1, ]), max(X[1, ]), length = 2 ), 2 ), sort(rep(seq(min(X[2, ]), max(X[2, ]), length = 2 ), 2 ))), nrow = 2, byrow = TRUE), testN = c(11 , 12 , 13 , 14 , 15 ), p = 2, B = 5 , tildeB = 2 , same_N = TRUE, ncores = 1) > args(QuantifQuantile.d) The R Journal Vol. XX/YY, AAAA ISSN 2073-4859 function(X, Y, x, alpha = c( . 5, .25, .5, .75, .95), testN = c(35, 4 , 45, 5 , 55), p = 2, B = 5 , tildeB = 2 , same_N = TRUE, ncores = 1)

Figure 1 :

 1 Figure 1: For the sample considered in Example 1, this figure provides (a) the plot of N → ISE ¯(N) with N ∈ {10, 15, 20, 25, 30}, and (b) the resulting reference curves. The panels (c)-(d) provide the corresponding plots when taking N ∈ {10, 11, 12, . . . , 19, 20, 25, 30}.

Figure 3 :

 3 Figure 3: For the sample considered in Example 2, this figure plots N → ISE ¯(N) (a) for N ∈ {40, 50, 60, 70, 80} and (b) for N ∈ {80, 90, . . . , 120, 130}.

Figure 4 :

 4 Figure 4: For the sample considered in Example 2, this figure plots (with two different views) the estimated conditional quantile surfaces obtained with the plot function for α =0.05, 0.25, 0.50, 0.75 and 0.95.

Figure 5 :

 5 Figure 5: Estimated conditional quantile curves obtained with QuantifQuantile (left) and rqss (right), for regression R 1 (top) and regression R 2 (bottom). In each case, the quantile orders considered are α = 0.05, 0.25, 0.50, 0.75 and 0.95.

Figure 7 :

 7 Figure7: For n = 2, 000 (left) and n = 20, 000 (right), an initial grid of size 15 obtained by sampling without replacement among a random sample of size n from the uniform distribution over (-2, 2) (top), and the corresponding optimal grid returned by choice.grid (bottom).

Table 1 :

 1 Averages of the computing times (in seconds) to obtain 50 times the conditional quantile curves in Example 3 for QuantifQuantile and rqss, respectively; standard deviations are reported in parentheses.

	The R Journal Vol. XX/YY, AAAA	ISSN 2073-4859

Our method actually requires assumptions on the link function between Y and X and on the error term of the model. These assumptions in particular guarantee that, for any x, the conditional distribution of Y given X = x is absolutely continuous; we refer toCharlier et al. (2015a) for details.

Since the number J of points in a default value of x obtained in this fashion would increase exponentially with the dimension d, we did not adopt the same approach for d ≥

Acknowledgments

The authors are grateful to the Editor and two anonymous referees for their careful reading and insightful comments that led to a significant improvement of the original manuscript. The first author's research is supported by a Bourse FRIA of the Fonds National de la Recherche Scientifique, Communauté française de Belgique. The second author's research is supported by an A.R.C. contract from the Communauté Française de Belgique and by the IAP research network grant P7/06 of the Belgian government (Belgian Science Policy).

Appendix

Illustration of choice.grid

We here put to work the function choice.grid in the univariate and bivariate cases. This function provides the "optimal" grid generated by the stochastic gradient algorithm described earlier. As above mentioned, quantization was extensively used in many other fields as numerical integration, cluster analysis, numerical probability or finance [START_REF] Pagès | A space quantization method for numerical integration[END_REF][START_REF] Pagès | Optimal quantization methods and applications to numerical problems in finance[END_REF]. Therefore, this function can be of interest outside the regression setup considered here.

We start with the univariate case and generate a random sample of size n = 500 from the uniform distribution over (-2, 2). With N=15 and ng=1, this function provides a single initial grid (obtained by sampling without replacement among the uniform sample) and the corresponding optimal grid returned by the algorithm. Figure 6 represents the observations (in grey), the initial grid (in red), and the optimal grid (in green). The same exercise is repeated with sample size n = 5, 000, and the results are also given in Figure 6.

set.seed(643625) n <-5 X <-runif(n, -2, 2) N <-15 ng <-1 res <-choice.grid(X, N, ng) # Plots of the initial and optimal grids plot(X, rep(1, n), col = "grey", cex = .5, ylim = c(-.1, 1.1), yaxt = "n", ylab = "") points(res$init_grid, rep( .5, N), col = "red", pch = 16, cex = 1.2) points(res$opti_grid, rep( , N), col = "forestgreen", pch = 16, cex = 1.2)

Since the parent distribution is uniform over (-2, 2), the population optimal grid is the equispaced grid on that interval [START_REF] Pagès | A space quantization method for numerical integration[END_REF]. For both sample sizes considered, the optimal grid provided by the choice.grid function is much closer to the population optimal grid than the initial one. Recalling that the stochastic gradient algorithm in choice.grid performs as many iterations as observations in the original sample, it is not surprising that the optimal grid associated with the sample of size 5, 000 better approximates the population optimal grid than the optimal grid associated with the sample of size 500.

The R Journal Vol. XX/YY, AAAA ISSN 2073-4859