hide long namesshow long names
hide short namesshow short names
Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

NAG Toolbox: nag_correg_robustm_corr_huber (g02hk)

Purpose

nag_correg_robustm_corr_huber (g02hk) computes a robust estimate of the covariance matrix for an expected fraction of gross errors.

Syntax

[cov, theta, nit, ifail] = g02hk(x, eps, 'n', n, 'm', m, 'maxit', maxit, 'nitmon', nitmon, 'tol', tol)
[cov, theta, nit, ifail] = nag_correg_robustm_corr_huber(x, eps, 'n', n, 'm', m, 'maxit', maxit, 'nitmon', nitmon, 'tol', tol)
Note: the interface to this routine has changed since earlier releases of the toolbox:
Mark 22: n has been made optional
Mark 23: nitmon, tol now optional
.

Description

For a set of nn observations on mm variables in a matrix XX, a robust estimate of the covariance matrix, CC, and a robust estimate of location, θθ, are given by
C = τ2(ATA)1,
C=τ2(ATA)-1,
where τ2τ2 is a correction factor and AA is a lower triangular matrix found as the solution to the following equations:
zi = A(xiθ),
zi=A(xi-θ),
n
1/nw(zi2)zi = 0,
i = 1
1n i= 1nw(zi2)zi=0,
and
n
1/nu(zi2)ziziTI = 0,
i = 1
1ni=1nu(zi2)zi ziT -I=0,
where xixi is a vector of length mm containing the elements of the iith row of x,
zizi is a vector of length mm,
II is the identity matrix and 00 is the zero matrix,
and ww and uu are suitable functions.
nag_correg_robustm_corr_huber (g02hk) uses weight functions:
u(t) = (au)/(t2), if ​t < au2
u(t) = 1, if ​au2tbu2
u(t) = (bu)/(t2), if ​t > bu2
u(t)= aut2, if ​t<au2 u(t)=1, if ​au2tbu2 u(t)= but2, if ​t>bu2
and
w(t) = 1, if ​tcw
w(t) = (cw)/t, if ​t > cw
w(t)= 1, if ​tcw w(t)= cwt, if ​t>cw
for constants auau, bubu and cwcw.
These functions solve a minimax problem considered by Huber (see Huber (1981)). The values of auau, bubu and cwcw are calculated from the expected fraction of gross errors, εε (see Huber (1981) and Marazzi (1987)). The expected fraction of gross errors is the estimated proportion of outliers in the sample.
In order to make the estimate asymptotically unbiased under a Normal model a correction factor, τ2τ2, is calculated, (see Huber (1981) and Marazzi (1987)).
The matrix CC is calculated using nag_correg_robustm_corr_user_deriv (g02hl). Initial estimates of θjθj, for j = 1,2,,mj=1,2,,m, are given by the median of the jjth column of XX and the initial value of AA is based on the median absolute deviation (see Marazzi (1987)). nag_correg_robustm_corr_huber (g02hk) is based on routines in ROBETH; see Marazzi (1987).

References

Huber P J (1981) Robust Statistics Wiley
Marazzi A (1987) Weights for bounded influence regression in ROBETH Cah. Rech. Doc. IUMSP, No. 3 ROB 3 Institut Universitaire de Médecine Sociale et Préventive, Lausanne

Parameters

Compulsory Input Parameters

1:     x(ldx,m) – double array
ldx, the first dimension of the array, must satisfy the constraint ldxnldxn.
x(i,j)xij must contain the iith observation for the jjth variable, for i = 1,2,,ni=1,2,,n and j = 1,2,,mj=1,2,,m.
2:     eps – double scalar
εε, the expected fraction of gross errors expected in the sample.
Constraint: 0.0eps < 1.00.0eps<1.0.

Optional Input Parameters

1:     n – int64int32nag_int scalar
Default: The first dimension of the array x.
nn, the number of observations.
Constraint: n > 1n>1.
2:     m – int64int32nag_int scalar
Default: The second dimension of the array x.
mm, the number of columns of the matrix XX, i.e., number of independent variables.
Constraint: 1mn1mn.
3:     maxit – int64int32nag_int scalar
The maximum number of iterations that will be used during the calculation of the covariance matrix.
Default: 150150
Constraint: maxit > 0maxit>0.
4:     nitmon – int64int32nag_int scalar
Indicates the amount of information on the iteration that is printed.
nitmon > 0nitmon>0
The value of AA, θθ and δδ (see Section [Accuracy]) will be printed at the first and every nitmon iterations.
nitmon0nitmon0
No iteration monitoring is printed.
When printing occurs the output is directed to the current advisory message unit (see nag_file_set_unit_advisory (x04ab)).
Default: 00
5:     tol – double scalar
The relative precision for the final estimates of the covariance matrix.
Default: 5e-55e-5
Constraint: tol > 0.0tol>0.0.

Input Parameters Omitted from the MATLAB Interface

ldx wk

Output Parameters

1:     cov(m × (m + 1) / 2m×(m+1)/2) – double array
A robust estimate of the covariance matrix, CC. The upper triangular part of the matrix CC is stored packed by columns. CijCij is returned in cov((j × (j1) / 2 + i))cov(j×(j-1)/2+i), ijij.
2:     theta(m) – double array
The robust estimate of the location parameters θjθj, for j = 1,2,,mj=1,2,,m.
3:     nit – int64int32nag_int scalar
The number of iterations performed.
4:     ifail – int64int32nag_int scalar
ifail = 0ifail=0 unless the function detects an error (see [Error Indicators and Warnings]).

Error Indicators and Warnings

Errors or warnings detected by the function:
  ifail = 1ifail=1
On entry,n1n1,
orm < 1m<1,
orn < mn<m,
orldx < nldx<n,
oreps < 0.0eps<0.0,
oreps1.0eps1.0,
ortol0.0tol0.0,
ormaxit0maxit0.
  ifail = 2ifail=2
On entry,a variable has a constant value, i.e., all elements in a column of XX are identical.
  ifail = 3ifail=3
The iterative procedure to find CC has failed to converge in maxit iterations.
  ifail = 4ifail=4
The iterative procedure to find CC has become unstable. This may happen if the value of eps is too large for the sample.

Accuracy

On successful exit the accuracy of the results is related to the value of tol; see Section [Parameters]. At an iteration let
(i) d1 = d1= the maximum value of the absolute relative change in AA
(ii) d2 = d2= the maximum absolute change in u(zi2)u(zi2)
(iii) d3 = d3= the maximum absolute relative change in θjθj
and let δ = max (d1,d2,d3)δ=max(d1,d2,d3). Then the iterative procedure is assumed to have converged when δ < tolδ<tol

Further Comments

The existence of AA, and hence CC, will depend upon the function uu (see Marazzi (1987)); also if XX is not of full rank a value of AA will not be found. If the columns of XX are almost linearly related, then convergence will be slow.

Example

function nag_correg_robustm_corr_huber_example
x = [3.4, 6.9, 12.2;
     6.4, 2.5, 15.1;
     4.9, 5.5, 14.2;
     7.3, 1.9, 18.2;
     8.8, 3.6, 11.7;
     8.4, 1.3, 17.9;
     5.3, 3.1, 15;
     2.7, 8.1, 7.7;
     6.1, 3, 21.9;
     5.3, 2.2, 13.9];
epsilon = 0.1;
[covar, theta, nit, ifail] = nag_correg_robustm_corr_huber(x, epsilon)
 

covar =

    3.4611
   -3.6806
    5.3477
    4.6818
   -6.6445
   14.4389


theta =

    5.8178
    3.6813
   15.0369


nit =

                   23


ifail =

                    0


function g02hk_example
x = [3.4, 6.9, 12.2;
     6.4, 2.5, 15.1;
     4.9, 5.5, 14.2;
     7.3, 1.9, 18.2;
     8.8, 3.6, 11.7;
     8.4, 1.3, 17.9;
     5.3, 3.1, 15;
     2.7, 8.1, 7.7;
     6.1, 3, 21.9;
     5.3, 2.2, 13.9];
epsilon = 0.1;
[covar, theta, nit, ifail] = g02hk(x, epsilon)
 

covar =

    3.4611
   -3.6806
    5.3477
    4.6818
   -6.6445
   14.4389


theta =

    5.8178
    3.6813
   15.0369


nit =

                   23


ifail =

                    0



PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2013