domain: scco shs info

text:
Attigodu C Ganesh

Frédéric Berthommier

Jean-Luc Schwartz
email: jean-luc.schwartz@gipsa-lab.grenoble-inp.fr

Audiovisual speech scene analysis in the context of competing sources

Keywords: audiovisual fusion, McGurk effect, scene analysis, attention, seniors

Audiovisual fusion in speech perception is generally conceived as a process independent from scene analysis, which is supposed to occur separately in the auditory and visual domain. On the contrary, we have been proposing in the last years that scene analysis such as what takes place in the cocktail party effect was an audiovisual process. We review here a series of experiments illustrating how audiovisual speech scene analysis occurs in the context of competing sources. Indeed, we show that a short contextual audiovisual stimulus made of competing auditory and visual sources modifies the perception of a following McGurk target. We interpret this in terms of binding, unbinding and rebinding processes, and we show how these processes depend on audiovisual correlations in time, attentional processes and differences between junior and senior participants.

Introduction

The classical cocktail party effect and the problems it raises for speech perception in adverse conditions [START_REF] Cherry | Some Experiments on the Recognition of Speech, with One and with Two Ears[END_REF] has generated two series of theoretical and experimental developments, which remain surprisingly separate in the years. On the one hand, Auditory Scene Analysis (ASA) puts at its agenda the search for auditory mechanisms enabling to group together auditory cues into auditory primitives, based on principles such as temporal synchrony, correlations in time or space, and more globally common fate [START_REF] Bregman | Auditory scene analysis[END_REF]. On the other hand, Audiovisual Speech Perception (AVSP) is focused on the cognitive processes enabling to fuse together the auditory and the visual input corresponding to the speaker's utterances (e.g. [START_REF] Massaro | Speech Perception by Ear and Eye[END_REF][START_REF] Summerfield | Some preliminaries to a comprehensive account of audiovisual speech perception[END_REF][START_REF] Schwartz | Ten years after Summerfield. A taxonomy of models for audiovisual fusion in speech perception[END_REF]), while assuming implicitly or explicitly that scene analysis occurs independently in the auditory and the visual domain.

As a matter of fact, some experimental studies display evidence where unimodal perceptual grouping precedes multisensory integration (e.g. [START_REF] Sanabria | Intramodal perceptual grouping modulates multisensory integration: evidence from the crossmodal dynamic capture task[END_REF][START_REF] Keetels | Auditory grouping occurs prior to intersensory pairing: evidence from temporal ventriloquism[END_REF][START_REF] Munhall | Integration of Speech in a Bistable Illusion[END_REF]). However, several recent behavioral and neurophysiological studies have suggested that the presentation of a visual stream can affect primary auditory streaming by enhancing segregation or integration ( [START_REF] Rahne | Visual cues can modulate integration and segregation of objects in auditory scene analysis[END_REF][START_REF] Marozeau | The Effect of Visual Cues on Auditory Stream Segregation in Musicians and Non-Musicians[END_REF][START_REF] Berthommier | Audiovisual streaming in voicing perception: new evidence for a low-level interaction between audio and visual modalities[END_REF][START_REF] Devergie | The effect of lip-reading on primary stream segregation[END_REF][START_REF] Maddox | Auditory selective attention is enhanced by a task-irrelevant temporally coherent visual stimulus in human listeners[END_REF]). This suggests that audiovisual speech perception likely incorporates a stage of audiovisual scene analysis before fusion might operate [START_REF] Andersen | The role of visual spatial attention in audiovisual speech perception[END_REF][START_REF] Alsius | Searching for audiovisual correspondence in multiple speaker scenarios[END_REF]. This led our group propose since a number of years that the ASA and AVSP frameworks should be combined within a single framework, that we tentatively called "audiovisual speech scene analysis" (AVSSA). In this framework, it is assumed that audiovisual speech perception is a two-stage process (as is auditory perception in the ASA framework), beginning by a scene analysis process where auditory and visual cues are bound within "audiovisual primitives", before audiovisual fusion for decision at a later stage. Fusion would hence operate on a set of audiovisual cues -an audiovisual source -selected at the first stage, and the result of the fusion process would depend on the coherence of the auditory and visual inputs at the first stage [START_REF] Berthommier | A phonetically neutral model of the low-level audio-visual interaction[END_REF].

The two-stage model of audiovisual speech perception in the framework of AVSSA has been elaborated in the last years, taking advantage of a specific paradigm that we developed to demonstrate that audiovisual fusion does indeed depend on the coherence of the unisensory inputs. This paradigm is based on the classical "McGurk effect", in which a visual input incongruent with an auditory input may change the perception of the sound: typically a visual "ga" dubbed on an auditory "ba" leads to the perception of "da" or "tha" [START_REF] Mcgurk | Hearing lips and seeing voices[END_REF]. In two series of experiments [START_REF] Nahorna | Binding and unbinding the auditory and visual streams in the McGurk effect[END_REF][START_REF] Nahorna | Audiovisual speech scene analysis: characterization of the dynamics of unbinding and rebinding the McGurk effect[END_REF], we were able to show that if such a McGurk stimulus is preceded by a period of audiovisual material displaying incoherent variations, the participants are lead to consider that the auditory and visual inputs should not be bound together, hence the McGurk effect decreases. This is what we called "unbinding", which can be possibly followed by a period of audiovisual coherence leading to "rebinding".

The binding/unbinding/rebinding processes should in our view play a crucial role in audiovisual speech perception in adverse but ecological conditions where various speakers discuss altogether: typically the cocktail party situation, where the listener should be able to attend to a source and adequately select the corresponding pieces of audio and video information before being able to fuse and understand the message. This is why we recently attempted to evaluate the binding/unbinding/rebinding paradigm with mixtures of sources, to assess whether the AVSSA framework would be able to account for such kinds of configurations. In this paper we will report the major results of this series of experiments about audiovisual speech scene analysis in the context of competing sources.

In Section 2 we present the binding/unbinding/rebinding paradigm together with the previously obtained results [START_REF] Nahorna | Binding and unbinding the auditory and visual streams in the McGurk effect[END_REF][START_REF] Nahorna | Audiovisual speech scene analysis: characterization of the dynamics of unbinding and rebinding the McGurk effect[END_REF]. Then, Section 3 reports the results of two series of experiments dealing with AVSSA with a mixture of audiovisual sources. Section 4 will be devoted to results obtained on a population of older subjects, showing how the binding abilities vary with age.

Binding, unbinding, rebinding

The experimental paradigm developed previously [START_REF] Nahorna | Binding and unbinding the auditory and visual streams in the McGurk effect[END_REF][START_REF] Nahorna | Audiovisual speech scene analysis: characterization of the dynamics of unbinding and rebinding the McGurk effect[END_REF] and which will be adapted in some of the present experiments is displayed in Fig. 1. The principle is the following. A pure McGurk "target" made of an auditory "ba" and a visual "ga" is not 100% perceived as "ba" because of the McGurk effect leading a number of participants to perceive it as "da": this is "binding" (Fig. 1, left). The assumption is that listeners have a "default binding stage", in which they consider that in lack of any further evidence, the auditory and visual inputs should be bound together for audiovisual fusion.

If the target is preceded by a given duration of incoherent audiovisual material (context), the amount of fusion decreases, hence the score of "ba" responses increases: this is unbinding (Fig. 1,middle). The material presented in all the previous experiments [START_REF] Nahorna | Binding and unbinding the auditory and visual streams in the McGurk effect[END_REF][START_REF] Nahorna | Audiovisual speech scene analysis: characterization of the dynamics of unbinding and rebinding the McGurk effect[END_REF] and in the present paper were prepared from audiovisual material produced by a French male speaker, JLS, with lips painted in blue to allow precise video analysis of lip movements [START_REF] Lallouache | Un poste 'visage-parole.' Acquisition et traitement de contours labiaux (A 'face-speech' workstation. Acquisition and processing of labial contours)[END_REF]. The videos consisted of the entire speaker's face, keeping natural colors apart from the blue make-up. Recordings were digitized at an acoustic sampling frequency of 44.1 kHz and a video sampling frequency of 50 Hz (25 images per second with two frames per image).

The target was either a congruent audiovisual "ba" syllable, or an incongruent McGurk stimulus with an audio "ba" dubbed on a video "ga" with precise temporal synchronization at the plosive burst. The focus was actually on McGurk targets and the congruent "ba" targets were only presented as controls.

The context and reset stimuli in all the experiments were prepared from various kinds of mixtures (that will be described later) of two types of audiovisual material. The first type, called "syllables", was made of random sequences of audiovisual syllables in the set: "pa", "ta", "va", "fa", "za", "sa", "ka", "ra", "la", "ja", "cha", "ma" or "na". The syllable rhythm was about 1.5 Hz, hence for sequences of 2 or 4 syllables that were used in the experiments presented here, the context duration varied between 1.3 and 2.7 s depending on the number of uttered syllables. The second type, called "sentences", consisted of excerpts from sentences invented online by the speaker at the recording stage.

Unbinding and rebinding with competing sources

The princeps experiments [START_REF] Nahorna | Binding and unbinding the auditory and visual streams in the McGurk effect[END_REF][START_REF] Nahorna | Audiovisual speech scene analysis: characterization of the dynamics of unbinding and rebinding the McGurk effect[END_REF] involved either coherent audiovisual contexts made of the coherent auditory and visual content of syllable sequences (A and V syllables in the following), or incoherent contexts in which A syllables were mixed with V sentences with the adequate duration.

Of course, the incoherent context material actually corresponded to two differing sources, one syllabic source presented in the auditory modality and one sentence source presented in the visual modality. However, there was no competition of sources in individual modalities.

The objective in present experiments was to go towards more realistic situations involving a competition of sources inside the auditory modality -apart from possible incoherence between modalities.

Adding noise in the audiovisual context (Exp. 1)

In a first step (Exp. 1), we considered acoustic noise added to the acoustic source in the context part. Therefore, we exploited the "unbinding/rebinding" paradigm presented in Fig. 1, though with an important variant, that is the addition of acoustic noise in the context before the McGurk target [START_REF] Ganesh | Effect of context, re-binding and noise, on audiovisual speech fusion[END_REF].

The experiment involved two types of context before the two targets (see Fig. 2). Firstly, a coherent context made of coherent A and V syllables (Fig. 2, bottom), provided a baseline. Secondly, an incoherent context made of 2 or 4 A syllables dubbed on excerpts from V sentences with the adequate duration was followed by a coherent reset made of a variable number (0, 1, 2 or 3) of coherent AV syllables (Fig. 2, top). This should produce unbinding (by the incoherent context) and a given amount of rebinding (depending on the number of syllables in the reset). The (context + reset) portion was presented either in clear or in acoustic noise (Gaussian white noise at 0 dB SNR), with no noise on the target. Thirty-one participants (22 women and 9 men; 30 righthanded and 1 left-handed; mean age=31.7 years; SD=11.7 years) took part in this experiment. The task consisted in monitoring online for the perception of either "ba" or "da" syllables -reminding the reader that there were in fact no "da" syllables all along the experiment. The participants had to signal the perception of either "ba" or "da" by pressing as soon as possible an adequate response button. Responses were only considered within a given amount of time after the target onset defined as the burst onset (between 200 ms and 1200 ms, see [START_REF] Nahorna | Binding and unbinding the auditory and visual streams in the McGurk effect[END_REF] for more explanations on the experimental procedure). "

" "

!

The results are displayed on Fig. 3 (averaged over the two context durations, 2 and 4 syllables). Fig. 3a displays the results obtained with no noise in the context period. They basically replicate the results in the original study [START_REF] Nahorna | Audiovisual speech scene analysis: characterization of the dynamics of unbinding and rebinding the McGurk effect[END_REF], displayed in a schematic way on Fig. 1: while the amount of McGurk responses is large with a coherent context (context "1" in the figure, with less than 30% "ba" responses, hence more than 70% McGurk responses), this amount decreases with an incoherent context (context "2", with more than 50% "ba" responses) and comes back to its original value with a coherent reset made of 3 coherent audiovisual syllables (context "3"). This is the unbinding-rebinding phenomenon described earlier. The interpretation we suggest is that noise in the context period is used by the subjects as an indication that the auditory sensor is less reliable and should hence be considered with a smaller weight in the fusion process. This is in line with various proposals about a "weighted fusion process", according to which fusion would depend globally on the weight each subject would attribute to the auditory and visual channels, possibly varying with subject, language, noise and a number of other factors (e.g. [START_REF] Heckmann | Noise adaptive stream weighting in audio-visual speech recognition[END_REF][START_REF] Schwartz | A reanalysis of McGurk data suggests that audiovisual fusion in speech perception is subject-dependent[END_REF][START_REF] Huyse | Degradation of labial information modifies audiovisual speech perception in cochlear-implanted children[END_REF]).

Adding a competing source in the audiovisual context (Exp. 2)

In a second step, we aimed to further explore the possibility that a scene analysis process would take place in the course of AV fusion. For this aim, we presented contexts made of a mixture of sources. More precisely, we prepared mixtures of two audio sources associated with one video input, being coherent with the one of the two audio sources [START_REF] Ganesh | Audio Visual integration with competing sources in the framework of Audio Visual Speech Scene Analysis[END_REF].

We predicted that the binding stage in AVSSA would now serve two roles: 1) compute partial correlations, which could enable the system to select the audio source coherent with the video input, and 2) assess the binding state modulating AV fusion. To test these predictions, we mixed two audio sources which have very different properties over time, and which are hence likely to lead to very different correlations with their corresponding video counterpart: a syllable stream and a sentence stream (see Fig. 4). Indeed, syllables correspond to stronger AV modulations in time and hence stronger AV coherence than sentences. Therefore, the association between the visual input and the corresponding auditory input should be stronger for syllables than for sentences. Hence, we predicted that the coherence of the AV context would be stronger for visual syllables, and would lead to a larger visual weight and more McGurk effect than with visual sentences. Experiment 2 actually consisted of two parts. Firstly (Exp. 2a) the subjects performed the same task as before -that is monitored for the presence of "ba" or "da" targets" -with one or the other context, to assess whether the "video syllables" context would produce more "da" responses than the "video sentences" context. Secondly (Exp. 2b), the same subjects performed the same task though with an additional demand: to attend selectively either to the syllables or to the sentences (with a counterbalanced order of the "attention to syllables" and "attention to sentences" conditions). The experiment was passed by twenty-nine French participants without hearing or vision problems (22 women and 7 men; 27 right-handed and 2 left-handed; mean age= 29.2 years; SD=10. 4 years). In the first part with no specific attentional demand (Exp. 2a), the results were actually in line with our prediction, with a 10% higher amount of McGurk responses (10% lower amount of "ba" scores) with video syllables. This is due in our view to a larger audiovisual correlation with syllables, hence more binding and more McGurk responses. The results obtained in the second part with an attentional focus on syllables or sentences (Exp. 2b) are displayed on Fig. 5. For the attentional focus on syllables ("Attention syllables") we replicate the results with no attention: there is unbinding from the video syllables (V. syll.) to the video sentences (V. sent.) context, hence more "ba" responses in the second case. The fact that attention is put on syllables does not modify the " ! ! !

results with no attentional focus, possibly because video syllables are so correlated that they "pop out" automatically as coherent in the binding process. However, attention put on sentences ("Attention sentences") seems to "rebind" the auditory and visual streams for the visual sentences context (V. sent.), hence removing the difference between the two contexts "V. syll." and "V. sent." This could be due to topdown processes based on audiovisual schemas, according to which the intrinsically low binding associated with audiovisual sentences would be enhanced thanks to the attentional process.

Assessing binding in seniors

Our last set of experiments consisted in assessing whether binding could vary with age. A series of recent papers have attempted to establish whether audiovisual integration was different in senior (typically above 65 years of age) contrasted to junior adults. The conclusions are not fully convergent. Altogether, while seniors display a degradation in their auditory perception but less so in their visual perception of speech, their integration abilities seem to stay stable, with actually a trend for a higher role of the visual input in speech perception in seniors (see a review in [START_REF] Mozolic | Multisensory Integration and Aging[END_REF], and recent evidence in [START_REF] Sekiyama | Enhanced audiovisual integration with aging in speech perception: a heightened McGurk effect in older adults[END_REF] that the visual influence is greater in older adults compared with younger ones not only with equal SNRs but also with SNRs calibrated to equalize unisensory performance).

The question we decided to explore in this final study was hence if binding would operate in a different way in seniors compared to juniors [START_REF] Ganesh | Audiovisual Binding in Elderly Population[END_REF]. Indeed, since the visual modality is of particular importance for seniors, it could be envisioned that the visual input would be exploited and fused even in case of incoherence, hence a decrease in unbinding in older adults. However, since the audiovisual scenes used in our experimental paradigm are rather complex, it could be suggested quite in the contrary that unbinding would be larger with a larger modulation by context and attention.

To test these two inverse predictions, we exploited two of the previous experimental paradigms in this study that are Experiment 1 without noise (the primary "unbindingrebinding" paradigm without noise developed in [START_REF] Nahorna | Audiovisual speech scene analysis: characterization of the dynamics of unbinding and rebinding the McGurk effect[END_REF]) and Experiment 2b (unbinding with context including two competing audio streams, and with attentional focus either on syllables or on sentences). The results are displayed on Fig. 6. They compare the results obtained with our senior group with those obtained with the younger group reported previously. They display two very interesting results. Firstly, in Experiment 1, while the average McGurk score was typically the same in our two groups of subjects with a coherent context (no significant difference in the percentage of "ba" responses), there was a strong and significant increase in unbinding in seniors (more than 15% increase in "ba" responses from juniors to seniors with the incoherent context). Rebinding resulted in recovering the same amount of fusion in the two populations. Secondly, in Experiment 2b, the results were qualitatively similar to those obtained with juniors, with effects of both context and attention as described in Section 3.2, without obvious differences in the size of the effects. The lessons of this set of experiments on seniors are rather clear. Firstly, globally, they provide a replication of the results with juniors, confirming (i) that context matters in fusion; (ii) that unbinding-rebinding processes modulate the McGurk process; (iii) and that within competing sources syllables produce more binding but attention may rebind in the case of audiovisual sentences.

Secondly, they suggest that seniors unbind more than juniors. The larger effect of unbinding in older subjects could be related to the fact that under cognitive load, integration reduces (see [START_REF] Alsius | Audiovisual integration of speech falters under high attention demands[END_REF][START_REF] Alsius | Attention to touch weakens audiovisual speech integration[END_REF]). Indeed, it could be assumed that in the case of incoherence, a certain amount of attention is required for keeping audition and vision bound together and hence produce binding. If the ability to maintain this amount of attention is decreased in seniors, this would result in less fusion and more unbinding, which is actually what happens here.

Conclusions

All the data presented in the present paper confirm that audiovisual context may modulate audiovisual fusion displayed by the McGurk effect. They are all in good agreement with the two-stage "binding-and-fusion" model that we are developing in our team, in the framework of Audio-Visual Speech Scene Analysis (AVSSA). They converge towards a model in which the binding stage has actually various roles, such as audiovisual estimation of the reliability of the sources (Exp. 1), and selection of the adequate unisensory sources or pieces of sources for construction of the audiovisual source to process (Exp. 2; see also [START_REF] Berthommier | Audiovisual streaming in voicing perception: new evidence for a low-level interaction between audio and visual modalities[END_REF]). Interestingly, senior data also suggest that binding could be a computationally costly process, which may suffer some limitations in case of increased cognitive load (see also [START_REF] Alsius | Audiovisual integration of speech falters under high attention demands[END_REF][START_REF] Alsius | Attention to touch weakens audiovisual speech integration[END_REF]). The binding paradigm exploited in the present study once again appears as a very efficient tool for studying the cognitive processes likely to take place in Audio-Visual Speech Scene Analysis.

Acknowledgements

This project has been supported by Academic Research Community "Quality of life and ageing" (ARC 2) of the Rhône-Alpes Region, which provided a doctoral funding for Ganesh Attigodu Chandrashekara

  If a "reset" stimulus made of audiovisual coherent material is presented after the incoherent context and before the McGurk target, subjects recover the original McGurk effect, hence the percentage of "ba" responses decreases down to its level for a pure McGurk target: this is rebinding (Fig 1, right).

Figure 1 :

 1 Figure 1: The unbinding/rebinding paradigm

Figure 2 :

 2 Figure 2: Material for Experiment 1

Fig. 3b displays

  Fig. 3b displays how results change when there is noise on the context -though NOT in the target. The amount of McGurk effect largely increases for all contexts -hence the percentage of "ba" responses decreases by 13% to 30% depending on context. There remains unbinding and rebinding, since the scores still vary from the coherent to the incoherent context and back to the incoherent + coherent reset context. However, the fluctuations are largely decreased in noise.

Figure 3 :

 3 Figure 3: Results for Experiment 1

Figure 4 :

 4 Figure 4: Material for Experiment 2

Figure 5 :

 5 Figure 5: Results for Experiment 2b

  Twenty-five native French speaking older adults participated in the experiments (2 women and 23 men; from 60 to 75 years, 21 right-handed and 4 left-handed, mean age= 65.3 years; SD=3.9 years). None of them reported any hearing, vision (after correction) or neurological disorders.

Figure 6 :

 6 Figure 6: Results for Experiments 1 and 2 on seniors

  . The research leading to these results has received funding from the European Research Council under the European Community's Seventh Framework Programme (FP7/2007-2013 Grant Agreement no. 339152, "Speech Unit(e)s", J.-L. Schwartz PI).


Julyan Arbel
email: julyan.arbel@inria.fr

Jean-Bernard Salomond
email: jean-bernard.salomond@u-pec.fr

Sequential Quasi Monte Carlo for Dirichlet Process Mixture Models

Keywords: Bayesian Nonparametrics, Density estimation, Quasi random variables, Monte Carlo methods, Sequential sampling

In mixture models, latent variables known as allocation variables play an essential role by indicating, at each iteration, to which component of the mixture observations are linked. In sequential algorithms, these latent variables take on the interpretation of particles. We investigate the use of quasi Monte Carlo within sequential Monte Carlo methods (a technique known as sequential quasi Monte Carlo) in nonparametric mixtures for density estimation. We compare them to sequential and non sequential Monte Carlo algorithms. We highlight a critical distinction of the allocation variables exploration of the latent space under each of the three sampling approaches.

Introduction

Sequential quasi Monte Carlo (SQMC) sampling is a novel sampling scheme proposed by [START_REF] Gerber | Sequential quasi Monte Carlo[END_REF] in a paper read before The Royal Statistical Society. It bridges ideas from quasi Monte Carlo (QMC) methods and from sequential Monte Carlo (SMC) sampling. The latter, also called particle filtering, focuses on sampling in state space models based on sequential data, while the former's initial goal is to enhance Monte Carlo efficiency by using evenly spread vectors instead of (unconstrained) random vectors. As Gerber & Chopin put it, we can only think that the full potential of QMC in statistics remains underexplored.

In this note, we explore applications of SQMC to the field of Bayesian nonparametrics. More specifically, we focus on nonparametric mixtures for density estimation. The method of nonparametric mixtures of kernels appeared in [START_REF] Ferguson | A Bayesian analysis of some nonparametric problems[END_REF] and [START_REF] Lo | On a class of Bayesian nonparametric estimates: I. Density estimates[END_REF] with Dirichlet process mixtures (DPM). Posterior sampling for DPM is usually carried out by Markov chain Monte Carlo (MCMC, see [START_REF] Jain | A split-merge Markov chain Monte Carlo procedure for the Dirichlet process mixture model[END_REF]. See also [START_REF] Blei | Variational inference for Dirichlet process mixtures[END_REF] for a variational Bayes approach. A known drawback of MCMC is its difficulty to cope with multimodality, ie the Markov chain can be trapped in local modes. Another limitation of MCMC is that it is not suited for online data acquisition, that is when data points are observed one at a time (see [START_REF] Caron | Bayesian Inference for Linear Dynamic Models With Dirichlet Process Mixtures[END_REF]. Sequential Monte Carlo methods have proved successful in addressing these issues, in general and also in the particular context of nonparametric mixture models. These models were first cast as SMC samplers by [START_REF] Liu | Nonparametric hierarchical Bayes via sequential imputation[END_REF][START_REF] Maceachern | Sequential importance sampling for nonparametric Bayes models: The next generation[END_REF] as follows. Observations are spread out into unobserved clusters. Hence the cluster labels, or allocation variables, are latent variables which act as the states of the observations in the context of SMC. Two features of the states in this setting are unusual in a QMC setting and will be discussed in this note: they are discrete and non Markovian. The original algorithm was later extended to particle filters by [START_REF] Fearnhead | Particle filters for mixture models with an unknown number of components[END_REF], and further variations and applications have been proposed by Ülker et al. (2010); [START_REF] Carvalho | Particle learning for general mixtures[END_REF]; [START_REF] Griffin | Sequential Monte Carlo methods for mixtures with normalized random measures with independent increments priors[END_REF]; [START_REF] Tsiligkaridis | Adaptive Low-Complexity Sequential Inference for Dirichlet Process Mixture Models[END_REF], among others.

Our goal here is to extend these existing SMC samplers for DPM to the SQMC paradigm. How does SQMC fare in this nonparametric Bayes framework? The extension that we pursue is not trivial and challenging due to the two aforementioned peculiarities: (i) the state-space is discrete, unlike the examples developed by [START_REF] Gerber | Sequential quasi Monte Carlo[END_REF] with continuous state-spaces; (ii) the allocation variables (the states) are not Markovian. The challenges can be formulated as follows: are QMC methods still advantageous in discrete state-space settings? Is the non-Markovianity of the latent variables an impediment to practical implementation, or further to computational gains? For our investigation, we confine attention to Dirichlet process mixtures. Nevertheless, the ideas developed here reach far beyond this model.

2 Background on sampling: SMC, QMC and SQMC Sequential inference refers to the ability to update the estimation as new observations arrive. Sequential Monte Carlo (SMC), or Particle filtering, is a principled technique which sequentially approximates the full posterior using particles [START_REF] Doucet | Sequential Monte Carlo methods in Practice[END_REF][START_REF] Del Moral | Sequential Monte Carlo samplers[END_REF][START_REF] Andrieu | Particle Markov chain Monte Carlo methods[END_REF]. It focuses on sequential state-space models: the density of the observations y t conditionally on Markov states x t in X ⊆ R d is given by y

t |x t ∼ f Y (y t |x t ), with kernel x 0 ∼ f X 0 (x 0 ), x t |x t-1 ∼ f X (x t |x t-1
).

(1) The output of an SMC sampler with N particles on (1 : T ) time period is a collection of weighted particles W

(n) T , x (n) 1:T for n = 1, . . . , N , W (n) T > 0, N n=1 W (n) T
= 1 which approximate the posterior distribution of the states, p(x 1:T |y 1:T ), by a discrete distribution whose support points are the particles. Thus the posterior expectation of any function g of the states can be approximated by

E[g(x 1:T |y 1:T )] ≈ N n=1 W (n) T g x (n) 1:T
with an error rate of order O P (N -1/2 ).

The initial motivation of quasi Monte Carlo (QMC) is to use low discrepancy vectors instead of unconstrained random vectors in order to improve the calculation of integrals via Monte Carlo. The discrepancy D(u) of a sample of vectors u = u (1:N ) is more precisely defined as the total variation distance (TV)

D(u) = d T V (P u , λ)
between the empirical measure of the sample u,

P u (•) = 1 N N n=1 1(u (n) ∈ •)
and the Lebesgue measure λ. Recall that the TV distance between two probability measures P and Q is defined by

d T V (P, Q) = sup A∈A |P (A) -Q(A)|
where A is a set of measurable sets, for instance Borel sets.

The reason why discrepancy is of interest to Monte Carlo sampling is perhaps best exemplified through the Koksma-Hlawka inequality which bounds the Monte Carlo error by a constant (depending on the integrand) times the discrepancy of the sample used. Hence the discrepancy alone leads the asymptotic error rate, hence the obvious aim of seeking lowest discrepancy. The best methods produce discrepancy of order O(N -1+ ) for any > 0, which is of course better than the Monte Carlo rate O P (N -1/2 ).

In order to utilize the potential of QMC in the context of SMC, Gerber & Chopin introduce a sequential quasi Monte Carlo (SQMC) methodology. This assumes the existence of transforms Γ t mapping uniform random variables to the state variables. More specifically, it is assumed that (1) can be rewritten as x

(n) 0 = Γ 0 (u (n) 0 ) to generate x (n) 0 ∼ f X 0 (dx (n) 0 )
, and as

x (n) 1:t = Γ t (x (n) 1:t-1 , u (n) t ) to generate x (n) 1:t |x (n) 1:t-1 ∼ f X t (dx (n) 1:t |x (n) 1:t-1 ), where u (n) t ∼ U([0, 1) d ).
We show in the next section that such a deterministic function Γ t that is easy to evaluate is available in the context of Dirichlet process mixture models, however at the cost of Markovianity.

Sequential quasi Monte Carlo for Infinite Mixture Models

Nonparametric mixtures are commonly used for density estimation and can be thought of as an extension of finite mixture models when the number of clusters is unknown. Observations y 1:T follow a DPM model with kernel ψ, ie y → ψ(y; θ) are probability density functions for any θ in a parameter set Θ,

y t |G i.i.d. ∼ ψ(y; θ)dG(θ), t ∈ (1 : T ),
where G is endowed with a Dirichlet process prior DP(α, G 0 ) with precision parameter α and base measure G 0 . The discrete structure of the Dirichlet process induces partitions on the observations y 1:t into clusters identified by the cluster labels x t , or allocation variables. They link each observation y t to a mixture component θ * xt . For any t, y 1:t are clustered into a set of k t ≤ t clusters of sizes n t,j , for j ∈ (1 : k t ). The allocation variables can be interpreted as the states of y 1:t in the context of SMC. Note that they are discrete and elements of {1, . . . , k t }. The generalized Pólya urn scheme (GPUS) is essentially a sequential formulation of this model: it describes the joint distribution of the states x 1:t by its sequence of prior predictive distributions pt,j = P(x t = j|x 1:t-1 ) for t ∈ (1 : T ). Interestingly, the posterior predictive distribution also takes the form of a GPUS p t,j = P(x t = j|x 1:t-1 , y 1:t ). Both prior and posterior predictives, as well as details on the model, are provided in Appendix, Section A.

Our SMC sampler follows similar lines as [START_REF] Fearnhead | Particle filters for mixture models with an unknown number of components[END_REF]. The specific resampling technique, known to be unavoidable in sequential samplers, is detailed in Appendix, Section B. As stressed in Section 2, the key ingredient of SQMC is the availability of a deterministic transform Γ t to sample x t given x t-1 with the use of an additional quasi uniform variable. Here, such a function is available, however at the cost of Markovianity, thus using x 1:t-1 instead of x t-1 . This transform utilizes the posterior predictive distribution of Equation 3 and takes the following form in our setting

Γ t (x (n) 1:t-1 , u (n) t ) = min j ∈ {1, . . . , k (n) t-1 + 1} : j i=1 p (n) t,i > u (n) t
for any particle n, where p (n) t,i is the associated GPUS weight (3) and with u t ∼ U([0, 1)) (see [START_REF] Arbel | Discussion of "Sequential Quasi-Monte Carlo" by Gerber and Chopin[END_REF]. Although the particles state-space dimension is growing with t, notice that this second argument of the Γ t transform remains a univariate uniform random variable. This is of crucial importance for our methodology since SQMC samplers are known to break for large dimensional u t [START_REF] Gerber | Sequential quasi Monte Carlo[END_REF]. Thus, Γ t relies on quasi uniform vectors u

(1:N ) t ∈ [0, 1) N . Generating such a vector is straightforward, for instance by adding a random uniform variable to a regular grid of [0, 1) modulo 1, such as u (1:N ) t = u + (0, 1/N, . . . , (N -1)/N ) mod.1, with u ∼ U[0, 1). This is also faster than generating N uniform random variables required in the SMC setting, since it only needs one such variable.

Impact of quasi Monte Carlo on diversity of the allocation variables

We compare some properties of the allocation variables trajectories x (n) 1:T , n = 1, . . . , N of DP normal mixture models under three samplers: the non sequential Monte Carlo method (Gibbs sampler) of [START_REF] Neal | Markov chain sampling methods for Dirichlet process mixture models[END_REF] (MC), a sequential Monte Carlo method of [START_REF] Griffin | Sequential Monte Carlo methods for mixtures with normalized random measures with independent increments priors[END_REF] (SMC) and our proposed sequential quasi Monte Carlo method (SQMC)3 . We run the three samplers on three simulated datasets of size T = 200, respectively with heavy-tailed, skewed, and multimodal distributions (see densities on Figure 1). The number of iterations for MC (after a burn-in period) and the number of particles in both sequential samplers is set to the same value of N = 1000.

For the three distributions, the posterior mean obtained for the three methods have equivalent fit (see left panel of Figure 1). However, we observe that the SMC method does not satisfactorily explore the allocation space. In particular, we observe that a single path for the allocation vector x 1:t is strongly favored by this approach up to some t ≤ T . Using quasi random sequences is a simple way to get around this problem by insuring at each step a good spread of the random sequence. In practice, this makes a significant difference and we observe allocations x (n) 1:T that present closer behavior in terms of dispersion to the one obtained with non sequential method. To measure that dispersion, we run a principal component analysis (PCA) on the allocation vector x (n) 1:T , n = 1, . . . , N . The PCA gives information on how different the trajectories are from one another. We represent on the right panel of Figure 1 the proportion of variance explained by number of components in the PCA. For instance, if the samples have similar allocation vectors, then most of the variance will be on few principal components, yielding a high curve.

The first two examples, heavy-tailed and skewed distributions, cannot be written straightforwardly as a mixture of normal kernels. Thus in these cases, we do not expect a clear clustering but some heterogeneity in the sampled allocation vectors. On the contrary, the third example is exactly a mixture of normal kernels, so there should be much less variability in the allocation vector here. From the experiments, we observe that the diversity of the allocation trajectories under SQMC behaves as one would expect: it is high in the first two examples, and moderate in the third. This feature is not as clear for SMC. Additionally, having more diversity in the allocation vector is a sign of a better exploration of the latent space. This diversity for SQMC is always higher than the one for SMC, and is closer to the one obtained from a non sequential approach (MC), indicating that SQMC outperforms SMC in terms of latent space exploration.

Discussion

We introduce sequential quasi Monte Carlo sampling for DPM models, and discuss the applicability of the proposed methodology to a broader class of BNP mixture models. For these models, the structure of the latent space is non standard for sequential methods. We show that favoring diversity by using quasi random sequences improves the exploration of the allocation space. leads to an increasing variance and so-called particle degeneracy. A popular criterion for measuring degeneracy is the effective sample size (ESS) defined by

N n=1 W (n) t 2 -1 .
which takes on values in (1, N ), and it is common practice to resample the particles when the ESS falls below a threshold. Once resampled, the particles are all given the same weight (1/N ). It is common to require the resampling to be unbiased, meaning that the expected value of any function of the particles is unchanged by the resampling. The typical resampling technique consists in sampling particles from x

(1:N ) 1:t using the multinomial distribution Multinomial N ; W , which leads to duplication of particle values with positive probability (see for instance [START_REF] Gordon | Novel approach to nonlinear/non-Gaussian Bayesian state estimation[END_REF]. In the case of discrete (finite) state-space, specific resampling schemes have been tailored in order to take advantage of the finiteness of the states with the aim of reducing wasteful particle duplication. One defines as 'putative particles' the set of M t = N n=1 k (n) t-1 + 1 possible particles at step t.

Computing their weights is feasible, and sampling from them instead of from the N < M t particles of step t -1 improves efficiency from a Rao-Blackwell point of view. However, keeping track of all possible trajectories till time t is a task beyond computational capabilities even for limited t, since the trajectories number is equal to the number of partitions of a set of size t which explodes essentially as t t . Hence resampling appears to be unavoidable. The approach that we follow consists in the two steps

• Computing the normalized weights of the M t putative particles, x (n) 1:t-1 , j , for n ∈ (1 : N ) and j ∈ (1 : (k

(n) t-1 + 1)), by W (n) t,j ∝ W (n) t-1 p x (n) 1:t-1 , j|y (n) 1:t p x (n) 1:t-1 |y (n) 1:t-1
where the ratio in the right hand side, also known as the incremental weight, can be evaluated up to a normalizing constant in conjugate models.

• Choosing N among the M t by the multinomial Multinomial N ; W (n) t,j , n ∈ (1 : N ), j ∈ (1 : (k (n) t-1 + 1)) .

Fig. 1 :

 1 Fig. 1: Left: Density fit; Right: Diversity index (PCA) for three samplers, non sequential Monte Carlo (MC), sequential Monte Carlo (SMC) and sequential quasi Monte Carlo (SQMC). Top row: Heavy tailed distribution (student with degree of freedom two); Middle row: Skewed distribution (log-Gamma); Bottom row: Multimodal distribution with three well separated modes (mixture of normals).

The code is available upon request.

A Generalized Pólya urn schemes for infinite mixture models

The generalized Pólya urn scheme (GPUS) is essentially a sequential formulation of this model: it describes the joint distribution of the states x 1:t by its sequence of prior predictive distributions pt,j = P(x t = j|x 1:t-1 ) for t ∈ (1 : T ). For a Dirichlet process prior DP (α, G 0 ) with precision parameter α and base measure G 0 pt,j ∝ n t-1,t for j ∈ (1 :

The posterior predictive distribution also takes the form of a GPUS p t,j = P(x t = j|x 1:t-1 , y 1:t ) with

where

and

where one denotes the allocation sets by s t,j = {l ≤ t : x l = j}. The weights ψ j (y t |x 1:t-1 ) and ψ 0 (y t ) can be calculated in conjugate cases and approximated otherwise. Note that posterior predictive distributions in the form of GPUS are also available for broader classes of random probability measures including the two-parameter Poisson-Dirichlet process [START_REF] Pitman | The two-parameter Poisson-Dirichlet distribution derived from a stable subordinator[END_REF] and normalized random measures with independent increments [START_REF] Regazzini | Distributional results for means of normalized random measures with independent increments[END_REF][START_REF] James | Posterior analysis for normalized random measures with independent increments[END_REF]. Algorithms are detailed by [START_REF] Griffin | Sequential Monte Carlo methods for mixtures with normalized random measures with independent increments priors[END_REF].

The most widespread setting [START_REF] Escobar | Bayesian density estimation and inference using mixtures[END_REF][START_REF] Fearnhead | Particle filters for mixture models with an unknown number of components[END_REF], that we also pursue here, considers Gaussian mixtures with a Gaussian base measure G 0 . More specifically, we consider location scale mixtures where the hierarchical representation of the model assumes that observations in each cluster are draw iid from Gaussian distributions N (µ, σ 2 ). The conjugate prior for θ = (µ, s) where s = σ -2 has the conditional representation s ∼ Gamma(a, b) and µ|s ∼ N (η, τ /s), with fixed hyperparameters a, b, η and τ . Posterior distributions of µ j and s j for cluster j parameters, as well as the posterior probability of the states x 1:t are given by [START_REF] Fearnhead | Particle filters for mixture models with an unknown number of components[END_REF].

Marginalising out the location scale parameter θ has been shown to improve substantially the efficiency of sequential samplers [START_REF] Chen | Mixture Kalman filters[END_REF]. Indeed, it allows us to consider only the allocations x 1:t as the states, thus removing the part on θ which would make our QMC methodology more cumbersome.

B Resampling algorithm

A typical feature of sequential methods is that the weights W (n) t can become very skewed as t increases, ie most of the particles weights tend to zero while only some remain active, which


Julyan Arbel
email: julyan.arbel@inria.fr

Igor Prünster

Truncation error of a superposed gamma process in a decreasing order representation

Keywords: Approximate inference, Bayesian Nonparametrics, Completely random measures, Ferguson & Klass algorithm, Normalized random measures, Truncation

), we study a moment-matching criterion which consists in evaluating a measure of discrepancy between actual moments of the CRM and moments based on the simulation output. This note focuses on a general class of CRMs, namely the superposed gamma process, which suitably transformed have already been successfully implemented in Bayesian Nonparametrics. To this end, we show that the moments of this class of processes can be obtained analytically.

Introduction

Completely random measures (CRMs) represent a fundamental building block of countless popular stochastic models and play a prominent role within Bayesian nonparametric modeling (see [START_REF] Lijoi | Models beyond the Dirichlet process[END_REF][START_REF] Jordan | Hierarchical models, nested models and completely random measures[END_REF]. For instance, the popular Dirichlet process [START_REF] Ferguson | A Bayesian analysis of some nonparametric problems[END_REF] can be obtained as normalization or exponentiation of suitable CRMs. Implementation of CRM-based models usually requires to simulate the CRMs trajectories. As infinite dimensional objects, representable as infinite random series (1), truncation of the series is needed leading to an approximation error. Different strategies for representing the series and choosing the truncation threshold M are available in the literature. A thorough account on CRM size-biased representations has been recently provided in [START_REF] Campbell | Truncated Random Measures[END_REF], where the truncation error is examined from the angle of the L 1 distance between marginal data densities under truncated and untruncated models, an approach first considered by [START_REF] Ishwaran | Gibbs sampling methods for stick-breaking priors[END_REF] for stick-breaking processes.

30th Conference on Neural Information Processing Systems (NIPS 2016), Barcelona, Spain.

Another important type of representation due to [START_REF] Ferguson | A representation of independent increment processes without Gaussian components[END_REF] (see also [START_REF] Walker | Representations of Lévy processes without Gaussian components[END_REF] is arguably one of the most useful ones in that it displays the weights in decreasing order. This yields the very interesting feature that the approximation error is minimized over the whole sample space for a given truncation level. This competitive advantage was exploited in many works, including [START_REF] Argiento | A blocked Gibbs sampler for NGGmixture models via a priori truncation[END_REF]; [START_REF] Barrios | Modeling with normalized random measure mixture models[END_REF][START_REF] De Blasi | A class of neutral to the right priors induced by superposition of beta processes[END_REF]; [START_REF] Epifani | Exponential functionals and means of neutralto-the-right priors[END_REF]; [START_REF] Griffin | Posterior simulation of normalized random measure mixtures[END_REF]Nieto-Barajas et al. (2004); Nieto- Barajas and Walker (2004) to cite just a few in Bayesian Nonparametrics. The quality of the approximation, addressed only heuristically in those previous works, is the focus of this note.

In [START_REF] Arbel | A moment-matching Ferguson & Klass algorithm[END_REF] it shown how moments of the CRMs can be used in order to assess the quality of approximation due to the truncation. Since moments of CRMs are simple to compute, the quality of the approximation can be quantified by evaluating a measure of discrepancy between the actual moments of the CRM at issue and the moments computed based on the sampled realizations of the CRM. The truncation level is selected so that measure of discrepancy does not exceed a given threshold, say 5%. In [START_REF] Arbel | A moment-matching Ferguson & Klass algorithm[END_REF] the methodology is illustrated on two classes of CRMs, namely the generalized gamma process and the stable-beta process.

In the proposed talk we review the results of [START_REF] Arbel | A moment-matching Ferguson & Klass algorithm[END_REF] and, in addition, analyze another broad class called the superposed gamma process (see [START_REF] Regazzini | Distributional results for means of normalized random measures with independent increments[END_REF][START_REF] Lijoi | Bayesian nonparametric analysis for a generalized Dirichlet process prior[END_REF]. After defining CRMs and the superposed gamma process (Section 2), we display a bound in probability on the truncation error in Section 3.1 and then show the applicability of the moment-matching criterion by deriving analytically the moments of the superposed gamma process in Section 3.2. An illustration of the Ferguson & Klass algorithm and a proof are deferred to the Appendix.

Superposed gamma process

Completely random measures (CRM) have been successfully employed in a wide spectrum of modern applications, including survival analysis, random sparse networks, biology, to cite just a few. A CRM μ on X is a random measure which spreads out mass independently in the space, which means that the random measure μ, when evaluated on disjoint sets A 1 , . . . , A n , leads to independent random variables μ(A 1 ), . . . , μ(A n ). Importantly, [START_REF] Kingman | Completely random measures[END_REF] showed that the only way to spread out mass in a completely random fashion (without deterministic components) is by randomly scattering point masses in the space. In other words, CRM present the interesting feature that they select (almost surely) discrete measures and hence can be represented as μ =

i≥1 J i δ Z i (1)
where both the jumps J i and the locations Z i are random and are controlled by the so-called Lévy intensity which characterizes the CRM. It is a measure on R + × X which can be written as ν(dv, dx) = ρ(dv)α(dx) for so-called homogeneous CRM, which are considered here and correspond to the case of jumps independent of the locations. The function ρ controls the intensity of the jumps. The measure α, if the CRM is (almost surely) finite, which is assumed throughout, splits up in α = aP 0 where a > 0 is called the total mass parameter and the probability distribution P 0 tunes the locations.

Ever-popular CRM include the generalized gamma process introduced by Brix (1999) and the stable-beta process, or three-parameter beta process, defined by [START_REF] Teh | Indian buffet processes with power-law behavior[END_REF] as an extension of the beta process [START_REF] Hjort | Nonparametric bayes estimators based on beta processes in models for life history data[END_REF]). Here we consider another large class of completely random measures called superposed gamma process, introduced by [START_REF] Regazzini | Distributional results for means of normalized random measures with independent increments[END_REF]. It is identified by the jump intensity

ρ(dv) = 1 -e -ηv 1 -e -v e -v v dv, η > 0.
(2)

As noted by [START_REF] Lijoi | Bayesian nonparametric analysis for a generalized Dirichlet process prior[END_REF], one usually restricts attention to the case of positive integer η. Under this assumption, the superposed gamma process takes the form of a genuine superposition of independent gamma processes with increasing integer-valued scale parameter, with jump intensity ρ(dv) = 1 v e -v + e -2v + . . . + e -ηv dv/v. The specification of integer values for η has also the advantage to lead to analytic computation of the moments. Note that the special case η = 1 reduces to the gamma process, which gives rise to the Dirichlet process by normalization. Alternatively, the normalization of the superposed gamma process for unspecified η provides the so-called generalized Dirichlet process [START_REF] Lijoi | Bayesian nonparametric analysis for a generalized Dirichlet process prior[END_REF]. [START_REF] Ferguson | A representation of independent increment processes without Gaussian components[END_REF] devise a constructive representation of a CRM which produces the jumps in decreasing order. This corresponds to the (almost surely unique) ordering of the sum elements in (1) where J 1 > J 2 > • • • . Indeed, the jumps are obtained as andξ 1 , ξ 2 , ∼ Exp(1). Figure 2 in Appendix illustrates the function N ( • ) which takes the following form in the superposed gamma process case

ξ i = N (J i ), where N (v) = ν([v, ∞), R) is a decreasing function,
N (v) = aE η (v), where E η (v) = η l=1 E 1 (kv) and E 1 (v) = ∞ v u -1 e -u du, (3)
and where the function E 1 denotes the exponential integral function. Since it is impossible to sample an infinite number of jumps, approximate simulation of μ is in order. This becomes a question of determining the number M of jumps to sample leading to the truncation μM and truncation error T M as follow

μM = M i=1 J i δ Z i , T M = ∞ i=M +1 J i . ( 4
)
The Ferguson & Klass representation has the key advantage of generating the jumps in decreasing order implicitly minimizing such an approximation error. However, a precise evaluation of T M , for example in expectation, is a daunting task due to the non independence of the jumps in the Ferguson & Klass representation. The algorithm is summarized in Algorithm 1.

Algorithm 1: Ferguson & Klass algorithm 1: sample ξ i ∼ PP for i = 1, . . . , M 2: define J i = N -1 (ξ i ) for i = 1, . . . , M 3: sample Z i ∼ P * for i = 1, . . . , M 4: approximate μ by M i=1 J i δ Z i 3 Truncation error of the superposed gamma process

Bound in probability

We provide an evaluation in probability of the truncation error T M in (4). Proposition 1. Let (ξ j ) j≥1 be the jump times for a homogeneous Poisson process on R + with unit intensity. Then for any ∈ (0, 1), the tail sum of the superposed gamma process (4) satisfies

P T M ≤ t M ≥ 1 -, for t M = C (η!) 1/η e 1-M C , where C = 2eaη .
The simple proof is deferred to Appendix B. It is interesting to note that the bound t M for the superposed gamma process is equal to its counterpart for the beta process with concentration parameter c set to η, all else things being equal (total mass parameter a and threshold ). See Proposition 1 in [START_REF] Arbel | A moment-matching Ferguson & Klass algorithm[END_REF]. This finding provides a nice connection between both processes otherwise seemingly unrelated.

The bound t M obtained in Proposition 1 is exponentially decreasing with M , which is reminiscent of the results obtained by [START_REF] Brix | Generalized gamma measures and shot-noise Cox processes[END_REF] and [START_REF] Arbel | A moment-matching Ferguson & Klass algorithm[END_REF], respectively, for the generalized gamma process and the stable-beta process with no stable component. As already pointed out by these authors, the bound t M is very conservative due to a crude lower bound on the quantiles q j (notation of the proof). The left panel of Figure 1 displays this bound t M , while the right panel illustrates the truncation level M (in log-scale) required in order to guarantee with 95% probability an upper bound on T M of t max ∈ {1, 10, 100}, for varying values of η. Inspection of the plots demonstrates the rapid increase with η of the number of jumps needed in order to assess a given bound in probability.

Moment-matching criterion

Given the limited practical usefulness of the bound provided in Proposition 1, we propose to use the alternative route of the moment-matching methodology of [START_REF] Arbel | A moment-matching Ferguson & Klass algorithm[END_REF] and derive the relevant quantities in order to implement it for the superposed gamma process. Let us consider the n-th (raw) moment of the (random) total mass μ(X) defined by m n = E μn (X) . Given the Laplace transform of CRM, the moments (see, e.g., Proposition 1 in [START_REF] Arbel | A moment-matching Ferguson & Klass algorithm[END_REF] take on the form

m n = ( * ) ( n k 1 •••kn ) n i=1 κ i /i! k i , (5)
where the sum ( * ) is over all n-tuples of nonnegative integers (k 1 , . . . , k n ) satisfying the constraint dv). Simple algebra leads to the following expression for the cumulants of the superposed gamma process

k 1 + 2k 2 + • • • + nk n = n and κ i is the ith cumulant defined by κ i = a ∞ 0 v i ρ(
κ i = a(i -1)!ζ η (i) (6)
which displays the incomplete Euler-Riemann zeta function ζ η (i) = η l=1 1 l i . Hence the moment-matching methodology introduced by Arbel and Prünster (2016) can be readily applied by making use of ( 5) and (6).

A Ferguson and Klass illustration

From Figure 2 it is apparent that increasing η leads to larger jumps which in turn leads to the need of a higher truncation level in order to match a given precision level. This is not surprising given the CRM at hand can be thought of as a superposition of η gamma CMRs. Such an intuition is made precise in Section 3 in the main text by (i) deriving a bound in probability on T M and (ii) obtaining the moments of the superposed gamma process, thus showing the applicability of the moment-criterion introduced by Arbel and Prünster (2016).

0 J 5 J 4 J 3 J 2 J 1 0 ξ 1 ξ 2 ξ 3 ξ 4 ξ 5 0

B Proof

Proof of Proposition 1. The proof follows along the same lines as the proof of Theorem A.1. by [START_REF] Brix | Generalized gamma measures and shot-noise Cox processes[END_REF] for the generalized gamma process and Proposition 4 by [START_REF] Arbel | A moment-matching Ferguson & Klass algorithm[END_REF] for the stable-beta process. Let q j denote the 2 M -j quantile, for j = M + 1, M + 2, . . ., of a gamma distribution with mean and variance equal to j. Then

P ∞ j=M +1 N -1 (ξ j ) ≤ ∞ j=M +1
N -1 (q j ) ≥ 1 -.

Denote t M = ∞ j=M +1 N -1 (q j ) = ∞ j=M +1 E -1 η (q j /a), and let us upper bound E -1 η . By using E 1 (u) ≤ 1 -log(u), one gets

E η (u) = η l=1 E 1 (lu) ≤ η - η l=1 log(lu) = η -log η!u η ,
which can be inverted to obtain

E -1 η (x) ≤ 1 (η!) 1/η e 1-x η .
Additionally, since the quantiles satisfy q j ≥ 2e j, we can conclude that

t M ≤ 1 (η!) 1/η ∞ j=M +1 e 1- q j aη ≤ 1 (η!) 1/η ∞ j=M +1
e 1-j 2eaη ≤ 2eaη (η!) 1/η e 1-M 2eaη .

  . . . are jump times of a standard Poisson process (PP) of unit rate: ξ 1 , ξ 2 -ξ 1 , . . . i.i.d.

Fig. 1 :

 1 Fig. 1: Left: variation of M → t M for η ∈ {1, 2, 5, 10}. Right: variation of the threshold function η → M needed to match an error bound of t max ∈ {1, 10, 100} with η ∈ {1, . . . , 20}, log scale on y-axis.

Fig. 2 :

 2 Fig. 2: Left: illustration of Ferguson and Klass representation through the inversion of the jumps times ξ 1 , . . . , ξ 5 for a homogeneous Poisson process on R + to the jumps J 1 , . . . , J 5 of a CRM. Right: Tail of the Lévy measure N ( • ) of the superposed gamma process with η ∈ {1, . . . , 10}, η = 1 for the lowest curve, η = 10 for the highest curve.


Guillaume Kon

Kam King

Julyan Arbel
email: julyan.arbel@inria.fr

Igor Prünster

A Bayesian nonparametric approach to ecological risk assessment

Keywords: Bayesian Nonparametrics, Ecotoxicology, HC 5, Mixture models, Normalized random measures, Species Sensitivity Distribution

, where the method is generalised to censored data and a systematic comparison on simulated data is also presented, along with a study of the clustering induced by the mixture model to examine patterns in species sensitivity.

Introduction

Assessing the response of a community of species to an environmental stress is critical for ecological risk assessment. Methods for this purpose vary in levels of complexity and realism. Species Sensitivity Distribution (SSD) represents an intermediate tier, more refined than rudimentary assessment factors [START_REF] Posthuma | Species sensitivity distributions in ecotoxicology[END_REF] but practical enough for routine use by environmental managers and regulators in most developed countries (Australia, Canada, China, EU, South Africa, USA. . . ). The SSD approach is intended to provide, for a given contaminant, a description of the tolerance of all species possibly exposed using information collected on a sample of those species. This information consists of Critical Effect Concentrations (CECs), a concentration specific to a species which marks a limit over which the species suffers a critical level of effect. This is for instance the concentration at which 50% of the tested organisms died (Lethal Concentration 50% (LC 50 )), or the concentration which inhibited growth or reproduction by 50% compared to the control experiment (Effect Concentration 50% (EC 50 )). Each CEC is the summary of long and costly bioassay experiments for a single species, so they are rarely available in large number. Typical sample sizes range from 10 to 15 (ECHA, 2008).

To describe the tolerance of all species to be protected, the distribution of the CECs is then estimated from the sample. In practice, a parametric distributional assumption is often adopted [START_REF] Forbes | Species Sensitivity Distributions Revisited: A Critical Appraisal[END_REF]: the CECs are assumed to follow a log-normal [START_REF] Aldenberg | Uncertainty of the hazardous concentration and fraction affected for normal species sensitivity distributions[END_REF], log-logistic [START_REF] Kooijman | A safety factor for LC 50 values allowing for differences in sensitivity among species[END_REF]), triangular (Van Straalen, 2002;[START_REF] Zhao | Species sensitivity distribution for chlorpyrifos to aquatic organisms: Model choice and sample size[END_REF] or BurrIII [START_REF] Shao | Estimation for hazardous concentrations based on NOEC toxicity data: An alternative approach[END_REF] distribution.

Once the response of the community is characterised by the distribution, the goal of risk assessment is to define a safe concentration protecting all or most of the species. In the case of distributions without a lower threshold strictly above 0, a cut-off value is often chosen as the safe concentration. Typically, this is the Hazardous Concentration for 5% of the Species (HC 5 ), which is the 5th percentile of the distribution. Reasonings behind this choice include: that the lowest bound of the confidence interval around the 5th percentile will be used instead of the estimate, that a safety factor will be subsequently applied to that value and that ecosystems have a certain resilience to perturbations.

The lack of justification for the choice of any given parametric distribution has sparked several research directions. Some authors [START_REF] Xu | Key issues for the development and application of the species sensitivity distribution (SSD) model for ecological risk assessment[END_REF][START_REF] He | Ecological risk assessment and priority setting for typical toxic pollutants in the water from Beijing-Tianjin-Bohai area using Bayesian matbugs calculator (BMC)[END_REF][START_REF] Jagoe | Bootstrap estimation of community NOEC values[END_REF][START_REF] Van Straalen | Threshold models for species sensitivity distributions applied to aquatic risk assessment for zinc[END_REF][START_REF] Xing | A comparison of statistical methods for deriving freshwater quality criteria for the protection of aquatic organisms[END_REF][START_REF] Zhao | Species sensitivity distribution for chlorpyrifos to aquatic organisms: Model choice and sample size[END_REF]) have sought to find the best parametric distribution by model comparison using goodness-of-fit measures. The general understanding is that no single distribution seems to provide a superior fit and that the answer is dataset dependent [START_REF] Forbes | Species Sensitivity Distributions Revisited: A Critical Appraisal[END_REF]. Therefore, the log-normal distribution has become the customary choice, notably because it readily provides confidence intervals on the HC 5 , and because model com-parison and goodness of fit tests have relatively low power on small datasets, precluding the emergence of a definite answer to the question. Another research direction consisted in seeking to avoid any reference to a distribution, using so-called nonparametric or distribution-free approaches. Those efforts included using the empirical distribution function [START_REF] Suter | Ecological Risk Assessment in a Large River-Reservoir: 2. Fish Community[END_REF][START_REF] Jones | Ecological risk assessment in a large river-reservoir: 3. Benthic invertebrates[END_REF], methods based on ranks [START_REF] Van Der Hoeven | Estimating the 5-percentile of the species sensitivity distributions without any assumptions about the distribution[END_REF][START_REF] Chen | A conservative, nonparametric estimator for the 5th percentile of the species sensitivity distributions[END_REF], bootstrap resampling [START_REF] Jagoe | Bootstrap estimation of community NOEC values[END_REF][START_REF] Wang | Development of species sensitivity distributions and estimation of HC(5) of organochlorine pesticides with five statistical approaches[END_REF] or nonparametric kernel density estimation [START_REF] Wang | Non-parametric kernel density estimation of species sensitivity distributions in developing water quality criteria of metals[END_REF]. All these approaches have in common that they require large sample sizes to be effectively applicable. Finally, authors have considered the possibility that the distribution of the CECs might not be a single distribution but rather a mixture of distributions [START_REF] Zajdlik | Estimating Water Quality Guidelines for Environmental Contaminants Using Multimodal Species Sensitivity Distributions: A Case Study with Atrazine[END_REF], datasets being an assemblage of several log-normally distributed subgroups [START_REF] Kefford | Global scale variation in the salinity sensitivity of riverine macroinvertebrates: Eastern Australia, France, Israel and South Africa[END_REF][START_REF] Craig | Exploring novel ways of using species sensitivity distributions to establish PNECs for industrial chemicals: Final report to Project Steering Group[END_REF]. This is more realistic from an ecological point of view because several factors influence the tolerance of a species to a contaminant such as the taxonomic group or the mode of action, and contaminant such as pesticides might even target specific species groups. Therefore, there is strong evidence in favour of the presence of groups of CECs, although the CECs within a group might remain log-normally distributed.

Ignorance of the group structure is a strong motivation for a nonparametric approach. However, the method must remain applicable to small datasets, which suggests trying to improve on the existing frequentist nonparametric methods. Bayesian nonparametric mixture models offer an interesting solution for both large and small datasets, because the complexity of the mixture model adapts to the size of the dataset. It offers a good compromise between a simplistic one-component parametric model and a kernel density method which in a certain sense lacks flexibility and might cause overfitting. Moreover, the low amount of information available in small datasets to estimate the groups parameters can be complemented via the prior, as some a priori degree of information is generally available from other species or contaminants [START_REF] Awkerman | Development of species sensitivity distributions for wildlife using interspecies toxicity correlation models[END_REF][START_REF] Craig | Exploring novel ways of using species sensitivity distributions to establish PNECs for industrial chemicals: Final report to Project Steering Group[END_REF][START_REF] Craig | Species nonexchangeability in probabilistic ecotoxicological risk assessment[END_REF]. This paper summarises the results of the complete study in Kon Kam [START_REF] Kam King | Species Sensitivity Distribution revisited: a Bayesian nonparametric approach[END_REF].

The rest of the article is organised as follows. In Section 2 we present the Bayesian nonparametric (BNP) model and existing frequentist models for SSD and explain how to obtain a density estimate. Then in Section 3 we compare the different methods on a real dataset, illustrating the benefits of the BNP SSD. We conclude with a final discussion in Section 4.

Models for SSD

Given that concentrations vary on a wide range, it is common practice to work on log-transformed concentrations. Consider a sample of n log-concentrations denoted by X = (X 1 , . . . , X n ). We propose to carry out density estimation for the SSD based on sample X by use of nonparametric mixtures. Bayesian nonparametric mixtures were introduced in [START_REF] Lo | On a Class of Bayesian Nonparametric Estimates: I. Density Estimates[END_REF] with Dirichlet process mixtures (DPM). Generalizations of the DPM correspond to allowing the mixing distribution to be any discrete nonparametric prior. A large class of such prior distributions is obtained by normalizing increasing additive processes [START_REF] Sato | Lévy Processes and Infinitely Divisible Distributions[END_REF]. The normalization step, under suitable conditions, gives rise to so-called normalized measures with independent increments (NRMI) as defined by [START_REF] Regazzini | Distributional results for means of normalized random measures with independent increments[END_REF], see also [START_REF] Barrios | Modeling with Normalized Random Measure Mixture Models[END_REF] for a recent review. An NRMI mixture model is defined hierarchically as:

X i |µ i , σ ind ∼ k(•|µ i , σ), µ i | P i.i.d. ∼ P , i = 1, . . . , n, (1) P ∼ NRMI, σ ∼ Ga(a σ , b σ ).
where k is a kernel, which we assume parametrized by some θ = (µ, σ) ∈ R × R + , and P is a random probability on R whose distribution is an NRMI. In our model, all clusters have a common variance. This is easier to fit on a small dataset, because information about the variance is pooled across clusters. Similar mixture SSD models described in Craig ( 2013) also assume common variance. As described in the Introduction, concentrations are commonly fitted with a log-normal distribution. Our aim is to move from this parametric model to the nonparametric one in (1). In order to allow comparisons to be made, we stick to the normal specification for k on the log-concentrations X by letting: k(x|µ, σ) = N (x|µ, σ). Under this framework, density estimation is carried out by evaluating the posterior predictive density along the lines of [START_REF] Barrios | Modeling with Normalized Random Measure Mixture Models[END_REF]:

f (x| P , X) = k(x|µ, σ)dπ(σ)d P (µ) (2)
for any x in R, where π denotes the posterior distribution of σ.

To specify the prior, we choose as mixing random measure the normalized stable process [START_REF] Kingman | Random discrete distributions[END_REF] with: i a stability parameter γ = 0.4, which controls the flatness of the prior on the number of clusters. The parameter γ can take values in (0, 1). Taking the limit γ → 0 reduces the model to a Dirichlet process, larger values of γ lead to less informative priors on the number of clusters. The parameter γ was chosen as a good compromise between model flexibility and numerical stability. The total mass parameter is, without loss of generality, set equal to 1. ii a base measure (which corresponds to the mean of the random probability measure)

P 0 ( • ) = N ( • | ϕ 1 , ϕ 2 )
with mean ϕ 1 and standard deviation ϕ 2 , hyperparameters fixed a priori to specify a certain knowledge in the degree of smoothness iii a common variance for all the clusters with a vaguely informative prior distribution Ga(0.5, 0.5).

Recent years have witnessed the appearance of a wealth of softwares dedicated to implement Bayesian nonparametric models and sample from their posterior. To cite a few, the R package DPpackage [START_REF] Jara | DPpackage: Bayesian Semi-and Nonparametric Modeling in R[END_REF], is a rather comprehensive bundle of functions for Bayesian nonparametric models, while Bayesian Regression [START_REF] Karabatsos | A menu-driven software package of Bayesian nonparametric (and parametric) mixed models for regression analysis and density estimation[END_REF]) is a software for Bayesian nonparametric regression. For posterior sampling, we use the R package BN-Pdensity and the function MixNRMI1 which implements BNP density models under a general specification of normalized random measures based on the generalised gamma processes (see [START_REF] Barrios | Modeling with Normalized Random Measure Mixture Models[END_REF]. The package is available from the Comprehensive R Archive Network (CRAN).

To illustrate the interest of the Bayesian nonparametric SSD, we compare our proposed BNP model to two commonly used frequentist models: the normal distribution [START_REF] Aldenberg | Uncertainty of the hazardous concentration and fraction affected for normal species sensitivity distributions[END_REF] and the nonparametric Kernel Density Estimate (KDE) recently proposed by [START_REF] Wang | Non-parametric kernel density estimation of species sensitivity distributions in developing water quality criteria of metals[END_REF]. For both frequentist approaches, the data is assumed to be iid. Density estimates take on respectively the following form (μ and σ are MLE)

fN (x) = N (x | μ, σ) and fKDE (x) = 1 n n i=1 N (x | X i , 1.06σn -1 5 ). (3)

Model comparison and cross-validation

For the purpose of comparing the predictive performance of the model, we resort to Leave-One-Out (LOO) cross-validation. We compute the LOOs for each of the methods as ∀i,

LOO i = f (X i | X -i ) where f (x | X -i )
is the density for one of the three methods estimated from X with X i left out. The LOOs for the BNP model correspond to the conditional predictive ordinates (CPOs) statistics which are commonly used in applications, see [START_REF] Gelfand | Model determination using sampling-based methods[END_REF]. A CPO statistic is defined for each log-concentration X i as follows:

CPO i = f (X i |X -i ) = k(X i |θ)dπ(θ|X -i ) (4)
where X -i denotes the all sample X but X i , dπ(θ|X -i ) is the posterior distribution associated to X -i and f is the (cross-validated) posterior predictive distribution of Equation ( 2). As shown by [START_REF] Barrios | Modeling with Normalized Random Measure Mixture Models[END_REF], CPOs can be easily approximated by Monte Carlo as

CPO i = 1 T T t=1 1 k(X i |θ (t) ) -1 (5)
where {θ (t) , t = 1, 2, . . . , T } is an MCMC sample from the posterior distribution.

Quantile estimation and HC 5

The quantity of interest for ecological risk assessment is the HC 5 , which corresponds to the 5th percentile of the SSD distribution. We choose as an estimator the median of the posterior distribution of the 5th percentile, while the 95% credible bands are formed by the 2.5% and 97.5% quantiles of the posterior distribution of the 5th percentile. The 5th percentile of the KDE is obtained by numerical inversion of the cumulative distribution function, and the confidence intervals using nonparametric bootstrap. The 5th percentile of the normal SSD and its confidence intervals are obtained following the classical method of [START_REF] Aldenberg | Uncertainty of the hazardous concentration and fraction affected for normal species sensitivity distributions[END_REF].

Application to real data

We applied this model to a selection of contaminants extracted from a large database collected by the National Institute for Public Health and the Environment (RIVM). This database was prepared, studied and published by [START_REF] Hickey | On the quantification of intertest variability in ecotoxicity data with application to species sensitivity distributions[END_REF]. We only considered non censored data, left or right censored data were discarded, while interval censored data were replaced by the centre of the interval. Kon Kam [START_REF] Kam King | Species Sensitivity Distribution revisited: a Bayesian nonparametric approach[END_REF] will describe how the method can be adapted to include censored data. Using a continuous distribution for the CECs implies that the model does not support ties (or, in other words, observing ties has zero probability). However, ties may appear in the dataset due to the rounding of concentrations. Hence, we used a small jittering of the data.

We selected three example datasets which feature three typical sample sizes: a relatively large carbaryl dataset (CAS: 63-25-2, insecticide, 55 species), a medium-sized temephos dataset (CAS: 3383-96-8, mosquito larvicide, 21 species), and a small captan dataset (CAS: 133-06-2, fungicide, 13 species). Datasets for new contaminants are always small, the minimum requirement set by the European Chemical Agency being 10 species. The datasets can be visualised on the histograms of Figure 1 (left panel).

These datasets illustrate different features of the three approaches: when there is a clear multimodality in the data, the BNP SSD is more flexible than the fixed bandwidth KDE SSD (Figure 1, carbaryl and captan). When the data do not exhibit strong multimodality, as for temephos, the BNP reduces to the normal SSD model, whereas the KDE remains by construction a mixture of many normal components.

One might think to increase the flexibility of the KDE by simply decreasing the bandwidth. However, that would also decrease the robustness of the method. On the middle panel of Figure 1, the LOOs give an indication of the robustness to over-fitting of the three methods. For carbaryl and captan, they show that the superior flexibility of the BNP SSD compared to the KDE SSD does not come at the expense of robustness, because the median CPO 0.0 0.2 0.4 0.6 0.8 0.0 2.5 5.0 7.5 q q q q q q 0.01 0.10 q q q q q q 0.4 0.8 1.2 0.0 0.2 0.4 0.6 -3 0 3 6 q q q q q q q q q q 0.01 0.10 q q q q q q -1 0 1 2 0.0 0.5 1.0 1.5 0 2 4 6

(a) Density q q q q q q q q 0.01 1.00 (b) LOO and CPO q q q q q q 0.75 1.00 1.25 1.50 1.75

(c) log HC 5 The box hinges extend to the inner quartiles. The whiskers extend to cover points up to one and a half times the inter-quartile distance away from the hinges. For both frequentist methods, the n LOOs are obtained by fitting the model n times, while an analytical expression is available for the BNP method (Equation 5). Right: log HC 5 and associated confidence/credible intervals (for Normal, KDE and BNP).

of the BNP SSD is higher than the other two. In the case of temephos, the median LOO likelihood estimate of the normal model is very similar to the median CPO for the BNP SSD, sign that there is little over-fitting. This generally illustrates the fact that model complexity in a BNP model scales with the amount and structure of the data. On the right panel of Figure 1, the credible intervals of the HC 5 for the BNP SSD are generally larger than the confidence interval of the normal SSD, which is coherent with the model uncertainty of the nonparametric approach.

Discussion

The BNP SSD seems to perform well when the dataset deviates from a normal distribution. Its great flexibility is an asset to describe the variability of the data, while it does not seem prone to over-fitting. It can be thought of as an intermediate model between the normal SSD with a single component on the one hand, and the KDE which counts as many components as there are species on the other hand. We chose to base the BNP SSD on NRMI rather than on the more common Dirichlet Process, because it is more robust in case of misspecification of the number of clusters [START_REF] Lijoi | Controlling the reinforcement in Bayesian non-parametric mixture models[END_REF][START_REF] Barrios | Modeling with Normalized Random Measure Mixture Models[END_REF]. The BNP SSD provides several benefits for risk assessment: it is an effective and robust standard model which adapts to many datasets. Moreover, it readily provides credible intervals. While it is always possible to obtain confidence intervals for a frequentist method using bootstrap, it can be difficult to stabilise the interval for small datasets even with a large number of bootstrap samples. As such, the BNP SSD represents a safe tool to remove one of the arbitrary parametric assumptions of SSD [START_REF] Forbes | Species Sensitivity Distributions Revisited: A Critical Appraisal[END_REF]. The extended paper supporting the BNP SSD [START_REF] Kam King | Species Sensitivity Distribution revisited: a Bayesian nonparametric approach[END_REF] will include a comparison of methods on simulated data, an extension to the case of censored data and an emphasis on the potential benefits of the approach from a biological point of view.

Fig. 1 :

 1 Fig.1: The top panel represents large-size carbaryl dataset, the middle panel represents the medium-sized temephos dataset, the bottom panel represents small-sized captan dataset. Fits of the Normal (in green), KDE (in red) and BNP (in blue) models. Concentrations are log transformed. Left: Histogram and density estimates. Centre: Boxplot for the LOOs (for Normal and KDE) and the CPO (for BNP) on logarithmic scale. The horizontal line corresponds to the median. The box hinges extend to the inner quartiles. The whiskers extend to cover points up to one and a half times the inter-quartile distance away from the hinges. For both frequentist methods, the n LOOs are obtained by fitting the model n times, while an analytical expression is available for the BNP method (Equation5). Right: log HC 5 and associated confidence/credible intervals (for Normal, KDE and BNP).

Acknowledgements

The authors would like to thank two anonymous reviewers for suggesting several clarifications to the article. J. Arbel was a postdoc at Bocconi University and Collegio Carlo Alberto, Italy, when this article was submitted. J. Arbel and I. Prünster are supported by the European Research Council (ERC) through StG "N-BNP" 306406.


Rafa L Kozik
email: rafal.kozik@utp.edu.pl

Micha L Choraś
email: mchoras@itti.com.pl

Rafa L Renk
email: renk@amu.edu.pl

Witold Ho Lubowicz

A proposal of algorithm for web applications cyber attack detection

Keywords: web attacks detection, web applications firewall, machine learning, data mining

Injection attacks (e.g. XSS or SQL) are ranked at the first place in world-wide lists (e.g. MITRE and OWASP). These types of attacks can be easily obfuscated. Therefore it is difficult or even impossible to provide a reliable signature for firewalls that will detect such attacks. In this paper, we have proposed an innovative method for modelling the normal behaviour of web applications. The model is based on information obtained from HTTP requests generated by a client to a web server. We have evaluated our method on CSIC 2010 HTTP Dataset achieving satisfactory results.

Introduction

Currently, providing effective cyber-security solutions for web applications is very challenging. This happens due to the fact that the commonly used IDS (Intrusion Detection System) and IPS (Intrusion Prevention System) systems have problems in recognising new attacks (0-day exploits), since these systems are based on the signature-based approach. In such a mode, when the system does not have an attack signature in its database, such attack is not detected. Therefore, there is a need to develop more sophisticated methods that are both capable of adapting domain expert knowledge [START_REF] Choraś | Ontology Applied in Decision Support System for Critical Infrastructures Protection[END_REF] [START_REF] Choraś | Network Events Correlation for Federated Networks Protection System[END_REF] as well as emerging cyber security solutions (e.g. event correlation [START_REF] Choraś | Network Event Correlation and Semantic Reasoning for Federated Networks Protection System[END_REF] and data mining).

The list of top 10 most critical risks related to web applications security, provided by OWASP (Open Web Application Security Project [START_REF]The Ten Most Critical Web Application Security Risks[END_REF]) indicates "Injection" vulnerabilities as the major vulnerability. The Injection flaws, such as SQL, OS, and LDAP injection occur when improperly validated data containing malicious code is sent to an interpreter as part of a command or query.

According to the OWASP ranking, the second on the list are attacks related to Broken Authentication and Session Management. Incorrectly implemented authentication usually allows attackers to compromise passwords, keys, or session tokens, or to exploit other implementation flaws to assume other users identities. The XSS (Cross Site Scripting) take the third place on the OWASP list. The XSS flaws occur whenever an application takes untrusted data and sends it to a web browser without a proper validation or escaping. XSS allows attackers to execute scripts in the victim's browser which can hijack user sessions, deface web sites, or redirect the user to malicious sites.

Therefore, in order to counter those top ranked problems,in this paper we have proposed an innovative method for modelling the normal behaviour of web application. The model is based on information obtained from HTTP requests generated by client to a web server.

The remainder of this paper is structured as follows. First, we give an overview of methods for we application attacks detection. Next, the detailed method description is presented. The experiments set-ups as well as the results are presented in following section. Conclusions are given thereafter.

Overview of methods for web application attacks detection

There are several tools and methods for detecting the cyber attacks targeting web applications. Some of the frequently used tools use static code analysis approaches in order to find the vulnerabilities that may be exploited by any cyber attack. Some examples of such tools include PhpMiner II [START_REF] Shar | Predicting common web application vulnerabilities from input validation and sanitization code patterns[END_REF], STRANGER [START_REF] Yu | Stranger: An automata-based string analysis tool for PHP[END_REF], AMNESIA [START_REF] Halfond | AMNESIA: analysis and monitoring for NEutralizing SQL-injection attacks[END_REF]. However, as it is stated in [START_REF]Source Code Analysis Tools[END_REF], the difficulty relates to the fact that many kinds of security vulnerabilities are hard to find automatically (e.g. access control issues, authentication problems). Therefore, currently such tools are only able to automatically find a relatively small fraction of application security flaws.

One of the most commonly used and popular class of tools for web application cyber attacks detection adapts signature-based approach to describe (and detect) cyber attacks. Some examples include PHP-IDS [10], SCALP [START_REF]Apache Scalp Project homepage[END_REF], Snort [START_REF]Snort project homepage[END_REF]. The biggest advantage of such tools is their ability to process huge amounts of data. This is due to the fact that there are efficient algorithms that are able to check a given piece of text against a pattern (usually expressed as PCRE [START_REF]Perl-compatible regular expressions[END_REF] regular expressions) in a short time. However, the common drawback is that an expert knowledge is required to build such patters describing cyber attack. Moreover, such attacks like SQL injection are easy to obfuscate (e.g. using URL encoding). This makes the problem of providing a reliable pattern of an attack very difficult.

Proposed method overview

The proposed method overview is shown in Fig. 2. It adapts a machine-learning paradigm, therefore two distinct phases are presented on the diagram. During the learning phase, the labelled data is required in order to establish the model parameters of normal application behaviour. As mentioned before, only HTTP request headers are used for model training and these need to be labelled as either normal or anomalous.

As it is shown in Fig. 2 A, the HTTP requests need to be parsed in order to extract significant parts. In this approach, first the URL (e.g. https://host/users) is extracted and concatenated with the HTTP method name (e.g. GET, POST, PUT, etc.). Those two parameters are used to create a key (address) to entry in hashmap. During the learning phase, the hashmap is populated only with a normal (legitimate) HTTP request. This allows us to build a whitelist of resources that are usually requested by users via HTTP protocol. Whenever the HTTP request has parameters (e.g. parameter1=value1&parameter2=value2), it is encoded using the method described in section 3.1. However, the method produces vectors that are of different length. Therefore, to make it possible to learn a classifier, we transform this vector to histograms of constant length. The final feature vector is extended with information whenever a given request is on the whitelist or not.

During the testing phase (or when the algorithm operates in a production environment), the key is established using the same procedure as before (URL is concatenated with the HTTP method). If the HTTP request contains the parameters, it is encoded in order to produce the feature vector that is extended with information whenever a given request is on the whitelist or not. Such feature vector is recognized as normal or anomalous with the classifier learnt before.

Encoding HTTP parameters

In order to encode the HTTP parameters as a feature vector, we use dictionary D that maps parts of text to a set of natural numbers (see equation ( 1)).

D : word → {i : i ∈ N} (1)
The dictionary is established on a learning set using the algorithm (1) that adapts a modification of LZW compression method (Lempel-Ziv-Welch [START_REF]LZW algorithm[END_REF]). In contrast to the original LZW method we first establish the dictionary D during the learning phase and then encode the text (we do not extend the dictionary further when the method operates in the production environment). The algorithm scans through the input request parameters S for successively longer substrings until it finds one that is not in the dictionary. If a given substring in not in the dictionary, then it is added and the whole procedure is repeated until the whole learning data set is processed.

Data: Set of HTTP request parameters S Result: Dictionary D s = empty string while there is still data to be read in S do ch ← read a character;

if (s + ch) ∈ D then s ← s+ch; else D ← D ∪ (s + ch); s ← ch; end end
Algorithm 1: Algorithm for establishing dictionary D Once the dictionary D is established, the HTTP request parameters are encoded. The algorithm scans again through the input request parameters S for successively longer substrings until it finds one that is not in the dictionary. The longest substring is encoded with the natural number that indicates its position in the dictionary D. The procedure is described by algorithm [START_REF] Choraś | Network Events Correlation for Federated Networks Protection System[END_REF].

Data: Set of HTTP request parameters S, dictionary D Result: Vector V of natural numbers s = empty string while there is still data to be read in S do ch ← read a character; Moreover, the requests targeting hidden (or unavailable) resources are also considered as anomalies. Some examples from this group of anomalies include client requests for: configuration files, default files or session ID in URL (symptoms of HTTP session take over attempt). What is more, the requests not having the appropriate format (e.g. telephone number composed of letters) are also considered anomalous. As the authors of the dataset explained, such requests may not have a malicious intention, but they do not follow the normal behaviour of the web application.

V ← 0 ; if (s + ch) ∈ D then s ← s+ch; else V ← V ∪ D(
According to the authors knowledge, there is no other publicly available dataset for the web attack detection problem. The datasets like DARPA or KDD'99 are outdated and do not cover many of the current attacks.

Results

For evaluation purposes we have adapted the 10-fold cross-validation technique.

For that approach, the data obtained for learning and evaluation purposes is divided randomly into 10 parts (sets). For each part it is intended to preserve the proportions of labels (number of anomalies and normal feature vectors) in the full dataset. One part (10% of full dataset) is used for evaluation while the remaining 90% is used for training (e.g. establishing model parameters).

When the classifier is learnt the evaluation data set is used to calculate the error rates. The whole procedure is repeated 10 times, so each time different part is used for evaluation and different part of data set is used for training.

The result for all 10 runs (10-folds) are averaged to yield an overall error estimate.

In these experiments we have evaluated such classifiers as: J48, PART, Ad-aBoost, and NaiveBayes. In Tab.1, a comparison of different classifiers is presented. Moreover, we reported the effectiveness of our method using as a baseline the approach proposed by Nguyen et al. in [START_REF] Nguyen | Application of the Generic Feature Selection Measure in Detection of Web Attacks[END_REF] (the authors of CSIC 2010 HTTP Dataset).

The ROC curve for different classifiers has been presented in Fig. 3. It was generated with the WEKA tool [START_REF]ROC curve generation[END_REF], which varies the threshold on the class probability estimates. The best results have been observed for J48 tree classifier. It was possible to achieve a better detection rate while having lower false positive rate in comparison to the method proposed by Nguyen et al. in [START_REF] Nguyen | Application of the Generic Feature Selection Measure in Detection of Web Attacks[END_REF]. Additional, we have used t-test to evaluate the statistical significance of the obtained results. Test showed that the results are statistically significant at 0.95 level.

Conclusions

In this paper we have proposed an innovative method for detecting current cyber attacks targeting web applications.

We have particularly focused on solutions that are using HTTP protocol to communicate clients with the servers. We have shown that recent cyber incidents reports prove that there is an increasing number of attacks targeting these webbased applications.

The analysis show that the attacks exploiting injection vulnerabilities are still one of the most dangerous and frequently reported by institutions gathering statistics about the network incidents.

The proposed algorithm for detecting the cyber attacks targeting the web applications relies on fact that it is more effective to model normal behaviour of an application (observing the HTTP traffic) than to produce the reliable attack signature.

We have evaluated the proposed method using CSIC 2010 HTTP Dataset. The experiments have shown that the proposed method achieves satisfactory results. Moreover, we have compared our method with method proposed by CSIC Dataset authors. We report that our method that it is able to achieve the higher detection rate while having lower false positive rate.

Fig. 1 .

 1 Fig. 1. The 2013 vs. 2012 and 2011 Web Application Vulnerabilities Trends(source [5])

Fig. 2 .

 2 Fig. 2. The proposed algorithm (A indicates learning phase, while B indicates testing/classification phase)

Fig. 3 .

 3 Fig. 3. ROC curve for different classifiers.

  For the experiment the CSIC'10 dataset[15] was used. It contains several thousands of HTTP protocol requests which are organised in a form similar to Apache Access Log. The dataset was developed at the Information Security Institute of CSIC (Spanish Research National Council) and it contains the generated traffic targeted to an e-Commerce web application. For convenience the data was split into anomalous, training, and normal sets. There are over 36000 normal and 25000 anomalous requests. The anomalous requests refer to a wide range of application layer attacks, such as: SQL injection, buffer overflow, information gathering, files disclosure, CRLF injection, XSS, and parameter tampering.

	s);
	s ← ch;
	end
	end
	Algorithm 2: Algorithm for encoding text with dictionary D
	4 Experimental set-up

Table 1 .

 1 Effectiveness for CSIC 2010 HTTP Dataset.

		Detection Rate (True Positives) False Positive Rate
	Nguyen et al. [16]	93.65%	6.9%
	NaiveBayes	88.89%	6.26%
	AdaBoost	83.23%	15.45%
	PART	93.35%	2.79%
	J48	95.97%	3.54%