g02 Chapter Contents
g02 Chapter Introduction
NAG Library Manual

# NAG Library Function Documentnag_robust_m_corr_user_fn_no_derr (g02hmc)

## 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 #include
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, $\theta$, are given by
 $C=τ2ATA-1,$
where ${\tau }^{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
 $1n∑i=1nuzi2zi ziT -vzi2I=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\left(t\right)=1$ for all $t$, (ii) $v\left(t\right)=u\left(t\right)$.
The robust covariance matrix may be calculated from a weighted sum of squares and cross-products matrix about $\theta$ using weights ${\mathit{wt}}_{i}=u\left(‖{z}_{i}‖\right)$. In case (i) a divisor of $n$ is used and in case (ii) a divisor of $\sum _{i=1}^{n}{\mathit{wt}}_{i}$ is used. If $w\left(.\right)=\sqrt{u\left(.\right)}$, then the robust covariance matrix can be calculated by scaling each row of $X$ by $\sqrt{{\mathit{wt}}_{i}}$ and calculating an unweighted covariance matrix about $\theta$.
In order to make the estimate asymptotically unbiased under a Normal model a correction factor, ${\tau }^{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 ${S}_{k}=\left({s}_{jl}\right)$, for $\mathit{j}=1,2,\dots ,m$ and $\mathit{l}=1,2,\dots ,m$ is a lower triangular matrix such that
 $sjl= -minmaxhjl/D2,-BL,BL, j>l -minmax12hjj/D2-1,-BD,BD, j=l ,$
where
• ${D}_{1}=\sum _{i=1}^{n}w\left({‖{z}_{i}‖}_{2}\right)$
• ${D}_{2}=\sum _{i=1}^{n}u\left({‖{z}_{i}‖}_{2}\right)$
• ${h}_{jl}=\sum _{i=1}^{n}u\left({‖{z}_{i}‖}_{2}\right){z}_{ij}{z}_{il}$, for $j\ge l$
• ${b}_{j}=\sum _{i=1}^{n}w\left({‖{z}_{i}‖}_{2}\right)\left({x}_{ij}-{b}_{j}\right)$
and $\mathit{BD}$ and $\mathit{BL}$ are suitable bounds.
The value of $\tau$ 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 ${\mathbf{order}}=\mathrm{Nag_RowMajor}$. See Section 3.2.1.3 in the Essential Introduction for a more detailed explanation of the use of this argument.
Constraint: ${\mathbf{order}}=\mathrm{Nag_RowMajor}$ or $\mathrm{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: ${\mathbf{u}}\ge 0.0$.
3:     wdouble *Output
On exit: the value of the $w$ function at the point t.
Constraint: ${\mathbf{w}}\ge 0.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.1 in the Essential Introduction).
3:     indmIntegerInput
On entry: indicates which form of the function $v$ will be used.
${\mathbf{indm}}=1$
$v=1$.
${\mathbf{indm}}\ne 1$
$v=u$.
4:     nIntegerInput
On entry: $n$, the number of observations.
Constraint: ${\mathbf{n}}>1$.
5:     mIntegerInput
On entry: $m$, the number of columns of the matrix $X$, i.e., number of independent variables.
Constraint: $1\le {\mathbf{m}}\le {\mathbf{n}}$.
6:     x[$\mathit{dim}$]const doubleInput
Note: the dimension, dim, of the array x must be at least
• $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{pdx}}×{\mathbf{m}}\right)$ when ${\mathbf{order}}=\mathrm{Nag_ColMajor}$;
• $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}×{\mathbf{pdx}}\right)$ when ${\mathbf{order}}=\mathrm{Nag_RowMajor}$.
Where ${\mathbf{X}}\left(i,j\right)$ appears in this document, it refers to the array element
• ${\mathbf{x}}\left[\left(j-1\right)×{\mathbf{pdx}}+i-1\right]$ when ${\mathbf{order}}=\mathrm{Nag_ColMajor}$;
• ${\mathbf{x}}\left[\left(i-1\right)×{\mathbf{pdx}}+j-1\right]$ when ${\mathbf{order}}=\mathrm{Nag_RowMajor}$.
On entry: ${\mathbf{X}}\left(\mathit{i},\mathit{j}\right)$ must contain the $\mathit{i}$th observation on the $\mathit{j}$th variable, for $\mathit{i}=1,2,\dots ,n$ and $\mathit{j}=1,2,\dots ,m$.
7:     pdxIntegerInput
On entry: the stride separating row or column elements (depending on the value of order) in the array x.
Constraints:
• if ${\mathbf{order}}=\mathrm{Nag_ColMajor}$, ${\mathbf{pdx}}\ge {\mathbf{n}}$;
• if ${\mathbf{order}}=\mathrm{Nag_RowMajor}$, ${\mathbf{pdx}}\ge {\mathbf{m}}$.
8:     cov[${\mathbf{m}}×\left({\mathbf{m}}+1\right)/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}_{ij}$ is returned in ${\mathbf{cov}}\left[j×\left(j-1\right)/2+i-1\right]$, $i\le j$.
9:     a[${\mathbf{m}}×\left({\mathbf{m}}+1\right)/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 $\text{}\ne 0$, and in practice will usually be $\text{}>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: ${\mathbf{a}}\left[\mathit{j}×\left(\mathit{j}-1\right)/2+\mathit{j}\right]\ne 0.0$, for $\mathit{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: ${\mathbf{wt}}\left[\mathit{i}-1\right]$ contains the weights, ${\mathit{wt}}_{\mathit{i}}=u\left({‖{z}_{\mathit{i}}‖}_{2}\right)$, for $\mathit{i}=1,2,\dots ,n$.
11:   theta[m]doubleInput/Output
On entry: an initial estimate of the location argument, ${\theta }_{\mathit{j}}$, for $\mathit{j}=1,2,\dots ,m$.
In many cases an initial estimate of ${\theta }_{\mathit{j}}=0$, for $\mathit{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, ${\theta }_{\mathit{j}}$, for $\mathit{j}=1,2,\dots ,m$.
12:   bldoubleInput
On entry: the magnitude of the bound for the off-diagonal elements of ${S}_{k}$, $BL$.
Suggested value: ${\mathbf{bl}}=0.9$.
Constraint: ${\mathbf{bl}}>0.0$.
13:   bddoubleInput
On entry: the magnitude of the bound for the diagonal elements of ${S}_{k}$, $BD$.
Suggested value: ${\mathbf{bd}}=0.9$.
Constraint: ${\mathbf{bd}}>0.0$.
14:   maxitIntegerInput
On entry: the maximum number of iterations that will be used during the calculation of $A$.
Suggested value: ${\mathbf{maxit}}=150$.
Constraint: ${\mathbf{maxit}}>0$.
15:   nitmonIntegerInput
On entry: indicates the amount of information on the iteration that is printed.
${\mathbf{nitmon}}>0$
The value of $A$, $\theta$ and $\delta$ (see Section 7) will be printed at the first and every nitmon iterations.
${\mathbf{nitmon}}\le 0$
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 ${\mathbf{outfile}}=\mathbf{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 $\delta$ (see Section 7) is less than tol.
Constraint: ${\mathbf{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.
On entry, argument $⟨\mathit{\text{value}}⟩$ had an illegal value.
NE_CONST_COL
Column $⟨\mathit{\text{value}}⟩$ of x has constant value.
NE_CONVERGENCE
Iterations to calculate weights failed to converge.
NE_FUN_RET_VAL
$u$ value returned by ${\mathbf{ucv}}<0.0$: $u\left(⟨\mathit{\text{value}}⟩\right)=⟨\mathit{\text{value}}⟩$.
$w$ value returned by ${\mathbf{ucv}}<0.0$: $w\left(⟨\mathit{\text{value}}⟩\right)=⟨\mathit{\text{value}}⟩$.
NE_INT
On entry, ${\mathbf{m}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{m}}\ge 1$.
On entry, ${\mathbf{maxit}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{maxit}}>0$.
On entry, ${\mathbf{n}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{n}}>1$.
On entry, ${\mathbf{pdx}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{pdx}}>0$.
NE_INT_2
On entry, ${\mathbf{m}}=⟨\mathit{\text{value}}⟩$ and ${\mathbf{n}}=⟨\mathit{\text{value}}⟩$.
Constraint: $1\le {\mathbf{m}}\le {\mathbf{n}}$.
On entry, ${\mathbf{n}}=⟨\mathit{\text{value}}⟩$ and ${\mathbf{m}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{n}}\ge {\mathbf{m}}$.
On entry, ${\mathbf{pdx}}=⟨\mathit{\text{value}}⟩$ and ${\mathbf{m}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{pdx}}\ge {\mathbf{m}}$.
On entry, ${\mathbf{pdx}}=⟨\mathit{\text{value}}⟩$ and ${\mathbf{n}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{pdx}}\ge {\mathbf{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 $⟨\mathit{\text{value}}⟩$.
NE_NOT_WRITE_FILE
Cannot open file $⟨\mathit{\text{value}}⟩$ for writing.
NE_REAL
On entry, ${\mathbf{bd}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{bd}}>0.0$.
On entry, ${\mathbf{bl}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{bl}}>0.0$.
On entry, ${\mathbf{tol}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{tol}}>0.0$.
NE_ZERO_DIAGONAL
On entry, diagonal element $⟨\mathit{\text{value}}⟩$ of a is $0.0$.
NE_ZERO_SUM
Sum of $u$'s ($\mathit{D2}$) is zero.
Sum of $w$'s ($\mathit{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=\text{}$ the maximum value of $\left|{s}_{jl}\right|$ (ii) $d2=\text{}$ the maximum absolute change in $wt\left(i\right)$ (iii) $d3=\text{}$ the maximum absolute relative change in ${\theta }_{j}$
and let $\delta =\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(d1,d2,d3\right)$. Then the iterative procedure is assumed to have converged when $\delta <{\mathbf{tol}}$.

## 8  Parallelism and Performance

nag_robust_m_corr_user_fn_no_derr (g02hmc) is not threaded by NAG in any implementation.
nag_robust_m_corr_user_fn_no_derr (g02hmc) makes calls to BLAS and/or LAPACK routines, which may be threaded within the vendor library used by this implementation. Consult the documentation for the vendor library for further information.

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.

## 10  Example

A sample of $10$ observations on three variables is read in along with initial values for $A$ and $\theta$ 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 $\theta$.
ucv computes the Huber's weight functions:
 $ut=1, if t≤cu2 ut= cut2, if t>cu2$
and
 $wt= 1, if t≤cw wt= cwt, if t>cw.$

### 10.1  Program Text

Program Text (g02hmce.c)

### 10.2  Program Data

Program Data (g02hmce.d)

### 10.3  Program Results

Program Results (g02hmce.r)