NAG Library Function Document

nag_robust_m_corr_user_fn_no_derr (g02hmc)

+− Contents

1 Purpose

2 Specification

3 Description

4 References

5 Arguments

6 Error Indicators and Warnings

7 Accuracy

8 Further Comments

+− 9 Example

9.1 Program Text

9.2 Program Data

9.3 Program Results

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 = τ^{2} {(A^{T} A)}^{- 1},

where

τ^{2}

is a correction factor and

A

is a lower triangular matrix found as the solution to the following equations.

z_{i} = A (x_{i} - θ)

\frac{1}{n} \sum_{i = 1}^{n} w ({‖z_{i}‖}_{2}) z_{i} = 0

and

\frac{1}{n} \sum_{i = 1}^{n} u ({‖z_{i}‖}_{2}) z_{i} z_{i}^{T} - v ({‖z_{i}‖}_{2}) I = 0,

where	$x_{i}$ is a vector of length $m$ containing the elements of the $i$ th row of $X$ ,
	$z_{i}$ 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)	$v (t) = 1$ for all $t$ ,
(ii)	$v (t) = u (t)$ .

The robust covariance matrix may be calculated from a weighted sum of squares and cross-products matrix about

θ

using weights

{wt}_{i} = u (‖z_{i}‖)

. In case (i) a divisor of

n

is used and in case (ii) a divisor of

\sum_{i = 1}^{n} {wt}_{i}

is used. If

w (.) = \sqrt{u (.)}

, then the robust covariance matrix can be calculated by scaling each row of

X

\sqrt{{wt}_{i}}

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).

A_{k} = (S_{k} + I) A_{k - 1}

and

θ_{j_{k}} = \frac{b_{j}}{D_{1}} + θ_{j_{k - 1}},

where

S_{k} = (s_{j l})

, for

j = 1, 2, \dots, m

and

l = 1, 2, \dots, m

is a lower triangular matrix such that

s_{j l} = \{\begin{cases} - \min [\max (h_{j l} / D_{2}, - BL), BL], & j > l \\ - \min [\max (\frac{1}{2} (h_{j j} / D_{2} - 1), - BD), BD], & j = l \end{cases},

where

$D_{1} = \sum_{i = 1}^{n} w ({‖z_{i}‖}_{2})$
$D_{2} = \sum_{i = 1}^{n} u ({‖z_{i}‖}_{2})$
$h_{j l} = \sum_{i = 1}^{n} u ({‖z_{i}‖}_{2}) z_{i j} z_{i l}$ , for $j \geq l$
$b_{j} = \sum_{i = 1}^{n} w ({‖z_{i}‖}_{2}) (x_{i j} - b_{j})$

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: order – Nag_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: ucv – function, 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: t – doubleInput

On entry: the argument for which the functions

u

and

w

must be evaluated.

2: u – double *Output

On exit: the value of the

u

function at the point t.

Constraint:

u \geq 0.0

3: w – double *Output

On exit: the value of the

w

function at the point t.

Constraint:

w \geq 0.0

4: comm – Nag_Comm *

Pointer to structure of type Nag_Comm; the following members are relevant to ucv.

user – double *
iuser – Integer *
p – Pointer: 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: indm – IntegerInput

On entry: indicates which form of the function

v

will be used.

$indm = 1$: $v = 1$ .
$indm \neq 1$: $v = u$ .

4: n – IntegerInput

On entry:

n

, the number of observations.

Constraint:

n > 1

5: m – IntegerInput

On entry:

m

, the number of columns of the matrix

X

, i.e., number of independent variables.

Constraint:

1 \leq m \leq n

6: x[ $\dim$ ] – const doubleInput

Note: the dimension, dim, of the array x must be at least

$\max (1, pdx \times m)$ when $order = Nag_ColMajor$ ;
$\max (1, n \times pdx)$ when $order = Nag_RowMajor$ .

Where

X (i, j)

appears in this document, it refers to the array element

$x [(j - 1) \times pdx + i - 1]$ when $order = Nag_ColMajor$ ;
$x [(i - 1) \times pdx + j - 1]$ when $order = Nag_RowMajor$ .

On entry:

X (i, j)

must contain the

i

th observation on the

j

th variable, for

i = 1, 2, \dots, n

and

j = 1, 2, \dots, m

7: pdx – IntegerInput

On entry: the stride separating row or column elements (depending on the value of order) in the array x.

Constraints:

if $order = Nag_ColMajor$ , $pdx \geq n$ ;
if $order = Nag_RowMajor$ , $pdx \geq m$ .

8: cov[ $m \times (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

C_{i j}

is returned in

cov [j \times (j - 1) / 2 + i - 1]

i \leq j

9: a[ $m \times (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

\neq 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 \times (j - 1) / 2 + j] \neq 0.0

, for

j = 0, 1, \dots, 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,

{wt}_{i} = u ({‖z_{i}‖}_{2})

, for

i = 1, 2, \dots, n

11: theta[m] – doubleInput/Output

On entry: an initial estimate of the location argument,

θ_{j}

, for

j = 1, 2, \dots, m

In many cases an initial estimate of

θ_{j} = 0

, for

j = 1, 2, \dots, 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, \dots, m

12: bl – doubleInput

On entry: the magnitude of the bound for the off-diagonal elements of

S_{k}

B L

Suggested value:

bl = 0.9

Constraint:

bl > 0.0

13: bd – doubleInput

On entry: the magnitude of the bound for the diagonal elements of

S_{k}

B D

Suggested value:

bd = 0.9

Constraint:

bd > 0.0

14: maxit – IntegerInput

On entry: the maximum number of iterations that will be used during the calculation of

A

Suggested value:

maxit = 150

Constraint:

maxit > 0

15: nitmon – IntegerInput

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.
$nitmon \leq 0$: No iteration monitoring is printed.

16: outfile – const 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: tol – doubleInput

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: nit – Integer *Output

On exit: the number of iterations performed.

19: comm – Nag_Comm *Communication Structure

The NAG communication argument (see Section 3.2.1.1 in the Essential Introduction).

20: fail – NagError *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$ : $u (〈value〉) = 〈value〉$ .; $w$ value returned by $ucv < 0.0$ : $w (〈value〉) = 〈value〉$ .
NE_INT: On entry, $m = 〈value〉$ .
Constraint: $m \geq 1$ .; On entry, $maxit = 〈value〉$ .
Constraint: $maxit > 0$ .; On entry, $n = 〈value〉$ .
Constraint: $n > 1$ .; On entry, $n = 〈value〉$ .
Constraint: $n \geq 2$ .; On entry, $pdx = 〈value〉$ .
Constraint: $pdx > 0$ .
NE_INT_2: On entry, $m = 〈value〉$ and $n = 〈value〉$ .
Constraint: $1 \leq m \leq n$ .; On entry, $n = 〈value〉$ and $m = 〈value〉$ .
Constraint: $n \geq m$ .; On entry, $pdx = 〈value〉$ and $m = 〈value〉$ .
Constraint: $pdx \geq m$ .; On entry, $pdx = 〈value〉$ and $n = 〈value〉$ .
Constraint: $pdx \geq n$ .
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)	$d 1 =$ the maximum value of $\|s_{j l}\|$
(ii)	$d 2 =$ the maximum absolute change in $w t (i)$
(iii)	$d 3 =$ the maximum absolute relative change in $θ_{j}$

and let

δ = \max (d 1, d 2, d 3)

. 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,

c_{u}

and

c_{w}

. 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:

\begin{array}{l} u (t) = 1, & if t \leq c_{u}^{2} \\ u (t) = \frac{c_{u}}{t^{2}}, & if t > c_{u}^{2} \end{array}

and

\begin{array}{l} w (t) = 1, & if t \leq c_{w} \\ w (t) = \frac{c_{w}}{t}, & if t > c_{w} . \end{array}

NAG Library Function Documentnag_robust_m_corr_user_fn_no_derr (g02hmc)

+− Contents

1 Purpose

2 Specification

3 Description

4 References

5 Arguments

6 Error Indicators and Warnings

7 Accuracy

8 Further Comments

9 Example

9.1 Program Text

9.2 Program Data

9.3 Program Results

NAG Library Function Document

nag_robust_m_corr_user_fn_no_derr (g02hmc)