Monte Carlo Estimation Research Articles

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.

Interrater reliability for multilevel data: A generalizability theory approach.

Current interrater reliability (IRR) coefficients ignore the nested structure of multilevel observational data, resulting in biased estimates of both subject- and cluster-level IRR. We used generalizability theory to provide a conceptualization and estimation method for IRR of continuous multilevel observational data. We explain how generalizability theory decomposes the variance of multilevel observational data into subject-, cluster-, and rater-related components, which can be estimated using Markov chain Monte Carlo (MCMC) estimation. We explain how IRR coefficients for each level can be derived from these variance components, and how they can be estimated as intraclass correlation coefficients (ICC). We assessed the quality of MCMC point and interval estimates with a simulation study, and showed that small numbers of raters were the main source of bias and inefficiency of the ICCs. In a follow-up simulation, we showed that a planned missing data design can diminish most estimation difficulties in these conditions, yielding a useful approach to estimating multilevel interrater reliability for most social and behavioral research. We illustrated the method using data on student-teacher relationships. All software code and data used for this article is available on the Open Science Framework: https://osf.io/bwk5t/. (PsycInfo Database Record (c) 2022 APA, all rights reserved).

Toise: a framework to describe the performance of high-energy neutrino detectors

Neutrinos offer a unique window to the distant, high-energy universe. Several next-generation instruments are being designed and proposed to characterize the flux of TeV–EeV neutrinos. The projected physics reach of the detectors is often quantified with simulation studies. However, a complete Monte Carlo estimate of detector performance is costly from a computational perspective, restricting the number of detector configurations considered when designing the instruments. In this paper, we present a new Python-based software framework, toise, which forecasts the performance of a high-energy neutrino detector using parameterizations of the detector performance, such as the effective areas, angular and energy resolutions, etc. The framework can be used to forecast performance of a variety of physics analyses, including sensitivities to diffuse fluxes of neutrinos and sensitivity to both transient and steady state point sources. This parameterized approach reduces the need for extensive simulation studies in order to estimate detector performance, and allows the user to study the influence of single performance metrics, like the angular resolution, in isolation. The framework is designed to allow for multiple detector components, each with different responses and exposure times, and supports paramterization of both optical- and radio-Cherenkov (Askaryan) neutrino telescopes. In the paper, we describe the mathematical concepts behind toise and introduce the reader to the use of the framework.

Multi-fidelity uncertainty quantification of particle deposition in turbulent pipe flow

Particle deposition in fully-developed turbulent pipe flow is quantified taking into account uncertainty in electric charge, van der Waals strength, and temperature effects. A framework is presented for obtaining variance-based sensitivity in multiphase flow systems via a multi-fidelity Monte Carlo approach that optimally manages model evaluations for a given computational budget. The approach combines a high-fidelity model based on direct numerical simulation and a lower-order model based on a one-dimensional Eulerian description of the two-phase flow. Significant speedup is obtained compared to classical Monte Carlo estimation. Deposition is found to be most sensitive to electrostatic interactions and exhibits largest uncertainty for mid-sized (i.e., moderate Stokes number) particles.