[Robast-commits] r132 - in branches/robast-0.6/pkg: RandVar/inst/doc RobAStBase RobAStBase/R RobAStBase/man RobLox/R RobLox/man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Jul 31 17:07:02 CEST 2008
Author: stamats
Date: 2008-07-31 17:07:01 +0200 (Thu, 31 Jul 2008)
New Revision: 132
Added:
branches/robast-0.6/pkg/RobAStBase/R/RobAStBaseOptions.R
branches/robast-0.6/pkg/RobAStBase/R/kStepEstimator.R
branches/robast-0.6/pkg/RobAStBase/man/RobAStBaseOptions.Rd
branches/robast-0.6/pkg/RobAStBase/man/kStepEstimator.Rd
Modified:
branches/robast-0.6/pkg/RandVar/inst/doc/RandVar.pdf
branches/robast-0.6/pkg/RobAStBase/DESCRIPTION
branches/robast-0.6/pkg/RobAStBase/NAMESPACE
branches/robast-0.6/pkg/RobAStBase/R/AllClass.R
branches/robast-0.6/pkg/RobAStBase/R/AllGeneric.R
branches/robast-0.6/pkg/RobAStBase/R/ContIC.R
branches/robast-0.6/pkg/RobAStBase/R/IC.R
branches/robast-0.6/pkg/RobAStBase/R/TotalVarIC.R
branches/robast-0.6/pkg/RobAStBase/R/bALEstimate.R
branches/robast-0.6/pkg/RobAStBase/R/getBiasIC.R
branches/robast-0.6/pkg/RobAStBase/R/getRiskIC.R
branches/robast-0.6/pkg/RobAStBase/R/locMEstimator.R
branches/robast-0.6/pkg/RobAStBase/R/oneStepEstimator.R
branches/robast-0.6/pkg/RobAStBase/R/optIC.R
branches/robast-0.6/pkg/RobAStBase/man/ALEstimate-class.Rd
branches/robast-0.6/pkg/RobAStBase/man/BdStWeight-class.Rd
branches/robast-0.6/pkg/RobAStBase/man/BoundedWeight-class.Rd
branches/robast-0.6/pkg/RobAStBase/man/ContIC-class.Rd
branches/robast-0.6/pkg/RobAStBase/man/ContIC.Rd
branches/robast-0.6/pkg/RobAStBase/man/HampIC-class.Rd
branches/robast-0.6/pkg/RobAStBase/man/HampelWeight-class.Rd
branches/robast-0.6/pkg/RobAStBase/man/IC-class.Rd
branches/robast-0.6/pkg/RobAStBase/man/IC.Rd
branches/robast-0.6/pkg/RobAStBase/man/InfluenceCurve-class.Rd
branches/robast-0.6/pkg/RobAStBase/man/InfluenceCurve.Rd
branches/robast-0.6/pkg/RobAStBase/man/MEstimate-class.Rd
branches/robast-0.6/pkg/RobAStBase/man/RobAStControl-class.Rd
branches/robast-0.6/pkg/RobAStBase/man/RobWeight-class.Rd
branches/robast-0.6/pkg/RobAStBase/man/TotalVarIC-class.Rd
branches/robast-0.6/pkg/RobAStBase/man/TotalVarIC.Rd
branches/robast-0.6/pkg/RobAStBase/man/checkIC.Rd
branches/robast-0.6/pkg/RobAStBase/man/comparePlot.Rd
branches/robast-0.6/pkg/RobAStBase/man/evalIC.Rd
branches/robast-0.6/pkg/RobAStBase/man/generateIC.Rd
branches/robast-0.6/pkg/RobAStBase/man/generateICfct.Rd
branches/robast-0.6/pkg/RobAStBase/man/getBiasIC.Rd
branches/robast-0.6/pkg/RobAStBase/man/getRiskIC.Rd
branches/robast-0.6/pkg/RobAStBase/man/getweight.Rd
branches/robast-0.6/pkg/RobAStBase/man/infoPlot.Rd
branches/robast-0.6/pkg/RobAStBase/man/kStepEstimate-class.Rd
branches/robast-0.6/pkg/RobAStBase/man/locMEstimator.Rd
branches/robast-0.6/pkg/RobAStBase/man/makeIC-methods.Rd
branches/robast-0.6/pkg/RobAStBase/man/oneStepEstimator.Rd
branches/robast-0.6/pkg/RobAStBase/man/optIC.Rd
branches/robast-0.6/pkg/RobLox/R/rlOptIC.R
branches/robast-0.6/pkg/RobLox/R/rlsOptIC_AL.R
branches/robast-0.6/pkg/RobLox/R/roblox.R
branches/robast-0.6/pkg/RobLox/R/rowRoblox.R
branches/robast-0.6/pkg/RobLox/R/rsOptIC.R
branches/robast-0.6/pkg/RobLox/man/rlOptIC.Rd
branches/robast-0.6/pkg/RobLox/man/rlsOptIC.AL.Rd
branches/robast-0.6/pkg/RobLox/man/rsOptIC.Rd
Log:
major updates towards k-step estimation
Modified: branches/robast-0.6/pkg/RandVar/inst/doc/RandVar.pdf
===================================================================
(Binary files differ)
Modified: branches/robast-0.6/pkg/RobAStBase/DESCRIPTION
===================================================================
--- branches/robast-0.6/pkg/RobAStBase/DESCRIPTION 2008-07-28 12:27:14 UTC (rev 131)
+++ branches/robast-0.6/pkg/RobAStBase/DESCRIPTION 2008-07-31 15:07:01 UTC (rev 132)
@@ -1,6 +1,6 @@
Package: RobAStBase
Version: 0.1.0
-Date: 2008-07-21
+Date: 2008-07-30
Title: Robust Asymptotic Statistics
Description: Base S4-classes and functions for robust asymptotic statistics.
Depends: R(>= 2.6.0), methods, distr(>= 2.0), distrEx(>= 2.0), distrMod(>= 2.0), RandVar(>= 0.6.2)
@@ -9,4 +9,3 @@
LazyLoad: yes
License: GPL version 2 or later
URL: http://robast.r-forge.r-project.org/
-Packaged: Thu Jan 3 20:00:08 2008; btm722
Modified: branches/robast-0.6/pkg/RobAStBase/NAMESPACE
===================================================================
--- branches/robast-0.6/pkg/RobAStBase/NAMESPACE 2008-07-28 12:27:14 UTC (rev 131)
+++ branches/robast-0.6/pkg/RobAStBase/NAMESPACE 2008-07-31 15:07:01 UTC (rev 132)
@@ -27,7 +27,8 @@
exportMethods("Curve",
"Risks", "Risks<-", "addRisk<-",
"Infos", "Infos<-", "addInfo<-",
- "CallL2Fam", "CallL2Fam<-",
+ "CallL2Fam", "CallL2Fam<-",
+ "modifyIC",
"generateIC",
"checkIC",
"evalIC",
@@ -39,7 +40,8 @@
"clipLo", "clipLo<-",
"clipUp", "clipUp<-",
"optIC")
-exportMethods("oneStepEstimator",
+exportMethods("oneStepEstimator",
+ "kStepEstimator",
"locMEstimator")
exportMethods("weight", "weight<-",
"getweight",
@@ -49,8 +51,11 @@
exportMethods("getRiskIC")
exportMethods("getBiasIC")
exportMethods("comparePlot")
-exportMethods("pIC", "steps", "Mroot")
+exportMethods("pIC", "asbias",
+ "steps",
+ "Mroot")
export("ContNeighborhood", "TotalVarNeighborhood")
export("FixRobModel", "InfRobModel")
export("InfluenceCurve", "IC", "ContIC", "TotalVarIC")
export(".eq", ".getDistr")
+export("RobAStBaseOptions", "getRobAStBaseOption")
Modified: branches/robast-0.6/pkg/RobAStBase/R/AllClass.R
===================================================================
--- branches/robast-0.6/pkg/RobAStBase/R/AllClass.R 2008-07-28 12:27:14 UTC (rev 131)
+++ branches/robast-0.6/pkg/RobAStBase/R/AllClass.R 2008-07-31 15:07:01 UTC (rev 132)
@@ -5,6 +5,10 @@
require("distrMod", character = TRUE, quietly = TRUE)
require("RandVar", character = TRUE, quietly = TRUE)
}
+.onAttach <- function(library, pkg){
+ unlockBinding(".RobAStBaseOptions", asNamespace("RobAStBase"))
+ invisible()
+}
## neighborhood
setClass("Neighborhood",
@@ -85,14 +89,16 @@
else TRUE
})
## partial incluence curve
-setClass("IC", representation(CallL2Fam = "call"),
+setClass("IC", representation(CallL2Fam = "call",
+ modifyIC = "OptionalFunction"),
prototype(name = "square integrable (partial) influence curve",
Curve = EuclRandVarList(RealRandVariable(Map = list(function(x){x}),
Domain = Reals())),
Risks = list(),
Infos = matrix(c(character(0),character(0)), ncol=2,
dimnames=list(character(0), c("method", "message"))),
- CallL2Fam = call("L2ParamFamily")),
+ CallL2Fam = call("L2ParamFamily"),
+ modifyIC = NULL),
contains = "InfluenceCurve",
validity = function(object){
L2Fam <- eval(object at CallL2Fam)
@@ -119,6 +125,7 @@
Infos = matrix(c(character(0),character(0)), ncol=2,
dimnames=list(character(0), c("method", "message"))),
CallL2Fam = call("L2ParamFamily"),
+ modifyIC = NULL,
stand = as.matrix(1),
lowerCase = NULL,
neighborRadius = 0,
@@ -147,6 +154,7 @@
Infos = matrix(c(character(0),character(0)), ncol=2,
dimnames=list(character(0), c("method", "message"))),
CallL2Fam = call("L2ParamFamily"),
+ modifyIC = NULL,
clip = Inf, cent = 0, stand = as.matrix(1),
lowerCase = NULL,
neighborRadius = 0, weight = new("HampelWeight"),
@@ -180,9 +188,11 @@
Infos = matrix(c(character(0),character(0)), ncol=2,
dimnames=list(character(0), c("method", "message"))),
CallL2Fam = call("L2ParamFamily"),
+ modifyIC = NULL,
clipLo = -Inf, clipUp = Inf, stand = as.matrix(1),
lowerCase = NULL,
- neighborRadius = 0, weight = new("BdStWeight")),
+ neighborRadius = 0, weight = new("BdStWeight"),
+ biastype = symmetricBias(), NormType = NormType()),
contains = "HampIC",
validity = function(object){
if(any(object at neighborRadius < 0)) # radius vector?!
@@ -202,10 +212,15 @@
## ALEstimate
setClassUnion("OptionalInfluenceCurve", c("InfluenceCurve", "NULL"))
setClass("ALEstimate",
- representation(pIC = "OptionalInfluenceCurve"),
+ representation(pIC = "OptionalInfluenceCurve",
+ asbias = "OptionalNumeric"),
prototype(name = "Asymptotically linear estimate",
estimate = numeric(0),
+ samplesize = numeric(0),
+ asvar = NULL,
+ asbias = NULL,
pIC = NULL,
+ nuis.idx = NULL,
Infos = matrix(c(character(0),character(0)), ncol=2,
dimnames=list(character(0), c("method", "message")))),
contains = "Estimate")
@@ -213,7 +228,11 @@
representation(steps = "integer"),
prototype(name = "k-step estimate",
estimate = numeric(0),
+ samplesize = numeric(0),
+ asvar = NULL,
+ asbias = NULL,
pIC = NULL,
+ nuis.idx = NULL,
steps = integer(0),
Infos = matrix(c(character(0),character(0)), ncol=2,
dimnames=list(character(0), c("method", "message")))),
@@ -222,7 +241,11 @@
representation(Mroot = "numeric"),
prototype(name = "M estimate",
estimate = numeric(0),
+ samplesize = numeric(0),
+ asvar = NULL,
+ asbias = NULL,
pIC = NULL,
+ nuis.idx = NULL,
Mroot = numeric(0),
Infos = matrix(c(character(0),character(0)), ncol=2,
dimnames=list(character(0), c("method", "message")))),
Modified: branches/robast-0.6/pkg/RobAStBase/R/AllGeneric.R
===================================================================
--- branches/robast-0.6/pkg/RobAStBase/R/AllGeneric.R 2008-07-28 12:27:14 UTC (rev 131)
+++ branches/robast-0.6/pkg/RobAStBase/R/AllGeneric.R 2008-07-31 15:07:01 UTC (rev 132)
@@ -87,8 +87,12 @@
}
if(!isGeneric("oneStepEstimator")){
setGeneric("oneStepEstimator",
- function(x, IC, start) standardGeneric("oneStepEstimator"))
+ function(x, IC, start, ...) standardGeneric("oneStepEstimator"))
}
+if(!isGeneric("kStepEstimator")){
+ setGeneric("kStepEstimator",
+ function(x, IC, start, ...) standardGeneric("kStepEstimator"))
+}
if(!isGeneric("locMEstimator")){
setGeneric("locMEstimator", function(x, IC, ...) standardGeneric("locMEstimator"))
}
@@ -106,7 +110,7 @@
}
if(!isGeneric("weight<-")){
setGeneric("weight<-",
- function(object, value, ...) standardGeneric("weight<-"))
+ function(object, value) standardGeneric("weight<-"))
}
if(!isGeneric("clip")){
setGeneric("clip",
@@ -163,9 +167,15 @@
if(!isGeneric("pIC")){
setGeneric("pIC", function(object) standardGeneric("pIC"))
}
+if(!isGeneric("asbias")){
+ setGeneric("asbias", function(object) standardGeneric("asbias"))
+}
if(!isGeneric("steps")){
setGeneric("steps", function(object) standardGeneric("steps"))
}
if(!isGeneric("Mroot")){
setGeneric("Mroot", function(object) standardGeneric("Mroot"))
}
+if(!isGeneric("modifyIC")){
+ setGeneric("modifyIC", function(object) standardGeneric("modifyIC"))
+}
Modified: branches/robast-0.6/pkg/RobAStBase/R/ContIC.R
===================================================================
--- branches/robast-0.6/pkg/RobAStBase/R/ContIC.R 2008-07-28 12:27:14 UTC (rev 131)
+++ branches/robast-0.6/pkg/RobAStBase/R/ContIC.R 2008-07-31 15:07:01 UTC (rev 132)
@@ -4,7 +4,8 @@
Domain = Reals())),
Risks, Infos, clip = Inf, cent = 0, stand = as.matrix(1),
lowerCase = NULL, neighborRadius = 0, w = new("HampelWeight"),
- normtype = NormType(), biastype = symmetricBias()){
+ normtype = NormType(), biastype = symmetricBias(),
+ modifyIC = NULL){
if(missing(name))
name <- "IC of contamination type"
if(missing(Risks))
@@ -40,6 +41,7 @@
contIC at weight <- w
contIC at biastype <- biastype
contIC at normtype <- normtype
+ contIC at modifyIC <- modifyIC
return(contIC)
# return(new("ContIC", name = name, Curve = Curve, Risks = Risks, Infos = Infos,
@@ -56,9 +58,8 @@
a <- res$a
b <- res$b
d <- res$d
- normtype <- res$normtype
- biastype <- res$biastype
- if(is.null(res$w)) res$w <- new("HampelWeight")
+ normtype <- res$normtype
+ biastype <- res$biastype
w <- res$w
return(ContIC(
name = "IC of contamination type",
@@ -83,6 +84,7 @@
lowerCase = d,
w = w,
neighborRadius = neighbor at radius,
+ modifyIC = res$modifyIC,
normtype = normtype,
biastype = biastype,
Risks = res$risk,
@@ -94,6 +96,7 @@
setMethod("clip", "ContIC", function(object) object at clip)
setMethod("cent", "ContIC", function(object) object at cent)
+setMethod("neighbor", "ContIC", function(object) ContNeighborhood(radius = object at neighborRadius) )
## replace methods
setReplaceMethod("clip", "ContIC",
@@ -107,7 +110,8 @@
normW = object at normtype)
res <- list(A = object at stand, a = object at cent, b = value, d = object at lowerCase,
risk = object at Risks, info = object at Infos, w = w,
- normtype = object at normtype, biastype = object at biastype)
+ normtype = object at normtype, biastype = object at biastype,
+ modifyIC = object at modifyIC)
object <- generateIC(neighbor = ContNeighborhood(radius = object at neighborRadius),
L2Fam = L2Fam, res = res)
addInfo(object) <- c("clip<-", "The clipping bound has been changed")
@@ -125,7 +129,8 @@
normW = object at normtype)
res <- list(A = object at stand, a = value, b = object at clip, d = object at lowerCase,
risk = object at Risks, info = object at Infos, w = w,
- normtype = object at normtype, biastype = object at biastype)
+ normtype = object at normtype, biastype = object at biastype,
+ modifyIC = object at modifyIC)
object <- generateIC(neighbor = ContNeighborhood(radius = object at neighborRadius),
L2Fam = L2Fam, res = res)
addInfo(object) <- c("cent<-", "The centering constant has been changed")
@@ -143,7 +148,8 @@
normW = object at normtype)
res <- list(A = value, a = object at cent, b = object at clip, d = object at lowerCase,
risk = object at Risks, info = object at Infos, w = w,
- normtype = object at normtype, biastype = object at biastype)
+ normtype = object at normtype, biastype = object at biastype,
+ modifyIC = object at modifyIC)
object <- generateIC(neighbor = ContNeighborhood(radius = object at neighborRadius),
L2Fam = L2Fam, res = res)
addInfo(object) <- c("stand<-", "The standardizing matrix has been changed")
@@ -156,7 +162,8 @@
L2Fam <- eval(object at CallL2Fam)
res <- list(A = object at stand, a = object at cent, b = object at clip, d = value,
risk = object at Risks, info = object at Infos, w = object at weight,
- normtype = object at normtype, biastype = object at biastype)
+ normtype = object at normtype, biastype = object at biastype,
+ modifyIC = object at modifyIC)
object <- generateIC(neighbor = ContNeighborhood(radius = object at neighborRadius),
L2Fam = L2Fam, res = res)
addInfo(object) <- c("lowerCase<-", "The slot 'lowerCase' has been changed")
@@ -168,7 +175,8 @@
L2Fam <- eval(value)
res <- list(A = object at stand, a = object at cent, b = object at clip, d = object at lowerCase,
risk = object at Risks, info = object at Infos, w = object at weight,
- normtype = object at normtype, biastype = object at biastype)
+ normtype = object at normtype, biastype = object at biastype,
+ modifyIC = object at modifyIC)
object <- generateIC(neighbor = ContNeighborhood(radius = object at neighborRadius),
L2Fam = L2Fam, res = res)
addInfo(object) <- c("CallL2Fam<-", "The slot 'CallL2Fam' has been changed")
Modified: branches/robast-0.6/pkg/RobAStBase/R/IC.R
===================================================================
--- branches/robast-0.6/pkg/RobAStBase/R/IC.R 2008-07-28 12:27:14 UTC (rev 131)
+++ branches/robast-0.6/pkg/RobAStBase/R/IC.R 2008-07-31 15:07:01 UTC (rev 132)
@@ -1,6 +1,7 @@
## generating function
IC <- function(name, Curve = EuclRandVarList(RealRandVariable(Map = list(function(x){x}),
- Domain = Reals())), Risks, Infos, CallL2Fam = call("L2ParamFamily")){
+ Domain = Reals())), Risks, Infos, CallL2Fam = call("L2ParamFamily"),
+ modifyIC = NULL){
if(missing(name))
name <- "square integrable (partial) influence curve"
if(missing(Risks))
@@ -32,12 +33,14 @@
IC1 at Risks <- Risks
IC1 at Infos <- Infos
IC1 at CallL2Fam <- CallL2Fam
+ IC1 at modifyIC <- modifyIC
return(IC1)
}
## access methods
setMethod("CallL2Fam", "IC", function(object) object at CallL2Fam)
+setMethod("modifyIC", "IC", function(object) object at modifyIC)
## replace methods
setReplaceMethod("CallL2Fam", "IC",
@@ -126,25 +129,42 @@
trafo <- trafo(L2Fam at param)
IC1 <- as(diag(dimension(IC at Curve)) %*% IC at Curve, "EuclRandVariable")
cent <- E(D1, IC1)
+ IC1 <- IC1 - cent
dims <- length(L2Fam at param)
if(dimension(Domain(IC at Curve[[1]])) != dims)
stop("Dimension of IC and parameter must be the equal")
-
+
L2deriv <- as(diag(dims) %*% L2Fam at L2deriv, "EuclRandVariable")
E1 <- matrix(E(L2Fam, IC1 %*% t(L2deriv)), dims, dims)
-
stand <- trafo %*% solve(E1)
- Y <- as(stand %*% L2Fam at L2deriv - cent, "EuclRandVariable")
- ICfct <- vector(mode = "list", length = dims)
- ICfct[[1]] <- function(x){Y(x)}
- return(IC(name = name(IC),
- Curve = EuclRandVarList(EuclRandVariable(Map = ICfct,
- Domain = Y at Domain,Range = Y at Range)),
- Risks=list(), Infos=matrix(c("IC<-",
- "generated by affine linear trafo to enforce consistency"), ncol=2,
- dimnames=list(character(0), c("method", "message"))),
- CallL2Fam = IC at CallL2Fam))
+ Y <- as(stand %*% IC1, "EuclRandVariable")
+ #ICfct <- vector(mode = "list", length = dims)
+ #ICfct[[1]] <- function(x){Y(x)}
+
+ modifyIC <- function(L2Fam, IC){ makeIC(IC, L2Fam) }
+
+ CallL2Fam <- call("L2ParamFamily",
+ name = L2Fam at name,
+ distribution = L2Fam at distribution,
+ distrSymm = L2Fam at distrSymm,
+ param = L2Fam at param,
+ props = L2Fam at props,
+ modifyParam = L2Fam at modifyParam,
+ L2deriv.fct = L2Fam at L2deriv.fct,
+ L2derivSymm = L2Fam at L2derivSymm,
+ L2derivDistr = L2Fam at L2derivDistr,
+ L2derivDistrSymm = L2Fam at L2derivDistrSymm,
+ FisherInfo.fct = L2Fam at FisherInfo.fct,
+ FisherInfo = L2Fam at FisherInfo)
+
+ return(IC(name = name(IC),
+ Curve = EuclRandVarList(Y),
+ Risks = list(),
+ Infos=matrix(c("IC<-",
+ "generated by affine linear trafo to enforce consistency"),
+ ncol=2, dimnames=list(character(0), c("method", "message"))),
+ CallL2Fam = CallL2Fam,
+ modifyIC = modifyIC))
})
-
Added: branches/robast-0.6/pkg/RobAStBase/R/RobAStBaseOptions.R
===================================================================
--- branches/robast-0.6/pkg/RobAStBase/R/RobAStBaseOptions.R (rev 0)
+++ branches/robast-0.6/pkg/RobAStBase/R/RobAStBaseOptions.R 2008-07-31 15:07:01 UTC (rev 132)
@@ -0,0 +1,26 @@
+.RobAStBaseOptions <- list(
+ kStepUseLast = FALSE
+)
+
+RobAStBaseOptions <- function(...) {
+ if (nargs() == 0) return(.RobAStBaseOptions)
+ current <- .RobAStBaseOptions
+ temp <- list(...)
+ if (length(temp) == 1 && is.null(names(temp))) {
+ arg <- temp[[1]]
+ switch(mode(arg),
+ list = temp <- arg,
+ character = return(.RobAStBaseOptions[arg]),
+ stop("invalid argument: ", sQuote(arg)))
+ }
+ if (length(temp) == 0) return(current)
+ n <- names(temp)
+ if (is.null(n)) stop("options must be given by name")
+ changed <- current[n]
+ current[n] <- temp
+ env <- if (sys.parent() == 0) asNamespace("RobAStBase") else parent.frame()
+ assign(".RobAStBaseOptions", current, envir = env)
+ invisible(current)
+}
+
+getRobAStBaseOption <- function(x) RobAStBaseOptions(x)[[1]]
Modified: branches/robast-0.6/pkg/RobAStBase/R/TotalVarIC.R
===================================================================
--- branches/robast-0.6/pkg/RobAStBase/R/TotalVarIC.R 2008-07-28 12:27:14 UTC (rev 131)
+++ branches/robast-0.6/pkg/RobAStBase/R/TotalVarIC.R 2008-07-31 15:07:01 UTC (rev 132)
@@ -2,7 +2,9 @@
TotalVarIC <- function(name, CallL2Fam = call("L2ParamFamily"),
Curve = EuclRandVarList(RealRandVariable(Map = c(function(x){x}), Domain = Reals())),
Risks, Infos, clipLo = -Inf, clipUp = Inf, stand = as.matrix(1),
- lowerCase = NULL, neighborRadius = 0, w = new("BdStWeight")){
+ lowerCase = NULL, neighborRadius = 0, w = new("BdStWeight"),
+ normtype = NormType(), biastype = symmetricBias(),
+ modifyIC = NULL){
if(missing(name))
name <- "IC of total variation type"
@@ -32,6 +34,9 @@
IC1 at lowerCase <- lowerCase
IC1 at neighborRadius <- neighborRadius
IC1 at weight <- w
+ IC1 at biastype <- biastype
+ IC1 at normtype <- normtype
+ IC1 at modifyIC <- modifyIC
return(IC1)
}
@@ -44,7 +49,6 @@
A <- res$A
clipLo <- sign(as.vector(A))*res$a
b <- res$b
- if(is.null(res$w)) res$w <- new("BdStWeight")
w <- res$w
ICfct <- vector(mode = "list", length = 1)
Y <- as(A %*% L2Fam at L2deriv, "EuclRandVariable")
@@ -75,7 +79,10 @@
stand = A,
lowerCase = res$d,
w = w,
- neighborRadius = neighbor at radius,
+ modifyIC = res$modifyIC,
+ normtype = res$normtype,
+ biastype = res$biastype,
+ neighborRadius = neighbor at radius,
Risks = res$risk,
Infos = matrix(res$info, ncol = 2,
dimnames = list(character(0), c("method", "message")))))
@@ -84,6 +91,7 @@
## Access methods
setMethod("clipLo", "TotalVarIC", function(object) object at clipLo)
setMethod("clipUp", "TotalVarIC", function(object) object at clipUp)
+setMethod("neighbor", "TotalVarIC", function(object) TotalVarNeighborhood(radius = object at neighborRadius) )
## Replace methods
setReplaceMethod("clipLo", "TotalVarIC",
@@ -97,7 +105,8 @@
normW = object at normtype)
res <- list(A = object at stand, a = value, b = object at clipUp-value,
d = object at lowerCase, risk = object at Risks, info = object at Infos,
- w = w, normtype = object at normtype, biastype = object at biastype)
+ w = w, normtype = object at normtype, biastype = object at biastype,
+ modifyIC = object at modifyIC)
object <- generateIC(neighbor = TotalVarNeighborhood(radius = object at neighborRadius),
L2Fam = L2Fam, res = res)
addInfo(object) <- c("clipLo<-", "The lower clipping bound has been changed")
@@ -115,7 +124,8 @@
normW = object at normtype)
res <- list(A = object at stand, a = object at clipLo, b = value-object at clipLo,
d = object at lowerCase, risk = object at Risks, info = object at Infos,
- w = w, normtype = object at normtype, biastype = object at biastype)
+ w = w, normtype = object at normtype, biastype = object at biastype,
+ modifyIC = object at modifyIC)
object <- generateIC(neighbor = TotalVarNeighborhood(radius = object at neighborRadius),
L2Fam = L2Fam, res = res)
addInfo(object) <- c("clipUp<-", "The upper clipping bound has been changed")
@@ -133,7 +143,8 @@
normW = object at normtype)
res <- list(A = value, a = object at clipLo, b = object at clipUp-object@clipLo,
d = object at lowerCase, risk = object at Risks, info = object at Infos,
- w = w, normtype = object at normtype, biastype = object at biastype)
+ w = w, normtype = object at normtype, biastype = object at biastype,
+ modifyIC = object at modifyIC)
object <- generateIC(neighbor = TotalVarNeighborhood(radius = object at neighborRadius),
L2Fam = L2Fam, res = res)
addInfo(object) <- c("stand<-", "The standardizing matrix has been changed")
@@ -146,7 +157,8 @@
L2Fam <- eval(object at CallL2Fam)
res <- list(A = object at stand, a = object at clipLo, b = object at clipUp-object@clipLo,
d = value, risk = object at Risks, info = object at Infos,
- w = object at weight, normtype = object at normtype, biastype = object at biastype)
+ w = object at weight, normtype = object at normtype, biastype = object at biastype,
+ modifyIC = object at modifyIC)
object <- generateIC(neighbor = TotalVarNeighborhood(radius = object at neighborRadius),
L2Fam = L2Fam, res = res)
addInfo(object) <- c("lowerCase<-", "The slot 'lowerCase' has been changed")
@@ -158,7 +170,8 @@
L2Fam <- eval(value)
res <- list(A = object at stand, a = object at clipLo, b = object at clipUp-object@clipLo,
d = object at lowerCase, risk = object at Risks, info = object at Infos,
- w = object at weight, normtype = object at normtype, biastype = object at biastype)
+ w = object at weight, normtype = object at normtype, biastype = object at biastype,
+ modifyIC = object at modifyIC)
object <- generateIC(neighbor = TotalVarNeighborhood(radius = object at neighborRadius),
L2Fam = L2Fam, res = res)
addInfo(object) <- c("CallL2Fam<-", "The slot 'CallL2Fam' has been changed")
Modified: branches/robast-0.6/pkg/RobAStBase/R/bALEstimate.R
===================================================================
--- branches/robast-0.6/pkg/RobAStBase/R/bALEstimate.R 2008-07-28 12:27:14 UTC (rev 131)
+++ branches/robast-0.6/pkg/RobAStBase/R/bALEstimate.R 2008-07-31 15:07:01 UTC (rev 132)
@@ -3,5 +3,6 @@
###############################################################################
setMethod("pIC", "ALEstimate", function(object) object at pIC)
+setMethod("asbias", "ALEstimate", function(object) object at asbias)
setMethod("steps", "kStepEstimate", function(object) object at steps)
setMethod("Mroot", "MEstimate", function(object) object at Mroot)
Modified: branches/robast-0.6/pkg/RobAStBase/R/getBiasIC.R
===================================================================
--- branches/robast-0.6/pkg/RobAStBase/R/getBiasIC.R 2008-07-28 12:27:14 UTC (rev 131)
+++ branches/robast-0.6/pkg/RobAStBase/R/getBiasIC.R 2008-07-31 15:07:01 UTC (rev 132)
@@ -8,8 +8,10 @@
numbeval = 1e5){
misF <- FALSE
- if(missing(L2Fam))
- {misF <- TRUE; L2Fam <- eval(IC at CallL2Fam)}
+ if(missing(L2Fam)){
+ misF <- TRUE
+ L2Fam <- eval(IC at CallL2Fam)
+ }
D1 <- L2Fam at distribution
if(dimension(Domain(IC at Curve[[1]])) != dimension(img(D1)))
stop("dimension of 'Domain' of 'Curve' != dimension of 'img' of 'distribution' of 'L2Fam'")
@@ -23,7 +25,7 @@
prec <- if(misF) checkIC(IC, out = FALSE) else
checkIC(IC, L2Fam, out = FALSE)
if(prec > tol)
- warning("The maximum deviation from the exact IC properties is", prec,
+ warning("The maximum deviation from the exact IC properties is ", prec,
"\nThis is larger than the specified 'tol' ",
"=> the result may be wrong")
return(list(asBias = list(distribution = .getDistr(L2Fam),
Modified: branches/robast-0.6/pkg/RobAStBase/R/getRiskIC.R
===================================================================
--- branches/robast-0.6/pkg/RobAStBase/R/getRiskIC.R 2008-07-28 12:27:14 UTC (rev 131)
+++ branches/robast-0.6/pkg/RobAStBase/R/getRiskIC.R 2008-07-31 15:07:01 UTC (rev 132)
@@ -25,7 +25,7 @@
prec <- checkIC(IC, L2Fam, out = FALSE)
if(prec > tol)
- warning("The maximum deviation from the exact IC properties is", prec,
+ warning("The maximum deviation from the exact IC properties is ", prec,
"\nThis is larger than the specified 'tol' ",
"=> the result may be wrong")
@@ -57,7 +57,7 @@
prec <- checkIC(IC, L2Fam, out = FALSE)
if(prec > tol)
- warning("The maximum deviation from the exact IC properties is", prec,
+ warning("The maximum deviation from the exact IC properties is ", prec,
"\nThis is larger than the specified 'tol' ",
"=> the result may be wrong")
@@ -113,7 +113,7 @@
prec <- checkIC(IC, L2Fam, out = FALSE)
if(prec > tol)
- warning("The maximum deviation from the exact IC properties is", prec,
+ warning("The maximum deviation from the exact IC properties is ", prec,
"\nThis is larger than the specified 'tol' ",
"=> the result may be wrong")
Added: branches/robast-0.6/pkg/RobAStBase/R/kStepEstimator.R
===================================================================
--- branches/robast-0.6/pkg/RobAStBase/R/kStepEstimator.R (rev 0)
+++ branches/robast-0.6/pkg/RobAStBase/R/kStepEstimator.R 2008-07-31 15:07:01 UTC (rev 132)
@@ -0,0 +1,430 @@
+###############################################################################
+## one-step estimator
+###############################################################################
+setMethod("kStepEstimator", signature(x = "numeric",
+ IC = "IC",
+ start = "numeric"),
+ function(x, IC, start, steps = 1L, useLast = getRobAStBaseOption("kStepUseLast")){
+ if(!is.integer(steps))
+ steps <- as.integer(steps)
+ if(steps < 1)
+ stop("steps needs to be a positive integer")
+
+ nrvalues <- dimension(IC at Curve)
+ if(is.list(start)) start <- unlist(start)
+ if(nrvalues != length(start))
+ stop("dimension of 'start' != dimension of 'Curve'")
+
+ res <- start + rowMeans(evalIC(IC, as.matrix(x)), na.rm = TRUE)
+
+ L2Fam <- eval(CallL2Fam(IC))
+ Infos <- matrix(c("kStepEstimator",
+ paste(steps, "-step estimate for ", name(L2Fam), sep = "")),
+ ncol = 2)
+ colnames(Infos) <- c("method", "message")
+ if(is(L2Fam, "L2GroupParamFamily")) useLast <- TRUE
+
+ if(steps == 1){
+ if(useLast && !is(modifyIC(IC), "NULL") ){
+ newParam <- param(L2Fam)
+ main(newParam) <- res
+ newL2Fam <- modifyModel(L2Fam, newParam)
+ IC <- modifyIC(IC)(newL2Fam, IC)
+ Infos <- rbind(Infos, c("oneStepEstimator",
+ "computation of IC, asvar and asbias via useLast = TRUE"))
+ }else{
+ if(useLast && is(modifyIC(IC), "NULL")){
+ warning("'useLast = TRUE' only possible if slot 'modifyIC' of 'IC'
+ is filled with some function!")
+ Infos <- rbind(Infos, c("oneStepEstimator",
+ "slot 'modifyIC' of 'IC' was not filled!"))
+ }
+ Infos <- rbind(Infos, c("oneStepEstimator",
+ "computation of IC, asvar and asbias via useLast = FALSE"))
+ }
+ if("asCov" %in% names(Risks(IC)))
+ asVar <- Risks(IC)$asCov
+ else
+ asVar <- getRiskIC(IC, risk = asCov())$asCov$value
+
+ if("asBias" %in% names(Risks(IC))){
+ asBias <- neighborRadius(IC)*Risks(IC)$asBias
+ }else{
+ if(is(IC, "HampIC")){
+ r <- neighborRadius(IC)
+ asBias <- r*getRiskIC(IC, risk = asBias(), neighbor = neighbor(IC))$asBias$value
+ }else{
+ asBias <- NULL
+ }
+ }
+ return(new("kStepEstimator",
+ name = paste(steps, "-step estimate", sep = ""),
+ estimate = res, samplesize = length(x), asvar = asVar,
+ asbias = asBias, pIC = IC, steps = 1L, Infos = Infos))
+ }else{
+ if(is(modifyIC(IC), "NULL"))
+ stop("slot 'modifyIC' of 'IC' is 'NULL'!")
+ for(i in 2:steps){
+ start <- res
+ newL2Fam <- eval(CallL2Fam(IC))
+ newParam <- param(newL2Fam)
+ main(newParam) <- start
+ newL2Fam <- modifyModel(newL2Fam, newParam)
+ IC <- modifyIC(IC)(newL2Fam, IC)
+ res <- start + rowMeans(evalIC(IC, as.matrix(x)), na.rm = TRUE)
+ }
+ if(useLast){
+ newL2Fam <- eval(CallL2Fam(IC))
+ newParam <- param(newL2Fam)
+ main(newParam) <- res
+ newL2Fam <- modifyModel(newL2Fam, newParam)
+ IC <- modifyIC(IC)(newL2Fam, IC)
+ Infos <- rbind(Infos, c("oneStepEstimator",
+ "computation of IC, asvar and asbias via useLast = TRUE"))
+ }else{
+ Infos <- rbind(Infos, c("oneStepEstimator",
+ "computation of IC, asvar and asbias via useLast = FALSE"))
+ }
+ if("asCov" %in% names(Risks(IC)))
+ asVar <- Risks(IC)$asCov
+ else
+ asVar <- getRiskIC(IC, risk = asCov())$asCov$value
+
+ if("asBias" %in% names(Risks(IC))){
+ asBias <- neighborRadius(IC)*Risks(IC)$asBias
+ }else{
+ if(is(IC, "HampIC")){
+ r <- neighborRadius(IC)
+ asBias <- r*getRiskIC(IC, risk = asBias(), neighbor = neighbor(IC))$asBias$value
+ }else{
+ asBias <- NULL
+ }
+ }
+ return(new("kStepEstimate",
+ name = paste(steps, "-step estimate", sep = ""),
+ estimate = res, samplesize = length(x), asvar = asVar,
+ asbias = asBias, pIC = IC, steps = steps, Infos = Infos))
+ }
+ })
+
+setMethod("kStepEstimator", signature(x = "matrix",
+ IC = "IC",
+ start = "numeric"),
+ function(x, IC, start, steps = 1, useLast = getRobAStBaseOption("kStepUseLast")){
+ if(!is.integer(steps))
+ steps <- as.integer(steps)
+ if(steps < 1)
+ stop("steps needs to be a positive integer")
+
+ nrvalues <- dimension(IC at Curve)
+ if(is.list(start)) start <- unlist(start)
+ if(nrvalues != length(start))
+ stop("dimension of 'start' != dimension of 'Curve'")
+ if(ncol(x) != IC at Curve[[1]]@Domain at dimension)
+ stop("'x' has wrong dimension")
+
+ res <- start + rowMeans(evalIC(IC, x), na.rm = TRUE)
+
+ L2Fam <- eval(CallL2Fam(IC))
+ Infos <- matrix(c("kStepEstimator",
+ paste(steps, "-step estimate for ", name(L2Fam), sep = "")),
+ ncol = 2)
+ colnames(Infos) <- c("method", "message")
+ if(is(L2Fam, "L2GroupParamFamily")) useLast <- TRUE
+
+ if(steps == 1){
+ if(useLast && !is(modifyIC(IC), "NULL") ){
+ newParam <- param(L2Fam)
+ main(newParam) <- res
+ newL2Fam <- modifyModel(L2Fam, newParam)
+ IC <- modifyIC(IC)(newL2Fam, IC)
+ Infos <- rbind(Infos, c("oneStepEstimator",
+ "computation of IC, asvar and asbias via useLast = TRUE"))
+ }else{
+ if(useLast && is(modifyIC(IC), "NULL")){
+ warning("'useLast = TRUE' only possible if slot 'modifyIC' of 'IC'
+ is filled with some function!")
+ Infos <- rbind(Infos, c("oneStepEstimator",
+ "slot 'modifyIC' of 'IC' was not filled!"))
+ }
+ Infos <- rbind(Infos, c("oneStepEstimator",
+ "computation of IC, asvar and asbias via useLast = FALSE"))
+ }
+ if("asCov" %in% names(Risks(IC)))
+ asVar <- Risks(IC)$asCov
+ else
+ asVar <- getRiskIC(IC, risk = asCov())$asCov$value
+
+ if("asBias" %in% names(Risks(IC))){
+ asBias <- neighborRadius(IC)*Risks(IC)$asBias
+ }else{
+ if(is(IC, "HampIC")){
+ r <- neighborRadius(IC)
+ asBias <- r*getRiskIC(IC, risk = asBias(), neighbor = neighbor(IC))$asBias$value
+ }else{
+ asBias <- NULL
+ }
+ }
+ return(new("kStepEstimator",
+ name = paste(steps, "-step estimate", sep = ""),
+ estimate = res, samplesize = ncol(x), asvar = asVar,
+ asbias = asBias, pIC = IC, steps = 1L, Infos = Infos))
+ }else{
+ if(is(modifyIC(IC), "NULL"))
+ stop("slot 'modifyIC' of 'IC' is 'NULL'!")
+ for(i in 2:steps){
+ start <- res
+ newL2Fam <- eval(CallL2Fam(IC))
+ newParam <- param(newL2Fam)
+ main(newParam) <- start
+ newL2Fam <- modifyModel(newL2Fam, newParam)
+ IC <- modifyIC(IC)(newL2Fam, IC)
+ res <- start + rowMeans(evalIC(IC, x), na.rm = TRUE)
+ }
+ if(useLast){
+ newL2Fam <- eval(CallL2Fam(IC))
+ newParam <- param(newL2Fam)
+ main(newParam) <- res
+ newL2Fam <- modifyModel(newL2Fam, newParam)
+ IC <- modifyIC(IC)(newL2Fam, IC)
+ Infos <- rbind(Infos, c("oneStepEstimator",
+ "computation of IC, asvar and asbias via useLast = TRUE"))
+ }else{
+ Infos <- rbind(Infos, c("oneStepEstimator",
+ "computation of IC, asvar and asbias via useLast = FALSE"))
+ }
+ if("asCov" %in% names(Risks(IC)))
+ asVar <- Risks(IC)$asCov
+ else
+ asVar <- getRiskIC(IC, risk = asCov())$asCov$value
+
+ if("asBias" %in% names(Risks(IC))){
+ asBias <- neighborRadius(IC)*Risks(IC)$asBias
+ }else{
+ if(is(IC, "HampIC")){
+ r <- neighborRadius(IC)
+ asBias <- r*getRiskIC(IC, risk = asBias(), neighbor = neighbor(IC))$asBias$value
+ }else{
+ asBias <- NULL
+ }
+ }
+ return(new("kStepEstimate",
+ name = paste(steps, "-step estimate", sep = ""),
+ estimate = res, samplesize = ncol(x), asvar = asVar,
+ asbias = asBias, pIC = IC, steps = steps, Infos = Infos))
+ }
+ })
+setMethod("kStepEstimator", signature(x = "numeric",
+ IC = "IC",
+ start = "Estimate"),
+ function(x, IC, start, steps = 1, useLast = getRobAStBaseOption("kStepUseLast")){
+ if(!is.integer(steps))
+ steps <- as.integer(steps)
+ if(steps < 1)
+ stop("steps needs to be a positive integer")
+
+ nrvalues <- dimension(IC at Curve)
+ start0 <- estimate(start)
+ if(is.list(start0)) start0 <- unlist(start0)
+ if(nrvalues != length(start0))
+ stop("dimension of slot 'estimate' of 'start' != dimension of 'Curve'")
+
+ res <- start0 + rowMeans(evalIC(IC, as.matrix(x)), na.rm = TRUE)
+
+ L2Fam <- eval(CallL2Fam(IC))
+ Infos <- matrix(c("kStepEstimator",
+ paste(steps, "-step estimate for ", name(L2Fam), sep = "")),
+ ncol = 2)
+ colnames(Infos) <- c("method", "message")
+ if(is(L2Fam, "L2GroupParamFamily")) useLast <- TRUE
+
+ if(steps == 1){
+ if(useLast && !is(modifyIC(IC), "NULL") ){
+ newParam <- param(L2Fam)
+ main(newParam) <- res
+ newL2Fam <- modifyModel(L2Fam, newParam)
+ IC <- modifyIC(IC)(newL2Fam, IC)
+ Infos <- rbind(Infos, c("oneStepEstimator",
+ "computation of IC, asvar and asbias via useLast = TRUE"))
+ }else{
+ if(useLast && is(modifyIC(IC), "NULL")){
+ warning("'useLast = TRUE' only possible if slot 'modifyIC' of 'IC'
+ is filled with some function!")
+ Infos <- rbind(Infos, c("oneStepEstimator",
+ "slot 'modifyIC' of 'IC' was not filled!"))
+ }
+ Infos <- rbind(Infos, c("oneStepEstimator",
+ "computation of IC, asvar and asbias via useLast = FALSE"))
+ }
+ if("asCov" %in% names(Risks(IC)))
+ asVar <- Risks(IC)$asCov
+ else
+ asVar <- getRiskIC(IC, risk = asCov())$asCov$value
+
+ if("asBias" %in% names(Risks(IC))){
+ asBias <- neighborRadius(IC)*Risks(IC)$asBias
+ }else{
+ if(is(IC, "HampIC")){
+ r <- neighborRadius(IC)
+ asBias <- r*getRiskIC(IC, risk = asBias(), neighbor = neighbor(IC))$asBias$value
+ }else{
+ asBias <- NULL
+ }
+ }
+ return(new("kStepEstimator",
+ name = paste(steps, "-step estimate", sep = ""),
+ estimate = res, samplesize = length(x), asvar = asVar,
+ asbias = asBias, pIC = IC, steps = 1L, Infos = Infos))
+ }else{
+ if(is(modifyIC(IC), "NULL"))
+ stop("slot 'modifyIC' of 'IC' is 'NULL'!")
+ for(i in 2:steps){
+ start0 <- res
+ newL2Fam <- eval(CallL2Fam(IC))
+ newParam <- param(newL2Fam)
+ main(newParam) <- start0
+ newL2Fam <- modifyModel(newL2Fam, newParam)
+ IC <- modifyIC(IC)(newL2Fam, IC)
+ res <- start0 + rowMeans(evalIC(IC, as.matrix(x)), na.rm = TRUE)
+ }
+ if(useLast){
+ newL2Fam <- eval(CallL2Fam(IC))
+ newParam <- param(newL2Fam)
+ main(newParam) <- res
+ newL2Fam <- modifyModel(newL2Fam, newParam)
+ IC <- modifyIC(IC)(newL2Fam, IC)
+ Infos <- rbind(Infos, c("oneStepEstimator",
+ "computation of IC, asvar and asbias via useLast = TRUE"))
+ }else{
+ Infos <- rbind(Infos, c("oneStepEstimator",
+ "computation of IC, asvar and asbias via useLast = FALSE"))
+ }
+ if("asCov" %in% names(Risks(IC)))
+ asVar <- Risks(IC)$asCov
+ else
+ asVar <- getRiskIC(IC, risk = asCov())$asCov$value
+
+ if("asBias" %in% names(Risks(IC))){
+ asBias <- neighborRadius(IC)*Risks(IC)$asBias
+ }else{
+ if(is(IC, "HampIC")){
+ r <- neighborRadius(IC)
+ asBias <- r*getRiskIC(IC, risk = asBias(), neighbor = neighbor(IC))$asBias$value
+ }else{
+ asBias <- NULL
+ }
+ }
+ return(new("kStepEstimate",
+ name = paste(steps, "-step estimate", sep = ""),
+ estimate = res, samplesize = length(x), asvar = asVar,
+ asbias = asBias, pIC = IC, steps = steps, Infos = Infos))
+ }
+ })
+setMethod("kStepEstimator", signature(x = "matrix",
+ IC = "IC",
+ start = "Estimate"),
+ function(x, IC, start, steps = 1, useLast = getRobAStBaseOption("kStepUseLast")){
+ if(!is.integer(steps))
+ steps <- as.integer(steps)
+ if(steps < 1)
+ stop("steps needs to be a positive integer")
+
+ nrvalues <- dimension(IC at Curve)
+ start0 <- estimate(start)
+ if(is.list(start0)) start0 <- unlist(start0)
+ if(nrvalues != length(start0))
+ stop("dimension of slot 'estimate' of 'start' != dimension of 'Curve'")
+ if(ncol(x) != IC at Curve[[1]]@Domain at dimension)
+ stop("'x' has wrong dimension")
+
+ res <- start0 + rowMeans(evalIC(IC, x), na.rm = TRUE)
+
+ L2Fam <- eval(CallL2Fam(IC))
+ Infos <- matrix(c("kStepEstimator",
+ paste(steps, "-step estimate for ", name(L2Fam), sep = "")),
+ ncol = 2)
+ colnames(Infos) <- c("method", "message")
+ if(is(L2Fam, "L2GroupParamFamily")) useLast <- TRUE
+
+ if(steps == 1){
+ if(useLast && !is(modifyIC(IC), "NULL") ){
+ newParam <- param(L2Fam)
+ main(newParam) <- res
+ newL2Fam <- modifyModel(L2Fam, newParam)
+ IC <- modifyIC(IC)(newL2Fam, IC)
+ Infos <- rbind(Infos, c("oneStepEstimator",
+ "computation of IC, asvar and asbias via useLast = TRUE"))
+ }else{
+ if(useLast && is(modifyIC(IC), "NULL")){
+ warning("'useLast = TRUE' only possible if slot 'modifyIC' of 'IC'
+ is filled with some function!")
+ Infos <- rbind(Infos, c("oneStepEstimator",
+ "slot 'modifyIC' of 'IC' was not filled!"))
+ }
+ Infos <- rbind(Infos, c("oneStepEstimator",
+ "computation of IC, asvar and asbias via useLast = FALSE"))
+ }
+ if("asCov" %in% names(Risks(IC)))
+ asVar <- Risks(IC)$asCov
+ else
+ asVar <- getRiskIC(IC, risk = asCov())$asCov$value
+
+ if("asBias" %in% names(Risks(IC))){
+ asBias <- neighborRadius(IC)*Risks(IC)$asBias
+ }else{
+ if(is(IC, "HampIC")){
+ r <- neighborRadius(IC)
+ asBias <- r*getRiskIC(IC, risk = asBias(), neighbor = neighbor(IC))$asBias$value
+ }else{
+ asBias <- NULL
+ }
+ }
+ return(new("kStepEstimator",
+ name = paste(steps, "-step estimate", sep = ""),
+ estimate = res, samplesize = ncol(x), asvar = asVar,
+ asbias = asBias, pIC = IC, steps = 1L, Infos = Infos))
+ }else{
+ if(is(modifyIC(IC), "NULL"))
+ stop("slot 'modifyIC' of 'IC' is 'NULL'!")
+ for(i in 2:steps){
+ start0 <- res
+ newL2Fam <- eval(CallL2Fam(IC))
+ newParam <- param(newL2Fam)
+ main(newParam) <- start0
+ newL2Fam <- modifyModel(newL2Fam, newParam)
+ IC <- modifyIC(IC)(newL2Fam, IC)
+ res <- start0 + rowMeans(evalIC(IC, x), na.rm = TRUE)
+ }
+ if(useLast){
+ newL2Fam <- eval(CallL2Fam(IC))
+ newParam <- param(newL2Fam)
+ main(newParam) <- res
+ newL2Fam <- modifyModel(newL2Fam, newParam)
+ IC <- modifyIC(IC)(newL2Fam, IC)
+ Infos <- rbind(Infos, c("oneStepEstimator",
+ "computation of IC, asvar and asbias via useLast = TRUE"))
+ }else{
+ Infos <- rbind(Infos, c("oneStepEstimator",
+ "computation of IC, asvar and asbias via useLast = FALSE"))
+ }
+ if("asCov" %in% names(Risks(IC)))
+ asVar <- Risks(IC)$asCov
+ else
+ asVar <- getRiskIC(IC, risk = asCov())$asCov$value
+
+ if("asBias" %in% names(Risks(IC))){
+ asBias <- neighborRadius(IC)*Risks(IC)$asBias
+ }else{
+ if(is(IC, "HampIC")){
+ r <- neighborRadius(IC)
+ asBias <- r*getRiskIC(IC, risk = asBias(), neighbor = neighbor(IC))$asBias$value
+ }else{
+ asBias <- NULL
+ }
+ }
+ return(new("kStepEstimate",
+ name = paste(steps, "-step estimate", sep = ""),
+ estimate = res, samplesize = ncol(x), asvar = asVar,
+ asbias = asBias, pIC = IC, steps = steps, Infos = Infos))
+ }
+ })
Modified: branches/robast-0.6/pkg/RobAStBase/R/locMEstimator.R
===================================================================
--- branches/robast-0.6/pkg/RobAStBase/R/locMEstimator.R 2008-07-28 12:27:14 UTC (rev 131)
+++ branches/robast-0.6/pkg/RobAStBase/R/locMEstimator.R 2008-07-31 15:07:01 UTC (rev 132)
@@ -11,16 +11,22 @@
}
res <- uniroot(f = mest, interval = c(min(x), max(x)),
tol = eps, x = x, IC = IC)
+
if(is(IC, "IC")){
L2Fam <- eval(CallL2Fam(IC))
Infos <- matrix(c("locMEstimator",
- paste("Location M estimate for", name(L2Fam))))
+ paste("Location M estimate for", name(L2Fam))),
+ ncol = 2)
colnames(Infos) <- c("method", "message")
+ asVar <- getRiskIC(IC, risk = asCov(), L2Fam = L2Fam)$asCov$value
+ asBias <- getRiskIC(IC, risk = asBias(), L2Fam = L2Fam)$asBias$value
}else{
- Infos <- matrix(c(character(0),character(0)), ncol=2,
- dimnames=list(character(0), c("method", "message")))
+ Infos <- matrix(c("locMEstimator", "Location M estimate"), ncol = 2)
+ colnames(Infos) <- c("method", "message")
+ asvar <- NULL
}
new("MEstimate", name = "Location M estimate", estimate = res$root,
- pIC = IC, Mroot = res$f.root, Infos = Infos)
+ pIC = IC, Mroot = res$f.root, Infos = Infos, samplesize = length(x),
+ asvar = asVar, asbias = asBias)
})
Modified: branches/robast-0.6/pkg/RobAStBase/R/oneStepEstimator.R
===================================================================
--- branches/robast-0.6/pkg/RobAStBase/R/oneStepEstimator.R 2008-07-28 12:27:14 UTC (rev 131)
+++ branches/robast-0.6/pkg/RobAStBase/R/oneStepEstimator.R 2008-07-31 15:07:01 UTC (rev 132)
@@ -4,7 +4,7 @@
setMethod("oneStepEstimator", signature(x = "numeric",
IC = "InfluenceCurve",
start = "numeric"),
- function(x, IC, start){
+ function(x, IC, start, useLast = getRobAStBaseOption("kStepUseLast")){
nrvalues <- dimension(IC at Curve)
if(is.list(start)) start <- unlist(start)
if(nrvalues != length(start))
@@ -12,12 +12,60 @@
res <- start + rowMeans(evalIC(IC, as.matrix(x)), na.rm = TRUE)
- return(res)
+ if(is(IC, "IC")){
+ L2Fam <- eval(CallL2Fam(IC))
+ Infos <- matrix(c("oneStepEstimator",
+ paste("1-step estimate for", name(L2Fam))),
+ ncol = 2)
+ colnames(Infos) <- c("method", "message")
+ if(is(L2Fam, "L2GroupParamFamily")) useLast <- TRUE
+ if(useLast && !is(modifyIC(IC), "NULL") ){
+ newParam <- param(L2Fam)
+ main(newParam) <- res
+ newL2Fam <- modifyModel(L2Fam, newParam)
+ IC <- modifyIC(IC)(newL2Fam, IC)
+ Infos <- rbind(Infos, c("oneStepEstimator",
+ "computation of IC, asvar and asbias via useLast = TRUE"))
+ }else{
+ if(useLast && is(modifyIC(IC), "NULL")){
+ warning("'useLast = TRUE' only possible if slot 'modifyIC' of 'IC'
+ is filled with some function!")
+ Infos <- rbind(Infos, c("oneStepEstimator",
+ "slot 'modifyIC' of 'IC' was not filled!"))
+ }
+ Infos <- rbind(Infos, c("oneStepEstimator",
+ "computation of IC, asvar and asbias via useLast = FALSE"))
+ }
+ if("asCov" %in% names(Risks(IC)))
+ asVar <- Risks(IC)$asCov
+ else
+ asVar <- getRiskIC(IC, risk = asCov())$asCov$value
+
+ if("asBias" %in% names(Risks(IC))){
+ asBias <- neighborRadius(IC)*Risks(IC)$asBias
+ }else{
+ if(is(IC, "HampIC")){
+ r <- neighborRadius(IC)
+ asBias <- r*getRiskIC(IC, risk = asBias(), neighbor = neighbor(IC))$asBias$value
+ }else{
+ asBias <- NULL
+ }
+ }
+ }else{
+ Infos <- matrix(c("oneStepEstimator", "1-step estimate"), ncol = 2)
+ colnames(Infos) <- c("method", "message")
+ asVar <- NULL
+ asBias <- NULL
+ }
+
+ new("kStepEstimate", name = "1-step estimate", estimate = res,
+ samplesize = length(x), asvar = asVar, asbias = asBias, pIC = IC,
+ steps = 1L, Infos = Infos)
})
setMethod("oneStepEstimator", signature(x = "matrix",
IC = "InfluenceCurve",
start = "numeric"),
- function(x, IC, start){
+ function(x, IC, start, useLast = getRobAStBaseOption("kStepUseLast")){
nrvalues <- dimension(IC at Curve)
if(is.list(start)) start <- unlist(start)
if(nrvalues != length(start))
@@ -27,12 +75,60 @@
res <- start + rowMeans(evalIC(IC, x), na.rm = TRUE)
- return(res)
+ if(is(IC, "IC")){
+ L2Fam <- eval(CallL2Fam(IC))
+ Infos <- matrix(c("oneStepEstimator",
+ paste("1-step estimate for", name(L2Fam))),
+ ncol = 2)
+ colnames(Infos) <- c("method", "message")
+ if(is(L2Fam, "L2GroupParamFamily")) useLast <- TRUE
+ if(useLast && !is(modifyIC(IC), "NULL") ){
+ newParam <- param(L2Fam)
+ main(newParam) <- res
+ newL2Fam <- modifyModel(L2Fam, newParam)
+ IC <- modifyIC(IC)(newL2Fam, IC)
+ Infos <- rbind(Infos, c("oneStepEstimator",
+ "computation of IC, asvar and asbias via useLast = TRUE"))
+ }else{
+ if(useLast && is(modifyIC(IC), "NULL")){
+ warning("'useLast = TRUE' only possible if slot 'modifyIC' of 'IC'
+ is filled with some function!")
+ Infos <- rbind(Infos, c("oneStepEstimator",
+ "slot 'modifyIC' of 'IC' was not filled!"))
+ }
+ Infos <- rbind(Infos, c("oneStepEstimator",
+ "computation of IC, asvar and asbias via useLast = FALSE"))
+ }
+ if("asCov" %in% names(Risks(IC)))
+ asVar <- Risks(IC)$asCov
+ else
+ asVar <- getRiskIC(IC, risk = asCov())$asCov$value
+
+ if("asBias" %in% names(Risks(IC))){
+ asBias <- neighborRadius(IC)*Risks(IC)$asBias
+ }else{
+ if(is(IC, "HampIC")){
+ r <- neighborRadius(IC)
+ asBias <- r*getRiskIC(IC, risk = asBias(), neighbor = neighbor(IC))$asBias$value
+ }else{
+ asBias <- NULL
+ }
+ }
+ }else{
+ Infos <- matrix(c("oneStepEstimator", "1-step estimate"), ncol = 2)
+ colnames(Infos) <- c("method", "message")
+ asVar <- NULL
+ asBias <- NULL
+ }
+
+ new("kStepEstimate", name = "1-step estimate", estimate = res,
+ samplesize = ncol(x), asvar = asVar, asbias = asBias, pIC = IC, steps = 1L,
+ Infos = Infos)
})
setMethod("oneStepEstimator", signature(x = "numeric",
IC = "InfluenceCurve",
start = "Estimate"),
- function(x, IC, start){
+ function(x, IC, start, useLast = getRobAStBaseOption("kStepUseLast")){
nrvalues <- dimension(IC at Curve)
start0 <- estimate(start)
if(is.list(start0)) start0 <- unlist(start0)
@@ -41,12 +137,60 @@
res <- start0 + rowMeans(evalIC(IC, as.matrix(x)), na.rm = TRUE)
- return(res)
+ if(is(IC, "IC")){
+ L2Fam <- eval(CallL2Fam(IC))
+ Infos <- matrix(c("oneStepEstimator",
+ paste("1-step estimate for", name(L2Fam))),
+ ncol = 2)
+ colnames(Infos) <- c("method", "message")
+ if(is(L2Fam, "L2GroupParamFamily")) useLast <- TRUE
+ if(useLast && !is(modifyIC(IC), "NULL") ){
+ newParam <- param(L2Fam)
+ main(newParam) <- res
+ newL2Fam <- modifyModel(L2Fam, newParam)
+ IC <- modifyIC(IC)(newL2Fam, IC)
+ Infos <- rbind(Infos, c("oneStepEstimator",
+ "computation of IC, asvar and asbias via useLast = TRUE"))
+ }else{
+ if(useLast && is(modifyIC(IC), "NULL")){
+ warning("'useLast = TRUE' only possible if slot 'modifyIC' of 'IC'
+ is filled with some function!")
+ Infos <- rbind(Infos, c("oneStepEstimator",
+ "slot 'modifyIC' of 'IC' was not filled!"))
+ }
+ Infos <- rbind(Infos, c("oneStepEstimator",
+ "computation of IC, asvar and asbias via useLast = FALSE"))
+ }
+ if("asCov" %in% names(Risks(IC)))
+ asVar <- Risks(IC)$asCov
+ else
+ asVar <- getRiskIC(IC, risk = asCov())$asCov$value
+
+ if("asBias" %in% names(Risks(IC))){
+ asBias <- neighborRadius(IC)*Risks(IC)$asBias
+ }else{
+ if(is(IC, "HampIC")){
+ r <- neighborRadius(IC)
+ asBias <- r*getRiskIC(IC, risk = asBias(), neighbor = neighbor(IC))$asBias$value
+ }else{
+ asBias <- NULL
+ }
+ }
+ }else{
+ Infos <- matrix(c("oneStepEstimator", "1-step estimate"), ncol = 2)
+ colnames(Infos) <- c("method", "message")
+ asVar <- NULL
+ asBias <- NULL
+ }
+
+ new("kStepEstimate", name = "1-step estimate", estimate = res,
+ samplesize = length(x), asvar = asVar, asbias = asBias, pIC = IC, steps = 1L,
+ Infos = Infos)
})
setMethod("oneStepEstimator", signature(x = "matrix",
IC = "InfluenceCurve",
start = "Estimate"),
- function(x, IC, start){
+ function(x, IC, start, useLast = getRobAStBaseOption("kStepUseLast")){
nrvalues <- dimension(IC at Curve)
start0 <- estimate(start)
if(is.list(start0)) start0 <- unlist(start0)
@@ -57,5 +201,53 @@
res <- start0 + rowMeans(evalIC(IC, x), na.rm = TRUE)
- return(res)
+ if(is(IC, "IC")){
+ L2Fam <- eval(CallL2Fam(IC))
+ Infos <- matrix(c("oneStepEstimator",
+ paste("1-step estimate for", name(L2Fam))),
+ ncol = 2)
+ colnames(Infos) <- c("method", "message")
+ if(is(L2Fam, "L2GroupParamFamily")) useLast <- TRUE
+ if(useLast && !is(modifyIC(IC), "NULL") ){
+ newParam <- param(L2Fam)
+ main(newParam) <- res
+ newL2Fam <- modifyModel(L2Fam, newParam)
+ IC <- modifyIC(IC)(newL2Fam, IC)
+ Infos <- rbind(Infos, c("oneStepEstimator",
+ "computation of IC, asvar and asbias via useLast = TRUE"))
+ }else{
+ if(useLast && is(modifyIC(IC), "NULL")){
+ warning("'useLast = TRUE' only possible if slot 'modifyIC' of 'IC'
+ is filled with some function!")
+ Infos <- rbind(Infos, c("oneStepEstimator",
+ "slot 'modifyIC' of 'IC' was not filled!"))
+ }
+ Infos <- rbind(Infos, c("oneStepEstimator",
+ "computation of IC, asvar and asbias via useLast = FALSE"))
+ }
+ if("asCov" %in% names(Risks(IC)))
+ asVar <- Risks(IC)$asCov
+ else
+ asVar <- getRiskIC(IC, risk = asCov())$asCov$value
+
+ if("asBias" %in% names(Risks(IC))){
+ asBias <- neighborRadius(IC)*Risks(IC)$asBias
+ }else{
+ if(is(IC, "HampIC")){
+ r <- neighborRadius(IC)
+ asBias <- r*getRiskIC(IC, risk = asBias(), neighbor = neighbor(IC))$asBias$value
+ }else{
+ asBias <- NULL
+ }
+ }
+ }else{
+ Infos <- matrix(c("oneStepEstimator", "1-step estimate"), ncol = 2)
+ colnames(Infos) <- c("method", "message")
+ asVar <- NULL
+ asBias <- NULL
+ }
+
+ new("kStepEstimate", name = "1-step estimate", estimate = res,
+ samplesize = ncol(x), asvar = asVar, asbias = asBias, pIC = IC, steps = 1L,
+ Infos = Infos)
})
Modified: branches/robast-0.6/pkg/RobAStBase/R/optIC.R
===================================================================
--- branches/robast-0.6/pkg/RobAStBase/R/optIC.R 2008-07-28 12:27:14 UTC (rev 131)
+++ branches/robast-0.6/pkg/RobAStBase/R/optIC.R 2008-07-31 15:07:01 UTC (rev 132)
@@ -6,6 +6,8 @@
Curve <- as((model at param@trafo %*% solve(model at FisherInfo)) %*% model at L2deriv, "EuclRandVariable")
asCov <- model at param@trafo %*% solve(model at FisherInfo) %*% t(model at param@trafo)
+ modifyIC <- function(L2Fam, IC){ optIC(L2Fam, asCov()) }
+
return(IC(
name = paste("Classical optimal influence curve for", model at name),
CallL2Fam = call("L2ParamFamily",
@@ -22,7 +24,8 @@
L2derivDistrSymm = model at L2derivDistrSymm,
FisherInfo = model at FisherInfo,
FisherInfo.fct = model at FisherInfo.fct),
- Curve = EuclRandVarList(Curve),
+ Curve = EuclRandVarList(Curve),
+ modifyIC = modifyIC,
Risks = list(asCov = asCov, trAsCov = sum(diag(asCov))),
Infos = matrix(c("optIC", "optimal IC in sense of Cramer-Rao bound"),
ncol = 2, dimnames = list(character(0), c("method", "message")))))
Modified: branches/robast-0.6/pkg/RobAStBase/man/ALEstimate-class.Rd
===================================================================
--- branches/robast-0.6/pkg/RobAStBase/man/ALEstimate-class.Rd 2008-07-28 12:27:14 UTC (rev 131)
+++ branches/robast-0.6/pkg/RobAStBase/man/ALEstimate-class.Rd 2008-07-31 15:07:01 UTC (rev 132)
@@ -3,6 +3,8 @@
\alias{ALEstimate-class}
\alias{pIC}
\alias{pIC,ALEstimate-method}
+\alias{asbias}
+\alias{asbias,ALEstimate-method}
\title{ALEstimate-class.}
\description{Class of asymptotically linear estimates.}
@@ -15,8 +17,16 @@
name of the estimator. }
\item{\code{estimate}:}{Object of class \code{"ANY"}:
estimate. }
+ \item{\code{samplesize}:}{Object of class \code{"numeric"}:
+ sample size. }
+ \item{\code{asvar}:}{Optional object of class \code{"matrix"}:
+ asymptotic variance. }
+ \item{\code{asbias}:}{Optional object of class \code{"numeric"}:
+ asymptotic bias. }
\item{\code{pIC}:}{Optional object of class \code{InfluenceCurve}:
influence curve. }
+ \item{\code{nuis.idx}:}{ object of class \code{"OptionalNumeric"}:
+ indices of \code{estimate} belonging to the nuisance part. }
\item{\code{Infos}:}{ object of class \code{"matrix"}
with two columns named \code{method} and \code{message}:
additional informations. }
@@ -35,6 +45,9 @@
\author{Matthias Kohl \email{Matthias.Kohl at stamats.de}}
%\note{}
\seealso{\code{\link[distrMod]{Estimate-class}}}
-%\examples{}
+\examples{
+## prototype
+new("ALEstimate")
+}
\concept{estimate}
\keyword{classes}
Modified: branches/robast-0.6/pkg/RobAStBase/man/BdStWeight-class.Rd
===================================================================
--- branches/robast-0.6/pkg/RobAStBase/man/BdStWeight-class.Rd 2008-07-28 12:27:14 UTC (rev 131)
+++ branches/robast-0.6/pkg/RobAStBase/man/BdStWeight-class.Rd 2008-07-31 15:07:01 UTC (rev 132)
@@ -49,7 +49,8 @@
\seealso{\code{\link{BoundedWeight-class}}, \code{\link{RobWeight-class}},
\code{\link{IC}}, \code{\link{InfluenceCurve-class}}}
\examples{
-%
+## prototype
+new("BdStWeight")
}
\concept{influence curve}
\keyword{classes}
Modified: branches/robast-0.6/pkg/RobAStBase/man/BoundedWeight-class.Rd
===================================================================
--- branches/robast-0.6/pkg/RobAStBase/man/BoundedWeight-class.Rd 2008-07-28 12:27:14 UTC (rev 131)
+++ branches/robast-0.6/pkg/RobAStBase/man/BoundedWeight-class.Rd 2008-07-31 15:07:01 UTC (rev 132)
@@ -44,7 +44,8 @@
%\note{}
\seealso{\code{\link{RobWeight-class}}, \code{\link{IC}}, \code{\link{InfluenceCurve-class}}}
\examples{
-%
+## prototype
+new("BoundedWeight")
}
\concept{influence curve}
\keyword{classes}
Modified: branches/robast-0.6/pkg/RobAStBase/man/ContIC-class.Rd
===================================================================
--- branches/robast-0.6/pkg/RobAStBase/man/ContIC-class.Rd 2008-07-28 12:27:14 UTC (rev 131)
+++ branches/robast-0.6/pkg/RobAStBase/man/ContIC-class.Rd 2008-07-31 15:07:01 UTC (rev 132)
@@ -14,6 +14,7 @@
\alias{lowerCase<-,ContIC-method}
\alias{stand<-}
\alias{stand<-,ContIC-method}
+\alias{neighbor,ContIC-method}
\alias{generateIC,ContNeighborhood,L2ParamFamily-method}
\alias{show,ContIC-method}
@@ -41,6 +42,12 @@
\item{\code{Curve}:}{ object of class \code{"EuclRandVarList"}}
+ \item{\code{modifyIC}:}{ Object of class \code{"OptionalFunction"}:
+ function of two arguments, which are an L2 parametric family
+ and an optional influence curve. Returns an object of
+ class \code{"IC"}. This slot is mainly used for internal
+ computations! }
+
\item{\code{Risks}:}{ object of class \code{"list"}:
list of risks; cf. \code{\link[distrMod]{RiskType-class}}. }
@@ -59,10 +66,10 @@
\item{\code{weight}:}{ object of class \code{"HampelWeight"}:
weight function }
-
+
\item{\code{biastype}:}{ object of class \code{"BiasType"}:
bias type (symmetric/onsided/asymmetric) }
-
+
\item{\code{normtype}:}{ object of class \code{"NormType"}:
norm type (Euclidean, information/self-standardized)}
@@ -102,6 +109,10 @@
\item{lowerCase<-}{\code{signature(object = "ContIC")}:
replacement function for slot \code{lowerCase}. }
+ \item{neighbor}{\code{signature(object = "ContIC")}:
+ generates an object of class \code{"ContNeighborhood"} with
+ radius given in slot \code{neighborRadius}. }
+
\item{generateIC}{\code{signature(neighbor = "ContNeighborhood", L2Fam = "L2ParamFamily")}:
generate an object of class \code{"ContIC"}. Rarely called directly. }
Modified: branches/robast-0.6/pkg/RobAStBase/man/ContIC.Rd
===================================================================
--- branches/robast-0.6/pkg/RobAStBase/man/ContIC.Rd 2008-07-28 12:27:14 UTC (rev 131)
+++ branches/robast-0.6/pkg/RobAStBase/man/ContIC.Rd 2008-07-31 15:07:01 UTC (rev 132)
@@ -17,7 +17,8 @@
Domain = Reals())),
Risks, Infos, clip = Inf, cent = 0, stand = as.matrix(1),
lowerCase = NULL, neighborRadius = 0, w = new("HampelWeight"),
- normtype = NormType(), biastype = symmetricBias())
+ normtype = NormType(), biastype = symmetricBias(),
+ modifyIC = NULL)
}
\arguments{
\item{name}{ object of class \code{"character"}. }
@@ -38,6 +39,11 @@
contamination neighborhood. }
\item{biastype}{ BiasType: type of the bias}
\item{normtype}{ NormType: type of the norm}
+ \item{modifyIC}{ object of class \code{"OptionalFunction"}:
+ function of two arguments, which are an L2 parametric family
+ and an optional influence curve. Returns an object of
+ class \code{"IC"}. This function is mainly used for internal
+ computations! }
}
%\details{}
\value{Object of class \code{"ContIC"}}
@@ -55,4 +61,4 @@
plot(IC1)
}
\concept{influence curve}
-\keyword{}
+\keyword{robust}
Modified: branches/robast-0.6/pkg/RobAStBase/man/HampIC-class.Rd
===================================================================
--- branches/robast-0.6/pkg/RobAStBase/man/HampIC-class.Rd 2008-07-28 12:27:14 UTC (rev 131)
+++ branches/robast-0.6/pkg/RobAStBase/man/HampIC-class.Rd 2008-07-31 15:07:01 UTC (rev 132)
@@ -30,6 +30,12 @@
\item{\code{Curve}:}{ object of class \code{"EuclRandVarList"}}
+ \item{\code{modifyIC}:}{ Object of class \code{"OptionalFunction"}:
+ function of two arguments, which are an L2 parametric family
+ and an optional influence curve. Returns an object of
+ class \code{"IC"}. This slot is mainly used for internal
+ computations! }
+
\item{\code{Risks}:}{ object of class \code{"list"}:
list of risks; cf. \code{\link[distrMod]{RiskType-class}}. }
Modified: branches/robast-0.6/pkg/RobAStBase/man/HampelWeight-class.Rd
===================================================================
--- branches/robast-0.6/pkg/RobAStBase/man/HampelWeight-class.Rd 2008-07-28 12:27:14 UTC (rev 131)
+++ branches/robast-0.6/pkg/RobAStBase/man/HampelWeight-class.Rd 2008-07-31 15:07:01 UTC (rev 132)
@@ -52,7 +52,8 @@
\code{\link{BoundedWeight-class}}, \code{\link{RobWeight-class}},
\code{\link{IC}}, \code{\link{InfluenceCurve-class}}}
\examples{
-%
+## prototype
+new("HampelWeight")
}
\concept{influence curve}
\keyword{classes}
Modified: branches/robast-0.6/pkg/RobAStBase/man/IC-class.Rd
===================================================================
--- branches/robast-0.6/pkg/RobAStBase/man/IC-class.Rd 2008-07-28 12:27:14 UTC (rev 131)
+++ branches/robast-0.6/pkg/RobAStBase/man/IC-class.Rd 2008-07-31 15:07:01 UTC (rev 132)
@@ -5,6 +5,8 @@
\alias{CallL2Fam,IC-method}
\alias{CallL2Fam<-}
\alias{CallL2Fam<-,IC-method}
+\alias{modifyIC}
+\alias{modifyIC,IC-method}
\alias{checkIC,IC,missing-method}
\alias{checkIC,IC,L2ParamFamily-method}
\alias{evalIC,IC,numeric-method}
@@ -25,6 +27,11 @@
\item{\code{CallL2Fam}:}{Object of class \code{"call"}:
creates an object of the underlying L2-differentiable
parametric family. }
+ \item{\code{modifyIC}:}{ Object of class \code{"OptionalFunction"}:
+ function of two arguments, which are an L2 parametric family
+ and an optional influence curve. Returns an object of
+ class \code{"IC"}. This slot is mainly used for internal
+ computations! }
\item{\code{name}:}{Object of class \code{"character"}. }
\item{\code{Curve}:}{Object of class \code{"EuclRandVarList"}.}
\item{\code{Risks}:}{Object of class \code{"list"}:
@@ -45,6 +52,9 @@
\item{CallL2Fam<-}{\code{signature(object = "IC")}:
replacement function for slot \code{CallL2Fam}. }
+ \item{modifyIC}{\code{signature(object = "IC")}:
+ accessor function for slot \code{modifyIC}. }
+
\item{checkIC}{\code{signature(IC = "IC", L2Fam = "missing")}:
check centering and Fisher consistency of \code{IC} assuming
the L2-differentiable parametric family which can
@@ -86,3 +96,4 @@
}
\concept{influence curve}
\keyword{classes}
+\keyword{robust}
Modified: branches/robast-0.6/pkg/RobAStBase/man/IC.Rd
===================================================================
--- branches/robast-0.6/pkg/RobAStBase/man/IC.Rd 2008-07-28 12:27:14 UTC (rev 131)
+++ branches/robast-0.6/pkg/RobAStBase/man/IC.Rd 2008-07-31 15:07:01 UTC (rev 132)
@@ -8,7 +8,7 @@
\usage{
IC(name, Curve = EuclRandVarList(RealRandVariable(Map = list(function(x){x}),
Domain = Reals())),
- Risks, Infos, CallL2Fam = call("L2ParamFamily"))
+ Risks, Infos, CallL2Fam = call("L2ParamFamily"), modifyIC = NULL)
}
\arguments{
\item{name}{ Object of class \code{"character"}. }
@@ -20,6 +20,11 @@
list of risks; cf. \code{\link[distrMod]{RiskType-class}}. }
\item{Infos}{ matrix of characters with two columns
named \code{method} and \code{message}: additional informations. }
+ \item{modifyIC}{ Object of class \code{"OptionalFunction"}:
+ function of two arguments, which are an L2 parametric family
+ and an optional influence curve. Returns an object of
+ class \code{"IC"}. This function is mainly used for internal
+ computations! }
}
%\details{}
\value{Object of class \code{"IC"}}
@@ -38,20 +43,6 @@
\examples{
IC1 <- IC()
plot(IC1)
-
-## The function is currently defined as
-IC <- function(name, Curve = EuclRandVarList(RealRandVariable(Map = list(function(x){x})),
- Domain = Reals()), Risks, Infos, CallL2Fam = call("L2ParamFamily")){
- if(missing(name))
- name <- "square integrable (partial) influence curve"
- if(missing(Risks))
- Risks <- list()
- if(missing(Infos))
- Infos <- matrix(c(character(0),character(0)), ncol=2,
- dimnames=list(character(0), c("method", "message")))
- return(new("IC", name = name, Curve = Curve, Risks = Risks,
- Infos = Infos, CallL2Fam = CallL2Fam))
}
-}
\concept{influence curve}
-\keyword{}
+\keyword{robust}
Modified: branches/robast-0.6/pkg/RobAStBase/man/InfluenceCurve-class.Rd
===================================================================
--- branches/robast-0.6/pkg/RobAStBase/man/InfluenceCurve-class.Rd 2008-07-28 12:27:14 UTC (rev 131)
+++ branches/robast-0.6/pkg/RobAStBase/man/InfluenceCurve-class.Rd 2008-07-31 15:07:01 UTC (rev 132)
@@ -98,3 +98,4 @@
}
\concept{influence curve}
\keyword{classes}
+\keyword{robust}
Modified: branches/robast-0.6/pkg/RobAStBase/man/InfluenceCurve.Rd
===================================================================
--- branches/robast-0.6/pkg/RobAStBase/man/InfluenceCurve.Rd 2008-07-28 12:27:14 UTC (rev 131)
+++ branches/robast-0.6/pkg/RobAStBase/man/InfluenceCurve.Rd 2008-07-31 15:07:01 UTC (rev 132)
@@ -49,4 +49,4 @@
}
}
\concept{influence curve}
-\keyword{}
+\keyword{robust}
Modified: branches/robast-0.6/pkg/RobAStBase/man/MEstimate-class.Rd
===================================================================
--- branches/robast-0.6/pkg/RobAStBase/man/MEstimate-class.Rd 2008-07-28 12:27:14 UTC (rev 131)
+++ branches/robast-0.6/pkg/RobAStBase/man/MEstimate-class.Rd 2008-07-31 15:07:01 UTC (rev 132)
@@ -17,8 +17,16 @@
name of the estimator. }
\item{\code{estimate}:}{Object of class \code{"ANY"}:
estimate. }
+ \item{\code{samplesize}:}{Object of class \code{"numeric"}:
+ sample size. }
+ \item{\code{asvar}:}{Optional object of class \code{"matrix"}:
+ asymptotic variance. }
+ \item{\code{asbias}:}{Optional object of class \code{"numeric"}:
+ asymptotic bias. }
\item{\code{pIC}:}{Optional object of class \code{InfluenceCurve}:
influence curve. }
+ \item{\code{nuis.idx}:}{ object of class \code{"OptionalNumeric"}:
+ indices of \code{estimate} belonging to the nuisance part. }
\item{\code{Mroot}:}{Object of class \code{"numeric"}: value of
the M equation at the estimate. }
\item{\code{Infos}:}{ object of class \code{"matrix"}
Added: branches/robast-0.6/pkg/RobAStBase/man/RobAStBaseOptions.Rd
===================================================================
--- branches/robast-0.6/pkg/RobAStBase/man/RobAStBaseOptions.Rd (rev 0)
+++ branches/robast-0.6/pkg/RobAStBase/man/RobAStBaseOptions.Rd 2008-07-31 15:07:01 UTC (rev 132)
@@ -0,0 +1,49 @@
+\name{RobAStBaseOptions}
+\alias{RobAStBaseOptions}
+\alias{getRobAStBaseOption}
+\alias{kStepUseLast}
+
+\title{Function to change the global variables of the package `RobAStBase' }
+\description{With \code{RobAStBaseOptions} you can inspect and change
+ the global variables of the package \pkg{RobAStBase}. }
+\usage{
+RobAStBaseOptions(...)
+getRobAStBaseOption(x)
+}
+\arguments{
+ \item{\dots}{ any options can be defined, using name = value or by passing a list of such tagged values. }
+ \item{x}{ a character string holding an option name.}
+}
+%\details{}
+\value{
+ \code{RobAStBaseOptions()} returns a list of the global variables.\newline
+ \code{RobAStBaseOptions(x)} returns the global variable \var{x}.\newline
+ \code{getRobAStBaseOption(x)} returns the global variable \var{x}.\newline
+ \code{RobAStBaseOptions(x=y)} sets the value of the global variable \var{x} to \var{y}.
+}
+\section{Global Options}{
+\describe{
+ \item{kStepUseLast:}{ The default value of argument \code{kStepUseLast} is
+ \code{FALSE}. Explicitly setting \code{kStepUseLast} to \code{TRUE} should
+ be done with care as in this situation the influence curve in case of
+ \code{oneStepEstimator} and \code{kStepEstimator} is re-computed using
+ the value of the one- resp. k-step estimate which may take quite a long
+ time depending on the model. }
+}
+}
+%\references{}
+\author{Matthias Kohl \email{Matthias.Kohl at stamats.de}}
+%\note{}
+\seealso{\code{\link[base]{options}}, \code{\link[base]{getOption}}}
+\examples{
+RobAStBaseOptions()
+RobAStBaseOptions("kStepUseLast")
+RobAStBaseOptions("kStepUseLast" = TRUE)
+# or
+RobAStBaseOptions(kStepUseLast = 1e-6)
+getRobAStBaseOption("kStepUseLast")
+}
+\keyword{misc}
+\keyword{robust}
+\concept{global options}
+\concept{options}
Modified: branches/robast-0.6/pkg/RobAStBase/man/RobAStControl-class.Rd
===================================================================
--- branches/robast-0.6/pkg/RobAStBase/man/RobAStControl-class.Rd 2008-07-28 12:27:14 UTC (rev 131)
+++ branches/robast-0.6/pkg/RobAStBase/man/RobAStControl-class.Rd 2008-07-31 15:07:01 UTC (rev 132)
@@ -7,17 +7,14 @@
\title{Control classes in package RobAStBase}
\description{Control classes in package \pkg{RobAStBase}.}
\section{Objects from the Class}{
+ This class is virtual; that is no objects may be created.
}
\section{Slots}{
-% \describe{
-% \item{\code{CallL2Fam}:}{Object of class \code{"call"}:
-% creates an object of the underlying L2-differentiable
-% parametric family. }
-% }
+ \describe{
+ \item{\code{name}:}{Object of class \code{"character"}:
+ name of the control object. }
+ }
}
-%\section{Extends}{
-%Class \code{"InfluenceCurve"}, directly.
-%}
\section{Methods}{
\describe{
\item{name}{\code{signature(object = "RobAStControl")}:
@@ -39,11 +36,7 @@
}
\author{Peter Ruckdeschel \email{Peter.Ruckdeschel at itwm.fraunhofer.de}}
%\note{}
-\seealso{
-%
-}
-\examples{
-%
-}
+%\seealso{}
+%\examples{}
\concept{influence curve}
\keyword{classes}
Modified: branches/robast-0.6/pkg/RobAStBase/man/RobWeight-class.Rd
===================================================================
--- branches/robast-0.6/pkg/RobAStBase/man/RobWeight-class.Rd 2008-07-28 12:27:14 UTC (rev 131)
+++ branches/robast-0.6/pkg/RobAStBase/man/RobWeight-class.Rd 2008-07-31 15:07:01 UTC (rev 132)
@@ -48,7 +48,8 @@
%\note{}
\seealso{\code{\link{InfluenceCurve-class}}, \code{\link{IC}}}
\examples{
-%
+## prototype
+new("RobWeight")
}
\concept{influence curve}
\keyword{classes}
Modified: branches/robast-0.6/pkg/RobAStBase/man/TotalVarIC-class.Rd
===================================================================
--- branches/robast-0.6/pkg/RobAStBase/man/TotalVarIC-class.Rd 2008-07-28 12:27:14 UTC (rev 131)
+++ branches/robast-0.6/pkg/RobAStBase/man/TotalVarIC-class.Rd 2008-07-31 15:07:01 UTC (rev 132)
@@ -11,6 +11,7 @@
\alias{clipUp<-}
\alias{clipUp<-,TotalVarIC-method}
\alias{lowerCase<-,TotalVarIC-method}
+\alias{neighbor,TotalVarIC-method}
\alias{show,TotalVarIC-method}
\alias{stand<-,TotalVarIC-method}
\alias{generateIC,TotalVarNeighborhood,L2ParamFamily-method}
@@ -39,6 +40,12 @@
\item{\code{Curve}:}{ object of class \code{"EuclRandVarList"}.}
+ \item{\code{modifyIC}:}{ Object of class \code{"OptionalFunction"}:
+ function of two arguments, which are an L2 parametric family
+ and an optional influence curve. Returns an object of
+ class \code{"IC"}. This slot is mainly used for internal
+ computations! }
+
\item{\code{Risks}:}{ object of class \code{"list"}:
list of risks; cf. \code{\link[distrMod]{RiskType-class}}. }
@@ -58,6 +65,12 @@
\item{\code{weight}:}{ object of class \code{"BdStWeight"}:
weight function }
+ \item{\code{biastype}:}{ object of class \code{"BiasType"}:
+ bias type (symmetric/onsided/asymmetric) }
+
+ \item{\code{normtype}:}{ object of class \code{"NormType"}:
+ norm type (Euclidean, information/self-standardized)}
+
\item{\code{neighborRadius}:}{ object of class \code{"numeric"}:
radius of the corresponding (unconditional) contamination
neighborhood. }
@@ -88,6 +101,13 @@
\item{stand<-}{\code{signature(object = "TotalVarIC")}:
replacement function for slot \code{stand}. }
+ \item{lowerCase<-}{\code{signature(object = "TotalVarIC")}:
+ replacement function for slot \code{lowerCase}. }
+
+ \item{neighbor}{\code{signature(object = "TotalVarIC")}:
+ generates an object of class \code{"TotalVarNeighborhood"} with
+ radius given in slot \code{neighborRadius}. }
+
\item{generateIC}{\code{signature(neighbor = "TotalVarNeighborhood", L2Fam = "L2ParamFamily")}:
generate an object of class \code{"TotalVarIC"}. Rarely called directly. }
@@ -109,3 +129,4 @@
}
\concept{influence curve}
\keyword{classes}
+\keyword{robust}
Modified: branches/robast-0.6/pkg/RobAStBase/man/TotalVarIC.Rd
===================================================================
--- branches/robast-0.6/pkg/RobAStBase/man/TotalVarIC.Rd 2008-07-28 12:27:14 UTC (rev 131)
+++ branches/robast-0.6/pkg/RobAStBase/man/TotalVarIC.Rd 2008-07-31 15:07:01 UTC (rev 132)
@@ -16,7 +16,9 @@
Curve = EuclRandVarList(RealRandVariable(Map = c(function(x) {x}),
Domain = Reals())),
Risks, Infos, clipLo = -Inf, clipUp = Inf, stand = as.matrix(1),
- lowerCase = NULL, neighborRadius = 0, w = new("BdStWeight"))
+ lowerCase = NULL, neighborRadius = 0, w = new("BdStWeight"),
+ normtype = NormType(), biastype = symmetricBias(),
+ modifyIC = NULL)
}
\arguments{
\item{name}{ object of class \code{"character"}. }
@@ -35,6 +37,13 @@
\item{lowerCase}{ optional constant for lower case solution. }
\item{neighborRadius}{ radius of the corresponding (unconditional)
contamination neighborhood. }
+ \item{biastype}{ BiasType: type of the bias}
+ \item{normtype}{ NormType: type of the norm}
+ \item{modifyIC}{ object of class \code{"OptionalFunction"}:
+ function of two arguments, which are an L2 parametric family
+ and an optional influence curve. Returns an object of
+ class \code{"IC"}. This function is mainly used for internal
+ computations! }
}
%\details{}
\value{Object of class \code{"TotalVarIC"}}
@@ -52,4 +61,4 @@
plot(IC1)
}
\concept{influence curve}
-\keyword{}
+\keyword{robust}
Modified: branches/robast-0.6/pkg/RobAStBase/man/checkIC.Rd
===================================================================
--- branches/robast-0.6/pkg/RobAStBase/man/checkIC.Rd 2008-07-28 12:27:14 UTC (rev 131)
+++ branches/robast-0.6/pkg/RobAStBase/man/checkIC.Rd 2008-07-31 15:07:01 UTC (rev 132)
@@ -33,4 +33,4 @@
checkIC(IC1)
}
\concept{influence curve}
-\keyword{}
+\keyword{robust}
Modified: branches/robast-0.6/pkg/RobAStBase/man/comparePlot.Rd
===================================================================
--- branches/robast-0.6/pkg/RobAStBase/man/comparePlot.Rd 2008-07-28 12:27:14 UTC (rev 131)
+++ branches/robast-0.6/pkg/RobAStBase/man/comparePlot.Rd 2008-07-31 15:07:01 UTC (rev 132)
@@ -39,4 +39,4 @@
comparePlot(IC1,IC2)
}
}
-\keyword{}
+\keyword{robust}
Modified: branches/robast-0.6/pkg/RobAStBase/man/evalIC.Rd
===================================================================
--- branches/robast-0.6/pkg/RobAStBase/man/evalIC.Rd 2008-07-28 12:27:14 UTC (rev 131)
+++ branches/robast-0.6/pkg/RobAStBase/man/evalIC.Rd 2008-07-31 15:07:01 UTC (rev 132)
@@ -27,4 +27,4 @@
\seealso{\code{\link{IC-class}}}
%\examples{}
\concept{influence curve}
-\keyword{}
+\keyword{robust}
Modified: branches/robast-0.6/pkg/RobAStBase/man/generateIC.Rd
===================================================================
--- branches/robast-0.6/pkg/RobAStBase/man/generateIC.Rd 2008-07-28 12:27:14 UTC (rev 131)
+++ branches/robast-0.6/pkg/RobAStBase/man/generateIC.Rd 2008-07-31 15:07:01 UTC (rev 132)
@@ -27,4 +27,4 @@
\seealso{\code{\link{IC-class}}, \code{\link{ContIC-class}}, \code{\link{TotalVarIC-class}}}
%\examples{}
\concept{influence curve}
-\keyword{}
+\keyword{robust}
Modified: branches/robast-0.6/pkg/RobAStBase/man/generateICfct.Rd
===================================================================
--- branches/robast-0.6/pkg/RobAStBase/man/generateICfct.Rd 2008-07-28 12:27:14 UTC (rev 131)
+++ branches/robast-0.6/pkg/RobAStBase/man/generateICfct.Rd 2008-07-31 15:07:01 UTC (rev 132)
@@ -30,8 +30,6 @@
\author{Peter Ruckdeschel \email{Peter.Ruckdeschel at itwm.fraunhofer.de}}
%\note{}
\seealso{\code{\link[distrMod]{L2ParamFamily-class}}, \code{\link{IC-class}}}
-\examples{
-%
-}
+%\examples{}
\concept{influence curve}
-\keyword{}
+\keyword{robust}
Modified: branches/robast-0.6/pkg/RobAStBase/man/getBiasIC.Rd
===================================================================
--- branches/robast-0.6/pkg/RobAStBase/man/getBiasIC.Rd 2008-07-28 12:27:14 UTC (rev 131)
+++ branches/robast-0.6/pkg/RobAStBase/man/getBiasIC.Rd 2008-07-31 15:07:01 UTC (rev 132)
@@ -23,7 +23,7 @@
\item{tol}{ the desired accuracy (convergence tolerance).}
\item{numbeval}{number of evalation points.}
}
-\details{}
+%\details{}
\value{The bias of the IC is computed.}
\section{Methods}{
\describe{
@@ -63,4 +63,4 @@
\seealso{\code{\link{getRiskIC-methods}}, \code{\link{InfRobModel-class}}}
%\examples{}
\concept{influence curve}
-\keyword{}
+\keyword{robust}
Modified: branches/robast-0.6/pkg/RobAStBase/man/getRiskIC.Rd
===================================================================
--- branches/robast-0.6/pkg/RobAStBase/man/getRiskIC.Rd 2008-07-28 12:27:14 UTC (rev 131)
+++ branches/robast-0.6/pkg/RobAStBase/man/getRiskIC.Rd 2008-07-31 15:07:01 UTC (rev 132)
@@ -119,4 +119,4 @@
\seealso{\code{\link[ROptEst]{getRiskIC-methods}}, \code{\link{InfRobModel-class}}}
%\examples{}
\concept{influence curve}
-\keyword{}
+\keyword{robust}
Modified: branches/robast-0.6/pkg/RobAStBase/man/getweight.Rd
===================================================================
--- branches/robast-0.6/pkg/RobAStBase/man/getweight.Rd 2008-07-28 12:27:14 UTC (rev 131)
+++ branches/robast-0.6/pkg/RobAStBase/man/getweight.Rd 2008-07-31 15:07:01 UTC (rev 132)
@@ -75,9 +75,6 @@
\seealso{\code{\link{BdStWeight-class}},
\code{\link{HampelWeight-class}},
\code{\link{IC-class}}}
-\examples{
-%
-}
-
+%\examples{}
\concept{influence curve}
-\keyword{}
+\keyword{robust}
Modified: branches/robast-0.6/pkg/RobAStBase/man/infoPlot.Rd
===================================================================
--- branches/robast-0.6/pkg/RobAStBase/man/infoPlot.Rd 2008-07-28 12:27:14 UTC (rev 131)
+++ branches/robast-0.6/pkg/RobAStBase/man/infoPlot.Rd 2008-07-31 15:07:01 UTC (rev 132)
@@ -33,4 +33,4 @@
}
\concept{absolute information}
\concept{relative information}
-\keyword{}
+\keyword{robust}
Modified: branches/robast-0.6/pkg/RobAStBase/man/kStepEstimate-class.Rd
===================================================================
--- branches/robast-0.6/pkg/RobAStBase/man/kStepEstimate-class.Rd 2008-07-28 12:27:14 UTC (rev 131)
+++ branches/robast-0.6/pkg/RobAStBase/man/kStepEstimate-class.Rd 2008-07-31 15:07:01 UTC (rev 132)
@@ -17,8 +17,16 @@
name of the estimator. }
\item{\code{estimate}:}{Object of class \code{"ANY"}:
estimate. }
+ \item{\code{samplesize}:}{Object of class \code{"numeric"}:
+ sample size. }
+ \item{\code{asvar}:}{Optional object of class \code{"matrix"}:
+ asymptotic variance. }
+ \item{\code{asbias}:}{Optional object of class \code{"numeric"}:
+ asymptotic bias. }
\item{\code{pIC}:}{Optional object of class \code{InfluenceCurve}:
influence curve. }
+ \item{\code{nuis.idx}:}{ object of class \code{"OptionalNumeric"}:
+ indices of \code{estimate} belonging to the nuisance part. }
\item{\code{steps}:}{Object of class \code{"integer"}: number
of steps. }
\item{\code{Infos}:}{ object of class \code{"matrix"}
Added: branches/robast-0.6/pkg/RobAStBase/man/kStepEstimator.Rd
===================================================================
--- branches/robast-0.6/pkg/RobAStBase/man/kStepEstimator.Rd (rev 0)
+++ branches/robast-0.6/pkg/RobAStBase/man/kStepEstimator.Rd 2008-07-31 15:07:01 UTC (rev 132)
@@ -0,0 +1,83 @@
+\name{kStepEstimator}
+\alias{kStepEstimator}
+\alias{kStepEstimator-methods}
+\alias{kStepEstimator,numeric,IC,numeric-method}
+\alias{kStepEstimator,matrix,IC,numeric-method}
+\alias{kStepEstimator,numeric,IC,Estimate-method}
+\alias{kStepEstimator,matrix,IC,Estimate-method}
+
+\title{Generic function for the computation of k-step estimates}
+\description{
+ Generic function for the computation of k-step estimates.
+}
+\usage{
+kStepEstimator(x, IC, start, ...)
+
+\S4method{kStepEstimator}{numeric,IC,numeric}(x, IC, start, steps = 1L,
+ useLast = getRobAStBaseOption("kStepUseLast"))
+\S4method{kStepEstimator}{matrix,IC,numeric}(x, IC, start, steps = 1L,
+ useLast = getRobAStBaseOption("kStepUseLast"))
+\S4method{kStepEstimator}{numeric,IC,Estimate}(x, IC, start, steps = 1L,
+ useLast = getRobAStBaseOption("kStepUseLast"))
+\S4method{kStepEstimator}{matrix,IC,Estimate}(x, IC, start, steps = 1L,
+ useLast = getRobAStBaseOption("kStepUseLast"))
+}
+\arguments{
+ \item{x}{ sample }
+ \item{IC}{ object of class \code{"IC"} }
+ \item{start}{ initial estimate }
+ \item{steps}{ integer: number of steps }
+ \item{useLast}{ which parameter estimate (initial estimate or
+ k-step estimate) shall be used to fill the slots \code{pIC},
+ \code{asvar} and \code{asbias} of the return value. }
+ \item{...}{ additional parameters }
+}
+\details{
+ Given an initial estimation \code{start}, a sample \code{x}
+ and an influence curve \code{IC} the corresponding k-step
+ estimator is computed.
+
+ The default value of argument \code{useLast} is set by the
+ global option \code{kStepUseLast} which by default is set to
+ \code{FALSE}. In case of general models \code{useLast}
+ remains unchanged during the computations. However, if
+ slot \code{CallL2Fam} of \code{IC} generates an object of
+ class \code{"L2GroupParamFamily"} the value of \code{useLast}
+ is changed to \code{TRUE}.
+ Explicitly setting \code{useLast} to \code{TRUE} should
+ be done with care as in this situation the influence curve
+ is re-computed using the value of the one-step estimate
+ which may take quite a long time depending on the model.
+
+ If \code{useLast} is set to \code{TRUE} and slot \code{modifyIC}
+ of \code{IC} is filled with some function (which can be
+ used to re-compute the IC for a different parameter), the
+ computation of \code{asvar}, \code{asbias} and \code{IC} is
+ based on the k-step estimate.
+}
+\value{Object of class \code{"kStepEstimate"}.}
+\section{Methods}{
+\describe{
+ \item{x = "numeric", IC = "IC", start = "numeric"}{
+ univariate samples. }
+ \item{x = "matrix", IC = "IC", start = "numeric"}{
+ multivariate samples. }
+ \item{x = "matrix", IC = "IC", start = "Estimate"}{
+ multivariate samples. }
+ \item{x = "matrix", IC = "IC", start = "Estimate"}{
+ multivariate samples. }
+}}
+\references{
+ Rieder, H. (1994) \emph{Robust Asymptotic Statistics}. New York: Springer.
+
+ Kohl, M. (2005) \emph{Numerical Contributions to the Asymptotic Theory of Robustness}.
+ Bayreuth: Dissertation.
+}
+\author{Matthias Kohl \email{Matthias.Kohl at stamats.de}}
+%\note{}
+\seealso{\code{\link{IC-class}}, \code{\link{kStepEstimate-class}} }
+%\examples{}
+\concept{k-step estimator}
+\concept{estimator}
+\keyword{univar}
+\keyword{robust}
Modified: branches/robast-0.6/pkg/RobAStBase/man/locMEstimator.Rd
===================================================================
--- branches/robast-0.6/pkg/RobAStBase/man/locMEstimator.Rd 2008-07-28 12:27:14 UTC (rev 131)
+++ branches/robast-0.6/pkg/RobAStBase/man/locMEstimator.Rd 2008-07-31 15:07:01 UTC (rev 132)
@@ -18,9 +18,10 @@
\item{\dots}{ additional parameters }
\item{eps}{ the desired accuracy (convergence tolerance). }
}
-%\details{}
+\details{ Given some sample \code{x} and some influence curve \code{IC}
+ an M estimate is computed by solving the corresponding
+ M equation. }
\value{Object of class \code{"MEstimate"}}
-}
\section{Methods}{
\describe{
\item{x = "numeric", IC = "InfluenceCurve"}{ univariate location. }
Modified: branches/robast-0.6/pkg/RobAStBase/man/makeIC-methods.Rd
===================================================================
--- branches/robast-0.6/pkg/RobAStBase/man/makeIC-methods.Rd 2008-07-28 12:27:14 UTC (rev 131)
+++ branches/robast-0.6/pkg/RobAStBase/man/makeIC-methods.Rd 2008-07-31 15:07:01 UTC (rev 132)
@@ -33,9 +33,26 @@
%\note{}
\seealso{\code{\link[distrMod]{L2ParamFamily-class}}, \code{\link{IC-class}}}
\examples{
+## default IC
IC1 <- new("IC")
+
+## L2-differentiable parametric family
B <- BinomFamily(13, 0.3)
-makeIC(IC1,B)
+
+## check IC properties
+checkIC(IC1, B)
+
+## make IC
+IC2 <- makeIC(IC1, B)
+
+## check IC properties
+checkIC(IC2)
+
+## slot modifyIC is filled in case of IC2
+IC3 <- modifyIC(IC2)(BinomFamily(13, 0.2), IC2)
+checkIC(IC3)
+## identical to
+checkIC(IC3, BinomFamily(13, 0.2))
}
\concept{influence curve}
-\keyword{}
+\keyword{robust}
Modified: branches/robast-0.6/pkg/RobAStBase/man/oneStepEstimator.Rd
===================================================================
--- branches/robast-0.6/pkg/RobAStBase/man/oneStepEstimator.Rd 2008-07-28 12:27:14 UTC (rev 131)
+++ branches/robast-0.6/pkg/RobAStBase/man/oneStepEstimator.Rd 2008-07-31 15:07:01 UTC (rev 132)
@@ -11,19 +11,54 @@
Generic function for the computation of one-step estimates.
}
\usage{
-oneStepEstimator(x, IC, start)
+oneStepEstimator(x, IC, start, ...)
+
+\S4method{oneStepEstimator}{numeric,InfluenceCurve,numeric}(x, IC, start,
+ useLast = getRobAStBaseOption("kStepUseLast"))
+\S4method{oneStepEstimator}{matrix,InfluenceCurve,numeric}(x, IC, start,
+ useLast = getRobAStBaseOption("kStepUseLast"))
+\S4method{oneStepEstimator}{numeric,InfluenceCurve,Estimate}(x, IC, start,
+ useLast = getRobAStBaseOption("kStepUseLast"))
+\S4method{oneStepEstimator}{matrix,InfluenceCurve,Estimate}(x, IC, start,
+ useLast = getRobAStBaseOption("kStepUseLast"))
}
\arguments{
\item{x}{ sample }
\item{IC}{ object of class \code{"InfluenceCurve"} }
\item{start}{ initial estimate }
+ \item{useLast}{ which parameter estimate (initial estimate or
+ one-step estimate) shall be used to fill the slots \code{pIC},
+ \code{asvar} and \code{asbias} of the return value. }
+ \item{...}{ additional arguments }
}
\details{
Given an initial estimation \code{start}, a sample \code{x}
and an influence curve \code{IC} the corresponding one-step
- estimator is computed
+ estimator is computed.
+
+ In case \code{IC} is an object of class \code{"IC"}
+ the slots \code{asvar} and \code{asbias} of the return
+ value are filled (based on the initial estimate).
+
+ The default value of argument \code{useLast} is set by the
+ global option \code{kStepUseLast} which by default is set to
+ \code{FALSE}. In case of general models \code{useLast}
+ remains unchanged during the computations. However, if
+ slot \code{CallL2Fam} of \code{IC} generates an object of
+ class \code{"L2GroupParamFamily"} the value of \code{useLast}
+ is changed to \code{TRUE}.
+ Explicitly setting \code{useLast} to \code{TRUE} should
+ be done with care as in this situation the influence curve
+ is re-computed using the value of the one-step estimate
+ which may take quite a long time depending on the model.
+
+ If \code{useLast} is set to \code{TRUE} and slot \code{modifyIC}
+ of \code{IC} is filled with some function (which can be
+ used to re-compute the IC for a different parameter), the
+ computation of \code{asvar}, \code{asbias} and \code{IC} is
+ based on the one-step estimate.
}
-\value{The one-step estimation is computed.}
+\value{Object of class \code{"kStepEstimate"}}
\section{Methods}{
\describe{
\item{x = "numeric", IC = "InfluenceCurve", start = "numeric"}{
@@ -43,7 +78,7 @@
}
\author{Matthias Kohl \email{Matthias.Kohl at stamats.de}}
%\note{}
-\seealso{\code{\link{InfluenceCurve-class}}}
+\seealso{\code{\link{InfluenceCurve-class}}, \code{\link{kStepEstimate-class}} }
%\examples{}
\concept{one-step estimator}
\concept{estimator}
Modified: branches/robast-0.6/pkg/RobAStBase/man/optIC.Rd
===================================================================
--- branches/robast-0.6/pkg/RobAStBase/man/optIC.Rd 2008-07-28 12:27:14 UTC (rev 131)
+++ branches/robast-0.6/pkg/RobAStBase/man/optIC.Rd 2008-07-31 15:07:01 UTC (rev 132)
@@ -47,4 +47,4 @@
}
\concept{robust influence curve}
\concept{influence curve}
-\keyword{}
+\keyword{robust}
Modified: branches/robast-0.6/pkg/RobLox/R/rlOptIC.R
===================================================================
--- branches/robast-0.6/pkg/RobLox/R/rlOptIC.R 2008-07-28 12:27:14 UTC (rev 131)
+++ branches/robast-0.6/pkg/RobLox/R/rlOptIC.R 2008-07-31 15:07:01 UTC (rev 132)
@@ -18,12 +18,25 @@
biastype = symmetricBias(),
normW = NormType())
+ modIC <- function(L2Fam, IC){
+# if(is(L2Fam, "L2LocationFamily") && is(distribution(L2Fam), "Norm")){
+ if(is(distribution(L2Fam), "Norm")){
+ CallL2New <- call("NormLocationFamily",
+ mean = main(L2Fam))
+ CallL2Fam(IC) <- CallL2New
+ return(IC)
+ }else{
+ stop("'L2Fam' is not compatible with 'CallL2Fam' of 'IC'!")
+ }
+ }
+
return(generateIC(neighbor = ContNeighborhood(radius = r),
L2Fam = NormLocationFamily(mean = mean, sd = sd),
res = list(A = as.matrix(A), a = 0, b = b, d = NULL,
risk = list(asMSE = A, asBias = b, asCov = A - r^2*b^2),
info = c("rlOptIC", "optimally robust IC for AL estimators and 'asMSE'"),
- w = w, biastype = symmetricBias(), normtype = NormType())))
+ w = w, biastype = symmetricBias(), normtype = NormType(),
+ modifyIC = modIC)))
}else{
return(list(A = A, a = 0, b = b))
}
Modified: branches/robast-0.6/pkg/RobLox/R/rlsOptIC_AL.R
===================================================================
--- branches/robast-0.6/pkg/RobLox/R/rlsOptIC_AL.R 2008-07-28 12:27:14 UTC (rev 131)
+++ branches/robast-0.6/pkg/RobLox/R/rlsOptIC_AL.R 2008-07-31 15:07:01 UTC (rev 132)
@@ -27,7 +27,7 @@
###############################################################################
## computation of a1, a2 and a3
###############################################################################
-.ALrlsGeta1a2a3 <- function(b, a1, a2, a3, inteps=1e-12){
+.ALrlsGeta1a2a3 <- function(b, a1, a2, a3){
integrand1 <- function(x, b, a1, a2, a3){
x^2*.ALrlsGetw(x, b, a1, a2, a3)*dnorm(x)
}
@@ -54,8 +54,27 @@
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
###############################################################################
@@ -78,7 +97,9 @@
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
+ 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,
@@ -120,6 +141,7 @@
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
@@ -129,18 +151,57 @@
if(computeIC){
w <- new("HampelWeight")
clip(w) <- b
- cent(w) <- c(0, a2-1)
+ 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){
+ ICL2Fam <- eval(CallL2Fam(IC))
+# if(is(L2Fam, "L2LocationScaleFamily") && is(distribution(L2Fam), "Norm")){
+ if(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
+ weight(w) <- getweight(w, neighbor = ContNeighborhood(radius = neighborRadius(IC)),
+ biastype = biastype(IC),
+ normW = normtype(IC))
+ A <- sdneu^2*stand(IC)/sdalt^2
+ b <- sdneu*clip(IC)/sdalt
+ mse <- sum(diag(A))
+ r <- neighborRadius(IC)
+ a1 <- A[1, 1]/sdneu^2
+ a3 <- A[2, 2]/sdneu^2
+ a2 <- a[1]/sd/a3 + 1
+ asVar <- sd^2*.ALrlsVar(b = b/sd, 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)
+ addInfo(IC) <- c("modifyIC", "The IC has been modified")
+ addInfo(IC) <- c("modifyIC", "The entries in 'Infos' may be wrong")
+ return(IC)
+ }else{
+ stop("'L2Fam' is not compatible with 'CallL2Fam' of 'IC'!")
+ }
+ }
+
return(generateIC(neighbor = ContNeighborhood(radius = r),
L2Fam = NormLocationScaleFamily(mean = mean, sd = sd),
res = list(A = as.matrix(A), a = a, b = b, d = NULL,
- risk = list(asMSE = mse, asBias = b, trAsCov = mse - r^2*b^2),
+ 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())))
+ w = w, biastype = symmetricBias(), normtype = NormType(),
+ modifyIC = modIC)))
}else{
return(list(A = A, a = a, b = b))
}
Modified: branches/robast-0.6/pkg/RobLox/R/roblox.R
===================================================================
--- branches/robast-0.6/pkg/RobLox/R/roblox.R 2008-07-28 12:27:14 UTC (rev 131)
+++ branches/robast-0.6/pkg/RobLox/R/roblox.R 2008-07-31 15:07:01 UTC (rev 132)
@@ -113,9 +113,7 @@
}
.kstep.sc <- function(x, initial.est, A, a, b, mean, k){
est <- .onestep.sc(x = x, initial.est = initial.est, A = A, a = a, b = b, mean = mean)
- if(k == 1){
- return(est)
- }else{
+ if(k > 1){
for(i in 2:k){
A <- est^2*A/initial.est^2
a <- est*a/initial.est
@@ -123,8 +121,12 @@
initial.est <- est
est <- .onestep.sc(x = x, initial.est = est, A = A, a = a, b = b, mean = mean)
}
- return(est)
}
+ A <- est^2*A/initial.est^2
+ a <- est*a/initial.est
+ b <- est*b/initial.est
+
+ return(list(est = est, A = A, a = a, b = b))
}
.onestep.locsc <- function(x, initial.est, A1, A2, a, b){
mean <- initial.est[1]
@@ -137,9 +139,7 @@
}
.kstep.locsc <- function(x, initial.est, A1, A2, a, b, mean, k){
est <- .onestep.locsc(x = x, initial.est = initial.est, A1 = A1, A2 = A2, a = a, b = b)
- if(k == 1){
- return(est)
- }else{
+ if(k > 1){
for(i in 2:k){
A1 <- est[2]^2*A1/initial.est[2]^2
A2 <- est[2]^2*A2/initial.est[2]^2
@@ -148,8 +148,17 @@
initial.est <- est
est <- .onestep.locsc(x = x, initial.est = est, A1 = A1, A2 = A2, a = a, b = b)
}
- return(est)
}
+ A1 <- est[2]^2*A1/initial.est[2]^2
+ A2 <- est[2]^2*A2/initial.est[2]^2
+ a <- est[2]*a/initial.est[2]
+ b <- est[2]*b/initial.est[2]
+ a1 <- A1/est[2]^2
+ a3 <- A2/est[2]^2
+ a2 <- a[2]/est[2]/a3 + 1
+ asVar <- est[2]^2*.ALrlsVar(b = b, a1 = a1, a2 = a2, a3 = a3)
+
+ return(list(est = est, A1 = A1, A2 = A2, a = a, b = b, asvar = asVar))
}
@@ -236,30 +245,74 @@
mse <- A1 + A2
}
robEst <- .kstep.locsc(x = x, initial.est = c(mean, sd), A1 = A1, A2 = A2, a = a, b = b, k = k)
- names(robEst) <- c("mean", "sd")
+ names(robEst$est) <- c("mean", "sd")
Info.matrix <- matrix(c("roblox",
paste("optimally robust estimate for contamination 'eps' =", round(eps, 3),
"and 'asMSE'")),
ncol = 2, dimnames = list(NULL, c("method", "message")))
if(returnIC){
w <- new("HampelWeight")
- clip(w) <- b
- cent(w) <- a/A2
- stand(w) <- diag(c(A1, A2))
+ clip(w) <- robEst$b
+ cent(w) <- robEst$a/robEst$A2
+ stand(w) <- diag(c(robEst$A1, robEst$A2))
weight(w) <- getweight(w, neighbor = ContNeighborhood(radius = r),
biastype = symmetricBias(),
normW = NormType())
+ mse <- robEst$A1 + robEst$A2
+ modIC <- function(L2Fam, IC){
+ ICL2Fam <- eval(CallL2Fam(IC))
+# if(is(L2Fam, "L2LocationScaleFamily") && is(distribution(L2Fam), "Norm")){
+ if(is(distribution(L2Fam), "Norm")){
+ sdneu <- main(L2Fam)[2]
+ sdalt <- main(ICL2Fam)[2]
+ r <- neighborRadius(IC)
+ w <- weight(IC)
+ clip(w) <- sdneu*clip(w)/sdalt
+ cent(w) <- sdalt*cent(w)/sdneu
+ stand(w) <- sdneu^2*stand(w)/sdalt^2
+ weight(w) <- getweight(w, neighbor = ContNeighborhood(radius = r),
+ biastype = biastype(IC),
+ normW = normtype(IC))
+ A <- sdneu^2*stand(IC)/sdalt^2
+ b <- sdneu*clip(IC)/sdalt
+ a <- sdneu*cent(IC)/sdalt
+ mse <- sum(diag(A))
+ a1 <- A[1, 1]/sdneu^2
+ a3 <- A[2, 2]/sdneu^2
+ a2 <- a[2]/sdneu/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)
+ addInfo(IC) <- c("modifyIC", "The IC has been modified")
+ addInfo(IC) <- c("modifyIC", "The entries in 'Infos' may be wrong")
+ return(IC)
+ }else{
+ stop("'L2Fam' is not compatible with 'CallL2Fam' of 'IC'!")
+ }
+ }
IC1 <- generateIC(neighbor = ContNeighborhood(radius = r),
- L2Fam = NormLocationScaleFamily(mean = mean, sd = sd),
- res = list(A = diag(c(A1, A2)), a = a, b = b, d = NULL,
- risk = list(asMSE = mse, asBias = b, asCov = mse - r^2*b^2),
+ L2Fam = NormLocationScaleFamily(mean = robEst$est[1], sd = robEst$est[2]),
+ res = list(A = diag(c(robEst$A1, robEst$A2)), a = robEst$a,
+ b = robEst$b, d = NULL,
+ risk = list(asMSE = mse, asBias = robEst$b,
+ trAsCov = mse - r^2*robEst$b^2,
+ asCov = robEst$asVar),
info = c("roblox", "optimally robust IC for AL estimators and 'asMSE'"),
- w = w, biastype = symmetricBias(), normtype = NormType()))
+ w = w, biastype = symmetricBias(), normtype = NormType(),
+ modifyIC = modIC))
return(new("kStepEstimate", name = "Optimally robust estimate",
- estimate = robEst, steps = k, pIC = IC1, Infos = Info.matrix))
+ estimate = robEst$est, samplesize = length(x), asvar = robEst$asvar,
+ asbias = r*robEst$b, steps = k, pIC = IC1, Infos = Info.matrix))
}else
return(new("kStepEstimate", name = "Optimally robust estimate",
- estimate = robEst, steps = k, pIC = NULL, Infos = Info.matrix))
+ estimate = robEst$est, samplesize = length(x), asvar = robEst$asvar,
+ asbias = r*robEst$b, steps = k, pIC = NULL, Infos = Info.matrix))
}else{
sqrtn <- sqrt(length(x))
rlo <- sqrtn*eps.lower
@@ -296,7 +349,7 @@
}
}
robEst <- .kstep.locsc(x = x, initial.est = c(mean, sd), A1 = A1, A2 = A2, a = a, b = b, k = k)
- names(robEst) <- c("mean", "sd")
+ names(robEst$est) <- c("mean", "sd")
Info.matrix <- matrix(c(rep("roblox", 3),
paste("radius-minimax estimate for contamination interval [",
round(eps.lower, 3), ", ", round(eps.upper, 3), "]", sep = ""),
@@ -305,18 +358,60 @@
ncol = 2, dimnames = list(NULL, c("method", "message")))
if(returnIC){
w <- new("HampelWeight")
- clip(w) <- b
- cent(w) <- a/A2
- stand(w) <- diag(c(A1, A2))
+ clip(w) <- robEst$b
+ cent(w) <- robEst$a/robEst$A2
+ stand(w) <- diag(c(robEst$A1, robEst$A2))
weight(w) <- getweight(w, neighbor = ContNeighborhood(radius = r),
biastype = symmetricBias(),
normW = NormType())
+ mse <- robEst$A1 + robEst$A2
+ modIC <- function(L2Fam, IC){
+ ICL2Fam <- eval(CallL2Fam(IC))
+# if(is(L2Fam, "L2LocationScaleFamily") && is(distribution(L2Fam), "Norm")){
+ if(is(distribution(L2Fam), "Norm")){
+ sdneu <- main(L2Fam)[2]
+ sdalt <- main(ICL2Fam)[2]
+ r <- neighborRadius(IC)
+ w <- weight(IC)
+ clip(w) <- sdneu*clip(w)/sdalt
+ cent(w) <- sdalt*cent(w)/sdneu
+ stand(w) <- sdneu^2*stand(w)/sdalt^2
+ weight(w) <- getweight(w, neighbor = ContNeighborhood(radius = r),
+ biastype = biastype(IC),
+ normW = normtype(IC))
+ A <- sdneu^2*stand(IC)/sdalt^2
+ b <- sdneu*clip(IC)/sdalt
+ a <- sdneu*cent(IC)/sdalt
+ mse <- sum(diag(A))
+ a1 <- A[1, 1]/sdneu^2
+ a3 <- A[2, 2]/sdneu^2
+ a2 <- a[2]/sdneu/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)
+ addInfo(IC) <- c("modifyIC", "The IC has been modified")
+ addInfo(IC) <- c("modifyIC", "The entries in 'Infos' may be wrong")
+ return(IC)
+ }else{
+ stop("'L2Fam' is not compatible with 'CallL2Fam' of 'IC'!")
+ }
+ }
IC1 <- generateIC(neighbor = ContNeighborhood(radius = r),
- L2Fam = NormLocationScaleFamily(mean = mean, sd = sd),
- res = list(A = diag(c(A1, A2)), a = a, b = b, d = NULL,
- risk = list(asMSE = mse, asBias = b, asCov = mse - r^2*b^2),
+ L2Fam = NormLocationScaleFamily(mean = robEst$est[1], sd = robEst$est[2]),
+ res = list(A = diag(c(robEst$A1, robEst$A2)), a = robEst$a,
+ b = robEst$b, d = NULL,
+ risk = list(asMSE = mse, asBias = robEst$b,
+ trAsCov = mse - r^2*robEst$b^2,
+ asCov = robEst$asvar),
info = c("roblox", "optimally robust IC for AL estimators and 'asMSE'"),
- w = w, biastype = symmetricBias(), normtype = NormType()))
+ w = w, biastype = symmetricBias(), normtype = NormType(),
+ modifyIC = modIC))
Infos(IC1) <- matrix(c(rep("roblox", 3),
paste("radius-minimax IC for contamination interval [",
round(eps.lower, 3), ", ", round(eps.upper, 3), "]", sep = ""),
@@ -324,10 +419,12 @@
paste("maximum MSE-inefficiency: ", round(ineff, 3), sep = "")),
ncol = 2, dimnames = list(NULL, c("method", "message")))
return(new("kStepEstimate", name = "Optimally robust estimate",
- estimate = robEst, steps = k, pIC = IC1, Infos = Info.matrix))
+ estimate = robEst$est, samplesize = length(x), asvar = robEst$asvar,
+ asbias = r*robEst$b, steps = k, pIC = IC1, Infos = Info.matrix))
}else
return(new("kStepEstimate", name = "Optimally robust estimate",
- estimate = robEst, steps = k, pIC = NULL, Infos = Info.matrix))
+ estimate = robEst$est, samplesize = length(x), asvar = robEst$asvar,
+ asbias = r*robEst$b, steps = k, pIC = NULL, Infos = Info.matrix))
}
}else{
if(missing(mean)){
@@ -364,19 +461,35 @@
cent(w) <- 0
stand(w) <- as.matrix(A)
weight(w) <- getweight(w, neighbor = ContNeighborhood(radius = r),
- biastype = symmetricBias(),
- normW = NormType())
+ biastype = symmetricBias(),
+ normW = NormType())
+ modIC <- function(L2Fam, IC){
+# if(is(L2Fam, "L2LocationFamily") && is(distribution(L2Fam), "Norm")){
+ if(is(distribution(L2Fam), "Norm")){
+ CallL2New <- call("NormLocationFamily",
+ mean = main(L2Fam))
+ CallL2Fam(IC) <- CallL2New
+ return(IC)
+ }else{
+ stop("'L2Fam' is not compatible with 'CallL2Fam' of 'IC'!")
+ }
+ }
IC1 <- generateIC(neighbor = ContNeighborhood(radius = r),
- L2Fam = NormLocationFamily(mean = mean, sd = sd),
+ L2Fam = NormLocationFamily(mean = robEst, sd = sd),
res = list(A = as.matrix(A), a = 0, b = b, d = NULL,
- risk = list(asMSE = A, asBias = b, asCov = b^2),
+ risk = list(asMSE = A, asBias = b, asCov = A-r^2*b^2),
info = c("roblox", "optimally robust IC for AL estimators and 'asMSE'"),
- w = w, biastype = symmetricBias(), normtype = NormType()))
+ w = w, biastype = symmetricBias(), normtype = NormType(),
+ modifyIC = modIC))
return(new("kStepEstimate", name = "Optimally robust estimate",
- estimate = robEst, steps = k, pIC = IC1, Infos = Info.matrix))
+ estimate = robEst, samplesize = length(x),
+ asvar = as.matrix(A-r^2*b^2),
+ asbias = r*b, steps = k, pIC = IC1, Infos = Info.matrix))
}else
return(new("kStepEstimate", name = "Optimally robust estimate",
- estimate = robEst, steps = k, pIC = NULL, Infos = Info.matrix))
+ estimate = robEst, samplesize = length(x),
+ asvar = as.matrix(A-r^2*b^2),
+ asbias = r*b, steps = k, pIC = NULL, Infos = Info.matrix))
}else{
sqrtn <- sqrt(length(x))
rlo <- sqrtn*eps.lower
@@ -420,12 +533,24 @@
weight(w) <- getweight(w, neighbor = ContNeighborhood(radius = r),
biastype = symmetricBias(),
normW = NormType())
+ modIC <- function(L2Fam, IC){
+# if(is(L2Fam, "L2LocationFamily") && is(distribution(L2Fam), "Norm")){
+ if(is(distribution(L2Fam), "Norm")){
+ CallL2New <- call("NormLocationFamily",
+ mean = main(L2Fam))
+ CallL2Fam(IC) <- CallL2New
+ return(IC)
+ }else{
+ stop("'L2Fam' is not compatible with 'CallL2Fam' of 'IC'!")
+ }
+ }
IC1 <- generateIC(neighbor = ContNeighborhood(radius = r),
- L2Fam = NormLocationFamily(mean = mean, sd = sd),
+ L2Fam = NormLocationFamily(mean = robEst, sd = sd),
res = list(A = as.matrix(A), a = 0, b = b, d = NULL,
- risk = list(asMSE = A, asBias = b, asCov = b^2),
+ risk = list(asMSE = A, asBias = b, asCov = A-r^2*b^2),
info = c("roblox", "optimally robust IC for AL estimators and 'asMSE'"),
- w = w, biastype = symmetricBias(), normtype = NormType()))
+ w = w, biastype = symmetricBias(), normtype = NormType(),
+ modifyIC = modIC))
Infos(IC1) <- matrix(c(rep("roblox", 3),
paste("radius-minimax IC for contamination interval [",
round(eps.lower, 3), ", ", round(eps.upper, 3), "]", sep = ""),
@@ -433,10 +558,14 @@
paste("maximum MSE-inefficiency: ", round(ineff, 3), sep = "")),
ncol = 2, dimnames = list(NULL, c("method", "message")))
return(new("kStepEstimate", name = "Optimally robust estimate",
- estimate = robEst, steps = k, pIC = IC1, Infos = Info.matrix))
+ estimate = robEst, samplesize = length(x),
+ asvar = as.matrix(A-r^2*b^2),
+ asbias = r*b, steps = k, pIC = IC1, Infos = Info.matrix))
}else
return(new("kStepEstimate", name = "Optimally robust estimate",
- estimate = robEst, steps = k, pIC = NULL, Infos = Info.matrix))
+ estimate = robEst, samplesize = length(x),
+ asvar = as.matrix(A-r^2*b^2),
+ asbias = r*b, steps = k, pIC = NULL, Infos = Info.matrix))
}
}
if(missing(sd)){
@@ -467,30 +596,66 @@
b <- sd*.getb.sc(r)
}
robEst <- .kstep.sc(x = x, initial.est = sd, A = A, a = a, b = b, mean = mean, k = k)
- names(robEst) <- "sd"
+ names(robEst$est) <- "sd"
Info.matrix <- matrix(c("roblox",
paste("optimally robust estimate for contamination 'eps' =", round(eps, 3),
"and 'asMSE'")),
ncol = 2, dimnames = list(NULL, c("method", "message")))
if(returnIC){
w <- new("HampelWeight")
- clip(w) <- b
- cent(w) <- a/A
- stand(w) <- as.matrix(A)
+ clip(w) <- robEst$b
+ cent(w) <- robEst$a/robEst$A
+ stand(w) <- as.matrix(robEst$A)
weight(w) <- getweight(w, neighbor = ContNeighborhood(radius = r),
biastype = symmetricBias(),
normW = NormType())
+ modIC <- function(L2Fam, IC){
+ ICL2Fam <- eval(CallL2Fam(IC))
+# if(is(L2Fam, "L2ScaleFamily") && is(distribution(L2Fam), "Norm")){
+ if(is(distribution(L2Fam), "Norm")){
+ sdneu <- main(L2Fam)
+ sdalt <- main(ICL2Fam)
+ r <- neighborRadius(IC)
+ w <- weight(IC)
+ clip(w) <- sdneu*clip(w)/sdalt
+ cent(w) <- sdalt*cent(w)/sdneu
+ stand(w) <- sdneu^2*stand(w)/sdalt^2
+ weight(w) <- getweight(w, neighbor = ContNeighborhood(radius = r),
+ biastype = biastype(IC),
+ normW = normtype(IC))
+ A <- sdneu^2*stand(IC)/sdalt^2
+ b <- sdneu*clip(IC)/sdalt
+ res <- list(A = A, a = sdneu*cent(IC)/sdalt, b = b, d = NULL,
+ risk = list(asMSE = A, asBias = b, asCov = A-r^2*b^2),
+ info = Infos(IC), w = w,
+ normtype = normtype(IC), biastype = biastype(IC),
+ modifyIC = modifyIC(IC))
+ IC <- generateIC(neighbor = ContNeighborhood(radius = r),
+ L2Fam = L2Fam, res = res)
+ addInfo(IC) <- c("modifyIC", "The IC has been modified")
+ addInfo(IC) <- c("modifyIC", "The entries in 'Infos' may be wrong")
+ return(IC)
+ }else{
+ stop("'L2Fam' is not compatible with 'CallL2Fam' of 'IC'!")
+ }
+ }
IC1 <- generateIC(neighbor = ContNeighborhood(radius = r),
- L2Fam = NormScaleFamily(mean = mean, sd = sd),
- res = list(A = as.matrix(A), a = a, b = b, d = NULL,
- risk = list(asMSE = A, asBias = b, asCov = b^2),
+ L2Fam = NormScaleFamily(mean = mean, sd = robEst$est),
+ res = list(A = as.matrix(robEst$A), a = robEst$a, b = robEst$b, d = NULL,
+ risk = list(asMSE = robEst$A, asBias = robEst$b,
+ asCov = robEst$A-r^2*robEst$b^2),
info = c("roblox", "optimally robust IC for AL estimators and 'asMSE'"),
- w = w, biastype = symmetricBias(), normtype = NormType()))
+ w = w, biastype = symmetricBias(), normtype = NormType(),
+ modifyIC = modIC))
return(new("kStepEstimate", name = "Optimally robust estimate",
- estimate = robEst, steps = k, pIC = IC1, Infos = Info.matrix))
+ estimate = robEst$est, samplesize = length(x),
+ asvar = as.matrix(robEst$A-r^2*robEst$b^2),
+ asbias = r*robEst$b, steps = k, pIC = IC1, Infos = Info.matrix))
}else
return(new("kStepEstimate", name = "Optimally robust estimate",
- estimate = robEst, steps = k, pIC = NULL, Infos = Info.matrix))
+ estimate = robEst$est, samplesize = length(x),
+ asvar = as.matrix(robEst$A-r^2*robEst$b^2),
+ asbias = r*robEst$b, steps = k, pIC = NULL, Infos = Info.matrix))
}else{
sqrtn <- sqrt(length(x))
rlo <- sqrtn*eps.lower
@@ -521,7 +686,7 @@
}
}
robEst <- .kstep.sc(x = x, initial.est = sd, A = A, a = a, b = b, mean = mean, k = k)
- names(robEst) <- "sd"
+ names(robEst$est) <- "sd"
Info.matrix <- matrix(c(rep("roblox", 3),
paste("radius-minimax estimate for contamination interval [",
round(eps.lower, 3), ", ", round(eps.upper, 3), "]", sep = ""),
@@ -530,18 +695,50 @@
ncol = 2, dimnames = list(NULL, c("method", "message")))
if(returnIC){
w <- new("HampelWeight")
- clip(w) <- b
- cent(w) <- a/A
- stand(w) <- as.matrix(A)
+ clip(w) <- robEst$b
+ cent(w) <- robEst$a/robEst$A
+ stand(w) <- as.matrix(robEst$A)
weight(w) <- getweight(w, neighbor = ContNeighborhood(radius = r),
biastype = symmetricBias(),
normW = NormType())
+ modIC <- function(L2Fam, IC){
+ ICL2Fam <- eval(CallL2Fam(IC))
+# if(is(L2Fam, "L2ScaleFamily") && is(distribution(L2Fam), "Norm")){
+ if(is(distribution(L2Fam), "Norm")){
+ sdneu <- main(L2Fam)
+ sdalt <- main(ICL2Fam)
+ r <- neighborRadius(IC)
+ w <- weight(IC)
+ clip(w) <- sdneu*clip(w)/sdalt
+ cent(w) <- sdalt*cent(w)/sdneu
+ stand(w) <- sdneu^2*stand(w)/sdalt^2
+ weight(w) <- getweight(w, neighbor = ContNeighborhood(radius = r),
+ biastype = biastype(IC),
+ normW = normtype(IC))
+ A <- sdneu^2*stand(IC)/sdalt^2
+ b <- sdneu*clip(IC)/sdalt
+ res <- list(A = A, a = sdneu*cent(IC)/sdalt, b = b, d = NULL,
+ risk = list(asMSE = A, asBias = b, asCov = A-r^2*b^2),
+ info = Infos(IC), w = w,
+ normtype = normtype(IC), biastype = biastype(IC),
+ modifyIC = modifyIC(IC))
+ IC <- generateIC(neighbor = ContNeighborhood(radius = r),
+ L2Fam = L2Fam, res = res)
+ addInfo(IC) <- c("modifyIC", "The IC has been modified")
+ addInfo(IC) <- c("modifyIC", "The entries in 'Infos' may be wrong")
+ return(IC)
+ }else{
+ stop("'L2Fam' is not compatible with 'CallL2Fam' of 'IC'!")
+ }
+ }
IC1 <- generateIC(neighbor = ContNeighborhood(radius = r),
- L2Fam = NormScaleFamily(mean = mean, sd = sd),
- res = list(A = as.matrix(A), a = a, b = b, d = NULL,
- risk = list(asMSE = A, asBias = b, asCov = b^2),
+ L2Fam = NormScaleFamily(mean = mean, sd = robEst$est),
+ res = list(A = as.matrix(robEst$A), a = robEst$a, b = robEst$b, d = NULL,
+ risk = list(asMSE = robEst$A, asBias = robEst$b,
+ asCov = robEst$A-r^2*robEst$b^2),
info = c("roblox", "optimally robust IC for AL estimators and 'asMSE'"),
- w = w, biastype = symmetricBias(), normtype = NormType()))
+ w = w, biastype = symmetricBias(), normtype = NormType(),
+ modifyIC = modIC))
Infos(IC1) <- matrix(c(rep("roblox", 3),
paste("radius-minimax IC for contamination interval [",
round(eps.lower, 3), ", ", round(eps.upper, 3), "]", sep = ""),
@@ -549,10 +746,14 @@
paste("maximum MSE-inefficiency: ", round(ineff, 3), sep = "")),
ncol = 2, dimnames = list(NULL, c("method", "message")))
return(new("kStepEstimate", name = "Optimally robust estimate",
- estimate = robEst, steps = k, pIC = IC1, Infos = Info.matrix))
+ estimate = robEst$est, samplesize = length(x),
+ asvar = as.matrix(robEst$A-r^2*robEst$b^2),
+ asbias = r*robEst$b, steps = k, pIC = IC1, Infos = Info.matrix))
}else
return(new("kStepEstimate", name = "Optimally robust estimate",
- estimate = robEst, steps = k, pIC = NULL, Infos = Info.matrix))
+ estimate = robEst$est, samplesize = length(x),
+ asvar = as.matrix(robEst$A-r^2*robEst$b^2),
+ asbias = r*robEst$b, steps = k, pIC = NULL, Infos = Info.matrix))
}
}
}
Modified: branches/robast-0.6/pkg/RobLox/R/rowRoblox.R
===================================================================
--- branches/robast-0.6/pkg/RobLox/R/rowRoblox.R 2008-07-28 12:27:14 UTC (rev 131)
+++ branches/robast-0.6/pkg/RobLox/R/rowRoblox.R 2008-07-31 15:07:01 UTC (rev 132)
@@ -22,9 +22,7 @@
}
.kstep.sc.matrix <- function(x, initial.est, A, a, b, mean, k){
est <- .onestep.sc.matrix(x = x, initial.est = initial.est, A = A, a = a, b = b, mean = mean)
- if(k == 1){
- return(est)
- }else{
+ if(k > 1){
for(i in 2:k){
A <- est^2*A/initial.est^2
a <- est*a/initial.est
@@ -32,8 +30,11 @@
initial.est <- est
est <- .onestep.sc.matrix(x = x, initial.est = est, A = A, a = a, b = b, mean = mean)
}
- return(est)
}
+ A <- est^2*A/initial.est^2
+ a <- est*a/initial.est
+ b <- est*b/initial.est
+ return(list(est = est, A = A, a = a, b = b))
}
.onestep.locsc.matrix <- function(x, initial.est, A1, A2, a, b){
mean <- initial.est[,1]
@@ -48,9 +49,7 @@
}
.kstep.locsc.matrix <- function(x, initial.est, A1, A2, a, b, mean, k){
est <- .onestep.locsc.matrix(x = x, initial.est = initial.est, A1 = A1, A2 = A2, a = a, b = b)
- if(k == 1){
- return(est)
- }else{
+ if(k > 1){
for(i in 2:k){
A1 <- est[,2]^2*A1/initial.est[,2]^2
A2 <- est[,2]^2*A2/initial.est[,2]^2
@@ -59,8 +58,13 @@
initial.est <- est
est <- .onestep.locsc.matrix(x = x, initial.est = est, A1 = A1, A2 = A2, a = a, b = b)
}
- return(est)
}
+ A1 <- est[,2]^2*A1/initial.est[,2]^2
+ A2 <- est[,2]^2*A2/initial.est[,2]^2
+ a <- est[,2]*a/initial.est[,2]
+ b <- est[,2]*b/initial.est[,2]
+
+ return(list(est = est, A1 = A1, A2 = A2, a = a, b = b))
}
@@ -148,13 +152,18 @@
}
robEst <- .kstep.locsc.matrix(x = x, initial.est = cbind(mean, sd),
A1 = A1, A2 = A2, a = a, b = b, k = k)
- colnames(robEst) <- c("mean", "sd")
+ colnames(robEst$est) <- c("mean", "sd")
Info.matrix <- matrix(c("roblox",
paste("optimally robust estimates for contamination 'eps' =", round(eps, 3),
"and 'asMSE'")),
ncol = 2, dimnames = list(NULL, c("method", "message")))
return(new("kStepEstimate", name = "Optimally robust estimate",
- estimate = robEst, steps = k, pIC = NULL, Infos = Info.matrix))
+ estimate = robEst$est, samplesize = ncol(x), steps = k,
+ pIC = NULL, Infos = Info.matrix))
+## we need a class like "list of estimates" to set asvar and asbias consistently ...
+# return(new("kStepEstimate", name = "Optimally robust estimate",
+# estimate = robEst$est, samplesize = ncol(x), asvar = NULL,
+# asbias = r*robEst$b, steps = k, pIC = NULL, Infos = Info.matrix))
}else{
sqrtn <- sqrt(ncol(x))
rlo <- sqrtn*eps.lower
@@ -192,7 +201,7 @@
}
robEst <- .kstep.locsc.matrix(x = x, initial.est = cbind(mean, sd),
A1 = A1, A2 = A2, a = a, b = b, k = k)
- colnames(robEst) <- c("mean", "sd")
+ colnames(robEst$est) <- c("mean", "sd")
Info.matrix <- matrix(c(rep("roblox", 3),
paste("radius-minimax estimates for contamination interval [",
round(eps.lower, 3), ", ", round(eps.upper, 3), "]", sep = ""),
@@ -200,7 +209,12 @@
paste("maximum MSE-inefficiency: ", round(ineff[1], 3), sep = "")),
ncol = 2, dimnames = list(NULL, c("method", "message")))
return(new("kStepEstimate", name = "Optimally robust estimate",
- estimate = robEst, steps = k, pIC = NULL, Infos = Info.matrix))
+ estimate = robEst$est, samplesize = ncol(x), steps = k,
+ pIC = NULL, Infos = Info.matrix))
+## we need a class like "list of estimates" to set asvar and asbias consistently ...
+# return(new("kStepEstimate", name = "Optimally robust estimate",
+# estimate = robEst$est, samplesize = ncol(x), asvar = NULL,
+# asbias = r*robEst$b, steps = k, pIC = NULL, Infos = Info.matrix))
}
}else{
if(missing(mean)){
@@ -219,6 +233,8 @@
stop("'initial.est' has wrong dimension")
mean <- initial.est
}
+ if(length(sd) == 1)
+ sd <- rep(sd, length(mean))
if(!missing(eps)){
r <- sqrt(ncol(x))*eps
@@ -236,7 +252,12 @@
"and 'asMSE'")),
ncol = 2, dimnames = list(NULL, c("method", "message")))
return(new("kStepEstimate", name = "Optimally robust estimate",
- estimate = robEst, steps = k, pIC = NULL, Infos = Info.matrix))
+ estimate = robEst, samplesize = ncol(x), steps = k,
+ pIC = NULL, Infos = Info.matrix))
+## we need a class like "list of estimates" to set asvar and asbias consistently ...
+# return(new("kStepEstimate", name = "Optimally robust estimate",
+# estimate = robEst$est, samplesize = ncol(x), asvar = as.matrix(A - r^2*b^2),
+# asbias = r*b, steps = k, pIC = NULL, Infos = Info.matrix))
}else{
sqrtn <- sqrt(ncol(x))
rlo <- sqrtn*eps.lower
@@ -273,7 +294,12 @@
paste("maximum MSE-inefficiency: ", round(ineff[1], 3), sep = "")),
ncol = 2, dimnames = list(NULL, c("method", "message")))
return(new("kStepEstimate", name = "Optimally robust estimate",
- estimate = robEst, steps = k, pIC = NULL, Infos = Info.matrix))
+ estimate = robEst, samplesize = ncol(x), steps = k,
+ pIC = NULL, Infos = Info.matrix))
+## we need a class like "list of estimates" to set asvar and asbias consistently ...
+# return(new("kStepEstimate", name = "Optimally robust estimate",
+# estimate = robEst$est, samplesize = ncol(x), asvar = as.matrix(A - r^2*b^2),
+# asbias = r*b, steps = k, pIC = NULL, Infos = Info.matrix))
}
}
if(missing(sd)){
@@ -308,14 +334,19 @@
a <- sd*.geta.sc(r)
b <- sd*.getb.sc(r)
}
- robEst <- as.matrix(.kstep.sc.matrix(x = x, initial.est = sd, A = A, a = a, b = b, mean = mean, k = k))
- colnames(robEst) <- "sd"
+ robEst <- .kstep.sc.matrix(x = x, initial.est = sd, A = A, a = a, b = b, mean = mean, k = k)
+ names(robEst$est) <- "sd"
Info.matrix <- matrix(c("roblox",
paste("optimally robust estimates for contamination 'eps' =", round(eps, 3),
"and 'asMSE'")),
ncol = 2, dimnames = list(NULL, c("method", "message")))
return(new("kStepEstimate", name = "Optimally robust estimate",
- estimate = robEst, steps = k, pIC = NULL, Infos = Info.matrix))
+ estimate = as.matrix(robEst$est), samplesize = ncol(x), steps = k,
+ pIC = NULL, Infos = Info.matrix))
+## we need a class like "list of estimates" to set asvar and asbias consistently ...
+# return(new("kStepEstimate", name = "Optimally robust estimate",
+# estimate = robEst$est, samplesize = ncol(x), asvar = as.matrix(robEst$A - r^2*robEst$b^2),
+# asbias = r*robEst$b, steps = k, pIC = NULL, Infos = Info.matrix))
}else{
sqrtn <- sqrt(ncol(x))
rlo <- sqrtn*eps.lower
@@ -345,8 +376,8 @@
ineff <- (A - b^2*(r^2 - rlo^2))/Alo
}
}
- robEst <- as.matrix(.kstep.sc.matrix(x = x, initial.est = sd, A = A, a = a, b = b, mean = mean, k = k))
- colnames(robEst) <- "sd"
+ robEst <- .kstep.sc.matrix(x = x, initial.est = sd, A = A, a = a, b = b, mean = mean, k = k)
+ names(robEst$est) <- "sd"
Info.matrix <- matrix(c(rep("roblox", 3),
paste("radius-minimax estimates for contamination interval [",
round(eps.lower, 3), ", ", round(eps.upper, 3), "]", sep = ""),
@@ -354,7 +385,12 @@
paste("maximum MSE-inefficiency: ", round(ineff[1], 3), sep = "")),
ncol = 2, dimnames = list(NULL, c("method", "message")))
return(new("kStepEstimate", name = "Optimally robust estimate",
- estimate = robEst, steps = k, pIC = NULL, Infos = Info.matrix))
+ estimate = as.matrix(robEst$est), samplesize = ncol(x), steps = k,
+ pIC = NULL, Infos = Info.matrix))
+## we need a class like "list of estimates" to set asvar and asbias consistently ...
+# return(new("kStepEstimate", name = "Optimally robust estimate",
+# estimate = robEst$est, samplesize = ncol(x), asvar = as.matrix(robEst$A - r^2*robEst$b^2),
+# asbias = r*robEst$b, steps = k, pIC = NULL, Infos = Info.matrix))
}
}
}
Modified: branches/robast-0.6/pkg/RobLox/R/rsOptIC.R
===================================================================
--- branches/robast-0.6/pkg/RobLox/R/rsOptIC.R 2008-07-28 12:27:14 UTC (rev 131)
+++ branches/robast-0.6/pkg/RobLox/R/rsOptIC.R 2008-07-31 15:07:01 UTC (rev 132)
@@ -31,7 +31,7 @@
tol = .Machine$double.eps^0.5, r = r, z = z)$root
iter <- 0
- repeat{
+ repeat{
iter <- iter + 1
if(iter > itmax)
stop("Algorithm did not converge => increase 'itmax'!")
@@ -64,18 +64,49 @@
if(computeIC){
w <- new("HampelWeight")
clip(w) <- b
- cent(w) <- z-1
+ cent(w) <- (z-1)/sd
stand(w) <- as.matrix(A)
weight(w) <- getweight(w, neighbor = ContNeighborhood(radius = r),
biastype = symmetricBias(),
normW = NormType())
+ modIC <- function(L2Fam, IC){
+ ICL2Fam <- eval(CallL2Fam(IC))
+# if(is(L2Fam, "L2ScaleFamily") && is(distribution(L2Fam), "Norm")){
+ if(is(distribution(L2Fam), "Norm")){
+ sdneu <- main(L2Fam)
+ sdalt <- main(ICL2Fam)
+ w <- weight(IC)
+ clip(w) <- sdneu*clip(w)/sdalt
+ cent(w) <- sdalt*cent(w)/sdneu
+ stand(w) <- sdneu^2*stand(w)/sdalt^2
+ weight(w) <- getweight(w, neighbor = ContNeighborhood(radius = neighborRadius(IC)),
+ biastype = biastype(IC),
+ normW = normtype(IC))
+ A <- sdneu^2*stand(IC)/sdalt^2
+ b <- sdneu*clip(IC)/sdalt
+ res <- list(A = as.matrix(A), a = sdneu*cent(IC)/sdalt, b = b, d = NULL,
+ risk = list(asMSE = A, asBias = b, asCov = A - r^2*b^2),
+ info = Infos(IC), w = w,
+ normtype = normtype(IC), biastype = biastype(IC),
+ modifyIC = modifyIC(IC))
+ IC <- generateIC(neighbor = ContNeighborhood(radius = neighborRadius(IC)),
+ L2Fam = L2Fam, res = res)
+ addInfo(IC) <- c("modifyIC", "The IC has been modified")
+ addInfo(IC) <- c("modifyIC", "The entries in 'Infos' may be wrong")
+ return(IC)
+ }else{
+ stop("'L2Fam' is not compatible with 'CallL2Fam' of 'IC'!")
+ }
+ }
+
return(generateIC(neighbor = ContNeighborhood(radius = r),
L2Fam = NormScaleFamily(sd = sd, mean = mean),
res = list(A = as.matrix(A), a = a, b = b, d = NULL,
risk = list(asMSE = A, asBias = b, asCov = A - r^2*b^2),
info = c("rlOptIC", "optimally robust IC for AL estimators and 'asMSE'"),
- w = w, biastype = symmetricBias(), normtype = NormType())))
+ w = w, biastype = symmetricBias(), normtype = NormType(),
+ modifyIC = modIC)))
}else{
return(list(A = A, a = a, b = b))
}
Modified: branches/robast-0.6/pkg/RobLox/man/rlOptIC.Rd
===================================================================
--- branches/robast-0.6/pkg/RobLox/man/rlOptIC.Rd 2008-07-28 12:27:14 UTC (rev 131)
+++ branches/robast-0.6/pkg/RobLox/man/rlOptIC.Rd 2008-07-31 15:07:01 UTC (rev 132)
@@ -24,7 +24,7 @@
}
\value{
If 'computeIC' is 'TRUE' an object of class \code{"ContIC"} is returned,
- otherwise a list of Lagrane multipliers
+ otherwise a list of Lagrange multipliers
\item{A}{ standardizing constant }
\item{a}{ centering constant; always '= 0' is this symmetric setup }
\item{b}{ optimal clipping bound }
Modified: branches/robast-0.6/pkg/RobLox/man/rlsOptIC.AL.Rd
===================================================================
--- branches/robast-0.6/pkg/RobLox/man/rlsOptIC.AL.Rd 2008-07-28 12:27:14 UTC (rev 131)
+++ branches/robast-0.6/pkg/RobLox/man/rlsOptIC.AL.Rd 2008-07-31 15:07:01 UTC (rev 132)
@@ -37,7 +37,7 @@
}
\value{
If 'computeIC' is 'TRUE' an object of class \code{"ContIC"} is returned,
- otherwise a list of Lagrane multipliers
+ otherwise a list of Lagrange multipliers
\item{A}{ standardizing matrix }
\item{a}{ centering vector }
\item{b}{ optimal clipping bound }
@@ -63,8 +63,8 @@
plot(IC1)
infoPlot(IC1)
-## one-step estimation
-## see also: ?roblox
+## k-step estimation
+## better use function roblox (see ?roblox)
## 1. data: random sample
ind <- rbinom(100, size=1, prob=0.05)
x <- rnorm(100, mean=0, sd=(1-ind) + ind*9)
@@ -73,19 +73,29 @@
median(x)
mad(x)
-## 2. Kolmogorov(-Smirnov) minimum distance estimator
+## 2. Kolmogorov(-Smirnov) minimum distance estimator (default)
## -> we use it as initial estimate for one-step construction
-(est0 <- MDEstimator(x, ParamFamily = NormLocationScaleFamily(), distance = KolmogorovDist))
+(est0 <- MDEstimator(x, ParamFamily = NormLocationScaleFamily()))
-## 3. one-step estimation: radius known
+## 3.1 one-step estimation: radius known
IC1 <- rlsOptIC.AL(r = 0.5, mean = estimate(est0)[1], sd = estimate(est0)[2])
(est1 <- oneStepEstimator(x, IC1, est0))
-## 4. one-step estimation: radius unknown
+## 3.2 k-step estimation: radius known
+## Choose k = 3
+(est2 <- kStepEstimator(x, IC1, est0, steps = 3L))
+
+## 4.1 one-step estimation: radius unknown
## take least favorable radius r = 0.579
## cf. Table 8.1 in Kohl(2005)
IC2 <- rlsOptIC.AL(r = 0.579, mean = estimate(est0)[1], sd = estimate(est0)[2])
-(est2 <- oneStepEstimator(x, IC2, est0))
+(est3 <- oneStepEstimator(x, IC2, est0))
+
+## 4.2 k-step estimation: radius unknown
+## take least favorable radius r = 0.579
+## cf. Table 8.1 in Kohl(2005)
+## choose k = 3
+(est4 <- kStepEstimator(x, IC2, est0, steps = 3L))
}
\concept{normal location and scale}
\concept{influence curve}
Modified: branches/robast-0.6/pkg/RobLox/man/rsOptIC.Rd
===================================================================
--- branches/robast-0.6/pkg/RobLox/man/rsOptIC.Rd 2008-07-28 12:27:14 UTC (rev 131)
+++ branches/robast-0.6/pkg/RobLox/man/rsOptIC.Rd 2008-07-31 15:07:01 UTC (rev 132)
@@ -26,7 +26,7 @@
}
\value{
If 'computeIC' is 'TRUE' an object of class \code{"ContIC"} is returned,
- otherwise a list of Lagrane multipliers
+ otherwise a list of Lagrange multipliers
\item{A}{ standardizing constant }
\item{a}{ centering constant }
\item{b}{ optimal clipping bound }
More information about the Robast-commits
mailing list