Monte Carlo Errors Research Articles

Sequential multilevel assimilation of inverted seismic data

We consider estimation of absolute permeability from inverted seismic data. Large amounts of simultaneous data, such as inverted seismic data, enhance the negative effects of Monte Carlo errors in ensemble-based Data Assimilation (DA). Multilevel (ML) models consist of a selection of models with different fidelities. Multilevel Data Assimilation (MLDA) attempts to obtain a better statistical accuracy with a small sacrifice of the numerical accuracy. Spatial grid coarsening is one way of generating an ML model. It has been shown that coarsening the spatial grid results in a problem with weaker nonlinearity, and hence, in a less challenging problem than the problem on the original fine grid. Accordingly, formulating a sequential MLDA algorithm which uses the coarser models in the first steps of the DA, followed by the finer models, helps to find an approximation to the solution of the inverse problem at the first steps and gradually converge to the solution. We present two variants of a sequential MLDA algorithm and compare their performance with both conventional DA algorithms and a simultaneous (i.e., using all the models on the different grids simultaneously) MLDA algorithm using numerical experiments. Both posterior parameters and posterior model forecasts are compared qualitatively and quantitatively. The results from numerical experiments suggest that all MLDA algorithms generally perform better than the conventional DA algorithms. In estimation of the posterior parameter fields, the simultaneous MLDA algorithm and one of the variants of sequential MLDA (SMLES-H) perform similarly and slightly better than the other variant (SMLES-S). While in estimation of the posterior model forecasts, SMLES-S clearly performs better than both the simultaneous MLDA algorithm and SMLES-H.

A stochastic gradient algorithm with momentum terms for optimal control problems governed by a convection–diffusion equation with random diffusivity

In this paper, we focus on a numerical investigation of a strongly convex and smooth optimization problem subject to a convection–diffusion equation with uncertain terms. Our approach is based on stochastic approximation where true gradient is replaced by a stochastic ones with suitable momentum term to minimize the objective functional containing random terms. A full error analysis including Monte Carlo, finite element, and stochastic momentum gradient iteration errors is done. Numerical examples are presented to illustrate the performance of the proposed stochastic approximations in the PDE-constrained optimization setting.

Bridging the gap between Monte Carlo simulations and measurements of the LISA Pathfinder test-mass charging for LISA

Context. Cubic gold-platinum free-falling test masses (TMs) constitute the mirrors of future LISA and LISA-like interferometers for low-frequency gravitational wave detection in space. High-energy particles of Galactic and solar origin charge the TMs and thus induce spurious electrostatic and magnetic forces that limit the sensitivity of these interferometers. Prelaunch Monte Carlo simulations of the TM charging were carried out for the LISA Pathfinder (LPF) mission, that was planned to test the LISA instrumentation. Measurements and simulations were compared during the mission operations. The measured net TM charging agreed with simulation estimates, while the charging noise was three to four times higher. Aims. We aim to bridge the gap between LPF TM charging noise simulations and observations. Methods. New Monte Carlo simulations of the LPF TM charging due to both Galactic and solar particles were carried out with the FLUKA/LEI toolkit. This allowed propagating low-energy electrons down to a few electronvolt. Results. These improved FLUKA/LEI simulations agree with observations gathered during the mission operations within statistical and Monte Carlo errors. The charging noise induced by Galactic cosmic rays is about one thousand charges per second. This value increases to tens of thousands charges per second during solar energetic particle events. Similar results are expected for the LISA TM charging.

Multivariate Bayesian Arm-Based Network Meta-Analysis of Pharmacological Interventions for the Treatment of Acute Bipolar Mania in Adults.

In a network meta-analysis (NMA), multiple treatments can be compared simultaneously by aggregating pieces of evidence from direct as well as indirect treatment comparisons in different randomized controlled trials (RCTs). Conventional NMA are performed using a normal approximation approach and can be applied for arm-level binary outcome data as well. This study aimed to estimate the treatment effects within a Bayesian framework using a binomial likelihood for a multivariate NMA model. The dataset consists of 57 RCTs comparing the effect of ten pharmacological drugs and a placebo for acute bipolar mania in adults. The binary outcomes of interest were treatment response and all-cause dropouts measured three weeks from the baseline. Binomial distribution was adopted for the number of events and the probability of event occurrence modeled on the logit scale. Jeffrey's Beta prior was considered for the heterogeneity and inconsistency of standard deviation (SD) parameters. Cholesky and spherical decomposition strategies were adopted for the between-study variance-covariance matrix. Deviance information criterion (DIC) indices were computed to determine the model fit. All results pertaining to Markov chain Monte Carlo simulations and all analyses were carried out in WinBUGS software. The estimated common heterogeneity SDs were similar, and the DIC values did not provide any evidence for superiority between the two decomposition strategies. The correlation (95% credible interval) between the outcomes was estimated as -0.31 (-0.71, -0.02) and -0.37 (-0.73, -0.03) for the Cholesky and spherical decompositions, respectively. Gelman-Rubin convergence statistics were stable, and Monte Carlo errors for all the parameters were around 0.005. Overall, olanzapine, paliperidone, and quetiapine were both significantly more effective and acceptable than a placebo when both the study outcomes were considered simultaneously. The findings favoring olanzapine, paliperidone, and quetiapine possess an excellent concordance with the one adopted in clinical practice, and the Canadian Network for Mood and Anxiety Treatments and Royal Australian and New Zealand College of Psychiatrists guidelines recommend these as first-line drugs for treating bipolar disorder.

Connecting the Dots: Numerical Randomized Hamiltonian Monte Carlo with State-Dependent Event Rates

Numerical generalized randomized Hamiltonian Monte Carlo is introduced, as a robust, easy to use and computationally fast alternative to conventional Markov chain Monte Carlo methods for continuous target distributions. A wide class of piecewise deterministic Markov processes generalizing Randomized HMC (Bou-Rabee and Sanz-Serna) by allowing for state-dependent event rates is defined. Under very mild restrictions, such processes will have the desired target distribution as an invariant distribution. Second, the numerical implementation of such processes, based on adaptive numerical integration of second order ordinary differential equations (ODEs) is considered. The numerical implementation yields an approximate, yet highly robust algorithm that, unlike conventional Hamiltonian Monte Carlo, enables the exploitation of the complete Hamiltonian trajectories (hence, the title). The proposed algorithm may yield large speedups and improvements in stability relative to relevant benchmarks, while incurring numerical biases that are negligible relative to the overall Monte Carlo errors. Granted access to a high-quality ODE code, the proposed methodology is both easy to implement and use, even for highly challenging and high-dimensional target distributions. Supplementary materials for this article are available online.

From offshore to onshore probabilistic tsunami hazard assessment via efficient Monte Carlo sampling.

SUMMARYOffshore Probabilistic Tsunami Hazard Assessments (offshore PTHAs) provide large-scale analyses of earthquake-tsunami frequencies and uncertainties in the deep ocean, but do not provide high-resolution onshore tsunami hazard information as required for many risk-management applications. To understand the implications of an offshore PTHA for the onshore hazard at any site, in principle the tsunami inundation should be simulated locally for every earthquake scenario in the offshore PTHA. In practice this is rarely feasible due to the computational expense of inundation models, and the large number of scenarios in offshore PTHAs. Monte Carlo methods offer a practical and rigorous alternative for approximating the onshore hazard, using a random subset of scenarios. The resulting Monte Carlo errors can be quantified and controlled, enabling high-resolution onshore PTHAs to be implemented at a fraction of the computational cost. This study develops efficient Monte Carlo approaches for offshore-to-onshore PTHA. Modelled offshore PTHA wave heights are used to preferentially sample scenarios that have large offshore waves near an onshore site of interest. By appropriately weighting the scenarios, the Monte Carlo errors are reduced without introducing bias. The techniques are demonstrated in a high-resolution onshore PTHA for the island of Tongatapu in Tonga, using the 2018 Australian PTHA as the offshore PTHA, while considering only thrust earthquake sources on the Kermadec-Tonga trench. The efficiency improvements are equivalent to using 4–18 times more random scenarios, as compared with stratified-sampling by magnitude, which is commonly used for onshore PTHA. The greatest efficiency improvements are for rare, large tsunamis, and for calculations that represent epistemic uncertainties in the tsunami hazard. To facilitate the control of Monte Carlo errors in practical applications, this study also provides analytical techniques for estimating the errors both before and after inundation simulations are conducted. Before inundation simulation, this enables a proposed Monte Carlo sampling scheme to be checked, and potentially improved, at minimal computational cost. After inundation simulation, it enables the remaining Monte Carlo errors to be quantified at onshore sites, without additional inundation simulations. In combination these techniques enable offshore PTHAs to be rigorously transformed into onshore PTHAs, with quantification of epistemic uncertainties, while controlling Monte Carlo errors.

Iterative multilevel assimilation of inverted seismic data

In ensemble-based data assimilation (DA), the ensemble size is usually limited to around one hundred. Straightforward application of ensemble-based DA can therefore result in significant Monte Carlo errors, often manifesting themselves as severe underestimation of parameter uncertainties. Localization is the conventional remedy for this problem. Assimilation of large amounts of simultaneous data enhances the negative effects of Monte Carlo errors. Use of lower-fidelity models reduces the computational cost per ensemble member and therefore renders the possibility to reduce Monte Carlo errors by increasing the ensemble size, but it also adds to the modeling error. Multilevel data assimilation (MLDA) uses a selection of models forming hierarchies of both computational cost and computational accuracy, and tries to balance between Monte Carlo errors and modeling errors. In this work, we assess a recently developed MLDA algorithm, the Multilevel Hybrid Ensemble Smoother (MLHES), and introduce and assess an iterative version of this algorithm, the Iterative Multilevel Hybrid Ensemble Smoother (IMLHES). In our assessments, we compare these algorithms with conventional single-level DA algorithms with localization. To this end, a typical example of large amount of spatially distributed data, i.e. inverted seismic data, is considered and three data sets of this kind are assimilated in three different petroleum reservoir models. Qualitatively evaluating the DA outcomes, it is found that multilevel algorithms outperform their conventional single-level counterparts in obtaining the posterior statistics of both uncertain parameters and model forecasts. Additionally, it is observed that IMLHES performs better than MLHES in the same regard, and also successfully converges to the proximity of solution in a case where the considered iterative single-level algorithm did not converge to the global optimum.

Multilevel Assimilation of Inverted Seismic Data With Correction for Multilevel Modeling Error

With large amounts of simultaneous data, like inverted seismic data in reservoir modeling, negative effects of Monte Carlo errors in straightforward ensemble-based data assimilation (DA) are enhanced, typically resulting in underestimation of parameter uncertainties. Utilization of lower fidelity reservoir simulations reduces the computational cost per ensemble member, thereby rendering the possibility of increasing the ensemble size without increasing the total computational cost. Increasing the ensemble size will reduce Monte Carlo errors and therefore benefit DA results. The use of lower fidelity reservoir models will however introduce modeling errors in addition to those already present in conventional fidelity simulation results. Multilevel simulations utilize a selection of models for the same entity that constitute hierarchies both in fidelities and computational costs. In this work, we estimate and approximately account for the multilevel modeling error (MLME), that is, the part of the total modeling error that is caused by using a multilevel model hierarchy, instead of a single conventional model to calculate model forecasts. To this end, four computationally inexpensive approximate MLME correction schemes are considered, and their abilities to correct the multilevel model forecasts for reservoir models with different types of MLME are assessed. The numerical results show a consistent ranking of the MLME correction schemes. Additionally, we assess the performances of the different MLME-corrected model forecasts in assimilation of inverted seismic data. The posterior parameter estimates from multilevel DA with and without MLME correction are compared to results obtained from conventional single-level DA with localization. It is found that multilevel DA (MLDA) with and without MLME correction outperforms conventional DA with localization. The use of all four MLME correction schemes results in posterior parameter estimates with similar quality. Results obtained with MLDA without any MLME correction were also of similar quality, indicating some robustness of MLDA toward MLME.

Distributing Monte Carlo Errors as a Blue Noise in Screen Space by Permuting Pixel Seeds Between Frames

AbstractRecent work has shown that distributing Monte Carlo errors as a blue noise in screen space improves the perceptual quality of rendered images. However, obtaining such distributions remains an open problem with high sample counts and high‐dimensional rendering integrals. In this paper, we introduce a temporal algorithm that aims at overcoming these limitations. Our algorithm is applicable whenever multiple frames are rendered, typically for animated sequences or interactive applications. Our algorithm locally permutes the pixel sequences (represented by their seeds) to improve the error distribution across frames. Our approach works regardless of the sample count or the dimensionality and significantly improves the images in low‐varying screen‐space regions under coherent motion. Furthermore, it adds negligible overhead compared to the rendering times. Note: our supplemental material provides more results with interactive comparisons against previous work.

Reduced basis method for the adapted mesh and Monte Carlo methods applied to an elliptic stochastic problem

In this paper, we consider a stochastic elliptic partial differential system and we aim to approximate the solution using the Monte Carlo method based on the finite elements method. To speed up the resolution and reduce the CPU time of computation, we propose to couple the reduced basis method with the adapted mesh method based on an a posteriori error estimate. Balancing the discretization and the Monte Carlo errors is very important to avoid performing an excessive number of iterations. Numerical experiments show and confirm the efficiency of our proposed algorithm.

Evaluating Factors Influencing the Severity of Three-Plus Multiple-Vehicle Crashes using Real-Time Traffic Data

Multiple-vehicle crashes involving at least two vehicles constitute over 70% of fatal and injury crashes in the U.S. Moreover, multiple-vehicle crashes involving three or more vehicles (3+) are usually more severe compared with the crashes involving only two vehicles. This study focuses on developing 3+ multiple-vehicle crash severity models for a freeway section using real-time traffic data and crash data for the years 2014–2016. The study corridor is a 111-mile section on I-4 in Orlando, Florida. Crash injury severity was classified as a binary outcome (fatal/severe injury and minor/no injury crashes). For the purpose of identifying the reliable relationship between the 3+ severe multiple-vehicle crashes and the identified explanatory variables, a binary probit model with Dirichlet random effect parameter was used. More specifically, Dirichlet random effect model was introduced to account for unobserved heterogeneity in the crash data. The probit model was implemented using a Bayesian framework and the ratios of the Monte Carlo errors were monitored to achieve parameter estimation convergence. The following variables were found significant at the 95% Bayesian credible interval: logarithm of average vehicle speed, logarithm of average equivalent 10-minute hourly volume, alcohol involvement, lighting condition, and number of vehicles involved (3, or >3) in multiple-vehicle crashes. Further analysis involved analyzing the posterior probability distributions of these significant variables. The study findings can be used to associate certain traffic conditions with severe injury crashes involving 3+ multiple vehicles, and can help develop effective crash injury reduction strategies based on real-time traffic data.

Melting Sensitivities - Exact and Approximate Margin Valuation Adjustments

In this paper we present various methods for fast and accurate calculation of margin valuation adjustments (MVA). We consider the calculation of an MVA using sensitivities for the determination of the initial margin and consider the IDSA SIMM (ISDA Standard Initial Margin Model) as a benchmark. To assess the quality of the various approximations we calculate an - up to Monte-Carlo and discretization errors - exact MVA applying adjoint algorithmic differentiation (AAD). This exact MVA and exact forward initial margin serves as a benchmark for the various approximation methods. It turns out that the conversion from model sensitivities to market rate sensitivities - as required by SIMM - constitutes a computationally expensive calculation step. The step requires the calculation of a (pseudo-) inverse of stochastic matrices. We therefore introduce approximation methods allowing for more efficient MVA calculation.

Local Kernel Dimension Reduction in Approximate Bayesian Computation

Approximate Bayesian Computation (ABC) is a popular sampling method in applications involving intractable likelihood functions. Without evaluating the likelihood function, ABC approximates the posterior distribution by the set of accepted samples which are simulated with parameters drawn from the prior distribution, where acceptance is determined by the distance between the summary statistics of the sample and the observation. The sufficiency and dimensionality of the summary statistics play a central role in the application of ABC. This paper proposes Local Gradient Kernel Dimension Reduction (LGKDR) to construct low dimensional summary statistics for ABC. The proposed method identifies a sufficient subspace of the original summary statistics by implicitly considers all nonlinear transforms therein, and a weighting kernel is used for the concentration of the projections. No strong assumptions are made on the marginal distributions nor the regression model, permitting usage in a wide range of applications. Experiments are done with both simple rejection ABC and sequential Monte Carlo ABC methods. Results are reported as competitive in the former and substantially better in the latter cases in which Monte Carlo errors are compressed as much as possible.

Spatio-Temporal Variation of HIV Infection in Kenya

Disease mapping is the study of the distribution of disease relative risks or rates in space and time, and normally uses generalized linear mixed models (GLMMs) which includes fixed effects and spatial, temporal, and spatio-temporal random effects. Model fitting and statistical inference are commonly accomplished through the empirical Bayes (EB) and fully Bayes (FB) approaches. The EB approach usually relies on the penalized quasi-likelihood (PQL), while the FB approach, which has increasingly become more popular in the recent past, usually uses Markov chain Monte Carlo (McMC) techniques. However, there are many challenges in conventional use of posterior sampling via McMC for inference. This includes the need to evaluate convergence of posterior samples, which often requires extensive simulation and can be very time consuming. Spatio-temporal models used in disease mapping are often very complex and McMC methods may lead to large Monte Carlo errors if the dimension of the data at hand is large. To address these challenges, a new strategy based on integrated nested Laplace approximations (INLA) has recently been recently developed as a promising alternative to the McMC. This technique is now becoming more popular in disease mapping because of its ability to fit fairly complex space-time models much more quickly than the McMC. In this paper, we show how to fit different spatio-temporal models for disease mapping with INLA using the Leroux CAR prior for the spatial component, and we compare it with McMC using Kenya HIV incidence data during the period 2013-2016.

Maximum Likelihood Estimation of Latent Variable Models by SMC with Marginalization and Data Cloning

A data-cloning SMC² method is proposed as a general purpose optimization routine for estimating latent variable models by maximum likelihood. The latent variables are first marginalized out by SMC at any fixed parameter value, and the model parameters are then estimated by density tempered SMC. The data-cloning step is employed to efficiently reduce Monte Carlo errors inherent in the SMC² algorithm and also to effectively address multi-modality present in typical objective functions. This new method has wide applicability and can be massively parallelized to take advantage of typical computers today.

Maximum Likelihood Estimation of Latent Variable Models by SMC with Marginalization and Data Cloning

A data-cloning SMC2 maximum likelihood estimation algorithm is proposed as a general-purpose optimization routine for models with latent variables. Our algorithm first marginalizes out latent variables by applying one layer of SMC at a fixed parameter value, and then estimates the model parameters by another layer of SMC utilizing density-tempering. Data-cloning is employed to effectively reduce Monte Carlo errors inherent in the SMC2 algorithm, and also to address multi-modality present in typical latent variable models. This new method has wide applicability and can be massively parallelized to take advantage of typical computers today.

The characterization of Monte Carlo errors for the quantification of the value of forensic evidence

ABSTRACTRecent developments in forensic science have lead to a proliferation of methods for quantifying the probative value of evidence by constructing a Bayes Factor that allows a decision-maker to select between the prosecution and defense models. Unfortunately, the analytical form of a Bayes Factor is often computationally intractable. A typical approach in statistics uses Monte Carlo integration to numerically approximate the marginal likelihoods composing the Bayes Factor. This article focuses on developing a generally applicable method for characterizing the numerical error associated with Monte Carlo integration techniques used in constructing the Bayes Factor. The derivation of an asymptotic Monte Carlo standard error (MCSE) for the Bayes Factor will be presented and its applicability to quantifying the value of evidence will be explored using a simulation-based example involving a benchmark data set. The simulation will also explore the effect of prior choice on the Bayes Factor approximations and corresponding MCSEs.

Circulant Embedding of Approximate Covariances for Inference From Gaussian Data on Large Lattices

ABSTRACTRecently proposed computationally efficient Markov chain Monte Carlo (MCMC) and Monte Carlo expectation–maximization (EM) methods for estimating covariance parameters from lattice data rely on successive imputations of values on an embedding lattice that is at least two times larger in each dimension. These methods can be considered exact in some sense, but we demonstrate that using such a large number of imputed values leads to slowly converging Markov chains and EM algorithms. We propose instead the use of a discrete spectral approximation to allow for the implementation of these methods on smaller embedding lattices. While our methods are approximate, our examples indicate that the error introduced by this approximation is small compared to the Monte Carlo errors present in long Markov chains or many iterations of Monte Carlo EM algorithms. Our results are demonstrated in simulation studies, as well as in numerical studies that explore both increasing domain and fixed domain asymptotics. We compare the exact methods to our approximate methods on a large satellite dataset, and show that the approximate methods are also faster to compute, especially when the aliased spectral density is modeled directly. Supplementary materials for this article are available online.

Reducing Monte Carlo error in the Bayesian estimation of risk ratios using log-binomial regression models.

In cohort studies, binary outcomes are very often analyzed by logistic regression. However, it is well known that when the goal is to estimate a risk ratio, the logistic regression is inappropriate if the outcome is common. In these cases, a log-binomial regression model is preferable. On the other hand, the estimation of the regression coefficients of the log-binomial model is difficult owing to the constraints that must be imposed on these coefficients. Bayesian methods allow a straightforward approach for log-binomial regression models and produce smaller mean squared errors in the estimation of risk ratios than the frequentist methods, and the posterior inferences can be obtained using the software WinBUGS. However, Markov chain Monte Carlo methods implemented in WinBUGS can lead to large Monte Carlo errors in the approximations to the posterior inferences because they produce correlated simulations, and the accuracy of the approximations are inversely related to this correlation. To reduce correlation and to improve accuracy, we propose a reparameterization based on a Poisson model and a sampling algorithm coded in R.

The τ leptons theory and experimental data: Monte Carlo, fits, software and systematic errors

Status of τ lepton decay Monte Carlo generator TAUOLA is reviewed. Recent efforts on development of new hadronic currents are presented. Multitude new channels for anomalous τ decay modes and parameterization based on defaults used by BaBar collaboration are introduced. Also parameterization based on theoretical considerations are presented as an alternative. Lesson from comparison and fits to the BaBar and Belle data is recalled. It was found that as in the past, in particular at a time of comparisons with CLEO and ALEPH data, proper fitting, to as detailed as possible representation of the experimental data, is essential for appropriate developments of models of τ decays.In the later part of the presentation, use of the TAUOLA program for phenomenology of W, Z, H decays at LHC is adressed. Some new results, relevant for QED bremsstrahlung in such decays are presented as well.