Mathias Bourel

Jairo Cugliari

Bagging of Density Estimators

Keywords:

In this work we give new density estimators by averaging classical density estimators such as the histogram, the frequency polygon and the kernel density estimators obtained over different bootstrap samples of the original data. We prove the L 2consistency of these new estimators and compare them to several similar approaches by simulations. Based on them, we give also a way to construct non parametric pointwise confidence intervals for the target density.

Introduction

Ubiquitous in data analysis, density estimation techniques are certainly the most used unsupervised learning technique on low dimension. Whether for studying asymmetry, normality, residual diagnostic or bump hunting among others, one usually relies on a visual inspection of a plot of the density to take a primary decision, mostly in one or two dimensions.

The general aim is to gather basic information about the unobserved data generation mechanism out of a sample of n observations say x 1 , . . . , x n . One usually supposes that the observations are realizations of a random variable X that admits a probability density function f (i.e. f is non negative and integrates 1). Then, the learning task is to estimate f as accurately as possible. First, by obtaining a point wise estimate f (x) of f (x) for all x ∈ R d , and second, by assessing the uncertainty of the point wise estimate through the construction of a confidence interval for f . Of course both problems implies different difficulties and involves specific techniques. In what follows, we focus on nonparametric approaches for the first objective and give a possible way to construct a pointwise confidence interval for the true density. In particular we center our attention on three classes of base or individual estimators for the density: histograms, frequency polygons and the kernel density estimator. We postpone the formal definition up to the next section but we provide a discussion guided by intuitive descriptions here.

Histograms are undoubtedly the most popular construction for density estimation. They rely on the best constant-wise approximation of f given by the data using a binning argument. The power of this simple construction combined with the ease of its interpretation makes them accessible to non technical users. Besides, theoretical properties can be derived showing that histograms are consistent estimators. Some of the lacks of the histograms are inherent to the constant level for at each partition, on one hand their discontinuities and on the other hand they have null derivative everywhere. Frequency polygons are constructed on top of histograms (and so take profit of the binning advantages) providing a linear piecewise estimator. Although the added regularity was once pointed out as a flaw (see [START_REF] Fisher | Statistical Methods for Research Workers[END_REF]) it is now well known that it increases the quality of the estimator. Theoretically this is shown by a faster rate of convergence (cf. section 2). Further regularity can be gained using kernel density estimators. Essentially one first picks a kernel function with the desired degree of regularity for the final estimate. Then, the empirical measure is convolved to produce the kernel density estimation of f . It has been proved in [START_REF] Scott | Frequency polygons: Theory and application[END_REF] that the rates of convergence of the frequency polygon are similar to the kernel density estimator.All the three approaches are exhaustively studied both from on practice and theory as individual estimators and a general reference for the subject is [START_REF] Scott | Multivariate Density Estimation: Theory, Practice, and Visualization[END_REF]. However, a reasonable question is to ask whether further improvement can be achieved from these construction by means of aggregation schemes.

Ensemble learning or aggregation methods are increasingly used in the supervised framework: these methods combine intermediate predictors to obtain an aggregated model with the aim to obtain a better estimator. Bagging [START_REF] Breiman | Bagging predictors[END_REF] (Bootstrap and AGGregatING), Boosting [START_REF] Freund | A decision-theoretic generalization of on-line learning and an application to boosting[END_REF], Stacking [START_REF] Wolpert | Stacked Generalization[END_REF], and Random Forests [START_REF] Breiman | Random forests[END_REF] have been broadly studied in the case of classification (principally binary classification) or regression from the theoretical viewpoint and have very high performances when tested over tens of various datasets selected from the machine learning benchmark. Several extensions are still under study: multivariate regression, multiclass classification, and adaptation to functional data or time series. Very few developments exist for ensemble learning for unsupervised techniques such as clustering analysis and density estimation. Only on few works several authors look at the adaptation of the aggregation procedure to estimate a density under somehow restrictive conditions. One of the first is the mean of the most simplest density estimator, the histograms, each one constructed over several different deterministic grids in Average Shifted Histograms, ASH [START_REF] Scott | Averaged shifted histogram: Effective nonparametric density estimators inseveral dimensions[END_REF]. With a combination of several kernel density estimators with different bandwidths, often in a normal context, and varying the form of the aggregation we can cite [START_REF] Ridgeway | Looking for lumps: Boosting and bagging for density estimation[END_REF], [START_REF] Glodek | Ensemble gaussian mixture models for probability density estimation[END_REF], [START_REF] Song | Density boosting for gaussian mixtures[END_REF], [START_REF] Rosset | Boosting density estimation[END_REF], [START_REF] Smyth | Linearly combining density estimators via stacking[END_REF] and [START_REF] Rigollet | Linear and convex aggregation of density estimators[END_REF].

Another kind of aggregation can be obtained by introducing randomness in the individual estimators. In [START_REF] Bourel | Aggregating density estimators: an empirical study[END_REF] the authors include randomness in the construction of the intermediate histograms using then different aggregation schemes: AggregHist, using simple aggregation, BagHist using Bagging or StackHist using Stacking. Mathematically well sound, these approaches were explored thorough empirical simulation without theoretical framework. The Random Average Shifted Histogram, RASH [START_REF] Bourel | Random average shifted histograms[END_REF] is constructed as the mean average of different histograms, each one constructed over a random translated grid of the initial breakpoints of the initial histogram. RASH is show to be consistent and perform well in comparison to the precedent aggregation schemes.

In this work, we explore the contribution of Bagging to the density estimation task (like for BagHist) where the intermediate estimators are either histograms, frequency polygons or kernel density estimators. This article is organized as follow. Section 2 introduces notation and reviews the framework context we need about density estimators. In Section 3 we present our three methods: BagHist, BagFP and BagKDE and establish a theoretical result about their consistency. The results of comprehensive simulations are the object of Section 4. For this, we use several target densities available on literature to evaluate the performance of our estimators and comparing them to classical density estimators. We also explore the construction of pointwise confidence intervals as an approximation of a confidence band using the bootstrap procedure. The work concludes with a discussion in Section 5.

Some density estimators

We look at the principal results using different density estimators. While for the histogram and the kernel density estimator these results are quite popular we give some detail in order to fix notation. For a detailed exposition see [START_REF] Scott | Multivariate Density Estimation: Theory, Practice, and Visualization[END_REF].

In all cases, the starting point is an independent and identically distributed sample x 1 , . . . , x n of a real random variable with density f .

• Histogram. Let B j = [jh, (j + 1)h] be a set of intervals defined over the support of f , and h → 0 as n → +∞.

The ordinary histogram is defined as:

f h n (x) = 1 nh n i=1 j∈Z 1 B j (x i )1 B j (x) = ν j nh (1)
where ν j ∼ Bin(n, p j ), with p j = B j f (t)dt, is the number of observations of the sample that fall in bin B j . If x ∈ B j , we have that:

E f h n (x) = 1 nh n i=1 j∈Z P(x i ∈ B j )1 B j (x) = 1 nh nP(x i ∈ B j ) = p j h = f (ξ j ) Var f h n (x) = 1 n 2 h 2 Var(ν j ) = p j (1 -p j ) nh 2 ≤ p j nh 2 = f (ξ j ) nh
for some ξ j ∈ B j . For a x ∈ B j a fixed point, when h → 0, n → ∞ and nh → ∞, we get the classical properties for the histogram:

E f h n (x) → f (x), Var f h n (x) → 0.
and, moreover, if f is locally Lipschitz the histogram is mean square consistent, i.e MSE f h n (x) = Bias 2 f h n (x) + Var f h n (x) → 0. It is widely used in many fields, because of its computational simplicity. The histogram depends on two parameters: the bin width h and an origin x 0 to fix the grid. There is a huge literature that proposes several optimal choices for h using different criteria. If we suppose that the underlying density f is Gaussian, it can be shown (see [START_REF] Scott | On optimal and data-based histograms[END_REF]) that an optimal choice for h is of order n -1/3 . With this value the histogram has a rate of convergence of order n -2/3 with respect to the Mean Integrated Squared Error (MISE).

• Frequency Polygon. Frequency polygons are constructed on top of histograms connecting with straight lines the midpoint of two consecutive bin values. The expression of the frequency polygon for an

x ∈ B j = [(j -1/2)h, (j + 1/2)h] is f fp n (x) = 1 2 + j - x h ν j nh + 1 2 -j + x h ν j+1 nh (2)
The frequency polygon was deeply studied in [START_REF] Scott | Frequency polygons: Theory and application[END_REF]. With respect to the histogram, it presents the the advantages of being continuous and smooth. Under weak conditions, an optimal choice for h is of order n -1/5 and it achieves a rate of convergence of order order n -4/5 with respect to the MISE ( [START_REF] Scott | Frequency polygons: Theory and application[END_REF]).

• Kernel Density Estimators. A Kernel Density Estimator, KDE, is a function defined by

f kde n (x) = 1 nh n i=1 K x -x i h (3)
for all x ∈ R where K is a kernel function, i.e a non-negative, symmetric and unimodal function such that K(u) du = 1. Parameter h is called the bandwidth of the estimator f kde n , who inherits all the mathematical properties of K. The function K indicates the weight that observation x i has in the estimation of x: observations close to x are weighted more important. It can be shown that

E f kde n (x) = f (x) + f (x) 2 µ 2 (K)h 2 + o(h 2 ) Var f kde n (x) = ||K|| 2 2 f (x) nh + o 1 nh
where µ r (K) = u r K(u) du. When h → 0 and nh → ∞, we get the classical properties as for the histogram:

E f kde n (x) → f (x), Var f kde n (x) → 0.
As for the histogram, if h is large the variance decreases but the bias is large. On the other hand, if h is small, the bias is small but the variance is large. The optimal rate of convergence of KDE is of order n -4/5 as for the frequency polygon.

Bagging of estimators

The bootstrap method was introduced in [START_REF] Efron | Bootstrap methods: Another look at the jackknife[END_REF] and have the purpose of doing statistical inference using resamples of the original set of data. More precisely, if we have a data set L = {x 1 , x 2 , . . . , x n } with distribution F , the non parametric bootstrap procedure consists to draw, with replacement, a new sample L * = {x * 1 , x * 2 , . . . , x * n } from L of the same size. Then, the sample L * has a distribution F n , the empirical distribution of L. The main idea is that the sample L * is to the original sample L what the sample L is to the population, so the method treats the empirical of a distribution of sample data as the true distribution.

Reiterating this procedure several times and obtaining many bootstrap samples is a cornerstone to the construction of several bootstrap based approaches.

The bootstrap has three big applications: bias correction, construction of confidence interval and hypothesis testing [START_REF] Efron | An Introduction to the Bootstrap[END_REF].

However, the use of bootstrap in nonparametric density estimation requires some caution, particularly concerning the bias of the estimators. In our setting, let f be a nonparametric density estimator for f obtained from the sample L.

Now we draw a bootstrap sample L * of L.

In the case of the kernel density estimator, the estimation over L * is without bias,

E f * (x)|L = 1 nh n i=1 E K x -x * i h = 1 h E K x -x * i h = 1 h y∈R x * i K x -y h P(x * i = y) = 1 h n j=1 K x -x j h 1 n = 1 nh n j=1 K x -x j h = f (x).
This simple results holds also for the histogram and for the frequency polygon too (see demonstration of theorem 1), that is E( f * (x)|L) = f . The fact is treated in [START_REF] Hall | The Bootstrap and Edgeworth Expansion[END_REF] in this terms: "far from accurately estimating the substantial bias of f , the bootstrap sets the bias of this kind of density estimator equal to zero". In supervised learning, the main application of bootstrap is definitely the Bagging. It is a parallel aggregation method of individual entities that at each step b draw a bootstrap sample L * b of the original sample L and compute an estimator (a classifier in classification or a predictor in regression) over L * b . For an input x, the output of the Bagging method is the average in regression or the majority rule in classification of the intermediate estimators at x. We follow this aggregation strategy to construct new density estimators of a density function f . Our procedures run as follows (Figure 1): Let L = {x 1 , . . . , x n } be a sample with unknown distribution F admitting a density f . Also, considerer f n a density estimator evaluated in L. For b ∈ 1, . . . , B:

1. obtain L * b = {x * 1 , . . . , x * n } a bootstrap sample from L;
2. construct f * b the density estimator obtained over this bootstrap sample. In particular the bandwidth h is calculated over L * b Output: The final estimator is the simple pointwise average of the individual estimators i.e.

f * (x) = 1 B B b=1 f * b (x)

L 2 consistency of the Bagging of estimators

Here we prove L 2 consistency, for three estimators BagHist, BagFP and BagKDE which correspond to bagging of histograms, bagging of polygon frequencies and bagging of kernel density estimators proving that, since h → 0, n → +∞ and nh → +∞, for all x in bin B j (in case of histogram or frequency polygon) or for all x ∈ R (in case of kernel density estimator):

E[( f baghist n (x) -f (x)) 2 ] → 0, E[( f bagfp n (x) -f (x)) 2 ] → 0, E[( f bagkde n (x) -f (x)) 2 ] → 0.
Theorem 1 BagHist, BagFP and BagKDE are L 2 -consistent.

Proof 1 We will give a global proof, inspired by [START_REF] Hall | The Bootstrap and Edgeworth Expansion[END_REF] and [START_REF] Scott | Multivariate Density Estimation: Theory, Practice, and Visualization[END_REF] to encompass the different methods, because the demonstration for all these estimators follows the same steps. We have to compute:

(1) for the expectation

E( f (x)) = E E[ f (x)|L]
(2) and to calculate the variance we use the decomposition

Var( f (x)) = E(Var( f (x)|L)) (A) + Var(E( f (x)|L)) (B)
Without loss of generality, in some calculation for histogram or frequency polygon, and with the aim to simplify notations, we will assume that x ∈ B 0 .

(1)

• BagHist. If x ∈ B j = [jh, (j + 1)h] , then we have

E[ f baghist n (x)|L] =E 1 B B b=1 f * b (x)|L = E f * b (x)|L = E ν * j nh = ν j nh = f h n (x)
and if f is locally Lipschitz

E(E[ f (x)|L]) -f (x) = p j h -f (x) = hf (ξ j ) h -f (x) ≤γ j |ξ j -x| ≤ γ j h → 0 • BagFP. If x ∈ B j = (j -1 2 )h, (j + 1 2 )h , then we have E[ f bagfp n (x)|L] =E 1 B B b=1 f * b (x)|L = E f * b (x)|L =E 1 2 + j - x h ν * j nh + 1 2 -j + x h ν * j+1 nh |L = 1 2 + j - x h ν j nh + 1 2 -j + x h ν j+1 nh = f fp n (x)
and, if f has second derivative, if x ∈ B 0 :

E E[ f bagfp n (x)|L] -f (x) ≈ f (ξ 0 )(h 2 -3x 2 ) ≤ 4h 2 f (ξ 0 ) → 0 • BagKDE.
For BagKDE we have:

E[ f bagkde n (x)|L] =E 1 B B b=1 f * b (x)|L = E f * b (x)|L =E 1 nh n i=1 K x -x * i h = f kde n
and this is well known that

E E[ f bagkde n (x)|L] -f (x) → 0
(2) For variance we use the well known formula defined above.

• BagHist.

(A) Because of the independence and identical distribution of the bootstrap samples, if x ∈ B 0 :

Var[ f baghist n (x)|L] =Var 1 B B b=1 f * b (x)|L = 1 B Var f * b (x)|L = 1 B np * 0 (1 -p * 0 ) (nh) 2
where p * 0 is equal to ν 0 n . Taking expectation over L we have

1 Bnh 2 E (p * 0 (1 -p * 0 )) = 1 Bnh 2 E ν 0 n -Var ν 0 n -E ν 0 n 2 = 1 Bnh 2 p 0 - 1 n p 0 - 1 n p 2 0 -p 2 0 = 1 Bnh 2 hf (ξ 0 ) - h n f (ξ 0 ) - h 2 n f (ξ 0 ) 2 -h 2 f (ξ 0 ) 2 = 1 Bnh f (ξ 0 ) - 1 Bn 2 h f (ξ 0 ) - 1 Bn 2 f (ξ 0 ) 2 - 1 Bn f (ξ 0 ) 2
which tends to 0 as n → ∞, nh → ∞ and h → 0.

(B) Var E( f baghist n (x)|L) = Var( f h n (x)) ≤ p 0 nh 2 = f (ξ 0 ) nh → 0 • BagFP.
(A) Because of the independence and identical distribution of the bootstrap samples:

Var[ f bagfp n (x)|L] = Var 1 B B b=1 f * b (x)|L = 1 B Var f * b (x)|L = 1 B 1 2 - x h Var( f * 0 ) + 1 2 + x h Var( f * 1 ) + 2 B 1 4 - x 2 h 2 Cov( f * 0 , f * 1 )
where f * 0 and f * 1 are the histogram estimations over [-h, 0] and [0, h] respectively. As Var(

f * 0 ) = np * o (1-p * o ) n 2 h 2
, then taking expectation:

E(Var( f * 0 )) =E n ν 0 n (1 -ν 0 n ) n 2 h 2 = 1 n 2 h 2 E ν 0 1 - ν 0 n = 1 n 2 h 2 E(ν 0 ) - 1 n E(ν 2 0 ) = 1 n 2 h 2 E(ν 0 ) - 1 n (Var(ν 0 ) + [E(ν 0 )] 2 ) = 1 n 2 h 2 np 0 - 1 n (np 0 (1 -p 0 ) + (np 0 ) 2 ) = nhf (ξ 0 ) n 2 h 2 - hf (ξ 0 ) n 2 h 2 + h 2 f (ξ 0 ) 2 n 2 h 2 - nh 2 f (ξ 0 ) 2 n 2 h 2 → 0
In the same way E(Var(

f * 1 )) → 0. As Cov( f * 0 , f * 1 ) = -np * 0 p * 1 n 2 h 2 ,
then taking expectation we have:

E Cov( f * 0 , f * 1 ) = 1 nh 2 E ν 0 n ν 1 n ≤ 1 n 3 h 2 E(ν 2 0 )E(ν 2 1 ) = 1 n 3 h 2 [np 0 (1 -p 0 ) + p 2 0 ][np 1 (1 -p 1 ) + p 2 1 ] = 1 n 3 h 2 (nhf (ξ 0 ) -nh 2 f (ξ 0 ) 2 + h 2 f (ξ 0 ) 2 )(nhf (ξ 1 ) -nh 2 f (ξ 1 ) 2 + h 2 f (ξ 1 ) 2 ) = 1 n 3 h 2 (n 2 h 2 f (ξ 0 )f (ξ 1 ) -n 2 h 3 f (ξ 0 )f (ξ 1 ) 2 + nh 3 f (ξ 0 )f (ξ 1 ) 2 + 1 n 3 h 2 (n 2 h 4 f (ξ 0 ) 2 f (ξ 1 ) 2 -n 2 h 3 f (ξ 0 ) 2 f (ξ 1 ) -nh 4 f (ξ 0 ) 2 f (ξ 1 ) 2 ) + 1 n 3 h 2 (nh 3 f (ξ 0 ) 2 f (ξ 1 ) -nh 4 f (ξ 0 ) 2 f (ξ 1 ) 2 + h 4 f (ξ 0 ) 2 f (ξ 1 ) 2 ) = 1 n f (ξ 0 )f (ξ 1 ) - h n f (ξ 0 )f (ξ 1 ) 2 + h n 2 f (ξ 0 )f (ξ 1 ) 2 + h 2 n f (ξ 0 ) 2 f (ξ 1 ) 2 - h n f (ξ 0 ) 2 f (ξ 1 ) - h 2 n 2 f (ξ 0 ) 2 f (ξ 1 ) 2 + h n 2 f (ξ 0 ) 2 f (ξ 1 ) - h 2 n 2 f (ξ 0 ) 2 f (ξ 1 ) 2 + h 2 n 3 f (ξ 0 ) 2 f (ξ 1 ) 2 → 0 So we conclude that E Var[ f bagfp n (x)|L] → 0 (B)
We recall from [START_REF] Scott | Multivariate Density Estimation: Theory, Practice, and Visualization[END_REF] pag. 103 that

Var(E( f bagfp n (x)|L)) = Var( f fp n (x)) = 2x 2 nh 3 + 1 2nh f (ξ 0 )- f (ξ 0 ) 2 n +o 1 n
Then if nh → ∞ and n → +∞:

|Var(E( f bagfp n (x)|L))| → 0 • BagKDE.
(A) Because of the independence and identical distribution of the bootstrap samples:

Var[ f bagkde n (x)|L] = Var 1 B B b=1 f * b (x)|L = 1 B Var f * b (x)|L So we compute Var f * b (x)|L : Var f * b (x)|L = 1 n n i=1 1 (nh) 2 K 2 x -x * i h (a) - 1 n 2 1 nh n i=1 K x -x * i h 2 (b)
and therefore Var f * b (x)|L ≤ (a) + (b). Taking expectation:

(a)

E 1 n 3 h 2 n i=1 K 2 x -x * i h = 1 n 3 h 2 n i=1 E K 2 x -x * i h = 1 n 2 h 2 n j=1 K 2 x -x ij h P(x * i = x ij ) = 1 n 3 h 2 n j=1 K 2 x -x ij h ≤ 1 n 3 h 2 n j=1 C = C (nh) 2 → 0 because since K is bounded, K 2 also. (b) E   1 n 2 1 nh n i=1 K x -x * i h 2   = 1 n 4 h 2 E n i=1 K x -x * i h 2 = 1 n 4 h 2 n j=1 n i=1 K x -x * i h 2 P(x * i = x ij ) ≤ 1 n 5 (nh) 2 n j=1 n i=1 C 2 = C 2 n 2 (nh) 2 → 0 So we conclude that E Var[ f bagkde n (x)|L] → 0 (B)

Experiments

We describe in this section a series of numerical experiments aiming to show the practical performance of the bagged versions of the classical density estimators. First, we obtain a numerical estimate of the MISE on simulated data sets created following baseline densities.

The impact of the aggregation is analyzed. We also use the bootstrapped version of the density estimators to construct a confidence interval for it and study its performance.

Simulations

Among the numerous possibilities of univariate densities, we choose eight simulation models partially following the work of [START_REF] Bourel | Random average shifted histograms[END_REF]. This choice presents a different degree of difficulty related to the number of modes, asymmetry, tail behavior and regularity. We denote them by M1 to M8. Their definition is the object of Table 1 and Figure 2 shows a graphical display of the densities. The notation N (µ, σ Table 1: Simulated univariate densities.

At each replication we draw two datasets following each density. The first one is used for estimation purposes while the second one is left-out for evaluation (either MISE or empirical covering). All the simulations are done with the R software ([? ]).

Quality of the estimation

We compare density estimators of different nature. On one hand side we use three individual estimators: histograms (H), frequency polygons (FP) and kernel density estimators (KDE), on the other hand, their bagged versions, respectively BagHist, BagFP and BagKDE. Also we include the RASH estimator. We use cross validation to calibrate the bandwidth h at each step of all the intermediate estimation methods. An alternative would be to use maximum likelihood as in [START_REF] Bourel | Random average shifted histograms[END_REF]. In our framework cross validation has, in general, a better computation behavior. Also it is more general and may be used for example with dependent data as in time series. 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 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 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 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 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 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 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 Triangular Uniform 0-1 Mix2 Mix3 Figure 3 represents the dependence of MISE on the sample size n for the different combinations of densities and estimators. Each point represents the average of M = 100 times the MISE×100 of the method using B = 200 intermediate estimators for the aggregating methods. Notice that these plots are in log-log scale which is useful to highlight the convergence behaviour.Individual values of these plots are presented in Appendix A.

Normal

Let us comment these plots. First, the adjusted lines are of relative good quality since the points corresponding to each combination density-estimator are almost aligned. Remember that each point is the mean average of M = 100 replicates and so the inner replicate variability is reduced even for a few points determining each line. However, in some few cases the quality of the fit is quite poor. Now, for each panel most of the adjusted lines are almost parallels which means the methods share a similar convergence behavior. Comparing the four simple densities (leftmost panels) and the four mixtures (rightmost panels) a difference in the behaviour seems to appear, at least if one looks at some bagged version as for instance the BagFP. Second, if one compares each individual estimator (H, FP and KDE) with their bagged versions, the latter success to reduce the MISE in most of the situations. An important exception if BagHist which produces almost always worst results than the intermediate estimator, i.e. the histogram. However, the asymptotic behaviour is such that with large sample sizes it is able to catch the quality of the histogram and even overwhelm it on the four mixtures. On regular targets, KDE (or at least its bagged version) shows a very competitive performance. However, in presence of multi-modality they loss relative competitiveness with large sample sizes. The fact that the results are not entirely satisfactory for the bagged version of KDE may be because KDE it is already a good and stable density estimator (more stable in any case than histogram) and, according with [START_REF] Breiman | Bagging predictors[END_REF], bagging kernel density estimators may be degrade the performance of this stable procedure.

Reduction of MISE due to aggregation

We concentrate now on aggregating methods. A natural matter to look at is the quality of the aggregation as the number of bootstrap samples increases. For this, we examine the MISE of the bagged versions for a range of increasing bootstrap samples. We replicate M = 100 times each combination of density simulation to construct the different curves. The result of experiments are presented in Figure 4 in a log-log scale with n = 500 observations.

Globally we observe that MISE decreases monotonically with increasing values of B until some point between 20 and 50 bootstrap samples after which more samples does not produce further enhancement.

Variability bands

A natural by-product of bootstrap samples is the construction of confidence bands. For some level α, one wants to estimate the quantities ln (x) and ûn (x) that verify

P{ ln (x) ≤ f (x) ≤ ûn (x)} ≥ 1 -α, ∀ x
that is, the quantities are the borders of an interval that covers at the true density f (x) at some confidence level (1 -α) × 100%. We tackle here its construction for the density estimator. Generally, a confidence band for f is centered over an estimator fn of f and has the form fn (x) ± c σn (x) for all x, with c > 0. However, since nonparametric density estimators are biased, the usual construction does not yields on a really a confidence band for f . Indeed, for a fixed x, due to the bias E fn (x) -f (x), it is not easy to derive a confidence interval using the pivotal quantity fn(x)-f (x) σ(x)

. So, the interval is usually centered at f n (x) = E fn (x) instead of being around f (x). For this reason, these confidence bands are often called variability bands.

We describe two popular constructions to compare with our procedure.

1. Variability Band for histograms. Under mild conditions [START_REF] Wasserman | All of Nonparametric Statistics[END_REF] the histogram estimator f h n (x) is approximately unbiased for the target density f (x). But the approximate variance is f (x)/(nh) where h = 1/m is the inverse of the number of bins m. Its dependence on the unknown target is an obstacle. To circumvent it, [START_REF] Scott | Multivariate Density Estimation: Theory, Practice, and Visualization[END_REF] looks at Var f h n (x) which is approximately 1/(4nh) and thus independent of f (x). We

define f h n = E[ f h n (x)
] as the target and as we say before the confidence band will not take account of the bias but only of the variability of the estimator. Then, using a normal approximation it is easy to show that [22, p. 130]:

l n (x) = max f h n (x) -c, 0 2 , u n (x) = f h n (x) + c 2 where c = z α/(2m) 2
m n give an approximate variability band for f h n at (1 -α) × 100% of confidence.

2. Variability Band for KDE. As we have shown with the histogram, since the variance of f kde n (x) also involves the true density f , it is more suitable to use the square root (see [START_REF] Bowman | Applied Smoothing Techniques for Data Analysis: The Kernel Approach with S-Plus Illustrations[END_REF]). In the case of the kernel density estimator Var

f kde n (x) ≈ ||K|| 2 2
4nhn and again does not depend on the true unknown density f . On this square root scale, for a fixed point x we consider the interval that back to the original scale is given by

l n (x) = f kde n (x) - ||K|| 2 √ 4nh n 2 , u n (x) = f kde n (x) + ||K|| 2 √ 4nh n 2 .
As we said before this is not a confidence band for the true density f , because of the bias so we will talk about a variability band.

Bootstrap based confidence interval and resulting tube

The bootstrapped sample induces a distribution that can be used to asses the variability of the estimator. Indeed, the simple superposition of the individual estimators (histogram, frequency polygon, kernel density estimator) gives a coarse idea of the uncertainty around the aggregate estimator. More the scatter of individual individual density estimators is dispersed, higher is the variance of the estimator. For the concrete construction of the confidence band we first fix the abscissa x ∈ R. Then we consider the set of bootstrapped density estimators evaluated at that point, i.e. { f * 1 (x), . . . , f * B (x)}. This set is a collection of B univariate measures. Note that the bagged estimation is the mean average of this collection. Then, a (1-α)×100% confidence interval can be obtained by considering the empirical quantiles at α/2 and 1 -α/2 of this ensemble. In Figure 5 we show the different constructions generated by these methods for a standard normal target. Each panel correspond to one of the three intermediate density estimators (Histogram, frequency polygons and KDE). At each time, the bootstraped construction for the confidence interval is represented. For the histogram and KDE we also draw the constructions detailed above. The confidence interval for the histogram succeeds to cover the true density but produces an arguably too large band. For KDE the confidence for the intermediate estimator is very good, and the one obtained by bootstrap sampling presents a relatively high variability.

We compare the alternative constructions of the confidence band using two metrics. The aim is to obtain the narrowest band that warranties a given nominal coverage. For this, we consider the empirical coverage of the bands and its mean width. Let us call ln (t i ) and ûn (t i ) the lower and upper bounds of the confidence bands, evaluated at points t i , i = 1, 2, . . . , N . Then, we call the empirical mean coverage of the target f (x) the quantity 1

N N i=1 I { ln(ti)≤f (t i )≤ûn(t i )} ,
where I A is the indicator function of the set A. The mean width of the interval is defined by

1 N N i=1 ûn (t i ) -ln (t i ) .
In this experiment we set the confidence level at 95% to construct the variability bands. We give as reference the variability band constructed through the kernel density estimator as explained before (we denote this method as KDE-sm). The construction for the individual histograms produces too larges bands which always cover the true density. For this reason they are not presented in Table 2.

In general, the bands cover reasonably well the simplest densities (on the top rows) having more difficulties with the more exotic density models (on the bottom rows). Particularly, the last three densities in the table are too difficult targets producing very low empirical coverings. If we look at the methods, the construction using bootstraped histograms is the best one among the boostrap based ones, and give a fair competitor to KDE-sm.

Coverage

Mean

Conclusions

In this work we present three univariate density estimators obtained by aggregation such as in Bagging. For each method, the intermediate estimators are histograms, frequency polygons or kernel density estimators. We prove the L 2 consistency of the three estimators and do several simulations over densities with different characteristics. Also, we bring a way to compute a kind of confidence band, which is more close to a point wise variability band in the sense that the authors who studied on this subject give. This construction needs a deeper study to be able to draw more conclusive conclusions about it. Another clue for future work is to investigate the natural extension of considering the bagged construction over multivariate densities.

A Additional results

Quality of the estimation

For sake of completeness we present in this appendix the individual values of Figure 3. In the following tables (Tables 3 to 7

Figure 1 :

 1 Figure 1: Bagging of density estimators To obtain the bagged histogram (BagHist) estimator f baghist n we simply use at each step b, histograms f h n defined in (1) as f * b . Analogously for bagged frequency polygons (BagFP) f bagfp n

  It is a well known result that Var E[ f bagkde n (x)|L] = Var( f kde n ) → 0 So, with the usual assumption of n → ∞, h → 0, nh → ∞ this implies L 2 convergence for f baghist n , f bagfp n and f bagkde n .

  2 ) is used to refer to a normal distribution with mean equal to µ and variance equal to σ 2 , U[a, b] is the uniform density over the support [a, b], and χ 2 ν is a Chi squared density with ν degrees of freedom. Models 3, 4, 7 and 8 are mixtures of densities. Models 2 and 7 are asymmetrical. Model Description (M1) : Normal Standard Standard Gaussian density N (0

Figure 2 :

 2 Figure 2: Densities used for the simulations.

Figure 3 :

 3 Figure 3: MISE by estimation method for the six simulation data sets in scale log.

Figure 4 :

 4 Figure 4: MISE error vs number of aggregates, n=500, M=100, B=200 in log-log scale.

Figure 5 :

 5 Figure 5: Variability bands for the estimation of a standard normal target (red) based on Histograms (left), Frequency polygons (middle) and KDE (right). Confidence intervals are either constructed for individual estimators (in blue) when available or constructed using bootstrap sampling (in green).

Table 2 :

 2 Mean empirical coverage and mean interval widths for the densities and estimators considered.

	width

Table 3 :

 3 ), values are 100×MISE obtained as mean average over 100 replicates. At each line, best results are shown in blue. MISE for sample size n = 50

		Hist	FP	KDE BagHist BagFP BagKDE RASH
	Normal	1.1447	0.6171 0.3181	10.3243 0.2487	1.5195 0.6025
	Chi10	0.0696	0.0374 0.0254	0.5635 0.0181	0.0919 0.0334
	Mix1	3.8260	2.1742 1.3801	27.4794 1.8937	3.0423 1.9464
	Claw	4.9526	3.6531 2.7909	17.8066 2.2902	4.1301 3.4380
	Triangular	7.1220	4.4168 2.2008	60.1150 1.7936	8.6147 4.0222
	Uniform 0-1 17.2112 11.4841 6.8795 175.0100 6.9725	21.2978 8.5623
	Mix2	6.2811	6.0856 6.3647	16.3470 5.8546	5.7913 5.9907
	Mix3	4.9186	3.0408 1.8700	22.3239 2.0799	2.5847 2.3357
		Hist	FP	KDE BagHist BagFP BagKDE RASH
	Normal 0.7098 0.4016 0.2665	4.9699 0.1665	0.9510 0.3908
	Chi10 0.0422 0.0215 0.0144	0.2689 0.0089	0.0549 0.0198
	Mix1 2.5451 1.3417 0.8085 12.7401 2.5901	1.9115 1.1964
	Claw 2.6535 1.9604 1.6735	8.1400 2.1407	2.3149 1.8267
	Triangular 5.0242 3.1358 1.2952 33.9450 0.8356	5.8185 2.8512
	Uniform 0-1 7.9732 5.2470 5.2962 87.8440 3.8460	14.8139 4.1718
	Mix2 5.8523 5.5056 6.0074	8.9022 5.9432	4.5841 5.2529
	Mix3 2.4888 1.6908 1.4754 11.2250 2.1891	1.8702 1.1925

Table 4 :

 4 MISE for sample size n = 100

	Hist	FP	KDE BagHist BagFP BagKDE RASH
	Normal 0.4077 0.2059 0.1394	2.4724 0.0827	0.5840 0.1937
	Chi10 0.0242 0.0121 0.0084	0.1246 0.0075	0.0305 0.0111
	Mix1 1.6231 0.8703 0.5597	6.1511 3.2618	1.2028 0.8753
	Claw 1.7414 1.1899 0.9338	3.7563 2.0185	1.4601 1.0360
	Triangular 2.4070 1.4362 0.7807 16.6818 0.5730	3.4612 1.3875
	Uniform 0-1 3.7493 2.5354 3.9036 44.4120 2.7891	10.1004 2.1922
	Mix2 4.2920 3.2531 5.7958	4.3796 5.9281	3.5878 2.6659
	Mix3 1.6959 1.2517 1.0518	5.5502 2.3213	1.2049 0.8032

Table 5 :

 5 MISE for sample size n = 200

Acknowledgements

We would like to thank project ECOS-2014 Aprendizaje Automático para la Modelización y el Análisis de Recursos Naturales, n o U14E02, the LIA-IFUM and the ANII -Uruguay for their financial support.


Robert E Alvarez

Invertibility of spectral x-ray data with pileup-two dimension-two spectrum case

Keywords: Jacobian, photon counting, dead time, pileup, spectral x-ray, dual energy, energy selective

In the Alvarez-Macovski[1] method, the line integrals of the x-ray basis set coefficients are computed from measurements with multiple spectra. An important question is whether the transformation from measurements to line integrals is invertible. This paper presents a proof that for a system with two spectra and a photon counting detector, pileup does not affect the invertibility of the system. If the system is invertible with no pileup, it will remain invertible with pileup although the reduced Jacobian may lead to increased noise.

I. INTRODUCTION

In the Alvarez-Macovski [START_REF] Alvarez | Energy-selective reconstructions in X-ray computerized tomography[END_REF] method, the line integrals of the x-ray basis set coefficients are computed from measurements with multiple spectra. The introduction of photon counting detectors into medical x-ray imaging [START_REF] Taguchi | Vision 20/20: Single photon counting x-ray detectors in medical imaging[END_REF] gives the possibility of providing the spectral data by pulse height analysis (PHA). These detectors, however, have multiple defects [START_REF] Overdick | Towards direct conversion detectors for medical imaging with X-rays[END_REF] that affect the information they provide. In some cases, the defects have been found to produce sharply increased noise in the estimates of the line integrals at specific object attenuation [START_REF] Alvarez | Non-invertibility of spectral x-ray photon counting data with pileup[END_REF]. The increased noise may indicate potential non-invertibility of the transformation between the spectral measurements and the line integrals. Therefore, it is important to develop mathematical descriptions of the invertibility of the transformation.

This paper is a step towards this mathematical description. It presents a proof that measurements with two spectra and a photon counting detector with pileup do not affect the invertibility of the system. If the system is invertible with no pileup, it will remain invertible with pileup although the reduced Jacobian may lead to increased noise. An example of the system analyzed is measurements with a photon counting detector of the transmitted flux with two different x-ray tube voltages. Note that the results of this paper are not applicable to pulse height analysis (PHA) since pileup changes the effective spectra of the bins [START_REF] Alvarez | Signal to noise ratio of energy selective x-ray photon counting systems with pileup[END_REF]. My recent paper [START_REF] Alvarez | Non-invertibility of spectral x-ray photon counting data with pileup[END_REF] gives an example of a three bin PHA system that becomes non-invertible for high pileup.

This paper addresses only the invertibility with deterministic, non-noisy measurements although it describes the conditioning of the system, which also affects noise.

aprendtech.com, ralvarez@aprendtech.com

II. THE ALVAREZ-MACOVSKI METHOD

For biological materials we can approximate the xray attenuation coefficient µ(r, E) accurately with a two function basis set [START_REF]Dimensionality and noise in energy selective x-ray imaging[END_REF]

µ(r, E) = a 1 (r)f 1 (E) + a 2 (r)f 2 (E). (1)
In this equation, a i (r) are the basis set coefficients and f i (E) are the basis functions, i = 1 . . . 2. As implied by the notation, the coefficients a i (r) are functions only of the position r within the object and the functions f i (E) depend only on the x-ray energy E. If there is a high atomic number contrast agent, we need to extend the basis set to higher dimensions. Neglecting scatter, the expected value of the number of transmitted photons λ k with an effective measurement spectrum S k (E) is

λ k = ˆSk (E)e -´L µ(r,E)dr dE (2)
where the line integral in the exponent is on a line L from the x-ray source to the detector. Using the decomposition, Eq. 1, the line integral in Eq.

is

ˆL µ (r, E) dr

= A 1 f 1 (E) + A 2 f 2 (E). (3)
where A i = ´ai (r) dr, i = 1 . . . 2 are the line integrals of the basis set coefficients. Summarizing the A i as the components of the A-vector, A, and the basis functions at energy E as a vector

f (E) = [f 1 (E), f 2 (E)],
we can write the line integral as the inner product of A and f (E),

ˆL µ (r, E) dr = A • f (E). (4)
The measurements are summarized by a vector, N, whose components are the expected photon counts with each effective spectrum. Since the body transmission is exponential in A, we can approximately linearize the measurements by taking logarithms. The results is the log measurement vector L = -log(N/N 0 ), where N 0 is the expected value of the measurements with no object in the beam and the division means that corresponding members of the vectors are divided. Equations 2 define a relationship between A and the expected value measurement vector, L(A). The invertibility of this transformation is the subject of this paper.

III. INVERTIBILITY WITHOUT PILEUP-SIMPLE CASES

Before discussing the invertibility in general, two simple but important cases will be discussed.

A. Delta function spectra

The first case is two monoenergetic spectra with energies E 1 and E 2 . In this case, L i = A • f (E i ), i = 1, 2 and the transformation can be written L(A) = MA where M is a matrix with coefficients

M = f 1 (E 1 ) f 2 (E 1 ) f 1 (E 2 ) f 2 (E 2 )
The transformation is invertible if the determinant of M is not equal to zero:

f 1 (E 1 ) f 2 (E 1 ) f 1 (E 2 ) f 2 (E 2 ) = 0 (5)
That is if

f 1 (E 1 ) f 2 (E 1 ) = f 1 (E 2 ) f 2 (E 2 ) (6)
The bottom graph of Fig. 1 shows the ratio of the Compton scattering/photoelectric cross-sections basis set functions [START_REF] Alvarez | Energy-selective reconstructions in X-ray computerized tomography[END_REF] versus energy in the medical diagnostic region. Note that the ratio is monotonically decreasing. Thus, the condition in equation ( 6) will be true if the two energies are different. The M matrix for any other basis set that spans the space of attenuation coefficients will be M = TM where T is an invertible matrix so it has a nonzero determinant. Since these are square matrices, det (M ) = det (T) det (M) and if det (M) = 0 then det (M ) = 0 so this result is true for any valid basis set.

B. Single material object

Another useful special case is when the object is known to consist of a single material with basis set coefficient vector a. In this case, a single photon count measurement suffices to determine the material thickness. With a single material of thickness t, the A vector is A = at and the expected number of photons is a function of only

t λ(t) = ˆS1 (E)exp[-µ(E)t]dE.
The derivative of λ with respect to t is always negative

λ = dλ dt = -´µ(E)S 1 (E)exp[-µ(E)t]dE < 0. ( 7
)
because the integrand is always greater than zero. Since λ(t) is monotonically decreasing, it can be inverted to compute t.

IV. INVERTIBILITY FOR THE GENERAL CASE WITH NO

PILEUP

For the general case, the following theorem is useful (Fulks 1978[7] page 284):

Let F be a continuously differentiable mapping defined on an open region D in E2, with range R in E2 , and lets its Jacobian be never zero in D. Suppose further that C is a simple closed curve that, together with its interior (recall the Jordan curve theorem), lies in D, and that F is one-to-one on C. Then the image T of C is a simple closed curve that, together with its interior, lies in R. Furthermore, F is one-to-one on the closed region consisting of C and its interior, so that the inverse transformation can be defined on the closed region consisting of T and its interior. which I paraphrase as

If the Jacobian of a continuously differentiable two dimensional mapping is nonzero throughout an open region D and if the mapping is one to one on a simple closed curve C which lies in D, then the mapping is one to one on C and its interior. The first quadrant will be used as the region with the closed curve C consisting of segments along the positive axes and a circle joining the ends of the segments, as shown in Figure 2. This is a region of theoretical and practical importance because a basis set consisting of the attenuation coefficients of the calibration materials is usually used, Since only positive equivalent thicknesses of the calibration materials can be used, this region must contain all the measured values. .

In the theorem, the Jacobian is the determinant of the matrix of all the partial derivatives of the transformation. Since L is the logarithm of the measurements, the components of its Jacobian matrix are

M ij = ∂Li ∂Aj = -´fj (E)Si(E)e -A•f (E) dE ´Si(E)e -A•f (E) dE i, j = 1, 2 (8)
Note that by defining the normalized spectra

ŝi (E) = S i (E)e -A•f (E) ´Si (E)e -A•f (E) dE i = 1, 2 (9)
the M ij are

M ij = ˆfj (E)ŝ i (E)e -A•f (E) dE
That is, they are the effective values of the basis functions in the spectra transmitted through the object.

A. Application of the theorem to the dual energy transformation

To apply the first part of the theorem, the transformation must be shown to be invertible on a closed curve in the domain. The simple cases discussed Sec. III may be used for this proof.

The parts of the curve along the axes are special cases of the single material case. Each axis corresponds to different thicknesses of one of the basis materials if the attenuation coefficient of the calibration material is used as a basis function. If other basis functions used, the coordinates can be transformed to the attenuation coefficients of real materials and the proof applies in the transformed coordinates.

The circle of large radius is an approximation of the single energy case. For large radius, there will be high attenuation. With beam hardening, the transmitted spectrum with large attenuation therefore approaches the two monoenergetic spectra case where the energies are the maximum energies in the spectra. If the maximum energies are different, we can approach the known invertible monoenergetic case arbitrarily closely by making the radius larger and larger.

The remaining part of the theorem requires us to show that J (A) is non-zero inside C. This must be tested with individual spectra. The following sections show that, for the system studied, pileup does not affect this condition. If the Jacobian is non-zero without pileup, it will also be non-zero with pileup.

V. EXPECTED NUMBER OF PHOTONS RECORDED WITH

PILEUP

The response time of a photon counting detector is modeled using the dead time, τ , which is defined to be the minimum time between two photons that are recorded as separate events [START_REF] Knoll | Radiation Detection and Measurement[END_REF]. The dead time is an abstraction that combines the contributions to the response time of all the physical effects in the detector. In this model, the detector is assumed to start in a "live" state. With the arrival of a photon, the detector enters a separate state where it does not count additional photons. The nonparalyzable model will be used where the time in the separate state is assumed to be fixed and independent of the arrival of any other photons during the dead time. There is a second model commonly used, called paralyzable, where the arrival of photons extends the time in the non-counting state. Both models give similar recorded counts at low interaction rates but give different results at high rates where the probability of multiple interactions during the dead time becomes significant. Measurements by Taguchi et al. [START_REF] Taguchi | Modeling the performance of a photon counting x-ray detector for CT: energy response and pulse pileup effects[END_REF] indicate that for the detectors they studied the non-paralyzable model is more accurate at higher count rates. It also leads to simpler analytical results [START_REF] Yu | Mean and variance of single photon counting with deadtime[END_REF].

In a previous paper [START_REF] Alvarez | Signal to noise ratio of energy selective x-ray photon counting systems with pileup[END_REF], I used the central limit theorem of renewal processes [START_REF] Ross | Stochastic Processes[END_REF] to show that, for recording times much greater than the mean inter-event time, the expected value of recorded counts with pileup approaches

N rec = λ 1 + η . (10)
In this equation, λ is the expected number of photons incident on the detector during the measurement time and η is the expected number of photons arriving during the dead time τ . If ρ is the average rate of photon arrivals then η = ρτ . If T rec is the measurement time then ρ = λ /Trec. Defining b = τ /Trec, η = bλ and the expected recorded counts are

N rec = λ 1 + bλ . (11)

VI. INVERTIBILITY OF TWO BASIS FUNCTIONS-TWO SPECTRA CASE WITH PILEUP

In this section I give a proof that with photon counting detector measurements of the total number of transmitted photons with two spectra, pileup does not affect invertibility. An example would be making two sequential measurements of an object using an x-ray tube with different voltages. Note that this does not prove that any two spectrum measurement with pileup is invertible.

For example with two bin photon counting with PHA, pileup causes the recorded spectrum to change so the assumptions of this section would not be met.

The proof is analogous to the proof for the measurements with no pileup described in the previous Sec. III.

A. Invertibility with pileup on the contour C

First, I will show that if the transformation is invertible on the path in the first quadrant shown in Fig. 2 without pileup it is also invertible with pileup.

The proof for invertibility on the circular segment joining the segments on the axes is also applicable with pileup since, for large thicknesses, the count rate is very low so the pileup parameter η is essentially equal to zero and the pileup counts are the same as those without pileup.

Next we need to show that the data are invertible on the paths from the origin along the coordinates axes. The equation for the expected value of the recorded number of photons with pileup is

N rec = λ 1 + bλ . ( 12
)
Differentiating this equation along the axes

N rec = λ 1 + bλ - b (λ ) λ (1 + bλ) 2
where a prime denotes a derivative with respect to object thickness. Factoring the equation

N rec = λ 1+bλ 1 -bλ 1+bλ = λ (1+bλ) 2 . ( 13
)
From Eq. 7, for a single material λ is always less than zero. Since the denominator of Eq. 13 is always positive, the derivative of the recorded counts with pileup with is never equal to zero so the transformation of the recorded counts with pileup is invertible along the coordinate axes. This shows that if the transformation without pileup is invertible on C then the transformation with pileup is also invertible on the curve.

B. Jacobian inside C

The log measurement with pileup is

L rec,i = -log N rec,i N rec,i,0
and the Jacobian matrix with pileup has elements

M rec,ij = ∂Lrec,i ∂Aj = -1 Nrec,i ∂Nrec,i ∂Aj . (14)
From Eq. 11

N rec,i = λ i 1 + bλ i and from Eq. 13

∂N rec,i ∂A j = λ i (1 + bλ i ) 2 .
Substituting in Eq. 14

M rec,ij = 1 1 + bλ i λ i λ i (15)
From the definition in Eq. 8, the term in brackets in Eq. 15 is the element of the Jacobian matrix without pileup. Therefore,

M rec,ij = M 0,ij 1 + bλ i
Since for the dual energy case the matrices are 2 × 2, the determinant is .

J rec = M rec
where J 0 is the Jacobian without pileup. Since the terms in the denominator are always positive, the zero values of Jacobian determinant with pileup, if any, will occur at the same points as the Jacobian without pileup.

VII. DISCUSSION

The data acquisition model used is unrealistic-most systems with photon counting detectors would also use PHA. Nevertheless, it provides an example of invertibility with pileup.

As discussed in my previous paper [START_REF] Alvarez | Signal to noise ratio of energy selective x-ray photon counting systems with pileup[END_REF], due to the exponential probability distribution of x-ray photon inter-arrival times, pulse pileup is a fundamental effect. It will always be present in photon counting systems no matter how small the response time. Pileup has accurate analytical models and its analysis may lead to an understanding of invertibility with the other detector defects such as charge trapping and sharing, polarization, and incomplete photon energy deposition due to Compton scattering and K radiation escape [START_REF] Overdick | Towards direct conversion detectors for medical imaging with X-rays[END_REF].

VIII. CONCLUSION

A proof is given that pileup does not affect the invertibility of a two spectrum, photon counting data with pileup system.

Fig. 1 .

 1 Fig.1. Ratio of the Compton scattering/photoelectric cross-sections basis set functions (bottom graph). Since the ratio is monotonically decreasing, the ratios in Eq. 6 are different if the delta function energies are different.

Fig. 2 .

 2 Fig. 2. Closed contour used in proof of invertibility. The contour consists of segments of the two positive axes and a quadrant of a circle joining the ends of the line segments.

, 11 M

 11 rec,22 -M rec,12 M rec,21 = J0 (1+bλ1)(1+bλ2)


Adnane Saoud

Antoine Girard

Optimal Multirate Sampling in Symbolic Models for Incrementally Stable Switched Systems

Keywords: Approximate bisimulation, switched systems, symbolic control, multirate sampling, incremental stability

Methods for computing approximately bisimilar symbolic models for incrementally stable switched systems are often based on discretization of time and space, where the value of time and space sampling parameters must be carefully chosen in order to achieve a desired precision. These approaches can result in symbolic models that have a very large number of transitions, especially when the time sampling, and thus the space sampling parameters are small. In this paper, we present an approach to the computation of symbolic models for switched systems with dwell-time constraints using multirate time sampling, where the period of symbolic transitions is a multiple of the control (i.e. switching) period. We show that all the multirate symbolic models, resulting from the proposed construction, are approximately bisimilar to the original incrementally stable switched system with the precision depending on the sampling parameters, and the sampling factor between transition and control periods. The main contribution of the paper is the explicit determination of the optimal sampling factor, which minimizes the number of transitions in the class of proposed symbolic models for a prescribed precision. Interestingly, we prove that this optimal sampling factor is mainly determined by the state space dimension and the number of modes of the switched system. Finally, an illustration of the proposed approach is shown on an example, which shows the benefit of multirate symbolic models in reducing the computational cost of abstraction-based controller synthesis.

Introduction

A switched system is a dynamical system consisting of a finite number of subsystems and a law that controls the switching among them [START_REF] Liberzon | Switching in systems and control, ser. systems & control: Foundations & applications[END_REF][START_REF] Sun | Stability theory of switched dynamical systems[END_REF][START_REF] Lin | Stability and stabilizability of switched linear systems: a survey of recent results[END_REF]. The literature on switched systems principally focuses on the stability and stabilization problems. However, other objectives need also to be considered such as safety, reachability or more A preliminary version of this paper was presented at the 20th IFAC World Congress. This work was partially supported by the Labex DigiCosme (project ANR-11-LABEX-0045-DIGICOSME) operated by ANR as part of the program "Investissement d'Avenir" Idex Paris Saclay (ANR-11-IDEX-0003-02). This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 725144). Corresponding author A. Girard; e-mail: Adnane.Saoud@l2s.centralesupelec.fr, Antoine.Girard@l2s.centralesupelec.fr. complex objectives such as those expressed in linear temporal logic. For this reason, over recent years, several studies focused on the use of discrete abstractions and symbolic control techniques. The area of symbolic control is concerned with the use of algorithmic discrete synthesis techniques for designing controllers for continuous and hybrid dynamical systems (see e.g. [START_REF] Tabuada | Verification and control of hybrid systems: a symbolic approach[END_REF][START_REF] Rungger | Symbolic synthesis for cyber-physical systems[END_REF] and the references therein). A key concept in symbolic control is that of symbolic models, which consist in discrete abstractions of the continuous dynamics, and which are amenable to automata theoretic techniques for the synthesis of controllers enforcing a broad range of specifications [START_REF] Bloem | Synthesis of reactive (1) designs[END_REF][START_REF] Belta | Formal Methods for Discrete-Time Dynamical Systems[END_REF]. Controllers for the original system, with strong formal guarantees, can then be obtained through dedicated refinement procedures [START_REF] Tabuada | Verification and control of hybrid systems: a symbolic approach[END_REF][START_REF] Girard | Controller synthesis for safety and reachability via approximate bisimulation[END_REF][START_REF] Reissig | Feedback refinement relations for the synthesis of symbolic controllers[END_REF]. This latter step requires the original system and the symbolic model to be related by some formal behavioral relationship such as simulation, bisimulation relations or their approximate and alternated versions [START_REF] Tabuada | Verification and control of hybrid systems: a symbolic approach[END_REF][START_REF] Girard | Approximation metrics for discrete and continuous systems[END_REF].

Numerous works have been dedicated to the computation of symbolic models for various classes of dynamical systems. Focusing on approximately bisimilar abstractions, existing approaches make it possible to deal with nonlinear systems [START_REF] Pola | Approximately bisimilar symbolic models for nonlinear control systems[END_REF][START_REF] Pola | Symbolic models for nonlinear control systems: Alternating approximate bisimulations[END_REF], switched systems [START_REF] Girard | Approximately bisimilar symbolic models for incrementally stable switched systems[END_REF], timedelay systems [START_REF] Pola | Symbolic models for nonlinear time-delay systems using approximate bisimulations[END_REF][START_REF] Pola | Symbolic models for time-varying time-delay systems via alternating approximate bisimulation[END_REF], networked control systems [START_REF] Borri | Symbolic control design of nonlinear networked control systems[END_REF][START_REF] Zamani | Finite abstractions of networked control systems[END_REF], stochastic systems [START_REF] Zamani | Approximately bisimilar symbolic models for randomly switched stochastic systems[END_REF][START_REF] Zamani | Symbolic control of stochastic systems via approximately bisimilar finite abstractions[END_REF]... All these approaches are essentially based on discretization of time and space and require the considered system to satisfy some kind of incremental stability property [START_REF] Angeli | A Lyapunov approach to incremental stability properties[END_REF]. However, incremental stability can be dropped if one seeks symbolic models related only by one-sided approximate simulation relations [START_REF] Tabuada | An approximate simulation approach to symbolic control[END_REF][START_REF] Zamani | Symbolic models for nonlinear control systems without stability assumptions[END_REF]. In most cases, symbolic models of arbitrary precision can be obtained by carefully choosing time and space sampling parameters. However, for a given precision, the choice of a small time sampling parameter imposes to choose a small space sampling parameter resulting in symbolic models with a prohibitively large number of transitions. This constitutes a limiting factor of the approach because the size of the symbolic models is crucial for computational efficiency of discrete controller synthesis algorithms. Several studies have been conducted in order to address this issue by enabling the computation of more parsimonious symbolic models with smaller numbers of transitions. Such approaches include compositional abstraction schemes where symbolic models of a system are built from symbolic models of its components [START_REF] Tazaki | Bisimilar finite abstractions of interconnected systems[END_REF][START_REF] Pola | Symbolic models for networks of control systems[END_REF][START_REF] Majumdar | Compositional synthesis of finite state abstractions[END_REF]; multi-resolution or multi-scale symbolic models computed using non-uniform adaptive space discretizations [START_REF] Tazaki | Discrete-state abstractions of nonlinear systems using multi-resolution quantizer[END_REF][START_REF] Girard | Safety controller synthesis for incrementally stable switched systems using multiscale symbolic models[END_REF]; symbolic models where the set of symbolic states is not given by a discretization of the state-space but by input sequences [START_REF] Le Corronc | Mode sequences as symbolic states in abstractions of incrementally stable switched systems[END_REF][START_REF] Zamani | Symbolic models for stochastic switched systems: A discretization and a discretization-free approach[END_REF].

In this paper, we show how the size of symbolic models can be reduced using multirate sampling. Multirate sampling has been introduced in the area of sampled-data systems to face some of the sampling processes disadvantages such as the loss of relative degree and changes in the properties of the zero dynamics (see e.g. [START_REF] Monaco | An introduction to motion planning under multirate digital control[END_REF][START_REF] Grizzle | Feedback linearization of sampled-data systems[END_REF][START_REF] Monaco | Issues on nonlinear digital control[END_REF]). In this paper, we present an approach to the computation of multirate symbolic models for incrementally stable switched systems, where the period of symbolic transitions is a multiple of the control (i.e. switching) period. A similar approach has been explored in the symbolic control literature in the context of nonlinear digital control systems [START_REF] Majumdar | Approximately bisimilar symbolic models for digital control systems[END_REF]. The first contribution of the paper is to extend this approach to the class of switched systems, with dwell-time constraints. We show that the obtained multirate symbolic models are approximately bisimilar to the original switched system. Then, the second and main contribution of the paper lies in the explicit determination of the optimal sampling factor between transition and control periods, which minimizes the number of transitions in the class of proposed symbolic models for a prescribed precision; this problem is not considered in [START_REF] Majumdar | Approximately bisimilar symbolic models for digital control systems[END_REF]. Interestingly, we show that the optimal sampling factor is mainly determined by the state space dimension and the number of modes of the switched system. This paper is organized as follows. In Section 2, we introduce the class of incrementally stable switched systems under study and we present the abstraction framework used in the paper. In Section 3, we present the construction of symbolic models for incrementally stable switched systems with dwell-time constraints, using multirate sampling. In Section 4, we establish the optimal sampling factor between control and transition periods which minimizes the number of transitions in the symbolic model. Finally, in Section 5, we illustrate our approach using an example taken from [START_REF] Girard | Approximately bisimilar symbolic models for incrementally stable switched systems[END_REF], which shows the benefits of the proposed multirate symbolic models.

A preliminary version of this work has been presented in the conference paper [START_REF] Saoud | Multirate symbolic models for incrementally stable switched systems[END_REF] where switched systems without dwell-time constraints are considered. The current paper extends the approach to consider dwell-time constraints; results of [START_REF] Saoud | Multirate symbolic models for incrementally stable switched systems[END_REF] being recovered as particular cases. We also provide novel numerical experiments.

Notations: Z, N and N + denote the sets of integers, of non-negative integers and of positive integers, respectively. R, R + 0 and R + denote the sets of real numbers, of non-negative real numbers, and of positive real numbers, respectively. For s ∈ R + 0 , s denote its integer part, i.e. the largest nonnegative integer r ∈ N such that r ≤ s. For x ∈ R n , x denotes the Euclidean norm (i.e. the 2-norm) of x. A continuous function γ : R + 0 → R + 0 is said to belong to class K if it is strictly increasing and γ(0) = 0; γ is said to belong to class K ∞ if γ is of class K and γ(s) → ∞ as s → ∞. A continuous function β : R + 0 × R + 0 → R + 0 is said to belong to class KL if, for all fixed t ∈ R + 0 , the map β(•, t) belongs to class K, and for all fixed s ∈ R + , the map β(s, •) is strictly decreasing and β(s, t) → 0 as t → ∞.

Preliminaries

Incrementally stable switched systems

We introduce the class of switched systems: Definition 1 A switched system is a quadruple Σ = (R n , P, P, F ), consisting of:

• a state space R n ; • a finite set of modes P = {1, . . . , m}; • a set of switching signals P ⊆ S(R + 0 , P ), where S(R + 0 , P ) denotes the set of piecewise constant functions from R + 0 to P , continuous from the right and with a finite number of discontinuities on every bounded interval of R + 0 ; • a collection of vector fields F = {f 1 , . . . , f m }, indexed by P .

The discontinuities 0 < t 1 < t 2 < . . . of a switching signal are called switching times; by definition of S(R + 0 , P ), there are only a finite number of switching times on every bounded interval of R + 0 and thus Zeno behaviors are avoided. A switching signal p ∈ S(R + 0 , P ) has dwelltime τ d ∈ R + if the sequence of switching times satisfies t k+1 -t k ≥ τ d , for all k ≥ 1. The set of switching signals with dwell-time τ d is denoted S τ d (R + 0 , P ).

A piecewise C 1 function x : R + 0 → R n is said to be a trajectory of Σ if it is continuous and there exists a switching signal p ∈ P such that, at each t ∈ R + 0 where the function p is continuous, x is continuously differentiable and satisfies: ẋ(t) = f p(t) (x(t)).

(1) We make the assumption that the vector fields f p , p ∈ P , are locally Lipschtiz and forward complete (see e.g. [START_REF] Angeli | Forward completeness, unboundedness observability, and their Lyapunov characterizations[END_REF] for necessary and sufficient conditions), so that for all switching signals p ∈ P and all initial states x ∈ R n , there exists a unique trajectory, solution to (1) with x(0) = x, denoted x(., x, p). We will denote by φ p t the flow associated to the vector field f p . Then, for a constant switching signal given by p(t) = p, for all t ∈ R + 0 , we have x(t, x, p) = φ p t (x), for all t ∈ R + 0 .

In the following, we consider incrementally globally uniformly asymptotically stable (δ-GUAS) switched systems, see [START_REF] Girard | Approximately bisimilar symbolic models for incrementally stable switched systems[END_REF] for a formal definition. Intuitively, incremental stability means that all trajectories associated to the same switching signal converge to the same trajectory, independently of their initial conditions. Sufficient conditions for incremental stability are given in [START_REF] Girard | Approximately bisimilar symbolic models for incrementally stable switched systems[END_REF] in terms of existence of multiple Lyapunov functions.

Definition 2 : Smooth functions V p : R n × R n → R + 0 , p ∈ P , are multiple δ-GUAS Lyapunov functions for Σ if there exist K ∞ functions α, α, κ ∈ R + and µ ≥ 1 such that for all x, y ∈ R n , and p, p ∈ P ,

α( x -y ) ≤ V p (x, y) ≤ α( x -y );
(2)

∂V p ∂x (x, y)f p (x) + ∂V p ∂y (x, y)f p (y) ≤ -κV p (x, y); (3) V p (x, y) ≤ µV p (x, y). (4)
In [START_REF] Girard | Approximately bisimilar symbolic models for incrementally stable switched systems[END_REF], it is proved that Σ = (R n , P, P, F ) is δ-GUAS if there exist multiple δ-GUAS Lyapunov functions for Σ and the set of switching signals

P ⊆ S τ d (R + 0 , P ) with dwell-time τ d > ln(µ) κ .
In this paper, we assume that the previous condition holds and in order to construct symbolic models for the switched systems, we shall make the supplementary assumption that there exists a K ∞ function γ such that for multiple δ-GUAS Lyapunov functions we have:

∀x, y, z ∈ R n , p ∈ P, |V p (x, y) -V p (x, z)| ≤ γ( y -z ). ( 5
)
Remark 1 In [START_REF] Girard | Approximately bisimilar symbolic models for incrementally stable switched systems[END_REF], it is shown that if we are interested in the dynamics of the switched system on a compact set C ⊆ R n and V p , p ∈ P , are C 1 on C, then, (5) holds with the linear K ∞ function given by γ(s) = c γ s where

c γ = max x,y∈C,p∈P ∂V p ∂y (x, y) ,
Remark 2 For all x ∈ R n , (2) implies that for all p ∈ P , V p (x, x) = 0, then for all x, y ∈ R n , p ∈ P we have from (5) that:

V p (x, y) ≤ |V p (x, y) -V p (x, x)| ≤ γ( x -y ).
Then, there is no loss of generality in assuming that the second inequality in (2) holds with α = γ.

Remark 3 The smoothness assumption on functions V p , p ∈ P , in Definition 2 can be relaxed, and condition (3) can be replaced by

V p (φ p t (x), φ p t (y)) ≤ e -κt V p (x, y) (6)
for all x, y ∈ R n , p ∈ P and t ∈ R + 0 .

Transition systems

We present the notion of transition systems, which allows us to describe, in a unified framework, switched systems and their symbolic models.

Definition 3 A transition system is a tuple T = (X, U, Y, ∆, X 0 ) consisting of:

• a set of states X;

• a set of inputs U ;

• a set of outputs Y ;

• a transition relation ∆ ⊆ X × U × X × Y ; • a set of initial states X 0 ⊆ X.
T is said to be metric if the set of outputs Y is equipped with a metric d, symbolic if X and U are finite or countable sets.

The transition (x, u, x , y) ∈ ∆ will be denoted (x , y) ∈ ∆(x, u) and means that the system can evolve from state x to state x under the action of input u, while producing output y. An input u ∈ U belongs to the set of enabled inputs at state x ∈ X, denoted enab

∆ (x), if ∆(x, u) = ∅.
T is said to be deterministic if for all x ∈ X and for all u ∈ enab ∆ (x), ∆(x, u) consists of a unique element. State x ∈ X is said to be blocking if enab ∆ (x) = ∅, otherwise it is said to be non-blocking.

A trajectory of the transition system is a finite or infinite sequence of transitions σ = (x 0 , u 0 , y

0 )(x 1 , u 1 , y 1 ) (x 2 , u 2 , y 2 ) . . . where (x i+1 , y i ) ∈ ∆(x i , u i ), for i ≥ 0. It is initialized if x 0 ∈ X 0 . A state x ∈ X is reach- able if
there exists an initialized trajectory such that x i = x, for some i ≥ 0. The transition system is said to be non-blocking if all reachable states are non-blocking.

The output behavior associated to the trajectory σ is the sequence of outputs y 0 y 1 y 2 . . .

In this paper, we consider the approximation relationship for transition systems based on the notion of approximate bisimulation [START_REF] Girard | Approximation metrics for discrete and continuous systems[END_REF], which requires that the distance between the output behaviors of two transition systems remains bounded by some specified precision.

The following definition is taken from [START_REF] Girard | Safety controller synthesis for incrementally stable switched systems using multiscale symbolic models[END_REF] and generalizes that of [START_REF] Girard | Approximation metrics for discrete and continuous systems[END_REF] to accommodate the encoding of the output map within the transition relation.

Definition 4 Let T 1 = (X 1 , U, Y, ∆ 1 , X 0 1 ) and T 2 = (X 2 , U, Y, ∆ 2 , X 0
2 ) be two metric transition systems with the same input set U and the same output set Y equipped with a metric d. Let ε ≥ 0 be a given precision. A relation R ⊆ X 1 × X 2 is said to be an ε-approximate bisimulation relation between T 1 and T 2 if for all (x 1 , x 2 ) ∈ R, u ∈ U :

∀(x 1 , y 1 ) ∈ ∆ 1 (x 1 , u), ∃(x 2 , y 2 ) ∈ ∆ 2 (x 2 , u), d(y 1 , y 2 ) ≤ ε and (x 1 , x 2 ) ∈ R; ∀(x 2 , y 2 ) ∈ ∆ 2 (x 2 , u), ∃(x 1 , y 1 ) ∈ ∆ 1 (x 1 , u), d(y 1 , y 2 ) ≤ ε and (x 1 , x 2 ) ∈ R.
The transition systems T 1 and T 2 are said to be εapproximately bisimilar, denoted T 1 ∼ ε T 2 , if:

• ∀x 1 ∈ X 0 1 , ∃x 2 ∈ X 0 2 , such that (x 1 , x 2 ) ∈ R; • ∀x 2 ∈ X 0 2 , ∃x 1 ∈ X 0 1 , such that (x 1 , x 2 ) ∈ R.
The approximate bisimulation relation guarantees that for each output behavior of T 1 (respectively of T 2 ), there exists an output behavior of T 2 (respectively of T 1 ) such that the distance between these output behaviors is uniformly bounded by ε (see [START_REF] Girard | Approximation metrics for discrete and continuous systems[END_REF]).

Symbolic models with multirate sampling

In this section, we extend the results of [START_REF] Saoud | Multirate symbolic models for incrementally stable switched systems[END_REF] to the case of switched systems with dwell-time constraints. Let us consider a switched system Σ τ d = (R n , P, P τ d , F ), in which the switching is periodically controlled with control period τ ∈ R + and in which a dwell-time τ d ∈ R + is imposed on switching signals. For simplicity, we assume that τ = τ d /k where k ∈ N + .

The sampled dynamics of Σ τ d can then be described by the transition system T τ (Σ τ d ) = (X, U, Y, ∆ τ , X 0 ) as follows:

• the set of states is X = R n × P ;

• the set of inputs is U = P ;

• the set of outputs is Y = R n ∪ R k×n ;

• the transition relation is given for (x, p), (x , p ) ∈ X, u ∈ U , y ∈ Y , by ((x , p ), y) ∈ ∆ τ ((x, p), u) if and only if

x = φ u τ (x), p = u y = x if u = p x = φ u kτ (x), p = u, y = (x, φ u τ (x), . . . , φ u (k-1)τ (x)) if u = p • the set of initial states X 0 = R n × P .
We should emphasize that transitions in T τ (Σ τ d ) have either duration τ or τ d = kτ . The state (x, p) ∈ X indicates that the state of the switched system is x ∈ R n and that the active mode is p ∈ P . Then, one can either go on with mode p, which corresponds to the first type of transitions of duration τ ; or switch to another mode p = p, which corresponds to the second type of transitions where the new mode p is held for duration τ d . It is easy to see that the dwell-time constraint is fulfilled by construction. It is noteworthy that this construction differs from, and is more compact than, that of [START_REF] Girard | Approximately bisimilar symbolic models for incrementally stable switched systems[END_REF].

T τ (Σ τ d ) is non-blocking (enab ∆τ ((x, p)) = U for all (x, p) ∈ X),

Multirate sampling of switched systems

In the previous transition system, one transition coincides with the control period (of duration τ or kτ if the transition system keeps the same mode or changes it, respectively). In this paper, we consider multirate sampling where a transition corresponds to a sequence of r control periods, where the sampling factor r ∈ N + . Thus, multirate transition systems are obtained by concatenating r successive transitions of T τ (Σ τ d ).

Let us define T r τ (Σ τ d ) = (X, U r , Y r , ∆ r τ , X 0 ) where:

• the set of states is X = R n × P ;

• the set of inputs is U r = P r ;

• the set of outputs Y r = (R n ∪ R k×n ) r ;

• the transition relation is given for (x, p), (x , p ) ∈ X, u ∈ U r , with u = (u 1 , . . . , u r ), and y ∈ Y r , with y = (y 1 , . . . , y r ) by ((x , p ), y) ∈ ∆ r τ ((x, p), u) if and only if (x, p) = (x 1 , p 1 ), (x , p ) = (x r+1 , p r+1 ), with ((x i+1 , p i+1 ), y i ) ∈ ∆ τ ((x i , p i ), u i ), i = 1, . . . , r

• the set of initial states X 0 = R n × P .

T r τ (Σ τ d ) is non-blocking (enab ∆ r τ ((x, p)) = U r for all (x, p) ∈ X), deterministic, and metric when the set of outputs Y r is equipped with the following metric d Y r :

∀y = (y 1 , . . . , y r ), y = (y 1 , . . . , y r ) ∈ Y r , d Y r (y, y ) = r max i=1 d Y (y i , y i ). ( 7
)
Let us remark that for r = 1, T r τ (Σ τ d ) coincides with T τ (Σ τ d ). When r = 1, the following result shows that T τ (Σ τ d ) and T r τ (Σ τ d ) produce equivalent infinite output behaviors. This result is straightforward and is stated without proof.

Claim 1 For any infinite output behavior (y 0 , y 1 , y 2 , . . . ) of T τ (Σ τ d ), there exists an infinite output behavior

(z 0 , z 1 , z 2 , . . . ) of T r τ (Σ τ d ) with z i = (z i 1 , . . . , z i r ) such that ∀i ∈ N, j = 1, . . . , r, z i j = y ir+j-1 . (8)
Conversely, for any infinite output behavior (z 0 , z 1 , z 2 , . . . ) of T r τ (Σ τ d ) with z i = (z i 1 , . . . , z i r ), there exists an infinite output behavior (y 0 , y 1 , y 2 , . . . ) of T τ (Σ τ d ) such that (8) holds.

Remark 4 Using T τ (Σ τ d ) or T r τ (Σ τ d ) for the purpose of synthesis provides identical guarantees on the sampled behavior of the switched system, since the infinite output behaviors of both transition systems are equivalent. However, it leads to different implementations of switching controllers. For controllers synthesized using T τ (Σ τ d ), the sensing and actuation periods (of duration τ or kτ ) are equal; while for controllers synthesized using T r τ (Σ τ d ), the sensing period consists of r actuation periods. In the latter case, at sensing instants, the controller selects a sequence of r modes, each of which is actuated for a duration τ or kτ .

Construction of symbolic models

For an incrementally stable switched system Σ with multiple δ-GUAS Lyapunov functions, a construction of symbolic models that are approximately bisimilar to T τ (Σ τ d ) has been proposed in [START_REF] Girard | Approximately bisimilar symbolic models for incrementally stable switched systems[END_REF], based on a discretization of the state-space R n . Theorem 4.2 in that paper, shows that symbolic models of arbitrary precision can be computed by using a sufficiently fine discretization of the state-space. However, this usually results in symbolic models that have a very large number of transitions, especially when the control period τ is small.

In this section, we establish a similar result for the multirate transition system T r τ (Σ τ d ). This idea is inspired by the work presented in [START_REF] Majumdar | Approximately bisimilar symbolic models for digital control systems[END_REF], in which symbolic models are computed for digital control systems using multirate sampling. Our results can be seen as an extension to the class of switched systems with dwell-time. In addition, in the following sections, we will provide a theoretical analysis allowing us to choose the optimal sampling factor r, minimizing the number of transitions in the symbolic model, which is not available in [START_REF] Majumdar | Approximately bisimilar symbolic models for digital control systems[END_REF].

Let η ∈ R + be a space sampling parameter, the set of states R n is approximated by the lattice:

[R n ] η = q ∈ R n |q i = k i 2η √ n , k i ∈ Z, i = 1, . . . , n .
We associate a quantizer

Q η : R n -→ [R n ] η given by Q η (x) = q if and only if ∀i = 1, . . . , n, q i - η √ n ≤ x i < q i + η √ n .
where x i and q i denote the i-th coordinates, i = 1, . . . , n of x and q, respectively. We can easily show that for all

x ∈ R n , Q η (x) -x ≤ η.
Let us define the transition system T r τ,η (Σ τ d ) = (X η , U r , Y r , ∆ r τ,η , X 0 η ) as follows:

• the set of states is

X η = [R n ] η × P ; • the set of inputs is U r = P r ; • the set of outputs Y r = (R n ∪ R k×n ) r ;
• the transition relation is given for (q, p), (q , p ) ∈ X η , u ∈ U r , y ∈ Y r , by ((q , p ), y) ∈ ∆ r τ,η ((q, p), u) if and only if q = Q η (x ) and ((x , p ), y) ∈ ∆ r τ ((q, p), u);

• the set of initial states is

X 0 η = [R n ] η × P .
T r τ,η (Σ τ d ) is symbolic, non-blocking (enab ∆ r τ,η ((q, p)) = U r for all (q, p) ∈ X η ), deterministic and metric when the set of outputs Y r is equipped with the metric d Y r given by [START_REF] Girard | Controller synthesis for safety and reachability via approximate bisimulation[END_REF].

Theorem 1 Consider a switched system Σ τ d , and let us assume that there exist multiple δ-GUAS Lyapunov functions V p , p ∈ P , for Σ τ d such that (5) holds for some K ∞ function γ, let the dwell-time τ d > ln(µ) κ . Let time and space sampling parameters τ, η ∈ R + , sampling factor r ∈ N + and precision ε ∈ R + satisfy:

η ≤ γ -1 1 µ (1 -λ(τ ) r )α(ε) (9)
where λ(τ ) = max(e -κτ , µe -κτ d ), then, the transition systems T r τ (Σ τ d ) and T r τ,η (Σ τ d ) are ε-approximately bisimilar.

PROOF. Let us prove that the relation R defined by: R = ((x, p 1 ), (q, p

2 )) ∈ X × X η p 1 = p 2 = p V p (x, q) ≤ 1 µ α(ε)
is an ε-approximate bisimulation relation between T r τ (Σ τ d ) and T r τ,η (Σ τ d ).

Let ((x, p 1 ), (q, p 2 )) ∈ R, then we have p 1 = p 2 = p and V p (x, q) ≤ 1 µ α(ε). Let u = (u 1 , . . . , u r ) ∈ U r and ((x , p ), y) ∈ ∆ r τ ((x, p), u), where y = (y 1 , . . . , y r ) ∈ Y r , then by definition of ∆ r τ :

(x, p) = (x 1 , p 1 ), (x , p ) = (x r+1 , p r+1 ), with ((x i+1 , p i+1 ), y i ) ∈ ∆ τ ((x i , p i ), u i ), i = 1, . . . , r.
Similarly, let ((q , p ), z) ∈ ∆ r τ,η ((q, p), u), where z = (z 1 , . . . , z r ) ∈ Y r , then by definition of ∆ r τ,η :

(q, p) = (q 1 , p 1 ), (q , p ) = (Q η (q r+1 ), p r+1 ), with ((q i+1 , p i+1 ),

y i ) ∈ ∆ τ ((q i , p i ), u i ), i = 1, . . . , r.
By the definition of ∆ τ , we have for all i = 1, . . . , r, p i+1 = u i , and

V pi+1 (x i+1 , q i+1 ) ≤ e -κτ V pi (x i , q i ), if p i = p i+1 , V pi+1 (x i+1 , q i+1 ) ≤ µ e -κτ d V pi (x i , q i ), if p i = p i+1 .
where the two inequalities are obtained by ( 3) and ( 4). Then, it follows that for all i = 1, . . . , r + 1,

V pi (x i , q i ) ≤ λ(τ ) i-1 V p1 (x 1 , q 1 ) ≤ λ(τ ) i-1 1 µ α(ε). ( 10
)
Then, from ( 5), ( 10) and ( 9), we have

V p (x , q ) = V pr+1 (x r+1 , Q η (q r+1 )) ≤ V pr+1 (x r+1 , q r+1 ) + γ(η) ≤ λ(τ ) r 1 µ α(ε) + γ(η) ≤ 1 µ α(ε).
Thus, ((x , p ), (q , p )) ∈ R. 2), [START_REF] Girard | Approximately bisimilar symbolic models for incrementally stable switched systems[END_REF] and since λ(τ ) ≤ 1 and

Let i = 1, . . . , r, if u i = p i , we have y i = x i , z i = q i , then from (
1 µ ≤ 1, d Y (y i , z i ) = x i -q i ≤ α -1 (V pi (x i , q i )) ≤ ε. ( 11
)
If u i = p i , we have y i = (y i,1 , . . . , y i,k ), z i = (z i,1 , . . . , z i,k ) where y i,j = φ ui (j-1)τ (x i ) and z i,j = φ ui (j-1)τ (q i ), j = 1, . . . , k. Then, from (3), ( 4), [START_REF] Girard | Approximately bisimilar symbolic models for incrementally stable switched systems[END_REF] and since λ(τ ) ≤ 1, we have for all j = 1, . . . , k, V ui (y i,j , z i,j ) ≤ V ui (x i , q i ) ≤ µV pi (x i , q i ) ≤ α(ε).

Then, by (2), we have for all j = 1, . . . , k,

y i,j -z i,j ≤ α -1 (V ui (y i,j , z i,j )) ≤ ε. Hence, d Y (y i , z i ) = k max j=1 y i,j -z i,j ≤ ε. ( 12
)
It then follows from ( 11), ( 12) that

d Y r (y, z) = r max i=1 d Y (y i , z i ) ≤ ε.
Hence, the first condition of Definition 4 holds.

In a similar way, we prove that for all ((q , p ), z) ∈ ∆ r τ,η ((q, p), u) there exists ((x , p ), y) ∈ ∆ r τ ((x, p), u) such that ((x , p ), (q , p )) ∈ R and d Y r (y, z) ≤ ε. Hence, R is an ε-approximate bisimulation relation between T r τ (Σ τ d ) and T r τ,η (Σ τ d ).

Now, let (x, p) ∈ X 0 = R n × P , and (q, p) ∈ X 0 η = [R n ] η × P , given by q = Q η (x), then x -q ≤ η. Following Remark 2, we have that the second inequality of (2) holds with α = γ. It follows that

V p (x, q) ≤ γ( x -q ) ≤ γ(η) ≤ 1 µ α(ε)
where the last inequality comes from [START_REF] Girard | Approximation metrics for discrete and continuous systems[END_REF]. Hence ((x, p), (q, p)) ∈ R. Conversely, for all (q, p) ∈ X 0 η = [R n ] η ×P , let (x, p) ∈ X 0 = R n ×P , given by x = q, then V p (x, q) = 0 and ((x, p), (q, p)) ∈ R. Hence, T r τ (Σ τ d ) and T r τ,η (Σ τ d ) are ε-approximately bisimilar. 2 Some remarks regarding the size of the symbolic models are in order. It appears from ( 9) that, for a given precision ε ∈ R + and control period τ ∈ R + , using larger sampling factor r ∈ N + allows us to use larger values of η ∈ R + and thus coarser discretizations of the state space. This results in symbolic models with fewer symbolic states. However, the number of transitions initiating from a symbolic state is m r and thus grows exponentially with the sampling factor. Hence, the advantage of using multirate symbolic models in terms of number of transitions in the symbolic model is still unclear. This issue is addressed in the following section, where we determine the optimal value of the sampling factor.

Remark 5

The results of [START_REF] Saoud | Multirate symbolic models for incrementally stable switched systems[END_REF] on multirate symbolic models for incrementally stable switched systems (without dwell-time) can be recovered easily from the particular case of the present study when µ = 1 (i.e. there exists common δ-GUAS Lyapunov function), and k = 1 (i.e. the dwell-time coincides with the control period, τ d = τ ).

Optimal sampling factor

In the following, we consider multirate symbolic models T r τ,η (Σ τ d ) computed using the approach described above, where we restrict the set of states to some compact set C ⊆ R n with nonempty interior. The number of symbolic states in X η ∩ (C × P ) can be accurately estimated by v C η n × m, where v C ∈ R + is a positive constant proportional to the volume of C. Then the number of symbolic transitions initiating from states in

X η ∩ (C × P ) is v C m r+1
η n . We assume that the number of modes m ≥ 2.

Problem formulation

In this section, given a desired precision ε ∈ R + , and a control period τ ∈ R + , we establish the optimal values r * ∈ N + and η * ∈ R + , which characterizes the multirate symbolic model T r τ,η (Σ τ d ) of precision ε (as guaranteed by Theorem 1) with the minimal number of symbolic transitions initiating from states in X η ∩ (C × P ).

Since, C is a compact set, following Remark 1, we assume that (5) holds for a linear K ∞ function γ given by γ(s) = c γ s where c γ ∈ R + . Thus, we aim at solving the following mixed integer nonlinear program:

Minimize v C m r+1 η n over r ∈ N + , η ∈ R + under η ≤ (1 -λ(τ ) r ) α(ε) µcγ (13)
where the inequality constraints comes from (9) to guarantee a precision ε.

Let us first remark that for a given r ∈ N + , the optimal value η ∈ R + is obviously obtained as

η = (1 - λ(τ ) r ) α(ε)
µcγ . It follows that ( 13) is equivalent to the following integer program:

Minimize v C (µcγ ) n (α(ε)) n m r+1 (1-λ(τ ) r ) n over r ∈ N + (14)
The value v C (µcγ ) n (α(ε)) n ∈ R + does not depend on r and thus does not affect the solution of ( 14), which can finally be equivalently formulated as:

Minimize g(r) = m r+1 (1-λ(τ ) r ) n over r ∈ N + (15)
A first interesting information that can be inferred from [START_REF] Majumdar | Compositional synthesis of finite state abstractions[END_REF] is that the optimal sampling factor only depends on the control period τ ∈ R + , the dimension of the statespace n ∈ N + , the number of modes m ∈ N + , the dwell time τ d , the constant µ ≥ 1 given in ( 4) and the decay rate κ ∈ R + of the multiple δ-GUAS Lyapunov functions. In particular, it is noteworthy that it is independent of the desired precision ε ∈ R + and of the compact set C.

Explicit solution

In this section, we show that the previous optimization problems can be solved explicitly. We first consider the relaxation of the integer program [START_REF] Majumdar | Compositional synthesis of finite state abstractions[END_REF] over the positive real numbers:

Lemma 1 Let g : R + → R + be given as in [START_REF] Majumdar | Compositional synthesis of finite state abstractions[END_REF]. Then, g has a unique minimizer r * ∈ R + given by

r * = 1 -ln(λ(τ )) ln 1 - n ln(λ(τ )) ln(m) (16)
with λ(τ ) = max(e -κτ , µe -κτ d ).

Moreover, g is strictly decreasing on (0, r * ] and strictly increasing on [r * , +∞).

PROOF. Let us compute the first order derivative of g:

g (r) = m (1 -λ(τ ) r ) 2n ln(m)m r (1 -λ(τ ) r ) n + m r n ln(λ(τ ))λ(τ ) r (1 -λ(τ ) r ) n-1 = m r+1 (1 -λ(τ ) r ) n+1 ln(m)(1 -λ(τ ) r ) + n ln(λ(τ ))λ(τ ) r = ln(m)m r+1 (1 -λ(τ ) r ) n+1 1 -λ(τ ) r 1 - n ln(λ(τ ))
ln(m) .

By remarking that ln(m)m r+1 (1-λ(τ ) r ) n+1 > 0 for all r ∈ R + , it is easy to see that g (r) has the same sign as

1 -λ(τ ) r 1 -n ln(λ(τ )) ln(m)
, and that it is negative on (0, r * ), zero at r * and positive on (r * , +∞). The result stated in Lemma 1 follows immediately. 2

We can now state the main result of the section: Theorem 2 For any desired precision ε ∈ R + , and any control period τ ∈ R + , the optimal parameters r * ∈ N + and η * ∈ R + , solutions of [START_REF] Liberzon | Switching in systems and control, ser. systems & control: Foundations & applications[END_REF], which minimize the number of symbolic transitions of T r τ,η (Σ τ d ), initiating from states in X η ∩ (C × P ), while satisfying [START_REF] Girard | Approximation metrics for discrete and continuous systems[END_REF], are given by r * = r * or r * = r * + 1 [START_REF] Monaco | An introduction to motion planning under multirate digital control[END_REF] and

η * = (1 -λ(τ ) r * ) α(ε) µcγ (18)
where r * is given by [START_REF] Majumdar | Approximately bisimilar symbolic models for digital control systems[END_REF].

PROOF.

From Lemma 1, it follows that ∀r ∈ N + , with r < r * , g(r) > g( r * ) and ∀r ∈ N + , with r > r * + 1, g(r) > g( r * + 1).

Then, it follows that the minimal value of g over N + is obtained for r * = r * or r * = r * + 1. Then, from the discussions in Section 4.1, it follows that the solution of ( 13) is given by r * and η * = (1 -λ(τ ) r * )

α(ε) µcγ . 2
In practice, we compute the optimal parameters of the multirate symbolic models by evaluating the function g at r * and r * + 1. We then pick the one, out of two possible values of r * , which minimizes g and compute η * using [START_REF] Monaco | Issues on nonlinear digital control[END_REF].

We would like to point out that the previous result can be applied to either linear or nonlinear switched systems.

The only requirement is that we restrict the analysis to a compact subset of R n . Finally, it is interesting to remark that for small values of the control period τ ∈ R + , the optimal sampling factor r * is mainly determined by the state space dimension and the number of modes.

Corollary 1 There exists τ ∈ R + , such that for any desired precision ε ∈ R + , and any control period τ ∈ (0, τ ], the optimal parameters r * ∈ N + and η * ∈ R + , solutions of [START_REF] Liberzon | Switching in systems and control, ser. systems & control: Foundations & applications[END_REF], which minimize the number of symbolic transitions of T r τ,η (Σ τ d ), initiating from states in X η ∩ (C × P ), while satisfying [START_REF] Girard | Approximation metrics for discrete and continuous systems[END_REF], are given by

r * = n ln(m) or r * = n ln(m) + 1 and η * = (1 -e -r * κτ ) α(ε) µc γ .
PROOF. First let us remark that for τ ≤ τ d -ln(µ) κ , we have λ(τ ) = e -κτ . Then, let τ be given by Now, let us remark that for all θ ∈ R + , we have that θ(1 -θ 2 ) ≤ ln(1 + θ) ≤ θ. Let r * be given by ( 16), then it follows from the previous inequalities that for all τ ∈ R + .

τ = min   τ d - ln(µ) κ , 2 ln(m) nκ   1 - n ln(m) n ln(m)     .
n ln(m) 1 - nκτ 2 ln(m) ≤ r * ≤ n ln(m)
.

Then, using [START_REF] Pola | Approximately bisimilar symbolic models for nonlinear control systems[END_REF], it follows that for all τ ∈ (0, τ ],

n ln(m) ≤ r * ≤ n ln(m)
which implies that r * = n ln(m) . The stated result is then a consequence of Theorem 2. 2

Illustrating example

In this section, we illustrate our main results and demonstrate the benefits of the proposed approach by considering the same example as in [START_REF] Girard | Approximately bisimilar symbolic models for incrementally stable switched systems[END_REF]. We consider a twodimensional switched affine system with two modes (i.e. n = 2, m = 2) and given by ẋ

(t) = A p(t) x(t) + b p(t) with b 1 = [-0.25 -2] T , b 2 = [0.25 1] T and A 1 = -0.25 1 -2 -0.25 , A 2 = -0.25 2 -1 -0.25 .
The system does not have a common δ-GUAS Lyapunov function but admits multiple δ-GUAS Lyapunov functions of the form V p (x, y) = (x -y) T M p (x -y), with

M 1 = [ 2 0 0 1 ] , M 2 = [ 1 0 0 2 ] .
Let us first remark that even if the Lyapunov function V i , i ∈ {1, 2} is not smooth, it satisfies condition (6) with κ = 0.25. The equations ( 2), ( 4) and ( 5) hold with α

(s) = s, α(s) = √ 2s, µ = √ 2 and γ(s) = √ 2s. Imposing a dwell-time τ d = 2 > ln(µ)
κ , the switched system is incrementally stable.

We compute multirate symbolic models using the approach described in Section 5. We set the control period τ = 0.2 (i.e. k = 10) and the desired precision ε = 0.25 × √ 2. We restrict the dynamics to a compact subset of R 2 given by C = [-6, 6] × [-4, 4]. We compute the symbolic models for several sampling factors r = 1, . . . , 9, the space sampling parameter is then chosen as η = (1 -λ(τ ) r ) α(ε) µcγ . Figure 1 shows the number of symbolic transitions as a function of r and we can see that this number is minimal for r = 3.

Using [START_REF] Pola | Approximately bisimilar symbolic models for nonlinear control systems[END_REF], we compute τ = 0.61. Thus, τ ∈ (0, τ ] and the assumptions of Corollary 1 hold. In particular, since n ln(m) = 2.89, the optimal sampling factor is either 2 or 3. We can then check numerically that the optimal sampling factor is r * = 3, which is consistent with the experimental data.

We now synthesize safety controllers (see e.g. [START_REF] Tabuada | Verification and control of hybrid systems: a symbolic approach[END_REF]), which keep the output of the symbolic models inside the compact region C while avoiding C = [-1.5, 1.5] × [-1, 1]. Figure 1 reports the computation times for generating symbolic models and synthesizing controllers for r = 1, . . . , 9. The time for generating the symbolic model is linear with respect to the number of transitions and thus perfectly correlated with the number of transitions. However, the time for synthesizing the controller depends on the fixed point algorithm (see [START_REF] Tabuada | Verification and control of hybrid systems: a symbolic approach[END_REF]) for which the worst case complexity is polynomial in the number of transitions which explains why the number of transitions and the CPU time are not perfectly correlated (a higher number of iterations is needed to reach the fixed point for r = 1 and r = 2 than for the other values). We can check that using the optimal sampling factor r = 3 allows us to reduce, for that example the computation times by about 70% in comparison to the approach corresponding to r = 1. For r = 3, Figure 2 shows a trajectory of the switched system and the associated switching signal controlled with the symbolic controller for the initial state x 0 = [0 3] T .

Conclusion

In this paper, we have proposed the use of multirate sampling for the computation of symbolic models for incrementally stable switched systems, with dwell-time constraints. We have demonstrated that our technique makes it possible to use more compact abstractions (i.e. with fewer transitions) than the standard existing approach presented in [START_REF] Girard | Approximately bisimilar symbolic models for incrementally stable switched systems[END_REF]. Moreover, the optimal sampling factor has been determined theoretically and we provided a simple expression depending solely on the number of modes and on the dimension of the state space, which makes it possible to use this result as a rule of thumb when computing symbolic models of switched systems. Our approach has been validated experimentally on a numerical example, which showed that multirate symbolic models indeed enable controller synthesis at a reduced computational cost. We are confident that similar results can be established for other classes of incrementally stable dynamical systems.

  deterministic, and metric when the set of outputs Y is equipped with the metric given by d Y (y, y ) = +∞ if y, y do not have the same dimension and d Y (y, y ) = y -y if y, y ∈ R n d Y (y, y ) = k max j=1 y j -y j if y, y ∈ R k×n with y = (y 1 , . . . , y k ), y = (y 1 , . . . , y k ).

( 19 )

 19 From Theorem 2.2 in[START_REF] Baker | Transcendental number theory[END_REF], we have that for all n, m ∈ N + with m ≥ 2, n ln(m) ∈ R + \ N + . Then, it follows that n ln(m)< n ln(m) and that 2 ln(m) τ d -ln(µ) κ > 0. Hence, τ > 0.

Fig. 1 .

 1 Fig. 1. Number of symbolic transitions in the multirate symbolic models T r τ,η (Στ d ) of the switched system with dwelltime for different values of the sampling factor r; Computation times for generating symbolic models and synthesizing safety controllers (MATLAB implementation, processor 2.2 GHz Intel Core i7, memory 16 GB 1600 MHz DDR3).

Fig. 2 .

 2 Fig.2. Top: Trajectory of the switched system with dwelltime and the associated switching signal controlled with the symbolic controller for the initial state x 0 = [0 3] T . The control period is τ = 0.2 while the sampling factor r = 3 (instants of transitions are indicated with circles). Bottom: Same trajectory in the state-space.