naginterfaces.library.univar.robust_​1var_​mestim

naginterfaces.library.univar.robust_1var_mestim(isigma, x, ipsi, c, h1, h2, h3, dchi, theta, sigma, tol, maxit=50)

robust_1var_mestim computes an $M$-estimate of location with (optional) simultaneous estimation of the scale using Huber’s algorithm.

For full information please refer to the NAG Library document for g07db

https://support.nag.com/numeric/nl/nagdoc_30.3/flhtml/g07/g07dbf.html

Parameters

isigmaint

The value assigned to $isigma$ determines whether $\hat{\sigma }$ is to be simultaneously estimated.

$isigma=0$

The estimation of $\hat{\sigma }$ is bypassed and $sigma$ is set equal to ${\sigma }_c$.

$isigma=1$

$\hat{\sigma }$ is estimated simultaneously.

xfloat, array-like, shape $\left(n\right)$

The vector of observations, $x_1,x_2,\dots ,x_n$.

ipsiint

Which $\psi$ function is to be used.

$ipsi=0$

$\psi \left(t\right)=t$.

$ipsi=1$

Huber’s function.

$ipsi=2$

Hampel’s piecewise linear function.

$ipsi=3$

Andrew’s sine wave,

$ipsi=4$

Tukey’s bi-weight.

cfloat

If $ipsi=1$, $c$ must specify the parameter, $c$, of Huber’s $\psi$ function. $c$ is not referenced if $ipsi\ne 1$.

h1float

If $ipsi=2$, $h1$, $h2$ and $h3$ must specify the parameters, $h_1$, $h_2$, and $h_3$, of Hampel’s piecewise linear $\psi$ function. $h1$, $h2$ and $h3$ are not referenced if $ipsi\ne 2$.

h2float

If $ipsi=2$, $h1$, $h2$ and $h3$ must specify the parameters, $h_1$, $h_2$, and $h_3$, of Hampel’s piecewise linear $\psi$ function. $h1$, $h2$ and $h3$ are not referenced if $ipsi\ne 2$.

h3float

If $ipsi=2$, $h1$, $h2$ and $h3$ must specify the parameters, $h_1$, $h_2$, and $h_3$, of Hampel’s piecewise linear $\psi$ function. $h1$, $h2$ and $h3$ are not referenced if $ipsi\ne 2$.

dchifloat

$d$, the parameter of the $\chi$ function. $dchi$ is not referenced if $ipsi=0$.

thetafloat

If $sigma>0$ then $theta$ must be set to the required starting value of the estimation of the location parameter $\hat{\theta }$. A reasonable initial value for $\hat{\theta }$ will often be the sample mean or median.

sigmafloat

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 calculations of $\hat{\theta }$ and $\hat{\sigma }$. If $sigma\le 0.0$ then robust_1var_mestim will determine the starting points of $\hat{\theta }$ and $\hat{\sigma }$. Otherwise the value assigned to $sigma$ will be taken as the starting point for $\hat{\sigma }$, and $theta$ must be assigned a value before entry, see above;

if $isigma=0$, $sigma$ must be assigned a value which determines the value of ${\sigma }_c$, which is held fixed during the iterations, and the starting value for the calculation of $\hat{\theta }$. If $sigma\le 0$, robust_1var_mestim will determine the value of ${\sigma }_c$ as the median absolute deviation adjusted to reduce bias (see robust_1var_median()) and the starting point for $\hat{\theta }$. Otherwise, the value assigned to $sigma$ will be taken as the value of ${\sigma }_c$ and $theta$ must be assigned a relevant value before entry, see above.

tolfloat

The relative precision for the final estimates. Convergence is assumed when the increments for $theta$, and $sigma$ are less than $tol\times max(1.0,{\sigma }_{k-1})$.

maxitint, optional

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

Returns

thetafloat

The $M$-estimate of the location parameter, $\hat{\theta }$.

sigmafloat

Contains the $M$-estimate of the scale parameter, $\hat{\sigma }$, if $isigma$ was assigned the value $1$ on entry, otherwise $sigma$ will contain the initial fixed value ${\sigma }_c$.

rsfloat, ndarray, shape $\left(n\right)$

The Winsorized residuals.

nitint

The number of iterations that were used during the estimation.

wrkfloat, ndarray, shape $\left(n\right)$

If $sigma\le 0.0$ on entry, $wrk$ will contain the $n$ observations in ascending order.

Raises

NagValueError

(errno $1$)

On entry, $ipsi=⟨value⟩$.

Constraint: $ipsi=0$, $1$, $2$, $3$ or $4$.

(errno $1$)

On entry, $isigma=⟨value⟩$.

Constraint: $isigma=0$ or $1$.

(errno $1$)

On entry, $maxit=⟨value⟩$.

Constraint: $maxit>0$.

(errno $1$)

On entry, $tol=⟨value⟩$.

Constraint: $tol>0.0$.

(errno $1$)

On entry, $n=⟨value⟩$.

Constraint: $n>1$.

(errno $2$)

On entry, $dchi=⟨value⟩$.

Constraint: $dchi>0.0$.

(errno $2$)

On entry, $h1=⟨value⟩$, $h2=⟨value⟩$ and $h3=⟨value⟩$.

Constraint: $0\le h1\le h2\le h3$ and $h3>0.0$.

(errno $2$)

On entry, $c=⟨value⟩$.

Constraint: $c>0.0$.

(errno $3$)

All elements of $x$ are equal.

(errno $4$)

Current estimate of $sigma$ is zero or negative: $sigma=⟨value⟩$.

(errno $5$)

Number of iterations required exceeds $maxit$: $maxit=⟨value⟩$.

(errno $6$)

All winsorized residuals are zero.

Notes

In the NAG Library the traditional C interface for this routine uses a different algorithmic base. Please contact NAG if you have any questions about compatibility.

The data consists of a sample of size $n$, denoted by $x_1,x_2,\dots ,x_n$, drawn from a random variable $X$.

The $x_i$ are assumed to be independent with an unknown distribution function of the form

$$F\left((x_i-\theta )/\sigma \right)$$

where $\theta$ is a location parameter, and $\sigma$ is a scale parameter. $M$-estimators of $\theta$ and $\sigma$ are given by the solution to the following system of equations:

$$\sum _{i=1}^n\psi \left((x_i-\hat{\theta })/\hat{\sigma }\right)=0$$

$$\sum _{i=1}^n\chi \left((x_i-\hat{\theta })/\hat{\sigma }\right)=(n-1)\beta$$

where $\psi$ and $\chi$ are given functions, and $\beta$ is a constant, such that $\hat{\sigma }$ is an unbiased estimator when $x_{\text{i}}$, for $\text{i}=1,2,\dots ,n$ has a Normal distribution. Optionally, the second equation can be omitted and the first equation is solved for $\hat{\theta }$ using an assigned value of $\sigma ={\sigma }_c$.

The values of $\psi \left(\frac{x_i-\hat{\theta }}{\hat{\sigma }}\right)\hat{\sigma }$ are known as the Winsorized residuals.

The following functions are available for $\psi$ and $\chi$ in robust_1var_mestim:

1. Null Weights

   $\psi \left(t\right)=t$

   $\chi \left(t\right)=\frac{t^2}{2}$

   Use of these null functions leads to the mean and standard deviation of the data.

2. Huber’s Function

   $\psi \left(t\right)=max(-c,min(c,t))$
   $\chi \left(t\right)=\frac{{\left|t\right|}^2}{2}$	$\left|t\right|\le d$
   $\chi \left(t\right)=\frac{d^2}{2}$	$\left|t\right|>d$

3. Hampel’s Piecewise Linear Function

   ${\psi }_{h_1,h_2,h_3}\left(t\right)=-{\psi }_{h_1,h_2,h_3}(-t)$
   ${\psi }_{h_1,h_2,h_3}\left(t\right)=t$	$0\le t\le h_1$
   ${\psi }_{h_1,h_2,h_3}\left(t\right)=h_1$	$h_1\le t\le h_2$
   ${\psi }_{h_1,h_2,h_3}\left(t\right)=h_1(h_3-t)/(h_3-h_2)$	$h_2\le t\le h_3$
   ${\psi }_{h_1,h_2,h_3}\left(t\right)=0$	$t>h_3$
   $\chi \left(t\right)=\frac{{\left|t\right|}^2}{2}$	$\left|t\right|\le d$
   $\chi \left(t\right)=\frac{d^2}{2}$	$\left|t\right|>d$

4. Andrew’s Sine Wave Function

   $\psi \left(t\right)=\sin \left(t\right)$	$-\pi \le t\le \pi$
   $\psi \left(t\right)=0$	otherwise
   $\chi \left(t\right)=\frac{{\left|t\right|}^2}{2}$	$\left|t\right|\le d$
   $\chi \left(t\right)=\frac{d^2}{2}$	$\left|t\right|>d$

5. Tukey’s Bi-weight

   $\psi \left(t\right)=t{(1-t^2)}^2$	$\left|t\right|\le 1$
   $\psi \left(t\right)=0$	otherwise
   $\chi \left(t\right)=\frac{{\left|t\right|}^2}{2}$	$\left|t\right|\le d$
   $\chi \left(t\right)=\frac{d^2}{2}$	$\left|t\right|>d$

where $c$, $h_1$, $h_2$, $h_3$ and $d$ are constants.

Equations (1) and (2) are solved by a simple iterative procedure suggested by Huber:

$${\hat{\sigma }}_k=\sqrt{\frac{1}{\beta (n-1)}\left(\sum _{i=1}^n\chi \left(\frac{x_i-{\hat{\theta }}_{k-1}}{{\hat{\sigma }}_{k-1}}\right)\right){\hat{\sigma }}_{k-1}^2}$$

and

$${\hat{\theta }}_k={\hat{\theta }}_{k-1}+\frac{1}{n}\sum _{i=1}^n\psi \left(\frac{x_i-{\hat{\theta }}_{k-1}}{{\hat{\sigma }}_k}\right){\hat{\sigma }}_k$$

or

$${\hat{\sigma }}_k={\sigma }_c\text{, if}\sigma \text{is fixed.}$$

The initial values for $\hat{\theta }$ and $\hat{\sigma }$ may either be user-supplied or calculated within robust_1var_mestim as the sample median and an estimate of $\sigma$ based on the median absolute deviation respectively.

robust_1var_mestim 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