Search

Stochastic simulation models are generative models that mimic complex systems to help with decision-making. The reliability of these models heavily depends on well-calibrated input model parameters. However, in many practical scenarios, only output-level data are available to learn the input model parameters, which is challenging due to the often intractable likelihood of the stochastic simulation model. Moreover, stochastic simulation models are frequently inexact, with discrepancies between the model and the target system. No existing methods can effectively learn and quantify the uncertainties of input parameters using only output-level data. In this paper, we propose to learn differentiable input parameters of stochastic simulation models using output-level data via kernel score minimization with stochastic gradient descent. We quantify the uncertainties of the learned input parameters using a frequentist confidence set procedure based on a new asymptotic normality result that accounts for model inexactness. The proposed method is evaluated on exact and inexact G/G/1 queueing models., Comment: 31 pages, 12 tables, 4 figures

Agent-based models capture heterogeneity among individuals in a population and are widely used in studies of multi-cellular systems, disease, epidemics and demography to name a few. However, existing frameworks consider discrete time-step simulation or assume that agents' states only change as a result of discrete events. In this note, we present AgentBasedModeling$.$jl, a Julia package for simulating stochastic agent-based population models in continuous time. The tool allows to easily specify and simulate agents evolving through generic continuous-time jump-diffusions and interacting via continuous-rate processes. AgentBasedModeling$.$jl provides a powerful methodology for studying the effects of stochasticity on structured population dynamics., Comment: 6 pages, 1 figure

Reproducibility is a foundational standard for validating scientific claims in computational research. Stochastic computational models are employed across diverse fields such as systems biology, financial modelling and environmental sciences. Existing infrastructure and software tools support various aspects of reproducible model development, application, and dissemination, but do not adequately address independently reproducing simulation results that form the basis of scientific conclusions. To bridge this gap, we introduce the Empirical Characteristic Function Equality Convergence Test (EFECT), a data-driven method to quantify the reproducibility of stochastic simulation results. EFECT employs empirical characteristic functions to compare reported results with those independently generated by assessing distributional inequality, termed EFECT error, a metric to quantify the likelihood of equality. Additionally, we establish the EFECT convergence point, a metric for determining the required number of simulation runs to achieve an EFECT error value of a priori statistical significance, setting a reproducibility benchmark. EFECT supports all real-valued and bounded results irrespective of the model or method that produced them, and accommodates stochasticity from intrinsic model variability and random sampling of model inputs. We tested EFECT with stochastic differential equations, agent-based models, and Boolean networks, demonstrating its broad applicability and effectiveness. EFECT standardizes stochastic simulation reproducibility, establishing a workflow that guarantees reliable results, supporting a wide range of stakeholders, and thereby enhancing validation of stochastic simulation studies, across a model's lifecycle. To promote future standardization efforts, we are developing open source software library libSSR in diverse programming languages for easy integration of EFECT., Comment: 25 pages, 4 figures

In the infectious disease literature, significant effort has been devoted to studying dynamics at a single scale. For example, compartmental models describing population-level dynamics are often formulated using differential equations. In cases where small numbers or noise play a crucial role, these differential equations are replaced with memoryless Markovian models, where discrete individuals can be members of a compartment and transition stochastically. Classic stochastic simulation algorithms, such as Gillespie's algorithm and the next reaction method, can be employed to solve these Markovian models exactly. The intricate coupling between models at different scales underscores the importance of multiscale modelling in infectious diseases. However, several computational challenges arise when the multiscale model becomes non-Markovian. In this paper, we address these challenges by developing a novel exact stochastic simulation algorithm. We apply it to a showcase multiscale system where all individuals share the same deterministic within-host model while the population-level dynamics are governed by a stochastic formulation. We demonstrate that as long as the within-host information is harvested at a reasonable resolution, the novel algorithm we develop will always be accurate. Moreover, the novel algorithm we develop is general and can be easily applied to other multiscale models in (or outside) the realm of infectious diseases., Comment: 23 pages, 7 figures

Risk management is particularly concerned with extreme events, but analysing these events is often hindered by the scarcity of data, especially in a multivariate context. This data scarcity complicates risk management efforts. Various tools can assess the risk posed by extreme events, even under extraordinary circumstances. This paper studies the evaluation of univariate risk for a given risk factor using metrics that account for its asymptotic dependence on other risk factors. Data availability is crucial, particularly for extreme events where it is often limited by the nature of the phenomenon itself, making estimation challenging. To address this issue, two non-parametric simulation algorithms based on multivariate extreme theory are developed. These algorithms aim to extend a sample of extremes jointly and conditionally for asymptotically dependent variables using stochastic simulation and multivariate Generalised Pareto Distributions. The approach is illustrated with numerical analyses of both simulated and real data to assess the accuracy of extreme risk metric estimations.

The effective planning and allocation of resources in modern breeding programs is a complex task. Breeding program design and operational management have a major impact on the success of a breeding program and changing parameters such as the number of selected/phenotyped/genotyped individuals will impact genetic gain, genetic diversity, and costs. As a result, careful assessment and balancing of design parameters is crucial, considering the trade-offs between different breeding goals and associated costs. In a previous study, we optimized the resource allocation strategy in a dairy cattle breeding scheme via the combination of stochastic simulations and kernel regression, aiming to maximize a target function containing genetic gain and the inbreeding rate under a given budget. However, the high number of simulations required when using the proposed kernel regression method to optimize a breeding program with many parameters weakens the effectiveness of such a method. In this work, we are proposing an optimization framework that builds on the concepts of kernel regression but additionally makes use of an evolutionary algorithm to allow for a more effective and general optimization. The key idea is to consider a set of potential parameterizations of the breeding program, evaluate their performance based on stochastic simulations, and use these outputs to derive new parametrization to test in an iterative procedure. The evolutionary algorithm was implemented in a Snakemake pipeline to allow for efficient scaling on large distributed computing platforms. The algorithm achieved convergence to the same optimum with a massively reduced number of simulations. Thereby, the incorporation of class variables and accounting for a higher number of parameters in the optimization pipeline leads to substantially reduced computing time and better scaling for the desired optimization of a breeding program.

Generic open quantum systems are notoriously difficult to simulate unless one looks at specific regimes. In contrast, classical dissipative systems can often be effectively described by stochastic processes, which are generally less computationally expensive. Here, we use the paradigmatic case of a dissipative quantum oscillator to give a pedagogic introduction into the modelling of open quantum systems using quasiclassical methods, i.e. classical stochastic methods that use a 'quantum' noise spectrum to capture the influence of the environment on the system. Such quasiclassical methods have the potential to offer insights into the impact of the quantum nature of the environment on the dynamics of the system of interest whilst still being computationally tractable., Comment: Comments welcome

Although Gillespie's algorithm is justified under a set of axioms based on the assumption of homogeneity of the system, many chemical systems deviate from this assumption, as is the case for reactions taking place in low-mobility media. Using instead the generalized q formalism, we propose a new stochastic simulation algorithm by redefining the probability with which a mu reaction occurs between t + tao and t + tao + delta tao as P(tao q, mu) = a mu exp q(-a0 tao q), taking into account the separation of the natural exponential by the q parameter. Our algorithm has been implemented in the study of binary annihilation reactions, demonstrating a wider amplitude within the range of established physicochemical reactions, being the stochastic Gillespie scheme and the classical deterministic approach, particular cases of this new proposal. The effect of the nonextensivity parameter, q, on the reaction rate is analyzed and its relationship with the reaction order, n, and the heterogeneity parameter, h, is determined for two different reactant concentrations in the annihilation reaction. Different behaviors of these parameters are observed for the two types of samples, especially as q moves away from 1, confirming that quasi-second order reactions occur when reactant concentrations are similar and quasi-first order reactions when they are different. Empirical expressions between the classical reaction order and the degree of heterogeneity are proposed. The results obtained allow us to associate the behavior of sub and supra Arrhenius kinetics reported in other different reactions with the degree of non-extensiveness of the reaction., Comment: 20 pages, 10 figures, physics.chem-ph (Chemical Physics)

Our study focuses on fractional order compartment models derived from underlying physical stochastic processes, providing a more physically grounded approach compared to models that use the dynamical system approach by simply replacing integer-order derivatives with fractional order derivatives. In these models, inherent stochasticity becomes important, particularly when dealing with the dynamics of small populations far from the continuum limit of large particle numbers. The necessity for stochastic simulations arises from deviations of the mean states from those obtained from the governing equations in these scenarios. To address this, we introduce an exact stochastic simulation algorithm designed for fractional order compartment models, based on a semi-Markov process. We have considered a fractional order resusceptibility SIS model and a fractional order recovery SIR model as illustrative examples, highlighting significant disparities between deterministic and stochastic dynamics when the total population is small. Beyond its modeling applications, the algorithm presented serves as a versatile tool for solving fractional order differential equations via Monte Carlo simulations., Comment: 26 pages, 7 figures

We provide finite-sample performance guarantees for control policies executed on stochastic robotic systems. Given an open- or closed-loop policy and a finite set of trajectory rollouts under the policy, we bound the expected value, value-at-risk, and conditional-value-at-risk of the trajectory cost, and the probability of failure in a sparse cost setting. The bounds hold, with user-specified probability, for any policy synthesis technique and can be seen as a post-design safety certification. Generating the bounds only requires sampling simulation rollouts, without assumptions on the distribution or complexity of the underlying stochastic system. We adapt these bounds to also give a constraint satisfaction test to verify safety of the robot system. We provide a thorough analysis of the bound sensitivity to sim-to-real distribution shifts and provide results for constructing robust bounds that can tolerate some specified amount of distribution shift. Furthermore, we extend our method to apply when selecting the best policy from a set of candidates, requiring a multi-hypothesis correction. We show the statistical validity of our bounds in the Ant, Half-cheetah, and Swimmer MuJoCo environments and demonstrate our constraint satisfaction test with the Ant. Finally, using the 20 degree-of-freedom MuJoCo Shadow Hand, we show the necessity of the multi-hypothesis correction., Comment: Submitted to IEEE-TRO

Gene regulatory networks, as a powerful abstraction for describing complex biological interactions between genes through their expression products within a cell, are often regarded as virtually deterministic dynamical systems. However, this view is now being challenged by the fundamentally stochastic, 'bursty' nature of gene expression revealed at the single cell level. We present a Python package called Harissa which is dedicated to simulation and inference of such networks, based upon an underlying stochastic dynamical model driven by the transcriptional bursting phenomenon. As part of this tool, network inference can be interpreted as a calibration procedure for a mechanistic model: once calibrated, the model is able to capture the typical variability of single-cell data without requiring ad hoc external noise, unlike ordinary or even stochastic differential equations frequently used in this context. Therefore, Harissa can be used both as an inference tool, to reconstruct biologically relevant networks from time-course scRNA-seq data, and as a simulation tool, to generate quantitative gene expression profiles in a non-trivial way through gene interactions.

There is a large professional literature on the correct measurement of the funded status of and indicated employer contributions to government employee pension plans. But static measures do not provide a quantification of the risk that plans could represent in the future in various possible investment environments. This is better done through stochastic simulation projections. Using a 2022 data base of 187 large pension plans, basic plan features and conditions, actuarial relationships, simple economic projections, and varied bond and stock investment returns based on past ten-year historical periods in the US and the current asset allocation of the plans, I create a risk model. This risk assessment is denominated as the range of reasonably possible future actuarial funded ratios and amounts of indicated employer contributions across plans, when investment performance is poor but within historical experience. Several plans would be at or near insolvency, about 13 percent of plans would see very large drops in funded status, while employer contributions would, on average, double, and for about one in eight plans, triple. Highlights * At the median of historical investment experience, the mean plan funded ratio is projected to increase from about 75 percent in 2022 to around 82 percent in 2032. But at the 25th and 10th percentile simulations, the average plan funded ratio is 66 and 56 percent, respectively, in 2032. * At the fiftieth percentile, employer contributions decline by $17 billion in aggregate from 2022 to 2032 or 13 percent. But if investment conditions do not turn out so well, at the 75th percentile (the obverse of the 25th percentile for funded ratios), aggregate contributions increase by $73 billion, or 57 percent, and at the 90th percentile, in a rough investment (and presumably economic) environment, aggregate contributions nearly double, increasing by $127 billion. * Another measure of risk--the ratio of the standard deviation to the mean of indicated contributions projected in 2032--is high, at 0.81, for the median plan, and above 0.97 for about a quarter of the plans. Keywords and JEL Codes Public pension plans, funding, costs, investment risk, stochastic projections G17, G23, H75, J32, Introduction There is a large professional literature on the correct measurement of the funded status of or indicated employer contributions to defined benefit plans, both in the private and public [...]