Yves Rozenholc

Markus Reiß

Remco Duits

Erik Franken

Bernhard Burgeth

Al

Markus Van Almsick

Yaniv Gur

Hanno Scharr

Kai 2 Krajsek

Luc Florack

Michael Felsberg

Jörg Polzehl

Karsten 3 Tabelow

Michal 4 Haindl

Markus Van

Preserving time structures while denoising a dynamical image

Keywords: Statistical Hypothesis Testing, Time Series, Denoising, Adaptive Estimation, DCE-CT, Microcirculation, Angiogenesis

In restoration or denoising of a movie, the classical procedures often do not take into account the full information provided by the movie. These procedures are either applied spatially "image per image" or locally on some neighborhood combining both closeness in one image and closeness in time. The local information is then combined homogeneously in order to realize the treatment. There is one type of movie where both approaches fail to provide a relevant treatment. Such a movie, called dynamical image, represents the same scene along the time with only variations, not in the positions of the objects in the scene but, in their gray levels, colors or contrasts. Hence, at each point of the image, one observes some regular temporal dynamics. This is the typical output using Dynamic Contrast Enhanced Computed Tomography (DCE-CT) or Dynamic Contrast Enhanced Magnetic Resonance Imaging (DCE-MRI) where at each spatial location (called voxel) of the image a time series is recorded.

In such a case, in order to preserve the full temporal information in the image, a proper denoising procedure should be based on averaging over spatial neighborhoods, but using the full dynamics of all pixels (or voxels) within the neighborhood. It is thus necessary to search homogeneous spatial neighborhoods with respect to the full dynamical information.

We introduce a new algorithm which meets these requirements. It is based on two main tools : a multiple hypothesis testing for the zero mean of a vector and an adaptive construction of homogeneous neighborhoods. For both tools, results of mathematical statistics ensure the quality of the global procedure. Illustrations from medical image sequences show a very good practical performance.

Introduction

Classically, the restoration or the denoising of a movie involve 2D techniques for smoothing [START_REF] Lepski | Optimal pointwise adaptive methods in nonparametric estimation[END_REF]Spokoiny, 1997, Buades et al., 2008], enhancing (see contributions on structure enhancement and use of line diffusion), filtering (see contribution of Michael Felsberg and references within), partitioning (see [START_REF] Polzehl | Adaptive weights smoothing with applications to image restoration[END_REF] or contributions of J. Polzehl and K. Tabelow), etc. applied image per image. If some methods use instead of one image a small number of images to take into account the time domain, all methods ignore the full time structure. A natural mathematical reason, outside the computational cost, is that the information contained in one image is considered to be less and less connected to the next images as the time increases. In some sense, the time regular structures existing in a movie do not have a long range effect and disappear quickly as the time distance between two images increases. One exception to this last rule is given by movies made from one dynamical image where at each 2D location (pixel or voxel) one observes not only a gray level or a color, but a full dynamical process represented by a time series or a vector. This is typically the case when the objects in the picture do not have changes in their position, but have changes only their brightness, contrast or color. Medical examples of interest are sequences like DCE-CT or DCE-MRI which register the variations of a contrast medium in the internal tissues of the body in response to a dynamical injection.

Our aim is to develop a new technique, adapted to these dynamical images, where each pixel presents a time series. In order to denoise by local averaging, we combine two approaches: statistical multiple hypothesis testing in the spirit of the works [Baraud et al., 2003, Durot andRozenholc, 2006] to compare the time series between two pixels, together with a spatial growth of locally homogeneous neighborhoods closely related to [Lepskiȋ, 1990, Lepski and[START_REF] Lepski | Optimal pointwise adaptive methods in nonparametric estimation[END_REF]. The originality of our method is not only to use the full dynamics to compare pixels, but also the way of sequentially growing the neighborhood. While usual paradigms to grow such neighborhoods compare the estimates built on the nested neighborhoods as in [START_REF] Lepski | Optimal pointwise adaptive methods in nonparametric estimation[END_REF], our method compares, at each step, independent estimates to stop or not the recursive growth. In the context of noise with heavy tails like Laplace noise for example, this comparison using independent estimates and robust statistics like the median offers a clear benefit in reducing the error.

Because the test procedure is based on the comparison to zero of the difference between two noisy vectors, i.e. two noisy time series, no modeling assumption is used. As we use an adaptive multiple hypothesis testing procedure, neither regularization in time is needed nor made. From a theoretical point of view, we only need the dynamics to be not too wild in the sense that the enhancement differences -viewed as a function of time -should remain in some Hölder regularity ball.

The article is organized as follows: In Section 1.2, we introduce the statistical framework of dynamical images and summarize our method. The two main statistical tools used to construct this method, multiple testing and neighborhood/ring growth, are respectively studied in 1.3 and in 1.4. In these sections, results from mathematical statistics are provided to guarantee the behavior of our method. Section 1.5 is devoted to the comparison of our estimation technique with estimates introduced in [START_REF] Lepski | Optimal pointwise adaptive methods in nonparametric estimation[END_REF] and in [START_REF] Polzehl | Adaptive weights smoothing with applications to image restoration[END_REF]. Finally in Section 1.6, we briefly describe an application to DCE-CT of renal metastasis in the liver, an extract from the specific work [Reiß et al., 2009a].

Denoising a dynamical image

We consider the following model called dynamical image indexed by the two quantities t and x representing time and spatial location, respectively:

I = {I x (t), x ∈ X , t ∈ {t 1 , t 2 , . . . , t K }} .
Here I x (.) denotes a noisy time series at location x. The notation I x will denote the vector (I x (t 1 ), . . . , I x (t K )), of the discrete observation at times t 1 ,. . . , t K of this time series. In our setting, X denotes the finite grid of pixels.

We assume that the noisy observations I x (t k ) may be written as

I x (t k ) = i x (t k ) + σε x k , (1.1)
where i x (t) denotes an unobservable true gray level at time t, where ε x k denotes a standardized noise with fixed distribution and where σ denotes the level of noise. We assume that the noise variables ε x k are independent and identically distributed with respect to both, space location x and (time index) k. It is not necessary to assume that σ and the distribution of ε x are known, but we will illustrate our construction using the simple setting σ known and ε x having a Gaussian distribution.

In order to compare the full time series, we use multiple hypothesis tests to decide whenever the vectors I x and I y at spatial locations x and y are statistically similar or not. For this purpose, we ask if the expectation of difference vector I x -I y is the zero vector or not at a given level α which could be considered as a natural tuning parameter. Such a test is presented in Section 1.3.

Given a set V of pixels, let us define the estimated dynamics using V as the coordinate-wise empirical mean

ÎV = 1 |V | y∈V I y , (1.2)
where |A| denotes the cardinality of a set A. Our aim is to construct at each spatial location x, a spatial neighborhood V x of x made from voxels y such that the difference I x -I y does not deviate significantly from the zero vector and such that the bias of the estimated dynamics obtained by replacing V in (1.2) by V x remains under control. Controlling the statistical bias means that we aim at using a procedure which ensures that the principal goal of nonparametric statistics, which is the trade-off between bias and variance, is achieved. This method is introduced in Section 1.4 and can be summarized as follows:

Let a geometrically increasing sequence (n i ) of positive integers be given with n 1 = 1. For a fixed location x ∈ X , we consider the set W x of all pixels y such that, at a level α ∈ (0, 1), I x -I y is accepted to have a zero mean by the previous test. We set i = 0 as well as V 0 = ∅ and W 0 = {x}. Our method repeats the following steps until one test, at least, rejects the null hypothesis:

• Increment i. • Define V i = V i-1 ∪W i-1 and build the set W i of the n i closest y in W x \V i . • For j = 1, . . . , i, test the i -1 null hypotheses E( ÎWi -ÎVj ) = 0 for j = 1, . . . , i -1 against a non-zero expectation.
The final estimate is ÎVi-1 . The algorithm realizes a pointwise adaptive selection of homogeneous dynamics. Because the last step involves a multiple hypothesis test, in order to ensure that its level is α, each individual test has to be calibrated with respect to α and i by a Bonferroni correction.

Multiple testing

Given two spatial locations x and y, we present the statistical test used to compare to the zero vector the expectation of the difference vector Z with components

Z k := I y (t k ) -I x (t k ), k = 1 . . . K.
Such a test of comparison to the zero vector derives from the theoretical works [Baraud et al., 2003, Durot andRozenholc, 2006] which consider a more general framework where σ is unknown and/or where ε is not necessarily assumed to be Gaussian. The theoretical study relies on large K asymptotics, which in our case means that we observe the continuous time signal more and more frequently. For the sake of simplicity, we write Z = f + √ 2σε and we introduce these tests in their simplified version which is the Gaussian case with σ known. We aim to test that the mean vector f is zero or not and hence consider the hypotheses

H 0 : "f = 0" versus H 1 : "f = 0".
To enlighten the notations, we suppose that K, the number of time indices, is of the form K = 2 d and we consider the regular dyadic decomposition of the observation times t 1 . . . t 2 d . For j = 0 . . . d -1, we denote by T j 1 . . . T j 2 j the 2 j intervals of time indexes

T j l = {t k : k = 2 d-j (l -1) + (1 . . . 2 d-j )}.
The set T j l contains 2 d-j time indices. Given j in 0 . . . d -1, let us introduce the projection of the observation Z onto the space spanned by the vectors with same coordinates on each time index T j l :

Π j Z = (m j 1 , . . . , m j 1 2 d-j times , . . . , m j 2 j , . . . , m j 2 j 2 d-j times ) with the local mean over T j l m j l = 1 2 d-j t∈T j l Z t .
The test is based on the size of the squared Euclidean norm Π j Z 2 K which is:

2 d-j 2 j l=1 m j l 2 = 1 2 d-j 2 j l=1   t∈T j l Z t   2 .
Under H 0 , the difference vector Z = I y -I x is a centered Gaussian vector with covariance matrix 2σ 2 Id K , where Id K denotes the identity matrix in R K . Hence, under H 0 , Π j Z 2 n /2σ 2 ∼ χ 2 (2 j ) which defines a distribution-free test statistic. Let us denote by Ψ -1 D the inverse cumulative distribution function of a χ 2 with D degrees of freedom, using a Bonferroni correction, our test procedure works as follows:

Reject H 0 at level α if for any j = 0 . . . d -1, Π j Z 2 K /2σ 2 > Ψ -1 2 j (α/d).
If σ is unknown, instead of a χ 2 -test based procedure, [START_REF] Baraud | Adaptive tests of linear hypotheses by model selection[END_REF]] proposes a Fisher like construction and considers the following procedure:

Reject H 0 at level α if for any j = 0 . . . d -1,

Π j Z 2 K /2 j Z -Π j Z 2 K /2 K-j > F -1 2 j ,2 K-j (α/d),
where F -1 2 j ,2 K-j denotes the inverse cumulative distribution function of a Fisher distribution with 2 j and 2 K-j degrees of freedom. Moreover, if the distribution of the noise variables ε x k is unknown, as ε := (ε xε y )/ √ 2 is by construction standardized and symmetrical, under the weak condition that the common distribution of the ε does not put mass in 0, we can follow the construction proposed in [START_REF] Durot | An adaptive test for zero mean[END_REF]:

Consider p Rademacher vectors R 1 , . . . , R p of length K. Rademacher vectors have independent coordinates R i with P (R i = +1) = P (R i = -1) = 1/2. Introduce the vectors N j = Z ⊗ R j where ⊗ denotes the coordinate-wise product. The previous constructions can be generalized by considering the following procedure:

Reject H 0 at level α if for any j = 0 . . . d -1,

Π j Z 2 K > Q -1 Π j Z ⊗ R 1 2 K , . . . , Π j Z ⊗ R p 2 K (α/d),
where Q -1 [X 1 , . . . , X K ] denotes the inverse cumulative distribution of the empirical distribution of the X i , i = 1, . . . , K.

Introducing the continuous-time signal F defined by F (t) := i x (t)i y (t), the difference vector Z = I x -I y follows the fixed-design regression model

Z k = f k + √ 2ε k = F (t k ) + √ 2σε k .
In this framework, the introduced multiple hypothesis tests are known to have a power which may be controlled on balls of the form f K ≥ ρ(K) where ρ(K) is a radius decreasing with K, the number of components of f (= number of observation times in the dynamical image). Moreover, these multiple hypothesis tests are adaptive with respect to the unknown Hölder regularity s of the unknown function F : at a fixed level α and for a given fixed power 1-β, this test automatically achieves (up to logarithmic factors) the minimax rate of testing ρ s (K) obtained in [START_REF] Gayraud | Adaptive minimax testing in the discrete regression scheme[END_REF] for all regularities s > 1/4. More precisely, we can derive following [Baraud et al., 2003, Th 1] and [Durot and Rozenholc, 2006, Th 4]) the following control of the power :

Theorem 1 : If the distribution of ε 1 does not put mass on 0, assuming that : (Bernstein's conditions) max 1≤i≤n E(ε 2p i ) ≤ γp!µ p-2 holds for all integers p ≥ 1 and some positive γ and µ , then there are A 1 , . . . , A 4 > 0 such that the multi-test procedure is of power 1β outside the ball defined by f 2 K ≤ inf j ∆(J j , f, A 1 , . . . , A 5 ) where

∆(J j , f, A 1 , . . . , A 5 ) = A 1 f -Π Jj f 2 K + A 2 γ β 2 (j-1)/2 + + A 3 ( √ γ + µ + max k=1..K (f 2 k )) 1 + j log 2 + log(1/β) 2 log 2d βα + + A 4 ( √ γ + µ + max k=1..K (f 2 k ))2 (j-1)/2 1 + j log 2 + log(1/β) 2 log 2d βα Corollary 1 If f k = F (t k
) then the test is adaptative w.r.t. the Hölder regularity s of F and achieves the optimal minimax rate for s > 1/4. This theorem is the exact translation of [Durot and Rozenholc, 2006, Thm. 4], using our simple setting with K = 2 d and dyadic partitions. It ensures a good power of the test (up to logarithmic terms) if the true unknown function to be tested is, at least for one dyadic partition, at a squared distance from 0 large enough to ensure that a noisy version has also its squared distance from 0 large enough. Let us recall that the squared distance of a noisy version is of the order of the sum of the variance term given by the number of pieces of the partition and the bias term f -Π Jj f 2 K which measures how well f can be approximated by a piecewise constant vector on the index partition.

We have now at hand the tool to compare dynamics at two locations x and y with respect to the known noise level σ. We denote "I y ≡ α σ 2 I x " when the enhancement difference average vector i xi y is accepted to be the zero vector at level α with respect to a noise variance 2σ 2 following the above construction.

Spatial neighborhood growth

We now present the construction of the spatial neighborhood V x for a fixed location x. In all that follows, the distance between spatial locations x in X can be measured by any metric ; usually the Euclidean metric will be taken.

The main idea of the sequential neighborhood growth can be explained in terms of testing. Suppose that neighborhoods

V 1 ⊂ V 2 ⊂ • • • ⊂ V M of x are given.
Then the stochastic fluctuations in the mean ÎVi of the observed signal over V i , as measured by the variances, decrease with i. On the other hand, the bias terms usually increase with i (we smooth with a larger bandwidth). We can assume that the bias in ÎV1 is negligible (or we just set V 1 = {x}). Considering next the empirical mean ÎV2\V1 over the observations in the set difference V 2 \ V 1 , a test is constructed to decide on the hypothesis that E( ÎV2\V1 -ÎV1 ) = 0 (recall that EZ denotes the coordinate-wise expectation of a vector Z). If the test rejects, we conclude that the bias of the new observations is already too large and we take V 1 as neighborhood, thus keeping ÎV1 as estimator. Otherwise, we accept V 2 as new neighborhood and consider V 3 \ V 2 . If ÎV3\V2 does not deviate significantly from ÎV1 and from ÎV2 , we continue. Otherwise, we stop and select V 2 as neighborhood, ÎV2 as estimator. This way we continue growing the neighborhoods until the first time the new observations differ statistically significantly from any of the estimators on the previously built neighborhoods or until we reach the last index M . This simple idea is reminiscent of Lepski's method [Lepskiȋ, 1991, Lepski and[START_REF] Lepski | Optimal pointwise adaptive methods in nonparametric estimation[END_REF]] in mathematical statistics with the important difference that at the ith step we do not test whether the new estimator ÎVi+1 is statistically homogeneous with all previous ÎVj , j ≤ i, but we compare only the average ÎVi+1\Vi over the new observations with the previous estimates. This has the clear advantage that the underlying test statistics are based on independent quantities and that in practice the number of too large grown neighborhoods is reduced because changes in the signal are detected earlier.

In order to correctly define our method, it is, of course, of great importance how to define precisely the tests and how to preselect the potential neighborhoods (V i ). In the static case (i.e. we only observe a scalar, not a time series at each point) this is described in full detail in [START_REF] Reiß | Pointwise adaptive estimation for robust and quantile regression[END_REF]], here we describe the main points. Let us first indicate how two averages ÎV and ÎV ′ over disjoint sets V and V ′ are compared. The test statistic is always based on the size of the difference Z V,V ′ = ÎV -ÎV ′ . In the static case and for Gaussian noise, we can standardize Z 2 V,V ′ by its expectation under the null, which is ρ(

V, V ′ ) = |V | -1 + |V ′ | -1
, and compare it with some critical value z V,V ′ . This critical value can be specified by a Bonferroni correction of the corresponding χ 2 -quantile, taking into account the multiple testing. A more refined and powerful, but also computationally more demanding calibration of the quantiles in this case is proposed in [START_REF] Reiß | Pointwise adaptive estimation for robust and quantile regression[END_REF]. Together with the prescription of geometrically growing neighborhoods (V i ), it is proved there that this results in estimators that optimally adapt to the local regularity of the signal (in the asymptotic minimax sense over Hölder balls). We thus achieve an unsupervised estimation procedure which is pointwise adaptive.

Here, the focus is on dynamic images and we need to compare the mean of time series over different spatial locations V and V ′ . For this we use exactly the adaptive multiple testing procedure presented in the preceding section. A second feature is that an a priori prescription of potentially good neighborhoods (V i ) around each pixel or voxel x is not feasible. The potential neighborhoods should rather be preselected based on the data observed in order to adapt to possibly heterogeneous geometric structures like borders or bands where a priori defined balls or squares are obviously not adequate. In a first preprocessing step we therefore exclude pixels that carry with high probability a significantly different time series. On the pixels kept after preprocessing we apply our neighborhood growth, which eventually results in the following algorithm.

We consider the set W x = {y ∈ X s.t. y = x and y ≡ α σ 2 x} of the spatial locations for which the time series are statistically similar to those of x with respect to the multi-test introduced in Section 1.3 above. Let us emphasize that W x is neither necessarily contiguous with {x} nor connected.

Our algorithm works as follows:

1. Start with i = 0 and set V 0 = {x} and Î0 = I x as the first neighborhood and its associated denoised time series. 2. Find W i the subset of the C 0 2 i closest points to V i in W x \ V i . We call this set the "ring" around V i , alluding to the case W x = X with the Euclidean distance. 3. Compute Ĵi the estimated time series using locations in the ring W i .

Test if

Ĵi is statistically not to be distinguished from all previously constructed estimates Îj for j = 0 . . . i. See hereafter for precision on this test. 5. If equality is accepted, consider

V i+1 = V i ∪ W i
and define the estimate Îi+1 using locations in V i+1 . Return to 2) with i = i + 1. 6. Else stop, define V x = V i and use Îx = Îi as denoised time series in x.

To compare Ĵi with the previously constructed estimates Î0 , . . . , Îi , we test the hypothesis H 0 : "E( Ĵi -Îj ) = 0, for all j = 0 . . . i", against H 1 : "E( Ĵi -Îj ) = 0, for at least one j", using a generalization of our test introduced in Section 1.3. We accept

H 0 if Ĵi ≡ α/(i+1) σ 2 ρ(|Vj |,|Wi|) Îj .
The use of α/(i + 1) is a Bonferroni correction due to multiple hypothesis tests. It would be possible to use more clever corrections as proposed in [START_REF] Benjamini | Controlling the false discovery rate: a practical and powerful approach to multiple testing[END_REF] for example. The correction factor ρ in the noise level takes into account that the estimates come from independent samples with respective size |V j | and |W i |, which in the Gaussian case is just the variance ρ(

|V j |, |W i |) = |V j | -1 + |W i | -1 .

Static toy examples

Not yet exploring the full scope of the method for treating dynamical images, we present first how the denoising is accomplished for one-and twodimensional signals (images) without an attached time series structure.

1D example

We propose a simple 1D example to compare the use of rings combined with balls with the algorithm introduced in [START_REF] Lepski | Optimal pointwise adaptive methods in nonparametric estimation[END_REF] based only on balls. Given an increasing sequence of neighborhoods V j , indexed by j = 1, . . . , j max , the selected index in [START_REF] Lepski | Optimal pointwise adaptive methods in nonparametric estimation[END_REF]] is defined by

 = inf j ≥ 0 s.t. ∃ℓ ≤ j, ÎVj -ÎV ℓ 2 E[ ÎVj -ÎV ℓ 2 ] > z 2 ℓ ∧ j max
while for our method the index is selected as follow

 = inf j ≥ 0 s.t. ∃ℓ ≤ j, ÎVj\V ℓ -ÎV ℓ 2 E[ ÎVj\V ℓ -ÎV ℓ 2 ] > z 2 ℓ ∧ j max
Both methods achieve the same rate of convergence for α-Hölder classes of signals in R d :

E[ ÎV  -i 2 ] ≤ c log n n 2α/(2α+d)
,

where n is the size of the neighborhood. This rate is known to be optimal in a minimax sense over classes of Hölder functions, already in the classical mean regression model with known errors (see [Lepskiȋ, 1991]). The logarithmic term has to be paid in order to be adaptive. Note that the dimension d, driven by the dimension of the neighborhood, is two for classical images, the length K of the time series itself does not change the rate, but will certainly influence the finite sample performance: more data in time increases the signal to noise ratio. In the left part of Figure 1.1 the true signal is presented together with a typical sample of 200 data points obtained by adding iid Laplace distributed noise. The estimation point is at x = 0 in the center. As possible balls, i.e. intervals, around x we take V k := [-a j , +a j ] with a j = 4 • (5/4) j /200 with j = 1, 2, . . . We compare the following methods : (1) balls as in Lepski [START_REF] Lepski | Optimal pointwise adaptive methods in nonparametric estimation[END_REF]] using the local mean; (2) balls/ring using the local mean; (3) balls as in Lepski [START_REF] Lepski | Optimal pointwise adaptive methods in nonparametric estimation[END_REF]] using the local median; (4) balls/ring using the local median; all these methods are compared to the (5) oracle where the median is taken over the (a posteriori optimal) interval [-0.39, +0.39], which is slightly larger than V 10 . We have used the empirical median in addition to the empirical mean (1.2) since for Laplace distributed noise the median is more efficient and robust. The calibration of the methods follows always the same algorithm as recommended in [START_REF] Reiß | Pointwise adaptive estimation for robust and quantile regression[END_REF]. The boxplots present the errors computed for the same 1000 samples with each method. The boxplots in Figure 1.1 show a comparable performance for our method in the mean case with the benchmark Lepski approach while in the median case a clear advantage for our ball/ring method can be seen. It is remarkable that the oracle error is almost achieved. The left subfigure shows a noisy image obtained using a pattern with two values and a Laplace noise with noise level equal to the difference of these two values. The middle subfigure shows the denoised version using AWS, a typical benchmark procedure for off-the-shelf spatially adaptive denoising. We have used the R package AWS version 1.3-3.1 developed by the authors of [START_REF] Polzehl | Adaptive weights smoothing with applications to image restoration[END_REF] available on the web site of the CRAN using the default 2D settings to produce this subfigure. Finally, the right subfigure shows the denoised image obtained by our method. While, at the used default setting, AWS offers a very clean picture in the right lower corner, all details are lost in the upper left corner. Using our method, the signal is less smoothed out, but allows to recover details even in the last part of the upper corner. Clearly more details can be found in this last image even if one could complain about the remaining noise due the noise level used. In Figure 1.3, four examples of selected neighborhoods are presented. One can remark the special effect of the first selection of the set W x in the presented examples: the neighborhood follows the borders and can even be not connected, allowing for rich geometrical structures in the signal.

A 2D example

Practical results

A complete study, using our technique in the framework of DCE-CT, has been conducted with radiologists at the European George Pompidou Hospital (HEGP), Paris, and the complete results are reported in [Reiß et al., 2009a]. In order to provide a flavor of how our method works on real data, we have extracted from this study a few illustrations. The precise setting of the medical experiments as well as the way the tuning parameters like σ and α are fixed is left to that work.

Due to the limited irradiation dose used during the sequential acquisition, the dynamical images in this study suffer from a poor signal-to-noise ratio. Figure 1.4 shows typical time series obtained from voxels in manually selected Regions of Interest (of size ≃100 voxels) within the Aorta (left curves) and a tumor (right curves).

The Figures 1.5 and 1.6 present the result of our method applied at two different voxels, one inside the aorta and the other inside the tumor.

Each figure is divided into 3 sub-figures. Left: the denoised slice at a specific time, with the selected voxel x (black cross) and the voxels from its neighborhood V x (pink dots). Center: the original enhancement vector I x (red thick curve) and the associated enhancement vectors I y for y ∈ V x (background curves). Right: the estimated enhancement using a generalized median as a robust center of the selected neighborhood (see Section 1.4).

Fig. 1

 1 Fig. 1.1 A simple example with a parabolic signal. Estimation point is x = 0 and noise is Laplace distributed. Errors from 1000 Monte Carlo replications

Figure 1 .

 1 Figure1.2 compares our method to AWS introduced in[START_REF] Polzehl | Adaptive weights smoothing with applications to image restoration[END_REF]. The left subfigure shows a noisy image obtained using a pattern with two values and a Laplace noise with noise level equal to the difference of these two values. The middle subfigure shows the denoised version using AWS, a typical benchmark procedure for off-the-shelf spatially adaptive denoising. We have used the R package AWS version 1.3-3.1 developed by the authors of[START_REF] Polzehl | Adaptive weights smoothing with applications to image restoration[END_REF] available on the web site of the CRAN using the default 2D settings to produce this subfigure. Finally, the right subfigure shows the denoised image obtained by our method. While, at the used default setting, AWS offers a very clean picture in the right lower corner, all details are lost in the upper left corner. Using our method, the signal is less smoothed out, but allows to recover details even in the last part of the upper corner. Clearly more details can be found in this last image even if one could complain about the remaining noise due the noise level used. In Figure1.3, four examples of selected neighborhoods are presented. One can remark the special effect of the first selection of the set W x in the presented examples: the neighborhood follows the borders and can even be not connected, allowing for rich geometrical structures in the signal.

Fig. 1

 1 Fig. 1.2 (left) Noisy data -(center) AWS -(right) Balls/Rings

  Fig. 1.4 Typical curves (time, enhancement) of individual voxels in ROI manually drawn within the Aorta (left) and the tumor (right)

Preserving time structures while denoising a dynamical image

Acknowledgment

This work was supported in 2008 by a grant "Bonus Qualité Recherche" from the University Paris Descartes for the project "Cancer Angiogénèse et Outils Mathématiques".


Yves Rozenholc
email: yves.rozenholc@parisdescartes.fr

Thoralf Mildenberger

Ursula Gather
email: gather@statistik.tu-dortmund.de

Combining Regular and Irregular Histograms by Penalized Likelihood

Keywords: irregular histogram, density estimation, penalized likelihood, dynamic programming

A new fully automatic procedure for the construction of histograms is proposed. It consists of constructing both a regular and an irregular histogram and then choosing between the two. To choose the number of bins in the irregular histogram, two different penalties motivated by recent work in model selection are proposed. A description of the algorithm and a proper tuning of the penalties is given. Finally, different versions of the procedure are compared to other existing proposals for a wide range of densities and sample sizes. In the simulations, the squared Hellinger risk of the new procedure is always at most twice as large as the risk of the best of the other methods.

The procedure is implemented in an R-Package.

Introduction

For a sample (X 1 , X 2 , . . . , X n ) of a real random variable X with an unknown density f w.r.t. Lebesgue measure, we denote the realizations by (x 1 , x 2 , . . . , x n ) and the realizations of the order statistics by x (1) ≤ x (2) ≤ • • • ≤ x (n) . The goal in nonparametric density estimation is to construct an estimate f of f from the sample. In this work, we focus on estimation by histograms, which are defined as piecewise constant densities. The procedure we propose consists of constructing both a regular and an irregular histogram (both to be defined below) and then choosing between the two. Although other types of nonparametric density estimators are known to be superior to histograms according to several optimality criteria, histograms still play an important role in practice. The main reason is their simplicity and hence their interpretability [START_REF] Birgé | How many bins should be put in a regular histogram?[END_REF]. Often, the histogram is the only density estimator taught to future researchers in non-mathematical subject areas, usually introduced in an exploratory context without reference to optimality criteria.

We first introduce histograms and describe the connection to Maximum Likelihood estimation: Given (x 1 , x 2 , . . . , x n ) and a set of densities F, the maximum likelihood estimate -if it exists -is given by an element f ∈ F that maximizes the likelihood n i=1 f (x i ) or equivalently its logarithm, the log-likelihood L(f, x 1 , . . . , x n ):

f := argmax f ∈F L(f, x 1 , . . . , x n ) := argmax f ∈F n i=1 log(f (x i )).
Without further restrictions on the class F, the log-likelihood is unbounded, and hence, no maximum likelihood estimate exists. One possibility is to restrict F to a set of histograms. Consider a partition I := {I 1 , . . . , I D } of a compact interval K ⊂ R into D intervals I 1 , . . . , I D , such that I i ∩ I j = ∅ for i = j and I i = K. Now consider the set F I of all histograms that are piecewise constant on I and zero outside I:

F I := f f = D j=1
h j 1 I j , h j ≥ 0, j = 1, . . . , D and

:= argmax f ∈F I L(f, x 1 , . . . , x n ) = 1 n D j=1 N j |I j | 1 I j , (1)
with N j = n i=1 1 I j (x i ). Its log-likelihood is

L( fI , x 1 , . . . , x n ) = D j=1 N j log N j n|I j | . (2)
In the following, we consider partitions I := I D := (I 1 , . . . , I D ) of the interval I := [x (1) , x (n) ], consisting of D intervals of the form

I j :=    [t 0 , t 1 ] j = 1 (t j-1 , t j ] j = 2, . . . , D , with breakpoints x (1) =: t 0 < t 1 < • • • < t D := x (n) . A histogram is called
regular if all intervals have the same length and irregular otherwise. The intervals are also referred to as bins.

We will only consider ML histograms in this work, and use the term "histogram" synonymously with "ML histogram" unless explicitly stated otherwise. We focus on finding a data-driven construction of a histogram with good risk behavior. Given a distance measure d between densities, the risk is defined as the expected distance between the true and the estimated density:

R n (f, fI , d) := E f [d(f, fI (X 1 , . . . , X n ))].
We consider the risks w.r.t. (normalized) squared Hellinger distance

d H (f, g) = 1 2 ( f (t) -g(t)) 2 dt, (3)
and w.r.t. powers of the L p -norms (p = 1, 2) defined by

d p := f -g p p = |f (t) -g(t)| p dt. (4)
For a detailed discussion on the choice of loss functions in histogram density estimation, see Birgé and Rozenholc (2006, sec. 2.2) and the references given there.

Given the sample, the histogram fI depends only on the partition I = (I 1 , . . . , I D ) as the values on the intervals of the partition are given by (1). In order to achieve good performance in terms of risk, the crucial point is thus choosing the partition. Since partitions with too many bins will result in a large likelihood without yielding a good estimate of f , a naïve comparison of the likelihoods of histograms with different numbers of bins is misleading. But also without any further restrictions on the allowed partitions the likelihood can be made arbitrarily large for a fixed number of bins.

Many approaches exist for the special case of regular histograms where I is divided into D equal sized bins; the problem is then reduced to the choice of D, cf. [START_REF] Birgé | How many bins should be put in a regular histogram?[END_REF], [START_REF] Davies | A comparison of automatic histogram constructions[END_REF] and the references given there. Using irregular partitions can reduce bias and therefore can improve performance for spatially inhomogenous densities, but the increased difficulty of choosing a good partition may lead to an increase in risk for more well-behaved densities. The idea of constructing both a regular and an irregular histogram and then choosing between the two is briefly discussed in [START_REF] Birgé | How many bins should be put in a regular histogram?[END_REF]. To our knowledge, this approach has not yet been put into practice.

Our recommendation is to construct a regular histogram that maximizes the penalized log-likelihood

L( fI , x 1 , . . . , x n ) -(D -1) -log 2.5 D (5)
among all regular partitions of [x (1) , x (n) ] with D = 1, . . . , ⌊n/ log n⌋ bins (where ⌊x⌋ denotes the largest integer not larger than x) and an irregular histogram that maximizes the penalized log-likelihood

L( fI , x 1 , . . . , x n ) -log n -1 D -1 -(D -1) -log 2.5 D (6)
among a set of partitions of [x (1) , x (n) ] with breakpoints equal to the sample points, where again D is the number of bins in a given partition. The final estimate is then the one with the larger penalized log-likelihood. The penalty in (5) for the regular case was proposed in [START_REF] Birgé | How many bins should be put in a regular histogram?[END_REF], while the motivation for ( 6) is developed later in this paper, where we consider different penalty forms for the irregular case. Note that the difference between the penalties is the term log n-1 D-1 which is needed because in the irregular case, the best partition with D bins has to be chosen from a set of n-1 D-1 partitions, while there is only one partition with D bins in the regular case. The necessity of taking into account in the penalty not only the number of parameters in a model but also the number of candidate models with same number of parameters is one of the key points in [START_REF] Barron | Risk bounds for model selection via penalization[END_REF]. Specific penalty forms for histogram estimators were derived in [START_REF] Castellan | Modified Akaike's criterion for histogram density estimation[END_REF] and for more general situations in [START_REF] Massart | Concentration inequalities and model selection[END_REF]. Note that both penalties (and hence the penalized log-likelihoods) coincide for D = 1. The penalties are both designed to achieve a good control on the Hellinger risk in term of Oracle Inequality as derived for example in (Massart, 2007, Th.??? eq(???)).Thoralf, I don't have the Massart book with me, could you extract an Oracle Inequality from a theorem connected to histogram Several methods for choosing a good irregular histogram have been developed previously. [START_REF] Kogure | Asymptotically optimal cells for a histogram[END_REF] gives asymptotic results for the optimal choice of bins. His approach is based on using blocks of equisized bins, and he explores the dependence on tuning parameters via simulations [START_REF] Kogure | Optimal cells for a histogram[END_REF]. It does not result in a fully automatic procedure. [START_REF] Kanazawa | An optimal variable cell histogram[END_REF] proposes to control the Hellinger distance between the unknown true density and the estimated histogram and introduces a dynamic programming algorithm to find the best partition with a given number of bins. [START_REF] Kanazawa | An optimal variable cell histogram based on the sample spacings[END_REF] derives the asymptotically optimal choice of the number of bins which depends on derivatives of the unknown density, making the procedure inapplicable in practice. [START_REF] Celisse | Nonparametric density estimation by exact leave-p-out cross-validation[END_REF] give explicit formulas for L 2 leave-p-out cross-validation for regular and irregular histograms. They only briefly comment on the case of irregular histograms and only show simulations with ad-hoc choices of the set of partitions. In our simulations, we use their explicit formula to compare risk behavior of cross-validation and our penalized likelihood approach when both are used to choose an irregular histogram from the same set of partitions. The multiresolution histogram by [START_REF] Engel | The multiresolution histogram[END_REF] is based on a tree of dyadic partitions. Its performance crucially depends on the finest resolution level, for which no universally usable recommendation is given. Some other tree-based procedures have been suggested for the multivariate case. They can be used for the univariate case, but they either perform a complete search over a restricted set of partitions [START_REF] Blanchard | Optimal dyadic decision trees[END_REF][START_REF] Klemelä | Multivariate histograms with data-dependent partitions[END_REF] or a greedy search over a full set of partitions [START_REF] Klemelä | Density estimation with stagewise optimization of the empirical risk[END_REF] to deal with computational problems that do not occur in the univariate case. Conditions for consistency of histogram estimates with data-driven and possibly irregular partitions are derived in [START_REF] Chen | Almost sure L 1 -norm convergence for databased histogram density estimators[END_REF]; [START_REF] Lugosi | Consistency of data-driven histogram methods for density estimation and classification[END_REF] and [START_REF] Zhao | Almost sure L r -norm convergence for data-based histogram estimates[END_REF]. [START_REF] Devroye | Bin width selection in multivariate histograms by the combinatorial method[END_REF] give a construction of histograms where bin widths are allowed to vary according to a pre-specified function. [START_REF] Hartigan | Bayesian histograms[END_REF] considers regular and irregular histogram construction from a Bayesian point of view. However, we are not aware of any fully tuned automatic Bayesian procedure for irregular histogram construction. [START_REF] Rissanen | Density estimation by stochastic complexity[END_REF] give a construction based on the Minimum Description Length (MDL) paradigm, which leads to a penalized likelihood estimator. A choice of several discretization parameters is needed, and the recommendation given by the authors is to perform an exhaustive search over all possible combinations of values, which is computationally expensive. A more recent MDL-based proposal by [START_REF] Kontkanen | MDL histogram density estimation[END_REF] involves a discretization which results in the estimate not being a proper density. [START_REF] Catoni | Data compression and adaptive histograms[END_REF] suggests a multi-stage procedure based on coding ideas that computes a density estimate by aggregating histograms.

The taut string procedure introduced by [START_REF] Davies | Densities, spectral densities and modality[END_REF] can also be used to generate an irregular histogram as described in [START_REF] Davies | A comparison of automatic histogram constructions[END_REF]. Regularization is achieved not by controlling the number of bins but by constructing an estimate that has a minimum number of modes subject to a constraint on the distance between the empirical and the estimated distribution function. The main idea is to construct a piecewise linear spline of minimal length (the taut string) in a tube around the empirical cdf and then take its derivative, which is piecewise constant. With some modifications this gives a histogram that fulfills definition (1). The main tuning parameter is the tube width, and an automatic choice is suggested by the authors.

Although not designed to minimize risk, the procedure has performed well w.r.t. classical loss functions [START_REF] Davies | A comparison of automatic histogram constructions[END_REF], and therefore is included in our simulations.

For our construction of irregular histograms, we will focus on penalized likelihood maximization techniques. For a good data-driven histogram one needs an appropriate penalization to provide an automatic choice of D as well as of the partition I = (I 1 , . . . , I D ). Since Akaike's Information Criterion (AIC) introduced by [START_REF] Akaike | A new look at the statistical model identification[END_REF], penalized likelihood has been used with many different penalty terms. AIC aims at ensuring a good risk behavior of the resulting estimate. Another widely used criterion is the Bayesian Information Criterion (BIC) introduced by [START_REF] Schwarz | Estimating the dimension of a model[END_REF]. It is constructed in such a way as to consistently estimate the smallest true model order, which in histogram density estimation would lead to very large models unless the true density is piecewise constant. In practice, criteria like AIC and BIC are routinely applied in many different statistical models, often without reference to their different conceptual backgrounds and without appropriate modifications for the model under consideration. In their original forms, both AIC and BIC do not account for multiple partitions with the same number of bins. See Chapter 7.3 of [START_REF] Massart | Concentration inequalities and model selection[END_REF] for a critique of the use of AIC in histogram density estimation. Since both are widely used, we include them in our comparisons. Our penalties are motivated by recent model selection works by [START_REF] Barron | Risk bounds for model selection via penalization[END_REF], [START_REF] Castellan | Modified Akaike's criterion for histogram density estimation[END_REF][START_REF] Castellan | Sélection d'histogrammes à l'aide d'un critère de type Akaike[END_REF] and [START_REF] Massart | Concentration inequalities and model selection[END_REF].

The regular histogram construction proposed in [START_REF] Birgé | How many bins should be put in a regular histogram?[END_REF] with which we combine our irregular histogram is based on the same ideas.

Our paper is structured as follows: In Section 2, we review the problem of constructing an irregular histogram using penalized likelihood. Section 3 gives a description of the practical implementation including calibration of the penalty. Section 4 gives the results of a simulation study and conclusions.

Penalized likelihood construction of irregular histograms

Constructing an irregular histogram by penalized likelihood means maximizing

L( fI , x 1 , . . . , x n ) -pen n (I), (7) w
.r.t. partitions I = (I 1 , . . . , I |I| ) of [x (1) , x (n) ]
, where pen n (I) is a penalty term depending only on the partition I and possibly on the sample (datadriven). We will introduce a new choice here motivated by work of [START_REF] Barron | Risk bounds for model selection via penalization[END_REF], [START_REF] Castellan | Modified Akaike's criterion for histogram density estimation[END_REF][START_REF] Castellan | Sélection d'histogrammes à l'aide d'un critère de type Akaike[END_REF] and [START_REF] Massart | Concentration inequalities and model selection[END_REF].

Optimizing w.r.t. the partition I with |I| fixed in (7) leaves us with a continuous optimization problem. Without further restrictions, for |I| ≥ 2 the likelihood is unbounded. One possibility is to restrict to all partitions that are built with endpoints on the observations; the optimization problem (7) can then be solved using a dynamic programming algorithm first used for histogram construction by [START_REF] Kanazawa | An optimal variable cell histogram[END_REF]. More details are given in Section 3.

With D = |I|, we propose the following families of penalties parametrized by two constants c and α:

pen A n (I) = c log n -1 D -1 + α(D -1) + ε (1) c,α (D), ( 8
) pen B n (I) = c log n -1 D -1 + α(D -1) + ε (2) (D) (9) and pen R n (I) = c log n -1 D -1 + α n D j=1 N j |I j | + ε (2) (D), (10)
where

ε (1) c,α (D) = c k log D + 2 cα(D -1)(log n -1 D -1 + k log D) (11)
and

ε (2) (D) = log 2.5 D. (12)
The precise choices for c and α obtained by simulations are described in Section 3. Note that, while the penalties given in ( 8) and ( 9) depend only on the number of bins of the partition, the penalty in formula ( 10) is a random penalty in the sense that it also depends on the data.

We now give arguments to explain the origins of these penalties. The penalty defined by ( 8) is derived from Theorem 3.2 in [START_REF] Castellan | Modified Akaike's criterion for histogram density estimation[END_REF], which is also stated as Theorem 7.9 in Massart (2007, p. 232) and from eq. (7.32) in Theorem 7.7 in Massart (2007, p.219). From the penalty form in Theorem 7.9 in [START_REF] Massart | Concentration inequalities and model selection[END_REF] we derive ε (1) :

pen n (I) = c 1 ( √ D -1 + √ c 2 x I ) 2 , ( 13
)
where the weights x I are chosen such that

D |I|=D e -x I ≤ Σ (14)
for an absolute constant Σ. Because the endpoints of our partitions are fixed, there are n-1 D-1 different partitions with cardinality D, and we assign equal weights to every partition I with |I| = D:

x I = log n -1 D -1 + k log D.
Choosing k > 1 ensures that the sum in ( 14) is converging and that Σ is finite. Substitution into (13) gives

pen n (I) = c 1 D -1 + c 2 log n -1 D -1 + k log D +2 c 2 (D -1) log n -1 D -1 + k log D . ( 15
)
Let us emphasize that Theorem 7.9 in Massart (2007, p. 232) requires

c 1 > 1/2 and c 2 = 2(1 + 1/c 1 ). Coming back to our notations, with α = c 1 , c = c 1 c 2 we obtain Equation (11).
We now use Theorem 7.7 in Massart (2007, p. 219) to justify the random penalty in (10). The orthonormal basis considered in this theorem for a given partition I consists of all 1 I / |I| for all I in I. The least squares contrast used in this theorem in our framework is -n -2 I∈I N 2 I /|I|. To link the minimization of the least squares contrast and the maximization of the log-likelihood, we consider the following approximation:

L( fI , x 1 , . . . , x n ) = D j=1 N j log N j n|I j | ≈ D j=1 N j N j n|I j | -1 = 1 n D j=1 N 2 j |I j | -n.
From the penalty form (7.32) and the use of M = 1 and ε = 0 in Theorem 7.7 in Massart (2007, p. 219), following the same derivation for ε (1) , we find the penalty in ( 15) with c 1 = 1 and c 2 = 2.

Using the least squares approximation, we can use the random penalty (7.33) in Theorem 7.7 in [START_REF] Massart | Concentration inequalities and model selection[END_REF]. Let us emphasize that V m defined by Massart is in our framework I∈I N I /n|I| with m = I. To derive ε (2) in (10) we start from the penalty defined in (7.33) in [START_REF] Massart | Concentration inequalities and model selection[END_REF]:

pen n (I) = (1 + ε) 5 V I + 2M L I D 2 .
Following the same derivations as for the penalty (13), setting M = 1, ε = 0 and L I = D -1 (log n-1 D-1 + k log D) we obtain:

pen n (I) = V I + 2 log n -1 D -1 + 2k log D +2 2 V I log n -1 D -1 + k log D .
Let us emphasize that, because of terms of the form ϕ(D) V I , the expression in the square root above prevents the use of dynamic programming to compute the maximum of the penalized log-likelihood defined in (7). To avoid this problem we propose, following penalty forms proposed in [START_REF] Birgé | How many bins should be put in a regular histogram?[END_REF] and [START_REF] Comte | A new algorithm for fixed design regression and denoising[END_REF], to replace the remainder expression 2k log D + 2 2 V I log n-1 D-1 + k log D by a power of log D. We have tried several values of the power and found that formula (12) leads to a good choice. Finally, we also replaced ε

(1) c,α in formula (8) by ε (2) , leading to the penalty given in (9).

Practical Implementation

We will describe briefly the implementation of our method. For a more detailed description, see [START_REF] Rozenholc | Combining regular and irregular histograms by penalized likelihood[END_REF]. To calibrate the penalties, histograms with the endpoints of the partitions placed on the observations and with different choices of the constants α and c in the penalties given in ( 8), ( 9) and ( 10) were evaluated by means of simulations. We used the same densities for calibration as in the simulations described in Section 4 but different samples and a smaller number of replications. The loss functions d = d H , d 1 , d 2 were evaluated by numerical integration using a trapeze rule.

We focused on the Hellinger risk to obtain good choices of the penalties, but the behavior w.r.t. L 1 loss is very similar. For minimizing L 2 risk, other choices may be preferable. Since no single penalty is best in all cases, we describe in the following what we consider to be a good compromise.

In formula (8) we tried different combinations of α ∈ {0.5, 1} and c between 1 and 4, some of which were motivated by Theorems 7.7 (eq. 7.32) and 7.9 in [START_REF] Massart | Concentration inequalities and model selection[END_REF]. We always set k = 2. From these experiments, the most satisfactory choice is c = 1 and α = 0.5. We also ran experiments replacing ε

(1) c,α by ε (2) , leading to the penalty given in (9). In this case, the most satisfactory choice is c = 1 and α = 1, and this choice is even better than ε

(1) 2,1 . Note that the resulting penalty, given in (6), exactly corresponds to the penalty in (5) proposed in [START_REF] Birgé | How many bins should be put in a regular histogram?[END_REF] for the regu-lar case, except for the additional term log n-1 D-1 that is needed to account for multiple partitions with the same number of bins. This term vanishes for D = 1, making the maxima of penalized likelihoods directly comparable. Because ( 8) and ( 9) are very similar, we only use this version in our simulations in Section 4.

For the random penalty in formula ( 10) we tried all combinations of c ∈ {0.5, 1, 2} and α ∈ {0.5, 1}. Let us emphasize that c = 2 and α = 1 correspond to formula (7.33) in [START_REF] Massart | Concentration inequalities and model selection[END_REF] up to our choice of ε (2) defined in (12). From our point of view, the most satisfactory choice is c = 1 and α = 0.5. In order to make the maximum of the log-likelihood penalized by (10) comparable to the maximum of (5), we add the constant α(x (n)x (1) ), which does not change the maximizer but makes the penalized log-likelihoods coincide for D = 1.

We now briefly describe the algorithm for constructing the irregular histogram. We consider partitions I built with endpoints on the observations:

I = ([x (1) , x (k 1 ) ], (x (k 1 ) , x (k 2 ) ], (x (k 2 ) , x (k 3 ) ], . . . , (x (k D-2 ) , x (k D-1 ) ], (x (k D-1 ) , x (n) ]),
where 1 < k 1 < . . . < k D-1 < n. We start from a "finest" partition I max defined by D max < n and the choice 1 < k 1 < . . . < k Dmax-1 < n. Our aim is to construct a sub-partition I of I max that maximizes (7). For all penalties (8)-( 10) this can be achieved by a dynamic programming algorithm [START_REF] Kanazawa | An optimal variable cell histogram[END_REF][START_REF] Comte | A new algorithm for fixed design regression and denoising[END_REF]. The total complexity of the algorithm is of order D 3 max . We reduce this to the order n by first using a greedy algorithm to construct a partition with ⌊max{n 1/3 , 100}⌋ bins if this number is smaller than n. Starting with a partition with one bin, in each step we recursively add to the current partition an endpoint (at an observation) in an existing bin in order to achieve the best maximization (in one step) of the likelihood. The resulting partition is used as the finest partition for the dynamic programming algorithm.

Let us remark that the theoretical results by [START_REF] Castellan | Modified Akaike's criterion for histogram density estimation[END_REF][START_REF] Castellan | Sélection d'histogrammes à l'aide d'un critère de type Akaike[END_REF] and Massart (2007, ch. 7), are derived for the case of a finest regular grid (depending on n but not on the sample) with bin sizes not smaller than a constant times log 2 (n)/n. However, we found that, in practice, we can improve performance drastically for some densities by using a data-dependent finest grid imposing no restrictions on the smallest bins without losing much at other densities. More comments on this are given in Section 4.

Simulation Study and Conclusions

The performance of our proposals given in Section 3 is compared to the performance of other available methods for 12 of the 28 test-bed densities introduced by Berlinet and [START_REF] Berlinet | A comparison of kernel density estimates[END_REF] which include standard distributions like the uniform and the normal as well as more special ones. We denote these by f 1 , . . . , f 12 . We also add 4 histogram densities f 13 , . . . , f 16 which are defined in [START_REF] Rozenholc | Combining regular and irregular histograms by penalized likelihood[END_REF]. Some of the densities do not fulfill the conditions in Theorem 3.2 in [START_REF] Castellan | Modified Akaike's criterion for histogram density estimation[END_REF] for example if they are not bounded away from zero. They are included in order to explore the behavior of the procedure in cases not covered by theory. All densities are implemented in the R-package benchden [START_REF] Mildenberger | benchden: 28 benchmark densities from Berlinet/Devroye[END_REF] and are depicted in Figure 1.

The sample sizes are 50,100,500,1000,5000 and 10000. We use 500 replications for each scenario. The methods compared in the simulations are: • B Penalized maximum likelihood using penalty (9) with c = 1 and

α = 1.
R Penalized maximum likelihood using random penalty (10) with c = 1 and α = 0.5.

CV Leave-one-out cross-validation using formula (11) given in [START_REF] Celisse | Nonparametric density estimation by exact leave-p-out cross-validation[END_REF]. We also tried formula ( 12) of [START_REF] Celisse | Nonparametric density estimation by exact leave-p-out cross-validation[END_REF] for different values of p without finding a big difference.

Maximization is performed over a data-driven finest grid as described in Section 3 without restrictions on the minimum bin width.

• Methods B, R, CV using the same data-driven grid but with the additional constraint that the minimum bin length allowed is (x (n) -

x (1) ) log 1.5 (n)/n. These are denoted by B ′ , R ′ and CV ′ .

• Methods B, R, CV but using a full optimization over a finest regular partition with bin width (x (n)x (1) )⌊n/ log 1.5 (n)⌋ -1 . This is the grid considered in [START_REF] Castellan | Modified Akaike's criterion for histogram density estimation[END_REF], except that we slightly relax log 2 (n) to log 1.5 (n). These are denoted B ′′ , R ′′ and CV ′′ .

• Penalized maximum likelihood using Information Criteria: AIC Akaike's Information Criterion [START_REF] Akaike | A new look at the statistical model identification[END_REF]. The penalty is

pen AIC n (D) = (D -1).
BIC Bayesian Information Criterion [START_REF] Schwarz | Estimating the dimension of a model[END_REF]. The penalty is pen BIC n (D) = 0.5 log(n)(D -1).

• The taut-string method TS introduced by [START_REF] Davies | Densities, spectral densities and modality[END_REF].

We use the function pmden() implemented in the R-package ftnonpar [START_REF] Davies | ftnonpar: Features and strings for nonparametric regression[END_REF] with the default values except that we set localsq=FALSE as local squeezing of the tube does not give a ML histogram. The histogram is then constructed using the knots of the string as the boundaries of the bins. This coincides with the derivative of the string except on intervals where the string switches from the upper to the lower boundary of the tube or vice versa.

• Regular histogram construction BR due to [START_REF] Birgé | How many bins should be put in a regular histogram?[END_REF].

The penalty is pen BR n (D) = (D-1)+log(D) 2.5 , where the log-likelihood is maximized over all regular partitions with 1, . . . , ⌊n/ log n⌋ bins.

• Combined approach: A regular histogram and an irregular histogram are constructed. After adding a constant such that all penalties are 0 for a histogram with 1 bin, the penalized log-likelihoods are compared. Such a procedure mixing several collections of models was proposed in [START_REF] Barron | Risk bounds for model selection via penalization[END_REF] but, to our knowledge, it has never been applied from a practical point of view. The aim is to achieve the best trade-off between bias and variance. Indeed, very large collections like those consisting of irregular histograms may lead to a small bias but may increase the stochastic part of the risk. The necessary price are terms of order log n-1 D-1 in the penalty which lead to larger upper bounds on the risk in the oracle inequality. It is then of interest to compare more than one strategy: one giving good results when the density is well-behaved (using the regular histogram collection) and one for less well-behaved (e.g. less smooth) densities (using the irregular histogram collection).

Because the density is unknown and so is its regularity, it is proposed to directly compare the penalized likelihoods in each collection. The histogram which has the larger penalized likelihood is chosen as the final estimate: B * Chooses between a regular histogram using BR and an irregular histogram using B R * Chooses between a regular histogram using BR and an irregular histogram using R.

Except for TS, all methods have been implemented in the R-package histogram (Mildenberger et al., 2009a). The program code to reproduce all simulation results is available from the corresponding author's website. For the discussion of the results, we focus on squared Hellinger risk. Table 1 shows the binary logarithms of relative squared Hellinger risks w.r.t. the best method for any given n and density: log 2 ( R method n / R best n ) for all 15 methods in the simulation study. Thus, a value of 0 means that the method was best in this particular setting and a value of 1 means that the risk of the method is twice as large as the risk of the best method. We omitted n = 5000 for space reasons; complete tables giving the absolute risk values, tables for other risks and a more detailed discussion of the results can be found in Rozenholc et al. (2009). Generally, the L 1 and Hellinger risks behave similarly, while the results for L 2 are much less conclusive and no clear recommendation for using one of the methods under consideration can be given when minimizing the L 2 risk is the main aim.

As was to be expected, the table shows no clear overall best method for all scenarios. The relative risk w.r.t. the best method is always smaller than 2 for our proposal B * . In this sense, it can be seen as the best method. Let us remark that other methods could be better in particular situations. In many cases, TS or one of our proposals B and R is either the best or the binary logarithm of relative risk w.r.t. the best method is close to zero. These three methods are -except for B * -also the only ones in the simulation study for which this quantity is always strictly smaller than log 2 3 ≈ 1.58, meaning that the empirical risk is never greater than three times the risk achieved by the best method. The random penalty R seems to be slightly better than B in many cases, the most notable exception being the trimodal uniform density f 10 for n = 50. Cross-validation using the same set of partitions performs rather poorly.

n B R CV B ′ R ′ CV ′ B ′′ R ′′ CV ′′ AIC BIC TS BR B * R * f 1 50 0.
Note that it is particularly bad for the uniform density f 1 . Relative performance of CV w.r.t. the best method becomes generally worse when sample size increases. If we compare our penalties and cross-validation for the case of full dynamic programming optimization over a finest regular partition with bin length (x (n)x (1) ) log 1.5 (n)/n (B ′′ , R ′′ ,CV ′′ ), the picture changes. Overall, the performance of all three methods is not bad, in particular CV ′′ often outperforms B ′′ and R ′′ , which behave very similarly. But CV ′′ performs badly for histogram densities, especially the uniform. Putting a constraint on the minimum bin size causes a problem for all three methods when the density has very sharp peaks (especially f 10 ). For the intermediate case,

i.e. using both penalties and cross-validation for a data-driven finest grid but adding a constraint on the bin widths, we see that B ′ and R ′ share the catastrophic behavior of B ′′ and R ′′ at f 10 without offering a real improvement over B and R at the other densities. On the other hand, CV ′ is a good compromise between CV and CV ′′ , but still shows bad behavior for f 1 and f 10 .

The taut string method TS shows a particularly good behavior in terms of Hellinger risk, although it was derived for different aims than performing well w.r.t. a given loss function [START_REF] Davies | Densities, spectral densities and modality[END_REF], and many questions regarding behavior in a more classical framework remain open. It does not control the number of bins but the modality of the estimate, thereby avoiding overfitting while still being able to chose a large number of bins to give sufficient detail.

Using AIC leads to very bad results. It is known to underpenalize even for regular histograms [START_REF] Birgé | How many bins should be put in a regular histogram?[END_REF][START_REF] Castellan | Modified Akaike's criterion for histogram density estimation[END_REF][START_REF] Massart | Concentration inequalities and model selection[END_REF]. This tendency becomes even worse when used for irregular histograms, since now there are many partitions with the same number of bins. This leads to problems similar to those arising in multiple testing. In this case, an additional penalty is needed [START_REF] Castellan | Modified Akaike's criterion for histogram density estimation[END_REF][START_REF] Massart | Concentration inequalities and model selection[END_REF]. BIC does not aim for a good control of risk but at asymptotically identifying the "smallest true model", if it exists. It shows some good behavior in particular for small sample sizes that deteriorates when samples become larger. Particularly noteworthy is the bad performance for "simple" models like the uniform and the 5 bin regular histogram density f 13 .

The regular histogram method BR, which improves on Akaike's penalization, is the best method for f 3 , f 11 and f 13 , at least when the sample size is not very small. This shows that the greater flexibility of an irregular

histogram over a regular one may be outweighed by the greater difficulty in choosing a good partition, as was already remarked by [START_REF] Birgé | How many bins should be put in a regular histogram?[END_REF]. Regular histograms are inferior for spatially inhomogeneous densities like f 4 and f 10 .

The simulations show that one can successfully combine the advantages of regular and irregular histograms (B * and R * ). The Hellinger risk for B * is always within twice the risk of the best method for a given situation. While R * is even better in most cases, it shares the bad behaviour of BR for f 10 and n = 50. Table 2 shows that for larger sample sizes both B * and R * almost always choose an irregular histogram for densities where this is advantageous (the spatially inhomogenous densities f 4 , f 5 , f 10 , f 14 and f 16 ) and a regular partition in most of the other cases.

To summarize, we propose a practical method of irregular histogram construction inspired by theoretical works by [START_REF] Barron | Risk bounds for model selection via penalization[END_REF][START_REF] Castellan | Modified Akaike's criterion for histogram density estimation[END_REF][START_REF] Castellan | Sélection d'histogrammes à l'aide d'un critère de type Akaike[END_REF] and [START_REF] Massart | Concentration inequalities and model selection[END_REF]. It can be easily implemented using a 0.02 0.13 0.10 0.68 1.00 0.11 0.27 0.06 0.12 0.00 0.02 0.16 0.03 0.07 0.22 0.43 100 B * 0.04 0.07 0.04 0.74 1.00 0.11 0.23 0.02 0.19 1.00 0.12 0.43 0.01 0.12 0.24 0.65 R * 0.00 0.10 0.10 0.64 1.00 0.24 0.40 0.05 0.36 1.00 0.05 0.28 0.01 0.11 0.33 0.65 500 B * 0.03 0.00 0.00 0.73 1.00 0.01 0.03 0.00 0.42 1.00 0.00 0.47 0.00 0.69 0.36 0.99 R * 0.00 0.00 0.00 0.47 1.00 0.04 0.07 0.00 0.53 1.00 0.00 0.45 0.00 0.75 0.48 0.99 1000 B * 0.02 0.00 0.00 0.89 1.00 0.00 0.00 0.00 0.58 1.00 0.00 0.41 0.00 0.96 0.49 1.00 R * 0.01 0.00 0.00 0.63 1.00 0.01 0.00 0.00 0.64 1.00 0.00 0.41 0.00 0.98 0.57 1.00 5000 B * 0.01 0.00 0.00 1.00 1.00 0.00 0.00 0.00 0.65 1.00 0.00 0.26 0.00 1.00 0.43 1.00 R * 0.00 0.00 0.00 1.00 1.00 0.00 0.00 0.00 0.70 1.00 0.00 0.28 0.00 1.00 0.47 1.00 10000 B * 0.00 0.00 0.00 1.00 1.00 0.00 0.00 0.00 0.66 1.00 0.00 0.18 0.00 1.00 0.35 1.00 R * 0.00 0.00 0.00 1.00 0.98 0.00 0.00 0.00 0.71 1.00 0.00 0.20 0.00 1.00 0.39 1.00

Table 2: Frequency of choosing an irregular partition.

dynamic programming algorithm and it performs well for a wide range of different densities and sample sizes, even for some cases not covered by the underlying theory. Performance is shown to be improved when combined with the regular histogram approach proposed in [START_REF] Birgé | How many bins should be put in a regular histogram?[END_REF]. All procedures proposed here are available in the R-package histogram (Mildenberger et al., 2009a).

  j | = 1 , where 1 A denotes the indicator function of a set A and |I| the length of the interval I. If K contains [x (1) , x (n) ], the Maximum Likelihood Histogram (ML histogram) is defined as the maximizer of the log-likelihood in F I and is given by fI

Figure 1 :

 1 Figure 1: The densities used in the simulation study.

Table 1 :

 1 Binary logarithms of relative squared Hellinger Risks w.r.t the best method.

			13 0.05 1.53 0.03 0.03 0.76 0	0	0.44 2.07 1.27 0.07 0.12 0.18 0.13
		100	0.08 0.02 2.31 0.02 0.02 1.25 0	0.01 0.72 2.91 1.56 0.1	0.11 0.17 0.12
		500	0.06 0.03 4.2	0.01 0.02 2.56 0	0.01 1.69 4.64 2.04 0.19 0.13 0.17 0.14
		1000	0.04 0.02 4.64 0.01 0.01 3.05 0	0.01 2.3	4.98 2.1	0.12 0.11 0.14 0.12
		10000	0	0	5.39 0	0	3.98 0	0	4.71 5.59 1.85 0.2	0.18 0.19 0.18
	f 2	50 100	0.42 0.34 0.5 0.57 0.51 0.81 0.56 0.51 0.05 0.56 0.54 0 0.38 0.34 0.09 0.52 0.47 0	0.91 0.53 0.13 0.14 0.17 0.16 1.34 0.61 0.02 0.2 0.24 0.25
		500	0.71 0.67 1.79 0.7	0.66 0.04 0.69 0.67 0.07 2.2	0.75 0	0.2	0.2	0.2
		1000	0.74 0.7	1.94 0.72 0.7	0.08 0.72 0.7	0.15 2.25 0.73 0	0.18 0.18 0.18
		10000	0.69 0.71 1.59 0.69 0.71 0.21 0.67 0.68 0.68 1.7	0.56 0	0.07 0.07 0.07
	f 3	50 100	0.51 0.42 0.68 0.49 0.43 0.14 0.51 0.48 0 0.47 0.4 1.07 0.45 0.39 0.17 0.48 0.43 0	1.15 0.64 0.51 0.23 0.24 0.23 1.62 0.72 0.25 0.15 0.16 0.17
		500	0.65 0.6	2.08 0.64 0.6	0.36 0.65 0.59 0.1	2.49 0.81 0.13 0	0	0
		1000	0.73 0.69 2.28 0.73 0.69 0.54 0.73 0.66 0.25 2.59 0.84 0.19 0	0	0
		10000	1.03 0.99 2.19 1.02 0.99 0.94 1	0.96 1.34 2.32 0.93 0.57 0	0	0
	f 4	50 100	0.56 0.5 0.63 0.6	0.68 0.74 0.72 0.45 0.68 0.65 0.41 1.12 0.67 0 1.08 0.72 0.72 0.38 0.7 0.69 0.38 1.62 0.79 0	0.55 0.56 0.52 0.61 0.62 0.59
		500	0.56 0.61 1.79 0.71 0.72 0.36 0.66 0.69 0.41 2.2	0.62 0	0.67 0.55 0.62
		1000	0.31 0.42 1.74 0.72 0.74 0.4	0.7	0.73 0.46 2.05 0.39 0	0.75 0.34 0.52
		10000	0.55 0.77 1.62 1.48 1.56 1.27 1.27 1.37 1.25 1.73 0.45 0	2.05 0.55 0.77
	f 5	50 100	0.25 0.19 0.5 0.47 0.44 0.76 2.71 2.71 2.73 1.81 1.82 1.77 1.07 0.48 0 2.18 2.18 2.18 1.6 1.59 1.57 0.64 0.33 0	1.12 0.25 0.19 0.41 0.47 0.44
		500	0.79 0.77 1.73 1.69 1.69 1.92 1.89 1.9	1.68 2.05 0.81 0	1.18 0.79 0.77
		1000	0.89 0.91 1.92 1.56 1.56 1.78 1.29 1.3	0.98 2.16 0.85 0	0.95 0.89 0.91
		10000	0.98 1	1.63 0.98 1	0.47 0.96 1	0.22 1.71 0.82 0	0.68 0.98 1
	f 6	50 100	0.4 0.36 0.29 1.14 0.35 0.3 0.31 0.75 0.38 0.32 0.24 0.33 0.3 0.31 0.45 0.39 0 0	1.22 0.66 0.51 0.28 0.3 1.7 0.8 0.44 0.31 0.31 0.29 0.27
		500	0.5	0.45 1.98 0.5	0.45 0.31 0.48 0.42 0	2.38 0.71 0.1	0.09 0.09 0.1
		1000	0.55 0.49 2.11 0.54 0.49 0.42 0.52 0.47 0.11 2.42 0.67 0.06 0	0	0.01
		10000	0.92 0.87 2.1	0.91 0.87 0.85 0.91 0.86 1.26 2.23 0.83 0	0.11 0.11 0.11
	f 7	50 100	0.28 0.16 0.48 0.25 0.17 0.04 0.39 0.35 0 0.05 0.01 0.66 0.02 0 0.08 0.2 0.18 0.05 1.15 0.53 0.06 0.14 0.12 0.08 0.86 0.45 0.39 0.19 0.2 0.15
		500	0.68 0.64 1.44 0.66 0.63 0.14 0.69 0.66 0.15 1.84 0.68 0	0.27 0.28 0.3
		1000	0.89 0.87 1.8	0.88 0.86 0.44 1.15 1.14 0.29 2.1	0.8	0	0.45 0.45 0.45
		10000	1.25 1.24 1.85 1.25 1.24 0.76 1.26 1.26 1.11 1.94 1.11 0	0.77 0.77 0.77
	f 8	50 100	0.3 0.37 0.28 0.27 0.37 0.3 0.24 0.08 0.29 0.26 0.05 0.34 0.32 0 0 0.63 0.55 0.09 0.67 0.27 0.17 0.09 0.09 0.09 0.4 0.21 0.22 0.14 0.16 0.14
		500	0.44 0.36 0.9	0.51 0.46 0.24 0.66 0.59 0	1.25 0.35 0.03 0.1	0.1	0.1
		1000	0.47 0.39 1.05 0.54 0.49 0.27 0.64 0.58 0.09 1.33 0.37 0	0.19 0.19 0.19
		10000	0.99 0.93 1.17 1.01 0.96 0.61 0.99 0.93 0.66 1.25 0.82 0	0.7	0.7	0.7
	f 9	50 100	0.44 0.42 0 0.46 0.28 0.72 0.49 0.28 0 0.45 0.45 0.16 0.44 0.44 0.36 0.34 0.03 0.26 0.36 0.38 0.37 1.37 1.34 1.09 1.13 0.47 0.17 0.6 0.51 0.37
		500	0.55 0.5	1.81 0.52 0.48 0.36 1.64 1.62 1.36 2.17 0.68 0	0.85 0.63 0.55
		1000	0.61 0.56 2.02 0.59 0.55 0.5	1.65 1.62 1.37 2.3	0.72 0	1.08 0.66 0.61
		10000	0.73 0.69 1.87 0.72 0.69 0.74 1.61 1.59 1.69 1.98 0.65 0	1.31 0.77 0.74
	f 10	50 100	0 0.06 0 1.31 0.55 3.15 3.15 3.14 3.12 3.12 3.12 0.68 0.31 0.15 2.68 0 1.02 4.2 4.2 4.2 4.02 4.02 4.02 1.32 0.53 0.32 3.29 0.06 0 2.68
		500	0.02 0	2.4	5.82 5.82 5.82 5.64 5.64 5.64 2.74 0.79 0.51 5.08 0.02 0
		1000	0.01 0	2.8	6.8	6.8	6.8	6.34 6.34 6.34 3.07 0.82 0.56 5.63 0.01 0
		10000	0	0	3.48 10.1 10.1 10.1 6.46 6.46 6.46 3.59 0.65 0.81 0.17 0	0
	f 11	50 100	0.06 0.02 0.54 0.01 0.01 0.22 0 0.02 0.01 0.39 0 0 0.24 0	0 0	0.13 0.81 0.51 0.19 0.03 0.08 0.04 0.13 0.74 0.43 0.25 0.03 0.05 0.03
		500	0.86 0.61 0.8	0.89 0.7	0.27 1.06 1.06 0.15 1.14 0.48 0.11 0	0	0
		1000	0.63 0.58 0.98 0.62 0.59 0.4	0.86 0.69 0.08 1.24 0.63 0.09 0	0	0
		10000	1.06 0.99 0.93 1.05 0.98 0.57 1.09 1.01 0.37 0.97 0.91 0.17 0	0	0
	f 12	50 100	0.21 0.11 0.88 0.01 0.02 0.25 0 0.23 0.17 1.17 0.13 0.14 0.27 0.12 0.13 0 0.01 0.14 1.36 0.79 0.03 0.15 0.27 0.18 1.75 0.73 0.15 0.16 0.24 0.18
		500	0.37 0.33 2.41 0.3	0.3	0.88 0.29 0.28 0.28 2.83 0.9	0	0.29 0.37 0.34
		1000	0.52 0.49 2.63 0.48 0.47 1.09 0.52 0.51 0.5	2.95 0.89 0	0.36 0.45 0.42
		10000	0.68 0.66 2.48 0.67 0.66 1.18 0.67 0.65 1.61 2.62 0.74 0	0.57 0.6	0.59
	f 13	50 100	0.1 0.43 0.4 0.07 0.69 0.04 0.05 0.02 0.01 0.02 0 1.21 0.39 0.38 0.25 0.38 0.38 0.04 1.77 0.87 0.4 1.16 0.65 0.06 0.06 0.09 0.06 0 0.01 0.01
		500	1.41 1.34 3.43 1.39 1.33 1.86 1.53 1.5	1.41 3.85 1.88 1.1	0	0	0
		1000	1.06 0.98 3.98 1.02 0.98 2.44 0.83 0.78 1.7	4.29 1.97 1.11 0	0	0
		10000	0.92 0.92 4.83 0.9	0.91 3.57 1.39 1.39 4.1	5	1.93 1.25 0	0	0
	f 14	50 100	0.18 0.13 0.65 0.12 0.11 0.14 0.12 0.13 0 0.26 0.23 1 0.24 0.23 0.3 0.36 0.33 0.19 1.54 0.65 0.29 0 1.1 0.59 0.19 0.02 0.08 0.04 0.03 0.02
		500	0.08 0.03 2.41 0.02 0	0.85 0.58 0.56 0.67 2.83 0.83 0.09 0.88 0.35 0.28
		1000	0.14 0.11 2.98 0.12 0.11 1.45 0.83 0.82 1.25 3.3	0.98 0	0.94 0.21 0.15
		10000	0.01 0.01 3.96 0	0	2.67 0.92 0.93 3.31 4.13 1.04 0.57 2.64 0.01 0.01
	f 15	50 100	0.1 0.27 0.18 0.65 0.24 0.21 0 0 0.06 0.35 0.29 0.06 0.46 0.42 0.06 0.43 0.09 0.11 0.03 0.04 0.01 0.69 0.65 0.35 1.11 0.46 0.1 0.27 0.25 0.2
		500	0.15 0	1.79 0.15 0.01 0.31 1.15 1.09 0.87 2.17 0.48 0.21 0.48 0.27 0.13
		1000	0.06 0.01 2.39 0.03 0	0.87 1.26 1.25 1.32 2.69 0.67 0.34 0.87 0.26 0.16
		10000	0.01 0	2.8	0.01 0	1.69 1.5	1.48 2.59 2.95 0.51 0.1	0.85 0.14 0.1
	f 16	50 100	0.11 0 0.1 0.03 0.72 0 0.29 0.96 0.93 0.74 0.74 0.71 0.48 0.65 0.29 0 0.01 0 0.9 0.91 0.78 1.2 0.46 0.02 0.55 0.25 0.21 0.34 0.26 0.2
		500	0.07 0.02 1.8	0.58 0.57 0.79 1.53 1.52 1.32 2.19 0.5	0	1.21 0.09 0.03
		1000	0	0	2.18 0.27 0.28 0.84 1.21 1.22 1.26 2.49 0.55 0.03 1.24 0	0
		10000	0.04 0	2.97 0.04 0	1.83 1.17 1.15 2.58 3.12 0.59 0.8	1.65 0.04 0

Acknowledgments

This work has been supported in part by the Collaborative Research Centers "Reduction of Complexity in Multivariate Data Structures" (SFB 475) and "Statistical Modelling of Nonlinear Dynamic Processes" (SFB 823) of the German Research Foundation (DFG). The authors also wish to thank Henrike Weinert for discussions and programming in earlier stages of the work as well as two anonymous referees for giving comments that greatly improved the paper.