Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
43,881 characters · 35 sections · 22 citation commands
critband: A Python Package for Critical Bandwidth Analysis of Multimodal Distributions
\newgeometry{ a4paper, left=0.9in, right=0.9in, top=0.9in, bottom=0.8in, } \footnotetext[1]{The critband Python package can be installed via: \colorbox{verylightgray}{pip install critband}. For documentation and tutorials, visit: \url{https://pypi.org/project/critband/}.}
\restoregeometry
The question of whether a distribution is unimodal or multimodal is a fundamental inferential problem in scientific computing that cuts across nearly every empirical discipline. In ecology, bimodality in body size distributions signals niche partitioning holling1992; in economics, distributional separation is central to the measurement of polarization esteban1994; in genomics, bimodal gene-expression patterns can mark biologically meaningful patient or cellular subgroups bessarabova2010; and in astronomy, the bimodal color distribution of galaxies separates the red sequence from the blue cloud baldry2004. Yet despite the centrality of this question, the computational tools for multimodal assessment---particularly the critical bandwidth method and its associated numerical machinery---have remained fragmented in the Python scientific computing ecosystem.
This paper introduces critband, a Python package for critical bandwidth bimodality detection based on Silverman's kernel density approach silverman1981, wand1995. The critical bandwidth is defined as the smallest bandwidth $h$ for which a Gaussian kernel density estimate (KDE) of a sample has at most $k$ modes. When $k = 2$, the critical bandwidth $h_{\text{crit}}$ quantifies the separation between two components: a large $h_{\text{crit}}$ indicates strong, well-separated bimodality, while a small $h_{\text{crit}}$ suggests borderline or weak bimodality.
The R ecosystem has long dominated this analytical space. The multimode package ameijeirasalonso2021 provides the modetest() function implementing Silverman's test with bootstrap calibration. The multimode package also provides nmodes() for mode counting at a given bandwidth. These tools are mature, well-cited, and reliable. However, Python---despite its growing dominance in scientific computing---has lacked an equivalent. SciPy's gaussian_kde performs density estimation but provides no multimodality testing. Libraries such as sklearn.mixture focus on parametric mixture models rather than nonparametric mode assessment. The user who wishes to compute a critical bandwidth, run Silverman's test, or decompose a bimodal distribution into components has had no cohesive Python solution.
critband fills this gap. It implements a bracketed mode-count critical bandwidth solver with FFT-accelerated KDE, and provides additional features including $k$-mode detection, component decomposition, bimodality strength quantification, and excess mass estimation.
In the following sections, we review the theoretical background (Section (ref)), describe the package architecture and API (Section (ref)), validate against twelve benchmark cases with R cross-comparison (Section (ref)), report performance benchmarks (Section (ref)), highlight additional features (Section (ref)), demonstrate the package on three real-world-inspired application examples (Section (ref)), and conclude with discussion (Section (ref)).
Given an independent and identically distributed sample $X_1, \dots, X_n$ from an unknown density $f$, the Gaussian kernel density estimate is \[ \hat{f}_h(x) = \frac{1}{n h \sqrt{2\pi}} \sum_{i=1}^{n} \exp\left(-\frac{(x - X_i)^2}{2h^2}\right), \] where $h > 0$ is the bandwidth (smoothing parameter). The bandwidth governs the bias-variance tradeoff: small $h$ produces wiggly estimates that may capture spurious structure, while large $h$ oversmooths and may obscure genuine features. A natural starting point is the Silverman-style robust rule-of-thumb bandwidth used by the package silverman1986, jones1996: \[ \hat{h}_{\text{silverman}} = 1.06 \cdot \min\left(\hat{\sigma}, \frac{\text{IQR}}{1.34}\right) \cdot n^{-1/5}, \] where $\hat{\sigma}$ is the sample standard deviation and IQR is the interquartile range. This implementation combines the $1.06$ normal-reference constant with the robust scale estimate $\min(\hat{\sigma}, \mathrm{IQR}/1.34)$; the alternative $0.9$ constant is often used with the robust scale in density-estimation texts. We retain the implemented convention here because it defines the public silverman_bandwidth(x) API used throughout the validation.
Silverman silverman1981 proposed using the KDE as a basis for multimodality testing. The critical bandwidth $h_{\text{crit}}(k)$ is defined as \[ h_{\text{crit}}(k) = \inf\{h > 0 : \hat{f}_h \text{ has at most } k \text{ modes}\}. \]
When $k = 1$, the critical bandwidth is the smallest $h$ such that $\hat{f}_h$ is unimodal; for $k = 2$, it is the smallest $h$ such that $\hat{f}_h$ has at most two modes. The null hypothesis of $k$-modality is rejected when $h_{\text{crit}}(k)$ is larger than what one would expect under the null. Because the null distribution of $h_{\text{crit}}(k)$ depends on the unknown true density, Silverman proposed a bootstrap procedure: resample from a $k$-modal density estimate with bandwidth $h_{\text{crit}}(k)$, recompute the critical bandwidth on each bootstrap sample, and compute the $p$-value as the proportion of bootstrap replicates where $h_{\text{crit}}^* \geq h_{\text{crit}}(k)$. The calibration properties of this procedure and related excess-mass alternatives have been studied extensively hall2001, fisher2001.
Hartigan and Hartigan hartigan1985 introduced the dip test, which measures unimodality by computing the maximum difference between the empirical cumulative distribution function (ECDF) and the unimodal distribution function that best approximates it. The dip statistic is \[ D = \inf_{F \in \mathcal{U}} \sup_x |F_n(x) - F(x)|, \] where $\mathcal{U}$ is the class of unimodal distribution functions and $F_n$ is the ECDF. Large values of $D$ indicate departure from unimodality. The dip test is complementary to Silverman's test: it assesses unimodality directly rather than through the KDE smoothing parameter, and it requires no tuning of bandwidth choices.
M\"uller and Sawitzki muller1991 proposed the excess mass approach, which measures the amount of probability mass that lies above a threshold $p$ for intervals where the density exceeds that threshold. For a $k$-modal density, the excess mass functional is \[ E(p) = \sum_{j=1}^{k} \int_{I_j(p)} (f(x) - p) \, dx, \] where $I_j(p)$ are the intervals above the threshold $p$. The slope of the excess mass function provides information about the number and separation of modes. The critband package implements this as excess_mass().
The R package multimode ameijeirasalonso2021 is the current gold standard. Its modetest() function performs Silverman's test with several bandwidth calibration methods. The package also provides nmodes() for mode counting at a specified bandwidth and locmodes() for locating mode positions. These tools are well-maintained and widely cited.
In Python, no equivalent package existed prior to critband. Some functionality is distributed across libraries: scipy.stats.gaussian_kde provides KDE evaluation virtanen2020; scipy.signal.find_peaks can count modes on an evaluated density; and sklearn.mixture.GaussianMixture pedregosa2011 fits parametric mixtures but does not test for the number of components in a nonparametric sense. There is no cohesive implementation of critical bandwidth search, Silverman's bootstrap test, or component decomposition.
critband is designed around three principles:
The critical bandwidth computation poses several algorithmic challenges that shaped the design of critband.
Mode-count search and trough-ratio diagnostics. Silverman's original formulation defines $h_{\text{crit}}$ as the infimum over bandwidths yielding at most $k$ modes---a discrete property. critband therefore uses the KDE mode count itself as the default search criterion. The implementation also exposes a trough-ratio diagnostic $r(h) = f_{\text{valley}} / \max(f_{\text{peak1}}, f_{\text{peak2}})$, where $f_{\text{valley}}$ is the KDE height at the lowest point between the two largest peaks. This diagnostic is useful for describing bimodality strength and for the optional Brent-based solver path, but the reported validation below uses the robust mode-count search.
Bracketed binary solver. Pure binary search on the mode-count criterion is robust ($O(\log(1/\varepsilon))$ iterations) and directly targets the definition of $h_{\text{crit}}$. The default “auto” strategy first narrows the search interval with a short bracketing pass and then applies binary search inside that bracket. The explicit method="brent" option is kept for exploratory solver comparisons and falls back to binary search if the continuous objective does not produce a verified transition.
FFT-accelerated KDE. The dominant computational cost is the inner KDE evaluation loop. Direct evaluation of $\hat{f}_h$ at $g$ grid points from $n$ samples is $O(n \cdot g)$. For $n > 5000$, critband switches to an FFT-based algorithm that discretizes the data onto the grid and convolves with a Gaussian filter in the frequency domain, reducing complexity to $O(g \log g)$. The grid size scales adaptively as $g = \max(800, \min(5000, n/2))$, ensuring fine resolution for small-bandwidth searches without wasting computation for large samples.
The package exports the following primary functions:
The following example demonstrates a typical workflow: load data from a CSV file, compute the critical bandwidth, run both Silverman's test and Hartigan's dip test, and decompose the distribution.
The critband.io module provides two entry points:
The first reads from a file path with automatic format detection based on extension. The second reads from a BytesIO buffer, making it suitable for web applications. Supported formats include CSV, TSV/TXT, JSON, Markdown pipe tables, HTML tables, XLSX, XLS, DOCX tables, and PDF tables. The DOCX and PDF adapters target tabular numeric extraction rather than arbitrary document understanding.
We evaluated critband against twelve benchmark cases spanning a range of bimodality scenarios. Each case is a Gaussian mixture with known parameters. For each case, we report the mean critical bandwidth $\bar{h}_{\text{crit}}$ and its standard deviation across 10 independent random seeds, the coefficient of variation (CV$\%$), and the number of modes $\hat{k}$ detected by find_modes at Silverman's bandwidth.
The critical bandwidth values range from $0.226$ (barely separated) to $4.682$ (extreme separation), correctly reflecting the degree of bimodal separation in each case. The ten-seed Monte Carlo experiment confirms that the $h_{\text{crit}}$ estimates are stable for clearly separated cases: the coefficient of variation is below $5\%$ for nine of the twelve cases. The three high-CV cases are all boundary or overlap cases near the bimodal--unimodal transition, where small sampling changes can flip the detected mode count.
We compared critband's results against the R multimode package. Each case was generated with a unique random seed to ensure reproducibility across environments.
$^\dagger$R's nmodes(bw = bw.nrd0()) returns 772 spurious modes on the unequal-weights mixture due to severe undersmoothing of the dominant component's tail. The alternative Sheather--Jones bandwidth selector (bw.SJ) produces even more undersmoothing (744--1113 modes across five random seeds), confirming that this is not a bw.nrd0-specific pathology but a fundamental limitation of automatic bandwidth selection for unequal-weight mixtures. critband's find_modes() correctly returns 2 modes by using an explicit Silverman-rule bandwidth with user control.
Key observations:
For most benchmark cases, critband's find_modes(x, h=silverman_bandwidth(x)) returns the expected number of modes (2 for bimodal cases and 3 for the trimodal case). The boundary and overlap cases occasionally fluctuate between one and two modes across seeds, which is the expected behavior near the unimodal--bimodal transition rather than a stable failure of the mode counter.
To assess how critband's nonparametric multimodality summary compares with parametric model selection, we evaluated eight benchmark cases using both critband's critical bandwidth calculation and scikit-learn's GaussianMixture (GMM) with Bayesian Information Criterion (BIC) and Akaike Information Criterion (AIC). For each case, GMM was fitted with $k = 1, 2, 3$ components and the model minimizing each criterion was selected. We report the continuous critical-bandwidth value and the mode count, while leaving hypothesis-test $p$-values to Table (ref); this avoids mixing R modetest calibration with a separate Python bootstrap configuration.
Several observations emerge. First, BIC and AIC select $k = 2$ for all seven two-component cases and $k = 3$ for the trimodal case. This reflects a fundamental difference in inferential target: BIC selects the model with the best information-theoretic fit under a Gaussian-mixture family, whereas critband's nonparametric workflow reports a critical bandwidth, a mode count, and a bootstrap test of unimodality. The two approaches should therefore be read as complementary rather than interchangeable.
Second, both methods identify the trimodal case, with BIC strongly favoring $k = 3$ ($\Delta\text{BIC} \approx 560$ over $k = 2$). This confirms that critband's nonparametric mode detection and GMM's parametric selection converge on the correct structure when modes are well separated.
Third, critband offers information that GMM does not. The continuous $h_{\text{crit}}$ metric ranges from $0.28$ (barely separated) to $1.88$ (well-separated), quantifying the degree of smoothing required before modes merge. GMM's selected component count is discrete. For the near-unimodal case ($h_{\text{crit}} = 0.42$, only ${\sim}1.5\times$ the Silverman bandwidth), critband reveals that the distribution is weakly separated and close to the unimodal threshold---information lost by BIC's unambiguous $k = 2$ selection.
critband and GMM are complementary rather than competing approaches. critband's nonparametric workflow makes no distributional assumptions about component shape and provides a direct mode-based summary. GMM provides statistically efficient parameter estimation when the Gaussian assumption holds. We recommend critband for the initial exploratory question “does this distribution show modal separation?” and GMM for the follow-up “what are the component parameters?” when Gaussian components are substantively plausible.
We benchmarked critband against the relevant R package on the twelve benchmark cases. All timings were collected on a standard workstation (Apple M2, 16 GB RAM). R timings used modetest(B = 199) and nmodes(bw = bw.nrd0()). Python timings used critical_bandwidth() (default settings) and find_modes().
The critical bandwidth search in critband is 3--10 times faster than R's modetest on individual benchmark cases. This speedup arises from three design choices:
We measured the scaling behavior of critical_bandwidth() across sample sizes $n = \{100, 1000, 10000\}$ using the well-separated benchmark case:
The transition to FFT-based KDE at $n > 5000$ kept the runtime sub-second at $n = 10000$ in this local benchmark. The R comparison at this scale should be read as a machine- and configuration-specific timing reference rather than as a portable performance guarantee.
Beyond replicating existing R functionality, critband provides several additional features not directly available in the multimode R package:
The critical_bandwidth(x, k=n) function accepts any integer $k \geq 1$, enabling detection of trimodality ($k = 3$), quadmodality ($k = 4$), and higher-order multimodality. The trimodal case in our benchmark suite confirms this: critical_bandwidth(x, k=3) returns $h_{\text{crit}} = 1.3810$ for the $\mathcal{N}(-3,0.3) \cup \mathcal{N}(0,0.3) \cup \mathcal{N}(3,0.3)$ mixture. R's modetest supports $k$-mode testing as well; critband provides a convenient interface for the same functionality.
The detect_components() function decomposes a bimodal distribution into two Gaussian components by splitting the data at the KDE trough. It returns a BimodalDecomposition object containing:
This decomposition is valuable for exploratory downstream analysis: in genomics, the component means and weights can summarize active vs.\ silenced expression states; in economics, they can summarize the apparent gap between distributional groups. Because the method is a KDE trough split rather than a generative model, these summaries should not be treated as confirmatory component estimates without additional checks.
The return_ci option on critical_bandwidth() provides a bootstrap estimate of the sampling distribution of $h_{\text{crit}}$ efron1993:
The bootstrap procedure resamples the observed data, computes $h_{\text{crit}}^*$ on each resample, and returns a bootstrap interval along with the standard error of $h_{\text{crit}}$. For the well-separated case with $n = 400$, we obtain $h_{\text{crit}} = 1.8650$ with a 95% bootstrap interval of $[1.72, 2.01]$, reflecting moderate sampling variability. We treat this interval as provisional uncertainty support rather than as a final inferential guarantee; broader coverage studies would be needed before using it as a standalone calibration target hall2001.
The bimodality_strength() function computes a quantitative, continuous measure of bimodality strength based on the ratio $h_{\text{crit}} / \hat{h}_{\text{silverman}}$. In the current benchmarks, higher ratios correspond to stronger separation and ratios below 1.0 correspond to weaker or borderline bimodality. We treat the weak/moderate/strong bands as heuristic summaries of the observed cases rather than calibrated decision thresholds. This provides a more nuanced assessment than a binary reject/fail-to-reject decision from a hypothesis test.
The excess_mass() function implements the M\"uller--Sawitzki excess mass approach muller1991, providing an alternative nonparametric test for multimodality. The excess mass functional is estimated from the KDE, and its slope at small thresholds provides diagnostic information about the number of modes. This is particularly useful for validating results from Silverman's test with an independent methodology.
In this section we illustrate the practical utility of critband on three synthetic datasets that mimic canonical scientific problems. Each example demonstrates a complete analysis workflow: computing the critical bandwidth, testing for bimodality, and interpreting the results in the context of the original scientific question.
A central question in community ecology is whether body size distributions of coexisting species exhibit multimodal structure, indicating niche partitioning holling1992. We simulate a community of 300 species drawn from a mixture of two log-normal components: small-bodied invertebrates ($\log\mathcal{N}(1.5, 0.3)$, 60% weight) and large-bodied vertebrates ($\log\mathcal{N}(3.0, 0.4)$, 40% weight), mimicking the size distribution of a forest-floor arthropod and small-mammal assemblage.
In this synthetic log-normal example, critband estimates $h_{\text{crit}} \approx 2.64$ and the default Silverman bootstrap path does not reject unimodality at conventional levels. The trough-split component summary remains descriptively close to the two simulated size groups, with means near $4.9$ and $23.2$ on the original scale. This example is therefore best read as an exploratory decomposition of a visibly structured sample, not as a strong hypothesis-test rejection.
Galaxy color distributions are known to separate galaxies into a “red sequence” of passively evolving galaxies and a “blue cloud” of star-forming galaxies baldry2004. We simulate 500 galaxies with a bimodal $(g-r)$ color distribution: a red sequence at $g-r = 0.8$ ($\sigma = 0.15$, 55% weight) and a blue cloud at $g-r = 0.3$ ($\sigma = 0.12$, 45% weight), with mild overlap reflecting the “green valley” transition.
The critical bandwidth is approximately $h_{\text{crit}} = 0.18$ in this synthetic example, and a 99-resample smoke run of Silverman's bootstrap path rejects unimodality. The bimodality strength ratio $h_{\text{crit}} / \hat{h}_{\text{silverman}}$ is about 2.1, indicating visible but not threshold-calibrated modal separation. Component decomposition yields means of about $0.30$ and $0.80$, accurately recovering the input parameters.
In genomics, bimodal expression can indicate biologically meaningful switching or subgroup structure rather than a single homogeneous expression regime bessarabova2010. We simulate nonnegative log$_2(\text{counts}+1)$ expression values for 200 cells: an active state centered at $4.0$ ($\sigma = 0.4$, 70% weight) and a silenced state centered at $0.5$ ($\sigma = 0.3$, 30% weight), truncated at zero to respect the lower bound of transformed count data.
Silverman's test rejects unimodality in this synthetic example. Component decomposition identifies groups near $\sim0.5$ (30% weight) and $\sim4.0$ (70% weight), matching the generative parameters up to component ordering. The component weights provide an interpretable exploratory summary of the simulated regulatory states---information that would be lost under a unimodal analysis that simply reports mean expression across all cells.
critband provides the Python ecosystem with a comprehensive, well-tested, and performant implementation of critical bandwidth bimodality detection. The package closes a long-standing gap: whereas R users have had access to multimode for over a decade, Python users were limited to piecemeal solutions involving manual KDE evaluation, ad-hoc mode counting, and bespoke bootstrap code.
Our validation against twelve benchmark cases demonstrates that critband produces stable critical bandwidth values on clearly separated synthetic cases and correctly flags boundary cases as unstable. The R cross-validation supports the same broad modality conclusions through modetest $p$-values and nmodes checks. It also revealed a practical limitation of R's default bandwidth selector on unequal-weight mixtures---a problem that critband avoids through explicit bandwidth specification and robust mode counting.
The performance benchmarks show that critband's critical bandwidth search delivers per-case speedups of 3--10$\times$ over R's modetest in the benchmark timing setup. This speedup is achieved through a bracketed mode-count search, FFT-accelerated KDE for large samples, and adaptive grid sizing. For interactive analysis, the core critical-bandwidth and component-decomposition steps complete quickly for typical dataset sizes; bootstrap tests remain governed by the requested number of resamples.
Beyond replication, critband provides several additional features: $k$-mode detection for arbitrary $k$, component decomposition via KDE trough splitting, bootstrap intervals for $h_{\text{crit}}$, a continuous bimodality strength metric, and excess mass estimation. The multi-format I/O subsystem makes the package suitable for production pipelines that ingest data from diverse sources.
critband is designed as a native component of the Python scientific computing ecosystem, building on NumPy harris2020 and SciPy virtanen2020. A user can pass a numpy.ndarray into any critband function, chain the output into scipy.stats for further analysis, or visualize results with Matplotlib---all without format conversion or data copying. The package occupies a specific niche between SciPy's density estimation (scipy.stats.gaussian_kde) and scikit-learn's parametric mixture modeling (sklearn.mixture.GaussianMixture): it provides nonparametric multimodality testing that neither library offers, filling a gap that has required R interop or bespoke implementation.
The package is distributed under the Apache 2.0 license and is intended for installation via pip install critband. In the local repository state used for this manuscript, the full test suite contains 289 tests. This provides useful regression coverage, though it should not be read as a guarantee of reliability across all platforms and Python versions.
The current implementation has several limitations that suggest directions for future work.
Optional Brent solver. The optional Brent path uses the trough-ratio objective $r(h) = f_{\text{valley}} / \max(f_{\text{peak1}}, f_{\text{peak2}})$ and therefore assumes useful continuity near the modal transition. In practice, discontinuities can occur when the relative heights of the two largest peaks change with $h$, causing the $\max$ function to switch between them. For this reason, the default solver remains the bracketed binary mode-count search, and the Brent path is treated as an exploratory option with binary fallback rather than as the evidentiary backbone of the paper.
Component decomposition. The detect_components() function splits the data at the KDE trough and computes per-side statistics. This is a fast and interpretable heuristic, but it lacks the statistical rigor of parametric approaches such as Gaussian mixture models (GMM) fitted via expectation-maximization with BIC selection (\S(ref)). The trough-split method can be biased when the true components overlap substantially or have unequal variances, and it does not provide uncertainty estimates for the component parameters. We recommend it as an exploratory tool, with GMM as the preferred option for downstream parameter estimation in hard cases.
Bimodality strength thresholds. The bimodality_strength() metric $h_{\text{crit}} / \hat{h}_{\text{silverman}}$ provides a useful continuous measure, but the weak/moderate/strong labels are heuristic summaries of the current benchmark set rather than calibrated decision rules. They are useful for descriptive reporting, but they should not be read as formal cutoffs.
Bootstrap confidence intervals. The bootstrap intervals for $h_{\text{crit}}$ (\S(ref)) are useful for exploratory uncertainty support, but they should still be treated cautiously until broader coverage studies are available. The current implementation uses a manual percentile bootstrap with failure accounting and explicit provenance; that is an improvement over the earlier percentile-only version, but it does not by itself justify strong inferential language.
Large-sample testing. Exact nonparametric tests can become expensive for datasets with $n > 5000$. For large-sample analyses, users may prefer Silverman's bootstrap test or a smaller n_boot setting, and the implementation now warns when the default bootstrap count is used on large samples.
Additional limitations. The current implementation defaults to the Gaussian kernel, the standard choice in much of the multimodality literature, while also exposing additional kernel options for robustness checks. The bootstrap implementation does not yet support parallel execution across multiple cores, a natural extension for large-scale simulation studies. Finally, an adaptive bandwidth selection method that performs well on unequal-weight mixtures would be a valuable addition, potentially combining Silverman's rule with a pilot estimate of component structure.
We gratefully acknowledge the developers of the multimode R package, whose work provided the validation benchmark for this project. We thank the open-source Python community for the foundational libraries (NumPy, SciPy) upon which critband is built.