Computes a robust estimate of a covariance matrix and mean vector.

The type double function is imsls_d_robust_covariances.

Matrix of size n_variables by n_variables containing the matrix of covariances.

Function imsls_f_robust_covariances computes robust M-estimates of the mean and covariance matrix from a matrix of observations. A pooled estimate of the covariance matrix is computed when multiple groups are present in the input data. M-estimate weights are obtained using the “minimax” weights of Huber (1981, pp. 231-235), with percentage expected gross errors. Huber’s (1981) weighting equations are given by:

User specified observation weights and frequencies may be given for each row in x. Listwise deletion of missing values is assumed so that all observations used are “complete”.

Let f (x;μi, Σ) denote the density of an observation p-vector x in population (group) i with mean vector μi, for i = 1, …, τ. Let the covariance matrix Σ be such that Σ = RTR. If

then

It is assumed that g(y) is a spherically symmetric density in p-dimensions.

In imsls_f_robust_covariances, Σ and μi are estimated as the solutions

of the estimation equations

and

where i indexes the τ groups, ni, is the number of observations in group i, fij is the frequency for the j‑th observation in group i, wij is the observation weight specified in column iwt of x, Ip is a p ×  p identity matrix,

w(r) and u(r) are the weighting functions, and where β is a constant computed by the program to make the expected weighted Mahalanobis distance (yTy) equal the expected Mahalanobis distance from a multivariate normal distribution (see Marazzi 1985). The constant β is described more fully below.

Function imsls_f_robust_covariances uses one of two algorithms for solving the estimation equations. The first algorithm is discussed in detail in Huber (1981) and is a variant of the conjugate gradient method. The second algorithm is due to Stahel (1981) and is discussed in detail by Marazzi (1985). In both algorithms, correction vectors Tki for the group i means and correction matrix Wk = Ip + Uk for the Cholesky factorization of Σ are found such that the updated mean vectors are given by

and the updated matrix R is given as

where k is the iteration number and

When all elements of Uk and Tki are less than ɛ = tolerance, convergence is assumed.

Three methods for obtaining estimates are allowed. In the first method, the sample weighted estimate of Σ is computed. In the second method, estimates based upon the median and the interquartile range are used. Finally, in the last method, the user inputs initial estimates.

Function imsls_f_robust_covariances computes estimates based on the “minimax” weights discussed above. The constant β is chosen such that E (u(r)r2) = ρβ where the expectation is with respect to a standard p‑variate multivariate normal distribution. This yields estimates with the correct expectation for the multivariate normal distribution (for given mean vector). The expectation is computed via integration of estimated spline function. 200 knots are used on an equally spaced grid from 0.0 to the 99.999 percentile of

distribution. An error estimate is computed based upon 100 of these knots. If the estimated relative error is greater than 0.0001, a warning message is issued. If β is not computed accurately (i.e., if the warning message is issued), the computed estimates are still optimal, but the scale of the estimated covariance matrix may need to be multiplied by a constant in order for

to have the correct multivariate normal covariance expectation.

The following example computes a robust variance-covariance matrix. The last column of the data set is the group indicator.

The following example computes estimates of the pooled covariance matrix for the Fisher’s iris data. For comparison, the estimates are first computed via function imsls_f_pooled_covariances. Function imsls_f_robust_covariances with percentage = 2.0 is then used to compute the robust estimates. As can be seen from the output, the resulting estimates are quite similar.

Next, three observations are made into outliers, and again, estimates are computed using functions imsls_f_pooled_covariances and imsls_f_robust_covariances. When outliers are present, the estimates of imsls_f_pooled_covariances are adversely affected, while the estimates produced by imsls_f_robust_covariances are close the estimates produced when no outliers are present.