Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

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 θ$\theta$ be the regression parameters and let C$C$ be the asymptotic variance-covariance matrix of θ̂$\stackrel{^}{\theta }$. Then for Huber type regression
 C = fH(XTX) − 1σ̂2, $C=fH(XTX)-1σ^2,$
where
 fH = 1/(n − m) ( ∑ 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/σ^)-1n∑i=1nψ′ (ri/σ^) ) 2 (1n ∑i=1nψ′ (riσ^) ) 2 ,$
see Huber (1981) and Marazzi (1987).
For Mallows and Schweppe type regressions, C$C$ is of the form
 (σ̂)/n2S1 − 1S2S1 − 1, $σ^n2S1-1S2S1-1,$
where S1 = 1/nXTDX${S}_{1}=\frac{1}{n}{X}^{\mathrm{T}}DX$ and S2 = 1/nXTPX${S}_{2}=\frac{1}{n}{X}^{\mathrm{T}}PX$.
D$D$ is a diagonal matrix such that the i$i$th element approximates E(ψ(ri / (σwi)))$E\left({\psi }^{\prime }\left({r}_{i}/\left(\sigma {w}_{i}\right)\right)\right)$ in the Schweppe case and E(ψ(ri / σ)wi)$E\left({\psi }^{\prime }\left({r}_{i}/\sigma \right){w}_{i}\right)$ in the Mallows case.
P$P$ is a diagonal matrix such that the i$i$th element approximates E(ψ2(ri / (σwi))wi2)$E\left({\psi }^{2}\left({r}_{i}/\left(\sigma {w}_{i}\right)\right){w}_{i}^{2}\right)$ in the Schweppe case and E(ψ2(ri / σ)wi2)$E\left({\psi }^{2}\left({r}_{i}/\sigma \right){w}_{i}^{2}\right)$ in the Mallows case.
Two approximations are available in nag_correg_robustm_user_varmat (g02hf):
1. Average over the ri${r}_{i}$
$Schweppe Mallows Di=(1n∑j=1nψ′ (rjσ^wi ) ) wi Di=(1n∑j=1nψ′ (rjσ^) ) wi Pi=(1n∑j=1nψ2 (rjσ^wi ) ) wi2 Pi=(1n∑j=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$
In all cases σ̂$\stackrel{^}{\sigma }$ is a robust estimate of σ$\sigma$.
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 ψ$\psi$ 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)${\psi }^{\prime }\left(t\right)=\frac{d}{dt}\psi \left(t\right)$ 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 < 0${\mathbf{indw}}<0$
Mallows type regression.
indw = 0${\mathbf{indw}}=0$
Huber type regression.
indw > 0${\mathbf{indw}}>0$
Schweppe type regression.
4:     indc – int64int32nag_int scalar
If indw0${\mathbf{indw}}\ne 0$, indc must specify the approximation to be used.
If indc = 1${\mathbf{indc}}=1$, averaging over residuals.
If indc1${\mathbf{indc}}\ne 1$ , replacing expected by observed.
If indw = 0${\mathbf{indw}}=0$, indc is not referenced.
5:     sigma – double scalar
The value of σ̂$\stackrel{^}{\sigma }$, as given by nag_correg_robustm_user (g02hd).
Constraint: sigma > 0.0${\mathbf{sigma}}>0.0$.
6:     x(ldx,m) – double array
ldx, the first dimension of the array, must satisfy the constraint ldxn$\mathit{ldx}\ge {\mathbf{n}}$.
The values of the X$X$ matrix, i.e., the independent variables. x(i,j)${\mathbf{x}}\left(\mathit{i},\mathit{j}\right)$ must contain the ij$\mathit{i}\mathit{j}$th element of X$X$, for i = 1,2,,n$\mathit{i}=1,2,\dots ,n$ and j = 1,2,,m$\mathit{j}=1,2,\dots ,m$.
7:     rs(n) – double array
n, the dimension of the array, must satisfy the constraint n > 1${\mathbf{n}}>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 > 1${\mathbf{n}}>1$.
If indw0${\mathbf{indw}}\ne 0$, 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 = 0${\mathbf{indw}}=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.)
n$n$, the number of observations.
Constraint: n > 1${\mathbf{n}}>1$.
2:     m – int64int32nag_int scalar
Default: The second dimension of the array x.
m$m$, the number of independent variables.
Constraint: 1m < n$1\le {\mathbf{m}}<{\mathbf{n}}$.

ldx ldc

### Output Parameters

1:     c(ldc,m) – double array
ldcm$\mathit{ldc}\ge {\mathbf{m}}$.
The estimate of the variance-covariance matrix.
2:     wk(m × (n + m + 1) + 2 × n${\mathbf{m}}×\left({\mathbf{n}}+{\mathbf{m}}+1\right)+2×{\mathbf{n}}$) – double array
If indw0${\mathbf{indw}}\ne 0$, wk(i)${\mathbf{wk}}\left(\mathit{i}\right)$, for i = 1,2,,n$\mathit{i}=1,2,\dots ,n$, will contain the diagonal elements of the matrix D$D$ and wk(i)${\mathbf{wk}}\left(\mathit{i}\right)$, for i = n + 1,,2n$\mathit{i}=n+1,\dots ,2n$, will contain the diagonal elements of matrix P$P$.
The rest of the array is used as workspace.
3:     ifail – int64int32nag_int scalar
${\mathrm{ifail}}={\mathbf{0}}$ unless the function detects an error (see [Error Indicators and Warnings]).

## Error Indicators and Warnings

Errors or warnings detected by the function:
ifail = 1${\mathbf{ifail}}=1$
 On entry, n ≤ 1${\mathbf{n}}\le 1$, or m < 1${\mathbf{m}}<1$, or n ≤ m${\mathbf{n}}\le {\mathbf{m}}$, or ldc < m$\mathit{ldc}<{\mathbf{m}}$, or ldx < n$\mathit{ldx}<{\mathbf{n}}$.
ifail = 2${\mathbf{ifail}}=2$
 On entry, sigma ≤ 0.0${\mathbf{sigma}}\le 0.0$.
ifail = 3${\mathbf{ifail}}=3$
If indw = 0${\mathbf{indw}}=0$ then the matrix XTX${X}^{\mathrm{T}}X$ is either not positive definite, possibly due to rounding errors, or is ill-conditioned.
If indw0${\mathbf{indw}}\ne 0$ then the matrix S1${S}_{1}$ is singular or almost singular. This may be due to many elements of D$D$ being zero.
ifail = 4${\mathbf{ifail}}=4$
Either the value of 1/ni = 1nψ ((ri)/(σ̂)) = 0$\frac{1}{n}\sum _{i=1}^{n}{\psi }^{\prime }\left(\frac{{r}_{i}}{\stackrel{^}{\sigma }}\right)=0$,
or κ = 0$\kappa =0$,
or i = 1nψ2 ((ri)/(σ̂)) = 0$\sum _{i=1}^{n}{\psi }^{2}\left(\frac{{r}_{i}}{\stackrel{^}{\sigma }}\right)=0$.
In this situation nag_correg_robustm_user_varmat (g02hf) returns C$C$ as (XTX)1${\left({X}^{\mathrm{T}}X\right)}^{-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).

nag_correg_robustm_user_varmat (g02hf) is only for situations in which X$X$ has full column rank.
Care has to be taken in the choice of the ψ$\psi$ function since if ψ(t) = 0${\psi }^{\prime }\left(t\right)=0$ for too wide a range then either the value of fH${f}_{H}$ will not exist or too many values of Di${D}_{i}$ will be zero and it will not be possible to calculate C$C$.

## 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

```