[Robast-commits] r1111 - in branches/robast-1.2/pkg/RobAStBase: . R inst inst/chkTimeCode man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Aug 10 01:43:02 CEST 2018
Author: ruckdeschel
Date: 2018-08-10 01:43:02 +0200 (Fri, 10 Aug 2018)
New Revision: 1111
Added:
branches/robast-1.2/pkg/RobAStBase/inst/chkTimeCode/
branches/robast-1.2/pkg/RobAStBase/inst/chkTimeCode/TimingChecks.R
Modified:
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/AllShow.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/bALEstimate.R
branches/robast-1.2/pkg/RobAStBase/R/combinedICs.R
branches/robast-1.2/pkg/RobAStBase/R/generateICfct.R
branches/robast-1.2/pkg/RobAStBase/R/getPIC.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/ALEstimate-class.Rd
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/internals.Rd
Log:
[RobAStBase] branch 1.1
we are finally there...
I reverted most changes from yesterday and now we have decent timings again...
The clue was that to force optimal ICs to respect the support of the model
distribution, during evaluation of kStepEstimator it is prohibitive to
put line liesInSupport in each of the coordinate functions as this blows
up the integration time for covariances;
instead, we use helper .fixInLiesInSupport in file generateICfct.R
which after computation of variances inserts this in the Maps of
the IC
+ for time checking use file TimingChecks.R (with the preparation that
the lines commented out by ##-t-## in kStepEstimator.R have to be activated;
this uses helper function .addTime to produce a matrix with detailed timing
information which can be read out as argument ) -- it is in package
system folder "chkTimeCode" (in inst/chkTimeCode in r-forge)
BTW: by preallocating memory for the matrix with the timings we do not
gain speed unfortunately, so commenting in and out seems to be the best
option
Modified: branches/robast-1.2/pkg/RobAStBase/NAMESPACE
===================================================================
--- branches/robast-1.2/pkg/RobAStBase/NAMESPACE 2018-08-08 21:49:35 UTC (rev 1110)
+++ branches/robast-1.2/pkg/RobAStBase/NAMESPACE 2018-08-09 23:43:02 UTC (rev 1111)
@@ -20,12 +20,13 @@
"FixRobModel",
"InfRobModel")
exportClasses("InfluenceCurve",
- "IC", "HampIC", ".fastIC",
+ "IC", "HampIC",
"ContIC",
"TotalVarIC")
-exportClasses("RobAStControl", "RobWeight", "BoundedWeight",
+exportClasses("RobAStControl", "RobWeight", "BoundedWeight",
"BdStWeight", "HampelWeight")
-exportClasses("ALEstimate", "MCALEstimate", "kStepEstimate", "MEstimate")
+exportClasses("ALEstimate", "MCALEstimate", "kStepEstimate", "MEstimate",
+ "CvMMD.ALEstimate", "ML.ALEstimate" )
exportClasses("cutoff")
exportClasses("interpolRisk", "OMSRRisk","MBRRisk","RMXRRisk")
exportClasses("StartClass", "pICList", "OptionalpICList", "OptionalCall",
@@ -44,7 +45,7 @@
"modifyIC",
"generateIC",
"checkIC",
- "evalIC", "evalIC.v",
+ "evalIC",
"clip", "clip<-",
"cent", "cent<-",
"stand", "stand<-",
@@ -88,3 +89,4 @@
export(".rescalefct",".plotRescaledAxis",".makedotsP",".makedotsLowLevel",".SelectOrderData")
export(".merge.lists")
export("InfoPlot", "ComparePlot", "PlotIC")
+export(".fixInLiesInSupport")
\ No newline at end of file
Modified: branches/robast-1.2/pkg/RobAStBase/R/AllClass.R
===================================================================
--- branches/robast-1.2/pkg/RobAStBase/R/AllClass.R 2018-08-08 21:49:35 UTC (rev 1110)
+++ branches/robast-1.2/pkg/RobAStBase/R/AllClass.R 2018-08-09 23:43:02 UTC (rev 1111)
@@ -105,6 +105,7 @@
stop("'Infos' must have two columns")
else TRUE
})
+## comment 20180809: reverted changes in rev 1110
## partial incluence curve
setClass("IC", representation(CallL2Fam = "call",
modifyIC = "OptionalFunction"),
@@ -127,11 +128,6 @@
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",
@@ -139,7 +135,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())),
@@ -153,7 +149,7 @@
neighborRadius = 0,
biastype = symmetricBias(),
NormType = NormType()),
- contains = ".fastIC",
+ contains = "IC",
validity = function(object){
if(any(object at neighborRadius < 0)) # radius vector?!
stop("'neighborRadius' has to be in [0, Inf]")
@@ -265,6 +261,8 @@
pIC = NULL),
contains = c("ALEstimate","MCEstimate")
)
+setClass("CvMMD.ALEstimate",contains = c("MCALEstimate","CvMMDEstimate"))
+setClass("ML.ALEstimate",contains = c("MCALEstimate","MLEstimate"))
setClass("kStepEstimate",
representation(steps = "integer",
Modified: branches/robast-1.2/pkg/RobAStBase/R/AllGeneric.R
===================================================================
--- branches/robast-1.2/pkg/RobAStBase/R/AllGeneric.R 2018-08-08 21:49:35 UTC (rev 1110)
+++ branches/robast-1.2/pkg/RobAStBase/R/AllGeneric.R 2018-08-09 23:43:02 UTC (rev 1111)
@@ -43,12 +43,10 @@
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"))
}
+## comment 20180809: reverted changes in rev 1110
if(!isGeneric("clip")){
setGeneric("clip", function(x1, ...) standardGeneric("clip"))
}
Modified: branches/robast-1.2/pkg/RobAStBase/R/AllShow.R
===================================================================
--- branches/robast-1.2/pkg/RobAStBase/R/AllShow.R 2018-08-08 21:49:35 UTC (rev 1110)
+++ branches/robast-1.2/pkg/RobAStBase/R/AllShow.R 2018-08-09 23:43:02 UTC (rev 1111)
@@ -105,6 +105,20 @@
show(pIC(object))
}
})
+setMethod("show", "MCALEstimate",
+ function(object){
+ digits <- getOption("digits")
+ getMethod("show","MCEstimate")(object)
+ if(getdistrModOption("show.details") != "minimal"){
+ cat("asymptotic bias:\n")
+ print(asbias(object), quote = FALSE)
+ }
+ if(getdistrModOption("show.details") == "maximal" && !is.null(pIC(object))){
+ cat("(partial) influence curve:\n")
+ show(pIC(object))
+ }
+ })
+
setMethod("show", "kStepEstimate",
function(object){
digits <- getOption("digits")
Modified: branches/robast-1.2/pkg/RobAStBase/R/ContIC.R
===================================================================
--- branches/robast-1.2/pkg/RobAStBase/R/ContIC.R 2018-08-08 21:49:35 UTC (rev 1110)
+++ branches/robast-1.2/pkg/RobAStBase/R/ContIC.R 2018-08-09 23:43:02 UTC (rev 1111)
@@ -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, .fastFct = NULL){
+ modifyIC = NULL){
if(missing(name))
name <- "IC of contamination type"
if(missing(Risks))
@@ -42,7 +42,6 @@
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,
@@ -67,7 +66,6 @@
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,
@@ -172,5 +170,4 @@
addInfo(object) <- c("CallL2Fam<-", "The entries in 'Risks' and 'Infos' may be wrong")
object
})
-
-
+## comment 20180809: reverted changes in rev 1110
\ No newline at end of file
Modified: branches/robast-1.2/pkg/RobAStBase/R/HampIC.R
===================================================================
--- branches/robast-1.2/pkg/RobAStBase/R/HampIC.R 2018-08-08 21:49:35 UTC (rev 1110)
+++ branches/robast-1.2/pkg/RobAStBase/R/HampIC.R 2018-08-09 23:43:02 UTC (rev 1111)
@@ -18,22 +18,4 @@
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)
- }
- })
+## comment 20180809: reverted changes in rev 1110
\ No newline at end of file
Modified: branches/robast-1.2/pkg/RobAStBase/R/IC.R
===================================================================
--- branches/robast-1.2/pkg/RobAStBase/R/IC.R 2018-08-08 21:49:35 UTC (rev 1110)
+++ branches/robast-1.2/pkg/RobAStBase/R/IC.R 2018-08-09 23:43:02 UTC (rev 1111)
@@ -85,8 +85,6 @@
return(prec)
})
-
-
## evaluate IC
setMethod("evalIC", signature(IC = "IC", x = "numeric"),
function(IC, x){
@@ -115,12 +113,7 @@
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){
@@ -129,8 +122,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)
@@ -149,27 +142,40 @@
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
@@ -214,3 +220,4 @@
if(forceIC) IC.0 <- makeIC(IC.0, L2Fam)
return(IC.0)
})
+## comment 20180809: reverted changes in rev 1110
\ No newline at end of file
Modified: branches/robast-1.2/pkg/RobAStBase/R/TotalVarIC.R
===================================================================
--- branches/robast-1.2/pkg/RobAStBase/R/TotalVarIC.R 2018-08-08 21:49:35 UTC (rev 1110)
+++ branches/robast-1.2/pkg/RobAStBase/R/TotalVarIC.R 2018-08-09 23:43:02 UTC (rev 1111)
@@ -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, .fastFct = NULL){
+ modifyIC = NULL){
if(missing(name))
name <- "IC of total variation type"
@@ -37,7 +37,6 @@
IC1 at biastype <- biastype
IC1 at normtype <- normtype
IC1 at modifyIC <- modifyIC
- IC1 at .fastFct <- .fastFct
return(IC1)
}
@@ -66,7 +65,6 @@
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,
@@ -172,3 +170,4 @@
addInfo(object) <- c("CallL2Fam<-", "The entries in 'Risks' and 'Infos' may be wrong")
object
})
+## comment 20180809: reverted changes in rev 1110
\ No newline at end of file
Modified: branches/robast-1.2/pkg/RobAStBase/R/bALEstimate.R
===================================================================
--- branches/robast-1.2/pkg/RobAStBase/R/bALEstimate.R 2018-08-08 21:49:35 UTC (rev 1110)
+++ branches/robast-1.2/pkg/RobAStBase/R/bALEstimate.R 2018-08-09 23:43:02 UTC (rev 1111)
@@ -5,20 +5,22 @@
setMethod("pIC", "ALEstimate", function(object){
pIC0 <- .getPIC(object)
- eval.parent(substitute(object at pIC <- pIC0))
+ if(is(pIC0,"IC")) eval.parent(substitute(object at pIC <- pIC0))
return(pIC0)
})
setMethod("pIC", "MCEstimate", function(object){
if("pIC" %in% slotNames(class(object))){
pIC0 <- .getPIC(object)
- eval.parent(substitute(object at pIC <- pIC0))
+ if(is(pIC0,"IC")) eval.parent(substitute(object at pIC <- pIC0))
return(pIC0)
}else{
return(getPIC(object))
}})
setMethod("pIC", "MCALEstimate", getMethod("pIC", "ALEstimate"))
+setMethod("pIC", "ML.ALEstimate", getMethod("pIC", "ALEstimate"))
+setMethod("pIC", "CvMMD.ALEstimate", getMethod("pIC", "ALEstimate"))
setMethod("asbias", "ALEstimate", function(object) object at asbias)
setMethod("steps", "kStepEstimate", function(object) object at steps)
Modified: branches/robast-1.2/pkg/RobAStBase/R/combinedICs.R
===================================================================
--- branches/robast-1.2/pkg/RobAStBase/R/combinedICs.R 2018-08-08 21:49:35 UTC (rev 1110)
+++ branches/robast-1.2/pkg/RobAStBase/R/combinedICs.R 2018-08-09 23:43:02 UTC (rev 1111)
@@ -1,3 +1,9 @@
+################################################################################
+if(FALSE){
+################################################################################
+## 20180809: reverted changes from rev 1110
+################################################################################
+
combineOrthPICs <- function(pIC1, pIC2, combinedName = "combined IC", dim){
## adds to complementary pICs to give one IC
## the orthogonality is not checked here
@@ -51,3 +57,7 @@
fastIC at .fastFct <- .fastFct
return(fastIC)
}
+################################################################################
+## end if(FALSE)
+################################################################################
+}
\ No newline at end of file
Modified: branches/robast-1.2/pkg/RobAStBase/R/generateICfct.R
===================================================================
--- branches/robast-1.2/pkg/RobAStBase/R/generateICfct.R 2018-08-08 21:49:35 UTC (rev 1110)
+++ branches/robast-1.2/pkg/RobAStBase/R/generateICfct.R 2018-08-09 23:43:02 UTC (rev 1111)
@@ -8,70 +8,68 @@
d <- if(!is.null(res$d)) res$d else 0
w <- weight(res$w)
nrvalues <- nrow(A)
- dim <- ncol(A)
+ dims <- ncol(A)
ICfct <- vector(mode = "list", length = nrvalues)
- L <- as(diag(dim)%*%L2Fam at L2deriv, "EuclRandVariable")
+ L <- as(diag(dims)%*%L2Fam at L2deriv, "EuclRandVariable")
distr <- distribution(L2Fam)
+
L.fct <- function(x) evalRandVar(L,x)
if(nrvalues == 1){
if(!is.null(res$d)){
ICfct[[1]] <- function(x){}
- if(all(dim(trafo(L2Fam at param)) == c(1, 1))){
+ if(dims==1L){
body(ICfct[[1]]) <- substitute(
- { indS <- liesInSupport(Di,x,checkFin=TRUE)
- Lx <- L(x)
+ { Lx <- L(x)
Yx <- A %*% Lx - a
ind <- 1-.eq(Yx)
- (Yx*w(Lx) + zi*(1-ind)*d*b)*indS },
+ (Yx*w(Lx) + zi*(1-ind)*d*b) },
list(L = L.fct, w = w, b = b, d = d, A = A, a = a,
- zi = sign(trafo(L2Fam at param)), .eq = .eq, Di = distr))
+ zi = sign(trafo(L2Fam at param)), .eq = .eq))
}else{
body(ICfct[[1]]) <- substitute(
- { indS <- liesInSupport(Di,x,checkFin=TRUE)
- Lx <- L(x)
+ { Lx <- L(x)
Yx <- A %*% Lx - a
ind <- 1-.eq(Yx)
- ifelse(ind, Yx*w(Lx), NA)*indS },
+ ifelse(ind, Yx*w(Lx), NA) },
list(L = L.fct, w = w, b = b, d = d, A = A, a = a,
- .eq = .eq, Di = distr))
+ .eq = .eq))
}
}else{
ICfct[[1]] <- function(x){}
- body(ICfct[[1]]) <- substitute({ indS <- liesInSupport(Di,x,checkFin=TRUE)
- Lx <- L(x)
+ body(ICfct[[1]]) <- substitute({ Lx <- L(x)
Yx <- A %*% Lx - a
- Yx*w(Lx)*indS },
- list(L = L.fct, A = A, a = a, w = w, Di = distr))
+ Yx*w(Lx) },
+ list(L = L.fct, A = A, a = a, w = w))
}
}else{
if(!is.null(res$d))
for(i in 1:nrvalues){
ICfct[[i]] <- function(x){}
- body(ICfct[[i]]) <- substitute({indS <- liesInSupport(Di,x,checkFin=TRUE)
- Lx <- L(x)
+ body(ICfct[[i]]) <- substitute({Lx <- L(x)
Yix <- Ai %*% Lx - ai
ind <- 1-.eq(Yix)
- (ind*Yix*w(Lx) + (1-ind)*di)*indS
+ (ind*Yix*w(Lx) + (1-ind)*di)
},
list(L = L.fct, Ai = A[i,,drop=FALSE], ai = a[i], w = w,
- di = d[i], Di = distr))#, .eq = .eq))
+ di = d[i]))#, .eq = .eq))
}
else
for(i in 1:nrvalues){
ICfct[[i]] <- function(x){}
- body(ICfct[[i]]) <- substitute({indS <- liesInSupport(Di,x,checkFin=TRUE)
- Lx <- L(x)
+ body(ICfct[[i]]) <- substitute({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))
+ Yix*w(Lx) },
+ list(L = L.fct, Ai = A[i,,drop=FALSE], ai = a[i], w = w))
}
}
return(EuclRandVarList(EuclRandVariable(Map = ICfct, Domain = L at Domain,
Range = Reals()))) # EuclideanSpace(dimension = nrvalues))))
})
+## comment 20180809: reverted changes in rev 1110 as to generate.fast.fc:
## generate fast IC fct
## for internal use only!
+if(FALSE){
generateIC.fast.fct <- function(neighbor, L2Fam, res){
A <- as.matrix(res$A)
a <- if(is(neighbor,"TotalVarNeighborhood")) 0 else res$a
@@ -102,4 +100,13 @@
.eq = .eq, Di = distr))
return(fastFct)
}
+}
+.fixInLiesInSupport<- function(IC, distr){
+ MapL <- IC at Curve[[1]]@Map
+ for(i in 1:length(MapL))
+ body(IC at Curve[[1]]@Map[[i]]) <- substitute({
+ liesInSupport(distr,x,checkFin=TRUE)*fct(x)
+ }, list(fct = MapL[[i]], distr=distr))
+ return(IC)
+}
Modified: branches/robast-1.2/pkg/RobAStBase/R/getPIC.R
===================================================================
--- branches/robast-1.2/pkg/RobAStBase/R/getPIC.R 2018-08-08 21:49:35 UTC (rev 1110)
+++ branches/robast-1.2/pkg/RobAStBase/R/getPIC.R 2018-08-09 23:43:02 UTC (rev 1111)
@@ -5,13 +5,21 @@
setMethod(".checkEstClassForParamFamily",
signature=signature(PFam="ANY",estimator="MCEstimate"),
- function(PFam, estimator){
+ function(PFam, estimator) .extendbyPIC(PFam, estimator, "MCALEstimate"))
+setMethod(".checkEstClassForParamFamily",
+ signature=signature(PFam="ANY",estimator="MLEstimate"),
+ function(PFam, estimator) .extendbyPIC(PFam, estimator, "ML.ALEstimate"))
+setMethod(".checkEstClassForParamFamily",
+ signature=signature(PFam="ANY",estimator="CvMMDEstimate"),
+ function(PFam, estimator) .extendbyPIC(PFam, estimator, "CvMMD.ALEstimate"))
+
+.extendbyPIC <- function(PFam, estimator, toClass){
fromSlotNames <- slotNames(class(estimator))
- to <- new("MCALEstimate")
+ to <- new(toClass)
for(item in fromSlotNames) slot(to, item) <- slot(estimator,item)
to at pIC <- substitute(getPIC(estimator0), list(estimator0=estimator))
to
- } )
+ }
.getPIC <- function(object){
if(is.null(object at pIC)) return(NULL)
Modified: branches/robast-1.2/pkg/RobAStBase/R/kStepEstimator.R
===================================================================
--- branches/robast-1.2/pkg/RobAStBase/R/kStepEstimator.R 2018-08-08 21:49:35 UTC (rev 1110)
+++ branches/robast-1.2/pkg/RobAStBase/R/kStepEstimator.R 2018-08-09 23:43:02 UTC (rev 1111)
@@ -24,6 +24,18 @@
setMethod("neighborRadius","ANY",function(object)NA)
+.addTime <- function(timold,timnew,namenew){
+ nameold <- rownames(timold)
+ tim <- rbind(timold,timnew)
+ rownames(tim) <- c(nameold,namenew)
+ return(tim)
+}
+
+.ensureDim2 <- function(x){
+ d <- dim(x)
+ if(length(d)==3L && d[3]==1L) dim(x) <- d[1:2]
+ x }
+
### no dispatch on top layer -> keep product structure of dependence
kStepEstimator <- function(x, IC, start = NULL, steps = 1L,
useLast = getRobAStBaseOption("kStepUseLast"),
@@ -41,8 +53,12 @@
es.call[[1]] <- as.name("kStepEstimator")
## get some dimensions
- CallL2Fam <- CallL2Fam(IC)
- L2Fam <- eval(CallL2Fam)
+##-t-## syt <- system.time({
+ L2Fam <- eval(CallL2Fam(IC))
+##-t-## })
+##-t-## sytm <- matrix(syt,nrow=1)
+##-t-## rownames(sytm) <- "eval(CallL2Fam(IC))"
+##-t-## colnames(sytm) <- names(syt)
Param <- param(L2Fam)
tf <- trafo(L2Fam,Param)
@@ -54,8 +70,6 @@
p <- nrow(Dtau)
k <- ncol(Dtau)
- CallL2FamK <- CallL2Fam
- if(p!=k) CallL2FamK$trafo <- diag(k)
lmx <- length(main(L2Fam))
lnx <- length(nuisance(L2Fam))
@@ -83,9 +97,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)
@@ -94,11 +108,20 @@
### use dispatch here (dispatch only on start)
#a.var <- if( is(start, "Estimate")) asvar(start) else NULL
+##-t-## syt <- system.time({
IC.UpdateInKer.0 <- if(is(start,"ALEstimate")) pIC(start) else NULL
+##-t-## })
+##-t-## sytm <- .addTime(sytm,syt,"pIC(start)")
+ ## pIC(start) instead of start at pIC to potentially eval a call
+
force(startArgList)
+
+##-t-## syt <- system.time({
start.val <- kStepEstimator.start(start, x=x0, nrvalues = k,
na.rm = na.rm, L2Fam = L2Fam,
startList = startArgList)
+##-t-## })
+##-t-## sytm <- .addTime(sytm,syt,"kStepEstimator.start")
### use Logtransform here in scale models
sclname <- ""
@@ -106,6 +129,7 @@
logtrf <- is(L2Fam, "L2ScaleUnion") &
withLogScale & sclname %in% names(start.val)
### a starting value in k-space
+# print(start.val)
u.theta <- start.val
theta <- if(is(start.val,"Estimate")) estimate(start.val)
else trafoF(u.theta[idx])$fval
@@ -125,40 +149,65 @@
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.C %*% t(IC.C)),dim,dim))
+ return(matrix(E(L2, IC %*% t(IC)),dim,dim))
}else{
- return(matrix(E(L2, IC.C %*% t(IC.C)),dim,dim, dimnames = dimn))
+ return(matrix(E(L2, IC %*% t(IC)),dim,dim, dimnames = dimn))
}
}
+##-t-## updStp <- 0
### update - function
updateStep <- function(u.theta, theta, IC, L2Fam, Param,
withPreModif = FALSE,
- withPostModif = TRUE, with.u.var = FALSE
+ withPostModif = TRUE, with.u.var = FALSE,
+ withEvalAsVar.0 = FALSE
){
+##-t-## updStp <<- updStp + 1
if(withPreModif){
main(Param)[] <- .deleteDim(u.theta[idx])
+# print(Param)
if (lnx) nuisance(Param)[] <- .deleteDim(u.theta[nuis.idx])
+# print(Param)
+# print(L2Fam)
+##-t-## syt <- system.time({
L2Fam <- modifyModel(L2Fam, Param,
.withL2derivDistr = L2Fam at .withEvalL2derivDistr)
+##-t-## })
+##-t-## sytm <<- .addTime(sytm,syt,paste("modifyModel-PreModif-",updStp))
+# print(L2Fam)
+##-t-## syt <- system.time({
IC <- modifyIC(IC)(L2Fam, IC, withMakeIC = FALSE)
- CallL2Fam <- IC at CallL2Fam
- if(steps==1L && withMakeIC) IC <- makeIC(IC, L2Fam)
+##-t-## })
+##-t-## sytm <<- .addTime(sytm,syt,paste("modifyIC-PreModif-",updStp))
+ if(steps==1L && withMakeIC){
+##-t-## syt <- system.time({
+ IC <- makeIC(IC, L2Fam)
+##-t-## })
+##-t-## sytm <<- .addTime(sytm,syt,paste("modifyIC-makeIC-",updStp))
+# IC at modifyIC <- oldmodifIC
+ }
+ # print(IC)
}
- IC.c <- .fastIC(Curve=EuclRandVarList(as(diag(p) %*% IC at Curve, "EuclRandVariable")), CallL2Fam = CallL2Fam)
+##-t-## syt <- system.time({
+ IC.c <- as(diag(p) %*% IC at Curve, "EuclRandVariable")
+##-t-## })
+##-t-## sytm <<- .addTime(sytm,syt,paste("IC.c <- as(diag(p) %*%-",updStp))
+# 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 <- .fastIC(Curve=EuclRandVarList(as(Dminus %*% IC.c at Curve, "EuclRandVariable")), CallL2Fam = CallL2FamK)
+ IC.tot1 <- Dminus %*% IC.c
+# IC.tot2 <- 0 * IC.tot1
IC.tot2.isnull <- TRUE
if(sum(diag(projker))>0.5 && ### is EM-D^-D != 0 (i.e. rk D<p)
@@ -166,33 +215,41 @@
if(!is.null(IC.UpdateInKer)&&!is(IC.UpdateInKer,"IC"))
warning("'IC.UpdateInKer' is not of class 'IC'; we use default instead.")
if(is.null(IC.UpdateInKer)){
- IC.tot2 <- .fastIC(Curve=EuclRandVarList(getBoundedIC(L2Fam, D = projker)),
- CallL2Fam = CallL2FamK)
+##-t-## syt <- system.time({
+ IC.tot2 <- getBoundedIC(L2Fam, D = projker)
+##-t-## })
+##-t-## sytm <<- .addTime(sytm,syt,paste("getBoundedIC-",updStp))
}else{
- IC.tot2 <- .fastIC(Curve=EuclRandVarList(as(projker %*% IC.UpdateInKer at Curve,
- "EuclRandVariable")),
- CallL2Fam = CallL2FamK)
+##-t-## syt <- system.time({
+ IC.tot2 <- as(projker %*% IC.UpdateInKer at Curve, "EuclRandVariable")
+##-t-## })
+##-t-## sytm <<- .addTime(sytm,syt,paste("IC.tot2<-as(projker...-",updStp))
}
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.0 <- IC.tot1 + IC.tot2
+ }else{ if(is.null(IC.UpdateInKer.0)){
+ IC.tot.0 <- NULL
+ }else{
+##-t-## syt <- system.time({
+ if(is.call(IC.UpdateInKer.0))
+ IC.UpdateInKer.0 <- eval(IC.UpdateInKer.0)
+##-t-## })
+##-t-## sytm <<- .addTime(sytm,syt,paste("eval(IC.UpdateInKer.0)-",updStp))
+##-t-## syt <- system.time({
+ IC.tot.0 <- IC.tot1 + as(projker %*%
+ IC.UpdateInKer.0 at Curve,
+ "EuclRandVariable")
+##-t-## })
+##-t-## sytm <<- .addTime(sytm,syt,paste("IC.tot.0 <- IC.tot1 + as(proj-",updStp))
+ }
}
-
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
- }
+ if(!IC.tot2.isnull) IC.tot <- IC.tot1 + IC.tot2
+##-t-## syt <- system.time({
+ indS <- liesInSupport(distribution(L2Fam),x0,checkFin=TRUE)
+ correct <- rowMeans(t(t(.ensureDim2(evalRandVar(IC.tot, x0)))*indS), na.rm = na.rm)
+##-t-## })
+##-t-## sytm <<- .addTime(sytm,syt,paste("Dtau-not-Unit:correct <- rowMeans-",updStp))
iM <- is.matrix(u.theta)
names(correct) <- if(iM) rownames(u.theta) else names(u.theta)
if(logtrf){
@@ -204,8 +261,16 @@
theta <- (tf$fct(u.theta[idx]))$fval
}else{
- correct <- rowMeans(evalIC.v(IC.c, x0), na.rm = na.rm )
+# print("HU2!")
+##-t-## syt <- system.time({
+ indS <- liesInSupport(distribution(L2Fam),x0,checkFin=TRUE)
+ correct <- rowMeans(t(t(.ensureDim2(evalRandVar(IC.c, x0)))*indS), na.rm = na.rm)
+##-t-## })
+##-t-## sytm <<- .addTime(sytm,syt,paste("Dtau=Unit:correct <- rowMeans-",updStp))
iM <- is.matrix(theta)
+# print(sclname)
+# print(names(theta))
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/robast -r 1111
More information about the Robast-commits
mailing list