g07dc:: Univariate Estimation (NAG Toolbox)

The

x_{i}

are assumed to be independent with an unknown distribution function of the form,

F ((x_{i} - θ) / σ)

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;

\begin{array}{lcl} \sum_{i = 1}^{n} ψ ((x_{i} - \hat{θ}) / \hat{σ}) & = & 0 \\ \sum_{i = 1}^{n} χ ((x_{i} - \hat{θ}) / \hat{σ}) & = & (n - 1) β \end{array}

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

\hat{θ}

using an assigned value of

σ = σ_{c}

The constant

β

should be chosen so that

\hat{σ}

is an unbiased estimator when

x_{i}

, for

i = 1, 2, \dots, n

has a Normal distribution. To achieve this the value of

β

is calculated as:

β = E (χ) = \int_{- \infty}^{\infty} χ (z) \frac{1}{\sqrt{2 π}} \exp \{\frac{- z^{2}}{2}\} d z

The values of

ψ (\frac{x_{i} - \hat{θ}}{\hat{σ}}) \hat{σ}

are known as the Winsorized residuals.

The equations are solved by a simple iterative procedure, suggested by Huber:

{\hat{σ}}_{k} = \sqrt{\frac{1}{β (n - 1)} (\sum_{i = 1}^{n} χ (\frac{x_{i} - {\hat{θ}}_{k - 1}}{{\hat{σ}}_{k - 1}})) {\hat{σ}}_{k - 1}^{2}}

and

{\hat{θ}}_{k} = {\hat{θ}}_{k - 1} + \frac{1}{n} \sum_{i = 1}^{n} ψ (\frac{x_{i} - {\hat{θ}}_{k - 1}}{{\hat{σ}}_{k}}) {\hat{σ}}_{k}

{\hat{σ}}_{k} = σ_{c}

σ

is fixed.

The initial values for

\hat{θ}

and

\hat{σ}

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.

References

Parameters

Compulsory Input Parameters

Optional Input Parameters

Output Parameters

Error Indicators and Warnings

Accuracy

Further Comments

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

\frac{x_{i} - {\hat{θ}}_{k}}{σ}

will be large. If the redescending

ψ

functions are used, i.e.,

ψ = 0

|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

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

On entry,	$n \leq 1$ ,
or	$maxit \leq 0$ ,
or	$tol \leq 0.0$ ,
or	$isigma \neq 0$ or $1$ .

NAG Toolbox: nag_univar_robust_1var_mestim_wgt (g07dc)

▸▿ Contents

Purpose

Syntax

Description