Standard Monte Carlo Techniques Research Articles

Differentiable stochastic halo occupation distribution

ABSTRACT In this work, we demonstrate how differentiable stochastic sampling techniques developed in the context of deep reinforcement learning can be used to perform efficient parameter inference over stochastic, simulation-based, forward models. As a particular example, we focus on the problem of estimating parameters of halo occupation distribution (HOD) models that are used to connect galaxies with their dark matter haloes. Using a combination of continuous relaxation and gradient re-parametrization techniques, we can obtain well-defined gradients with respect to HOD parameters through discrete galaxy catalogue realizations. Having access to these gradients allows us to leverage efficient sampling schemes, such as Hamiltonian Monte Carlo, and greatly speed up parameter inference. We demonstrate our technique on a mock galaxy catalogue generated from the Bolshoi simulation using a standard HOD model and find near-identical posteriors as standard Markov chain Monte Carlo techniques with an increase of ∼8× in convergence efficiency. Our differentiable HOD model also has broad applications in full forward model approaches to cosmic structure and cosmological analysis.

Complex Langevin approach to interacting Bose gases

Quantitative numerical analyses of interacting dilute Bose-Einstein condensates are most often based on semiclassical approximations. Since the complex-valued field-theoretic action of the Bose gas does not offer itself to standard Monte Carlo techniques, simulations beyond their scope almost exclusively rely on quantum-mechanical techniques, which are inherently limited to small particle numbers. Here we explore an alternative approach based on a Langevin-type sampling in an extended state space, known as complex Langevin (CL) algorithm. While the use of the CL technique has a long-standing history in high-energy physics, in particular in the simulation of QCD at finite baryon density, applications to ultracold atoms are still in their infancy. Here we examine the applicability of the CL approach for a one- and two-component, three-dimensional non-relativistic Bose gas in thermal equilibrium, below and above the Bose-Einstein phase transition. By comparison with analytic descriptions at the Gaussian level, including Bogoliubov and Hartree-Fock theory, we find that the method allows computing reliably and efficiently observables in the regime of experimentally accessible parameters. Close to the transition, quantum corrections lead to a shift of the critical temperature which we reproduce within the numerical range known in the literature. With this work, we aim to provide a first test of CL as a potential out-of-the-box tool for the simulation of experimentally realistic situations, including trapping geometries and multicomponent/-species models.

Sharp-edge diffraction under Bessel beam illumination: a catastrophe optics perspective.

Boundary diffraction wave theory and catastrophe optics have already proved to be a formidable combination in computational optics. In the present paper, a general paraxial theory aimed at dealing paraxial diffraction of Bessel beams by arbitrarily shaped sharp-edge apertures is developed. A key ingredient of our analysis is the δ-like nature of the angular spectrum of nondiffracting beams. This allows the diffracted wavefield to be effectively represented through two-dimensional integrals defined onto rectangular domains, whose numerical evaluation is easily achievable via standard Monte Carlo techniques. As a byproduct of the present analysis, a simple explanation of a recently observed property of some "heart-like" apertures to flatten the axial intensity of apodized Bessel beams is also provided.

Comparison of continuity equation and Gaussian mixture model for long-term density propagation using semi-analytical methods

This paper compares the continuum evolution for density equation modelling and the Gaussian mixture model on the 2D phase space long-term density propagation problem in the context of high-altitude and high area-to-mass ratio satellite long-term propagation. The density evolution equation, a pure numerical and pointwise method for the density propagation, is formulated under the influence of solar radiation pressure and Earth’s oblateness using semi-analytical methods. Different from the density evolution equation and Monte Carlo techniques, for the Gaussian mixture model, the analytical calculation of the density is accessible from the first two statistical moments (i.e. the mean and the covariance matrix) corresponding to each sub-Gaussian distribution for an initial Gaussian density distribution. An insight is given into the phase space long-term density propagation problem subject to nonlinear dynamics. The efficiency and validity of the density propagation are demonstrated and compared between the density evolution equation and the Gaussian mixture model with respect to standard Monte Carlo techniques.

Data-driven low-fidelity models for multi-fidelity Monte Carlo sampling in plasma micro-turbulence analysis

The linear micro-instabilities driving turbulent transport in magnetized fusion plasmas (as well as the respective nonlinear saturation mechanisms) are known to be sensitive with respect to various physical parameters characterizing the background plasma and the magnetic equilibrium. Therefore, uncertainty quantification is essential for achieving predictive numerical simulations of plasma turbulence. However, the high computational costs of the required gyrokinetic simulations and the large number of parameters render standard Monte Carlo techniques intractable. To address this problem, we propose a multi-fidelity Monte Carlo approach in which we employ data-driven low-fidelity models that exploit the structure of the underlying problem such as low intrinsic dimension and anisotropic coupling of the stochastic inputs. The low-fidelity models are efficiently constructed via sensitivity-driven dimension-adaptive sparse grid interpolation using both the full set of uncertain inputs and subsets comprising only selected, important parameters. We illustrate the power of this method by applying it to two plasma turbulence problems with up to 14 stochastic parameters, demonstrating that it is up to four orders of magnitude more efficient than standard Monte Carlo methods measured in single-core performance, which translates into a runtime reduction from around eight days to one hour on 240 cores on parallel machines.

A Discrete Quaternion Particle Filter Based on Deterministic Sampling for IMU Attitude Estimation

A discrete quaternion filtering algorithm following high-density uncertainty matching-based deterministic sampling scheme is proposed for IMU attitude estimation. Existing quaternion filters rely on the specific distributions directly or the on-tangent-plane Gaussian distribution indirectly to model the underlying uncertainty. The key idea of this paper is to apply a Dirac mixture at few deterministic weighted samples for non-parametric interpreting the high-density uncertainty of unit quaternion manifold. In contrast to standard Monte Carlo techniques based on random sampling, the proposed approach is deterministic. Furthermore, the deterministic quaternion particles are obtained by projecting on-sphere uniform points to the unit quaternion manifold via an exponential map. The on-manifold quaternions are then exploited in a discrete Bayesian filtering procedure. We compare the new filter with the genetic algorithm embedded quaternion particle filter (GA-QPF) and the unscented quaternion estimator (USQUE) for accuracy and computation cost by using simulated and practical IMU data. The results illustrate that the proposed estimator possesses lower computational cost than the GA-QPF while providing a comparable accuracy, higher accuracy than the USQUE, as well as more effective against large initialization error.

Computational analysis of conductive heat transfer in a rectangular slab of stable boundary using Monte Carlo method

Heat transfer is of immense importance in many engineering studies. Monte Carlo technique are broadly utilized in the operation research fields and atomic physics in which difficult problems above the existing tools of theoretical mathematics were resolved. The target of this study is to confirm the short duration in which heat transfer occurs in a quadrilateral slab where the temperature is provided throughout the borderline. An effort was made towards providing appropriate condition for the explanation of the heat transfer in a substance. The Shrinking Periphery Monte Carlo technique was utilized to obtain heat transfer in a helical pattern, upward and downward movement, which was compared with the standard Monte Carlo technique. The result of the study showed that the dimension of the quadrilateral slab determines the duration to calculate temperature dissemination in the system. The study revealed that the helical pattern is the shortest route in computing run time for temperature dispersal in a slab. The helical pattern is paramount in determining temperature distribution in a quadrilateral slab of stable state. The application of this technique to examine the conduction of heat in quadrilateral slabs produced good outcomes.

Bayesian Inference, Model Selection and Likelihood Estimation using Fast Rejection Sampling: The Conway-Maxwell-Poisson Distribution

Bayesian inference for models with intractable likelihood functions represents a challenging suite of problems in modern statistics. In this work we analyse the Conway-Maxwell-Poisson (COM-Poisson) distribution, a two parameter generalisation of the Poisson distribution. COM-Poisson regression modelling allows the flexibility to model dispersed count data as part of a generalised linear model (GLM) with a COM-Poisson response, where exogenous covariates control the mean and dispersion level of the response. The major difficulty with COM-Poisson regression is that the likelihood function contains multiple intractable normalising constants and is not amenable to standard inference and Markov Chain Monte Carlo (MCMC) techniques. Recent work by Chanialidis et al. (2018) has seen the development of a sampler to draw random variates from the COM-Poisson likelihood using a rejection sampling algorithm. We provide a new rejection sampler for the COM-Poisson distribution which significantly reduces the central processing unit (CPU) time required to perform inference for COM-Poisson regression models. An extension of this work shows that for any intractable likelihood function with an associated rejection sampler it is possible to construct unbiased estimators of the intractable likelihood which proves useful for model selection or for use within pseudo-marginal MCMC algorithms (Andrieu and Roberts, 2009). We demonstrate all of these methods on a real-world dataset of takeover bids.

Impact probability computation of near-Earth objects using Monte Carlo line sampling and subset simulation

This work introduces two Monte Carlo (MC)-based sampling methods, known as line sampling and subset simulation, to improve the performance of standard MC analyses in the context of asteroid impact risk assessment. Both techniques sample the initial uncertainty region in different ways, with the result of either providing a more accurate estimate of the impact probability or reducing the number of required samples during the simulation with respect to standard MC techniques. The two methods are first described and then applied to some test cases, providing evidence of the increased accuracy or the reduced computational burden with respect to a standard MC simulation. Finally, a sensitivity analysis is carried out to show how parameter setting affects the accuracy of the results and the numerical efficiency of the two methods.

A new deviational asymptotic preserving Monte Carlo method for the homogeneous Boltzmann equation

In this work, we introduce a new Monte Carlo method for solving the Boltzmann model of rarefied gas dynamics. The method works by reformulating the original problem through a micro-macro decomposition and successively in solving a suitable equation for the perturbation from the local thermodynamic equilibrium. This equation is then discretized by using unconditionally stable exponential schemes in time which project the solution over the corresponding equilibrium state when the time step is sent to infinity. The Monte Carlo method is designed on this time integration method and it only describes the perturbation from the final state. In this way, the number of samples diminishes during the time evolution of the solution and when the final equilibrium state is reached, the number of statistical samples becomes automatically zero. The resulting method is computationally less expensive as the solution approaches the equilibrium state as opposite to standard methods for kinetic equations which computational cost increases with the number of interactions. At the same time, the statistical error decreases as the system approaches the equilibrium state. In a last part, we show the behaviors of this new approach in comparison with standard Monte Carlo techniques and in comparison with spectral methods on different prototype problems. Mathematics Subject Classification: 76P05, 82C80, 82D05, 65C05, 35B25

Convergent series for polynomial lattice models with complex actions

Lattice models with complex actions are important for the understanding of matter at finite densities, but not accessible by the standard Monte Carlo techniques due to the sign problem. Here, we propose a new approach aiming to avoid the complex action/sign problem, by extending the method of convergent series (CS) with a non-Gaussian initial approximation. The main properties of the new series are demonstrated on the example of the two-dimensional oscillating integral.

Typical representatives of free homotopy classes in multi-punctured plane

We show that a uniform probability measure supported on a specific set of piecewise linear loops in a nontrivial free homotopy class in a multi-punctured plane is overwhelmingly concentrated around loops of minimal lengths. Our approach is based on extending Mogulskii’s theorem to closed paths, which is a useful result of independent interest. In addition, we show that the above measure can be sampled using standard Markov Chain Monte Carlo techniques, thus providing a simple method for approximating shortest loops.

English

The gamma ray log has over the years provided the conventional means for shale volume (Vsh) estimation. Knowledge of Vsh is used in the prediction of petrophysical parameters like effective porosity and water saturation, which are the input parameters for the calculation of oil in place. Currently, many studies have been conducted on the Bornu Basin of Nigeria, to access its hydrocarbon potential. Unfortunately, the practice of using best gamma ray log value for the computation of gamma ray index, IGR, and subsequently Vsh estimation incorporates huge uncertainty in the estimated volumes. Uncertainty is best captured when estimates are represented in a possible range rather than single value measurements. To the best of our knowledge, this is the first time shale volume has been estimated from the gamma ray log using sampling techniques. The gamma ray log data of the two upper shaly intervals of the NGAMMAEAST_1 well, which penetrates the Gombe formation of the basin, were utilized for this study. The gamma ray log response of the zone of interest is the uncertain parameter in Vsh estimation. A histogram plot of the uncertain log data was used to assume the probability distribution of the data. In the MATLAB platform, Standard Monte Carlo (MC) and Latin Hypercube sampling techniques were used to model the uncertain log response using random numbers. Possible input log data generated from the distribution of the uncertain log data were used in the linear and non-linear models for shale volume estimation to run a series of simulations to determine the possible range of estimates with their probabilities. The Latin hypercube method has shown to be a quick and accurate alternative method to the standard MC method. The approach presented here sets a guideline for the implementation of a probabilistic approach for the volume of shale estimation. Key words: Shale volume, Monte Carlo, Latin hypercube, sampling techniques, gamma ray log.

Asymptotic expansion for some local volatility models arising in finance

In this paper, we study the small noise asymptotic expansions for certain classes of local volatility models arising in finance. We provide explicit expressions for the involved coefficients as well as accurate estimates on the remainders. Moreover, we perform a detailed numerical analysis, with accuracy comparisons, of the obtained results by means of the standard Monte Carlo technique as well as exploiting the Polynomial Chaos Expansion approach.

Multi-scale control variate methods for uncertainty quantification in kinetic equations

Kinetic equations play a major rule in modeling large systems of interacting particles. Uncertainties may be due to various reasons, like lack of knowledge on the microscopic interaction details or incomplete informations at the boundaries. These uncertainties, however, contribute to the curse of dimensionality and the development of efficient numerical methods is a challenge. In this paper we consider the construction of novel multi-scale methods for such problems which, thanks to a control variate approach, are capable to reduce the variance of standard Monte Carlo techniques.

Configurational-Bias Monte Carlo Back-Mapping Algorithm for Efficient and Rapid Conversion of Coarse-Grained Water Structures into Atomistic Models

Coarse-grained molecular dynamics (MD) simulations represent a powerful approach to simulate longer time scale and larger length scale phenomena than those accessible to all-atom models. The gain in efficiency, however, comes at the cost of atomistic details. The reverse transformation, also known as back mapping, of coarse-grained beads into their atomistic constituents represents a major challenge. Most existing approaches are limited to specific molecules or specific force fields and often rely on running a long-time atomistic MD of the back-mapped configuration to arrive at an optimal solution. Such approaches are problematic when dealing with systems with high diffusion barriers. Here, we introduce a new extension of the configurational-bias Monte Carlo (CBMC) algorithm, which we term the crystalline-configurational-bias Monte Carlo (C-CBMC) algorithm, which allows rapid and efficient conversion of a coarse-grained model back into its atomistic representation. Although the method is generic, we use a coarse-grained water model as a representative example and demonstrate the back mapping or reverse transformation for model systems ranging from the ice-liquid water interface to amorphous and crystalline ice configurations. A series of simulations using the TIP4P/Ice model are performed to compare the new CBMC method to several other standard Monte Carlo and molecular dynamics-based back-mapping techniques. In all of the cases, the C-CBMC algorithm is able to find optimal hydrogen-bonded configuration many thousand evaluations/steps sooner than the other methods compared within this paper. For crystalline ice structures, such as a hexagonal, cubic, and cubic-hexagonal stacking disorder structures, the C-CBMC was able to find structures that were between 0.05 and 0.1 eV/water molecule lower in energy than the ground-state energies predicted by the other methods. Detailed analysis of the atomistic structures shows a significantly better global hydrogen positioning when contrasted with the existing simpler back-mapping methods. The errors in the radial distribution functions (RDFs) of back-mapped configuration relative to reference configuration for the C-CBMC, MD, and MC were found to be 6.9, 8.7, and 12.9, respectively, for the hexagonal system. For the cubic system, the relative errors of the RDFs for the C-CBMC, MD, and MC were found to be 18.2, 34.6, and 39.0, respectively. Our results demonstrate the efficiency and efficacy of our new back-mapping approach, especially for crystalline systems where simple force-field-based relaxations have a tendency to get trapped in local minima.

A Bayesian Genomic Regression Model with Skew Normal Random Errors.

Genomic selection (GS) has become a tool for selecting candidates in plant and animal breeding programs. In the case of quantitative traits, it is common to assume that the distribution of the response variable can be approximated by a normal distribution. However, it is known that the selection process leads to skewed distributions. There is vast statistical literature on skewed distributions, but the skew normal distribution is of particular interest in this research. This distribution includes a third parameter that drives the skewness, so that it generalizes the normal distribution. We propose an extension of the Bayesian whole-genome regression to skew normal distribution data in the context of GS applications, where usually the number of predictors vastly exceeds the sample size. However, it can also be applied when the number of predictors is smaller than the sample size. We used a stochastic representation of a skew normal random variable, which allows the implementation of standard Markov Chain Monte Carlo (MCMC) techniques to efficiently fit the proposed model. The predictive ability and goodness of fit of the proposed model were evaluated using simulated and real data, and the results were compared to those obtained by the Bayesian Ridge Regression model. Results indicate that the proposed model has a better fit and is as good as the conventional Bayesian Ridge Regression model for prediction, based on the DIC criterion and cross-validation, respectively. A computing program coded in the R statistical package and C programming language to fit the proposed model is available as supplementary material.

Performance analysis of waveguide-mode resonant optical filters with stochastic design parameters.

The performance of optical filters with resonant waveguide gratings is investigated numerically in a stochastic context, assuming random fluctuations of various design variables. Specifically, we derive stochastic models based on polynomial chaos expansions, whose involved coefficients are obtained by computing spectral projections via sparse-grid quadrature. The latter exploits purely deterministic results from a rigorous coupled-wave analysis solver and requires less simulation data than standard Monte Carlo (MC) techniques. The statistical moments of the filter's spectral response are calculated reliably, as the comparison against reference results from MC analysis verifies, and the extraction of the Sobol indices reveals the structure's sensitivity with respect to specific design parameters. Moreover, the present analysis clearly points out that neglecting even small geometric variations in the filter design may produce misleading conclusions regarding the corresponding performance, with undesirable consequences in real-life applications.

An efficient alternating direction method of multipliers for optimal control problems constrained by random Helmholtz equations

Based on the alternating direction method of multipliers (ADMM), we develop three numerical algorithms incrementally for solving the optimal control problems constrained by random Helmholtz equations. First, we apply the standard Monte Carlo technique and finite element method for the random and spatial discretization, respectively, and then ADMM is used to solve the resulting system. Next, combining the multi-modes expansion, Monte Carlo technique, finite element method, and ADMM, we propose the second algorithm. In the third algorithm, we preprocess certain quantities before the ADMM iteration, so that nearly no random variable is in the inner iteration. This algorithm is the most efficient one and is easy to implement. The error estimates of these three algorithms are established. The numerical experiments verify the efficiency of our algorithms.

Large Sample Properties of the Adaptive Gaussian Mixture Filter

In high-dimensional dynamic systems, standard Monte Carlo techniques that asymptotically reproduce the posterior distribution are computationally too expensive. Alternative sampling strategies are usually applied and among these the ensemble Kalman filter (EnKF) is perhaps the most popular. However, the EnKF suffers from severe bias if the model under consideration is far from linear. Another class of sequential Monte Carlo methods is kernel-based Gaussian mixture filters, which reduce the bias but maintain the robustness of the EnKF. Although many hybrid methods have been introduced in recent years, not many have been analyzed theoretically. Here it is shown that the recently proposed adaptive Gaussian mixture filter can be formulated in a rigorous Bayesian framework and that the algorithm can be generalized to a broader class of interpolated kernel filters. Two parameters—the bandwidth of the kernel and a weight interpolation factor—determine the filter performance. The new formulation of the filter includes particle filters, EnKF, and kernel-based Gaussian mixture filters as special cases. Techniques from particle filter literature are used to calculate the asymptotic bias of the filter as a function of the parameters and to derive a central limit theorem. The asymptotic theory is then used to determine the parameters as a function of the sample size in a robust way such that the error norm vanishes asymptotically, whereas the normalized error is sample independent and bounded. The parameter choice is tested on the Lorenz 63 model, where it is shown that the error is smaller or equal to the EnKF and the optimal particle filter for a varying sample size.