naginterfaces.library.correg.robustm_​user

naginterfaces.library.correg.robustm_user(psi, psip0, beta, indw, isigma, x, y, wgt, theta, sigma, chi=None, tol=5e-05, eps=5e-06, maxit=50, nitmon=0, data=None, io_manager=None)

robustm_user performs bounded influence regression ($M$-estimates) using an iterative weighted least squares algorithm.

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

https://support.nag.com/numeric/nl/nagdoc_30.3/flhtml/g02/g02hdf.html

Parameters

psicallable retval = psi(t, data=None)

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

Parameters

tfloat

The argument for which $psi$ must be evaluated.

dataarbitrary, optional, modifiable in place

User-communication data for callback functions.

Returns

retvalfloat

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

psip0float

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

betafloat

If $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 stat.inv_cdf_normal).

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 specfun.cdf_normal).

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

$$\begin{array}{c}\begin{array}{ccc}{\beta }_2=& {\int }_{-\infty }^{\infty }\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{\int }_{-\infty }^{\infty }\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{\int }_{-\infty }^{\infty }\chi \left(z/w_i\right)\varphi \left(z\right)dz\text{,}& \text{in the Schweppe case;}\end{array}\end{array}$$

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

If $isigma=0$, $beta$ is not referenced.

indwint

Determines the type of regression to be performed.

$indw=0$

Huber type regression.

$indw<0$

Mallows type regression.

$indw>0$

Schweppe type regression.

isigmaint

Determines how $\sigma$ is to be estimated.

$isigma=0$

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

$isigma<0$

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

$isigma>0$

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

xfloat, array-like, shape $(n,m)$

The values of the $X$ matrix, i.e., the independent variables. $x[\text{i}-1,\text{j}-1]$ must contain the $\text{i}\text{j}$th element of $x$, for $\text{j}=1,2,\dots ,m$, for $\text{i}=1,2,\dots ,n$.

If $indw<0$, during calculations the elements of $x$ will be transformed as described in Notes.

Before exit the inverse transformation will be applied.

As a result there may be slight differences between the input $x$ and the output $x$.

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

The data values of the dependent variable.

$y[\text{i}-1]$ must contain the value of $y$ for the $\text{i}$th observation, for $\text{i}=1,2,\dots ,n$.

If $indw<0$, during calculations the elements of $y$ will be transformed as described in Notes.

Before exit the inverse transformation will be applied.

As a result there may be slight differences between the input $y$ and the output $y$.

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

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

If $indw<0$, during calculations elements of $wgt$ will be transformed as described in Notes.

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 $wgt[i-1]\le 0$, the $i$th observation is not included in the analysis.

If $indw=0$, $wgt$ is not referenced.

thetafloat, array-like, shape $\left(m\right)$

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

sigmafloat

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.

chiNone or callable retval = chi(t, data=None), optional

Note: if this argument is None then a NAG-supplied facility will be used.

If $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.

Parameters

tfloat

The argument for which $chi$ must be evaluated.

dataarbitrary, optional, modifiable in place

User-communication data for callback functions.

Returns

retvalfloat

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

tolfloat, optional

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 \text{machine precision}$.

epsfloat, optional

A relative tolerance to be used to determine the rank of $X$. See linsys.real_gen_solve for further details.

If $eps<\text{machine precision}$ or $eps>1.0$, 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.

maxitint, optional

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

A value of $maxit=50$ should be adequate for most uses.

nitmonint, optional

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

$nitmon\le 0$

No information is printed.

$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 file object associated with the advisory I/O unit (see FileObjManager).

dataarbitrary, optional

User-communication data for callback functions.

io_managerFileObjManager, optional

Manager for I/O in this routine.

Returns

xfloat, ndarray, shape $(n,m)$

Unchanged, except as described above.

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

Unchanged, except as described above.

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

Unchanged, except as described above.

thetafloat, ndarray, shape $\left(m\right)$

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

kint

The column rank of the matrix $X$.

sigmafloat

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

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

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

nitint

The number of iterations that were used during the estimation.

Raises

NagValueError

(errno $1$)

On entry, $n=⟨value⟩$ and $\text{ldx}=⟨value⟩$.

Constraint: $\text{ldx}\ge n$.

(errno $1$)

On entry, $n=⟨value⟩$ and $m=⟨value⟩$.

Constraint: $n>m$.

(errno $1$)

On entry, $m=⟨value⟩$.

Constraint: $m\ge 1$.

(errno $1$)

On entry, $n=⟨value⟩$.

Constraint: $n\ge 2$.

(errno $2$)

On entry, $beta=⟨value⟩$.

Constraint: $beta>0.0$.

(errno $2$)

On entry, $sigma=⟨value⟩$.

Constraint: $sigma>0.0$.

(errno $3$)

On entry, $maxit=⟨value⟩$.

Constraint: $maxit>0$.

(errno $3$)

On entry, $tol=⟨value⟩$.

Constraint: $tol>0.0$.

(errno $4$)

Value given by $chi$ function $<0$: $chi(⟨value⟩)=⟨value⟩$.

The value of $chi$ must be non-negative.

(errno $5$)

Estimated value of $sigma$ is zero.

(errno $6$)

Iterations to solve the weighted least squares equations failed to converge.

(errno $8$)

The function has failed to converge in $maxit$ iterations.

(errno $9$)

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

Warns

NagAlgorithmicWarning

(errno $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 $errno$ = 7 and $k=m$ is possible (see linsys.real_gen_solve).

Notes

For the linear regression model

$$y=X\theta +ϵ\text{,}$$

where	$y$ is a vector of length $n$ of the dependent variable,
	$X$ is an $n\times m$ matrix of independent variables of column rank $k$,
	$\theta$ is a vector of length $m$ of unknown parameters,
and	$ϵ$ is a vector of length $n$ of unknown errors with var $\left({ϵ}_i\right)={\sigma }^2$,

robustm_user calculates the M-estimates given by the solution, $\hat{\theta }$, to the equation

$$\sum _{i=1}^n\psi \left(r_i/\left(\sigma w_i\right)\right)w_i{x}_{ij}=0\text{,}j=1,2,\dots ,m\text{,}$$

where	$r_i$ is the $i$th residual, i.e., the $i$th element of the vector $r=y-X\hat{\theta }$,
	$\psi$ is a suitable weight function,
	$w_i$ are suitable weights such as those that can be calculated by using output from robustm_wts(),
and	$\sigma$ may be estimated at each iteration by the median absolute deviation of the residuals $\hat{\sigma }={med}_i\left(\left[\left|r_i\right|\right]/{\beta }_1\right)$

or as the solution to

$$\sum _{i=1}^n\chi \left(r_i/\left(\hat{\sigma }{w}_i\right)\right)w_i^2=(n-k){\beta }_2$$

for a suitable weight function $\chi$, where ${\beta }_1$ and ${\beta }_2$ are constants, chosen so that the estimator of $\sigma$ is asymptotically unbiased if the errors, ${ϵ}_i$, have a Normal distribution. Alternatively $\sigma$ may be held at a constant value.

The above describes the Schweppe type regression. If the $w_i$ are assumed to equal $1$ for all $i$, then Huber type regression is obtained. A third type, due to Mallows, replaces (1) by

$$\sum _{i=1}^n\psi \left(r_i/\sigma \right)w_i{x}_{ij}=0\text{,}j=1,2,\dots ,m\text{.}$$

This may be obtained by use of the transformations

$$\begin{array}{c}\begin{array}{cc}{w}_i^{\ast }& \leftarrow \sqrt{w_i}\\ y_i^{\ast }& \leftarrow y_i\sqrt{w_i}\\ x_{ij}^{\ast }& \leftarrow x_{ij}\sqrt{w_i}\text{,}j=1,2,\dots ,m\end{array}\end{array}$$

(see Marazzi (1987)).

The calculation of the estimates of $\theta$ can be formulated as an iterative weighted least squares problem with a diagonal weight matrix $G$ given by

$$\begin{array}{c}{G}_{ii}=\{\begin{array}{cc}\frac{\psi \left(r_i/\left(\sigma w_i\right)\right)}{\left(r_i/\left(\sigma w_i\right)\right)}\text{,}& r_i\ne 0\\ & \\ {\psi }^{\prime }\left(0\right)\text{,}& r_i=0\text{.}\end{array}\text{.}\end{array}$$

The value of $\theta$ at each iteration is given by the weighted least squares regression of $y$ on $X$. This is carried out by first transforming the $y$ and $X$ by

$$\begin{array}{c}\begin{array}{cc}{\tilde{y}}_i& =y_i\sqrt{G_{ii}}\\ {\tilde{x}}_{ij}& =x_{ij}\sqrt{G_{ii}}\text{,}j=1,2,\dots ,m\end{array}\end{array}$$

and then using linsys.real_gen_solve. If $X$ is of full column rank then an orthogonal-triangular ($QR$) decomposition is used; if not, a singular value decomposition is used.

Observations with zero or negative weights are not included in the solution.

Note: there is no explicit provision in the function for a constant term in the regression model. However, the addition of a dummy variable whose value is $1.0$ for all observations will produce a value of $\hat{\theta }$ corresponding to the usual constant term.

robustm_user is based on routines in ROBETH, 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 and bounded influence regression in ROBETH, Cah. Rech. Doc. IUMSP, No. 3 ROB 2, Institut Universitaire de Médecine Sociale et Préventive, Lausanne