[Robast-commits] r1233 - in branches/robast-1.3/pkg/ROptEst: R inst inst/scripts
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Jul 28 16:39:58 CEST 2021
Author: ruckdeschel
Date: 2021-07-28 16:39:58 +0200 (Wed, 28 Jul 2021)
New Revision: 1233
Modified:
branches/robast-1.3/pkg/ROptEst/R/CheckMakeContIC.R
branches/robast-1.3/pkg/ROptEst/R/L1L2normL2deriv.R
branches/robast-1.3/pkg/ROptEst/R/LowerCaseMultivariate.R
branches/robast-1.3/pkg/ROptEst/R/getInfCent.R
branches/robast-1.3/pkg/ROptEst/R/getInfGamma.R
branches/robast-1.3/pkg/ROptEst/R/getInfStand.R
branches/robast-1.3/pkg/ROptEst/R/getInfV.R
branches/robast-1.3/pkg/ROptEst/R/roptest.new.R
branches/robast-1.3/pkg/ROptEst/R/updateNorm.R
branches/robast-1.3/pkg/ROptEst/inst/NEWS
branches/robast-1.3/pkg/ROptEst/inst/scripts/NormalLocationScaleModel.R
Log:
[ROptEst] branch 1.3:
Bug-Fix under the hood:
+ in calls of form do.call(E, .....) we only use functions with one argument;
(in calls where a function was passed on as argument, this threw errors...)
-- to check this try:
#----------------------------------------------------------------
library(RobExtremes)
## robust starting estimator
GP <- GParetoFamily()
## original data set AM1 missing !!!
### instead artificially generate it
set.seed(20210713)
AM1 <- c(1.7 + r(GPareto(loc=0, scale=14, shape=1.5))(97),
10^c(10,14,15))
est0 <- medkMAD(AM1-1.6, GP, k=10)
## opt-robust estimator
roptest(x = AM1 - 1.6, GP)
## Geschwindigkeit:
RMXEstimator(x = AM1 - 1.6, GP)
#----------------------------------------------------------------
Modified: branches/robast-1.3/pkg/ROptEst/R/CheckMakeContIC.R
===================================================================
--- branches/robast-1.3/pkg/ROptEst/R/CheckMakeContIC.R 2021-04-01 16:47:39 UTC (rev 1232)
+++ branches/robast-1.3/pkg/ROptEst/R/CheckMakeContIC.R 2021-07-28 14:39:58 UTC (rev 1233)
@@ -208,8 +208,8 @@
res2 <- numeric(nrvalues)
for(i in 1:nrvalues){
if(z.comp[i]){
- Eargs <- c(list(object = Distr, fun = integrand2,
- L2.i = L2deriv at Map[[i]]), dotsI)
+ integrand2i <- function(x) integrand2(x,L2deriv at Map[[i]])
+ Eargs <- c(list(object = Distr, fun = integrand2i), dotsI)
res2[i] <- buf <- do.call(E,Eargs)
if(diagnostic){k <- k + 1; diagn[[k]] <- attr(buf,"diagnostic")}
}else{
@@ -229,9 +229,9 @@
for(i in 1:nrvalues){
for(j in i:nrvalues){
if(A.comp[i,j]){
- Eargs <- c(list(object = Distr, fun = integrandA,
- L2.i = L2deriv at Map[[i]],
- L2.j = L2deriv at Map[[j]], i = i, j = j), dotsI)
+ integrandAij <- function(x) integrandA(x, L2.i = L2deriv at Map[[i]],
+ L2.j = L2deriv at Map[[j]], i = i, j = j)
+ Eargs <- c(list(object = Distr, fun = integrandAij), dotsI)
erg[i, j] <- buf <- do.call(E,Eargs)
if(diagnostic){k <- k + 1; diagn[[k]] <- attr(buf,"diagnostic")}
}
Modified: branches/robast-1.3/pkg/ROptEst/R/L1L2normL2deriv.R
===================================================================
--- branches/robast-1.3/pkg/ROptEst/R/L1L2normL2deriv.R 2021-04-01 16:47:39 UTC (rev 1232)
+++ branches/robast-1.3/pkg/ROptEst/R/L1L2normL2deriv.R 2021-07-28 14:39:58 UTC (rev 1233)
@@ -12,14 +12,13 @@
dotsI <- .filterEargsWEargList(list(...))
if(is.null(dotsI$useApply)) dotsI$useApply <- FALSE
- integrandG <- function(x, L2, stand, cent){
- X <- evalRandVar(L2, as.matrix(x))[,,1] - cent
+ integrandG <- function(x){
+ X <- evalRandVar(L2deriv, as.matrix(x))[,,1] - cent
Y <- apply(X, 2, "%*%", t(stand))
res <- fct(normtype)(Y)
return((res > 0)*res)
}
-
- return(do.call(E, c(list(object = Distr, fun = integrandG, L2 = L2deriv,
- stand = stand, cent = cent),dotsI)))
+ retval <- do.call(E, c(list(object = Distr, fun = integrandG),dotsI))
+ return(retval)
})
Modified: branches/robast-1.3/pkg/ROptEst/R/LowerCaseMultivariate.R
===================================================================
--- branches/robast-1.3/pkg/ROptEst/R/LowerCaseMultivariate.R 2021-04-01 16:47:39 UTC (rev 1232)
+++ branches/robast-1.3/pkg/ROptEst/R/LowerCaseMultivariate.R 2021-07-28 14:39:58 UTC (rev 1233)
@@ -62,8 +62,9 @@
w <<- w0
}
- E1 <- do.call(E,c(list(object = Distr, fun = abs.fct, L2 = L2deriv, stand = A,
- cent = z, normtype.0 = normtype), dotsI))
+ abs.fct.0 <- function(x) abs.fct(x, L2deriv, A, z, normtype)
+
+ E1 <- do.call(E,c(list(object = Distr, fun = abs.fct.0), dotsI))
stA <- if (is(normtype,"QFNorm"))
QuadForm(normtype)%*%A else A
# erg <- E1/sum(diag(stA %*% t(trafo)))
@@ -130,8 +131,10 @@
p <- 1
A <- matrix(param, ncol = k, nrow = 1)
# print(A)
- E1 <- do.call(E, c(list( object = Distr, fun = pos.fct,
- L2 = L2deriv, stand = A), dotsI))
+
+ pos.fct.0 <- function(x) pos.fct(x, L2deriv, A)
+
+ E1 <- do.call(E, c(list( object = Distr, fun = pos.fct.0), dotsI))
erg <- E1/sum(diag(A %*% t(trafo)))
return(erg)
}
@@ -145,14 +148,14 @@
b <- 1/erg$value
stand(w) <- A
- pr.fct <- function(x, L2, pr.sign=1){
- X <- evalRandVar(L2, as.matrix(x)) [,,1]
+ pr.fct <- function(x, pr.sign=1){
+ X <- evalRandVar(L2deriv, as.matrix(x)) [,,1]
Y <- as.numeric(A %*% X)
return(as.numeric(pr.sign*Y>0))
}
- p.p <- do.call(E, c(list( object = Distr, fun = pr.fct, L2 = L2deriv,
+ p.p <- do.call(E, c(list( object = Distr, fun = pr.fct,
pr.sign = 1), dotsI))
- m.p <- do.call(E, c(list( object = Distr, fun = pr.fct, L2 = L2deriv,
+ m.p <- do.call(E, c(list( object = Distr, fun = pr.fct,
pr.sign = -1), dotsI))
Modified: branches/robast-1.3/pkg/ROptEst/R/getInfCent.R
===================================================================
--- branches/robast-1.3/pkg/ROptEst/R/getInfCent.R 2021-04-01 16:47:39 UTC (rev 1232)
+++ branches/robast-1.3/pkg/ROptEst/R/getInfCent.R 2021-07-28 14:39:58 UTC (rev 1233)
@@ -85,8 +85,9 @@
res2 <- numeric(nrvalues)
for(i in 1:nrvalues){
if(z.comp[i]){
- res2[i] <- do.call(E, c(list(object = Distr, fun = integrand2,
- L2.i = L2deriv at Map[[i]]), dotsI))
+ integrand2i <- function(x) integrand2(x, L2.i = L2deriv at Map[[i]])
+ res2[i] <- do.call(E, c(list(object = Distr, fun = integrand2i),
+ dotsI))
}else{
res2[i] <- 0
}
Modified: branches/robast-1.3/pkg/ROptEst/R/getInfGamma.R
===================================================================
--- branches/robast-1.3/pkg/ROptEst/R/getInfGamma.R 2021-04-01 16:47:39 UTC (rev 1232)
+++ branches/robast-1.3/pkg/ROptEst/R/getInfGamma.R 2021-07-28 14:39:58 UTC (rev 1233)
@@ -34,8 +34,8 @@
dotsI <- .filterEargsWEargList(list(...))
if(is.null(dotsI$useApply)) dotsI$useApply <- FALSE
- integrandG <- function(x, L2, stand, cent, clip){
- X <- evalRandVar(L2, as.matrix(x))[,,1] - cent
+ integrandG <- function(x){
+ X <- evalRandVar(L2deriv, as.matrix(x))[,,1] - cent
Y <- stand %*% X
res <- norm(risk)(Y) - clip
@@ -42,8 +42,7 @@
return((res > 0)*res^power)
}
- res <- do.call(E, c(list(object = Distr, fun = integrandG, L2 = L2deriv,
- stand = stand, cent = cent, clip = clip),dotsI))
+ res <- do.call(E, c(list(object = Distr, fun = integrandG),dotsI))
return(-res)
})
@@ -57,8 +56,8 @@
dotsI <- .filterEargsWEargList(list(...))
if(is.null(dotsI$useApply)) dotsI$useApply <- FALSE
- integrandG <- function(x, L2, stand, cent, clip){
- X <- evalRandVar(L2, as.matrix(x))[,,1] - cent
+ integrandG <- function(x){
+ X <- evalRandVar(L2deriv, as.matrix(x))[,,1] - cent
Y <- stand %*% X
res <- Y - clip
@@ -65,8 +64,7 @@
return((res > 0)*res^power)
}
- res <- do.call(E, c(list(object = Distr, fun = integrandG, L2 = L2deriv,
- stand = stand, cent = cent, clip = clip),dotsI))
+ res <- do.call(E, c(list(object = Distr, fun = integrandG),dotsI))
return(-res)
})
###############################################################################
Modified: branches/robast-1.3/pkg/ROptEst/R/getInfStand.R
===================================================================
--- branches/robast-1.3/pkg/ROptEst/R/getInfStand.R 2021-04-01 16:47:39 UTC (rev 1232)
+++ branches/robast-1.3/pkg/ROptEst/R/getInfStand.R 2021-07-28 14:39:58 UTC (rev 1233)
@@ -39,9 +39,10 @@
for(i in 1:nrvalues)
for(j in i:nrvalues)
if(A.comp[i,j]){
- erg[i, j] <- do.call(E, c(list(object = Distr, fun = integrandA,
- L2.i = L2deriv at Map[[i]],
- L2.j = L2deriv at Map[[j]], i = i, j = j), dotsI))
+ integrandAij <- function(x) integrandA(x,L2.i = L2deriv at Map[[i]],
+ L2.j = L2deriv at Map[[j]], i = i, j = j)
+ erg[i, j] <- do.call(E, c(list(object = Distr, fun = integrandAij),
+ dotsI))
}
erg[col(erg) < row(erg)] <- t(erg)[col(erg) < row(erg)]
Modified: branches/robast-1.3/pkg/ROptEst/R/getInfV.R
===================================================================
--- branches/robast-1.3/pkg/ROptEst/R/getInfV.R 2021-04-01 16:47:39 UTC (rev 1232)
+++ branches/robast-1.3/pkg/ROptEst/R/getInfV.R 2021-07-28 14:39:58 UTC (rev 1233)
@@ -56,9 +56,10 @@
for(i in 1:nrvalues)
for(j in i:nrvalues)
if(V.comp[i,j]){
- eArgs <- c(list(object = Distr, fun = integrandV,
+ integrandVij <- function(x) integrandV(x,
L2.i = L2deriv at Map[[i]],
- L2.j = L2deriv at Map[[j]], i = i, j = j), dotsI)
+ L2.j = L2deriv at Map[[j]], i = i, j = j)
+ eArgs <- c(list(object = Distr, fun = integrandVij), dotsI)
erg[i, j] <- do.call(E,eArgs)
}
erg[col(erg) < row(erg)] <- t(erg)[col(erg) < row(erg)]
Modified: branches/robast-1.3/pkg/ROptEst/R/roptest.new.R
===================================================================
--- branches/robast-1.3/pkg/ROptEst/R/roptest.new.R 2021-04-01 16:47:39 UTC (rev 1232)
+++ branches/robast-1.3/pkg/ROptEst/R/roptest.new.R 2021-07-28 14:39:58 UTC (rev 1233)
@@ -390,7 +390,8 @@
withLogScale = kStepCtrl$withLogScale,
withEvalAsVar = withEvalAsVarkStep,
withMakeIC = withMakeICkStep)
- print(argList) }
+ print(argList)
+ }
sy.kStep <- system.time({
kStepArgList <- list(x, IC = ICstart, start = initial.est,
steps = steps, useLast = kStepCtrl$useLast,
@@ -406,8 +407,14 @@
nms <- names(kStepCtrl$E.arglist)
for(nmi in nms) kStepArgList[[nmi]] <- kStepCtrl$E.arglist[[nmi]]
}
- res <- do.call(kStepEstimator, kStepArgList)
- })
+ if(debug){
+ print(substitute({
+ res <- do.call(kStepEstimator, kStepArgList0)
+ }, list(kStepArgList0=kStepEstimator)))
+ }else{
+ res <- do.call(kStepEstimator, kStepArgList)
+ }
+ })
sy.OnlykStep <- attr(res,"timings")
kStepDiagn <- attr(res,"diagnostic")
if (withTimings) print(sy.kStep)
Modified: branches/robast-1.3/pkg/ROptEst/R/updateNorm.R
===================================================================
--- branches/robast-1.3/pkg/ROptEst/R/updateNorm.R 2021-04-01 16:47:39 UTC (rev 1232)
+++ branches/robast-1.3/pkg/ROptEst/R/updateNorm.R 2021-07-28 14:39:58 UTC (rev 1233)
@@ -1,9 +1,10 @@
-#setMethod("updateNorm", "NormType", function(normtype, ...) normtype)
-#setMethod("updateNorm", "InfoNorm", function(normtype, FI, ...)
+# setMethod("updateNorm", "NormType", function(normtype, ...) normtype)
+# setMethod("updateNorm", "InfoNorm", function(normtype, FI, ...)
# {QuadForm(normtype) <- PosSemDefSymmMatrix(FI); normtype})
+
setMethod("updateNorm", "SelfNorm", function(normtype, L2, neighbor, biastype,
- Distr, V.comp, cent, stand, w)
- {Cv <- getInfV(L2deriv = L2, neighbor = neighbor,
+ Distr, V.comp, cent, stand, w){
+ Cv <- getInfV(L2deriv = L2, neighbor = neighbor,
biastype = biastype, Distr = Distr,
V.comp = V.comp, cent = cent, stand = stand, w = w)
QuadForm(normtype) <- PosSemDefSymmMatrix(distr::solve(Cv))
Modified: branches/robast-1.3/pkg/ROptEst/inst/NEWS
===================================================================
--- branches/robast-1.3/pkg/ROptEst/inst/NEWS 2021-04-01 16:47:39 UTC (rev 1232)
+++ branches/robast-1.3/pkg/ROptEst/inst/NEWS 2021-07-28 14:39:58 UTC (rev 1233)
@@ -11,6 +11,10 @@
version 1.3
#######################################
+under the hood:
++ in calls of form do.call(E, .....) we only use functions with one argument;
+ (in calls where a function was passed on as argument, this threw errors...)
+
#######################################
version 1.2.1
#######################################
Modified: branches/robast-1.3/pkg/ROptEst/inst/scripts/NormalLocationScaleModel.R
===================================================================
--- branches/robast-1.3/pkg/ROptEst/inst/scripts/NormalLocationScaleModel.R 2021-04-01 16:47:39 UTC (rev 1232)
+++ branches/robast-1.3/pkg/ROptEst/inst/scripts/NormalLocationScaleModel.R 2021-07-28 14:39:58 UTC (rev 1233)
@@ -21,7 +21,7 @@
N0.Rob1 <- InfRobModel(center = N0, neighbor = ContNeighborhood(radius = 0.5))
N0.Rob1 # show N0.Rob1
-## OBRE solution (ARE = .95)
+## OBRE solution (ARE = .95 acc. to Anscombe criterion)
system.time(N0.ICA <- optIC(model = N0.Rob1, risk = asAnscombe(), upper=NULL,lower=NULL, verbose=TRUE))
checkIC(N0.ICA)
Risks(N0.ICA)
@@ -28,6 +28,12 @@
plot(N0.ICA)
infoPlot(N0.ICA)
+## in two dimensions, other norms are possible
+## acc. to Rieder[94], we consider norms of type |x|^2 = x' C^{-1} x for some positive definite C
+## in particular, for C = FisherInfo^{-1}, we obtain information-standardization
+## in particular, for C = Cov(IC), we obtain self-standardization
+## notationally, we abbreviate this with index .i or index .s
+
system.time(N0.ICA.i <- optIC(model = N0.Rob1, risk = asAnscombe(eff=0.95, normtype=InfoNorm()), upper=NULL,lower=NULL, verbose=TRUE))
## MSE solution
@@ -37,6 +43,10 @@
plot(N0.IC1)
infoPlot(N0.IC1)
+## the infoPlot plots
+## + for IC psi the map x -> |psi(x)|^2 (where |y| denotes Euclidean norm of the two coordinates for location and scale)
+## + for each coordinate i of the IC psi the map x -> psi^{(i)}(x)^2/|psi(x)|^2
+
system.time(N0.IC1.i <- optIC(model = N0.Rob1, risk = asMSE(normtype=InfoNorm())))
checkIC(N0.IC1.i)
Risks(N0.IC1.i)
@@ -49,6 +59,7 @@
plot(N0.IC1.s)
infoPlot(N0.IC1.s)
+## the comparePlot compares the coordinates of the IFs
comparePlot(N0.IC1,N0.IC1.i,N0.IC1.s)
@@ -107,10 +118,13 @@
# a non-trivial trafo:
###############################################################################
+## here we use x[1]+qnorm(.95)*x[2] which, for x = (mu, sigma)
+## gives a parametric quantile in the normal distribution
+
tfct <- function(x){
nms0 <- c("mean","sd")
- nms <- "comb"
- fval0 <- x[1]+2*x[2]
+ nms <- "q95"
+ fval0 <- x[1]+qnorm(.95)*x[2]
names(fval0) <- nms
mat0 <- matrix(c(1,2), nrow = 1, dimnames = list(nms,nms0))
return(list(fval = fval0, mat = mat0))
@@ -217,6 +231,46 @@
estimate(ROest1)
estimate(ROest2)
+### plotting
+qKI <- qnorm(.975)
+n.chem <- length(chem)
+
+m.chem <- mean(chem)
+s.chem <- sd(chem)
+m.chem.rob <- estimate(ROest2)[1]
+s.chem.rob <- estimate(ROest2)[2]
+
+plot(chem, ylim=c(-10,35))
+
+## classical, non-robust estimate:
+abline(h=m.chem, col = "red")
+## confidence bands (for mu)
+abline(h=m.chem + qKI*s.chem/sqrt(n.chem), col = "red", lty = 2)
+abline(h=m.chem - qKI*s.chem/sqrt(n.chem), col = "red", lty = 2)
+## prediction bands (for all observations)
+abline(h=m.chem + qKI*s.chem, col = "red", lty = 3)
+abline(h=m.chem - qKI*s.chem, col = "red", lty = 3)
+
+## RMX-based estimate:
+abline(h=m.chem.rob, col = "darkblue")
+## confidence bands (for mu)
+abline(h=m.chem.rob + qKI*s.chem.rob/sqrt(n.chem), col = "darkblue", lty = 2)
+abline(h=m.chem.rob - qKI*s.chem.rob/sqrt(n.chem), col = "darkblue", lty = 2)
+## prediction bands (for all observations)
+abline(h=m.chem.rob + qKI*s.chem.rob, col = "darkblue", lty = 3)
+abline(h=m.chem.rob - qKI*s.chem.rob, col = "darkblue", lty = 3)
+
+legend("topright", legend = c("mu.MLE", "mu.rob", "95%-confidence band for mu - classic",
+ "95%-confidence band for mu - robust",
+ "95%-prediction bands for obs - classic",
+ "95%-prediction bands for obs - robust"), cex=0.7,
+ col=c("red", "darkblue"), lty=c(1,1,2,2,3,3))
+
+text(17,chem[17],"out", col="red",adj=1.05)
+text(13,chem[13],"out", col="blue",adj=1.05)
+### observation 17 sticks out as outlier in both classic and robust fit
+### observation 13 is only identified as outlier in robust fit
+
## confidence intervals
confint(ROest1, symmetricBias())
confint(ROest2, symmetricBias())
More information about the Robast-commits
mailing list