[Robast-commits] r528 - in branches/robast-0.9/pkg: ROptEst/R ROptEst/man RobAStBase RobAStBase/R RobAStBase/man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Jan 11 02:30:50 CET 2013
Author: ruckdeschel
Date: 2013-01-11 02:30:50 +0100 (Fri, 11 Jan 2013)
New Revision: 528
Added:
branches/robast-0.9/pkg/ROptEst/R/makedots.R
branches/robast-0.9/pkg/RobAStBase/R/getRiskBV.R
branches/robast-0.9/pkg/RobAStBase/R/plotRescaledAxis.R
branches/robast-0.9/pkg/RobAStBase/R/selectorder.R
branches/robast-0.9/pkg/RobAStBase/man/getRiskFctBV-methods.Rd
Modified:
branches/robast-0.9/pkg/ROptEst/R/cniperCont.R
branches/robast-0.9/pkg/ROptEst/man/cniperCont.Rd
branches/robast-0.9/pkg/RobAStBase/NAMESPACE
branches/robast-0.9/pkg/RobAStBase/R/AllGeneric.R
branches/robast-0.9/pkg/RobAStBase/R/AllPlot.R
branches/robast-0.9/pkg/RobAStBase/R/ddPlot_utils.R
branches/robast-0.9/pkg/RobAStBase/R/infoPlot.R
branches/robast-0.9/pkg/RobAStBase/man/infoPlot.Rd
branches/robast-0.9/pkg/RobAStBase/man/plot-methods.Rd
Log:
modularized diagnostic plots; still checking + debugging necessary
Modified: branches/robast-0.9/pkg/ROptEst/R/cniperCont.R
===================================================================
--- branches/robast-0.9/pkg/ROptEst/R/cniperCont.R 2013-01-09 13:01:17 UTC (rev 527)
+++ branches/robast-0.9/pkg/ROptEst/R/cniperCont.R 2013-01-11 01:30:50 UTC (rev 528)
@@ -1,4 +1,5 @@
-setMethod("cniperCont", signature(IC1 = "IC",
+PETER <- FALSE
+setMethod("cniperCont", signature(IC1 = "IC",
IC2 = "IC",
L2Fam = "L2ParamFamily",
neighbor = "ContNeighborhood",
@@ -28,7 +29,8 @@
abline(h = 0)
invisible()
})
-setMethod("cniperPoint", signature(L2Fam = "L2ParamFamily",
+
+setMethod("cniperPoint", signature(L2Fam = "L2ParamFamily",
neighbor = "ContNeighborhood",
risk = "asMSE"),
function(L2Fam, neighbor, risk, lower, upper){
@@ -40,18 +42,22 @@
maxMSE <- Risks(eta)$asMSE$value
Delta <- sqrt(maxMSE - tr.invF)/neighbor at radius
fun <- function(x){
- y <- evalIC(psi, 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",
+
+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)),
+ function(L2Fam, neighbor, risk, lower, upper, n = 101, ...,
+ scaleX = FALSE, scaleX.fct, scaleX.inv,
+ scaleY = FALSE, scaleY.fct = pnorm, scaleY.inv=qnorm,
+ ){
+ 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)))
@@ -60,7 +66,7 @@
eta <- optIC(model = robMod, risk = asMSE())
maxMSE <- Risks(eta)$asMSE$value
fun <- function(x){
- y <- evalIC(psi, 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)
@@ -68,7 +74,7 @@
colSet <- ltySet <- lwdSet <- FALSE
if(!is.null(dots$col)) {colSet <- TRUE; colo <- eval(dots$col)}
if(colSet) {
- colo <- rep(colo,length.out=2)
+ colo <- rep(colo,length.out=2)
dots$col <- colo[1]
}
if(!is.null(dots$lwd)) {lwdSet <- TRUE; lwdo <- eval(dots$lwd)}
@@ -78,7 +84,7 @@
}
if(!is.null(dots$lty)) {ltySet <- TRUE; ltyo <- eval(dots$lty)}
if(ltySet && ((!is.numeric(ltyo) && length(ltyo)==1)||
- is.numeric(ltyo))){
+ is.numeric(ltyo))){
ltyo <- list(ltyo,ltyo)
dots$lty <- ltyo[[1]]
}else{ if (ltySet && !is.numeric(ltyo) && length(ltyo)==2){
@@ -87,9 +93,13 @@
}
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))
+ 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)
@@ -101,6 +111,284 @@
do.call(abline, dots)
invisible()
})
-#
-#cniperPointPlot(L2Fam=N0, neighbor=ContNeighborhood(radius = 0.5), risk=asMSE(),lower=-12, n =30, upper=8, lwd=c(2,4),lty=list(c(5,1),3),col=c(2,4))
-
\ No newline at end of file
+
+if(PETER){
+.plotData <- function(data, dots, origCl, fun, L2Fam, IC){
+ dotsP <- .makedotsP(dots)
+ dotsP$col <- rep(origCl$col.pts, length.out=n)
+ dotsP$pch <- rep(origCl$pch.pts, length.out=n)
+
+ 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)
+ 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, origCl$scaleY.inv,
+ dots$xlim, dots$ylim, dots)
+
+ dotsP$x <- resc.dat$X
+ dotsP$y <- resc.dat$Y
+
+ 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))}
+
+ 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
+ absy <- absInfoEval(x.d, absy.f)
+
+ dotsP$cex <- log(absy+1)*3*rep(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)
+ return(invisible(NULL))
+ }
+
+cniperContPlot <- function(IC1, IC2, data = NULL, ...,
+ neighbor, risk, lower=getdistrOption("DistrResolution"),
+ upper=1-getdistrOption("DistrResolution"), n = 101,
+ scaleX = FALSE, scaleX.fct, scaleX.inv,
+ scaleY = FALSE, scaleY.fct = pnorm, scaleY.inv=qnorm,
+ cex.pts = 1, col.pts = par("col"),
+ pch.pts = 1, jitter.fac = 1, with.lab = FALSE,
+ lab.pts = NULL, lab.font = NULL,
+ which.lbs = NULL, which.Order = NULL,
+ return.Order = FALSE){
+
+ mc <- match.call(call = sys.call(sys.parent(1)),
+ expand.dots = FALSE)
+ dots <- 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))
+
+
+ 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))
+ }
+ x <- q(L2Fam)(seq(lower,upper,length=n))
+ resc <- RobAStBase:::.rescalefct(x, function(u) sapply(u,fun),
+ scaleX, scaleX.fct, scaleX.inv,
+ scaleY, scaleY.fct, scaleY.inv,
+ dots$xlim, dots$ylim, dots)
+ 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)"
+
+ 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]]
+ }
+ }
+
+ do.call(plot,dots)
+
+ 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)
+
+ RobAStBase:::.plotRescaledAxis(scaleX, scaleX.fct, scaleX.inv,
+ scaleY,scaleY.fct, scaleY.inv,
+ dots$xlim, dots$ylim, x, ypts = 400)
+ if(!is.null(data))
+ return(.plotData(data, dots, mc, fun, L2Fam, IC1))
+ invisible(NULL)
+}
+
+cniperPoint <- function(L2Fam, neighbor, risk= asMSE(),
+ 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
+
+ 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))
+ }
+
+ res <- uniroot(fun, lower = lower, upper = upper)$root
+ names(res) <- "cniper point"
+ res
+}
+
+cniperPointPlot <- function(L2Fam, data=NULL, ..., neighbor, risk= asMSE(),
+ lower=getdistrOption("DistrResolution"),
+ upper=1-getdistrOption("DistrResolution"), n = 101,
+ scaleX = FALSE, scaleX.fct, scaleX.inv,
+ scaleY = FALSE, scaleY.fct = pnorm, scaleY.inv=qnorm,
+ cex.pts = 1, col.pts = par("col"),
+ pch.pts = 1, jitter.fac = 1, with.lab = FALSE,
+ lab.pts = NULL, lab.font = NULL,
+ which.lbs = NULL, which.Order = NULL,
+ return.Order = FALSE){
+
+ mc <- match.call(call = sys.call(sys.parent(1)),
+ expand.dots = FALSE)
+ dots <- 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 <- RobAStBase:::.rescalefct(x, function(u) sapply(u,fun),
+ scaleX, scaleX.fct, scaleX.inv,
+ scaleY, scaleY.fct, scaleY.inv,
+ 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")
+ if(is.null(dots$ylab))
+ dots$ylab <- gettext("Asymptotic Risk difference (classic - robust)")
+ dots$type <- "l"
+ do.call(plot, dots)
+
+ 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)
+ RobAStBase:::.plotRescaledAxis(scaleX, scaleX.fct, scaleX.inv,
+ scaleY,scaleY.fct, scaleY.inv,
+ dots$xlim, dots$ylim, x, ypts = 400)
+
+ if(!is.null(data))
+ return(.plotData(data, dots, mc, fun, L2Fam, eta))
+ return(invisible(NULL))
+}
+
+}
+
Added: branches/robast-0.9/pkg/ROptEst/R/makedots.R
===================================================================
--- branches/robast-0.9/pkg/ROptEst/R/makedots.R (rev 0)
+++ branches/robast-0.9/pkg/ROptEst/R/makedots.R 2013-01-11 01:30:50 UTC (rev 528)
@@ -0,0 +1,41 @@
+.makedotsLowLevel <- function(dots){
+ dots$sub <- dots$xlab <- dots$ylab <- dots$main <- dots$type <- NULL
+ dots$xlim <- dots$ylim <- dots$yaxt <- dots$axes <- dots$xaxt <- NULL
+ dots$panel.last <- dots$panel.first <- dots$frame.plot <- dots$ann <-NULL
+ dots$log <- dots$asp <- NULL
+ return(dots)
+}
+.deleteDotsABLINE <- function(dots){
+ dots$reg <- dots$a <- dots$b <- NULL
+ dots$untf <- dots$h <- dots$v <- NULL
+ dots
+}
+.deleteDotsTEXT <- function(dots){
+ dots$labels <- dots$offset <- dots$vfont <- dots$pos <- dots$font <- NULL
+ dots
+}
+.makedotsL <- function(dots){
+ dots <- .makedotsLowLevel(dots)
+ dots$pch <- dots$cex <- NULL
+ .deleteDotsABLINE(.deleteDotsTEXT(dots))
+}
+.makedotsP <- function(dots){
+ dots <- .makedotsLowLevel(dots)
+ dots$lwd <- NULL
+ .deleteDotsABLINE(.deleteDotsTEXT(dots))
+}
+.makedotsPt <- function(dots){
+ dots <- dots[names(dots) %in% c("bg", "lwd", "lty")]
+ if (length(dots) == 0 ) dots <- NULL
+ return(dots)
+}
+.makedotsAB <- function(dots){
+ dots <- .makedotsLowLevel(dots)
+ dots <- .deleteDotsTEXT(dots)
+ dots$pch <- dots$cex <- NULL
+}
+.makedotsT <- function(dots){
+ dots <- .makedotsLowLevel(dots)
+ dots <- .deleteDotsABLINE(dots)
+ dots
+}
Modified: branches/robast-0.9/pkg/ROptEst/man/cniperCont.Rd
===================================================================
--- branches/robast-0.9/pkg/ROptEst/man/cniperCont.Rd 2013-01-09 13:01:17 UTC (rev 527)
+++ branches/robast-0.9/pkg/ROptEst/man/cniperCont.Rd 2013-01-11 01:30:50 UTC (rev 528)
@@ -17,6 +17,29 @@
other procedure.
}
\usage{
+cniperContPlot(IC1, IC2, data = NULL, ...,
+ neighbor, risk, lower = getdistrOption("DistrResolution"),
+ upper = 1-getdistrOption("DistrResolution"), n = 101,
+ scaleX = FALSE, scaleX.fct, scaleX.inv,
+ scaleY = FALSE, scaleY.fct = pnorm, scaleY.inv=qnorm,
+ cex.pts = 1, col.pts = par("col"),
+ pch.pts = 1, jitter.fac = 1, with.lab = FALSE,
+ lab.pts = NULL, lab.font = NULL,
+ which.lbs = NULL, which.Order = NULL,
+ return.Order = FALSE)
+
+cniperPointPlot(L2Fam, data=NULL, ..., neighbor, risk= asMSE(),
+ lower=getdistrOption("DistrResolution"),
+ upper=1-getdistrOption("DistrResolution"), n = 101,
+ scaleX = FALSE, scaleX.fct, scaleX.inv,
+ scaleY = FALSE, scaleY.fct = pnorm, scaleY.inv=qnorm,
+ cex.pts = 1, col.pts = par("col"),
+ pch.pts = 1, jitter.fac = 1, with.lab = FALSE,
+ lab.pts = NULL, lab.font = NULL,
+ which.lbs = NULL, which.Order = NULL,
+ return.Order = FALSE)
+
+
cniperCont(IC1, IC2, L2Fam, neighbor, risk, ...)
\S4method{cniperCont}{IC,IC,L2ParamFamily,ContNeighborhood,asMSE}(IC1,
IC2, L2Fam, neighbor, risk, lower, upper, n = 101)
@@ -37,8 +60,49 @@
\item{risk}{ object of class \code{RiskType} }
\item{\dots}{ additional parameters (in particular to be passed on to \code{plot}). }
\item{lower, upper}{ the lower and upper end points of the
- contamination interval. }
+ contamination interval (in prob-scale). }
\item{n}{ number of points between \code{lower} and \code{upper}}
+ \item{scaleX}{logical; shall X-axis be rescaled (by default according to the cdf of
+ the underlying distribution)?}
+ \item{scaleY}{logical; shall Y-axis be rescaled (by default according to a probit scale)?}
+ \item{scaleX.fct}{an isotone, vectorized function mapping the domain of the IC(s)
+ to [0,1]; if \code{scaleX} is \code{TRUE} and \code{scaleX.fct} is
+ missing, the cdf of the underlying observation distribution.}
+ \item{scaleX.inv}{the inverse function to \code{scale.fct}, i.e., an isotone,
+ vectorized function mapping [0,1] to the domain of the IC(s)
+ such that for any \code{x} in the domain,
+ \code{scaleX.inv(scaleX.fct(x))==x}; if \code{scaleX} is \code{TRUE}
+ and \code{scaleX.inv} is
+ missing, the quantile function of the underlying observation distribution.}
+ \item{scaleY.fct}{an isotone, vectorized function mapping for each coordinate the
+ range of the respective coordinate of the IC(s)
+ to [0,1]; defaulting to the cdf of \eqn{{\cal N}(0,1)}{N(0,1)}.}
+ \item{scaleY.inv}{an isotone, vectorized function mapping for each coordinate
+ the range [0,1] into the range of the respective coordinate of the IC(s);
+ defaulting to the quantile function of \eqn{{\cal N}(0,1)}{N(0,1)}.}
+ \item{cex.pts}{size of the points of the second argument plotted}
+ \item{col.pts}{color of the points of the second argument plotted}
+ \item{pch.pts}{symbol of the points of the second argument plotted}
+ \item{with.lab}{logical; shall labels be plotted to the observations?}
+ \item{lab.pts}{character or NULL; labels to be plotted to the observations; if \code{NULL}
+ observation indices;}
+ \item{lab.font}{font to be used for labels}
+ \item{jitter.fac}{jittering factor used in case of a \code{DiscreteDistribution}
+ for plotting points of the second argument in a jittered fashion.}
+ \item{which.lbs}{either an integer vector with the indices of the observations
+ to be plotted into graph or \code{NULL} --- then no observation is excluded}
+ \item{which.Order}{we order the observations (descending) according to the norm given by
+ \code{normtype(object)}; then \code{which.Order}
+ either is an integer vector with the indices of the \emph{ordered}
+ observations (remaining after a possible reduction by argument \code{which.lbs})
+ to be plotted into graph or \code{NULL} --- then no (further) observation
+ is excluded.}
+ \item{return.Order}{logical; if \code{TRUE}, an order vector
+ is returned; more specifically, the order of the (remaining) observations
+ given by their original index is returned (remaining means: after a possible
+ reduction by argument \code{which.lbs}, and ordering is according to the norm given by
+ \code{normtype(object)});
+ otherwise we return \code{invisible()} as usual.}
}
\details{
In case of \code{cniperCont} the difference between the risks of two ICs
Modified: branches/robast-0.9/pkg/RobAStBase/NAMESPACE
===================================================================
--- branches/robast-0.9/pkg/RobAStBase/NAMESPACE 2013-01-09 13:01:17 UTC (rev 527)
+++ branches/robast-0.9/pkg/RobAStBase/NAMESPACE 2013-01-11 01:30:50 UTC (rev 528)
@@ -62,6 +62,7 @@
exportMethods("ddPlot", "qqplot")
exportMethods("cutoff.quantile", "cutoff.quantile<-")
exportMethods("samplesize<-", "samplesize")
+exportMethods("getRiskFctBV")
export("oneStepEstimator", "kStepEstimator")
export("ContNeighborhood", "TotalVarNeighborhood")
export("FixRobModel", "InfRobModel")
@@ -71,3 +72,4 @@
export("cutoff","cutoff.chisq","cutoff.sememp")
export("outlyingPlotIC", "RobAStBaseMASK")
export("OMSRRisk","MBRRisk","RMXRRisk")
+export("getRiskFctBV")
\ No newline at end of file
Modified: branches/robast-0.9/pkg/RobAStBase/R/AllGeneric.R
===================================================================
--- branches/robast-0.9/pkg/RobAStBase/R/AllGeneric.R 2013-01-09 13:01:17 UTC (rev 527)
+++ branches/robast-0.9/pkg/RobAStBase/R/AllGeneric.R 2013-01-11 01:30:50 UTC (rev 528)
@@ -225,3 +225,7 @@
setGeneric("samplesize<-",
function(object, value) standardGeneric("samplesize<-"))
}
+if(!isGeneric("getRiskFctBV")){
+ setGeneric("getRiskFctBV", function(risk, biastype) standardGeneric("getRiskFctBV"))
+}
+
Modified: branches/robast-0.9/pkg/RobAStBase/R/AllPlot.R
===================================================================
--- branches/robast-0.9/pkg/RobAStBase/R/AllPlot.R 2013-01-09 13:01:17 UTC (rev 527)
+++ branches/robast-0.9/pkg/RobAStBase/R/AllPlot.R 2013-01-11 01:30:50 UTC (rev 528)
@@ -1,15 +1,21 @@
setMethod("plot", signature(x = "IC", y = "missing"),
- function(x,...,withSweave = getdistrOption("withSweave"),
+ function(x, ...,withSweave = getdistrOption("withSweave"),
main = FALSE, inner = TRUE, sub = FALSE,
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,
mfColRow = TRUE, to.draw.arg = NULL){
xc <- match.call(call = sys.call(sys.parent(1)))$x
dots <- match.call(call = sys.call(sys.parent(1)),
expand.dots = FALSE)$"..."
+ dotsLeg <- dotsT <- dotsL <- .makedotsLowLevel(dots)
-
if(!is.logical(inner)){
if(!is.list(inner))
inner <- as.list(inner)
@@ -21,6 +27,10 @@
L2Fam <- eval(x at CallL2Fam)
+ if(missing(scaleX.fct)){
+ scaleX.fct <- p(L2Fam)
+ scaleX.inv <- q(L2Fam)
+ }
trafO <- trafo(L2Fam at param)
dims <- nrow(trafO)
@@ -39,6 +49,16 @@
nrows <- trunc(sqrt(dims0))
ncols <- ceiling(dims0/nrows)
+ MBRB <- matrix(rep(t(MBRB), length.out=dims0*2),ncol=2, byrow=T)
+ if(withMBR && all(is.na(MBRB))){
+ 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)
+ else withMBR <- FALSE
+ }
+ MBRB <- MBRB * MBR.fac
e1 <- L2Fam at distribution
if(!is(e1, "UnivariateDistribution")) stop("not yet implemented")
@@ -169,7 +189,20 @@
}
}
+ if(with.legend){
+ if(missing(legend.location)){
+ legend.location <- distr:::.fillList(list("bottomright"), dims0)
+ }else{
+ legend.location <- as.list(legend.location)
+ legend.location <- distr:::.fillList(legend.location, dims0)
+ }
+ if(is.null(legend)){
+ legend <- vector("list",dims0)
+ legend <- distr:::.fillList(as.list(xc),dims0)
+ }
+ }
+
w0 <- getOption("warn")
options(warn = -1)
on.exit(options(warn = w0))
@@ -188,28 +221,66 @@
do.call(par,args=parArgs)
- dotsT <- dots
- dotsT["main"] <- NULL
- dotsT["cex.main"] <- NULL
- dotsT["col.main"] <- NULL
- dotsT["line"] <- NULL
+ dotsT["pch"] <- dotsT["cex"] <- NULL
+ dotsT["col"] <- dotsT["lwd"] <- NULL
+ dotsL["cex"] <- dotsLeg["bg"] <- dotsLeg["cex"] <- NULL
+ dots$ylim <- NULL
- dots$ylim <- NULL
for(i in 1:dims0){
indi <- to.draw[i]
if(!is.null(ylim)) dots$ylim <- ylim[,i]
- do.call(plot, args=c(list(x.vec, sapply(x.vec, IC1 at Map[[indi]]),
+ resc <-.rescalefct(x.vec, IC1 at Map[[indi]], scaleX, scaleX.fct,
+ scaleX.inv, scaleY, scaleY.fct, scaleY.inv,
+ xlim[,i], ylim[,i], dots)
+ dots <- resc$dots
+ x.vec1 <- resc$X
+ y.vec1 <- resc$Y
+ do.call(plot, args=c(list(x.vec1, y.vec1,
type = plty, lty = lty,
xlab = "x", ylab = "(partial) IC"),
dots))
+
+ .plotRescaledAxis(scaleX, scaleX.fct, scaleX.inv,
+ scaleY,scaleY.fct, scaleY.inv,
+ xlim[,i], ylim[,i], x.vec1, ypts = 400)
+ if(withMBR){
+ MBR.i <- MBRB[i,]
+ if(scaleY) MBR.i <- scaleY.fct(MBR.i)
+ abline(h=MBR.i, col=col.MBR, lty=lty.MBR, lwd = lwd.MBR)
+ }
if(is(e1, "DiscreteDistribution")){
- x.vec1 <- seq(from = min(x.vec), to = max(x.vec), length = 1000)
- do.call(lines,args=c(list(x.vec1, sapply(x.vec1, IC1 at Map[[indi]]),
- lty = "dotted"), dots))
+ x.vec1a <- seq(from = min(x.vec), to = max(x.vec), length = 1000)
+ if(scaleX){
+ if(!is.null(xlim)){
+ dots$xlim <- scaleX.fct(xlim[,i])
+ x.vec10 <- x.vec1a[xvec1a>=xlim[1,i] & xvec1a<=xlim[2,i]]
+ }
+ x.vec1 <- scaleX.fct(x.vec10)
+ x.vec1 <- distr:::.DistrCollapse(x.vec1, 0*x.vec1+1/length(x.vec1))
+ dots$axes <- NULL
+ dots$xaxt <- "n"
+ }
+ y.vec1 <- sapply(x.vec1, IC1 at Map[[indi]])
+ if(scaleY){
+ y.vec1 <- scaleY.fct(y.vec)
+ if(!is.null(ylim)) dots$ylim <- scaleY.fct(ylim[,i])
+ dots$axes <- NULL
+ dots$yaxt <- "n"
+ }
+
+ dotsL$lty <- NULL
+ do.call(lines,args=c(list(x.vec1, y.vec1,
+ lty = "dotted"), dotsL))
+
}
do.call(title,args=c(list(main = innerT[indi]), dotsT, line = lineT,
cex.main = cex.inner, col.main = col.inner))
+ if(with.legend)
+ legend(.legendCoord(legend.location[[i]], scaleX, scaleX.fct,
+ scaleY, scaleY.fct), bg = legend.bg,
+ legend = legend[[i]], dotsLeg, cex = legend.cex*fac.leg)
+
}
if(!hasArg(cex.main)) cex.main <- par("cex.main") else cex.main <- dots$"cex.main"
if(!hasArg(col.main)) col.main <- par("col.main") else col.main <- dots$"col.main"
@@ -236,16 +307,9 @@
expand.dots = FALSE)$"..."
n <- if(!is.null(dim(y))) nrow(y) else length(y)
- oN0 <- NULL
- if(is.null(which.lbs))
- which.lbs <- 1:n
- which.lbs0 <- (1:n) %in% which.lbs
- which.lbx <- rep(which.lbs0, length.out=length(y))
- y0 <- y[which.lbx]
- n <- if(!is.null(dim(y0))) nrow(y0) else length(y0)
- oN <- (1:n)[which.lbs0]
+ pch.pts <- rep(pch.pts, length.out=n)
+ lab.pts <- if(is.null(lab.pts)) paste(1:n) else rep(lab.pts,n)
-
L2Fam <- eval(x at CallL2Fam)
trafO <- trafo(L2Fam at param)
dims <- nrow(trafO)
@@ -260,22 +324,13 @@
ICMap <- IC1 at Map
absInfo <- sapply(y, absInfo at Map[[1]])
- absInfo0 <- absInfo[which.lbs]/max(absInfo)
- if (n==length(y0)) {
- oN <- order(absInfo0)
- oN0 <- order(absInfo)
- oN0 <- oN0[oN0 %in% which.lbs]
- y0 <- y0[oN]
- if(!is.null(which.Order)){
- oN <- oN0[(n+1)-which.Order]
- y0 <- y[oN]
- absInfo0 <- absInfo[oN]/max(absInfo[oN])
- }
- }
- if(is.null(lab.pts)) lab.pts <- paste(oN)
- else {lab.pts <- rep(lab.pts, length.out=length(y))
- lab.pts <- lab.pts[oN]}
+ sel <- .SelectOrderData(y, function(x)sapply(x, absInfo at Map[[1]]),
+ which.lbs, which.Order)
+ i.d <- sel$ind
+ i0.d <- sel$ind1
+ x.d <- sel$data
+ n <- length(i.d)
dots.without <- dots
dots.without$col <- dots.without$cex <- dots.without$pch <- NULL
@@ -286,23 +341,41 @@
dots$panel.last <- NULL
pL <- substitute({
+ y1 <- y0s
ICy <- sapply(y0s,ICMap0[[indi]])
+ resc.dat <-.rescalefct(y0s, function(x) sapply(x,ICMap0[[indi]]),
+ scaleX, scaleX.fct, scaleX.inv,
+ scaleY, scaleY.fct, scaleY.inv,
+ dwo0$xlim, dwo0$ylim, dwo0)
+ y1 <- resc.dat$X
+ ICy <- resc.dat$Y
+
if(is(e1, "DiscreteDistribution"))
ICy <- jitter(ICy, factor = jitter.fac0)
- do.call(points, args=c(list(y0s, ICy, cex = log(absy0+1)*3*cex0,
+
+
+ do.call(points, args=c(list(y1, ICy, cex = log(absy0+1)*3*cex0,
col = col0, pch = pch0), dwo0))
if(with.lab0){
text(x = y0s, y = ICy, labels = lab.pts0,
cex = log(absy0+1)*1.5*cex0, col = col0)
}
pL0
- }, list(pL0 = pL, ICMap0 = ICMap, y0s = y0, absy0 = absInfo0,
- dwo0 = dots.without, cex0 = cex.pts, pch0 = pch.pts,
- col0 = col.pts, with.lab0 = with.lab, lab.pts0 = lab.pts,
+ }, list(pL0 = pL, ICMap0 = ICMap, y0s = x.d, absy0 = absInfo0,
+ dwo0 = dots.without, cex0 = cex.pts, pch0 = pch.pts[i.d],
+ col0 = col.pts, with.lab0 = with.lab, lab.pts0 = lab.pts[i.d],
jitter.fac0 = jitter.fac
))
do.call("plot", args = c(list(x = x, panel.last = pL), dots))
- if(return.Order) return(oN0)
+ if(return.Order) return(i0.d)
invisible()
})
+
+.getExtremeCoordIC <- function(IC, D, indi, n = 50000){
+ 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)))
+}
\ No newline at end of file
Modified: branches/robast-0.9/pkg/RobAStBase/R/ddPlot_utils.R
===================================================================
--- branches/robast-0.9/pkg/RobAStBase/R/ddPlot_utils.R 2013-01-09 13:01:17 UTC (rev 527)
+++ branches/robast-0.9/pkg/RobAStBase/R/ddPlot_utils.R 2013-01-11 01:30:50 UTC (rev 528)
@@ -98,43 +98,43 @@
if(is.null(dots$lty)) dots$lty <- par("lty")
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/robast -r 528
More information about the Robast-commits
mailing list