[START_REF] Lonsway | Rape myths. in review[END_REF]defined rape myths as "attitudes and beliefs which are generally false but which are widespread and persistent and which serve to deny and justify male sexual assault against women" (p. 134). Yet despite awareness campaigns, these myths persist in the collective imagination. Women are still more blamed when they are close to their attackers [START_REF] Ben-David | Rape perceptions, gender role attitudes, and victimperpetrator acquaintance[END_REF]. While the perpetrator of a rape is -logically -the one and only person responsible, the attribution of responsibility is shared between the aggressor and the victim.

However, there is little research to explain the blame against women. In this research, we investigated the degree of blame of both victim and/or perpetrator with the attitudes towards victims of sexual assault (PCRS, [START_REF] Cowan | Cognitive style and attitudinal correlates of the perceived causes of rape scale[END_REF], the beliefs that we have about the events that happen to us or that happen to others (Just World Beliefs), dark personality traits (Machiavellianism, psychopathy and narcissism) and various socio-demographic variables (age, sex, etc.).

Conclusion

Unlike the few existing studies, we have not found a link between psychopathy and the factors investigated regarding perpetrator blame. It therefore seems difficult to identify stable factors explaining the blame of the aggressors. Conversely, the victim's blame depends on several factors, some of which have been previously identified in the literature, such as Machiavellianism. However, we have shown that this impacts the blame of the victim through male sexuality and female provocation. So, if the Machiavellians blame the victims, it is because they consider both that the women provoked their attackers, but also that the men -because the men -fail to control their sexual outbursts. These same mediators also allow Belief in a Just World to be expressed: the more people are considered responsible for their behavior, the more female provocation is perceived as a factor increasing the blame of the victim.

MethodResultsFigure 1 .

 1 Figure 1. Antecedent Factors and Mediators of the Attribution of Responsibility to the Perpetrator and the Victim.

Factors Modulating Responsibility (Victim and Aggressor) in Rape

		University of Tours
		International Conference
	"From Seduction to Sexual Assault? Questioning the harassment."
		November 2021, 29 and 30
	Annick Lemoine (PhD candidate)	Samuel Demarchi (Assistant Professor, PhD)
	Université Paris 8	Université Paris 8
	Laboratoire CHArt (EA 224)	Laboratoire CHArt (EA 224)
	annicklemoine7595@gmail.com


Emilie Lebarbier

Nicolas Marie

Amélie Rosier

TRENDTM: A R PACKAGE FOR THE TREND OF HIGH-DIMENSIONAL TIME SERIES ESTIMATION

In this paper, we present the R package TrendTM dedicated to the trend estimation of high dimensional time series matrices. The main features of this package is the possibility to take into account different types of temporal structures on the data (none, smooth or periodic) and to select, according to the user's wishes, either the rank k of the matrix, the parameter τ linked to the temporal structure or both. Then, a two-stage heuristic is provided for the joint selection of the rank k and of the tuning parameter τ . Moreover, the TrendTM function is applied to time series clustering and principal component analysis on real datasets. The package is available on the CRAN.

Introduction

Since the 1970's, it is usual to model a one-dimensional time series by a process (X t ) t∈Z satisfying [START_REF] Alquier | Matrix Factorization for Multivariate Time Series Analysis[END_REF] F θ (X t+q , . . . , X t-q , η t+p , . . . , η t-p ) = 0 ; t ∈ Z, where p, q ∈ N, η = (η t ) t∈Z is a second order stationary process, often a white noise, and F = (F θ ) θ∈Θ is a family of continuous maps from R 2(p+q+1) into R indexed in a set Θ. For instance, ARMA models, GARCH models, and all their extensions are defined this way. An advantage of Model [START_REF] Alquier | Matrix Factorization for Multivariate Time Series Analysis[END_REF] is that F can be chosen in order to take into account properties known on the dynamics of the modeled phenomenon regardless to the datas. However, except in simple cases, Model (1) is difficult to extend to the high dimensional framework. It is also difficult to bypass the stationary condition on η. For a good reference on the classic time series models, see Gouriéroux and Monfort [START_REF] Gouriéroux | Time Series and Dynamic Models[END_REF]. For an introduction to GARCH models, the reader can refer to Brockwell and Davis [START_REF] Brockwell | Introduction to Time Series and Forecasting[END_REF], Chapter 7. Finally, for an introduction to time series analysis with R, see Cowpertwait and Metcalfe [START_REF] Cowpertwait | Introductory Time Series with R[END_REF]. Independently, for almost two decades, in particular thanks to the Netflix challenge on movies recommendations, the matrix factorization for the completion and denoising of high dimensional matrices with i.i.d. entries has been deeply investigated on the theoretical side (see [START_REF] Cai | ROP: Matrix Recovery via Rank-One Projections[END_REF][START_REF] Klopp | Robust Matrix Completion[END_REF][START_REF] Klopp | Structured Matrix Estimation and Completion[END_REF][START_REF] Koltchinskii | Nuclear-Norm Penalization and Optimal Rates for Noisy Low-Rank Matrix Completion[END_REF][START_REF] Moridomi | Tighter Generalization Bounds for Matrix Completion via Factorization into Constrained Matrices[END_REF]). Recently, Alquier and Marie [START_REF] Alquier | Matrix Factorization for Multivariate Time Series Analysis[END_REF] (resp. Alquier, Marie and Rosier [START_REF] Alquier | Tight Risk Bound For High Dimensional Time Series Completion[END_REF]) has extended this method in order to take into account time series trends properties as periodicity in the denoising (resp. completion) process of time series matrices. On the theoretical side, to take into account trends properties in the definition of the denoising and completion estimators allowed the authors to improve existing risk bounds.

The TrendTM package, for "Trend of High-Dimensional Time Series Matrix Estimation", implements the denoising method investigated on the theoretical side in Alquier and Marie [START_REF] Alquier | Matrix Factorization for Multivariate Time Series Analysis[END_REF]. The selection of both the rank k of the trends matrix and the parameter τ controlling the temporal structure (the period for seasonal time series) of the dataset is a model selection issue, that can be difficult to perform in practice when some unknown constants exist in the associated criterion and in the particular context of two parameters to be selected. This issue is well-managed by the TrendTM package by using the strategy of Devijver et al. [START_REF] Devijver | Nonlinear network-based quantitative trait prediction from transcriptomic data[END_REF]. Moreover, since the denoising method implemented in the TrendTM package is based on matrix factorization algorithms, natural applications to time series clustering and principal component analysis (PCA) are investigated in this paper. To illustrate the efficiency of the package in these application fields, we present an application to ECG profiles and children growth profiles clustering.

The paper is organized as follows. Section 2 recalls the estimation procedure proposed by Alquier and Marie [START_REF] Alquier | Matrix Factorization for Multivariate Time Series Analysis[END_REF]. Section 3 discusses the model selection issue in practice. More precisely, data-driven existing 1 heuristics are recalled and presented for our joint selection on the rank and the trend parameter whose performances are studied in Section 4 on simulated data. Section 5 gives some precisions and guidelines of the proposed method in the package TrendTM and shows an application on real data. In Section 6, we show that this method can be used for two classical statistical problems which are the time series clustering and PCA. As mentioned above, illustrations on real datasets are given.

The trend estimation procedure

In this section, we recall the model and the estimation procedure proposed by Alquier and Marie [START_REF] Alquier | Matrix Factorization for Multivariate Time Series Analysis[END_REF]. When we want to take into account for a possible trend property in the time series, we assume that M = MΛ, where M is a d × τ matrix of low rank k (thus with τ > k) and Λ is a known τ × n full-rank matrix. Two trend's property cases are considered:

• periodic series: if the trend of X is τ -periodic then Λ = (I τ | • • • |I τ ) where I τ is the identity matrix in M τ,τ (R), • smooth series: if the form of the trend is t ∈ {1, . . . , n} → f (t/n) with f ∈ L 2 ([0, 1]; R d ), then Λ = ϕ t n ( ,t)∈{1,...,τ }×{1,...,n}
, where τ is odd and (ϕ 1 , . . . , ϕ τ ) is the τ -dimensional trigonometric basis defined by

ϕ (x) :=    1 if n = 1 √ 2 cos(2πmx) if = 2m √ 2 sin(2πmx) if = 2m + 1
for every x ∈ [0, 1] and m ∈ {1, . . . , (τ -1)/2}.

We define an auxiliary model that will be used for the estimation of M according to model [START_REF] Alquier | Tight Risk Bound For High Dimensional Time Series Completion[END_REF] and the definition of M that is

(3) X = M + ε,
where X := XΛ + , ε := εΛ + and Λ + = Λ * (ΛΛ * ) -1 is the Moore-Penrose inverse of Λ. This model doesn't embed some trend's property anymore. Note that the case with no trend corresponds to take τ = n and Λ = I n .

2.2. Estimation procedure. The estimation procedure consists in two steps: estimate M for fixed k and τ , then choose these parameters which is a model selection issue.

Step 1: Estimation of M, k and τ being fixed. Consider the least-square estimator of the matrix M, that is

(4) M k,τ ∈ arg min A∈S k,τ X -A 2 F ,
where . F is the Frobenius norm (for a matrix A, A F := trace(AA * ) 1/2 ) and S k,τ ⊂ {UV ; U ∈ M d,k (R) and V ∈ M k,τ (R)}. So, a natural estimator of M is given by

M k,τ := M k,τ Λ.
Step 2: Choice of k and τ . Alquier and Marie [START_REF] Alquier | Matrix Factorization for Multivariate Time Series Analysis[END_REF] (see Section 4) provides a data-driven selection criterion. Precisely, for a fixed s > 0, the final estimator of M is M s := M k(s), τ (s) where

( k(s), τ (s)) ∈ arg min (k,τ )∈K×T { X -M k,τ 2
F + pen s (k, τ )}
with K ⊂ {1, . . . , d ∧ n}, T ⊂ {1, . . . , n}, and for (k, τ ) ∈ K × T ,

(5)

pen s (k, τ ) := c pen Σ ε op k(d + τ + s)
where c pen > 0 is a deterministic constant and the operator norm is given by A op = sup x =1 Ax with . the Euclidean norm on R n . Alquier and Marie [START_REF] Alquier | Matrix Factorization for Multivariate Time Series Analysis[END_REF] (Theorem 4.1) ensures that for every θ ∈ (0, 1), with probability larger than 1 -2e -s ,

M s -M 2 F min (k,τ )∈K×T min A∈S k,τ 1 + θ 1 -θ 2 AΛ -M 2 F + 4 θ(1 -θ) 2 pen s (k, τ ) .

Heuristic for the calibration of the penalty constants

The penalty function defined by (5) depends on some constants s and c pen that must to be chosen or calibrated in practice. If the constant s can be easily chosen, the choice of the fixed c pen is a difficult task. It is now usual to calibrate this constant in a data-driven manner. In the literature, different and now well known heuristics have been developed to this purpose. However, they are all dedicated to the selection of one dimension or parameter, that is our case if we aim at selecting k for a fixed τ or τ for a fixed k. In the case of the selection of two parameters, Devijver et al. [START_REF] Devijver | Nonlinear network-based quantitative trait prediction from transcriptomic data[END_REF] proposed a two-stage heuristic using for each step one of the previous heuristic. We first recall some heuristics for the selection of one parameter, then we present the two-stage heuristic for the joint selection of (k, τ ).

First, in practice, we could take s = -log((1 -α)/2) with α fixed to 99%, 95% or 90%. Here we choose to fix s = 4. The penalty function is thus reduced to

pen(k, τ ) := c cal k(d + τ + s) ; ∀(k, τ ) ∈ K × T ,
where c cal > 0 is a penalty constant that needs to be calibrated. The resulting adaptive estimator is denoted by

M s = M k, τ .
Up to our knowledge, there exit the three following heuristics dedicated to the constant calibration question in the model selection frameworks of one parameter:

• the one proposed by Lavielle [START_REF] Lavielle | Using penalized contrasts for the change-point problem[END_REF], denoted here ML, that involves a threshold S which is fixed to S = 0.75 as suggested by the author. • the two proposed by Birgé and Massart [START_REF] Birgé | Gaussian model selection[END_REF] (see the more recent versions [START_REF] Arlot | Data-driven calibration of penalties for least-squares regression[END_REF] and [START_REF] Arlot | Minimal penalties and the slope heuristics: a survey[END_REF]) that are two versions of the well known slope heuristic: the 'dimension jump' and the 'slope', denoted here BJ and Slope respectively. The both heuristics have been implemented in the R package capushe described in [START_REF] Baudry | Slope heuristics: overview and implementation[END_REF].

For the joint selection of (k, τ ), following the strategy proposed by Devijver et al. [START_REF] Devijver | Nonlinear network-based quantitative trait prediction from transcriptomic data[END_REF], we propose the following two-stage heuristic: first, we choose the best τ for each k ∈ K via the criterion

τ (k) ∈ arg min τ ∈T { X -M k,τ 2
F + c cal,τ k(d + τ + s)},
where the penalty constant c cal,τ is calibrated using one of the previous heuristic, then we select the best k among them via the criterion

k ∈ arg min k∈K { X -M k, τ (k) 2 F + c cal,k k(d + τ (k) + s)},
where the penalty constant c cal,k is calibrated using the same heuristic to be consistant, and τ = τ ( k).

Note that, in practice K = {1, . . . , k max } and T = {k + 1, . . . , τ max } where k max is the maximal rank and τ max is the maximal value of τ where these two quantites need to be specified.

Simulation study

In this section, we first study the behavior of the three heuristics for the model selection issue for k and τ separatly (Study 1), then we illustrate the importance of accounting for the trend in the estimation procedure when it exists in the series (Study 2) and finally we study the perfomance of the proposed method when both k and τ are selected (Study 3) on simulated data. In this framework, we only consider the smooth trend case. (1) we generate a matrix M = UV by simulating

U ∈ M d,k (R) and V ∈ M k,τ (R)
for which the entries of U and V are assumed to be i.i.d. and follows a centered Gaussian distribution with same standard deviation σ uv fixed to 0.5; (2) two cases are considered according to the presence or not of a trend in the simulated series: if there is no trend, then τ = n and M = M, and otherwise M = MΛ with the matrix Λ of the smooth case. To distinguish between these two cases in the sequel, we call them dataset NoTrend and dataset Trend respectively; (3) the rows of the error matrix ε are assumed to be i.i.d. and to follow a centered Gaussian distribution of variance σ 2 (i.e. Σ ε = σ 2 I n ). We fix k = 3 and τ = 25. We consider different values for the residual standard deviation σ in order to have different levels of difficulty for the estimation problem. First, according to the previous considerations, Var(M ij ) = kσ 4 uv for dataset NoTrend and Var(M ij ) = τ kσ 4 uv for dataset Trend . Let us consider s ∈ {0.1, 0.5, 1.5, 2}. In order to have the same estimation difficulty (same ratio between σ and the standard deviation of M ij ) for the two datasets, we set σ = s for dataset NoTrend and σ = √ τ s for dataset Trend . The obtained four cases are judged as 'Easy', 'Medium', 'Difficult' and 'Hard' respectively. For each combination of parameters, we've simulated 200 datasets.

Let us precise that when the trend is not considered in the estimation procedure, the resulting estimator is M k or k,n (if k is selected or not) and when it is considered the resulting estimator is M k or k,τ or τ (if both k and τ are selected or one of them or none). 4.1.2. Quality criteria. The performance of our procedure is assessed via:

• the estimated k and/or τ ; and

• the squared Frobenius norm between M and its estimate M k, τ . Moreover we also consider the Frobenius norm between M and

• the estimator of M for the true k and/or τ , that is M k,τ ; and

• the trajectorial oracle that is M k,τ where ( k, τ ) = arg min (k,τ )∈K×T M -M k,τ 2
F when both k and τ are selected, M k,τ where τ = arg min τ ∈T M -M k,τ 2
F when k is fixed, and M k,n where k = arg min k∈K M -M k,n 2
F when no trend is considered. 4.2. Study 1: behavior of the three heuristics for the selection of k or τ . We first study the selection of k for dataset NoTrend when no trend is considered in the estimation procedure. We consider two different values of the maximal rank k max ∈ {15, 35}. The results are presented in Figure 1. When the noise is small, i.e. the estimation problem is easy (cases 'Easy' and 'Medium'), all the heuristics recover the true rank, and therefore the obtained estimators perform as well as M k,n (the estimator of M for k fixed to its true value). When the estimation problem gets more difficult (cases 'Difficult' and 'Hard'), the heuristics have a tendency to underestimate the rank. This underestimation behavior seems to be logical and even desirable in the particular 'Hard' case. Indeed, we observe that in terms of Frobenius norm, the obtained estimators perform better compared to the one with the true rank. Moreover, they have performance close to the oracle. Comparing the three heuristics, the Slope heuristic shows better performance compared to the two other heuristics. This is particularly marked for the 'Medium' case and k max = 15. We can note that the behavior of the three heuristics can be affected by the choice of k max . This problem is well known for both the BJ and Slope heuristics (see Arlot [START_REF] Arlot | Minimal penalties and the slope heuristics: a survey[END_REF] for more explanations in the case univariate series analysis).

Then, we study the selection of τ for dataset Trend for k fixed to the true value. We fix τ max = 55. The results are presented in Figure 2. Except with BJ that is more unstable, the heuristics retrieve the true value of τ whatever the estimation difficulty with same performance as the oracle.

From this study, we choose the Slope heuristic for the model selection issue for both k and τ in the sequel and in the developed package. 4.3. Study 2: accounting for the smooth structure in the trend. We compare the performance of the procedure on the dataset Trend when the trend is considered (τ = τ ) or not (τ = n) for k fixed to the true value. We choose τ max = 55. The results are represented in Figure 3. Whatever the difficulty of the estimation problem (different values of σ), accounting for the trend increases the precision of the estimation. This is more marked for high values of σ. Note that, similarly as Study 1, the estimation naturally degrades with the increasing of σ.

4.4. Study 3: selection of k and τ . Figure 4 shows that the joint heuristic retrieves the true values of k and τ whatever the difficulty of the estimation problem, except very few times. Thus, the perfomance of the estimator of M is comparable to the one of the estimator M k,τ and moreover it has performance close to the oracle (see Figure 5). Compared to Study 1 where τ = n, here for difficult estimation problems, k is not underestimated.

Using package TrendTM

5.1.

Comments on the package. The package is organized around the main and unique function TrendTM. In this section, we present the arguments used in a call to this function:

TrendTM(X, k.select=FALSE, k.max=20, struct.temp="none", tau.select=FALSE, tau.max=floor(n/2), type.soft="als") This function returns a list containing six elements:

• k.est, the estimated k or the true k when no selection is chosen;

• tau.est, the estimated τ or the true τ when no selection is chosen;

• U.est, the component U of the decomposition of the final estimator of M;

• V.est, the component V of the decomposition of the final estimator of M;

• M.est, the estimator of M;

• contrast, the squared Frobenius norm of X-M.est. If k and τ are fixed, the contrast is an unique value; if k is selected and τ is fixed or if τ is selected and k is fixed, the contrast is a vector containing the norms for each visiting values of k or τ respectively; and if k and τ are selected, the contrast is a matrix with k max rows and τ max colums such that contrast k,τ = X -M k,τ

Model selection.

The selection of k or/and τ is requested using the options k.select or/and tau.select that are booleans. When there is no selection, the option is set to FALSE and k.max= k or/and tau.max= τ . Note that if no trend is considered in the estimation procedure, τ = n, otherwise tau.max must be a smaller than n and larger than k.max+2 in order to ensure that the rank of M is k.

5.1.3.

Accounting for a trend. Let us give more details about the different arguments of TrendTM that need to be specified when accounting for a temporal structure in the estimation procedure.

q q q q q q q q q q q q q q q q q q 0 200 400 600 q q q q q q q q q q q q q q q q q q q q q q q 3 4 5 6

Easy

Medium Difficult Hard

CaseDiff (select_k,select_tau) (true_k,true_tau) oracle q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 25.0 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 100 200 Two temporal structures are considered: periodic trend and smooth trend. This can be specified using the option struct.temp, struct.temp="periodic" or struct.temp="smooth" respectively. Recall that the selection of τ is only possible when a smooth trend is considered. Thus, when

• struct.temp="periodic", then tau.select=FALSE and tau.max= τ . In this case, τ must be such that n is a multiple of τ ; • struct.temp="smooth", then tau.select is either FALSE or TRUE. Whatever this choice, tau.max must be an odd number. When no trend is taken into account, struct.temp="none" and tau.max= n.

The function TrendTM returns M.est=U.est%*%V.est%*%Lambda, where the matrix Λ depends on the considered temporal structure. If needed, the code to compute this matrix and its Moore-Penrose inverse Λ + (denoted below LambdaP) for both temporal structures and tau.max= τ is if (struct.temp=="periodic"){ p <-n/tau.max A <-matrix(1, ncol=p,nrow=1) B <-diag(tau) Lambda <-kronecker(A,B) LambdaP <-kronecker(t(A),B)*(1/p) } if (struct.temp=="smooth"){ times.eval <-seq(1/n,1,by=1/n) fbasis_obj= fda::create.fourier.basis(rangeval=c(0,1),nbasis=tau,period = 1) fbasis_evals <-fda::eval.basis(times.eval, fbasis_obj) Lambda <-t(fbasis_evals) LambdaP <-t(Lambda)/n } 5.1.4. Estimation of M, k and τ being fixed (Step 1). The least-square estimator of M, M k,τ given by ( 4), is obtained by using the softImpute function from the R package of the same name developed by Hastie and Mazumder for matrix completion [START_REF] Hastie | softImpute: Matrix Completion via Iterative Soft-Thresholded SVD[END_REF]. In this package, two algorithms are implemented: 'svd' and 'als'. In a simulation study, we observed that they have both provided the same accuracy of the estimator (results not shown). We decide to use the 'als' algorithm by default but the choice is left free to the user in our package TrendTM. In this package, this choice is specified using the option type.soft.

5.1.5.

The Slope heuristic. Let us now focus on the selection problem of the rank k, and write the penalty as pen(k) = c cal,k ϕ(k). The Slope heuristic, proposed by [START_REF] Birgé | Gaussian model selection[END_REF], consists in estimating the slope s of the contrast X -M k,n 2 F as a function of ϕ(k) with k 'large enough' and defining c cal,k = -2 s. The implementation of this heuristic requires the choice of the dimensions on which to perform the regression, that can be difficult in practice. To deal with this problem, Baudry et. al [START_REF] Baudry | Slope heuristics: overview and implementation[END_REF] proposed to make robust regressions for dimensions between k and k max for k = 1, 2, . . ., resulting in different selected k. The choice of the final dimension is the maximal value k such that the length of successive same k is greater than the option point of the function DDSE. In order to avoid some implementation problems as such condition is not reached and no k is selected, we decide to take the value k associated to the maximal length of successive same k.

Application to pollution dataset.

Let us use the package on a real dataset1 . The dataset contains the amounts of d = 13 toxic gases in the air recorded n = 9357 times during one year. We do a first step of data imputation using the function complete of the package softImpute since missing values (coded with -200) exist in this dataset (see [START_REF] Alquier | Tight Risk Bound For High Dimensional Time Series Completion[END_REF] for more details on high dimensional time series completion). Our function trendTM is then applied on the completed matrix with a selection of both k and τ . The code is the following: ## Used package library('softImpute') library('TrendTM') ## Data importation AirPollution <-read.table('AirQualityUCI.csv',sep=";",header=TRUE)[,-c(1:2)] AirPollution <-t(as.matrix(AirPollution)) ## Imputation AirPollution[AirPollution==-200] <-NA AirPollutionSoft=softImpute::softImpute(AirPollution,rank=1,lambda=0) AirPollutionImp=softImpute::complete(AirPollution,AirPollutionSoft) ## Trend estimation Trend.AirPollution = TrendTM(AirPollutionImp, k.select=TRUE, k.max=13, struct.temp="smooth", tau.select=TRUE, tau.max=101) Trend.AirPollution$k.est [START_REF] Alquier | Matrix Factorization for Multivariate Time Series Analysis[END_REF] 7 Trend.AirPollution$tau.est [START_REF] Alquier | Matrix Factorization for Multivariate Time Series Analysis[END_REF] 13

The procedure selects k = 7 and τ = 13. Figure 6 shows the obtained trend estimation for 4 toxic gases among the 13. The denoising process seems to have been well applied to the data.

Method used for usual problems in statistics

This section deals with an application of Model (2), first to times series clustering, and then to PCA of times series. These applications were not developed in Alquier and Marie [START_REF] Alquier | Matrix Factorization for Multivariate Time Series Analysis[END_REF].

6.1. Application to time series clustering. To assume that M = UL with L = ( j (t)) j,t := VΛ means that for any i ∈ {1, . . . , d}, the trend m i (.) of the i-th row of X satisfies

m i (t) = k j=1
U i,j j (t) ; ∀t ∈ {1, . . . , T }.

In other words, the trend of the i-th time series stored in X is a linear combination of the trends 1 (.), . . . , k (.) of k latent time series. Moreover, since L = VΛ, 1 (.), . . . , k (.) have the same usual time series trend's property, characterized by Λ, than m 1 (.), . . . , m d (.). The identification of the aforementioned latent time series is part of curves or time series clustering framework, where k corresponds to the number of clusters. To this aim, we propose here to apply the well-known Hierarchical Agglomerative Clustering with the Ward's linkage on the rows of U k , where

M k,τ = U k V τ or M k,τ = U k V τ Λ.
In [START_REF] Jacques | Functional data clustering: a survey[END_REF], Jacques and Preda applied and compared different curves clustering methods proposed in the literature on real datasets. We propose to apply our clustering strategy on two datasets: the ECG and Growth datasets. The ECG dataset (taken from the UCR Time Series Classification and Clustering website) consists in 200 electrocardiogram from 2 groups of patients sampled at 96 time instants in which 133 are classified as normal and 67 as abnormal, and the Growth dataset (available in the R package fda) contains the heights of 54 girls and 39 boys measured at 31 stages from 1 to 18 years. In both datasets, we have two groups. The times series of the two datasets are plotted in Figure 7.

We apply our proposed clustering strategy based on the procedure without and with taking into account for a trend (periodic for ECG and smooth for Growth) with the known number of groups k = 2. The Correct Classification Rates (CCR) according to the known partitions are given in Table 1 for the ECG dataset and in Table 2 for the Growth dataset. We also report the CCR obtained for the best method among the ones tested in [START_REF] Jacques | Functional data clustering: a survey[END_REF]. In addition, we indicate the time taken by the different methods on a laptop 1.6 GHz CPU (note that the time for the method HDDC on FPCA scores is not given since this method is not available). For the ECG dataset, accounting for the trend improves significantly the clustering performance. This is not the case with the Growth dataset which already has a very high CCR without trend. Our clustering performances are the same compared to the best clustering methods but it is much more faster. Note that among the compared methods in [START_REF] Jacques | Functional data clustering: a survey[END_REF], the best ones for the two datasets are not the same.

The R code of our clustering strategy for the Growth dataset without accounting the trend is the following library(fda) data("growth") Growth <-t(cbind(matrix(growth$hgtm,31,39),matrix(growth$hgtf,31,54))) cls <-c(rep(1,39),rep(0,54)) k.max <-2 res.Growth <-TrendTM::TrendTM(Growth,k.max=k.max) Uf.cr <-scale(res.Growth$U.est,scale=TRUE, center=TRUE) Uf.ward<-hclust(dist(Uf.cr, method = "euclidean")^2, method = "ward.D") cluster.Uf <-cutree(Uf.ward, k =k.max) table(cls,cluster.Uf) Note that for the Growth dataset without trend, the Slope heuristic selects k = 6 clusters and for the ECG dataset with a 32-periodic trend, it selects k = 7. In Figure 8, the series of the ECG dataset are 6.2. Application to time series PCA. The PCA problem can be rephrased as a k-rank approximation of a matrix X (information captured by the first k axes). More precisely, if the matrix X is of dimension d × n, the simple PCA solution of rank k is given by

M k,n = U k V k ∈ arg min U∈M d,k ,V∈M k,n X c -UV 2 F ,
where X c is the centered matrix (if A j denotes the j-th column of the matrix A, X c j = X j -Xj with Xj = d i=1 X ij /d). The matrix of principal components is thus X c V * k = U k . Each line of this d × k matrix contains the coordinates of the projection of the associated time series on the first k axis. Two projected time series are close if they share globally the same trend's property. However, in high-dimensional spaces, the Euclidean distance used in PCA can lose its meaning and a local trend similarity could be preferred. We thus propose to perform the PCA on the transformed matrix X c = X c Λ + . The solution is M k,τ = M k,τ Λ = U k V τ Λ and the matrix of principal components is

X c V * k = X c Λ + V * k .
The projections of the times series on the principal plan (coordinates given by the two first columns of the matrix of principal components) of X c and X c are plotted in Figure 9 on the left and the right respectively for the ECG dataset colored according to the true partition. Recall that for the ECG dataset, we consider that the trend is periodic with period τ = 32.

To illustrate the trend reduction using Λ + on a period with length τ , the 28th and the 121th time series of X c and X c are plotted in Figure 10 on left and right respectively. As expected according to their features in both matrices, they are not close in the PCA of X c whereas they are close in the PCA of X c . That also explains the difference when the custering is performed without trend or with trend (see the previous section). For example, these two time series are clustered in the same group when a trend property is taking into account, whereas they are not in the same group when no trend property is considered.

4. 1 .

 1 Simulation design and quality criteria.

4. 1 . 1 .

 11 Simulation design. We simulated datasets with d = 100 and n = 600 as follows:

2 F . 5 . 1 . 1 .

 2511 Data structure. The structure of the data is a matrix with d rows and n columns containing the d time series with size n each. A simulated dataset without temporal structure is available in the package with d = 30 and n = 100. The user can load the data using the following command:

kFigure 1 .Figure 2 .

 12 Figure 1. Comparison of the three heuristics for the selection of k for dataset NoTrend . Left: estimated number of the rank k and right: boxplot of M -M k,n F for two values of k max = 15 (first line) and k max = 35 (second line), and different values of σ. On each graph and for each value of σ, from left to right, we have the result from ML, BJ, Slope ( k), the true rank (k) and the oracle ( k).

Figure 3 .

 3 Figure 3. Boxplot of M -M k,τ F with τ = τ (select_tau) and τ = n (tau=n) for different values of σ.

Figure 4 .

 4 Figure 4. Left: estimated k. Right: estimated τ for different values of σ.

Figure 5 .

 5 Figure 5. Boxplot of M -M k,τ F with (k, τ ) = ( k, τ ) the selected k and τ , (k, τ ) = (k, τ ) the true values and (k, τ ) = ( k, τ ) the oracle for different values of σ.

Figure 6 .

 6 Figure 6. Estimation of the trend of 4 toxic gases (red) compared to the initial series (grey).

Figure 7 .

 7 Figure 7. Left: Growth dataset (black: boys, red: girls).Right: ECG dataset (black: normal, red: abnormal).

  plotted and colored according to the k = 7 obtained clusters with a 32-periodic trend on the left, and the associated obtained trend m i (.) for i = 1, . . . , 6 are plotted on the right.

Figure 8 .

 8 Figure 8. Left: ECG dataset colored according to the k = 7 obtained clusters with a 32-periodic trend. Right: the associated trends m i (.) for i = 1, . . . , 6.

Figure 9 .

 9 Figure 9. PCA on the ECG dataset. Left: on X c . Right: on X c with a periodic (τ = 32) trend.

Figure 10 .

 10 Figure 10. The 28th (in black) and the 121th (in red) time series. Left: on X c . Right: on X c with a periodic (τ = 32) trend.

  2.1. A high dimensional matrix based time series model. Consider the model

	(2)	X = M + ε,

where X is an observed d × n matrix which rows are time series, M is a deterministic matrix of low rank k ∈ N * (i.e. k d ∧ n), and ε is a random matrix which rows are i.i.d. samples of a centered Gaussian stochastic process with a covariance matrix Σ ε .

Table 1 .

 1 Correct classification rates (CCR) in percentage accounting or not for a trend on the ECG dataset. Mean times in second obtained on 30 runs.

		proc NoTrend	proc Trend	Best method in [14]
			τ = 32-periodic	Funclust
	CCR	74.5	83	84 (reported from [14])
	mean times in second (on 30 runs)	0.015	0.008	19.2
	proc NoTrend	proc Trend	Best method in [14]
		smooth with τ = 13	HDDC on FPCA scores
	CCR	97.85	97.85	97.85 (reported from [14])
	mean times in second (on 30 runs)	0.003	0.0028	.

Table 2 .

 2 Correct

classification rates (CCR) in percentage accounting or not for a trend on the Growth dataset. Mean times in second obtained on 30 runs.

available at https://archive.ics.uci.edu/ml/datasets/Air+quality

† Laboratoire Modal'X, Université Paris Nanterre, Nanterre, France Email address: emilie.lebarbier@parisnanterre.fr Email address: nmarie@parisnanterre.fr

Acknowledgments

This work has been conducted as part of the project Labex MME-DII (ANR11-LBX-0023-01) and within the FP2M federation (CNRS FR 2036).

ESME Sudria, Ivry-sur-Seine, France Email address: amelie.rosier@esme.fr