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)

Purpose

nag_univar_robust_1var_mestim_wgt (g07dc) computes an MM-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 nn, denoted by x1,x2,,xnx1,x2,,xn, drawn from a random variable XX.
The xixi are assumed to be independent with an unknown distribution function of the form,
F((xiθ) / σ)
F((xi-θ)/σ)
where θθ is a location parameter, and σσ is a scale parameter. MM-estimators of θθ and σσ are given by the solution to the following system of equations;
n
ψ((xiθ̂) / σ̂)
i = 1
= 0
n
χ((xiθ̂) / σ̂)
i = 1
= (n1)β
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σ=σc.
The constant ββ should be chosen so that σ̂σ^ is an unbiased estimator when xixi, for i = 1,2,,ni=1,2,,n has a Normal distribution. To achieve this the value of ββ is calculated as:
β = E(χ) = χ(z)1/(sqrt(2π))exp{(z2)/2}dz
β=E(χ)=-χ(z)12πexp{-z22}dz
The values of ψ ((xiθ̂)/(σ̂)) σ̂ψ ( xi-θ^σ^) σ^ are known as the Winsorized residuals.
The equations are solved by a simple iterative procedure, suggested by Huber:
σ̂k =
sqrt(1/(β(n − 1))(n ) ∑ χ((xi − θ̂k − 1)/(σ̂k − 1))i = 1 σ̂k − 12)
σ^k=1β(n-1) (i=1nχ (xi-θ^k-1σ^k-1) ) σ^k-12
and
n
θ̂k = θ̂k 1 + 1/nψ((xiθ̂k 1)/(σ̂k))σ̂k
i = 1
θ^k=θ^k- 1+1n i= 1nψ (xi-θ^k- 1σ^k) σ^k
or
σ̂k = σc
σ^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 result of the function.
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 result of the function.
3:     isigma – int64int32nag_int scalar
The value assigned to isigma determines whether σ̂σ^ is to be simultaneously estimated.
isigma = 0isigma=0
The estimation of σ̂σ^ is bypassed and sigma is set equal to σcσc.
isigma = 1isigma=1
σ̂σ^ is estimated simultaneously.
4:     x(n) – double array
n, the dimension of the array, must satisfy the constraint n > 1n>1.
The vector of observations, x1,x2,,xnx1,x2,,xn.
5:     beta – double scalar
The value of the constant ββ of the chosen chi function.
Constraint: beta > 0.0beta>0.0.
6:     theta – double scalar
If sigma > 0sigma>0, then theta must be set to the required starting value of the estimate of the location parameter θ̂θ^. 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 = 1isigma=1, sigma must be assigned a value which determines the values of the starting points for the calculation of θ̂θ^ and σ̂σ^. If sigma0.0sigma0.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 = 0isigma=0, sigma must be assigned a value which determines the values of σcσc, which is held fixed during the iterations, and the starting value for the calculation of θ̂θ^. If sigma0sigma0, then nag_univar_robust_1var_mestim_wgt (g07dc) will determine the value of σcσ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σ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 × max (1.0,σk1)tol×max(1.0,σk-1).
Constraint: tol > 0.0tol>0.0.

Optional Input Parameters

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

Input Parameters Omitted from the MATLAB Interface

None.

Output Parameters

1:     theta – double scalar
The MM-estimate of the location parameter θ̂θ^.
2:     sigma – double scalar
The MM-estimate of the scale parameter σ̂σ^, if isigma was assigned the value 11 on entry, otherwise sigma will contain the initial fixed value σcσc.
3:     rs(n) – double array
The Winsorized residuals.
4:     nit – int64int32nag_int scalar
The number of iterations that were used during the estimation.
5:     wrk(n) – double array
If sigma0.0sigma0.0 on entry, wrk will contain the nn observations in ascending order.
6:     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,
ormaxit0maxit0,
ortol0.0tol0.0,
orisigma0isigma0 or 11.
  ifail = 2ifail=2
On entry,beta0.0beta0.0.
  ifail = 3ifail=3
On entry,all elements of the input array x are equal.
  ifail = 4ifail=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 = 5ifail=5
The number of iterations required exceeds maxit.
  ifail = 6ifail=6
On completion of the iterations, the Winsorized residuals were all zero. This may occur when using the isigma = 0isigma=0 option with a redescending ψψ function, i.e., ψ = 0ψ=0 if |t| > τ|t|>τ, for some positive constant ττ.
If the given value of σσ is too small, then the standardized residuals (xiθ̂k)/(σc) xi-θ^kσc , will be large and all the residuals may fall into the region for which ψ(t) = 0ψ(t)=0. This may incorrectly terminate the iterations thus making theta and sigma invalid.
Re-enter the function with a larger value of σcσc or with isigma = 1isigma=1.
  ifail = 7ifail=7
The value returned by the chi function is negative.

Accuracy

On successful exit the accuracy of the results is related to the value of tol, see Section [Parameters].

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 MM-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)/σ xi-θ^kσ  will be large. If the redescending ψψ functions are used, i.e., ψ = 0ψ=0 if |t| > τ|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

function nag_univar_robust_1var_mestim_wgt_example
isigma = int64(1);
x = [13;
     11;
     16;
     5;
     3;
     18;
     9;
     8;
     6;
     27;
     7];
beta = 0.3892326;
theta = 0;
sigma = -1;
tol = 0.0001;
[thetaOut, sigmaOut, rs, nit, wrk, ifail] = ...
    nag_univar_robust_1var_mestim_wgt(@chi, @psi, isigma, x, beta, theta, sigma, tol)

function [result] = chi(t)
  ps = min(1.5, abs(t));
  result = ps*ps/2;
function [result] = psi(t)

  if abs(t) < 4.5
    if abs(t) < 3
      result=min(1.5, abs(t));
    else
      result=1.5*(4.5-abs(t))/1.5;
    end
    if t < 0
      result = -result;
    end
  else
    result=0;
  end
 

thetaOut =

   10.5487


sigmaOut =

    6.3247


rs =

    2.4513
    0.4513
    5.4513
   -5.5487
   -7.5487
    7.4513
   -1.5487
   -2.5487
   -4.5487
   16.4513
   -3.5487


nit =

                    8


wrk =

     3
     5
     6
     7
     8
     9
    11
    13
    16
    18
    27


ifail =

                    0


function g07dc_example
isigma = int64(1);
x = [13;
     11;
     16;
     5;
     3;
     18;
     9;
     8;
     6;
     27;
     7];
beta = 0.3892326;
theta = 0;
sigma = -1;
tol = 0.0001;
[thetaOut, sigmaOut, rs, nit, wrk, ifail] = ...
    g07dc(@chi, @psi, isigma, x, beta, theta, sigma, tol)

function [result] = chi(t)
  ps = min(1.5, abs(t));
  result = ps*ps/2;
function [result] = psi(t)

  if abs(t) < 4.5
    if abs(t) < 3
      result=min(1.5, abs(t));
    else
      result=1.5*(4.5-abs(t))/1.5;
    end
    if t < 0
      result = -result;
    end
  else
    result=0;
  end
 

thetaOut =

   10.5487


sigmaOut =

    6.3247


rs =

    2.4513
    0.4513
    5.4513
   -5.5487
   -7.5487
    7.4513
   -1.5487
   -2.5487
   -4.5487
   16.4513
   -3.5487


nit =

                    8


wrk =

     3
     5
     6
     7
     8
     9
    11
    13
    16
    18
    27


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