NAG Toolbox: nag_correg_robustm_user (g02hd)

Purpose

Syntax

Description

References

Parameters

Compulsory Input Parameters

1:     $\mathrm{chi}$ – function handle or string containing name of m-file

If $\mathbf{isigma}>0$, chi must return the value of the weight function $\chi$ for a given value of its argument. The value of $\chi$ must be non-negative.

[result] = chi(t)

Input Parameters

1:     $\mathrm{t}$ – double scalar

The argument for which chi must be evaluated.

Output Parameters

1:     $\mathrm{result}$ – double scalar

The value of the weight function $\chi$ evaluated at t.

If $\mathbf{isigma}\le 0$, the actual argument chi may be the string nag_correg_robustm_user_dummy_chi (g02hdz). (nag_correg_robustm_user_dummy_chi (g02hdz) is included in the NAG Toolbox.)

2:     $\mathrm{psi}$ – function handle or string containing name of m-file

psi must return the value of the weight function $\psi$ for a given value of its argument.

[result] = psi(t)

Input Parameters

1:     $\mathrm{t}$ – double scalar

The argument for which psi must be evaluated.

Output Parameters

1:     $\mathrm{result}$ – double scalar

The value of the weight function $\psi$ evaluated at t.

3:     $\mathrm{psip0}$ – double scalar

The value of $\psi \prime \left(0\right)$.

4:     $\mathrm{beta}$ – double scalar

If $\mathbf{isigma}<0$, beta must specify the value of ${\beta }_1$.

For Huber and Schweppe type regressions, ${\beta }_1$ is the $75$th percentile of the standard Normal distribution (see nag_stat_inv_cdf_normal (g01fa)). For Mallows type regression ${\beta }_1$ is the solution to

$$\frac{1}{n}\sum _{i=1}^n\Phi \left({\beta }_1/\sqrt{w_i}\right)=0.75\text{,}$$

where $\Phi$ is the standard Normal cumulative distribution function (see nag_specfun_cdf_normal (s15ab)).

If $\mathbf{isigma}>0$, beta must specify the value of ${\beta }_2$.

$$\begin{array}{lll}{\beta }_2=& \underset{-\infty }{\overset{\infty }{\int }}\chi \left(z\right)\varphi \left(z\right)dz\text{,}& \text{in the Huber case;}\\ & & \\ & & \\ {\beta }_2=& \frac{1}{n}\sum _{i=1}^n{w}_i\underset{-\infty }{\overset{\infty }{\int }}\chi \left(z\right)\varphi \left(z\right)dz\text{,}& \text{in the Mallows case;}\\ & & \\ & & \\ {\beta }_2=& \frac{1}{n}\sum _{i=1}^n{w}_i^2\underset{-\infty }{\overset{\infty }{\int }}\chi \left(z/w_i\right)\varphi \left(z\right)dz\text{,}& \text{in the Schweppe case;}\end{array}$$

where $\varphi$ is the standard normal density, i.e., $\frac{1}{\sqrt{2\pi }}\exp \left(-\frac{1}{2}{x}^2\right)$.

If $\mathbf{isigma}=0$, beta is not referenced.

Constraint: if $\mathbf{isigma}\ne 0$, $\mathbf{beta}>0.0$.

5:     $\mathrm{indw}$ – int64int32nag_int scalar

Determines the type of regression to be performed.

$\mathbf{indw}=0$

Huber type regression.

$\mathbf{indw}<0$

Mallows type regression.

$\mathbf{indw}>0$

Schweppe type regression.

6:     $\mathrm{isigma}$ – int64int32nag_int scalar

Determines how $\sigma$ is to be estimated.

$\mathbf{isigma}=0$

$\sigma$ is held constant at its initial value.

$\mathbf{isigma}<0$

$\sigma$ is estimated by median absolute deviation of residuals.

$\mathbf{isigma}>0$

$\sigma$ is estimated using the $\chi$ function.

7:     $\mathrm{x}\left(\mathit{ldx},\mathbf{m}\right)$ – double array

ldx, the first dimension of the array, must satisfy the constraint $\mathit{ldx}\ge \mathbf{n}$.

The values of the $X$ matrix, i.e., the independent variables. $\mathbf{x}\left(\mathit{i},\mathit{j}\right)$ must contain the $\mathit{i}\mathit{j}$th element of $\mathbf{x}$, for $\mathit{i}=1,2,\dots ,n$ and $\mathit{j}=1,2,\dots ,m$.

If $\mathbf{indw}<0$, during calculations the elements of x will be transformed as described in Description. Before exit the inverse transformation will be applied. As a result there may be slight differences between the input x and the output x.

8:     $\mathrm{y}\left(\mathbf{n}\right)$ – double array

The data values of the dependent variable.

$\mathbf{y}\left(\mathit{i}\right)$ must contain the value of $y$ for the $\mathit{i}$th observation, for $\mathit{i}=1,2,\dots ,n$.

If $\mathbf{indw}<0$, during calculations the elements of y will be transformed as described in Description. Before exit the inverse transformation will be applied. As a result there may be slight differences between the input y and the output y.

9:     $\mathrm{wgt}\left(\mathbf{n}\right)$ – double array

The weight for the $\mathit{i}$th observation, for $\mathit{i}=1,2,\dots ,n$.

If $\mathbf{indw}<0$, during calculations elements of wgt will be transformed as described in Description. Before exit the inverse transformation will be applied. As a result there may be slight differences between the input wgt and the output wgt.

If $\mathbf{wgt}\left(i\right)\le 0$, the $i$th observation is not included in the analysis.

If $\mathbf{indw}=0$, wgt is not referenced.

10:   $\mathrm{theta}\left(\mathbf{m}\right)$ – double array

Starting values of the argument vector $\theta$. These may be obtained from least squares regression. Alternatively if $\mathbf{isigma}<0$ and $\mathbf{sigma}=1$ or if $\mathbf{isigma}>0$ and sigma approximately equals the standard deviation of the dependent variable, $y$, then $\mathbf{theta}\left(\mathit{i}\right)=0.0$, for $\mathit{i}=1,2,\dots ,m$ may provide reasonable starting values.

11:   $\mathrm{sigma}$ – double scalar

A starting value for the estimation of $\sigma$. sigma should be approximately the standard deviation of the residuals from the model evaluated at the value of $\theta$ given by theta on entry.

Constraint: $\mathbf{sigma}>0.0$.

Optional Input Parameters

1:     $\mathrm{n}$ – int64int32nag_int scalar

Default: the dimension of the arrays y, wgt and the first dimension of the array x. (An error is raised if these dimensions are not equal.)

$n$, the number of observations.

Constraint: $\mathbf{n}>1$.

2:     $\mathrm{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 independent variables.

Constraint: $1\le \mathbf{m}<\mathbf{n}$.

3:     $\mathrm{tol}$ – double scalar

Default: $5e-5$

The relative precision for the final estimates. Convergence is assumed when both the relative change in the value of sigma and the relative change in the value of each element of theta are less than tol.

It is advisable for tol to be greater than $100\times \mathit{machine precision}$.

Constraint: $\mathbf{tol}>0.0$.

4:     $\mathrm{eps}$ – double scalar

Default: $5e-6$

A relative tolerance to be used to determine the rank of $X$. See nag_linsys_real_gen_solve (f04jg) for further details.

If $\mathbf{eps}<\mathit{machine precision}$ or $\mathbf{eps}>1.0$ then machine precision will be used in place of tol.

A reasonable value for eps is $5.0\times {10}^{-6}$ where this value is possible.

5:     $\mathrm{maxit}$ – int64int32nag_int scalar

Default: $50$

The maximum number of iterations that should be used during the estimation.

A value of $\mathbf{maxit}=50$ should be adequate for most uses.

Constraint: $\mathbf{maxit}>0$.

6:     $\mathrm{nitmon}$ – int64int32nag_int scalar

Default: $0$

Determines the amount of information that is printed on each iteration.

$\mathbf{nitmon}\le 0$

No information is printed.

$\mathbf{nitmon}>0$

On the first and every nitmon iterations the values of sigma, theta and the change in theta during the iteration are printed.

When printing occurs the output is directed to the current advisory message unit (see nag_file_set_unit_advisory (x04ab)).

Output Parameters

1:     $\mathrm{x}\left(\mathit{ldx},\mathbf{m}\right)$ – double array

Unchanged, except as described above.

2:     $\mathrm{y}\left(\mathbf{n}\right)$ – double array

Unchanged, except as described above.

3:     $\mathrm{wgt}\left(\mathbf{n}\right)$ – double array

Unchanged, except as described above.

4:     $\mathrm{theta}\left(\mathbf{m}\right)$ – double array

The M-estimate of ${\theta }_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,m$.

5:     $\mathrm{k}$ – int64int32nag_int scalar

The column rank of the matrix $X$.

6:     $\mathrm{sigma}$ – double scalar

The final estimate of $\sigma$ if $\mathbf{isigma}\ne 0$ or the value assigned on entry if $\mathbf{isigma}=0$.

7:     $\mathrm{rs}\left(\mathbf{n}\right)$ – double array

The residuals from the model evaluated at final value of theta, i.e., rs contains the vector $\left(y-X\hat{\theta }\right)$.

8:     $\mathrm{nit}$ – int64int32nag_int scalar

The number of iterations that were used during the estimation.

9:     $\mathrm{ifail}$ – int64int32nag_int scalar

$\mathbf{ifail}=\mathbf{0}$ unless the function detects an error (see Error Indicators and Warnings).

Error Indicators and Warnings

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

   $\mathbf{ifail}=1$

On entry,	$\mathbf{n}\le 1$,
or	$\mathbf{m}<1$,
or	$\mathbf{n}\le \mathbf{m}$,
or	$\mathit{ldx}<\mathbf{n}$.

   $\mathbf{ifail}=2$

On entry,	$\mathbf{beta}\le 0.0$, and $\mathbf{isigma}\ne 0$,
or	$\mathbf{sigma}\le 0.0$.

   $\mathbf{ifail}=3$

On entry,	$\mathbf{tol}\le 0.0$,
or	$\mathbf{maxit}\le 0$.

   $\mathbf{ifail}=4$

A value returned by the chi function is negative.

   $\mathbf{ifail}=5$

During iterations a value of $\mathbf{sigma}\le 0.0$ was encountered.

   $\mathbf{ifail}=6$

A failure occurred in nag_linsys_real_gen_solve (f04jg) . This is an extremely unlikely error. If it occurs, please contact NAG.

W  $\mathbf{ifail}=7$

The weighted least squares equations are not of full rank. This may be due to the $X$ matrix not being of full rank, in which case the results will be valid. It may also occur if some of the $G_{ii}$ values become very small or zero, see Further Comments. The rank of the equations is given by k. If the matrix just fails the test for nonsingularity then the result $\mathbf{ifail}=\mathbf{7}$ and $\mathbf{k}=\mathbf{m}$ is possible (see nag_linsys_real_gen_solve (f04jg)).

   $\mathbf{ifail}=8$

The function has failed to converge in maxit iterations.

   $\mathbf{ifail}=9$

Having removed cases with zero weight, the value of $\mathbf{n}-\mathbf{k}\le 0$, i.e., no degree of freedom for error. This error will only occur if $\mathbf{isigma}>0$.

   $\mathbf{ifail}=-99$

An unexpected error has been triggered by this routine. Please contact NAG.

   $\mathbf{ifail}=-399$

Your licence key may have expired or may not have been installed correctly.

   $\mathbf{ifail}=-999$

Dynamic memory allocation failed.

Accuracy

Further Comments

Example

Open in the MATLAB editor: g02hd_example

function g02hd_example


fprintf('g02hd example results\n\n');

global dchi
dchi = 1.5;

x = [1, -1, -1;
     1, -1,  1;
     1,  1, -1;
     1,  1,  1;
     1,  0,  3];
y   = [10.5;    11.3;     12.6;     13.4;     17.1   ];
wgt = [0.4039;   0.5012;   0.4039;   0.5012;   0.3862];

[n,m] = size(x);

% Calculate beta
[beta] = betcal(wgt);

% Control Parameters
psip0 = 1;
indw = int64(1);
isigma = int64(1);

% Intial values
sigma = 1;
theta = zeros(m,1);

% Perform bounded influence regression
[x, y, wgt, theta, k, sigma, rs, nit, ifail] = ...
  g02hd( ...
         @chi, @psi, psip0, beta, indw, isigma, x, y, wgt, theta, sigma);

fprintf(' iterations to convergence = %4d\n', nit);
fprintf('                         k = %4d\n',k);
fprintf('                     sigma = %9.4f\n',sigma);
fprintf('Theta:\n');
disp(theta');
fprintf('\n  Weights  Residuals\n');
fprintf('%9.4f%9.4f\n', [wgt rs]');



function [result] = chi(t)
  global dchi

  if (abs(t) < dchi)
    ps=t;
  else
    ps=dchi;
  end
  result = ps*ps/2;


function [result] = psi(t)
  global dchi

  if t < -dchi
    result = -dchi;
  elseif abs(t) < dchi
    result = t;
  else
    result = dchi;
  end;

function [beta] = betcal(wgt)
  %  Calculate beta for Schweppe type regression
  global dchi

  n = numel(wgt);
  amaxex = -log(x02ak);
  anormc = sqrt(2*pi);
  d2 = dchi*dchi;
  beta = 0;
  for i = 1:n
    w2 = wgt(i)*wgt(i);
    dw = wgt(i)*dchi;
    [pc, ifail] = s15ab(dw);
    dw2 = dw*dw;
    if dw2<amaxex
      dc = exp(-dw2/2)/anormc;
    else
      dc = 0;
    end
    b = (-dw*dc+pc-0.5)/w2 + (1-pc)*d2;
    beta = b*w2/n + beta;
  end

g02hd example results

 iterations to convergence =    5
                         k =    3
                     sigma =    2.7783
Theta:
   12.2321    1.0500    1.2464


  Weights  Residuals
   0.4039   0.5643
   0.5012  -1.1286
   0.4039   0.5643
   0.5012  -1.1286
   0.3862   1.1286