Monte Carlo Estimation Research Articles

Two multifidelity kriging-based strategies to control discretization error in reliability analysis exploiting a priori and a posteriori error estimators

This paper presents two approaches to tackle the issue of discretization error in the reliability assessment of structures. The first method (AGSK-MCS for Adaptive Guaranteed State Kriging Monte Carlo Sampling) uses discretization error bounds to guarantee the state safe or failed of the points used to build the Kriging metamodel of the limit state function. Two kriging metamodels interpolating lower and upper bounds can be constructed. These metamodels allow to compute discretization error bounds on the probability of failure through Monte Carlo sampling, which can then be used to validate the choice of the mesh. However, discretization error bounds are not available for any solver and any mechanical problem. In that case, a Mesh Size parameterized Kriging (MSK) metamodel can be used to check mesh convergence of the probability of failure. First, finite element simulations are spread on different mesh sizes. Second, the metamodel is used to compute the probability of failure for a given set of mesh sizes using Monte Carlo estimation. The mesh convergence of the probability of failure can be checked and may guide the user toward remeshing. These two strategies are illustrated on two 2-D mechanical problems.

Tissue equivalence of 3D printing materials with respect to attenuation and absorption of X-rays used for diagnostic and interventional imaging.

Three-dimensional printing is a promising technology to produce phantoms for quality assurance and dosimetry in X-ray imaging. Crucial to this, however, is the use of tissue equivalent printing materials. It was thus the aim of this study to evaluate the properties of a larger number of commercially available printing filaments with respect to their attenuation and absorption of X-rays. Apparent kerma attenuation coefficients (AKACs) and absorbed doses for different X-ray spectra (tube voltages, 70-140kV) were measured and simulated by Monte-Carlo computations for a larger number of fused-deposition-modeling (FDM) materials. The results were compared with the respective values simulated for reference body tissues. In addition, the properties of polylactide acid samples printed with reduced infill densities were investigated. Measured and simulated AKACs and absorbed doses agreed well with each other and in case of AKACs also with attenuation coefficients derived from the reference database of the National Institute of Standards and Technology (NIST). For lung, adipose, muscle, and bulk soft tissue as well as for spongiosa (cancellous bone), printed materials with equivalent attenuation as well as absorption properties could be identified. In contrast, none of the considered printed materials was equivalent to cortical bone. Several FDM materials have been identified as well-suited substitutes for body tissues in terms of the investigated material characteristics. They can therefore be used for in-house production of individualized and task-specific phantoms for image quality assessment and dose measurements in X-ray imaging.

Bayesian inference in spatial GARCH models: an application to US house price returns

ABSTRACT In this paper we consider a high-order spatial generalized autoregressive conditional heteroskedasticity (GARCH) model to account for the volatility clustering patterns observed over space. The model consists of a log-volatility equation that includes the high-order spatial lags of the log-volatility term and the squared outcome variable. We use a transformation approach to turn the model into a mixture of normals model, and then introduce a Bayesian Markov chain Monte Carlo (MCMC) estimation approach coupled with a data-augmentation technique. Our simulation results show that the Bayesian estimator has good finite sample properties. We apply a first-order version of the spatial GARCH model to US house price returns at the metropolitan statistical area level over the period 2006Q1–2013Q4 and show that there is significant variation in the log-volatility estimates over space in each period.

Level-crossing distributions of kinematic variables in multidirectional second-order ocean waves

The conditional value of a stationary random process, given the level-upcrossing of another dependent stationary random process, is considered. Assuming that both processes are weakly non-Gaussian, an analytical approximation for the related conditional distribution is derived. It is based on a trivariate Edgeworth expansion truncated to non-Gaussian terms of lowest order, to which Rice’s formula is then applied. As an application, the effect of level-upcrossing conditioning in second-order ocean waves is investigated. Upcrossing events are monitored for the sea surface elevation. The conditional distributions of different kinematic variables, given upcrossing, are considered for different sea-state configurations. Predictions from the analytical model are compared with numerical data obtained from Monte Carlo experiments. It is found that the analytical approximation provides conditional mean and variance in good agreement with numerical data, although moderate discrepancies appear for the sea states with the most severe wave steepnesses. Regarding the conditional skewness, given upcrossing, results are mixed, with significant discrepancies between the analytical approximation and numerical estimates, in a number of cases. An Edgeworth-type approximation is also provided for the upcrossing frequency and compared with Monte Carlo estimates; this analytical estimate is found to be accurate over a wide range of crossing levels.

Adsorption of atomic and molecular monolayers on Pt-supported graphene

Adsorption energies of monolayers of noble gases (Ar, Kr, Xe) and small molecules (N2, O2, CO, CH4, C2H6, and C3H8) on graphene supported by Pt(111) surface were computed by DFT, corrected for temperature and zero-point effects, and, compared to the recent temperature-programmed desorption experiments (Smith et al., 2016). The results indicate that some of the considered DFT exchange–correlation functionals reasonably well describe molecule–graphene/Pt(111) noncovalent interactions. For selected models, fixed-node diffusion Monte Carlo computations were used for ultimate cross-validation of interaction energies between the molecule and free-standing graphene.

Radar analysis algorithm for determining meteor head echo parameter probability distributions

ABSTRACT We present an automated radar data analysis algorithm developed to calculate probability distributions of meteor- and meteoroid parameters for head echoes detected with the Middle and Upper atmosphere (MU) radar in Shigaraki, Japan. The algorithm utilizes direct Monte Carlo simulations of uncertainties, with Bayesian Markov-chain Monte Carlo estimation of meteor model parameters and N-body propagation of distributions to perform orbit determination. The implementation has been validated using raw data simulations and a comparison with previous analysis methods. The concepts are applicable on a wide range of possible head echo measurements with other radar systems. The generated probability distributions provide quantitative reliability, which enables improved statistical studies and investigating the origins of detected meteoroids. The methodology section is highly detailed in order for the methods to be reproducible and provide a solid reference foundation for future studies. One such study is presented in a companion paper called ‘High-altitude meteors detected by the interferometric MU radar’.

Mixed vine copula flows for flexible modeling of neural dependencies.

Recordings of complex neural population responses provide a unique opportunity for advancing our understanding of neural information processing at multiple scales and improving performance of brain computer interfaces. However, most existing analytical techniques fall short of capturing the complexity of interactions within the concerted population activity. Vine copula-based approaches have shown to be successful at addressing complex high-order dependencies within the population, disentangled from the single-neuron statistics. However, most applications have focused on parametric copulas which bear the risk of misspecifying dependence structures. In order to avoid this risk, we adopted a fully non-parametric approach for the single-neuron margins and copulas by using Neural Spline Flows (NSF). We validated the NSF framework on simulated data of continuous and discrete types with various forms of dependency structures and with different dimensionality. Overall, NSFs performed similarly to existing non-parametric estimators, while allowing for considerably faster and more flexible sampling which also enables faster Monte Carlo estimation of copula entropy. Moreover, our framework was able to capture low and higher order heavy tail dependencies in neuronal responses recorded in the mouse primary visual cortex during a visual learning task while the animal was navigating a virtual reality environment. These findings highlight an often ignored aspect of complexity in coordinated neuronal activity which can be important for understanding and deciphering collective neural dynamics for neurotechnological applications.

Constructing unbiased gradient estimators with finite variance for conditional stochastic optimization

We study stochastic gradient descent for solving conditional stochastic optimization problems, in which an objective to be minimized is given by a parametric nested expectation with an outer expectation taken with respect to one random variable and an inner conditional expectation with respect to the other random variable. The gradient of such a parametric nested expectation is again expressed as a nested expectation, which makes it hard for the standard nested Monte Carlo estimator to be unbiased. In this paper, we show under some conditions that a multilevel Monte Carlo gradient estimator is unbiased and has finite variance and finite expected computational cost, so that the standard theory from stochastic optimization for a parametric (non-nested) expectation directly applies. We also discuss a special case for which yet another unbiased gradient estimator with finite variance and cost can be constructed.

A Simulation Study Comparing the Performance of Time-Varying Inverse Probability Weighting and G-Computation in Survival Analysis.

Inverse probability weighting (IPW) and g-computation are commonly used in time-varying analyses. To inform decisions on which to use, we compared these methods using a plasmode simulation based on data from the Effects of Aspirin in Gestation and Reproduction (EAGeR) Trial (June 15, 2007-July 15, 2011). In our main analysis, we simulated a cohort study of 1,226 individuals followed for up to 10 weeks. The exposure was weekly exercise, and the outcome was time to pregnancy. We controlled for 6 confounding factors: 4 baseline confounders (race, ever smoking, age, and body mass index) and 2 time-varying confounders (compliance with assigned treatment and nausea). We sought to estimate the average causal risk difference by 10 weeks, using IPW and g-computation implemented using a Monte Carlo estimator and iterated conditional expectations (ICE). Across 500 simulations, we compared the bias, empirical standard error (ESE), average standard error, standard error ratio, and 95% confidence interval coverage of each approach. IPW (bias = 0.02; ESE = 0.04; coverage = 92.6%) and Monte Carlo g-computation (bias = -0.01; ESE = 0.03; coverage = 94.2%) performed similarly. ICE g-computation was the least biased but least precise estimator (bias = 0.01; ESE = 0.06; coverage = 93.4%). When choosing an estimator, one should consider factors like the research question, the prevalences of the exposure and outcome, and the number of time points being analyzed.

Monte Carlo fPINNs: Deep learning method for forward and inverse problems involving high dimensional fractional partial differential equations

We introduce a sampling-based machine learning approach, Monte Carlo fractional physics-informed neural networks (MC-fPINNs), for solving forward and inverse fractional partial differential equations (FPDEs). As a generalization of the physics-informed neural networks (PINNs), MC-fPINNs utilize a Monte Carlo approximation strategy to compute the fractional derivatives of the DNN outputs, and construct an unbiased estimation of the physical soft constraints in the loss function. Our sampling approach can yield lower overall computational cost compared to fPINNs proposed in Pang et al.(2019), hence it can solve high dimensional FPDEs at reasonable cost. We validate the performance of MC-fPINNs via several examples, including high dimensional integral fractional Laplacian equations, parametric identification of time–space fractional PDEs, and fractional diffusion equation with random inputs. The results show that MC-fPINNs are flexible and quite effective in tackling high dimensional FPDEs.

Simulation of the drawdown and its duration in Lévy models via stick-breaking Gaussian approximation

We develop a computational method for expected functionals of the drawdown and its duration in exponential Lévy models. It is based on a novel simulation algorithm for the joint law of the state, supremum and time the supremum is attained of the Gaussian approximation for a general Lévy process. We bound the bias for various locally Lipschitz and discontinuous payoffs arising in applications and analyse the computational complexities of the corresponding Monte Carlo and multilevel Monte Carlo estimators. Monte Carlo methods for Lévy processes (using Gaussian approximation) have been analysed for Lipschitz payoffs, in which case the computational complexity of our algorithm is up to two orders of magnitude smaller when the jump activity is high. At the core of our approach are bounds on certain Wasserstein distances, obtained via the novel stick-breaking Gaussian (SBG) coupling between a Lévy process and its Gaussian approximation. Numerical performance, based on the implementation in Cázares and Mijatović (SBG approximation. GitHub repository. Available online at https://github.com/jorgeignaciogc/SBG.jl (2020)), exhibits a good agreement with our theoretical bounds. Numerical evidence suggests that our algorithm remains stable and accurate when estimating Greeks for barrier options and outperforms the “obvious” algorithm for finite-jump-activity Lévy processes.

Expectation-Maximization-Based Simultaneous Localization and Mapping for Millimeter-Wave Communication Systems.

In this paper, we proposed a novel expectation–maximization-based simultaneous localization and mapping (SLAM) algorithm for millimeter-wave (mmW) communication systems. By fully exploiting the geometric relationship among the access point (AP) positions, the angle difference of arrival (ADOA) from the APs and the mobile terminal (MT) position, and regarding the MT positions as the latent variable of the AP positions, the proposed algorithm first reformulates the SLAM problem as the maximum likelihood joint estimation over both the AP positions and the MT positions in a latent variable model. Then, it employs a feasible stochastic approximation expectation–maximization (EM) method to estimate the AP positions. Specifically, the stochastic Monte Carlo approximation is employed to obtain the intractable expectation of the MT positions’ posterior probability in the E-step, and the gradient descent-based optimization is used as a viable substitute for estimating the high-dimensional AP positions in the M-step. Further, it estimates the MT positions and constructs the indoor map based on the estimated AP topology. Due to the efficient processing capability of the stochastic approximation EM method and taking full advantage of the abundant spatial information in the crowd-sourcing ADOA data, the proposed method can achieve a better positioning and mapping performance than the existing geometry-based mmW SLAM method, which usually has to compromise between the computation complexity and the estimation performance. The simulation results confirm the effectiveness of the proposed algorithm.

Modeling Rapid Guessing Behaviors in Computer-Based Testlet Items.

In traditional test models, test items are independent, and test-takers slowly and thoughtfully respond to each test item. However, some test items have a common stimulus (dependent test items in a testlet), and sometimes test-takers lack motivation, knowledge, or time (speededness), so they perform rapid guessing (RG). Ignoring the dependence in responses to testlet items can negatively bias standard errors of measurement, and ignoring RG by fitting a simpler item response theory (IRT) model can bias the results. Because computer-based testing captures response times on testlet responses, we propose a mixture testlet IRT model with item responses and response time to model RG behaviors in computer-based testlet items. Two simulation studies with Markov chain Monte Carlo estimation using the JAGS program showed (a) good recovery of the item and person parameters in this new model and (b) the harmful consequences of ignoring RG (biased parameter estimates: overestimated item difficulties, underestimated time intensities, underestimated respondent latent speed parameters, and overestimated precision of respondent latent estimates). The application of IRT models with and without RG to data from a computer-based language test showed parameter differences resembling those in the simulations.

Cross-sector comovements and policy impact in the COVID-19 stock market: A dynamic factor approach

This study investigates the remarkable comovements in U.S. equity returns during the COVID-19 pandemic. It constructs a dynamic factor model (DFM) to illuminate the sources of the comovements and their implications. Using the Markov Chain Monte Carlo (MCMC) estimation method, the study finds that the comovements had a weak daily oscillation pattern during the pandemic. With that pattern, the study also finds significant monetary policy effects on the equity returns of several key sectors. In addition, it estimates the impact of news shocks, including monetary policy news, fiscal stimulus news, and unemployment news, on cross-sector equity returns. For any given sector, the conventional and unconventional monetary policy news shocked the sector in opposite directions. Among the positive monetary news shocks, the strongest were from interest rate policy surprises. Conversely, fiscal stimulus news had the most substantial positive impact and triggered all sectors to rebound from the bear market at the end of March 2020. Furthermore, by applying Natural Language Processing (NLP) sentiment analysis, this study sheds light on the positive correlation between comovements and news sentiment.

Perfect Sampling of the Posterior in the Hierarchical Pitman-Yor Process.

The predictive probabilities of the hierarchical Pitman-Yor process are approximated through Monte Carlo algorithms that exploits the Chinese Restaurant Franchise (CRF) representation. However, in order to simulate the posterior distribution of the hierarchical Pitman-Yor process, a set of auxiliary variables representing the arrangement of customers in tables of the CRF must be sampled through Markov chain Monte Carlo. This paper develops a perfect sampler for these latent variables employing ideas from the Propp-Wilson algorithm and evaluates its average running time by extensive simulations. The simulations reveal a significant dependence of running time on the parameters of the model, which exhibits sharp transitions. The algorithm is compared to simpler Gibbs sampling procedures, as well as a procedure for unbiased Monte Carlo estimation proposed by Glynn and Rhee. We illustrate its use with an example in microbial genomics studies.

How do Monte Carlo estimates affect stochastic geometric numerical integration?

In this work, we investigate the numerical conservation of characteristic properties of stochastic Hamiltonian problems driven by additive noise under time discretizations, in the spirit of stochastic geometric numerical integration. In particular, we aim to understand how Monte Carlo and multilevel Monte Carlo estimators for the expected values eventually affect the conservation properties of the numerical scheme. Specifically, our analysis is focused on the role of random number generation in the conservation of invariant laws for stochastic Hamiltonian problems under time discretization by the drift-preserving numerical scheme introduced in C. Chen, D. Cohen, R. D'Ambrosio and A. Lang, [Drift-preserving numerical integrators for stochastic Hamiltonian systems, Adv. Comput. Math. 46(2) (2020), pp. 27.] and the stochastic perturbation of a symplectic Runge–Kutta method, introduced in K. Burrage and P.M. Burrage, [Low rank Runge–Kutta methods, symplecticity and stochastic Hamiltonian problems with additive noise, J. Comput. Appl. Math. 236(16) (2012), pp. 3920–3930.]. The numerical evidence confirms the aforementioned theoretical analysis.

The nonzero gain coefficients of Sobol's sequences are always powers of two

When a plain Monte Carlo estimate on n samples has variance σ2/n, then scrambled digital nets attain a variance that is o(1/n) as n→∞. For finite n and an adversarially selected integrand, the variance of a scrambled (t,m,s)-net can be at most Γσ2/n for a maximal gain coefficient Γ<∞. The most widely used digital nets and sequences are those of Sobol'. It was previously known that Γ⩽2t3s for any nets in base 2. For digital nets, Dick and Pillichshammer (2010) obtained the bound 2t+s. In this paper we study digital nets in base 2 and show that Γ⩽2t+s−1 for such nets. This bound is a simple, but apparently unnoticed, consequence of a microstructure analysis by Niederreiter and Pirsic in 2001. We obtain a sharper bound that is smaller than this for some digital nets. Our main finding is that all nonzero gain coefficients must be powers of two. A consequence of this latter fact is a simplified algorithm for computing gain coefficients of digital nets in base 2.

Computational scatter correction in near real-time with a fast Monte Carlo photon transport model for high-resolution flat-panel CT

In computed tomography (CT), scattering causes server quality degradation of the reconstructed CT images by introducing streaks and cupping artifacts which reduce the detectability of low contrast objects. Monte Carlo (MC) simulation is considered the most accurate approach for scatter estimation. However, the existing MC estimators are computationally expensive, especially for high-resolution flat-panel CT. In this paper, we propose a fast and accurate MC photon transport model which describes the physics within the 1 keV to 1 MeV range using multiple controllable key parameters. Based on this model, scatter computation for a single projection can be completed within a range of a few seconds under well-defined model parameters. Smoothing and interpolation are performed on the estimated scatter to accelerate the scatter calculation without compromising accuracy too much compared to measured near scatter-free projection images. Combining the fast scatter estimation with the filtered backprojection (FBP), scatter correction is performed effectively in an iterative manner. To evaluate the proposed MC model, we have conducted extensive experiments on the simulated data and real-world high-resolution flat-panel CT. Compared to the state-of-the-art MC simulators, the proposed MC model achieved a 15times acceleration on a single-GPU compared to the GPU implementation of the Penelope simulator (MCGPU) utilizing several acceleration techniques, and a 202 times speed-up on a multi-GPU system compared to the multi-threaded state-of-the-art EGSnrc MC simulator. Furthermore, it is shown that for high-resolution images, scatter correction with sufficient accuracy is accomplished within one to three iterations using a FBP and the proposed fast MC photon transport model.

Rater characteristics, response content, and scoring contexts: Decomposing the determinates of scoring accuracy.

Raters may introduce construct-irrelevant variance when evaluating written responses to performance assessments, threatening the validity of students’ scores. Numerous factors in the rating process, including the content of students’ responses, the characteristics of raters, and the context in which the scoring occurs, are thought to influence the quality of raters’ scores. Despite considerable study of rater effects, little research has examined the relative impacts of the factors that influence rater accuracy. In practice, such integrated examinations are needed to afford evidence-based decisions of rater selection, training, and feedback. This study provides the first naturalistic, integrated examination of rater accuracy in a large-scale assessment program. Leveraging rater monitoring data from an English language arts (ELA) summative assessment program, I specified cross-classified, multilevel models via Bayesian (i.e., Markov chain Monte Carlo) estimation to decompose the impact of response content, rater characteristics, and scoring contexts on rater accuracy. Results showed relatively little variation in accuracy attributable to teams, items, and raters. Raters did not collectively exhibit differential accuracy over time, though there was significant variation in individual rater’s scoring accuracy from response to response and day to day. I found considerable variation in accuracy across responses, which was in part explained by text features and other measures of response content that influenced scoring difficulty. Some text features differentially influenced the difficulty of scoring research and writing content. Multiple measures of raters’ qualification performance predicted their scoring accuracy, but general rater background characteristics including experience and education did not. Site-based and remote raters demonstrated comparable accuracy, while evening-shift raters were slightly less accurate, on average, than day-shift raters. This naturalistic, integrated examination of rater accuracy extends previous research and provides implications for rater recruitment, training, monitoring, and feedback to improve human evaluation of written responses.

Adaptive Multilevel Monte Carlo for Probabilities

We consider the numerical approximation of $\mathbb{P}[G\in \Omega]$, where the $d$-dimensional random variable $G$ cannot be sampled directly, but there is a hierarchy of increasingly accurate approximations $\{G_\ell\}_{\ell\in\mathbb{N}}$ which can be sampled. The cost of standard Monte Carlo estimation scales poorly with accuracy in this setup since it compounds the approximation and sampling cost. A direct application of multilevel Monte Carlo (MLMC) improves this cost scaling slightly, but returns suboptimal computational complexities since estimation of the probability involves a discontinuous functional of $G_\ell$. We propose a general adaptive framework which is able to return the MLMC complexities seen for smooth or Lipschitz functionals of $G_\ell$. Our assumptions and numerical analysis are kept general allowing the methods to be used for a wide class of problems. We present numerical experiments on nested simulation for risk estimation, where $G = \mathbb{E}[X|Y]$ is approximated by an inner Monte Carlo estimate. Further experiments are given for digital option pricing, involving an approximation of a $d$-dimensional SDE.