Geostatistical Modelling Using Non‐Gaussian Matérn Fields

This work provides a class of non‐Gaussian spatial Matérn fields which are useful for analysing geostatistical data. The models are constructed as solutions to stochastic partial differential equations driven by generalized hyperbolic noise and are incorporated in a standard geostatistical setting with irregularly spaced observations, measurement errors and covariates. A maximum likelihood estimation technique based on the Monte Carlo expectation‐maximization algorithm is presented, and a Monte Carlo method for spatial prediction is derived. Finally, an application to precipitation data is presented, and the performance of the non‐Gaussian models is compared with standard Gaussian and transformed Gaussian models through cross‐validation.     

Scalable statistical inference for averaged implicit stochastic gradient descent

Stochastic gradient descent (SGD) provides a scalable way to compute parameter estimates in applications involving large‐scale data or streaming data. As an alternative version, averaged implicit SGD (AI‐SGD) has been shown to be more stable and more efficient. Although the asymptotic properties of AI‐SGD have been well established, statistical inferences based on it such as interval estimation remain unexplored. The bootstrap method is not computationally feasible because it requires to repeatedly resample from the entire data set. In addition, the plug‐in method is not applicable when there is no explicit covariance matrix formula. In this paper, we propose a scalable statistical inference procedure, which can be used for conducting inferences based on the AI‐SGD estimator. The proposed procedure updates the AI‐SGD estimate as well as many randomly perturbed AI‐SGD estimates, upon the arrival of each observation. We derive some large‐sample theoretical properties of the proposed procedure and examine its performance via simulation studies.     

Predictive Inference for Big,Spatial, Non‐Gaussian Data: MODIS Cloud Data and its Change‐of‐Support

Remote sensing of the earth with satellites yields datasets that can be massive in size, nonstationary in space, and non‐Gaussian in distribution. To overcome computational challenges, we use the reduced‐rank spatial random effects (SRE) model in a statistical analysis of cloud‐mask data from NASA's Moderate Resolution Imaging Spectroradiometer (MODIS) instrument on board NASA's Terra satellite. Parameterisations of cloud processes are the biggest source of uncertainty and sensitivity in different climate models’ future projections of Earth's climate. An accurate quantification of the spatial distribution of clouds, as well as a rigorously estimated pixel‐scale clear‐sky‐probability process, is needed to establish reliable estimates of cloud‐distributional changes and trends caused by climate change. Here we give a hierarchical spatial‐statistical modelling approach for a very large spatial dataset of 2.75 million pixels, corresponding to a granule of MODIS cloud‐mask data, and we use spatial change‐of‐Support relationships to estimate cloud fraction at coarser resolutions. Our model is non‐Gaussian; it postulates a hidden process for the clear‐sky probability that makes use of the SRE model, EM‐estimation, and optimal (empirical Bayes) spatial prediction of the clear‐sky‐probability process. Measures of prediction uncertainty are also given.     

Efficient Bayesian Multivariate Surface Regression

Methods for choosing a fixed set of knot locations in additive spline models are fairly well established in the statistical literature. The curse of dimensionality makes it nontrivial to extend these methods to nonadditive surface models, especially when there are more than a couple of covariates. We propose a multivariate Gaussian surface regression model that combines both additive splines and interactive splines, and a highly efficient Markov chain Monte Carlo algorithm that updates all the knot locations jointly. We use shrinkage prior to avoid overfitting with different estimated shrinkage factors for the additive and surface part of the model, and also different shrinkage parameters for the different response variables. Simulated data and an application to firm leverage data show that the approach is computationally efficient, and that allowing for freely estimated knot locations can offer a substantial improvement in out‐of‐sample predictive performance.     

A calibrated imputation method for secondary data analysis of survey data

In practical survey sampling, missing data are unavoidable due to nonresponse, rejected observations by editing, disclosure control, or outlier suppression. We propose a calibrated imputation approach so that valid point and variance estimates of the population (or domain) totals can be computed by the secondary users using simple complete‐sample formulae. This is especially helpful for variance estimation, which generally require additional information and tools that are unavailable to the secondary users. Our approach is natural for continuous variables, where the estimation may be either based on reweighting or imputation, including possibly their outlier‐robust extensions. We also propose a multivariate procedure to accommodate the estimation of the covariance matrix between estimated population totals, which facilitates variance estimation of the ratios or differences among the estimated totals. We illustrate the proposed approach using simulation data in supplementary materials that are available online.     

Non‐Gaussian geostatistical modeling using (skew) t processes

We propose a new model for regression and dependence analysis when addressing spatial data with possibly heavy tails and an asymmetric marginal distribution. We first propose a stationary process with t marginals obtained through scale mixing of a Gaussian process with an inverse square root process with Gamma marginals. We then generalize this construction by considering a skew‐Gaussian process, thus obtaining a process with skew‐t marginal distributions. For the proposed (skew) t process, we study the second‐order and geometrical properties and in the t case, we provide analytic expressions for the bivariate distribution. In an extensive simulation study, we investigate the use of the weighted pairwise likelihood as a method of estimation for the t process. Moreover we compare the performance of the optimal linear predictor of the t process versus the optimal Gaussian predictor. Finally, the effectiveness of our methodology is illustrated by analyzing a georeferenced dataset on maximum temperatures in Australia.     

M-Estimates of regression when the scale is unknown and the error distribution is possibly asymmetric: A minimax result

Huber (1964) found the minimax-variance M-estimate of location under the assumption that the scale parameter is known; Li and Zamar (1991) extended this result to the case when the scale is unknown. We consider the robust estimation of the regression coefficients (β₁,…,β_p) when the scale and the intercept parameters are unknown. The minimax-variance estimates of (β₁,…,β_p) with respect to the trace of their asymptotic covariance matrix are derived. The maximum is taken over ?-contamination neighbourhoods of a central regression model with Gaussian errors (asymmetric contamination is allowed), and the minimum is taken over a large class of generalized M-estimates of regression of the Mallow type. The optimal choice of estimates for the nuisance parameters (scale and intercept) is also considered.     

A Non‐Parametric Estimator of the Spectral Density of a Continuous‐Time Gaussian Process Observed at Random Times

Abstract. In numerous applications data are observed at random times and an estimated graph of the spectral density may be relevant for characterizing and explaining phenomena. By using a wavelet analysis, one derives a non‐parametric estimator of the spectral density of a Gaussian process with stationary increments (or a stationary Gaussian process) from the observation of one path at random discrete times. For every positive frequency, this estimator is proved to satisfy a central limit theorem with a convergence rate depending on the roughness of the process and the moment of random durations between successive observations. In the case of stationary Gaussian processes, one can compare this estimator with estimators based on the empirical periodogram. Both estimators reach the same optimal rate of convergence, but the estimator based on wavelet analysis converges for a different class of random times. Simulation examples and an application to biological data are also provided.     

Estimation of functions of the parameters of a normal distribution

Several authors have considered the problem of estimating parameters of a distribution after some fixed Gaussian inducing transformation has been applied to the observations. This paper extends this work to the situation where the observations represent a noisy version of a true process, the parameters of the latter requiring estimation     

M‐estimation for general ARMA Processes with Infinite Variance

Abstract. General autoregressive moving average (ARMA) models extend the traditional ARMA models by removing the assumptions of causality and invertibility. The assumptions are not required under a non‐Gaussian setting for the identifiability of the model parameters in contrast to the Gaussian setting. We study M‐estimation for general ARMA processes with infinite variance, where the distribution of innovations is in the domain of attraction of a non‐Gaussian stable law. Following the approach taken by Davis et al. (1992) and Davis (1996) , we derive a functional limit theorem for random processes based on the objective function, and establish asymptotic properties of the M‐estimator. We also consider bootstrapping the M‐estimator and extend the results of Davis & Wu (1997) to the present setting so that statistical inferences are readily implemented. Simulation studies are conducted to evaluate the finite sample performance of the M‐estimation and bootstrap procedures. An empirical example of financial time series is also provided.     

Bayesian analysis of the scatterometer wind retrieval inverse problem: some new approaches

Summary.  The retrieval of wind vectors from satellite scatterometer observations is a non-linear inverse problem. A common approach to solving inverse problems is to adopt a Bayesian framework and to infer the posterior distribution of the parameters of interest given the observations by using a likelihood model relating the observations to the parameters, and a prior distribution over the parameters. We show how Gaussian process priors can be used efficiently with a variety of likelihood models, using local forward (observation) models and direct inverse models for the scatterometer. We present an enhanced Markov chain Monte Carlo method to sample from the resulting multimodal posterior distribution. We go on to show how the computational complexity of the inference can be controlled by using a sparse, sequential Bayes algorithm for estimation with Gaussian processes. This helps to overcome the most serious barrier to the use of probabilistic, Gaussian process methods in remote sensing inverse problems, which is the prohibitively large size of the data sets. We contrast the sampling results with the approximations that are found by using the sparse, sequential Bayes algorithm.     

Bayesian parameter estimation via variational methods

We consider a logistic regression model with a Gaussian prior distribution over the parameters. We show that an accurate variational transformation can be used to obtain a closed form approximation to the posterior distribution of the parameters thereby yielding an approximate posterior predictive model. This approach is readily extended to binary graphical model with complete observations. For graphical models with incomplete observations we utilize an additional variational transformation and again obtain a closed form approximation to the posterior. Finally, we show that the dual of the regression problem gives a latent variable density model, the variational formulation of which leads to exactly solvable EM updates.     

Improved robust Bayes estimators of the error variance in linear models

We consider the problem of estimating the error variance in a general linear model when the error distribution is assumed to be spherically symmetric, but not necessary Gaussian. In particular we study the case of a scale mixture of Gaussians including the particularly important case of the multivariate-t distribution. Under Stein's loss, we construct a class of estimators that improve on the usual best unbiased (and best equivariant) estimator. Our class has the interesting double robustness property of being simultaneously generalized Bayes (for the same generalized prior) and minimax over the entire class of scale mixture of Gaussian distributions.     

The Self‐Similar and Multifractal Nature of a Network Traffic Model

Abstract

We look at a family of models for Internet traffic with increasing input rates and consider approximation models which exhibit self‐similarity at large time scales and multifractality at small time scales. Depending on whether the input rate is fast or slow, the total cumulative input traffic can be approximated by a self‐similar stable Lévy motion or a self‐similar Gaussian process. The stable Lévy limit does not depend on the behavior of the individual transmission schedules but the Gaussian limit does. Also, the models and their approximations show multifractal behavior at small time scales.     

Use of likelihood ratio tests to detect outliers under the variance shift outlier model

In this paper, we revisit the alternative outlier model of Thompson [A note on restricted maximum likelihood estimation with an alternative outlier model, J. Roy. Stat. Soc. Ser. B 47 (1985), pp. 53–55] for detecting outliers in the linear model. Gumedze et al. [A variance shift model for detection of outliers in the linear mixed model, Comput. Statist. Data Anal. 54 (2010), pp. 2128–2144] called this model the variance shift outlier model (VSOM). The basic idea behind the VSOM is to detect observations with inflated variance and isolate them for further investigation. The VSOM is appealing because it downweights an outlier in the analysis, with the weighting determined automatically as part of the estimation procedure. We set up the VSOM as a linear mixed model and then use the likelihood ratio test (LRT) statistic as an objective measure for determining whether the weighting is required, i.e. whether the observation is an outlier. We also derived one-step updates of the variance parameter estimates based on observed, expected and average information matrices to obtain one-step LRT statistics which usually require less computation. Both the fully iterated and one-step LRTs are functions of the squared standard residuals from the null model and therefore can be computed directly without the need to fit the VSOM. We investigated the properties of the likelihood ratio tests and compare them. An extension of the model to detect a group of outliers is also given. We illustrate the proposed methodology using simulated datasets and a real dataset.     

The auto-regression and the moving-average

We explore some relationships in the second-order properties of a causal auto-regression and an invertible moving-average process with the same polynomial. We reveal that the inverse variance matrix for random variables from the auto-regression is equal to a conditional variance matrix of Gaussian random variables from the moving-average and vice versa. While the inverse variance matrix for the auto-regression can be written explicitly, we manage to write down the exact Gaussian likelihood of consecutive observations from the moving-average process, by using the properties of the auto-regression.     

Diagnosing explainable heterogeneity of variance in random‐effects models

Data‐analytic tools for models other than the normal linear regression model are relatively rare. Here we develop plots and diagnostic statistics for nonconstant variance for the random‐effects model (REM). REMs for longitudinal data include both within‐ and between‐subject variances. A basic assumption is that the two variance terms are constant across subjects. However, we often find that these variances are functions of covariates, and the data set has what we call explainable heterogeneity, which needs to be allowed for in the model. We characterize several types of heterogeneity of variance in REMs and develop three diagnostic tests using the score statistic: one for each of the two variance terms, and the third for a form of multivariate nonconstant variance. For each test we present an adjusted residual plot which can identify cases that are unusually influential on the outcome of the test.     

Bayesian Robustness Modelling of Location and Scale Parameters

Abstract. The modelling process in Bayesian Statistics constitutes the fundamental stage of the analysis, since depending on the chosen probability laws the inferences may vary considerably. This is particularly true when conflicts arise between two or more sources of information. For instance, inference in the presence of an outlier (which conflicts with the information provided by the other observations) can be highly dependent on the assumed sampling distribution. When heavy‐tailed (e.g. t) distributions are used, outliers may be rejected whereas this kind of robust inference is not available when we use light‐tailed (e.g. normal) distributions. A long literature has established sufficient conditions on location‐parameter models to resolve conflict in various ways. In this work, we consider a location–scale parameter structure, which is more complex than the single parameter cases because conflicts can arise between three sources of information, namely the likelihood, the prior distribution for the location parameter and the prior for the scale parameter. We establish sufficient conditions on the distributions in a location–scale model to resolve conflicts in different ways as a single observation tends to infinity. In addition, for each case, we explicitly give the limiting posterior distributions as the conflict becomes more extreme.     

Goodness-of-fit testing under long memory

This paper discusses the problem of fitting a distribution function to the marginal distribution of a long memory moving average process. Because of the uniform reduction principle, unlike in the i.i.d. set up, classical tests based on empirical process are relatively easy to implement. More importantly, we discuss fitting the marginal distribution of the error process in location, scale, location–scale and linear regression models. An interesting observation is that in the location model, location–scale model, or more generally in the linear regression models with non-zero intercept parameter, the null weak limit of the first order difference between the residual empirical process and the null model is degenerate at zero, and hence it cannot be used to fit an error distribution in these models for the large samples. This finding is in sharp contrast to a recent claim of Chan and Ling (2008) that the null weak limit of such a process is a continuous Gaussian process. This note also proposes some tests based on the second order difference for the location case. Another finding is that residual empirical process tests in the scale problem are robust against not knowing the scale parameter.     

Gaussian mixture analysis of covariance

ABSTRACT

In many real-world applications, the traditional theory of analysis of covariance (ANCOVA) leads to inadequate and unreliable results because of violation of the response variable observations from the essential Gaussian assumption that may be due to the heterogeneity of population, the presence of outlier or both of them. In this paper, we develop a Gaussian mixture ANCOVA model for modelling heterogeneous populations with a finite number of subpopulation. We provide the maximum likelihood estimates of the model parameters via an EM algorithm. We also drive the adjusted effects estimators for treatments and covariates. The Fisher information matrix of the model and asymptotic confidence intervals for the parameter are also discussed. We performed a simulation study to assess the performance of the proposed model. A real-world example is also worked out to explained the methodology.