Valentina Zelaya Mendizábal
email: valentina.zelaya@etu.univ-paris1.fr

Marc Boullé
email: marc.boulle@orange.com

Fabrice Rossi
email: fabrice.rossi@dauphine.psl.eu

Fast and fully-automated histograms for large-scale data sets

Keywords: Density estimation, Histograms, Model selection, Minimum description length

G-Enum histograms are a new fast and fully automated method for irregular histogram construction. By framing histogram construction as a density estimation problem and its automation as a model selection task, these histograms leverage the Minimum Description Length principle (MDL) to derive two different model selection criteria. Several proven theoretical results about these criteria give insights about their asymptotic behavior and are used to speed up their optimisation. These insights, combined to a greedy search heuristic, are used to construct histograms in linearithmic time rather than the polynomial time incurred by previous works. The capabilities of the proposed MDL density estimation method are illustrated with reference to other fully automated methods in the literature, both on synthetic and large real-world data sets.

Introduction

Histograms are popular non-parametric density estimators available in all statistical computing packages. They are particularly adapted for univariate visualisation (see for instance [START_REF] Zubiaga | Graphical perception of value distributions: An evaluation of non-expert viewers' data literacy[END_REF]) and can serve as the starting point for more complex analyses. Histograms are also quite useful because they provide a compressed representation of the distribution of a random variable. This property has been used for instance in database query optimization, as illustrated in [START_REF] Oommen | The Efficiency of Histogram-like Techniques for Database Query Optimization[END_REF]; [START_REF] Ioannidis | The history of histograms (abridged)[END_REF].

In theory, histograms require very few parameters and can adapt to any kind of distribution given enough bins. In practice however, the choice of the binning can have an unforeseeable effect on the accuracy of data estimation, on the readability of the visualisation and on the size of the representation.

Numerous strategies exist to automatically select the number of bins in a regular histogram with equal width bins (for instance the pioneering Sturges' formula [START_REF] Sturges | The choice of a class interval[END_REF] and the Freedman-Diaconis' rule [START_REF] Freedman | On the histogram as a density estimator:l2 theory[END_REF]). A major problem with regular histograms is that regions with higher or lower densities are treated the same, which gives histograms their limited reputation. For complex distributions, irregular histograms with bins of different widths are more adapted [START_REF] Rissanen | Density estimation by stochastic complexity[END_REF], but fully automated construction methods are scarce or limited, as pointed out in [START_REF] Davies | A comparison of automatic histogram constructions[END_REF]; [START_REF] Rozenholc | Combining regular and irregular histograms by penalized likelihood[END_REF]: many methods rely on user adjustable parameters, while others have high computational loads that hinder their use on large scale data. Scalable and fully automated approaches that specify the location, number and widths of histogram bins, based only on the observed properties of the data, are still rare overall.

Although alternative methods of density estimation, such as kernel density estimators, are usually recommended as a fitting solution to this drawback, we argue that histograms can stay relevant for density estimation if the choice of binning is done properly, especially because of the ease with which they are interpreted [START_REF] Zubiaga | Graphical perception of value distributions: An evaluation of non-expert viewers' data literacy[END_REF]. Most fully automated methods in the literature view histogram construction as a model selection task and implement it via some form of penalized quality criterion that balances the quality of the histogram as a representation of the data with its complexity. Among these, the Minimum Description Length principle (MDL [START_REF] Rissanen | Modeling by shortest data description[END_REF][START_REF] Grunwald | The minimum description length principle[END_REF]) provides a sound general framework to implement model selection. The key idea of MDL is that any regularity in the data can be used to "compress the data", i.e. to describe it in a shorter manner. More formally, the best choice for a model and its parameters is the one that minimizes the coding length of the model parameters and of the data given the model.

MDL has been applied to histogram construction, using different derivations of the coding length. For instance, the Bayesian Mixture criterion and a uniform prior were used to construct irregular histograms in [START_REF] Rissanen | Density estimation by stochastic complexity[END_REF]. Another work formalized a MDL criterion via the Normalized Maximum likelihood (NML) distribution -which, unlike the Bayesian Mixture criterion, has several desirable optimality properties (Kontkanen and Myllymäki, 2007). Notice however that the NML criterion for histograms relies on a single user-defined parameter, the accuracy at which the data is to be approximated, . In addition, while the experimental results reported on Gaussian distributions are very good, this approach has some scalability limitations. These issues are consequences of the computational complexities of the NML criteria evaluation and of the search for the optimal model in the space of all histograms.

We introduce a new MDL based histogram construction method that aims to reduce the computational burden of previous solutions and to enable automatic tuning of the accuracy parameter . Our formulation is based on a Bayesian maximum a posteriori interpretation of MDL. The criterion derived from this formulation does not contain the NML parametric complexity term from Kontkanen and Myllymäki (2007) which enables both a simpler analysis and a faster evaluation. Leveraging the theoretical properties of the proposed criterion, we derive a simple search heuristic over the space of histograms with a linearithmic time complexity (O(n log n)) rather than a polynomial one. As previous MDL based criteria, our criterion uses an accuracy parameter. We study the effect of this parameter on the histograms obtained by optimising our criterion. Then we introduce a new granulated criterion, derived from the base one, that automates the choice of .

The rest of this paper is organized as follows. Section 2 provide a short overview of other fully automated histogram construction methods. Section 3 specifies in details the histogram construction problem. Section 4 describes the NML approach from Kontkanen and Myllymäki (2007) and discusses its limitations. Section 5 introduces two enumerative criteria for histogram construction, while their theoretical properties and consistency are discussed in section 6. An efficient optimisation algorithm that benefits from these properties is then presented in section 7. In section 8, experiments demonstrate the performance of all three criteria and other state-of-the-art methods on synthetic and real large-scale datasets. Concluding remarks are given in section 9.

Related work

As pointed out in [START_REF] Birge | How many bins should be put in a regular histogram[END_REF] and [START_REF] Davies | A comparison of automatic histogram constructions[END_REF], among others, almost all non-heuristic automatic histogram construction is based on a notion of risk. The goal is to minimise

R n ( f, fθ , l) = E f l f, fθ (x 1 , . . . , x n ) , (1)
where f is the true unknown density of the data, x 1 , . . . , x n a sample from f , l a loss function, and fθ the estimation procedure with its parameters θ (e.g. the number of bins in a regular histogram). An ideal choice of θ would minimise the risk. Most statistical computing software still include simple methods, such as the Sturges rule [START_REF] Sturges | The choice of a class interval[END_REF], which divides the data domain into K * = 1 + log 2 n intervals. Other more principled approaches derived from asymptotic analyses of the risk are also included. For instance, Scott's rule [START_REF] Scott | On optimal and data based histograms[END_REF] is derived from the asymptotic behaviour of the risk for the squared L2-loss l( f, g) = fg 2 2 . Notice that this type of analysis relies on smoothness assumptions on the true density f , which are used to derive the optimal bin width from characteristics of this unknown density. Those characteristics are in turn estimated from the data either using the Normal distribution as a reference as in Scott's rule, using heuristic consideration as in Freedman-Diaconis' rule [START_REF] Freedman | On the histogram as a density estimator:l2 theory[END_REF] or based on plug-in methods as in [START_REF] Wand | Data-based choice of histogram bin width[END_REF]. See [START_REF] Sulewski | Equal-bin-width histogram versus equal-bin-count histogram[END_REF] and the references therein for other examples of such approaches. Notice however that none of these approaches can be used to build irregular histograms. As pointed out earlier, being able to have intervals of varying widths is advantageous to describe denser and narrower ranges of values.

An alternative to asymptotic analyses is to leverage the cross-validation (CV) principle to directly estimate the risk. For instance [START_REF] Rudemo | Empirical choice of histograms and kernel density estimators[END_REF] proposes to use leave-one-out estimates of the risk. Those can be computed in closed-form for the squared L2-loss and the Kullback-Leibler loss (see [START_REF] Hall | Akaike's information criterion and Kullback-Leibler loss for histogram density estimation[END_REF] for the latter). Other versions of the CV principle have been applied to risk estimation and model selection, in particular the leave-p-out cross-validation which generalizes leave-one-out and the V-fold cross-validation. While leave-p-out CV is generally too computationally intensive as it needs evaluating n p models, Celisse and Robin provide in [START_REF] Celisse | Nonparametric density estimation by exact leave-p-out cross-validation[END_REF] an efficient closed-form formula for the squared L2-loss of histograms, including irregular ones. The same paper shows that the V-fold CV has a larger variance than the equivalent leave-p-out CV (i.e. when p = n V ) and is therefore less reliable. Additional results on leave-p-out CV are proved in [START_REF] Celisse | Optimal cross-validation in density estimation with the L 2 -loss[END_REF]. It is shown that the leave-one-out CV is optimal in terms of risk estimation (for the squared L2-loss) but not in terms of model selection, i.e. when the risk estimate is used to select θ in (1). In this latter task a leave-p-out CV is preferable and p should be adapted both to the size of the data and to the complexity of the model collection. Unfortunately, no rule for choosing an optimal value is currently known and experiments on simulated data confirm a strong link between p n and the quality of the associated histogram, see [START_REF] Celisse | Optimal cross-validation in density estimation with the L 2 -loss[END_REF]. Another way to address limitations of the simple techniques consists in using a penalized likelihood approach, a very common approach in model selection problems. The classical penalties are Akaike's information criterion AIC [START_REF] Akaike | Information theory and an extension of the maximum likelihood principle[END_REF] (used by [START_REF] Taylor | Akaike's information criterion and the histogram[END_REF] for histograms) and Schwarz' Bayesian information criterion BIC [START_REF] Schwarz | Estimating the Dimension of a Model[END_REF]. As shown in [START_REF] Rozenholc | Combining regular and irregular histograms by penalized likelihood[END_REF] AIC systematically underpenalises complex histograms (even regular ones) which makes it unsuitable for model section in this particular task (BIC does not suffer from this issue). Like simpler rules, AIC and BIC are based on asymptotic considerations only. This motivated Birgé and Rozenholc to use non asymptotic results on penalized maximum likelihood from [START_REF] Castellan | Modified Akaike's criterion for histogram density estimation[END_REF] to derive a modified AIC criterion for regular histograms [START_REF] Birge | How many bins should be put in a regular histogram[END_REF]). An even stronger penalty is needed for irregular histograms, as shown in [START_REF] Rozenholc | Combining regular and irregular histograms by penalized likelihood[END_REF].

Penalized estimators can also be constructed using Rissanen's approaches to minimal complexity. Hall and Hannan derive in [START_REF] Hall | On stochastic complexity and nonparametric density estimation[END_REF] two such criteria: one is based on Rissanen's stochastic complexity ideas [START_REF] Rissanen | Stochastic complexity and modeling[END_REF]) and the other one on Rissanen's minimum description length (MDL, Rissanen (1978)). An extension of the stochastic complexity based model to irregular histograms is proposed in [START_REF] Rissanen | Density estimation by stochastic complexity[END_REF]. As pointed out in the introduction, another MDL approach is proposed (Kontkanen and Myllymäki, 2007), using Normalized Maximum Likelihood (NML). We will describe in more detail this approach in Section 4.

Penalized likelihood approaches can generally be seen from a Bayesian point of view as instances of the maximum a posteriori (MAP) principle. Several authors have argued the advantages of using a MAP approach for histogram construction. For instance, the Bayesian block method [START_REF] Scargle | Studies in Astronomical Time Series Analysis. VI. Bayesian Block Representations[END_REF] addresses optimal segmentation with a MAP approach and proposes to apply its general framework to irregular histogram construction. A similar solution that uses a different likelihood and a different prior on parameters is proposed in [START_REF] Knuth | Optimal data-based binning for histograms and histogram-based probability density models[END_REF] for regular histograms.

Most of other methods introduced are for regular histograms (see for instance references in [START_REF] Birge | How many bins should be put in a regular histogram[END_REF]; [START_REF] Davies | A comparison of automatic histogram constructions[END_REF]; [START_REF] Rozenholc | Combining regular and irregular histograms by penalized likelihood[END_REF]). Among methods adapted to irregular histograms, the taut string approach [START_REF] Davies | Densities, spectral densities and modality[END_REF][START_REF] Davies | A comparison of automatic histogram constructions[END_REF] is interesting as it provides one of the few fully automated construction methods. It has been shown to give good results in extensive benchmarks, especially in terms of identifying the modes of a distribution [START_REF] Davies | A comparison of automatic histogram constructions[END_REF][START_REF] Rozenholc | Combining regular and irregular histograms by penalized likelihood[END_REF]. The approach is based on constructing a piecewise linear spline of minimal length, the taut string, that lies in a tube around the empirical cumulative distribution function. A histogram is constructed from the derivative of the taut string, see [START_REF] Davies | Densities, spectral densities and modality[END_REF]; [START_REF] Davies | A comparison of automatic histogram constructions[END_REF] for details. While the method is fully automated, it leverages the so-called κ-order Kuiper metrics whose parameters are set to some default values based on the data size via a heuristic. As stated by the authors, the performance of the taut string approach depends on those parameters. In addition, histograms built by this approach are not maximum likelihood histograms: densities associated to intervals are not directly derived from data frequencies in those intervals.

Another approach that does not fit directly into the penalized likelihood framework is the essential histogram proposed by [START_REF] Li | The essential histogram[END_REF]. The main idea of the method is to build a set of histograms that are close enough to the empirical distribution with respect to a collection of statistical tests. The essential histogram is the simplest of the collection. While appealing, the method uses a user defined confidence region specified by a unique parameter whose influence on the result cannot be discarded. The method is therefore not fully automated.

In the present section, we focused on the criteria and heuristics proposed to automatically select a possibly irregular histogram from the data. We have deliberately postponed the discussion of the algorithmic considerations to the next section.

Irregular histograms and dynamic programming

A histogram on an interval [a, b] is associated to a partition of [a, b] into K intervals {[a, c 1 ], ]c 1 , c 2 ], ..., ]c K-1 , b]}.
Optimising a quality criterion for a histogram implies to optimise over a subset of such partitions. The case of regular histograms is very simple, as there is only one partition associated to a number of sub-intervals K. On the contrary, irregular histograms are associated to an infinite set of partitions, which is not tractable in general.

Restricted irregular histograms

There are essentially two ways to address this problem. The data driven approach consists in restricting the subinterval endpoints c k to be observations, as in [START_REF] Davies | A comparison of automatic histogram constructions[END_REF]; [START_REF] Rozenholc | Combining regular and irregular histograms by penalized likelihood[END_REF]. Additional regularisation can be implemented by forbidding too short intervals but experiments in [START_REF] Rozenholc | Combining regular and irregular histograms by penalized likelihood[END_REF] show this has a limited and mostly negative impact on performance. The second approach consists in using a regular grid from which the endpoints are selected [START_REF] Celisse | Nonparametric density estimation by exact leave-p-out cross-validation[END_REF]Kontkanen and Myllymäki, 2007;[START_REF] Rissanen | Density estimation by stochastic complexity[END_REF]. This can be seen as enforcing an approximation accuracy on the observations, as detailed below. In both cases, the optimisation problem is now restricted to a finite set of partitions.

In this paper we use the fixed grid approach. We assume (for now) a given approximation accuracy . An observation x in the interval [x min , x max ] is approximated by the closest value x in X = {x min + t ; t = 0, ..., E -1}, E = 1 + L and L = x maxx min is the 'domain length' of the data. We restrict the possible accuracies such that E ∈ N. Following Kontkanen and Myllymäki (2007), we define the set of possible endpoints for sub-intervals as

C = x min - 2 + t ; t = 0, ..., E .
These endpoints define E elementary bins of length , which we call -bins (see figure 1). They are the building Notice that such a histogram is not bound to any given data set, despite the use of the term "counts" for the h k . Most of the methods studied in this paper, including our proposals, follow a maximum likelihood principle and therefore given K and C, the h k will be directly computed from the observations. Considering arbitrary histograms is nevertheless useful to include non maximum likelihood methods such as the taut string approach from [START_REF] Davies | Densities, spectral densities and modality[END_REF]; [START_REF] Davies | A comparison of automatic histogram constructions[END_REF], and to define easily prior distributions on histograms (see Section 5.1.1).

The piecewise constant density associated to M = (K, C, (h k ) 1≤k≤K ) is given by

f M (x) = 1 n K k=1 h k c k -c k-1 I ]c k-1 ,c k ] (x), (2) = 1 n K k=1 h k E k I ]c k-1 ,c k ] (x),
where E k is the number of -bins in interval [c k-1 , c k ], n = K k=1 h k and I A is the indicator function of the set A. Notice that the intervals are exclusive of their lower bound, even for the first interval ]c 0 , c 1 ]. This is by virtue of the definition of c 0 which is outside of the range of the data: no observation can be missed despite the exclusion of c 0 from the formal definition of the density. This enables us to keep a simple notation but it has no effect on the rest of the analysis.

Dynamic programming

While the two classical restrictions over the partitions we highlighted (data-driven grids and fixed grids) make the optimisation space finite, it remains large. There are E K-1 different partitions of E elementary bins in K intervals, and E may be big. However, most criteria used in the literature are additive and can therefore be optimised using the dynamic programming principle [START_REF] Bellman | On the approximation of curves by line segments using dynamic programming[END_REF]. This was proposed originally in [START_REF] Kanazawa | An optimal variable cell histogram[END_REF] for irregular histograms with endpoints constrained to observations. Since then, dynamic programming has been used systematically for all irregular histogram estimation techniques we are aware of.

Provided the chosen criterion can be computed efficiently, dynamic programming has a computational complexity of O(E 2 • K max ) when the histograms have at most K max bins (and thus a maximal complexity of O(E 3 )).

Data driven grid versus fixed grid

Both finite grids have advantages and drawbacks. Data driven grids do not introduce an additional precision parameter but this comes with a potential large computational cost. Dynamic programming has a computational cost of O(n 2 • K max ) for a data driven grid with n observations and a limit of K max bins in the histograms. If one sets K max = n to avoid introducing yet another parameter, data driven grids lead to a cubic cost O(n 3 ) and thus do not scale to large data sets. To circumvent this problem [START_REF] Rozenholc | Combining regular and irregular histograms by penalized likelihood[END_REF] use a greedy clustering of the data driven grid into max(n 1 3 , 100) bins prior to applying a dynamic programming approach. This blurs the distinction between data driven grids and fixed ones, and introduces a new parameter (the cut-off value of 100 used in [START_REF] Rozenholc | Combining regular and irregular histograms by penalized likelihood[END_REF]).

Another limitation of data driven grids is that they assume a perfect representation of real numbers. This is acceptable for small data sets but for large scale data sets, the probability of getting identical observations increases and the limits of real number representation start to play a role. This problem is discussed in [START_REF] Knuth | Optimal data-based binning for histograms and histogram-based probability density models[END_REF]; [START_REF] Knuth | Identifying excessively rounded or truncated data[END_REF] and in Section 6.2.

Finally in a Bayesian perspective, endpoints prior distributions have to be specified. With a data driven grid, continuous densities are needed while discrete distributions are needed when using a fixed grid. As a consequence, specifying priors that favor simple histograms is easier in the discrete case than in the continuous one.

Fixed grids bypass the above issues, at the cost of introducing a new accuracy parameter. We propose in Section 5.2 a way to automate the choice of this parameter. This mitigates the drawbacks of the fixed grid, while keeping all its advantages.

An NML criterion for histogram density estimation

Among minimum description length approaches, Kontkanen and Myllymäki's solution (Kontkanen and Myllymäki, 2007) derived using the Normalized Maximum Likelihood (NML) is the closest to our proposal. This approach consists in minimising a criterion c NML over all possible histograms defined on a fixed grid as in Section 3.1. The criterion is evaluated for a data set D = (x i ) 1≤i≤n of n observations.

To ease the comparison of this criterion to our own, we report the following simplified expression for c NML , using notations of Section 3.1 (the derivation is provided in Appendix A):

c NML (M|D) = log E K -1 + log R n M (3) + log n n h 1 h 1 • • • h K h K + K k=1 h k log E k ,
where log R n M is called the NML parametric complexity. Notice that the h k are computed from the data set D based on the histogram specification (see Section 5.1.2 for a discussion on this aspect).

The NML parametric complexity can be interpreted as the logarithm of the number of different distributions for a given model class [START_REF] Grunwald | The minimum description length principle[END_REF]. The term log E K-1 represents the code length of the choice of K -1 cut-points among the E possible candidates. Notice that this term is not directly derived from the NML distribution but was added by the authors as a way to index the different model choices, as recommended in [START_REF] Grunwald | The minimum description length principle[END_REF]. The last two terms represent the log-likelihood of the partition of n entries into K parts of size (h k ) 1≤k≤K and the log-likelihood of data distribution within each interval respectively.

The NML density has several important theoretical optimality properties: it is the minimax optimal universal model [START_REF] Shtarkov | Universal sequential coding of individual messages[END_REF]. Additionally, codes that are derived from it minimize the expected regret when choosing the worst-case generating density model for the data [START_REF] Rissanen | Strong optimality of the normalized ml models as universal codes and information in data[END_REF]. In Kontkanen and Myllymäki (2007), the authors show how to find the NML-optimal cut points via a dynamic programming in a polynomial time with respect to E, the total number of cut-points available for selection. Furthermore, their experiments showed that the minimization of the NML criterion produced good quality histograms of several Normal distributions of n = 10, 000 samples, all for a fixed approximation accuracy ( = 0.1).

We argue however that the NML approach has the following limitations.

NML computation: the NML approach results in a criterion whose exact computation remains costly. The cost of computing the NML parametric complexity was reduced from O(n 2 ) [START_REF] Kontkanen | Efficient computing of stochastic complexity[END_REF] to O(n + K) algorithm in Kontkanen and Myllymäki (2007).

Approximations of the NML parametric complexity have been introduced to reduce its computation time. For instance a O( √ dn + K) algorithm (where d is the precision in digits) is given in [START_REF] Mononen | Computing the multinomial stochastic complexity in sub-linear time[END_REF]. Szpankowski proposes in [START_REF] Szpankowski | On asymptotics of certain recurrences arising in universal coding[END_REF] an approximate computation in O(1) time but with an asymptotical error in O(n -3/2 ). The trade-off of using these faster methods is either a loss in precision or an error that is hard to evaluate non asymptotically w.r.t. n and K.

Notice that the main issue of the computational cost of the NML complexity is induced by the need to evaluate it at least O(n) times to build an optimal histogram, leading to minimal total cost of O(n 2 ) for an exact calculation.

Choice of : Although the authors argue that the sole user-defined parameter does not play a role in the model selection process, we show on the contrary in Section 6.2 that its impact cannot be overlooked. An automated choice is therefore needed.

Notice also that (Kontkanen and Myllymäki, 2007) uses a dynamic programming approach to optimise c NML which has a O(E 2 • K max ) or a O(E 3 ) computational cost, as already explained in Section 3.2. We improve over the NML approach on those three issues by introducing in Section 5 an enumerative criterion that avoids the NML computation cost, in Section 5.2 a granulated version of the our criterion that automates the choice of the precision parameter and finally in Section 7 greedy search heuristics that together reduce the global computational cost from O(E 3 ) to O(n log n).

G-Enum: a granulated enumerative criterion

This paper's first contribution is the introduction of a granulated enumerative two-part code for histogram density estimation, which we call G-Enum. It builds on a simpler enumerative approach, for which we introduce granularity later.

Enumerative criterion

Our enumerative criterion for histogram model selection exploits a maximum a posteriori Bayesian interpretation of the MDL principle which has the same optimality properties as the NML distribution [START_REF] Boullé | Revisiting enumerative two-part crude MDL for Bernoulli and multinomial distributions (Extended version)[END_REF]. In this framework, the best model is the one that maximises the probability of the model given the data, p(M|D) where D = (x i ) 1≤i≤n is a data set of n observations (we allow non unique values in D). This formulation is well known to be equivalent to a penalised likelihood approach. Indeed we have log p(M|D) = log p(M) + log p(D|M)log p(D).

As log p(D) does not depend on M, maximising p(M|D) is equivalent to maximising a compromise between a large likelihood log p(D|M) and a small prior probability log p(M) for a complex model and the reverse for a simple one. In practice, we minimise an enumerative criterion c Enum (M|D) =log p(M)log p(D|M).

(4)

Prior distribution

Using notations from Section 3.1, where K is the number of intervals, C is the set of endpoints of the histograms invertals and {h k } are the counts of values in each of these intervals, we factor p(M) as

p(M) = p(K) • p(C|K) • p((h k ) 1≤k≤K |K, C),
without loss of generality. Then, we use a hierarchical prior on the parameters and we assume conditional independence of C and {h k } 1≤k≤K given K. Thus, we have

p(M) = p(K) • p(C|K) • p((h k ) 1≤k≤K |K),
with the following components:

• Number of intervals

For the number of intervals K, we could use a uniform prior distribution leading to p(K) = 1 E , but Rissanen's universal prior for integers [START_REF] Rissanen | A universal prior for integers and estimation by minimum description length[END_REF] is more adapted as it favors small integers, i.e. simpler histograms. We have

p(K) = exp(-log * K), (5)
with log * K = log 2 (κ 0 ) + j>1 max log ( j) 2 (K), 0 , where κ 0 = p>1 2 -log * 2 (p) = 2.865 and log ( j) 2 (K) is the j-th composition of log 2 , i.e. log (1) 2 (K) = log 2 (K) and log ( j) 2 (K) = log 2 (log ( j-1)

(K)).

For the rest of the parameters, we use uniform priors.

• Interval composition

Specifying interval endpoints C is equivalent to specifying the number of elementary bins gathered into each interval. We use a uniform distribution over all the ordered partitions of the E elementary bins into K contiguous subsets of size

E 1 , E 2 , . . . E K such that E 1 + E 2 + • • • + E K = E.
We allow empty subsets, which leads to

p(C|K) = 1 E+K-1 K-1 . ( 6
)
Notice that we could forbid empty subsets (i.e. identical endpoints) which would lead to p

(C) = 1 ( E K-1 )
but this later solution is non monotone with respect to K and favors values around E 2 which is not a desirable property.

• Prior on (h k ) 1≤k≤K
For a given number of intervals K, the (h k ) 1≤k≤K specify the number of observations in each interval. A value for {h k } 1≤k≤K is therefore a vector of K non negative integers which sum to n. There are of n+K-1 K-1 such vectors and we use a uniform distribution over them, leading to

p((h k ) 1≤k≤K |K) = 1 n+K-1 K-1 (7)

Likelihood

The likelihood p(D|M) is obtained using a generative model for histograms specified by M which is also based on uniform distributions. Importantly, we do not assume the observations to be independent, contrarily to most of the histogram models.

Notice also that the {h k } 1≤k≤K in M are not computed from the data but parameters. This means that p(D|M) is non zero if and only if M and D are compatible.

Definition 1. A data set D = (x i ) 1≤i≤n and a histogram M = (K, C, (h k ) 1≤k≤K ) are compatible if and only if ∀k ∈ {1, . . . , K}, |{i ∈ {1, . . . , n} | x i ∈]c k-1 , c k ]}| = h k ,
where |S | denotes the cardinal of set S .

Notice that in practice the data "set" D can contain multiple times the same value and thus we take into account the index set when counting observations in this definition (see Section 6.2 for additional details on this aspect).

Given M, a data set D is generated by a hierarchical distribution based on uniform distributions. 1. Observations are generated by first choosing which interval is responsible for generating each individual observation. This corresponds to building a mapping I map from 1, . . . , n to 1, . . . , K specifying that x i will be generated in interval

I map (i). I map is compatible with M if ∀k ∈ {1, . . . , K}, i ∈ {1, . . . , n} | I map (i) = k = h k .
There are n! h 1 !...h K ! such mappings and we assume a uniform distribution on them, that is

p(I map |M) = h 1 !...h K ! n! (8)
2. Then given I map we generate observations assigned independently to distinct intervals, i.e.

p(D|I map , M) = K k=1 p((x i ) I map (i)=k |M). (9)
3. In a given interval, observations are generated by choosing their elementary bins independently and uniformly.

That is we assume

p((x i ) I map (i)=k |M) = 1 E k h k , (10)
as there are E k elementary bins in interval k and we have h k data points to generate. Notice than if h k = 0, then we can have E k = 0 as in this case ∀i, I map(i) k by compatibility.

Then when M and D are compatible, we have

p(D|M) = h 1 !...h K ! n! • K k=1 1 E k h k . ( 11
)
The simple enumerative criterion for histogram density estimation is thus given by

c Enum (M|D) = log * K + log E + K -1 K -1 + log n + K -1 K -1 + log n! h 1 !...h K ! + K k=1 h k log E k . ( 12
)
We use the convention that 0 × log 0 = 0 to avoid the case where h k = 0. Notice that equation ( 12) is valid only when M and D are compatible. When M and D are not compatible, p(D|M) = 0 and we set accordingly c Enum (M|D) = +∞.

As a consequence, histograms that are not compatible with the data are excluded from the solution space of our approach. This ensure that histograms obtained by minimizing the proposed criterion use a classical maximum likelihood estimation of the density in each of their intervals.

Comparison to the NML criterion

A side by side comparison of the Enum and NML criteria terms is presented in table 1. 2007)

log E K -1 log R n M + log n n h 1 h 1 • • • h K h K K k=1 h k log E k Enum log * K + log E + K -1 K -1 log n + K -1 K -1 + log n! h 1 !...h K ! K k=1 h k log E k
The terms that represent the indexing of -bins within each interval are exactly the same for both criteria. The differences lie in model indexing and the encoding of the multinomial distributions.

Model indexing terms. As in the NML criterion, we too have a term that serves as an index for model families. In the Enum criterion, model indexing is done through a log E+K-1 K-1 term. This term is monotone with K: it penalizes models with many intervals. The corresponding term in the NML criterion is not monotone and might thus favour choices with many cut-points (especially when K ∼ E).

Multinomial terms.

These terms encode the multinomial distributions of n observations into K intervals. Both versions are universal codes [START_REF] Grunwald | The minimum description length principle[END_REF] and have been proved to have the same optimality properties [START_REF] Boullé | Revisiting enumerative two-part crude MDL for Bernoulli and multinomial distributions (Extended version)[END_REF]. What sets them apart is that the enumerative version is easier to compute and to analyze.

The classical NML parametric complexity (log R n M ) is replaced by a simpler term in the Enum approach (log n+K-1 K-1 ). As shown in Kontkanen and Myllymäki (2007), the NML complexity can be computed exactly in O(n + K) time. In stark contrast, the equivalent term in our enumerative criterion can be computed in O(1) time.

As we can see, both criteria are very similar, though the Enum criterion is clearly simpler to compute and to analyse theoretically, using its closed form expression (see Section 6).

G-Enum criterion

The Enum criterion solves the computational issues induced by the NML parametric complexity without introducing approximations, but we still have to choose the approximation accuracy associated to the fixed grid. We introduce granularity to mitigate this issue.

Granularity and

The main issue with is that it plays two roles. It serves as a way to set a precision limit for the representation of real numbers, as discussed in Section 3.3. It also plays a crucial role in setting the resolution of the histograms themselves through E. The first role is more related to the data collection process and to the way the fixed grid is specified than to modelling, while the second appears directly in the model selection process via E.

We propose to explicitly separate those roles by introducing an intermediate parameter G, the granularity. For a given E, we assume that the numerical domain can be split into G bins (1 ≤ G ≤ E) of equal width. Each of these new elementary bins, that we call g-bins, is composed of g = E/G -bins. Each of the intervals of any constructed histogram has a length that is a multiple of these g-bins. In other words, each interval is no longer composed of a multiple of -bins but rather composed of G k g-bins. To better grasp this new setting, Figure 2 shows the case where three intervals were made from g-bins. In this illustration, we have E = 20 and G = 5, so that each g-bin is composed of 4 -bins.

By keeping , we still have a fixed grid and user defined resolution for both the data and the grid itself, while the introduction of G enables to chose the overall maximal complexity of the histogram.

Granulated criterion

Each granulated bin can either contain a whole number or a fraction of -bins. To simplify our reasoning and keep this paper self-contained, we assume that the number of elementary -bins E is a multiple of G, i.e. g = E/G ∈ N * .

Using this new parameter, the criterion changes both in model indexing terms and in the bin indexing terms:

• change in model indexing terms If we assume that G is distributed according to the universal prior for integers, the prior term for choosing the length of the intervals turns from log

E + K -1 K -1 into log * G + log G + K -1 K -1 .
• change in bin indexing terms Interval lengths are no longer described by E k -bins but rather by G k g-bins.

Given that each g-bin is actually an aggregation of g = E/G -bins, the likelihood term of the data distribution within an interval is now given by:

K k=1 h k log G k • E G = K k=1 h k log G k + K k=1 h k log E G = K k=1 h k log G k + n log E G (13)
This new criterion, which we call G-Enum is very similar to the original criterion introduced in this paper, as shown in table 2. The accuracy remains under the control of the user and can be fixed using standard considerations on real number computer representation (see Section 6.2), all the while G is optimised as part of model selection.

Enum log * K + log E + K -1 K -1 log n + K -1 K -1 + log n! h 1 !...h K ! K k=1 h k log E k G-Enum log * K + log * G + log G + K -1 K -1 log n + K -1 K -1 + log n! h 1 !...h K ! K k=1 h k log G k + n log E G

Theoretical properties

An analytical study of both the Enum and G-Enum criteria provides insights into how the enumerative MDL criterion works. Most of the propositions that follow cannot be extended to the NML criterion because of the non-monotone nature of some of its terms and because of the lack of closed-form expression for the parametric complexity. We focus solely on Enum and G-Enum histograms in what follows, as only these criteria are fit for analysis.

Compelling properties of enumerative histograms

We study the behaviour of the criterion in certain configurations in order to identify the most important factors that determine the overall shape of an optimal histogram (proofs of these properties are provided in Appendix B). Notice that although our analysis focuses on the encoding vision of the criterion, these properties can also be understood as statements of the most likely outcomes in terms of probability. This is due to the natural link between a coding length (the criterion value) and a posteriori probability at the basis of information theory: a more probable model will have a shorter coding length. Note also that in the properties, models are always assumed to be compatible with the data.

We show first an elementary property of optimal histograms.

Proposition 1. Let M = (K, (c k ) 0≤k≤K , (h k ) 1≤k≤K
) be an optimal histogram for the data set D. Then

∀k, 1 ≤ k ≤ K c k-1 < c k .
In other words, an optimal histogram cannot contain zero-length intervals.

Then we show that too complex histograms will not be selected by our criterion.

Proposition 2. Let D be a data set with n observations. Let us denote M K=n a histogram compatible with D such that there is one observation per interval and M K=1 a histogram compatible with D with only one interval. Then the coding length of M K=1 is shorter than the one of M K=n :

c Enum (M K=n |D) > c Enum (M K=1 |D).
This can be extended to more complex solutions with empty intervals.

Proposition 3. Let D be a data set with n observations. Let us denote M K>n a histogram compatible with D consisting of either singleton or empty intervals, one interval for each observation and empty intervals in-between, and let M K=1 be as in Proposition 2. Then the coding length of M K=1 is shorter than the one of M K>n :

c Enum (M K>n |D) > c Enum (M K=1 |D).
In addition, optimal histograms will never exhibit two consecutive empty intervals.

Proposition 4. For any data set D, the coding length of a histogram with two adjacent empty intervals is always longer than the coding length of a histogram with no consecutive empty intervals.

c Enum (M K,(h A ,h B )=0 |D) > c Enum (M K-1,h A∪B =0 |D)
Despite those general rules, local optimisation remains possible in the sense that empty intervals or intervals with only a single observation can be present in optimal histograms. Proposition 5. There exist data sets such that the optimal histogram has at least one interval which contains only a single observation. Proposition 6. There exist data sets such that the optimal histogram has at least one interval that does not contain any observation.

We can also characterise the structure of optimal histograms in different ways.

Proposition 7. An optimal histogram has at most 2n -2 intervals (K * ≤ 2n -2).

Proposition 8. In an optimal histogram, each interval endpoint is at most away from one of the values of the data set.

Propositions 4, 7 and 8 will serve as the basis for further improvements on the search of the best histogram. For instance, as Proposition 4 states that an optimal histogram will not have consecutive empty intervals, our search heuristic will systematically favour the merge of two consecutive empty bins. See Section 7 for additional details.

Role of the approximation accuracy

The authors of the NML criterion argue that the effect of the approximation accuracy can be ignored during the model selection process (Kontkanen and Myllymäki, 2007) and set to the accuracy at which the data have been recorded.

In their experiments, is fixed to a rather large value (0.1) which is inconsistent with current recording practices.

In this section, we conduct a theoretical investigation of the behaviour of the Enum criterion when → 0 (or equivalently, when E → ∞). Experimental results illustrate this behaviour in Section Appendix F.

Behaviour when → 0

To study the behaviour of our enumerative criterion as tends to 0, we first introduce the definition of singular intervals.

Definition 2. Let M = (K, C, (h k ) 1≤k≤K ) be a histogram compatible with a data set D and built on a grid. An interval

I k =]c k-1 , c k ] of M is singular if the following conditions are met 1. c k -c k-1 = ; 2. h k > 0;
3. all values of D that belong to I k are identical.

We denote by S (M, D) the number of data points from D that belong to singular intervals. In other words

S (M, D) = I k is singular h k .
Then we have the following theorem (see Appendix C for the proof). (

) 14
As c Enum (M |D) is dominated by the log 1 term when → 0, equation ( 14) shows that there is an inherent trade off between K and S (M , D). With a higher value of K, a histogram has more opportunities to have singular intervals and thus to have a larger value of S (M , D). If D contains multiple instances of the same value (i.e. there are i and j such that i j and x i = x j ), S (M , D) can grow fast enough to compensate the increase of K and it is therefore not possible to characterise further the asymptotic behaviour of c Enum (M |D).

However in the particular case where all values are distinct in D, we have the following corollary (see Appendix C for the proof).

Corollary 1. Let D be a data set with n distinct observations. Then for sufficiently small, the optimal histogram build on bins for the Enum criterion is the trivial one with a single interval

M = (1, {x min - 2 , x max + 2 }, n).
Enumerative histograms have a detrimental asymptotic behaviour: if we strive to be more precise in terms of cut point positions, the quality of the model will be poorer.

Illustration of the asymptotic behaviour

To complement the asymptotic analysis provided by Theorem 1 and its corollary, we study a simple example in the present section. Experiments on synthetic data are provided in Appendix F.

Let An optimal histogram should contain two intervals with a cut point close to α, provided that the mixture density θU[0, α] + (1 -θ)U[α, 1], is distinct enough from U[0, 1]. However Corollary 1 applies and if the precision is small enough, a single interval will be preferred. We can study on this simple example the relationship between , n and the mixture density to get a better understanding of the limitations proved by our results.

We study the optimal histogram with a unique interval, M 1, as well as the ideal histogram with two intervals M 2, with a cut point in α. Let ∆(n, , α, θ) be

∆(n, , α, θ) = c Enum (M 2, |D n (θ, α)) -c Enum (M 1, |D n (θ, α)).
We show in Appendix D that

∆(n, , α, θ) = log 1 + 1 -nD KL (θ α) + O(log n),
where D KL (θ α) is Kullback-Leibler divergence between the Bernoulli distribution with parameter θ (B(θ)) and the one with parameter α (B(α)). The histogram with two intervals is the optimal one if ∆(n, , α, θ) < 0.

For a fixed value of , the histogram with two intervals is preferred for larger values of n. The value of n needed to induce this preference is reduced by an increased KL divergence between B(θ) and B(α). An important aspect is that E = 1 influences the criterion in a logarithmic way while the influence of n is linear. This means that, in this example, large values of E can be compensated by relatively small values of n.

To illustrate further the behaviour, we computed exactly ∆(n, , α, θ) for θ = 1/2 and α = 1/10, and for different values of n and E. Table 3 summarizes those results by showing for several n the minimum value of E above which the single interval solution is preferred over the two interval ones. In this simple and extreme example, the transition appears only for very large values of E. However, as shown in Appendix F, in more realistic settings the effects of E manifests in a more reasonable range of values justifying the need for selecting it (or equivalently ) carefully.

Behaviour of the G-Enum criterion

While the decoupling between the precision of the grid and the resolution of the histograms provided by the G-Enum enables us to optimise the latter, it does not change the properties of the criterion for a fixed value of G.

This can be seen by interpreting a histogram M computed on G g-bins constructed from -bins as if it was constructed directly on g-bins. In other words, we can compute its quality according to the G-Enum criterion, taking into account the underlying -bins or according to the Enum criterion in which E is set to G. Then we have

c G-Enum (M, D) G,E = c Enum (M, D) G + log * G + n log E G ,
where the dependency to G and E has been made explicit using subscripts. This means that all propositions established for the Enum criterion still hold for the G-Enum criterion, for any fixed G parameter. In addition, the best histogram for G-Enum when G is fixed to a given value is exactly the best histogram for Enum with E = G. Thus, we can use any optimisation strategy designed to find the best Enum histogram to obtain the best G-Enum histogram for any fixed value of G.

Moreover, if we consider two different values of G, G 1 and G 2 , both divisors of E, we have

c G-Enum (M, D) G 1 ,E -c G-Enum (M, D) G 2 ,E = c Enum (M, D) G 1 -c Enum (M, D) G 2 + log * G 1 -log * G 2 + n log G 2 G 1 .
This shows that the optimal model for G-Enum depends only in a limited way on E. The choice of E simply constrains the possible values of G to its divisors, which in turn constrains the space of possible histograms. Importantly, in this space, the choice of the optimal histogram does not depend on E.

Notice also that while the G-Enum criterion does no longer suffer from the increase of E, it will behave similarly as G increases: for distinct data points, when G is larger than a data dependent bound G max (D) the optimal histogram for a fixed G consists in a single interval.

As a counter measure, we propose to set E to a very large value using the limits of numerical representation precision as a guideline. Given that integer values are encoded using four bytes and have their values between INT min ≈ -2.10 9 and INT max ≈ 2.10 9 , one can set E up to E max = 10 9 ≈ INT max /2.

Then we can select a set of values of G among the divisors of E, compute the optimal histogram for each of those values according to Enum and select the best overall one using G-Enum. Details about this procedure are given in Section 7.3.

Efficient search for MDL-optimal histograms

As recalled in Section 3.2, most irregular histogram construction methods leverage the dynamic programming principle to obtain an optimal histogram (with respect to the chosen criterion). This is in particular the case of the NML approach which has as a consequence a O(E 3 ) complexity (Kontkanen and Myllymäki, 2007).

We describe in this Section several techniques used to reduce this complexity. First, we introduce a greedy approach that can be applied to any additive criterion to obtain in O(E log E) time a sub-optimal histogram of good quality. Then, we show how to reduce this complexity to O(n log n) leveraging properties of enumerative histograms. Finally, we detail and discuss the particularities of the optimisation of granular models.

Greedy search

We propose a greedy search heuristic that combines a classic bottom-up approach and post-optimisations. It is based on a similar greedy algorithm to minimize additive criteria in the case of supervised discretisation [START_REF] Boullé | MODL: a Bayes optimal discretization method for continuous attributes[END_REF]. It applies to any additive criterion, defined as follows.

Definition 3. Let M = (K, C, (h k ) 1≤k≤K ) be a histogram compatible with a data set D of size n. The histogram is based on the intervals ]c k-1 , c k ] defined by C = (c k ) 0≤k≤K . A criterion that evaluate the quality of M with respect to D, q(M|D), is additive if it can be written q(M|D) = q 1 (n, K)

+ K k=1 q 2 (c k-1 , c k , h k ).
Algorithm 1 Optimisation with a greedy search Require: , D = (x i ) 1≤i≤n , additive criterion q Ensure: Histogram model M = (K, C, (h k ) 1≤k≤K )

Initialisation: 1: Sort(D) in O(n log n) 2: L bins = Create grid(D, )
with E bins 3: for consecutive bins A and B do in O(E)

4: ∆q ← q(M K-1,h A +h B |D) -q(M K,(h A ,h B ) |D) 5:
L merges ← ∆q 6: end for 7: Sort(L merges )

in O(E log E)

Optimisation: 8: ∆ ← Head(L merges ) 9: while ∆ < 0 do end for 16: end while

The main consequence of using an additive criterion is it locality. For instance when two adjacent intervals are merged, the value of q(M|D) for the resulting histogram can be computed from the value of the criterion of the original model considering only the two merged intervals. We exploit this property in the greedy search approach described by Algorithm 1.

The algorithm is rather simple: we start with the most refined histogram based of E -bins. A priority queue is used to store the effects on the quality criterion of merging two adjacent intervals (here -bins). Then we coarsen the histogram by implementing the best merge until the criterion cannot be improved (∆ < 0). The key point is that the merge qualities can be updated efficiently as a consequence of the additivity of the criterion.

Like most regularized histogram quality criteria, the Enum and NML criteria are additive: the overall computation time of the search for the best model can thus go from O(E 3 ) to O(E log E).

Although the heuristic has the advantage of being time and memory efficient, it may fall into a local optimum. We use two heuristics suggested in [START_REF] Boullé | MODL: a Bayes optimal discretization method for continuous attributes[END_REF] in to improve the quality of the final model. Rather than stopping the merge process as soon as ∆ < 0, we merge intervals up to the final histogram with a single interval. Then the best histogram is selected among all the histograms explored by this greedy merging approach. We post-optimise this histogram using local modifications of the intervals, as detailed in [START_REF] Boullé | MODL: a Bayes optimal discretization method for continuous attributes[END_REF]. This consists in choosing a set of simple operation on contiguous intervals (split, merge, merge and split, etc.) and in applying those operations in a greedy way until no improvement is possible. These operations as chosen in a way that do not modify the overall complexity of the algorithm.

Both heuristics are an important addition: extensive experiments reported in [START_REF] Boullé | MODL: a Bayes optimal discretization method for continuous attributes[END_REF] show that the greedy search alone produces an optimal solution only in roughly 50 % of the test cases, while the success rate increases to 95 % when heuristics are included.

Optimisation gains for enumerative histograms

The reduction from O(E 3 ) to O(E log E) is very important but our recommendation for setting E (see Section 6.2.2) would yield unacceptable running time even with the greedy algorithm. However theoretical properties of the Enum criterion can be leveraged to reduce further the complexity.

Proposition 8, in particular, shows that interval end points must be close to data points: regardless of the precision of the grid, i.e. of E, we need only to consider 2n -1 candidate splits which are the approximations of the data points. Thus, regardless of E, the exact optimisation of the Enum criterion can be done in O(n 3 ) while the greedy search has a complexity of O(n log n).

Optimisation of granular models

To optimise the G-Enum criterion, we use a large precision parameter E = 2 30 ≈ 10 9 , as large as possible w.r.t. the limits of numerical accuracy on computers (see section 6.3). We then exploit a loop on power of two granularities G = 2 i , 0 ≤ i ≤ 30 to retrieve the best model per granularity. As E is a multiple of G for each optimised granularity, the exact criterion of Table 2 holds. We exploit only the power of two granularities as a trade-off to both keep the computation time tractable and to explore a wide set of granularities. Assuming that n ≤ E, each step has a O(G i log G i ) time complexity for G i ≤ n and O(n log n) otherwise, since at most 2n intervals need to be considered for the optimisation of Enum histograms. Overall, the total number of operations t(n, E) is defined as follows.

t(n, E) = log 2 n i=1 2 i log 2 i + log 2 E i=1+log 2 n n log n, ≤ log 2 n i=1 2 i log n + (log 2 E -log 2 n)n log n, ≤ 2n log n + (log 2 E -log 2 n)n log n, ≤ (2 + log 2 E -log 2 n)n log n.
As E = 10 9 is a constant, the time complexity w.r.t. n is O(n log n). For a given G parameter, the optimisation of the G-Enum criterion can thus be performed in O(n log n).

Efficient search for other methods

As pointed out in Section 3.3, methods that use data driven grid are generally based on dynamic programming and have a O(n 3 ) complexity. Thus, they do not face the computational issues associated to a very fine fixed grid with a cost of O(E 3 ), contrarily to e.g. the NML criterion (Kontkanen and Myllymäki, 2007). An additional heuristic for fixed grid methods consists in restricting the cut points between intervals to approximations of the data points, i.e. to apply Proposition 8 even when its applicability has not been proved.

In the irregular histogram context, the applicability of dynamic programming is directly linked to the additivity of the criterion. Therefore, most of the methods reviewed earlier could benefit from the greedy search algorithm proposed above. Note however that its complexity is tied to the fact that the criterion can be updated in O(1). This is the case for instance for the penalized likelihood variants studied in [START_REF] Rozenholc | Combining regular and irregular histograms by penalized likelihood[END_REF] but not for the NML criterion (Kontkanen and Myllymäki, 2007). As recalled in Section 4, the parametric complexity term must be evaluated at least n times with an overall O(n 2 ) cost for an exact calculation, which would mitigate any advantage the heuristic could provide.

Experimental evaluation

The experimental evaluation is divided in three parts. First, we compare the performance of the three MDL criteria analysed in this paper. Secondly, we compare our MDL methods to the state-of-the-art algorithms identified as fully automated approaches to histogram building in the related work section. Both of these benchmarking comparisons are done on synthetic data sets of varying sizes, for which the performance of the different estimators can be objectively measured. Finally, we showcase the performance and practical relevance of our method for exploring real-world data of large scale.

A binary standalone implementation of the G-Enum based histogram construction is available here: http:// marc-boulle.fr/genum/. Hellinger distance.

The same can be said for HD values: NML, Enum and G-Enum histograms are very similar estimators. A steady decrease in HD values can be observed for the 5 variants, which means that all methods produce better estimations as the sample size increases. We highlight however that across all data sets, the G-Enum method has slightly better HD values, especially for large sample sizes (see table 5 for an overview).

Computation time. The main difference between the methods lies in the computation time. The results obtained across the different distributions showcase the advantage of using the greedy search heuristic to optimise MDL criteria. A particularly striking difference is observed for the Cauchy distribution (figure E.8), where NML histograms are computed up to 10 times faster with the heuristic than with the optimal algorithm for sample sizes from n = 100 to n = 10 4 . From n = 10 5 , the optimal algorithm took over an hour long. Additionally, we observed that for the uniform, triangle and triangle mixture distributions, the same amount of time is needed to compute NML histograms with either algorithm as the sample sizes increase. This shows that, for some data sets, the cost of computing the NML criterion itself can outweigh the benefits of using the search heuristic. For the Enum histograms the gap in computation time is not as easily closed. The Enum histograms computed with the optimal algorithm take slightly less time than the NML histograms for the Normal, Cauchy and Gaussian mixture distributions. For the uniform, triangle and triangle mixtures however, using the Enum even with the optimal algorithm can cut computation time by a factor of 100 as sample sizes increase. Enumerative histograms produced with the heuristic remain the fastest of all MDL-based methods throughout all experiments.

It is worth noting that, although the G-Enum method takes slightly longer, it remains consistently at a mid-point between the methods. Overall conclusion. NML and Enum histograms are interchangeable in terms of interval count and Hellinger distance, and this regardless of what algorithm is chosen to optimise the criteria. In terms of computation time however, there is a clear advantage in preferring the search heuristic and the simpler enumerative criteria. Finally, the use of G-Enum histograms increases slightly the computational time and improves also slightly the estimation quality in most cases in terms of Hellinger distance.

Comparison with other fully automated histogram methods

We compare in this section our MDL based methods with state-of-the-art fully automated methods. More precisely, the comparison includes:

• NML + heuristic, the NML criterion optimised with our search heuristic;

• G-Enum, our fully automated enumerative criterion;

• Taut string histograms [START_REF] Davies | Densities, spectral densities and modality[END_REF][START_REF] Davies | A comparison of automatic histogram constructions[END_REF], as implemented in the ftnonpar R package [START_REF] Davies | ftnonpar: Features and Strings for Nonparametric Regression[END_REF];

• RMG histograms [START_REF] Rozenholc | Combining regular and irregular histograms by penalized likelihood[END_REF], as implemented in the histogram R package [START_REF] Mildenberger | histogram: Construction of Regular and Irregular Histograms with Different Options for Automatic Choice of Bins[END_REF];

• Sturges rule histograms, as implemented in Python's numpy library [START_REF] Harris | Array programming with NumPy[END_REF];

• Freedman-Diaconis rule histograms [START_REF] Freedman | On the histogram as a density estimator:l2 theory[END_REF], as implemented in Python's numpy library;

• Bayesian blocks [START_REF] Scargle | Studies in Astronomical Time Series Analysis. VI. Bayesian Block Representations[END_REF], as implemented in Python's AstroPy library (Astropy Collaboration et al., 2022).

A very complete evaluation provided in [START_REF] Rozenholc | Combining regular and irregular histograms by penalized likelihood[END_REF] concluded that cross-validation based estimators were not among the best performing solutions and we excluded them from the comparison. In addition [START_REF] Rozenholc | Combining regular and irregular histograms by penalized likelihood[END_REF] showed that RMG histograms tend to provide the best overall performances.

For each studied distribution, we provide a visual example of histograms produced by the 7 methods, as well as the variations in interval counts, computation times and HD values as sample sizes increase. Detailed results are provided in Appendix G.

Given the extensive nature of these experiments, we selected a subset of the results to provide an overview for a n = 10 4 sample size in tables 6, 7 and 8.

Both the tables and the more detailed plots show that the other fully automated methods compared here work their best in the specific cases they were designed for. Although rarely the best for each distribution type, G-Enum histograms are consistently among the best estimators, and this without the high variability of the other methods. Focusing on irregular histograms, G-Enum is certainly among the most parsimonious in number of intervals. For exploratory analysis, this is an important quality because it makes the interpretation of the results easier and more reliable. G-Enum is also by far the fastest of irregular methods, making it suitable to large data sets.

Illustration on a large-scale real-world data set

We focus now on the performance and practical relevance of the G-Enum method for modelling an unknown distribution from a real-world large scale data set.

The Lunar Crater Database 1 contains 1.3 million entries on lunar impact craters larger than 1 to 2 km in diameter [START_REF] Robbins | A new global database of lunar impact craters >1-2 km: 1. crater locations and sizes, comparisons with published databases, and global analysis[END_REF]. Craters were manually identified and measured from images of NASA's Lunar Reconnaissance Orbiter (LRO), taken from 2011 until 2018. The run time and histogram sizes obtained on this data set are given in Table 9.

Although we ignore which law governs craters' diameter distribution, the granulated criterion produces a smoothlooking histogram with only K * = 75 intervals for the 1.3 million entries (Figure 3, with both axis in log scale). This fairly compact representation is easy to interpret. We notice for example, that the first 10 intervals are the most dense ones; they account for about 40% of all entries. Intervals become less dense for higher crater diameters. The last interval spans over about 1500 km of crater diameters and only accounts for 3 data entries.

Overall, the shape of our irregular histogram also shows that our approach can capture a power law decrease of the densities. This is in line with astrophysics literature: power or multiple power laws are often used to fit the crater size distribution [START_REF] Wang | Determining proportions of lunar crater populations by fitting crater size distribution[END_REF][START_REF] Minton | The equilibrium size-frequency distribution of small craters reveals the effects of distal ejecta on lunar landscape morphology[END_REF]. The same information would be hard to convey with an 1 https://astrogeology.usgs.gov/search/map/Moon/Research/Craters/lunar_crater_database_robbins_2018 The density estimation experiment on this data set also shows the scalability of our method: the granular histogram was computed in about 1.3 seconds. The taut-string method and the RMG approach remain usable with respectively about 54 seconds and 37 seconds, but both NML histogram (≈ 17 minutes of calculation) and the Bayesian Blocks method (almost 1 hour and 15 minutes of run time) use an unreasonable quantity of computational resources.

Figure 4 displays the histograms obtained by those methods while Table 10 gives the Hellinger distances between the densities estimated by the different methods. The NML histogram is the most unique one. The log log display used on the figures hides to some extent the main source of disagreement between the histograms which is the shorter first interval used by the NML histogram compared to the others. All other histograms are fairly close to one another. This emphasizes the non-parsimonious nature of the Taut string method which uses 478 intervals to produce an histogram that is fairly close to the one obtained by the Bayesian Blocks method with only 86 intervals. In addition, while this does not play a major role in the Hellinger distance because of its quite low support, the estimation of the tail of the distribution by the Taut string histogram seems to be more crude than other solutions. Overall, this experiment confirms the computational efficiency of G-Enum and its ability to produce compact representations of complex and unknown distributions. Importantly, the resulting models are better or comparable to the ones obtained by less efficient solutions. G-Enum histograms provide interesting insights on the previously unknown distribution laws of large data sets without needing much computation time.

Conclusion

We presented a simple yet robust and very efficient enumerative criterion for histogram model selection that produces indistinguishable results to those obtained with much more compute-costly NML approach presented in a previous work.

By pairing our criterion with a search heuristic rather than the optimal but costly original optimisation algorithm, we achieve substantial gains in computation time. By introducing granularity to alleviate our dependency to the sole user parameter of this problem, the approximation accuracy , we achieve substantial gains in robustness.

With our theoretical and experimental evaluation of these criteria we show that our granulated MDL criterion fills a gap in the current histogram model selection landscape : it's a resilient, efficient and fully automated approach to histogram density estimation that can scale to explore known or unknown distribution laws in large data sets. and

h k =        h k when k < l h k+1 when k ≥ l.
Obviously, M defines the same density as M. In addition, as h l = 0

K k=1 h k log E k = K k=1,k l h k log E k , = K k=1,k l h k log c k -c k-1 , = K-1 k=1 h k log c k -c k-1 , = K-1 k=1 h k log E k ., and
log n! K k=1 h k ! = log n! K k=1,k l h k ! , = log n! K-1 k=1 h k ! . Therefore ∆ = c Enum (M|D) -c Enum (M |D) is given by ∆ = log * K -log * (K -1) + log E + K -1 K -1 -log E + K -2 K -2 + log n + K -1 K -1 -log n + K -2 K -2 ,
which shows that ∆ > 0 and thus that M cannot be optimal.

Appendix B.2. Proof of proposition 2 Proposition 2. Let D be a data set with n observations. Let us denote M K=n a histogram compatible with D such that there is one observation per interval and M K=1 a histogram compatible with D with only one interval. Then the coding length of M K=1 is shorter than the one of M K=n :

c Enum (M K=n |D) > c Enum (M K=1 |D).
Proof. If M K=n is compatible with D and there is one observation per interval, then ∀k, h k = 1 and obviously, E > n. Also ∀k, E k > 0.

Let us define δc Enum (n, 1) = c Enum (M K=n )c Enum (M K=1 ).

Then we have

δc Enum (n, 1) = log E + n -1 n -1 + log 2n -1 n -1 + log n! + n k=1 log E k + log * n -log * 1 -n log E. Since log E+n-1 n-1 = n-1 k=1 log(E + k) -log(n -1)!, we get δc Enum (n, 1) = n-1 k=1 log(E + k) + log 2n -1 n -1 + log n + n k=1 log E k + log * n -log * 1 -n log E.
As function f (x) = log(1 + x) is strictly concave and f (0) = 0, it is sub-additive on [0, +∞[. We obtain

n k=1 log E k = n k=1 f (E k -1) ≥ f        n k=1 (E k -1)        ≥ log(E -n + 1).
Back to δc Enum (n, 1), we get

δc Enum (n, 1) ≥ n-1 k=1 log E(1 + k E ) + log 2n -1 n -1 log n + log E(1 - n -1 E ) + log * n -log * 1 -n log E ≥ log * n -log * 1 + log n + log 2n -1 n -1 + n-1 k=1 log(1 + k E ) + log(1 - n -1 E ).
As n ≤ E, we get

log(1 - n -1 E ) + log n = log(n - n E (n -1)) ≥ log(n -(n -1)) ≥ 0.
Finally, δc Enum (n, 1) > 0, which proves the claim that the histogram with one single interval is always more probable than the one with n singleton intervals.

Appendix B.3. Proof of proposition 3 Proposition 3. Let D be a data set with n observations. Let us denote M K>n a histogram compatible with D consisting of either singleton or empty intervals, one interval for each observation and empty intervals in-between, and let M K=1 be as in Proposition 2. Then the coding length of M K=1 is shorter than the one of M K>n :

c Enum (M K>n |D) > c Enum (M K=1 |D).
Proof. Let us consider a histogram consists with K > n intervals composed of E k -bins, with either h k = 1 for the singleton intervals and h k = 0 for the empty intervals.

δc Enum (K, 1) = log E + K -1 K -1 + log n + K -1 K -1 + log * K -log * 1 + log n! + n k=1 log E k -n log E. We have log E + K -1 K -1 = K-1 k=1 log(E + k) log k = n k=1 log(E + k) log k + K-1 k=n+1 log(E + k) log k Given that log(E + k) > log E, we have n k=1 log(E+k) log k > n log E -log n!. Likewise, we know that log(E + k) > log k, so the sum K-1 k=n+1 log(E+k) log k is greater than K-1 k=n+1 1 = K -n -1.
Given these two lower bounds, we can write that

log E + K -1 K -1 > n log E -log n! + (K -n -1)
We then obtain δc Enum (K, 1) > (Kn -1)

+ log n + K -1 K -1 + log * K -log * 1 + n k=1 log E k .
Finally, δc Enum (K, 1) > 0, which proves the claim that the histogram with one single interval is more probable that any histogram where entries are in singleton intervals.

Note 1. Similar results to Proposition 2 and Proposition 3 could not be obtained for the NML criterion given by equation (3) for two reasons. First, the analysis of the criterion for histograms with K singleton intervals (n ≤ K ≤ E) is complex, since the prior term log E K-1 is not monotonous and decreases for K > E/2. Second, the parametric complexity log R n M K term has no closed-form formula and its best known approximations [START_REF] Kontkanen | Computationally efficient methods for MDL-optimal density estimation and data clustering[END_REF][START_REF] Szpankowski | On asymptotics of certain recurrences arising in universal coding[END_REF]) are given for fixed K as n goes to infinity. These approximations cannot be used when n ≤ K ≤ E.

Appendix B.4. Proof of proposition 4

We begin by evaluating the expression of the cost variation ∆c after the merge of two intervals. Let A and B be two adjacent intervals composed of E A , E B -length elementary bins and of data counts h A , h B in a K-bin histogram. They are merged into a single interval composed of E A∪B = E A + E B -length elementary bins with h A∪B = h A + h B elements, creating a histogram with K -1 bins (we consider only histograms compatible with the data).

For this fusion, the variation for the Enum criterion, ∆c, is given by

∆c Enum =c Enum (M K-1 ) -c Enum (M K ) (B.1) = log K -1 E + K + log K -1 n + K -1 + log * (K -1) -log * K + log h A !h B ! (h A + h B )! + log (E A + E B ) (h A +h B ) E A h A • E B h B
The fusion of two intervals is advantageous if the code length for the histogram with K -1 intervals is smaller than the code length for the histogram with K intervals, that is if ∆c Enum < 0. Proposition 4 corresponds to an interesting particular case.

c Enum (M K (x)|D) = log E + K -1 K -1 + log n + K -1 K -1 + log * K + log n! h 1 ! . . . h i !h j ! . . . h K ! + k<i h k log E k + h i log(E i + x) + h j log( E -(E i + x)) + k> j h k log E k
Seeing as we focus on the frontier between just two intervals, the model index term does not change. Since the data counts of each interval are preserved regardless of the position of the interval's endpoints, the corresponding part of the likelihood term will not change. In this setting, the sole term responsible for the cost variation between models is:

f h i ,h j (x) = h i log(E i + x) + h j log( E -(E i + x))
From Proposition 4, an optimal histogram cannot contain two consecutive empty intervals, therefore at least one of the two intervals is non empty. We study the three possible cases.

1. Let us first consider the case where both intervals are non empty (h i > 0 and h j > 0). As f h i ,h j is a concave function on the closed interval [0, E -E i -E j ], its minimum values are only reached at its extremities, that is for x = 0 and x = E -E i -E j .

We have thus shown that the best cut-point between the two intervals is placed at at most , either from last value of the first interval or from the first value of the second interval. 2. Let us now consider the case where the first interval is non empty, that is h i > 0 and h j = 0. Then f h i ,h j (x) = h i log(E i + x). In this case, the minimum value of f h i ,h j (x) is obtained for x = 0, that is with a cut-point placed at at most from last value of the first interval. 3. Similarly, if the second interval is non empty, the cut-point is placed at at most from the first value of the second interval. Let us notice that in both cases, the empty interval is of maximal length.

Overall, the endpoints of the optimal histogram are always at at most from one value of the data set.

Note 2. Intuitively, the choice of which one of these two placements is best will also depend on the data count of each interval. As we have highlighted before, a part of the likelihood term will favour shorter dense intervals over large ones. See the proof of proposition 5 for an illustration in a simple case. Finally, as in the proof of proposition 5, the optimal histogram with K = 3 intervals consists in two non-empty intervals, each with E A = 1, E B = 1 -bins, surrounding an empty interval composed of E -2 -bins. Therefore

c Enum (M K=3 ) = log * 3 + log (E + 1)(E + 2) 2 + log (n + 1)(n + 2) 2 + log n! h A !h B ! = n log 2 + 2 log E + 3 2 log n + O(1).
Therefore, even for rather small n and E, the model cost is minimal for K = 3, for the histogram where the second interval is empty and has (E -2) • . (C.1)

Proof. Let D be a fixed data set with n observations. For any such that L is an integer, M = (K, (c k ) 0≤k≤K , (h k ) 1≤k≤K )

be a optimal histogram constructed on the E = L + 1 regular grid of -bins. The Enum criterion for M is given by:

c Enum (M |D) = log * K + log n + K -1 K -1 + log E + K -1 K -1 + log n! h 1 !...h K ! + K k=1 h k log E k, ,
where the E k, are the sizes of the intervals expressed in terms of -bins. More precisely, let c 0, , c 1, , . . . , c K, be the end points of the K intervals of M . By definition of the grid, we have

L k, = c k, -c k-1, = E k, × , 1 ≤ k ≤ K and L + = c K , -c 0, = E × . Using log E + K -1 K -1 = K-1 k=1 log(E + k) -log(K -1)!,
and the notations above, we have

c Enum (M |D) = log * K + log n + K -1 K -1 + log n! h 1 !...h K ! -log(K -1)! + K-1 k=1 log L + k + 1 + K k=1 h k log L k, .
Obviously, the main influence of on the criterion is through the last two terms. Let us first notice that ∀k,

1 ≤ k ≤ K -1 log L + k + 1 = log 1 + log(L + (k + 1) ).
and thus By proposition 8 we know that ∀0

K-1 k=1 log L + k + 1 = (K -1) log 1 + K k=2 log(L + k ). (C.2) To analyse the last term K k=1 h k log L k, ,
≤ k ≤ K, min x∈D |c k, -x| ≤ , as M is optimal. This shows that |L k, -D k, | ≤ 2 . Using D k, , we have K k=1 h k log L k, = (C.3) D k, =0 h k log L k, + D k, >0 h k log 1 + log L k, D k, + log D k, .
Using equations (C.2) and (C.3), we can decompose the Enum criterion c Enum (M |D) into the sum of the following three terms:

c 1 (M |D) = log * K + log n + K -1 K -1 + log n! h 1 !...h K ! -log(K -1)! + D k, =0 h k log L k, + D k, >0 h k log D k, , c 2 (M |D) = K k=2 log(L + k ) + D k, >0 log L k, D k, . c 3 (M |D) = (K -1) log 1 + D k, >0 h k log 1 ,
As M is optimal, K ≤ 2n -2 (proposition 7) and as a consequence c 1 (M |D) is upper and lower bounded by constants that depend only on the data set. This is obviously the case of the terms log K, log n+K-1 K-1 andlog(K -1)!. The term log n! h 1 !...h K ! can take only a finite number of positive values for a fixed n and any K ≤ 2n -2 and is therefore bounded by the largest of those values. Moreover, when D k, = 0, ≤ L k, ≤ 2 and thus 0

≤ D k, =0 h k log L k, ≤ n log 2.
It is also obvious that D k, ≤ L and thus

D k, >0 h k log D k, ≤ n log L.
To lower bound this quantity, we introduce

D min = min{|x i -x j | | x i ∈ D, x j ∈ D, x i x j }, (C.4)
which depends only on the data set D and is strictly positive (by hypothesis on D). By definition for all k such that D k, > 0, D k, ≥ D min , and therefore

D k, >0 h k log D k, ≥ min(log D min , 0). Thus there exists C 1 (D) ≥ 0 such that for all > 0, -C 1 (D) ≤ c 1 (M |D) ≤ C 1 (D).
The second term c 2 (M |D) can also be upper and lower bounded with a minimal condition on . Indeed we have

(K -1) log L ≤ K k=2 log(L + k ) ≤ (K -1) log(L + K ).
Assuming that ≤ L and with K ≤ 2n -2, we have

(2n -3) log L ≤ K k=2 log(L + k ) ≤ (2n -3) log L + log(2n -1) . 36
We also have for D k, > 0

-2 ≤ L k, -D k, ≤ 2 -2 D k, ≤ L k, D k, -1 ≤ 2 D k, 1 -2 D k, ≤ L k, D k, ≤ 1 + 2 D k, 1 -2 D min ≤ L k, D k, ≤ 1 + 2 D min
Thus if we assume ≤ D min 4 < L, we have

-log 2 ≤ log L k, D k, ≤ log 3 2 .
Thus if ≤ D min 4 , we have

-K log 2 ≤ D k, >0 log L k, D k, ≤ K log 3 2 -(2n -2) log 2 ≤ D k, >0 log L k, D k, ≤ (2n -2) log 3 2
Thus there exists C 2 (D) ≥ 0 such that for all K ≤ 2|D| -2 and for all 0

< ≤ D min 4 , -C 2 (D) ≤ c 2 (M |D) ≤ C 2 (D)
. The last term c 3 is the only one that depends on in a non bounded way. It can be interpreted in a more direct way using

D k, >0 h k = n - D k, =0 h k ,
and characterising intervals such that h k > 0 and D k, = 0. Let I k be such an interval and let x be a point from D in I k . When D k, = 0, L k, = or L k, = 2 . Combined with the assumption that ≤ D min 4 , all data points in I k take the same value x.

Using proposition 8, we can then rule out the case L k, = 2 . Indeed we have

I k =]c k-1 , c k ].
According to proposition 8 both c k-1 and c k must be at most distant of from a data point. Considering D k, = 0 this must be the same value x . If c kc k-1 = 2 , then the only possibility is that c k = x + and c k-1 = x -, and thus x = x, the only value taken by data points in I k . However based on Remark 3 above, we should then have c k = x. Thus, L k, = and I k is a singular interval as per definition 2. Therefore D k, =0 h k = S (M * , D), which concludes the proof. is 1. Therefore when < e -2C(D) , the sign of

c Enum (M K, |D) -c Enum (M K , |D),
is the one of K -K + S (M * K , , D) -S (M * K, , D) as long as this quantity is not zero. In other words, histograms can be compared based on this difference as long as is small enough. Our goal is to show that M 1, is optimal. This histogram is preferred over M * K, when

1 -K + S (M * K, , D) -S (M * 1, , D) < 0, S (M * K, , D) < K -1,
as M 1, does not contain singular intervals.

Let us now focus on the assumption that the n values in the data set are distinct. This implies that a singular interval I k is associated to h k = 1 and thus S (M * K, , D) is exactly the number of singular intervals. This number is controlled by the fact that when is small enough, a histogram cannot have adjacent singular intervals. Indeed singular intervals have a maximal width of . Let I k and I k+1 be two such adjacent intervals, with the associated values x i ∈ I k and x j ∈ I j . Then D min ≤ |x jx i | ≤ 2 , with D min has defined in equation (C.4). The lower bound is induced by the fact that I k and I k+1 contain only a single value each (by definition) and are adjacent: there is no value from the data set in ]x i , x j [. The inequalities can be fulfilled simultaneously only if ≥ D min 2 and thus when < D min 2 , singular intervals cannot be adjacent. We assume in the rest of the proof that < min D min 2 , e -2C(D) . Then This shows that M * 1, is always preferred to by the Enum criterion to histograms with S (M * K, , D) ≥ 3. We need now to discuss the remaining cases, i.e., situations when S (M * K, , D) ≤ 2. If K ≥ S (M * K, , D) + 2, M * 1, is preferred and thus we have two interesting particular cases to handle:

1. K = 2 and S (M * 2, , D) = 1: this is a particular case with a single singular interval completed by a single large non singular one; 2. K = 3 and S (M * 3, , D) = 2: this is again a particular case where we have two singular intervals separated by a single large non singular interval.

In both cases, the dominating term of the Enum criterion is n log 1 , exactly as for M * 1, and we need to compare with more precision the values of the criterion to choose the optimal histogram.

For a single interval, we have In stark contrast, G-Enum remains at 3 intervals throughout. The same consistency can be observed in regards to computation times and Hellinger distance. While the NML and Enum methods take more time and produce lower quality results as the approximation accuracy increases, the G-Enum steadily produces quality results in less time.

Note also that the best value for is different for both sample sizes. For a distribution of n = 100 samples, having = 0.1 guarantees the lowest Hellinger distance for the NML and Enum methods. For n = 1000 samples, it is preferable to set = 0.5. When we do not know anything about the nature of the data, there is no guarantee we will make the right choice for NML and Enum. On the other hand, our granularised approach automatically selects * ≈ 1.324 and * ≈ 0.178 for samples sizes n = 100 and n = 1000 respectively. While these choices certainly do not achieve the lowest HD technically possible, they have the merit of not being too far off. G-Enum histograms seem to select fair enough accuracy values for a fully automated estimation.

These sensitivity experiments have shown that the value of cannot be overlooked as easily as thought: for some distributions and middle-range sample sizes, the approximation accuracy can play an important role in computation time and estimation quality. In an exploratory analysis context, where little is known about the data, a truly fully automated approach such as G-Enum is preferable than the other MDL methods.

Figure 1 :

 1 Figure 1: A possible choice of intervals, with their respective data-counts (h k ) 1≤k≤K and lengths expressed in number of -bins (E k ) 1≤k≤K . blocks of histogram intervals: possible grouping of consecutive -bins into K intervals, with K ranging from 1 to E, defines a histogram. A histogram M is fully specified by M = (K, C, (h k ) 1≤k≤K ) where K is the number of intervals, C = (c k ) 0≤k≤K ⊂ C are the endpoints and h k are the counts in each bin (we have c 0 = x min -/2 and c K = x max + /2).Notice that such a histogram is not bound to any given data set, despite the use of the term "counts" for the h k . Most of the methods studied in this paper, including our proposals, follow a maximum likelihood principle and therefore given K and C, the h k will be directly computed from the observations. Considering arbitrary histograms is nevertheless useful to include non maximum likelihood methods such as the taut string approach from[START_REF] Davies | Densities, spectral densities and modality[END_REF];[START_REF] Davies | A comparison of automatic histogram constructions[END_REF], and to define easily prior distributions on histograms (see Section 5.1.1).

Figure 2 :

 2 Figure 2: Replacing -bins by g-bins in the model building

Theorem 1 .

 1 Let D be a data set with n observations. There exists two positive values C(D) and E(D) that depends only on D such that for all ≤ E(D) for any optimal histogram M have c Enum (M |D) -K -1 + n -S (M , D) log 1 ≤ C(D).

  U[a, b] denote the uniform distribution on the interval [a, b]. Let α and θ be two real numbers with α ∈]0, 1 2 ] and θ ∈] 1 2 , 1]. Let us consider a data set D n (θ, α) consisting of nθ independent observations generated by U[0, α] and n(1 -θ) independent observations generated by U[α, 1].

Figure 4 :

 4 Figure 4: Different histograms obtained for the Lunar data set

Method

Figure B. 5 :

 5 Figure B.5: Illustration and notations of the optimal cut-point between two values problem

Figure B. 6 :

 6 Figure B.6: Graphical representation of the different cases involved in proposition 5.

  Appendix C. Asymptotic behaviour of the MODL criterion when → 0 Appendix C.1. proof of Theorem 1 Theorem 1. Let D be a data set with n observations. There exists two positive values C(D) and E(D) that depends only on D such that for all ≤ E(D) for any optimal histogram M have c Enum (M |D) -K -1 + n -S (M , D) log 1 ≤ C(D).

  let us define D k, by D k, = arg min x∈D |c k, -x|arg min x∈D |c k-1, -x|, that is the distance between the data points that are the closest to the boundaries of interval ]c k-1, , c k, ]. By definition, we have |L k, -D k, | ≤ min x∈D |c k, -x| + min x∈D |c k-1, -x|.

K

  ≥ 2S (M * K, , D) -1,as we need "in between" non singular intervals to separate singular ones. Therefore we haveK -S (M * K, , D) ≥ S (M * K, , D) -1, and if S (M * K, , D) ≥ 3, S (M * K, , D) ≤ K -2.

c

  Enum (M 1, |D) = log * 1 + n log L + 1 . (C.5)We use the fact that if λ > 0, ρ > 0 and κ do not depend on , for large enough such that ρ + κ > 0,λ log ρ + κ = λ log ρ + λ log 1 + o( ). (C.6) Thus we have c Enum (M 1, |D) = n log 1 + log * 1 + n log L + o( ).

Figure E. 8 :

 8 Figure E.8: Comparison between MDL methods over a Cauchy distribution of different sample size

Figure F. 14 :

 14 Figure F.14: Comparison of MDL methods over a Normal distribution of size n = 1000 for different values of

Table 1 :

 1 Term comparison of the NML and enumerative criteria

	Crit.	Indexing terms	Multinomial terms	Bin index terms
	NML
	Kontkanen and
	Myllymäki
	(

Table 2 :

 2 Term comparison of the Enum and G-Enum criteria

	Crit.	Indexing terms	Multinomial terms	Bin indexing terms

Table 3 :

 3 Transition from two intervals to one single interval.

	n	E
	10	30
	12	80
	16	530
	20	3700
	30 5.05 × 10 5
	40 7.28 × 10 7
	50 1.08 × 10 10

Table 4 :

 4 All the densities used for benchmarking

Table 5 :

 5 Hellinger distance values over different distributions of sample size n = 10 4 , for MDL methods

	Distribution	NML + optimal	NML + heuristic	Enum + optimal	Enum + heuristic	G-Enum
	Normal	0.046 ± 0.002	0.047 ± 0.002	0.046 ± 0.002	0.047 ± 0.002	0.045 ± 8 • 10 -4
	Cauchy	0.074 ± 0.003	0.073 ± 0.003	0.075 ± 0.004	0.075 ± 0.004	0.060 ± 0.004
	Uniform	0.051 ± 0.005	0.051 ± 0.005	0.052 ± 0.006	0.052 ± 0.006	0.024 ± 0.001
	Triangle	0.038 ± 0.002	0.038 ± 0.002	0.038 ± 0.002	0.039 ± 0.002	0.039 ± 0.002
	Triangle mixture	0.038 ± 0.003	0.038 ± 0.003	0.039 ± 0.003	0.038 ± 0.003	0.039 ± 0.002
	Gaussian mixture	0.059 ± 0.002	0.059 ± 0.002	0.058 ± 0.002	0.058 ± 0.002	0.057 ± 0.002

Table 6 :

 6 Hellinger distance values over different distributions of sample size n = 10 4

	Distribution	G-Enum	NML	BB	TS	RMG	FD	Sturges
			Kontkanen and	Scargle et al.	Davies et al.	Rozenholc et al.	Freedman
			Myllymäki	(2013)	(2009)	(2010)	and Diaconis
			(2007)				(1981)
	Normal	0.045 ± 6.10 -4	0.046 ± 0.002	0.047 ± 0.002	0.040 ± 0.002	0.034 ± 0.002	0.033 ± 0.002	0.055 ± 0.002
	Cauchy	0.061 ± 0.004	0.074 ± 0.003	0.064 ± 0.002	0.045 ± 0.005	0.064 ± 0.001	0.138 ± 0.002	0.862 ± 0.036
	Uniform	0.024 ± 0.001	0.050 ± 0.005	0.025 ± 0.004	0.031 ± 0.015	0.029 ± 0.011	0.082 ± 0.011	0.028 ± 0.002
	Triangle	0.039 ± 0.002	0.038 ± 0.0025	0.039 ± 0.001	0.084 ± 0.024	0.084 ± 0.029	0.032 ± 0.002	0.049 ± 9.10 -4
	Triangle mixture	0.037 ± 0.002	0.038 ± 0.003	0.040 ± 0.003	0.078 ± 0.029	0.069 ± 0.026	0.032 ± 0.002	0.043 ± 4.10 -4
	Gaussian mixture	0.057 ± 0.002	0.059 ± 0.002	0.060 ± 0.002	0.040 ± 0.001	0.052 ± 0.002	0.060 ± 0.001	0.142 ± 0.013

Table 7 :

 7 Computation times (in seconds) over different distributions of sample size n = 10 4

	Distribution	G-Enum	NML	BB	TS	RMG	FD	Sturges
			Kontkanen and	Scargle et al.	Davies et al.	Rozenholc et al.	Freedman
			Myllymäki	(2013)	(2009)	(2010)	and Diaconis
			(2007)				(1981)
	Normal	0.014 ± 0.003	2.724 ± 0.283	5.785 ± 0.479	0.014 ± 0.002	1.239 ± 0.085	0.002 ± 2.10 -4	6.10 -4 ± 2.10 -4
	Cauchy	0.028 ± 0.006	121.60 ± 4.99	3.250 ± 0.112	0.116 ± 0.205	0.906 ± 0.107	0.009 ± 0.014	0.001 ± 0.003
	Uniform	0.015 ± 0.002	0.168 ± 0.012	5.989 ± 0.167	0.011 ± 0.002	1.387 ± 0.139	0.002 ± 3.10 -4	6.10 -4 ± 2.10 -4
	Triangle	0.014 ± 0.005	0.169 ± 0.005	5.962 ± 0.113	0.015 ± 0.002	1.291 ± 0.091	0.002 ± 2.10 -4	6.10 -4 ± 2.10 -4
	Triangle mixture	0.012 ± 0.006	0.103 ± 0.027	3.004 ± 0.245	0.013 ± 0.006	0.954 ± 0.138	0.002 ± 0.005	0.0 ± 0.0
	Gaussian mixture	0.017 ± 0.002	1.91 ± 0.085	4.165 ± 0.369	0.048 ± 0.005	1.056 ± 0.077	0.006 ± 0.008	0.002 ± 0.005

Table 8 :

 8 Optimal bin counts over different distributions of sample size n = 10 4

	Distribution	G-Enum	NML	BB	TS	RMG	FD	Sturges
			Kontkanen and	Scargle et al.	Davies et al.	Rozenholc et al.	Freedman
			Myllymäki	(2013)	(2009)	(2010)	and Diaconis
			(2007)				(1981)
	Normal	16.30 ± 0.46	15.60 ± 1.02	15.90 ± 1.04	72.50 ± 4.41	39.70 ± 7.34	62.40 ± 2.29	15.0 ± 0.0
	Cauchy	30.90 ± 2.43	23.60 ± 1.02	29.40 ± 1.56	144.90 ± 9.26	29.50 ± 1.57	110711.90 ±	15.0 ± 0.0
							132580.43
	Uniform	1.0 ± 0.0	2.80 ± 0.60	1.30 ± 0.90	3.70 ± 5.44	1.70 ± 1.80	22.0 ± 0.0	15.0 ± 0.0
	Triangle	12.50 ± 0.92	13.60 ± 0.66	12.60 ± 0.66	48.0 ± 5.85	33.70 ± 8.25	32.10 ± 0.30	15.0 ± 0.0
	Triangle mixture	11.20 ± 0.75	12.00 ± 0.77	10.90 ± 0.70	42.30 ± 4.90	27.60 ± 8.0	30.80 ± 0.40	15.0 ± 0.0
	Gaussian mixture	28.90 ± 1.22	27.30 ± 1.49	27.00 ± 2.00	134.90 ± 9.37	100.40 ± 11.60	66.40 ± 2.42	15.0 ± 0.0

Table 9 :

 9 Results for the Lunar data set

		G-Enum	NML	BB	TS	RMG
			Kontkanen and	Scargle et al.	Davies et al.	Rozenholc et al.
			Myllymäki	(2013)	(2009)	(2010)
			(2007)
	Number of intervals	75	65	86	478	73
	Computation time (seconds)	1.3	1025	4489.26	53.59	37.02

Table 10 :

 10 Hellinger distances between the densities estimated by the methods on the Lunar Crater data set.

Preprint submitted to Computational Statistics & Data AnalysisDecember 22, 2022

Acknowledgments

We thanks the two anonymous reviewers and the associated editor for their insightful comments which help us to improve the quality of the manuscript.

Experimental protocol

The methods are evaluated on a collection of data sets. For artificial data set, we use samples of increasing size from n = 10 to n = 100, 000 or n = 1, 000, 000 samples when possible.

For each evaluation metric, we report the mean and standard deviation of results obtained over 10 independent samples (for each distribution and sample size).

Implementation

All experiments are carried out on a Windows 10 machine equipped with an AMD Ryzen 2-core processor and 6 GB of RAM memory.

The MDL criteria presented in this paper as well as the optimisation algorithms are implemented in C++. As proposed in Section 7.4, we restrict the search of the cut points of for NML criterion to approximations of the data points, even if Proposition 8 does not apply. This brings some computational efficiency to this method even when E is large.

Metrics

The comparison between histograms models is based on the number of intervals each histogram has, as well as the time it took to compute it. Additionally, the relevance of each histogram is evaluated by computing the Hellinger distance to the original model density.

The Hellinger Distance (HD) H(p, q) for p, q being probability density functions, is defined as

A HD close to 0 indicates a strong similarity between probability distributions. It is worth noting that most authors used the squared Hellinger distance, which may produce an impression of a better convergence of the estimations (Kontkanen and Myllymäki, 2007;[START_REF] Davies | A comparison of automatic histogram constructions[END_REF][START_REF] Luosto | Construction of irregular histograms by penalized maximum likelihood: a comparative study[END_REF]. The Hellinger distances presented in this paper are not squared. HD measures reported are obtained via numerical integration to estimate the probability distributions of the model density and of the histogram that models it.

Reference distributions

For the benchmarking experiments, we compare all methods over the 6 distributions described in table 4.

Comparison of MDL-based methods

In this subset of experiments, we focus on illustrating the similarities and differences between histograms produced by the NML, Enum and G-Enum criteria. The NML and Enum criteria were optimised either using the optimal dynamic programming algorithm presented in Kontkanen and Myllymäki (2007) or with the search heuristic presented in this paper.

We tested the NML, Enum and G-Enum criteria over the synthetic data sets described in table 4. We set the approximation accuracy to = 0.01 across all experiments on synthetic data for the NML and Enum methods. The accuracy of G-Enum is optimised as part of the estimation process. Detailed results are provided in the appendix for each distribution (figures E.7 to E.12). Section Appendix F provides in addition a study of the impact of on the results on a synthetic data set. In this section, we provide only the main insights we could gather from this benchmark.

Intervals. Across all data sets and all sample sizes, NML, Enum and G-Enum histograms have the same number of intervals. This confirms that the criteria are very similar. Moreover, we can see that choosing either the dynamic programming algorithm or the search heuristic for the optimisation results in almost identical interval counts. The final form of K&M's NML criterion in their paper is as follows.

To ease the comparison of both criteria, we simplified the likelihood term as described heareafter.

Proofs are not given in the same order as in the main text. Indeed proofs of propositions 5 and 6 are based on proposition 8 but for the clarity of exposure, we preferred to present this latter proposition before the former ones. Of course, the proof of proposition 8 uses neither proposition 5 nor proposition 6, which is emphasized by the natural mathematical order used in the present section.

In other words, an optimal histogram cannot contain zero-length intervals.

Proof. Let M = (K, (c k ) 0≤k≤K , (h k ) 1≤k≤K ) be a histogram compatible with D. Let us assume that c l-1 = c l for some l > 0 and let M = (K -1, (c k ) 0≤k≤K-1 , (h k ) 1≤k>K ) be the histogram derived from M by removing the zero-length interval. It is defined formally by

Proposition 4. The coding length of a histogram with two adjacent empty intervals is always longer than the coding length of a histogram with no consecutive empty intervals.

Proof. We analyse the cost variation when two adjacent intervals are empty for the Enum criterion (h A = 0, h B = 0, so

We have thus shown that, for the Enum criterion, the cost variation is always strictly negative after the merge of two adjacent empty intervals. This means that keeping a model with two adjacent empty bins is always more costly in terms of code length, so our search algorithm will systematically favor the merge of two consecutive empty bins.

Appendix B.5. Proof of proposition 7 Proposition 5. An optimal histogram has at most 2n -2 intervals (K * ≤ 2n -2).

Proof. The maximal number of non empty intervals is obviously bounded by the number of observations n. Provided the data set has n distinct observations and that is small enough, we can consider a histogram with n intervals containing each a single observation, separated by empty intervals. According to proposition 6 an optimal histogram cannot have consecutive empty intervals and thus the maximal number of separating intervals is n -1. Therefore an optimal histogram cannot have more than 2n -1 intervals.

But proposition 3 shows that a histogram with a single interval is preferred by the criterion over a histogram consisting only of empty intervals and intervals containing each a single observation. Therefore the maximal model with 2n -1 intervals cannot be optimal. Then the optimal number of intervals K * is bounded by 2n -2.

Appendix B.6. Proof of proposition 8 Proposition 8. In an optimal histogram, each endpoint is at most away from one of the values of the data set.

Proof. Let M = (K, {c k } 0≤k≤K , {h k } 1≤k≤K ) be an optimal histogram built from a data set D.

As by definition, c 0 = x min -/2 and c K = x max + /2, the property is valid for c 0 and c K .

Let us now assume that K > 1 and focus on endpoints c k , 0 < k < K, that is on cut-points between adjacent intervals.

Let i, j = i + 1 be the index of two such adjacent intervals. We study the location of the cut-point between intervals i and j, c i , while considering the exterior boundaries of intervals i and j (c i-1 and c j ) to be fixed. We denote L i, j = c jc i-1 . We consider also the data counts of the intervals h i and h j to be fixed. In other words, we study to what extent the cut-point between intervals i and j can be freely set at grid positions in the empty space between the last data point in i and the first data point in j, as illustrated by Figure B.5.

Let D i (resp D j ) be the data point in interval i (resp. j). We define

and

i.e. L i and L j are the minimum lengths of intervals that contain respectively D i and D j .

In terms of elementary -length bins, let E i = L i / , E j = L j / , E = L i, j / . We note x = l/ ∈ [0, E -E i -E j ] the number of -length bins at which the cut-point c i is placed.

Appendix B.7. Proof of proposition 5

Proposition 5. There exist data sets such that the optimal histogram has at least one interval which contains only a single observation.

Proof. Let us consider a data set with n observations and only two distinct values

We choose v A and v B such that v Bv A ≥ 3 . According to proposition 4 an optimal histogram cannot contain successive empty intervals and thus we need only to consider three cases: histograms with one, two or three intervals. Those cases are illustrated on Figure B.6.

For K = 1, we have

For K = 2, the optimal split is necessarily within of v A or v B according to proposition 8. If the split is close to v A ,

we have

If the split is close to v B , we have

For K = 3, the optimal split is necessarily right next to each value and we thus have two non-empty intervals composed E A = 1, E B = 1 -bins, surrounding an empty interval of E -2 -bins (this is again a consequence of proposition 8). We get

Therefore, even for rather small n and E, the model cost is minimal for K = 2 with the split within from v B , i.e. for the histogram where the first interval contains the singleton value v A .

Appendix B.8. Proof of proposition 6

Proposition 6. There exist data sets such that the optimal histogram has at least one interval that does not contain any observation.

Proof. Let us consider a data set with n observations and only two values v A , v B of frequencies h A = h B = n 2 . We choose v A and v B such that v Bv A ≥ 3 . For the same reasons as in the proof of proposition 5 we need only to consider three cases: histograms with one, two or three intervals.

For K = 1, we have

For K = 2, the optimal split is necessarily within of v A or v B according to proposition 8. Because of the symmetry of the setting, we can set the split to be near v B . Thus we have two intervals composed of

We get

To analyse this quantity, we use the approximation given in Grunwald (2007) (formula 4.36) which states that for

proof of Corollary 1

Corollary 1. Let D be a data set with n distinct observations. Then for sufficiently small, the optimal histogram build on bins for the Enum criterion is the trivial one with a single interval

Proof. Let us consider two optimal histograms M K, and M K , with K K . We have

As K (resp. K ) and S (M * K, , D) (resp. S (M * K , , D)) are integers, the minimum non zero value of

Using equation (C.6) to handle the two logarithmic terms with , we have c Enum (M 2, |D) =n log 1

Finally we have for K = 3 and

Using again equation (C.6), we have c Enum (M 3, |D) =n log 1

When converges to 0, c Enum (M k, |D)n log 1n log L reduces to log k plus some positive terms when k > 1. As log is an increasing function, when is mall enough, c Enum (M 1, |D) becomes smaller than c Enum (M k, |D) for k > 1 and thus M 1, is the optimal histogram.

Appendix D. Illustration of the asymptotic behaviour

We have:

39 where we used the same approximation from [START_REF] Grunwald | The minimum description length principle[END_REF] as in Section Appendix B.8. Using equation (C.5), we obtain

This can be simplified into

Finally we conclude using

Appendix E. Benchmarking and comparison of MDL methods

The first series of figures of the appendix compares the NML, Enum and G-Enum criteria. The three methods are evaluated on different sample sizes of a Normal (figure E.7, Cauchy (figure E.8), uniform (figure E.9), triangle (figure E.10), triangle mixture (figure E.11) and Gaussia mixture (figure E.12) distributions.

We highlight that all results are means over 10 experiments for each sample size. All results are shown in log scale. The standard deviations of the metrics are represented asymmetrically in the graphs to because of the log scale.

Appendix F. Role of the approximation accuracy,

As a complement to the theoretical analysis conducted in Section 6.2, we conducted a series of experiments to show the effect in practice of using a too small value for . We compare all MDL methods over a Normal distribution of small size (n = 100 and n = 1000), for ranging from 2.0 to 10 -8 . Figures F.13 and F.14 summarise the results. For both sample sizes we can observe that, as → 0, the number of intervals in NML and Enum histograms decreases. Particularly, for the n = 100 distribution, the estimations models start at 3 intervals and NML and Enum histograms slowly decrease to having a single one. This observation is in line with Corollary 1. We highlight that all results are means over 10 experiments for each sample size. All results are shown in log scale. The standard deviations of the metrics are represented asymmetrically in the graphs to because of the log scale.