High-dimensional Integrals Research Articles

Two-part models with stochastic processes for modelling longitudinal semicontinuous data: Computationally efficient inference and modelling the overall marginal mean.

Several researchers have described two-part models with patient-specific stochastic processes for analysing longitudinal semicontinuous data. In theory, such models can offer greater flexibility than the standard two-part model with patient-specific random effects. However, in practice, the high dimensional integrations involved in the marginal likelihood (i.e. integrated over the stochastic processes) significantly complicates model fitting. Thus, non-standard computationally intensive procedures based on simulating the marginal likelihood have so far only been proposed. In this paper, we describe an efficient method of implementation by demonstrating how the high dimensional integrations involved in the marginal likelihood can be computed efficiently. Specifically, by using a property of the multivariate normal distribution and the standard marginal cumulative distribution function identity, we transform the marginal likelihood so that the high dimensional integrations are contained in the cumulative distribution function of a multivariate normal distribution, which can then be efficiently evaluated. Hence, maximum likelihood estimation can be used to obtain parameter estimates and asymptotic standard errors (from the observed information matrix) of model parameters. We describe our proposed efficient implementation procedure for the standard two-part model parameterisation and when it is of interest to directly model the overall marginal mean. The methodology is applied on a psoriatic arthritis data set concerning functional disability.

Different numerical estimators for main effect global sensitivity indices

The variance-based method of global sensitivity indices based on Sobol' sensitivity indices became very popular among practitioners due to its easiness of interpretation. For complex practical problems computation of Sobol' indices generally requires a large number of function evaluations to achieve reasonable convergence. Four different direct formulas for computing Sobol’ main effect sensitivity indices are compared on a set of test models for which there are analytical results. Considered test functions represent various types of models that are found in practice. Formulas are based on high-dimensional integrals which are evaluated using Monte Carlo (MC) and Quasi Monte Carlo (QMC) techniques. Direct formulas are also compared with a different approach based on the so-called “double loop reordering” formula. It is found that the “double loop reordering” (DLR) approach shows a superior performance among all methods both for models with independent and dependent variables.

Copula-Based Random Effects Models for Clustered Data

Sorting and spillovers can create correlation in individual outcomes. In this situation, standard discrete choice estimators cannot consistently estimate the probability of joint and conditional events, and alternative estimators can yield incoherent statistical models or intractable estimators. I propose a random effects estimator that models the dependence among the unobserved heterogeneity of individuals in the same cluster using a parametric copula. This estimator makes it possible to compute joint and conditional probabilities of the outcome variable, and is statistically coherent. I describe its properties, establishing its efficiency relative to standard random effects estimators, and propose a specification test for the copula. The likelihood function for each cluster is an integral whose dimension equals the size of the cluster, which may require high-dimensional numerical integration. To overcome the curse of dimensionality from which methods like Monte Carlo integration suffer, I propose an algorithm that works for Archimedean copulas. I illustrate this approach by analysing labour supply in married couples.

Maximum likelihood estimation for stochastic volatility in mean models with heavy-tailed distributions.

In this article, we introduce a likelihood-based estimation method for the stochastic volatility in mean (SVM) model with scale mixtures of normal (SMN) distributions (Abanto-Valle et al., 2012). Our estimation method is based on the fact that the powerful hidden Markov model (HMM) machinery can be applied in order to evaluate an arbitrarily accurate approximation of the likelihood of an SVM model with SMN distributions. The method is based on the proposal of Langrock et al. (2012) and makes explicit the useful link between HMMs and SVM models with SMN distributions. Likelihood-based estimation of the parameters of stochastic volatility models in general, and SVM models with SMN distributions in particular, is usually regarded as challenging as the likelihood is a high-dimensional multiple integral. However, the HMM approximation, which is very easy to implement, makes numerical maximum of the likelihood feasible and leads to simple formulae for forecast distributions, for computing appropriately defined residuals, and for decoding, i.e., estimating the volatility of the process.

An offline approach for output-only Bayesian identification of stochastic nonlinear systems using unscented Kalman filtering

In this paper an offline approach for output-only Bayesian identification of stochastic nonlinear systems is presented. The approach is based on a re-parameterization of the joint posterior distribution of the parameters that define a postulated state-space stochastic model class. In the re-parameterization the state predictive distribution is included, marginalized, and estimated recursively in a state estimation step using an unscented Kalman filter, bypassing state augmentation as required by existing online methods. In applications expectations of functions of the parameters are of interest, which requires the evaluation of potentially high-dimensional integrals; Markov chain Monte Carlo is adopted to sample the posterior distribution and estimate the expectations. The proposed approach is suitable for nonlinear systems subjected to non-stationary inputs whose realization is unknown, and that are modeled as stochastic processes. Numerical verification and experimental validation examples illustrate the effectiveness and advantages of the approach, including: (i) an increased numerical stability with respect to augmented-state unscented Kalman filtering, avoiding divergence of the estimates when the forcing input is unmeasured; (ii) the ability to handle arbitrary prior and posterior distributions. The experimental validation of the approach is conducted using data from a large-scale structure tested on a shake table. It is shown that the approach is robust to inherent modeling errors in the description of the system and forcing input, providing accurate prediction of the dynamic response when the excitation history is unknown.

Low-rank canonical-tensor decomposition of potential energy surfaces: application to grid-based diagrammatic vibrational Green's function theory

ABSTRACTA new method is proposed for a fast evaluation of high-dimensional integrals of potential energy surfaces (PES) that arise in many areas of quantum dynamics. It decomposes a PES into a canonical low-rank tensor format, reducing its integral into a relatively short sum of products of low-dimensional integrals. The decomposition is achieved by the alternating least squares (ALS) algorithm, requiring only a small number of single-point energy evaluations. Therefore, it eradicates a force-constant evaluation as the hotspot of many quantum dynamics simulations and also possibly lifts the curse of dimensionality. This general method is applied to the anharmonic vibrational zero-point and transition energy calculations of molecules using the second-order diagrammatic vibrational many-body Green's function (XVH2) theory with a harmonic-approximation reference. In this application, high dimensional PES and Green's functions are both subjected to a low-rank decomposition. Evaluating the molecular integrals over a low-rank PES and Green's functions as sums of low-dimensional integrals using the Gauss–Hermite quadrature, this canonical-tensor-decomposition-based XVH2 (CT-XVH2) achieves an accuracy of 0.1 cm−1 or higher and nearly an order of magnitude speedup as compared with the original algorithm using force constants for water and formaldehyde.

Bayesian estimation of acoustic emissions source in plate structures using particle-based stochastic filtering

The application of particle-based stochastic filters to acoustic emission source localization in plate structures is presented. The approach employs time-of-flight measurements of guided waves using triangulation to estimate the acoustic emission source coordinates in a probabilistic framework using Bayesian inference, incorporating uncertainties related to material properties, measurement noise, and geometry of the system of interest. The estimate of the source location is given by a probability density function conditional on the guided wave measurements, found using particle-based stochastic simulation algorithms; in this setting, a set of particles is used to explore the space of possible source locations and efficiently estimate the posterior. The use of 2 filters is explored: the ensemble Kalman filter and the particle filter. The former filter assumes that the posterior distribution can be approximated by a Gaussian distribution, although the latter provides a nonparametric estimate of the posterior in the form of a weighted set of samples, overcoming the challenges related to the evaluation of high-dimensional integrals in an efficient way. Results of an experimental validation study conducted in a laboratory environment demonstrate the accuracy and efficiency of the particle filter-based approach. In particular, it is shown that the proposed particle filter-based approach has the capability to locate the emission source under minimal instrumentation, providing confidence intervals as a quantitative measure of the uncertainty in the estimates.

A Variational Maximization-Maximization Algorithm for Generalized Linear Mixed Models with Crossed Random Effects.

We present a variational maximization-maximization algorithm for approximate maximum likelihood estimation of generalized linear mixed models with crossed random effects (e.g., item response models with random items, random raters, or random occasion-specific effects). The method is based on a factorized variational approximation of the latent variable distribution given observed variables, which creates a lower bound of the log marginal likelihood. The lower bound is maximized with respect to the factorized distributions as well as model parameters. With the proposed algorithm, a high-dimensional intractable integration is translated into a two-dimensional integration problem. We incorporate an adaptive Gauss-Hermite quadrature method in conjunction with the variational method in order to increase computational efficiency. Numerical studies show that under the small sample size conditions that are considered the proposed algorithm outperforms the Laplace approximation.

Bayesian Computing with INLA: A Review

The key operation in Bayesian inference is to compute high-dimensional integrals. An old approximate technique is the Laplace method or approximation, which dates back to Pierre-Simon Laplace (1774). This simple idea approximates the integrand with a second-order Taylor expansion around the mode and computes the integral analytically. By developing a nested version of this classical idea, combined with modern numerical techniques for sparse matrices, we obtain the approach of integrated nested Laplace approximations (INLA) to do approximate Bayesian inference for latent Gaussian models (LGMs). LGMs represent an important model abstraction for Bayesian inference and include a large proportion of the statistical models used today. In this review, we discuss the reasons for the success of the INLA approach, the R-INLA package, why it is so accurate, why the approximations are very quick to compute, and why LGMs make such a useful concept for Bayesian computing.

Digital inversive vectors can achieve polynomial tractability for the weighted star discrepancy and for multivariate integration

We study high-dimensional numerical integration in the worst-case setting. The subject of tractability is concerned with the dependence of the worst-case integration error on the dimension. Roughly speaking, an integration problem is tractable if the worst-case error does not grow exponentially fast with the dimension. Many classical problems are known to be intractable. However, sometimes tractability can be shown. Often such proofs are based on randomly selected integration nodes. Of course, in applications, true random numbers are not available and hence one mimics them with pseudorandom number generators. This motivates us to propose the use of pseudorandom vectors as underlying integration nodes in order to achieve tractability. In particular, we consider digital inverse vectors and present two examples of problems, the weighted star discrepancy and integration of Hölder continuous, absolute convergent Fourier and cosine series, where the proposed method is successful.

Comparing the performance of linear and nonlinear principal components in the context of high-dimensional genomic data integration.

Linear principal component analysis (PCA) is a widely used approach to reduce the dimension of gene or miRNA expression data sets. This method relies on the linearity assumption, which often fails to capture the patterns and relationships inherent in the data. Thus, a nonlinear approach such as kernel PCA might be optimal. We develop a copula-based simulation algorithm that takes into account the degree of dependence and nonlinearity observed in these data sets. Using this algorithm, we conduct an extensive simulation to compare the performance of linear and kernel principal component analysis methods towards data integration and death classification. We also compare these methods using a real data set with gene and miRNA expression of lung cancer patients. First few kernel principal components show poor performance compared to the linear principal components in this occasion. Reducing dimensions using linear PCA and a logistic regression model for classification seems to be adequate for this purpose. Integrating information from multiple data sets using either of these two approaches leads to an improved classification accuracy for the outcome.

A new method to explore the structure of effective dimensions for functions

Effective dimension, an indicator for the difficulty of high-dimensional integration, describes whether a function can be well approximated by low-dimensional terms or sums of low-order terms. Some problems in option pricing are believed to have low effective dimensions, which help explain the success of quasi-Monte Carlo (QMC) methods recently observed in financial engineering. This paper provides a way of studying the structure of effective dimensions by finding a proper space the function of interest belongs to and then determining the effective dimension of that space. To this end, we extend the definitions of effective dimensions to weighted function spaces with product-order-dependent weights and give bounds on norms and variances. Furthermore, we show that the proposed method is applicable to functions arising in option pricing and consequently offers some hints on the performance of QMC methods.

Numerical Methods for Discrete Double Barrier Option Pricing Based on Merton Jump Diffusion Model

As a kind of weak-path dependent options, barrier options are an important kind of exotic options. Because the pricing formula for pricing barrier options with discrete observations cannot avoid computing a high dimensional integral, numerical calculation is time-consuming. In the current studies, some scholars just obtained theoretical derivation, or gave some simulation calculations. Others impose underlying assets on some strong assumptions, for example, a lot of calculations are based on the Black-Scholes model. This thesis considers Merton jump diffusion model as the basic model to derive the pricing formula of discrete double barrier option; numerical calculation method is used to approximate the continuous convolution by calculating discrete convolution. Then we compare the results of theoretical calculation with simulation results by Monte Carlo method, to verify their efficiency and accuracy. By comparing the results of degeneration constant parameter model with the results of previous models we verified the calculation method is correct indirectly. Compared with the Monte Carlo simulation method, the numerical results are stable. Even if we assume the simulation results are accurate, the time consumed by the numerical method to achieve the same accuracy is much less than the Monte Carlo simulation method.

A New Sparse Gauss-Hermite Cubature Rule Based on Relative-Weight-Ratios for Bearing-Ranging Target Tracking

A new sparse Gauss-Hermite cubature rule is designed to avoid dimension explosion caused by the traditional full tensor-product based Gauss-Hermite cubature rule. Although Smolyak’s quadrature rule can successfully generate sparse cubature points for high dimensional integral, it has a potential drawback that some cubature points generated by Smolyak’s rule have negative weights, which may result in instability for the computation. A relative-weight-ratio criterion based sparse Gauss-Hermite rule is presented in this paper, in which cubature points are kept symmetric in the input space and corresponding weights are guaranteed to be positive. The generation of the new sparse cubature points set is simple and meaningful for practice. The difference between our new sparse Gauss-Hermite cubature rule and other cubature rules is analysed. Simulation results show that, compared with Kalman filter with those types of full tensor-product based Gauss-Hermite rules, our new sparse Gauss-Hermite cubature rule based Kalman filter can lead to a substantially reduced number of cubature points, more stable computation capability, and maintaining the accuracy of integration at the same time.

Bayesian Option Pricing Framework with Stochastic Volatility for FX Data

The application of stochastic volatility (SV) models in the option pricing literature usually assumes that the market has sufficient option data to calibrate the model’s risk-neutral parameters. When option data are insufficient or unavailable, market practitioners must estimate the model from the historical returns of the underlying asset and then transform the resulting model into its risk-neutral equivalent. However, the likelihood function of an SV model can only be expressed in a high-dimensional integration, which makes the estimation a highly challenging task. The Bayesian approach has been the classical way to estimate SV models under the data-generating (physical) probability measure, but the transformation from the estimated physical dynamic into its risk-neutral counterpart has not been addressed. Inspired by the generalized autoregressive conditional heteroskedasticity (GARCH) option pricing approach by Duan in 1995, we propose an SV model that enables us to simultaneously and conveniently perform Bayesian inference and transformation into risk-neutral dynamics. Our model relaxes the normality assumption on innovations of both return and volatility processes, and our empirical study shows that the estimated option prices generate realistic implied volatility smile shapes. In addition, the volatility premium is almost flat across strike prices, so adding a few option data to the historical time series of the underlying asset can greatly improve the estimation of option prices.

Reducing the two-loop large-scale structure power spectrum to low-dimensional, radial integrals

Modeling the large-scale structure of the universe on nonlinear scales has the potential to substantially increase the science return of upcoming surveys by increasing the number of modes available for model comparisons. One way to achieve this is to model nonlinear scales perturbatively. Unfortunately, this involves high-dimensional loop integrals that are cumbersome to evaluate. Trying to simplify this, we show how two-loop (next-to-next-to-leading order) corrections to the density power spectrum can be reduced to low-dimensional, radial integrals. Many of those can be evaluated with a one-dimensional Fast Fourier Transform, which is significantly faster than the five-dimensional Monte-Carlo integrals that are needed otherwise. The general idea of this FFT-PT method is to switch between Fourier and position space to avoid convolutions and integrate over orientations, leaving only radial integrals. This reformulation is independent of the underlying shape of the initial linear density power spectrum and should easily accommodate features such as those from baryonic acoustic oscillations. We also discuss how to account for halo bias and redshift space distortions.

Monte Carlo approximation of likelihood function in spatial GLMMs through an empirical Bayes method

ABSTRACTIn spatial generalized linear mixed models (SGLMMs), statistical inference encounters problems, since random effects in the model imply high-dimensional integrals to calculate the marginal likelihood function. In this article, we temporarily treat parameters as random variables and express the marginal likelihood function as a posterior expectation. Hence, the marginal likelihood function is approximated using the obtained samples from the posterior density of the latent variables and parameters given the data. However, in this setting, misspecification of prior distribution of correlation function parameter and problems associated with convergence of Markov chain Monte Carlo (MCMC) methods could have an unpleasant influence on the likelihood approximation. To avoid these challenges, we utilize an empirical Bayes approach to estimate prior hyperparameters. We also use a computationally efficient hybrid algorithm by combining inverse Bayes formula (IBF) and Gibbs sampler procedures. A simulation study is conducted to assess the performance of our method. Finally, we illustrate the method applying a dataset of standard penetration test of soil in an area in south of Iran.

Bayesian zero-failure reliability modeling and assessment method for multiple numerical control (NC) machine tools

A new problem that classical statistical methods are incapable of solving is reliability modeling and assessment when multiple numerical control machine tools (NCMTs) reveal zero failures after a reliability test. Thus, the zero-failure data form and corresponding Bayesian model are developed to solve the zero-failure problem of NCMTs, for which no previous suitable statistical model has been developed. An expert-judgment process that incorporates prior information is presented to solve the difficulty in obtaining reliable prior distributions of Weibull parameters. The equations for the posterior distribution of the parameter vector and the Markov chain Monte Carlo (MCMC) algorithm are derived to solve the difficulty of calculating high-dimensional integration and to obtain parameter estimators. The proposed method is applied to a real case; a corresponding programming code and trick are developed to implement an MCMC simulation in WinBUGS, and a mean time between failures (MTBF) of 1057.9 h is obtained. Given its ability to combine expert judgment, prior information, and data, the proposed reliability modeling and assessment method under the zero failure of NCMTs is validated.

COMPUTING HIGH DIMENSIONAL INTEGRALS WITH MONTE CARLO METHODS

High dimensional integrals are usually solved with Monte Carlo algorithms and quasi Monte Carlo algorithms. We are doing numerical testing which compare low discrepancy and Monte Carlo algorithms. It is well known that Sobol algorithm has some advantageous over the other low discrepancy sequences, that’s why we use this algorithm for our numerical example. The obtained relative error confirms this superiority of the presented Monte Carlo and quasi Monte Carlo algorithms even when small number of sample points are used. It is very interesting that the presented high dimensional integral gives very low relative error even for computational time less than one second which shows the great importance of the developed algorithms.

Monte Carlo MP2 on Many Graphical Processing Units.

In the Monte Carlo second-order many-body perturbation (MC-MP2) method, the long sum-of-product matrix expression of the MP2 energy, whose literal evaluation may be poorly scalable, is recast into a single high-dimensional integral of functions of electron pair coordinates, which is evaluated by the scalable method of Monte Carlo integration. The sampling efficiency is further accelerated by the redundant-walker algorithm, which allows a maximal reuse of electron pairs. Here, a multitude of graphical processing units (GPUs) offers a uniquely ideal platform to expose multilevel parallelism: fine-grain data-parallelism for the redundant-walker algorithm in which millions of threads compute and share orbital amplitudes on each GPU; coarse-grain instruction-parallelism for near-independent Monte Carlo integrations on many GPUs with few and infrequent interprocessor communications. While the efficiency boost by the redundant-walker algorithm on central processing units (CPUs) grows linearly with the number of electron pairs and tends to saturate when the latter exceeds the number of orbitals, on a GPU it grows quadratically before it increases linearly and then eventually saturates at a much larger number of pairs. This is because the orbital constructions are nearly perfectly parallelized on a GPU and thus completed in a near-constant time regardless of the number of pairs. In consequence, an MC-MP2/cc-pVDZ calculation of a benzene dimer is 2700 times faster on 256 GPUs (using 2048 electron pairs) than on two CPUs, each with 8 cores (which can use only up to 256 pairs effectively). We also numerically determine that the cost to achieve a given relative statistical uncertainty in an MC-MP2 energy increases as O(n3) or better with system size n, which may be compared with the O(n5) scaling of the conventional implementation of deterministic MP2. We thus establish the scalability of MC-MP2 with both system and computer sizes.