1
mscale <- function(u, na.rm=FALSE)
2
# Scale M-estimator with 50% breakdown
3
# Yohai (1987) Annals, Stromberg (1993) JASA.
6
# Revised 17 April 2010
19
if(mean(u==0) >= 0.5) return(0)
21
s <- median(U)/0.6744898
26
d1 <- mean(.rho.hampel(z))-3.75
27
d2 <- mean(z*.psi.hampel(z))
30
warning("Max iterations exceeded")
33
if(abs(d1/d2) < 1e-13) break
38
.rho.hampel <- function(u, a = 1.5, b = 3.5, c = 8)
40
# Integral of Hampel's redescending psi function (Hampel, Ronchetti,
41
# Rousseeuw and Stahel, 1986, Robust Statistics, Wiley, page 150).
42
# Default values are as in Stromberg (1993) JASA.
47
A <- (U <= a) #increasing
48
B <- (U > a) & (U <= b) #flat
49
C <- (U > b) & (U <= c) #descending
52
rho[A] <- (U[A] * U[A])/2
53
rho[B] <- a * (U[B] - a/2)
54
rho[C] <- a * (b - a/2) + a * (U[C] - b) * (1 - (U[C] - b)/(c - b)/2)
55
rho[D] <- (a * (b - a + c))/2
59
.psi.hampel <- function(u, a = 1.5, b = 3.5, c = 8)
61
# Hampel's redescending psi function (Hampel, Ronchetti,
62
# Rousseeuw and Stahel, 1986, Robust Statistics, Wiley, page 150).
63
# Default values are as in Stromberg (1993) JASA.
68
B <- (U > a) & (U <= b) #flat
69
C <- (U > b) & (U <= c) #descending
72
psi[B] <- sign(u[B]) * a
73
psi[C] <- sign(u[C]) * a * (c - U[C])/(c - b)
b'\\ No newline at end of file'
added
removed
R/mscale.R
