Simulation-based Inference Research Articles

Simulation-based inference for parameter estimation of complex watershed simulators

Abstract. High-resolution, spatially distributed process-based (PB) simulators are widely employed in the study of complex catchment processes and their responses to a changing climate. However, calibrating these PB simulators using observed data remains a significant challenge due to several persistent issues, including the following: (1) intractability stemming from the computational demands and complex responses of simulators, which renders infeasible calculation of the conditional probability of parameters and data, and (2) uncertainty stemming from the choice of simplified representations of complex natural hydrologic processes. Here, we demonstrate how simulation-based inference (SBI) can help address both of these challenges with respect to parameter estimation. SBI uses a learned mapping between the parameter space and observed data to estimate parameters for the generation of calibrated simulations. To demonstrate the potential of SBI in hydrologic modeling, we conduct a set of synthetic experiments to infer two common physical parameters – Manning's coefficient and hydraulic conductivity – using a representation of a snowmelt-dominated catchment in Colorado, USA. We introduce novel deep-learning (DL) components to the SBI approach, including an “emulator” as a surrogate for the PB simulator to rapidly explore parameter responses. We also employ a density-based neural network to represent the joint probability of parameters and data without strong assumptions about its functional form. While addressing intractability, we also show that, if the simulator does not represent the system under study well enough, SBI can yield unreliable parameter estimates. Approaches to adopting the SBI framework for cases in which multiple simulator(s) may be adequate are introduced using a performance-weighting approach. The synthetic experiments presented here test the performance of SBI, using the relationship between the surrogate and PB simulators as a proxy for the real case.

Inferring Cosmological Parameters on SDSS via Domain-generalized Neural Networks and Light-cone Simulations

We present a proof-of-concept simulation-based inference on Ωm and σ 8 from the Sloan Digital Sky Survey (SDSS) Baryon Oscillation Spectroscopic Survey (BOSS) LOWZ Northern Galactic Cap (NGC) catalog using neural networks and domain generalization techniques without the need of summary statistics. Using rapid light-cone simulations L-picola, mock galaxy catalogs are produced that fully incorporate the observational effects. The collection of galaxies is fed as input to a point cloud-based network, Minkowski-PointNet . We also add relatively more accurate Gadget mocks to obtain robust and generalizable neural networks. By explicitly learning the representations that reduce the discrepancies between the two different data sets via the semantic alignment loss term, we show that the latent space configuration aligns into a single plane in which the two cosmological parameters form clear axes. Consequently, during inference, the SDSS BOSS LOWZ NGC catalog maps onto the plane, demonstrating effective generalization and improving prediction accuracy compared to non-generalized models. Results from the ensemble of 25 independently trained machines find Ωm = 0.339 ± 0.056 and σ 8 = 0.801 ± 0.061, inferred only from the distribution of galaxies in the light-cone slices without relying on any indirect summary statistics. A single machine that best adapts to the Gadget mocks yields a tighter prediction of Ωm = 0.282 ± 0.014 and σ 8 = 0.786 ± 0.036. We emphasize that adaptation across multiple domains can enhance the robustness of the neural networks in observational data.

From Simulations to Inference: Using Machine Learning to Tune Patient-Specific Finite-Element Models of the Middle Ear Towards Objective Diagnosis.

Computational models, particularly finite-element (FE) models, are essential for interpreting experimental data and predicting system behavior, especially when direct measurements are limited. A major challenge in tuning these models is the large number of parameters involved. Traditional methods, such as one-by-one sensitivity analyses, are time-consuming, subjective, and often return only a single set of parameter values, focusing on reproducing averaged data rather than capturing the full variability of experimental measurements. In this study, we applied simulation-based inference (SBI) using neural posterior estimation (NPE) to tune an FE model of the human middle ear. The training dataset consisted of 10,000 FE simulations of stapes velocity, ear-canal (EC) input impedance, and absorbance, paired with seven FE parameter values randomly sampled within plausible ranges. The neural network learned the association between parameters and simulation outcomes, returning the probability distribution of parameter values that can reproduce experimental data. Our approach successfully identified parameter sets that reproduced three experimental datasets simultaneously. By accounting for experimental noise and variability during training, the method provided a probability distribution of parameters, representing all valid combinations that could fit the data, rather than tuning to averaged values. The network demonstrated robustness to noise and exhibited an efficient learning curve due to the large training dataset. SBI offers an objective alternative to laborious sensitivity analyses, providing probability distributions for each parameter and uncovering interactions between them. This method can be applied to any biological FE model, and we demonstrated its effectiveness using a middle-ear model. Importantly, it holds promise for objective differential diagnosis of conductive hearing loss by providing insight into the mechanical properties of the middle ear.

Template switching during DNA replication is a prevalent source of adaptive gene amplification.

Copy number variants (CNVs)-gains and losses of genomic sequences-are an important source of genetic variation underlying rapid adaptation and genome evolution. However, despite their central role in evolution little is known about the factors that contribute to the structure, size, formation rate, and fitness effects of adaptive CNVs. Local genomic sequences are likely to be an important determinant of these properties. Whereas it is known that point mutation rates vary with genomic location and local DNA sequence features, the role of genome architecture in the formation, selection, and the resulting evolutionary dynamics of CNVs is poorly understood. Previously, we have found that the GAP1 gene in Saccharomyces cerevisiae undergoes frequent and repeated amplification and selection under long-term experimental evolution in glutamine-limiting conditions. The GAP1 gene has a unique genomic architecture consisting of two flanking long terminal repeats (LTRs) and a proximate origin of DNA replication (autonomously replicating sequence, ARS), which are likely to promote rapid GAP1 CNV formation. To test the role of these genomic elements on CNV-mediated adaptive evolution, we performed experimental evolution in glutamine-limited chemostats using engineered strains lacking either the adjacent LTRs, ARS, or all elements. Using a CNV reporter system and neural network simulation-based inference (nnSBI) we quantified the formation rate and fitness effect of CNVs for each strain. We find that although GAP1 CNVs repeatedly form and sweep to high frequency in strains with modified genome architecture, removal of local DNA elements significantly impacts the rate and fitness effect of CNVs and the rate of adaptation. We performed genome sequence analysis to define the molecular mechanisms of CNV formation for 177 CNV lineages. We find that across all four strain backgrounds, between 26% and 80% of all GAP1 CNVs are mediated by Origin Dependent Inverted Repeat Amplification (ODIRA) which results from template switching between the leading and lagging strand during DNA synthesis. In the absence of the local ARS, a distal ARS can mediate CNV formation via ODIRA. In the absence of local LTRs, homologous recombination mechanisms still mediate gene amplification following de novo insertion of retrotransposon elements at the locus. Our study demonstrates the remarkable plasticity of the genome and reveals that template switching during DNA replication is a frequent source of adaptive CNVs.

Bayesian inference of frequency-specific functional connectivity in MEG imaging using a spectral graph model

Abstract Understanding the relationship between structural connectivity (SC) and functional connectivity (FC) of the human brain is an important goal of neuroscience. Highly detailed mathematical models of neural masses exist that can simulate the interactions between functional activity and structural wiring. These models are often complex and require intensive computation. Most importantly, they do not provide a direct or intuitive interpretation of this structure–function relationship. In this study, we employ the emerging concepts of spectral graph theory to obtain this mapping in terms of graph harmonics, which are eigenvectors of the structural graph’s Laplacian matrix. In order to imbue these harmonics with biophysical underpinnings, we leverage recent advances in parsimonious spectral graph modeling (SGM) of brain activity. Here, we show that such a model can indeed be cast in terms of graph harmonics, and can provide a closed-form prediction of FC in an arbitrary frequency band. The model requires only three global, spatially invariant parameters, yet is capable of generating rich FC patterns in different frequency bands. Only a few harmonics are sufficient to reproduce realistic FC patterns. We applied the method to predict FC obtained from pairwise magnitude coherence of source-reconstructed resting-state magnetoencephalography (MEG) recordings of 36 healthy subjects. To enable efficient model inference, we adopted a deep neural network-based Bayesian procedure called simulation-based inference. Using this tool, we were able to speedily infer not only the single most likely model parameters, but also their full posterior distributions. We also implemented several other benchmark methods relating SC to FC, including graph diffusion and coupled neural mass models. The present method was shown to give the best performance overall. Notably, we discovered that a single biophysical parameterization is capable of fitting FCs from all relevant frequency bands simultaneously, an aspect that did not receive adequate attention in prior computational studies.

Simulation-based inference of black hole ringdowns in the time domain

Simulation-based inference of black hole ringdowns in the time domain

Ensemble surrogate modeling of advective-dispersive transport with intraparticle diffusion model for column-leaching test

Column-leaching tests are a common approach for assessing the leaching behavior and resulting environmental risks of contaminated soils and waste materials, which are frequently reused for various construction purposes. The observed breakthrough curves of the contaminants are influenced by the complex dynamics of solute transport and kinetic inter-phase mass transfer. Disentangling these interactions necessitates numerical models. However, inverse modeling and sensitivity analysis can be time-consuming, especially when sorption kinetics are explicitly described by intraparticle diffusion, which requires discretizing the domain both in the flow direction along the column axis and inside the grains. To circumvent the need for such computationally intensive models, we have developed two different ensemble surrogate models. These models employ two separate ensemble methods: random forest stacking and inverse-distance weighted interpolation. Each method is applied to base surrogate models that cover different parts of the parameter space. The base surrogate models use the method of Extremely randomized Trees (ExtraTrees). The defined parameter range is based on the German standard for column-leaching tests. To optimize the base surrogate models, we utilized adaptive-sampling methods based on three distinct infill criteria: maximizing the expected improvement, staying within a certain Mahalanobis distance to the best estimate (both for exploitation), and maximizing the standard deviation (for exploration). The ensemble surrogate model demonstrates excellent performance in emulating the behavior of the original numerical model, with a relative root mean squared error of 0.09. We applied our proposed ensemble surrogate model to estimate the complete posterior parameter distribution using Simulation-Based Inference, specifically Neural Posterior Estimation, to determine the full parameter distribution conditioned on copper-leaching data from two different soils. Samples drawn from the posterior distribution align perfectly with the observed data for both the surrogate and original models.

Constraining the Higgs potential with neural simulation-based inference for di-Higgs production

Determining the form of the Higgs potential is one of the most exciting challenges of modern particle physics. Higgs pair production directly probes the Higgs self-coupling and should be observed in the near future at the High-Luminosity LHC. We explore how to improve the sensitivity to physics beyond the Standard Model through per-event kinematics for di-Higgs events. In particular, we employ machine learning through simulation-based inference to estimate per-event likelihood ratios and gauge potential sensitivity gains from including this kinematic information. In terms of the Standard Model Effective Field Theory, we find that adding a limited number of observables can help to remove degeneracies in Wilson coefficient likelihoods and significantly improve the experimental sensitivity. Published by the American Physical Society 2024

Inferring galaxy cluster masses from cosmic microwave background lensing with neural simulation based inference

Gravitational lensing by massive galaxy clusters distorts the observed cosmic microwave background (CMB) on arcminute scales, and these distortions carry information about cluster masses. Standard approaches to extracting cluster mass constraints from the CMB cluster lensing signal are either sub-optimal, ignore important physical or observational effects, are computationally intractable, or require additional work to turn the lensing measurements into constraints on cluster masses. We apply simulation based inference (SBI) using neural likelihood models to the problem. We show that in circumstances where the exact likelihood can be computed, the SBI constraints on cluster masses are in agreement with the exact likelihood, demonstrating that the SBI constraints are close to optimal. In scenarios where the exact likelihood cannot be feasibly computed, SBI still recovers unbiased estimates of individual cluster masses and combined constraints from multiple clusters. SBI will be a powerful tool for constraining the masses of galaxy clusters detected by future cosmic surveys. Code to run the analyses presented here will be made publicly available.

Neural simulation-based inference of the neutron star equation of state directly from telescope spectra

Neutron stars provide a unique opportunity to study strongly interacting matter under extreme density conditions. The intricacies of matter inside neutron stars and their equation of state are not directly visible, but determine bulk properties, such as mass and radius, which affect the star's thermal X-ray emissions. However, the telescope spectra of these emissions are also affected by the stellar distance, hydrogen column, and effective surface temperature, which are not always well-constrained. Uncertainties on these nuisance parameters must be accounted for when making a robust estimation of the equation of state. In this study, we develop a novel methodology that, for the first time, can infer the full posterior distribution of both the equation of state and nuisance parameters directly from telescope observations. This method relies on the use of neural likelihood estimation, in which normalizing flows use samples of simulated telescope data to learn the likelihood of the neutron star spectra as a function of these parameters, coupled with Hamiltonian Monte Carlo methods to efficiently sample from the corresponding posterior distribution. Our approach surpasses the accuracy of previous methods, improves the interpretability of the results by providing access to the full posterior distribution, and naturally scales to a growing number of neutron star observations expected in the coming years.

Simulation-based Inference of Reionization Parameters from 3D Tomographic 21 cm Light-cone Images. II. Application of Solid Harmonic Wavelet Scattering Transform

The information regarding how the intergalactic medium is reionized by astrophysical sources is contained in the tomographic three-dimensional 21 cm images from the epoch of reionization. In Zhao et al. (“Paper I”), we demonstrated for the first time that density estimation likelihood-free inference (DELFI) can be applied efficiently to perform a Bayesian inference of the reionization parameters from the 21 cm images. Nevertheless, the 3D image data needs to be compressed into informative summaries as the input of DELFI by, e.g., a trained 3D convolutional neural network (CNN) as in Paper I (DELFI-3D CNN). Here in this paper, we introduce an alternative data compressor, the solid harmonic wavelet scattering transform (WST), which has a similar, yet fixed (i.e., no training), architecture to CNN, but we show that this approach (i.e., solid harmonic WST with DELFI) outperforms earlier analyses based on 3D 21 cm images using DELFI-3D CNN in terms of credible regions of parameters. Realistic effects, including thermal noise and residual foreground after removal, are also applied to the mock observations from the Square Kilometre Array. We show that under the same inference strategy using DELFI, the 21 cm image analysis with solid harmonic WST outperforms the 21 cm power spectrum analysis. This research serves as a proof of concept, demonstrating the potential to harness the strengths of WST and simulation-based inference to derive insights from future 21 cm light-cone image data.

An introduction to Bayesian simulation-based inference for quantum machine learning with examples

An introduction to Bayesian simulation-based inference for quantum machine learning with examples

Inferring astrophysical parameters using the 2D cylindrical power spectrum from reionisation

Abstract Enlightening our understanding of the first galaxies responsible for driving reionisation requires detecting the 21-cm signal from neutral hydrogen. Interpreting the wealth of information embedded in this signal requires Bayesian inference. Parameter inference from the 21-cm signal is primarily restricted to the spherically averaged power spectrum (1D PS) owing to its relatively straightforward derivation of an analytic likelihood function enabling traditional Monte-Carlo Markov-Chain (MCMC) approaches. However, in recent years, simulation-based inference (SBI) has become feasible which removes the necessity of having an analytic likelihood, enabling more complex summary statistics of the 21-cm signal to be used for Bayesian inference. In this work, we use SBI, specifically marginal neural ratio estimation to learn the likelihood-to-evidence ratio with Swyft, to explore parameter inference using the cylindrically averaged 2D PS. Since the 21-cm signal is anisotropic, the 2D PS should yield more constraining information compared to the 1D PS which isotropically averages the signal. For this, we consider a mock 1000 hr observation of the 21-cm signal using the SKA and compare the performance of the 2D PS relative to the 1D PS. Additionally, we explore two separate foreground mitigation strategies, perfect foreground removal and wedge avoidance. We find the 2D PS outperforms the 1D PS by improving the marginalised uncertainties on individual astrophysical parameters by up to ∼30 − 40 per cent irrespective of the foreground mitigation strategy. Primarily, these improvements stem from how the 2D PS distinguishes between the transverse, k⊥, and redshift dependent, k∥ information which enables greater sensitivity to the complex reionisation morphology.

Exploring the role of the halo mass function for inferring astrophysical parameters during reionisation

Abstract Detecting the 21-cm signal at z ≳ 6 will reveal insights into the properties of the first galaxies responsible for driving reionisation. To extract this information, we perform parameter inference with 3D simulations of the 21-cm signal embedded within a Bayesian inference pipeline. Presently, when performing inference, we must choose which sources of uncertainty to sample and which to hold fixed. Since the astrophysics of galaxies is much more uncertain than that of the underlying halo-mass function (HMF), we typically parameterise and model the former while fixing the latter. However, doing so may bias our inference of the galaxy properties. In this work, we explore the consequences of assuming an incorrect HMF and quantify the relative biases on our inferred astrophysical model parameters when considering the wrong HMF. We then relax this assumption by constructing a generalised five parameter HMF model and simultaneously recover it with our underlying astrophysical model. For this, we use 21CMFAST and perform Simulation-Based Inference using marginal neural ratio estimation to learn the likelihood-to-evidence ratio with Swyft. Using a mock 1000-hour observation of the 21-cm power spectrum from the forthcoming Square Kilometre Array, conservatively assuming foreground wedge avoidance, we find that assuming the incorrect HMF can bias the recovered astrophysical parameters by up to ∼3 − 4σ even when including independent information from observed luminosity functions. Using our generalised HMF model, although we recover our astrophysical parameters with a factor of ∼2 − 4 larger marginalised uncertainties, the constraints are unbiased, agnostic to the underlying HMF and therefore more conservative.

Kilonova Spectral Inverse Modelling with Simulation-based Inference: An Amortized Neural Posterior Estimation Analysis

Kilonovae represent a category of astrophysical transients, identifiable as the electromagnetic (EM) counterparts associated with the coalescence events of binary systems comprising neutron stars and neutron star–black hole pairs. They act as probes for heavy-element nucleosynthesis in astrophysical environments. These studies rely on an inference of the physical parameters (e.g., ejecta mass, velocity, composition) that describe kilonovae-based on EM observations. This is a complex inverse problem typically addressed with sampling-based methods such as Markov Chain Monte Carlo or nested sampling algorithms. However, repeated inferences can be computationally expensive, due to the sequential nature of these methods. This poses a significant challenge to ensuring the reliability and statistical validity of the posterior approximations and, thus, the inferred kilonova parameters themselves. We present a novel approach: simulation-based inference using simulations produced by KilonovaNet. Our method employs an ensemble of amortized neural posterior estimation (ANPE) with an embedding network to directly predict posterior distributions from simulated spectral energy distributions. We take advantage of the quasi-instantaneous inference time of ANPE to demonstrate the reliability of our posterior approximations using diagnostics tools, including coverage diagnostic and posterior predictive checks. We further test our model with real observations from AT 2017gfo, the only kilonova with multimessenger data, demonstrating agreement with previous likelihood-based methods while reducing inference time down to a few seconds. The inference results produced by ANPE appear to be conservative and reliable, paving the way for testable and more efficient kilonova parameter inference.

Simulation-based inference of developmental EEG maturation with the spectral graph model

The spectral content of macroscopic neural activity evolves throughout development, yet how this maturation relates to underlying brain network formation and dynamics remains unknown. Here, we assess the developmental maturation of electroencephalogram spectra via Bayesian model inversion of the spectral graph model, a parsimonious whole-brain model of spatiospectral neural activity derived from linearized neural field models coupled by the structural connectome. Simulation-based inference was used to estimate age-varying spectral graph model parameter posterior distributions from electroencephalogram spectra spanning the developmental period. This model-fitting approach accurately captures observed developmental electroencephalogram spectral maturation via a neurobiologically consistent progression of key neural parameters: long-range coupling, axonal conduction speed, and excitatory:inhibitory balance. These results suggest that the spectral maturation of macroscopic neural activity observed during typical development is supported by age-dependent functional adaptations in localized neural dynamics and their long-range coupling across the macroscopic structural network.

Simulation-based inference on virtual brain models of disorders

Connectome-based models, also known as virtual brain models (VBMs), have been well established in network neuroscience to investigate pathophysiological causes underlying a large range of brain diseases. The integration of an individual’s brain imaging data in VBMs has improved patient-specific predictivity, although Bayesian estimation of spatially distributed parameters remains challenging even with state-of-the-art Monte Carlo sampling. VBMs imply latent nonlinear state space models driven by noise and network input, necessitating advanced probabilistic machine learning techniques for widely applicable Bayesian estimation. Here we present simulation-based inference on VBMs (SBI-VBMs), and demonstrate that training deep neural networks on both spatio-temporal and functional features allows for accurate estimation of generative parameters in brain disorders. The systematic use of brain stimulation provides an effective remedy for the non-identifiability issue in estimating the degradation limited to smaller subset of connections. By prioritizing model structure over data, we show that the hierarchical structure in SBI-VBMs renders the inference more effective, precise and biologically plausible. This approach could broadly advance precision medicine by enabling fast and reliable prediction of patient-specific brain disorders.

Deriving the star formation histories of galaxies from spectra with simulation-based inference

High-resolution galaxy spectra encode information about the stellar populations within galaxies. The properties of the stars, such as their ages, masses, and metallicities, provide insights into the underlying physical processes that drive the growth and transformation of galaxies over cosmic time. We explore a simulation-based inference (SBI) workflow to infer from optical absorption spectra the posterior distributions of metallicities and the star formation histories (SFHs) of galaxies (i.e. the star formation rate as a function of time). We generated a dataset of synthetic spectra to train and test our model using the spectroscopic predictions of the MILES stellar population library and non-parametric SFHs. We reliably estimate the mass assembly of an integrated stellar population with well-calibrated uncertainties. Specifically, we reach a score of $0.97\,R^2$ for the time at which a given galaxy from the test set formed $50<!PCT!>$ of its stellar mass, obtaining samples of the posteriors in only $\,s. We then applied the pipeline to real observations of massive elliptical galaxies, recovering the well-known relationship between the age and the velocity dispersion, and show that the most massive galaxies ($ km/s) built up to 90<!PCT!> of their total stellar masses within $1$\,Gyr of the Big Bang. The inferred properties also agree with the state-of-the-art inversion codes, but the inference is performed up to five orders of magnitude faster. This SBI approach coupled with machine learning and applied to full spectral fitting makes it possible to address large numbers of galaxies while performing a thick sampling of the posteriors. It will allow both the deterministic trends and the inherent uncertainties of the highly degenerated inversion problem to be estimated for large and complex upcoming spectroscopic surveys, such as DESI, WEAVE, or 4MOST.

Fast Parameter Inference on Pulsar Timing Arrays with Normalizing Flows.

Pulsar timing arrays perform Bayesian posterior inference with expensive Markov chain MonteCarlo (MCMC) methods. Given a dataset of ∼10-100 pulsars and O(10^{3}) timing residuals each, producing a posterior distribution for the stochastic gravitational wave background (SGWB) can take days to a week. The computational bottleneck arises because the likelihood evaluation required for MCMC is extremely costly when considering the dimensionality of the search space. Fortunately, generating simulated data is fast, so modern simulation-based inference techniques can be brought to bear on the problem. In this Letter, we demonstrate how conditional normalizing flows trained on simulated data can be used for extremely fast and accurate estimation of the SGWB posteriors, reducing the sampling time from weeks to a matter of seconds.

Applying simulation-based inference to spectral and spatial information from the Galactic Center gamma-ray excess

The two most favored explanations of the Fermi Galactic Center gamma-ray excess (GCE) are millisecond pulsars and self annihilation of the smooth dark matter halo of the galaxy. In order to distinguish between these possibilities, we would like to optimally use all information in the available data, including photon direction and energy information.To date, analyses of the GCE have generally treated directional and energy information separately, or have ignored one or the other completely.Here, we develop a method for analyzing the GCE that relies on simulation-basedinference with neural posterior models to jointly analyze photon directional and spectral information while correctly accounting for the spatial and energy resolution of the telescope, here assumed to be the Fermi Large Area Telescope (LAT).Our results also have implications for analyses of the diffuse gamma-ray background, which we discuss.