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_univar_robust_1var_mestim_wgt (g07dc)

 Contents

    1  Purpose
    2  Syntax
    7  Accuracy
    9  Example

Purpose

nag_univar_robust_1var_mestim_wgt (g07dc) computes an M-estimate of location with (optional) simultaneous estimation of scale, where you provide the weight functions.

Syntax

[theta, sigma, rs, nit, wrk, ifail] = g07dc(chi, psi, isigma, x, beta, theta, sigma, tol, 'n', n, 'maxit', maxit)
[theta, sigma, rs, nit, wrk, ifail] = nag_univar_robust_1var_mestim_wgt(chi, psi, isigma, x, beta, theta, sigma, tol, 'n', n, 'maxit', maxit)

Description

The data consists of a sample of size n, denoted by x1,x2,,xn, drawn from a random variable X.
The xi are assumed to be independent with an unknown distribution function of the form,
Fxi-θ/σ  
where θ is a location argument, and σ is a scale argument. M-estimators of θ and σ are given by the solution to the following system of equations;
i=1nψxi-θ^/σ^ = 0 i=1nχxi-θ^/σ^ = n-1β  
where ψ and χ are user-supplied weight functions, and β is a constant. Optionally the second equation can be omitted and the first equation is solved for θ^ using an assigned value of σ=σc.
The constant β should be chosen so that σ^ is an unbiased estimator when xi, for i=1,2,,n has a Normal distribution. To achieve this the value of β is calculated as:
β=Eχ=-χz12πexp-z22dz 
The values of ψ xi-θ^σ^ σ^ are known as the Winsorized residuals.
The equations are solved by a simple iterative procedure, suggested by Huber:
σ^k=1βn-1 i=1nχ xi-θ^k-1σ^k-1 σ^k-12  
and
θ^k=θ^k- 1+1n i= 1nψ xi-θ^k- 1σ^k σ^k  
or
σ^k=σc  
if σ is fixed.
The initial values for θ^ and σ^ may be user-supplied or calculated within nag_univar_robust_1var_mestim (g07db) as the sample median and an estimate of σ based on the median absolute deviation respectively.
nag_univar_robust_1var_mestim_wgt (g07dc) is based upon function LYHALG within the ROBETH library, 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 estimation of location and scale in ROBETH Cah. Rech. Doc. IUMSP, No. 3 ROB 1 Institut Universitaire de Médecine Sociale et Préventive, Lausanne

Parameters

Compulsory Input Parameters

1:     chi – function handle or string containing name of m-file
chi must return the value of the weight function χ for a given value of its argument. The value of χ must be non-negative.
[result] = chi(t)

Input Parameters

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

Output Parameters

1:     result – double scalar
The value of the weight function ϕ evaluated at t.
2:     psi – function handle or string containing name of m-file
psi must return the value of the weight 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 value of the weight function ψ evaluated at t.
3:     isigma int64int32nag_int scalar
The value assigned to isigma determines whether σ^ is to be simultaneously estimated.
isigma=0
The estimation of σ^ is bypassed and sigma is set equal to σc.
isigma=1
σ^ is estimated simultaneously.
4:     xn – double array
The vector of observations, x1,x2,,xn.
5:     beta – double scalar
The value of the constant β of the chosen chi function.
Constraint: beta>0.0.
6:     theta – double scalar
If sigma>0, then theta must be set to the required starting value of the estimate of the location argument θ^. A reasonable initial value for θ^ will often be the sample mean or median.
7:     sigma – double scalar
The role of sigma depends on the value assigned to isigma as follows.
If isigma=1, sigma must be assigned a value which determines the values of the starting points for the calculation of θ^ and σ^. If sigma0.0, then nag_univar_robust_1var_mestim_wgt (g07dc) will determine the starting points of θ^ and σ^. Otherwise, the value assigned to sigma will be taken as the starting point for σ^, and theta must be assigned a relevant value before entry, see above.
If isigma=0, sigma must be assigned a value which determines the values of σc, which is held fixed during the iterations, and the starting value for the calculation of θ^. If sigma0, then nag_univar_robust_1var_mestim_wgt (g07dc) will determine the value of σc as the median absolute deviation adjusted to reduce bias (see nag_univar_robust_1var_median (g07da)) and the starting point for θ. Otherwise, the value assigned to sigma will be taken as the value of σc and theta must be assigned a relevant value before entry, see above.
8:     tol – double scalar
The relative precision for the final estimates. Convergence is assumed when the increments for theta, and sigma are less than tol×max1.0,σk-1.
Constraint: tol>0.0.

Optional Input Parameters

1:     n int64int32nag_int scalar
Default: the dimension of the array x.
n, the number of observations.
Constraint: n>1.
2:     maxit int64int32nag_int scalar
Default: 50
The maximum number of iterations that should be used during the estimation.
Constraint: maxit>0.

Output Parameters

1:     theta – double scalar
The M-estimate of the location argument θ^.
2:     sigma – double scalar
The M-estimate of the scale argument σ^, if isigma was assigned the value 1 on entry, otherwise sigma will contain the initial fixed value σc.
3:     rsn – double array
The Winsorized residuals.
4:     nit int64int32nag_int scalar
The number of iterations that were used during the estimation.
5:     wrkn – double array
If sigma0.0 on entry, wrk will contain the n observations in ascending order.
6:     ifail int64int32nag_int scalar
ifail=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
On entry,n1,
ormaxit0,
ortol0.0,
orisigma0 or 1.
   ifail=2
On entry,beta0.0.
   ifail=3
On entry,all elements of the input array x are equal.
   ifail=4
sigma, the current estimate of σ, is zero or negative. This error exit is very unlikely, although it may be caused by too large an initial value of sigma.
   ifail=5
The number of iterations required exceeds maxit.
   ifail=6
On completion of the iterations, the Winsorized residuals were all zero. This may occur when using the isigma=0 option with a redescending ψ function, i.e., ψ=0 if t>τ, for some positive constant τ.
If the given value of σ is too small, then the standardized residuals xi-θ^kσc , will be large and all the residuals may fall into the region for which ψt=0. This may incorrectly terminate the iterations thus making theta and sigma invalid.
Re-enter the function with a larger value of σc or with isigma=1.
   ifail=7
The value returned by the chi function is negative.
   ifail=-99
An unexpected error has been triggered by this routine. Please contact NAG.
   ifail=-399
Your licence key may have expired or may not have been installed correctly.
   ifail=-999
Dynamic memory allocation failed.

Accuracy

On successful exit the accuracy of the results is related to the value of tol, see Arguments.

Further Comments

Standard forms of the functions ψ and χ are given in Hampel et al. (1986), Huber (1981) and Marazzi (1987). nag_univar_robust_1var_mestim (g07db) calculates M-estimates using some standard forms for ψ and χ.
When you supply the initial values, care has to be taken over the choice of the initial value of σ. If too small a value is chosen then initial values of the standardized residuals xi-θ^kσ  will be large. If the redescending ψ functions are used, i.e., ψ=0 if t>τ, for some positive constant τ, then these large values are Winsorized as zero. If a sufficient number of the residuals fall into this category then a false solution may be returned, see page 152 of Hampel et al. (1986).

Example

The following program reads in a set of data consisting of eleven observations of a variable X.
The psi and chi functions used are Hampel's Piecewise Linear Function and Hubers chi function respectively.
Using the following starting values various estimates of θ and σ are calculated and printed along with the number of iterations used:
(a) nag_univar_robust_1var_mestim_wgt (g07dc) determined the starting values, σ is estimated simultaneously.
(b) You must supply the starting values, σ is estimated simultaneously.
(c) nag_univar_robust_1var_mestim_wgt (g07dc) determined the starting values, σ is fixed.
(d) You must supply the starting values, σ is fixed.
function g07dc_example


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

global dchi h1 h2 h3;

dchi = 1.5;
h1   = 1.5;
h2   = 3.0;
h3   = 4.5;

x = [13; 11; 16;  5;  3; 18;  9;  8;  6; 27;  7];

% Controll parameter
beta = 0.3892326;
tol  = 0.0001;

% Loop over combinations of isigma sigma and theta
isigma = int64([ 1  1  0  0]);
sigma  =         [-1  7 -1  7];
theta  =         [ 0  2  0  2];

fprintf('           Input parameters     Output parameters\n');
fprintf(' isigma   sigma   theta   tol    sigma  theta\n');

for j = 1:numel(theta)

  fprintf('%3d   %8.4f%8.4f%8.4f', isigma(j), sigma(j), theta(j), tol);

  [thetaOut, sigmaOut, rs, nit, wrk, ifail] = ...
  g07dc( ...
         @chi, @psi, isigma(j), x, beta, theta(j), sigma(j), tol);

  fprintf(' %8.4f%8.4f\n', sigmaOut, thetaOut);

end



function [result] = chi(t)
  % Hubers Chi function
  global dchi;

  ps = min(dchi, abs(t));
  result = ps*ps/2;

function [result] = psi(t)
  % Hampels piecewise linear function
  global h1 h2 h3;

  if abs(t) < h3
    if abs(t) < h2
      result=min(h1, abs(t));
    else
      result=h1*(h3-abs(t))/(h3-h2);
    end
    if t < 0
      result = -result;
    end
  else
    result=0;
  end
g07dc example results

           Input parameters     Output parameters
 isigma   sigma   theta   tol    sigma  theta
  1    -1.0000  0.0000  0.0001   6.3247 10.5487
  1     7.0000  2.0000  0.0001   6.3249 10.5487
  0    -1.0000  0.0000  0.0001   5.9304 10.4896
  0     7.0000  2.0000  0.0001   7.0000 10.6500

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–2015