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.
21,157 characters · 11 sections · 23 citation commands
Estimating Sloppy Directions via KDE: The Case of Kirman's Ants
\email{[email removed]}
Confronted with highly parameterized models and a limited set of observations, scientists can nonetheless make accurate predictions. The reason is that many models are sloppy: their predictions depend only on a handful of stiff parameter combinations, while the remaining sloppy combinations have little to no effect on the predicted outcomes BrownSethna2003StatisticalMechanicalApproaches,BrownEtAl2004StatisticalMechanicsComplex. The notion of stiff and sloppy parameter directions was first introduced by Brown and Sethna BrownSethna2003StatisticalMechanicalApproaches while studying models of biochemical regulation. For each of the individual parameter estimates, “95% confidence intervals each spanned more than a factor of 50” GutenkunstEtAl2007UniversallySloppyParameter, yet the models retained accurate predictive power because the stiff directions were well-constrained by data. Sloppiness has since been identified across cell signaling, radioactive decay, neural networks, quantum wave functions, the Ising model, and macroeconomic agent-based models QuinnEtAl2023InformationGeometryMultiparameter,WaterfallEtAl2006SloppyModelUniversalityClass,MachtaEtAl2013ParameterSpaceCompression.
The standard framework identifies stiff and sloppy directions via the eigendecomposition of the Hessian matrix of a loss function that measures the change in model predictions when parameters are perturbed. Two loss functions dominate the literature: the mean-squared error (MSE) on predicted time series, and the symmetrized Kullback--Leibler (sKL) divergence between predicted distributions KullbackLeibler1951InformationSufficiency,Jeffreys1948TheoryProbability. The Hessian of the sKL divergence equals the FIM QuinnEtAl2023InformationGeometryMultiparameter, a natural Riemannian metric on the manifold of model predictions TranstrumEtAl2011GeometryNonlinearLeast,TranstrumEtAl2010WhyAreNonlinear. Comparing distributions rather than time series also reduces the sensitivity to stochastic path noise and makes the distributional approach preferable for nonlinear stochastic models.
Yet existing worked examples have a common feature: the distribution being compared is either known analytically MachtaEtAl2013ParameterSpaceCompression or fit to a known parametric form GutenkunstEtAl2007UniversallySloppyParameter. In practice, for complex agent-based models (ABMs) the stationary distribution is unknown and must be estimated nonparametrically from simulation data. Prior work on sloppiness in ABMs Naumann-WoleskeEtAl2025ExplorationParameterSpace applied the sloppy-models framework to a high-dimensional macroeconomic ABM but relied on the MSE loss for stochastic time-series observables rather than a distributional loss. This leads to high uncertainty in the estimated stiff and sloppy directions, which in turns requires a large simulation budget, suggesting that an alternative distributional approach may be preferable but in this case there is no analytical prior on the distribution of outcomes. To our knowledge, whether a nonparametric estimate of the distribution recovers the true sloppy directions has not been demonstrated.
I address this gap using Kirman's ant recruitment model Kirman1993AntsRationalityRecruitment, a simple non-trivial stochastic model with a known analytical stationary distribution. With only two parameters $(\rho, \mu)$, the model has a closed-form FIM, providing an exact ground truth. I show that a Gaussian KDE with an Improved Sheather--Jones bandwidth BotevEtAl2010KernelDensityEstimation and a logit transform to handle the bounded state variable recovers the analytical eigenvectors and eigenvalues, and we characterize the convergence as a function of the simulation budget. Beyond validation, we demonstrate how the stiff direction enables efficient exploration of the model's phase diagram.
Consider a set of $N$ ants facing two equivalent food sources $A$ and $B$, where $k$ is the number of ants at source $A$. Each period, a randomly drawn ant may switch its food source with probability $\rho \in [0,1]$ (random switching), or alternatively recruit another randomly drawn ant with probability $\mu \in [0,1]$ (herding). Let $x_t = k/N$ denote the fraction of ants at source $A$ at time $t$.
Following Kirman1993AntsRationalityRecruitment and MoranEtAl2020SchrodingersAntsContinuous, in the $N\to \infty$ limit the dynamics of $x$ obey the It\^o stochastic differential equation
where $W_t$ is a standard Wiener process. This discretizes via Euler--Maruyama with step $\Delta t$ as
with $\xi_t \sim \mathcal{N}(0,1)$ i.i.d. For $\rho > 0$, the normalized stationary distribution is a symmetric Beta distribution Kirman1993AntsRationalityRecruitment,MoranEtAl2020SchrodingersAntsContinuous:
The ratio $\frac{\rho}{\mu}$ determines the qualitative behavior of the model: for $\frac{\rho}{\mu} < 1$ the distribution is bimodal (density diverging at the boundaries), for $\frac{\rho}{\mu} > 1$ it is unimodal (peaked at $x = \tfrac{1}{2}$), and $\frac{\rho}{\mu} = 1$ yields a uniform distribution. The model thus displays a phase transition at $\frac{\rho}{\mu} = 1$, controlled entirely by this single combination of the two parameters.
Figure (ref) shows representative dynamics and stationary distributions for the three regimes.
Throughout, we write $\Phi \equiv (\rho,\mu)^\top$ for the parameter vector and $\delta\Phi$ for an infinitesimal perturbation in log-parameter space. To identify stiff and sloppy directions, we decompose the Hessian matrix of a loss function $\mathcal{L}(\Phi, \Phi+\delta\Phi)$ that measures the change in model predictions under such a perturbation. Working in log-parameter space (to make derivatives scale-free and parameter-unit-agnostic), the $(i,j)$ element of the Hessian is
The second equality uses the fact that $\delta\Phi = 0$ is a stationary point of $\mathcal{L}^{sKL}(\Phi,\Phi+\delta\Phi)$, eliminating the chain-rule cross-term between log- and bare-parameter derivatives. By a second-order Taylor expansion, the loss is approximated as $\mathcal{L}(\Phi, \Phi+\delta\Phi) \approx \tfrac{1}{2}\delta\Phi^\top H \delta\Phi$, so the eigenvalue-eigenvector pairs $(\lambda_i, v_i)$ of $H$ identify the principal directions and their stiffness. The stiff direction (largest $\lambda$) induces the largest loss for a unit step; the sloppy direction (smallest $\lambda$) induces the smallest.
I use the symmetrized Kullback--Leibler divergence $\mathcal{L}^{sKL}(\Phi,\Phi+\delta\Phi) = \tfrac{1}{2}[KL\!\left(P\!\left(x|\Phi\right)\|P(x|\Phi^\star)\right) + KL\!\left(P(x|\Phi^\star)\|P\!\left(x|\Phi\right)\right)]$ as the loss function given its prominence in prior work. Its Hessian equals the FIM QuinnEtAl2023InformationGeometryMultiparameter:
For the Beta distribution in Eq. (ref) the FIM evaluates to
where $\psi^{(1)}(\cdot)$ is the trigamma function. The eigendecomposition yields
The stiff direction $v_1 \propto (-1,1)^\top$ corresponds to changes in the log-parameter ratio $\log(\rho/\mu)$, which is precisely the only parameter combination that affects the stationary distribution. The sloppy direction $v_2 \propto (1,1)^\top$ changes both $\log\rho$ and $\log\mu$ equally, leaving $\rho/\mu$ unchanged and hence all predictions invariant; its eigenvalue is exactly zero. I note here too that the MSE-based Hessian, computed on the log-variance of the Beta distribution as a representative moment, yields the same eigenvectors because the variance of $\mathrm{Beta}(\alpha,\alpha)$ depends only on $\alpha = \rho/\mu$, the same parameter combination that determines the entire distribution; the eigenvalue magnitude differs: $\lambda_1^{MSE} = 4(\rho/\mu)^2(2\rho/\mu+1)^{-2}$.
Figure (ref) shows $\lambda_1(\frac{\rho}{\mu})$ for both loss functions across the parameter space.
For general stochastic models the stationary distribution is unknown. I estimate $P(x|\Phi)$ nonparametrically using kernel density estimation (KDE) on the simulated time series $\mathcal{X} = \{x_{s,t}\}_{s=1,\dots,S;\;t=1,\dots,T}$, this being a common method for getting an estimate of the distribution (in recent work, Mottes et al. MottesEtAl2026GradientbasedOptimizationExact alternatively use differentiable histograms. However, this would then depend on the binning choice).
I apply a Gaussian KDE with bandwidth matrix $\mathbf{\Sigma}$:
While the kernel choice is not important, the bandwidth choice is critical Epanechnikov1969NonparametricEstimationMultivariate. I use the Improved Sheather--Jones (ISJ) bandwidth estimator BotevEtAl2010KernelDensityEstimation, which is data-adaptive and avoids the over-smoothing of rules-of-thumb.
The state variable $x \in [0,1]$ is bounded, causing standard KDE to incur boundary bias. Following KoekemoerSwanepoel2008TransformationKernelDensity and WandEtAl1991TransformationsDensityEstimation, we apply a logit transformation $g(x) = \log\frac{x}{1-x}$, mapping $x$ to $\mathbb{R}$, estimate the KDE on the transformed sample, and recover the density on $[0,1]$ via
Since the KL divergence is invariant under invertible transformations CsiszarShields2004InformationTheoryStatistics, the FIM computed in the transformed space equals that in the original space.
Given the KDE estimate $\hat{P}(x|\Phi)$, we compute the log-parameter derivatives of $\log\hat{P}$ by central finite differences with step size $\varepsilon$ in log-space:
where $e_i$ is the $i$-th unit vector. The numerical FIM element $\hat{H}^{sKL}_{ij}$ is then estimated by integrating over the discretized support,
averaged across the $S$ independent simulation seeds. The eigenvectors and eigenvalues are then obtained by SVD of the assembled $2\times 2$ matrix.
Figure (ref) shows the convergence of the estimated first eigenvector $\hat{v}_1$ and eigenvalue $\hat{\lambda}_1$ to their analytical counterparts as a function of the number of simulation seeds $S$ for a fixed simulation length $T=100$.
The KDE-based estimate converges to the analytical eigenvector $v_1 = (-1,1)^\top/\sqrt{2}$ for both the unimodal and bimodal regimes. Convergence in the bimodal case is somewhat slower, reflecting the heavier tails of the Beta distribution near the boundaries, which require more data for the logit-transformed KDE to resolve accurately. For simulation lengths $T=100$ and seed counts $S \geq 10$, the angular error of the first eigenvector falls below $5^\circ$ across the parameter space. The eigenvalue $\hat{\lambda}_1$ converges on a similar budget, though with larger relative error due to the additional integration step.
These results confirm that a standard nonparametric KDE pipeline recovers the theoretical sloppy structure of Kirman's ant model without requiring knowledge of the analytical distribution.
The stiff direction $v_1$ also serves as a tool for exploring the model's phase diagram. A second-order Taylor expansion shows that a step of size $\varepsilon$ in the direction of eigenvector $v_i$ produces a loss $\mathcal{L} \approx \tfrac{1}{2}\varepsilon^2\lambda_i$. Equivalently, to achieve a fixed loss magnitude $l$, the required step size is
This scaling law links the geometry of the loss surface to the “distance” between model predictions.
Figure (ref) demonstrates this relationship in its dual form: at fixed step $\varepsilon$ along $v_1$, the loss tracks $\tfrac{1}{2}\varepsilon^2\lambda_1$ across the parameter space. Equivalently, achieving a fixed loss requires step sizes $\varepsilon \propto 1/\sqrt{\lambda_1(\frac{\rho}{\mu})}$, smoothly traversing the model manifold from the unimodal phase through the critical point at $\frac{\rho}{\mu}=1$ into the bimodal phase, recovering all qualitatively distinct dynamics in a single directed walk through parameter space.
By following the stiff direction, one discovers the phase transition at $\frac{\rho}{\mu} = 1$ without any prior knowledge of the model's qualitative behavior. This illustrates a potential advantage of the sloppy approach in stochastic models more broadly: the exploratory path is derived from the model structure and simulation outputs alone, without requiring ex-ante knowledge of the modeler.
I have demonstrated that a standard KDE pipeline is sufficient to recover the analytical FIM eigenvectors and eigenvalues of Kirman's ant model. The sKL Hessian (FIM) has a single non-zero eigenvalue, with eigenvector $v_1 = (-1,1)^\top/\sqrt{2}$ corresponding to the log-ratio $\log(\rho/\mu)$, the only parameter combination affecting the stationary distribution. The KDE-based estimate converges to the analytical ground truth with modest simulation budgets ($S \geq 10$, $T \geq 10^6$).
This worked example provides a proof of concept for KDE-based FIM estimation when the stationary distribution is not known analytically. In contrast, the MSE approach on time series would instead remain sensitive to realization noise. Extension to the higher-dimensional parameter spaces and joint observables typical of full agent-based models is left to future work, also to test further the geometric decay of eigenvalues. The step-size scaling law further enables systematic phase exploration once the FIM eigenvectors are estimated, offering a principled alternative to exhaustive grid search or surrogate modeling tenBroekeEtAl2021UseSurrogateModels,LampertiEtAl2018AgentBasedModelCalibration.