Polynomial Stein discrepancies (PSD) provide a scalable alternative to kernel Stein methods for measuring sample quality and goodness-of-fit testing, but their statistical properties remain poorly understood. We show that increasing polynomial degree primarily amplifies signal without adequately controlling variance, rather than directly optimising the signal-to-noise ratio (SNR). Under suitable assumptions, this might lead to a failure mode in which the SNR^2 can provably decay exponentially with polynomial degree. Motivated by this observation, we reformulate Stein discrepancy construction as an explicit SNR^2 maximisation problem, yielding a Rayleigh quotient over Stein features. This perspective motivates λ-PSD, an approximate scalable covariance-aware reweighting scheme defined in a low-dimensional subspace. Under Gaussian settings, we show that λ-PSD avoids the exponential SNR^2 collapse and achieves a stable SNR^2. Empirically, λ-PSD substantially improves test power while retaining linear-time complexity in the number of samples, highlighting the importance of SNR-aware design for scalable Stein discrepancies.

Ordinary differential equation (ODE) models are widely used to describe systems in many areas of science. To ensure these models provide accurate and interpretable representations of real-world dynamics, it is often necessary to infer parameters from data, which involves specifying the form of the ODE system as well as a statistical model describing the observational process. A popular and convenient choice for the error model is a Gaussian distribution with constant variance. However, the choice may not be realistic in many systems, since the variance of the observational error may vary over time or have some dependence on the system state (heteroscedastic), reflecting changes in measurement conditions, environmental fluctuations, or intrinsic system variability. Misspecification of the error model can lead to substantial inaccuracies of the posterior estimates of the ODE model parameters and predictions. More elaborate parametric error models could be specified, but this would increase computational cost because additional parameters would need to be estimated within the MCMC procedure and may still be misspecified. In this work we propose a two-step semi-parametric framework for Bayesian parameter estimation of ODE model parameters when there exists heteroscedasticity in the error process. The first step applies a heteroscedastic Gaussian process to estimate the time-dependent error, and the second step performs Bayesian inference for the ODE model parameters using the estimated time-dependent error estimated from step one in the likelihood function. Through a simulation study and two real-world applications, we demonstrate that the proposed approach yields more reliable posterior inference and predictive uncertainty compared to the standard homoscedastic models. Although our focus is on heteroscedasticity, the framework could be applied to handle more complex error processes.

Stochastic models can be highly computationally expensive. This limits the range of parameters and scenarios that can be realistically explored. Previously, a queuing network model was developed for the insulin-stimulated intracellular translocation of the glucose transporter GLUT4. Whilst one hypothesis of insulin action was tested, alternative hypotheses were too computationally expensive for parameter inference. In this study, a deterministic surrogate model is developed for the queuing network. The surrogate model uses feedback terms in a system of differential equations to approximate the blocking mechanisms seen in the queuing network. A sensitivity analysis of the surrogate model was performed and its correspondence to the queuing network assessed. This surrogate model may be useful in a parameter inference recalibration process, allowing posteriors for the queuing network to be acquired with lower computational cost.

Mixture-of-Experts (MoE) architectures combine specialized predictors through a learned gate and are effective across regression and classification, but for classification with softmax multinomial-logistic gating, rigorous guarantees for stable maximum-likelihood training and principled model selection remain limited. We address both issues in the full-data (batch) regime. First, we derive a batch minorization-maximization (MM) algorithm for softmax-gated multinomial-logistic MoE using an explicit quadratic minorizer, yielding coordinate-wise closed-form updates that guarantee monotone ascent of the objective and global convergence to a stationary point (in the standard MM sense), avoiding approximate M-steps common in EM-type implementations. Second, we prove finite-sample rates for conditional density estimation and parameter recovery, and we adapt dendrograms of mixing measures to the classification setting to obtain a sweep-free selector of the number of experts that achieves near-parametric optimal rates after merging redundant fitted atoms. Experiments on biological protein--protein interaction prediction validate the full pipeline, delivering improved accuracy and better-calibrated probabilities than strong statistical and machine-learning baselines.

Decision-making in hierarchical systems requires probabilistic forecasts at all cross-sectional levels. Current hierarchical forecasting methods typically generate independent forecasts at each level and reconcile them post hoc to ensure coherence between upper and lower levels. Such post hoc corrections do not incorporate hierarchical structure or decision goals into the underlying parameter estimation. We propose a fully Bayesian hierarchical forecasting framework that shares information more effectively between and across levels than reconciliation alone. Our approach has the flexibility to softly penalise incoherence, subject to model specification, and to focus the global model and coherence update on hierarchical levels most relevant to decision outcomes. This yields parameter estimates that are focused towards the forecasting goals and capture the requirement for coherency, removing the need to estimate covariance matrices for multi-step forecasting horizons. We demonstrate improvements in predictive accuracy metrics on both simulated data and Australian domestic tourism forecasting.

Processing high-volume, streaming data is increasingly common in modern statistics and machine learning, where batch-mode algorithms are often impractical because they require repeated passes over the full dataset. This has motivated incremental stochastic estimation methods, including the incremental stochastic Expectation-Maximization (EM) algorithm formulated via stochastic approximation. In this work, we revisit and analyze an incremental stochastic variant of the Majorization-Minimization (MM) algorithm, which generalizes incremental stochastic EM as a special case. Our approach relaxes key EM requirements, such as explicit latent-variable representations, enabling broader applicability and greater algorithmic flexibility. We establish theoretical guarantees for the incremental stochastic MM algorithm, proving consistency in the sense that the iterates converge to a stationary point characterized by a vanishing gradient of the objective. We demonstrate these advantages on a softmax-gated mixture of experts (MoE) regression problem, for which no stochastic EM algorithm is available. Empirically, our method consistently outperforms widely used stochastic optimizers, including stochastic gradient descent, root mean square propagation, adaptive moment estimation, and second-order clipped stochastic optimization. These results support the development of new incremental stochastic algorithms, given the central role of softmax-gated MoE architectures in contemporary deep neural networks for heterogeneous data modeling. Beyond synthetic experiments, we also validate practical effectiveness on two real-world datasets, including a bioinformatics study of dent maize genotypes under drought stress that integrates high-dimensional proteomics with ecophysiological traits, where incremental stochastic MM yields stable gains in predictive performance.

Background Health care systems face growing fiscal pressure while AI reaches clinical parity in several domains. UK National Health Service expenditure rose by 52%, while Australia's health expenditure grew by 29% between 2019 and 2023. Yet large-scale cost savings from AI remain limited, largely because implementation constraints continue to outweigh technical capability. Objective This study develops a Bayesian budget-impact framework to estimate AI-driven gross cost savings in radiology, workflow optimization, and workforce optimization in the United Kingdom and Australia, explicitly accounting for adoption uncertainty, effectiveness, and implementation risk. Methods We used a sequential Monte Carlo simulation with 1000 particles to estimate annual gross cost savings. The core savings function combined expenditure base, sector weight, adoption, effectiveness, and implementation risk. Priors were informed by a structured review of multiple studies per domain. The savings likelihood used a heteroscedastic noise specification in which the SD followed an exponential prior, with the mean set to 15% of each sector’s observed savings estimate, ranging from US $12.0 million for Australian radiology to US $87.2 million for UK workforce optimization. The likelihood was also augmented with sector-specific observations that anchored effectiveness to published cost-reduction estimates and implementation risk to observed deployment failure rates. Scenario analysis applied multipliers to 1000 bootstrap posterior draws across optimistic, conservative, and pessimistic settings. Sensitivity analyses varied the σ scaling factor from 0.10 to 0.20 and perturbed prior means for adoption, effectiveness, and implementation risk by ±20%. Results Posterior annual savings were US $949 million (95% credible interval [CrI] US $720.6-US $1173.5 million) for the United Kingdom and US $737 million (95% CrI US $526.1-US $953.6 million) for Australia. Workforce optimization generated the largest share of savings in both countries, contributing 62.4% of UK savings (US $591.9 million, 95% CrI US $398.8-US $854.8 million) and 76.2% of Australian savings (US $561.4 million, 95% CrI US $331.9-US $815.3 million). Posterior implementation risk estimates ranged from 35.5% to 47.4%, below prior means of 49.1% to 57.2%, reflecting the empirical anchoring introduced through deployment-failure data. Across scenarios, projected savings ranged from US $357 million to US $1.845 billion in the UK and from US $267.2 million to US $1.454 billion in Australia. Baseline cumulative projections for 2024-2030 were US $10.1 billion for the United Kingdom and US $8.0 billion for Australia. The sensitivity analyses confirmed the robustness of the posterior savings estimates. Conclusions AI could generate substantial expenditure reductions in both health systems, but implementation risk remains the main constraint on realizing those gains. Workforce savings should be interpreted primarily as capacity gains that can be redeployed to higher-value care, not as automatic cash savings. The Bayesian framework offers probabilistic planning ranges rather than point forecasts and provides a practical basis for policy planning under uncertainty.

When effective vaccines are available, vaccination programs are typically one of the best defences against the spread of an infectious disease. Such vaccination programs become particularly important during severe epidemics or pandemics to ensure sufficient vaccination coverage is achieved to increase protection or reduce transmission to a level that enables relaxation of non-pharmaceutical interventions, such as lockdowns, travel restrictions, or social distancing. Unfortunately, vaccination uptake in the community may be slow if there are substantial levels of vaccine hesitancy in the population. As a result, it is important to identify when these hesitancy behaviours are present in the community. Furthermore, understanding the main drivers of such behaviour can inform adjustments to public health strategies to improve community uptake. In this study, we consider the problem of identifying vaccination hesitancy behaviour during a vaccination roll-out that occurs in response to a severe epidemic. Specifically, our aim is to explore the extent to which mathematical modelling of reported case, death, and vaccination counts can be used to detect vaccine hesitancy and possible drivers. To do this, we develop a novel susceptible-exposed-infectious-recovered (SEIR) epidemiological model of disease transmission that incorporates changes in population behaviour relating to non-pharmaceutical interventions and vaccine uptake that are influenced by information reported through media or data dashboards about cases, deaths, and vaccination rates. We then use a Bayesian approach to analyse simulated data representing various hesitancy scenarios. Through this simulation study, our key findings are that individual parameters values related to drivers of vaccine hesitancy often cannot be identified. However, posterior correlation structures between these parameters enable the presence of vaccine hesitancy in the community to be detected and provide some insight into the relative influence of key factors, such as vaccine safety concerns or complacency. While our simulation study is inspired by the public health response to the COVID-19 pandemic, our tools and techniques are general and could enable vaccination programs of various infectious diseases to be adapted rapidly in response to community behaviours in the future.

Abstract Motivation Kinetic models are central to systems biology, but enzyme-kinetic parameters compiled from the literature and databases are often incomplete, inconsistent, and measured under heterogeneous conditions. Classical parameter balancing helps infer missing parameters, yet it often lacks calibrated uncertainty, robustness to misspecification, and explicit treatment of source-level heterogeneity. Results We develop a formal Bayesian parameter balancing framework that enforces thermodynamic constraints, estimates full posterior uncertainty, and validates calibration using leave-one-out cross-validation and posterior-predictive coverage. Beyond the classical Gaussian formulation, we introduce robust Student- t and skewed error models to improve reliability under outliers and model misspecification, and incorporate random effects to account for source-level or group-level variability across studies. The resulting approach yields thermodynamically consistent parameter sets with well-calibrated credible intervals on held-out data, offering a Bayesian parameter balancing approach useful to systems biology researchers. Availability and implementation Source code, data, workflows, a Julia package and command-line usage are available at the project GitHub repository . Graphical Abstract

Likelihood-free inference (LFI) methods, such as approximate Bayesian computation, have become commonplace for conducting inference in complex models. Many approaches are based on summary statistics or discrepancies derived from synthetic data. However, determining which summary statistics or discrepancies to use for constructing the posterior remains a challenging question, both practically and theoretically. Instead of relying on a single vector of summaries for inference, we propose a new pooled posterior that optimally combines inferences from multiple LFI posteriors. This pooled approach eliminates the need to select a single vector of summaries or even a specific LFI algorithm. Our approach is straightforward to implement and avoids performing a high-dimensional LFI analysis involving all summary statistics. We give theoretical guarantees for the improved performance of the pooled posterior mean in terms of asymptotic frequentist risk and demonstrate the effectiveness of the approach in a number of benchmark examples.

Forecasting ecosystem changes due to disturbances or conservation interventions is essential to improve ecosystem management and anticipate unintended consequences of conservation decisions. Mathematical models allow practitioners to understand the potential effects and unintended consequences via simulation. However, calibrating these models is often challenging due to a paucity of appropriate ecological data. Ensemble ecosystem modelling (EEM) is a quantitative method used to parameterize models from theoretical ecosystem features rather than data. Two approaches have been considered to find parameter values satisfying those features: a standard accept–reject algorithm, appropriate for small ecosystem networks, and a sequential Monte Carlo (SMC) algorithm that is more computationally efficient for larger ecosystem networks. In practice, using SMC for EEM generation requires advanced statistical and mathematical knowledge, as well as strong programming skills, which might limit its uptake. In addition, current EEM approaches have been developed for only one model structure (generalised Lotka–Volterra). To facilitate the usage of EEM methods, we introduce EEMtoolbox, an R package for calibrating quantitative ecosystem models. Our package allows the generation of parameter sets satisfying ecosystem features by using either the standard accept–reject algorithm or the novel SMC procedure. Our package extends the existing EEM methodology, originally developed for the generalised Lotka–Volterra model, to two additional model structures (the multispecies Gompertz and the Bimler–Baker model) and additionally allows users to define their own model structures. We demonstrate the usage of EEMtoolbox by simulating changes in species abundance immediately after the release of the sihek ( Todiramphus cinnamominus , extinct‐in‐the‐wild species) on Palmyra Atoll in the Pacific Ocean. With its simple interface, our package facilitates straightforward generation of EEM parameter sets, thus unlocking advanced statistical methods supporting conservation decisions using ecosystem network models.

Mixture of Experts (MoE) models constitute a widely utilized class of ensemble learning approaches in statistics and machine learning, known for their flexibility and computational efficiency. They have become integral components in numerous state-of-the-art deep neural network architectures, particularly for analyzing heterogeneous data across diverse domains. Despite their practical success, the theoretical understanding of model selection, especially concerning the optimal number of mixture components or experts, remains limited and poses significant challenges. These challenges primarily stem from the inclusion of covariates in both the Gaussian gating functions and expert networks, which introduces intrinsic interactions governed by partial differential equations with respect to their parameters. In this paper, we revisit the concept of dendrograms of mixing measures and introduce a novel extension to Gaussian-gated Gaussian MoE models that enables consistent estimation of the true number of mixture components and achieves the pointwise optimal convergence rate for parameter estimation in overfitted scenarios. Notably, this approach circumvents the need to train and compare a range of models with varying numbers of components, thereby alleviating the computational burden, particularly in high-dimensional or deep neural network settings. Experimental results on synthetic data demonstrate the effectiveness of the proposed method in accurately recovering the number of experts. It outperforms common criteria such as the Akaike information criterion, the Bayesian information criterion, and the integrated completed likelihood, while achieving optimal convergence rates for parameter estimation and accurately approximating the regression function.

Mathematical models connect theory with the real world through data, enabling us to interpret, understand, and predict complex phenomena. However, scientific knowledge often extends beyond what can be empirically measured, offering valuable insights into complex and uncertain systems. Here, we introduce a statistical framework for calibrating mathematical models using non-empirical information. Through examples in ecology, biology, and medicine, we demonstrate how expert knowledge, scientific theory, and qualitative observations can meaningfully constrain models. In each case, these non-empirical insights guide models toward more realistic dynamics and more informed predictions than empirical data alone could achieve. Now, our understanding of the systems represented by mathematical models is not limited by the data that can be obtained; they instead sit at the edge of scientific understanding.

Predictions of animal movement are vital for understanding and managing wild populations. However, the fine-scale, complex decision-making of animals can pose challenges for the accurate prediction of trajectories. Integrated step selection functions (iSSFs), a common tool for inferring relationships between animal movement and the environment, are also increasingly used to simulate animal trajectories for prediction. Although admitting a lot of flexibility, the iSSF framework is limited to its reliance on pre-defined functional forms for fitting to data, and iSSFs that involve complex functional forms to model detailed processes can be prohibitively difficult to fit and interpret. Here, we present deepSSF, an approach to fit and predict animal movement data using deep learning. The deepSSF approach replaces the log-linear model of an iSSF with a neural network architecture that receives multiple environmental layers and scalar values as inputs and outputs a single layer representing the next-step probability. We demonstrate an example deepSSF model, built in PyTorch, consisting of distinct but interacting habitat selection and movement subnetworks. This allows for explicit representation of both selection and movement processes, thus giving interpretable intermediate outputs. We apply our model to GPS data of introduced water buffalo (Bubalus bubalis) in the tropical savannas of Northern Australia. Our deepSSF model was able to learn features that are present in the habitat covariate layers, such as linear features (rivers, forest edges) and the composition of certain habitat areas, without having to specify them pre-emptively within the model framework. It was able to capture complex interactions between the habitat covariates as well as temporal dynamics across time of day and year. Finally, our deepSSF model generally had better in- and out-of-sample predictive accuracy than the analogous iSSF model. We expect that the deepSSF approach will generate accurate and informative predictions about animal movement, which can be used for deepening our understanding of animal-environment systems and for the practical management of species. We discuss how the wide range of existing deep learning tools could enable the deepSSF approach to be extended to represent memory and social dynamic processes, with the potential for integrating non-spatial data sources such as accelerometers and physiological sensors.

Approximate Bayesian computation (ABC) is one of the most popular "likelihood-free" methods. These methods have been applied in a wide range of fields by providing solutions to intractable likelihood problems in which exact Bayesian approaches are either infeasible or computationally costly. However, the performance of ABC can be unreliable when dealing with model misspecification. To circumvent the poor behavior of ABC in these settings, we propose a novel ABC approach that is robust to model misspecification. This new method can deliver more accurate statistical inference under model misspecification than alternatives and also enables the detection of summary statistics that are incompatible with the assumed data-generating process. We demonstrate the effectiveness of our approach through several simulated examples, where it delivers more accurate point estimates and uncertainty quantification over standard ABC approaches when the model is misspecified. Additionally, we apply our approach to an empirical example, further showcasing its advantages over alternative methods.

Simulation-based Bayesian inference (SBI) methods are widely used for parameter estimation in complex models where evaluating the likelihood is challenging but generating simulations is relatively straightforward. However, these methods commonly assume that the simulation model accurately reflects the true data-generating process, an assumption that is frequently violated in realistic scenarios. In this paper, we focus on the challenges faced by SBI methods under model misspecification. We consolidate recent research aimed at mitigating the effects of misspecification, highlighting three key strategies: i) robust summary statistics, ii) generalised Bayesian inference, and iii) error modelling and adjustment parameters. To illustrate both the vulnerabilities of popular SBI methods and the effectiveness of misspecification-robust alternatives, we present empirical results on an illustrative example.

We develop a unified statistical framework for softmax-gated Gaussian mixture of experts (SGMoE) that addresses three long-standing obstacles in parameter estimation and model selection: (i) non-identifiability of gating parameters up to common translations, (ii) intrinsic gate-expert interactions that induce coupled differential relations in the likelihood, and (iii) the tight numerator-denominator coupling in the softmax-induced conditional density. Our approach introduces Voronoi-type loss functions aligned with the gate-partition geometry and establishes finite-sample convergence rates for the maximum likelihood estimator (MLE). In over-specified models, we reveal a link between the MLE's convergence rate and the solvability of an associated system of polynomial equations characterizing near-nonidentifiable directions. For model selection, we adapt dendrograms of mixing measures to SGMoE, yielding a consistent, sweep-free selector of the number of experts that attains pointwise-optimal parameter rates under overfitting while avoiding multi-size training. Simulations on synthetic data corroborate the theory, accurately recovering the expert count and achieving the predicted rates for parameter estimation while closely approximating the regression function. Under model misspecification (e.g., ε-contamination), the dendrogram selection criterion is robust, recovering the true number of mixture components, while the Akaike information criterion, the Bayesian information criterion, and the integrated completed likelihood tend to overselect as sample size grows. On a maize proteomics dataset of drought-responsive traits, our dendrogram-guided SGMoE selects two experts, exposes a clear mixing-measure hierarchy, stabilizes the likelihood early, and yields interpretable genotype-phenotype maps, outperforming standard criteria without multi-size training.

Updating a priori information given some observed data is the core tenet of Bayesian inference. Bayesian transfer learning extends this idea by incorporating information from a related dataset to improve the inference on the observed target dataset which may have been collected under slightly different settings. The use of related information can be useful when the target dataset is scarce, for example. There exist various Bayesian transfer learning methods that decide how to incorporate the related data in different ways. Unfortunately, there is no principled approach for comparing Bayesian transfer methods in real data settings. Additionally, some Bayesian transfer learning methods, such as the so-called power prior approaches, rely on conjugacy or costly specialised techniques. In this paper, we find an effective approach to compare Bayesian transfer learning methods is to apply leave-one-out cross validation on the target dataset. Further, we introduce a new framework, transfer sequential Monte Carlo, that efficiently implements power prior methods in an automated fashion. We demonstrate the performance of our proposed methods in two comprehensive simulation studies.

The Bayesian Synthetic Likelihood (BSL) method is a widely-used tool for likelihood-free Bayesian inference. This method assumes that some summary statistics are normally distributed, which can be incorrect in many applications. We propose a transformation, called the Wasserstein Gaussianization transformation, that uses a Wasserstein gradient flow to approximately transform the distribution of the summary statistics into a Gaussian distribution. BSL also implicitly requires compatibility between simulated summary statistics under the working model and the observed summary statistics. A robust BSL variant which achieves this has been developed in the recent literature. We combine the Wasserstein Gaussianization transformation with robust BSL, and an efficient Variational Bayes procedure for posterior approximation, to develop a highly efficient and reliable approximate Bayesian inference method for likelihood-free problems.