Marc G Genton

The supplement contains additional analyses and plots in support to the main findings in the paper. The code for this work is available at the following GitHub repository: github.com/Env-an-Stat-group/24.Zhang.AoAS.

The unified skew-t (SUT) is a flexible parametric multivariate distribution that accounts for skewness and heavy tails in the data. A few of its properties can be found scattered in the literature or in a parameterization that does not follow the original one for unified skew-normal (SUN) distributions, yet a systematic study is lacking. In this work, explicit properties of the multivariate SUT distribution are presented, such as its stochastic representations, moments, SUN-scale mixture representation, linear transformation, additivity, marginal distribution, canonical form, quadratic form, conditional distribution, change of latent dimensions, Mardia measures of multivariate skewness and kurtosis, and non-identifiability issue. These results are given in a parametrization that reduces to the original SUN distribution as a sub-model, hence facilitating the use of the SUT for applications. Several models based on the SUT distribution are provided for illustration.

The broad class of multivariate unified skew-normal (SUN) distributions has been recently shown to possess fundamental conjugacy properties. When used as priors for the vector of parameters in general probit, tobit, and multinomial probit models, these distributions yield posteriors that still belong to the SUN family. Although such a core result has led to important advancements in Bayesian inference and computation, its applicability beyond likelihoods associated with fully-observed, discretized, or censored realizations from multivariate Gaussian models remains yet unexplored. This article covers such an important gap by proving that the wider family of multivariate unified skew-elliptical (SUE) distributions, which extends SUNs to more general perturbations of elliptical densities, guarantees conjugacy for broader classes of models, beyond those relying on fully-observed, discretized or censored Gaussians. Such a result leverages the closure under linear combinations, conditioning and marginalization of SUE to prove that such a family is conjugate to the likelihood induced by general multivariate regression models for fully-observed, censored or dichotomized realizations from skew-elliptical distributions. This advancement substantially enlarges the set of models that enable conjugate Bayesian inference to general formulations arising from elliptical and skew-elliptical families, including the multivariate Student's t and skew-t, among others.

We introduce a multivariate version of the modified skew-normal distribution, which contains the multivariate normal distribution as a special case. Unlike the Azzalini multivariate skew-normal distribution, this new distribution has a nonsingular Fisher information matrix when the skewness parameters are all zero, and its profile log-likelihood of the skewness parameters is always a non-monotonic function. We study some basic properties of the proposed family of distributions and present an expectation-maximization (EM) algorithm for parameter estimation that we validate through simulation studies. Finally, we apply the proposed model to the univariate frontier data and to a trivariate wind speed data, and compare its performance with the Azzalini skew-normal model.

Spatial statistical modeling and prediction involve generating and manipulating an n*n symmetric positive definite covariance matrix, where n denotes the number of spatial locations. However, when n is large, processing this covariance matrix using traditional methods becomes prohibitive. Thus, coupling parallel processing with approximation can be an elegant solution to this challenge by relying on parallel solvers that deal with the matrix as a set of small tiles instead of the full structure. Each processing unit can process a single tile, allowing better performance. The approximation can also be performed at the tile level for better compression and faster execution. The Tile Low-Rank (TLR) approximation, a tile-based approximation algorithm, has recently been used in spatial statistics applications. However, the quality of TLR algorithms mainly relies on ordering the matrix elements. This order can impact the compression quality and, therefore, the efficiency of the underlying linear solvers, which highly depends on the individual ranks of each tile. Thus, herein, we aim to investigate the accuracy and performance of some existing ordering algorithms that are used to order the geospatial locations before generating the spatial covariance matrix. Furthermore, we highlight the pros and cons of each ordering algorithm in the context of spatial statistics applications and give hints to practitioners on how to choose the ordering algorithm carefully. We assess the quality of the compression and the accuracy of the statistical parameter estimates of the Mat\'ern covariance function using TLR approximation under various ordering algorithms and settings of …

Asymptotic methods for hypothesis testing in high-dimensional data usually require the dimension of the observations to increase to infinity, often with an additional condition on its rate of increase compared to the sample size. On the other hand, multivariate asymptotic methods are valid for fixed dimension only, and their practical implementations in hypothesis testing methodology typically require the sample size to be large compared to the dimension for yielding desirable results. However, in practical scenarios, it is usually not possible to determine whether the dimension of the data at hand conform to the conditions required for the validity of the high-dimensional asymptotic methods, or whether the sample size is large enough compared to the dimension of the data. In this work, a theory of asymptotic convergence is proposed, which holds uniformly over the dimension of the random vectors. This theory attempts to unify the asymptotic results for fixed-dimensional multivariate data and high-dimensional data, and accounts for the effect of the dimension of the data on the performance of the hypothesis testing procedures. The methodology developed based on this asymptotic theory can be applied to data of any dimension. An application of this theory is demonstrated in the two-sample test for the equality of locations. The test statistic proposed is unscaled by the sample covariance, similar to usual tests for high-dimensional data. Using simulated examples, it is demonstrated that the proposed test exhibits better performance compared to several popular tests in the literature for high-dimensional data. Further, it is demonstrated in simulated models …

We thank Florian Stijven, Ariel Alonso Abad, and Gökçe Deliorman for pointing out typos in the signs of the formula for the mutual information index for the Student's t case (p. 47, last formula of Section 2.4). The correct formula is:I XY T n+ m (Ω, ν)= I XY N n+ m (Ω)+ log Γ (ν/2) Γ {(ν+ n+ m)/2} Γ {(ν+ n)/2} Γ {(ν+ m)/2}+ ν+ m 2 ψ ν+ m 2+ ν+ n 2 ψ ν+ n 2− ν+ n+ m 2 ψ ν+ n+ m 2− ν 2 ψ ν 2. $${\displaystyle\begin {array}{ll}\hfill {I} _ {\mathbf {XY}}^{T_ {n+ m}}\left (\boldsymbol {\Omega},\nu\right) &={I} _ {\mathbf {XY}}^{N_ {n+ m}}\left (\boldsymbol {\Omega}\right)+\log\left [\frac {\Gamma\left (\nu/2\right)\Gamma\left\{\left (\nu+ n+ m\right)/2\right\}}{\Gamma\left\{\left (\nu+ n\right)/2\right\}\Gamma\left\{\left (\nu+ m\right)/2\right\}}\right]+\frac {\nu+ m}{2}\psi\left (\frac {\nu+ m}{2}\right)\\{}\hfill &\kern1em+\frac {\nu+ n}{2}\psi\left (\frac {\nu+ n}{2}\right)-\frac {\nu+ n+ m}{2}\psi\left (\frac {\nu+ n+ m}{2}\right)-\frac {\nu}{2}\psi\left (\frac {\nu}{2 …

Gaussian processes (GPs) are commonly used for geospatial analysis, but they suffer from high computational complexity when dealing with massive data. For instance, the log-likelihood function required in estimating the statistical model parameters for geospatial data is a computationally intensive procedure that involves computing the inverse of a covariance matrix with size n X n, where n represents the number of geographical locations. As a result, in the literature, studies have shifted towards approximation methods to handle larger values of n effectively while maintaining high accuracy. These methods encompass a range of techniques, including low-rank and sparse approximations. Vecchia approximation is one of the most promising methods to speed up evaluating the log-likelihood function. This study presents a parallel implementation of the Vecchia approximation, utilizing batched matrix computations on contemporary GPUs. The proposed implementation relies on batched linear algebra routines to efficiently execute individual conditional distributions in the Vecchia algorithm. We rely on the KBLAS linear algebra library to perform batched linear algebra operations, reducing the time to solution compared to the state-of-the-art parallel implementation of the likelihood estimation operation in the ExaGeoStat software by up to 700X, 833X, 1380X on 32GB GV100, 80GB A100, and 80GB H100 GPUs, respectively. We also successfully manage larger problem sizes on a single NVIDIA GPU, accommodating up to 1M locations with 80GB A100 and H100 GPUs while maintaining the necessary application accuracy. We further assess the …

We introduce a new family of multivariate distributions by taking the component-wise Tukey-h transformation of a random vector following a skew-normal distribution with an alternative parameterization. The proposed distribution is named the skew-normal-Tukey-h distribution and is an extension of the skew-normal distribution for handling heavy-tailed data. We compare this proposed distribution to the skew-t distribution, which is another extension of the skew-normal distribution for modeling tail-thickness, and demonstrate that when there are substantial differences in marginal kurtosis, the proposed distribution is more appropriate. Moreover, we derive many appealing stochastic properties of the proposed distribution and provide a methodology for the estimation of the parameters that can be applied to large dimensions. Using simulations, as well as a wine and a wind speed data application, we illustrate how to …

The enclosed datasets have been generated by the internal spatial data generator tool included in the ExaGeoStat software (https://github. com/ecrc/exageostat). It was used for the 2023 KAUST Competition on Spatial Statistics for Large Datasets (https://cemse. kaust. edu. sa/stsds/news/2023-kaust-competition-spatial-statistics-large-datasets). The competition had four parts (Sub-competition 1a, Sub-competition 1b, Sub-competition 2a, and Sub-competition 2b). The main purpose of the competition was to reassess existing approximation methods on large spatial datasets in a uniform way that guarantees a fair comparison. More information about the datasets can be found in the file" Competition Description and True Model".

There is a wide availability of methods for testing normality under the assumption of independent and identically distributed data. When data are dependent in space and/or time, however, assessing and testing the marginal behavior is considerably more challenging, as the marginal behavior is impacted by the degree of dependence. We propose a new approach to assess normality for dependent data by non-linearly incorporating existing statistics from normality tests as well as sample moments such as skewness and kurtosis through a neural network. We calibrate (deep) neural networks by simulated normal and non-normal data with a wide range of dependence structures and we determine the probability of rejecting the null hypothesis. We compare several approaches for normality tests and demonstrate the superiority of our method in terms of statistical power through an extensive simulation study. A real world application to global temperature data further demonstrates how the degree of spatio-temporal aggregation affects the marginal normality in the data.

Large-scale parallel computing is crucial in Gaussian regressions to reduce the time complexity of spatial statistics applications. The log-likelihood function is utilized to evaluate the Gaussian model for a set of measurements in ???? geographical locations. Several studies have shown an utilization of modern hardware to scale the log-likelihood function for handling large numbers of locations. ExaGeoStat is an example of software that allows parallel statistical parameter estimation from the log-likelihood function. However, generating a covariance matrix is mandatory and challenging when estimating the log-likelihood function. In ExaGeoStat, the generation process was performed on CPU hardware due to missing math functions in CUDA libraries, eg, the modified Bessel function of the second kind. This study aims to optimize the generation process using GPU with two proposed generation schemes: pure GPU and hybrid. Our implementations demonstrate up to 6X speedup with pure GPU and up to 1.5 X speedup with the hybrid scheme.

HPC-based applications often have complex workflows with many software dependencies that hinder their portability on contemporary HPC architectures. In addition, these applications often require extraordinary efforts to deploy and execute at performance potential on new HPC systems, while the users expert in these applications generally have less expertise in HPC and related technologies. This paper provides a dynamic solution that facilitates containerization for transferring HPC software onto diverse parallel systems. The study relies on the HPC Workflow as a Service (HPCWaaS) paradigm proposed by the EuroHPC eFlows4HPC project. It offers to deploy workflows through containers tailored for any of a number of specific HPC systems. Traditional container image creation tools rely on OS system packages compiled for generic architecture families (x86\_64, amd64, ppc64, ...) and specific MPI or GPU runtime library versions. The containerization solution proposed in this paper leverages HPC Builders such as Spack or Easybuild and multi-platform builders such as buildx to create a service for automating the creation of container images for the software specific to each hardware architecture, aiming to sustain the overall performance of the software. We assess the efficiency of our proposed solution for porting the geostatistics ExaGeoStat software on various parallel systems while preserving the computational performance. The results show that the performance of the generated images is comparable with the native execution of the software on the same architectures. On the distributed-memory system, the containerized version can …

The assumption of normality has underlain much of the development of statistics, including spatial statistics, and many tests have been proposed. In this work, we focus on the multivariate setting and first review the recent advances in multivariate normality tests for i.i.d. data, with emphasis on the skewness and kurtosis approaches. We show through simulation studies that some of these tests cannot be used directly for testing normality of spatial data. We further review briefly the few existing univariate tests under dependence (time or space), and then propose a new multivariate normality test for spatial data by accounting for the spatial dependence. The new test utilises the union‐intersection principle to decompose the null hypothesis into intersections of univariate normality hypotheses for projection data, and it rejects the multivariate normality if any individual hypothesis is rejected. The individual hypotheses for …

The prevalence of multivariate space-time data collected from monitoring networks and satellites, or generated from numerical models, has brought much attention to multivariate spatio-temporal statistical models, where the covariance function plays a key role in modeling, inference, and prediction. For multivariate space-time data, understanding the spatio-temporal variability, within and across variables, is essential in employing a realistic covariance model. Meanwhile, the complexity of generic covariances often makes model fitting very challenging, and simplified covariance structures, including symmetry and separability, can reduce the model complexity and facilitate the inference procedure. However, a careful examination of these properties is needed in real applications. In the work presented here, we formally define these properties for multivariate spatio-temporal random fields and use functional data …

Correlated binary response data with covariates are ubiquitous in longitudinal or spatial studies. Among the existing statistical models, the most well‐known one for this type of data is the multivariate probit model, which uses a Gaussian link to model dependence at the latent level. However, a symmetric link may not be appropriate if the data are highly imbalanced. Here, we propose a multivariate skew‐elliptical link model for correlated binary responses, which includes the multivariate probit model as a special case. Furthermore, we perform Bayesian inference for this new model and prove that the regression coefficients have a closed‐form unified skew‐elliptical posterior with an elliptical prior. The new methodology is illustrated by an application to COVID‐19 data from three different counties of the state of California, USA. By jointly modeling extreme spikes in weekly new cases, our results show that the spatial …

To establish a more sustainable future, Saudi Arabia is striving to reduce its dependency on fossil fuels and promote renewable energies. Solar energy is a major resource in Saudi Arabia because of the country's geographical location with year‐round clear skies and plentiful sunlight. However, although solar energy is a clean and safe renewable energy resource, it can be quite unpredictable. Therefore, solar irradiance needs to be forecasted and simulated as accurately as possible on both regional and national scales in the Kingdom to assist in planning for reserve usage, switching sources and short‐term power purchases. Based on an hourly solar diffuse horizontal irradiance (DHI) dataset from 45 different Saudi Arabian solar monitoring stations, this study proposes a novel spatio‐temporal model. We identify the temporal dependency of DHI as cyclostationary and incorporate this key observation in the model …

Given the computational challenges involved in calculating the maximum likelihood estimates for large spatial datasets, there has been significant interest in the research community regarding approximation methods for estimation and subsequent predictions. However, prior studies examining the evaluation of these methods have primarily focused on scenarios where the data are observed on a regular grid or originate from a uniform distribution of locations. Nevertheless, non-uniformly distributed locations are commonplace in fields like meteorology and ecology. Examples include gridded data with missing observations acquired through remote sensing techniques. To assess the reliability and effectiveness of cutting-edge approximation methods, we have initiated a competition focused on estimation and prediction for large spatial datasets with non-uniformly distributed locations. Participants were invited to …

A novel elastic time distance for sparse multivariate functional data is proposed and used to develop a robust distance-based two-layer partition clustering method. With this proposed distance, the new approach not only can detect correct clusters for sparse multivariate functional data under outlier settings but also can detect those outliers that do not belong to any clusters. Classical distance-based clustering methods such as density-based spatial clustering of applications with noise (DBSCAN), agglomerative hierarchical clustering, and K-medoids are extended to the sparse multivariate functional case based on the newly-proposed distance. Numerical experiments on simulated data highlight that the performance of the proposed algorithm is superior to the performances of existing model-based and extended distance-based methods. The effectiveness of the proposed approach is demonstrated using Northwest …

Spatial processes observed in various fields, such as climate and environmental science, often occur on a large scale and demonstrate spatial nonstationarity. Fitting a Gaussian process with a nonstationary Mat\'ern covariance is challenging. Previous studies in the literature have tackled this challenge by employing spatial partitioning techniques to estimate the parameters that vary spatially in the covariance function. The selection of partitions is an important consideration, but it is often subjective and lacks a data-driven approach. To address this issue, in this study, we utilize the power of Convolutional Neural Networks (ConvNets) to derive subregions from the nonstationary data. We employ a selection mechanism to identify subregions that exhibit similar behavior to stationary fields. In order to distinguish between stationary and nonstationary random fields, we conducted training on ConvNet using various simulated data. These simulations are generated from Gaussian processes with Mat\'ern covariance models under a wide range of parameter settings, ensuring adequate representation of both stationary and nonstationary spatial data. We assess the performance of the proposed method with synthetic and real datasets at a large scale. The results revealed enhanced accuracy in parameter estimations when relying on ConvNet-based partition compared to traditional user-defined approaches.