[Robast-commits] r497 - in branches/robast-0.9/pkg: ROptEst ROptEst/R ROptEst/man ROptReg/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Jul 1 23:39:06 CEST 2012
Author: ruckdeschel
Date: 2012-07-01 23:39:06 +0200 (Sun, 01 Jul 2012)
New Revision: 497
Added:
branches/robast-0.9/pkg/ROptEst/R/getStartIC.R
branches/robast-0.9/pkg/ROptEst/R/internal-getpsi.R
branches/robast-0.9/pkg/ROptEst/R/internal.roptest.R
branches/robast-0.9/pkg/ROptEst/R/roptest.new.R
Removed:
branches/robast-0.9/pkg/ROptReg/R/ContIC.R
branches/robast-0.9/pkg/ROptReg/R/FixRobRegrTypeModel.R
branches/robast-0.9/pkg/ROptReg/R/TotalVarIC.R
branches/robast-0.9/pkg/ROptReg/R/getFiRiskRegTS.R
branches/robast-0.9/pkg/ROptReg/R/getFixClipRegTS.R
branches/robast-0.9/pkg/ROptReg/R/getFixRobRegTypeIC_fiUnOvShoot.R
Modified:
branches/robast-0.9/pkg/ROptEst/DESCRIPTION
branches/robast-0.9/pkg/ROptEst/NAMESPACE
branches/robast-0.9/pkg/ROptEst/R/AllGeneric.R
branches/robast-0.9/pkg/ROptEst/R/getReq.R
branches/robast-0.9/pkg/ROptEst/man/getReq.Rd
Log:
somehow forgot to commit these changes
Modified: branches/robast-0.9/pkg/ROptEst/DESCRIPTION
===================================================================
--- branches/robast-0.9/pkg/ROptEst/DESCRIPTION 2012-07-01 13:54:28 UTC (rev 496)
+++ branches/robast-0.9/pkg/ROptEst/DESCRIPTION 2012-07-01 21:39:06 UTC (rev 497)
@@ -4,7 +4,7 @@
Title: Optimally robust estimation
Description: Optimally robust estimation in general smoothly parameterized models using S4
classes and methods.
-Depends: R(>= 2.7.0), methods, distr(>= 2.0), distrEx(>= 2.0), distrMod(>= 2.0), RandVar(>=
+Depends: R(>= 2.10.0), methods, distr(>= 2.0), distrEx(>= 2.0), distrMod(>= 2.0), RandVar(>=
0.6.4), RobAStBase
Suggests: RobLox, MASS
Author: Matthias Kohl, Peter Ruckdeschel
Modified: branches/robast-0.9/pkg/ROptEst/NAMESPACE
===================================================================
--- branches/robast-0.9/pkg/ROptEst/NAMESPACE 2012-07-01 13:54:28 UTC (rev 496)
+++ branches/robast-0.9/pkg/ROptEst/NAMESPACE 2012-07-01 21:39:06 UTC (rev 497)
@@ -28,7 +28,7 @@
"getModifyIC",
"cniperCont", "cniperPoint", "cniperPointPlot")
exportMethods("updateNorm", "scaleUpdateIC", "eff",
- "get.asGRisk.fct")
+ "get.asGRisk.fct", "getStartIC")
export("getL2normL2deriv",
"asAnscombe", "asL1", "asL4",
"getReq", "getMaxIneff")
Modified: branches/robast-0.9/pkg/ROptEst/R/AllGeneric.R
===================================================================
--- branches/robast-0.9/pkg/ROptEst/R/AllGeneric.R 2012-07-01 13:54:28 UTC (rev 496)
+++ branches/robast-0.9/pkg/ROptEst/R/AllGeneric.R 2012-07-01 21:39:06 UTC (rev 497)
@@ -92,3 +92,6 @@
if(!isGeneric("get.asGRisk.fct")){
setGeneric("get.asGRisk.fct", function(Risk) standardGeneric("get.asGRisk.fct"))
}
+if(!isGeneric("getStartIC")){
+ setGeneric("getStartIC", function(model, risk, ...) standardGeneric("getStartIC"))
+}
Modified: branches/robast-0.9/pkg/ROptEst/R/getReq.R
===================================================================
--- branches/robast-0.9/pkg/ROptEst/R/getReq.R 2012-07-01 13:54:28 UTC (rev 496)
+++ branches/robast-0.9/pkg/ROptEst/R/getReq.R 2012-07-01 21:39:06 UTC (rev 497)
@@ -2,8 +2,9 @@
list(s = getRiskIC(IC,risk=trAsCov())$trAsCov$value^.5,
b = getRiskIC(IC,risk=asBias(),neighbor=neighbor)$asBias$value)
-getReq <- function(Risk,neighbor,IC1,IC2,n=1,upper=15){
- if(!is(IC1,"IC")||!is(IC2,"IC"))
+getReq <- function(Risk,neighbor,IC1,IC2,n=1,upper=15, radOrOutl=c("radius","Outlier")){
+ radOrOutl <- match.arg(radOrOutl)
+ if(!is(IC1,"IC")||!is(IC2,"IC"))
stop("Arguments IC1, IC2 must be of class 'IC'.")
if(!identical(IC1 at CallL2Fam,IC2 at CallL2Fam))
stop("Arguments IC1, IC2 must be of defined for the same model.")
@@ -24,6 +25,8 @@
get.asGRisk.fct(Risk)(r,s=sb2$s,b=sb2$b)
}
r0 <- uniroot(dRisk,lower=0, upper=upper)$root/n^.5
+ if(radOrOutl=="Outlier")
+ r0 <- if(sb1$s<=sb2$s) floor(r0*n) else ceiling(r0*n)
if(sb1$s<=sb2$s)
return(c(0,r0))
else
Added: branches/robast-0.9/pkg/ROptEst/R/getStartIC.R
===================================================================
--- branches/robast-0.9/pkg/ROptEst/R/getStartIC.R (rev 0)
+++ branches/robast-0.9/pkg/ROptEst/R/getStartIC.R 2012-07-01 21:39:06 UTC (rev 497)
@@ -0,0 +1,89 @@
+setMethod("getStartIC",signature(model = "ANY", risk = "ANY"),
+ function(model, risk, ...) stop("not yet implemented"))
+
+setMethod("getStartIC",signature(model = "L2ParamFamily", risk = "asRisk"),
+ function(model, risk, ..., ..debug=FALSE){
+ mc <- match.call(expand.dots=FALSE)
+ dots <- mc$"..."
+ fsCor <- dots[["fsCor"]]
+ eps <- dots[["eps"]]
+ dots[["eps"]] <- NULL
+ dots[["fsCor"]] <- NULL
+ neighbor <- dots$neighbor
+ dots$neighbor <- NULL
+
+ sm.rmx <- selectMethod("radiusMinimaxIC", signature(
+ class(model),class(neighbor),class(risk)))
+ dots.rmx <- .fix.in.defaults(dots, sm.rmx)
+ dots.rmx$L2Fam <- NULL
+ dots.rmx$neighbor <- NULL
+ dots.rmx$risk <- NULL
+
+ infMod <- InfRobModel(center = model, neighbor = neighbor)
+ sm.optic <- selectMethod("optIC", signature(
+ class(infMod),class(risk)))
+ dots.optic <- .fix.in.defaults(dots, sm.optic)
+ dots.optic$model <- NULL
+ dots.optic$risk <- NULL
+ if(is.null(eps$e)){
+ dots.rmx$loRad <- eps$sqn * eps$lower
+ dots.rmx$upRad <- eps$sqn * eps$upper
+ arg.rmx <- c(list(L2Fam = model, neighbor = neighbor,
+ risk = risk), dots.rmx)
+ if(..debug) print(c(arg.rmx=arg.rmx))
+ ICstart <- do.call(radiusMinimaxIC, args=arg.rmx)
+ if(..debug) print(ICstart)
+ if(!isTRUE(all.equal(fsCor, 1, tol = 1e-3))){
+ neighbor at radius <- neighborRadius(ICstart)*fsCor
+ arg.optic <- c(list( model = infMod, risk = risk), dots.optic)
+ if(..debug) print(c(arg.optic=arg.optic))
+ ICstart <- do.call(optIC, args = arg.optic)
+ }
+ }else{
+ neighbor at radius <- eps$sqn * fsCor * eps$e
+# print(neighbor)
+ infMod <- InfRobModel(center = model, neighbor = neighbor)
+ arg.optic <- c(list(model = infMod, risk = risk),
+ dots.optic)
+ if(..debug) print(c(arg.optic=arg.optic))
+# print(arg.optic)
+# print("----------------------------------------------------")
+ ICstart <- do.call(optIC, args = arg.optic)
+ }
+ return(ICstart)
+ })
+
+
+
+setMethod("getStartIC",signature(model = "L2ScaleShapeUnion", risk = "interpolRisk"),
+ function(model, risk, ...){
+
+ mc <- match.call(expand.dots=TRUE)
+
+ gridn <- type(risk)
+ nam <- name(model)
+ xi <- main(param(model))["shape"] #[scaleshapename(model)["shape"]]
+ beta <- main(param(model))["scale"] #[scaleshapename(model)["scale"]]
+ nsng <- character(0)
+ sng <- try(getFromNamespace(gridn, ns = "RobExtremes"),silent=TRUE)
+ if(!is(sng,"try-error")) nsng <- names(sng)
+ if(length(nsng)){
+ if(nam %in% nsng){
+ interpolfct <- sng[[nam]]$fct
+ #print(xi)
+ #print(beta)
+ .modifyIC <- function(L2Fam, IC){
+ para <- param(L2Fam)
+ xi0 <- main(para)["shape"]#[scaleshapename(L2Fam)["scale"]]
+ beta0 <- main(para)["scale"]#[scaleshapename(L2Fam)["scale"]]
+ .getPsi(xi0,beta0, interpolfct, L2Fam, type(risk))}
+ IC0 <- .getPsi(xi, beta, interpolfct, model, type(risk))
+ IC0 at modifyIC <- .modifyIC
+ return(IC0)
+ }
+ }
+ mc$risk <- if(type(risk)==".MBRE") asMSE() else asBias()
+ mc$neighbor <- ContNeighborhood(radius=0.5)
+ return(do.call(getStartIC, as.list(mc[-1])))
+ })
+
Added: branches/robast-0.9/pkg/ROptEst/R/internal-getpsi.R
===================================================================
--- branches/robast-0.9/pkg/ROptEst/R/internal-getpsi.R (rev 0)
+++ branches/robast-0.9/pkg/ROptEst/R/internal-getpsi.R 2012-07-01 21:39:06 UTC (rev 497)
@@ -0,0 +1,62 @@
+.getPsi <- function(xi, beta, fct, L2Fam , type){
+
+ L2deriv <- L2Fam at L2deriv
+ dbeta <- diag(c(beta,1))
+ b <- fct[[1]](xi)
+ a <- c(dbeta%*%c(fct[[2]](xi),fct[[3]](xi)))
+ aw <- c(dbeta%*%c(fct[[4]](xi),fct[[5]](xi)))
+ am <- mean(c(fct[[7]](xi),fct[[8]](xi)))
+ A <- dbeta%*%matrix(c(fct[[6]](xi),am,am,fct[[9]](xi)),2,2)%*%dbeta
+ am <- mean(c(fct[[11]](xi),fct[[12]](xi)))
+ Aw <- dbeta%*%matrix(c(fct[[10]](xi),am,am,fct[[13]](xi)),2,2)%*%dbeta
+
+
+
+ normt <- NormType()
+ biast <- symmetricBias()
+ nb <- ContNeighborhood(radius=0.5)
+ ICT <- paste("optimally robust IC for", switch(type,
+ ".OMSE"="maxMSE",".RMXE"="RMX", ".MBRE"="maxBias"))
+ riskT <- if(type!=".MBRE") "asGRisk" else "asBias"
+
+ w <- new("HampelWeight")
+ stand(w) <- Aw
+ cent(w) <- aw
+ clip(w) <- b
+
+ if(type!=".MBRE"){
+ weight(w) <- getweight(w, neighbor = nb, biastype = biast,
+ normW = normt)
+ }else weight(w) <- minbiasweight(w, neighbor = nb, biastype = biast,
+ normW = normt)
+
+ res <- list(a = a, A = A, b = b, d = 0*a,
+ normtype = normt, biastype = biast, w = w,
+ info = c("optIC", ICT), risk = list(),
+ modifyIC = NULL)
+
+
+ IC <- generateIC(nb, L2Fam, res)
+ return(IC)
+}
+
+if(FALSE){
+ res <- list(a = a, A = A, b = b, d = 0*a, w = w)
+
+ IC <- ContIC(name = "interpolated IC of contamination type",
+ CallL2Fam = L2Fam at fam.call,
+ Curve = generateIC.fct(nb, L2Fam, res),
+ clip = b,
+ cent = a,
+ stand = A,
+ lowerCase =0*a,
+ w = w,
+ neighborRadius = nb at radius,
+ modifyIC = NULL,
+ normtype = normt,
+ biastype = biast,
+ Risks = list(),
+ Infos = matrix(c("optIC", ICT), ncol = 2,
+ dimnames = list(character(0), c("method", "message"))
+ ))
+}
\ No newline at end of file
Copied: branches/robast-0.9/pkg/ROptEst/R/internal.roptest.R (from rev 476, branches/robast-0.9/pkg/RobExtremesBuffer/internal.roptest.R)
===================================================================
--- branches/robast-0.9/pkg/ROptEst/R/internal.roptest.R (rev 0)
+++ branches/robast-0.9/pkg/ROptEst/R/internal.roptest.R 2012-07-01 21:39:06 UTC (rev 497)
@@ -0,0 +1,239 @@
+.constructArg.list <- function(fun,matchCall, onlyFormal=FALSE, debug =FALSE){
+ dprint <- function(...) if(debug) print(...)
+ form0 <- form <- formals(fun)
+ form0["..."] <- NULL
+ dprint("----------------------")
+ es.call.0 <- list(as.list(matchCall)[[1]])
+ if(! "..." %in% names(form) && onlyFormal){
+ match.i <- 2
+ while(match.i<= length(matchCall)){
+ nm0 <- names(matchCall)[match.i]
+ if(!nm0 %in% names(form0)){
+ es.call.0[[nm0]] <- if(is.call(matchCall[[match.i]])){
+ dprint("HU");eval(matchCall[[match.i]])
+ }else matchCall[[match.i]]
+ matchCall[[match.i]] <- NULL
+ }else match.i <- match.i + 1
+ }
+ }else{es.call.0 <- matchCall}
+ dprint("----------------------")
+ dprint(matchCall)
+ dprint("----------------------")
+ for(form.i in seq(along=form0)) {
+ if(!is.null(names(form0)[form.i])){
+ dprint(form.i)
+ dprint(names(form0)[form.i])
+ dprint(form0[[form.i]])
+ dprint("----------------------")
+ nam0 <- names(form0)[form.i]
+ if(!nam0 %in% names(matchCall) &&
+ !is.symbol(form0[[form.i]])
+ ){
+ dprint("Try assign")
+ dprint(paste(nam0,": ",sep=""))
+ dprint(matchCall[[nam0]])
+ if(!is.null(form0[[form.i]])){
+ dprint(dum <- str(form0[[form.i]]))
+ xu <- if(is.call(form0[[form.i]])){
+ dprint("HU1"); eval(form0[[form.i]])
+ }else form0[[form.i]]
+ if(!is.null(xu)) {dprint(xu)
+ matchCall[[nam0]] <- xu
+ dprint(matchCall[[nam0]])}
+ }else{
+ matchCall[nam0] <- NULL
+ dprint("NULL")
+ }
+ dprint("----------------------")
+
+ }
+ }
+ }
+ mc <- list(mc=matchCall)
+ if(length(es.call.0)) mc$esc <- as.call(es.call.0)
+ dprint(mc)
+ return(mc)
+}
+
+.fix.in.defaults <- function(call.list, fun, withEval=TRUE){
+ formals.fun <- formals(fun)
+ k <- length(call.list)
+ L <- length(formals.fun)
+ if("..." %in% names(formals.fun)) L <- L-1
+# cat("\nFORMALS\n");print(formals.fun); cat("\nCALL\n");print(call.list)
+ for(i in 1:L){
+ foP <- paste(formals.fun[[i]])
+ if(length(formals.fun[[i]]))
+ if(length(foP)||foP!=""){
+ if(!names(formals.fun)[i] %in% names(call.list)&&!is.null(formals.fun[[i]])){
+ k <- k + 1
+ call.list[[k]] <- formals.fun[[i]]
+ if(withEval && is.call(formals.fun[[i]])){
+ xu <- eval(formals.fun[[i]])
+ if(!is.null(xu)) call.list[[k]] <- xu
+ }
+ if(length(foP)){
+ if(withEval && is.name(formals.fun[[i]])){
+ xu <- if(foP!="") get(foP) else NULL
+ if(!is.null(xu)) call.list[[k]] <- xu
+ }
+ }
+ names(call.list)[k] <- names(formals.fun)[i]
+ }
+ }
+ }
+# cat("\nNEWLIST\n"); print(call.list)
+ return(call.list)
+
+}
+
+.pretreat <- function(x, na.rm = TRUE){
+ if(missing(x))
+ stop("'x' is missing with no default")
+ if(!is.numeric(x)){
+ if(is.data.frame(x))
+ x <- data.matrix(x)
+ else
+ x <- as.matrix(x)
+ if(!is.matrix(x))
+ stop("'x' has to be a numeric vector resp. a matrix or data.frame")
+ }
+ completecases <- complete.cases(x)
+ if(na.rm) x <- na.omit(x)
+ return(list(x=x,completecases=completecases))
+}
+
+.check.eps <- function(...){
+ mc <- match.call(expand.dots=TRUE)
+
+ eps <- eps.lower <- eps.upper <- NULL
+# lapply(mc, function(xx) print(str(xx)))
+
+ ine <- is.null(mc[["eps"]]) || is.symbol(mc[["eps"]])
+ inl <- is.null(mc[["eps.lower"]]) || is.symbol(mc[["eps.lower"]])
+ inu <- is.null(mc[["eps.upper"]]) || is.symbol(mc[["eps.upper"]])
+
+ if(ine && inl && inu){
+ eps.lower <- 0
+ eps.upper <- 0.5
+ }
+ if(ine){
+ if(!inl && !inu)
+ eps.lower <- mc[["eps.lower"]]
+ eps.upper <- mc[["eps.upper"]]
+ if(!inl && inu)
+ eps.lower <- mc[["eps.lower"]]
+ eps.upper <- 0.5
+ if(inl && !inu)
+ eps.lower <- 0
+ eps.upper <- mc[["eps.upper"]]
+ if(length(eps.lower) != 1 || length(eps.upper) != 1)
+ stop("'eps.lower' and 'eps.upper' have to be of length 1")
+ if(!is.numeric(eps.lower) || !is.numeric(eps.upper) || eps.lower >= eps.upper)
+ stop("'eps.lower' < 'eps.upper' is not fulfilled")
+ if((!inl &&(eps.lower < 0)) || (!inu&& (eps.upper > 0.5)))
+ stop("'eps.lower' and 'eps.upper' have to be in [0, 0.5]")
+ }else{
+ eps <- mc[["eps"]]
+ if(length(eps) != 1)
+ stop("'eps' has to be of length 1")
+ if(!ine&& (eps == 0))
+ stop("'eps = 0'! => use functions 'mean' and 'sd' for estimation")
+ if(!ine && ((eps < 0) || (eps > 0.5)))
+ stop("'eps' has to be in (0, 0.5]")
+ }
+
+ x <- mc[["x"]]
+ if(is.matrix(x))
+ sqrtn <- sqrt(ncol(x))
+ else
+ sqrtn <- sqrt(length(x))
+ erg <- list(sqn = sqrtn)
+ erg[["e"]] <- eps
+ erg[["lower"]] <- eps.lower
+ erg[["upper"]] <- eps.upper
+ return(erg)
+}
+
+.isOKsteps <- function(steps){
+ if(!is.integer(steps))
+ steps <- as.integer(steps)
+ if(steps < 1){
+ stop("'steps' has to be some positive integer value")
+ }
+ if(length(steps) != 1){
+ stop("'steps' has to be of length 1")
+ }
+ return(invisible(NULL))
+}
+
+.isOKfsCor <- function(fsCor){
+ if(fsCor <= 0)
+ stop("'fsCor' has to be positive")
+ if(length(fsCor) != 1)
+ stop("'fsCor' has to be of length 1")
+ return(invisible(NULL))
+}
+
+
+#.getROptICstart <- function(..., ..eval=FALSE){
+# mc <- match.call(expand.dots=TRUE)
+# eps <- mc$eps
+# dots <- mc$dots
+#
+# if(is.null(eps$e)){
+# r.lower <- eps$sqn * eps$lower
+# r.upper <- eps$sqn * eps$upper
+# ICstart <- substitute(do.call(radiusMinimaxIC,
+# c(list(L2Fam = mc$L2FamStart, neighbor = mc$neighbor,
+# risk = mc$risk,
+# loRad = r.lower, upRad = r.upper,
+# verbose = mc$verbose,
+# OptOrIter = mc$OptOrIter),dots)))
+# if(..eval) ICstart <- eval(ICstart)
+#
+# if(!isTRUE(all.equal(mc$fsCor, 1, tol = 1e-3))){
+# neighbor at radius <- neighborRadius(ICstart)*mc$fsCor
+# infMod <- InfRobModel(center = mc$L2FamStart, neighbor = mc$neighbor)
+# ICstart <- substitute(do.call(optIC, c(list( model = mc$infMod, risk = mc$risk,
+# verbose = mc$verbose, OptOrIter = mc$OptOrIter),
+# dots)))
+# if(..eval) ICstart <- eval(ICstart)
+# }
+# }else{
+# neighbor at radius <- eps$sqn*eps$e*mc$fsCor
+# infMod <- InfRobModel(center = mc$L2FamStart, neighbor = mc$neighbor)
+# ICstart <- substitute(do.call(optIC, c(list(model = mc$infMod, risk = mc$risk,
+# verbose = mc$verbose, OptOrIter = mc$OptOrIter),
+# dots)))
+# if(..eval) ICstart <- eval(ICstart)
+#
+# }
+# return(ICstart)
+#}
+
+genkStepCtrl <- function(useLast = getRobAStBaseOption("kStepUseLast"),
+ withUpdateInKer = getRobAStBaseOption("withUpdateInKer"),
+ IC.UpdateInKer = getRobAStBaseOption("IC.UpdateInKer"),
+ withICList = getRobAStBaseOption("withICList"),
+ withPICList = getRobAStBaseOption("withPICList"),
+ scalename = "scale", withLogScale = TRUE){
+ es.call <- match.call()
+ es.list <- as.list(es.call[-1])
+ es.list <- .fix.in.defaults(es.list,genkStepCtrl)
+ return(es.list)
+}
+genstartCtrl<- function(initial.est = NULL, initial.est.ArgList = NULL,
+ startPar = NULL, distance = CvMDist){
+ es.call <- match.call()
+ es.list <- as.list(es.call[-1])
+ es.list <- .fix.in.defaults(es.list,genstartCtrl)
+ return(es.list)
+}
+gennbCtrl <- function(neighbor = ContNeighborhood(),
+ eps, eps.lower, eps.upper){
+ es.call <- match.call()
+ es.list <- as.list(es.call[-1])
+ es.list <- .fix.in.defaults(es.list,gennbCtrl)
+ return(es.list)
+}
Copied: branches/robast-0.9/pkg/ROptEst/R/roptest.new.R (from rev 476, branches/robast-0.9/pkg/RobExtremesBuffer/roptest.new.R)
===================================================================
--- branches/robast-0.9/pkg/ROptEst/R/roptest.new.R (rev 0)
+++ branches/robast-0.9/pkg/ROptEst/R/roptest.new.R 2012-07-01 21:39:06 UTC (rev 497)
@@ -0,0 +1,231 @@
+###############################################################################
+## Optimally robust estimation
+###############################################################################
+
+roptest.n <- function(x, L2Fam, eps, eps.lower, eps.upper, fsCor = 1, initial.est,
+ neighbor = ContNeighborhood(), risk = asMSE(), steps = 1L,
+ distance = CvMDist, startPar = NULL, verbose = NULL,
+ OptOrIter = "iterate",
+ useLast = getRobAStBaseOption("kStepUseLast"),
+ withUpdateInKer = getRobAStBaseOption("withUpdateInKer"),
+ IC.UpdateInKer = getRobAStBaseOption("IC.UpdateInKer"),
+ withICList = getRobAStBaseOption("withICList"),
+ withPICList = getRobAStBaseOption("withPICList"),
+ na.rm = TRUE, initial.est.ArgList, ...,
+ withLogScale = TRUE,..withCheck=FALSE){
+ es.call <- match.call()
+ es.call0 <- match.call(expand.dots=FALSE)
+ dots <- es.call0[["..."]]
+ es.call0$"..." <- NULL
+ es.call1 <- .constructArg.list(roptest.n,es.call0, onlyFormal=FALSE,
+ debug = ..withCheck)$mc
+
+ res <- .constructArg.list(gennbCtrl,es.call1, onlyFormal=TRUE,
+ debug = ..withCheck)
+ nbCtrl0 <- as.list(res$mc[-1])
+ es.call1 <- res$esc
+ res <- .constructArg.list(genstartCtrl,es.call1, onlyFormal=TRUE,
+ debug = ..withCheck)
+ startCtrl0 <- as.list(res$mc[-1])
+ es.call1 <- res$esc
+ res <- .constructArg.list(genkStepCtrl,es.call1, onlyFormal=TRUE,
+ debug = ..withCheck)
+ kStepCtrl0 <- as.list(res$mc[-1])
+ es.call1 <- res$esc
+
+ list1 <- as.list(es.call1[-1])
+ res <- do.call(gennbCtrl,nbCtrl0)
+ if(length(res)) list1$nbCtrl <- res
+ res <- do.call(genstartCtrl,startCtrl0)
+ if(length(res)) list1$startCtrl <- res
+ res <- do.call(genkStepCtrl,kStepCtrl0)
+ if(length(res)) list1$kStepCtrl <- res
+ list1 <- c(list1,dots)
+ list1$..withCheck <- NULL
+ if(..withCheck) print(list1)
+ if(..withCheck) return(substitute(do.call(roptEst, lis), list(lis=list1)))
+ else{
+ res <- do.call(roptEst, list1)
+ res at estimate.call <- es.call
+ return(res)}
+}
+#roptest.n(x=1:10,L2Fam=GammaFamily(),also=3,..withCheck=TRUE)
+
+#roptest.n(x=1:10,L2Fam=GammaFamily(),also=3)
+
+roptEst <- function(x, L2Fam, fsCor = 1,
+ risk = asMSE(), steps = 1L,
+ verbose = NULL,
+ OptOrIter = "iterate",
+ nbCtrl = gennbCtrl(),
+ startCtrl = genstartCtrl(),
+ kStepCtrl = genkStepCtrl(),
+ na.rm = TRUE, ..., debug = FALSE){
+
+ es.call <- match.call()
+ es.call0 <- match.call(expand.dots=FALSE)
+ dots <- es.call0[["..."]]
+ es.call0$"..." <- NULL
+# if(!is.null(dots[["escall"]])) es.call <- dots[["escall"]]
+
+ if(missing(nbCtrl)||is.null(nbCtrl)) nbCtrl <- gennbCtrl()
+ if(missing(startCtrl)||is.null(startCtrl)) startCtrl <- genstartCtrl()
+ if(missing(kStepCtrl)||is.null(kStepCtrl)) kStepCtrl <- genkStepCtrl()
+ nbCtrl <- .fix.in.defaults(nbCtrl, gennbCtrl)
+ startCtrl <- .fix.in.defaults(startCtrl, genstartCtrl)
+ kStepCtrl <- .fix.in.defaults(kStepCtrl, genkStepCtrl)
+
+ es.list <- as.list(es.call0[-1])
+ es.list <- c(es.list,nbCtrl)
+ es.list$nbCtrl <- NULL
+ es.list0 <- c(es.list,dots)
+
+ if(missing(verbose)|| is.null(verbose))
+ es.list$verbose <- getRobAStBaseOption("all.verbose")
+
+ if(missing(L2Fam))
+ stop("'L2Fam' is missing with no default")
+ if(!is(L2Fam, "L2ParamFamily"))
+ stop("'L2Fam' must be of class 'L2ParamFamily'.")
+
+
+ res.x <- .pretreat(x,na.rm)
+ x <- res.x$x
+ completecases <- res.x$completecases
+
+ es.list[["x"]] <- x
+ if(debug){ cat("\nes.list:\n");print(es.list);cat("\n\n")}
+
+ .isOKfsCor(fsCor)
+
+ .isOKsteps(steps)
+
+ if(debug){
+ if(is.null(startCtrl$initial.est)){
+ print(substitute(MDEstimator(x = x0, ParamFamily = L2Fam0,
+ distance = distance0,
+ startPar = startPar0,dots0),
+ list(x0=x,L2Fam0=L2Fam, distance0=startCtrl$distance,
+ startPar0=startCtrl$startPar0, dots0=dots)))
+ startCtrl$initial.est <- "BLUB"
+ }
+ }else{
+ if(is.null(startCtrl$initial.est))
+ startCtrl$initial.est <- MDEstimator(x = x, ParamFamily = L2Fam,
+ distance = startCtrl$distance,
+ startPar = startCtrl$startPar, ...)
+ }
+ nrvalues <- length(L2Fam at param)
+
+ if(debug){
+ if(!is.null(startCtrl$initial.est.ArgList)){
+ initial.est <- print(substitute({kStepEstimator.start(initial.est0, x = x0,
+ nrvalues = nrvalues0, na.rm = na.rm0,
+ L2Fam = L2Fam0,
+ startList = startCtrl0)},
+ list(x0=x,L2Fam0=L2Fam,
+ initial.est0=startCtrl$initial.est,
+ nrvalues0=nrvalues,na.rm0=na.rm,
+ startCtrl0 = startCtrl$initial.est.ArgList)
+ ))
+ }else{
+ print(substitute(kStepEstimator.start(initial.est0, x = x0,
+ nrvalues = nrvalues0, na.rm = na.rm0,
+ L2Fam = L2Fam0),list(x0=x,L2Fam0=L2Fam,
+ initial.est0=startCtrl$initial.est,
+ nrvalues0=nrvalues,na.rm0=na.rm)
+ ))
+ }
+ }else{
+ sy <- system.time({
+ initial.est <- kStepEstimator.start(eval(startCtrl$initial.est), x = x,
+ nrvalues = nrvalues, na.rm = na.rm,
+ L2Fam = L2Fam)
+ })
+ print(sy)
+ }
+
+ newParam <- param(L2Fam)
+
+ if(!debug){
+ main(newParam)[] <- as.numeric(initial.est)
+ L2FamStart <- modifyModel(L2Fam, newParam)
+ }
+ if(debug) print(risk)
+
+ if(!is(risk,"interpolRisk"))
+ es.list0[["eps"]] <- do.call(.check.eps, args=c(nbCtrl,list("x"=x)))
+ es.list0$risk <- NULL
+ es.list0$L2Fam <- NULL
+ neighbor <- eval(es.list0$neighbor)
+ es.list0$neighbor <- NULL
+ es.list0$kStepCtrl <- NULL
+ es.list0$startCtrl <- NULL
+ es.list0$distance <- NULL
+ es.list0$x <- NULL
+ es.list0$steps <- NULL
+ es.list0$eps.lower <- NULL
+ es.list0$eps.upper <- NULL
+ es.list0$na.rm <- NULL
+ if(debug) {cat("\n\n\n::::\n\n")
+ argList <- c(list(model=L2Fam,risk=risk,neighbor=neighbor),
+ es.list0)
+ print(argList)
+ cat("\n\n\n")
+ }
+ if(!debug){
+ sy <- system.time({
+ ICstart <- do.call(getStartIC, args=c(list(model=L2FamStart,risk=risk,
+ neighbor=neighbor),es.list0))
+ })
+ print(sy)
+ }
+ if(debug){
+ argList <- list(x, IC = ICstart, start = initial.est, steps = steps,
+ useLast = kStepCtrl$useLast,
+ withUpdateInKer = kStepCtrl$withUpdateInKer,
+ IC.UpdateInKer = kStepCtrl$IC.UpdateInKer,
+ withICList = kStepCtrl$withICList,
+ withPICList = kStepCtrl$withPICList,
+ na.rm = na.rm,
+ scalename = kStepCtrl$scalename,
+ withLogScale = kStepCtrl$withLogScale)
+ print(argList) }
+ sy <- system.time({
+ res <- kStepEstimator(x, IC = ICstart, start = initial.est, steps = steps,
+ useLast = kStepCtrl$useLast,
+ withUpdateInKer = kStepCtrl$withUpdateInKer,
+ IC.UpdateInKer = kStepCtrl$IC.UpdateInKer,
+ withICList = kStepCtrl$withICList,
+ withPICList = kStepCtrl$withPICList,
+ na.rm = na.rm,
+ scalename = kStepCtrl$scalename,
+ withLogScale = kStepCtrl$withLogScale)
+ })
+ print(sy)
+
+ if(!debug) res at estimate.call <- es.call
+ Infos <- matrix(c("roptest",
+ paste(steps, "-step estimate for ", name(L2Fam), sep = "")),
+ ncol = 2)
+ colnames(Infos) <- c("method", "message")
+
+
+ if(! distrMod:::.isUnitMatrix(trafo(L2Fam)))
+ Infos <- rbind(Infos, c("method"="roptest",
+ "message"=paste("computation of IC",
+ ifelse(kStepCtrl$withUpdateInKer,"with","without") ,
+ "modification in ker(trafo)")))
+
+ Infos <- rbind(Infos, c("method"="roptest",
+ "message"=paste("computation of IC, asvar and asbias via useLast =",
+ kStepCtrl$useLast)))
+
+ if(debug) print(Infos)
+ if(debug) return(NULL)
+
+ Infos(res) <- Infos
+ res at completecases <- completecases
+ res at start <- initial.est
+ return(res)
+}
Modified: branches/robast-0.9/pkg/ROptEst/man/getReq.Rd
===================================================================
--- branches/robast-0.9/pkg/ROptEst/man/getReq.Rd 2012-07-01 13:54:28 UTC (rev 496)
+++ branches/robast-0.9/pkg/ROptEst/man/getReq.Rd 2012-07-01 21:39:06 UTC (rev 497)
@@ -1,11 +1,13 @@
\name{getReq}
\alias{getReq}
-\title{getReq -- computation of the radius interval where IC1 is better than IC2}
+\title{getReq -- computation of the radius interval where IC1 is better than IC2.}
\description{
- (tries to) compute a radius interval where IC1 is better than IC2
+ (tries to) compute a radius interval where IC1 is better than IC2,
+ respectively the number of (worst-case) outliers interval where IC1 is
+ better than IC2.
}
-\usage{getReq(Risk,neighbor,IC1,IC2,n=1,upper=15)}
+\usage{getReq(Risk,neighbor,IC1,IC2,n=1,upper=15, radOrOutl=c("radius","Outlier"))}
\arguments{
\item{Risk}{an object of class \code{"asGRisk"} -- the risk at which IC1 is better than IC2.}
\item{neighbor}{ object of class \code{"Neighborhood"}; the neighborhood at which to compute the bias. }
@@ -15,6 +17,9 @@
in the shrinking neighborhood setting of Rieder[94]. Otherwise the radius interval is scaled
down accordingly.}
\item{upper}{the upper bound of the radius interval in which to search}
+ \item{radOrOutl}{ a character string specifying whether an interval of
+ radii or a number of outliers is returned; must be one
+ of "radius" (default) and "Outlier". }
}
%\details{}
\value{The radius interval (given by its endpoints) where \code{IC1} is better than \code{IC2}
@@ -57,6 +62,7 @@
d.qn <- (2^.5*qnorm(5/8))^-1
IC.qn <- makeIC(function(x) d.qn*(1/4 - pnorm(x+1/d.qn) + pnorm(x-1/d.qn)), L2M)
getReq(asMSE(), neighbor, IC.mad, IC.qn)
+getReq(asMSE(), neighbor, IC.mad, IC.qn, radOrOutl = "Outlier", n = 30)
# => MAD is better once r > 0.5144 (i.e. for more than 2 outliers for n = 30)
}
\concept{Hampel risk}
Deleted: branches/robast-0.9/pkg/ROptReg/R/ContIC.R
===================================================================
--- branches/robast-0.9/pkg/ROptReg/R/ContIC.R 2012-07-01 13:54:28 UTC (rev 496)
+++ branches/robast-0.9/pkg/ROptReg/R/ContIC.R 2012-07-01 21:39:06 UTC (rev 497)
@@ -1,72 +0,0 @@
-## generate IC
-## for internal use only!
-setMethod("generateIC", signature(neighbor = "ContNeighborhood",
- L2Fam = "L2RegTypeFamily"),
- function(neighbor, L2Fam, res){
- A <- res$A
- a <- res$a
- b <- res$b
- d <- res$d
- nrvalues <- nrow(A)
- ICfct <- vector(mode = "list", length = nrvalues)
- Y <- as(A %*% L2Fam at L2deriv - a, "EuclRandVariable")
- if(nrvalues == 1){
- if(!is.null(d)){
- ICfct[[1]] <- function(x){
- ind <- (Y(x) != 0)
- b*(ind*Y(x)/(ind*absY(x) + (1-ind)*1) + zi*(1-ind)*d)
- }
- body(ICfct[[1]]) <- substitute(
- { ind <- (Y(x) != 0)
- b*(ind*Y(x)/(ind*absY(x) + (1-ind)*1) + zi*(1-ind)*d) },
- list(Y = Y at Map[[1]], absY = abs(Y)@Map[[1]], b = b, d = d,
- zi = sign(L2Fam at param@trafo)))
- }else{
- ICfct[[1]] <- function(x){ Y(x)*pmin(1, b/absY(x)) }
- body(ICfct[[1]]) <- substitute({ Y(x)*pmin(1, b/absY(x)) },
- list(Y = Y at Map[[1]], absY = abs(Y)@Map[[1]], b = b))
- }
- }
- else{
- absY <- sqrt(Y %*% Y)
- if(!is.null(d))
- for(i in 1:nrvalues){
- ICfct[[i]] <- function(x){ ind <- (Yi(x) != 0) ; ind*b*Yi(x)/absY(x) + (1-ind)*d }
- body(ICfct[[i]]) <- substitute({ ind <- (Yi(x) != 0) ; ind*b*Yi(x)/absY(x) + (1-ind)*d },
- list(Yi = Y at Map[[i]], absY = absY at Map[[1]], b = b, d = d[i]))
- }
- else
- for(i in 1:nrvalues){
- ICfct[[i]] <- function(x){ Yi(x)*pmin(1, b/absY(x)) }
- body(ICfct[[i]]) <- substitute({ Yi(x)*pmin(1, b/absY(x)) },
- list(Yi = Y at Map[[i]], absY = absY at Map[[1]], b = b))
- }
- }
- return(ContIC(
- name = "IC of contamination type",
- CallL2Fam = call("L2RegTypeFamily",
- name = L2Fam at name,
- distribution = L2Fam at distribution,
- param = L2Fam at param,
- props = L2Fam at props,
- ErrorDistr = L2Fam at ErrorDistr,
- ErrorSymm = L2Fam at ErrorSymm,
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/robast -r 497
More information about the Robast-commits
mailing list