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_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 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)ziziTv(zi2)I = 0,
i = 1
1ni=1nu(zi2)zi ziT -v(zi2)I=0,
where xixi is a vector of length mm containing the elements of the iith row of XX,
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_user_deriv (g02hl) covers two situations:
(i) v(t) = 1v(t)=1 for all tt,
(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 wti = u(zi)wti=u(zi). In case (i) a divisor of nn is used and in case (ii) a divisor of i = 1nwtii=1nwti is used. If w( . ) = sqrt(u( . ))w(.)=u(.), then the robust covariance matrix can be calculated by scaling each row of XX by sqrt(wti)wti and calculating an unweighted covariance matrix about θθ.
In order to make the estimate asymptotically unbiased under a Normal model a correction factor, τ2τ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 AA using the iterative procedure as given by Huber.
Ak = (Sk + I)Ak1
Ak=(Sk+I)Ak-1
and
θjk = (bj)/(D1) + θjk 1,
θjk=bjD1+θjk- 1,
where Sk = (sjl)Sk=(sjl), for j = 1,2,,mj=1,2,,m and l = 1,2,,ml=1,2,,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
,
sjl={ -min[max(hjl/D3,-BL),BL], j>l -min[max((hjj/(2D3-D4/D2)),-BD),BD], j=l ,
where 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 uu and ww 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 uu and ww 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 uu function at the point t.
3:     ud – double scalar
The value of the derivative of the uu function at the point t.
4:     w – double scalar
The value of the ww function at the point t.
5:     wd – double scalar
The value of the derivative of the ww function at the point t.
2:     indm – int64int32nag_int scalar
Indicates which form of the function vv will be used.
indm = 1indm=1
v = 1v=1.
indm1indm1
v = uv=u.
3:     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 on the jjth variable, for i = 1,2,,ni=1,2,,n and j = 1,2,,mj=1,2,,m.
4:     a(m × (m + 1) / 2m×(m+1)/2) – double array
An initial estimate of the lower triangular real matrix AA. Only the lower triangular elements must be given and these should be stored row-wise in the array.
The diagonal elements must be 00, and in practice will usually be > 0>0. If the magnitudes of the columns of XX are of the same order, the identity matrix will often provide a suitable initial value for AA. If the columns of XX are of different magnitudes, the diagonal elements of the initial value of AA should be approximately inversely proportional to the magnitude of the columns of XX.
Constraint: a(j × (j1) / 2 + j)0.0aj×(j-1)/2+j0.0, for j = 1,2,,mj=1,2,,m.
5:     theta(m) – double array
m, the dimension of the array, must satisfy the constraint 1mn1mn.
An initial estimate of the location parameter, θjθj, for j = 1,2,,mj=1,2,,m.
In many cases an initial estimate of θj = 0θj=0, for j = 1,2,,mj=1,2,,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.
nn, the number of observations.
Constraint: n > 1n>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.)
mm, the number of columns of the matrix XX, i.e., number of independent variables.
Constraint: 1mn1mn.
4:     bl – double scalar
The magnitude of the bound for the off-diagonal elements of SkSk, BLBL.
Default: 0.90.9
Constraint: bl > 0.0bl>0.0.
5:     bd – double scalar
The magnitude of the bound for the diagonal elements of SkSk, BDBD.
Default: 0.90.9
Constraint: bd > 0.0bd>0.0.
6:     maxit – int64int32nag_int scalar
The maximum number of iterations that will be used during the calculation of AA.
Default: 150150
Constraint: maxit > 0maxit>0.
7:     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
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-55e-5
Constraint: tol > 0.0tol>0.0.

Input Parameters Omitted from the MATLAB Interface

ruser ldx wk

Output Parameters

1:     user – Any MATLAB object
2:     cov(m × (m + 1) / 2m×(m+1)/2) – double array
Contains a robust estimate of the covariance matrix, CC. The upper triangular part of the matrix CC is stored packed by columns (lower triangular stored by rows), CijCij is returned in cov((j × (j1) / 2 + i))cov(j×(j-1)/2+i), ijij.
3:     a(m × (m + 1) / 2m×(m+1)/2) – double array
The lower triangular elements of the inverse of the matrix AA, stored row-wise.
4:     wt(n) – double array
wt(i)wti contains the weights, wti = u(zi2)wti=u(zi2), for i = 1,2,,ni=1,2,,n.
5:     theta(m) – double array
Contains the robust estimate of the location parameter, θjθj, for j = 1,2,,mj=1,2,,m.
6:     nit – int64int32nag_int scalar
The number of iterations performed.
7:     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:

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

  ifail = 1ifail=1
On entry,n1n1,
orm < 1m<1,
orn < mn<m,
orldx < nldx<n.
  ifail = 2ifail=2
On entry,tol0.0tol0.0,
ormaxit0maxit0,
ordiagonal element of a = 0.0a=0.0,
orbl0.0bl0.0,
orbd0.0bd0.0.
  ifail = 3ifail=3
A column of x has a constant value.
  ifail = 4ifail=4
Value of u or w returned by ucv < 0ucv<0.
  ifail = 5ifail=5
The function has failed to converge in maxit iterations.
W ifail = 6ifail=6
One of the following is zero: D1D1, D2D2 or D3D3.
This may be caused by the functions uu or ww being too strict for the current estimate of AA (or CC). You should try either a larger initial estimate of AA or make uu and ww 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 = d1= the maximum value of |sjl||sjl|
(ii) d2 = d2= the maximum absolute change in wt(i)wt(i)
(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 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_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
 

user =

     4     2


covar =

    3.2778
   -3.6918
    5.2841
    4.7391
   -6.4086
   11.8371


aOut =

    0.5523
    1.0614
    0.9424
   -0.1881
    0.4776
    0.5021


wt =

    1.0000
    1.0000
    1.0000
    1.0000
    0.2339
    1.0000
    1.0000
    0.9385
    0.4012
    0.7579


thetaOut =

    5.6998
    3.8636
   14.7036


nit =

                   25


ifail =

                    0


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
 

user =

     4     2


covar =

    3.2778
   -3.6918
    5.2841
    4.7391
   -6.4086
   11.8371


aOut =

    0.5523
    1.0614
    0.9424
   -0.1881
    0.4776
    0.5021


wt =

    1.0000
    1.0000
    1.0000
    1.0000
    0.2339
    1.0000
    1.0000
    0.9385
    0.4012
    0.7579


thetaOut =

    5.6998
    3.8636
   14.7036


nit =

                   25


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