[Robast-commits] r1110 - in branches/robast-1.2/pkg/RobAStBase: . R inst man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Aug 8 23:49:35 CEST 2018
Author: ruckdeschel
Date: 2018-08-08 23:49:35 +0200 (Wed, 08 Aug 2018)
New Revision: 1110
Added:
branches/robast-1.2/pkg/RobAStBase/R/combinedICs.R
Modified:
branches/robast-1.2/pkg/RobAStBase/DESCRIPTION
branches/robast-1.2/pkg/RobAStBase/NAMESPACE
branches/robast-1.2/pkg/RobAStBase/R/AllClass.R
branches/robast-1.2/pkg/RobAStBase/R/AllGeneric.R
branches/robast-1.2/pkg/RobAStBase/R/ContIC.R
branches/robast-1.2/pkg/RobAStBase/R/HampIC.R
branches/robast-1.2/pkg/RobAStBase/R/IC.R
branches/robast-1.2/pkg/RobAStBase/R/TotalVarIC.R
branches/robast-1.2/pkg/RobAStBase/R/ddPlot_utils.R
branches/robast-1.2/pkg/RobAStBase/R/generateICfct.R
branches/robast-1.2/pkg/RobAStBase/R/kStepEstimator.R
branches/robast-1.2/pkg/RobAStBase/inst/NEWS
branches/robast-1.2/pkg/RobAStBase/man/ContIC-class.Rd
branches/robast-1.2/pkg/RobAStBase/man/ContIC.Rd
branches/robast-1.2/pkg/RobAStBase/man/HampIC-class.Rd
branches/robast-1.2/pkg/RobAStBase/man/IC-class.Rd
branches/robast-1.2/pkg/RobAStBase/man/TotalVarIC-class.Rd
branches/robast-1.2/pkg/RobAStBase/man/TotalVarIC.Rd
branches/robast-1.2/pkg/RobAStBase/man/outlyingPlotIC.Rd
Log:
[RobAStBase] branch 1.2:
+ force optimal ICs to respect the support of the model distribution
+ and a forgotten no longer used instance of oldmodif in kStepEstimator
+ updated required package versions in DESCRIPTION
+ force optimal ICs to respect the support of the model distribution
+ in kStepEstimator got back from RandVar-evaluation to IC - evaluation
background: updates should be fast (I saw examples with 60s for 3step...
with fast LMs...) -> to this end: [so far things only got worse....]
(a) (for internal purposes) introduce new intermediate S4 class ".fastIC"
(with non-exported generator .fastIC in file combinedICs.R) which is
inbetween class IC and HampIC and has a new slot ".fastFct".
".fastFct" is an optional (= can be NULL) mere function in one argument
which returns the vector-valued IC; this way coordinatewise repeated
checking whether x is in support of distr (and evaluation of the weight)
can be avoided
(b) new slot ".fastFct" is filled automatically for our Hamepl-type
ICs in generators ContIC and TotalVarIC by analogue generateIC.fast.fct
to generateIC.fct in file generateICfct.R.
(c) class .fastIC is intermediate as we need it, too, for non-Hampel type ICs
as arise when either the covariance of our opt-rob IC is singular or
one works with pICs and has to reconstruct full ICs by filling the parts
in the orthogonal complement of Range IC;
(d) to this last issue instead of adding two random variables, as was done
beforehand in kStepEstimator, one uses the new helper function combineOrthPICs
in file combinedICs.R which combines (without checking orthogonality) two
pICs to one full IC by adding the curves (and the fast functions).
(e) in kStepEstimator, we now use evalIC.v, a (sapply-)vectorized version
of evalIC; this is an exported method and has a particular method for
class ".fastIC" which uses slot ".fastFct" instead of the evaluation
of the pIC through evalRandVar ...
(f) generateIC.fct has also been revised: it avoids using random variable
Y(x)/Yi(x) and instead computes them right away from Lambda;
this also has as background that checkIC/makeIC should be enhanced;
ultimately, this enhancement is passed to ROptEst -- idea is to
reuse infrastructure from getInfStand getInfCent which automatically
does symmetry checking ...
Modified: branches/robast-1.2/pkg/RobAStBase/DESCRIPTION
===================================================================
--- branches/robast-1.2/pkg/RobAStBase/DESCRIPTION 2018-08-07 23:53:41 UTC (rev 1109)
+++ branches/robast-1.2/pkg/RobAStBase/DESCRIPTION 2018-08-08 21:49:35 UTC (rev 1110)
@@ -3,7 +3,7 @@
Date: 2018-08-03
Title: Robust Asymptotic Statistics
Description: Base S4-classes and functions for robust asymptotic statistics.
-Depends: R(>= 2.14.0), methods, rrcov, distr(>= 2.5.2), distrEx(>= 2.8.0), distrMod(>= 2.8.0),
+Depends: R(>= 2.14.0), methods, rrcov, distr(>= 2.8.0), distrEx(>= 2.8.0), distrMod(>= 2.8.0),
RandVar(>= 1.1.0)
Suggests: ROptEst(>= 1.1.0), RUnit(>= 0.4.26)
Imports: startupmsg, graphics, grDevices, stats
Modified: branches/robast-1.2/pkg/RobAStBase/NAMESPACE
===================================================================
--- branches/robast-1.2/pkg/RobAStBase/NAMESPACE 2018-08-07 23:53:41 UTC (rev 1109)
+++ branches/robast-1.2/pkg/RobAStBase/NAMESPACE 2018-08-08 21:49:35 UTC (rev 1110)
@@ -20,7 +20,7 @@
"FixRobModel",
"InfRobModel")
exportClasses("InfluenceCurve",
- "IC", "HampIC",
+ "IC", "HampIC", ".fastIC",
"ContIC",
"TotalVarIC")
exportClasses("RobAStControl", "RobWeight", "BoundedWeight",
@@ -44,7 +44,7 @@
"modifyIC",
"generateIC",
"checkIC",
- "evalIC",
+ "evalIC", "evalIC.v",
"clip", "clip<-",
"cent", "cent<-",
"stand", "stand<-",
Modified: branches/robast-1.2/pkg/RobAStBase/R/AllClass.R
===================================================================
--- branches/robast-1.2/pkg/RobAStBase/R/AllClass.R 2018-08-07 23:53:41 UTC (rev 1109)
+++ branches/robast-1.2/pkg/RobAStBase/R/AllClass.R 2018-08-08 21:49:35 UTC (rev 1110)
@@ -127,6 +127,11 @@
return(TRUE)
})
+
+## internal class
+setClass(".fastIC", representation(.fastFct = "OptionalFunction"),
+ prototype(.fastFct = NULL), contains="IC")
+
## HampIC -- common mother class to ContIC and TotalVarIC
setClass("HampIC",
representation(stand = "matrix",
@@ -134,7 +139,7 @@
neighborRadius = "numeric",
weight = "RobWeight",
biastype = "BiasType",
- normtype = "NormType"),
+ normtype = "NormType"),
prototype(name = "IC of total-var or contamination type",
Curve = EuclRandVarList(RealRandVariable(Map = list(function(x){x}),
Domain = Reals())),
@@ -148,7 +153,7 @@
neighborRadius = 0,
biastype = symmetricBias(),
NormType = NormType()),
- contains = "IC",
+ contains = ".fastIC",
validity = function(object){
if(any(object at neighborRadius < 0)) # radius vector?!
stop("'neighborRadius' has to be in [0, Inf]")
Modified: branches/robast-1.2/pkg/RobAStBase/R/AllGeneric.R
===================================================================
--- branches/robast-1.2/pkg/RobAStBase/R/AllGeneric.R 2018-08-07 23:53:41 UTC (rev 1109)
+++ branches/robast-1.2/pkg/RobAStBase/R/AllGeneric.R 2018-08-08 21:49:35 UTC (rev 1110)
@@ -43,6 +43,9 @@
if(!isGeneric("evalIC")){
setGeneric("evalIC", function(IC, x) standardGeneric("evalIC"))
}
+if(!isGeneric("evalIC.v")){
+ setGeneric("evalIC.v", function(IC, x) standardGeneric("evalIC.v"))
+}
if(!isGeneric("makeIC")){
setGeneric("makeIC", function(IC, L2Fam, ...) standardGeneric("makeIC"))
}
Modified: branches/robast-1.2/pkg/RobAStBase/R/ContIC.R
===================================================================
--- branches/robast-1.2/pkg/RobAStBase/R/ContIC.R 2018-08-07 23:53:41 UTC (rev 1109)
+++ branches/robast-1.2/pkg/RobAStBase/R/ContIC.R 2018-08-08 21:49:35 UTC (rev 1110)
@@ -5,7 +5,7 @@
Risks, Infos, clip = Inf, cent = 0, stand = as.matrix(1),
lowerCase = NULL, neighborRadius = 0, w = new("HampelWeight"),
normtype = NormType(), biastype = symmetricBias(),
- modifyIC = NULL){
+ modifyIC = NULL, .fastFct = NULL){
if(missing(name))
name <- "IC of contamination type"
if(missing(Risks))
@@ -42,6 +42,7 @@
contIC at biastype <- biastype
contIC at normtype <- normtype
contIC at modifyIC <- modifyIC
+ contIC at .fastFct <- .fastFct
return(contIC)
# return(new("ContIC", name = name, Curve = Curve, Risks = Risks, Infos = Infos,
@@ -66,6 +67,7 @@
name = "IC of contamination type",
CallL2Fam = L2call,
Curve = generateIC.fct(neighbor, L2Fam, res),
+ .fastFct = generateIC.fast.fct(neighbor, L2Fam, res),
clip = b,
cent = a,
stand = A,
@@ -170,3 +172,5 @@
addInfo(object) <- c("CallL2Fam<-", "The entries in 'Risks' and 'Infos' may be wrong")
object
})
+
+
Modified: branches/robast-1.2/pkg/RobAStBase/R/HampIC.R
===================================================================
--- branches/robast-1.2/pkg/RobAStBase/R/HampIC.R 2018-08-07 23:53:41 UTC (rev 1109)
+++ branches/robast-1.2/pkg/RobAStBase/R/HampIC.R 2018-08-08 21:49:35 UTC (rev 1110)
@@ -18,3 +18,22 @@
object
})
+## evaluate IC
+setMethod("evalIC.v", signature(IC = ".fastIC", x = "numeric"),
+ function(IC, x){
+ if(is.null(IC at .fastFct)){
+ res <- setMethod("evalIC.v", signature(IC = "IC", x = "numeric"))(IC,x)
+ ## cast to matrix ICdim x nobs
+ }else{
+ res <- IC at .fastFct(x)
+ }
+ })
+setMethod("evalIC.v", signature(IC = ".fastIC", x = "matrix"),
+ function(IC, x){
+ if(is.null(IC at .fastFct)){
+ res <- setMethod("evalIC.v", signature(IC = "IC", x = "matrix"))(IC,x)
+ ## cast to matrix ICdim x nobs
+ }else{
+ res <- IC at .fastFct(x)
+ }
+ })
Modified: branches/robast-1.2/pkg/RobAStBase/R/IC.R
===================================================================
--- branches/robast-1.2/pkg/RobAStBase/R/IC.R 2018-08-07 23:53:41 UTC (rev 1109)
+++ branches/robast-1.2/pkg/RobAStBase/R/IC.R 2018-08-08 21:49:35 UTC (rev 1110)
@@ -85,6 +85,8 @@
return(prec)
})
+
+
## evaluate IC
setMethod("evalIC", signature(IC = "IC", x = "numeric"),
function(IC, x){
@@ -113,7 +115,12 @@
else
return(evalRandVar(Curve, x)[,,1])
})
+## evaluate IC
+setMethod("evalIC.v", signature(IC = "IC", x = "numeric"),
+ function(IC, x) sapply(x, function(x) evalIC(IC,x))
+ )
+
## make some L2function a pIC at a model
setMethod("makeIC", signature(IC = "IC", L2Fam = "missing"),
function(IC){
@@ -122,8 +129,8 @@
})
## make some L2function a pIC at a model
-setMethod("makeIC", signature(IC = "IC", L2Fam = "L2ParamFamily"),
- function(IC, L2Fam){
+setMethod("makeIC", signature(IC = "IC", L2Fam = "L2ParamFamily"),
+ function(IC, L2Fam){
dims <- length(L2Fam at param)
if(dimension(IC at Curve) != dims)
@@ -142,40 +149,27 @@
E10 <- E(L2Fam, IC1 %*% t(L2deriv))
E1 <- matrix(E10, dims, dims)
- stand <- trafo %*% solve(E1)
+ stand <- trafo %*% solve(E1)
Y <- as(stand %*% IC1, "EuclRandVariable")
- #ICfct <- vector(mode = "list", length = dims)
- #ICfct[[1]] <- function(x){Y(x)}
-
if(!is.function(IC at modifyIC))
IC at modifyIC <- function(L2Fam, IC, withMakeIC = FALSE) return(makeIC(IC,L2Fam))
-# modifyIC <- ..modifnew
-# }else{
-# .modifyIC <- IC at modifyIC
-# if(!is.null(attr(IC at modifyIC,"hasMakeICin.modifyIC"))){
-# modifyIC <- .modifyIC
-# }else{
-# modifyIC <- function(L2Fam, IC){ IC. <- .modifyIC(L2Fam, IC)
-# return(makeIC(IC., L2Fam)) }
-# }
-# }
-# }
-# attr(modifyIC,"hasMakeICin.modifyIC") <- TRUE
CallL2Fam <- L2Fam at fam.call
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"))),
+ 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 = IC at modifyIC))
})
+
+
# alias to IC needed here:
.IC <- IC
Modified: branches/robast-1.2/pkg/RobAStBase/R/TotalVarIC.R
===================================================================
--- branches/robast-1.2/pkg/RobAStBase/R/TotalVarIC.R 2018-08-07 23:53:41 UTC (rev 1109)
+++ branches/robast-1.2/pkg/RobAStBase/R/TotalVarIC.R 2018-08-08 21:49:35 UTC (rev 1110)
@@ -4,7 +4,7 @@
Risks, Infos, clipLo = -Inf, clipUp = Inf, stand = as.matrix(1),
lowerCase = NULL, neighborRadius = 0, w = new("BdStWeight"),
normtype = NormType(), biastype = symmetricBias(),
- modifyIC = NULL){
+ modifyIC = NULL, .fastFct = NULL){
if(missing(name))
name <- "IC of total variation type"
@@ -37,6 +37,7 @@
IC1 at biastype <- biastype
IC1 at normtype <- normtype
IC1 at modifyIC <- modifyIC
+ IC1 at .fastFct <- .fastFct
return(IC1)
}
@@ -65,6 +66,7 @@
name = "IC of total variation type",
CallL2Fam = L2call,
Curve = generateIC.fct(neighbor, L2Fam, res),
+ .fastFct = generateIC.fast.fct(neighbor, L2Fam, res),
clipUp = clipUp,
clipLo = clipLo,
stand = A,
Added: branches/robast-1.2/pkg/RobAStBase/R/combinedICs.R
===================================================================
--- branches/robast-1.2/pkg/RobAStBase/R/combinedICs.R (rev 0)
+++ branches/robast-1.2/pkg/RobAStBase/R/combinedICs.R 2018-08-08 21:49:35 UTC (rev 1110)
@@ -0,0 +1,53 @@
+combineOrthPICs <- function(pIC1, pIC2, combinedName = "combined IC", dim){
+ ## adds to complementary pICs to give one IC
+ ## the orthogonality is not checked here
+
+ IC <- new(".fastIC")
+ IC at name <- combinedName
+ pICC1 <- as(diag(dim)%*%pIC1 at Curve,"EuclRandVariable")
+ pICC2 <- as(diag(dim)%*%pIC2 at Curve,"EuclRandVariable")
+ IC at Curve <- EuclRandVarList(pICC1+pICC2)
+ IC at Risks <- pIC1 at Risks
+ if(length(pIC2 at Risks)) addRisk(IC) <- pIC2 at Risks
+ IC at Infos <- pIC1 at Infos
+ if(nrow(pIC2 at Infos)) addInfo(IC) <- pIC2 at Infos
+ IC at CallL2Fam <- pIC1 at CallL2Fam
+ .modifyIC.0 <- function(L2Fam, IC, withMakeIC = FALSE){
+ pic1 <- pic1 at modifyIC(L2Fam, pIC1, withMakeIC)
+ pic2 <- pic2 at modifyIC(L2Fam, pIC2, withMakeIC)
+ IC1 <- combineOrthPICs(pic1, pic2,combinedName)
+ return(IC1)
+ }
+ .modifyIC.1 <- function(L2Fam, IC, withMakeIC = FALSE){
+ IC1 <- .modifyIC.0(L2Fam, IC, withMakeIC)
+ IC1 at modifyIC <- .modifyIC.1
+ return(IC1)
+ }
+
+ IC at modifyIC <- .modifyIC.1
+ IC at .fastFct <- function(x){pIC1 at .fastFct(x)+pIC2 at .fastFct(x)}
+ return(IC)
+}
+
+
+.fastIC <- function(name ="", Curve = EuclRandVarList(RealRandVariable(Map = list(function(x){x}),
+ Domain = Reals())), Risks, Infos, CallL2Fam = call("L2ParamFamily"),
+ modifyIC = NULL, .fastFct = NULL){
+fastIC <- new(".fastIC")
+if(missing(Infos)) Infos <- fastIC at Infos
+if(missing(Risks)) Risks <- fastIC at Risks
+IC.0 <- IC(name, Curve, Risks, Infos, CallL2Fam, modifyIC)
+slotNms <- slotNames(class(IC.0))
+for(sN in slotNms) slot(fastIC, sN) <- slot(IC.0,sN)
+if(is.null(.fastFct)||missing(.fastFct)){
+ ICM <- IC.0 at Curve[[1]]@Map
+ .fastFct <- function(x){
+ if(is.null(dim(x)))
+ sapply(x, function(u) sapply(ICM, function(s)s(u)))
+ else
+ apply(x, 1,function(u) sapply(ICM, function(s)s(u)))
+ }
+}
+fastIC at .fastFct <- .fastFct
+return(fastIC)
+}
Modified: branches/robast-1.2/pkg/RobAStBase/R/ddPlot_utils.R
===================================================================
--- branches/robast-1.2/pkg/RobAStBase/R/ddPlot_utils.R 2018-08-07 23:53:41 UTC (rev 1109)
+++ branches/robast-1.2/pkg/RobAStBase/R/ddPlot_utils.R 2018-08-08 21:49:35 UTC (rev 1110)
@@ -266,9 +266,15 @@
ndata.y0[!isna] <- jitter(ndata.y0[!isna], factor=jitter.pts[2])
pdots$col <- col
+ inax <- is.na(ndata.x)
+ inay <- is.na(ndata.y)
+
+ nonina <- !inax&!inay
+
retV <- list(id.x=id0.x, id.y= id0.y, id.xy = id0.xy,
- qtx = quantile(ndata.x), qty = quantile(ndata.y),
- cutoff.x.v = co.x, cutoff.y.v = co.y)
+ qtx = quantile(ndata.x[nonina]),
+ qty = quantile(ndata.y[nonina]),
+ cutoff.x.v = co.x, cutoff.y.v = co.y)
if(doplot){
plotInfo<- list("plotArgs"=NULL)
Modified: branches/robast-1.2/pkg/RobAStBase/R/generateICfct.R
===================================================================
--- branches/robast-1.2/pkg/RobAStBase/R/generateICfct.R 2018-08-07 23:53:41 UTC (rev 1109)
+++ branches/robast-1.2/pkg/RobAStBase/R/generateICfct.R 2018-08-08 21:49:35 UTC (rev 1110)
@@ -5,53 +5,101 @@
A <- as.matrix(res$A)
a <- if(is(neighbor,"TotalVarNeighborhood")) 0 else res$a
b <- res$b
- d <- res$d
+ d <- if(!is.null(res$d)) res$d else 0
w <- weight(res$w)
nrvalues <- nrow(A)
dim <- ncol(A)
ICfct <- vector(mode = "list", length = nrvalues)
- Y <- as(A %*% L2Fam at L2deriv - a, "EuclRandVariable")
L <- as(diag(dim)%*%L2Fam at L2deriv, "EuclRandVariable")
+ distr <- distribution(L2Fam)
L.fct <- function(x) evalRandVar(L,x)
if(nrvalues == 1){
- if(!is.null(d)){
+ if(!is.null(res$d)){
ICfct[[1]] <- function(x){}
if(all(dim(trafo(L2Fam at param)) == c(1, 1))){
body(ICfct[[1]]) <- substitute(
- { ind <- 1-.eq(Y(x))
- Y(x)*w(L(x)) + zi*(1-ind)*d*b },
- list(Y = Y at Map[[1]], L = L.fct, w = w, b = b, d = d,
- zi = sign(trafo(L2Fam at param)), .eq = .eq))
+ { indS <- liesInSupport(Di,x,checkFin=TRUE)
+ Lx <- L(x)
+ Yx <- A %*% Lx - a
+ ind <- 1-.eq(Yx)
+ (Yx*w(Lx) + zi*(1-ind)*d*b)*indS },
+ list(L = L.fct, w = w, b = b, d = d, A = A, a = a,
+ zi = sign(trafo(L2Fam at param)), .eq = .eq, Di = distr))
}else{
body(ICfct[[1]]) <- substitute(
- { ind <- 1-.eq(Y(x))
- ifelse(ind, Y(x)*w(L(x)), NA) },
- list(Y = Y at Map[[1]], L = L.fct, w = w, b = b, d = d,
- .eq = .eq))
+ { indS <- liesInSupport(Di,x,checkFin=TRUE)
+ Lx <- L(x)
+ Yx <- A %*% Lx - a
+ ind <- 1-.eq(Yx)
+ ifelse(ind, Yx*w(Lx), NA)*indS },
+ list(L = L.fct, w = w, b = b, d = d, A = A, a = a,
+ .eq = .eq, Di = distr))
}
}else{
ICfct[[1]] <- function(x){}
- body(ICfct[[1]]) <- substitute({ Y(x)*w(L(x)) },
- list(Y = Y at Map[[1]], L = L.fct, w = w))
+ body(ICfct[[1]]) <- substitute({ indS <- liesInSupport(Di,x,checkFin=TRUE)
+ Lx <- L(x)
+ Yx <- A %*% Lx - a
+ Yx*w(Lx)*indS },
+ list(L = L.fct, A = A, a = a, w = w, Di = distr))
}
}else{
- if(!is.null(d))
+ if(!is.null(res$d))
for(i in 1:nrvalues){
ICfct[[i]] <- function(x){}
- body(ICfct[[i]]) <- substitute({ind <- 1-.eq(Yi(x))
- ind*Yi(x)*w(L(x)) + (1-ind)*d
+ body(ICfct[[i]]) <- substitute({indS <- liesInSupport(Di,x,checkFin=TRUE)
+ Lx <- L(x)
+ Yix <- Ai %*% Lx - ai
+ ind <- 1-.eq(Yix)
+ (ind*Yix*w(Lx) + (1-ind)*di)*indS
},
- list(Yi = Y at Map[[i]], L = L.fct, w = w,
- b = b, d = d[i]))#, .eq = .eq))
+ list(L = L.fct, Ai = A[i,,drop=FALSE], ai = a[i], w = w,
+ di = d[i], Di = distr))#, .eq = .eq))
}
else
for(i in 1:nrvalues){
ICfct[[i]] <- function(x){}
- body(ICfct[[i]]) <- substitute({ Yi(x)*w(L(x)) },
- list(Yi = Y at Map[[i]], L = L.fct, w = w))
+ body(ICfct[[i]]) <- substitute({indS <- liesInSupport(Di,x,checkFin=TRUE)
+ Lx <- L(x)
+ Yix <- Ai %*% Lx - ai
+ Yix*w(Lx)*indS },
+ list(L = L.fct, Ai = A[i,,drop=FALSE], ai = a[i], w = w, Di = distr))
}
}
- return(EuclRandVarList(EuclRandVariable(Map = ICfct, Domain = Y at Domain,
- Range = Y at Range)))
+ return(EuclRandVarList(EuclRandVariable(Map = ICfct, Domain = L at Domain,
+ Range = Reals()))) # EuclideanSpace(dimension = nrvalues))))
})
+## generate fast IC fct
+## for internal use only!
+generateIC.fast.fct <- function(neighbor, L2Fam, res){
+ A <- as.matrix(res$A)
+ a <- if(is(neighbor,"TotalVarNeighborhood")) 0 else res$a
+ b <- res$b
+ d <- res$d
+ w <- weight(res$w)
+ nrvalues <- nrow(A)
+ dims <- ncol(A)
+ L <- as(diag(dims)%*%L2Fam at L2deriv, "EuclRandVariable")
+ distr <- distribution(L2Fam)
+ L.fct <- function(x) evalRandVar(L,x)
+ fastFct <- function(x){}
+ if(nrvalues==1L){
+ d0 <- if(dims==1L) d else NA
+ }else{
+ d0 <- if(!is.null(d)) d else 0
+ }
+ zi0 <- if(nrvalues==1L && dims==1L) sign(trafo(L2Fam at param)) else 1
+ b0 <- if(nrvalues==1L) b else 1
+ body(fastFct) <- substitute({ indS <- liesInSupport(Di,x,checkFin=TRUE)
+ Lx <- L(x)
+ Yx <- A %*% Lx - a
+ ind <- 1-.eq(Yx)
+ ifelse(ind,Yx*w(Lx), zi*d*b)*indS
+ },
+ list(L = L.fct, w = w, b = b0,
+ d = d0 , A = A, a = a, zi = zi0,
+ .eq = .eq, Di = distr))
+ return(fastFct)
+ }
+
Modified: branches/robast-1.2/pkg/RobAStBase/R/kStepEstimator.R
===================================================================
--- branches/robast-1.2/pkg/RobAStBase/R/kStepEstimator.R 2018-08-07 23:53:41 UTC (rev 1109)
+++ branches/robast-1.2/pkg/RobAStBase/R/kStepEstimator.R 2018-08-08 21:49:35 UTC (rev 1110)
@@ -41,7 +41,8 @@
es.call[[1]] <- as.name("kStepEstimator")
## get some dimensions
- L2Fam <- eval(CallL2Fam(IC))
+ CallL2Fam <- CallL2Fam(IC)
+ L2Fam <- eval(CallL2Fam)
Param <- param(L2Fam)
tf <- trafo(L2Fam,Param)
@@ -53,6 +54,8 @@
p <- nrow(Dtau)
k <- ncol(Dtau)
+ CallL2FamK <- CallL2Fam
+ if(p!=k) CallL2FamK$trafo <- diag(k)
lmx <- length(main(L2Fam))
lnx <- length(nuisance(L2Fam))
@@ -80,9 +83,9 @@
### transform if necessary
x0 <- x
- x0 <- if(is.numeric(x) && ! is.matrix(x)) {
- x0 <- as.matrix(x)
- }
+ #x0 <- if(is.numeric(x) && ! is.matrix(x)) {
+ # x0 <- as.matrix(x)
+ # }
completecases <- complete.cases(x0)
if(na.rm) x0 <- na.omit(x0)
@@ -91,7 +94,7 @@
### use dispatch here (dispatch only on start)
#a.var <- if( is(start, "Estimate")) asvar(start) else NULL
- IC.UpdateInKer.0 <- if(is(start,"ALEstimate")) start at pIC else NULL
+ IC.UpdateInKer.0 <- if(is(start,"ALEstimate")) pIC(start) else NULL
force(startArgList)
start.val <- kStepEstimator.start(start, x=x0, nrvalues = k,
na.rm = na.rm, L2Fam = L2Fam,
@@ -122,10 +125,11 @@
ICList <- if(withICList) vector("list", steps) else NULL
cvar.fct <- function(L2, IC, dim, dimn =NULL){
+ IC.C <- as(diag(dim)%*%IC at Curve, "EuclRandVariable")
if(is.null(dimn)){
- return(matrix(E(L2, IC %*% t(IC)),dim,dim))
+ return(matrix(E(L2, IC.C %*% t(IC.C)),dim,dim))
}else{
- return(matrix(E(L2, IC %*% t(IC)),dim,dim, dimnames = dimn))
+ return(matrix(E(L2, IC.C %*% t(IC.C)),dim,dim, dimnames = dimn))
}
}
@@ -137,57 +141,58 @@
if(withPreModif){
main(Param)[] <- .deleteDim(u.theta[idx])
-# print(Param)
if (lnx) nuisance(Param)[] <- .deleteDim(u.theta[nuis.idx])
-# print(Param)
-# print(L2Fam)
L2Fam <- modifyModel(L2Fam, Param,
.withL2derivDistr = L2Fam at .withEvalL2derivDistr)
-# print(L2Fam)
IC <- modifyIC(IC)(L2Fam, IC, withMakeIC = FALSE)
- if(steps==1L && withMakeIC){
- IC <- makeIC(IC, L2Fam)
-# IC at modifyIC <- oldmodifIC
- }
- # print(IC)
+ CallL2Fam <- IC at CallL2Fam
+ if(steps==1L && withMakeIC) IC <- makeIC(IC, L2Fam)
}
- IC.c <- as(diag(p) %*% IC at Curve, "EuclRandVariable")
+ IC.c <- .fastIC(Curve=EuclRandVarList(as(diag(p) %*% IC at Curve, "EuclRandVariable")), CallL2Fam = CallL2Fam)
-# print(theta)
tf <- trafo(L2Fam, Param)
Dtau <- tf$mat
IC.tot.0 <- NULL
-# print(Dtau)
if(!.isUnitMatrix(Dtau)){
- # print("HU1!")
Dminus <- solve(Dtau, generalized = TRUE)
projker <- diag(k) - Dminus %*% Dtau
- IC.tot1 <- Dminus %*% IC.c
- IC.tot2 <- 0 * IC.tot1
+ IC.tot1 <- .fastIC(Curve=EuclRandVarList(as(Dminus %*% IC.c at Curve, "EuclRandVariable")), CallL2Fam = CallL2FamK)
+ IC.tot2.isnull <- TRUE
if(sum(diag(projker))>0.5 && ### is EM-D^-D != 0 (i.e. rk D<p)
withUpdateInKer){
if(!is.null(IC.UpdateInKer)&&!is(IC.UpdateInKer,"IC"))
warning("'IC.UpdateInKer' is not of class 'IC'; we use default instead.")
- IC.tot2 <- if(is.null(IC.UpdateInKer))
- getBoundedIC(L2Fam, D = projker) else
- as(projker %*% IC.UpdateInKer at Curve,
- "EuclRandVariable")
- IC.tot.0 <- IC.tot1 + IC.tot2
- }else{ if(!is.null(IC.UpdateInKer.0)){
- IC.tot.0 <- NULL
+ if(is.null(IC.UpdateInKer)){
+ IC.tot2 <- .fastIC(Curve=EuclRandVarList(getBoundedIC(L2Fam, D = projker)),
+ CallL2Fam = CallL2FamK)
}else{
- if(is.call(IC.UpdateInKer.0))
- IC.UpdateInKer.0 <- eval(IC.UpdateInKer.0)
- IC.tot.0 <- IC.tot1 + as(projker %*%
- IC.UpdateInKer.0 at Curve,
- "EuclRandVariable")
+ IC.tot2 <- .fastIC(Curve=EuclRandVarList(as(projker %*% IC.UpdateInKer at Curve,
+ "EuclRandVariable")),
+ CallL2Fam = CallL2FamK)
}
+ IC.tot2.isnull <- FALSE
+ IC.tot.0 <- combineOrthPICs(IC.tot1,IC.tot2,dim=k)
+ }else{if(is.null(IC.UpdateInKer.0)){
+ IC.tot.0 <- NULL
+ }else{
+ if(is.call(IC.UpdateInKer.0))
+ IC.UpdateInKer.0 <- eval(IC.UpdateInKer.0)
+ IC.tot.00 <- .fastIC(Curve= EuclRandVarList(as(projker %*% IC.UpdateInKer.0 at Curve,
+ "EuclRandVariable")),
+ CallL2Fam = CallL2FamK)
+ IC.tot.0 <- combineOrthPICs(IC.tot1,IC.tot.00,dim=k)
+ }
}
- IC.tot <- IC.tot1 + IC.tot2
- correct <- rowMeans(evalRandVar(IC.tot, x0), na.rm = na.rm)
+
+ IC.tot <- IC.tot1
+ correct <- rowMeans(evalIC.v(IC.tot1, x0), na.rm = na.rm)
+ if(!IC.tot2.isnull){
+ correct <- correct + rowMeans(evalIC.v(IC.tot2, x0), na.rm = na.rm)
+ IC.tot <- IC.tot.0
+ }
iM <- is.matrix(u.theta)
names(correct) <- if(iM) rownames(u.theta) else names(u.theta)
if(logtrf){
@@ -199,12 +204,8 @@
theta <- (tf$fct(u.theta[idx]))$fval
}else{
-# print("HU2!")
- correct <- rowMeans(evalRandVar(IC.c, x0), na.rm = na.rm )
+ correct <- rowMeans(evalIC.v(IC.c, x0), na.rm = na.rm )
iM <- is.matrix(theta)
-# print(sclname)
-# print(names(theta))
-# print(str(theta))
names(correct) <- if(iM) rownames(theta) else names(theta)
if(logtrf){
scl <- if(iM) theta[sclname,1] else theta[sclname]
@@ -217,7 +218,6 @@
IC.tot <- IC.c
u.theta <- theta
}
-# print("HU3!")
var0 <- u.var <- NULL
if(with.u.var){
@@ -228,8 +228,6 @@
L2F0 = L2Fam, IC0 = IC.tot.0, dim0 = k,
dimn0 = list(cnms,cnms)))
if(withEvalAsVar) u.var <- eval(u.var)
- # matrix(E(L2Fam, IC.tot.0 %*% t(IC.tot.0)),
- # k,k, dimnames = list(cnms,cnms))
}
if(!var.to.be.c){
var0 <- substitute(do.call(cfct, args = list(L2F0, IC0,
@@ -241,12 +239,9 @@
if(withPostModif){
main(Param)[] <- .deleteDim(u.theta[idx])
if (lnx) nuisance(Param)[] <- .deleteDim(u.theta[nuis.idx])
-# print(L2Fam)
L2Fam <- modifyModel(L2Fam, Param,
.withL2derivDistr = L2Fam at .withEvalL2derivDistr)
-# print(L2Fam)
IC <- modifyIC(IC)(L2Fam, IC, withMakeIC = withMakeIC)
-# print(IC)
}
return(list(IC = IC, Param = Param, L2Fam = L2Fam,
@@ -260,23 +255,17 @@
colnames(Infos) <- c("method", "message")
if(is(L2Fam, "L2GroupParamFamily")) useLast <- TRUE
- ### iteration
-# print(IC at Risks$asCov)
-# print(Risks(IC)$asCov)
-
ksteps <- matrix(0,ncol=steps, nrow = p)
uksteps <- matrix(0,ncol=steps, nrow = k)
rownames(ksteps) <- est.names
rownames(uksteps) <- u.est.names
if(!is(modifyIC(IC), "NULL") ){
for(i in 1:steps){
-# modif.old <- modifyIC(IC)
if(i>1){
IC <- upd$IC
L2Fam <- upd$L2Fam
if((i==steps)&&withMakeIC) IC <- makeIC(IC,L2Fam)
-# IC at modifyIC <- modif.old
Param <- upd$Param
tf <- trafo(L2Fam, Param)
@@ -285,17 +274,11 @@
upd <- updateStep(u.theta,theta,IC, L2Fam, Param,
withPreModif = withPre,
withPostModif = (steps>i) | useLast,
- with.u.var = (i==steps), oldmodifIC = modif.old)
+ with.u.var = (i==steps))
uksteps[,i] <- u.theta <- upd$u.theta
-# print(str(upd$theta))
-# print(nrow(ksteps))
ksteps[,i] <- theta <- upd$theta
if(withICList)
- ICList[[i]] <- new("InfluenceCurve",
- name = paste(gettext("(total) IC in step"),i),
- Risks = list(),
- Infos = matrix(c("",""),ncol=2),
- Curve = EuclRandVarList(upd$IC.tot))
+ ICList[[i]] <- upd$IC.tot
if(withPICList)
pICList[[i]] <- upd$IC.c
u.var <- upd$u.var
@@ -336,9 +319,6 @@
"computation of IC, asvar and asbias via useLast = FALSE"))
}
- ## if non-trivial trafo: info on how update was done
-# print(IC at Risks$asCov)
-# print(Risks(IC)$asCov)
if(! .isUnitMatrix(trafo(L2Fam)))
Infos <- rbind(Infos, c("kStepEstimator",
@@ -347,7 +327,6 @@
"modification in ker(trafo)")))
## some risks
-# print(list(u.theta=u.theta,theta=theta,u.var=u.var,var=var0))
if(var.to.be.c){
if("asCov" %in% names(Risks(IC)))
if(is.matrix(Risks(IC)$asCov) || length(Risks(IC)$asCov) == 1)
@@ -391,9 +370,11 @@
dimnames(asVar) <- list(nms.theta.idx, nms.theta.idx)
}
+ samplesize <- if(is.null(dim(x0))) length(x0) else nrow(x0)
+
estres <- new("kStepEstimate", estimate.call = es.call,
name = paste(steps, "-step estimate", sep = ""),
- estimate = theta, samplesize = nrow(x0), asvar = asVar,
+ estimate = theta, samplesize = samplesize, asvar = asVar,
trafo = tf, fixed = fixed, nuis.idx = nuis.idx,
untransformed.estimate = u.theta, completecases = completecases,
untransformed.asvar = u.var, asbias = asBias, pIC = IC,
Modified: branches/robast-1.2/pkg/RobAStBase/inst/NEWS
===================================================================
--- branches/robast-1.2/pkg/RobAStBase/inst/NEWS 2018-08-07 23:53:41 UTC (rev 1109)
+++ branches/robast-1.2/pkg/RobAStBase/inst/NEWS 2018-08-08 21:49:35 UTC (rev 1110)
@@ -15,7 +15,11 @@
+ slot function modifyIC of the different IC classes gains
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/robast -r 1110
More information about the Robast-commits
mailing list