Journal of Probability and Statistics

Volume 2017 (2017), Article ID 6579537, 24 pages

https://doi.org/10.1155/2017/6579537

## Numerical Reconstruction of the Covariance Matrix of a Spherically Truncated Multinormal Distribution

^{1}Istituto Nazionale di Statistica (ISTAT), Via Cesare Balbo 16, 00184 Rome, Italy^{2}Italian Agency for New Technologies, Energy and Sustainable Economic Development (ENEA), Via Enrico Fermi 45, 00044 Frascati, Italy

Correspondence should be addressed to Filippo Palombi; ti.aene@ibmolap.oppilif

Received 7 July 2016; Revised 7 October 2016; Accepted 12 October 2016; Published 10 January 2017

Academic Editor: Ramón M. Rodríguez-Dagnino

Copyright © 2017 Filippo Palombi et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

We relate the matrix of the second moments of a spherically truncated normal multivariate to its full covariance matrix and present an algorithm to invert the relation and reconstruct from . While the eigenvectors of are left invariant by the truncation, its eigenvalues are nonuniformly damped. We show that the eigenvalues of can be reconstructed from their truncated counterparts via a fixed point iteration, whose convergence we prove analytically. The procedure requires the computation of multidimensional Gaussian integrals over an Euclidean ball, for which we extend a numerical technique, originally proposed by Ruben in 1962, based on a series expansion in chi-square distributions. In order to study the feasibility of our approach, we examine the convergence rate of some iterative schemes on suitably chosen ensembles of Wishart matrices. We finally discuss the practical difficulties arising in sample space and outline a regularization of the problem based on perturbation theory.

#### 1. Introduction

It has been more than forty years since Tallis [1] worked out the moment-generating function of a normal multivariate , subject to the conditional eventThe perfect match between the symmetries of the ellipsoid and those of allows for an exact analytic result, from which the complete set of multivariate truncated moments can be extracted upon differentiation. Consider, for instance, the matrix of the second truncated moments, expressing the covariances among the components of within . From Tallis’ paper it turns out thatwith denoting the cumulative distribution function of a -variable with degrees of freedom. Inverting (2)—so as to express as a function of —is trivial, since is a scalar damping factor independent of . In this paper, we shall refer to such inverse relation as the* reconstruction* of from . Unfortunately, life is not always so easy. In general, the effects produced on the expectation of functions of by cutting off the probability density outside a generic domain can be hardly calculated in closed form, especially if the boundary of is shaped in a way that breaks the ellipsoidal symmetry of . Thus, for instance, unlike (2), the matrix of the second truncated moments is expected in general to display a nonlinear/nontrivial dependence upon .

In the present paper, we consider the case where is a -dimensional Euclidean ball with center in the origin and square radius . Specifically, we discuss the reconstruction of from the matrix of the spherically truncated second moments. To this aim, we need to mimic Tallis’ calculation, with (1) being replaced by the conditional eventThis is precisely an example of the situation described in the previous paragraph: although has a higher degree of symmetry than , still there is no possibility of obtaining a closed-form relation between and , since breaks the ellipsoidal symmetry of : technically speaking, in this case we cannot perform any change of variable under the defining integral of the moment-generating function, which may reduce the dimensionality of the problem, as in Tallis’ paper.

In spite of that, the residual symmetries characterizing the truncated distribution help simplify the problem in the following respects: (i) the reflection invariance of the whole setup still yields and (ii) the rotational invariance of preserves the possibility of defining the principal components of the distribution just like in the unconstrained case. In particular, the latter property means that and share the same orthonormal eigenvectors. In fact, the reconstruction of from amounts to solving a system of nonlinear integral equations, having the eigenvalues of as unknown variables and the eigenvalues of as input parameters. In a lack of analytic techniques to evaluate exactly the integrals involved, we resort to a numerical algorithm, of which we investigate feasibility, performance, and optimization.

The paper is organized as follows. In Section 2, we describe a few examples illustrating the occurrence of spherical truncations in practical situations. In Section 3, we show that the aforementioned integral equations have the analytic structure of a fixed point vector equation; that is to say, . This suggests achieving the reconstruction of numerically via suitably chosen iterative schemes. In Section 4, we prove the convergence of the simplest of them by inductive arguments, the validity of which relies upon the monotonicity properties of ratios of Gaussian integrals over . In Section 5, we review some numerical techniques for the computation of Gaussian integrals over with controlled systematic error. These are based on and extend a classic work by Ruben [2] on the distribution of quadratic forms of normal variables. For the sake of readability, we defer proofs of statements made in this section to the appendix. In Section 6, we report on our numerical experiences: since the simplest iterative scheme, namely, the Gauss–Jacobi iteration, is too slow for practical purposes, we investigate the performance of its improved version based on overrelaxation; as expected, we find that the latter has a higher convergence rate; yet it still slows down polynomially in as and exponentially in as ; in order to reduce the slowing down, we propose an acceleration technique, which boosts the higher components of the eigenvalue spectrum. A series of Monte Carlo simulations enables us to quantify the speedup. In Section 7 we discuss the problems arising when is affected by statistical uncertainty and propose a regularization technique based on perturbation theory. To conclude, we summarize our findings in Section 8.

#### 2. Motivating Examples

Spherical truncations of multinormal distributions may characterize different kinds of experimental devices and may occur in various problems of statistical and convex analysis. In this section, we discuss two motivating examples.

##### 2.1. A Two-Dimensional* Gedanken* Experiment in Classical Particle Physics

Consider the following ideal situation. An accelerator physicist prepares an elliptical beam of classical particles with Gaussian transversal profile. The experimenter knows* a priori* the spatial distribution of the beam, that is, the covariance matrix of the two-dimensional coordinates of the particles on a plane orthogonal to their flight direction. We can assume with no loss of generality that the transversal coordinate system has origin at the maximum of the beam intensity and axes along the principal components of the beam; thus it holds . The beam travels straightforward until it enters a linear coaxial pipeline with circular profile, schematically depicted in Figure 1, where the beam is longitudinally accelerated. While the outer part of the beam is stopped by an absorbing wall, the inner part propagates within the pipeline. At the end of the beam flight the physicist wants to know if the transversal distribution of the particles is changed, due to unknown disturbance factors arisen within the pipeline. Accordingly, he measures again the spatial covariance matrix of the beam. Unfortunately, the absorbing wall has cut off the Gaussian tail, thus damping the covariance matrix and making it no more comparable to the original one. To perform such a comparison in the general case , the covariance matrix of the truncated circular beam has to go through the reconstruction procedure described in next sections.