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_user_varmat (g02hf)

Purpose

nag_correg_robustm_user_varmat (g02hf) calculates an estimate of the asymptotic variance-covariance matrix for the bounded influence regression estimates (M-estimates). It is intended for use with nag_correg_robustm_user (g02hd).

Syntax

[c, wk, ifail] = g02hf(psi, psp, indw, indc, sigma, x, rs, wgt, 'n', n, 'm', m)
[c, wk, ifail] = nag_correg_robustm_user_varmat(psi, psp, indw, indc, sigma, x, rs, wgt, 'n', n, 'm', m)

Description

For a description of bounded influence regression see nag_correg_robustm_user (g02hd). Let θθ be the regression parameters and let CC be the asymptotic variance-covariance matrix of θ̂θ^. Then for Huber type regression
C = fH(XTX)1σ̂2,
C=fH(XTX)-1σ^2,
where
fH = 1/(nm) (i = 1nψ2 (ri / σ̂) )/((1/nψ((ri)/(σ̂)))2) κ2
fH=1n-m i= 1nψ2 (ri/σ^) (1nψ (riσ^) ) 2 κ2
κ2 = 1 + m/n (1/n ∑ i = 1n
( n )ψ(ri / σ̂) − 1/n ∑ ψ(ri / σ̂) i = 1 2
)/(
( n )1/n ∑ ψ((ri)/(σ̂)) i = 1 2
),
κ2=1+mn 1n i=1n (ψ (ri/σ^)-1ni=1nψ (ri/σ^) ) 2 (1n i=1nψ (riσ^) ) 2 ,
see Huber (1981) and Marazzi (1987).
For Mallows and Schweppe type regressions, CC is of the form
(σ̂)/n2S11S2S11,
σ^n2S1-1S2S1-1,
where S1 = 1/nXTDXS1=1nXTDX and S2 = 1/nXTPXS2=1nXTPX.
DD is a diagonal matrix such that the iith element approximates E(ψ(ri / (σwi)))E(ψ(ri/(σwi))) in the Schweppe case and E(ψ(ri / σ)wi)E(ψ(ri/σ)wi) in the Mallows case.
PP is a diagonal matrix such that the iith element approximates E(ψ2(ri / (σwi))wi2)E(ψ2(ri/(σwi))wi2) in the Schweppe case and E(ψ2(ri / σ)wi2)E(ψ2(ri/σ)wi2) in the Mallows case.
Two approximations are available in nag_correg_robustm_user_varmat (g02hf):
  1. Average over the riri 
    Schweppe Mallows Di=(1nj=1nψ (rjσ^wi ) ) wi Di=(1nj=1nψ (rjσ^) ) wi Pi=(1nj=1nψ2 (rjσ^wi ) ) wi2 Pi=(1nj=1nψ2 (rjσ^) ) wi2
    Schweppe Mallows
    Di = ( n )1/n ∑ ψ((rj)/(σ̂wi)) j = 1 wi 
    Di = ( n )1/n ∑ ψ((rj)/(σ̂)) j = 1 wi
    Pi = ( n )1/n ∑ ψ2((rj)/(σ̂wi)) j = 1 wi2
    Pi = ( n )1/n ∑ ψ2((rj)/(σ̂)) j = 1 wi2
  2. Replace expected value by observed
    Schweppe Mallows
    Di = ψ ((ri)/(σ̂wi)) wi Di = ψ ((ri)/(σ̂)) wi
    Pi = ψ2 ((ri)/(σ̂wi)) wi2 Pi = ψ2 ((ri)/(σ̂)) wi2
    Schweppe Mallows Di=ψ ( riσ ^wi ) wi Di=ψ ( riσ ^) wi Pi=ψ2 ( riσ ^wi ) wi2 Pi=ψ2 ( riσ ^) wi2
See Hampel et al. (1986) and Marazzi (1987).
In all cases σ̂σ^ is a robust estimate of σσ.
nag_correg_robustm_user_varmat (g02hf) is based on routines in ROBETH; see Marazzi (1987).

References

Hampel F R, Ronchetti E M, Rousseeuw P J and Stahel W A (1986) Robust Statistics. The Approach Based on Influence Functions Wiley
Huber P J (1981) Robust Statistics Wiley
Marazzi A (1987) Subroutines for robust and bounded influence regression in ROBETH Cah. Rech. Doc. IUMSP, No. 3 ROB 2 Institut Universitaire de Médecine Sociale et Préventive, Lausanne

Parameters

Compulsory Input Parameters

1:     psi – function handle or string containing name of m-file
psi must return the value of the ψψ function for a given value of its argument.
[result] = psi(t)

Input Parameters

1:     t – double scalar
The argument for which psi must be evaluated.

Output Parameters

1:     result – double scalar
The result of the function.
2:     psp – function handle or string containing name of m-file
psp must return the value of ψ(t) = d/(dt)ψ(t)ψ(t)=ddt ψ(t) for a given value of its argument.
[result] = psp(t)

Input Parameters

1:     t – double scalar
The argument for which psp must be evaluated.

Output Parameters

1:     result – double scalar
The result of the function.
3:     indw – int64int32nag_int scalar
The type of regression for which the asymptotic variance-covariance matrix is to be calculated.
indw < 0indw<0
Mallows type regression.
indw = 0indw=0
Huber type regression.
indw > 0indw>0
Schweppe type regression.
4:     indc – int64int32nag_int scalar
If indw0indw0, indc must specify the approximation to be used.
If indc = 1indc=1, averaging over residuals.
If indc1indc1 , replacing expected by observed.
If indw = 0indw=0, indc is not referenced.
5:     sigma – double scalar
The value of σ̂σ^, as given by nag_correg_robustm_user (g02hd).
Constraint: sigma > 0.0sigma>0.0.
6:     x(ldx,m) – double array
ldx, the first dimension of the array, must satisfy the constraint ldxnldxn.
The values of the XX matrix, i.e., the independent variables. x(i,j)xij must contain the ijijth element of XX, for i = 1,2,,ni=1,2,,n and j = 1,2,,mj=1,2,,m.
7:     rs(n) – double array
n, the dimension of the array, must satisfy the constraint n > 1n>1.
The residuals from the bounded influence regression. These are given by nag_correg_robustm_user (g02hd).
8:     wgt(n) – double array
n, the dimension of the array, must satisfy the constraint n > 1n>1.
If indw0indw0, wgt must contain the vector of weights used by the bounded influence regression. These should be used with nag_correg_robustm_user (g02hd).
If indw = 0indw=0, wgt is not referenced.

Optional Input Parameters

1:     n – int64int32nag_int scalar
Default: The dimension of the arrays rs, wgt and the first dimension of the array x. (An error is raised if these dimensions are not equal.)
nn, the number of observations.
Constraint: n > 1n>1.
2:     m – int64int32nag_int scalar
Default: The second dimension of the array x.
mm, the number of independent variables.
Constraint: 1m < n1m<n.

Input Parameters Omitted from the MATLAB Interface

ldx ldc

Output Parameters

1:     c(ldc,m) – double array
ldcmldcm.
The estimate of the variance-covariance matrix.
2:     wk(m × (n + m + 1) + 2 × nm×(n+m+1)+2×n) – double array
If indw0indw0, wk(i)wki, for i = 1,2,,ni=1,2,,n, will contain the diagonal elements of the matrix DD and wk(i)wki, for i = n + 1,,2ni=n+1,,2n, will contain the diagonal elements of matrix PP.
The rest of the array is used as workspace.
3:     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:
  ifail = 1ifail=1
On entry,n1n1,
orm < 1m<1,
ornmnm,
orldc < mldc<m,
orldx < nldx<n.
  ifail = 2ifail=2
On entry,sigma0.0sigma0.0.
  ifail = 3ifail=3
If indw = 0indw=0 then the matrix XTXXTX is either not positive definite, possibly due to rounding errors, or is ill-conditioned.
If indw0indw0 then the matrix S1S1 is singular or almost singular. This may be due to many elements of DD being zero.
  ifail = 4ifail=4
Either the value of 1/ni = 1nψ ((ri)/(σ̂)) = 0 1ni=1nψ ( riσ^)=0,
or κ = 0κ=0,
or i = 1nψ2 ((ri)/(σ̂)) = 0i=1nψ2 ( riσ^)=0.
In this situation nag_correg_robustm_user_varmat (g02hf) returns CC as (XTX)1(XTX)-1.

Accuracy

In general, the accuracy of the variance-covariance matrix will depend primarily on the accuracy of the results from nag_correg_robustm_user (g02hd).

Further Comments

nag_correg_robustm_user_varmat (g02hf) is only for situations in which XX has full column rank.
Care has to be taken in the choice of the ψψ function since if ψ(t) = 0ψ(t)=0 for too wide a range then either the value of fHfH will not exist or too many values of DiDi will be zero and it will not be possible to calculate CC.

Example

function nag_correg_robustm_user_varmat_example
indw = int64(1);
indc = int64(1);
sigma = 20.7783;
x = [1, -1, -1;
     1, -1, 1;
     1, 1, -1;
     1, 1, 1;
     1, 0, 3];
rs = [0.5643;
     -1.1286;
     0.5643;
     -1.1286;
     1.1286];
wgt = [0.4039;
     0.5012;
     0.4039;
     0.5012;
     0.3862];
[c, wk, ifail] = nag_correg_robustm_user_varmat(@psi, @psp, indw, indc, sigma, x, rs, wgt)

function [result] = psi(t)
  if t < -1.5
    result = -1.5;
  elseif abs(t) < 1.5
    result = t;
  else
    result = 1.5;
  end;
function [result] = psp(t)
  if (abs(t) < 1.5)
    result=1;
  else
    result=0;
  end
 

c =

    0.2070         0   -0.0478
         0    0.2229         0
   -0.0478         0    0.0796


wk =

    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    0.0021
    0.0021
    0.0021
    0.0021
    0.0021
    1.1607
         0
   -0.2679
         0
    1.2500
         0
   -0.2679
         0
    0.4464
         0
   -1.0000
    1.0000
    0.0021
    0.0021
         0
    0.0000
         0
    0.0021
         0
    0.0000
         0
    0.0021
   -0.0021
    0.0021
   -0.0021
    0.0021
    0.0062


ifail =

                    0


function g02hf_example
indw = int64(1);
indc = int64(1);
sigma = 20.7783;
x = [1, -1, -1;
     1, -1, 1;
     1, 1, -1;
     1, 1, 1;
     1, 0, 3];
rs = [0.5643;
     -1.1286;
     0.5643;
     -1.1286;
     1.1286];
wgt = [0.4039;
     0.5012;
     0.4039;
     0.5012;
     0.3862];
[c, wk, ifail] = g02hf(@psi, @psp, indw, indc, sigma, x, rs, wgt)

function [result] = psi(t)
  if t < -1.5
    result = -1.5;
  elseif abs(t) < 1.5
    result = t;
  else
    result = 1.5;
  end;
function [result] = psp(t)
  if (abs(t) < 1.5)
    result=1;
  else
    result=0;
  end
 

c =

    0.2070         0   -0.0478
         0    0.2229         0
   -0.0478         0    0.0796


wk =

    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    0.0021
    0.0021
    0.0021
    0.0021
    0.0021
    1.1607
         0
   -0.2679
         0
    1.2500
         0
   -0.2679
         0
    0.4464
         0
   -1.0000
    1.0000
    0.0021
    0.0021
         0
    0.0000
         0
    0.0021
         0
    0.0000
         0
    0.0021
   -0.0021
    0.0021
   -0.0021
    0.0021
    0.0062


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