Gildas Mazo

Yaroslav Averyanov

Constraining kernel estimators in semiparametric copula mixture models

Keywords: copula, kernel, semiparametric, nonparametric, mixture model, clustering

This paper presents a novel algorithm for performing inference and/or clustering in semiparametric copula-based mixture models. The algorithm replaces the standard kernel density estimator by a weighted version that permits to take into account the constraints put on the underlying marginal densities. Lower misclassification error rates and better estimates are obtained on simulations. The pointwise consistency of the weighted kernel density estimator is established under an assumption on the rate of convergence of the sample maximum.

Introduction

In modern data science, the observations of heterogeneous clusters is not uncommon. An example is given in [START_REF] Fujimaki | Online heterogeneous mixture modeling with marginal and copula selection[END_REF] where one can observe two heterogeneous clusters of data points described by blood pressure and medical costs. The first dimension has a skewed Gaussian distribution and the second a log-normal distribution. The first cluster has negative dependency and the second positive dependency. These data cannot be captured by the standard Gaussian mixture model. The Student-t mixture model [START_REF] Lee | Finite mixtures of multivariate skew tdistributions: some recent and new results[END_REF] [START_REF] Mclachlan | Finite mixture models[END_REF] is not able to deal with heterogeneous clusters either.

Recently more flexible models have been considered. On the one hand, there are copula-based methods [START_REF] Kim | Mixture of D-vine copulas for modeling dependence[END_REF][START_REF] Kosmidis | Model-based clustering using copulas with applications[END_REF]. Copula-based methods allow for a separate analysis of the marginals and the dependence structure. They have been successfully applied in Pattern Recognition [START_REF] Roy | Pair-copula based mixture models and their application in clustering[END_REF], Machine Learning [START_REF] Rey | Copula mixture model for dependency-seeking clustering[END_REF], Knowledge Discovery and Database Management [START_REF] Fujimaki | Online heterogeneous mixture modeling with marginal and copula selection[END_REF]. Copulas allow for concatenating discrete and continuous data, too [START_REF] Marbac | Model-based clustering of Gaussian copulas for mixed data[END_REF]. In this paper, we only consider continuous data.

On the other hand, there are nonparametric methods. Nonparametric methods do not need to pick parametric families for the component distributions (i.e., the distributions of the clusters) but at the cost of assuming independence within each component [START_REF] Benaglia | An EM-like algorithm for semi-and nonparametric estimation in multivariate mixtures[END_REF][START_REF] Mallapragada | Nonparametric mixture models for clustering[END_REF]. In nonparametric mixture models, the parameters are probability density functions, which are estimated by kernel density estimators embedded in pseudo-EM algorithms [START_REF] Bordes | A stochastic EM algorithm for a semiparametric mixture model[END_REF].

In this paper, following the work in [START_REF] Mazo | A semiparametric and location-shift copula-based mixture model[END_REF], we combine both the copula framework and nonparametric estimation into a single mixture model. This permits to capture a wide spectrum of dependence structures while avoiding the choice of setting up the parametric families for the marginals. However, there is an important difference between the model of [START_REF] Mazo | A semiparametric and location-shift copula-based mixture model[END_REF] and ours. In the former, the distributions in the clusters were not allowed to vary in scale. In the latter, change in scale is possible. This additional degree of freedom induces a structural constraint on the component marginal densities of the mixture. The constraint is not satisfied by the kernel density estimator used in the algorithm in [START_REF] Mazo | A semiparametric and location-shift copula-based mixture model[END_REF]. How can we take the constraint into account? Will the inference be improved? To answer the first question, we have built a random weighted kernel density estimator and proved its pointwise consistency. To answer the second, we compared the algorithms on simulated and real data.

The rest of this paper is as follows. We present the models in Section 2. The first part reviews the paradigms under which one can build mixture models (Gaussian, copula-based, nonparametric and semiparametric) and the second part presents the model of interest in this paper. We give the learning algorithms in Section 3. Section 4 contains the definition and the consistency result for the weighted kernel density estimator. This section is written in a generic framework and therefore can be read independently. Section 5 and Section 6 contain the simulation experiments and the real data analysis, respectively. A Summary closes the paper.

2 Four kinds of mixture models 2.1 A review of paradigms for mixture models

We consider mixture models of the form f (x 1 , . . . , x d ) = K z=1 π z f z (x 1 , . . . , x d ), [START_REF] Aas | Pair-copula constructions of multiple dependence[END_REF] where π 1 , . . . , π K are the proportions of the K components (or clusters) and f 1 , . . . , f K are the corresponding densities. The choice of the structure for the component densities f z specifies the kind of mixture model.

The first kind of mixture model is as follows. One picks a multivariate parametric family for the component densities and estimate their parameters by maximum likelihood through an EM algorithm. In the majority of cases one usually picks the multivariate Gaussian family, or, perhaps, the multivariate Student-t family. Note that all coordinates of a vector of variables are distributed according to the same distribution up to their parameters. For instance, all the coordinates of a vector distributed according to a Gaussian mixture model are Gaussian. This is an homogeneity assumption. We refer to a standard textbook [START_REF] Mclachlan | Finite mixture models[END_REF] for further details. We note that, to deal with the complex, high-dimensional and noisy data of modern science, reseachers build more sophisticated models [START_REF] Anderson | Rectangular latent markov models for time-specific clustering, with an analysis of the wellbeing of nations[END_REF][START_REF] Coretto | Robust improper maximum likelihood: tuning, computation, and a comparison with other methods for robust gaussian clustering[END_REF][START_REF] Fop | Variable selection methods for model-based clustering[END_REF][START_REF] Mcnicholas | A mixture of variance-gamma factor analyzers[END_REF][START_REF] Viroli | Deep gaussian mixture models[END_REF]. However, the Gaussian distribution remains at the core of statistics and is often used as an important building block of those [START_REF] Anderson | Rectangular latent markov models for time-specific clustering, with an analysis of the wellbeing of nations[END_REF][START_REF] Coretto | Robust improper maximum likelihood: tuning, computation, and a comparison with other methods for robust gaussian clustering[END_REF][START_REF] Fop | Variable selection methods for model-based clustering[END_REF][START_REF] Viroli | Deep gaussian mixture models[END_REF]. The copula method, presented next, is an interesting alternative.

The second kind of mixture model arises when one chooses to use the copula decomposition for each of the component densities, that is, one writes

f z (x 1 , . . . , x d ) = c z (F 1,z (x 1 ), . . . , F d,z (x d )) d j=1 f j,z (x j ), (2)
where c z is the copula density corresponding to f z and f 1,z , . . . , f d,z are the marginals. Here F 1,z , . . . , F d,z are the corresponding (cumulative) distribution functions. Sklar's theorem [START_REF] Nelsen | An introduction to copulas[END_REF][START_REF] Sklar | Fonction de répartition dont les marges sont données[END_REF] states that for any distribution function F z with continuous marginals F 1,z , . . . , F d,z , there exists a function

C z : [0, 1] d → [0, 1], called the copula, such that F z (x 1 , . . . , x d ) = C z (F 1,z (x 1 ), . . . , F d,z (x d )), (3)
for any (x 1 , . . . , x d ) in the domain of definition of F z . The decomposition (2) follows from Sklar's theorem by differentiation. The copula C z encodes the dependence structure of a random vector. One easily checks that C z is the distribution function of the random vector (F 1,z (X 1 ), . . . , F d,z (X d )) if F z is the distribution function of (X 1 , . . . , X d ). Copulas are typically parametrized by considering families of the form {C z (•, . . . , •; θ z ), θ z } for some parameters θ z . An example is given in Section

5. If in (3) C z (u 1 , . . . , u d ) = u 1 • • • u d , then c z = 1 in (2).
This means that the variables are independent conditionally on belonging to the cluster z. In copula-based models, one can choose different parametric families for the marginals within the same cluster but this heterogeneity property comes at a price. Indeed, the specification of all the parametric families (there are dK marginals) can be a daunting task. Estimation of copula-based mixture models can be performed by EM or EM-like algorithms [START_REF] Kosmidis | Model-based clustering using copulas with applications[END_REF].

The third kind of mixture model is of nonparametric flavor. In nonparametric mixture models, one assumes

f (x 1 , . . . , x d ) = K z=1 π z d j=1 f j,z (x j ).
That is, conditionally on the labels (i.e. conditionally on being in a certain cluster), the variables are assumed to be independent. But, in contrast to copulabased mixture models, one does not assume parametric marginals. Nonparametric estimation can be performed with kernel density estimators embedded in EM-like algorithms [START_REF] Benaglia | An EM-like algorithm for semi-and nonparametric estimation in multivariate mixtures[END_REF]. In [START_REF] Benaglia | An EM-like algorithm for semi-and nonparametric estimation in multivariate mixtures[END_REF], marginals of the form

f j,z (x j ) = 1 σ j,z g j x j -µ j,z σ j,z , (4)
where µ j,z and σ j,z are location and scale parameters, respectively, are also considered. The case σ j,z = 1 and d = 1 was considered in [START_REF] Bordes | A stochastic EM algorithm for a semiparametric mixture model[END_REF]. This work largely inspired further work on nonparametric mixture models from the kernel density estimation viewpoint. But nonparametric maximum likelihood estimation is also possible if one assumes log-concavity of the component densities [START_REF] Hu | Maximum likelihood estimation of the mixture of log-concave densities[END_REF].

The fourth kind of mixture model combines nonparametric estimation and copula modeling [START_REF] Mazo | A semiparametric and location-shift copula-based mixture model[END_REF]. It is of the form (1), ( 2) and (4). In (2), the distribution functions F j,z are given by F j,z (x j ) = G j ((x j -µ j,z )/σ j,z ) and

G j (x j ) = xj -∞ g j (t) dt. ( 5
)
The model [START_REF] Mazo | A semiparametric and location-shift copula-based mixture model[END_REF] is a particular case where σ j,z = 1. The g j (hereafter called the generators) are estimated in a nonparametric way but the copula is entirely parametric, thus the term semiparametric used for this kind of models. Inference can be performed with essentially the same algorithms as in [START_REF] Benaglia | An EM-like algorithm for semi-and nonparametric estimation in multivariate mixtures[END_REF][START_REF] Bordes | A stochastic EM algorithm for a semiparametric mixture model[END_REF] but with an additional step for estimating the copula parameters. Algorithm 1 in Section 3 is an example of such algorithms.

The model of interest

We consider a model of the fourth kind, a so called location-scale semiparametric copula-based mixture model of the form

f (x 1 , . . . , x d ; π, µ, σ, g, θ) = K z=1 π z c z G 1 x 1 -µ 1,z σ 1,z , . . . , G d x d -µ d,z σ d,z ; θ z × d j=1 1 σ j,z g j x j -µ j,z σ j,z ,
that is, of the form (1), ( 2) and ( 4) where the generators g j , j = 1, . . . , d, satisfy

x j g j (x j ) dx j = 0 (6) and

x 2 j g j (x j ) dx j = 1. (7)
Here π = {π z }, µ = {µ j,z }, σ = {σ j,z }, g = {g j }, θ = {θ z }, j = 1, . . . , d, z = 1, . . . , K, are the parameters of the model (note that θ z can be a multivariate parameter). Note that there is no loss of generality in assuming a unit variance in [START_REF] Bordes | A stochastic EM algorithm for a semiparametric mixture model[END_REF]. Indeed, if the variance would be σ 2 j , say, then we could find a unique reparametrization (given by gj (x j ) = σ j g j (σ j x j ) and σj,z = σ j σ j,z ) so that (7) would be true. The copulas are parametrized by vectors θ z . No specific parametric families are assumed for the generators.

Estimation

Given the model of interest in Section 2.2, one needs to estimate the proportions π z , locations µ j,z , scales σ j,z , generators g j and copulas parameters θ z for z = 1, . . . , K and j = 1, . . . , d. Note that the estimates of the distribution functions G j can be computed through [START_REF] Benaglia | An EM-like algorithm for semi-and nonparametric estimation in multivariate mixtures[END_REF]. The sample is denoted by (x

(i) 1 , . . . , x (i) d ), i = 1, . . . , n.
Two learning algorithms are presented in this section. Algorithm 1, is essentially the same as that in [START_REF] Mazo | A semiparametric and location-shift copula-based mixture model[END_REF], which itself is inspired from the algorithms in [START_REF] Benaglia | An EM-like algorithm for semi-and nonparametric estimation in multivariate mixtures[END_REF][START_REF] Bordes | A stochastic EM algorithm for a semiparametric mixture model[END_REF]. Hence we do not consider that Algorithm 1 is a contribution of the paper. The contribution is Algorithm 2.

Description of the learning algorithms

Building upon the work of [START_REF] Benaglia | An EM-like algorithm for semi-and nonparametric estimation in multivariate mixtures[END_REF][START_REF] Bordes | A stochastic EM algorithm for a semiparametric mixture model[END_REF][START_REF] Mazo | A semiparametric and location-shift copula-based mixture model[END_REF], the most natural algorithm one can build is Algorithm 1. Algorithm 1 requires initial estimates π 0 z , µ 0 j,z , σ 0 j,z , g 0 j , θ 0 z and then

Algorithm 1

Given initial estimates π 0 z , µ 0 j,z , σ 0 j,z , g 0 j , θ 0 z and for t = 1, 2, . . . (until some stopping criterion has been reached), follow the steps below. 1. Compute (for i = 1, . . . , n and z = 1, . . . , K)

w t i,z = π t z c z G t 1 x (i) 1 -µ t 1,z σ t 1,z , . . . , G t d x (i) d -µ t d,z σ t d,z ; θ t z d j=1 1 σ t j,z g t j x (i) j -µ t j,z σ t j,z K z=1 π t z c z G t 1 x (i) 1 -µ t 1,z σ t 1,z , . . . , G t d x (i) d -µ t d,z σ t d,z ; θ t z d j=1 1 σ t j,z g t j x (i) j -µ t j,z σ t j,z
2. Process through the following steps (j = 1, . . . , d, z = 1, . . . , K).

(a) Update the cluster proportions

π t+1 z = 1 n n i=1 w t i,z .
(b) Update the location parameters

µ t+1 j,z = n i=1 x (i) j w t i,z n i=1 w t i,z . (c) Update the scale parameters (σ t+1 j,z ) 2 = n i=1 (x (i) j -µ t+1 j,z ) 2 w t i,z n i=1 w t i,z .
3. To update the generators, proceed through the following steps (j = 1, . . . , d).

(a) Generate a random variable Z (i) from Multi(w t i,1 , . . . , w t i,K ), (b) Define xi j = (x i j -µ t j,Z (i) )/σ t j,Z (i) . (c) Update the generators

g t+1 j (x j ) = 1 nh j n i=1 K x j - x(i) j h j (8)
4. Update the copula parameters (z = 1, . . . , K)

θ t+1 z = arg max θz i w t i,z log c z    G t+1 1 x (i) 1 -µ t+1 1z σ t+1 1,z , . . . , G t+1 d   x (i) d -µ t+1 dz σ d,z t+1   ; θ z   
produces a sequence π t z , µ t j,z , σ t j,z , g t j , θ t z , for t = 1, 2, . . . until some stopping criterion has been reached. The first step is similar to the E step of any EM algorithm. The second step is also similar to the EM algorithm for Gaussian mixture models: the parameters are updated by computing weighted means where the weights w t i,z relate the observations to their probabilities of belonging to the given clusters. The third step is similar to the computations undertaken in [START_REF] Bordes | A stochastic EM algorithm for a semiparametric mixture model[END_REF]. Given the data x (i) j and given the weights computed at the t-th iteration, one generates a random label Z i ∈ {1, . . . , K} according to a multinomial distribution Multi(w t i,1 , . . . , w t i,d ). One then standardizes the data according to these simulated labels, that is, builds a pseudo-sample

x(1) j , . . . , x(n) j and constructs a kernel density estimator on the top of it for updating the generators. The kernel density estimator can be constructed by following the guidelines as those in a standard textbook [START_REF] Silverman | Density estimation for statistics and data analysis[END_REF]. In [START_REF] Coretto | Robust improper maximum likelihood: tuning, computation, and a comparison with other methods for robust gaussian clustering[END_REF], the kernel is denoted by K and the bandwidth by h j . Thanks to a straightforward extension of Lemma 1 in [START_REF] Bordes | A stochastic EM algorithm for a semiparametric mixture model[END_REF], one has that, at each iteration t of the algorithm,

x(1) j , . . . , x(n) j is a sample from g t j and therefore the choice of the bandwidth can be based on that sample. Finally in the last step, one maximizes a pseudo-likelihood for the copula parameters. See [START_REF] Mazo | A semiparametric and location-shift copula-based mixture model[END_REF] for more details about this step. Algorithm 1 empirically has been found to perform well on simulations (see Section 5) whenever one is concerned with the estimation of the parameters for their own sake. However, when one is interested in the task of clustering instead, Algorithm 1 appears to have no greater value than a standard Gaussian mixture model. See Figure 1 and Section 5.

Interestingly, one can improve on Algorithm 1 by taking the inherent structure of the model into account. Note that in Algorithm 1 the estimator of the generators is not a generator itself. That is, ( 6) and ( 7) hold true but in general x j g t+1 j (x j ) dx j = 0 and x 2 j g t+1 j (x j ) dx j = 1 [START_REF] Embrechts | Modelling Extremal Events for Insurance and Finance[END_REF] do not. By letting the estimators g t unconstrained in spite of ( 6) and ( 7), information may be lost. To overcome this problem, we propose to base inference on Algorithm 2. Algorithm 2 takes into account the inherent constraints of the model by replacing the standard kernel density estimator (8) by a weighted version [START_REF] Fop | Variable selection methods for model-based clustering[END_REF] satisfying the constraints at each iteration of the algorithm. The proof of pointwise consistency of the weighted kernel density estimator are postponed to Section 4.

Algorithm 2 proceeds as follows. First one follows the instructions of Algorithm 1 till the construction of the pseudo-samples x(i) j . Then one solves an optimization problem for each marginal to get the weights of an adaptive kernel density estimator which, at each iteration of the algorithm, satisfies the constraints (9) (see Section 4). The optimization problem is convex and easy to solve. Consistency of the resulting estimator is studied in Section 4. Finally, once the marginals have been updated, a last step is added to estimate the copula parameters, as in Algorithm 1.

Heuristics underlying the learning algorithms

Increase of the log-likelihood

The two learning algorithms of Section 3.1 are designed to increase the loglikelihood of the data. They start as the standard EM algorithm. In particular, 1. Follow the steps 1 and 2 in Algorithm 1.

Generate the random labels

Z (i) ∼ Multi(w t i,1 , . . . , w t i,K ) and build the pseudo-sample xi j = (x i j -µ t j,Z (i) )/σ t j,Z (i) as in Algorithm 1. 3. Compute M n,j =    1 • • • 1 x(1) j • • • x(n) j [x (1)
j ] 2 • • • [x (n) j ] 2    , and b n,j =   1 0 1 -h 2 j   .

Solve the optimization problems

min p∈R n p 2 2 such that M n,j p = b n,j p ≥ 0,
and denote the solutions by pj = (p

j , . . . , p(n) j ).

5. Follow step 3 of Algorithm 1 but substitute (8) for

g t+1 j (x j ) = 1 h j n i=1 p(i) j K x j - x(i) j h j (10)
6. Follow step 4 of Algorithm 1 to update the copula parameters. the E step of the EM algorithm and Step 1 of Algorithm 1 are the same. In the M step of the EM algorithm, at each iteration t, one can always write the objective function, given by z i

w t i,z log c z G 1 x (i) 1 -µ 1z σ 1,z , . . . , G d x (i) d -µ dz σ d,z ; θ z (11) + z i w t i,z    d j=1 log 1 σ j,z g j x (i) j -µ j,z σ j,z    + z i w t i,z log π z ,
where w i,z are the weights calculated at the E step. The heuristic is to find values of π z , µ j,z , σ j,z , g j , θ z at which the objective function is high. For π z , the solution is standard and independent of the other parameters: it is the formula in Step 1 of Algorithm 1. For the other parameters, we use the following tricks.

If updates µ t+1 j,z , σ t+1 j,z , g t+1 j were available, then they could be plugged in to the first term in [START_REF] Fujimaki | Online heterogeneous mixture modeling with marginal and copula selection[END_REF], which in turn could be maximized over θ z . To get the updates µ t+1 j,z , σ t+1 j,z , g t+1 j , one tries to "maximize" the second term in [START_REF] Fujimaki | Online heterogeneous mixture modeling with marginal and copula selection[END_REF]. But what does it mean to "maximize" the second term? A method of [START_REF] Bordes | A stochastic EM algorithm for a semiparametric mixture model[END_REF] is used. By the decoupling of the marginals, we are left with d univariate problems similar to those of [START_REF] Bordes | A stochastic EM algorithm for a semiparametric mixture model[END_REF]. Their idea consists of updating the location parameters µ j,z (the scale parameters σ j,z were assumed to be one) by proceeding as if the generators g j were Gaussian, leading to the formula in Step 2 (b) in Algorithm 1. Once updates µ t+1 j,z have been obtained, they go on by proposing an elegant way of updating the generators g j : this is Step 3 of Algorithm 1. Note that step 2 (c) of Algorithm 1 simply incorporates the estimation of the scale parameters σ j,z in a straightforward way. The novelty of Algorithm 2 is to replace the standard kernel density estimator of Algorithm 1 by a new one that satisfies identifiability constraints and is consistent, as shown in Section 4.

Let us come back to the "as if the g j were Gaussian" argument mentioned above. Although the original authors [START_REF] Bordes | A stochastic EM algorithm for a semiparametric mixture model[END_REF] did not mention the following argument, this "Gaussian trick" can be supported for densities symmetric about zero. Indeed, maximization over µ j,z of the second term in [START_REF] Fujimaki | Online heterogeneous mixture modeling with marginal and copula selection[END_REF] is tantamount to solving i,z

g j (x (i) j -µ) g j (x (i) j -µ) w t i,z = 0.
Now, if g j is symmetric and hence an even function, then g j /g j is an odd function and odd functions are quite close to the identity function around zero because the terms of even order in Taylor expansions are zero themselves exactly. Thus, up to quadratic approximation, and assuming that the data observations are tightly concentrated around their mean value within each cluster, the formula in Step 2 (b) is approximately correct.

Initialization of the algorithms

To get initial values π 0 z , µ 0 j,z , σ 0 j,z , g 0 j , θ 0 z , one can proceed as follows. A k-means algorithm is run and the returned centers provide values for the location parameters µ 0 j,z . The returned partition of the data is used to estimate the remaining parameters. The proportions of the returned clusters provide values for the proportion parameters π 0 z . The standard errors of the clusters provide values for the scale parameters σ 0 j,z . The estimation of the generators g 0 j is based on the shuffle of all the K standarized samples. That is, one builds the sample {(x (i) j -µ 0 j,z )/σ 0 j,z } i,z and performs kernel density estimation on it to get an estimate for g 0 j . Finally, copula parameters θ 0 z are estimated by standard methods [START_REF] Genest | Everything you always wanted to know about copula modeling but were afraid to ask[END_REF][START_REF] Genest | A semiparametric estimation procedure of dependence parameters in multivariate families of distributions[END_REF], cluster by cluster. Note that, as with the standard EM algorithm, Algorithm 1 and Algorithm 2 may depend on the initialization. It is therefore advisable to test as many starting points as possible.

Check of convergence

The convergence of Algorithm 1 and Algorithm 2 is checked visually. Since the likelihood contains all the information about the parameters, we check that the log-likelihood has increased until stabilization, in average. "In average" means that one must take into account the inherent randomness of the algorithms, which is why the check is done visually. Practically, we set a large number of iterations, let the algorithms run, and inspect the plot afterwards. If the loglikelihood has not stabilized around a mean value, we increase the number of iterations. Examples are given in Section 5 and Section 6.

Kernel density estimation under moment constraints

We consider the problem of estimating the common density g of independent random variables X 1 , . . . , X n . We assume that g verifies the regularity conditions in Assumption 1

Assumption 1. The density g is continuous on R, symmetric about zero and obeys

x 2 g(x) dx = 1 = x 4 g(x) dx < ∞.
Note that the assumed symmetry implies

xg(x)dx = 0.
Continuity is a standard assumption to ensure pointwise consistency of the standard kernel density estimator [START_REF] Parzen | On estimation of a probability density function and mode[END_REF] and the Nadaraya-Watson estimator [START_REF] Watson | Smooth regression analysis[END_REF].

The condition on the moment of second order stems from the structure of the model in Section 2.2. The moment of fourth order must have a different value than that of the moment of second order to ensure the convergence of a certain quantity (see the proof of Theorem 1 for details). We view this rather as a technical condition. For instance if g were the Gaussian density, its variance would have to be not equal to 1/3. As explained in Section 3, our aim is to construct an estimator ĝ that obeys xĝ(x) dx = 0, and

x 2 ĝ(x) dx = 1. ( 12
)
We define the estimator

ĝ(x) = n i=1 pn,i K hn (X i -x) (13)
where K hn (y) = K(y/h n )/h n is a kernel depending on a positive sequence h n and where pn = (p n,1 , . . . , pn,n ) (throughout stands for the transpose operation) is the unique solution of the random optimization problem

min p∈Rn p 2 2 (14) such that M n p = b n p ≥ 0, (15)
where p = (p 1 , . . . , p n ) and

M n =   1 • • • 1 X 1 • • • X n X 2 1 • • • X 2 n   , and b n =   1 0 1 -h 2 n   .
Each pn,i is a function of the random sample. For each realization of the sample, the optimization problem ( 14) is convex and hence admits a unique solution which is denoted by pn . The constraint [START_REF] Grønneberg | The copula information criteria[END_REF] ensures that ĝ satisfies [START_REF] Gallopin | Classification et inférence de réseaux pour les donneées RNAseq[END_REF]. Indeed, elementary calculations show that (12) holds if and only if

n i=1 pn,i X i = 0 and n i=1 pn,i X 2 i = 1 -h 2 n ,
respectively. The constraints i pn,i = 1 and pn,i ≥ 0, i = 1, . . . , n, must always hold to ensure that ĝ is a density.

As soon as n > 3 the system M n p = b has infinitely many solutions and hence there are infinitely many estimators that satisfy [START_REF] Gallopin | Classification et inférence de réseaux pour les donneées RNAseq[END_REF]. We chose to pick the closest one to the standard kernel density estimator. The standard kernel density estimator is an estimator of the form [START_REF] Genest | Everything you always wanted to know about copula modeling but were afraid to ask[END_REF] where pn,i = 1/n, and the solution of min (p1,...,pn)

E ( n i=1 p i K hn (X i -x) -g(x)) 2 dx.
In our case, we cannot set pn,i = 1/n because the constraint [START_REF] Grønneberg | The copula information criteria[END_REF] would not be satisfied. But we can project (1/n, . . . , 1/n) onto the feasible space given in [START_REF] Grønneberg | The copula information criteria[END_REF], which amounts to solve the optimization problem [START_REF] Genest | A semiparametric estimation procedure of dependence parameters in multivariate families of distributions[END_REF] because minimizing p 2 is the same as minimizing pe 2 , where e = (1, . . . , 1) . Thus, the minimization of p 2 is a heuristic justified by an analogy. Moreover, even though one can imagine other criteria [START_REF] Hall | Reducing bias in curve estimation by use of weights[END_REF] for choosing p, the choice of the euclidean norm is the easiest from a theoretical and computational point of view.

Having defined the estimator in [START_REF] Genest | Everything you always wanted to know about copula modeling but were afraid to ask[END_REF], it is natural to require at least pointwise consistency. The issue resides in the constraint p ≥ 0. Without such a constraint, Lemma 1 states that the solution of the optimization problem is explicit and yields a consistent estimate. In the presence of the constraint, Theorem 1 states that consistency can be achieved under a condition on the tail of the underlying density.

Theorem 1. Suppose Assumption 1 holds. If h n → 0, nh n → ∞ and there exist constants

a n > 0, b n ∈ R such that n -1/4 a n → 0, h n a n → 0, n -1/4 b n → 0, h n b n → 0 and a -1 n (max{X 1 , . . . , X n } -b n ) (16)
converges in distribution, then the estimator ( 13) is pointwise consistent.

The conditions h n → 0 and nh n → ∞ are necessary to ensure pointwise convergence of the standard kernel density estimator [START_REF] Silverman | Density estimation for statistics and data analysis[END_REF]. The condition ( 16) is standard in extreme value theory [START_REF] Resnick | Extreme values, regular variation and point processes[END_REF]. The conditions n -1/4 a n → 0 and n -1/4 b n → 0 state that the rate at which the sample maxima grows to infinity must not be too fast. The conditions h n a n → 0 and h n b n → 0 state that the rate at which the sample maxima grows to infinity must be smaller than the rate at which the bandwidth h n vanishes. If h n is the optimal bandwidth, that is if h n ∝ n -1/5 , then the conditions n -1/4 a n → 0 and n -1/4 b n → 0 are automatically satisfied. Example 1 and Example 3 give distributions which satisfy these conditions. Example 2 is a counter-example. Example 1 and Example 2 are drawn from [START_REF] Embrechts | Modelling Extremal Events for Insurance and Finance[END_REF], p. 153-157. The computation of the normalizing constants in Example 3 is given in the Appendix.

Example 1. Let h n ∝ n -1/5 . The Gaussian distribution (2π) -1/2 exp(-x 2 /2), x ∈ R, satisfies the conditions in Theorem 1 with

a n = (2 log n) -1/2 , b n = 2 log n - log(4π) + log log n 2(2 log n) 1/2
Example 2 (Counter-example). The Cauchy distribution g(x) = [π(1 + x 2 )] -1 , x ∈ R, does not satisfy the conditions in Theorem 1. Indeed, in addition to have infinite variance, the normalization constants are given by a n = n/π and b n = 0. The sequence (a n ) does not verifies n -1/4 a n → 0.

Computer experiments

In this section, we wish to compare Algorithm 1 (hereafter called cKDE for convenience) and Algorithm 2 (fKDE) in terms of the quality of the obtained estimates. The standard Gaussian Mixture Model (GMM) was also implemented as a benchmark. We generated 500 datasets of size n = 300, 500, 700, 900 according to the following data generating process. The number of clusters was set to K = 3 and their proportion parameters were all set of equal value. The Frank family of bivariate copulas, given by

C θz (u, v) = -1 θz log 1 + (e -θz u -1)(e -θz v -1) (e -θz -1) , if θ z ∈ (-∞, ∞) \ {0}, uv, if θ z = 0.
was chosen for all of the three copulas. The parameters were θ 1 = -3.45, θ 2 = 3.45 and θ 3 = 0, corresponding to negative, positive and null dependence levels, respectively. The generators for the marginals along the first, resp. second, axis (g 1 , resp. g 2 ), were a normal, resp. a Laplace, distribution with zero mean and unit variance. The three clusters had means (µ 1,1 = -3, µ 2,1 = 0), (µ 1,2 = 0, µ 2,2 = 3) and (µ 1,3 = 3, µ 2,3 = 0) respectively. The scale parameters along the first, resp. second, axis were set to σ 1,1 = 2, σ 1,2 = 0.7 and σ 1,3 = 1.4, resp. σ 2,1 = 0.7, σ 2,2 = 1.4 and σ 2,3 = 2.8. All the three algorithms were run with 100 iterations and initialized according to Section 3.1. The kernel and the bandwidth selection method used for building the kernel density estimators were the Gaussian kernel and the method given by (3.30) in p. 47 of [START_REF] Silverman | Density estimation for statistics and data analysis[END_REF]. That is, the bandwidth is 1.06An -1/5 , where A is the minimum between the standard deviation of the data and, the interquartile range divided by 1. [START_REF] Resnick | Extreme values, regular variation and point processes[END_REF].

In order to compare the algorithms, we computed the mean absolute errors, that is, the differences in absolute value between the true parameters and the estimates. These were averaged over the clusters and the coordinates (if any). For the generators, the L 1 norm was used instead. Only the errors for the location, scale and proportion parameters were computed for GMM. The misclassification rate was computed, too. All these error measures can be computed at each iteration of the algorithms and averaged over the replications. The results are shown in Fig. 1.

From a clustering point of view, the three learning algorithms can be compared on the basis of the misclassification error rate in Fig. 1 (f). Both the semiparametric algorithms perform better than GMM, especially for large sample sizes, where nonparametric modeling can express its potential (see Fig. 2). When the sample size is not so large, cKDE performs much better than fKDE. This indicates that the new kernel density estimator implemented in cKDE was a good idea. For larger sample sizes, the difference diminishes.

With regard to estimation, cKDE performs better than fKDE but the difference is approximately constant across the sample sizes. Both the semiparametric algorithms perform better than GMM for the proportion (Fig. 1 (b)), location (Fig. 1 (c)) and scale (Fig. 1 (d)) parameters in a way that is similar to the missclassification error. That is, the gain is much more important as the sample size gets larger.

The stability of the algorithms was checked by plotting the log-likelihood and the pointwise averaged error trajectories. Specifically, we focused on the case n = 300 and inspected the behaviors of the error trajectories across the iterations. These are displayed in Fig. 3. In Fig. 3 (f), we see that, in average, the log-likelihood increases and stabilizes after 20-30 iterations. This is in agreement with the heuristics discussed in Section 3.2. The log-likelihood of cKDE is higher than that of fKDE, which is in agreement with the results found in Figure 1. In the other panels, note that the method with the lowest trajectory cannot formally be claimed the best because the point at which would converge the algorithms is unknown. Of course if the sample size is large enough then the true parameter would be close to this point. Thus the error trajectories (except the log-likelihood) are here only to inspect convergence of the algorithms.

In view of the panels (a), (c) and (d) of Fig. 3, convergence seems to have been reached after 30 iterations. For the proportion parametres in panel (b), cKDE seems to have reached convergence while it is less clear for fKDE and even GMM. The density generators in panel (e) exhibit a high variability.

To get an idea of the variability of individual trajectories, Fig. 4 trajectories for the first replication of the experiment under focus. Finally, we repeated the computer experiment described at the beginning of this section with two modifications, one at a time. First, in five successive experiments, the proportion parameters π were set to (1/4, 3/8, 3/8), (1/5, 2/5, 2/5), (1/6, 5/12, 5/12), (1/7, 6/14, 6/14), (1/8, 7/16, 7/16). Second, a fourth component with location parameter (µ 1,4 , µ 2,4 ) given by (0, 8), (0, 6), (0, 4), (0, 2), (0, 0) with independent marginals and σ 1,4 = σ 2,4 = 1 was added. The cluster proportions were π = (1/4, 1/4, 1/4, 1/4). Note that, since the centers of the other clusters are (-3, 0), (0, 3) and (3, 0), the added fourth component gets closer to the other clusters.

The computed missclassification errors averaged over 100 replications are shown in Table 1. We see that the missclassification error seems to remain stable (it exhibits only a very slight increase) as the proportion of observations in the first cluster decreases from π 1 = 1/4 to 1/6. Thus, there is some robustness in the clustering in spite of nonparametric estimation.

In general, clustering methods typically perform worst as the amount of separation between clusters decreases. This is also illustrated in Table 1, where the error goes from 7-8% to 21-25%. Compared to GMM, the increase of the error is not worst for semiparametric copula models. 6 Illustrations on real data

π 1 = 1/4 1/

Illustration on the Iris data

We performed two data analyses of the well-known iris data available in the R software. In the first analysis, we took the first (sepal length) and the third (petal length) variables and ran Algorithm 2 for four families of copulas, namely, the Frank, Clayton, Gaussian copulas and the copula representing independence between the variables. The convergence of the algorithms was checked by inspection of the log-likelihood values, which stabilized after 50 iterations around -250, -245, -246, -288, for the four copula families, respectively. The missclassification errors are given in Table 2. The lowest errors correspond to the highest likelihoods. In terms of the estimated densities, the difference between the copulas is depicted in Figure 5. The Frank and Gaussian copulas seem to best fit the data, in agreement with the results in Table 2.

Frank Clayton Gaussian independence log-likelihood -250 -245 -246 -288 missclassification error (%) 16 10 8 14

Table 2: Missclassification errors and estimated log-likelihood for the different copula families

In the second analysis, we ran Algorithm 2 with the complete data, that is, we fitted a semiparametric copula-based mixture model of dimension d = 4. Only the Gaussian copula was tested because the other copulas do not generalize easily in higher dimensions. We let the algorithm run 50 iterations and observed that the log-likelihood seem to have stabilized, see Figure 6(a). The isocontours, depicted in Figure 6(b) do not differ much from those in Figure 5(c).

Illustration on RNA-seq data

The use of high-throughput sequencing technologies to sequence ribonucleic acid content results in the production of RNA-seq data. From a statistical point of view, the observations are (realizations of) random variables Y i,j , i = 1, . . . , n, j = 1, . . . , d, each of which is a measure of the digital gene expression (DGE) of the biological entity i (e.g., a gene) for the experimental condition j. For instance, Y i,j may be the number of reads of the ith gene for the jth condition aligned to a reference genome sequence. One question of interest deals with the clustering of DGE profiles [START_REF] Rau | Co-expression analysis of high-throughput transcriptome sequencing data with poisson mixture models[END_REF]. For instance, one may want to discover groups of co-expressed genes. In recent years several clustering methods have been proposed. Poisson mixture models can be applied but they need to assume that, within a cluster, the DGE measures are independent, a very strong assumption. More precisely, they are of the form [START_REF] Rau | Co-expression analysis of high-throughput transcriptome sequencing data with poisson mixture models[END_REF] f (y; ψ) = n i=1 K z=1 π z f z (y i ; χ iz ) where f z (y i ; χ iz ) = d j=1 rj l=1 P(y ijl ; µ ijlz ) and χ iz = {µ ijlz } jl , ψ = {π z , χ iz } i,z . Here r j = 1, j = 1, . . . , d and P denotes the Poisson density. Another approach consists of applying a transformation Y i,j → Ỹi,j , i = 1, . . . , n, j = 1, . . . , d, so that the transformed data, or pseudo data, are more appropriate for Gaussian mixture models [START_REF] Rau | Transformation and model choice for rna-seq co-expression analysis[END_REF]. One such transformation [START_REF] Gallopin | Classification et inférence de réseaux pour les donneées RNAseq[END_REF] is given by Ỹi,j = log Y i,j /N j + 1

m i + 1 ,
where

N j = n i=1 Y i,j /10 6 and m i = d -1 d j=1 N -1 j Y i,j
. This approach essentially amounts to assuming that the data are Gaussian on a log-scale. The semiparametric copula-based mixture models permit to relax this assumption.

In this section, we compare the Poisson mixture model of [START_REF] Rau | Co-expression analysis of high-throughput transcriptome sequencing data with poisson mixture models[END_REF], the Gaussian mixture model and the semiparametric copula-based mixture models with Gaussian and Frank copulas. The data are high-throughput transcriptome RNA-seq data [START_REF] Sultan | A global view of gene activity and alternative splicing by deep sequencing of the human transcriptome[END_REF] downloaded from the companion R package HTSCluster of [START_REF] Rau | Co-expression analysis of high-throughput transcriptome sequencing data with poisson mixture models[END_REF]. We removed the biological replicates so that d = 2. Estimation in the semiparametric copula-based models was performed with Algorithm 2. Estimation in the Poisson mixture model was performed with the function PoisMixClus of the package HTSCluster. All the algorithms were run with a fixed number of clusters, set to K = 10, corresponding to the number of clusters selected by the integrated completed likelihood criterion in the analysis performed in [START_REF] Rau | Co-expression analysis of RNA-seq data with the HTSCluster package[END_REF].

In order the compare the models, we reproduced Fig. 2 of [START_REF] Rau | Co-expression analysis of high-throughput transcriptome sequencing data with poisson mixture models[END_REF]. The bar heights in Fig. 7 stand for the quantities

n i=1 ŵi,z Y i,j n i=1 ŵi,z d j=1 Y i,j
, each of which, according to [START_REF] Rau | Co-expression analysis of high-throughput transcriptome sequencing data with poisson mixture models[END_REF], can be interpreted as the proportion of reads that are attributed to condition j in cluster z. The quantities ŵi,z are estimates of the probability that the i-th observation belongs to the z-th cluster, estimate of which depends on the fitted model (Poisson, GMM, or semiparametric copulabased). Bar widths are proportional to πz , the estimated cluster proportions. Each bar represents a cluster and each color represents a mean normalized expression profile, the value of which is given by the bar length of a given color. In Figure 7, the results for the Poisson model, the only one which does not take into account the dependence structure within the clusters, differ from all the other models. We note that the copula-based semiparametric models are both similar (compared to the Poisson model) and different from GMM. We take this as an encouragement for copula-based semiparametric models: there are not absurd since similar to GMM; there are potentially of practical interest since they differ from GMM.

Open problems and challenges

Despite the good results of the learning algorithms in Section 5 and Section 6, there are open problems that need to be addressed to make these algorithms fully applicable in practice.

High dimensions

Are Algorithm 1 and Algorithm 2 applicable to a high-dimensional setting? In principle, the algorithms are written for any dimension d. In practice, however, two issues make the problem challenging. First, as the dimension increases, the number of flexible copula families drops rapidly. Still, there are certain families such as Gaussian copulas, Vines copulas [START_REF] Aas | Pair-copula constructions of multiple dependence[END_REF] or factor copulas [START_REF] Krupskii | Factor copula models for multivariate data[END_REF] that might be appropriate. But then one must be able to solve the optimization problem in Step 4 of the algorithms.

Number of clusters

An important problem is that of the choice of the number of clusters. To this end, many criteria write as the observed log-likelihood minus a penalization term [START_REF] Akaike | Information theory and an extension of the maximum likelihood principle[END_REF][START_REF] Biernacki | Assessing a mixture model for clustering with the integrated completed likelihood[END_REF][START_REF] Grønneberg | The copula information criteria[END_REF][START_REF] Schwarz | Estimating the dimension of a model[END_REF]. But it is unclear what the penalization term should be in semiparametric copula-based mixture models. That said, let us notice that the observed log-likelihoods between different semiparametric copula-based mixture models with the same number of parameters may be compared to at least select the appropriate copula family.

Identifiability

As of today, identifiability of semiparametric copula-based mixture models has not been proved. Proving identifiability in general is known to be a difficult problem. In fact, even purely parametric copula-based mixture models, such as those in [START_REF] Kosmidis | Model-based clustering using copulas with applications[END_REF], have not been proved to be identifiable. Identifiability of elliptical mixtures, hence including Student-t mixtures, was not proved until 2006 [START_REF] Holzmann | Identifiability of finite mixtures of elliptical distributions[END_REF]. Identifiability, in a weaker sense, of nonparametric mixtures with independence components under certain assumptions was proved in 2009 [START_REF] Allman | Identifiability of parameters in latent structure models with many observed variables[END_REF]. Identifiability in the weak sense means that non-identifiable parameters in the strong sense belong to a subset of Lebesgue measure zero.

From a statistical perspective, that is, from an estimation point of view, identifiability is an important problem. In this respect, the simulations in Section 5 are reassuring: the true parameters with which the data were simulated could be recovered.

From a learning perspective, that is, from a clustering point of view, identifiability is less important. For example, it is well known that neural networks are not identifiable and still have been enjoying success throughout the sciences.

Convergence

Another important problem is that of checking the convergence of the algorithms given in Section 3.1. The current method is to check visually that the log-likelihood has increased and stabilized, taking into account the inherent stochasticity. The problem is that the standard criterion in EM algorithms, of the form |∇ φ i log f (X

(i) 1 , . . . , X (i)
d ; φ)| < , where ∇ φ denotes the gradient operator with respect to the parameters φ, is not applicable because of the inherent randomness of the algorithms. To address this issue, a smooth summary of the log-likelihood could be considered, as, for instance, its least concave majorant.

Summary

We proposed a novel algorithm which permitted to improve the inference in semiparametric copula-based mixture models in which the marginals have a location-scale structure. We did this by replacing the standard kernel density estimator by a weighted one in order to satisfy the inherent constraints of the model. Pointwise consistency of the estimator was proved under mild assumptions. An application to RNA-seq data and a benchmark dataset (the iris data) confirmed the ability of the models to fit real data.

Research on copula-based (and hence genuinely multivariate) semiparametric models has started only recently, and, therefore, many challenges still remain. A list of important open problems is given in Section 7. Among them stands the identifiability problem and the convergence of the algorithms. In fact, with regard to the last point, even for simpler algorithms such as those in [START_REF] Benaglia | An EM-like algorithm for semi-and nonparametric estimation in multivariate mixtures[END_REF][START_REF] Bordes | A stochastic EM algorithm for a semiparametric mixture model[END_REF][START_REF] Mazo | A semiparametric and location-shift copula-based mixture model[END_REF], the convergence properties are still unknown, even though a first step has been achieved in [START_REF] Levine | Maximum smoothed likelihood for multivariate mixtures[END_REF]. Addressing these problems should open the gate for designing sound convergence check methods and performing model selection (including selection of the correct number of clusters) through pseudo-AIC criteria.

A Appendix

(I -H n )X 2 X 2 (I -H n )X 2 (X 2 pn -1 + h 2 ) (19) pn = X 2 e -XX n(X 2 -X 2 ) ( 20
)
where H n = M n ( M n M n ) -1 M n is the projection matrix on the space spanned by e, X = (X 1 , . . . , X n ), X = n -1 n i X i , and X 2 = n -1 n i X 2 i . Moreover, the estimator [START_REF] Genest | Everything you always wanted to know about copula modeling but were afraid to ask[END_REF] with pn as in ( 17)-( 18) is pointwise consistent.

Proof of Lemma 1. Since the distribution of X i has no atom at zero, one has P (∀y ∈ R 3 , M n y = 0 or y = 0) = 1, The bound does not depend on i and vanishes asymptotically in probability by assumption on the sequences a n and b n . This is enough to conclude that (i) holds with probability tending to one.

We finally show (ii). It is convenient to introduce Lemma 3 the proof of which is deferred to the end of this Section. ), E n tends to a nonzero constant in probability, P (D n X i + E n > 0, i = 1, . . . , n) → 1 and

pn,i -pn,i pn,i = A n X 2 i + B n X i + C n D n X i + E n (22)
In view of 22, one has pn,i -pn,i pn,i

≤ |A n | ∨ n i=1 X 2 i + |B n | ∨ n i=1 X i + |C n | E n -∨ n i=1 -D n X i (23)
(we used the fact that min{y 1 , . . . , y n } = -max{-y 1 , . . . , -y n } for the denominator). By assumption and by symmetry, both ∨ n i=1 X i and ∨ n i=1 -X i are O P (a n ) + b n and by assumption on v n ,

v -2 n ∨ n i=1 X 2 i = max(v -1 n ∨ n i=1 X i , v -1 n ∨ n i=1 -X i ) 2 P → 0.
Hence the numerator in ( 23) is o P (1). The denominator equals E n + D n ∨ n i=1 X i if D n < 0 and equals E n -D n ∨ n i=1 -X i if D n > 0. Either way, the denominator tends to a constant in probability and pn,i -pn,i pn,i ≤

max |A n | ∨ n i=1 X 2 i + |B n | ∨ n i=1 X i + |C n | E n + D n ∨ n i=1 X i , |A n | ∨ n i=1 X 2 i + |B n | ∨ n i=1 X i + |C n | E n -D n ∨ n i=1 -X i .
This upper bound does not depend on i and vanishes asymptotically in probability. This proves (ii). It only remains to prove Lemma 3.

Proof of Lemma 3. Let δ i,j = 1 whenever i = j and δ i,j = 0 whenever i = j. Let H i,j denote the element at the i-th row and j-th column of H n . We have pn,i -pn,i pn,i =

- n j=1 (δ i,j -H i,j )X 2 j X 2 pn-1+h 2 X 2 (I-Hn)X 2 X 2 -XXi n(X 2 -X 2 )
.

Standard calculations yield n j=1

(δ i,j -H i,j )X 2 j = (X 2 -X 2 )X 2 i + (XX 2 -X 3 )X i + XX 3 -X 2 2 X 2 -X 2
and hence we can rewrite pn,i -pn,i pn

,i = [-(X 2 -X 2 )X 2 i -(XX 2 -X 3 )X i -XX 3 + X 2 2 ][X 2 pn -1 + h 2 ] [X 2 -XX i ][n -1 X 2 (I -H n )X 2 ] .
This is a ratio of polynomials in X i that can be identified with [START_REF] Lee | Finite mixtures of multivariate skew tdistributions: some recent and new results[END_REF]. One easily sees that X 2 pn -1 + h 2 is O P (n -1/2 ) + O P (h 2 ) and hence all the coefficients of the polynomial in the numerator are (at least) O P (n -1/2 ) + O P (h 2 ). By Lemma 2, n -1 X 2 (I -H n )X 2 tends to EX 4 1 -1 which nonzero by assumption. Therefore the desired equation ( 22) is satisfied with

v -2 n = n -1/2 + h 2 , A n = -(X 2 -X 2 )[X 2 pn -1 + h 2 ] B n = -(XX 2 -X 3 )[X 2 pn -1 + h 2 ] C n = (-XX 3 + X 2 2 )[X 2 pn -1 + h 2 ] E n = X 2 n -1 X 2 (I -H n )X 2 .
Indeed, v -2 n a 2 n = n -1/2 a 2 n + h 2 a 2 n → 0 by the assumptions in Theorem 1. Let us show that A n is O P (v -2 n ). We have

v 2 n A n =O p (v 2 n n -1/2 ) + O p (v 2 n h 2 ) =O p 1 1 + n 1/2 h 2 + O p 1 1 + n -1/2 h -2 =O p (1),
the last equality holding because the sequence (1 + n 1/2 h 2 ) -1 is bounded. The remaining conditions in Lemma 3 are checked in the same way. The proof of Lemma 3 is complete. Hence the proof of Theorem 1 is complete, too.

Example 3 .

 3 Let h n ∝ n -1/5 . The Laplace distribution g(x) = exp(-|x|/b)/(2b), b > 0, x ∈ R, satisfies the conditions in Theorem 1 with a n = b and b n = b log(n/2).

Figure 1 :

 1 Figure 1: Averaged error values for the various parameters and algorithms in terms of the sample size. 13

Figure 2 :

 2 Figure 2: Pointwise averaged marginal density estimates along the second axis in the first cluster for GMM and fKDE. The true underlying density is added for comparison.

  observed log-likelihood

Figure 3 :Figure 4 :

 34 Figure 3: Trajectories of the errors (from (a) to (e)) and log-likelihoods averaged over the replications. The x-line is the number of steps and the y-line the value of the error or observed log-likelihood. 15

Figure 5 :

 5 Figure 5: Isocontours of estimated densities under different copula assumptions: Frank copula (a), Clayton copula(b), Gaussian copula(c) and independence(d).

Figure 6 :

 6 Figure 6: Results of the semiparametric model for the fit in four dimensions: (a)Log-likelihood value of the complete Iris data; (b)isocontours of the estimated density, component by component (the estimated density is built on the four variables of the iris data and this is the coordinatewise projection on the first and third variables).

Figure 7 :

 7 Figure 7: Cluster profiles for the Poisson mixture model, the Gaussian mixture model and the semiparametric copula-based mixture models with Frank and Gauss copulas. Each bar represents a cluster and each color represents a mean normalized expression profile, the value of which is given by the bar length of a given color. The bar widths are proportional to the estimated cluster proportions.

A. 1 3 From [ 9 ]Lemma 1 .

 1391 Computation of the normalizing constants in Example , p. 155, we know that [E(c n x + d n ; 1/b)] n → Λ(x), n → ∞, x > 0, where E(x; 1/b) = 1 -exp(-x/b), b > 0 is the distribution function of the exponential distribution, Λ(x) = exp(-e -x ) is the distribution function of the Gumbel distribution and c n = b, d n = b log n. Let L(x; b) = exp(x/b)/2, x > 0, be the distribution function of the Laplace distribution on the positive real line. Let a n = c n , b n = d n -b log 2 and x > 0. By identification of the binomial coefficients in the binomial theorem, we have[L(a n x + b n )] n = [E(c n x + d n )] n → Λ(x),meaning that a n = b and b n = b log(n/2) are the appropriate constants. If x < 0, the same formula applies because a n x + b n → ∞.A.2 Proof of Theorem 1Theorem 1 shall be proved by first considering the optimization problem (14)-(15) without the constraint p ≥ 0. (This shall be called the simplified optimization problem.) Throughout the proofs, the bandwidth sequence h n is simply denoted by h. Let n ≥ 3. If h → 0 and nh → 0 then the solution pn of the simplified problem such that M n p = b n (18) obeys pn =p n -

Lemma 3 .

 3 Let v n be a positive sequence satisfying v -1 n → 0, v -1 n a n → 0, v -1 n b n → 0. There exist random quantities A n , B n , C n , D n , E n such that, as n → ∞, A n is O P (v -2 n ), B n , C n , D n are O P (v -1 n

  shows the

	error	0.4 0.5 0.6 0.7 0.8			fKDE cKDE	error	0.020 0.025 0.030 0.035 0.040		GMM fKDE cKDE
		300	500	700	900		300	500	700	900
			sample size				sample size
			(a) copula parameters			(b) proportion parameters
		0.22			GMM fKDE		0.22
					cKDE
		0.18					0.18
	error					error
		0.14					0.14
								GMM
		0.10					0.10	fKDE cKDE
		300	500	700	900		300	500	700	900
			sample size				sample size
			(c) location parameters				(d) scale parameters
	error	0.16 0.17 0.18 0.19 0.20	fKDE cKDE			error	0.055 0.060 0.065 0.070 0.075 0.080	GMM fKDE cKDE
		300	500	700	900		300	500	700	900
			sample size				sample size
			(e) density generators				(f) misclassification errors

Table 1 :

 1 Averaged missclassification error of the computer experiment under eight situations, in %.

			5 1/6 µ 2,4 = 8	6	4	2	0
	algo: cKDE	7	7	8	7 10 12 16 21
	fKDE	8	8	9	7 11 11 15 25
	GMM	9	10	11	8 11 14 18 23

Acknowledgment. The author thanks two anonymous referees and an associate editor who made suggestions that helped to improve this paper.

meaning that M n has full rank with probability one. Since n ≥ 3 this rank must be three. Hence M n M n has full rank equal to three and therfore is invertible. The optimization problem is convex hence there is a unique solution pn whose expression is easily found: the Lagrangian writes p p -λ( M n -b n ) for some λ > 0 and by equating its gradient to zero we get

In order to obtain the desired formulas [START_REF] Kim | Mixture of D-vine copulas for modeling dependence[END_REF] and [START_REF] Kosmidis | Model-based clustering using copulas with applications[END_REF] it is convenient to introduce

so that we have the decompositions by blocks:

M n be the projection matrix onto the linear space spanned by the rows of M n . With this notation, we have

) and applying formula [START_REF] Krupskii | Factor copula models for multivariate data[END_REF] then yields [START_REF] Kim | Mixture of D-vine copulas for modeling dependence[END_REF] with

bn , this last equality being equivalent to [START_REF] Kosmidis | Model-based clustering using copulas with applications[END_REF]. We now introduce an intermediate lemma in order to facilitate the study of remainder terms which shall appear in the proof of consistency. Lemma 2. Let (Z n,1 , . . . , Z n,n ) be i.i.d. random variables defined on the same probability space as X 1 , . . . , X n . They are assumed to obey n

To see why the limit holds, note that n( M n M n ) -1 converges elementwise to the identity matrix. We now prove the consistency statement of Lemma 1. We have ĝ(x) = g(x)+ ĝ(x) -g(x) with g(x) =

we easily get that g(x) → g(x) in probability. Now using ( 20)-( 21) and Lemma 2 we also get

where K = (K h (x -X 1 ), . . . , K h (x -X n )) . The proof of Lemma 1 is complete.

Proof of Theorem 1. In this proof, the symbol pn stands for the solution of the optimization problem ( 17)- [START_REF] Hu | Maximum likelihood estimation of the mixture of log-concave densities[END_REF], that is, without the positivity constraint, and the symbol p+ n stands for the solution of the optimization problem ( 14)- [START_REF] Grønneberg | The copula information criteria[END_REF], that is, with the positivity constraint. In view of Lemma 1, it is sufficient to show that

because, by definition of the optimization problems, this implies that P (p n,i = p+ n,i , i = 1, . . . , n) → 1 and therefore that the estimators are equal with probability tending to one.

We write pn,i = pn,i 1 + pn,i -pn,i pn,i and the proof will be complete if (i) P (p n,i ≥ 0, i = 1, . . . , n) → 1 and (ii) |(p n,i -pn,i )/p n,i | can be bounded above by a quantity which would not depend on i and would vanish asymptotically. We first show (i). We have

The first term in the right hand side is a O P (n -1 ) and does not depend on i. Now

where ∨ i X i is a compact notation for max{X 1 , . . . , X n }. By assumption, ∨ i a -1 n (X ib n ) converges in distribution. By symmetry, so does ∨ i -a -1 n (X i + b n ). Hence, by the continuous mapping theorem, the maximum of ∨ i a -1 n (X i -b n ) and ∨ i -a -1 n (X i + b n ) converges in distribution. Thus

=O P (n -1 ) + O P (n -1/2 a n )(O P (1) + a -1 n b n ).