# R/rlsOptIC_AL.R In RobLox: Optimally Robust Influence Curves and Estimators for Location and Scale

#### Documented in rlsOptIC.AL

```###############################################################################
## weight function
###############################################################################
.ALrlsGetw <- function(x, b, a1, a2, a3){
hvkt <- sqrt(a3^2*x^4 + (a1^2 - 2*a2*a3^2)*x^2 + a2^2*a3^2)
ind1 <- (hvkt < b)

return(ind1 + (1-ind1)*b/hvkt)
}

###############################################################################
## computation of r
###############################################################################
.ALrlsGetr <- function(b, r, a1, a2, a3){
integrandr <- function(x, b, a1, a2, a3){
hvkt <- sqrt(a3^2*x^4 + (a1^2 - 2*a2*a3^2)*x^2 + a2^2*a3^2)/b - 1
return((hvkt > 0)*hvkt*dnorm(x))
}
Int <- integrate(integrandr, lower = 0, upper = Inf,
rel.tol = .Machine\$double.eps^0.5, a1 = a1, a2 = a2,
a3 = a3, b = b)\$value

return(r-sqrt(2*Int))
}

###############################################################################
## computation of a1, a2 and a3
###############################################################################
.ALrlsGeta1a2a3 <- function(b, a1, a2, a3){
integrand1 <- function(x, b, a1, a2, a3){
x^2*.ALrlsGetw(x, b, a1, a2, a3)*dnorm(x)
}
Int1 <- 2*integrate(integrand1, lower = 0, upper = Inf,
rel.tol = .Machine\$double.eps^0.5, b = b, a1 = a1,
a2 = a2, a3 = a3)\$value
a1 <- 1/Int1

integrand2 <- function(x, b, a1, a2, a3){
.ALrlsGetw(x, b, a1, a2, a3)*dnorm(x)
}
Int2 <- 2*integrate(integrand2, lower = 0, upper = Inf,
rel.tol = .Machine\$double.eps^0.5, b = b, a1 = a1,
a2 = a2, a3 = a3)\$value
a2 <- Int1/Int2

integrand3 <- function(x, b, a1, a2, a3){
(x^2 - a2)^2*.ALrlsGetw(x, b, a1, a2, a3)*dnorm(x)
}
Int3 <- 2*integrate(integrand3, lower = 0, upper = Inf,
rel.tol = .Machine\$double.eps^0.5, b = b, a1 = a1,
a2 = a2, a3 = a3)\$value
a3 <- 1/Int3

return(list(a1=a1, a2=a2, a3=a3))
}
.ALrlsVar <- function(b, a1, a2, a3){
integrand1 <- function(x, b, a1, a2, a3){
x^2*.ALrlsGetw(x, b, a1, a2, a3)^2*dnorm(x)
}
Int1 <- 2*integrate(integrand1, lower = 0, upper = Inf,
rel.tol = .Machine\$double.eps^0.5, b = b, a1 = a1,
a2 = a2, a3 = a3)\$value
V1 <- a1^2*Int1

integrand2 <- function(x, b, a1, a2, a3){
(x^2 - a2)^2*.ALrlsGetw(x, b, a1, a2, a3)^2*dnorm(x)
}
Int2 <- 2*integrate(integrand2, lower = 0, upper = Inf,
rel.tol = .Machine\$double.eps^0.5, b = b, a1 = a1,
a2 = a2, a3 = a3)\$value
V2 <- a3^2*Int2

return(diag(c(V1, V2)))
}

###############################################################################
## optimal IC
###############################################################################
rlsOptIC.AL <- function(r, mean = 0, sd = 1, A.loc.start = 1, a.sc.start = 0,
A.sc.start = 0.5, bUp = 1000, delta = 1e-6, itmax = 100,
check = FALSE, computeIC = TRUE){
a1 <- A.loc.start; a2 <- 1+a.sc.start; a3 <- A.sc.start
b <- uniroot(.ALrlsGetr, lower = 1e-4, upper = bUp,
tol = .Machine\$double.eps^0.5, r = r, a1 = a1, a2 = a2,
a3 = a3)\$root

iter <- 0
repeat{
iter <- iter + 1
if(iter > itmax){
stop("Algorithm did not converge!\n",
"=> increase itmax or try different starting values",
"'A.loc.start', 'a.sc.start' and 'A.sc.start'\n")
}
a1.old <- a1; a2.old <- a2; a3.old <- a3; b.old <- b

a1a2a3 <- .ALrlsGeta1a2a3(b = b, a1 = a1, a2 = a2, a3 = a3)
a1 <- a1a2a3\$a1
a2 <- a1a2a3\$a2
a3 <- a1a2a3\$a3

b <- uniroot(.ALrlsGetr, lower = 1e-4, upper = bUp,
tol = .Machine\$double.eps^0.5, r = r, a1 = a1, a2 = a2,
a3 = a3)\$root
if(max(abs(a1.old-a1), abs(a2.old-a2), abs(a3.old-a3), abs(b.old-b))<delta)
break
}

if(check){
integrand1 <- function(x, b, a1, a2, a3){
x^2*.ALrlsGetw(x, b, a1, a2, a3)*dnorm(x)
}
Int1 <- 2*integrate(integrand1, lower = 0, upper = Inf,
rel.tol = .Machine\$double.eps^0.5, b = b, a1 = a1,
a2 = a2, a3 = a3)\$value
ch1 <- a1*Int1

integrand2 <- function(x, b, a1, a2, a3){
(x^2 - a2)^2*.ALrlsGetw(x, b, a1, a2, a3)*dnorm(x)
}
Int2 <- 2*integrate(integrand2, lower = 0, upper = Inf,
rel.tol = .Machine\$double.eps^0.5, b = b, a1 = a1,
a2 = a2, a3 = a3)\$value
ch2 <- a3*Int2

integrand3 <- function(x, b, a1, a2, a3){
(x^2 - a2)*.ALrlsGetw(x, b, a1, a2, a3)*dnorm(x)
}
Int3 <- 2*integrate(integrand3, lower=0, upper=Inf,
rel.tol = .Machine\$double.eps^0.5, b = b, a1 = a1,
a2 = a2, a3 = a3)\$value
ch3 <- a3*Int3

ch4 <- .ALrlsGetr(b = b, r = r, a1 = a1, a2 = a2, a3 = a3)

cat("Fisher consistency of eta.loc:\t", ch1-1, "\n")
cat("centering of eta.sc:\t", ch3, "\n")
cat("Fisher consistency of eta.sc:\t", ch2-1, "\n")
cat("MSE equation:\t", ch4, "\n")
}

asVar <- sd^2*.ALrlsVar(b = b, a1 = a1, a2 = a2, a3 = a3)
A <- sd^2*diag(c(a1, a3))
a <- sd*c(0, a3*(a2-1))
b <- sd*b
mse <- sd^2*(a1 + a3)

if(computeIC){
w <- new("HampelWeight")
clip(w) <- b
cent(w) <- c(0, a2-1)/sd
stand(w) <- A
weight(w) <- getweight(w, neighbor = ContNeighborhood(radius = r),
biastype = symmetricBias(),
normW = NormType())

modIC <- function(L2Fam, IC, withMakeIC, ...){
ICL2Fam <- eval(CallL2Fam(IC))
if(is(L2Fam, "L2LocationScaleFamily") && is(distribution(L2Fam), "Norm")){
sdneu <- main(L2Fam)[2]
sdalt <- main(ICL2Fam)[2]
w <- weight(IC)
clip(w) <- sdneu*clip(w)/sdalt
cent(w) <- sdalt*cent(w)/sdneu
stand(w) <- sdneu^2*stand(w)/sdalt^2
biastype = biastype(IC),
normW = normtype(IC))
A <- sdneu^2*stand(IC)/sdalt^2
b <- sdneu*clip(IC)/sdalt
mse <- sum(diag(A))
a1 <- A[1, 1]/sdneu^2
a3 <- A[2, 2]/sdneu^2
a2 <- a[1]/sd/a3 + 1
asVar <- sdneu^2*.ALrlsVar(b = b/sdneu, a1 = a1, a2 = a2, a3 = a3)
res <- list(A = A, a = sdneu*cent(IC)/sdalt, b = b, d = NULL,
risk = list(asMSE = mse, asBias = b, trAsCov = mse - r^2*b^2,
asCov = asVar),
info = Infos(IC), w = w,
normtype = normtype(IC), biastype = biastype(IC),
modifyIC = modifyIC(IC))
IC <- generateIC(neighbor = ContNeighborhood(radius = r),
L2Fam = L2Fam, res = res)
if(!any(grepl("Some entries in 'Infos' may be wrong", Infos(IC)[,2]))){
addInfo(IC) <- c("modifyIC", "The IC has been modified")
addInfo(IC) <- c("modifyIC", "Some entries in 'Infos' may be wrong")
}
return(IC)
}else{
makeIC(IC, L2Fam, ...)
}
}

L2Fam <- substitute(NormLocationScaleFamily(mean = m1, sd = s1),
list(m1 = mean, s1 = sd))
L2Fam = eval(L2Fam),
res = list(A = as.matrix(A), a = a, b = b, d = NULL,
risk = list(asMSE = mse, asBias = b, trAsCov = mse - r^2*b^2,
asCov = asVar),
info = c("rlOptIC", "optimally robust IC for AL estimators and 'asMSE'"),
w = w, biastype = symmetricBias(), normtype = NormType(),
modifyIC = modIC)))
}else{
return(list(A = A, a = a, b = b))
}
}
```

## Try the RobLox package in your browser

Any scripts or data that you put into this service are public.

RobLox documentation built on April 6, 2019, 3:04 a.m.