[Robast-commits] r476 - in branches/robast-0.9/pkg: . RobExtremesBuffer
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon May 21 17:49:11 CEST 2012
Author: ruckdeschel
Date: 2012-05-21 17:49:11 +0200 (Mon, 21 May 2012)
New Revision: 476
Added:
branches/robast-0.9/pkg/RobExtremesBuffer/
branches/robast-0.9/pkg/RobExtremesBuffer/FctWithCtrl.R
branches/robast-0.9/pkg/RobExtremesBuffer/README.txt
branches/robast-0.9/pkg/RobExtremesBuffer/getOptICGPD.R
branches/robast-0.9/pkg/RobExtremesBuffer/internal.roptest.R
branches/robast-0.9/pkg/RobExtremesBuffer/interpolGrids.R
branches/robast-0.9/pkg/RobExtremesBuffer/roptest.new.R
Log:
to share some ideas: some test files in connection with RobExtremes
Added: branches/robast-0.9/pkg/RobExtremesBuffer/FctWithCtrl.R
===================================================================
--- branches/robast-0.9/pkg/RobExtremesBuffer/FctWithCtrl.R (rev 0)
+++ branches/robast-0.9/pkg/RobExtremesBuffer/FctWithCtrl.R 2012-05-21 15:49:11 UTC (rev 476)
@@ -0,0 +1,229 @@
+###########################################################
+## function with Control as S4 class PR 21.5.12
+###########################################################
+
+
+## utilities:
+
+## we use a slightly different semantics from R here:
+## if formals are endowed with default values
+## and if they are not matched by exact or partial matching
+## they get set to default __before__ positional matching
+
+## .fix.in.defaults is to perform the default setting
+## (before positional matching)
+
+.fix.in.defaults <- function(call.list, fun){
+
+ formals.fun <- formals(fun)
+ k <- length(call.list)
+ L <- length(formals.fun)
+ if("..." %in% names(formals.fun)) L <- L-1
+ for(i in 1:L){
+ if(!is(formals.fun[[i]],"name")){
+ if(!names(formals.fun)[i] %in% names(call.list)&&
+ !is.null(formals.fun[[i]])){
+ k <- k + 1
+ call.list[[k]] <- formals.fun[[i]]
+ names(call.list)[k] <- names(formals.fun)[i]
+ }
+ }
+ }
+ return(call.list)
+
+}
+
+## as our Control argument is coming as a (NAMED!) list
+## and we want to be able to use these names, we attach the
+## list; at the end of the call (with on.exit(), so that this
+## is done even if an error was thrown), we call .clean.search
+## to remove no-longer needed positions in the search path.
+
+.clean.search<-function(what="structure\\(list\\("){
+
+ se <- search()
+ se0 <- grep(what,search(), value=TRUE)
+ du <- sapply(se0, detach, character.only=TRUE)
+ return(invisible(NULL))
+}
+
+## our S4 class FctWithCtrl: the function is _not_ a slot but rather
+## the class it self ( "is" instead of "has" )
+## and has additional slot Ctrl, a list
+## to be able to have access to Ctrl, one would need something
+## like a "self" which is not there in S4 classes
+## so instead we build something around it (see below)
+
+setClass("FctWithCtrl", representation(Ctrl="list"), contains="function")
+
+### generating function generates an object of class FctWithCtrl
+## (for simplicity of discussion assumed assigned to symbol "fu")
+## takes two arguments f (the function to be extended) and Ctrl, the
+## control list; this achieves
+## (a) [almost] arbitrary functions f can be used as "function" (to be enhanced
+## by Ctrl [restriction: no name of the Ctrl list may appear
+## in the formals nor in the actual arguments (passed through "...")
+## (b) an arbitrary list as Ctrl argument; no checking so far...
+## (c) a call to fu then behaves exactly as a call to argument f, except that
+## in the body of fu one has access to the items of Ctrl
+## -> no clever handling of masking so far
+
+## trick: we write to layers of wrapping functions which are:
+## + an anonymous top layer function doing the management of
+## call, names etc
+## + an inner function f1 which is a variation of f
+## with a modified arg list (the formals of f1 are the union of formals
+## of f and items of Ctrl)
+## with the same return value as f
+## with the body of f1 being the body of f with additional access
+## to the items of Ctrl
+### + the real trick is that Ctrl is an attribute to the generated
+## object and that we can access these arguments from an inner
+## function with sys.function(1)
+
+FctWithCtrl <- function(f, Ctrl){
+
+ g <- new("FctWithCtrl",
+ Ctrl = Ctrl,
+ function(...){
+
+ ## do exact & partial matching by names
+ mc <- match.call(expand=TRUE)
+
+ dots <- mc$...
+ mc$... <- NULL
+
+ Ctrl<- attributes(sys.function(1))$Ctrl
+ nmsall <- unique(c(names(Ctrl),names(formals(f))))
+
+ ## remove items from dots -> needed?
+ if(length(dots)&&!is.null(names(dots))){
+ for(i in 1: length(dots))
+ if(names(dots)[i] %in% nmsall) dots[[i]] <- NULL
+ }
+
+ ## remove items from Ctrl if as named arguments in mc
+ if(length(mc)&&!is.null(names(mc))){
+ for(i in 1: length(Ctrl))
+ if(names(Ctrl)[i] %in% names(mc)) Ctrl[[i]] <- NULL
+ }
+
+ ## combine to new call
+ mc <- as.call(c(as.list(mc),Ctrl,dots))
+
+ ## after exact and partial matching and before positional
+ ## matching set formals to defaults:
+ mc <- as.call(.fix.in.defaults(as.list(mc),f))
+
+ ## positional matching
+ mcfit <- as.list(mc[-1])
+ nmcfit <- names(mcfit)
+ unfit <- setdiff(names(formals(f)),nmcfit)
+ for(i in 1:length(nmcfit)){
+ if(nmcfit[i]=="" && length(unfit)>0){
+ names(mcfit)[i] <- if(unfit[1]!="...") unfit[1] else ""
+ unfit <- unfit[-1]
+ }
+ }
+
+ ## produce new call from the manipulated list
+ mc0 <- as.list(mc)
+ for(i in 1: length(mcfit))
+ {mc0[[i+1]] <- mcfit[[i]]
+ names(mc0)[i+1] <- names(mcfit)[i]
+ }
+ mc <- as.call(mc0)
+
+ ## generate variant f1 from f
+ ##
+ f1 <- f
+ ## combine formals
+ formals(f1) <- c(formals(f),Ctrl)
+ ## manipulate body:
+
+ body(f1) <- substitute({
+ on.exit(.clean.search())
+ attach(ctrl)
+ ..z <- body.f
+ return(..z)
+ }, list(lu=attributes(sys.function(1)),
+ ctrl=attributes(sys.function(1))$Ctrl,
+ body.f = body(f))
+ )
+ ## fall this f1 with the manipulated call
+ erg <- do.call(f1, as.list(mc[-1]))
+ return(erg)
+ })
+ return(g)
+}
+
+
+### EXAMPLE:
+
+## we want to produce a function f where all access to Ctrl arguments
+## is demonstrated: (attention
+##
+
+f.example <- function(x, z=2, ...){
+ ## show the manipulated matched call with items from Ctrl attached
+ mc0 <- match.call()
+ ## show that we have access to items of Ctrl
+ print(c(a=a,A=A))
+ ## show the ... argument:
+ if(length(list(...)))
+ cat(" ...", paste(names(list(...)),collapse=", "), "\n",
+ "...", paste(..., collapse=", "),"\n");
+ ##show how formals are matched:
+ print(c(x=x,z=z));
+ ## return a value
+ sin(x)}
+
+### attention does not work so far as a and A are not yet bound:
+# > f.example(2,3)
+# Error in print(c(a = a, A = A)) : object 'a' not found
+
+## fu as S4 object
+fu <- FctWithCtrl( f=f.example, Ctrl=list(a=2,A=5))
+
+## calls to fu
+fu(2,3)
+fu(x=22)
+fu(AA=3,3,.a=3)
+#> ## calls to fu
+#> fu(2,3)
+#a A
+#2 5
+# ...
+# ... 3
+#x z
+#2 2
+#[1] 0.9092974
+#> fu(x=22)
+#a A
+#2 5
+# x z
+#22 2
+#[1] -0.008851309
+#> fu(AA=3,3,.a=3)
+#a A
+#2 5
+# ... AA, .a
+# ... 3 3
+#x z
+#3 2
+#[1] 0.14112
+#>
+
+##manipulate Ctrl
+fu at Ctrl <- list(A=4,a=3)
+fu(AA=3,3,.a=3)
+
+#> fu at Ctrl <- list(A=4,a=3)
+#> fu(AA=3,3,.a=3)
+#a A
+#3 4
+# ... AA, .a
+# ... 3 3
+#x z
+#3 2
+#[1] 0.14112
Added: branches/robast-0.9/pkg/RobExtremesBuffer/README.txt
===================================================================
--- branches/robast-0.9/pkg/RobExtremesBuffer/README.txt (rev 0)
+++ branches/robast-0.9/pkg/RobExtremesBuffer/README.txt 2012-05-21 15:49:11 UTC (rev 476)
@@ -0,0 +1 @@
+some try and error R files ...
\ No newline at end of file
Added: branches/robast-0.9/pkg/RobExtremesBuffer/getOptICGPD.R
===================================================================
--- branches/robast-0.9/pkg/RobExtremesBuffer/getOptICGPD.R (rev 0)
+++ branches/robast-0.9/pkg/RobExtremesBuffer/getOptICGPD.R 2012-05-21 15:49:11 UTC (rev 476)
@@ -0,0 +1,211 @@
+.getpsi <- function(nameInSysdata =".OMSE",
+ PFam = GParetoFamily(),
+ with.correct=FALSE, r=.5){
+
+ sng <- try(getFromNamespace(nameInSysdata, ns = "ROptEst"),silent=TRUE)
+ if(is(sng,"try-error"))
+ stop("Grid of L.M.s for scale-shape family not yet availabe.")
+ sngr <- sng[[name(PFam)]][["grid"]]
+ snf <- sng[[name(PFam)]][["fct"]]
+ ncgr <- ncol(sngr)
+ b <- function(xi0) sng(xi0,1)
+ a <- function(xi0) c(sng(xi0,2),sng(xi0,3))
+ aw <- function(xi0) c(sng(xi0,4),sng(xi0,5))
+ A <- function(xi0) {am <- mean(c(sng(xi0,7),sng(xi0,8)))
+ matrix(c(sng(xi0,6),am,am,sng(xi0,9)),2,2)}
+ Aw <- function(xi0) {am <- mean(c(sng(xi0,11),sng(xi0,12)))
+ matrix(c(sng(xi0,10),am,am,sng(xi0,13)),2,2)}
+
+ normt <- NormType()
+ biast <- symmetricBias()
+ ICT <- paste("optimally robust IC for", switch(nameInSysdata,
+ c(".OMSE"="maxMSE",".RMXE"="RMX", ".MBRE"="maxBias")))
+ riskT <- if(nameInSysdata!=".MBRE") "asGRisk" else "asBias"
+
+ res.xi <- function(xi0){
+ neighbor <- ContNeighborhood(radius = r.0(xi0))
+ w <- new("HampelWeight")
+ stand(w) <- Aw(xi0)
+ cent(w) <- aw(xi0)
+ clip(w) <- b(xi0)
+ if(nameInSysdata!=".MBRE")
+ weight(w) <- getweight(w, neighbor = neighbor, biastype = biast,
+ normW = normt)
+ else weight(w) <- minbiasweight(w, neighbor = neighbor, biastype = biast,
+ normW = normt)
+ res = c(a = a(xi0), A = A(xi0), b = b(xi0), d = 0,
+ normtype = normt, biastype = biast, w = w,
+ info = c("optIC", ICT), risk = riskT,
+ modifyIC = NULL
+ )
+ return(res)
+ }
+ IC.fct.xi.0 <- function(xi0){
+ param0 <- ParamFamParameter(name = "theta", main = c(1,xi0),
+ fixed = 0,
+ trafo = param(PFam)@trafo)
+ PFam at param <- param0
+ Lx <- PFam at L2deriv.fct(param0)
+ IC.fctx <- function(x){
+ Lx.v <- t(cbind(Lx[[1]](x),Lx[[2]](x)))
+ Y <- A(xi0)%*%(Lx.v-a(xi0))*weight(w)(x)
+ return(Y)
+ }
+ return(IC.fctx)
+ }
+ IC.xi.0 <- function(xi0){
+ param0 <- ParamFamParameter(name = "theta", main = c(1,xi0),
+ fixed = 0,
+ trafo = param(PFam)@trafo)
+ PFam at param <- param0
+ res0 <- res(xi0)
+ res0$modifyIC <- function(L2Fam, IC){
+ param1 <- param(L2Fam)
+ xi <- main(param1)["shape"]
+ beta <- main(param1)["scale"]
+ IC <- IC.xi(xi)
+ IC1M <- IC[[1]]@Map
+ IC2M <- list(function(x)IC[[1]]@Map[[1]](x/beta),
+ function(x)IC[[1]]@Map[[2]](x/beta))
+ IC[[1]]@Map <- IC2M
+ return(IC)
+ }
+
+ IC <- generateIC(neighbor= neighbor,
+ L2Fam = PFam, res = res0)
+ return(IC)
+ }
+ IC.fct.xi <- IC.fct.xi.0
+ IC.xi <- IC.xi.0
+ if(with.correct){
+ IC.fct.xi <- function(xi0){
+ res0 <- IC.fct.xi.0(xi0)
+ distr <- distribution(PFam)
+
+ param1 <- param(distribution)
+ param1["shape"] <- xi0
+ param1["scale"] <- 1
+ distr at param <- param1
+
+ I.w <- E(distr, fun=weight(w))
+ I.1 <- E(distr, fun=function(x) weight(w)(x)*Lx[[1](x))
+ I.2 <- E(distr, fun=function(x) weight(w)(x)*Lx[[2](x))
+ z.c <- c(I.1,I.2)/I.w
+
+ L1.0 <- function(x) Lx[[1]](x)-z.c[1]
+ L2.0 <- function(x) Lx[[2]](x)-z.c[2]
+ I.11 <- E(distr, fun=function(x) weight(w)(x)*L1.0(x)^2)
+ I.21 <- E(distr, fun=function(x) weight(w)(x)*L1.0(x)*L2.0(x))
+ I.22 <- E(distr, fun=function(x) weight(w)(x)*L2.0(x)^2)
+
+ A.c <- solve(matrix(c(I.11,I.21,I.21,I.22),2,2))
+ IC.fct <- function(x) A.c%*%(res0$IC.fct(x)-z.c)
+ return(IC.fct)
+ }
+ IC.xi <- function(xi0){
+ res0 <- IC.xi.0(xi0)
+ param0 <- ParamFamParameter(name = "theta", main = c(1,xi0),
+ fixed = 0,
+ trafo = param(PFam)@trafo)
+ PFam at param <- param0
+ IC <- makeIC(res0$IC,PFam)
+ return(IC)
+ }
+ }
+ return(list(IC.fct.xi=IC.fct.xi,IC.xi=IC.xi))
+}
+
+roptestQuick <- function(x, PFam=GParetoFamily(), type=c("RMXE","OMSE","MBRE")){
+ theta0 <- estimate(medkMADhybr(x,PFam=PFam))
+ type <- match.arg(type)
+ nametype <- paste(".",type,sep="")
+ ICfct <- .getpsi(nameInSysdata = nametype,
+ PFam = PFam)$IC.fct.xi
+ IC <- ICfct(xi=theta0["shape"])
+ return(theta0+ rowMeans(IC(x)))
+}
+
+OMSE.5 <- function(x, PFam = GParetoFamily()) roptestQuick(x,PFam,"OMSE")
+RMXE <- function(x, PFam = GParetoFamily()) roptestQuick(x,PFam,"RMXE")
+MBRE <- function(x, PFam = GParetoFamily()) roptestQuick(x,PFam,"MBRE")
+
+roptestScSh <- function(x, PFam=GParetoFamily(), type=c("RMXE","OMSE","MBRE"),
+ initial.est, steps = 1L, verbose = NULL,
+ useLast = getRobAStBaseOption("kStepUseLast"),
+ withUpdateInKer = getRobAStBaseOption("withUpdateInKer"),
+ IC.UpdateInKer = getRobAStBaseOption("IC.UpdateInKer"),
+ withICList = getRobAStBaseOption("withICList"),
+ withPICList = getRobAStBaseOption("withPICList"),
+ na.rm = TRUE, initial.est.ArgList, ...,
+ scalename = "scale", withLogScale = TRUE){
+ if(missing(verbose)|| is.null(verbose))
+ verbose <- getRobAStBaseOption("all.verbose")
+ es.call <- match.call()
+
+ type <- match.arg(type)
+ nametype <- paste(".",type,sep="")
+
+ if(missing(x))
+ stop("'x' is missing with no default")
+ if(missing(PFam))
+ stop("'PFam' 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)
+
+ 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")
+ }
+
+ if(missing(initial.est))
+ initial.est <- medkMADhybr(x,PFam=PFam)
+
+ nrvalues <- length(PFam at param)
+ initial.est <- kStepEstimator.start(initial.est, x = x,
+ nrvalues = nrvalues, na.rm = na.rm,
+ L2Fam = PFam,
+ startList = initial.est.ArgList)
+
+ newParam <- param(PFam)
+ main(newParam)[] <- as.numeric(initial.est)
+ L2FamStart <- modifyModel(PFam, newParam)
+ ICStart <- .getpsi(nameInSysdata = nametype,
+ PFam = PFam, with.correct=TRUE)$IC.xi(main(newParam)[2])
+
+ res <- kStepEstimator(x, IC = ICstart, start = initial.est, steps = steps,
+ useLast = useLast, withUpdateInKer = withUpdateInKer,
+ IC.UpdateInKer = IC.UpdateInKer,
+ withICList = withICList, withPICList = withPICList,
+ na.rm = na.rm, scalename = "scale", withLogScale = TRUE)
+
+ res at estimate.call <- es.call
+ Infos <- matrix(c("roptestScSh",
+ paste(steps, "-step estimate for ", name(PFam), sep = "")),
+ ncol = 2)
+ colnames(Infos) <- c("method", "message")
+
+ if(! distrMod:::.isUnitMatrix(trafo(PFam)))
+ Infos <- rbind(Infos, c("roptestScSh",
+ paste("computation of IC",
+ ifelse(withUpdateInKer,"with","without") ,
+ "modification in ker(trafo)")))
+
+ Infos <- rbind(Infos, c("roptestScSh",
+ paste("computation of IC, asvar and asbias via useLast =", useLast)))
+ Infos(res) <- Infos
+ res at completecases <- completecases
+ res at start <- initial.est
+ return(res)
+}
Added: branches/robast-0.9/pkg/RobExtremesBuffer/internal.roptest.R
===================================================================
--- branches/robast-0.9/pkg/RobExtremesBuffer/internal.roptest.R (rev 0)
+++ branches/robast-0.9/pkg/RobExtremesBuffer/internal.roptest.R 2012-05-21 15:49:11 UTC (rev 476)
@@ -0,0 +1,145 @@
+.fix.in.defaults <- function(call.list, fun){
+ formals.fun <- formals(fun)
+ k <- length(call.list)
+ L <- length(formals.fun)
+ if("..." %in% names(formals.fun)) L <- L-1
+ for(i in 1:L){
+ if(!is(formals.fun[[i]],"name")){
+ if(!names(formals.fun)[i] %in% names(call.list)&&!is.null(formals.fun[[i]])){
+ k <- k + 1
+ call.list[[k]] <- formals.fun[[i]]
+ names(call.list)[k] <- names(formals.fun)[i]
+ }
+ }
+ }
+ 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)
+}
+.check.eps <- function(...){
+ mc <- match.call(expand=TRUE)
+
+ eps <- eps.lower <- eps.upper <- NULL
+ if(is.null(mc$eps) && is.null(mc$eps.lower) && is.null(mc$eps.upper)){
+ eps.lower <- 0
+ eps.upper <- 0.5
+ }
+ if(is.null(mc$eps)){
+ if(!is.null(mc$eps.lower) && is.null(mc$eps.upper))
+ eps.upper <- 0.5
+ if(is.null(mc$eps.lower) && !is.null(mc$eps.upper))
+ eps.lower <- 0
+ 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((eps.lower < 0) || (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(eps == 0)
+ stop("'eps = 0'! => use functions 'mean' and 'sd' for estimation")
+ if((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))
+
+ return(list(e=eps,lower=eps.lower, upper=eps.upper, sqn = sqrtn))
+}
+
+.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(...){
+ mc <- match.call(expand=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 <- 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(!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 <- do.call(optIC, c(list( model = mc$infMod, risk = mc$risk,
+ verbose = mc$verbose, OptOrIter = mc$OptOrIter),
+ dots))
+ }
+ }else{
+ neighbor at radius <- eps$sqn*eps$e*mc$fsCor
+ infMod <- InfRobModel(center = mc$L2FamStart, neighbor = mc$neighbor)
+ ICstart <- do.call(optIC, c(list(model = mc$infMod, risk = mc$risk,
+ verbose = mc$verbose, OptOrIter = mc$OptOrIter),
+ dots))
+ }
+ 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,genstartCtrl)
+ return(es.list)
+}
\ No newline at end of file
Added: branches/robast-0.9/pkg/RobExtremesBuffer/interpolGrids.R
===================================================================
--- branches/robast-0.9/pkg/RobExtremesBuffer/interpolGrids.R (rev 0)
+++ branches/robast-0.9/pkg/RobExtremesBuffer/interpolGrids.R 2012-05-21 15:49:11 UTC (rev 476)
@@ -0,0 +1,68 @@
+###generate first InterpolGrids
+if(FALSE){
+### if grid were dense enough....
+.myFolder <- "C:/rtest/RobASt/branches/robast-0.9/pkg/ROptEst/R"
+.myfolData <- "D:/Documents/Arbeit/Projekte/NataliyaDiss/Paper R-files/Extrapolation"
+newEnv1 <- new.env()
+load(file.path(.myfolData,"ExtrapolationMC_RMX1.RData"),envir=newEnv1)
+whatIn <- ls(envir=newEnv1, all.names=TRUE)
+xi <- get("xi.g",envir=newEnv1)
+Y <- get("Y",envir=newEnv1)
+rm(list=whatIn,envir=newEnv1)
+gc()
+.saveInterpGrid(xiGrid = xi,
+ PFam = GParetoFamily(),
+ sysRdaFolder = .myFolder,
+ getFun = getLMGrid, optFct = .RMXE.xi, nameInSysdata = ".RMXE",
+ withSmooth = TRUE, withPrint = TRUE,
+ Y=Y, elseFun= .MakeGridList)
+
+newEnv1 <- new.env()
+.myfolData <- "D:/Documents/Arbeit/Projekte/NataliyaDiss/Paper JOP"
+load(file.path(.myfolData,"ExtrapolationMC_OMSE.RData"),envir=newEnv1)
+whatIn <- ls(envir=newEnv1, all.names=TRUE)
+xi0 <- get("xi.g",envir=newEnv1)
+Y0 <- get("Y",envir=newEnv1)
+rm(list=whatIn,envir=newEnv1)
+gc()
+load(file.path(.myfolData,"ExtrapolationMC_OMSE_infmean.RData"),envir=newEnv1)
+whatIn <- ls(envir=newEnv1, all.names=TRUE)
+xi1 <- get("xi.g",envir=newEnv1)
+Y1 <- get("Y",envir=newEnv1)
+rm(list=whatIn,envir=newEnv1)
+gc()
+xi <- rbind(xi0,xi1)
+on <- order(xi)
+xi <- xi[on]
+Y <- (rbind(Y0[,1,],Y1[,1,]))[on,]
+.saveInterpGrid(xiGrid = xi,
+ PFam = GParetoFamily(),
+ sysRdaFolder = .myFolder,
+ getFun = getLMGrid, optFct = .OMSE.xi, nameInSysdata = ".OMSE",
+ withSmooth = TRUE, withPrint = TRUE,
+ Y=Y, elseFun= .MakeGridList)
+
+newEnv1 <- new.env()
+load(file.path(.myfolData,"ExtrapolationMC_OBRE.RData"),envir=newEnv1)
+whatIn <- ls(envir=newEnv1, all.names=TRUE)
+xi0 <- get("xi.g",envir=newEnv1)
+Y0 <- get("Y",envir=newEnv1)
+rm(list=whatIn,envir=newEnv1)
+gc()
+load(file.path(.myfolData,"ExtrapolationMC_OBRE_infmean.RData"),envir=newEnv1)
+whatIn <- ls(envir=newEnv1, all.names=TRUE)
+xi1 <- get("xi.g",envir=newEnv1)
+Y1 <- get("Y",envir=newEnv1)
+rm(list=whatIn,envir=newEnv1)
+gc()
+xi <- rbind(xi0,xi1)
+on <- order(xi)
+xi <- xi[on]
+Y <- (rbind(Y0[,1,],Y1[,1,]))[on,]
+.saveInterpGrid(xiGrid = xi,
+ PFam = GParetoFamily(),
+ sysRdaFolder = .myFolder,
+ getFun = getLMGrid, optFct = .MBRE.xi, nameInSysdata = ".MBRE",
+ withSmooth = TRUE, withPrint = TRUE,
+ Y=Y, elseFun= .MakeGridList)
+}
Added: branches/robast-0.9/pkg/RobExtremesBuffer/roptest.new.R
===================================================================
--- branches/robast-0.9/pkg/RobExtremesBuffer/roptest.new.R (rev 0)
+++ branches/robast-0.9/pkg/RobExtremesBuffer/roptest.new.R 2012-05-21 15:49:11 UTC (rev 476)
@@ -0,0 +1,89 @@
+###############################################################################
+## Optimally robust estimation
+###############################################################################
+roptest <- function(x, L2Fam, fsCor = 1,
+ risk = asMSE(), steps = 1L,
+ verbose = NULL,
+ OptOrIter = "iterate",
+ nbCtrl = gennbCtrl(neighbor = ContNeighborhood(),
+ eps, eps.lower, eps.upper)
+ startCtrl = genstartCtrl(initial.est = NULL,
+ initial.est.ArgList = NULL,
+ startPar = NULL, distance = CvMDist),
+ kstepCtrl = genkstepCtrl(
+ useLast = getRobAStBaseOption("kStepUseLast"),
+ withUpdateInKer = getRobAStBaseOption("withUpdateInKer"),
+ IC.UpdateInKer = getRobAStBaseOption("IC.UpdateInKer"),
+ withICList = getRobAStBaseOption("withICList"),
+ withPICList = getRobAStBaseOption("withPICList"),
+ withLogScale = TRUE),
+ na.rm = TRUE, ...){
+
+ es.call <- match.call()
+ dots <- match.call(expand=FALSE)$dots
+ es.list <- as.list(es.call[-1])
+ es.list <- .fix.in.defaults(es.list,roptest)
+ es.list <- c(es.list,nbCtrl)
+ es.list$dots <- dots
+
+ if(missing(verbose)|| is.null(verbose))
+ es.list$verbose <- getRobAStBaseOption("all.verbose")
+
+ if(missing(L2Fam))
+ stop("'L2Fam' is missing with no default")
+
+ x <- .pretreat(x,na.rm)
+
+ es.list$eps <- do.call(.check.eps, args=es.list)
+
+ .isOKfsCor(fsCor)
+
+ .isOKsteps(steps)
+
+ 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)
+ initial.est <- kStepEstimator.start(initial.est, x = x,
+ nrvalues = nrvalues, na.rm = na.rm,
+ L2Fam = L2Fam,
+ startList = startCtrl$initial.est.ArgList)
+
+
+ newParam <- param(L2Fam)
+ main(newParam)[] <- as.numeric(initial.est)
+ L2FamStart <- modifyModel(L2Fam, newParam)
+
+ ICstart <- do.call(.getROptICstart, args=es.list)
+
+ 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)
+
+
+ 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("roptest",
+ paste("computation of IC",
+ ifelse(withUpdateInKer,"with","without") ,
+ "modification in ker(trafo)")))
+
+ Infos <- rbind(Infos, c("roptest",
+ paste("computation of IC, asvar and asbias via useLast =", useLast)))
+ Infos(res) <- Infos
+ res at completecases <- completecases
+ res at start <- initial.est
+ return(res)
+}
More information about the Robast-commits
mailing list