[Robast-commits] r533 - in branches/robast-0.9/pkg/ROptEst: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Jan 15 01:31:11 CET 2013
Author: ruckdeschel
Date: 2013-01-15 01:31:10 +0100 (Tue, 15 Jan 2013)
New Revision: 533
Added:
branches/robast-0.9/pkg/ROptEst/R/versionSuff.R
branches/robast-0.9/pkg/ROptEst/man/getStartIC-methods.Rd
branches/robast-0.9/pkg/ROptEst/man/inputGenerator.Rd
branches/robast-0.9/pkg/ROptEst/man/internalInterpolate.Rd
branches/robast-0.9/pkg/ROptEst/man/internalMBRE.Rd
branches/robast-0.9/pkg/ROptEst/man/internalRobestHelpers.Rd
branches/robast-0.9/pkg/ROptEst/man/internal_Cniperplots.Rd
branches/robast-0.9/pkg/ROptEst/man/robest.Rd
Modified:
branches/robast-0.9/pkg/ROptEst/NAMESPACE
branches/robast-0.9/pkg/ROptEst/R/AllGeneric.R
branches/robast-0.9/pkg/ROptEst/R/AllPlot.R
branches/robast-0.9/pkg/ROptEst/R/cniperCont.R
branches/robast-0.9/pkg/ROptEst/R/comparePlot.R
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/roptest.R
branches/robast-0.9/pkg/ROptEst/R/roptest.new.R
branches/robast-0.9/pkg/ROptEst/man/0ROptEst-package.Rd
branches/robast-0.9/pkg/ROptEst/man/cniperCont.Rd
branches/robast-0.9/pkg/ROptEst/man/comparePlot.Rd
branches/robast-0.9/pkg/ROptEst/man/plot-methods.Rd
branches/robast-0.9/pkg/ROptEst/man/roptest.Rd
Log:
ROptEst: first warningless version with new robest (and everything documented)
Modified: branches/robast-0.9/pkg/ROptEst/NAMESPACE
===================================================================
--- branches/robast-0.9/pkg/ROptEst/NAMESPACE 2013-01-14 15:23:03 UTC (rev 532)
+++ branches/robast-0.9/pkg/ROptEst/NAMESPACE 2013-01-15 00:31:10 UTC (rev 533)
@@ -25,11 +25,14 @@
"lowerCaseRadius",
"minmaxBias", "getBiasIC",
"getL1normL2deriv",
- "getModifyIC",
- "cniperCont", "cniperPoint", "cniperPointPlot")
+ "getModifyIC")
exportMethods("updateNorm", "scaleUpdateIC", "eff",
- "get.asGRisk.fct", "getStartIC", "plot")
+ "get.asGRisk.fct", "getStartIC", "plot",
+ "comparePlot")
export("getL2normL2deriv",
"asAnscombe", "asL1", "asL4",
"getReq", "getMaxIneff")
-export("roptest","getLagrangeMultByOptim","getLagrangeMultByIter")
+export("roptest","roptest.old", "robest",
+ "getLagrangeMultByOptim","getLagrangeMultByIter")
+export("genkStepCtrl", "genstartCtrl", "gennbCtrl")
+export("cniperCont", "cniperPoint", "cniperPointPlot")
\ No newline at end of file
Modified: branches/robast-0.9/pkg/ROptEst/R/AllGeneric.R
===================================================================
--- branches/robast-0.9/pkg/ROptEst/R/AllGeneric.R 2013-01-14 15:23:03 UTC (rev 532)
+++ branches/robast-0.9/pkg/ROptEst/R/AllGeneric.R 2013-01-15 00:31:10 UTC (rev 533)
@@ -77,15 +77,6 @@
if(!isGeneric("scaleUpdateIC")){
setGeneric("scaleUpdateIC", function(neighbor, ...) standardGeneric("scaleUpdateIC"))
}
-if(!isGeneric("cniperCont")){
- setGeneric("cniperCont", function(IC1, IC2, L2Fam, neighbor, risk, ...) standardGeneric("cniperCont"))
-}
-if(!isGeneric("cniperPoint")){
- setGeneric("cniperPoint", function(L2Fam, neighbor, risk, ...) standardGeneric("cniperPoint"))
-}
-if(!isGeneric("cniperPointPlot")){
- setGeneric("cniperPointPlot", function(L2Fam, neighbor, risk, ...) standardGeneric("cniperPointPlot"))
-}
if(!isGeneric("eff")){
setGeneric("eff", function(object) standardGeneric("eff"))
}
Modified: branches/robast-0.9/pkg/ROptEst/R/AllPlot.R
===================================================================
--- branches/robast-0.9/pkg/ROptEst/R/AllPlot.R 2013-01-14 15:23:03 UTC (rev 532)
+++ branches/robast-0.9/pkg/ROptEst/R/AllPlot.R 2013-01-15 00:31:10 UTC (rev 533)
@@ -6,13 +6,13 @@
with.legend = FALSE, legend = NULL, legend.bg = "white",
legend.location = "bottomright", legend.cex = 0.8,
withMBR = FALSE, MBRB = NA, MBR.fac = 2, col.MBR = par("col"),
- lty.MBR = "dashed", lwd.MBR = 0.8,
+ lty.MBR = "dashed", lwd.MBR = 0.8, n.MBR = 10000,
scaleX = FALSE, scaleX.fct, scaleX.inv,
scaleY = FALSE, scaleY.fct = pnorm, scaleY.inv=qnorm,
scaleN = 9, x.ticks = NULL, y.ticks = NULL,
mfColRow = TRUE, to.draw.arg = NULL){
- mcl <- match.call(call = sys.call(sys.parent(1)))
+ mcl <- match.call(call = sys.call(sys.parent(1)), expand.dots = TRUE)
L2Fam <- eval(x at CallL2Fam); trafO <- trafo(L2Fam at param)
dims <- nrow(trafO); to.draw <- 1:dims
@@ -26,23 +26,24 @@
MBRB <- matrix(rep(t(MBRB), length.out=dims0*2),ncol=2, byrow=T)
if(withMBR && all(is.na(MBRB))){
- robModel <- InfRobModel(center = L2fam, neighbor =
+ robModel <- InfRobModel(center = L2Fam, neighbor =
ContNeighborhood(radius = 0.5))
ICmbr <- try(optIC(model = robModel, risk = asBias()), silent=TRUE)
if(!is(ICmbr,"try-error"))
- MBRB <- .getExtremeCoordIC(ICmbr, distribution(L2Fam), todraw)
+ MBRB <- .getExtremeCoordIC(ICmbr, distribution(L2Fam), to.draw,
+ n = n.MBR)
else withMBR <- FALSE
}
mcl$MBRB <- MBRB
mcl$withMBR <- withMBR
- do.call(getMethod("plot", signature(x = "IC", y = "missing"),
- where="RobAStBase"), mcl)
- })
+ plm <- getMethod("plot", signature(x = "IC", y = "missing"),
+ where="RobAStBase")
+ do.call(plm, as.list(mcl[-1]), envir=parent.frame(2))
+ return(invisible())
+ })
-.getExtremeCoordIC <- function(IC, D, indi, n = 50000){
+.getExtremeCoordIC <- function(IC, D, indi, n = 10000){
x <- q(D)(seq(1/2/n,1-1/2/n, length=n))
- li <- length(indi)
- ICx <- matrix(0,li,n)
- for( i in 1:li) ICx[i,] <- sapply(x, IC at Map[[indi[i]]])
- return(cbind(min=apply(ICx,1,min),max=apply(ICx,1,max)))
+ y <- (matrix(evalIC(IC,matrix(x,ncol=1)),ncol=n))[indi,]
+ return(cbind(min=apply(y,1,min),max=apply(y,1,max)))
}
\ No newline at end of file
Modified: branches/robast-0.9/pkg/ROptEst/R/cniperCont.R
===================================================================
--- branches/robast-0.9/pkg/ROptEst/R/cniperCont.R 2013-01-14 15:23:03 UTC (rev 532)
+++ branches/robast-0.9/pkg/ROptEst/R/cniperCont.R 2013-01-15 00:31:10 UTC (rev 533)
@@ -1,118 +1,7 @@
-PETER <- FALSE
-setMethod("cniperCont", signature(IC1 = "IC",
- IC2 = "IC",
- L2Fam = "L2ParamFamily",
- neighbor = "ContNeighborhood",
- risk = "asMSE"),
- function(IC1, IC2, L2Fam, neighbor, risk, lower, upper, n = 101){
- R1 <- Risks(IC1)[["trAsCov"]]
- if(is.null(R1)) R1 <- getRiskIC(IC1, risk = trAsCov(), L2Fam = L2Fam)
- if(length(R1) > 1) R1 <- R1$value
-
- R2 <- Risks(IC2)[["trAsCov"]]
- if(is.null(R2)) R2 <- getRiskIC(IC2, risk = trAsCov(), L2Fam = L2Fam)
- if(length(R2) > 1) R2 <- R2$value
-
- r <- neighbor at radius
-
- fun <- function(x){
- y1 <- evalIC(IC1, x)
- y2 <- evalIC(IC2, x)
- R1 - R2 + r^2*(as.vector(y1 %*% y1) - as.vector(y2 %*% y2))
- }
- x <- seq(from = lower, to = upper, length = n)
- y <- sapply(x, fun)
- plot(x, y, type = "l", main = "Cniper region plot",
- xlab = "Dirac point", ylab = "Asymptotic MSE difference (IC1 - IC2)")
-# text(min(x), max(y)/2, "IC2", pos = 4)
-# text(min(x), min(y)/2, "IC1", pos = 4)
- abline(h = 0)
- invisible()
- })
-
-setMethod("cniperPoint", signature(L2Fam = "L2ParamFamily",
- neighbor = "ContNeighborhood",
- risk = "asMSE"),
- function(L2Fam, neighbor, risk, lower, upper){
- D <- trafo(L2Fam at param)
- tr.invF <- sum(diag(D %*% solve(FisherInfo(L2Fam)) %*% t(D)))
- psi <- optIC(model = L2Fam, risk = asCov())
- robMod <- InfRobModel(center = L2Fam, neighbor = neighbor)
- eta <- optIC(model = robMod, risk = asMSE())
- maxMSE <- Risks(eta)$asMSE$value
- Delta <- sqrt(maxMSE - tr.invF)/neighbor at radius
- fun <- function(x){
- y <- evalIC(psi, x)
- sqrt(as.vector(y %*% y)) - Delta
- }
- res <- uniroot(fun, lower = lower, upper = upper)$root
- names(res) <- "cniper point"
- res
- })
-
-setMethod("cniperPointPlot", signature(L2Fam = "L2ParamFamily",
- neighbor = "ContNeighborhood",
- risk = "asMSE"),
- function(L2Fam, neighbor, risk, lower, upper, n = 101, ...,
- ){
- dots <- as.list(match.call(call = sys.call(sys.parent(1)),
- expand.dots = FALSE)$"...")
- D <- trafo(L2Fam at param)
- tr.invF <- sum(diag(D %*% solve(FisherInfo(L2Fam)) %*% t(D)))
- psi <- optIC(model = L2Fam, risk = asCov())
- robMod <- InfRobModel(center = L2Fam, neighbor = neighbor)
- eta <- optIC(model = robMod, risk = asMSE())
- maxMSE <- Risks(eta)$asMSE$value
- fun <- function(x){
- y <- evalIC(psi, x)
- tr.invF + as.vector(y %*% y)*neighbor at radius^2 - maxMSE
- }
- dots$x <- x <- seq(from = lower, to = upper, length = n)
- dots$y <- sapply(x, fun)
- colSet <- ltySet <- lwdSet <- FALSE
- if(!is.null(dots$col)) {colSet <- TRUE; colo <- eval(dots$col)}
- if(colSet) {
- colo <- rep(colo,length.out=2)
- dots$col <- colo[1]
- }
- if(!is.null(dots$lwd)) {lwdSet <- TRUE; lwdo <- eval(dots$lwd)}
- if(lwdSet) {
- lwdo <- rep(lwdo,length.out=2)
- dots$lwd <- lwdo[1]
- }
- if(!is.null(dots$lty)) {ltySet <- TRUE; ltyo <- eval(dots$lty)}
- if(ltySet && ((!is.numeric(ltyo) && length(ltyo)==1)||
- is.numeric(ltyo))){
- ltyo <- list(ltyo,ltyo)
- dots$lty <- ltyo[[1]]
- }else{ if (ltySet && !is.numeric(ltyo) && length(ltyo)==2){
- dots$lty <- ltyo[[1]]
- }
- }
- if(is.null(dots$main)) dots$main <- gettext("Cniper point plot")
- if(is.null(dots$xlab)) dots$xlab <- gettext("Dirac point")
- if(is.null(dots$ylab))
- dots$ylab <- gettext("Asymptotic MSE difference (classic - robust)")
- dots$type <- "l"
-
-
-
-
- do.call(plot, dots)
-# text(min(x), max(y)/2, "Robust", pos = 4)
-# text(min(x), min(y)/2, "Classic", pos = 4)
- dots$x <- dots$y <- dots$xlab <- dots$ylab <- dots$main <- dots$type <- NULL
- dots$h <- 0
- if(colSet) dots$col <- colo[2]
- if(lwdSet) dots$lwd <- lwdo[2]
- if(ltySet) dots$lty <- ltyo[[2]]
- do.call(abline, dots)
- invisible()
- })
-
-if(PETER){
.rescalefct <- RobAStBase:::.rescalefct
+.plotRescaledAxis <- RobAStBase:::.plotRescaledAxis
.makedotsP <- RobAStBase:::.makedotsP
+.makedotsLowLevel <- RobAStBase:::.makedotsLowLevel
.SelectOrderData <- RobAStBase:::.SelectOrderData
.plotData <- function(
@@ -125,25 +14,29 @@
IC # IC1 in cniperContPlot and eta in cniperPointPlot
){
dotsP <- .makedotsP(dots)
- dotsP$col <- rep(origCl$col.pts, length.out=n)
- dotsP$pch <- rep(origCl$pch.pts, length.out=n)
+ dotsP$col <- rep(eval(origCl$col.pts), length.out=n)
+ dotsP$pch <- rep(eval(origCl$pch.pts), length.out=n)
+ al <- eval(origCl$alpha.trsp)
+ if(!is.na(al))
+ dotsP$col <- sapply(dotsP$col, addAlphTrsp2col, alpha=al)
+
n <- if(!is.null(dim(data))) nrow(data) else length(data)
if(!is.null(lab.pts))
lab.pts <- rep(origCl$lab.pts, length.out=n)
sel <- .SelectOrderData(data, function(x)sapply(x,fun),
- origCl$which.lbs, origCl$which.Order)
+ eval(origCl$which.lbs),
+ eval(origCl$which.Order))
i.d <- sel$ind
i0.d <- sel$ind1
y.d <- sel$y
x.d <- sel$data
- x.d <- sel.C$data
n <- length(i.d)
resc.dat <- .rescalefct(x.d, function(x) sapply(x,fun),
- origCl$scaleX, origCl$scaleX.fct, origCl$scaleX.inv,
- origCl$scaleY, origCl$scaleY.fct,
+ eval(origCl$scaleX), origCl$scaleX.fct, origCl$scaleX.inv,
+ eval(origCl$scaleY), origCl$scaleY.fct,
dots$xlim, dots$ylim, dots)
dotsP$x <- resc.dat$X
@@ -152,28 +45,60 @@
trafo <- trafo(L2Fam at param)
dims <- nrow(trafo)
QF <- diag(dims)
- if(is(IC1,"ContIC") & dims>1 )
- {if (is(normtype(object),"QFNorm"))
- QF <- QuadForm(normtype(object))}
+ if(is(IC,"ContIC") & dims>1 )
+ {if (is(normtype(IC),"QFNorm"))
+ QF <- QuadForm(normtype(IC))}
absInfoEval <- function(x,y) sapply(x, y at Map[[1]])
- IC1.rv <- as(diag(dims) %*% IC at Curve, "EuclRandVariable")
- absy.f <- t(IC1.rv) %*% QF %*% IC1.rv
+ IC.rv <- as(diag(dims) %*% IC at Curve, "EuclRandVariable")
+ absy.f <- t(IC.rv) %*% QF %*% IC.rv
absy <- absInfoEval(x.d, absy.f)
- dotsP$cex <- log(absy+1)*3*rep(cex.pts, length.out=n)
+ if(is.null(origCl$cex.pts)) origCl$cex.pts <- par("cex")
+ dotsP$cex <- log(absy+1)*3*rep(origCl$cex.pts, length.out=n)
dotsT <- dotsP
dotsT$pch <- NULL
dotsT$cex <- dotsP$cex/2
dotsT$labels <- if(is.null(lab.pts)) i.d else lab.pts[i.d]
do.call(points,dotsP)
- if(with.lab) do.call(text,dotsT)
- if(return.Order) return(i0.d)
+ if(!is.null(origCl$with.lab))
+ if(origCl$with.lab) do.call(text,dotsT)
+ if(!is.null(origCl$return$order))
+ if(origCl$return.Order) return(i0.d)
return(invisible(NULL))
}
-cniperContPlot <- function(IC1, IC2, data = NULL, ...,
+
+.getFunCnip <- function(IC1,IC2, risk, L2Fam, r, b20=NULL){
+
+ riskfct <- getRiskFctBV(risk, biastype(risk))
+
+ .getTrVar <- function(IC){
+ R <- Risks(IC)[["trAsCov"]]
+ if(is.null(R)) R <- getRiskIC(IC, risk = trAsCov(), L2Fam = L2Fam)
+ if(length(R) > 1) R <- R$value
+ return(R)
+ }
+ R1 <- .getTrVar (IC1)
+ R2 <- .getTrVar (IC2)
+
+
+ fun <- function(x){
+ y1 <- evalIC(IC1,as.matrix(x,ncol=1))
+ r1 <- riskfct(var=R1,bias=r*fct(normtype(risk))(y1))
+ if(!is.null(b20))
+ r2 <- riskfct(var=R1,bias=b20) else{
+ y2 <- sapply(x,function(x0) evalIC(IC2,x0))
+ r2 <- riskfct(var=R2,bias=r*fct(normtype(risk))(y2))
+ }
+ r1 - r2
+ }
+
+ return(fun)
+}
+
+cniperCont <- function(IC1, IC2, data = NULL, ...,
neighbor, risk, lower=getdistrOption("DistrResolution"),
upper=1-getdistrOption("DistrResolution"), n = 101,
scaleX = FALSE, scaleX.fct, scaleX.inv,
@@ -181,53 +106,48 @@
scaleN = 9, x.ticks = NULL, y.ticks = NULL,
cex.pts = 1, col.pts = par("col"),
pch.pts = 1, jitter.fac = 1, with.lab = FALSE,
- lab.pts = NULL, lab.font = NULL,
+ lab.pts = NULL, lab.font = NULL, alpha.trsp = NA,
which.lbs = NULL, which.Order = NULL,
return.Order = FALSE){
- mc <- match.call(call = sys.call(sys.parent(1)),
- expand.dots = FALSE)
- dots <- mc$"..."
-
+ mc <- match.call(expand.dots = FALSE)
+ dots <- as.list(mc$"...")
if(!is(IC1,"IC")) stop ("IC1 must be of class 'IC'")
if(!is(IC2,"IC")) stop ("IC2 must be of class 'IC'")
if(!identical(IC1 at CallL2Fam, IC2 at CallL2Fam))
stop("IC1 and IC2 must be defined on the same model")
L2Fam <- eval(IC1 at CallL2Fam)
- riskfct <- getRiskBV(risk, biastype(risk))
+ b20 <- NULL
+ fCpl <- eval(dots$fromCniperPlot)
+ if(!is.null(fCpl))
+ if(fCpl) b20 <- neighbor at radius*Risks(IC2)$asBias$value
+ dots$fromCniperPlot <- NULL
+
+ fun <- .getFunCnip(IC1,IC2, risk, L2Fam, neighbor at radius, b20)
if(missing(scaleX.fct)){
scaleX.fct <- p(L2Fam)
scaleX.inv <- q(L2Fam)
}
- R1 <- Risks(IC1)[["trAsCov"]]
- if(is.null(R1)) R1 <- getRiskIC(IC1, risk = trAsCov(), L2Fam = L2Fam)
- if(length(R1) > 1) R1 <- R1$value
-
- R2 <- Risks(IC2)[["trAsCov"]]
- if(is.null(R2)) R2 <- getRiskIC(IC2, risk = trAsCov(), L2Fam = L2Fam)
- if(length(R2) > 1) R2 <- R2$value
-
- r <- neighbor at radius
-
- fun <- function(x){
- y1 <- evalIC(IC1, x)
- y2 <- evalIC(IC2, x)
- riskfct(R1,r*fct(normtype(risk))(y1))-
- riskfct(R2,r*fct(normtype(risk))(y2))
- }
+ if(!is.null(as.list(mc)$lower)) lower <- p(L2Fam)(lower)
+ if(!is.null(as.list(mc)$upper)) upper <- p(L2Fam)(upper)
x <- q(L2Fam)(seq(lower,upper,length=n))
- resc <- .rescalefct(x, function(u) sapply(u,fun), scaleX, scaleX.fct,
+ if(is(distribution(L2Fam), "DiscreteDistribution"))
+ x <- seq(q(L2Fam)(lower),q(L2Fam)(upper),length=n)
+ resc <- .rescalefct(x, fun, scaleX, scaleX.fct,
scaleX.inv, scaleY, scaleY.fct, dots$xlim, dots$ylim, dots)
- x <- dots$x <- resc$X
+ dots$x <- resc$X
dots$y <- resc$Y
+
+
dots$type <- "l"
- if(is.null(dots$main)) dots$main <- "Cniper region plot"
- if(is.null(dots$xlab)) dots$xlab <- "Dirac point"
- if(is.null(dots$ylab)) dots$ylab <- "Asymptotic Risk difference (IC1 - IC2)"
+ if(is.null(dots$main)) dots$main <- gettext("Cniper region plot")
+ if(is.null(dots$xlab)) dots$xlab <- gettext("Dirac point")
+ if(is.null(dots$ylab))
+ dots$ylab <- gettext("Asymptotic Risk difference (IC1 - IC2)")
colSet <- ltySet <- lwdSet <- FALSE
if(!is.null(dots$col)) {colSet <- TRUE; colo <- eval(dots$col)}
@@ -249,10 +169,9 @@
dots$lty <- ltyo[[1]]
}
}
-
do.call(plot,dots)
- dots <- makedotsLowLevel(dots)
+ dots <- .makedotsLowLevel(dots)
dots$x <- dots$y <- NULL
if(colSet) dots$col <- colo[2]
if(lwdSet) dots$lwd <- lwdo[2]
@@ -262,7 +181,7 @@
do.call(abline, dots)
.plotRescaledAxis(scaleX, scaleX.fct, scaleX.inv, scaleY,scaleY.fct,
- scaleY.inv, dots$xlim, dots$ylim, x, ypts = 400,
+ scaleY.inv, dots$xlim, dots$ylim, resc$X, ypts = 400,
n = scaleN, x.ticks = x.ticks, y.ticks = y.ticks)
if(!is.null(data))
return(.plotData(data, dots, mc, fun, L2Fam, IC1))
@@ -273,30 +192,21 @@
lower=getdistrOption("DistrResolution"),
upper=1-getdistrOption("DistrResolution")){
- lower.x <- q(L2Fam)(lower)
- upper.x <- q(L2Fam)(upper)
- riskfct <- getRiskBV(risk, biastype(risk))
- psi <- optIC(model = L2Fam, risk = asCov())
- R1 <- Risks(psi)[["trAsCov"]]
- if(is.null(R1)) R1 <- getRiskIC(psi, risk = trAsCov(), L2Fam = L2Fam)
- if(length(R1) > 1) R1 <- R1$value
+ mc <- match.call(expand.dots = FALSE)
+ if(!is.null(as.list(mc)$lower)) lower <- p(L2Fam)(lower)
+ if(!is.null(as.list(mc)$upper)) upper <- p(L2Fam)(upper)
+ lower <- q(L2Fam)(lower)
+ upper <- q(L2Fam)(upper)
+
robMod <- InfRobModel(center = L2Fam, neighbor = neighbor)
+
+ psi <- optIC(model = L2Fam, risk = asCov())
eta <- optIC(model = robMod, risk = risk)
- R2 <- Risks(eta)[["trAsCov"]]
- if(is.null(R2)) R2 <- getRiskIC(eta, risk = trAsCov(), L2Fam = L2Fam)
- if(length(R2) > 1) R2 <- R2$value
- r <- neighbor at radius
+ fun <- .getFunCnip(psi,eta, risk, L2Fam, neighbor at radius)
- fun <- function(x){
- y1 <- evalIC(psi, x)
- y2 <- evalIC(eta, x)
- riskfct(R1,r*fct(normtype(risk))(y1))-
- riskfct(R2,r*fct(normtype(risk))(y2))
- }
-
res <- uniroot(fun, lower = lower, upper = upper)$root
names(res) <- "cniper point"
res
@@ -305,97 +215,34 @@
cniperPointPlot <- function(L2Fam, data=NULL, ..., neighbor, risk= asMSE(),
lower=getdistrOption("DistrResolution"),
upper=1-getdistrOption("DistrResolution"), n = 101,
+ withMaxRisk = TRUE,
scaleX = FALSE, scaleX.fct, scaleX.inv,
scaleY = FALSE, scaleY.fct = pnorm, scaleY.inv=qnorm,
scaleN = 9, x.ticks = NULL, y.ticks = NULL,
cex.pts = 1, col.pts = par("col"),
pch.pts = 1, jitter.fac = 1, with.lab = FALSE,
- lab.pts = NULL, lab.font = NULL,
+ lab.pts = NULL, lab.font = NULL, alpha.trsp = NA,
which.lbs = NULL, which.Order = NULL,
return.Order = FALSE){
mc <- match.call(call = sys.call(sys.parent(1)),
expand.dots = FALSE)
- dots <- mc$"..."
+ mcl <- as.list(mc[-1])
+ dots <- as.list(mc$"...")
- if(missing(scaleX.fct)){
- scaleX.fct <- p(L2Fam)
- scaleX.inv <- q(L2Fam)
- }
-
- lower.x <- q(L2Fam)(lower)
- upper.x <- q(L2Fam)(upper)
- riskfct <- getRiskBV(risk, biastype(risk))
-
- psi <- optIC(model = L2Fam, risk = asCov())
- R1 <- Risks(psi)[["trAsCov"]]
- if(is.null(R1)) R1 <- getRiskIC(psi, risk = trAsCov(), L2Fam = L2Fam)
- if(length(R1) > 1) R1 <- R1$value
-
robMod <- InfRobModel(center = L2Fam, neighbor = neighbor)
- eta <- optIC(model = robMod, risk = risk)
- R2 <- Risks(eta)[["trAsCov"]]
- if(is.null(R2)) R2 <- getRiskIC(eta, risk = trAsCov(), L2Fam = L2Fam)
- if(length(R2) > 1) R2 <- R2$value
- r <- neighbor at radius
-
- fun <- function(x){
- y1 <- evalIC(psi, x)
- y2 <- evalIC(eta, x)
- riskfct(R1,r*fct(normtype(risk))(y1))-
- riskfct(R2,r*fct(normtype(risk))(y2))
- }
-
-
- x <- q(L2Fam)(seq(lower,upper,length=n))
- resc <- .rescalefct(x, function(u) sapply(u,fun), scaleX, scaleX.fct,
- scaleX.inv, scaleY, scaleY.fct, dots$xlim, dots$ylim, dots)
- x <- dots$x <- resc$X
- dots$y <- resc$Y
-
- colSet <- ltySet <- lwdSet <- FALSE
- if(!is.null(dots$col)) {colSet <- TRUE; colo <- eval(dots$col)}
- if(colSet) {
- colo <- rep(colo,length.out=2)
- dots$col <- colo[1]
- }
- if(!is.null(dots$lwd)) {lwdSet <- TRUE; lwdo <- eval(dots$lwd)}
- if(lwdSet) {
- lwdo <- rep(lwdo,length.out=2)
- dots$lwd <- lwdo[1]
- }
- if(!is.null(dots$lty)) {ltySet <- TRUE; ltyo <- eval(dots$lty)}
- if(ltySet && ((!is.numeric(ltyo) && length(ltyo)==1)||
- is.numeric(ltyo))){
- ltyo <- list(ltyo,ltyo)
- dots$lty <- ltyo[[1]]
- }else{ if (ltySet && !is.numeric(ltyo) && length(ltyo)==2){
- dots$lty <- ltyo[[1]]
- }
- }
- if(is.null(dots$main)) dots$main <- gettext("Cniper point plot")
- if(is.null(dots$xlab)) dots$xlab <- gettext("Dirac point")
+ mcl$IC1 <- optIC(model = L2Fam, risk = asCov())
+ mcl$IC2 <- optIC(model = robMod, risk = risk)
+ mcl$L2Fam <- NULL
if(is.null(dots$ylab))
- dots$ylab <- gettext("Asymptotic Risk difference (classic - robust)")
- dots$type <- "l"
- do.call(plot, dots)
+ mcl$ylab <- gettext("Asymptotic Risk difference (classic - robust)")
+ if(is.null(dots$main))
+ mcl$main <- gettext("Cniper point plot")
- dots <- makedotsLowLevel(dots)
- dots$x <- dots$y <- NULL
- if(colSet) dots$col <- colo[2]
- if(lwdSet) dots$lwd <- lwdo[2]
- if(ltySet) dots$lty <- ltyo[[2]]
-
- dots$h <- if(scaleY) scaleY.fct(0) else 0
- do.call(abline, dots)
- .plotRescaledAxis(scaleX, scaleX.fct, scaleX.inv, scaleY,scaleY.fct,
- scaleY.inv, dots$xlim, dots$ylim, x, ypts = 400,
- n = scaleN, x.ticks = x.ticks, y.ticks = y.ticks)
- if(!is.null(data))
- return(.plotData(data, dots, mc, fun, L2Fam, eta))
- return(invisible(NULL))
+ if(withMaxRisk) mcl$fromCniperPlot <- TRUE
+ do.call(cniperCont, mcl)
}
-}
+
Modified: branches/robast-0.9/pkg/ROptEst/R/comparePlot.R
===================================================================
--- branches/robast-0.9/pkg/ROptEst/R/comparePlot.R 2013-01-14 15:23:03 UTC (rev 532)
+++ branches/robast-0.9/pkg/ROptEst/R/comparePlot.R 2013-01-15 00:31:10 UTC (rev 533)
@@ -1,21 +1,26 @@
setMethod("comparePlot", signature("IC","IC"),
function(obj1,obj2, obj3 = NULL, obj4 = NULL, data = NULL,
- ..., withSweave = getdistrOption("withSweave"),
- main = FALSE, inner = TRUE, sub = FALSE,
- col = par("col"), lwd = par("lwd"), lty,
- col.inner = par("col.main"), cex.inner = 0.8,
- bmar = par("mar")[1], tmar = par("mar")[3],
- with.legend = TRUE, legend.bg = "white",
+ ..., withSweave = getdistrOption("withSweave"),
+ main = FALSE, inner = TRUE, sub = FALSE,
+ col = par("col"), lwd = par("lwd"), lty,
+ col.inner = par("col.main"), cex.inner = 0.8,
+ bmar = par("mar")[1], tmar = par("mar")[3],
+ with.legend = FALSE, legend = NULL, legend.bg = "white",
legend.location = "bottomright", legend.cex = 0.8,
+ withMBR = FALSE, MBRB = NA, MBR.fac = 2, col.MBR = par("col"),
+ lty.MBR = "dashed", lwd.MBR = 0.8,
+ scaleX = FALSE, scaleX.fct, scaleX.inv,
+ scaleY = FALSE, scaleY.fct = pnorm, scaleY.inv=qnorm,
+ scaleN = 9, x.ticks = NULL, y.ticks = NULL,
mfColRow = TRUE, to.draw.arg = NULL,
cex.pts = 1, col.pts = par("col"),
pch.pts = 1, jitter.fac = 1, with.lab = FALSE,
- lab.pts = NULL, lab.font = NULL,
+ lab.pts = NULL, lab.font = NULL, alpha.trsp = NA,
which.lbs = NULL, which.Order = NULL, return.Order = FALSE){
- mcl <- match.call(call = sys.call(sys.parent(1)))
+ mcl <- match.call(call = sys.call(sys.parent(1)), expand.dots = TRUE)
- L2Fam <- eval(x at CallL2Fam); trafO <- trafo(L2Fam at param)
+ L2Fam <- eval(obj1 at CallL2Fam); trafO <- trafo(L2Fam at param)
dims <- nrow(trafO); to.draw <- 1:dims
if(! is.null(to.draw.arg)){
if(is.character(to.draw.arg))
@@ -27,15 +32,17 @@
MBRB <- matrix(rep(t(MBRB), length.out=dims0*2),ncol=2, byrow=T)
if(withMBR && all(is.na(MBRB))){
- robModel <- InfRobModel(center = L2fam, neighbor =
+ robModel <- InfRobModel(center = L2Fam, neighbor =
ContNeighborhood(radius = 0.5))
ICmbr <- try(optIC(model = robModel, risk = asBias()), silent=TRUE)
if(!is(ICmbr,"try-error"))
- MBRB <- .getExtremeCoordIC(ICmbr, distribution(L2Fam), todraw)
+ MBRB <- .getExtremeCoordIC(ICmbr, distribution(L2Fam), to.draw)
else withMBR <- FALSE
}
mcl$MBRB <- MBRB
mcl$withMBR <- withMBR
do.call(getMethod("comparePlot", signature("IC","IC"),
- where="RobAStBase"), mcl)
- })
+ where="RobAStBase"), as.list(mcl[-1]),
+ envir=parent.frame(2))
+ return(invisible())
+ })
Modified: branches/robast-0.9/pkg/ROptEst/R/getStartIC.R
===================================================================
--- branches/robast-0.9/pkg/ROptEst/R/getStartIC.R 2013-01-14 15:23:03 UTC (rev 532)
+++ branches/robast-0.9/pkg/ROptEst/R/getStartIC.R 2013-01-15 00:31:10 UTC (rev 533)
@@ -3,14 +3,20 @@
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
+ mc <- match.call(expand.dots=FALSE, call = sys.call(sys.parent(1)))
+ dots <- as.list(mc$"...")
+ if("fsCor" %in% names(dots)){
+ fsCor <- dots[["fsCor"]]
+ dots$fsCor <- NULL
+ }else fsCor <- 1
+ if("eps" %in% names(dots)){
+ eps <- dots[["eps"]]
+ dots$eps <- NULL
+ }else eps <- NULL
+ if("neighbor" %in% names(dots)){
+ neighbor <- dots[["neighbor"]]
+ dots$neighbor <- NULL
+ }else neighbor <- ContNeighborhood()
sm.rmx <- selectMethod("radiusMinimaxIC", signature(
class(model),class(neighbor),class(risk)))
@@ -31,13 +37,13 @@
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)
+ ICstart <- do.call(radiusMinimaxIC, args=arg.rmx, envir=parent.frame(2))
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)
+ ICstart <- do.call(optIC, args = arg.optic, envir=parent.frame(2))
}
}else{
neighbor at radius <- eps$sqn * fsCor * eps$e
@@ -48,7 +54,7 @@
if(..debug) print(c(arg.optic=arg.optic))
# print(arg.optic)
# print("----------------------------------------------------")
- ICstart <- do.call(optIC, args = arg.optic)
+ ICstart <- do.call(optIC, args = arg.optic, envir=parent.frame(2))
}
return(ICstart)
})
@@ -65,7 +71,8 @@
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)
+ sng <- try(getFromNamespace(.versionSuff(gridn), ns = "RobExtremes"),
+ silent=TRUE)
if(!is(sng,"try-error")) nsng <- names(sng)
if(length(nsng)){
if(nam %in% nsng){
@@ -84,6 +91,6 @@
}
mc$risk <- if(type(risk)==".MBRE") asMSE() else asBias()
mc$neighbor <- ContNeighborhood(radius=0.5)
- return(do.call(getStartIC, as.list(mc[-1])))
+ return(do.call(getStartIC, as.list(mc[-1]), envir=parent.frame(2)))
})
Modified: branches/robast-0.9/pkg/ROptEst/R/internal-getpsi.R
===================================================================
--- branches/robast-0.9/pkg/ROptEst/R/internal-getpsi.R 2013-01-14 15:23:03 UTC (rev 532)
+++ branches/robast-0.9/pkg/ROptEst/R/internal-getpsi.R 2013-01-15 00:31:10 UTC (rev 533)
@@ -1,14 +1,14 @@
.getPsi <- function(xi, beta, fct, L2Fam , type){
L2deriv <- L2Fam at L2deriv
- dbeta <- diag(c(beta,1))
+ .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)))
+ 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
+ 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
+ Aw <- .dbeta%*%matrix(c(fct[[10]](xi),am,am,fct[[13]](xi)),2,2)%*%.dbeta
Modified: branches/robast-0.9/pkg/ROptEst/R/roptest.R
===================================================================
--- branches/robast-0.9/pkg/ROptEst/R/roptest.R 2013-01-14 15:23:03 UTC (rev 532)
+++ branches/robast-0.9/pkg/ROptEst/R/roptest.R 2013-01-15 00:31:10 UTC (rev 533)
@@ -1,7 +1,7 @@
###############################################################################
## Optimally robust estimation
###############################################################################
-roptest <- function(x, L2Fam, eps, eps.lower, eps.upper, fsCor = 1, initial.est,
+roptest.old <- 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",
Modified: branches/robast-0.9/pkg/ROptEst/R/roptest.new.R
===================================================================
--- branches/robast-0.9/pkg/ROptEst/R/roptest.new.R 2013-01-14 15:23:03 UTC (rev 532)
+++ branches/robast-0.9/pkg/ROptEst/R/roptest.new.R 2013-01-15 00:31:10 UTC (rev 533)
@@ -2,7 +2,7 @@
## Optimally robust estimation
###############################################################################
-roptest.n <- function(x, L2Fam, eps, eps.lower, eps.upper, fsCor = 1, initial.est,
+roptest <- function(x, L2Fam, eps, eps.lower, eps.upper, fsCor = 1, initial.est,
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/robast -r 533
More information about the Robast-commits
mailing list