{'fraction_non_alphanumeric': 0.02955339439146445, 'fraction_numerical': 0.01342649419318289, 'mean_word_length': 4.513378448724622, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Application of optimal data-based binning method to spatial analysis of ecological datasets

Investigation of highly structured data sets to unveil statistical regularities is of major importance in complex system research. The first step is to choose the scale at which to observe the process, the most informative scale being the one that includes the important features while disregarding noisy details in the data. In the investigation of spatial patterns, the optimal scale defines the optimal bin size of the histogram in which to visualize the empirical density of the pattern. In this paper we investigate a method proposed recently by K.~H.~Knuth to find the optimal bin size of an histogram as a tool for statistical analysis of spatial point processes. We test it through numerical simulations on various spatial processes which are of interest in ecology. We show that Knuth optimal bin size rule reducing noisy fluctuations performs better than standard kernel methods to infer the intensity of the underlying process. Moreover it can be used to highlight relevant spatial characteristics of the underlying distribution such as space anisotropy and clusterization. We apply these findings to analyse cluster-like structures in plants' arrangement of Barro Colorado Island rainforest.


Introduction
Nowadays, a huge quantity of data structured on different time and space scales are easily available. Analysis of these massive databases reveals that despite their diversity and complexity, natural phenomena are characterized by the emergence of regularities that are largely independent of biological and physiological details. One is the tendency, observed both in ecological communities and in human activities, to form spatial or temporal clusters (He et al. (1997); Condit et al. (2000); Plotkin et al. (2000); Adorisio et al. (2014)). However, classifying a spatial point pattern as clustered rather than regular can be a challenging task because establishing the main features of its spatial density function strongly depends on the scale through which we look at it.
More generally, it is well known that the form of a data-based density function may depend on the algorithm (binning rule) used for the binning of the data (Etienne and Haegeman (2010)). In our view the main flaw of many binning rules (Sturges (1926); Yule and Kendall (1950); Doane (1976); Freedman and Diaconis (1981); Stone (1984); Scott (2015)) is that they assume some knowledge on the data distribution. For example Sturges rule (Sturges (1926)) assumes that the data are normally distributed. This is a key point if you have in view applications to ecological datasets. In many cases it is not reasonable to assume such knowledge and the process generating the dataset must be considered unknown. Therefore any criteria based on some prior knowledge of the true density should not be applied as it often introduces a degree of arbitrariness that may produce biased conclusions.
In this paper we are concerned with the statistical analysis of spatial patterns describing the location of plants in a tropical forest study area. We intend to use a method proposed by K. H. Knuth (Knuth (2013)) to find the optimal bin size of a two-dimensional histogram. Knuth non parametric method selects the optimal scale from the data without any assumption on the underlying process that generated the data. We show that the Knuth method can be used to highlight relevant spatial characteristics on the underlying distri-bution such as space anisotropy and clusterization. We tested it against most currently used (Epanechnikov) kernel method for two-dimensional datasets and one-dimensional (Stone binning rule) and it results to be more efficient in detecting Complete Spatial Random processes and avoiding sample fluctuations.
Therefore our analysis validates it as a reliable method for determining the intensity function of a pattern. Additionally, it is not subject to the virtual aggregation phenomenon (see below Sec. 3.3). It correctly detects homogeneity or the presence of a gradient in the density function and the relative difference of the rectangular bin sides is a measure of the anisotropy of the pattern. It also allows to infer quantitative (cluster size) information on both first and secondorder statistics. Thus, it is not only a rule to choose the bin size in which to organize the data. Indeed our analysis proves that Knuth bin size is a good indicator of how finely structured is the dataset and that it can be used as a trusted tool for the preliminary statistical analysis of a spatial dataset. We show which are the relevant information contained in the size and the shape of the optimal bin and how they are related to the spatial features of the process/dataset. We test our findings on an ecological dataset consisting of the spatial coordinates of individuals belonging to 300 different species of plants located in a 50 ha rectangle of the Barro Colorado Island rainforest (BCI). In particular, we study cluster-like structures in plants' arrangement.
The present analysis should help inform future investigations of temporal or spatial features of different complex systems in ecology and human dynamics (Simini et al. (2012); Formentin et al. (2014); Sanli and Lambiotte (2015)).

On optimal binning rules
There exist diverse rules to determine the optimal number of bins of a histogram. Some of the most known (Scott (1979); Freedman and Diaconis (1981); Stone (1984); Scott (2015)) rely on the minimization of the L 2 norm between the histogram and the true underlying density on which they assume some prior knowledge. Assuming such prior information is not reasonable for ecolog-ical datasets and the process generating the data must be considered unknown.
Hence, any criteria determining the optimal bin size based on prior knowledge on the true density should not be applied. Moreover, some methods work well for unimodal densities while they are known to be suboptimal for multimodal ones.
In particular, Freedman and Diaconis rule (Freedman and Diaconis (1981)) is not valid for uniform or piece-wise constant density functions whereas Sturges rule (Sturges (1926)) is not suitable when the data exhibit skewness or any other non-normality. To overcome these difficulties, K. H. Knuth (Knuth (2013)) proposed a method based on Maximum a posteriori estimation and Bayes' Theorem. The main advantage is that no prior information about the density from which the data are sampled is assumed. Moreover, the Knuth method is more efficient in avoiding sample fluctuations and it best approximates the probability density in the L p norm sense, with p = 1, 2. See in Supporting Information the comparison between the Knuth method and the one developed by C. J. Stone (Stone (1984)) and in Sec. 2 below the comparison with standard kernel methods.

The Knuth method
We briefly present the Knuth method referring to Knuth (2013) for further details. Let us suppose to be working with a set of N 2-dimensional data sampled from an unknown probability density function we wish to estimate. In applications to spatial analysis of ecological datasets each 2-dimensional datum may represent the location of a tree.
Let us then cover our data span V with a rectangular grid of M x ×M y rectangular bins, where M x is the number of bins along the x-axis and M y the number along the y-one. We assume no measurement uncertainty about the data hence V is precisely known. We set M = (M x , M y ) and M = M x · M y . We denote with a x and a y the width of each bin along the x and y axis respectively and with a = a x · a y the bin area. We call h k , where k = 1, ..., M is the bin label, the correspondent histogram column's height. After the normalisation of the volume of the histogram, h k represents the constant value of the probability density function over the region of the bin. The volume of each histogram column π k = h k a is the probability mass of each bin, i.e. the probability of finding a datum in the range dictated by the k th bin. Therefore the piece-wise constant density function of the histogram with M bins is: where Π k (x, y) is the boxcar function, defined as Notice that due to the normalisation, only M − 1 probabilities masses are independent.
The likelihood function. When we arrange the data into a histogram, the probability that a datum falls within the k th bin is given by the probability mass of that bin. The associated probability density h k is called likelihood function.
Denoting with π = (π 1 , ..., π M −1 ) the independent probabilities masses, with d = (d 1 , ..., d N ) the vector of 2-dimensional data points, with n k the number of these latter contained in the k th bin and assuming that sampled data are independent, their joint likelihood reduces to the product of N factors Prior and Posterior Probability. The prior probability of the number of bins represents our knowledge of it a priori. As we do not know anything but the total range V of the data, it is reasonable to set where C is the maximum number of bins we wish to consider. Similarly, we take as the prior probability of the bin masses π 1 , ..., π M −1 the uniform probability on the simplex defined by the corners of an (M −1)-dimensional hypercube with unit side lengths due to the normalisation condition: where Γ(·) is the Gamma function and 0 ≤ π i ≤ 1, i = 1, ..., M − 1. Eq. (5) is called Jeffreys's non-informative prior (Box and Tiao (1992)) and it expresses complete ignorance about the form of the histogram.
To infer the posterior probability density for the number of bins M we use Bayes' Theorem, which leads to Inserting the expressions (3), (4) and (5) into the equation (6) above, we get the joint posterior probability for the bin parameters π and pair M : Notice that we disregarded p(M ) since it is constant.
By integrating equation 7 over all admissible values of π 1 , ..., π M −1 we get the posterior probability for the number of bins M : According to the method of Maximum a-posteriori Estimation, the optimal choice for M , the one that provides the best agreement of the model with the observed data, is the one which maximizes the logarithm of the posterior probability. Thus, the optimal number of bins iŝ A great advantage of this mathematical formalism, is that, once we have computed the optimal binning grid M , from (7) and (8) (see Knuth (2013) for details) we can analytically compute the mean and the variance of the bin probabilities which are given, respectively, by and and which allow to construct the optimal histogram with the proper error bars.

Knuth method for estimation of a Point Process's intensity
Given a point process φ defined on the plane, its intensity function λ or density is a first-order statistic measuring the mean number of points per unit area. Denoting with N φ (W ) the number of points falling within a region W ∈ R 2 , we have We are interested in estimating, from a sample of the process, its density λ.
Notice that if the process is assumed to be homogeneous, then its density is constant within the observation window W and it is therefore approximated by where A W is the area of the window.
Using the Knuth method, we can approximate the intensity of the possibly non homogeneous processes by a piece-wise constant function: the density in a point x is the height of the column over the bin containing x in the optimal histogram.
We compare Knuth answer with kernel methods which are widely used in literature (Illian et al. (2008); Shiffers et al. (2008); Wiegand and Moloney (2013)) (see Supporting Information for a comparison with Stone's non kernel method).
The idea of kernel methods is to estimate λ(x), x ∈ R 2 by looking at the number of points falling within the small disk B(x, R) of radius R (called bandwidth) centred in x and by dividing it by the area of the disk: where x i is a point of the process φ and k R (|x − x i |) is the kernel function. A popular choice for k R is the Epanechnikov kernel, defined as: This function weights the points x i ∈ φ within the disk B(x, R) according to their distance from the centre x: the smaller |x − x i |, the bigger the weight.
A weak point of kernel methods is that a general recipe for the choice of the bandwidth R, which gives the right smoothing of a rugged intensity function, does not exists (Illian et al. (2008) p. 115, Diggle (2013 p. 116). Some authors suggest the rough estimate R ∼ 1/ √ λ and a subsequent finer tuning of R using visual inspection (Wiegand and Moloney (2013), p. 97). On the contrary, the Knuth non parametric method select the optimal scale from the data without any assumption on the underlying process that generated the data. Notice that also simpler methods for assessing the CRS hypothesis like the quadrats count can be considered as parametric methods since their effectiveness depends on the choice of the size of the quadrats.
Below we compare kernel and the Knuth methods on datasets generated from both homogeneous and inhomogeneous processes to test their reliability at 1) reproducing the intensity of the process and 2) capturing the eventual presence of a gradient in the density functions.
Detection of CSR and inhomogeneity. We start applying the Knuth method and the Epanechnikov kernel estimation to a generated CSR pattern within a window W of area 500 × 500 units. On the left of Figure 1(a) we plot the CSR (Poisson) pattern generated with λ = 1/500, while in the two graphs on the right we can see the intensity function estimated by the Epanechnikov kernel and by the Knuth method. This latter arranges the data in a unique bin, so that the intensity function it returns is constant within the plot. Therefore it perfectly detects the underlying homogeneous structure of the process. By contrast, the kernel method results to be sensitive to sampling fluctuations.
Robustness of the Knuth method in detecting CSR processes is tested in Sup-porting Information. We then compare the two methods on an inhomogeneous Poisson point process where the true density function λ increases with the y coordinates (see Figure 1(b)). Although from the Epanechnikov kernel method it is evident that the number of points falling within the upper region of the window is bigger with respect to the lower region, it results to be still affected by sample fluctuations. By contrast, the Knuth method arranges the data into a 1 × 4 grid, perfectly detecting the homogeneity of the process along the x-axis and the density gradient along the y-axis.

Knuth method description of cluster features
As a preliminary analysis we investigated how the Knuth method reproduces the features of three different type of clusters: square, circular with constant density and circular with Gaussian density. We consider plots of area A 0 = 1000 × 500 units as is the BCI. For each type of cluster, we have generated 100 datasets setting the number of individual equal to N = 1000 and arranging them in a unique cluster positioned in the centre of the window. The characteristic size σ of the cluster (respectively the side, the radius and the standard deviation) varied from 1 to 100. For each dataset, we have computed Knuth optimal binning area a and the correlation between the characteristic size of the cluster σ (x-axis) and its optimal bin size representation a (y-axis). Results are displayed in Figure 2.
In all cases, the determination coefficient shows a strong correlation between the cluster size and the optimal bin area (R 2 > 0.9). In the first case (square clusters) the slope equals 1, meaning that the square cluster is arranged in a unique bin and that the density of points is correctly detected as constant.
In the other two cases, the slope is far from 1, meaning that the cluster is described using a higher number of bins. In particular, in the second case, the Knuth method aims at reproducing the circular boundary of the cluster while in the last case, the optimal bin serves at reproducing both the circular boundary of the cluster and the Gaussian shape of the density (see Supporting Information for the correspondent histograms and density plots). We also tested the Knuth method on uniform and Matérn clusters (see Figure 2 and Supporting Information).

Second-order statistics for estimation of cluster size and hard core radius
In spatial ecology, the CSR hypothesis mostly comes to fail due to several reasons: on one side, changes in environmental conditions, such as in physical features of the landscape or in chemical composition of the soil, may lead to inhomogeneous patterns. On the other side, different seed dispersal mechanisms may favour the formation of clumped structures as well as dispersed one. These pattern characteristics are revealed by second-order statistics, which take into account the correlations between pair of points due to possible interactions.
Below we briefly recall the most common second order statistics and their use.
The probably most used second-order statistics for homogeneous pattern is the Wiegand and Moloney (2013)), which is the expected number of points falling within a distance r from a point of the process, divided by the intensity function λ. Since under CSR hypothesis K(r) = π · r 2 scales quadratically with the distance, it is usually substituted with the L-function L(r) = K(r)/π − r, which takes constant zero value for the CSR model. Both Ripley's K and L functions are regarded in literature as cumulative statistics (Wiegand and Moloney (2013)), which do not permit to properly infer pattern characteristics at any specific scale.
This limitation can be avoided by taking the pair correlation function g(r) = K (r)/2πr (Adorisio et al. (2014); Azaele et al. (2015)). It is defined as the ratio between the density of points falling within a small ring of radius r centred at a point of the pattern and the constant density function λ of the process assumed to be homogeneous.
Since values of g(r) greater than 1 indicate clustering, whereas values less than 1 indicate dispersion, the point where the pair correlation function intercepts the line y = 1 gives a rough estimate of the average diameter of clusters for clumped patterns or the average hard core radius for overdispersed ones.
As pointed out in (Shiffers et al. (2008); Wiegand and Moloney (2013)), both Ripley's K and pair correlation function are window-dependent, in the sense that their ability in detecting the scale of clustering or dispersion depends on the window within which we arrange the data. For example, if the plot contains empty spaces, these statistics put in evidence a clustered structure which is not due to the pattern's generating process. This phenomenon is called virtual aggregation and it affects how these statistics perform in inhomogeneous pattern.
In this latter case, the apparent clumping nature of the pattern is due to the fact that both K and g are normalised with the empirical intensity λ, which is wrongly assumed as constant. Shiffers' K2-index (Shiffers et al. (2008)), using the derivative of the pair correlation function, is less sensitive to this phenomenon.
A main concern of this paper is to test the ability of the Knuth method to deal with clustered, dispersed and inhomogeneous patterns and to cope with the virtual aggregation phenomena which may arise when considering non homogeneous patterns. We have thus selected a model of clustered process, the modified Thomas process and a dispersed one with an hard core repulsion radius. They are briefly described below.
The Poisson Cluster Process or, more often, its simplest variant the modified Thomas process (mTp) (He et al. (1997) (2013)) is more efficient than others to capture important biological curves such as the Species-area curve (see, e. g. Plotkin et al. (2000)) or to model species occupancies at different spatial scales (Azaele et al. (2012)). Instead, as shown in Morlon et al. (2008), it inadequately reproduces the distance−decay relationship, indicating that some of the assumptions of mTp do not hold in nature.
The hard core process is generated from a uniform distribution with the additional constraint that if a point comes to fall within a fixed hard core distance from a pre-existing one, it is rejected. This model is of importance in ecology to describe reproductive mechanism where the seeds are shot apart from parents to avoid species-specific predators (Janzen-Connell effect, see Janzen (1970)

Interplay between second-order statistics and the Knuth method
Despite their ability to detect significant departures from CSR processes, the pair correlation function and the K2 -index lack in reliability at determining the cluster size when dealing with inhomogeneous patterns being subject to the virtual aggregation phenomenon.
To show this, we computed the pair correlation function and the K2 -index of Shiffers for three generated patterns, respectively: a) a modified Thomas process, b) an overdispersed pattern with fixed hard core radius and c) an inhomogeneous Poisson process.
To generate the virtual aggregation, they are firstly considered within a 500 × 500 units window (black curves of Figure 4), and then within a larger one, obtained from the previous by adding at its bottom an empty square box of the same area (grey curves). In Figure 4(a) we generate a point pattern according to a mTp with parameters ρ = 2 · 10 −4 , σ = 10, µ = 10, where ρ is the intensity of the Poisson process from which parents are generated, σ is the standard deviation of the Gaussian distribution of offspring around each parent and µ is the mean number of offspring per parent. The two second-order statistics g and K2 well capture the clustered structure of the pattern, either in the original window and in the expanded one. In particular, they reveal an average clump diameter around 37-40 units, which overestimate the true value 2σ · π/2 ≈ 25 units.
Looking at the grey curves, we can observe the phenomenon of the virtual aggregation. While the black pair correlation function intercepts the line y = 1, the grey one sees the species as clustered at any scale. By contrast, their corresponding K2 -indexes are much closer one another, meaning the addition of an empty box slightly affects this last statistic. By applying the Knuth method we get a 22 × 20 grid which results in a clump size of 26 units circa, which is therefore the closest to the real one. Notice that the Knuth method automatically restricts to the data span V (see Sec. 1) therefore it results to be insensitive to voids in the pattern and hence does not suffer from the virtual aggregation effect.
In Figure 4(b) we carry out the analysis for 500 points sampled from a uniform distribution with the additional constraint that if a point comes to fall within a fixed hard core distance from a pre-existing one, it is rejected. Once again all statistics are able to capture the overdispersion at small scale: they return an hard core radius around the true value of 10 units, with a slight difference between black and grey curves. Applying the Knuth method we obtain a 1 × 1 grid as in the CSR case: the estimated density is therefore correctly detected by the optimal histogram which sees the homogeneity of the pattern. Here we see a limitation of the Knuth method which is unable to reveal second-order information such as the dispersion of the process and the hard core radius.
The last pattern shown in Figure 4(c) is sampled from a Poisson process with intensity λ increasing with the x-coordinates, thus presenting a small gradient along the axis. Here the phenomenon of virtual aggregation is well visible looking at the graph of g: the pattern is detected as clustered at all scales by the statistic. Instead, the K2 index is not affected by this, although it cannot distinguish the process from a CSR one. By contrast, the Knuth method arrange data in a 4 × 1 grid, which permits us to capture the homogeneity along the y-axis and the gradient along the x-axis.
In conclusion, combining the Knuth method to the Shiffers' K2 index leads to a better understanding of the underlying process from which a sample is observed: the former permits to detect homogeneity against gradient in densities and gives a measure of how structured it is in the sense that small bins indicate clumping while large bin indicate Poisson-like distributions. The latter, on the other hand, allows to give quantitative information on the scale of both clumping and dispersion.

Detection of anisotropy
We tested the Knuth method's capability of detecting another relevant characteristic of a process: the anisotropy. By construction, Knuth algorithm is sensitive to the inversion of the orthogonal axes: if the two-dimensional pattern is tilted by a multiple of 90 • , also the resulting optimal grid is rotated of the same angle. In Supplementary Information we applied the Knuth method on a generated dataset aggregated in a unique anisotropic Gaussian cluster distributed according to a bivariate Gaussian with standard deviation along the x−axis twice than along the y−axis; it resulted sensitive to the inversion of variance's direction between 0 • -90 • and 45 • -135 • cases, which leads to the inversion of the optimal bin sizes.
These results confirm that the Knuth method is very efficient at detecting the anisotropy of the cluster structure. However, an obvious limitation of the method is that it returns a unique bin size for the whole plot. Therefore, when dealing with real patterns, where anisotropic clusters may be structured at different scales or oriented in many directions, the optimal bin size is the result of a compromise between these different sub-structures of the dataset. Notice that while different bin sizes denote that the Knuth method classifies the pattern as anisotropic, the contrary does not hold.
However, previous remarks suggest that the difference in bin width along two orthogonal axes is an index of how anisotropic is the spatial pattern. We define the anisotropy index of a pattern as where a x is the bin width along the x axis and a y the bin width along the y axis in order to be invariant with respect to inversion of the axes. To test the usefulness of the anisotropy index as a tool for detecting the anisotropy of a spatial pattern, we computed I an for a set of patterns generated from the mTp process with parameters N = 3000, ρ · A 0 = 5 and σ ∈ 0, ..., 100. The Azaele et al. (2015)), in many cases is not supported by the Knuth method.
Notice that checking the goodness-of-fit of a fitted model using a minimum contrast method when the theoretical form of the summary function (usually Ripley's K) is not known is a difficult task. For a mT p process the form of Ripley's K is known (see Supplementary Information) hence the fitted pattern is already optimal with respect to a minimum contrast goodness-of-fit criterion.
Nevertheless, the analysis of the difference of the optimal bin size for the real and the mT p-generated pattern may reveal a strong departure form the real data (see below).
To find the species that are suitable to be described by a clumped pattern, we have selected all species with abundance between 20 and 3000 individuals (204 species) and we have then computed the mTp parameters (ρ, σ, µ), for each species by fitting the empirical Ripley's K function (minimum contrast estimation, see Diggle (2013), Ch. 7, and Supporting Information). As in Morlon et al. (2008), we have discarded the species whose distribution is not quite distinguishable from a random one -in this case the Poisson cluster process cannot capture correctly the underlying structure of the data-according to the following criteria: i) the mean cluster diameter σ √ 2π is smaller than 500 m and ii) the number of clusters is smaller than the number of the individuals.
There are 183 of the 204 species satisfying both criteria. For each of the 183 species we have computed a/π which is the radius of the circle equivalent to a rectangular Knuth optimal bin of area a. We know this is a measure of the size of the underlying minimal structure of the dataset. The frequency histogram (see Supporting Information) indicates that there is no preferred choice for a common optimal bin area for all species, since the values of the optimal binning areas span from 100 m 2 to 5 · 10 5 m 2 circa. This is not surprising, since each species has its own distribution due to myriad of factors such as seed dispersal, gap recruitment or adaptation to the surrounding soil (Augspurger (1984); Plotkin et al. (2000)).

Difference index for BCI
We computed the Knuth optimal bin size a both for BCI real species with σ 2 π/2 < 10 4 m 2 and for their correspondent mTp generated counterparts and their difference ∆ = a(mT p) − a(real)). In Figure 5 we display the exploratory analysis (histogram, boxplot and QQ plot) of the difference data. We see that the histogram is right-skewed with fat tails. For 50% of the species the difference of the bin area ∆ is bigger than 2 times the smaller of the two. For these species which are in the tails of the distributions, the mT p process fails at reproducing the real pattern either because the generated clustered pattern has a coarser scale than the real one losing details of the original fine structure (a(real) a(mT P ), right tail) or because it introduces an artificial clustered structure on a more uniform real pattern (a(real) a(mT p), left tail).
To have a visual inspection of the difference between real and generated species, we have selected three species respectively i) in the right tail, ii) in the centre and iii) in the left tail. In Figure 6 we display the real and the mT p generated patterns with parameters fitted from the data. In the first case the mT p reproduces a coarser pattern than the real one while in the latter it introduces an artificial finer structure.
In Figure 7 the frequency histogram of the mT p clump area σ 2 · π/2 is superimposed to the one of a/π. In this cumulative plot the two histograms are respectively right and left-skewed showing that globally Knuth method assigns a finer structure to the real patterns with respect to the generated ones.

Anisotropy index for BCI
We computed the anisotropy index |a y − a x |/ max(a x , a y ) both for BCI real species with σ 2 π/2 < 10 4 m 2 and for their correspondent mTp generated counterparts. Figure 8 shows the frequency histogram of anisotropy index for real species superimposed on the one for mT p generated species. In the first case there is a high number of species whose index is far from 0, meaning that their underlying density function is not recognised as isotropic by the Knuth method.
In the second case, Knuth histogram is more shifted against the y−axis, meaning that the Knuth method sees the new generated species' process more isotropic than before. This fact suggests that a reason for which mTp process fails in capturing some important ecological curves is the fact that the hypothesis of isotropic clusters is too strong and is not supported by real data.

Relation with abundance
The relation between species abundance and degree of aggregation is still de- in Morlon et al. (2008), the correlation between these two quantities strongly depends on how they are measured. For example, for Pasoh forest, in (He et al. (1997)) the proposed Donnelly clumping index based on nearest-neighbour distance shows a slightly positive correlation between abundance and aggregation.
By contrast, the relative neighbourhood density Ω 0−10 of Condit et al. (2000) and the Cramer-von Mises-type k statistic of Plotkin et al. (2000) are negatively correlated to abundance.
To investigate if mTp and Knuth aggregation parameters are similarly correlated with species abundance or not, we have computed on the one hand the correlation between the abundance and the following mT p quantities: the mean number of parents ρ · A 0 , the mean clump radius σ · π/2, the mean number of offspring per parent µ and the relative neighbourhood density Ω 0−10 (see Morlon et al. (2008)) and on the other hand the correlation of species' abundance with Knuth optimal bin area and index of anisotropy. These latter, as we have seen, give us information on how structured is the data density function and how far it is from a uniform or isotropic one (see Supporting Information).
From the determination coefficients, only the relative neighbourhood density shows a negative correlation with abundance, while the other mTp parameters result to be slightly positive correlated with it. This is in accordance with the literature (Morlon et al. (2008)). Knuth optimal bin area ( Figure 17) and index of anisotropy are, instead, insignificantly correlated with the abundance with respect to the determination coefficient. This is quite reasonable because Knuth optimal grid depends only on data distribution, and not on their abundance.

Results and Discussion
Knuth optimal binning method, based on Bayes' Theorem and Maximum a-posteriori Estimation allows to infer the least biased estimate of the underlying density function of a point pattern. The optimal bin size sets the most informative scale at which to observe the data.
We showed how to use the Knuth optimal bin size and shape to estimate the intensity of a spatial process and infer characteristic spatial features as anisotropy and clusterization. We tested it against most currently used kernel method for two-dimensional datasets (Epanechnikov) and one-dimensional (Stone binning rule) ones and it resulted to be more efficient in detecting CSR processes and avoiding sample fluctuations in both cases. It is therefore a reliable method for determining the intensity function of a pattern. Moreover it does not need any a priori assumption about the phenomena that generated the data.
Since the Knuth method is based on a maximization procedure it does not contain adjustable parameters and it is not subject to the virtual aggregation phenomenon. When used in conjunction with the K2 -index, it allows to infer qualitative (it correctly detects homogeneity or the presence of a gradient in the density function, even if it is not strong) and quantitative (cluster size) information on both first and second-order statistics.
Our analysis indicates that Knuth bin size is a good indicator of how finely structured is the dataset and the relative difference of the rectangular bin sides is a measure of the anisotropy of the pattern. For the above reasons it can be used as a trusted tool for the preliminary statistical analysis of a spatial dataset.
Finally, we tested our findings on the BCI ecological dataset to have information about distribution of the size of cluster-like structure of plants, anisotropy of plant distribution and existence of uniformly distributed species. We found evidence that the choice of modelling a species' distribution through a modified Thomas process, which has been proven to be efficient in capturing some important biological curves but not others, is not always supported by the Knuth method. Moreover, we found that cluster size, measured by optimal bin area is insignificantly correlated with the abundance of a species. Globally, it provides a reliable method to test whether an hypothesis made on the underlying process of a pattern is justified or not. We are confident that this survey of Knuth algorithm' performance can be of help for the scientific community. consisted of 1000 points and for each case σ varies from 1 to 100. In the first case we plot a against the square of the cluster side σ 2 , getting a determination coefficient of R 2 = 1.00. In the second case we plot a/π against the squared radius of the cluster σ 2 and we get R 2 = 0.90 and in the last case, plotting a/π against the square of the mean distance of a point from the cluster centre σ 2 · π/2, we get R 2 = 0.94. The correlation between the two variables was very good in all cases.

Application of Knuth method on one-dimensional datasets
As in (Knuth, 2013) we test Knuth's method with four sets of data, each consisting of 1000 one-dimensional points sampled from known probability density functions. In order to do this we resort to M AT LAB software and Knuth's OPTBINS binning package v1.0. The results are shown in Figure 10.
For each set of data we plot the logarithm of the relative posterior given by eq.(8) of the main article and we then arrange data in an histogram with the found optimal bin numberM . It represents the piece-wise constant density function of our data we wished to obtain. To highlight the fact that the method is able to avoid sampling fluctuations while it captures the main characteristics of the underlying density function, we plot, under each optimal histogram, the one we would get using 100 bins instead ofM , which therefore better shows the spatial distributions of the data points.
Let us notice that the optimal binning numbers that we obtain in the four cases depend on the particular realisation of the process from which the data we are analysing are sampled. Anyway, as we will show later for CSR processes e. g., Knuth's method resulted to be very stable, so that the results it gives slightly differ from a sample to the other.
The first case we consider is the sampling of 1000 data points from a uniform density. From Figure 10(a), we see that the relative log posterior reaches its maximum atM = 1. Knuth's histogram therefore coincides perfectly with the underlying density function.
In the second example, the set of data is sampled from a four-steps density function. Again we can obtain it by using Knuth's algorithm, whose relative log posterior peaks at four bins.
In Figures 10(e) and (f) we show the results we obtain applying the method to a set of data sampled from a standard Gaussian density function N (0, 1). As we can see from Figure 10(e) the relative log posterior reaches its maximum at 12 bins, which capture very well the structure of the underlying density function.
Finally, in the last example, we sample 1000 data points from a density function consisting of three Gaussian peaks over a uniform background. In this case the relative log posterior has a maximum forM = 31. Looking at Figure 10(g) we can notice that the algorithm gives an optimal binning number much higher than the ones obtained in the previous examples. This is due to the fact that, since we are imposing that the bins have equal width v, we need an high number of bins in order to capture the real structure of the underlying density function and to identify the three peaks, in spite of the uniform background, usually arranged within only one bin, as we have seen in the first example.
This fact suggests that the number of bins given by the method could be a good parameter to distinguish clustered data from uniform ones, since the more structured is a species, the higher the number of bins we need to capture it correctly.
Moreover, we can notice that the relative log posterior decreases much slower than the previous ones, so that each value M ≥M could be taken as the number of bins of our histogram in order to faithfully approximate the underlying density function. This is a general feature that we can observe when we are dealing with data distributions characterised by the presence of clusters, such as the three Gaussian peaks of the last example.

Comparison with Stone
We now introduce a non-kernel method to select the optimal number of bins developed by Charles J. Stone in 1983(Stone, 1984. We will then compare the two methods to show how Knuth one is more efficient in avoiding sample fluctuations. Stone's optimal selection rule. Let us assume we have a set of N data (d i ) ∈ R n×N , sampled from a known probability density function p.
Let us denote with a = (a 1 , ..., a n ) ∈ R n the point in the space where we wish to start building our histogram and with b = (b 1 , ..., b n ) ∈ R n + the dimension of the n-dimensional bin we set. Our histogram will therefore depend on the pair (a, b). We wish now to compute the value of the piece-wise constant probability density of the histogram in the k th bin, whose integer coordinates, starting from a, are k = (k 1 , ...k n ) ∈ Z n and which therefore occupies, in R n , the place × (a n + (k n − 1)b n , a n + k n b n ).
We have obtained a partition of R n into identical bins of dimension v (a,b) = n i=1 b i . The probability mass of the k th bin, i.e. the volume of the column over it, is given by the empirical distribution From this we can compute the height of the column, which is simply The probability density function of the histogram is therefore given by where χ I (a,b),k is the characteristic function of I (a,b),k , defined as Similarly we can define the probability mass as follows: In order to best approximate the density function from which the data were sampled, Stone's idea relies on minimising the integrated squared error of h (a,b) Under some conditions, (see (Stone, 1984)), minimising (21) is asymptotically equivalent to minimising Comparison to one-dimensional datasets. Let us firstly remark that the condition above is satisfied, for example, if there is some non-empty open subset of R on which the derivative of p exists, it is continuous and non-zero (Stone, 1984).
As a consequence, if we are dealing with uniform or step density functions, such as we did in the previous section with Knuth, we cannot apply Stone's optimal selection rule, which has therefore a smaller range of applications. This is a drawback especially when working with a real dataset sampled from an unknown density function, since we cannot exclude a priori such densities.
We have compared the two methods in the one-dimensional case using L 1 and L 2 distance method (Gomes-Gonçalves et al., 2014).
We have generated 50 datasets each consisting of 1000 one-dimensional points sampled from a standard Gaussian density function N (0, 1). For each test, we have computed Knuth and Stone optimal binning number and then we have constructed the correspondent histograms. We finally have compared the L 1 and L 2 distances between these and the underlying density function.
In Figure 11(e) we plot the ratio between Stone's and Knuth's L 1 distance, while in Figure 11(f) we plot the logarithm of the ratio between Stone's and Knuth's As we can see in Figure 11(e), while the two methods are practically equivalent in the first case, Knuth's L 2 distance is about ten times smaller than Stone's one ( Figure 11(f)), due to the fact that the latter is more sensitive to sample fluctuations Knuth is able to avoid instead. The result may be surprising, since Stone's optimal selection rule relies on minimizing exactly the L 2 distance between the underlying density function and the one obtained by the histogram.
The matter is that the rule only assures that this minimum is reached when the number of data tends to infinity. In fact, it is in this asymptotic case that the sampling fluctuations diminish and thus the histogram recalls quite faithfully the distribution from which the data are sampled.
For the sake of completeness, we also insert the graphics we obtain testing the two methods on one dataset such as the ones described above (a Gaussian).
From Figure 11(d) and (c), it is clear that while Knuth's histogram correctly captures the main characteristics of the underlying probability density function, Stone is more sensitive to sampling fluctuations, which results in a higher optimal binning number.

Knuth sensitivity to uniform distribution
We firstly test how sensitive and robust is the algorithm with respect to the dataset in order to see if the answer that it gives depends from it rather than from the distribution from which it has been sampled.
With this goal in mind, we test Knuth algorithm's stability by seeing if, generating more times a dataset from a known probability density function, the result is always the same. In particular, we generate a dataset consisting of 1000 points within a 500 × 500 units sampled from a uniform density function. We therefore computed, for each generated dataset, Knuth optimal binning number.
From Section 5, we know that for a one-dimensional set of data uniformly distributed, the answer Knuth's algorithm gives isM = 1, which correctly captures the underlying distribution. We wish to see if the answer for the two-dimensional datasets is the same for almost every test, which would imply that the algorithm is stable with respect to the density function from which we sample the data.
We arrange our results in a 3 × 3 histogram, since the answers given by Knuth algorithm lie in this range.
From Figure 12 we can see that for almost all of the 200 tests we performed the optimal binning number results to be M = [1 1], so that the probability density function we get thanks to Knuth's histogram coincides perfectly with the one from which the data were sampled. This fact provides that the optimal binning number that we get with Knuth's algorithm for each species strictly depends on their distribution rather than on the data itself.

Poisson Cluster Process
Let us recall the definition of the modified Thomas process: • Parents are distributed according to a Poisson process with intensity ρ.
• To each parent a random number of offspring is assigned, drawn from a Poisson distribution of intensity µ.
• For each parent, its offspring is located according to a two-dimensional Gaussian distribution centred at the location of the parent in the plot and with standard deviation σ.
• The parents are removed from the plot, so that the clumped pattern is formed by offspring only.
Let us notice that with this model we are assuming that the formation of clusters is due to an isotropic local propagation of offspring from the parent, which is clearly an oversimplification of the complex natural mechanisms which actually determines their rise.
In order to compute the process parameters (ρ, σ, µ) from a dataset, we computed the empirical K-Ripley function, which, for a spatial Poisson point process with intensity λ, is defined by λ number of extra events within a distance d from an arbitrary event .
In our case, the arbitrary event is the location of any individual belonging to the species. If this latter is randomly distributed, then, since the mean number of stems within a circle of radius d is λπd 2 , we get that for a random placement Let us therefore consider a species with n individuals whose location in a plot of area A 0 are s 1 , ..., s n and let us define witĥ the estimator of the intensity λ for the species. Let us then denote with w(s i , s j ) the proportion of the circumference of the circle with centre s i passing through s j which lies in the plot. Lastly, let us define the indicator function where || · || indicates the euclidean distance between the individuals located in s i and s j respectively. The canonical edge-corrected estimator of Ripley's K function is then (Plotkin et al., 2000) K Defining with d max the larger distance at which we measure a cluster within the area A 0 of our rainforest, we evaluate the edge-corrected estimator between 0 to d max .
In order to get the empirical pairs (ρ, σ) for a generated dataset, we need to compute the Ripley's K function for the modified Thomas process with parameters (ρ, σ), which is given by (Cressie, 1993) We therefore have to choose the pair (ρ, σ) such that the above function best fits the empirical valuesK(0), ...,K(d max ), where we set d max = 300 m as in (Morlon et al., 2008).
We do this by the Method of Minimum Contrast (Diggle, 2013) Once obtained the parameters (ρ, σ) from the dataset, we generated a Poisson cluster model by the following four steps: • We simulate the Poisson cluster process by placing ρ × A 0 + 1/2 parents randomly distributed within the plot.
• We randomly assign each of the n individuals of the species to one of the previously generated parents.
• For each parent, we locate the associated stems according to a two-dimensional Gaussian distribution centred at the location of the parent in the plot and with variance σ. If the case that the offspring falls out of the plot, we impose toroidal boundary conditions.
• We remove the parents from the plot, so that only the offspring remains within it.

Knuth answer for different clustered structures
In the main article we found a strong correlation between Knuth optimal bin area and the characteristic size of three different type of clusters consisting of 1000 points: a square and a circular one with constant density and a circular one with Gaussian density. In Figure 13 we plot an example of three of these datasets with σ = 50.
We compare here Knuth method with Epanechnikov kernel method in estimating the intensity function of these different cluster-structures.
In the first case Knuth collects the points in a unique cluster, as we can see from the optimal grid we insert in the plot of the data distribution. Thus the resulting estimation of the intensity perfectly coincides with the underlying one, given by the constant value λ = 1000/50 2 = 0.4 in the 50 × 50 cluster area.
In the second and in the third case Knuth arranges data in a histogram with a higher number of bins in order to capture the circular boundary of the clusters.
In all cases, we remark that the bin sides are practically equal, meaning Knuth detects the isotropy of the cluster structures.
By contrast, looking at the density plots obtained by Epanechnikov kernel method, we can see that the three clusters are not quite distinguishable one from the other if not slightly for their sizes .

Anisotropic Gaussian clusters
We test the reliability of Knuth method on detecting the possible anisotropy of a point process. In particular, we apply the Knuth method on a dataset aggregated in a unique cluster distributed according to a bivariate Gaussian with standard deviation along the x−axis twice than along the y−axis. We then rotate the data points of different angles around the parent's coordinates and we apply the Knuth method to the new dataset. In Figure 14

Anisotropic uniform clusters
We test now on data generated form a Poisson cluster process where the offspring is distributed around its parent according to a uniform distribution on clusters of dimension 100 × 50 (Figure 15(a)), 50 × 100 (Figure 15(c)) and 50 × 50 (Figure 15(e)).
In Figure 15 we inserted the point distribution in the whole observation window (left column).
As we can see, once again Knuth well captures the uniform structure of the data: it returns bin sizes proportional to the length of the chosen dimensions of the clusters: 50 × 25 for the first case, 49 × 53 for the second and 25 × 26 for the last one.
On the right column we plot the data generated according to an mT p with the parameters that we get by fitting the original data. We get four clusters with σ = 4 for the first two cases, which therefore are considered equal in the mTp fitting and two clusters with σ = 18 for the last.

Application of Knuth to BCI
In (a)-(c) of Figure 16 we show the frequency histograms of the mTp parameters for BCI selected species (see the main article).
We have computed the optimal bin number with Knuth method for each species to investigate how similar or different are the species' distributions of the rainforest within the 50 ha plot and to compare the results with the previous ones.
In Figure 16(d), we see the frequency histogram of a/π, where a is the optimal bin area. This corresponds to the radius of the circle equivalent to the rectangular bin we get by Knuth and which therefore gives us measure of the size of the underlying minimal structure of the dataset.
As we can see from the histogram, there is no preferred choice for the common number of bins of all rainforest's species, since the values of the optimal bin area are quite different from each other spanning from 100 to 5 · 10 5 circa. This is not surprising, since each species has its own distribution due to myriad of factors such as seed dispersal, gap recruitment or adaptation to the surrounding soil [Augspurger, 1984;Plotkin et al., 2000].
In 17 we insert the plot the mT p parameters and Knuth ones against the abundance of species to see if they are linearly correlated. In agreement with the literature (Morlon et al., 2008), we find a slightly positive correlation for the mean clump radius, number of clumps and abundance per clump and a negative one for the relative neighbourhood density. Knuth optimal clump radius and index of anisotropy resulted completely uncorrelated with the abundance.  . On the right column we arrange the data in an histogram with that number of bins. To highlight the fact that the method is able to avoid sampling fluctuations we plot, under each optimal histogram, the one we would get using 100 bins instead ofM , which therefore better shows the spatial distributions of the data sets.