Spatio-temporal modeling of particulate matter concentration through the SPDE approach

In this work, we consider a hierarchical spatio-temporal model for particulate matter (PM) concentration in the North-Italian region Piemonte. The model involves a Gaussian Field (GF), affected by a measurement error, and a state process characterized by a first order autoregressive dynamic model and spatially correlated innovations. This kind of model is well discussed and widely used in the air quality literature thanks to its flexibility in modelling the effect of relevant covariates (i.e. meteorological and geographical variables) as well as time and space dependence. However, Bayesian inference—through Markov chain Monte Carlo (MCMC) techniques—can be a challenge due to convergence problems and heavy computational loads. In particular, the computational issue refers to the infeasibility of linear algebra operations involving the big dense covariance matrices which occur when large spatio-temporal datasets are present. The main goal of this work is to present an effective estimating and spatial prediction strategy for the considered spatio-temporal model. This proposal consists in representing a GF with Matérn covariance function as a Gaussian Markov Random Field (GMRF) through the Stochastic Partial Differential Equations (SPDE) approach. The main advantage of moving from a GF to a GMRF stems from the good computational properties that the latter enjoys. In fact, GMRFs are defined by sparse matrices that allow for computationally effective numerical methods. Moreover, when dealing with Bayesian inference for GMRFs, it is possible to adopt the Integrated Nested Laplace Approximation (INLA) algorithm as an alternative to MCMC methods giving rise to additional computational advantages. The implementation of the SPDE approach through the R-library INLA (www.r-inla.org) is illustrated with reference to the Piemonte PM data. In particular, providing the step-by-step R-code, we show how it is easy to get prediction and probability of exceedance maps in a reasonable computing time.     

Approximating hidden Gaussian Markov random fields

Summary.  Gaussian Markov random-field (GMRF) models are frequently used in a wide variety of applications. In most cases parts of the GMRF are observed through mutually independent data; hence the full conditional of the GMRF, a hidden GMRF (HGMRF), is of interest. We are concerned with the case where the likelihood is non-Gaussian, leading to non-Gaussian HGMRF models. Several researchers have constructed block sampling Markov chain Monte Carlo schemes based on approximations of the HGMRF by a GMRF, using a second-order expansion of the log-density at or near the mode. This is possible as the GMRF approximation can be sampled exactly with a known normalizing constant. The Markov property of the GMRF approximation yields computational efficiency.The main contribution in the paper is to go beyond the GMRF approximation and to construct a class of non-Gaussian approximations which adapt automatically to the particular HGMRF that is under study. The accuracy can be tuned by intuitive parameters to nearly any precision. These non-Gaussian approximations share the same computational complexity as those which are based on GMRFs and can be sampled exactly with computable normalizing constants. We apply our approximations in spatial disease mapping and model-based geostatistical models with different likelihoods, obtain procedures for block updating and construct Metropolized independence samplers.     

Bayesian analysis of multivariate stochastic volatility with skew return distribution

Multivariate stochastic volatility models with skew distributions are proposed. Exploiting Cholesky stochastic volatility modeling, univariate stochastic volatility processes with leverage effect and generalized hyperbolic skew t-distributions are embedded to multivariate analysis with time-varying correlations. Bayesian modeling allows this approach to provide parsimonious skew structure and to easily scale up for high-dimensional problem. Analyses of daily stock returns are illustrated. Empirical results show that the time-varying correlations and the sparse skew structure contribute to improved prediction performance and Value-at-Risk forecasts.     

[HDDA] sparse subspace constrained partial least squares

ABSTRACT

In this paper, we investigate the objective function and deflation process for sparse Partial Least Squares (PLS) regression with multiple components. While many have considered variations on the objective for sparse PLS, the deflation process for sparse PLS has not received as much attention. Our work highlights a flaw in the Statistically Inspired Modification of Partial Least Squares (SIMPLS) deflation method when applied in sparse PLS regression. We also consider the Nonlinear Iterative Partial Least Squares (NIPALS) deflation in sparse PLS regression. To remedy the flaw in the SIMPLS method, we propose a new sparse PLS method wherein the direction vectors are constrained to be sparse and lie in a chosen subspace. We give insight into this new PLS procedure and show through examples and simulation studies that the proposed technique can outperform alternative sparse PLS techniques in coefficient estimation. Moreover, our analysis reveals a simple renormalization step that can be used to improve the estimation of sparse PLS direction vectors generated using any convex relaxation method.     

Modeling material stress using integrated Gaussian Markov random fields

The equations of a physical constitutive model for material stress within tantalum grains were solved numerically using a tetrahedrally meshed volume. The resulting output included a scalar vonMises stress for each of the more than 94,000 tetrahedra within the finite element discretization. In this paper, we define an intricate statistical model for the spatial field of vonMises stress which uses the given grain geometry in a fundamental way. Our model relates the three-dimensional field to integrals of latent stochastic processes defined on the vertices of the one- and two-dimensional grain boundaries. An intuitive neighborhood structure of the said boundary nodes suggested the use of a latent Gaussian Markov random field (GMRF). However, despite the potential for computational gains afforded by GMRFs, the integral nature of our model and the sheer number of data points pose substantial challenges for a full Bayesian analysis. To overcome these problems and encourage efficient exploration of the posterior distribution, a number of techniques are now combined: parallel computing, sparse matrix methods, and a modification of a block update strategy within the sampling routine. In addition, we use an auxiliary variables approach to accommodate the presence of outliers in the data.     

Testing the homogeneity of risk differences with sparse count data

ABSTRACT

In this paper, we consider testing the homogeneity of risk differences in independent binomial distributions especially when data are sparse. We point out some drawback of existing tests in either controlling a nominal size or obtaining powers through theoretical and numerical studies. The proposed test is designed to avoid the drawbacks of existing tests. We present the asymptotic null distribution and asymptotic power function for the proposed test. We also provide numerical studies including simulations and real data examples showing the proposed test has reliable results compared to existing testing procedures.     

Bayesian P-splines and advanced computing in R for a changepoint analysis on spatio-temporal point processes

ABSTRACT

This work presents advanced computational aspects of a new method for changepoint detection on spatio-temporal point process data. We summarize the methodology, based on building a Bayesian hierarchical model for the data and declaring prior conjectures on the number and positions of the changepoints, and show how to take decisions regarding the acceptance of potential changepoints. The focus of this work is about choosing an approach that detects the correct changepoint and delivers smooth reliable estimates in a feasible computational time; we propose Bayesian P-splines as a suitable tool for managing spatial variation, both under a computational and a model fitting performance perspective. The main computational challenges are outlined and a solution involving parallel computing in R is proposed and tested on a simulation study. An application is also presented on a data set of seismic events in Italy over the last 20 years.     

Computing and approximating multivariate chi-square probabilities

We consider computational methods for evaluating and approximating multivariate chi-square probabilities in cases where the pertaining correlation matrix or blocks thereof have a low factorial representation. To this end, techniques from matrix factorization and probability theory are applied. We outline a variety of statistical applications of multivariate chi-square distributions and provide a system of MATLAB programs implementing the proposed algorithms. Computer simulations demonstrate the accuracy and the computational efficiency of our methods in comparison with Monte Carlo approximations, and a real data example from statistical genetics illustrates their usage in practice.     

An improved algorithm for reliability bounds of multistate networks

Abstract

Indirect approaches based on minimal path vectors (d-MPs) and/or minimal cut vectors (d-MCs) are reported to be efficient for the reliability evaluation of multistate networks. Given the need to find more efficient evaluation methods for exact reliability, such techniques may still be cumbersome when the size of the network and the states of component are relatively large. Alternatively, computing reliability bounds can provide approximated reliability with less computational effort. Based on Bai’s exact and indirect reliability evaluation algorithm, an improved algorithm is proposed in this study, which provides sequences of upper and lower reliability bounds of multistate networks. Novel heuristic rules with a pre-specified value to filter less important sets of unspecified states are then developed and incorporated into the algorithm. Computational experiments comparing the proposed methods with an existing direct bounding algorithm show that the new algorithms can provide tight reliability bounds with less computational effort, especially for the proposed algorithm with heuristic L1.     

Markov chain Monte Carlo exact tests for incomplete two-way contingency tables

We consider testing the quasi-independence hypothesis for two-way contingency tables which contain some structural zero cells. For sparse contingency tables where the large sample approximation is not adequate, the Markov chain Monte Carlo exact tests are powerful tools. To construct a connected chain over the two-way contingency tables with fixed sufficient statistics and an arbitrary configuration of structural zero cells, an algebraic algorithm proposed by Diaconis and Sturmfels [Diaconis, P. and Sturmfels, B. (1998). The Annals of statistics, 26, pp. 363–397.] can be used. However, their algorithm does not seem to be a satisfactory answer, because the Markov basis produced by the algorithm often contains many redundant elements and is hard to interpret. We derive an explicit characterization of a minimal Markov basis, prove its uniqueness, and present an algorithm for obtaining the unique minimal basis. A computational example and the discussion on further basis reduction for the case of positive sufficient statistics are also given.     

Iterative numerical methods for sampling from high dimensional Gaussian distributions

Many applications require efficient sampling from Gaussian distributions. The method of choice depends on the dimension of the problem as well as the structure of the covariance- (Σ) or precision matrix (Q). The most common black-box routine for computing a sample is based on Cholesky factorization. In high dimensions, computing the Cholesky factor of Σ or Q may be prohibitive due to accumulation of more non-zero entries in the factor than is possible to store in memory. We compare different methods for computing the samples iteratively adapting ideas from numerical linear algebra. These methods assume that matrix vector products, Qv, are fast to compute. We show that some of the methods are competitive and faster than Cholesky sampling and that a parallel version of one method on a Graphical Processing Unit (GPU) using CUDA can introduce a speed-up of up to 30x. Moreover, one method is used to sample from the posterior distribution of petroleum reservoir parameters in a North Sea field, given seismic reflection data on a large 3D grid.     

Regression with outlier shrinkage

We propose a robust regression method called regression with outlier shrinkage (ROS) for the traditional n>p$n>p$ cases. It improves over the other robust regression methods such as least trimmed squares (LTS) in the sense that it can achieve maximum breakdown value and full asymptotic efficiency simultaneously. Moreover, its computational complexity is no more than that of LTS. We also propose a sparse estimator, called sparse regression with outlier shrinkage (SROS), for robust variable selection and estimation. It is proven that SROS can not only give consistent selection but also estimate the nonzero coefficients with full asymptotic efficiency under the normal model. In addition, we introduce a concept of nearly regression equivariant estimator for understanding the breakdown properties of sparse estimators, and prove that SROS achieves the maximum breakdown value of nearly regression equivariant estimators. Numerical examples are presented to illustrate our methods.     

A new semiparametric Weibull cure rate model: fitting different behaviors within GAMLSS

ABSTRACT

We propose a new semiparametric Weibull cure rate model for fitting nonlinear effects of explanatory variables on the mean, scale and cure rate parameters. The regression model is based on the generalized additive models for location, scale and shape, for which any or all distribution parameters can be modeled as parametric linear and/or nonparametric smooth functions of explanatory variables. We present methods to select additive terms, model estimation and validation, where all computational codes are presented in a simple way such that any R user can fit the new model. Biases of the parameter estimates caused by models specified erroneously are investigated through Monte Carlo simulations. We illustrate the usefulness of the new model by means of two applications to real data. We provide computational codes to fit the new regression model in the R software.     

The construction of a partial least-squares biplot

Biplots are useful tools to explore the relationship among variables. In this paper, the specific regression relationship between a set of predictors X and set of response variables Y by means of partial least-squares (PLS) regression is represented. The PLS biplot provides a single graphical representation of the samples together with the predictor and response variables, as well as their interrelationships in terms of the matrix of regression coefficients.     

A representation for quadrivariate normal positive orthant probabilities

For four variables x₁,x₂, x₃ and x₄, which have a quadrivariate normal distribution with means equal to zero, the positive ortrhant probability is the probability that all of the x.'s are simultaneously positive. A representation for the quadrivariate normal positive orthant probability is obtained and it is a function of no more than three integrals over a single variable. Extensive testing has shown this representation to be very efficient on a computational basis.     

Randomized heuristic algorithms for orthogonal projection of a point onto a set

Abstract

The problem of orthogonal projection of a point onto a set is an essential problem of computational geometry. This problem has many practical applications in different areas such as robotics, computer graphics and so on. In the present paper three algorithms for solving this problem are proposed. This algorithms are based on the idea of heuristic random search. Numerical experiments illustrating the work of the proposed methods are presented.     

Integral Representation of a Distribution Associated with a Pure Birth Process

ABSTRACT

This article deals with a distribution associated with a pure birth process starting with no individuals, with birth rates λ _n  = λ for n = 0, 2,…, m ? 1 and λ _n  = μ for n ≥ m. The probability mass function is expressed in terms of an integral that is very convenient for computing probabilities, moments, generating functions, and others. Using this representation, the kth factorial moments of the distribution is obtained. Some other forms of this distribution are also given.     

Sparse Approximations of Fractional Matérn Fields

We consider fast lattice approximation methods for a solution of a certain stochastic non‐local pseudodifferential operator equation. This equation defines a Matérn class random field. We approximate the pseudodifferential operator with truncated Taylor expansion, spectral domain error functional minimization and rounding approximations. This allows us to construct Gaussian Markov random field approximations. We construct lattice approximations with finite‐difference methods. We show that the solutions can be constructed with overdetermined systems of stochastic matrix equations with sparse matrices, and we solve the system of equations with a sparse Cholesky decomposition. We consider convergence of the truncated Taylor approximation by studying band‐limited Matérn fields. We consider the convergence of the discrete approximations to the continuous limits. Finally, we study numerically the accuracy of different approximation methods with an interpolation problem.     

A Pointwise Estimator for the k-Fold Convolution of a Distribution Function

ABSTRACT

This article is concerned with some parametric and nonparametric estimators for the k-fold convolution of a distribution function. An alternative estimator is proposed and its unbiasedness, asymptotic unbiasedness, and consistency properties are investigated. The asymptotic normality of this estimator is established. Some applications of the estimator are given in renewal processes. Finally, the computational procedures are described and the relative performance of these estimators for small sample sizes is investigated by a simulation study.