nag_robust_m_corr_user_fn_no_derr (g02hmc) (PDF version)
g02 Chapter Contents
g02 Chapter Introduction
NAG C Library Manual

NAG Library Function Document

nag_robust_m_corr_user_fn_no_derr (g02hmc)

+ Contents

    1  Purpose
    7  Accuracy

1  Purpose

nag_robust_m_corr_user_fn_no_derr (g02hmc) computes a robust estimate of the covariance matrix for user-supplied weight functions. The derivatives of the weight functions are not required.

2  Specification

#include <nag.h>
#include <nagg02.h>
void  nag_robust_m_corr_user_fn_no_derr (Nag_OrderType order,
void (*ucv)(double t, double *u, double *w, Nag_Comm *comm),
Integer indm, Integer n, Integer m, const double x[], Integer pdx, double cov[], double a[], double wt[], double theta[], double bl, double bd, Integer maxit, Integer nitmon, const char *outfile, double tol, Integer *nit, Nag_Comm *comm, NagError *fail)

3  Description

For a set of n observations on m variables in a matrix X, a robust estimate of the covariance matrix, C, and a robust estimate of location, θ, are given by
C=τ2ATA-1,
where τ2 is a correction factor and A is a lower triangular matrix found as the solution to the following equations.
zi=Axi-θ
1n i= 1nwzi2zi=0
and
1ni=1nuzi2zi ziT -vzi2I=0,
where xi is a vector of length m containing the elements of the ith row of X,
zi is a vector of length m,
I is the identity matrix and 0 is the zero matrix.
and w and u are suitable functions.
nag_robust_m_corr_user_fn_no_derr (g02hmc) covers two situations:
(i) vt=1 for all t,
(ii) vt=ut.
The robust covariance matrix may be calculated from a weighted sum of squares and cross-products matrix about θ using weights wti=uzi. In case (i) a divisor of n is used and in case (ii) a divisor of i=1nwti is used. If w.=u., then the robust covariance matrix can be calculated by scaling each row of X by wti and calculating an unweighted covariance matrix about θ.
In order to make the estimate asymptotically unbiased under a Normal model a correction factor, τ2, is needed. The value of the correction factor will depend on the functions employed (see Huber (1981) and Marazzi (1987)).
nag_robust_m_corr_user_fn_no_derr (g02hmc) finds A using the iterative procedure as given by Huber; see Huber (1981).
Ak=Sk+IAk-1
and
θjk=bjD1+θjk- 1,
where Sk=sjl, for j=1,2,,m and l=1,2,,m is a lower triangular matrix such that
sjl= -minmaxhjl/D2,-BL,BL, j>l -minmax12hjj/D2-1,-BD,BD, j=l ,
where and BD and BL are suitable bounds.
The value of τ may be chosen so that C is unbiased if the observations are from a given distribution.
nag_robust_m_corr_user_fn_no_derr (g02hmc) is based on routines in ROBETH; see Marazzi (1987).

4  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

5  Arguments

1:     orderNag_OrderTypeInput
On entry: the order argument specifies the two-dimensional storage scheme being used, i.e., row-major ordering or column-major ordering. C language defined storage is specified by order=Nag_RowMajor. See Section 3.2.1.3 in the Essential Introduction for a more detailed explanation of the use of this argument.
Constraint: order=Nag_RowMajor or Nag_ColMajor.
2:     ucvfunction, supplied by the userExternal Function
ucv must return the values of the functions u and w for a given value of its argument.
The specification of ucv is:
void  ucv (double t, double *u, double *w, Nag_Comm *comm)
1:     tdoubleInput
On entry: the argument for which the functions u and w must be evaluated.
2:     udouble *Output
On exit: the value of the u function at the point t.
Constraint: u0.0.
3:     wdouble *Output
On exit: the value of the w function at the point t.
Constraint: w0.0.
4:     commNag_Comm *
Pointer to structure of type Nag_Comm; the following members are relevant to ucv.
userdouble *
iuserInteger *
pPointer 
The type Pointer will be void *. Before calling nag_robust_m_corr_user_fn_no_derr (g02hmc) you may allocate memory and initialize these pointers with various quantities for use by ucv when called from nag_robust_m_corr_user_fn_no_derr (g02hmc) (see Section 3.2.1 in the Essential Introduction).
3:     indmIntegerInput
On entry: indicates which form of the function v will be used.
indm=1
v=1.
indm1
v=u.
4:     nIntegerInput
On entry: n, the number of observations.
Constraint: n>1.
5:     mIntegerInput
On entry: m, the number of columns of the matrix X, i.e., number of independent variables.
Constraint: 1mn.
6:     x[dim]const doubleInput
Note: the dimension, dim, of the array x must be at least
  • max1,pdx×m when order=Nag_ColMajor;
  • max1,n×pdx when order=Nag_RowMajor.
Where Xi,j appears in this document, it refers to the array element
  • x[j-1×pdx+i-1] when order=Nag_ColMajor;
  • x[i-1×pdx+j-1] when order=Nag_RowMajor.
On entry: Xi,j must contain the ith observation on the jth variable, for i=1,2,,n and j=1,2,,m.
7:     pdxIntegerInput
On entry: the stride separating row or column elements (depending on the value of order) in the array x.
Constraints:
  • if order=Nag_ColMajor, pdxn;
  • if order=Nag_RowMajor, pdxm.
8:     cov[m×m+1/2]doubleOutput
On exit: a robust estimate of the covariance matrix, C. The upper triangular part of the matrix C is stored packed by columns (lower triangular stored by rows), that is Cij is returned in cov[j×j-1/2+i-1], ij.
9:     a[m×m+1/2]doubleInput/Output
On entry: an initial estimate of the lower triangular real matrix A. Only the lower triangular elements must be given and these should be stored row-wise in the array.
The diagonal elements must be 0, and in practice will usually be >0. If the magnitudes of the columns of X are of the same order, the identity matrix will often provide a suitable initial value for A. If the columns of X are of different magnitudes, the diagonal elements of the initial value of A should be approximately inversely proportional to the magnitude of the columns of X.
Constraint: a[j×j-1/2+j]0.0, for j=0,1,,m-1.
On exit: the lower triangular elements of the inverse of the matrix A, stored row-wise.
10:   wt[n]doubleOutput
On exit: wt[i-1] contains the weights, wti=uzi2, for i=1,2,,n.
11:   theta[m]doubleInput/Output
On entry: an initial estimate of the location argument, θj, for j=1,2,,m.
In many cases an initial estimate of θj=0, for j=1,2,,m, will be adequate. Alternatively medians may be used as given by nag_median_1var (g07dac).
On exit: contains the robust estimate of the location argument, θj, for j=1,2,,m.
12:   bldoubleInput
On entry: the magnitude of the bound for the off-diagonal elements of Sk, BL.
Suggested value: bl=0.9.
Constraint: bl>0.0.
13:   bddoubleInput
On entry: the magnitude of the bound for the diagonal elements of Sk, BD.
Suggested value: bd=0.9.
Constraint: bd>0.0.
14:   maxitIntegerInput
On entry: the maximum number of iterations that will be used during the calculation of A.
Suggested value: maxit=150.
Constraint: maxit>0.
15:   nitmonIntegerInput
On entry: indicates the amount of information on the iteration that is printed.
nitmon>0
The value of A, θ and δ (see Section 7) will be printed at the first and every nitmon iterations.
nitmon0
No iteration monitoring is printed.
16:   outfileconst char *Input
On entry: a null terminated character string giving the name of the file to which results should be printed. If outfile=NULL or an empty string then the stdout stream is used. Note that the file will be opened in the append mode.
17:   toldoubleInput
On entry: the relative precision for the final estimate of the covariance matrix. Iteration will stop when maximum δ (see Section 7) is less than tol.
Constraint: tol>0.0.
18:   nitInteger *Output
On exit: the number of iterations performed.
19:   commNag_Comm *Communication Structure
The NAG communication argument (see Section 3.2.1.1 in the Essential Introduction).
20:   failNagError *Input/Output
The NAG error argument (see Section 3.6 in the Essential Introduction).

6  Error Indicators and Warnings

NE_ALLOC_FAIL
Dynamic memory allocation failed.
NE_BAD_PARAM
On entry, argument value had an illegal value.
NE_CONST_COL
Column value of x has constant value.
NE_CONVERGENCE
Iterations to calculate weights failed to converge.
NE_FUN_RET_VAL
u value returned by ucv<0.0: uvalue=value.
w value returned by ucv<0.0: wvalue=value.
NE_INT
On entry, m=value.
Constraint: m1.
On entry, maxit=value.
Constraint: maxit>0.
On entry, n=value.
Constraint: n>1.
On entry, n=value.
Constraint: n2.
On entry, pdx=value.
Constraint: pdx>0.
NE_INT_2
On entry, m=value and n=value.
Constraint: 1mn.
On entry, n=value and m=value.
Constraint: nm.
On entry, pdx=value and m=value.
Constraint: pdxm.
On entry, pdx=value and n=value.
Constraint: pdxn.
NE_INTERNAL_ERROR
An internal error has occurred in this function. Check the function call and any array sizes. If the call is correct then please contact NAG for assistance.
NE_NOT_CLOSE_FILE
Cannot close file value.
NE_NOT_WRITE_FILE
Cannot open file value for writing.
NE_REAL
On entry, bd=value.
Constraint: bd>0.0.
On entry, bl=value.
Constraint: bl>0.0.
On entry, tol=value.
Constraint: tol>0.0.
NE_ZERO_DIAGONAL
On entry, diagonal element value of a is 0.0.
NE_ZERO_SUM
Sum of u's (D2) is zero.
Sum of w's (D1) is zero.

7  Accuracy

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

8  Further Comments

The existence of A will depend upon the function u (see Marazzi (1987)); also if X is not of full rank a value of A will not be found. If the columns of X are almost linearly related, then convergence will be slow.
If derivatives of the u and w functions are available then the method used in nag_robust_m_corr_user_fn (g02hlc) will usually give much faster convergence.

9  Example

A sample of 10 observations on three variables is read in along with initial values for A and θ and argument values for the u and w functions, cu and cw. The covariance matrix computed by nag_robust_m_corr_user_fn_no_derr (g02hmc) is printed along with the robust estimate of θ.
ucv computes the Huber's weight functions:
ut=1, if  tcu2 ut= cut2, if  t>cu2
and
wt= 1, if   tcw wt= cwt, if   t>cw.

9.1  Program Text

Program Text (g02hmce.c)

9.2  Program Data

Program Data (g02hmce.d)

9.3  Program Results

Program Results (g02hmce.r)


nag_robust_m_corr_user_fn_no_derr (g02hmc) (PDF version)
g02 Chapter Contents
g02 Chapter Introduction
NAG C Library Manual

© The Numerical Algorithms Group Ltd, Oxford, UK. 2012