Search

In this paper, we show the impact of Pascal on numerical software. Using recursion, a new and elegant algorithm for quadrature could be developed.

A method is introduced for approximate marginal likelihood inference via adaptive Gaussian quadrature in mixed models with a single grouping factor. The core technical contribution is an algorithm for computing the exact gradient of the approximate log-marginal likelihood. This leads to efficient maximum likelihood via quasi-Newton optimization that is demonstrated to be faster than existing approaches based on finite-differenced gradients or derivative-free optimization. The method is specialized to Bernoulli mixed models with multivariate, correlated Gaussian random effects; here computations are performed using an inverse log-Cholesky parameterization of the Gaussian density that involves no matrix decomposition during model fitting, while Wald confidence intervals are provided for variance parameters on the original scale. Simulations give evidence of these intervals attaining nominal coverage if enough quadrature points are used, for data comprised of a large number of very small groups exhibiting large between-group heterogeneity. The Laplace approximation is well-known to give especially poor coverage and high bias for data comprised of a large number of small groups. Adaptive quadrature mitigates this, and the methods in this paper improve the computational feasibility of this more accurate method. All results may be reproduced using code available at . [ABSTRACT FROM AUTHOR]

Physics-informed neural networks (PINNs) and their variants have recently emerged as alternatives to traditional partial differential equation (PDE) solvers, but little literature has focused on devising accurate numerical integration methods for neural networks (NNs), which is essential for getting accurate solutions. In this work, we propose adaptive quadratures for the accurate integration of neural networks and apply them to loss functions appearing in low-dimensional PDE discretisations. We show that at opposite ends of the spectrum, continuous piecewise linear (CPWL) activation functions enable one to bound the integration error, while smooth activations ease the convergence of the optimisation problem. We strike a balance by considering a CPWL approximation of a smooth activation function. The CPWL activation is used to obtain an adaptive decomposition of the domain into regions where the network is almost linear, and we derive an adaptive global quadrature from this mesh. The loss function is then obtained by evaluating the smooth network (together with other quantities, e.g., the forcing term) at the quadrature points. We propose a method to approximate a class of smooth activations by CPWL functions and show that it has a quadratic convergence rate. We then derive an upper bound for the overall integration error of our proposed adaptive quadrature. The benefits of our quadrature are evaluated on a strong and weak formulation of the Poisson equation in dimensions one and two. Our numerical experiments suggest that compared to Monte-Carlo integration, our adaptive quadrature makes the convergence of NNs quicker and more robust to parameter initialisation while needing significantly fewer integration points and keeping similar training times. [ABSTRACT FROM AUTHOR]

We provide the first stochastic convergence rates for a family of adaptive quadrature rules used to normalize the posterior distribution in Bayesian models. Our results apply to the uniform relative error in the approximate posterior density, the coverage probabilities of approximate credible sets, and approximate moments and quantiles, therefore, guaranteeing fast asymptotic convergence of approximate summary statistics used in practice. The family of quadrature rules includes adaptive Gauss-Hermite quadrature, and we apply this rule in two challenging low-dimensional examples. Further, we demonstrate how adaptive quadrature can be used as a crucial component of a modern approximate Bayesian inference procedure for high-dimensional additive models. The method is implemented and made publicly available in the aghq package for the R language, available on CRAN. for this article are available online. [ABSTRACT FROM AUTHOR]

The Cognitive Reflection Test (CRT) is a test designed to assess subjects' ability to override intuitively appealing but incorrect responses. Psychologists are concerned with whether subjects improve their scores on the test with repeated exposure, in which case, the test's predictive validity may be threatened. In this paper, we take a novel approach to modelling data recorded on subjects who took the CRT multiple times. We develop bivariate, longitudinal models to describe the responses, CRT score and time taken to complete the CRT. These responses serve as a proxy for the underlying latent variables "numeracy" and "reflectiveness", respectively---two components of "rationality". Our models allow for subpopulations of individuals whose responses exhibit similar patterns. We assess the reasonableness of our models via new visualizations of the data. We estimate their parameters by modifying the method of adaptive Gaussian quadrature. We then use our fitted models to address a range of subject-specific questions in a formal way. We find evidence of at least three subpopulations, which we interpret as representing individuals with differing combinations of numeracy and reflectiveness, and determine that, in some subpopulations, test exposure has a greater estimated effect on test scores than previously reported.

In this note, a rule for double integrals of precision-7 is investigated using precision-5 rules on different integrals with the help of adaptive technique for any vector function using Greens theorem for approximate evaluation of line integral. The rule is reasonable coincident with the single rule (5 point Clenshaw-Curtis rule). Convergence analysis as well as numerical experiments is also carried out through several examples. [ABSTRACT FROM AUTHOR]

In analysis of problems with parametric spline boundaries that are immersed or inserted into an underlying domain, the discretization on the underlying domain usually does not conform to the inserted boundaries. While the fixed underlying discretization is of great convenience as the immersed boundaries evolve, the field approximations near the inserted boundaries require refinement in the underlying domain, as do the quadrature cells. In this paper, a kd-tree data structure together with a sign-based and/or distance-based refinement strategy is proposed for local refinement near the inserted boundaries as well as for adaptive quadrature near the boundaries. The developed algorithms construct and utilize implicit forms of parametric Non-Uniform Rational B-Spline (NURBS) surfaces to algebraically (and non-iteratively) estimate distance as well as sign relative to the inserted boundary. The kd-tree local refinement is demonstrated to produce fewer sub-cells for the same accuracy of solution as compared to the classical quad/oct tree-based subdivision. Consistent with the kd-tree data structure, we describe a new a priori refinement algorithm based on the signed and unsigned distance from the inserted boundary. We first demonstrate the local refinement strategy coupled with the the kd-tree data structure by constructing Truncated Hierarchical B-spline (THB-spline) "meshes". We next demonstrate the accuracy and efficiency of the developed local refinement strategy through adaptive quadrature near NURBS boundaries inserted within volumetric three-dimensional NURBS discretizations. [ABSTRACT FROM AUTHOR]

The R package CFC performs cause-specific, competing-risk survival analysis by computing cumulative incidence functions from unadjusted, cause-specific survival functions. A high-level API in CFC enables end-to-end survival and competing-risk analysis, using a single-line function call, based on the parametric survival regression models in the survival package. A low-level API allows users to achieve more flexibility by supplying their custom survival functions, perhaps in a Bayesian setting. Utility methods for summarizing and plotting the output allow population-average cumulative incidence functions to be calculated, visualized and compared to unadjusted survival curves. Numerical and computational optimization strategies are employed for efficient and reliable computation of the coupled integrals involved. To address potential integrable singularities caused by infinite cause-specific hazards, particularly near time-from-index of zero, integrals are transformed to remove their dependency on hazard functions, making them solely functions of causespecific, unadjusted survival functions. This implicit variable transformation also provides for easier extensibility of CFC to handle custom survival models since it only requires the users to implement a maximum of one function per cause. The transformed integrals are numerically calculated using a generalization of Simpson's rule to handle the implicit change of variable from time to survival, while a generalized trapezoidal rule is used as reference for error calculation. An OpenMP-parallelized, efficient C++ implementation - using packages Rcpp and RcppArmadillo - makes the application of CFC in Bayesian settings practical, where a potentially large number of samples represent the posterior distribution of cause-specific survival functions.

In analysis of problems with parametric spline boundaries that are immersed or inserted into an underlying domain, the discretization on the underlying domain usually does not conform to the inserted boundaries. While the fixed underlying discretization is of great convenience as the immersed boundaries evolve, the field approximations near the inserted boundaries require refinement in the underlying domain, as do the quadrature cells. In this paper, a kd-tree data structure together with a sign-based and/or distance-based refinement strategy is proposed for local refinement near the inserted boundaries as well as for adaptive quadrature near the boundaries. The developed algorithms construct and utilize implicit forms of parametric Non-Uniform Rational B-Spline (NURBS) surfaces to algebraically (and non-iteratively) estimate distance as well as sign relative to the inserted boundary. The kd-tree local refinement is demonstrated to produce fewer sub-cells for the same accuracy of solution as compared to the classical quad/oct tree-based subdivision. Consistent with the kd-tree data structure, we describe a new a priori refinement algorithm based on the signed and unsigned distance from the inserted boundary. We first demonstrate the local refinement strategy coupled with the the kd-tree data structure by constructing Truncated Hierarchical B-spline (THB-spline) “meshes”. We next demonstrate the accuracy and efficiency of the developed local refinement strategy through adaptive quadrature near NURBS boundaries inserted within volumetric three-dimensional NURBS discretizations.

Generalized linear mixed models or latent variable models for categorical data are difficult to estimate if the random effects or latent variables vary at non-nested levels, such as persons and test items. Clayton and Rasbash (1999) suggested an Alternating Imputation Posterior (AIP) algorithm for approximate maximum likelihood estimation. For item response models with random item effects, the algorithm iterates between an item wing in which the item mean and variance are estimated for given person effects and a person wing in which the person mean and variance are estimated for given item effects. The person effects used for the item wing are sampled from the conditional posterior distribution estimated in the person wing and vice versa. Clayton and Rasbash (1999) used marginal quasi-likelihood (MQL) and penalized quasi-likelihood (PQL) estimation within the AIP algorithm, but this method has been shown to produce biased estimates in many situations, so we use maximum likelihood estimation with adaptive quadrature. We apply the proposed algorithm to the famous salamander mating data, comparing the estimates with many other methods, and to an educational testing dataset. We also present a simulation study to assess performance of the AIP algorithm and the Laplace approximation with different numbers of items and persons and a range of item and person variances. © 2010 Elsevier B.V. All rights reserved.

In mathematical finance, a process of calibrating stochastic volatility (SV) option pricing models to real market data involves a numerical calculation of integrals that depend on several model parameters. This optimization task consists of large number of integral evaluations with high precision and low computational time requirements. However, for some model parameters, many numerical quadrature algorithms fail to meet these requirements. We can observe an enormous increase in function evaluations, serious precision problems and a significant increase of computational time. In this paper, we numerically analyse these problems and show that they are especially caused by inaccurately evaluated integrands. We propose a fast regime switching algorithm that tells if it is sufficient to evaluate the integrand in standard double arithmetic or if a higher precision arithmetic has to be used. We compare and recommend numerical quadratures for typical SV models and different parameter values, especially for problematic cases. [ABSTRACT FROM AUTHOR]

A mixed quadrature rule of precision nine for approximate evaluation of real definite integrals has been constructed by blending Clenshaw-Curtis five point rule modified by Richardson Extrapolation and Gauss- Legendre four point rule. An error analysis for this mixed rule is provided. The efficiency of this rule is highlighted through numerical evaluation of some definite integrals at the end. [ABSTRACT FROM AUTHOR]

We discuss prediction of random effects and of expected responses in multilevel generalized linear models. Prediction of random effects is useful for instance in small area estimation and disease mapping, effectiveness studies and model diagnostics. Prediction of expected responses is useful for planning, model interpretation and diagnostics. For prediction of random effects, we concentrate on empirical Bayes prediction and discuss three different kinds of standard errors; the posterior standard deviation and the marginal prediction error standard deviation (comparative standard errors) and the marginal sampling standard deviation (diagnostic standard error). Analytical expressions are available only for linear models and are provided in an appendix. For other multilevel generalized linear models we present approximations and suggest using parametric bootstrapping to obtain standard errors. We also discuss prediction of expectations of responses or probabilities for a new unit in a hypothetical cluster, or in a new (randomly sampled) cluster or in an existing cluster. The methods are implemented in gllamm and illustrated by applying them to survey data on reading proficiency of children nested in schools. Simulations are used to assess the performance of various predictions and associated standard errors for logistic random-intercept models under a range of conditions. © 2009 Royal Statistical Society.

Studying behavior in economics, sociology, and statistics often involvesfitting models in which the response variable depends on a dummy variable- also known as a regime-switch variable- or in which the response variable is observed only if a particular selection condition is met. In either case, standard regression techniques deliver inconsistent estimators if unobserved factors that affect the re- sponse are correlated with unobserved factors that affect the switching or selection variable. Consistent estimators can be obtained by maximum likelihood estimation of a joint model of the outcome and switching or selection variable. This article describes a “wrapper” program, ssm, that calls gllamm (Rabe-Hesketh, Skrondal, and Pickles, GLLAMM Manual [University of California – Berkeley, Division of Bio- statistics, Working Paper Series, Paper No. 160]) to fit such models. The wrapper accepts data in a simple structure, has a straightforward syntax, and reports out- put that is easily interpretable. One important feature of ssm is that the log likelihood can be evaluated using adaptive quadrature (Rabe-Hesketh, Skrondal, and Pickles, Stata Journal 2: 1–21; Journal of Econometrics 128: 301–323). Copyright 2006 by StataCorp LP.

Gauss-Hermite quadrature is often used to evaluate and maximize the likelihood for random component probit models. Unfortunately, the estimates are biased for large cluster sizes and/or intraclass correlations. We show that adaptive quadrature largely overcomes these problems. We then extend the adaptive quadrature approach to general random coefficient models with limited and discrete dependent variables. The models can include several nested random effects (intercepts and coefficients) representing unobserved heterogeneity at different levels of a hierarchical dataset. The required multivariate integrals are evaluated efficiently using spherical quadrature rules. Simulations show that adaptive quadrature performs well in a wide range of situations. © 2004 Published by Elsevier B.V.

A unifying framework for generalized multilevel structural equation modeling is introduced. The models in the framework, called generalized linear latent and mixed models (GLLAMM), combine features of generalized linear mixed models (GLMM) and structural equation models (SEM) and consist of a response model and a structural model for the latent variables. The response model generalizes GLMMs to incorporate factor structures in addition to random intercepts and coefficients. As in GLMMs, the data can have an arbitrary number of levels and can be highly unbalanced with different numbers of lower-level units in the higher-level units and missing data. A wide range of response processes can be modeled including ordered and unordered categorical responses, counts, and responses of mixed types. The structural model is similar to the structural part of a SEM except that it may include latent and observed variables varying at different levels. For example, unit-level latent variables (factors or random coefficients) can be regressed on cluster-level latent variables. Special cases of this framework are explored and data from the British Social Attitudes Survey are used for illustration. Maximum likelihood estimation and empirical Bayes latent score prediction within the GLLAMM framework can be performed using adaptive quadrature in gllamm, a freely available program running in Stata.

Generalized linear mixed models or multilevel regression models have become increasingly popular. Several methods have been proposed for estimating such models. However, to date there is no single method that can be assumed to work well in all circumstances in terms of both parameter recovery and computational efficiency. Stata's xt commands for two-level generalized linear mixed models (e.g., xtlogit) employ Gauss–Hermite quadrature to evaluate and maximize the marginal log likelihood. The method generally works very well, and often better than common contenders such as MQL and PQL, but there are cases where quadrature performs poorly. Adaptive quadrature has been suggested to overcome these problems in the two-level case. We have recently implemented a multilevel versionofthismethodin gllamm, a program that fits a large class of multilevel latent variable models including multilevel generalized linear mixed models. As far as we know, this is the first time that adaptive quadrature has been proposed for multilevel models. We show that adaptive quadrature works well in problems where ordinary quadrature fails. Furthermore, even when ordinary quadrature works, adaptive quadrature is often computationally more efficient since it requires fewer quadrature points to achieve the same precision.

• hp-adaptive quadrature increases the accuracy of numerical integration. • Better error control for the population balance source terms. • Irregular integrands are efficiently dealt with. • Combination of piece wise Gaussian quadrature with weighted residual methods. The solution of the population balance equation requires the integration of several source terms. In the numerical weighted residuals methods, Gaussian quadrature is a natural candidate for numerical integration. Previous works using the weighted residuals methods for solving the population balance equation did use a fixed grid of quadrature points. This work shows that the use of adaptive quadrature points for the numerical integration can lead to more efficient and accurate solutions of the equation. For cases where the integrand shows a high degree of irregularity, the hp-optimization method distributes the quadrature points such that the method becomes more efficient than with a fixed grid. An additional improvement is that the amount of quadrature points changes to fit the need for each integral present, rather than having one set of quadrature points for all cases. A simple population balance model demonstrates the use of the adaptive quadrature approach. [ABSTRACT FROM AUTHOR]

In this paper, we study the deep Ritz method for solving the linear elasticity equation from a numerical analysis perspective. A modified Ritz formulation using the H 1 / 2 (Γ D) norm is introduced and analyzed for linear elasticity equation in order to deal with the (essential) Dirichlet boundary condition. We show that the resulting deep Ritz method provides the best approximation among the set of deep neural network (DNN) functions with respect to the "energy" norm. Furthermore, we demonstrate that the total error of the deep Ritz simulation is bounded by the sum of the network approximation error and the numerical integration error, disregarding the algebraic error. To effectively control the numerical integration error, we propose an adaptive quadrature-based numerical integration technique with a residual-based local error indicator. This approach enables efficient approximation of the modified energy functional. Through numerical experiments involving smooth and singular problems, as well as problems with stress concentration, we validate the effectiveness and efficiency of the proposed deep Ritz method with adaptive quadrature. [ABSTRACT FROM AUTHOR]

A lattice Boltzmann method for radiation transfer and Newton–Cotes formulas are used in this work to compute the propagation of polychromatic light in a biosuspension of phototrophic microorganisms. The polychromatic light field is obtained from monochromatic lattice Boltzmann simulations by integration across the visible spectrum. The effects of the spectral resolution, radiation characteristics and the chosen integration rule on the accuracy of the integration are investigated. It was found that reasonable results can be achieved on equidistant spectral grids with a grid spacing of Δλ ≤ 20 nm, although error compensation might be a serious issue if the trapezoidal rule is applied. Based on a priori information about the light field, an approach for the computation of adapted spectral grids is introduced, which aims at the efficient computation of polychromatic light fields. It was found that no significant increase of accuracy can be realized by usage of adapted spectral grids for spectral integration. It is presumed that this observation is caused by the changing shape of the light spectrum along the optical path. [ABSTRACT FROM AUTHOR]

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. [ABSTRACT FROM AUTHOR]

A novel collision source method has been developed to solve the Linear Boltzmann Equation (LBE) more efficiently by adaptation of the angular quadrature order. The angular adaptation method is unique in that the flux from each scattering source iteration is obtained, with potentially a different quadrature order used for each. Traditionally, the flux from every iteration is combined, with the same quadrature applied to the combined flux. Since the scattering process tends to distribute the radiation more evenly over angles (i.e., make it more isotropic), the quadrature requirements generally decrease with each iteration. This method allows for an optimal use of processing power, by using a high order quadrature for the first iterations that need it, before shifting to lower order quadratures for the remaining iterations. This is essentially an extension of the first collision source method, and is referred to as the adaptive collision source (ACS) method. The ACS methodology has been implemented in the 3-D, parallel, multigroup discrete ordinates code TITAN. This code was tested on a several simple and complex fixed-source problems. The ACS implementation in TITAN has shown a reduction in computation time by a factor of 1.5–4 on the fixed-source test problems, for the same desired level of accuracy, as compared to the standard TITAN code. [ABSTRACT FROM AUTHOR]

We propose a novel approach to the problem of polynomial approximation of rational Bézier triangular patches with prescribed boundary control points. The method is very efficient thanks to using recursive properties of the bivariate dual Bernstein polynomials and applying a smart algorithm for evaluating a collection of two-dimensional integrals. Some illustrative examples are given. [ABSTRACT FROM AUTHOR]

In this paper, we introduce a mixed quadrature rule using Fejer's second rule and Gaussian rule. This rule is taken as the base rule to develop an adaptive integration scheme. Using this scheme, some test integrals have been evaluated. The results are found to be more encouraging as compared to those obtained by using some other quadratures in this integration scheme. [ABSTRACT FROM AUTHOR]

The excess hazard regression model is an approach developed for the analysis of cancer registry data to estimate net survival, that is, the survival of cancer patients that would be observed if cancer was the only cause of death. Cancer registry data typically possess a hierarchical structure: individuals from the same geographical unit share common characteristics such as proximity to a large hospital that may influence access to and quality of health care, so that their survival times might be correlated. As a consequence, correct statistical inference regarding the estimation of net survival and the effect of covariates should take this hierarchical structure into account. It becomes particularly important as many studies in cancer epidemiology aim at studying the effect on the excess mortality hazard of variables, such as deprivation indexes, often available only at the ecological level rather than at the individual level. We developed here an approach to fit a flexible excess hazard model including a random effect to describe the unobserved heterogeneity existing between different clusters of individuals, and with the possibility to estimate non-linear and time-dependent effects of covariates. We demonstrated the overall good performance of the proposed approach in a simulation study that assessed the impact on parameter estimates of the number of clusters, their size and their level of unbalance. We then used this multilevel model to describe the effect of a deprivation index defined at the geographical level on the excess mortality hazard of patients diagnosed with cancer of the oral cavity. Copyright © 2016 John Wiley & Sons, Ltd. [ABSTRACT FROM AUTHOR]

Heap-based priority queues are very common dynamical data structures used in several fields, ranging from operating systems to scientific applications. However, the rise of new multicore CPUs introduced new challenges in the process of design of these data structures: in addition to traditional requirements like correctness and progress, the scalability is of paramount importance. It is a common opinion that these two demands are partially in conflict each other, so that in these computational environments it is necessary to relax the requirements of correctness and linearizability to achieve high performances. In this paper we introduce a loosely coordinated approach for the management of heap based priority queues on multicore CPUs, with the aim to realize a tradeoff between efficiency and sequential correctness. The approach is based on a sharing of information among only a small number of cores, so that to improve performance without completely losing the features of the data structure. The results obtained on a scientific problem show significant benefits both in terms of parallel efficiency, as well as in term of numerical accuracy. [ABSTRACT FROM AUTHOR]

We investigate the possibility of fast, accurate and reliable computation of the Cauchy principal value integrals ⨍ a b f ( x ) ( x − τ ) − 1 d x ( a < τ < b ) using a standard adaptive quadrature. In order to properly control the error tolerance for the adaptive quadrature and to obtain a reliable estimation of the approximation error, we research the possible influence of round-off errors on the computed result. As the numerical experiments confirm, the proposed method can successfully compete with other algorithms for computing such type integrals. Moreover, the presented method is very easy to implement on any system equipped with a reliable adaptive integration subroutine. [ABSTRACT FROM AUTHOR]

A common problem in the longitudinal data analysis is the missing data problem. Two types of missing patterns are generally considered in statistical literature: monotone and non-monotone missing data. Nonmonotone missing data occur when study participants intermittently miss scheduled visits, while monotone missing data can be from discontinued participation, loss to follow-up, and mortality. Although many novel statistical approaches have been developed to handle missing data in recent years, few methods are available to provide inferences to handle both types of missing data simultaneously. In this article, a latent random effects model is proposed to analyze longitudinal outcomes with both monotone and non-monotone missingness in the context of missing not at random. Another significant contribution of this article is to propose a new computational algorithm for latent random effects models. To reduce the computational burden of high-dimensional integration problem in latent random effects models, we develop a new computational algorithm that uses a new adaptive quadrature approach in conjunction with the Taylor series approximation for the likelihood function to simplify the E-step computation in the expectation-maximization algorithm. Simulation study is performed and the data from the scleroderma lung study are used to demonstrate the effectiveness of this method. [ABSTRACT FROM AUTHOR]