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)

NAG Toolbox: nag_correg_robustm_corr_user_deriv (g02hl)

Purpose

nag_correg_robustm_corr_user_deriv (g02hl) calculates a robust estimate of the covariance matrix for user-supplied weight functions and their derivatives.

Syntax

[user, cov, a, wt, theta, nit, ifail] = g02hl(ucv, indm, x, a, theta, 'user', user, 'n', n, 'm', m, 'bl', bl, 'bd', bd, 'maxit', maxit, 'nitmon', nitmon, 'tol', tol)

[user, cov, a, wt, theta, nit, ifail] = nag_correg_robustm_corr_user_deriv(ucv, indm, x, a, theta, 'user', user, 'n', n, 'm', m, 'bl', bl, 'bd', bd, '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 n

n

observations on m

m

variables in a matrix X

X

, a robust estimate of the covariance matrix, C

C

, and a robust estimate of location, θ

θ

, are given by:

C = τ²(A^TA)^{− 1},

C = τ^{2} {(A^{T} A)}^{- 1},

where τ²

τ^{2}

is a correction factor and A

A

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

z_i = A(x_i − θ)

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

	n
1/n	∑	w(‖z_i‖₂)z_i = 0
	i = 1

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

and

	n
1/n	∑	u(‖z_i‖₂)z_iz_i^T − v(‖z_i‖₂)I = 0,
	i = 1

\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 $x_{i}$ is a vector of length m $m$ containing the elements of the i $i$ th row of X $X$ ,
	z_i $z_{i}$ is a vector of length m $m$ ,
	I $I$ is the identity matrix and 0 $0$ is the zero matrix,
and	w $w$ and u $u$ are suitable functions.

nag_correg_robustm_corr_user_deriv (g02hl) covers two situations:

(i)	v(t) = 1 $v (t) = 1$ for all t $t$ ,
(ii)	v(t) = u(t) $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‖)

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

. In case (i) a divisor of n

n

is used and in case (ii) a divisor of ∑ _{i = 1}ⁿwt_i

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

is used. If w( . ) = sqrt(u( . ))

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

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

X

by sqrt(wt_i)

\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_correg_robustm_corr_user_deriv (g02hl) finds A

A

using the iterative procedure as given by Huber.

A_k = (S_k + I)A_{k − 1}

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

and

θ_{j_k} = (b_j)/(D₁) + θ_{j_{k − 1}},

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

where S_k = (s_jl)

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

, for j = 1,2, … ,m

j = 1, 2, \dots, m

and l = 1,2, … ,m

l = 1, 2, \dots, m

, is a lower triangular matrix such that:

sjl =

{	− min [max (hjl / D3, − BL),BL], j > l − min [max ((hjj / (2D3 − D4 / D2)), − BD),BD], j = l

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

where

D₁ = ∑ _{i = 1}ⁿ {w(‖z_i‖₂) + 1/mw^′(‖z_i‖₂)‖z_i‖₂} $D_{1} = \sum_{i = 1}^{n} {w ({‖ z_{i} ‖}_{2}) + \frac{1}{m} w^{'} ({‖ z_{i} ‖}_{2}) {‖ z_{i} ‖}_{2}}$
D₂ = ∑ _{i = 1}ⁿ {1/p(u^′(‖z_i‖₂)‖z_i‖₂ + 2u(‖z_i‖₂))‖z_i‖₂ − v^′(‖z_i‖₂)} ‖z_i‖₂ $D_{2} = \sum_{i = 1}^{n} {\frac{1}{p} (u^{'} ({‖ z_{i} ‖}_{2}) {‖ z_{i} ‖}_{2} + 2 u ({‖ z_{i} ‖}_{2})) {‖ z_{i} ‖}_{2} - v^{'} ({‖ z_{i} ‖}_{2})} {‖ z_{i} ‖}_{2}$
D₃ = 1/(m + 2) ∑ _{i = 1}ⁿ {1/m(u^′(‖z_i‖₂)‖z_i‖₂ + 2u(‖z_i‖₂)) + u(‖z_i‖₂)} ‖z_i‖₂² $D_{3} = \frac{1}{m + 2} \sum_{i = 1}^{n} {\frac{1}{m} (u^{'} ({‖ z_{i} ‖}_{2}) {‖ z_{i} ‖}_{2} + 2 u ({‖ z_{i} ‖}_{2})) + u ({‖ z_{i} ‖}_{2})} {‖ z_{i} ‖}_{2}^{2}$
D₄ = ∑ _{i = 1}ⁿ {1/mu(‖z_i‖₂)‖z_i‖₂² − v(‖z_i‖₂²)} $D_{4} = \sum_{i = 1}^{n} {\frac{1}{m} u ({‖ z_{i} ‖}_{2}) {‖ z_{i} ‖}_{2}^{2} - v ({‖ z_{i} ‖}_{2}^{2})}$
h_jl = ∑ _{i = 1}ⁿu(‖z_i‖₂)z_ijz_il $h_{j l} = \sum_{i = 1}^{n} u ({‖ z_{i} ‖}_{2}) z_{i j} z_{i l}$ , for j > l $j > l$
h_jj = ∑ _{i = 1}ⁿu(‖z_i‖₂)(z_ij² − ‖z_ij‖₂² / m) $h_{j j} = \sum_{i = 1}^{n} u ({‖ z_{i} ‖}_{2}) (z_{i j}^{2} - {‖ z_{i j} ‖}_{2}^{2} / m)$
b_j = ∑ _{i = 1}ⁿw(‖z_i‖₂)(x_ij − b_j) $b_{j} = \sum_{i = 1}^{n} w ({‖ z_{i} ‖}_{2}) (x_{i j} - b_{j})$
and BD $BD$ and BL $BL$ are suitable bounds.

nag_correg_robustm_corr_user_deriv (g02hl) 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: ucv – function handle or string containing name of m-file

ucv must return the values of the functions u

u

and w

w

and their derivatives for a given value of its argument.

[user, u, ud, w, wd] = ucv(t, user)

Input Parameters

1: t – double scalar: The argument for which the functions u $u$ and w $w$ must be evaluated.
2: user – Any MATLAB object: ucv is called from nag_correg_robustm_corr_user_deriv (g02hl) with the object supplied to nag_correg_robustm_corr_user_deriv (g02hl).

Output Parameters

1: user – Any MATLAB object
2: u – double scalar: The value of the u $u$ function at the point t.
3: ud – double scalar: The value of the derivative of the u $u$ function at the point t.
4: w – double scalar: The value of the w $w$ function at the point t.
5: wd – double scalar: The value of the derivative of the w $w$ function at the point t.

2: indm – int64int32nag_int scalar

Indicates which form of the function v

v

will be used.

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

3: x(ldx,m) – double array

ldx, the first dimension of the array, must satisfy the constraint ldx ≥ n

ldx \geq n

x(i,j)

x (i, j)

must contain the i

i

th observation on the j

j

th variable, for i = 1,2, … ,n

i = 1, 2, \dots, n

and j = 1,2, … ,m

j = 1, 2, \dots, m

4: a(m × (m + 1) / 2 $m \times (m + 1) / 2$ ) – double array

An initial estimate of the lower triangular real matrix A

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

\neq 0

, and in practice will usually be > 0

> 0

. If the magnitudes of the columns of X

X

are of the same order, the identity matrix will often provide a suitable initial value for A

A

. If the columns of X

X

are of different magnitudes, the diagonal elements of the initial value of A

A

should be approximately inversely proportional to the magnitude of the columns of X

X

Constraint: a(j × (j − 1) / 2 + j) ≠ 0.0

a (j \times (j - 1) / 2 + j) \neq 0.0

, for j = 1,2, … ,m

j = 1, 2, \dots, m

5: theta(m) – double array

m, the dimension of the array, must satisfy the constraint 1 ≤ m ≤ n

1 \leq m \leq n

An initial estimate of the location parameter, θ_j

θ_{j}

, for j = 1,2, … ,m

j = 1, 2, \dots, m

In many cases an initial estimate of θ_j = 0

θ_{j} = 0

, for j = 1,2, … ,m

j = 1, 2, \dots, m

, will be adequate. Alternatively medians may be used as given by nag_univar_robust_1var_median (g07da).

Optional Input Parameters

1: user – Any MATLAB object

user is not used by nag_correg_robustm_corr_user_deriv (g02hl), but is passed to ucv. Note that for large objects it may be more efficient to use a global variable which is accessible from the m-files than to use user.

2: n – int64int32nag_int scalar

Default: The first dimension of the array x.

n

, the number of observations.

Constraint: n > 1

n > 1

3: m – int64int32nag_int scalar

Default: The dimension of the array theta and the second dimension of the array x. (An error is raised if these dimensions are not equal.)

m

, the number of columns of the matrix X

X

, i.e., number of independent variables.

Constraint: 1 ≤ m ≤ n

1 \leq m \leq n

4: bl – double scalar

The magnitude of the bound for the off-diagonal elements of S_k

S_{k}

, BL

B L

Default: 0.9

0.9

Constraint: bl > 0.0

bl > 0.0

5: bd – double scalar

The magnitude of the bound for the diagonal elements of S_k

S_{k}

, BD

B D

Default: 0.9

0.9

Constraint: bd > 0.0

bd > 0.0

6: maxit – int64int32nag_int scalar

The maximum number of iterations that will be used during the calculation of A

A

Default: 150

150

Constraint: maxit > 0

maxit > 0

7: nitmon – int64int32nag_int scalar

Indicates the amount of information on the iteration that is printed.

nitmon > 0 $nitmon > 0$: The value of A $A$ , θ $θ$ and δ $δ$ (see Section [Accuracy]) will be printed at the first and every nitmon iterations.
nitmon ≤ 0 $nitmon \leq 0$: 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: 0

0

8: tol – double scalar

The relative precision for the final estimates of the covariance matrix. Iteration will stop when maximum δ

δ

(see Section [Accuracy]) is less than tol.

Default: 5e-5

5e-5

Constraint: tol > 0.0

tol > 0.0

Input Parameters Omitted from the MATLAB Interface

ruser ldx wk

Output Parameters

1: user – Any MATLAB object
2: cov(m × (m + 1) / 2 $m \times (m + 1) / 2$ ) – double array: Contains a robust estimate of the covariance matrix, C $C$ . The upper triangular part of the matrix C $C$ is stored packed by columns (lower triangular stored by rows), C_ij $C_{i j}$ is returned in cov((j × (j − 1) / 2 + i)) $cov ((j \times (j - 1) / 2 + i))$ , i ≤ j $i \leq j$ .
3: a(m × (m + 1) / 2 $m \times (m + 1) / 2$ ) – double array: The lower triangular elements of the inverse of the matrix A $A$ , stored row-wise.
4: wt(n) – double array: wt(i) $wt (i)$ contains the weights, wt_i = u(‖z_i‖₂) ${wt}_{i} = u ({‖ z_{i} ‖}_{2})$ , for i = 1,2, … ,n $i = 1, 2, \dots, n$ .
5: theta(m) – double array: Contains the robust estimate of the location parameter, θ_j $θ_{j}$ , for j = 1,2, … ,m $j = 1, 2, \dots, m$ .
6: nit – int64int32nag_int scalar: The number of iterations performed.
7: ifail – int64int32nag_int scalar: ifail = 0 $ifail = 0$ unless the function detects an error (see [Error Indicators and Warnings]).

Error Indicators and Warnings

Errors or warnings detected by the function:

Cases prefixed with W are classified as warnings and do not generate an error of type NAG:error_n. See nag_issue_warnings.

ifail = 1 $ifail = 1$

On entry,	n ≤ 1 $n \leq 1$ ,
or	m < 1 $m < 1$ ,
or	n < m $n < m$ ,
or	ldx < n $ldx < n$ .

ifail = 2 $ifail = 2$

On entry,	tol ≤ 0.0 $tol \leq 0.0$ ,
or	maxit ≤ 0 $maxit \leq 0$ ,
or	diagonal element of a = 0.0 $a = 0.0$ ,
or	bl ≤ 0.0 $bl \leq 0.0$ ,
or	bd ≤ 0.0 $bd \leq 0.0$ .

ifail = 3 $ifail = 3$: A column of x has a constant value.

ifail = 4 $ifail = 4$: Value of u or w returned by ucv < 0 $ucv < 0$ .

ifail = 5 $ifail = 5$: The function has failed to converge in maxit iterations.

W ifail = 6 $ifail = 6$: One of the following is zero: D₁ $D_{1}$ , D₂ $D_{2}$ or D₃ $D_{3}$ .

This may be caused by the functions u $u$ or w $w$ being too strict for the current estimate of A $A$ (or C $C$ ). You should try either a larger initial estimate of A $A$ or make u $u$ and w $w$ less strict.

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 = $d 1 =$ the maximum value of \|s_jl\| $\| s_{j l} \|$
(ii)	d2 = $d 2 =$ the maximum absolute change in wt(i) $w t (i)$
(iii)	d3 = $d 3 =$ the maximum absolute relative change in θ_j $θ_{j}$

and let δ = max (d1,d2,d3)

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

. Then the iterative procedure is assumed to have converged when δ < tol

δ < tol

Further Comments

The existence of A

A

will depend upon the function u

u

(see Marazzi (1987)); also if X

X

is not of full rank a value of A

A

will not be found. If the columns of X

X

are almost linearly related, then convergence will be slow.

Example

Open in the MATLAB editor: nag_correg_robustm_corr_user_deriv_example

function nag_correg_robustm_corr_user_deriv_example
indm = int64(1);
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];
a = [1;
     0;
     1;
     0;
     0;
     1];
theta = [0;
     0;
     0];
user = [4,2];
[user, covar, aOut, wt, thetaOut, nit, ifail] = ...
    nag_correg_robustm_corr_user_deriv(@ucv, indm, x, a, theta, 'user', user)

function [user, u, ud, w, wd] = ucv(t, user)
  cu = user(1);
  u = 1.0;
  ud = 0.0;
  if (t ~= 0)
    t2 = t*t;
    if (t2 > cu)
      u = cu/t2;
      ud = -2.0*u/t;
    end
  end
  % w function and derivative
  cw = user(2);
  if (t > cw)
    w = cw/t;
    wd = -w/t;
  else
    w = 1.0;
    wd = 0.0;
  end

Open in the MATLAB editor: g02hl_example

function g02hl_example
indm = int64(1);
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];
a = [1;
     0;
     1;
     0;
     0;
     1];
theta = [0;
     0;
     0];
user = [4,2];
[user, covar, aOut, wt, thetaOut, nit, ifail] = ...
    g02hl(@ucv, indm, x, a, theta, 'user', user)

function [user, u, ud, w, wd] = ucv(t, user)
  cu = user(1);
  u = 1.0;
  ud = 0.0;
  if (t ~= 0)
    t2 = t*t;
    if (t2 > cu)
      u = cu/t2;
      ud = -2.0*u/t;
    end
  end
  % w function and derivative
  cw = user(2);
  if (t > cw)
    w = cw/t;
    wd = -w/t;
  else
    w = 1.0;
    wd = 0.0;
  end

PDF version (NAG web site, 64-bit version, 64-bit version)

Chapter Contents

Chapter Introduction

NAG Toolbox