Statistics and probability are two separate branches of mathematics useful to analyze the relative frequency of events. Statistics involves the description and analysis of the frequency of past events attempting to make sense of observations in the real world, whereas probability deals with predictions concerning the likelihood of future events through the examination of the consequences of mathematical definitions that are issued independently of the real world (Skiena, 2003). When applied to spatial data (i.e., data that can be drawn in a map) statistics includes the generation of summaries describing spatial patterns and can be extended to comparisons between those summaries and expectations raised by theories of how the identified spatial patterns formed (Ripley, 1981). Any attempt to forecasting the spatial location and timing of future events should be considered part of a probabilistic approach. If the objects of interest are volcanic vents, studies involving their spatial distribution can therefore have either a statistical or a probabilistic orientation, depending on the emphasis made on the forecasting component.

There is a vast literature discussing several probability aspects that may produce variable and mutually inconsistent hazard and risk estimates in both seismic and volcanic contexts (e.g., Bernreuter et al., 1988; Marzocchi et al., 2008; Neri et al., 2008; Marzocchi and Bebbington, 2012; Marzocchi and Jordan, 2014; Ake et al., 2018; Bevilacqua et al., 2018; Selva et al., 2019; Marzocchi et al., 2021). Works with a more pronounced statistical approach, however, are not as common (Richardson et al., 2012; Cañón-Tapia, 2014; Cañón-Tapia and Mendoza-Borunda, 2014; Delcamp et al., 2019; Jacobo-Bojórquez and Cañón-Tapia, 2020; Cañón-Tapia, 2021b). It must be remarked that the difference made here between the probabilistic and statistic reports is based on the purpose of the study rather than on the tools used for the presentation of data. This might seem confusing at first because some studies use what would seem to be a statistical tool to assess the hazard of a future eruption (Connor, 1987; Connor and Hill, 1995; Condit and Connor, 1996; Conway et al., 1998; Weller et al., 2006; Jaquet et al., 2008; Connor and Connor, 2009; Kiyosugi et al., 2009; Bebbington and Cronin, 2011; Bebbington, 2015; Connor et al., 2019). Nevertheless, the use of a statistical description within a probabilistic framework is not uncommon, and sometimes is referred by the name of “statistical inference” (DeGroot and Schervish, 2012).

Distinction between the probabilistic and statistic approaches outlined above is important for several reasons, particularly when examining the literature aimed to characterize the spatial distribution of volcanic vents. Motivations for the study of vent distribution in volcanic fields aim to increase our understanding of a wide variety of aspects. To mention but a few those aspects include the relationship existing between polygenetic and monogenetic edifices, geochemical complexities that might arise within a field, the distinction between volcanic events and volcanic edifices, the role played by temporal gaps on volcanic activity or even issues related with the very definition of what is a volcanic field. Readers interested in those particular issues are referred to the works by Valentine and Gregg (2008), Szakács and Cañón-Tapia (2010), Kereszturi and Németh, 2013, Németh and Kereszturi, 2015, Cañón-Tapia (2016), and references therein. Underlying all of those issues at the most general level, the statistical approach aims to present a description of the data (e.g., already existing and identifiable volcanic vents) in such a form that some hypothesis concerning the physical structure leading to the observed distribution can be formulated. In contrast, the probabilistic approach focuses on the identification of the site and time with a larger likelihood for the occurrence of a future event. At a more subtle level, the assumptions made by either a statistically- or a probabilistically-oriented interpretation are also contrastingly different. The statistical interpretation should include assessment of the physical conditions favoring the formation of two or more vents during a single eruption, the orientation of tabular conduits that transport magma to the surface and their relationship with the prevailing stress orientation at the time, the possible shift of location of magma reservoirs, etc. A recent review of these topics, and relevant references are provided by Cañón-Tapia (2021a) and Rivalta et al. (2019). In contrast, a probabilistic or statistical-inference interpretation includes assumptions concerning the type of underlying distribution (Uniform, Clustered, Poisson, etc.), but most importantly it makes assumptions about the equality of likelihood of each outcome (past, present and future) and the continuity of the preferred mathematical distribution both in space and time to validate the interpretations of analysis with respect to a single underlying probabilistic model. Identification of the most basic underlying assumptions in many works, however, are not explicitly discussed. Consequently, it is not surprising that although much effort has been made to characterize vent distribution, it has not been possible to reach a consensus about the methodological approach that guarantees extraction of the larger amount of reliable information from a given set of vents, or even about the type of information that can be obtained under all circumstances.

Failure in understanding the differences between various methodological approaches might result in a negative influence on the scientific and technical acceptability of models which ultimately is detrimental for the advance of scientific knowledge in general. From a practical point of view, the mentioned failure might result in rejection of works submitted to scientific journals on the grounds of alleged errors, ignorance about the method, or indulgence on bespoken methods that produce what the author wants to see. Such statements are difficult to discuss in a formal context because most journals nowadays do not accept the inclusion of references to personal communications or to unpublished works. Nevertheless, those comments play an important role in shaping the evolution of scientific knowledge because acritical acceptance of such statements by some editors leads to rejections of works that are based on alternative methodologies and eliminate the possibility of having a truly open discussion on those subjects. At the end, acceptance of those criticisms followed by rejection of the papers that present alternative methodologies results in the propagation of unilateral points of view that may not have factual support in every circumstance. Although such extreme situations are not the rule across the scientific literature, it is in the best interest of science to reduce the occurrence of such events. This goal can be reached by keeping an open mind capable to work with multiple hypothesis and multiple methodologies to enrich the outcome of scientific enquiry in the most general terms. In addition, it is important to assist the geoscience community to grasp the fundamental aspects of different families of methods developed within theoretical statistics providing landmarks that can be used to assess the possible merits and faults of various works beyond personal opinions.

In this work I focus attention on several aspects of Kernel Density Estimation (KDE) in a volcanic context. Among others, KDE has been used with forecasting purposes in seismology by Frankel (1995) and Hiemer et al. (2014) and in a volcanic context by Connor et al. (2019). Nevertheless, aspects of KDE that have not been treated in detail include: a) the definition of the function used to estimate the probability-density of the data, b) fulfillment in the volcanic context (i.e., the physical world) of the conditions imposed when establishing that function, and c) the possible interpretations of that function in connection with the physical world. Although those issues have been discussed in the natural hazard literature (see references listed above) it would seem that the most fundamental aspects of those issues have not been assimilated by the community at large. Thus, in order to make accessible all of this information to the widest possible audience, the presentation in this work leans towards a more colloquial style. In addition, because the intention of this work is not to finger point any particular work that was produced by adopting alternative approaches than those described below (that judgement is better left to each reader), there are sections of the text that are presented without providing specific references, although the most relevant references are provided in masse. Despite these characteristics, by clarifying here many aspects of KDE that have not been addressed before in the context of the spatial distribution of volcanic vents, it is hoped that the community will have more elements than currently available to make better judgements concerning the conclusions that have been reached whenever this type of methodology has been employed. More importantly, by facilitating comprehension of the several aspects of this method discussed below, it is hoped that not only others might be encouraged to use this powerful tool, but a wider diversity of hypothesis can be generated and tested, ultimately contributing to the better understanding of volcanic systems in general.

This section presents arguments showing that the most common situations of volcanic interest are likely to involve observations that may not satisfy the assumptions for which KDE was conceived. Aspects examined include the definition of KDE, its interpretation in the context of the physical origin of volcanic vents, and the contrast between the idealized theoretical scenario with the conditions imposed by the real world in volcanic contexts.

After deciding which kernel function K( ) will be used to calculate f ^ , it is necessary to decide a value for the parameter h in Eq. 1. The specific choice of K( ) determines the form in which h influences the calculated f ^ as well as the numeric range within which h can be chosen. For example, Cañón-Tapia (2013) showed that Fisher and Gaussian kernels display the exact opposite effect of the calculated f ^ as h is increased; also it was shown that the numerical range of h is entirely different for each kernel. Despite those differences, it was also shown that identical f ^ s can be obtained with both kernels if suitable hs are used in each case. Thus, although the choice of K( ) may not be very large, in order to examine the influence of h it is important to restrict attention to only one type of kernel function.

Thus, to facilitate the discussion, hereafter all calculations will be made with reference to a univariate Gaussian (or normal) kernel. Upon this selection of K( ), Eq. 1 takes the form: f ^ ( x , y ) = 1 2 π N h 2 ∑ i = 1 N exp ( − 1 2 [ d i h ] 2 ) (2)

It must be noted that on Equation. 2, f ^ is evaluated at points with coordinates (x, y). Thus, the i-th observation (volcanic vent) has coordinates (Xi, Yi), and is at a distance di from the evaluation point. There are no prescribed units of h, but those units must coincide with the units used to express the separation between vents and evaluation points (i.e., d _i ). Coordinates for the vents and evaluation points can be latitude—longitude pairs, in which case distances must be calculated using as a reference a sphere of radius equal to 6,371 km, or another suitable ellipsoid to avoid the effect of curvature of the surface of a planet. If coordinates are given in other systems that already incorporate the effect of the curvature of the surface of the planet (e.g., UTM) the distances can be calculated as they would be in a Cartesian plane, and h can be expressed in either m or km, as deemed convenient.

It is remarked that even when the observations and evaluation points are represented as points in a plane, equation 2 is not strictly bivariate. Indeed, it requires only one variable (the distance between two points) for the calculation of f ^ . A truly bivariate form of Eq. 2 can be constructed, but in this case h is not longer a scalar. As explained by Wand and Jones (1993), in a bivariate version of Eq. 2 there can be between one and three independent smoothing parameters that are incorporated as elements of a 2 × 2 matrix. One element is a smoothing factor along the x direction, another element is the smoothing factor along the y direction and the off-diagonal elements on the matrix are equal to each other and allow the resulting f ^ to have axial orientations that are not parallel to the x or y directions. In those cases, if the h-matrix is visualized as representing sample variances, the off-diagonal elements can be visualized as proportional to the covariance between the x and y variances. Other changes to Eq. 2 are required to fully incorporate the effect of increasing the number of dimensions (variables) for which f ^ will be calculated. Nevertheless, the main principle behind Eq. 2 remains unaltered. Consequently, even when attention will be mostly focused on the remainder of this work on examining the effect of h in the univariate kernel represented by Eq. 2, the conclusions remain valid for the fully bivariate kernel, and actually for any kernel applied to higher dimension cases.

In those situations in which the real PDF from which the N observations have been drawn is known (that we should recall has been referred as F), it is reasonable to determine which of the various estimators of f ^ (each obtained with a different h) is the one that better approximates F. The most common approaches to address that issue involve the determination of the difference between the known F and f ^ , or the square of that difference. Arguments in favour and against the use of the Integrated Square Error (ISE), the Mean Integrated Square Error (MISE), or the Asymptotic Mean Squared Error (AMISE) as a measure of the difference between F and f ^ have been proposed (Jones, 1991; Jones et al., 1996). Other criteria examined also include the Mean Squared Error (MSE), the Asymptotic MSE (AMSE) and the Sum of the AMSE (SAMSE) (Wand and Jones, 1994; Duong and Hazelton, 2003). In addition to the criteria used to characterize the error, many different approaches have been proposed to minimize those measures of error. The most common methods include variants of Cross-Validation (Least squares, Biased, Smoothed, etc.), variants of Plug-in (Park and Marron, Implemented refined, Bootstrap, etc.), mixing and many other alternatives (Heidenreich et al., 2013). The list is so large that Heidenreich et al. (2013) counted more than 30 bandwidth selection methods, and probably some more have been added since that time. This count relates to bandwidth selectors of univariate kernels. If the methods extend to include bivariate kernels, the number of selectors is even larger. The problem then resides not so much in having one method for the selection of an “optimal”, “objective” bandwidth, but on deciding which method should be chosen (Duong, 2007).

Many comparisons of bandwidth selectors (hereafter referred only as selectors) have been made (Heidenreich et al., 2013; Schindler, 2011; Turlach, 1993 and references therein). Invariably those works have concluded that it has not been possible to find one selector that uniformly performs better than all alternatives when a variety of complex Fs are considered. In practice this indicates that the selector adopted in one situation may or may not have been the best choice, depending on the real F that was under study. In other words, the only possible form to be certain that the selector is the most adequate for the F that needs to be characterized is to know F before making that study. Thus, from a pragmatic point of view it is worth to question, if F is already known, therefore allowing us to calculate very precisely the error of a particular f ^ , what is the utility of the calculation of f ^ ? Furthermore, if F is already known, why should we embark in a series of complex calculations to select only one value of h, that may or may not produce a truly significant f ^ , depending on the selector of choice? Similarly, if F is not known with certainty, how can we be certain that a particular f ^ is not significant for the case at hand?

Evidently, it can be argued that any unbiased estimation of h is preferable over a subjective selection of a single h because it should lead to a more robust (stable, reproducible) estimation of f ^ , and even when that f ^ may not be an exact determination of the real F, its robustness justifies its selection. Another argument can also be raised in the sense that adoption of the most stable selector increases the possibility of adopting an h that is closer to the real F than the selection of a random h based on purely subjective biases. Although these arguments are not devoid of certain allure, a thorough scientist is likely to ask for more evidence that proves that indeed the real F is better approximated by what is considered to be at that moment the most stable selector. Unfortunately, if the real F is unknown, any argument claiming that a given f ^ is close or far from it results speculative in essence.

To avoid discussions concerning how close a given f ^ might be to an unknown F, it might be convenient to address the issue in a slightly different form. To begin with, it is necessary to realize that trying to decide whether a method for the selection of h based on the AMISE approach is better than one adopting a MISE or SAMSE approach might become a sterile discussion if all the selectors ultimately lead to similar values of h. On the other hand, if different methods lead to very different estimates of h, the discussion of which method should be adopted takes a greater relevance and can be approached from a different perspective. Thus, from a pragmatic point of view, it seems more illustrative to establish to what extent, if any, the selection of the method influences the value of h that is to be inserted in Eq. 2, and most importantly how different are the variously calculated f ^ s.

An indication that the method used to select an “objective” h exerts a large influence on the calculated f ^ is provided by Duong (2007). As shown in Figure 2, the range of possible f ^ s that can be calculated, each with a different h resulting from one specific selector, is particularly large. Although all the f ^ s coincide in identifying a distribution that is elongated in a NW-SE direction, some suggest an elliptical, unimodal distribution, while others indicate bimodality, three main modes, and even a fourth detached mode in the NW corner of the area of study. It is important to remark that the differences amongst the various panels of Figure 2 sometimes represent minute differences in the method employed to select the “optimal” h. For example, all three panels on the left column were obtained with a Plug-in selector that used a SAMSE approach, differing only on the form in which the data and the covariance matrix are used within the selector. Even so, the estimated f ^ s can variously have one, three or four modes.

The variability of possible results shown in Figure 2, all of which can be claimed to have been calculated in an “unbiased and objective form” indicates that deciding which selector is used turns out to be a subjective choice that includes several assumptions concerning the form in which errors should be handled (MISE, SAMSE, etc.), and even about the specific form of the distribution that underlies the observed group of vents. Thus, whether knowingly or not, the choice of a selector for the bandwidth, and consequently the resulting selection of just one h, are prone to bias. The cumbersome procedure followed and the involvement of a measure of error, however, contribute to disguise such subjective decisions behind an apparently objective choice.

Another important aspect of the results shown in Figure 2 is that there is no form that an outside observer could decide with 100% certainty which f ^ is more accurate (or has less error), because such observer does not know if the real F is unimodal, bimodal, etc. Thus, if a truly unbiased selection of h is to be made where a complex F is expected, as is probably the case in situations of volcanic interest, it would be necessary to test a variety of selectors, all of which claim to be entirely “objective”. The prospective of having to compare the results of a number of bandwidth selectors to decide which is the “more objective”, or the “less biased”, is not that different from the exploratory approach for the selection of h. This approach consists of choosing different values of h to visually inspect the resulting f ^ s (Silverman, 1981). Based on that inspection, a decision can be made about which h is deemed more appropriate, if any. The last remark is important because upon inspection of several estimates of F it might be concluded that more than one h is required to reveal various features of the underlying distribution. Thus, the exploratory method is more akin to the idea of multiple working hypothesis than it is to the uncompromised acceptance of a unique ruling theory (Chamberlin, 1897).

The exploratory approach has been deemed as very dated, also implying that it is undesirable because of its propensity to be affected by a subjective bias (Connor et al., 2019). Nevertheless, as discussed above, due to the diversity of alternative methods for the selection of h, and due to the striking differences that can be obtained from the adoption of even apparently similar methods for that “unbiased” selection, it cannot be asserted that use of a complicated selector is entirely devoid of biases. All things considered, perhaps the only difference is that when the exploratory approach is adopted the subjective choices are on the first plane whereas when an arbitrary selector is favoured the subjective choices are hidden behind the expectations (conscious or subconscious) of the analyst. In consequence, given the fact that subjective decisions are always involved in the selection of either a value of h directly, or through the selection of a method to obtain a supposedly “optimal” h, the merits of the exploration approach should not be dismissed without further examination.

In summary, in this section it has been shown that unless the real distribution from which the vents have been drawn, F, is known with 100% accuracy, it is not possible to decide in an unbiased form which value of h needs to be chosen to minimize the error between the estimated f ^ and the real F. If all things are considered, in a thorough study it is necessary to compare the results of different selectors to ensure that a really significant h is identified. From a pragmatic point of view, it seems much simpler to eliminate the middleman by examining directly several f ^ s obtained with a diversity of values of h. This approach is further explored in the following section.

The theory discussed in the previous section is applied to a few case studies in this section. Thus, in this Section I present the results of the various f ^ s calculated with various selectors available in the environment R available at https://www.r-project.org/. Vent locations from three different zones of volcanic interest were used as input to calculate the corresponding f ^ s. Each of the three zones was selected because there is enough information supporting identification of several key features. Previous geological knowledge about those zones is as close as we can get at this time to knowing the real F. In each of the following subsections several aspects concerning the volcanic systems that should ideally be identified or at least suggested by an analysis of volcanic vent distribution are described before presenting the corresponding f ^ s. Thus, evaluation of the predictive power of the analysis can be made in base of the features that can be identified on each of the calculated f ^ s. It is important to remark that those features have been established from the geologic studies of those zones, and therefore their validity does not depend in any form from the interpretation of the results of vent distributions. Actually, it is precisely because those features have been established from the geology that they can serve as the reference against which the predictions of the kernel method can be compared.

Although there is a large number of selectors that have been advanced in the literature, it is beyond the purpose of this work to make an exhaustive comparison of all of them. Furthermore, as many of those selectors might require access to very specialized programs or to algorithms that might not be easily accessible for the bulk of the geoscience community, the comparison of selectors presented below was restricted to only five. The five selectors included are those available in the generic function “density” of the R environment. These options include nrd0 and nrd which are the rule of thumb proposed by Silverman (1986) and its variation as proposed by Scott (1992), respectively; ucv and bcv, which are unbiased and biased cross-validation methods and sj corresponding to the selector presented by Sheather and Jones (1991). Despite the limited number of selectors presented here, the results serve to illustrate the variability that might characterize the use of different selectors, and therefore suffice to illustrate the main aspect highlighted in the discussion.

All the vent locations used are in Latitude-Longitude pairs, but these are not directly inserted into the functions available in R, which is designed to deal with Cartesian coordinates on a plane. Thus, to avoid issues related to the distortion associated with the shape of the Earth, Latitude-Longitude pairs were converted to UTM before using the density function in R. The resulting hs were then converted to km to produce the diagrams shown in the following sections.

As pointed out by Chamberlin (1965) investigations often proceed on the presumption that there is a definite process through which all results are of maximum excellence, and therefore the question of ‘what is the best method?’ is more often asked than ‘what are the special values of different methods, and what are their several advantageous applicabilities?‘. This frame of mind clearly promotes the point of view that problems often arise when assessment teams do not understand how to execute a specific method, when in reality there might not be errors, but only differences in opinion about the outcomes of different methods. Furthermore, in many cases it is implicitly assumed that the selected method of study is well conceived for the purposes of a specific application, which as shown in section 2, may not be always the case for KDE in many volcanic contexts. Nevertheless, as shown by the examples of section 4, the examination of a sequence of PDFs can provide enough clues to assess which of the conceptual situations depicted in Figure 1 seems better suited to describe a particular situation, or at least can be used to formulate better informed hypothesis that can be used to guide future studies in a region.

This approach has been well recognized in the context of inferential analysis on some branches of Earth sciences (e. g., Budnitz et al., 1997), but is not commonly appreciated in the context of statistical applications.

Thus, in parallel to this conceptual issue, throughout this work it has been shown that there is a large diversity of “unbiased” and/or “objective” estimators that can be used to select a single value of h that in turn can be used to produce a single estimator of the real PDF associated with a group of vents within a mainly statistical approach in mind. Because the real distribution of vents is a priori unknown, the range of “optimal” hs is not well constrained. Because some methods of estimation of h work better with some types of distribution than others, and also because we ignore the real distribution when examining the location of vents, the only form to be thorough in our analysis is to use more than one method to estimate h. In so doing, we are not far from adopting the exploratory approach in which a range of hs is used to produce a sequence of f ^ s, all of which are considered of some relevance. Examination of the complete sequence commonly provides clues to assess some aspects as probably spurious. Additionally, examination of a sequence of f ^ s is the only form to be sure that an analysis is not limited by the assumptions inherent to a theoretical-statistical point of view that may not correspond to the particular geologic-tectonic-volcanic situation under examination. Although in many cases it is probable that limiting the analysis to only one value of h may not be misleading, it is impossible to be sure if that is the case because, after all, we do not know for certain the real distribution. Consequently, despite the intentions of the researcher, any study that is based in only one h obtained by resorting to only one method of estimation, runs the risk of not only missing important information, but also it may even yield misleading clues concerning the physical characteristics of the volcanic system under examination.