# inst/scripts/NormalLocationModel.R In ROptEst: Optimally Robust Estimation

```###############################################################################
## Example: Normal Location
###############################################################################
require(ROptEst)
options("newDevice"=TRUE)

## generates Normal Location Family with scale = 3
N0 <- NormLocationFamily(mean=2, sd=3)
N0        # show G0
plot(N0)  # plot of Norm(mean = 2, sd = 3) and L_2 derivative
checkL2deriv(N0)

## classical optimal IC
N0.IC0 <- optIC(model = N0, risk = asCov())
N0.IC0       # show IC
plot(N0.IC0) # plot IC
checkIC(N0.IC0)
Risks(N0.IC0)

## L_2 family + infinitesimal neighborhood
N0.Rob1 <- InfRobModel(center = N0, neighbor = ContNeighborhood(radius = 0.5))
N0.Rob1     # show N0.Rob1
N0.Rob2 <- InfRobModel(center = N0, neighbor = TotalVarNeighborhood(radius = 0.5))

## OBRE solution (ARE 95%)

system.time(N0.ICA <- optIC(model = N0.Rob1, risk = asAnscombe(.95),upper=NULL,lower=NULL, verbose=TRUE))
checkIC(N0.ICA)
Risks(N0.ICA)
plot(N0.ICA)

system.time(N0.ICA2 <- optIC(model = N0.Rob2, risk = asAnscombe(.95),upper=NULL,lower=NULL, verbose=TRUE))
checkIC(N0.ICA2)
Risks(N0.ICA2)
plot(N0.ICA2)

## MSE solution
(N0.IC1 <- optIC(model=N0.Rob1, risk=asMSE()))
checkIC(N0.IC1)
Risks(N0.IC1)
plot(N0.IC1)

(N0.IC2 <- optIC(model=N0.Rob2, risk=asMSE()))
checkIC(N0.IC2)
Risks(N0.IC2)
plot(N0.IC2)

## L1 solution

(N0.IC1.L1 <- optIC(model=N0.Rob1, risk=asL1()))
checkIC(N0.IC1.L1)
Risks(N0.IC1.L1)
plot(N0.IC1.L1)

(N0.IC2.L1 <- optIC(model=N0.Rob2, risk=asL1()))
checkIC(N0.IC2.L1)
Risks(N0.IC2.L1)
plot(N0.IC2.L1)

## L4 solution

(N0.IC1.L4 <- optIC(model=N0.Rob1, risk=asL4()))
checkIC(N0.IC1.L4)
Risks(N0.IC1.L4)
plot(N0.IC1.L4)

(N0.IC2.L4 <- optIC(model=N0.Rob2, risk=asL4()))
checkIC(N0.IC2.L4)
Risks(N0.IC2.L4)
plot(N0.IC2.L4)

## semivar-solution

(N0.IC1.SemiV.p <- optIC(model=N0.Rob1, risk=asSemivar(sign=1)))
checkIC(N0.IC1.SemiV.p)
Risks(N0.IC1.SemiV.p)
plot(N0.IC1.SemiV.p)

(N0.IC1.SemiV.m <- optIC(model=N0.Rob1, risk=asSemivar(sign=-1)))
checkIC(N0.IC1.SemiV.m)
Risks(N0.IC1.SemiV.m)
plot(N0.IC1.SemiV.m)

## lower case solutions
(N0.IC3 <- optIC(model=N0.Rob1, risk=asBias()))
checkIC(N0.IC3)
Risks(N0.IC3)
plot(N0.IC3)
(N0.IC4 <- optIC(model=N0.Rob2, risk=asBias()))
checkIC(N0.IC4)
Risks(N0.IC4)
plot(N0.IC4)

## Hampel solution
(N0.IC5 <- optIC(model=N0.Rob1, risk=asHampel(bound=clip(N0.IC1))))
checkIC(N0.IC5)
Risks(N0.IC5)
plot(N0.IC5)
(N0.IC6 <- optIC(model=N0.Rob2, risk=asHampel(bound=Risks(N0.IC2)\$asBias\$value), maxiter = 200))
checkIC(N0.IC6)
Risks(N0.IC6)
plot(N0.IC6)

checkIC(N0.IC7)
Risks(N0.IC7)
plot(N0.IC7)
checkIC(N0.IC8)
Risks(N0.IC8)
plot(N0.IC8)

neighbor.0 <- ContNeighborhood()
getReq(asMSE(),neighbor.0,N0.ICA,N0.IC1,n=1)
getReq(asMSE(),neighbor.0,N0.ICA,N0.IC1,n=30)
getReq(asL1(),neighbor.0,N0.ICA,N0.IC1,n=30)
getReq(asL4(),neighbor.0,N0.ICA,N0.IC1,n=30)
getReq(asMSE(),neighbor.0,N0.ICA,N0.IC7,n=30)
getReq(asL1(),neighbor.0,N0.ICA,N0.IC7,n=30)
getReq(asL4(),neighbor.0,N0.ICA,N0.IC7,n=30)
getReq(asMSE(),neighbor.0,N0.IC1,N0.IC7,n=30)

getMaxIneff(N0.ICA,neighbor=ContNeighborhood())
getMaxIneff(N0.IC1,neighbor=ContNeighborhood())
getMaxIneff(N0.IC7,neighbor=ContNeighborhood())

## (may take quite some time!)
risk=asMSE(), rho=0.5))
risk=asMSE(), rho=1/3))

## For estimation use function roptest
ind <- rbinom(1e2, size=1, prob=0.05)
x <- rnorm(1e2, mean=1, sd = (1-ind)+ind*9)

## 1-step: contamination known
est1 <- roptest(x, eps = 0.05, L2Fam = NormLocationFamily(mean = 1))
est1v <- roptest(x, eps = 0.025, L2Fam = NormLocationFamily(mean = 1),
neighbor = TotalVarNeighborhood())

## k-step: contamination known
est2 <- roptest(x, eps = 0.05, L2Fam = NormLocationFamily(mean = 1), steps = 3)
est2v <- roptest(x, eps = 0.025, L2Fam = NormLocationFamily(mean = 1),
neighbor = TotalVarNeighborhood(), steps = 3)

## comparison
estimate(est1)
estimate(est2)
estimate(est1v)
estimate(est2v)

## confidence intervals
confint(est1, symmetricBias())
confint(est2, symmetricBias())
confint(est1v, symmetricBias())
confint(est2v, symmetricBias())
```

## Try the ROptEst package in your browser

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

ROptEst documentation built on April 6, 2019, 3:01 a.m.