[Robast-commits] r894 - in pkg/ROptEst: . R inst inst/scripts man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Sep 4 16:31:17 CEST 2016
Author: ruckdeschel
Date: 2016-09-04 16:31:17 +0200 (Sun, 04 Sep 2016)
New Revision: 894
Added:
pkg/ROptEst/R/internalutilsFromRobAStBase.R
Removed:
pkg/ROptEst/R/getFiRisk.R
Modified:
pkg/ROptEst/DESCRIPTION
pkg/ROptEst/NAMESPACE
pkg/ROptEst/R/00internal.R
pkg/ROptEst/R/AllPlot.R
pkg/ROptEst/R/cniperCont.R
pkg/ROptEst/R/comparePlot.R
pkg/ROptEst/R/getAsRisk.R
pkg/ROptEst/R/getInfRobIC_asAnscombe.R
pkg/ROptEst/R/getInfRobIC_asBias.R
pkg/ROptEst/R/getInfRobIC_asCov.R
pkg/ROptEst/R/getMaxIneff.R
pkg/ROptEst/R/getModifyIC.R
pkg/ROptEst/R/getReq.R
pkg/ROptEst/R/interpolLM.R
pkg/ROptEst/R/leastFavorableRadius.R
pkg/ROptEst/R/optIC.R
pkg/ROptEst/R/plotWrapper.R
pkg/ROptEst/R/roptest.R
pkg/ROptEst/R/roptest.new.R
pkg/ROptEst/R/updateNorm.R
pkg/ROptEst/inst/NEWS
pkg/ROptEst/inst/scripts/BinomialModel.R
pkg/ROptEst/inst/scripts/MBRE.R
pkg/ROptEst/man/0ROptEst-package.Rd
pkg/ROptEst/man/CniperPointPlotWrapper.Rd
pkg/ROptEst/man/cniperCont.Rd
pkg/ROptEst/man/getReq.Rd
pkg/ROptEst/man/internal-interpolate.Rd
pkg/ROptEst/man/internal_Cniperplots.Rd
pkg/ROptEst/man/minmaxBias.Rd
pkg/ROptEst/man/optIC.Rd
Log:
prepared ROptEst (trunk) for submission to CRAN
Modified: pkg/ROptEst/DESCRIPTION
===================================================================
--- pkg/ROptEst/DESCRIPTION 2016-09-04 14:30:41 UTC (rev 893)
+++ pkg/ROptEst/DESCRIPTION 2016-09-04 14:31:17 UTC (rev 894)
@@ -1,15 +1,18 @@
Package: ROptEst
-Version: 0.9.1
-Date: 2013-09-11
+Version: 1.0
+Date: 2016-09-04
Title: Optimally robust estimation
-Description: Optimally robust estimation in general smoothly parameterized models using S4 classes and methods.
-Depends: R(>= 2.14.0), methods, distr(>= 2.5.2), distrEx(>= 2.4), distrMod(>= 2.5.2), RandVar(>= 0.9.2),
- RobAStBase(>= 0.9)
+Description: Optimally robust estimation in general smoothly parameterized models using S4
+ classes and methods.
+Depends: R(>= 2.14.0), methods, distr(>= 2.5.2), distrEx(>= 2.4), distrMod(>= 2.5.2),
+ RandVar(>= 0.9.2), RobAStBase(>= 1.0)
+Imports: startupmsg
Suggests: RobLox, MASS
-Author: Matthias Kohl, Peter Ruckdeschel
-Maintainer: Matthias Kohl <Matthias.Kohl at stamats.de>
+Authors at R: c(person("Matthias", "Kohl", role=c("cre", "cph"), email="Matthias.Kohl at stamats.de"),
+ person("Mykhailo", "Pupashenko", role="ctb", comment="contributed wrapper functions for diagnostic plots"),
+ person("Gerald", "Kroisandt", role="ctb", comment="contributed testing routines"),
+ person("Peter", "Ruckdeschel", role=c("aut", "cph")))
ByteCompile: yes
-LazyLoad: yes
License: LGPL-3
URL: http://robast.r-forge.r-project.org/
Encoding: latin1
Modified: pkg/ROptEst/NAMESPACE
===================================================================
--- pkg/ROptEst/NAMESPACE 2016-09-04 14:30:41 UTC (rev 893)
+++ pkg/ROptEst/NAMESPACE 2016-09-04 14:31:17 UTC (rev 894)
@@ -4,13 +4,20 @@
import("RandVar")
import("distrMod")
import("RobAStBase")
+importFrom("startupmsg", "buildStartupMessage", "infoShow")
+importFrom("graphics", "abline", "axis", "box", "par",
+ "points", "text", "matlines", "matpoints", "mtext",
+ "title")
+importFrom("stats", "complete.cases", "dnorm", "na.omit", "optim",
+ "optimize", "pnorm", "qnorm", "uniroot")
+importFrom("utils", "read.csv", "read.table", "str", "write.table")
+importFrom("graphics", "abline")
exportClasses("asAnscombe", "asL1", "asL4")
exportMethods("optIC",
"getInfRobIC",
"getFixRobIC",
"getAsRisk",
- "getFiRisk",
"getInfRad",
"getInfClip",
"getFixClip",
Modified: pkg/ROptEst/R/00internal.R
===================================================================
--- pkg/ROptEst/R/00internal.R 2016-09-04 14:30:41 UTC (rev 893)
+++ pkg/ROptEst/R/00internal.R 2016-09-04 14:31:17 UTC (rev 894)
@@ -1,12 +1,3 @@
-.rescalefct <- RobAStBase:::.rescalefct
-.plotRescaledAxis <- RobAStBase:::.plotRescaledAxis
-.makedotsP <- RobAStBase:::.makedotsP
-.makedotsLowLevel <- RobAStBase:::.makedotsLowLevel
-.SelectOrderData <- RobAStBase:::.SelectOrderData
-### helper function to recursively evaluate list
-.evalListRec <- RobAStBase:::.evalListRec
-
-
if(packageVersion("distrMod")<"2.5"){
.isUnitMatrix <- function(m){
### checks whether m is unit matrix
Modified: pkg/ROptEst/R/AllPlot.R
===================================================================
--- pkg/ROptEst/R/AllPlot.R 2016-09-04 14:30:41 UTC (rev 893)
+++ pkg/ROptEst/R/AllPlot.R 2016-09-04 14:31:17 UTC (rev 894)
@@ -7,10 +7,10 @@
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, n.MBR = 10000,
- scaleX = FALSE, scaleX.fct, scaleX.inv,
+ x.vec = NULL, 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){
+ mfColRow = TRUE, to.draw.arg = NULL, withSubst = TRUE){
mcl <- match.call(call = sys.call(sys.parent(1)), expand.dots = TRUE)
@@ -47,6 +47,6 @@
.getExtremeCoordIC <- function(IC, D, indi, n = 10000){
x <- q(D)(seq(1/2/n,1-1/2/n, length=n))
- y <- (matrix(evalIC(IC,matrix(x,ncol=1)),ncol=n))[indi,]
+ y <- (matrix(evalIC(IC,matrix(x,ncol=1)),ncol=n))[indi,,drop=FALSE]
return(cbind(min=apply(y,1,min),max=apply(y,1,max)))
}
\ No newline at end of file
Modified: pkg/ROptEst/R/cniperCont.R
===================================================================
--- pkg/ROptEst/R/cniperCont.R 2016-09-04 14:30:41 UTC (rev 893)
+++ pkg/ROptEst/R/cniperCont.R 2016-09-04 14:31:17 UTC (rev 894)
@@ -1,243 +1,358 @@
-
-.plotData <- function(
- ## helper function for cniper-type plots to plot in data
- data, # data to be plot in
- dots, # dots from the calling function
- origCl, # call from the calling function
- fun, # function to determine risk difference
- L2Fam, # L2Family
- IC # IC1 in cniperContPlot and eta in cniperPointPlot
-){
- dotsP <- .makedotsP(dots)
- 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),
- eval(origCl$which.lbs),
- eval(origCl$which.Order))
- i.d <- sel$ind
- i0.d <- sel$ind1
- y.d <- sel$y
- x.d <- sel$data
- n <- length(i.d)
-
- resc.dat <- .rescalefct(x.d, function(x) sapply(x,fun),
- 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
- dotsP$y <- resc.dat$Y
-
- trafo <- trafo(L2Fam at param)
- dims <- nrow(trafo)
- QF <- diag(dims)
- 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]])
- IC.rv <- as(diag(dims) %*% IC at Curve, "EuclRandVariable")
- absy.f <- t(IC.rv) %*% QF %*% IC.rv
- absy <- absInfoEval(x.d, absy.f)
-
- 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(!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))
- }
-
-
-.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,
- 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 = 19, jitter.fac = 1, with.lab = FALSE,
- lab.pts = NULL, lab.font = NULL, alpha.trsp = NA,
- which.lbs = NULL, which.Order = NULL,
- return.Order = FALSE){
-
- 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)
-
- 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)
- }
-
- 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))
- 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)
- dots$x <- resc$X
- dots$y <- resc$Y
-
-
- dots$type <- "l"
- 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)}
- 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)
-
- .plotRescaledAxis(scaleX, scaleX.fct, scaleX.inv, scaleY,scaleY.fct,
- 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))
- invisible(NULL)
-}
-
-cniperPoint <- function(L2Fam, neighbor, risk= asMSE(),
- lower=getdistrOption("DistrResolution"),
- upper=1-getdistrOption("DistrResolution")){
-
-
- 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)
-
- fun <- .getFunCnip(psi,eta, risk, L2Fam, neighbor at radius)
-
- 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,
- 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 = 19, jitter.fac = 1, with.lab = FALSE,
- 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)
- mcl <- as.list(mc[-1])
- dots <- as.list(mc$"...")
-
- robMod <- InfRobModel(center = L2Fam, neighbor = neighbor)
-
- mcl$IC1 <- optIC(model = L2Fam, risk = asCov())
- mcl$IC2 <- optIC(model = robMod, risk = risk)
- mcl$L2Fam <- NULL
- if(is.null(dots$ylab))
- mcl$ylab <- gettext("Asymptotic Risk difference (classic - robust)")
- if(is.null(dots$main))
- mcl$main <- gettext("Cniper point plot")
-
- if(withMaxRisk) mcl$fromCniperPlot <- TRUE
- do.call(cniperCont, mcl)
-}
-
-
-
+.isReplicated <- function(x, tol=.Machine$double.eps){
+ tx <- table(x)
+ rx <- as.numeric(names(tx[tx>1]))
+ sapply(x, function(y) any(abs(y-rx)<tol))
+}
+
+.plotData <- function(
+ ## helper function for cniper-type plots to plot in data
+ data, # data to be plot in
+ dots, # dots from the calling function
+ fun, # function to determine risk difference
+ L2Fam, # L2Family
+ IC, # IC1 in cniperContPlot and eta in cniperPointPlot
+ jit.fac,
+ jit.tol
+){
+ dotsP <- .makedotsP(dots)
+
+ al <- dotsP$alpha.trsp
+ if(!is.null(al)) 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(dots$lab.pts))
+ lab.pts <- rep(lab.pts, length.out=n)
+
+
+ sel <- .SelectOrderData(data, function(x)sapply(x,fun),
+ dots$which.lbs, dots$which.Order)
+ i.d <- sel$ind
+ i0.d <- sel$ind1
+ y.d <- sel$y
+ x.d <- sel$data
+ n <- length(i.d)
+
+ resc.dat <- .rescalefct(x.d, function(x) sapply(x,fun),
+ dots$scaleX, dots$scaleX.fct, dots$scaleX.inv,
+ dots$scaleY, dots$scaleY.fct,
+ dots$xlim, dots$ylim, dots)
+
+ if(any(.isReplicated(resc.dat$X, jit.tol))&&jit.fac>0)
+ resc.dat$X <- jitter(resc.dat$X, factor = jit.fac)
+ if(any(.isReplicated(resc.dat$Y, jit.tol))&&jit.fac>0)
+ resc.dat$Y <- jitter(resc.dat$Y, factor = jit.fac)
+
+ dotsP$scaleX <- dotsP$scaleY <- dotsP$scaleN <-NULL
+ dotsP$scaleX.fct <- dotsP$scaleY.fct <- NULL
+ dotsP$scaleX.inv <- dotsP$scaleY.inv <- NULL
+ dotsP$cex.pts <- dotsP$col.pts <- dotsP$lab.pts <- dotsP$pch.pts <- NULL
+ dotsP$jit.fac <- dotsP$with.lab <- dotsP$alpha.trsp <- NULL
+ dotsP$return.Order <- dotsP$cex.pts.fun <- NULL
+ dotsP$x.ticks <- dotsP$y.ticks <- NULL
+ dotsP$lab.font <- dotsP$which.lbs <- dotsP$which.lbs <- NULL
+
+ dotsP$x <- resc.dat$X
+ dotsP$y <- resc.dat$Y
+
+ trafo <- trafo(L2Fam at param)
+ dims <- nrow(trafo)
+ QF <- diag(dims)
+ 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]])
+ IC.rv <- as(diag(dims) %*% IC at Curve, "EuclRandVariable")
+ absy.f <- t(IC.rv) %*% QF %*% IC.rv
+ absy <- absInfoEval(x.d, absy.f)
+
+ if(is.null(dots$cex.pts)) dots$cex.pts <- par("cex")
+
+ dotsT <- dotsP
+ dotsT$cex <- dotsP$cex/2
+ dotsP$cex <- .cexscale(absy,absy,cex=dots$cex.pts, fun = dots$cex.pts.fun)
+ dotsP$col <- dots$col.pts
+
+ dotsT$pch <- NULL
+ dotsT$labels <- if(is.null(dots$lab.pts)) i.d else dots$lab.pts[i.d]
+ do.call(points,dotsP)
+ if(!is.null(dots$with.lab))
+ if(dots$with.lab) do.call(text,dotsT)
+ if(!is.null(dots$return.Order))
+ if(dots$return.Order) return(i0.d)
+ return(invisible(NULL))
+ }
+
+
+.getFunCnip <- function(IC1,IC2, risk, L2Fam, r, b20=NULL){
+
+ bType <- biastype(risk)
+ nType <- normtype(risk)
+ fnorm <- fct(nType)
+ riskfct <- getRiskFctBV(risk, bType)
+
+ .getTrVar <- function(IC){
+ R <- Risks(IC)[["trAsCov"]]
+ if(is.null(R)) R <- getRiskIC(IC, risk = trAsCov(), L2Fam = L2Fam)
+ if("trAsCov" %in% names(R)) R <- R[["trAsCov"]]
+ if(length(R) > 1) R <- R$value
+ return(c(R))
+ }
+ R1 <- .getTrVar (IC1)
+ R2 <- .getTrVar (IC2)
+
+ fun <- function(x){
+ y1 <- sapply(x, function(x1)evalIC(IC1,as.matrix(x1,ncol=1)))
+ b1 <- r*fnorm(y1)
+ r1 <- riskfct(var=R1,bias=b1)
+ if(!is.null(b20)){
+ r2 <- riskfct(var=R2, bias=b20)
+ }else{
+ y2 <- sapply(x,function(x0) evalIC(IC2,x0))
+ b2 <- r*fnorm(y2)
+ r2 <- riskfct(var=R2,bias=b2)
+ }
+ r1 - r2
+ return(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,
+ scaleY = FALSE, scaleY.fct = pnorm, scaleY.inv=qnorm,
+ scaleN = 9, x.ticks = NULL, y.ticks = NULL,
+ cex.pts = 1, cex.pts.fun = NULL, col.pts = par("col"),
+ pch.pts = 19, jit.fac = 1, jit.tol = .Machine$double.eps, with.lab = FALSE,
+ lab.pts = NULL, lab.font = NULL, alpha.trsp = NA,
+ which.lbs = NULL, which.Order = NULL,
+ return.Order = FALSE,
+ draw.nonlbl = TRUE, ## should non-labelled observations also be drawn?
+ cex.nonlbl = 0.3, ## character expansion(s) for non-labelled observations
+ pch.nonlbl = ".", ## plotting symbol(s) for non-labelled observations
+ withSubst = TRUE){
+
+ mcD <- match.call(expand.dots = FALSE)
+ dots <- as.list(mcD$"...")
+ mc <- match.call(#call = sys.call(sys.parent(1)),
+ expand.dots = TRUE)
+ mcl <- as.list(mc[-1])
+ IC1c <- as.character(deparse(IC1))
+ IC2c <- as.character(deparse(IC2))
+
+ .mpresubs <- if(withSubst){
+ function(inx)
+ .presubs(inx, c("%C1", "%A1", "%C2", "%A2", "%D" ),
+ c(as.character(class(IC1)[1]),
+ IC1c,
+ as.character(class(IC2)[1]),
+ IC2c,
+ as.character(date())
+ ))
+ }else function(inx)inx
+
+ if(!is.null(dots$main)) dots$main <- .mpresubs(dots$main)
+ if(!is.null(dots$sub)) dots$sub <- .mpresubs(dots$sub)
+ if(!is.null(dots$xlab)) dots$xlab <- .mpresubs(dots$xlab)
+ if(!is.null(dots$ylab)) dots$ylab <- .mpresubs(dots$ylab)
+
+ 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)
+
+ b20 <- NULL
+ fCpl <- eval(dots$fromCniperPlot)
+ if(!is.null(fCpl)&&length(Risks(IC2)))
+ 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)
+ }
+
+ if("lower" %in% names(as.list(mc))) lower <- p(L2Fam)(lower)
+ if("upper" %in% names(as.list(mc))) upper <- p(L2Fam)(upper)
+
+ x <- q(L2Fam)(seq(lower,upper,length=n))
+ 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)
+
+ dotsPl <- dots
+ dotsPl$x <- resc$X
+ dotsPl$y <- resc$Y
+ dotsPl$type <- "l"
+ if(is.null(dotsPl$main)) dotsPl$main <- gettext("Cniper region plot")
+ if(is.null(dotsPl$xlab)) dotsPl$xlab <- gettext("Dirac point")
+ if(is.null(dotsPl$ylab))
+ dotsPl$ylab <- gettext("Asymptotic Risk difference (IC1 - IC2)")
+
+ colSet <- ltySet <- lwdSet <- FALSE
+ if(!is.null(dotsPl$col)) {colSet <- TRUE; colo <- eval(dotsPl$col)}
+ if(colSet) {
+ colo <- rep(colo,length.out=2)
+ dotsPl$col <- colo[1]
+ }
+ if(!is.null(dotsPl$lwd)) {lwdSet <- TRUE; lwdo <- eval(dotsPl$lwd)}
+ if(lwdSet) {
+ lwdo <- rep(lwdo,length.out=2)
+ dotsPl$lwd <- lwdo[1]
+ }
+ if(!is.null(dotsPl$lty)) {ltySet <- TRUE; ltyo <- eval(dotsPl$lty)}
+ if(ltySet && ((!is.numeric(ltyo) && length(ltyo)==1)||
+ is.numeric(ltyo))){
+ ltyo <- list(ltyo,ltyo)
+ dotsPl$lty <- ltyo[[1]]
+ }else{ if (ltySet && !is.numeric(ltyo) && length(ltyo)==2){
+ dotsPl$lty <- ltyo[[1]]
+ }
+ }
+ do.call(plot,dotsPl)
+
+ dots$x <- dots$y <- NULL
+ dotsl <- .makedotsLowLevel(dots)
+ if(colSet) dotsl$col <- colo[2]
+ if(lwdSet) dotsl$lwd <- lwdo[2]
+ if(ltySet) dotsl$lty <- ltyo[[2]]
+
+ dotsl$h <- if(scaleY) scaleY.fct(0) else 0
+ do.call(abline, dotsl)
+
+ .plotRescaledAxis(scaleX, scaleX.fct, scaleX.inv, scaleY,scaleY.fct,
+ scaleY.inv, dots$xlim, dots$ylim, resc$X, ypts = 400,
+ n = scaleN, x.ticks = x.ticks, y.ticks = y.ticks)
+ if(!is.null(data)){
+ dots$scaleX <- scaleX
+ dots$scaleX.fct <- scaleX.fct
+ if(!is.null(mcl$scaleX.inv)) dots$scaleX.inv <- scaleX.inv
+ dots$scaleY <- scaleY
+ dots$scaleY.fct <- scaleY.fct
+ dots$scaleY.inv <- scaleY.inv
+ dots$scaleN <- scaleN
+ dots$x.ticks <- x.ticks
+ dots$y.ticks <- y.ticks
+ dots$cex.pts <- cex.pts
+ dots$cex.pts.fun <- cex.pts.fun
+ dots$col.pts <- col.pts
+ dots$pch.pts <- pch.pts
+ dots$jit.fac <- jit.fac
+ dots$jit.tol <- jit.tol
+ dots$with.lab <- with.lab
+ dots$lab.pts <- lab.pts
+ dots$lab.font <- lab.font
+ dots$alpha.trsp <- alpha.trsp
+ dots$which.lbs <- which.lbs
+ dots$which.Order <- which.Order
+ dots$return.Order <- return.Order
+
+ return(.plotData(data=data, dots=dots, fun=fun, L2Fam=L2Fam, IC=IC1, jit.fac=jit.fac, jit.tol=jit.tol))
+ }
+ invisible(NULL)
+}
+
+cniperPoint <- function(L2Fam, neighbor, risk= asMSE(),
+ lower=getdistrOption("DistrResolution"),
+ upper=1-getdistrOption("DistrResolution")){
+
+
+ mc <- match.call(expand.dots = FALSE)
+
+ if(is.null(as.list(mc)$lower)) lower <- q(L2Fam)(lower)
+ if(is.null(as.list(mc)$upper)) upper <- q(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)
+
+ fun <- .getFunCnip(psi,eta, risk, L2Fam, neighbor at radius)
+ 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,
+ 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, cex.pts.fun = NULL, col.pts = par("col"),
+ pch.pts = 19, jit.fac = 1, jit.tol = .Machine$double.eps,
+ with.lab = FALSE,
+ lab.pts = NULL, lab.font = NULL, alpha.trsp = NA,
+ which.lbs = NULL, which.Order = NULL,
+ return.Order = FALSE,
+ draw.nonlbl = TRUE, ## should non-labelled observations also be drawn?
+ cex.nonlbl = 0.3, ## character expansion(s) for non-labelled observations
+ pch.nonlbl = ".", ## plotting symbol(s) for non-labelled observations
+ withSubst = TRUE){
+
+ mc0 <- match.call(#call = sys.call(sys.parent(1)),
+ expand.dots = FALSE)
+ mc <- match.call(#call = sys.call(sys.parent(1)),
+ expand.dots = TRUE)
+ mcl <- as.list(mc[-1])
+ dots <- as.list(mc0$"...")
+ L2Famc <- as.character(deparse(L2Fam))
+
+ .mpresubs <- if(withSubst){
+ function(inx)
+ .presubs(inx, c("%C", "%A", "%D" ),
+ c(as.character(class(L2Fam)[1]),
+ L2Famc,
+ as.character(date())
+ ))
+ }else function(inx)inx
+
+ if(!is.null(dots$main)) dots$main <- .mpresubs(dots$main)
+ if(!is.null(dots$sub)) dots$sub <- .mpresubs(dots$sub)
+ if(!is.null(dots$xlab)) dots$xlab <- .mpresubs(dots$xlab)
+ if(!is.null(dots$ylab)) dots$ylab <- .mpresubs(dots$ylab)
+ if(is.null(mcl$risk)) mcl$risk <- asMSE()
+
+ robMod <- InfRobModel(center = L2Fam, neighbor = neighbor)
+
+ mcl$IC1 <- optIC(model = L2Fam, risk = asCov())
+ mcl$IC2 <- if(is(risk,"interpolRisk")){
+ getStartIC(model=L2Fam, risk = risk)
+ }else optIC(model = robMod, risk = risk)
+ mcl$L2Fam <- NULL
+ if(is.null(dots$ylab))
+ mcl$ylab <- gettext("Asymptotic Risk difference (classic - robust)")
+ if(is.null(dots$main))
+ mcl$main <- gettext("Cniper point plot")
+
+ if(withMaxRisk) mcl$fromCniperPlot <- TRUE
+ mcl$withMaxRisk <- NULL
+ mcl$withSubst <- FALSE
+ do.call(cniperCont, mcl)
+}
+
+
+
+ .cexscale <- function(y, y1=y, maxcex=4,mincex=0.05,cex, fun=NULL){
+ if(is.null(fun)) fun <- function(x) log(1+abs(x))
+ ly <- fun(y)
+ ly1 <- fun(unique(c(y,y1)))
+ my <- min(ly1,na.rm=TRUE)
+ My <- max(ly1,na.rm=TRUE)
+ ly0 <- (ly-my)/My
+ ly1 <- ly0*(maxcex-mincex)+mincex
+ return(cex*ly1)
+ }
Modified: pkg/ROptEst/R/comparePlot.R
===================================================================
--- pkg/ROptEst/R/comparePlot.R 2016-09-04 14:30:41 UTC (rev 893)
+++ pkg/ROptEst/R/comparePlot.R 2016-09-04 14:31:17 UTC (rev 894)
@@ -1,3 +1,5 @@
+.oldcomparePlot <- getMethod("comparePlot", signature("IC","IC"))
+
setMethod("comparePlot", signature("IC","IC"),
function(obj1,obj2, obj3 = NULL, obj4 = NULL, data = NULL,
..., withSweave = getdistrOption("withSweave"),
@@ -45,8 +47,7 @@
}
mcl$MBRB <- MBRB
mcl$withMBR <- withMBR
- do.call(getMethod("comparePlot", signature("IC","IC"),
- where="RobAStBase"), as.list(mcl[-1]),
+ do.call(.oldcomparePlot, as.list(mcl[-1]),
envir=parent.frame(2))
return(invisible())
})
Modified: pkg/ROptEst/R/getAsRisk.R
===================================================================
--- pkg/ROptEst/R/getAsRisk.R 2016-09-04 14:30:41 UTC (rev 893)
+++ pkg/ROptEst/R/getAsRisk.R 2016-09-04 14:31:17 UTC (rev 894)
@@ -114,6 +114,7 @@
A.start = A.start, maxiter = maxiter,
tol = tol, verbose = verbose)
erg <- eerg$b
+ bias <- 1/erg$value
return(list(asBias = bias, normtype = eerg$normtype))
})
Deleted: pkg/ROptEst/R/getFiRisk.R
===================================================================
--- pkg/ROptEst/R/getFiRisk.R 2016-09-04 14:30:41 UTC (rev 893)
+++ pkg/ROptEst/R/getFiRisk.R 2016-09-04 14:31:17 UTC (rev 894)
@@ -1,198 +0,0 @@
-###############################################################################
-## finite-sample under-/overshoot risk
-###############################################################################
-
-# cdf of truncated normal distribution
-ptnorm <- function(x, mu, A, B){
- ((A <= x)*(x <= B)*(pnorm(x-mu)-pnorm(A-mu))/(pnorm(B-mu)-pnorm(A-mu))
- + (x > B))
-}
-
-# n-fold convolution for truncated normal distributions
-conv.tnorm <- function(z, A, B, mu, n, m){
- if(n == 1) return(ptnorm(z, mu = mu, A = A, B = B))
- if(z <= n*A) return(0)
- if(z >= n*B) return(1)
-
- M <- 2^m
- h <- (B-A)/M
- x <- seq(from = A, to = B, by = h)
- p1 <- ptnorm(x, mu = mu, A = A, B = B)
- p1 <- p1[2:(M + 1)] - p1[1:M]
-
- ## FFT
- pn <- c(p1, numeric((n-1)*M))
-
- ## convolution theorem for DFTs
- pn <- Re(fft(fft(pn)^n, inverse = TRUE)) / (n*M)
- pn <- (abs(pn) >= .Machine$double.eps)*pn
- i.max <- n*M-(n-2)
- pn <- c(0,pn[1:i.max])
- pn <- cumsum(pn)
-
- ## cdf with continuity correction h/2
- x <- c(n*A,seq(from = n*A+n/2*h, to = n*B-n/2*h, by=h),n*B)
- pnfun1 <- approxfun(x = x+0.5*h, y = pn, yleft = 0, yright = pn[i.max+1])
- pnfun2 <- function(x) pnfun1(x) / pn[i.max+1]
-
- return(pnfun2(z))
-}
-
-
-setMethod("getFiRisk", signature(risk = "fiUnOvShoot",
- Distr = "Norm",
- neighbor = "ContNeighborhood"
- ),
- function(risk, Distr, neighbor, clip, stand,
- sampleSize, Algo, cont){
- eps <- neighbor at radius
- tau <- risk at width
- n <- sampleSize
- m <- getdistrOption("DefaultNrFFTGridPointsExponent")
-
- if(Algo == "B"){
- if(cont == "left"){
- delta1 <- (1-eps)*(pnorm(-clip+tau) + pnorm(-clip-tau)) + eps
- K1 <- dbinom(0:n, size = n, prob = delta1)
- P1 <- (1-eps)*pnorm(-clip-tau) + eps
- p1 <- P1/delta1
-
- summe1 <- numeric(n+1)
- summe1[1] <- 1 - conv.tnorm(z = 0, A = -clip, B = clip, mu = -tau, n = n, m = m)
- summe1[n+1] <- (1 - 0.5*(pbinom(q = n/2, size = n, prob = p1)
- + pbinom(q = n/2-0.1, size = n, prob = p1)))
- for(k in 1:(n-1)){
- j <- 0:k
- z <- clip*(k-2*j)
- P1.ste <- sapply(z, conv.tnorm, A = -clip, B = clip, mu = -tau, n = n-k, m = m)
- summe1[k+1] <- sum((1-P1.ste)*dbinom(j, size = k, prob = p1))
- }
- erg <- sum(summe1*K1)
- }else{
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/robast -r 894
More information about the Robast-commits
mailing list