Monte Carlo Estimation Research Articles

Repeating caldera collapse events constrain fault friction at the kilometer scale

Fault friction is central to understanding earthquakes, yet laboratory rock mechanics experiments are restricted to, at most, meter scale. Questions thus remain as to the applicability of measured frictional properties to faulting in situ. In particular, the slip-weakening distance [Formula: see text] strongly influences precursory slip during earthquake nucleation, but scales with fault roughness and is challenging to extrapolate to nature. The 2018 eruption of K̄ılauea volcano, Hawaii, caused 62 repeatable collapse events in which the summit caldera dropped several meters, accompanied by [Formula: see text] 4.7 to 5.4 very long period (VLP) earthquakes. Collapses were exceptionally well recorded by global positioning system (GPS) and tilt instruments and represent unique natural kilometer-scale friction experiments. We model a piston collapsing into a magma reservoir. Pressure at the piston base and shear stress on its margin, governed by rate and state friction, balance its weight. Downward motion of the piston compresses the underlying magma, driving flow to the eruption. Monte Carlo estimation of unknowns validates laboratory friction parameters at the kilometer scale, including the magnitude of steady-state velocity weakening. The absence of accelerating precollapse deformation constrains [Formula: see text] to be [Formula: see text] mm, potentially much less. These results support the use of laboratory friction laws and parameters for modeling earthquakes. We identify initial conditions and material and magma-system parameters that lead to episodic caldera collapse, revealing that small differences in eruptive vent elevation can lead to major differences in eruption volume and duration. Most historical basaltic caldera collapses were, at least partly, episodic, implying that the conditions for stick-slip derived here are commonly met in nature.

Path-space differentiable rendering of participating media

Physics-based differentiable rendering---which focuses on estimating derivatives of radiometric detector responses with respect to arbitrary scene parameters---has a diverse array of applications from solving analysis-by-synthesis problems to training machine-learning pipelines incorporating forward-rendering processes. Unfortunately, existing general-purpose differentiable rendering techniques lack either the generality to handle volumetric light transport or the flexibility to devise Monte Carlo estimators capable of handling complex geometries and light transport effects. In this paper, we bridge this gap by showing how generalized path integrals can be differentiated with respect to arbitrary scene parameters. Specifically, we establish the mathematical formulation of generalized differential path integrals that capture both interfacial and volumetric light transport. Our formulation allows the development of advanced differentiable rendering algorithms capable of efficiently handling challenging geometric discontinuities and light transport phenomena such as volumetric caustics. We validate our method by comparing our derivative estimates to those generated using the finite differences. Further, to demonstrate the effectiveness of our technique, we compare both differentiable rendering and inverse rendering performance with state-of-the-art methods.

Monte Carlo estimators for differential light transport

Physically based differentiable rendering algorithms propagate derivatives through realistic light transport simulations and have applications in diverse areas including inverse reconstruction and machine learning. Recent progress has led to unbiased methods that can simultaneously compute derivatives with respect to millions of parameters. At the same time, elementary properties of these methods remain poorly understood. Current algorithms for differentiable rendering are constructed by mechanically differentiating a given primal algorithm. While convenient, such an approach is simplistic because it leaves no room for improvement. Differentiation produces major changes in the integrals that occur throughout the rendering process, which indicates that the primal and differential algorithms should be decoupled so that the latter can suitably adapt. This leads to a large space of possibilities: consider that even the most basic Monte Carlo path tracer already involves several design choices concerning the techniques for sampling materials and emitters, and their combination, e.g. via multiple importance sampling (MIS). Differentiation causes a veritable explosion of this decision tree: should we differentiate only the estimator, or also the sampling technique? Should MIS be applied before or after differentiation? Are specialized derivative sampling strategies of any use? How should visibility-related discontinuities be handled when millions of parameters are differentiated simultaneously? In this paper, we provide a taxonomy and analysis of different estimators for differential light transport to provide intuition about these and related questions.

Spatial modelling of Lexis mortality data

In this work we present a spatial approach to investigate mortality data referenced over a Lexis structure. We decompose the force of mortality into two interpretable components: a Markov random field, homogeneous with respect to age, period and cohort which explains the main pattern of mortality; and a secondary component of independent shocks, accounting for additional non structured patterns. Inference is based on a hierarchical Bayesian model with Markov chain Monte Carlo computations. We present an application to data from the Human Mortality Database, with respect to Italy and Sweden, two countries with very different histories in terms of demographic and epidemiological transitions. For each country the primary surface of spatially structured mortality and the secondary surface of additional mortality are estimated. The importance of each component is evaluated by the estimated value of the respective precision parameter. In both Italy and Sweden, we discovered an interesting band of extra mortality in the secondary surface across the time domain, in the age interval between 60 and 90 years, with a slightly positive slope. The band represents a significant amount of extra mortality for the elderly population, which is otherwise incompatible with a structured dynamics in age, period and cohort.

A Neural Network Monte Carlo Approximation for Expected Utility Theory

This paper proposes an approximation method to create an optimal continuous-time portfolio strategy based on a combination of neural networks and Monte Carlo, named NNMC. This work is motivated by the increasing complexity of continuous-time models and stylized facts reported in the literature. We work within expected utility theory for portfolio selection with constant relative risk aversion utility. The method extends a recursive polynomial exponential approximation framework by adopting neural networks to fit the portfolio value function. We developed two network architectures and explored several activation functions. The methodology was applied on four settings: a 4/2 stochastic volatility (SV) model with two types of market price of risk, a 4/2 model with jumps, and an Ornstein–Uhlenbeck 4/2 model. In only one case, the closed-form solution was available, which helps for comparisons. We report the accuracy of the various settings in terms of optimal strategy, portfolio performance and computational efficiency, highlighting the potential of NNMC to tackle complex dynamic models.

Multilevel and Quasi Monte Carlo Methods for the Calculation of the Expected Value of Partial Perfect Information.

The expected value of partial perfect information (EVPPI) provides an upper bound on the value of collecting further evidence on a set of inputs to a cost-effectiveness decision model. Standard Monte Carlo estimation of EVPPI is computationally expensive as it requires nested simulation. Alternatives based on regression approximations to the model have been developed but are not practicable when the number of uncertain parameters of interest is large and when parameter estimates are highly correlated. The error associated with the regression approximation is difficult to determine, while MC allows the bias and precision to be controlled. In this article, we explore the potential of quasi Monte Carlo (QMC) and multilevel Monte Carlo (MLMC) estimation to reduce the computational cost of estimating EVPPI by reducing the variance compared with MC while preserving accuracy. We also develop methods to apply QMC and MLMC to EVPPI, addressing particular challenges that arise where Markov chain Monte Carlo (MCMC) has been used to estimate input parameter distributions. We illustrate the methods using 2 examples: a simplified decision tree model for treatments for depression and a complex Markov model for treatments to prevent stroke in atrial fibrillation, both of which use MCMC inputs. We compare the performance of QMC and MLMC with MC and the approximation techniques of generalized additive model (GAM) regression, Gaussian process (GP) regression, and integrated nested Laplace approximations (INLA-GP). We found QMC and MLMC to offer substantial computational savings when parameter sets are large and correlated and when the EVPPI is large. We also found that GP and INLA-GP were biased in those situations, whereas GAM cannot estimate EVPPI for large parameter sets.

Intelligent sensitivity analysis framework based on machine learning for spacecraft thermal design

The variance-based method of global sensitivity analysis (GSA) has been widely applied in spacecraft thermal design, which is typically calculated using Monte Carlo estimations. However, such estimations require a large number of samples to ensure sufficient accuracy, which makes GSA expensive to perform when modeling is difficult. Moreover, multimodal or highly skewed output distributions may result in the use of variance as an uncertain agent that generates contradictory results. Therefore, an intelligent density-based GSA framework based on machine learning and multi-fidelity metamodels called IDGSA-3M is proposed. An intelligent batch processing system based on a real-time data interaction between MATLAB and NX/TMG was designed that uses many cheap low-fidelity sample points to reduce the cost of model evaluation while using a small number of expensive high-fidelity sample points to maintain high accuracy, thus achieving trade-offs between high accuracy and low computational cost. A radial basis function (RBF) neural network based on an improved mind evolutionary algorithm was applied to approximate the multi-fidelity metamodel of a spacecraft thermophysical model calculated using a batch processing system, which had a computational speed that was 1000+ times faster than that of the traditional thermophysical model and a high computational accuracy of 99%+. The output distributions of the RBF were then characterized by its cumulative distribution functions to obtain density-based sensitivity indices. Both the theoretical and experimental results of GSA for the thermal design parameters of the extreme ultraviolet radiation detector on the space-based Lyman-Alpha Solar Telescope, developed in China, demonstrated that the convergence rate of IDGSA-3M can be improved up to 10-fold for a fixed convergence level in comparison with two other GSA methods, thereby verifying its superiority.

Surrogate-based sequential Bayesian experimental design using non-stationary Gaussian Processes

Inferring arbitrary quantities of interest (QoI) using limited computational or, in realistic scenarios, financial budgets, is a challenging problem that requires sophisticated strategies for the optimal allocation of the available resources. Bayesian optimal experimental design identifies the optimal set of design locations for the purpose of solving a parameter inference problem and the optimality criterion is typically associated with maximizing the worth of information in the experimental measurements. Sequential design strategies further identify the optimal design in a sequential manner, starting from a initial budget and iteratively selecting new optimal points until either an accuracy threshold is reached, or a cost limit is exceeded. In this paper, we present a generic sequential Bayesian experimental design framework that relies on maximizing an information theoretic design criterion, namely the Expected Information Gain, in order to infer QoIs formed as nonlinear operators acting on black-box functions. Our framework relies on modeling the underlying response function using non-stationary Gaussian Processes, thus enabling efficient sampling from the QoI in order to provide Monte Carlo estimators for the design criterion. We demonstrate the performance of our method on an engineering problem of steel wire manufacturing and compare it with two classic approaches: uncertainty sampling and expected improvement.

On nonlinear Feynman–Kac formulas for viscosity solutions of semilinear parabolic partial differential equations

The classical Feynman–Kac identity builds a bridge between stochastic analysis and partial differential equations (PDEs) by providing stochastic representations for classical solutions of linear Kolmogorov PDEs. This opens the door for the derivation of sampling based Monte Carlo approximation methods, which can be meshfree and thereby stand a chance to approximate solutions of PDEs without suffering from the curse of dimensionality. In this paper, we extend the classical Feynman–Kac formula to certain semilinear Kolmogorov PDEs. More specifically, we identify suitable solutions of stochastic fixed point equations (SFPEs), which arise when the classical Feynman–Kac identity is formally applied to semilinear Kolmorogov PDEs, as viscosity solutions of the corresponding PDEs. This justifies, in particular, employing full-history recursive multilevel Picard (MLP) approximation algorithms, which have recently been shown to overcome the curse of dimensionality in the numerical approximation of solutions of SFPEs, in the numerical approximation of semilinear Kolmogorov PDEs.

Complexity Analysis of stochastic gradient methods for PDE-constrained optimal Control Problems with uncertain parameters

We consider the numerical approximation of an optimal control problem for an elliptic Partial Differential Equation (PDE) with random coefficients. Specifically, the control function is a deterministic, distributed forcing term that minimizes the expected squaredL2misfit between the state (i.e.solution to the PDE) and a target function, subject to a regularization for well posedness. For the numerical treatment of this risk-averse Optimal Control Problem (OCP) we consider a Finite Element discretization of the underlying PDE, a Monte Carlo sampling method, and gradient-type iterations to obtain the approximate optimal control. We provide full error and complexity analyses of the proposed numerical schemes. In particular we investigate the complexity of a conjugate gradient method applied to the fully discretized OCP (so called Sample Average Approximation), in which the Finite Element discretization and Monte Carlo sample are chosen in advance and kept fixed over the iterations. This is compared with aStochastic Gradientmethod on a fixed or varying Finite Element discretization, in which the expectation in the computation of the steepest descent direction is approximated by Monte Carlo estimators, independent across iterations, with small sample sizes. We show in particular that the second strategy results in an improved computational complexity. The theoretical error estimates and complexity results are confirmed by numerical experiments.

Primary-space Adaptive Control Variates Using Piecewise-polynomial Approximations

We present an unbiased numerical integration algorithm that handles both low-frequency regions and high-frequency details of multidimensional integrals. It combines quadrature and Monte Carlo integration by using a quadrature-based approximation as a control variate of the signal. We adaptively build the control variate constructed as a piecewise polynomial, which can be analytically integrated, and accurately reconstructs the low-frequency regions of the integrand. We then recover the high-frequency details missed by the control variate by using Monte Carlo integration of the residual. Our work leverages importance sampling techniques by working in primary space, allowing the combination of multiple mappings; this enables multiple importance sampling in quadrature-based integration. Our algorithm is generic and can be applied to any complex multidimensional integral. We demonstrate its effectiveness with four applications with low dimensionality: transmittance estimation in heterogeneous participating media, low-order scattering in homogeneous media, direct illumination computation, and rendering of distribution effects. Finally, we show how our technique is extensible to integrands of higher dimensionality by computing the control variate on Monte Carlo estimates of the high-dimensional signal, and accounting for such additional dimensionality on the residual as well. In all cases, we show accurate results and faster convergence compared to previous approaches.

Sequential Monte Carlo Filters with Parameters Learning for Commodity Pricing Models

In this article, an estimation methodology based on the sequential Monte Carlo algorithm is proposed, thatjointly estimate the states and parameters, the relationship between the prices of futures contracts and the spot prices of primary products is determined, the evolution of prices and the volatility of the historical data of the primary market (Gold and Soybean) are analyzed. Two stochastic models for an estimate the states and parameters are considered, the parameters and states describe physical measure (associated with the price) and risk-neutral measure (associated with the markets to futures), the price dynamics in the short-term through the reversion to the mean and volatility are determined, while that in the long term through markets to futures. Other characteristics such as seasonal patterns, price spikes, market dependent volatilities, and non-seasonality can also be observed. In the methodology, a parameter learning algorithm is used, specifically, three algorithms are proposed, that is the sequential Monte Carlo estimation (SMC) for state space modelswith unknown parameters: the first method is considered a particle filter that is based on the sampling algorithm of sequential importance with resampling (SISR). The second implemented method is the Storvik algorithm [19], the states and parameters of the posterior distribution are estimated that have supported in low-dimensional spaces, a sufficient statistics from the sample of the filtered distribution is considered. The third method is (PLS) Carvalho’s Particle Learning and Smoothing algorithm [31]. The cash prices of the contracts with future delivery dates are analyzed. The results indicate postponement of payment, the future prices on different maturity dates with the spot price are highly correlated. Likewise, the contracts with a delivery date for the last periods of the year 2017, the spot price lower than the prices of the contracts with expiration date for 12 and 24 months is found, opposite occurs in the contracts with expiration date for 1 and 6 months.

Solving high-dimensional Hamilton\u2013Jacobi\u2013Bellman PDEs using neural networks: perspectives from the theory of controlled diffusions and measures on path space

Optimal control of diffusion processes is intimately connected to the problem of solving certain Hamilton–Jacobi–Bellman equations. Building on recent machine learning inspired approaches towards high-dimensional PDEs, we investigate the potential of iterative diffusion optimisation techniques, in particular considering applications in importance sampling and rare event simulation, and focusing on problems without diffusion control, with linearly controlled drift and running costs that depend quadratically on the control. More generally, our methods apply to nonlinear parabolic PDEs with a certain shift invariance. The choice of an appropriate loss function being a central element in the algorithmic design, we develop a principled framework based on divergences between path measures, encompassing various existing methods. Motivated by connections to forward-backward SDEs, we propose and study the novel log-variance divergence, showing favourable properties of corresponding Monte Carlo estimators. The promise of the developed approach is exemplified by a range of high-dimensional and metastable numerical examples.

Bayesian inference of Unit Gompertz distribution based on dual generalized order statistics

In this article, we consider the estimation problem of Unit Gompertz distribution with parameters α and β under the framework of dual generalized order statistics. This article is purely devoted to present the Bayesian view of estimation of Unit Gompertz distribution. For this purpose, we consider two widely popular approximation methods called Markov chain Monte Carlo and Lindley approximation methods. The results are derived under the symmetric (squared error) and asymmetric (Linear exponential and General entropy) loss functions. Since the order statistics and lower record values are the particular cases of the dual generalized order statistics, a simulation study is provided for order statistics and lower record values to observe the behavior of estimators. The average lengths of highest posterior density intervals of α, β, and R(t) are calculated for 95% confidence coefficient. Finally, real data applications are reported for lower record values and order statistics, separately, to show the practical aspects of the derived results.

Multilevel Monte-Carlo for computing the SCR with the standard formula and other stress tests

This paper studies the multilevel Monte-Carlo estimator for the expectation of a maximum of conditional expectations. This problem arises naturally when considering many stress tests and appears in the calculation of the interest rate module of the standard formula for the SCR. We obtain theoretical convergence results that complement the recent work of Giles and Goda (2019) and give some additional tractability through a parameter that somehow describes regularity properties around the maximum. We then apply the MLMC estimator to the calculation of the SCR at future dates with the standard formula for an ALM savings business on life insurance. We compare it with estimators obtained with Least Squares Monte-Carlo or Neural Networks. We find that the MLMC estimator is computationally more efficient and has the main advantage to avoid regression issues, which is particularly significant in the context of projection of a balance sheet by an insurer due to the path dependency. Last, we discuss the potential of this numerical method and analyse in particular the effect of the portfolio allocation on the SCR at future dates.

Measuring Response Style Stability Across Constructs With Item Response Trees.

Individual response style behaviors, unrelated to the latent trait of interest, may influence responses to ordinal survey items. Response style can introduce bias in the total score with respect to the trait of interest, threatening valid interpretation of scores. Despite claims of response style stability across scales, there has been little research into stability across multiple scales from the beneficial perspective of item response trees. This study examines an extension of the IRTree methodology to include mixed item formats, providing an empirical example of responses to three scales measuring perceptions of social media, climate change, and medical marijuana use. Results show extreme and midpoint response styles were not stable across scales within a single administration and 5-point Likert-type items elicited higher levels of extreme response style than the 4-point items. Latent trait of interest estimation varied, particularly at the lower end of the score distribution, across response style models, demonstrating as appropriate response style model is important for adequate trait estimation using Bayesian Markov chain Monte Carlo estimation.

Invariant Procedure for Error Sensitivity Analysis Applied to Cyclic Creep Modelling

Cyclic creep, also known as ratcheting, is a progressive accumulation of inelastic strain under cyclic stress-controlled loading. This mechanical effect is of great importance for numerous applications. In engineering practice, phenomenological models of cyclic creep are calibrated against a limited set of macroscopic test data. Since the testing results are prone to systematic and non-systematic experimental errors, the impact of experimental errors on the quality of simulation has to be analysed. A simple inspection procedure is demonstrated and tested. Based on the Monte Carlo computations, it allows for analysis of error propagation through the simulation cycle. The focus of the paper is on the independence of the procedure from the chosen model parametrization. For demonstration purposes, cyclic creep of VT6 alloy is simulated. The corresponding macroscopic constitutive equations are based on the second Ohno-Wang model, combined with refined rule of isotropic hardening. Two different parametrizations are introduced to show that the procedure predicts the same results for both of them.

Antithetic multilevel sampling method for nonlinear functionals of measure

Let μ∈ P2(Rd), where P2(Rd) denotes the space of square integrable probability measures, and consider a Borel-measurable function Φ:P2(Rd)→R. In this paper we develop an antithetic Monte Carlo estimator (A-MLMC) for Φ(μ), which achieves sharp error bound under mild regularity assumptions. The estimator takes as input the empirical laws μN=1N∑ i=1NδX i, where (a) (Xi)i=1N is a sequence of i.i.d. samples from μ or (b) (Xi)i=1N is a system of interacting particles (diffusions) corresponding to a McKean–Vlasov stochastic differential equation (McKV-SDE). Each case requires a separate analysis. For a mean-field particle system, we also consider the empirical law induced by its Euler discretisation which gives a fully implementable algorithm. As by-products of our analysis, we establish a dimension-independent rate of uniform strong propagation of chaos, as well as an L2 estimate of the antithetic difference for i.i.d. random variables corresponding to general functionals defined on the space of probability measures.

Efficient Steady-State Simulation of High-Dimensional Stochastic Networks

We propose and study an asymptotically optimal Monte Carlo estimator for steady-state expectations of a d-dimensional reflected Brownian motion (RBM). Our estimator is asymptotically optimal in the sense that it requires [Formula: see text] (up to logarithmic factors in d) independent and identically distributed scalar Gaussian random variables in order to output an estimate with a controlled error. Our construction is based on the analysis of a suitable multilevel Monte Carlo strategy which, we believe, can be applied widely. This is the first algorithm with linear complexity (under suitable regularity conditions) for a steady-state estimation of RBM as the dimension increases.

Using global variance-based sensitivity analysis to prioritise bridge retrofits in a regional road network subject to seismic hazard

This paper presents a novel method for prioritising bridge retrofits within a regional road network subject to uncertain seismic hazard, using a technique that accounts for network performance while avoiding the combinatoric computational costs of exhaustive searches. Using global variance-based sensitivity analysis, a probabilistic ranking of bridges is determined according to how much their retrofit statuses influence the expected cost of the road network disruption. Bridges’ total-order sensitivity (Sobol’) indices are estimated with respect to the expected cost using the hybrid-point Monte Carlo approximation method. A bridge’s total-order Sobol’ index measures how much its retrofit status influences the variance of the expected cost of the road network performance and accounts for the effect of its interactions with other bridges’ retrofit states. For 71 highway bridges in San Francisco, a retrofit strategy based on bridges’ total-order Sobol’ indices outperforms other heuristic strategies. The proposed method remains computationally tractable while accounting for the probabilistic nature of the seismic hazard, the uniqueness of individual bridges, network effects, and decision-makers’ priorities. Because this method leverages existing risk assessment tools and models without imposing further assumptions, it should be extensible to other types of networks under different types of hazards and to other decision variables.