est destinée au dépôt

Using long time series of agricultural-derived nitrates for estimating catchment transit times


P Pernot

Z Peng

S Plessis

N Carrasco

Dissociative recombination exalts molecular growth in N 2 /CH 4 plasmas

We present a comparison of the outcomes of two models of N 2 /CH 4 VUV plasmas in two sets of conditions representative of a laboratory experiment and of Titan's ionosphere. The "FullDR" model incorporates all available experimental data on dissociative recombination od C x H y N + z ions thanks to the probabilistic tree technique recently developed in our group, whereas the "Hloss" model implements the H-loss pathways used by default in many models. We show that the FullDR model enhances globally the production of neutrals and ions, through the production of C 2 H x fragments. These species are efficient promoters of molecular growth, notably in the higher pressure conditions of the laboratory experiment.

Introduction

Complete sets of dissociative recombination branching ratios (BR) are known for only a handful of ions [1,[START_REF] Larsson | Dissociative Recombination of Molecular Ions[END_REF]. In most cases, notably for C x H y N + z ions, the distribution of H atoms amongst the fragments is not experimentally reachable. Without this information, the chemical identity of the fragments is unknown, and therefore their further reactivity. This situation has for long prevented the use of a large corpus of experimental BR data in detailed chemistry models, and many modelers have instead implemented a default "H-loss" scheme, where the parent ions looses a single H atom [START_REF] Keller | [END_REF][4][5][6][7][8]. The "H-loss" scheme is not supported by the experiments: many ions present major dissociation pathways involving heavy atoms bond breaking, for instance C 4 H + 8 [9] or HC 3 N + and other nitrile ions [10]. We have solved this problem by proposing a probabilistic representation enabling to handle incomplete BR sets [11,12]. For the first time, we have been able to take advantage of all available BR data for C x H y N + z ions in a model of the N 2 /CH 4 plasma in Titan's ionosphere [13]. This preliminary model considered a fixed neutral bath composition, and did not enable us to quantify direclly the impact of these new data on the predicted mole fractions of neutral species.

Recently, a new fully coupled ion-neutrals model has been developed and used to study the impact of different physical and chemical factors on the molecular growth in N 2 /CH 4 photochemical plasmas Table 1. Physical conditions of the two scenarii discussed in the text. P is the pressure in the reactor, T is the gas temperature and I 0 is the integrated photon flux between 50 and 200 nm. DISCO is a VUV beamline at the French SOLEIL synchrotron [20].

APSIS

"Titan" P /Pa 7 × 10 2 7 × 10 -4 T /K 300 150 Spectrum DISCO Solar I 0 /ph.m -2 .s -1 2 × 10 12 1 × 10 9 [START_REF] Peng | Experimental and Theoretical Simulations of Titan's VUV Photochemistry[END_REF][START_REF] Peng | [END_REF], with a specific focus on synchrotron-based experimental setups designed to mimic Titan ionospheric chemistry [17,18], and notably the APSIS reactor [19]. Of particular interest is the role of ion chemistry on the formation of heavy N-bearing molecules in ionospheric conditions [16].

The present article uses this model to show how dissociative recombination impacts molecular growth in these plasmas and to compare the predictions of our "FullDR" model with those of the corresponding "H-loss" model. The data needs resulting from this study are also underlined.

The model

A photochemical-transport model has been designed to represent different approximations of the APSIS reactor [19], from 0D (1-cell model) to 2D (cylindrical symmetry). In order to be able to manage the computational cost of the probabilistic tree representation of dissociative recombination BR sets [11], we use here a 1-cell geometry. The model assumes no transport and a uniform spatial distribution of all species, but takes into account the gas inflow (injection) and outflow (pumping) and the radiative transfer. Details of the model are presented elsewhere [START_REF] Peng | Experimental and Theoretical Simulations of Titan's VUV Photochemistry[END_REF][START_REF] Peng | [END_REF].

Simulations are performed in two sets of conditions representative of (1) the APSIS experiment; and (2) a crude proxy of Titan's ionosphere (noted "Titan" thereafter). In both cases a 90:10 N 2 /CH 4 mixture is considered. The two scenarii differ in several parameters, as reported in Table . 1: (i) the synchrotron DISCO flux at SOLEIL synchrotron has a rather flat wavelength distribution, in contrast with the solar spectrum mainly peaked at Lyman- [19]; (ii) APSIS experiments are led at room temperature, whereas a reference temperature of 150 K is used to represent Titan's ionospheric conditions; and (iii) a large (10 6 ) ratio between the APSIS and "Titan" pressure conditions, which is the major source of difference in the products distributions and abundances [START_REF] Peng | Experimental and Theoretical Simulations of Titan's VUV Photochemistry[END_REF][START_REF] Peng | [END_REF].

The chemical scheme incorporates the complete neutral chemistry by Hébrard et al. [21,22], with adapted photodissociation branching ratios following Gans et al. [23]. The model implements the ion chemistry database established by Carrasco et al. [24] and Plessis et al. [11], involving C, H and N. Oxygen chemistry has been excluded because it does not play a significant role in APSIS photochemistry [19].

The "FullDR" model considers 244 species (125 neutrals and 119 ions) reacting through 1708 reactions (33 photolysis, 415 bimolecular reactions, 82 termolecular reactions, 574 dissociative recombinations, 604 ion-molecule reactions). The rate parameters of these processes are all uncertain (some of them unknown). They are represented by random variables with prescribed distributions, following Hébrard et al. [21,22], Carrasco et al. [24] and Plessis et al. [11].

A derived "Hloss" model implements the "H-loss" scheme for dissociative recombination, where appropriate. In this scenario, the 448 dissociative recombination reactions coded as probabilistic trees for 58 ions are reduced to 63 reactions [START_REF] Peng | Experimental and Theoretical Simulations of Titan's VUV Photochemistry[END_REF][START_REF] Peng | [END_REF].

06003-p.2

Results and discussion

A sample of 500 realizations of the uncertain dissociative recombination rate parameters and branching ratios has been generated. All other parameters have been kept fixed. The photochemical code has been run for each realization, and the corresponding sample of stationary state mole fractions has been used to build the Figs. 1-3 discussed below.

To assess the global impact of the improved dissociative recombination description in the "FullDR" model when compared to the "Hloss" model, we evaluated for each species the ratio of the sample average mole fractions (y F ullDR /y H loss ). The empirical cumulative density function built from the values for all the species is presented in Fig. 1. In the APSIS setup at least 60% of the species have exalted mole fractions, whereas this ratio is about 30% in the "Titan" conditions. In both cases, the enhancement factors can reach very large values: more than 10 2 for 40% of the species in the APSIS conditions (15% in "Titan"). The prediction by Plessis et al. [13] of a global increased production mediated by the reactive fragments of dissociative recombination is therefore verified.

The mole fractions probability densities for selected species have been plotted for both BR models in the APSIS conditions in Fig. 2. The data for the "FullDR" case are systematically and significantly larger than for the "Hloss" scenario. A few species (C 2 H 2 , NH 3 , C 6 H 6 and HC 3 N) considered as emblematic in Titan's ionosphere are surrounded by red lines. The impact of dissociative recombination BR representation on most species is very strong (more than one order of magnitude). Flux analysis shows that the formation of C 2 H x unsaturated fragments through dissociative recombination of C 3 H + y ions is a major contribution to this production exaltation. In this model, 30% of C 2 H 2 comes from C 3 H + 7 , and C 2 H 2 is known to contribute very efficiently to molecular growth [START_REF] Peng | Experimental and Theoretical Simulations of Titan's VUV Photochemistry[END_REF][START_REF] Peng | [END_REF]26]. These pathways are absent from the "Hloss" model, although a full set of branching ratios has been measured [25]. This effect of an increased C 2 H 2 production is directly visible on the mole fractions of C 6 H 6 and HC 3 N.

One observes also larger uncertainty for the "FullDR" data, which is due to the assimilation of unknown BR sets in this model. The very large uncertainty of most N-bearing species can be traced 06003-p.3 1).

Figure 3.

Probability density plots for the stationary mole fractions of selected species in the "Titan" physical conditions (see Table 1). back to the dissociative recombination of CH 2 NH + 2 , for which rate constant and branching ratios are missing [11,27]. The enhanced NH 3 mole fraction is also linked to this process, through the NH 2 fragment.

A similar analysis has been done for "Titan" conditions in Fig. 3. The impact of the model is very strong for NH 3 , and remains significant for C 6 H 6 , but there is no visible impact on C 2 H 2 (not shown). In this setup, the contribution of dissociative recombination to the formation of C 2 H 2 and other unsaturated C 2 H x species is negligible. Conversely, the role of dissociative recombination is dominant for many of 06003-p.4

DR2013

the hydrocarbons (e.g. C 3 , C 4 , C 6 H 4 ). The large impact on the mole fraction of NH 3 comes here from the dissociative recombination of CH 3 CNH + , for which the BR set is incomplete [11,28]. For C 6 H 6 , many processes are involved with small contributions, and no dominant source has been identified yet.

Conclusion

This study evaluates directly the importance of implementing a full representation of dissociative recombination branching ratios in models of N 2 /CH 4 photochemical plasmas. A significant effect of this representation on the mole fractions of products is observed in two typical scenarii. The effect is enhanced in the experimental setup conditions at higher pressure. In this case, it has been shown that the "H-loss" model misses major pathways of C 2 H 2 production, with large impact on molecular growth promoted by this species.

The data request of Yelle et al. [27] for the branching ratios of CH 2 NH + 2 to predict correctly NH 3 abundance is confirmed, and we have shown that it has a more global impact on many N-bearing species. More generally, as the dissociation of ions into smaller neutrals, notably C 2 H x unsaturated species has a strong impact on molecular growth, the missing breakup patterns of "heavy" C x-y N y H + z (x ≥ 5) ions are urgently needed, even with unresolved H atoms.

DR2013Figure 1 .

 1 Figure 1. Impact of the dissociative recombination branching ratios representation on the mole fractions of the products of a N 2 /CH 4 photochemical plasma in two physical setups (see text). The shaded zone depicts the 95% confidence interval on the plotted values: the points outside of this zone can be considered as significantly different from zero.

Figure 2 .

 2 Figure2. Probability density plots for the stationary mole fractions of selected species in the APSIS physical conditions (see Table1).


Jean-Patrick Baudry
email: jean-patrick.baudry@upmc.fr

Gilles Celeux
email: gilles.celeux@inria.fr

J.-P Baudry

EM for mixtures Initialization requires special care

Keywords: Gaussian mixture models, EM algorithm, initialization strategies, recursive initialization, regularized likelihood, model selection criteria

Maximum likelihood through the EM algorithm is widely used to estimate the parameters in hidden structure models such as Gaussian mixture models. But the EM algorithm has well-documented drawbacks: its solution could be highly dependent from its initial position and it may fail as a result of degeneracies. We stress the practical dangers of theses limitations and how carefully they should be dealt with. Our main conclusion is that no method enables to address them satisfactory in all situations. But improvements are introduced by, first, using a penalized loglikelihood of Gaussian mixture models in a Bayesian regularization perspective and, second, choosing the best among several relevant initialisation strategies. In this perspective, we also propose new recursive initialization strategies which prove helpful. They are compared with standard initialization procedures through numerical experiments and their effects on model selection criteria are analyzed.

Introduction

The EM algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] is one of the most used algorithms in statistics (40 240 citations in Google Scholar, early February 2015). It is concerned with models with missing data structure and mixture models are certainly the favorite domain of application of the EM algorithm (see [START_REF] Mclachlan | Finite Mixture Models[END_REF]. Often the presence of missing data makes maximum likelihood (ml) inference difficult. The idea of the EM algorithm is to maximize at each iteration the conditional expectation of the complete data loglikelihood of the model at hand to derive the ml estimator of its parameter. The context is as follows: the observed data y = (y 1 , . . . , y n ) are assumed to arise from the probability distribution f (.; θ), θ being the model parameter. The complete data are x = (x 1 , . . . , x n ) with (x i = (y i , z i ), for i = 1, . . . , n), the z i s being missing. The density of a complete observation x = (y, z) is g(x; θ) and the two densities f and g are linked by the relation f (y; θ) = g(y, z; θ)dz.

Denoting L c (θ) = i g(x i ; θ) the complete likelihood, the two steps of the EM algorithm when θ r is the current parameter are

E step Compute Q(θ, θ r ) = E θ r [log L c (θ)|y]. M step θ r+1 = arg max θ Q(θ, θ r ).
The advantages and drawbacks of EM are well-documented (see for instance [START_REF] Mclachlan | The EM Algorithm and Extensions[END_REF]. Its advantages are that (i) the likelihood L(θ) = i f (y i ; θ) is increasing at each iteration, (ii) the E and M steps lead generally to closed form and simple formulas easy to program. Its main drawbacks are that (i) EM may converge painfully slowly, (ii) its solutions may be highly dependent on its initial position θ 0 . Its advantages made it the most popular algorithm to estimate missing structure models and many extensions have been proposed to answer its drawbacks (see for instance Chapters 5 and 6 in [START_REF] Mclachlan | The EM Algorithm and Extensions[END_REF].

This article is focused on the possible influence of the EM algorithm initialization for the estimation and the selection of a mixture model. To be specific, we restrict our presentation to Gaussian mixture models. But all our considerations are expected to be relevant for all latent structure models for which the EM algorithm encounters many local maxima, slow convergence situations, or degeneracies.

In a Gaussian mixture model, the observed data (y = y 1 , . . . , y n ), y i ∈ R d are assumed to arise from the density f (y; θ) = K k=1 p k φ(y; η k ), where p k ≥ 0 for k = 1, . . . , K with K k=1 p k = 1 are the mixing proportions, and where Φ(.; µ, Σ) is the density of a Gaussian distribution with mean µ and covariance matrix Σ. The vector parameter of the mixture model is denoted θ = (p 1 , . . . , p K-1 , µ 1 , . . . , µ K , Σ 1 , . . . , Σ K ). The missing data are the labels z i from which y i arises for i = 1, . . . , n. These labels are binary indicator vectors in {0, 1} K : z ik = 1 if y i arises from the kth component, otherwise z ik = 0. The EM algorithm for this model consists of E step For i = 1, . . . , n, k = 1, . . . , K compute the current conditional probabilities that y i arises from the Gaussian distribution φ(y

i ; µ r k , Σ r k ) τ r ik = p r k φ(y i ; µ r k , Σ r k ) K =1 p r φ(y i ; µ r , Σ r ) • M step p r+1 k = n i=1 τ r ik n µ r+1 k = n i=1 τ r ik y i n i=1 τ r ik Σ r+1 k = n i=1 τ r ik (y i -µ r+1 k )(y i -µ r+1 k ) n i=1 τ r ik •
Remark We gave the EM formulas for the most general multivariate Gaussian mixture model with no restrictive assumptions on the mixing proportions or the component covariance matrices. More parsimonious models could be assumed, based on the eigenvalue decomposition of the component covariance matrices (see [START_REF] Banfield | Model-based Gaussian and non-Gaussian clustering[END_REF]. Detailed formulas for EM with fourteen different models are in [START_REF] Celeux | Parsimonious Gaussian models in cluster analysis[END_REF]. Note that for some covariance matrix decompositions the M step is no more closed form and requires a numerical optimization. When using the EM algorithm, there is no guarantee that it provides the ml estimate of the model at hand since its solution depends on its initial position. Moreover, mixture ml estimation can be jeopardized with degeneracies. It can blur statistical inference for estimation and also for model selection. As a matter of fact, selecting a proper mixture model is an important and difficult problem. In particular, the choice of the number of components in a mixture model is quite influential. This choice depends on the purpose of the modeling, but it also depends on the quality of the assessed ml solutions. Actually, most model selection criteria and procedures are based on the likelihood of the consistent ml estimators. Using poor ml estimates at some place can lead to a misleading choice of mixture model.

All the traps of the EM algorithms (high dependence on starting values, slow convergence, convergence towards insensible or spurious maxima) arise more often for complex models. In particular, EM is jeopardized with mixture models with a great number of components or in high dimension settings or when many local maxima of the likelihood function are present. In such situations, this algorithm should be used with a great care to provide relevant and stable parameter estimates.

The aim of this article is to propose relevant procedures to get honest ml estimates derived with the EM algorithm. This is by the way a necessary condition to get a relevant assessing of the number of components in mixture models.

The article is organized as follows. Popular mixture model selection criteria are reviewed in Section 2. Section 3 is devoted to initialization strategies for EM to estimate a mixture model. First, standard strategies are reviewed. Then a recursive initialization procedure is proposed. The issue of degeneracies is considered in Section 5 and it is proposed to avoid them by replacing the likelihood with a regularized likelihood in a Bayesian perspective [START_REF] Fraley | Bayesian regularization for normal mixture estimation and model-based clustering[END_REF]. Section 6 is devoted to numerical experiments comparing different initialization procedures and their effect on mixture model selection.

The numerical experiments on both simulated and real data sets, highlight the performances of the compared initialization strategies to get sensible local maximizers of the penalized likelihood. A discussion section ends this article.

Mixture model selection criteria

This section is focused on the selection of the number of components in a mixture model. Thus, it is convenient to index the parameter θ with the number of components of the mixture at hand and denote it θ K . We only sketch some popular model selection criteria and refer the reader to the Chapter 6 in [START_REF] Mclachlan | Finite Mixture Models[END_REF] for an extensive presentation on assessing the order of a mixture model.

BIC criterion

If the purpose of the mixture modeling is to get a semi-parametric density estimation, the BIC criterion [START_REF] Schwarz | Estimating the dimension of a model[END_REF] has been proved to be practically adequate for choosing K [START_REF] Roeder | Practical Bayesian density estimation using mixtures of normals[END_REF][START_REF] Fraley | Model-based clustering, discriminant analysis and density estimation[END_REF]. Actually, under some regularity conditions, it is proved to asymptotically select the model minimizing the Kullback-Leibler distance to the true distribution, so the true model when it is available [START_REF] Keribin | Consistent estimation of the order of mixture models[END_REF]. This criterion has been derived within a Bayesian framework. It is an asymptotic approximation, as the number of observations goes to infinity, of the logarithm of the integrated likelihood

f K (y) = f (y; θ K )π(θ K )dθ K , π being a prior distribution for θ K . It is BIC(K) = log L( θK ) - D K 2 log n, (1)
D K being the dimension of θ K . Note that BIC is a pseudo-Bayesian criterion since it does not depend on the prior distribution π(θ K ).

ICL criterion If the purpose of the mixture modeling is model-based clustering the ICL criterion could be preferred to BIC (see [START_REF] Biernacki | Assessing a mixture model for clustering with the integrated completed likelihood[END_REF]. ICL is an asymptotic BIC-like approximation of the logarithm of the complete integrated likelihood

g K (y, z) = g(y, z; θ K )π(θ K )dθ K ,
where the missing labels z i are replaced with their most probable values ẑi for the ml estimate θK : for i = 1, . . . , n, ẑik = 1 if k = arg max τi and 0 otherwise. Here τik indicates the conditional probability of y i to arise from the kth component of the mixture with ml parameter θK . ICL is

ICL(K) = BIC(K) + ENT(K), (2)
where

ENT(K) = - n i=1 K k=1
τik log τik is the entropy of the fuzzy clustering associated to the (τ ik ) i=1,...,n;k=1,...,K .

The slope heuristic It is important to recall that BIC and ICL are asymptotic criteria which behavior can be jeopardized when the number d of descriptors is not small with respect to the number n of observations. In such cases the slope heuristics of [START_REF] Birgé | Minimal penalties for Gaussian model selection[END_REF], see also [START_REF] Baudry | Slope heuristics: overview and implementation[END_REF] can be a useful alternative. In the mixture context, it may be summarized as follows. Considering a criterion of the form

SH(K) = log L(K) + ρD K , (3)
it is assumed (i) that there exists a minimal ρ min such that any lighter penalty selects models with clearly too high complexities and such that heavier penalties select models with reasonable complexity; (ii) ρ opt = 2ρ min provides us with an optimal penalty. In [START_REF] Baudry | Slope heuristics: overview and implementation[END_REF], it is proposed to directly estimate ρ min by the slope of log L(K) with respect to D K via a robust regression method. Obviously, this procedure requires that log L(K) behaves linearly with respect to D K at least for large enough values of K. Thus, an efficient use of the slope heuristics requires to get stable maximum likelihood values for a large amount of numbers of mixture components. SH is expected to minimize the risk [START_REF] Birgé | Minimal penalties for Gaussian model selection[END_REF][START_REF] Baudry | Slope heuristics: overview and implementation[END_REF].

In order to lead to an honest selection of K, all these criteria require reasonable evaluations of the consistent ml estimates of θ K for a large amount of values of K. As a consequence, solid initialization procedures for EM are needed.

Initializing the EM algorithm

The choice of θ 0 is decisive for the EM algorithm. Several strategies have been proposed to initiate EM for estimating mixture parameters and are available in mixture softwares. Initialization strategies can be distinguished by the importance they give to randomness. Some procedures available in the mixmod software (http://www.mixmod.org) use intensively random initializations and have been proposed in [START_REF] Biernacki | Choosing starting values for the EM algorithm for getting the highest likelihood in multivariate gaussian mixture models[END_REF] (see also [START_REF] Berchtold | Optimisation of mixture models: Comparison of different strategies[END_REF].

Random This basic procedure consists of running the EM algorithm until convergence from several random positions and to keep the solution providing the largest likelihood. Small EM This procedure consists of using a large number of short runs of EM. By a short run of EM, we mean that we do not wait for convergence and that we stop the algorithm after a few iterations (say five). Then EM is run from the parameter value providing the largest likelihood from these short runs of EM. SEM The SEM algorithm is a stochastic algorithm incorporating between the E and M steps of EM a restoration of the unknown component labels z i , i = 1, . . . , n by drawing them at random from their current conditional distribution. The SEM algorithm is the same as the Monte Carlo EM algorithm with a single replication (see McLachlan and Krishnan, 2008, Chapter 6). This algorithm generates a Markov chain with a unique stationary distribution independent of its initial position. The SEM procedure consists of starting EM with the parameter value providing the largest likelihood from a long run of the SEM algorithm. CEM The CEM algorithm, see for instance [START_REF] Celeux | A classification EM algorithm for clustering and two stochastic versions[END_REF]) is a classification algorithm incorporating a classification step between the E and M steps of EM. This classification step consists of assigning each observation y i to the mixture component maximizing the current conditional probability (τ ik ) k∈{1,...,K} . This algorithm converges in a small number of iterations. The CEM procedure consists of repeating a large number of times the CEM algorithm from random initial positions and to start EM with the parameter value providing the largest likelihood from these runs of the CEM algorithm.

Moreover, it is recommended to involve several starts from the procedures Small EM and CEM. Among these procedures, Small EM is often preferred. It is the default intialization procedure in the mixmod software. But none of them has shown to outperform the other procedures. Moreover, all of them can appear to be disappointing in a large dimension setting because the domain parameter to be explored becomes very large or when the number of mixture components is large.

At the opposite, deterministic initialization procedures such as the hierarchical procedure proposed in the software Mclust do not use random starting solutions (http://www.stat.washington.edu/mclust). Such hierarchical initialization can be outperformed by the above mentioned random procedures. But, they provide more stable results which can be viewed as an advantage in a high dimensional setting or for large numbers of mixture components.

Herafter, we propose recursive intialisation procedures which can be expected to take profit of the advantages of both worlds.

Recursive Initialization

A problem which often occurs in a mixture analysis with different numbers of components is that some solutions may be suboptimal or are spurious. This problem is common with the random initialization procedures of Section 3. It could prevent from selecting a sensible number of mixture components with the criteria presented in Section 2. The recursive procedures we now present aim to avoid irrelevant parameter estimates.

Recursive algorithms

Assume that the user aims to choose a mixture component with a number of components belonging to {K min , . . . , K max }. Recursive initialization consists of splitting at random one of the K components into two components to get a (K + 1) solution.

-First the K min solution is thoroughly designed using for instance the Small EM procedure repeated a large number of times. -From K = K min , the initial position of the (K + 1)-component mixture is obtained by splitting one of the components of the K-component mixture into two components.

This strategy can be related to [START_REF] Pelleg | X-means: Extending k-means with efficient estimation of the number of clusters[END_REF], who introduce the X-means algorithm as an extension to the k-means. It enables to choose the number of clusters and to partially remedy to the initialization problem for k-means.

Three different strategies to choose the mixture component to be split will be considered.

-Random Choice: Choose at random the mixture component to be split [START_REF] Papastamoulis | On the estimation of mixtures of poisson regression models with large numbers of components[END_REF]. -Optimal Sequential Choice: Choose the mixture component to be split by optimizing a splitting criterion. -Complete Choice: Try to split all the K components and choose to split the component leading to the largest likelihood [START_REF] Baudry | Sélection de modèle pour la classification non supervisée[END_REF].

Other splitting strategies such as the strategy in [START_REF] Fraley | Incremental model-based clustering for large datasets with small clusters[END_REF] especially conceived to deal with large data sets are possible.

Splitting criteria

Possible criteria to choose the component to be split are now presented. They are derived from decompositions of the model likelihood

L(θ K ) log L(θ K ) = n i=1 log K k=1 p k φ(y i ; µ k , Σ k ).
Splitting the component with the weakest contribution to the likelihood Since

log L(θ K ) = K =1 n i=1 τ i log K k=1 p k φ(y i ; µ k , Σ k ),
it leads to split the component kL = arg min

1≤ ≤K n i=1 τ i log K k=1 p k φ(y i ; µ k , Σ k ). ( 4
)
An alternative criterion is the normalized contribution to the likelihood of a component kLn = arg min

1≤ ≤K n i=1 τ i log K k=1 p k φ(y i ; µ k , Σ k ) n i=1 τ i • (5)
Splitting the component with the weakest contribution to the complete likelihood We have

log L(θ K ) = K =1 n i=1 τ i log p φ(y i ; η ) - K =1 n i=1 τ i log τ i ,
τ i log τ i . (8)
If the normalized entropy is considered instead, it leads to split the component

kEn = arg max 1≤ ≤K - n i=1 τ i log τ i n i=1 τ i • (9)
4 Avoiding degeneracies of the EM algorithm

Deriving the maximum likelihood parameter estimate of a Gaussian mixture faces an important difficulty since the likelihood function of Gaussian mixtures with unrestricted component covariance matrices is unbounded. Thus, the EM algorithm is jeopardized by degenerated solutions. The EM algorithm can fail with some starting values, models, and numbers of components (see for instance McLachlan and Peel, 2000 Section 3.10). To get rid of this issue in our numerical experiments, we opt for a Bayesian regularization, as in [START_REF] Fraley | Bayesian regularization for normal mixture estimation and model-based clustering[END_REF] or [START_REF] Ciuperca | Penalized maximum likelihood estimator for normal mixtures[END_REF], and replace the ml estimate with the maximum a posteriori (MAP) estimate which maximizes the regularized loglikelihood log L(θ K ) + log π(θ K ), π being a prior distribution on the vector parameter θ K .

Since we make use of Bayesian inference merely to avoid degeneracies in the estimation of the component covariance matrices, we put no prior distribution on the mixing proportions and the component mean vectors. Following [START_REF] Fraley | Bayesian regularization for normal mixture estimation and model-based clustering[END_REF], we use an inverse Wishart(ν, Λ) exchangeable conjugate prior distribution for the component covariance matrices Σ k , k = 1, . . . , K. We choose the following prior hyperparameters: ν = d + 2 (as [START_REF] Fraley | Bayesian regularization for normal mixture estimation and model-based clustering[END_REF] and Λ depends on the model. Then, the MAP estimate of θ K is derived with the EM algorithm.

For the general model, Λ = σ 1/d 0 S |S| 1/d , S being the empirical covariance matrix of the data y and σ 0 a small positive number. We have |Λ| = σ 0 . The greater σ 0 the greater the regularization. Thus, this tuning parameter allows to control the regularization. It allows a weaker regularization than the hyperparameter [START_REF] Fraley | Bayesian regularization for normal mixture estimation and model-based clustering[END_REF]. The formulas of the EM algorithm remain the sames, except for the updating of the covariance matrices in the M step which becomes (at the (r + 1)th iteration) (see [START_REF] Fraley | Bayesian regularization for normal mixture estimation and model-based clustering[END_REF] for details)

Λ = S K 1/d , proposed in
Σ r+1 k = Λ + n i=1 τ r ik (y i -µ r+1 k )(y i -µ r+1 k ) ν + n i=1 τ r ik + d + 2

•

We also consider Gaussian mixture models with diagonal component covariance matrices B k , with B k = diag(B kj , j = 1, . . . , d) for k = 1, . . . , K. In this case, we use an inverseGamma( ν 2 , ζj 2 ) distribution for B kj with the prior hyperparameters ν = d + 2 and ζ j = (σ 0 ) 1/d sj (s1...s d ) 1/d , s j being the empirical variance of the variable j and σ 0 a small positive number. The updating of the B kj s at the (r + 1) iteration of the EM algorithm is

B r+1 kj = ζ j + n i=1 τ r ik (y ij -µ r+1 kj ) 2 ν + n i=1 τ r ik + 2 •
We also consider Gaussian mixture models with spherical covariance matrices λ k I(1 ≤ k ≤ K) proportional to the identity matrix. In this case, we use an inverseGamma( ν 2 , ζ 2 ) distribution for λ k with the prior hyperparameters ν = d + 2 and ζ = 2(σ 0 ) 1/d , σ 0 being a small positive number. The updating of the λ k s at the (r + 1) iteration of the EM algorithm is

λ r+1 k = ζ + n i=1 τ r ik (y i -µ r+1 k ) (y i -µ r+1 k ) ν + d n i=1 τ r ik + d + 2 •
The choice of σ 0 is obviously very important. A first idea is to make sure the chosen σ 0 does not hide the data structure. With simulated data, it is possible to monitor the evolution of a regularized EM starting from the true parameter and to exclude values of σ 0 such that the estimated parameter get too far away. For real data sets a non regularized ml parameter estimation can replace the true parameter...

Numerical experiments

The only way to assess the ability of EM initialization strategies to properly derive the ml estimate of a mixture model is to proceed to numerical experiments in practical situations. In this section, we present the results of numerical experiments on three different kinds of data sets. A six-component Gaussian mixture in R 2 is analyzed, first with "true" models, of which some contain the true distribution, and second with "false" models, of which none contains the true distribution. Then a simulated data set from a Gaussian mixture with a large number of components in R 3 is considered. Finally, numerical experiments are achieved on a real data set with a large number of observations to be clustered.

The initialization strategies considered are:

-Random. 13 runs of EM from random initializations.

-Small EM. Number of iterations of the small runs of EM: 5; number of small runs of EM: 50; repeat 10 times and keep the best. -SEM. Number of iterations in the SEM runs: 500; number of SEM runs: 8. -Mclust. Hierarchical initialization.

-KM1 Complete, L, Ln, etc. are KM1 strategies starting from the best Kcomponent solution among those obtained from the initialization strategies Random, Small EM, SEM, Mclust and KM1 (for K > K min ). -KM1 sComplete (s stands for sequential) is a KM1 strategy independent of the previous strategies (except for K = K min ).

The number of iterations of EM after the initialization step is the same for all: 1000. Moreover we consider the strategy "best without KM1" which consists of choosing the best solution after the final EM run, among those obtained from the strategies Random, Small EM, SEM, and Mclust. For all simulated numerical experiments (Sections 5.1 and 5.2), 30 data sets were simulated.

True and wrong models

We consider 600 observations from a six-component Gaussian mixture in R 2 where two components have a non diagonal covariance matrix. The mixture characteristics are given in Figure 1. At first, general Gaussian mixture models with free variance matrices were fitted to this data set. Secondly, the fitted Gaussian mixture models were restricted to diagonal covariance matrices (the family denoted by [p k λ k B k ] in [START_REF] Celeux | Parsimonious Gaussian models in cluster analysis[END_REF]. [START_REF] Celeux | Parsimonious Gaussian models in cluster analysis[END_REF]. R θ denotes the θ rotation matrix.

-2 0 2 4 6 8 10 12 -2 0 2 4 6 8 X 1 X 2 k Prop. Component Distribution 1 0.2 N 0 0 , R π 3 1 0 0 0.1 R -π 3 2 0.2 N 8 5 , R -π 3 1 0 0 0.1 R π 3 3 0.2 N 1 5 , 0.1 0 0 1 4 0.2 N 1 5 , 1 0 0 0.1 5 0.1 N 8 0 , 0.1 0 0 1 6 0.1 N 8 0 , 1 0 0 0.1
k λ k C k ] in
True models Despite the fact that it is a favorable situation for the Gaussian mixture model, the EM algorithms could be trapped when trying to maximize the non regularized likelihood as it is apparent from Figure 2. All the strategies encountered problems with the non regularized likelihood. Note that, by convention, when the EM algorithm is trapped, we set likelihood to the value obtained with K = 3 components. We do so to highlight these situations in the graph. In the following we use only regularized EM algorithms.

Note that according to Figure 3 (i) the MClust strategy becomes disappointing here for K > 4, (ii) the SEM strategy seems to outperform the other random strategies for large numbers of components and (iii) the KM1 Complete strategy provides the highest maximum for all K values.

Figure 4 shows that not surprisingly BIC provides most often the true number of mixture components except for the Mclust strategy. SH has a tendency here to overestimate the number of components. For the strategies which do not involve KM1, SH can select surprisingly high numbers of components. It is noteworthy that ICL prefers almost always the clustering solution in four clusters visible in Figure 1. This experiment illustrates clearly the fact that ICL is not aiming to recover the true number of components but to select a stable clustering solution.

- Wrong models First, for this more sensitive situation, we display the behavior of the KM1 strategies presented in Section 3.1.2 to get a sensible maximizer. It is apparent from Figure 5, that all the considered splitting criteria behave almost the same and are outperformed by the strategies KM1 Complete and KM1 sComplete. This behavior always occurs in all the situations of our numerical experiments. For simplicity, these splitting criteria will be omitted in the following. In practice, both "complete" KM1 strategies have analogous performances, but sometimes the KM1 sComplete strategy is outperformed by the KM1 Complete strategy. Since this strategy is also more natural (there is no reason not to choose the best solution among the K components), we only show the KM1 Complete strategy. Figure 6 confirms the superiority of the KM1 strategy and the good behavior of the SEM strategy. It is also important to remark that all strategies but KM1 can lead to a non increasing likelihood as a function of K. Figure 7 illustrates the ability of ICL to select a sensible number of clusters despite the fact that the family of models is not quite relevant. Not surprisingly, in this misspecified setting, BIC and SH have a tendency to overestimate the number of components (see [START_REF] Biernacki | Assessing a mixture model for clustering with the integrated completed likelihood[END_REF].

Bubbles data set

This simulated data set, plotted in Figure 8, is composed of 1000 observations in R 3 . It consists of an equiprobable mixture of three large bubble groups centered at ν 1 = (0, 0, 0), ν 2 = (6, 0, 0) and ν 3 = (0, 6, 0) respectively. Each bubble group j is simulated from a mixture of seven components according to From Figures 9 and 10, again KM1 appears to be the best strategy and the Best without KM1 solution produces a non increasing loglikehood in function of K. The SEM strategies behaves well. But, on the contrary to Mclust, the strategies Small EM and Random are amazingly suboptimal around K = 21, the true number of components. And actually, there is no sensitive difference between BIC and SH which have difficulties to recover the true number of clusters with these two strategies. Otherwise ICL does not select a stable number of clusters. terized the expression dynamics over 27 distinct stages of development during the life cycle of the fly using RNA-seq. [START_REF] Rau | Co-expression analysis of high-throughput transcriptome sequencing data with Poisson mixture models[END_REF] use a subset of these data from 12 embryonic samples that were collected at two-hour intervals for 24 hours, with one biological replicate for each time-point to illustrate a clustering method derived from a Poisson mixture. The phenotype tables and raw read counts for the 13,164 genes with at least one non-zero count among the 12 time-points were obtained from the ReCount online resource [START_REF] Frazee | ReCount: a multi-experiment resource of analysisready RNA-seq gene count datasets[END_REF]. In this paper, we consider this data set to illustrate initialization strategies of the EM algorithm for Gaussian mixture on a large data set. For this purpose, we transform the count data c with the log(c j + 1), j = 1, . . . , 12 transformation. Then we apply a PCA on the resulting data set in R 12 and apply a Gaussian mixture model to the first two coordinates of the PCA. The analyzed data set is represented in Figure 11.

For this data set, the superiority of KM1 is more marked. Moreover, the strategies without KM1 produce clearly non increasing likelihoods in function of K (Figure 12). ICL with KM1 chooses a sensible number of clusters (K = 3) depicted in Figure 11. The choice of K in a density estimation purpose seems more difficult. As a matter of fact, with the KM1 strategy, BIC chooses K = 12 components while SH chooses K = 17 components (see Table 1). Fig. 11 The Drosophila data set represented in the plane of its first principal components.

Discussion

We highlight some difficulties which can arise with the EM algorithm for maximizing the likelihood of a Gaussian mixture model. The degeneracies of the likelihood can be addressed by a Bayesian regularization. But the initialization issue remains and can be a most influencial factor.

To avoid the traps associated to degeneracies, we replace the unbounded loglikelihood with a regularized loglikelihood derived from a Bayesian regularization [START_REF] Fraley | Bayesian regularization for normal mixture estimation and model-based clustering[END_REF]. We modify their prior hyperparameter choices to get less regularized loglikelihoods, the strength of the regularization depending on a single scalar hyperparameter σ 0 . Thus the degeneracies in maximum likelihood methods are avoided. But it is important to choose carefully σ 0 to ensure stable results and meaningful estimates. [START_REF] Biernacki | Choosing starting values for the EM algorithm for getting the highest likelihood in multivariate gaussian mixture models[END_REF] proposed random Search/Run/Select initialization strategies as Small EM and SEM. They stated that no strategy outperforms the others and indicated that likelihood singularities can make the task more difficult. The present study confirms the difficulty to favor a strategy against the others. Actually, these authors remark that in their experiments the Small EM strategy is slightly better than the others. But in our numerical experiments, SEM seems to be slightly better than all the others, including Small EM... Beside random Search/Run/Select strategies, we propose recursive initialization strategies, so called KM1, which initialize the (K + 1) component mixture by the use of the K components solution. These recursive initialization strategies start from solutions which can be expected to be close to the ml estimate. Thus the EM algorithm is expected to converge faster with them than with fully random initializations. In this framework, we have defined several criteria to choose the component to be split. They all appear to be outperformed by the KM1 Complete procedure which consists of trying to split all the K components and splitting the component leading to the largest likelihood. Obviously this complete procedure is up to K times more expensive than the procedures using a criterion to split a component. But, this price remains often reasonable and ensures almost consistently better results. It is noteworthy that, on the contrary to the other initialization strategies, KM1 provides increasing loglikelihood graphs.

The influence of the initialization procedures on the choice of the complexity of a model can be important. It is not so important with the clustering criterion ICL which provides quite stable results. It is more important with the slope heuristics SH and BIC. In particular, Small EM and KM1 procedures can lead to different choices. And, in such cases, the solutions provided with the KM1 procedure can be preferred. On the other hand, it is rather difficult from these numerical experiments to compare the behavior of SH and BIC criteria, but it is beyond the scope of this article. BIC may be more able to recover the true number of components when it exists and provide stabler selections. Notice also that the KM1 strategy is helpful to produce a consistent estimation of the slope of the maximized likelihood.

Finally let us stress the importance of the choice of the parameter σ 0 to get an honest regularization. This is the topic of a current work.

  Splitting the component with the largest contribution to the mixture entropy It leads to split the component kE

Fig. 1

 1 Fig. 1 True and Wrong Models data set. 600 obs. from a six-component Gaussian mixture in R 2 (in the family denoted by [pk λ k C k ] in[START_REF] Celeux | Parsimonious Gaussian models in cluster analysis[END_REF]. R θ denotes the θ rotation matrix.

Fig. 3

 3 Fig. 3 True Model experiment. Graphs of the maximized likelihood for different strategies as a function of the number of components with the "true" mixture model [p k λ k C k ]; on the right figure, a zoom view

Fig. 4

 4 Fig. 4 True Model experiment. Frequencies of choices of the number of mixture components with the criteria SH, BIC, and ICL for the true mixture model [p k λ k C k ]

Fig. 8

 8 Fig. 8 Bubbles Data Set. 1000 obs. from a 21-component spherical Gaussian mixture in R 3

  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

Fig. 9 Fig. 10

 910 Fig. 9 Bubbles experiment. Behavior of the maximized likelihood for different strategies as a function of the number of components with the spherical mixture model [p k λ k I] for the Bubbles data set; on the right figure, a zoom view.

Fig. 12

 12 Fig. 12 Drosophila data set. Behavior of the maximized likelihood for different strategies as a function of the number of components ; on the right figure, a zoom view.

  True Model experiment. Graphs of the maximized likelihood for different strategies, including non regularized strategies (dashed), as a function of the number of components with the "true" mixture models [p k λ k C k ]

						-1800
						-2000
						-2200
					L(θK ^)	-2400
						-2600
																																Small EM
						-2800																										SEM Mclust
																																Small EM
						3000																										SEM Mclust
								1	2		3		4	5		6	7		8	9	10	11	12		13	14		15	16	17		18	19		20
																					K
	-2800 -2600 -2400 -2200 -2000 Fig. 2 * L(θK ^) -3000	*	*	*	* * * * * * * * * * * * * * * * Random Small EM SEM Mclust * Best without KM1 Best with KM1	L(θK ^)	-2150 -2100 -2050 -2000 -1950 -1900 -2200	*	*	*	*	*	*	*	*	*	*	*	*	*	* Random * * Small EM SEM Mclust Best without KM1 * * * Best with KM1
	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15	16	17	18	19	20			3	4	5	6	7	8	9	10	11	12	13		14	15	16	17	18	19	20
										K																					K

  Wrong Model experiment. Behavior of the maximized likelihood for the different KM1 strategies as a function of the number of components with the wrong diagonal mixture model [p k λ k B k ] Wrong Model experiment. Behavior of the maximized likelihood for different strategies as a function of the number of components with the wrong diagonal mixture model [p k λ k B k ]; on the right figure, a zoom view

		-2000 -1900														*		*		*		*	*			*		*		*		*		*			*		*	*
														*
	L(θK ^)	-2100							*		*																				*			Best without Km1 KM1 Complete
																																		KM1 Random
							*																											KM1 L KM1 Ln
		-2200																																KM1 Lc KM1 Lcn
																																		KM1 Ent
																																		KM1 Entn
				*																														KM1 sComplete
				3			4		5		6			7		8		9		10	11	12		13		14		15	16		17			18		19	20
																							K
	* Fig. 5 * -2800 -2600 -2400 -2200 -2000 L(θK ^) -3000	*	*	*	* * * * * * * * * * * * * * * Random Small EM SEM Mclust * Best without KM1 Best with KM1	L(θK ^)	-2200 -2100 -2000 -1900	*	*	*	*	*	*	*	*	*	*	*	*	*	* Random * * Small EM SEM Mclust Best without KM1 * * * Best with KM1
		1	2	3	4	5	6	7	8	9	10	11	12	13	14	15	16	17	18	19	20			3	4	5	6	7	8	9	10	11	12	13		14	15	16	17	18	19	20
											K																					K
	Fig. 6

Table 1

 1 The ellipses depict the three Gaussian components chosen by ICL with the KM1 strategy Drosophila data set. Number of mixture components selected by the model selection criteria SH, BIC and ICL.

	Crit.	Small SEM Random Mclust Best without KM1 Best with KM1
	SH	12	14	9	20	14	17
	BIC	12	14	12	18	14	12
	ICL	4	4	4	4	4	3

the following density distribution:

with µ 1 = (0, 0, 0), µ 2 = (0, 0, 1.5), µ 3 = (0, 1.5, 0), µ 4 = (1.5, 0, 0), µ 5 = (0, 0, -1.5), µ 6 = (0, -1.5, 0) and µ 7 = (-1.5, 0, 0). Thus the distribution of this data set is actually a 21-component Gaussian mixture. The reader is referred to [START_REF] Baudry | Sélection de modèle pour la classification non supervisée[END_REF], Chapter 5 for more details. This is a challenging data set because of the high number of components which overlap and have small sizes. A model collection (S m ) 1≤m≤Mmax of spherical Gaussian mixtures


Patrick Maillé
email: patrick.maille@telecom-bretagne.eu

Bruno Tuffin
email: bruno.tuffin@inria.fr

Users facing volume-based and flat-rate-based charging schemes at the same time

In the Internet, the data charging scheme has usually been flat rate. But more recently, especially for mobile data traffic, we have seen more diversity in the pricing offers, such as volume-based ones or cap-based ones. We propose in this paper to study the behavior of heterogeneous users facing two offers: a volume-based one and a flat-rate one. On top of that selection, we investigate 1) the relevance for an ISP to propose the two types of offers, and optimize the corresponding prices, and 2) the existence of a solution to the pricing game when the offers come from competing providers.

I. INTRODUCTION

The common pricing scheme in the Internet is the socalled flat rate pricing, where a user pays a fixed subscription fee, and can use the network as much as he wants. But there is currently a trend in telecommunications to move to usage-based pricing schemes, where the price you pay depends on your consumption pattern. The reasons invoked for those moves generally involve arguments of congestion and fairness between users. That trend has been observed with the introduction of cap-based pricing in broadband access, and the definition of the charging schemes for wireless networks. For an extensive description of flat-rate versus usage-based pricing, see [START_REF] Odlyzko | Know your limits: Considering the role of data caps and usage based billing in internet access service[END_REF], as well as [START_REF] Maillé | Telecommunication Network Economics: From Theory to Applications[END_REF] for a general presentation and discussion about Internet economics.

With the appearance of new wireless technologies (and also the deployment of fiber-to-the-home solutions), users are faced with several options, that are not all priced the same way. In this paper we investigate the co-existence of solutions based on flat-rate pricing with others based on usage-based pricing. Users being heterogeneous in their preferences, their preferred scheme may differ. Focusing on the case of two alternatives (one flat-rate and one usage-based), we analyze the impact of this diversity of users and choices in two settings: in the first one the two solutions are proposed by a monopolist provider trying to maximize revenues, while the second one considers competing providers who set their pricing parameters to attract users in order to make revenue (hence a noncooperative pricing game). Discussions on this can be found in [START_REF] Maillé | Telecommunication Network Economics: From Theory to Applications[END_REF] and the references therein, see also [START_REF] Songhurst | Charging Communication Networks[END_REF], [START_REF] Mo | Impacts of universal service regulation for broadband internet services[END_REF].

Our model is inspired by [START_REF] Cho | Pricing for mobile data services considering service evolution and change of user heterogeneity[END_REF], where the choice of which scheme to prefer is studied, but not the coexistence of both schemes. The authors of [START_REF] Cho | Pricing for mobile data services considering service evolution and change of user heterogeneity[END_REF] study the case of two schemes jointly proposed to users in [START_REF]Pricing mobile data services of different quality levels[END_REF] to analyze conditions for both schemes to be selected by users, but no revenue-maximization or competition problematics are considered. Our contributions can be summarized as follows:

• for a user with quadratic willingness-to-pay (in terms of the consumed quantity), we express analytically the choice made depending on the user-specific parameters and the parameters of both pricing solutions;

• then, assuming some distribution over user-specific parameters we numerically treat some examples of pricing decisions in the monopoly (revenue-maximization) and competition (non-cooperative game) cases.

• Our results suggest that , with our set of chosen parameters, offering the two options in the case of a monopoly increases very slightly the provider's revenue and has a limited impact on total demand. Also, when the options are offered by competitive providers, the pricing game will end up with a price war, i.e. null revenues, or very limited revenues if we introduce some data management costs in the model.

II. MODEL: PRICING SCHEMES AND USER PREFERENCES

We consider two ISPs (possibly controlled by the same entity) offering connectivity services to users. ISP 1 is implementing flat-rate pricing with access price p 1 per user (say, per month, like other values later on), after what the user can consume the amount of volume she wishes at no extra charge. ISP 2 on the other hand implements a volume-based pricing scheme, with still an access charge per user (and per month) p 2 , but also a charge per unit of volume q 2 , so that the (per month) total charge for a consumed volume v is p 2 + q 2 v.

Users are assumed heterogeneous. We index them by a (one-dimensional) type parameter θ characterizing their valuation (willingness-to-pay) function. We consider without loss of generality a total user mass of 1, and assume a density f (θ) for value θ with corresponding cdf F . Denote by V (θ, v) the amount (in monetary units) that a type-θ user is willing to pay for consuming a volume v per month: V (θ, v) is assumed nondecreasing in v, and constant after some value v max (θ) that a type θ user is interested in getting (or that the provider can serve per user; in this case it may not depend on θ) if the service were totally free.

We assume here that the networks are over-provisioned so that there is no QoS issue. This assumption is especially relevant for wired DSL or fiber access networks for instance.

III. USER BEHAVIOR

We consider in this section that price parameters p 1 , p 2 and q 2 are fixed, and study the choice for a particular user, together with her corresponding consumption level.

A user of type θ, if associated to ISP 1 (the ISP implementing flat-rate pricing), will consume a volume

v = v * 1 (θ) such that her net utility V (θ, v) -p 1 is maximized, i.e., v * 1 (θ) = argmax v≥0 (V (θ, v) -p 1 ). But since V (θ, •) is increasing, we get v * 1 (θ) = v max (θ). Similarly, if associated to ISP 2 (that implements usage- based pricing), she will consume v * 2 (θ) with v * 2 (θ) = argmax v≥0 (V (θ, v) -(p 2 + q 2 v
)) , (assuming that this argmax exists).

Overall, a type-θ user will select the option (among {ISP 1, ISP 2, no ISP}) maximizing her utility, considering that the "no ISP" choice yields a null utility. In case of equality between two options (marginal users), we assume that the user chooses the option yielding the largest V (θ, •), maximizing the generated valuation (like a company trying then to maximize its turnover).

In practice of course, p 2 < p 1 otherwise ISP 1 is always preferred to ISP 2. Assume as in our reference paper [START_REF] Cho | Pricing for mobile data services considering service evolution and change of user heterogeneity[END_REF] that valuation functions are quadratic, of the form

V (θ, v) = θv -v 2 2a 0 ≤ v ≤ θa a 2 θ 2 v ≥ aθ,
with a a constant and θ the (non-negative) user-specific parameter (her type). Such a valuation function corresponds to a demand D(p, θ) = a(θ -p) + with x + := max(0, x) (for more, see [START_REF] Cho | Pricing for mobile data services considering service evolution and change of user heterogeneity[END_REF]), and v max (θ) = aθ.

In that case v * 2 (θ) is obtained from the maximization of V (θ, v)-p 2 -q 2 v, hence from the relation ∂V ∂v (θ, v * 2 (θ))-q 2 = 0, leading to v * 2 (θ) = a(θ -q 2 ) + . In words, v * 2 (θ) is the level of demand at unit price q 2 .
The corresponding utilities for both choices are thus

V (θ, v max (θ)) -p 1 = aθ 2 2 -p 1 V (θ, v * 2 (θ)) -p 2 -q 2 v * 2 (θ) = a 2 (θ -q 2 ) 2 -p 2 .
Summarizing, a type-θ user will choose her service option according to the following threshold rule.

Proposition 1. Let θ -:= min 2p1 a , q 2 + 2p2 a and θ + := max 2p1 a , q2 2 + p1-p2 aq2 .
Then a type-θ user will i) prefer not to join any provider if θ ∈ [0, θ -);

ii) choose the volume-based provider if θ ∈ [θ -, θ + ); iii) choose the flat-rate provider if θ ∈ [θ + , ∞).
In particular, if prices are such that 2p1 a ≤ q 2 + 2p2 a then θ + = θ -and no user selects the volume-based provider.

Proof: The no-provider option is chosen when we jointly have aθ 2 2 -p 1 < 0 and a 2 ((θ -q 2 ) + ) 2 -p 2 < 0. Those two inequalities are equivalent to θ < 2p1 a and θ < q 2 + 2p2 a , respectively, giving i). We now assume that θ ≥ θ -, i.e., one provider is selected.

• If 2p1

a ≤ q 2 + 2p2 a , then the volume-based provider could be preferred to the flat-rate if a 2 ((θ

-q 2 ) + ) 2 -p 2 > aθ 2 2 -p 1 , that is, if θ < p1-p2 aq2 + q2 2 .
However, the condition

2p1 a ≤ q 2 + 2p2 a is equivalent to p 1 ≤ a 2 q 2 2 + p 2 + aq 2 2p2 a , which yields p1-p2 aq2 + q2 2 ≤ q 2 + 2p2 a . Hence θ ≥ 2p1 a implies that θ ≥ p1-p2
aq2 + q2 2 : the volume-based provider is therefore never chosen here.

• Now assume that q 2 + 2p2 a < 2p1 a . For θ ≥ θ -the volume-based provider is selected when θ < p1-p2 aq2 + q2 2 ; above that threshold, the flat-rate is preferred. Figure 1 summarizes the two situations for the repartition of users.

0 θ 2p 1 a q 2 + 2p 2 a
No provider Flat

0 θ q 2 + 2p 2 a 2p 1 a p 1 -p 2 aq 2 + q 2 2
No provider Volume Flat Fig. 1: The two situations for the repartition of users

As an illustration, we draw in Figure 2a the curves of

V (θ, v max (θ)) -p 1 and V (θ, v * 2 (θ)) -p 2 -q 2 v * 2 (θ), i.e.
, the utilities that type-θ users would get from choosing Provider 1 or 2, respectively, for the parameter values a = 2, p 2 = 0, p 1 = 0.002, q 2 = 0.02.

In this case, 2p1 a = 0.0447 > q 2 + 2p2 a = 0.02, and we indeed observe that under 0.02, no provider is chosen, the volume-based provider is the best option between 0.02 and p1-p2 aq2 + q2 2 = 0.06, and the flat-rate one should be chosen above this last value.

We also draw in Figure 2b the case when q 2 = 0.06. For those parameter values, 2p1 a = 0.0447 < q 2 + 2p2 a = 0.06, and we indeed observe that under 0.0447, no provider is chosen, and only the flat-rate one above this value.

In both cases, the curve for the volume-based provider starts at 0 because p 2 = 0, but any strictly positive value of p 2 would lead to negative utilities for low θ values. We however say that no ISP is chosen when the user is indifferent between "no ISP" and "ISP 2" here, because the consumed quantity with ISP 2 would be zero.

IV. DETERMINATION OF PRICES AT THE PROVIDER LEVEL

We now consider the provider point of view, and investigate the pricing decisions that revenue-driven providers would make. To do so we first express the revenue for each offer depending on the pricing parameters, then examine a numerical example of the expected outcome when both options are proposed by the same entity, and when they come from competing provider who may engage in a price war.

A. Expression of the revenue from the two pricing plans

Let m F (resp., m V ) be the mass of users associated to the flat-rate (resp., volume-based) charging plan. We have

m V = [θ -,θ + ) f (θ)dθ = F (θ + ) -F (θ -) m F = [θ + ,∞) f (θ)dθ = 1 -F (θ + )
and the corresponding revenues of providers are

R V = [θ -,θ + ) (p 2 + q 2 v * 2 (θ))f (θ)dθ R F = p 1 m F = p 1 (1 -F (θ + )).
Since analytical expressions for those revenues become very heavy even for simple distributions of θ, we will present a numerical study in the remainder of this section, to highlight some phenomena that can occur.

B. When the two pricing schemes come from a single provider

We assume in the subsection that the two options are proposed by the same provider. Offering more options may lead to higher revenues because it can allow the provider to segment the market by designing offers targeted at different types of users. The goal is then to find prices maximizing the total gain: max p1,p2,q2≥0

R F + R V .
We consider two different cases, with a = 2 in both cases.

1) Log-normal distribution for θ: In a first case, θ is lognormally distributed with µ = -1 and σ = 0.2. With that type of distribution, the density has a mode, and for those values, we numerically compute the maximum revenue for the provider:

max p1,p2,q2≥0 R F + R V = 8.373 × 10 -2 ,
obtained at p 1 = 0.137, p 2 = 4.43 × 10 -2 , q 2 = 9.36 × 10 -2 .

Under those values, we get from Proposition 1 that for θ < θ -= 0.304 users do not subscribe to the service, between θ -and θ + = 0.542 they choose the volume-based provider, and the flat-rate one above θ + . This results in masses m V = 0.803 and m F = 0.0264, and revenues R V = 0.0801 and R F = 0.00361. Hence most revenue comes from the volumebased offer, that is also selected by most users, the flat-rate alternative being used by the biggest users.

For comparison purposes, let us look at the provider revenue, had he offered only one option. If only flat rate was proposed, the optimal price would be p 1 = 0.1012 leading to a revenue 7.7547 × 10 -2 , while if we only had a volumebased option, the optimal prices would be p 2 = 4.55 × 10 -2 and q 2 = 9.0 × 10 -2 leading to a revenue 8.3713 × 10 -2 . The gain from proposing the two offers is minor here; we do not claim that it will systematically be the case though. The two offers nevertheless allow big users to consume as wished: those with the flat rate option consume a volume

∞ θ + v max (θ)f (θ)dθ = ∞ 0.54196 aθf (θ)dθ = 0.0309, which was only ∞ θ + v * 2 (θ)f (θ
)dθ = 0.02597 when there was only the volume-based option. Hence a 19% increase of aggregated volume for those "high-θ" users. Remark, to see the proportions, that the total volume for the volume-based option is

θ + θ -v *
2 (θ)f (θ)dθ = 0.4758. 2) Exponential distribution for θ: If on the other hand, θ is exponentially distributed with rate 1, then we numerically compute the maximum revenue max p1,p2,q2≥0

R F + R V = 7.358 × 10 -2
obtained at p 1 = 43.35, p 2 = 0, q 2 = 1.0.

Here, except for a very negligible proportion of the population, only the volume-based scheme is chose. More precisely, θ -= 1 and θ + = 22.175, which results in masses (proportions) m V = 0.3679 and m F = 2.34 × 10 -10 , and revenues R V = 0.7358 and R F = 1.015 × 10 -8 . We also observe here a negative aspect of monopolies, that is demand reduction with respect to competitive situations: the provider sets prices such that less than 37% of users finally subscribe to one plan.

Again, the gain from offering two schemes instead of one is very small: If we had a flat rate offer only, the optimal price would be p 1 = 3.9998, leading to a revenue 0.54134, while with volume-based only the optimal prices would have been p 2 = 0 and q 2 = 1.0 leading to the same revenue as in the two-scheme situation.

C. Pricing game when the two options are offered by providers in competition

We consider now that the two options are offered by two ISPs in competition, each one trying to maximize its own revenue. We then have a non-cooperative game where both providers play with their price (respectively p 1 and (p 2 , q 2 ) for the flat-rate and volume-based providers).

1) Without any cost: price war: Figure 3a displays the best responses of players when the flat-rate provider plays with p 1 and the volume-based one with q 2 (taking p 2 = 0) and F is the cdf of a log-normal distribution with µ = -1 and σ = 0.2. It can be readily seen that there is no intersection point, so that we end up with a price war: players reduce their price to attract customers, leading to zero prices.

The same result is obtained in the case where θ follows an exponential distribution with rate 1.

2) Game with management costs: Now, let us consider the case when we subtract a cost cm V and cm F to the revenues of respectively the volume-based and flat providers, corresponding to a management cost c per unit of mass of customers.

In Figure 3b, we draw the best-response curves of both ISPs when the flat-rate provider plays with p 1 and the volumebased one with q 2 (taking p 2 = 0) and the θ parameters are log-normally distributed with µ = -1 and σ = 0.2. In this A Nash equilibrium is then obtained at (q 2 = 0.018104, p 1 = 0.01319) giving revenues (taking into account the costs) R F = 0.00150139 and R V = 4.70877 × 10 -4 . Note that these revenues are very small when compared to the monopoly case for which, with this management cost, we obtain an optimal revenue R F + R V = 0.075583 at p 1 = 0.143, p 2 = 0.049 and q 2 = 0.0895. In this competition case, m V = 0.5293 and m F = 0.4707: all users subscribe to an offer (because p 2 = 0), and much more users than in the case of a monopoly choose the flat-rate option.

q 2 p 1 BR 1 (q 2 ) BR 2 (p 1 ) (a) No management cost q 2 p 1 BR 1 (q 2 ) BR 2 (p 1 ) (b) Management cost c = 0.01
The decreasing behavior of BR 2 (p 1 ) before starting increasing is intriguing. Figure 4 shows the revenue of the volume-based provider in terms of q 2 for three values of •10 -2 q 2 R V p 1 = 0.012 p 1 = 0.01319 p 1 = 0.016 Fig. 4: Revenue for different values of p 1 when when θs are log normally distributed and with a management cost c = 0.01. p 1 , and confirms that the value of q 2 optimizing R V is not monotonous in p 1 . To further understand what happens, we plot in Figure 5 the thresholds θ -(that is the same for the three values of p 1 selected) and θ + , recalling that users with θ ∈ [θ -, θ + ) select ISP 2. An explanation for this phenomenon comes from the fact that q 2 affects both θ -(positively) and θ + (negatively). Among the users with parameter close to θ -, there may be some "costly clients" with small consumption, whose impact on cost exceeds the additional revenue, but as illustrated in Figure 5, ISP 2 cannot deter them from entering the system (by increasing θ -) without losing "money-making" customers to the competitor (by reducing θ + ). But when p 1 decreases, θ + also decreases and the the savings on deterring costly clients from subscribing by increasing q 2 may now exceed the revenue from money-making customers that ISP 2 would lose: in that case ISP 2 is better off increasing q 2 . This is possible here even if the slope of θ -is very small, because the log-normal distribution has a mode, hence a small change in θ -may affect a significant proportion of users (the "costly" ones). This phenomenon justifies the introduction of a non-null fixed part p 2 in the volume-based pricing scheme, as a way for ISP 2 to avoid those costly clients.

V. CONCLUSIONS AND FUTURE WORK

This paper investigates the coexistence of flat-rate and volume-based pricing, where the choice is left to users, but the pricing parameter values are set by revenue-maximizing providers. When both options are controlled by a monopolist provider, our results indicate that the possibility of user segmentation offered by the variety of schemes does not bring large revenue improvements to the provider. On the other hand, if each option is offered by an independent provider, competition drives prices down, questioning the survivability of the providers.

Fig. 2 :

 2 Fig.2: Utility functions for a = 2, p 2 = 0, p 1 = 0.002.

Fig. 3 :

 3 Fig. 3: Best-responses when θs are log-normally distributed.

Fig. 5 :

 5 Fig. 5: Threshold values when θs are log normally distributed and with a management cost c = 0.01.

An aspect worth including in future models is the preference of users for flat-rate schemes when the expected price difference with volume-based is low. That preference, apparently contradicting the user rationality assumed in this paper, can be explained by users' aversion to risk and to having to monitor one's consumption. Other directions to investigate include the study of other types of schemes, such as the one proposed in [7] where a flat rate is applied but users can additionally purchase higher-quality service, charged on volume.