################################################################################
## R Code für Seminarvortrag vom 02.12.05
################################################################################

## load ROptEst - also loads: distr, evd, distrEx, RandVar
require(ROptEst)

###############
## Ideal Model: object of class "L2ParamFamily"
###############
(N <- NormLocationScaleFamily())

## L2-derivative
Map(L2deriv(N)[[1]])

## Fisher information
FisherInfo(N)

## marginal distributions of L2-derivative
L2derivDistr(N)

## further properties -> semi-symbolic calculus
L2derivSymm(N)
L2derivDistrSymm(N)

## check validity of N
checkL2deriv(N)

## let's have a look
plot(N)


###############
## Influence Curve: object of class "IC"
###############
## matrix multiplication for Euclidean random variables
(scores <- solve(N@FisherInfo) %*% N@L2deriv)

## classical scores
(IC.N <- IC(Curve = EuclRandVarList(scores[,1]),
           CallL2Fam = call("NormLocationScaleFamily")))

## check validity of N
checkIC(IC.N)

## or
checkIC(IC = IC.N, L2Fam = N)

## let's have a look
plot(IC.N)


###############
## Neighborhood: object of class "UncondNeighborhood"
###############
## convex contamination neighborhood (Tukey's gross error model)
(ngh.c <- ContNeighborhood(radius = 0.5))

## total variation neighborhood
(ngh.v <- TotalVarNeighborhood(radius = 0.5))


###############
## Robust Model: object of class "RobModel"
###############
## infinitesimal robust model: neighborhood radius shrinking at a rate of sqrt(n)
(RN <- InfRobModel(center = N, neighbor = ngh.c))

## fixed-size robust model: neighborhood radius/size fixed
FixRobModel(center = N, neighbor = ngh.c)


###############
## Risks: object of class "RiskType"
###############
## asymptotic mean square error
asMSE()

## asymptotic "Hampel risk"
asHampel()

## asymptotic bias
asBias()


###############
## Optimally robust IC
###############
## * = c: radius = 0.5
(IC1 <- optIC(model = RN, risk = asMSE()))

## check IC properties
checkIC(IC1)

## let's have a look
plot(IC1)

## absolute and relative Information
infoPlot(IC1)

## special case normal location and scale
## comparison of 18 different estimators in Chapter 8 of Kohl (2005)
require(RobLox)

## * = c: radius = 0.5
(IC2 <- rlsOptIC.AL(r=0.5)) # clearly faster than "optIC"

## Huber's (1964) proposal 2
(IC2 <- rlsOptIC.Hu1(r=0.5))
plot(IC2)
infoPlot(IC2)


###############
## Discrete Models and the gap condition
###############
(P <- PoisFamily(lambda = 5))

## lower case radius (*=c)
lowerCaseRadius(L2Fam = P, neighbor = ngh.c, risk = asMSE())

## lower case radius (*=v)
lowerCaseRadius(L2Fam = P, neighbor = ngh.v, risk = asMSE())


###############
## Radius--minimax Estimator
###############
radiusMinimaxIC(L2Fam=NormLocationFamily(), neighbor=ContNeighborhood(), 
                risk=asMSE(), loRad=0, upRad=Inf)
radiusMinimaxIC(L2Fam=NormScaleFamily(), neighbor=ContNeighborhood(), 
                risk=asMSE(), loRad=0, upRad=Inf)
radiusMinimaxIC(L2Fam=N, neighbor=ContNeighborhood(), 
                risk=asMSE(), loRad=0, upRad=Inf)


###############
## One-Step Estimation
###############
require(MASS)
plot(chem)

## initial estimate
est0 <- ksEstimator(chem, Norm())

## ideal Model
N1 <- NormLocationScaleFamily(mean = est0$mean, sd = est0$sd)

## robust Model - radius known
N1.Rob <- InfRobModel(center = N1, neighbor = ContNeighborhood(radius = 0.2))

## optimally robust IC
IC1 <- optIC(model = N1.Rob, risk = asMSE())

## optimally robust estimator
est1 <- oneStepEstimator(chem, IC1, est0)

## radius unknown r in [0.05, 0.5]
IC2 <- radiusMinimaxIC(L2Fam=N1, neighbor=ContNeighborhood(),risk=asMSE(), 
                       loRad=0.05, upRad=0.5) 

## radius-minimax estimator
est2 <- oneStepEstimator(chem, IC2, est0)

## special case: normal location and scale
## between 1 and 10 % gross errors
require(RobLox)
est3 <- roblox(chem, eps.lower = 0.01, eps.upper = 0.1, initial.est = "med")

## estimator comparison
c(mean(chem), sd(chem)) # classical optimal
c(median(chem), mad(chem))
est0 # Kolmogorov(-Smirnov) minimum distance estimator
est1 # radius r = 0.2
est2 # radius r in [0.05, 0.5]
est3 # gross error eps=r/sqrt(n) in [0.01, 0.1]

require(RColorBrewer)
mypalette <- brewer.pal(4,"Dark2")
palette(mypalette)

#par(mfrow=c(1,2))
#plot(chem)
plot(chem[-17])
abline(h = mean(chem), col = 1)
abline(h = mean(chem) - sd(chem), col = 1)
abline(h = mean(chem) + sd(chem), col = 1)
abline(h = median(chem), col = 2)
abline(h = median(chem) - mad(chem), col = 2, lty = 2)
abline(h = median(chem) + mad(chem), col = 2, lty = 2)
abline(h = est0$mean, col = 3)
abline(h = est0$mean - est0$sd, col = 3, lty = 2)
abline(h = est0$mean + est0$sd, col = 3, lty = 2)
abline(h = est3$mean, col = 4)
abline(h = est3$mean - est3$sd, col = 4, lty = 2)
abline(h = est3$mean + est3$sd, col = 4, lty = 2)
legend(1, max(chem[-17]), legend = c("mean", "median", "ksMD", "optRob"), fill = 1:4)