robest: Optimally robust estimation

In ROptEst: Optimally Robust Estimation

Optimally robust estimation

Description

Function to compute optimally robust estimates for L2-differentiable parametric families via k-step construction.

Usage

robest(x, L2Fam,  fsCor = 1, risk = asMSE(), steps = 1L,
        verbose = NULL, OptOrIter = "iterate", nbCtrl = gennbCtrl(),
        startCtrl = genstartCtrl(), startICCtrl = genstartICCtrl(),
        kStepCtrl = genkStepCtrl(), na.rm = TRUE, ..., debug = FALSE,
        withTimings = FALSE, diagnostic = FALSE)

Arguments

x	sample
L2Fam	object of class "L2ParamFamily"
fsCor	positive real: factor used to correct the neighborhood radius; see details.
risk	object of class "RiskType"
steps	positive integer: number of steps used for k-steps construction
verbose	logical: if TRUE, some messages are printed
OptOrIter	character; which method to be used for determining Lagrange multipliers A and a: if (partially) matched to "optimize", getLagrangeMultByOptim is used; otherwise: by default, or if matched to "iterate" or to "doubleiterate", getLagrangeMultByIter is used. More specifically, when using getLagrangeMultByIter, and if argument risk is of class "asGRisk", by default and if matched to "iterate" we use only one (inner) iteration, if matched to "doubleiterate" we use up to Maxiter (inner) iterations.
nbCtrl	a list specifying input concerning the used neighborhood; to be generated by a respective call to gennbCtrl.
startCtrl	a list specifying input concerning the used starting estimator; to be generated by a respective call to genstartCtrl.
startICCtrl	a list specifying input concerning the call to getStartIC which returns the starting influence curve; to be generated by a respective call to genstartICCtrl.
kStepCtrl	a list specifying input concerning the used variant of a kstepEstimator; to be generated by a respective call to genkStepCtrl.
na.rm	logical: if TRUE, the estimator is evaluated at complete.cases(x).
...	further arguments
debug	logical: if TRUE, only the respective calls within the function are generated for debugging purposes.
withTimings	logical: if TRUE, separate (and aggregate) timings for the three steps evaluating the starting value, finding the starting influence curve, and evaluating the k-step estimator is issued.
diagnostic	logical; if TRUE, diagnostic information on the performed integrations is gathered and shipped out as attributes kStepDiagnostic (for the kStepEstimator-step) and diagnostic for the remaining steps of the return value of robest.

Details

A new, more structured interface to the former function roptest. For details, see this function.

In some respects this functions allows for more granular arguments, in the sense that the different steps (a) computation of the inital estimator, resp. (a') in case initial.est is missing computation of the initial MDE, (b) computation of the optimal IC and (c) computation of the k-step estimator each can have individial arguments E.arglist to be passed on to calls to expectation operator E within each step.

These different arguments are passed through the input generating functions genstartCtrl, genstartICCtrl, and kStepCtrl

Diagnostics on the involved integrations are available if argument diagnostic is TRUE. Then there are attributes diagnostic and kStepDiagnostic attached to the return value, which may be inspected and assessed through showDiagnostic and getDiagnostic.

Value

Object of class "kStepEstimate". In addition, it has an attribute "timings" where computation time is stored.

Author(s)

Matthias Kohl [e-mail],
Peter Ruckdeschel [e-mail]

See Also

roblox, L2ParamFamily-class UncondNeighborhood-class, RiskType-class

Examples

## Don't test to reduce check time on CRAN

#############################
## 1. Binomial data
#############################
## generate a sample of contaminated data
set.seed(123)
ind <- rbinom(100, size=1, prob=0.05)
x <- rbinom(100, size=25, prob=(1-ind)*0.25 + ind*0.9)

## Family
BF <- BinomFamily(size = 25)
## ML-estimate
MLest <- MLEstimator(x, BF)
estimate(MLest)
confint(MLest)

## compute optimally robust estimator (known contamination)
nb <- gennbCtrl(eps=0.05)
robest1 <- robest(x, BF, nbCtrl = nb, steps = 3)
estimate(robest1)

confint(robest1, method = symmetricBias())
## neglecting bias
confint(robest1)
plot(pIC(robest1))
tmp <- qqplot(x, robest1, cex.pch=1.5, exp.cex2.pch = -.25,
              exp.fadcol.pch = .55, jit.fac=.9)

## compute optimally robust estimator (unknown contamination)
nb2 <- gennbCtrl(eps.lower = 0, eps.upper = 0.2)
robest2 <- robest(x, BF, nbCtrl = nb2, steps = 3)
estimate(robest2)
confint(robest2, method = symmetricBias())
plot(pIC(robest2))

## total variation neighborhoods (known deviation)
nb3 <- gennbCtrl(eps = 0.025, neighbor = TotalVarNeighborhood())
robest3 <- robest(x, BF, nbCtrl = nb3, steps = 3)
estimate(robest3)
confint(robest3, method = symmetricBias())
plot(pIC(robest3))

## total variation neighborhoods (unknown deviation)
nb4 <- gennbCtrl(eps.lower = 0, eps.upper = 0.1,
                 neighbor = TotalVarNeighborhood())
robest3 <- robest(x, BF, nbCtrl = nb4, steps = 3)
robest4 <- robest(x, BinomFamily(size = 25), nbCtrl = nb4, steps = 3)
estimate(robest4)
confint(robest4, method = symmetricBias())
plot(pIC(robest4))


#############################
## 2. Poisson data
#############################
## Example: Rutherford-Geiger (1910); cf. Feller~(1968), Section VI.7 (a)
x <- c(rep(0, 57), rep(1, 203), rep(2, 383), rep(3, 525), rep(4, 532), 
       rep(5, 408), rep(6, 273), rep(7, 139), rep(8, 45), rep(9, 27), 
       rep(10, 10), rep(11, 4), rep(12, 0), rep(13, 1), rep(14, 1))

## Family
PF <- PoisFamily()

## ML-estimate
MLest <- MLEstimator(x, PF)
estimate(MLest)
confint(MLest)

## compute optimally robust estimator (unknown contamination)
nb1 <- gennbCtrl(eps.upper = 0.1)
robest <- robest(x, PF, nbCtrl = nb1, steps = 3)
estimate(robest)

confint(robest, symmetricBias())
plot(pIC(robest))
tmp <- qqplot(x, robest, cex.pch=1.5, exp.cex2.pch = -.25,
              exp.fadcol.pch = .55, jit.fac=.9)
 
## total variation neighborhoods (unknown deviation)
nb2 <- gennbCtrl(eps.upper = 0.05, neighbor = TotalVarNeighborhood())
robest1 <- robest(x, PF, nbCtrl = nb2, steps = 3)
estimate(robest1)
confint(robest1, symmetricBias())
plot(pIC(robest1))


#############################
## 3. Normal (Gaussian) location and scale
#############################
## 24 determinations of copper in wholemeal flour
library(MASS)
data(chem)
plot(chem, main = "copper in wholemeal flour", pch = 20)

## Family
NF <- NormLocationScaleFamily()
## ML-estimate
MLest <- MLEstimator(chem, NF)
estimate(MLest)
confint(MLest)

## Don't run to reduce check time on CRAN
## Not run: 
## compute optimally robust estimator (known contamination)
## takes some time -> you can use package RobLox for normal 
## location and scale which is optimized for speed
nb1 <- gennbCtrl(eps = 0.05)
robEst <- robest(chem, NF, nbCtrl = nb1, steps = 3)
estimate.call(robEst)
attr(robEst,"timings")
estimate(robest)

confint(robest, symmetricBias())
plot(pIC(robest))
## plot of relative and absolute information; cf. Kohl (2005)
infoPlot(pIC(robest))

tmp <- qqplot(chem, robest, cex.pch=1.5, exp.cex2.pch = -.25,
              exp.fadcol.pch = .55, withLab = TRUE, which.Order=1:4,
              exp.cex2.lbl = .12,exp.fadcol.lbl = .45,
              nosym.pCI = TRUE, adj.lbl=c(1.7,.2),
              exact.pCI = FALSE, log ="xy")
             
## finite-sample correction
if(require(RobLox)){
    n <- length(chem)
    r <- 0.05*sqrt(n)
    r.fi <- finiteSampleCorrection(n = n, r = r)
    fsCor0 <- r.fi/r
    nb1 <- gennbCtrl(eps = 0.05)
    robest <- robest(chem, NF, nbCtrl = nb1, fsCor = fsCor0, steps = 3)
    estimate(robest)
}

## compute optimally robust estimator (unknown contamination)
## takes some time -> use package RobLox!
nb2 <- gennbCtrl(eps.lower = 0.05, eps.upper = 0.1)
robest1 <- robest(chem, NF, nbCtrl = nb2, steps = 3)
estimate(robest1)
confint(robest1, symmetricBias())
plot(pIC(robest1))
## plot of relative and absolute information; cf. Kohl (2005)
infoPlot(pIC(robest1))

## End(Not run)

ROptEst documentation built on Nov. 17, 2022, 1:06 a.m.