Mscale

Mscale finds the M estimator of the scale

expand all in page

Syntax

  • sc=Mscale(u, psifunc)example
  • sc=Mscale(u, psifunc, initialsc)example
  • sc=Mscale(u, psifunc, initialsc, tol)example
  • sc=Mscale(u, psifunc, initialsc, tol, maxiter)example

Description

example

sc =Mscale(u, psifunc) Example of M estimate of scale.

example

sc =Mscale(u, psifunc, initialsc) Estimate of scale using Hampel rho function.

example

sc =Mscale(u, psifunc, initialsc, tol) Use of options initialsc, tol and maxiter.

example

sc =Mscale(u, psifunc, initialsc, tol, maxiter) Compare scale estimate using two different link functions.

Examples

expand all

  • Example of M estimate of scale.

    • M estimate of the scale using Tukey biweight rho function with a value of c associated to a breakdown point of 0.5. 

    psifunc=struct;
psifunc.class='TB';
bdp=0.5;
c=TBbdp(bdp,1);
% kc = E(rho) = sup(rho)*bdp
kc=c^2/6*bdp;
psifunc.c1=c;
psifunc.kc1=kc;
n=10000;
shift=5;
u=2*randn(n,1);
u(1:10)=u(1:10)+shift;
s=Mscale(u,psifunc)

  • Estimate of scale using Hampel rho function.

    • M estimate of the scale using Hampel rho function with a value of c associated to a breakdown point of 0.5 

    psifunc=struct;
psifunc.class='HA'
abc=[1.5 3.5 8];
bdp=0.5;
c=HAbdp(bdp,1,abc);
% kc = E(rho) = sup(rho)*bdp
kc=HArho(c*abc(3),[c, abc])*bdp;
psifunc.c1=[c abc];
psifunc.kc1=kc;
n=10000;
shift=5;
u=3*randn(n,1);
u(1:10)=u(1:10)+shift;
s=Mscale(u,psifunc)

  • Use of options initialsc, tol and maxiter.

    • M estimate of the scale using Tukey biweight rho function with a value of c associated to a breakdown point of 0.5. 

    psifunc=struct;
psifunc.class='TB';
bdp=0.5;
c=TBbdp(bdp,1);
% kc = E(rho) = sup(rho)*bdp
kc=c^2/6*bdp;
psifunc.c1=c;
psifunc.kc1=kc;
n=10000;
shift=5;
u=2*randn(n,1);
u(1:10)=u(1:10)+shift;
s=Mscale(u,psifunc,3,1e-7,20)

  • Compare scale estimate using two different link functions.

    • psifunc=struct;
psifunc.class='HA'
abc=[1.5 3.5 8];
bdp=0.5;
c=HAbdp(bdp,1,abc);
% kc = E(rho) = sup(rho)*bdp
kc=HArho(c*abc(3),[c, abc])*bdp;
psifunc.c1=[c, abc];
psifunc.kc1=kc;
n=10000;
shift=100;
u=3*randn(n,1);
u(100:200)=u(100:200)+shift;
sHA=Mscale(u,psifunc)
% kc = E(rho) = sup(rho)*bdp
psifunc=struct;
psifunc.class='TB';
c=TBbdp(bdp,1);
psifunc.c1=c;
kc=c^2/6*bdp;
psifunc.kc1=kc;
sTB=Mscale(u,psifunc,3,1e-7,20);
sMLE=std(u);
cate=categorical({'Robust scale (Hampel)' 'Robust scale (TB)' 'Non robust scale'});
bar(cate,[sHA sTB sMLE])

    Input Arguments

    expand all

    u — residuals or Mahalanobis distances. Vector.

    n x 1 vector which contains the (unscaled) residuals or Mahalanobis distances

    Data Types: single | double

    psifunc — rho (psi) function. Structure.

    A structure specifying the class of rho (psi) function to use, the consistency factor, and the value associated with the Expectation of rho in correspondence of the consistency factor.

    psifunc must contain the following fields:

    Value Description
    class

    string identyfing the rho (psi) function to use.

    Admissible values for class are 'bisquare' (TB), 'optimal', (OPT) 'hyperbolic' (HYP), 'hampel' (HA) 'power divergence' (PD) and 'Andrew's sine function0 (AS).
    c1

    consistency factor (and other parameters) associated to required breakdown point or nominal efficiency.

    More precisely, psifunc.c1(1) contains consistency factor associated to required breakdown point or nominal efficiency psifunc.c1(2:end) contain other parameters associated with the rho (psi) function.

    For example, if psifunc.class='hampel', c1(2:4) must contain parameters (a, b and c) of Hampel rho function.
    kc1

    Expectation of rho associated with c1 to get a consistent estimator at the model distribution kc1=E(rho)

    Example: Example - psifunc.class='TB';psifunc.c1=1.5476;psifunc.kc1=0.1996

    Data Types: struct

    Optional Arguments

    initialsc — scalar. The initial estimate of the scale.

    If not defined, scaled MAD of vector |u| is used.

    Example: 'initialsc',0.34

    Data Types: double

    tol — scalar. The tolerance for controlling convergence.

    If not defined, tol is fixed to 1e-7.

    Example: 'tol',1e-10

    Data Types: double

    maxiter — scalar. Maximum number of iterations to find the scale.

    If not defined, maxiter is fixed to 200.

    Example: 'maxiter',100

    Data Types: double

    Output Arguments

    expand all

    sc —M-estimate of the scale. Scalar

    Robust M estimate of scale.

    This routine is called by Taureg and Sreg.m

    More About

    expand all

    Additional Details

    u = residuals or Mahalanobis distances (note that u is kept fixed in each iteration).

    Remark: the M estimator of scale must satisfy the following equation

    \[ (1/n) \sum_{i=1}^n \rho((u_i/c)/s) = kc \]

    This routine computes the value of $s$ which satisfies the above equation.

    References

    Huber, P.J. and Ronchetti, E.M. (2009), "Robust Statistics, 2nd Edition", Wiley. [equation 7.119, p. 176].

    See Also

    |

    This page has been automatically generated by our routine publishFS