[Robast-commits] r1057 - branches/robast-1.1/pkg/ROptEst branches/robast-1.1/pkg/ROptEst/R branches/robast-1.1/pkg/ROptEst/inst branches/robast-1.1/pkg/ROptEst/inst/scripts branches/robast-1.1/pkg/ROptEst/man branches/robast-1.2/pkg/ROptEst branches/robast-1.2/pkg/ROptEst/R branches/robast-1.2/pkg/ROptEst/inst branches/robast-1.2/pkg/ROptEst/inst/scripts branches/robast-1.2/pkg/ROptEst/man pkg/ROptEst pkg/ROptEst/R pkg/ROptEst/inst pkg/ROptEst/inst/scripts pkg/ROptEst/man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Jul 25 00:44:57 CEST 2018
Author: ruckdeschel
Date: 2018-07-25 00:44:57 +0200 (Wed, 25 Jul 2018)
New Revision: 1057
Modified:
branches/robast-1.1/pkg/ROptEst/DESCRIPTION
branches/robast-1.1/pkg/ROptEst/NAMESPACE
branches/robast-1.1/pkg/ROptEst/R/getIneffDiff.R
branches/robast-1.1/pkg/ROptEst/R/getInfV.R
branches/robast-1.1/pkg/ROptEst/R/getModifyIC.R
branches/robast-1.1/pkg/ROptEst/R/getRiskIC.R
branches/robast-1.1/pkg/ROptEst/R/getStartIC.R
branches/robast-1.1/pkg/ROptEst/R/roptest.new.R
branches/robast-1.1/pkg/ROptEst/inst/NEWS
branches/robast-1.1/pkg/ROptEst/inst/scripts/BinomialModel.R
branches/robast-1.1/pkg/ROptEst/inst/scripts/GammaModel.R
branches/robast-1.1/pkg/ROptEst/inst/scripts/MBRE.R
branches/robast-1.1/pkg/ROptEst/inst/scripts/NbinomModel.R
branches/robast-1.1/pkg/ROptEst/inst/scripts/NormalLocationScaleModel.R
branches/robast-1.1/pkg/ROptEst/man/getBiasIC.Rd
branches/robast-1.1/pkg/ROptEst/man/getRiskIC.Rd
branches/robast-1.2/pkg/ROptEst/DESCRIPTION
branches/robast-1.2/pkg/ROptEst/NAMESPACE
branches/robast-1.2/pkg/ROptEst/R/getIneffDiff.R
branches/robast-1.2/pkg/ROptEst/R/getInfV.R
branches/robast-1.2/pkg/ROptEst/R/getModifyIC.R
branches/robast-1.2/pkg/ROptEst/R/getRiskIC.R
branches/robast-1.2/pkg/ROptEst/R/getStartIC.R
branches/robast-1.2/pkg/ROptEst/R/roptest.new.R
branches/robast-1.2/pkg/ROptEst/inst/NEWS
branches/robast-1.2/pkg/ROptEst/inst/scripts/BinomialModel.R
branches/robast-1.2/pkg/ROptEst/inst/scripts/GammaModel.R
branches/robast-1.2/pkg/ROptEst/inst/scripts/MBRE.R
branches/robast-1.2/pkg/ROptEst/inst/scripts/NbinomModel.R
branches/robast-1.2/pkg/ROptEst/inst/scripts/NormalLocationScaleModel.R
branches/robast-1.2/pkg/ROptEst/man/getBiasIC.Rd
branches/robast-1.2/pkg/ROptEst/man/getRiskIC.Rd
pkg/ROptEst/DESCRIPTION
pkg/ROptEst/NAMESPACE
pkg/ROptEst/R/getIneffDiff.R
pkg/ROptEst/R/getInfV.R
pkg/ROptEst/R/getModifyIC.R
pkg/ROptEst/R/getRiskIC.R
pkg/ROptEst/R/getStartIC.R
pkg/ROptEst/R/roptest.new.R
pkg/ROptEst/inst/NEWS
pkg/ROptEst/inst/scripts/BinomialModel.R
pkg/ROptEst/inst/scripts/GammaModel.R
pkg/ROptEst/inst/scripts/MBRE.R
pkg/ROptEst/inst/scripts/NbinomModel.R
pkg/ROptEst/inst/scripts/NormalLocationScaleModel.R
pkg/ROptEst/man/getBiasIC.Rd
pkg/ROptEst/man/getRiskIC.Rd
Log:
[ROptEst] went through the scripts of ROptEst and where necessary fixed them... -> some bugfixes
+ DESCRIPTION now imports MASS, stats, graphics, utils
+ the return value of roptest (of class kStepEstimate) had in slot estimate.call the call to roptest,
but also, as an attribute the inner call to robest; subsequently, when printing the call this
cluttered the output -> changed this: the call to robest is now in an extra slot robestCall of class
OptionalCall which as a rule is NULL but is filled in roptest; it can be accessed by accessor
robestCall
+ some (non-standard)-scalenames were not correctly found in kStepEstimator; in addition, for nuisances, the restriction to main coordinates idx was missing in kStepEstimator
now have an internal function .fix.scalename to name the respective element as scale-coordinate
+ getBiasIC and getRiskIC now have argument withCheck in all methods
when checked, this is done through .checkICWithWarning (copied from RobAStBase)
+ in getStartIC when called from roptest / radiusMinimaxIC was missing argument debug
+ in getInvV.R we now need MASS::ginv to invert a non-square matrix (by using the generalized inverse)
+ in the self-standardized IC for radius minimax there was no mentioning of the radius to be searched for the lower BoundedWeight-class.R (getIneff.R)
+ in getModifyIC.R we matched ... wrongly
Modified: branches/robast-1.1/pkg/ROptEst/DESCRIPTION
===================================================================
--- branches/robast-1.1/pkg/ROptEst/DESCRIPTION 2018-07-24 22:31:37 UTC (rev 1056)
+++ branches/robast-1.1/pkg/ROptEst/DESCRIPTION 2018-07-24 22:44:57 UTC (rev 1057)
@@ -6,8 +6,8 @@
classes and methods.
Depends: R(>= 2.14.0), methods, distr(>= 2.5.2), distrEx(>= 2.5), distrMod(>= 2.5.2),
RandVar(>= 0.9.2), RobAStBase(>= 1.0)
-Imports: startupmsg
-Suggests: RobLox, MASS
+Imports: startupmsg, MASS, stats, graphics, utils
+Suggests: RobLox
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",
Modified: branches/robast-1.1/pkg/ROptEst/NAMESPACE
===================================================================
--- branches/robast-1.1/pkg/ROptEst/NAMESPACE 2018-07-24 22:31:37 UTC (rev 1056)
+++ branches/robast-1.1/pkg/ROptEst/NAMESPACE 2018-07-24 22:44:57 UTC (rev 1057)
@@ -12,6 +12,7 @@
"optimize", "pnorm", "qnorm", "uniroot")
importFrom("utils", "read.csv", "read.table", "str", "write.table")
importFrom("graphics", "abline")
+importFrom("MASS", "ginv")
exportClasses("asAnscombe", "asL1", "asL4")
exportMethods("optIC",
Modified: branches/robast-1.1/pkg/ROptEst/R/getIneffDiff.R
===================================================================
--- branches/robast-1.1/pkg/ROptEst/R/getIneffDiff.R 2018-07-24 22:31:37 UTC (rev 1056)
+++ branches/robast-1.1/pkg/ROptEst/R/getIneffDiff.R 2018-07-24 22:44:57 UTC (rev 1057)
@@ -89,7 +89,7 @@
normtype = upNorm, tol = eps,
numbeval = 1e4)
biasUp <- biasUpE$asBias$value
- ineffLo <- (p+biasLo^2*loRad^2)/loRisk
+ ineffLo <- (p+biasLo^2*radius^2)/loRisk
ineffUp <- if(upRad == Inf) biasUp^2/upRisk else
(p+biasUp^2*upRad^2)/upRisk
}else{
@@ -110,6 +110,7 @@
collapse = ""),"\n",sep="")
)
+# print(c(ineffLo,ineffUp))
if(withRetIneff) return(c(lo= ineffLo, up=ineffUp))
else return(ineffUp - ineffLo)
}else{
Modified: branches/robast-1.1/pkg/ROptEst/R/getInfV.R
===================================================================
--- branches/robast-1.1/pkg/ROptEst/R/getInfV.R 2018-07-24 22:31:37 UTC (rev 1056)
+++ branches/robast-1.1/pkg/ROptEst/R/getInfV.R 2018-07-24 22:44:57 UTC (rev 1057)
@@ -36,7 +36,12 @@
(weight(w)(evalRandVar(L2deriv, as.matrix(x)) [,,1]))^2
}
- cent0 <- solve(stand, cent)
+ .solve <- function(A0, b0) solve(A0,b0)
+ if(is.matrix(stand)){
+ if(nrow(stand)!=ncol(stand))
+ .solve <- function(A0,b0) MASS::ginv(A0)%*%b0
+ }
+ cent0 <- .solve(stand, cent)
integrandV <- function(x, L2.i, L2.j, i, j){
Modified: branches/robast-1.1/pkg/ROptEst/R/getModifyIC.R
===================================================================
--- branches/robast-1.1/pkg/ROptEst/R/getModifyIC.R 2018-07-24 22:31:37 UTC (rev 1056)
+++ branches/robast-1.1/pkg/ROptEst/R/getModifyIC.R 2018-07-24 22:44:57 UTC (rev 1057)
@@ -5,7 +5,7 @@
setMethod("getModifyIC", signature(L2FamIC = "L2ParamFamily",
neighbor = "Neighborhood", risk = "asRisk"),
function(L2FamIC, neighbor, risk, ...){
- dots <- list(...)
+ dots <- match.call(call = sys.call(sys.parent(1)), expand.dots=FALSE)[["..."]]
dots$verbose <- NULL
modIC <- function(L2Fam, IC){}
body(modIC) <- substitute({ verbose <- getRobAStBaseOption("all.verbose")
Modified: branches/robast-1.1/pkg/ROptEst/R/getRiskIC.R
===================================================================
--- branches/robast-1.1/pkg/ROptEst/R/getRiskIC.R 2018-07-24 22:31:37 UTC (rev 1056)
+++ branches/robast-1.1/pkg/ROptEst/R/getRiskIC.R 2018-07-24 22:44:57 UTC (rev 1057)
@@ -1,3 +1,15 @@
+.checkICWithWarning <- function(IC, L2Fam, tol){
+ if(!missing(L2Fam)){
+ prec <- checkIC(IC, L2Fam, out = FALSE)
+ }else{
+ prec <- checkIC(IC, out = FALSE)
+ }
+ if(prec > tol)
+ warning("The maximum deviation from the exact IC properties is ", prec,
+ "\nThis is larger than the specified 'tol' ",
+ "=> the result may be wrong")
+}
+
###############################################################################
## asymptotic covariance
###############################################################################
@@ -27,6 +39,7 @@
neighbor <- ContNeighborhood(1)
Cov <- getInfV(L2deriv = L2deriv, neighbor = neighbor,
biastype = biastype(IC), clip = c0, cent = z, stand = A)
+ if(withCheck) .checkICWithWarning(IC, L2Fam, tol=.Machine$double.eps^.25)
return(list(asCov = list(distribution = .getDistr(L2Fam),
value = Cov)))
}
@@ -39,15 +52,17 @@
risk = "asCov",
neighbor = "missing",
L2Fam = "L2ParamFamily"),
- function(IC, risk, L2Fam){
+ function(IC, risk, L2Fam, withCheck = TRUE){
Cov <- eval(IC at Risks[["asCov"]])
+ if(missing(withCheck)) withCheck <- TRUE
if (is.null(Cov)){
L2deriv <- L2Fam at L2derivDistr[[1]]
A <- as.vector(IC at stand)
c0 <- (IC at clipUp-IC@clipLo)/abs(A)
z <- IC at clipLo/abs(A)
neighbor <- TotalVarNeighborhood(1)
- Cov <- getInfV(L2deriv = L2deriv, neighbor = neighbor,
+ if(withCheck) .checkICWithWarning(IC, L2Fam, tol=.Machine$double.eps^.25)
+ Cov <- getInfV(L2deriv = L2deriv, neighbor = neighbor,
biastype = biastype(IC), clip = c0, cent = z, stand = A)
}
return(list(asCov = list(distribution = .getDistr(L2Fam), value = Cov)))
@@ -58,7 +73,7 @@
###############################################################################
setMethod("getBiasIC", signature(IC = "HampIC",
neighbor = "UncondNeighborhood"),
- function(IC, neighbor, L2Fam,...){
+ function(IC, neighbor, L2Fam,..., withCheck = TRUE){
if(missing(L2Fam))
L2Fam <- force(eval(IC at CallL2Fam))
@@ -70,7 +85,7 @@
setMethod("getBiasIC", signature(IC = "TotalVarIC",
neighbor = "UncondNeighborhood"),
- function(IC, neighbor, L2Fam,...){
+ function(IC, neighbor, L2Fam,..., withCheck = TRUE){
if(missing(L2Fam))
L2Fam <- force(eval(IC at CallL2Fam))
Modified: branches/robast-1.1/pkg/ROptEst/R/getStartIC.R
===================================================================
--- branches/robast-1.1/pkg/ROptEst/R/getStartIC.R 2018-07-24 22:31:37 UTC (rev 1056)
+++ branches/robast-1.1/pkg/ROptEst/R/getStartIC.R 2018-07-24 22:44:57 UTC (rev 1057)
@@ -6,6 +6,8 @@
..debug=FALSE){
mc <- match.call(expand.dots=FALSE, call = sys.call(sys.parent(1)))
dots <- as.list(mc$"...")
+# cat("HALLLO IN getstartIC\n"); print(..debug)
+ if(missing(..debug)||!is.logical(..debug)) ..debug <- FALSE
if("fsCor" %in% names(dots)){
fsCor <- eval(dots[["fsCor"]])
dots$fsCor <- NULL
Modified: branches/robast-1.1/pkg/ROptEst/R/roptest.new.R
===================================================================
--- branches/robast-1.1/pkg/ROptEst/R/roptest.new.R 2018-07-24 22:31:37 UTC (rev 1056)
+++ branches/robast-1.1/pkg/ROptEst/R/roptest.new.R 2018-07-24 22:44:57 UTC (rev 1057)
@@ -123,7 +123,7 @@
startCtrl = startCtrl, kStepCtrl = kStepCtrl,
na.rm = na.rm, ..., debug = ..withCheck,
withTimings = withTimings)
- attr(mc,"robest.call") <- retV at estimate.call
+ retV at robestCall <- retV at estimate.call
retV at estimate.call <- mc
return(retV)
}
@@ -231,6 +231,7 @@
startCtrl0 = startCtrl$initial.est.ArgList)
))
}else{
+ initial.est <- startCtrl$initial.est
print(substitute(kStepEstimator.start(initial.est0, x = x0,
nrvalues = nrvalues0, na.rm = na.rm0,
L2Fam = L2Fam0),list(x0=x,L2Fam0=L2Fam,
Modified: branches/robast-1.1/pkg/ROptEst/inst/NEWS
===================================================================
--- branches/robast-1.1/pkg/ROptEst/inst/NEWS 2018-07-24 22:31:37 UTC (rev 1056)
+++ branches/robast-1.1/pkg/ROptEst/inst/NEWS 2018-07-24 22:44:57 UTC (rev 1057)
@@ -26,12 +26,19 @@
+ DESCRIPTION tag SVNRevision changed to VCS/SVNRevision
+ getRiskIC and getBiasIC gain argument withCheck to speed up things if one does not want to call checkIC
+Return value of "roptest"
++ the return value of "roptest", an object of class "kStepEstimate" has a slot "estimate.call" which
+ contains the (matched) call to "roptest"; internally "roptest" calls "robest"; the call to "robest"
+ may be of interest, too, so we have a new slot "robestCall" of class "OptionalCall", ie a call
+ or NULL (default); it can be accessed via function robestCall()
+
Bugfix:
+ argument withMakeIC was not correctly used in roptest.new
under the hood:
+ wherever possible also use q.l internally instead of q to
provide functionality in IRKernel
++ fixed all scripts
#######################################
version 1.0.1
Modified: branches/robast-1.1/pkg/ROptEst/inst/scripts/BinomialModel.R
===================================================================
--- branches/robast-1.1/pkg/ROptEst/inst/scripts/BinomialModel.R 2018-07-24 22:31:37 UTC (rev 1056)
+++ branches/robast-1.1/pkg/ROptEst/inst/scripts/BinomialModel.R 2018-07-24 22:44:57 UTC (rev 1057)
@@ -1434,7 +1434,7 @@
#current radius: 0.6388282 inefficiency: 1.044571
#current radius: 0.6387762 inefficiency: 1.044584
# user system elapsed
-# 437.47 0.09 438.58
+# 125.41 1.20 143.03
r.rho1
@@ -1464,8 +1464,7 @@
#current radius: 0.3104966 inefficiency: 1.043323
#current radius: 0.3105373 inefficiency: 1.043323
# user system elapsed
-#2211.92 1.13 2235.70
-
+# 431.94 2.64 483.05
r.rho2
#$rho
Modified: branches/robast-1.1/pkg/ROptEst/inst/scripts/GammaModel.R
===================================================================
--- branches/robast-1.1/pkg/ROptEst/inst/scripts/GammaModel.R 2018-07-24 22:31:37 UTC (rev 1056)
+++ branches/robast-1.1/pkg/ROptEst/inst/scripts/GammaModel.R 2018-07-24 22:44:57 UTC (rev 1057)
@@ -74,6 +74,7 @@
RobG1 # show RobB1
## OBRE solution ARE 0.95 in ideal model
+## really need distrExOptions(ErelativeTolerance = 1e-8)
system.time(ICA <- optIC(model = RobG1, risk = asAnscombe(),
upper=NULL,lower=NULL, verbose=TRUE))
checkIC(ICA)
@@ -110,12 +111,12 @@
infoPlot(IC3)
## radius minimax IC
-## takes quite some time - about 30 min.
+## takes quite some time - about 180 sec.
system.time(IC4 <- radiusMinimaxIC(L2Fam=G, neighbor=ContNeighborhood(),
risk=asMSE(), loRad=0, upRad=Inf))
## least favorable radius
-## takes really long time - several hours!
+## takes really long time - 33 min!
#system.time(r.rho1 <- leastFavorableRadius(L2Fam=G, neighbor=ContNeighborhood(),
# risk=asMSE(), rho=0.5))
Modified: branches/robast-1.1/pkg/ROptEst/inst/scripts/MBRE.R
===================================================================
--- branches/robast-1.1/pkg/ROptEst/inst/scripts/MBRE.R 2018-07-24 22:31:37 UTC (rev 1056)
+++ branches/robast-1.1/pkg/ROptEst/inst/scripts/MBRE.R 2018-07-24 22:44:57 UTC (rev 1057)
@@ -6,16 +6,16 @@
### checks for lower case in various standardizations
N0 <- NormLocationScaleFamily(mean=-2, sd=3)
N0.Rob1<- InfRobModel(center = N0, neighbor = ContNeighborhood(radius = 15));
-N0.IC2 <- optIC(model = N0.Rob1, risk = asBias(), tol = 1e-10);print(stand(N0.IC2));print(cent(N0.IC2));print(stand(N0.IC2)/max(stand(N0.IC2)));print(cent(N0.IC2)/max(stand(N0.IC2)));print(clip(N0.IC2))
-plot(N0.IC2)
-N0.IC2.i <- optIC(model = N0.Rob1, risk = asMSE(), tol = 1e-10);print(stand(N0.IC2.i)/max(stand(N0.IC2.i)));print(cent(N0.IC2.i)/max(stand(N0.IC2.i)));print(clip(N0.IC2.i));print(stand(N0.IC2)/max(stand(N0.IC2)));print(cent(N0.IC2)/max(stand(N0.IC2)));print(clip(N0.IC2))
-plot(N0.IC2.i)
-N0.IC2.i <- optIC(model = N0.Rob1, risk = asMSE(normtype=InfoNorm()), tol = 1e-10);print(stand(N0.IC2.i)/max(stand(N0.IC2.i)));print(cent(N0.IC2.i)/max(stand(N0.IC2.i)));print(clip(N0.IC2.i));print(stand(N0.IC2)/max(stand(N0.IC2)));print(cent(N0.IC2)/max(stand(N0.IC2)));print(clip(N0.IC2))
-plot(N0.IC2.i)
-N0.IC2.i <- optIC(model = N0.Rob1, risk = asBias(normtype=InfoNorm()), tol = 1e-10);print(stand(N0.IC2.i));print(cent(N0.IC2.i));print(stand(N0.IC2.i)/max(stand(N0.IC2.i)));print(cent(N0.IC2.i)/max(stand(N0.IC2.i)));print(clip(N0.IC2.i));print(stand(N0.IC2)/max(stand(N0.IC2)));print(cent(N0.IC2)/max(stand(N0.IC2)));print(clip(N0.IC2))
-plot(N0.IC2.i)
-N0.IC2.i <- optIC(model = N0.Rob1, risk = asMSE(normtype=SelfNorm()), tol = 1e-10);print(stand(N0.IC2.i)/max(stand(N0.IC2.i)));print(cent(N0.IC2.i)/max(stand(N0.IC2.i)));print(clip(N0.IC2.i));print(stand(N0.IC2)/max(stand(N0.IC2)));print(cent(N0.IC2)/max(stand(N0.IC2)));print(clip(N0.IC2))
-plot(N0.IC2.i)
-N0.IC2.i <- optIC(model = N0.Rob1, risk = asBias(normtype=SelfNorm()), tol = 1e-10);print(stand(N0.IC2.i));print(cent(N0.IC2.i));print(stand(N0.IC2.i)/max(stand(N0.IC2.i)));print(cent(N0.IC2.i)/max(stand(N0.IC2.i)));print(clip(N0.IC2.i));print(stand(N0.IC2)/max(stand(N0.IC2)));print(cent(N0.IC2)/max(stand(N0.IC2)));print(clip(N0.IC2))
-plot(N0.IC2.i)
+N0.IC2.MBRE <- optIC(model = N0.Rob1, risk = asBias(), tol = 1e-10);print(stand(N0.IC2));print(cent(N0.IC2));print(stand(N0.IC2)/max(stand(N0.IC2)));print(cent(N0.IC2)/max(stand(N0.IC2)));print(clip(N0.IC2))
+plot(N0.IC2.MBRE)
+N0.IC2.OMSE <- optIC(model = N0.Rob1, risk = asMSE(), tol = 1e-10);print(stand(N0.IC2.i)/max(stand(N0.IC2.i)));print(cent(N0.IC2.i)/max(stand(N0.IC2.i)));print(clip(N0.IC2.i));print(stand(N0.IC2)/max(stand(N0.IC2)));print(cent(N0.IC2)/max(stand(N0.IC2)));print(clip(N0.IC2))
+plot(N0.IC2.OMSE)
+N0.IC2.MBRE.i <- optIC(model = N0.Rob1, risk = asBias(normtype=InfoNorm()), tol = 1e-10);print(stand(N0.IC2.i));print(cent(N0.IC2.i));print(stand(N0.IC2.i)/max(stand(N0.IC2.i)));print(cent(N0.IC2.i)/max(stand(N0.IC2.i)));print(clip(N0.IC2.i));print(stand(N0.IC2)/max(stand(N0.IC2)));print(cent(N0.IC2)/max(stand(N0.IC2)));print(clip(N0.IC2))
+plot(N0.IC2.MBRE.i)
+N0.IC2.OMSE.i <- optIC(model = N0.Rob1, risk = asMSE(normtype=InfoNorm()), tol = 1e-10);print(stand(N0.IC2.i)/max(stand(N0.IC2.i)));print(cent(N0.IC2.i)/max(stand(N0.IC2.i)));print(clip(N0.IC2.i));print(stand(N0.IC2)/max(stand(N0.IC2)));print(cent(N0.IC2)/max(stand(N0.IC2)));print(clip(N0.IC2))
+plot(N0.IC2.OMSE.i)
+N0.IC2.MBRE.s <- optIC(model = N0.Rob1, risk = asBias(normtype=SelfNorm()), tol = 1e-10);print(stand(N0.IC2.i));print(cent(N0.IC2.i));print(stand(N0.IC2.i)/max(stand(N0.IC2.i)));print(cent(N0.IC2.i)/max(stand(N0.IC2.i)));print(clip(N0.IC2.i));print(stand(N0.IC2)/max(stand(N0.IC2)));print(cent(N0.IC2)/max(stand(N0.IC2)));print(clip(N0.IC2))
+plot(N0.IC2.MBRE.s)
+N0.IC2.OMSE.s <- optIC(model = N0.Rob1, risk = asMSE(normtype=SelfNorm()), tol = 1e-10);print(stand(N0.IC2.i)/max(stand(N0.IC2.i)));print(cent(N0.IC2.i)/max(stand(N0.IC2.i)));print(clip(N0.IC2.i));print(stand(N0.IC2)/max(stand(N0.IC2)));print(cent(N0.IC2)/max(stand(N0.IC2)));print(clip(N0.IC2))
+plot(N0.IC2.OMSE.s)
}
Modified: branches/robast-1.1/pkg/ROptEst/inst/scripts/NbinomModel.R
===================================================================
--- branches/robast-1.1/pkg/ROptEst/inst/scripts/NbinomModel.R 2018-07-24 22:31:37 UTC (rev 1056)
+++ branches/robast-1.1/pkg/ROptEst/inst/scripts/NbinomModel.R 2018-07-24 22:44:57 UTC (rev 1057)
@@ -156,14 +156,15 @@
system.time(ICA <- optIC(model=RobN1, risk=asAnscombe(),
verbose=TRUE,lower=NULL,upper=10))
-
+# user system elapsed
+# 9.67 0.01 10.43
#-------------------------------------------------------------------------------
## MSE solution
#-------------------------------------------------------------------------------
system.time(IC1 <- optIC(model=RobN1, risk=asMSE()))
# user system elapsed
-# 10.53 0.02 11.21
+# 0.97 0.00 0.97
IC1
@@ -455,7 +456,7 @@
system.time(IC2 <- optIC(model=RobN2, risk=asMSE()))
# user system elapsed
-# 75.57 0.22 81.59
+# 2.46 0.00 2.46
IC2
@@ -1150,7 +1151,7 @@
system.time(IC7 <- radiusMinimaxIC(L2Fam=N, neighbor=ContNeighborhood(),
risk=asMSE(), loRad=0.01, upRad=3.9))
# user system elapsed
-# 33.26 0.02 33.64
+# 8.57 0.00 8.59
IC7
@@ -1267,7 +1268,7 @@
system.time(IC8 <- radiusMinimaxIC(L2Fam=N, neighbor=TotalVarNeighborhood(),
risk=asMSE(), loRad=0.01, upRad=1.8))
# user system elapsed
-# 565.58 0.21 586.05
+# 66.47 0.25 68.73
IC8
@@ -1405,7 +1406,7 @@
#current radius: 0.5626002 inefficiency: 1.044701
#current radius: 0.5625595 inefficiency: 1.044701
# user system elapsed
-# 361.25 0.12 369.05
+# 141.37 0.84 150.89
## same as for binomial????
@@ -1439,7 +1440,7 @@
#current radius: 0.2866482 inefficiency: 1.044425
#current radius: 0.2866889 inefficiency: 1.044456
# user system elapsed
-# 4891.07 1.90 5063.44
+# 707.48 3.17 760.09
r.rho2
@@ -2563,6 +2564,7 @@
IC12 <- radiusMinimaxIC(L2Fam=NbinomFamily(size=25, prob=estimate(est0)),
neighbor=TotalVarNeighborhood(), risk=asMSE(), loRad=0, upRad=Inf)
+(est4v <- kStepEstimator(x, IC=IC12, start=est0, steps = 3L))
#Evaluations of 3-step estimate:
#-------------------------------
Modified: branches/robast-1.1/pkg/ROptEst/inst/scripts/NormalLocationScaleModel.R
===================================================================
--- branches/robast-1.1/pkg/ROptEst/inst/scripts/NormalLocationScaleModel.R 2018-07-24 22:31:37 UTC (rev 1056)
+++ branches/robast-1.1/pkg/ROptEst/inst/scripts/NormalLocationScaleModel.R 2018-07-24 22:44:57 UTC (rev 1057)
@@ -81,16 +81,18 @@
plot(N0.IC4)
infoPlot(N0.IC4)
-(N0.IC4.i <- radiusMinimaxIC(L2Fam=N0, neighbor=ContNeighborhood(),
+system.time(N0.IC4.i <- radiusMinimaxIC(L2Fam=N0, neighbor=ContNeighborhood(),
risk=asMSE(normtype=InfoNorm()), loRad=0, upRad=Inf))
+print(N0.IC4.i)
checkIC(N0.IC4.i)
Risks(N0.IC4.i)
plot(N0.IC4.i)
infoPlot(N0.IC4.i)
## takes extremely long time:
-(N0.IC4.s <- radiusMinimaxIC(L2Fam=N0, neighbor=ContNeighborhood(),
+system.time(N0.IC4.s <- radiusMinimaxIC(L2Fam=N0, neighbor=ContNeighborhood(),
risk=asMSE(normtype=SelfNorm()), loRad=0, upRad=Inf))
+print(N0.IC4.s)
checkIC(N0.IC4.s)
Risks(N0.IC4.s)
plot(N0.IC4.s)
Modified: branches/robast-1.1/pkg/ROptEst/man/getBiasIC.Rd
===================================================================
--- branches/robast-1.1/pkg/ROptEst/man/getBiasIC.Rd 2018-07-24 22:31:37 UTC (rev 1056)
+++ branches/robast-1.1/pkg/ROptEst/man/getBiasIC.Rd 2018-07-24 22:44:57 UTC (rev 1057)
@@ -12,13 +12,16 @@
\usage{
getBiasIC(IC, neighbor, ...)
-\S4method{getBiasIC}{HampIC,UncondNeighborhood}(IC, neighbor, L2Fam, ...)
+\S4method{getBiasIC}{HampIC,UncondNeighborhood}(IC, neighbor, L2Fam, ..., withCheck = TRUE)
}
\arguments{
\item{IC}{ object of class \code{"InfluenceCurve"} }
\item{neighbor}{ object of class \code{"Neighborhood"}. }
\item{L2Fam}{ object of class \code{"L2ParamFamily"}. }
\item{\dots}{ additional parameters }
+ \item{withCheck}{logical: should a call to \code{checkIC} be done to
+ check accuracy (defaults to \code{TRUE}; ignored
+ if nothing is computed but simply a slot is read out).}
}
\details{ This function is rarely called directly. It is used by
other functions/methods. }
Modified: branches/robast-1.1/pkg/ROptEst/man/getRiskIC.Rd
===================================================================
--- branches/robast-1.1/pkg/ROptEst/man/getRiskIC.Rd 2018-07-24 22:31:37 UTC (rev 1056)
+++ branches/robast-1.1/pkg/ROptEst/man/getRiskIC.Rd 2018-07-24 22:44:57 UTC (rev 1057)
@@ -16,7 +16,7 @@
\S4method{getRiskIC}{HampIC,asCov,missing,missing}(IC, risk, withCheck= TRUE)
\S4method{getRiskIC}{HampIC,asCov,missing,L2ParamFamily}(IC, risk, L2Fam, withCheck= TRUE)
-\S4method{getRiskIC}{TotalVarIC,asCov,missing,L2ParamFamily}(IC, risk, L2Fam)
+\S4method{getRiskIC}{TotalVarIC,asCov,missing,L2ParamFamily}(IC, risk, L2Fam, withCheck = TRUE)
}
\arguments{
@@ -26,7 +26,8 @@
\item{\dots}{ additional parameters }
\item{L2Fam}{ object of class \code{"L2ParamFamily"}. }
\item{withCheck}{logical: should a call to \code{checkIC} be done to
- check accuracy (defaults to \code{TRUE}).}
+ check accuracy (defaults to \code{TRUE}; ignored
+ if nothing is computed but simply a slot is read out).}
}
\details{To make sure that the results are valid, it is recommended
to include an additional check of the IC properties of \code{IC}
Modified: branches/robast-1.2/pkg/ROptEst/DESCRIPTION
===================================================================
--- branches/robast-1.2/pkg/ROptEst/DESCRIPTION 2018-07-24 22:31:37 UTC (rev 1056)
+++ branches/robast-1.2/pkg/ROptEst/DESCRIPTION 2018-07-24 22:44:57 UTC (rev 1057)
@@ -1,20 +1,22 @@
Package: ROptEst
Version: 1.1.0
-Date: 2018-07-23
+Date: 2018-07-17
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.5), distrMod(>= 2.5.2), RandVar(>= 0.9.2), RobAStBase(>=
- 1.0)
-Imports: startupmsg
-Suggests: RobLox, MASS
-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")))
+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.5), distrMod(>= 2.5.2),
+ RandVar(>= 0.9.2), RobAStBase(>= 1.0)
+Imports: startupmsg, MASS, stats, graphics, utils
+Suggests: RobLox
+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
License: LGPL-3
URL: http://robast.r-forge.r-project.org/
Encoding: latin1
LastChangedDate: {$LastChangedDate$}
LastChangedRevision: {$LastChangedRevision$}
-VCS/SVNRevision: 1040
+VCS/SVNRevision: 940
Modified: branches/robast-1.2/pkg/ROptEst/NAMESPACE
===================================================================
--- branches/robast-1.2/pkg/ROptEst/NAMESPACE 2018-07-24 22:31:37 UTC (rev 1056)
+++ branches/robast-1.2/pkg/ROptEst/NAMESPACE 2018-07-24 22:44:57 UTC (rev 1057)
@@ -12,6 +12,7 @@
"optimize", "pnorm", "qnorm", "uniroot")
importFrom("utils", "read.csv", "read.table", "str", "write.table")
importFrom("graphics", "abline")
+importFrom("MASS", "ginv")
exportClasses("asAnscombe", "asL1", "asL4")
exportMethods("optIC",
Modified: branches/robast-1.2/pkg/ROptEst/R/getIneffDiff.R
===================================================================
--- branches/robast-1.2/pkg/ROptEst/R/getIneffDiff.R 2018-07-24 22:31:37 UTC (rev 1056)
+++ branches/robast-1.2/pkg/ROptEst/R/getIneffDiff.R 2018-07-24 22:44:57 UTC (rev 1057)
@@ -89,7 +89,7 @@
normtype = upNorm, tol = eps,
numbeval = 1e4)
biasUp <- biasUpE$asBias$value
- ineffLo <- (p+biasLo^2*loRad^2)/loRisk
+ ineffLo <- (p+biasLo^2*radius^2)/loRisk
ineffUp <- if(upRad == Inf) biasUp^2/upRisk else
(p+biasUp^2*upRad^2)/upRisk
}else{
@@ -110,6 +110,7 @@
collapse = ""),"\n",sep="")
)
+# print(c(ineffLo,ineffUp))
if(withRetIneff) return(c(lo= ineffLo, up=ineffUp))
else return(ineffUp - ineffLo)
}else{
Modified: branches/robast-1.2/pkg/ROptEst/R/getInfV.R
===================================================================
--- branches/robast-1.2/pkg/ROptEst/R/getInfV.R 2018-07-24 22:31:37 UTC (rev 1056)
+++ branches/robast-1.2/pkg/ROptEst/R/getInfV.R 2018-07-24 22:44:57 UTC (rev 1057)
@@ -36,7 +36,12 @@
(weight(w)(evalRandVar(L2deriv, as.matrix(x)) [,,1]))^2
}
- cent0 <- solve(stand, cent)
+ .solve <- function(A0, b0) solve(A0,b0)
+ if(is.matrix(stand)){
+ if(nrow(stand)!=ncol(stand))
+ .solve <- function(A0,b0) MASS::ginv(A0)%*%b0
+ }
+ cent0 <- .solve(stand, cent)
integrandV <- function(x, L2.i, L2.j, i, j){
Modified: branches/robast-1.2/pkg/ROptEst/R/getModifyIC.R
===================================================================
--- branches/robast-1.2/pkg/ROptEst/R/getModifyIC.R 2018-07-24 22:31:37 UTC (rev 1056)
+++ branches/robast-1.2/pkg/ROptEst/R/getModifyIC.R 2018-07-24 22:44:57 UTC (rev 1057)
@@ -5,7 +5,7 @@
setMethod("getModifyIC", signature(L2FamIC = "L2ParamFamily",
neighbor = "Neighborhood", risk = "asRisk"),
function(L2FamIC, neighbor, risk, ...){
- dots <- list(...)
+ dots <- match.call(call = sys.call(sys.parent(1)), expand.dots=FALSE)[["..."]]
dots$verbose <- NULL
modIC <- function(L2Fam, IC){}
body(modIC) <- substitute({ verbose <- getRobAStBaseOption("all.verbose")
Modified: branches/robast-1.2/pkg/ROptEst/R/getRiskIC.R
===================================================================
--- branches/robast-1.2/pkg/ROptEst/R/getRiskIC.R 2018-07-24 22:31:37 UTC (rev 1056)
+++ branches/robast-1.2/pkg/ROptEst/R/getRiskIC.R 2018-07-24 22:44:57 UTC (rev 1057)
@@ -1,3 +1,15 @@
+.checkICWithWarning <- function(IC, L2Fam, tol){
+ if(!missing(L2Fam)){
+ prec <- checkIC(IC, L2Fam, out = FALSE)
+ }else{
+ prec <- checkIC(IC, out = FALSE)
+ }
+ if(prec > tol)
+ warning("The maximum deviation from the exact IC properties is ", prec,
+ "\nThis is larger than the specified 'tol' ",
+ "=> the result may be wrong")
+}
+
###############################################################################
## asymptotic covariance
###############################################################################
@@ -27,6 +39,7 @@
neighbor <- ContNeighborhood(1)
Cov <- getInfV(L2deriv = L2deriv, neighbor = neighbor,
biastype = biastype(IC), clip = c0, cent = z, stand = A)
+ if(withCheck) .checkICWithWarning(IC, L2Fam, tol=.Machine$double.eps^.25)
return(list(asCov = list(distribution = .getDistr(L2Fam),
value = Cov)))
}
@@ -39,15 +52,17 @@
risk = "asCov",
neighbor = "missing",
L2Fam = "L2ParamFamily"),
- function(IC, risk, L2Fam){
+ function(IC, risk, L2Fam, withCheck = TRUE){
Cov <- eval(IC at Risks[["asCov"]])
+ if(missing(withCheck)) withCheck <- TRUE
if (is.null(Cov)){
L2deriv <- L2Fam at L2derivDistr[[1]]
A <- as.vector(IC at stand)
c0 <- (IC at clipUp-IC@clipLo)/abs(A)
z <- IC at clipLo/abs(A)
neighbor <- TotalVarNeighborhood(1)
- Cov <- getInfV(L2deriv = L2deriv, neighbor = neighbor,
+ if(withCheck) .checkICWithWarning(IC, L2Fam, tol=.Machine$double.eps^.25)
+ Cov <- getInfV(L2deriv = L2deriv, neighbor = neighbor,
biastype = biastype(IC), clip = c0, cent = z, stand = A)
}
return(list(asCov = list(distribution = .getDistr(L2Fam), value = Cov)))
@@ -58,7 +73,7 @@
###############################################################################
setMethod("getBiasIC", signature(IC = "HampIC",
neighbor = "UncondNeighborhood"),
- function(IC, neighbor, L2Fam,...){
+ function(IC, neighbor, L2Fam,..., withCheck = TRUE){
if(missing(L2Fam))
L2Fam <- force(eval(IC at CallL2Fam))
@@ -70,7 +85,7 @@
setMethod("getBiasIC", signature(IC = "TotalVarIC",
neighbor = "UncondNeighborhood"),
- function(IC, neighbor, L2Fam,...){
+ function(IC, neighbor, L2Fam,..., withCheck = TRUE){
if(missing(L2Fam))
L2Fam <- force(eval(IC at CallL2Fam))
Modified: branches/robast-1.2/pkg/ROptEst/R/getStartIC.R
===================================================================
--- branches/robast-1.2/pkg/ROptEst/R/getStartIC.R 2018-07-24 22:31:37 UTC (rev 1056)
+++ branches/robast-1.2/pkg/ROptEst/R/getStartIC.R 2018-07-24 22:44:57 UTC (rev 1057)
@@ -6,6 +6,8 @@
..debug=FALSE){
mc <- match.call(expand.dots=FALSE, call = sys.call(sys.parent(1)))
dots <- as.list(mc$"...")
+# cat("HALLLO IN getstartIC\n"); print(..debug)
+ if(missing(..debug)||!is.logical(..debug)) ..debug <- FALSE
if("fsCor" %in% names(dots)){
fsCor <- eval(dots[["fsCor"]])
dots$fsCor <- NULL
Modified: branches/robast-1.2/pkg/ROptEst/R/roptest.new.R
===================================================================
--- branches/robast-1.2/pkg/ROptEst/R/roptest.new.R 2018-07-24 22:31:37 UTC (rev 1056)
+++ branches/robast-1.2/pkg/ROptEst/R/roptest.new.R 2018-07-24 22:44:57 UTC (rev 1057)
@@ -123,7 +123,7 @@
startCtrl = startCtrl, kStepCtrl = kStepCtrl,
na.rm = na.rm, ..., debug = ..withCheck,
withTimings = withTimings)
- attr(mc,"robest.call") <- retV at estimate.call
+ retV at robestCall <- retV at estimate.call
retV at estimate.call <- mc
return(retV)
}
@@ -231,6 +231,7 @@
startCtrl0 = startCtrl$initial.est.ArgList)
))
}else{
+ initial.est <- startCtrl$initial.est
print(substitute(kStepEstimator.start(initial.est0, x = x0,
nrvalues = nrvalues0, na.rm = na.rm0,
L2Fam = L2Fam0),list(x0=x,L2Fam0=L2Fam,
Modified: branches/robast-1.2/pkg/ROptEst/inst/NEWS
===================================================================
--- branches/robast-1.2/pkg/ROptEst/inst/NEWS 2018-07-24 22:31:37 UTC (rev 1056)
+++ branches/robast-1.2/pkg/ROptEst/inst/NEWS 2018-07-24 22:44:57 UTC (rev 1057)
@@ -26,12 +26,19 @@
+ DESCRIPTION tag SVNRevision changed to VCS/SVNRevision
+ getRiskIC and getBiasIC gain argument withCheck to speed up things if one does not want to call checkIC
+Return value of "roptest"
++ the return value of "roptest", an object of class "kStepEstimate" has a slot "estimate.call" which
+ contains the (matched) call to "roptest"; internally "roptest" calls "robest"; the call to "robest"
+ may be of interest, too, so we have a new slot "robestCall" of class "OptionalCall", ie a call
+ or NULL (default); it can be accessed via function robestCall()
+
Bugfix:
+ argument withMakeIC was not correctly used in roptest.new
under the hood:
+ wherever possible also use q.l internally instead of q to
provide functionality in IRKernel
++ fixed all scripts
#######################################
version 1.0.1
Modified: branches/robast-1.2/pkg/ROptEst/inst/scripts/BinomialModel.R
===================================================================
--- branches/robast-1.2/pkg/ROptEst/inst/scripts/BinomialModel.R 2018-07-24 22:31:37 UTC (rev 1056)
+++ branches/robast-1.2/pkg/ROptEst/inst/scripts/BinomialModel.R 2018-07-24 22:44:57 UTC (rev 1057)
@@ -1434,7 +1434,7 @@
#current radius: 0.6388282 inefficiency: 1.044571
#current radius: 0.6387762 inefficiency: 1.044584
# user system elapsed
-# 437.47 0.09 438.58
+# 125.41 1.20 143.03
r.rho1
@@ -1464,8 +1464,7 @@
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/robast -r 1057
More information about the Robast-commits
mailing list