Monte Carlo Estimation Research Articles

Meta variance reduction for Monte Carlo estimation of energetic particle confinement during stellarator optimization

This work introduces meta estimators that combine multiple multifidelity techniques based on control variates, importance sampling, and information reuse to yield a quasi-multiplicative amount of variance reduction. The proposed meta estimators are particularly efficient within outer-loop applications when the input distribution of the uncertainties changes during the outer loop, which is often the case in reliability-based design and shape optimization. We derive asymptotic bounds of the variance reduction of the meta estimators in the limit of convergence of the outer-loop results. We demonstrate the meta estimators, using data-driven surrogate models and biasing densities, on a design problem under uncertainty motivated by magnetic confinement fusion, namely the optimization of stellarator coil designs to maximize the estimated confinement of energetic particles. The meta estimators outperform all of their constituent variance reduction techniques alone, ultimately yielding two orders of magnitude speedup compared to standard Monte Carlo estimation at the same computational budget.

Monte Carlo Estimation of Mass Attenuation Coefficients of Some ZnO-doped Glass Samples

Bu çalışma, Monte Carlo yöntemi kullanılarak ZnO katkılı bazı cam örneklerinin toplam kütle zayıflama katsayılarını hesaplamayı amaçlamaktadır. Simülasyonlar basitçe, silindir şeklinde soğurucu bir malzeme, NaI(Tl) detektörü,kolimatör ve soğurucuya doğru paralel bir ışın olarak yönlendirilen mono-enerjik fotonlar yayan nokta kaynaktan oluşmaktadır. Problem geometrisindeki tüm bileşenler, numune dışındaki malzemelerde herhangi bir etkileşimi önlemek için vakum bir küre ile çevrelenmiştir. Simülasyon düzeneği bu şekliyle, saçılan fotonların dedektördeki toplam akıya katkı yapması engellemiştir. Simülasyonlar, 10 keV-20 MeV enerji aralığında 39 farklı foton enerjisinde GAMOS 6.2 paketi kullanılarak gerçekleştirilmiştir. Bu çalışma .çinko oksitin cam örneklere ilave edilmesinin cam örneklerinin radyasyon soğurma özelliklerini arttırdığını göstermiştir. Çalışmanın sonuçları XCOM değerleri ile karşılaştırıldığında iyi bir uyum görülmektedir. Monte Carlo tekniğinin, basit bir model geometrisi kullanarak, geniş bir enerji aralığında, kütlesel zayıflama katsayılarının hesaplanması için bir alternatif olarak kullanılabileceğini ortaya koymuştur.

Scalable Bayesian optimization with randomized prior networks

Several fundamental problems in science and engineering consist of global optimization tasks involving unknown high-dimensional (black-box) functions that map a set of controllable variables to the outcomes of an expensive experiment. Bayesian Optimization (BO) techniques are known to be effective in tackling global optimization problems using a relatively small number objective function evaluations, but their performance suffers when dealing with high-dimensional outputs. To overcome the major challenge of dimensionality, here we propose a deep learning framework for BO and sequential decision making based on bootstrapped ensembles of neural architectures with randomized priors. Using appropriate architecture choices, we show that the proposed framework can approximate functional relationships between design variables and quantities of interest, even in cases where the latter take values in high-dimensional vector spaces or even infinite-dimensional function spaces. In the context of BO, we augmented the proposed probabilistic surrogates with re-parameterized Monte Carlo approximations of multiple-point (parallel) acquisition functions, as well as methodological extensions for accommodating black-box constraints and multi-fidelity information sources. We test the proposed framework against state-of-the-art methods for BO and demonstrate superior performance across several challenging tasks with high-dimensional outputs, including a constrained multi-fidelity optimization task involving shape optimization of rotor blades in turbo-machinery.

Research of nonparametric density estimation algorithms by applying clustering methods

One of the ways to improve the accuracy of probability density estimation is multi-mode density treating as the mixture of single-mode one. In this paper we offer to use data clustering in the first place and to estimate density in every cluster separately. To objectively compare the performance, Monte Carlo approximation is used. While using various methods to evaluate the accuracy of probability density estimations we tried to use clustered and not clustered data. In this paper we also tried to reveal the usefulness of using clustering for data generated by single-mode and multi-mode distributions.

A numerical method for solving snapping out Brownian motion in 2D bounded domains

Diffusion in heterogeneous media partitioned by semi-permeable interfaces has a wide range of applications in the physical and life sciences, including gas permeation in soils, diffusion magnetic resonance imaging (dMRI), drug delivery, thermal conduction in composite media, synaptic receptor trafficking, and intercellular gap junctions. At the single particle level, diffusion across a semi-permeable interface can be formulated in terms of so-called snapping out Brownian motion (SNOBM). The latter sews together successive rounds of reflected BM, each of which is restricted to one side of the interface. Each round of reflected BM is killed when the local time at the interface exceeds an independent, exponentially distributed random variable. (The local time specifies the amount of time a reflected Brownian particle spends in a neighborhood of the interface.) The particle then immediately resumes reflected BM on the same side or the other side of the interface according to a stochastic switch, and the process is iterated. In this paper, we develop a Monte Carlo algorithm for simulating a two-dimensional version of SNOBM, which is used to solve a first passage time (FPT) problem for diffusion in a domain with semi-permeable partially absorbing traps. Our method combines a walk-on-spheres (WOS) method with an efficient algorithm for computing the boundary local time that uses a Skorokhod integral representation of the latter. We validate our algorithm by comparing the Monte Carlo estimates of the mean first passage time (MFPT) to the exact solution for a single circular trap, and show that our simulations are consistent with asymptotic results obtained for the 2D narrow capture problem involving multiple small circular targets. We also discuss extensions to higher dimensions.

Path integrals formulations leading to propagator evaluation for coupled linear physics in large geometric models

Reformulating linear physics using second kind Fredholm equations is very standard practice. One of the straightforward consequences is that the resulting integrals can be expanded (when the Neumann expansion converges) and probabilized, leading to path statistics and Monte Carlo estimations. An essential feature of these algorithms is that they also allow to estimate propagators for all types of sources, including initial conditions. The resulting practice is a single Monte Carlo run, for one given set of sources, producing propagators that can later be used with any other set of sources for fast simulations, typically as parts of optimization, inversion, sensitivity analysis and command control algorithms. The present paper illustrates how this practice can be extended to problems involving several interacting physics, provided that their coupling is only at the boundary of the system or at interfaces between sub-parts, and may itself be given the form of a second kind Fredholm equation. A full practical implementation is described as part of the Stardis code, with the example of transferring heat via the coupling of radiation, reaction-diffusion and convection as typically expected in the multidisciplinary context of urban climate modeling. Besides, we show how recent advances in computer graphics indicate that these algorithms can be made numerically extremely efficient when facing large CAD geometries: computing the propagator becomes strictly independent of the geometry refinement, i.e. is identical whatever the number of triangles and tetraedra used to numerize the surface and volume descriptions. To the best of our knowledge this is the first report of propagator computations that remains practical for coupled physics in large CAD geometries. Program summaryProgram Title: Stardis 0.7.2 (built on stardis-solver 0.12.3)CPC Library link to program files:https://doi.org/10.17632/k76zrx4n6b.1Developer's repository link:https://www.meso-star.com/projects/stardis/stardis.htmlLicensing provisions: GPLv3Programming language: ANSI C and PythonNature of problem: Estimating temperatures in coupled heat transfer systems involving large CAD and/or large numbers of spatially distributed sources.Solution method: Stardis uses the Monte Carlo Method. Each temperature estimate constructs a propagator of each of the energy sources within the system: initial temperature, temperature boundary conditions, volume powers and surface fluxes. The propagator can be stored for further use outside the Monte Carlo code.Additional comments including restrictions and unusual features: The Stardis estimates of propagators are only reliable when radiative transfer can be linearized around a reference temperature. Stardis can deal with the nonlinearity of radiation when computing temperatures, but then nothing can be interpreted as a meaningfull propagator for any external usage.Python scripts were used to execute the different configurations, post-process the results and generate the different graphs.

A SAS Macro for Automated Stopping of Markov Chain Monte Carlo Estimation in Bayesian Modeling with PROC MCMC

A crucial challenge in Bayesian modeling using Markov chain Monte Carlo (MCMC) estimation is to diagnose the convergence of the chains so that the draws can be expected to closely approximate the posterior distribution on which inference is based. A close approximation guarantees that the MCMC error exhibits only a negligible impact on model estimates and inferences. However, determining whether convergence has been achieved can often be challenging and cumbersome when relying solely on inspecting the trace plots of the chain(s) or manually checking the stopping criteria. In this article, we present a SAS macro called %automcmc that is based on PROC MCMC and that automatically continues to add draws until a user-specified stopping criterion (i.e., a certain potential scale reduction and/or a certain effective sample size) is reached for the chain(s).

Variability Estimation in a Nonlinear Crack Growth Simulation Model With Controlled Parameters Using Designed Experiments Testing

Abstract Variability in multiple independent input parameters makes it difficult to estimate the resultant variability in a system's overall response. The propagation of errors (PE) and Monte Carlo (MC) techniques are two major methods to predict the variability of a system. However, the formalism of PE can lead to an inaccurate estimate for systems that have parameters varying over a wide range. For the latter, the results give a direct estimate of the variance of the response, but for complex systems with many parameters, the number of trials necessary to yield an accurate estimate can be sizeable to the point the technique becomes impractical. The effectiveness of a designed experiment (orthogonal array) methodology, as employed in Taguchi tolerance design (TD) method to estimate variability in complex systems is studied. We use a linear elastic three-point bending beam model and a nonlinear extended finite elements crack growth model to test and compare the PE and MC methods with the TD method. Our focus is on assessing their effectiveness under those specific circumstances. Results from an MC estimate, using 10,000 trials, serve as a referent to evaluate the result in both cases. We find that the PE method works suboptimal for a coefficient of variation above 5% in the input variables. In addition, we find that the TD method works very well with moderately sized trials of designed experiments for both models. Our results demonstrate how the variability estimation methods perform in the deterministic domain of numerical simulations and can assist in designing physical tests by providing a guideline performance measure.

Regularized Zero-Variance Control Variates

Zero-variance control variates (ZV-CV) are a post-processing method to reduce the variance of Monte Carlo estimators of expectations using the derivatives of the log target. Once the derivatives are available, the only additional computational effort lies in solving a linear regression problem. Significant variance reductions have been achieved with this method in low dimensional examples, but the number of covariates in the regression rapidly increases with the dimension of the target. In this paper, we present compelling empirical evidence that the use of penalized regression techniques in the selection of high-dimensional control variates provides performance gains over the classical least squares method. Another type of regularization based on using subsets of derivatives, or a priori regularization as we refer to it in this paper, is also proposed to reduce computational and storage requirements. Several examples showing the utility and limitations of regularized ZV-CV for Bayesian inference are given. The methods proposed in this paper are accessible through the R package ZVCV.

Integrative Bayesian tensor regression for imaging genetics applications.

Identifying biomarkers for Alzheimer's disease with a goal of early detection is a fundamental problem in clinical research. Both medical imaging and genetics have contributed informative biomarkers in literature. To further improve the performance, recently, there is an increasing interest in developing analytic approaches that combine data across modalities such as imaging and genetics. However, there are limited methods in literature that are able to systematically combine high-dimensional voxel-level imaging and genetic data for accurate prediction of clinical outcomes of interest. Existing prediction models that integrate imaging and genetic features often use region level imaging summaries, and they typically do not consider the spatial configurations of the voxels in the image or incorporate the dependence between genes that may compromise prediction ability. We propose a novel integrative Bayesian scalar-on-image regression model for predicting cognitive outcomes based on high-dimensional spatially distributed voxel-level imaging data, along with correlated transcriptomic features. We account for the spatial dependencies in the imaging voxels via a tensor approach that also enables massive dimension reduction to address the curse of dimensionality, and models the dependencies between the transcriptomic features via a Graph-Laplacian prior. We implement this approach via an efficient Markov chain Monte Carlo (MCMC) computation strategy. We apply the proposed method to the analysis of longitudinal ADNI data for predicting cognitive scores at different visits by integrating voxel-level cortical thickness measurements derived from T1w-MRI scans and transcriptomics data. We illustrate that the proposed imaging transcriptomics approach has significant improvements in prediction compared to prediction using a subset of features from only one modality (imaging or genetics), as well as when using imaging and transcriptomics features but ignoring the inherent dependencies between the features. Our analysis is one of the first to conclusively demonstrate the advantages of prediction based on combining voxel-level cortical thickness measurements along with transcriptomics features, while accounting for inherent structural information.

Generalized Method of Moments Estimation of Realized Stochastic Volatility Model

The purpose of this paper is to study the generalized method of moments (GMM) estimation procedures of the realized stochastic volatility model; we give the moment conditions for this model and then obtain the estimation of parameters. Then, we apply these moment conditions to the realized stochastic volatility model to improve the volatility prediction effect. This paper selects the Shanghai Composite Index (SSE) as the original data of model research and completes the volatility prediction under a realized stochastic volatility model. Markov chain Monte Carlo (MCMC) estimation and quasi-maximum likelihood (QML) estimation are applied to the parameter estimation of the realized stochastic volatility model to compare with the GMM method. And the volatility prediction accuracy of these three different methods is compared. The results of empirical research show that the effect of model prediction using the parameters obtained by the GMM method is close to that of the MCMC method, and the effect is obviously better than that of the quasi-maximum likelihood estimation method.

Inflatable aerodynamic decelerator for CubeSat reentry and recovery: IAD geometrical effects on the flowfield structure

In this investigation, Direct Simulation Monte Carlo computations were carried out to assess the influence of Inflatable Aerodynamic Decelerator (IAD) geometrical configurations on the flowfield during the atmospheric reentry. The non-reacting hypersonic rarefied flow over three IAD configurations coupled to a CubeSat was simulated considering a 0∘ angle of attack at an altitude of 105 km. According to the results, it was observed the formation of a strong and diffuse shock wave for all geometries considered in this investigation. However, a lower inflatable aeroshell angle is associated with a thinner shock wave and a maximum shock wave temperature closer to the shield's surface. These differences decrease in the flow expansion over the IAD shoulder. In the rear of the inflatable shields, a low-temperature and low-velocity region is observed, indicating that the IADs geometries successfully mitigate the harsh reentry conditions experienced by the payload. Finally, it was noticed that the wake region is larger for aerodynamic elongated shapes, in contrast to blunt geometries with the expense of slightly higher gas temperature closer to the front surface of the shield. No recirculation zone was observed in any of the simulated IAD configurations considered in this investigation.

BEATING A CONSTANT WEIGHT BENCHMARK: EASIER DONE THAN SAID

We determine a simple dynamic benchmark for asset allocation by solving an optimal stochastic control problem for outperforming the traditional constant proportion benchmark. An objective function based on a time averaged quadratic deviation from an elevated benchmark is proposed. We argue that this objective function combines the best features of tracking error and tracking difference. Assuming parametric models of the stock and bond processes, a closed form solution for the optimal control is obtained. The closed form optimal control is then clipped to prevent use of excessive leverage, and to prevent trading if insolvent. Monte Carlo computations using this clipped control are presented which show that for modest levels of outperformance (i.e. 80–170[Formula: see text]bps per year), this easily implementable strategy outperforms the traditional constant proportion benchmark with high probability. We advocate this clipped optimal strategy as a suitable benchmark for active asset allocation.

A point mass proposal method for Bayesian state-space model fitting

State-space models (SSMs) are commonly used to model time series data where the observations depend on an unobserved latent process. However, inference on the model parameters of an SSM can be challenging, especially when the likelihood of the data given the parameters is not available in closed-form. One approach is to jointly sample the latent states and model parameters via Markov chain Monte Carlo (MCMC) and/or sequential Monte Carlo approximation. These methods can be inefficient, mixing poorly when there are many highly correlated latent states or parameters, or when there is a high rate of sample impoverishment in the sequential Monte Carlo approximations. We propose a novel block proposal distribution for Metropolis-within-Gibbs sampling on the joint latent state and parameter space. The proposal distribution is informed by a deterministic hidden Markov model (HMM), defined such that the usual theoretical guarantees of MCMC algorithms apply. We discuss how the HMMs are constructed, the generality of the approach arising from the tuning parameters, and how these tuning parameters can be chosen efficiently in practice. We demonstrate that the proposed algorithm using HMM approximations provides an efficient alternative method for fitting state-space models, even for those that exhibit near-chaotic behavior.

Boundary Value Caching for Walk on Spheres

Grid-free Monte Carlo methods such as walk on spheres can be used to solve elliptic partial differential equations without mesh generation or global solves. However, such methods independently estimate the solution at every point, and hence do not take advantage of the high spatial regularity of solutions to elliptic problems. We propose a fast caching strategy which first estimates solution values and derivatives at randomly sampled points along the boundary of the domain (or a local region of interest). These cached values then provide cheap, output-sensitive evaluation of the solution (or its gradient) at interior points, via a boundary integral formulation. Unlike classic boundary integral methods, our caching scheme introduces zero statistical bias and does not require a dense global solve. Moreover we can handle imperfect geometry (e.g., with self-intersections) and detailed boundary/source terms without repairing or resampling the boundary representation. Overall, our scheme is similar in spirit to virtual point light methods from photorealistic rendering: it suppresses the typical salt-and-pepper noise characteristic of independent Monte Carlo estimates, while still retaining the many advantages of Monte Carlo solvers: progressive evaluation, trivial parallelization, geometric robustness, etc. We validate our approach using test problems from visual and geometric computing.

Insights into the effect of growth on the Ziff–Gulari–Barshad model and the film properties

We perform kinetic Monte Carlo computations with a modified Ziff–Gulari–Barshad (ZGB) model which considers the growth of a film. We show that the growth of the film significantly affects the conclusions that can be drawn from the ZGB model, even if the main mechanism, the surface reaction, remains the same. We compare the results of the growth model to the original ZGB and the phase transitions disappears; they are replaced by a smooth transition from 0 to full coverage. The latter observations qualitatively agree with experimental measurements for the CO2 formation. However, in the growth model the surface is always poisoned to a particular coverage values due to the local height differences of the lattice sites. Finally, a potential mechanism based only on surface phenomena which can lead to the decrement of the growth rate even if the amount of the precursor increases is explored.

Near-continuum, hypersonic oxygen flow over a double cone simulated by direct simulation Monte Carlo informed from quantum chemistry

A large-scale, fully resolved direct simulation Monte Carlo (DSMC) computation of a non-equilibrium, reactive flow of pure oxygen over a double cone is presented. Under the simulated near-continuum conditions, the computational demands are shown to be significant because of the wide range of length scales that must be resolved. Therefore, robust grid adaption capabilities and efficient parallelization of the Stochastic PArallel Rarefied-gas Time-accurate Analyzer (SPARTA) code that is utilized in this work are essential. The thermochemical and transport collision models were selected for efficiency and simplicity. First-principles data, obtained from the highly accurate direct molecular simulation method, were used to inform the collision models’ parameters. Importantly, because SPARTA implements molecular collision models using collision-specific energies, the resulting macroscopic relaxation rates were evaluateda posteriorivia zero-dimensional heat bath simulations. The comparisons of surface properties, namely heat flux and pressure, show very close agreement with previous computational fluid dynamics (CFD) results. Differences with the measurements were found to be similar to the CFD simulations. The unresolved discrepancy with the measurements could be due to inconsistent free stream conditions with the actual experimental data or missing physical phenomena altogether, for example atomic and molecular oxygen electronically excited states, three-dimensional effects, or more complex gas–surface interactions. As shown in this work, the advantages of obtaining a DSMC particle solution for these flows reside in the method's ability to be directly informed from first principles and to seamlessly describe internal energy non-equilibrium for all modes. With the advent of exascale computing and beyond, particle methods will be an increasingly important tool to verify the validity of physical assumptions in reduced-order models via fully resolved, experimental-scale simulations, down to the level of molecular-level distributions.

A deep learning approach to estimate x-ray scatter in digital breast tomosynthesis: From phantom models to clinical applications.

Digital breast tomosynthesis (DBT) has gained popularity as breast imaging modality due to its pseudo-3D reconstruction and improved accuracy compared to digital mammography. However, DBT faces challenges in image quality and quantitative accuracy due to scatter radiation. Recent advancements in deep learning (DL) have shown promise in using fast convolutional neural networks for scatter correction, achieving comparable results to Monte Carlo (MC) simulations. To predict the scatter radiation signal in DBT projections within clinically-acceptable times and using only clinically-available data, such as compressed breast thickness and acquisition angle. MC simulations to obtain scatter estimates were generated from two types of digital breast phantoms. One set consisted of 600 realistically-shaped homogeneous breast phantoms for initial DL training. The other set was composed of 80 anthropomorphic phantoms, containing realistic internal tissue texture, aimed at fine tuning the DL model for clinical applications. The MC simulations generated scatter and primary maps per projection angle for a wide-angle DBT system. Both datasets were used to train (using 7680 projections from homogeneous phantoms), validate (using 960 and 192 projections from the homogeneous and anthropomorphic phantoms, respectively), and test (using 960 and 48 projections from the homogeneous and anthropomorphic phantoms, respectively) the DL model. The DL output was compared to the corresponding MC ground truth using both quantitative and qualitative metrics, such as mean relative and mean absolute relative differences (MRD and MARD), and to previously-published scatter-to-primary (SPR) ratios for similar breast phantoms. The scatter corrected DBT reconstructions were evaluated by analyzing the obtained linear attenuation values and by visual assessment of corrected projections in a clinical dataset. The time required for training and prediction per projection, as well as the time it takes to produce scatter-corrected projection images, were also tracked. The quantitative comparison between DL scatter predictions and MC simulations showed a median MRD of 0.05% (interquartile range (IQR), -0.04% to 0.13%) and a median MARD of 1.32% (IQR, 0.98% to 1.85%) for homogeneous phantom projections and a median MRD of -0.21% (IQR, -0.35% to -0.07%) and a median MARD of 1.43% (IQR, 1.32% to 1.66%) for the anthropomorphic phantoms. The SPRs for different breast thicknesses and at different projection angles were within±15% of the previously-published ranges. The visual assessment showed good prediction capabilities of the DL model with a close match between MC and DL scatter estimates, as well as between DL-based scatter corrected and anti-scatter grid corrected cases. The scatter correction improved the accuracy of the reconstructed linear attenuation of adipose tissue, reducing the error from -16% and -11% to -2.3% and 4.4% for an anthropomorphic digital phantom and clinical case with similar breast thickness, respectively. The DL model training took 40 min and prediction of a single projection took less than 0.01s. Generating scatter corrected images took 0.03s per projection for clinical exams and 0.16s for one entire projection set. This DL-based method for estimating the scatter signal in DBT projections is fast and accurate, paving the way for future quantitative applications.

Efficient Path‐Space Differentiable Volume Rendering With Respect To Shapes

AbstractDifferentiable rendering of translucent objects with respect to their shapes has been a long‐standing problem. State‐of‐the‐art methods require detecting object silhouettes or specifying change rates inside translucent objects—both of which can be expensive for translucent objects with complex shapes.In this paper, we address this problem for translucent objects with no refractive or reflective boundaries. By reparameterizing interior components of differential path integrals, our new formulation does not require change rates to be specified in the interior of objects. Further, we introduce new Monte Carlo estimators based on this formulation that do not require explicit detection of object silhouettes.

A Monte Carlo algorithm for the extrema of tempered stable processes

AbstractWe develop a novel Monte Carlo algorithm for the vector consisting of the supremum, the time at which the supremum is attained, and the position at a given (constant) time of an exponentially tempered Lévy process. The algorithm, based on the increments of the process without tempering, converges geometrically fast (as a function of the computational cost) for discontinuous and locally Lipschitz functions of the vector. We prove that the corresponding multilevel Monte Carlo estimator has optimal computational complexity (i.e. of order $\varepsilon^{-2}$ if the mean squared error is at most $\varepsilon^2$ ) and provide its central limit theorem (CLT). Using the CLT we construct confidence intervals for barrier option prices and various risk measures based on drawdown under the tempered stable (CGMY) model calibrated/estimated on real-world data. We provide non-asymptotic and asymptotic comparisons of our algorithm with existing approximations, leading to rule-of-thumb principles guiding users to the best method for a given set of parameters. We illustrate the performance of the algorithm with numerical examples.