{'fraction_non_alphanumeric': 0.0465421891255544, 'fraction_numerical': 0.024986767416817244, 'mean_word_length': 4.095796130952381, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Mode testing, critical bandwidth and excess mass

The identification of peaks or maxima in probability densities, by mode testing or bump hunting, has become an important problem in applied fields. This task has been approached in the statistical literature from different perspectives, with the proposal of testing procedures which are based on kernel density estimators or on the quantification of excess mass. However, none of the existing proposals provides a satisfactory performance in practice. In this work, a new procedure which combines the previous approaches (smoothing and excess mass) is presented and compared with the existing methods, showing a superior behaviour. A real data example on philatelic data is also included for illustration purposes.


Introduction
Simple distribution models, such as the Gaussian density, may fail to capture the stochastic underlying structure driving certain mechanism in applied sciences. Complex measurements in geology, neurology, economics, ecology or astronomy exhibit some characteristics that cannot be reflected by unimodal densities. In addition, the identification of the (unknown) number of peaks or modes is quite common in these fields. Some examples include the study of the percentage of silica in chondrite meteors (Good and Gaskins, 1980), the analysis of the macaques neurons when performing an attention-demanding task (Mitchell et al., 2007), the distribution of household incomes of the United Kingdom (Marron and Schmitz, 1992), the study of the body-size in endangered fishes (Olden et al., 2007) or the analysis of the velocity at which galaxies are moving away from ours (Roeder, 1990).
In all these examples, identifying the number (and location) of local maxima of the density function (i.e. modes) is important per se, or as a previous step for applying other procedures.
An illustrative example which has been extensively considered in mode testing literature can be found in philately (the study of stamps and postal history and other related items).
Research in this field has been motivated by the use of stamps for investment purposes.
Before the 1940 decade, stamps images were printed on a variety of paper types and, in general, with a lack of quality control in manufactured paper, which led to important fluctuations in paper thickness, being thin stamps more likely to be produced than thick ones. Given that the price of any stamp depends on its scarcity, the thickness of the paper is crucial for determining its value. However, there is not a standard rule for classifying stamps according to their thickness (not being available such a classification in catalogues), becoming this problem even harder in stamp issues printed on a mixture of paper types with possible differences in their thickness. A stamp issue with this problematic is the Hidalgo catalogue, from the nineteenth century, and this particular example has been shown in several references in the literature as a paradigm of the problem of determining the number of modes/groups. This example will be revisited in this work, recalling previous analysis and comparing results with the ones provided by the new testing procedure presented in this paper.
A formal hypothesis test for a null hypothesis of a certain number of modes can be stated as follows. Let f be the density function of a real random variable X and denote by j the number of modes. For k ∈ Z + , the testing problem on the number of modes can be formulated as: (1) There have been quite a few proposals in the statistical literature for solving (1) and the different techniques can be classified in two groups: a first group of tests based on or using a critical bandwidth, introduced by Silverman (1981), further studied by Hall and York (2001) and also used by Fisher and Marron (2001); and a second group of tests based on the excess mass, such as those ones proposed by Hartigan and Hartigan (1985), Müller and Sawitzki (1991) and Cheng and Hall (1998). These methods are briefly revised and compared in this paper, where a new proposal gathering strength from both areas is also introduced, outperforming the existing procedures, in testing unimodality and more general hypotheses.
Apart from the formal testing procedures, and as a complementary tool for them, a first step when confronting the problem of identifying modes in a data distribution is the exploration of a nonparametric estimator of the underlying probability density, which can be done by kernel methods. Classical kernel density estimation (Wand and Jones, 1995, Ch. 2) allows for the reconstruction of the data density structure without imposing parametric restrictions (only subjected to mild regularity assumptions) but at the expense of choosing an appropriate bandwidth parameter, which controls the degree of smoothing. A direct observation of a kernel estimator may lead to inaccurate or even wrong conclusions about the mode density structure. This can be easily noticed from the plots shown in Figure   1, where the kernel density estimator for the stamp dataset with different bandwidths is shown. Based on this estimator, from an exploratory perspective, there are several alternatives for identifying modes such as the SiZer map (Chaudhuri and Marron, 1999), the mode tree and the random forest (Minnotte and Scott, 1993;Minnotte et al., 1998).
Although these tools are helpful in supporting the results of formal testing procedures on the number and location of modes, apart from giving some insight on the global mode structure, the interpretation of the outputs from these procedures requires an expert eye.
This paper presents a new testing procedure combining the use of a critical bandwidth and an excess mass statistic, which can be applied to solve (1) in a general setting. In Section 2, a review on mode testing methods is presented, considering both tests based on critical bandwidth and on excess mass, jointly with the new proposal. A simulation study comparing all the procedures, in terms of empirical size and power, is included in Section 3. Section 4 is devoted to data analysis, revising the stamp dataset and presenting new results. Some final comments and discussion are given in Section 5. Details on the simulated models, a modification of the proposal when the compact support of the variable is known, a more flexible testing scenario and the computation of the proposed test are provided in the Appendix.

A review on multimodality tests
Different proposals for multimodality tests will be briefly revised in this section. Section 2.1 includes a review on the methods using the critical bandwidth, and excess mass approaches are detailed in Section 2.2. A new proposal, borrowing strength from both alternatives, is presented in Section 2.3: an excess mass statistic will be calibrated from a modified nonparametric kernel density estimator using a critical bandwidth.

Tests based on the critical bandwidth
For a certain number of modes k ∈ Z + , the critical bandwidth (Silverman, 1981) is the smallest bandwidth such that the kernel density estimator has at most k modes: h k = min{h :f h has at most k modes}, wheref h denotes the kernel density estimator, computed from a random sample X = (X 1 , ..., X n ), with kernel K and bandwidth h: (2) Silverman (1981) proposed to use the critical bandwidth with the gaussian kernel as a statistic to test H 0 : j ≤ k vs. H a : j > k, being its use justified by the fact that, with a gaussian kernel, the number of modes off h is a nonincreasing function of h. Hence, H 0 is rejected for large values of h k , whose distribution is approached using bootstrap procedures. Specifically, the proposed methodology consists in generating B samples By computing the critical bandwidth, h * b k , from each sample X * b , given a significance level α, the null hypothesis is rejected if P(h * k ≤ h k |X ) ≥ 1 − α. When k = 1, the testing problem tackled by Silverman (1981) coincides with (1). However, for a general k, the null hypothesis in (1) is more restrictive than the one considered by Silverman (1981), but asymptotic consistency of the test is only derived for j = k (for a deeper insight see Hall and York, 2001, or Appendix C). Hall and York (2001) proved that the previous bootstrap algorithm does not provide a consistent approximation of the test statistic distribution under the null hypothesis and suggested a way for accurate calibration when k = 1. Given a closed interval I where the null hypothesis is tested (f has a single mode in I), if both the support of f an the interval I are unbounded then properties of h 1 (critical bandwidth when k = 1) are generally determined by extreme values in the sample, not by the modes of f . To avoid this issue, the testing problem (1) is reformulated as follows: H 0 : j = 1 in the interior of a given closed interval I and no local minimum in I, (3) and the critical bandwidth is redefined accordingly as: h HY = min{h :f h has exactly one mode in I}. (4) An issue that should be kept in mind in the computation of this critical bandwidth is that even if K is the gaussian kernel, the number of modes off h inside I is not necessarily a monotone function of h. But under relatively general conditions (see Hall and York, 2001), the probability that number of modes is monotone in h converges to 1 for such a kernel.
Hall and York (2001)  Unfortunately, the critical bandwidths for the bootstrap samples h * b HY , are smaller than h HY , so for an α-level test, a correction factor λ α to compute the p-value P(h * HY ≤ λ α h HY |X ) ≥ 1 − α must be considered. Two different methods were suggested for computing this factor λ α , the first one based on a polynomial approximation and a second one using Monte Carlo techniques considering a simple unimodal distribution.
The previous proposal could be extended, as mentioned by Hall and York (2001), to test that f has exactly k modes in I, against the alternative that it has (k + 1) or more modes there, extending the critical bandwidth in (4) for k modes, namely h HY,k . Nevertheless, in this scenario, the bootstrap test cannot be directly calibrated under the hypothesis that f has k modes and (k−1) antimodes, since it depends on the (2k−2) unknowns (c i /c 1 ), where , and t i being the ordered turning points of f in I with i = 1, · · · , (2k − 1); which notably complicates the computations.
Finally, it should be also commented that the use of the critical bandwidth for testing (1) is not limited to its use as a test statistic. Consider a Cramér-von Mises test statistic: where F 0 is a given continuous distribution function, {X (1) ≤ · · · ≤ X (n) } denotes the ordered sample and F n is the empirical distribution function. Fisher and Marron (2001) proposed the use of (5) for solving the general problem of testing k modes (H 0 : j ≤ k) by and derived the statistic: where the null hypothesis is rejected for large values of T k . To approximate the distribution of the test statistic (6) under the null hypothesis, a bootstrap procedure is also proposed.
It will be seen in Section 3 that the behaviour of the Fisher and Marron (2001) proposal is far from satisfactory. Müller and Sawitzki (1991) confront the testing problem (1), employing a different perspective, under the following premise: for k ∈ Z + , a mode is present where an excess of probability mass is concentrated. Specifically, given a continuous real density function f and a constant λ, the excess mass is defined as:

Tests based on excess mass
where C λ = {x : f (x) ≥ λ} and ||C λ || denotes the measure of C λ . If f has k modes, the excess mass can be defined as: where the supremum is taken over all families {C m : m = 1, · · · , k} of pairwise disjoint connected sets over the support of f . Under the assumption that f has k modes, E(f, λ) can be empirically estimated just replacing P(C m ) by P n (C m ) = (1/n) n i=1 I(X i ∈ C m ), being I the indicator function. The difference D n,k+1 (λ) = E n,k+1 (P n , λ) − E n,k (P n , λ) measures the plausibility of the null hypothesis, that is, large values of D n,k+1 (λ) would indicate that H 0 is false. Using these differences, Müller and Sawitzki (1991) proposed the following test statistic: rejecting the null hypothesis that f has k modes for large values of ∆ n,k+1 . Müller and Sawitzki (1991) showed that this statistic is an extension of the dip test introduced by Hartigan and Hartigan (1985), just valid for the unimodal case, since both quantities (dip and excess mass) coincide up to a factor, for the unimodality case. In addition, the proposal of Müller and Sawitzki (1991) for testing unimodality is the same as that one of Hartigan and Hartigan (1985) and considers a Monte Carlo calibration, generating resamples from the uniform distribution.
In view of the extremely conservative behaviour of the calibration of the previous proposals (see Section 3 for results), Cheng and Hall (1998) designed a calibration procedure based on the following result: for large samples and under the hypothesis that f is unimodal, the distribution of ∆ n,2 is independent of unknowns except for a factor c = (f 3 (x 0 )/|f (x 0 )|) 1/5 , where x 0 denotes the unique mode of f . Using this fact, for the case k = 1, Cheng and Hall (1998) approximated the distribution of ∆ n,2 employing the values of ∆ * n,2 obtained from the samples generated from a parametric calibration distribution Ψ(·, κ), being κ a certain parameter. Depending on the value of d = c −5 , different parametric distributions were suggested by the authors: a normal (d = 2π), a beta (d < 2π) or a rescaled Student t (d > 2π). For estimating d, if x 0 denotes the largest mode off h , then is used, weref andf are kernel estimators with a gaussian kernel and h and h are their respective asymptotically optimal global bandwidths, replacing the unknown quantities for the ones associated with a N (0,σ 2 ), beingσ 2 the sample variance.
The methodology proposed by Cheng and Hall (1998) consists in generating samples from Ψ(·,κ), whereκ and the distribution family are chosen usingd. The excess mass statistic given in (7) when k = 1, that is ∆ * n,2 , is computed from the resamples and, for a given significance level α, the null hypothesis is rejected if P(∆ * n,2 ≤ ∆ n,2 |X ) ≥ 1 − α.

A new proposal
The previous tests show some limitations for practical applications: first, just the proposals of Silverman (1981) and Fisher and Marron (2001) allow to test (1) for k > 1. Despite the efforts of Cheng and Hall (1998) and Hall and York (2001) for providing good calibration algorithms, it will be shown in Section 3 that the behaviour of all the proposals is far from satisfactory. Specifically, the test presented by Silverman (1981) is very conservative in general (although sometimes can show the opposite behaviour) and the proposal of Fisher and Marron (2001) does not have a good level accuracy. The new method proposed in this work overcomes these drawbacks by considering an excess mass statistic, as the one proposed by Müller and Sawitzki (1991) with bootstrap calibration. Unlike Cheng and Hall (1998), a completely data-driven procedure will be designed, using the critical bandwidth The proposal, in a nutshell. Consider then the testing problem (1) and take the excess mass statistic given in (7), under the null hypothesis. Given X , generate B resamples X * b (b = 1, . . . , B) of size n from (a modified version of)f h k . For a significance level α, the null hypothesis will be rejected if P(∆ * n,k+1 ≤ ∆ n,k+1 |X ) ≥ 1 − α, where ∆ * n,k+1 is the excess mass statistic obtained from the generated samples. It should be also noted that the procedure can be easily adapted to handle Hall and York (2001) scenario: to test the null hypothesis that f has almost k modes in the interior of a given closed interval I, if I is known, use (a modified version of)f h HY,k to generate the samples. From this brief description, two questions arise: How is this modified version off h k constructed? Does the procedure guarantee a correct calibration of the test? In fact, the construction of the modified kernel density estimator, namely the calibration function and subsequently denoted by g, is intended to ensure the correct calibration, under some regularity conditions. Under these premises, Cheng and Hall (1998) indicated that the distribution of the test statistic (7) is independent from the underlying density f , except for the values in the modes and antimodes of d i = |f (x i )|/f 3 (x i ) with i = 1, . . . , (2j − 1). It can be seen, assuming that f has k modes, that the distribution of ∆ n,k+1 can be approximated by ∆ * n,k+1 , just with bootstrap values obtained from samples coming from a calibration distribution with k modes. This calibration distribution must satisfy that the estimated values d i converge in probability to d i , as n → ∞, for i = 1, . . . , (2k − 1). As mentioned before, it is for that reason that Cheng and Hall (1998) proposed using a parametric family depending on the value of d when k = 1. Two issues appear related with their calibration procedure. First, it is not an easy task to construct a parametric family having the desired values of d i , for i = 1, . . . , (2k − 1), when k > 1. In addition, the second-order limiting properties of the test depend on the form of the density function. Then, a better behaviour is expected if the calibration function is "more similar" to the real density function. Our method deals with these two issues to get a test having a good performance in the finite-sample case and allowing to solve the general problem of testing k modes.
In order to estimate the location of modes and antimodes, as well as the value of the density in these points under H 0 , the kernel density estimator (2) with critical bandwidtĥ f h k is a good candidate. In addition, if f is also a bounded density with bounded support Mammen et al. (1992) holds and the critical bandwidth is of order n −1/5 , being then optimal for estimating f . Unfortunately, a suitable bandwidth for estimating f cannot be used for effectively estimating f , being, the optimal bandwidth for estimating the second derivative of the density is of order n −1/9 . For this reason, the kernel density estimator with critical bandwidth h k cannot be used as calibration function g, which needs to provide a good estimator of d i . Under the conditions of Mammen et al. (1992), an appropriate estimator can be calculated using the ratio where h PI is a plugin bandwidth. In this work, the plug-in rule for the second derivative will be obtained deriving the asymptotic mean integrated squared error and replacing f in its expression using a two-step procedure (see, for example, Wand and Jones, 1995, Ch. 3).
The calibration function g will be constructed by modifying f h k in a neighbourhood of {x : f h k (x) = 0}, being such values a finite collection (see Silverman, 1981), having positive estimated density. The expression of this calibration function g is given in (8) and explained below. Depending on the nature of these points, two modifications in their neighbourhood will be done. If the point x i is a mode or an antimode of f h k then the modification will preserve its location and its estimated density value, and its second derivative will satisfy g ( x i ) = f h PI ( x i ) 2 . For i = 1, . . . , (2k − 1), this will be the role of the function J, in the neighbourhood (r i , s i ) of x i for guaranteeing |g ( x i )|/g 3 ( x i ) = d i and satisfying (RC4). The second modification will remove the t saddle points of f h k , denoted as ζ p , with p = {1, . . . , t}.
In order to ensure that g satisfies (RC3), in a neighbourhood (z (2p−1) , z (2p) ) of ζ p , a function L will be introduced. Finally, as f h k verifies (RC1) and (RC2), and the modifications of the functions J and L will be done in bounded neighbourhoods and preserving condition (RC1), g will satisfy also these two conditions. Then, the calibration function g will be constructed as follows: where ς has k components ς i ∈ (0, 1/2), with i = 1, . . . , k, determining at which height of the kernel density estimation the modification is done. Values of ς i close to 0 imply a modification in an "small" neighbourhood around the mode or antimode. Note that a little abuse of notation was made as g will depend on the function f h k (not only on h k ) and on the values f h PI ( x i ), for i = 1, . . . , (2k − 1). An example of the effect of g can be seen in Figure 2 and it will be further described. Before defining functions J and L, to ensure that g has continuous derivative, a link function l must be introduced: where a 0 = a 1 and v > u. Two issues must be noticed in this function. are the same, then the first derivative of l will not be equal to 0 for any point inside [u, v].
The form of J is given in equation (10) and its construction guarantees that x i is the unique point in which the derivative is equal to 0 in the neighbourhood where it is defined.
The construction of J is achieved with the K function defined bellow and properly linked with the link function (9) to preserve (RC1). The K function is defined as follows being δ i a value indicating if x i is a mode (δ i = −1) or an antimode (δ i = 1). The value η i will be defined later and it will depend on ς i . The second derivative of this function exists and is Hölder continuous in ( The following equalities are also satisfied: , the J function can be defined as follows being v i = x i − η i /2 and w i = x i + η i /2. As it was mentioned, the function J described in (10) (and hence also the calibration function g) depends on the constant ς i ∈ (0, 1/2).
Ordering the modes and denoting as x 0 = −∞ and the remaining unknowns values in (10) will be obtained as follows.
First, it is necessary to decide at which height ϑ i the modification in f h k is done. For values of ς i close to 0, ϑ i will be close to f h k ( x i ); while for values close to 0.5, ϑ i will be in the middle point between f h k ( x i ) and the highest (or lowest if x i is an antimode) value of f h k in the two closest modes or antimodes ( x i−1 and x i+1 ). Second, once the height is decided, r i and s i will be the left and the right closest points to x i at which the density estimation is equal to ϑ i . Third, in order to link correctly the K function, it is necessary to define η i ensuring that K( x i ± η i /2; ρ i , η i ) will be higher (lower in the antimodes) than ϑ i . With this objective, η i is chosen in such a way that K( and as close as possible to the middle point between ϑ i and f h k ( x i ). Also, the value η i will ensure that the where the two link functions are employed. An example of the modifications achieved by the J function in the modes and antimodes of f h k is shown in Figure 2 and the complete characterization is provided bellow In order to proceed with the modification achieved with the L function, assume that this estimator has t saddle points ζ p , with p = 1, . . . , t. Define as ξ = min{|x − y| : x, y ∈ (ζ 1 , . . . , ζ t ) ∪ (r 1 , s 1 , . . . , r 2k−1 , s 2k−1 )}. Then, if ζ p is not inside the interval where the J funtions are defined, the neighbourhood used to remove the stationary and turning points will be delimited by z (2p−1) = ζ p − ξ and z (2p) = ζ p + ξ, with ∈ (0, 1/4). In the simulation study, the value of will be taken close enough to 0 to avoid an impact in the value of the integral associated to g. Once these points are calculated, the saddle points can be removed from g with the link function by taking L equal to To construct the calibration function, first, values of ς i ∈ (0, 1/2), for i ∈ {1, . . . , (2k − 1)} must be fixed. Then, using the J function (10) with the values given in (11) and the L function (12), the function g defined in (8) satisfies the specified regularity conditions and With this modification the calibration function also preserves the structure of the data under the hypothesis that f has k modes. However, this calibration function g may not be a density function since may not be equal to 1. To guarantee that g is indeed a density, a possible approach consists in proceeding with a search of values for ς such that q(ς) is equal to 1, which is guaranteed For convenience, in the simulation study, the employed approach will be followed using ς i close enough to 0 (∀i ∈ {1, . . . , 2k − 1}) in order to avoid an impact on the integral value.
To conclude this section it is important to mention the necessary condition of the bounded support in order to ensure the results of Mammen et al. (1992). Although, in general, it will not have an important effect as it can be seen in Section 3, in case of a known support, h HY,k can be used to create the calibration function. An alternative approach for known support is given in Appendix B. Once a conclusion about the number of modes is made, when both the estimation of the modes and antimodes is needed, this known support plays an special role. From the results of Romano (1988) and Mammen et al. (1992), it can be proved thatf h HY,k provides a good estimation of their localization.
Also, it should be reminded that, unlike Silverman (1981) and Fisher and Marron (2001) proposals, the one presented in this paper considers The aim of the following simulation study is to compare the different proposals presented in Section 2. Samples of size n = 50, n = 200 and n = 1000 (n = 100 instead of n = 1000 in power studies) were drawn from twenty five different distributions, ten of them unimodal (models M1-M10), ten bimodal (models M11-M20) and five trimodal (models M21-M25) (see Appendix A). For each choice of sampling distribution and sample size, 500 realizations of the sample were generated. Conditionally on each of those samples, for testing purposes, 500 resamples of size n were drawn from the population. Tables 1-5 report the percentages of rejections for significance levels α = 0.01, α = 0.05 and α = 0.10 under different scenarios: testing unimodality vs. multimodality (Tables 1 and 2); testing bimodality against more than two modes (Table 4) and power analysis (respectively Tables   3 and 5). The procedures considered include the proposals by Silverman (1981)  is used both for HY and for NP, when the interval containing the modes is assumed to be known (Table 6). Further computational details are included in Appendix E.
Testing unimodality vs. multimodality. From the results reported in Tables 1 and 2, it can be concluded that SI is quite conservative: even for high sample sizes, the percentage of rejections is below the significance level, and quite close to 0 even for α = 0.10.
An exception to this conservative behaviour can be found in models M5 (with outliers in tails) and M6 (with a flat tail ), finding percentages of rejections between 0.03 and 0.052 for n = 200 or n = 1000, when α = 0.10. Regarding FM, a systematic behaviour cannot be concluded: the percentage of rejections is above the significance level for models M1, M5, M7, M9 or M10, but it can be also below the true level for M2, M4 or M8.
The behaviour of HY is quite good when using I = [0, 1] for the different distributions and large sample sizes. For n = 1000, the percentage of rejections is quite close to α, except for model M5 (for α = 0.05, below level) and for models M6 and M7 (for α = 0.10, above level). However, the percentage of rejections is usually below the significance level for small sample sizes. Exceptions to this general pattern are found for model M1 (n = 200), M3 (n = 50) and M10 (n = 200), where percentage of rejections is close to α and models M3 (n = 200), M6 (n = 200) and M7 with percentages above α. Nevertheless, it should be kept in mind that the support where unimodality is tested must be known. Similarly to SI, the results obtained with HH are quite conservative. For instance, for n = 1000, even taking α = 0.10, the percentage of rejections is always below 0.002.
Calibration seems correct in simple models for CH, although slightly conservative in some cases, such as for models M4 (n = 1000), M5 (n = 200) and M8 (n = 200). As expected, the parametric calibration distributions do not capture, for example, the skewness and this affects the second-order properties in more complex models. This effect is reflected in the asymmetric M3, M7 and M10 (n = 1000), or model M9, where the percentage of rejections is below α, and for M6 where is considerably higher than the significance level.
Finally, regarding the new proposal NP, it can be concluded that the calibration is quite satisfactory, even for complicated models, with a slightly conservative performance for M3 (n = 200), M4 (n = 1000) or M7 (n = 200), being this effect more clear for model M9. The only scenario where the percentage of rejections is above α is for M6 with n = 200, but this behaviour is corrected when increasing the sample size. Although the performance is better for higher sample sizes (n = 200 or n = 1000), in some cases, such as M9 or M16 (in the bimodal case), it can be seen that even for n = 1000, a percentage of rejections close to α is hard to get. In this difficult cases, the knowledge of the support can be used for obtaining better results as it was reported in Table 6, where the percentage of rejections is close to α for the sample sizes n = 200 and n = 1000.
Regarding power behaviour (and just commenting on the three methods which exhibit a correct calibration), results are reported in For testing bimodality, NP presents good results. The percentage of rejections is close to the significance level, except for M12 (n = 200) and M13 (n = 200), slightly below α, and M11 (n = 50), M15 (n = 50, n = 200), M16 and M19 (n = 50), slightly above α. For n = 1000, all the results are good except for M16, but the calibration problem is corrected (as seen in Table 6) applying NP with known support, taking for that purpose I = [0, 1].
So, just the new proposal presents a correct calibration. Hence, power results reported in Table 5 are only judged for the new proposal: power increases with sample size, detecting that all the alternative distributions do not satisfy the null hypothesis, except for M21 (n = 50) and M24 (n = 50), but with small sample sizes.

Real data analysis
As explained in the Introduction, the value of stamps depends on its scarcity, and thickness is determinant in this sense. However, in general, the designation of thick, medium or thin stamps is relative and can only refer to a particular stamp issue. Otherwise, making uniform categories for all stamp issues may lead to inaccurate classifications. In addition, there is not such a differentiation between groups available in stamps catalogs, leaving this classification to a personal subjective judgment. The importance of establishing an objective criterion specially appears in stamp issues printed on a mixture of paper types, with possible differences in their thickness.
A stamp issue where the problem of determining the number of different groups of stamps appears is in the 1872 Hidalgo issue. First, for this particular issue, and in general in the Mexican ones, it is known that the handmade paper presents a high variability in the thickness of the paper. Second, since of scarcity of ordinary white wove paper, other types of paper were used to produce the Hidalgo issue. A small quantity of "vertically laid" paper, a fiscal type of white wove paper denominated Papel Sellado (some of them were watermarked vertically), other type of white wove, the La Croix-Freres of France (some of them with a watermark of LA+-F ) and also another unwatermarked white wove paper might also have been used. It is estimated than the watermark of Papel Sellado            watermark, it is estimated that the symbol appears in between 4 and 10 stamps if the sheet was watermarked, and some authors suggested that this watermark appears only once of every 4 sheets. In order to get more information about this particular problem and to obtain some further references, see Izenman and Sommer (1988).
This particular example has been explored in several references in the literature for determining the number of groups. From a non-parametric point of view, some examples of its utilization for mode testing can be shown in Efron and Tibshirani (1994, Ch. 16), Izenman and Sommer (1988) or in Fisher and Marron (2001). Also it was analysed using non-parametric exploratory tools in Wilson (1983), Minnotte and Scott (1993) and in Chaudhuri and Marron (1999). Some parametric studies of the 1872 Hidalgo issue can be found in Basford et al. (1997) or in McLachlan and Peel (2000, Ch. 6).
Taking a subsample of 437 stamps on white wove, Wilson (1983) made a histogram and the conclusion was that only two kinds of paper were used, the Papel Sellado and the La Croix-Freres, and that there was not a third kind of paper. Izenman and Sommer (1988) revisited the example considering a more complete collection, with 485 stamps. A histogram with the same parameters as those used by Wilson (1983) (same starting point and bin width) is shown in Figure 1 (top-left panel), revealing the same features as those noticed in the original reference. Two groups are also shown by a kernel density estimator, shown in the same plot, considering a gaussian kernel and a rule of thumb bandwidth (see Wand and Jones, 1995, Ch. 3.2). However, both approximations (histogram and kernel density estimator) depend heavily on the bin width and bandwidth, respectively.
Specifically, the use of an automatic rule for selecting the bandwidth value (focused on the global estimation of the entire density function) does not guarantee an appropriate recovery of the modes. In fact, using another automatic rule as the plug-in bandwidth (Figure 1, bottom-left panel), nine modes are observed. A histogram with a smaller bin width is also included in this plot, exhibiting apparently more modes than the initial one.
Given that the exploratory tools did not provide a formal way of determining if there are more than two groups, Papel Sellado and La Croix-Freres, Izenman and Sommer (1988) employed the multimodality test of Silverman (1981). Note that for this purpose, just FM, SI and the new proposal NP can be used, and the first two proposals present a poor calibration, as shown in the simulation study. Results from Izenman and Sommer (1988), applying SI with B = 100, are shown in Table 7 (note that with B = 500, different pvalues are obtained). For α = 0.05, the conclusions are the same, except in the crucial case of testing H 0 : j ≤ 2, where for B = 500, there are no evidences to reject the null hypothesis. These differences may be caused by the approximations implemented by Izenman and Sommer (1988) to obtain the critical bandwidth. Both Efron and Tibshirani (1994, Ch. 16) (using B = 500 bootstrap replicates) and Salgado-Ugarte et al. (1998) (employing B = 600) obtained similar results to ours. Hence, the null hypothesis must not be rejected when the hypothesis is that the distribution has at most two modes, but it has to be rejected when H 0 is that the distribution has at most six modes. This strange behaviour also happens in Izenman and Sommer (1988) analysis, when testing H 0 : j ≤ 3 and H 0 : j ≤ 6. Izenman and Sommer (1988) suggested non-rejecting the null hypothesis the first time that the p-value is higher than 0.4. The consideration of a flexible rule for rejecting the null hypothesis is justified by the fluctuations in the p-values of SI and, as Izenman and Sommer (1988) mentioned, by the "conservative" nature of this test. Under this premise, the result when applying SI would be that the null hypothesis is rejected until it is tested H 0 : j ≤ 7. Hence, Izenman and Sommer (1988) conclude that the number of groups in the 1872 Hidalgo Issue is seven.
As shown in Section 3, SI does not present a good calibration and sometimes it can be  Table 7: P-values obtained using different proposals for testing k-modality, with k between 1 and 9. Methods: SI, FM and NP. For SI method, B = 100 (first row; Izenman and Sommer, 1988) and B = 500.
also anticonservative. It is not surprising that SI behaves differently when testing H 0 : j ≤ k for k = 2, 3, with respect to the rest of cases until k = 7. Since NP has a good calibration behaviour, even with "small" sample sizes, this method is going to be used, first for testing the important case H 0 : j = 2 vs. H a : j > 2 and then to figure out how many groups are there in the 1872 Hidalgo Issue.
The computation of the excess mass statistic requires a non-discrete sample and the original data (denoted as X ) contained repeated values, the artificial sample Y = X + E will be employed for testing the number of modes, where E is a sample of size 485 from the U (−5 · 10 −4 , 5 · 10 −4 ) distribution. This modification of the data was also considered by Fisher and Marron (2001). The p-values obtained in their studio (using B = 200 bootstrap replicates) are shown in Table 7: it is not clear which conclusion has to be made. They mentioned that their results are consistent with the previous studies, detecting 7 modes.
But it should be noticed that, as shown in the simulation, FM does not present a good calibration behaviour.
Finally, the p-values obtained with NP are also shown in the Table 7, with B = 500.
Similar results can be obtained employing the interval I = [0.04, 0.15] in NP with known support, as Izenman and Sommer (1988) notice that the thickness of the stamps is always in this interval I. Employing a significance level α = 0.05 for testing H 0 : j = 2, it can be observed that the null hypothesis is rejected. It can be seen that the null hypothesis is rejected until k = 4, and then there is no evidences to reject H 0 employing greater values of k. Then, applying our new procedure, the conclusion is that the number of groups in the 1872 Hidalgo Issue is four.
In order to compare the results obtained by Izenman and Sommer (1988) and the ones derived applying the new proposal, two kernel density estimators, with gaussian kernel and critical bandwidths h 4 and h 7 are depicted in Figure 1 (bottom-right panel). Izenman and Sommer (1988) conclude that seven modes were present, and argued that the stamps could be divided in, first, three groups (pelure paper with mode at 0.072 mm, related with the forged stamps; the medium paper in the point 0.080 mm; and the thick paper at 0.090 mm). Given the efforts made in the new issue in 1872 to avoid forged stamps, it seems quite reasonable to assume that the group associated with the pelure paper had disappeared in this new issue. In that case, the asymmetry in the first mode using h 4 can be attributed to the modifications in the paper made by the manufacturers. Also, this first and asymmetric group, justifies the application of non-parametric techniques to determine the number of groups. It can be seen, in the Section 7 of Izenman and Sommer (1988) and in other references using mixtures of gaussian densities to model this data (see, for example, McLachlan and Peel, 2000, Ch. 6), that these parametric techniques have problems capturing this asymmetry, and they always determine that there are two modes in this first part of the density, one near the point 0.07 mm and another one near 0.08 mm. For the two modes near the points 0.10 and 0.11 mm, both corresponding to stamps produced in 1872. As Izenman and Sommer (1988) noticed, it seems that the stamps of 1872 were printed on two different paper types, one with the same characteristics as the unwatermarked white wove paper used in the 1868 issue, and a second much thicker paper that disappeared completely by the end of 1872. Using this explanation, it seems quite reasonable to think that the two final modes using h 4 , corresponds with the medium paper and the thick paper in this second block of stamps produced in 1872. Finally, for the two minor modes appearing near 0.12 and 0.13 mm, when h 7 is used, Izenman and Sommer (1988) do not find an explanation and they mention that probably they could be artefacts of the estimation procedure. This seems to confirm the conclusions obtained with our new procedure. The reason of determining more groups than the four obtained with our proposal, seems to be quite similar to that of the model M20 in our simulation study. This possible explanation is that the spurious data in the right tail of the last mode are causing the rejection of H 0 , when SI is used.

Discussion
Determining the number of modes in a distribution is a relevant practical problem in many applied sciences. The proposal presented in this paper provides a good perfomance for the testing problem (1), being in the case of a general number of modes k the only alternative with a reasonable behaviour. The totally nonparametric testing procedure can be extended to other contexts where a natural nonparametric estimator under the null hypothesis is available. For instance, the method can be adapted for dealing with periodic data, as it happens with the proposal by Fisher and Marron (2001).
In practical problems, where a large number of tests must be computed, obtaining a set of p-values is a crucial task. In this setting, such a computation should be accompanied by the application of FDR correction techniques. The proposal in this work, based on the use of critical bandwidth and excess mass ideas, and its combination with FDR, is computationally feasible. With the aim of making this procedure accessible for the scientific community, and therefore, enabling its use in large size practical problems, an R package has been developed.

A Models for simulation study
The specific formulas of those models considered in the simulation study carried out in Section 3 are given here with the notation l i=1 p i · ψ i , where each ψ i represents one of the component of the mixture and p i are the weights of these different components, with i = 1, . . . , l, satisfying l i=1 p i = 1. The unimodal density functions used as ψ i are the following models, as defined in Johnson et al. (1995): i ) and Weibull(δ i , c i ). All the models were created in such a way that f (0) ≈ f (1) ≈ 0.1 max x∈(0,1) f (x). The unimodal probability density functions are represented in Figure   3, the bimodal and trimodal models appear in Figure 4.

B New proposal when the support is known
When the bounded support is known, an alternative approach for the new proposal can be used in order to get better results in practice. This new proposal consists in replacing the critical bandwidth of Silverman (1981) for the one of Hall and York (2001) in the definition of the calibration function g. If the number of modes in the entire support is equal to k when testing H 0 : j = k (with (k − 1) antimodes) in [a, b], then no more changes are needed. If modes appear outside [a, b], then the link function (9) can be used in order to preserve the required regularity conditions. Denoting as a < x 1 < . . . < x 2k−1 < b, being x 1 and x 2k−1 modes, the points x 0 and x 2k , needed to obtain the values in (11), will be redefined to remove the modes outside [a, b]. If there are modes lower than x 1 , then x 0 = min{x : x ≥ a and f h HY,k (x) > 0} and if there are modes greater than x 2k−1 , then x 2k = max{x : x ≤ b and f h HY,k (x) < 0}. Once this change is done, two extra values, a < x 0 and b > x 2k , are needed in order to use the link function. The steps to obtain these two values will be defined later on and, from them, the g function defined in (8) x ≤ a and f h HY,k has modes lower than a, x 0 ) and f h HY,k has modes lower than a, (if it has modes greater than b) to obtain an integral equal to one. An option is to employ a and b satisfying It can be the case where the equality (14) is not satisfy for any pair (a, b), being a ∈ (−∞, x 0 ) and b ∈ ( x 2k , ∞). If this happens, the calibration function can be divided by the normalizing constant to correct the value of the integral. Another alternative can be to take other values of x 0 < x 1 and x 2k > x 2k−1 , such as f h HY,k (x) > 0, for all x ∈ [ x 0 , x 1 ), The approach considered in the simulation study (when the support is known) is, in the tails, try to find the value of a in the interval [ x 0 − b + a, x 0 ) and the value of b in If for all the possible values of a and b the integral is not equal to 1, then the solution is take a and b in such a way that q 2 is as close as possible to one and, then, employ the quotient g(·; h HY,k , h PI , ς, a, b)/q 2 as the calibration function.
C Testing k-modality when the true density has less than k modes As it has already been mentioned, the methods proposed by Silverman (1981) andMarron (2001) can be extended from unimodality to test a general null hypothesis as H 0 : j ≤ k. Nevertheless, the proposal presented in this work just allows to test The reason why the k-modal test should not be used when the true underlying density has less than k modes is that the test statistic in the bootstrap resamples converge in distribution to a random variable, depending only on the values d i with i = 1, . . . , (2k − 1) (see Cheng and Hall, 1998). When j < k, in the calibration function g, there exist (2k − 2j) turning points that they will not converge to any fixed value depending on the real density function. As the (asymptotic) distribution of the test statistic in the bootstrap resamples depends also on this (2k − 2j) values, one would expect that the sample distribution of the test statistic will not be correctly approximated with the bootstrap resamples.
Testing H 0 : j = k instead of H 0 : j ≤ k is not in general an important limitation for practical purposes. As it is done in the stamp example in Section 4, the usual procedure is to perform a stepwise algorithm starting with one mode and, if the null hypothesis is rejected, increasing the number of modes in the null hypothesis by one until there is no evidences for rejection. Despite this note of caution, it can be seen that, generally, testing H 0 : j = k, when j < k, reports also good calibration results.
In order to show the accuracy in practice when the bimodality test is employed in unimodal cases, Table 8 reports the percentages of rejections for significance levels α = 0.01, α = 0.05 and α = 0.10 for testing bimodality employing the proposals by Silverman (1981) (SI), Fisher and Marron (2001) (FM) and the new proposal (NP) presented in this paper.
With this goal, samples of size n = 50, n = 200 and n = 1000 were drawn from 10 unimodal distributions (models M1-M10). Again, for each choice of sampling distribution and sample size, 500 realizations were generated. Conditionally on each of those samples, for testing purposes, 500 resamples of size n were drawn from the population.
The conclusions from the results reported in Table 8 are quite similar to those given previously in Section 3 when H 0 : j = 1 was tested. First, although SI still reports a percentage of rejections below the significance level, it is less conservative than for the initial results (testing H 0 : j = 1). For all the available models and all sample sizes and significance levels, the percentage of rejections employing the bimodality test is greater or equal than the one obtained when the unimodality test was applied but lower than the significance level. In fact, taking into account the two largest sample sizes (n = 200 and n = 1000), the percentage of rejections is always strictly greater in the bimodal test with the exceptions of the models M3 (n = 50 and α = 0.05), M6 (n = 200 and α = 0.10) and M10 (n = 50 and α = 0.05). Regarding FM, again a systematic behaviour cannot be concluded, being the results similar to those ones reported in Section 3 when unimodality was tested on the same models. Finally, the results for the new proposal seem to be again quite satisfactory, with a slightly conservative performance in some models, such as M2, M3 (n = 1000 and α = 0.10), M4, M5 (n = 200), M7, M9 and M10 (n = 50). Observing these results, it seems that, in practice, NP can be used for testing H 0 : j ≤ k, but it should be kept in mind that a correct calibration is not guaranteed. An example of poor behaviour can be observed for (unimodal) model M26. Analysing the results reported in Table 9 for NP, it can be seen that, for n = 1000, when testing H 0 : j = 1 the percentage of rejections is close to the significance level, whereas when testing H 0 : j = 2, the percentage of rejections is bellow α, even employing the correction provided when the support is known.

D Some theoretical results
The complementary results needed to prove thatd i converge in probability to d i , as n → ∞, for i = 1, . . . , (2j − 1), where j is the number of modes of f , are given in this part. These proves are based on the results of Romano (1988). In order to provide the different results, some further conditions must be considered: (C1) f is a bounded and derivable density function.
(C2) The second derivative of f is continuous in a neighbourhood of the order points (C3) f has a bounded support [a, b].
(C6) Denoting as a n and b n the fixed sequence of numbers such as (nb 2 n )/(log(n)a n ) → ∞ and a n → 0. The bandwidth h (n) satisfies b n ≤ h (n) ≤ a n and f h (n) has j modes and (j − 1) antimodes.   Table 9: Percentages of rejections for testing H 0 : j = 1 (first column) and H 0 : j = 2 (second and third column), with 500 simulations (1.96 times their estimated standard deviation in parenthesis) and B = 500 bootstrap samples.
Although in this document the results are given under more general conditions, for the important case of the critical bandwidth h j , employing the results of Mammen et al. (1992), it can be seen that under (C1)-(C5), h j always satisfy (C6). Based on this fact, combined with previous results, it can be concluded that, under (C1)-(C5),d i converge in probability to d i , as n → ∞, for i = 1, . . . , (2j − 1). Note that for having this last result no more conditions that those ones given as Assumptions in Mammen et al. (1992) are needed.
To prove the results given in (R1) and (R2), in this document, the focus will be when x i is a mode, that is, i = 1, 3, . . . , 2j − 1. The proof for the antimodes is analogue to the given one. To begin the proof, first, it will be introduced the function Taking x 0 = 2a − x 1 and x 2j = 2b − x 2j−1 (note that condition (C3) is not needed in this case, if the support is unbounded these values can be replace conveniently with other values x 0 < x 1 and x 2j > x 2j−1 ), let δ i be the values such as δ i ≤ min{(x i − x i−1 )/2, (x i+1 − x i )/2} and also the following conditions are satisfied (see Romano, 1988 Combining (15) and (16), it is obtained that Using that by definition, fixed i ∈ {1, 3, . . . , 2j − 1}, for any δ ∈ (0, δ i ), the continuity of f , its Taylor expansion and Combining (16) and (18) Denoting asl i,h the point where the supreme off h is obtained, for a given value of h satisfying b n ≤ h ≤ a n in the interval (x i − δ, x i + δ) and then combining (17) and (19) Finally, using thatf h (n) has j modes and x i was defined as the closest mode off h (n) to Results (R1) and (R2) are derived from (20) and (21).

E Numerical approximations
Details of the numerical approaches used in this paper can be found bellow. All the functions were implemented in the statistical software R Core Team (2016).

Practical computation of the critical bandwidth
To obtain both the critical bandwidth of Silverman (1981) and that one of Hall and York (2001), a binary search procedure was used. In each step of the algorithm, denoting h a as the bandwidth in such a way than f ha has at most k modes and h b the bandwidth for which the kernel density estimation has more than k modes (both in the interval I if the critical bandwidth of Hall and York (2001) is being calculated). Then the dichotomy algorithm is stopped when (h b − h a ) < (h a 0 /2 10 ), where h a 0 is the initial value of h a . The last calculated value of h a is the one employed as the critical bandwidth.

Practical computation of the excess mass
To obtain the excess mass defined by Müller and Sawitzki (1991) when the null hypothesis of unimodality is being tested, the following result is used: the value of the excess mass statistic is exactly twice the value of the dip statistic introduced by Hartigan and Hartigan (1985). The value of the dip was obtained using the diptest package implemented by Maechler (2015).
For the general case, the following algorithm will be employed to obtain the excess mass statistic when the null hypothesis H 0 : j = k, with k > 1, is tested. First, assume that d k (p) is the minimum distance of the union of k intervals containing p data points. To get the possible values of λ corresponding to E n,k (P n , λ), and candidates to minimize the difference D n,k+1 (λ) = {E n,k+1 (P n , λ) − E n,k (P n , λ)}, the first step begins in q k (1) = n (where the number in parenthesis is the current iteration).
Then, it search from q k (1) − 1 until (k + 1) the integer q k (2) minimizing the following expression λ k (1) = min q k (2) q k (1) − q k (2) n(d k (q k (1)) − d k (q k (2))) , the value of λ k (1) is one of the possible values of λ minimizing D n,k+1 (λ). Then, if q k (2) = (k + 1) is the value minimizing the previous expression, the algorithm will be stopped, otherwise it is continued until q k (t k ) = (k+1) (where t k is the realised number of iterations).
After obtaining the vector λ k of possible values of λ minimizing D n,k+1 (λ), the algorithm is repeated for (k + 1) to get the vector λ k+1 . Then the excess mass statistic is easy to obtain using that ∆ n,k+1 = min λ k ∪λ k+1 D n,k+1 (λ).
To get the values of d k (p), one can obtain the exact result of the excess mass employing the algorithm provided by Müller and Sawitzki (1991). A similar algorithm was implemented and employed to get the exact results in Section 4. The problem of employing this algorithm is the high computing cost for an extensive study. For this reason an approximation was employed in Section 3 to get the excess mass statistic, ∆ n,3 . The new algorithm consist in, first, calculate d 1 (p) to obtain the exact values of λ 1 and secondly create a grid of l possible values between each λ 1 (j) and λ 1 (j + 1), with j and entire value between 1 and (t k − 1). Finally, to get the test statistic, that value of λ belonging the entire grid and minimising D n,3 (λ) is chosen. The employed size l was: l = 100 when n = 50, l = 40 for n = 100, l = 20 if n = 200 and l = 5 when n = 1000. This selection of points represents a balance between the accuracy and the computation time, as, in general, if n is large then the length of the vector t k is also large.