[Robast-commits] r158 - in branches/robast-0.6/pkg/ROptEst: R inst/scripts
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Aug 10 18:39:02 CEST 2008
Author: stamats
Date: 2008-08-10 18:39:02 +0200 (Sun, 10 Aug 2008)
New Revision: 158
Modified:
branches/robast-0.6/pkg/ROptEst/R/getIneffDiff.R
branches/robast-0.6/pkg/ROptEst/R/getInfRobIC_asGRisk.R
branches/robast-0.6/pkg/ROptEst/R/radiusMinimaxIC.R
branches/robast-0.6/pkg/ROptEst/inst/scripts/ExponentialScaleModel.R
branches/robast-0.6/pkg/ROptEst/inst/scripts/NormalLocationScaleModel.R
Log:
adaption to new implementations, some modifications to speed up computations
Modified: branches/robast-0.6/pkg/ROptEst/R/getIneffDiff.R
===================================================================
--- branches/robast-0.6/pkg/ROptEst/R/getIneffDiff.R 2008-08-08 17:53:31 UTC (rev 157)
+++ branches/robast-0.6/pkg/ROptEst/R/getIneffDiff.R 2008-08-10 16:39:02 UTC (rev 158)
@@ -75,8 +75,8 @@
ineffUp <- if(upRad == Inf) biasUp^2/upRisk else
(p+biasUp^2*upRad^2)/upRisk
}else{
- ineffLo <- (sum(diag(std%*%res$A%*%t(trafo))) -
- biasLo^2*(radius^2-loRad^2))/loRisk
+ ineffLo <- (sum(diag(std%*%res$A%*%t(trafo))) -
+ biasLo^2*(radius^2-loRad^2))/loRisk
if(upRad == Inf)
ineffUp <- biasUp^2/upRisk
else
Modified: branches/robast-0.6/pkg/ROptEst/R/getInfRobIC_asGRisk.R
===================================================================
--- branches/robast-0.6/pkg/ROptEst/R/getInfRobIC_asGRisk.R 2008-08-08 17:53:31 UTC (rev 157)
+++ branches/robast-0.6/pkg/ROptEst/R/getInfRobIC_asGRisk.R 2008-08-10 16:39:02 UTC (rev 158)
@@ -36,6 +36,7 @@
else
S <- FALSE
+ prec <- 1
repeat{
iter <- iter + 1
z.old <- z
@@ -76,9 +77,18 @@
clip = c0, cent = z, symm = S, trafo = trafo, tol.z = tol)
# cat("c0:\t", c0, "c0.old:\t", c0.old, "z:\t", z, "z.old:\t", z.old, "\n")
if(S) break
- if(max(abs(z - z.old), abs(c0-c0.old)) < tol) break
+
+ prec.old <- prec
+ prec <- max(abs(z - z.old), abs(c0-c0.old))
+ if(verbose)
+ cat("current precision in IC algo:\t", prec, "\n")
+ if(prec < tol) break
+ if(abs(prec.old - prec) < 1e-10){
+ cat("algorithm did not converge!\n", "achieved precision:\t", prec, "\n")
+ break
+ }
if(iter > maxiter){
- cat("maximum iterations reached!\n", "achieved precision:\t", abs(c0 - c0.old), "\n")
+ cat("maximum iterations reached!\n", "achieved precision:\t", prec, "\n")
break
}
}
@@ -177,6 +187,7 @@
A <- A.start
b <- 0
iter <- 0
+ prec <- 1
repeat{
iter <- iter + 1
z.old <- z
@@ -249,10 +260,15 @@
Distr = Distr, V.comp = A.comp, cent = as.vector(A %*% z),
stand = A, w = w)}
+ prec.old <- prec
prec <- max(abs(b-b.old), max(abs(A-A.old)), max(abs(z-z.old)))
if(verbose)
cat("current precision in IC algo:\t", prec, "\n")
if(prec < tol) break
+ if(abs(prec.old - prec) < 1e-10){
+ cat("algorithm did not converge!\n", "achieved precision:\t", prec, "\n")
+ break
+ }
if(iter > maxiter){
cat("maximum iterations reached!\n", "achieved precision:\t", prec, "\n")
break
@@ -274,7 +290,7 @@
biastype = biastype, Distr = Distr,
V.comp = A.comp, cent = a,
stand = A, w = w)
-
+
QF <- if(is(normtype,"QFNorm")) QuadForm(normtype) else diag(nrow(A))
trAsCov <- sum(diag(QF%*%Cov))
Modified: branches/robast-0.6/pkg/ROptEst/R/radiusMinimaxIC.R
===================================================================
--- branches/robast-0.6/pkg/ROptEst/R/radiusMinimaxIC.R 2008-08-08 17:53:31 UTC (rev 157)
+++ branches/robast-0.6/pkg/ROptEst/R/radiusMinimaxIC.R 2008-08-10 16:39:02 UTC (rev 158)
@@ -6,7 +6,7 @@
neighbor = "UncondNeighborhood",
risk = "asGRisk"),
function(L2Fam, neighbor, risk, loRad, upRad, z.start = NULL,
- A.start = NULL, upper = 1e5, maxiter = 100,
+ A.start = NULL, upper = 1e5, maxiter = 50,
tol = .Machine$double.eps^0.4, warn = FALSE, verbose = FALSE){
if(length(loRad) != 1)
stop("'loRad' is not of length == 1")
@@ -183,7 +183,7 @@
z.start = z.start, A.start = A.start, upper.b = upper.b, risk = risk,
loRad = loRad, upRad = upRad, loRisk = loRisk, upRisk = upRisk,
eps = tol, MaxIter = maxiter, warn = warn,
- loNorm = loNorm, upNorm = upNorm)$root
+ loNorm = loNorm, upNorm = upNorm, verbose = verbose)$root
neighbor at radius <- leastFavR
res <- getInfRobIC(L2deriv = L2deriv, neighbor = neighbor, risk = risk,
Distr = L2Fam at distribution, DistrSymm = L2Fam at distrSymm,
Modified: branches/robast-0.6/pkg/ROptEst/inst/scripts/ExponentialScaleModel.R
===================================================================
--- branches/robast-0.6/pkg/ROptEst/inst/scripts/ExponentialScaleModel.R 2008-08-08 17:53:31 UTC (rev 157)
+++ branches/robast-0.6/pkg/ROptEst/inst/scripts/ExponentialScaleModel.R 2008-08-10 16:39:02 UTC (rev 158)
@@ -3,8 +3,8 @@
###############################################################################
require(ROptEst)
-## generates Exponential Scale Family with rate = 2
-E1 <- ExpScaleFamily(rate = 2)
+## generates Exponential Scale Family with scale = 0.5 (rate = 2)
+E1 <- ExpScaleFamily(scale = 0.5)
E1 # show E1
plot(E1) # plot of Exp(rate = 1) and L_2 derivative
checkL2deriv(E1)
@@ -74,13 +74,13 @@
E1.x <- rexp(1e2, rate=(1-ind)*2+ind*10)
## 2. Kolmogorov(-Smirnov) minimum distance estimator
-(E1.est0 <- MDEstimator(x=E1.x, ExpScaleFamily(), interval = c(0, 10)))
+(E1.est0 <- MDEstimator(x=E1.x, ExpScaleFamily()))
## 3. one-step estimation: radius known
-E1.Rob3 <- InfRobModel(center=ExpScaleFamily(rate=1/E1.est0$estimate),
+E1.Rob3 <- InfRobModel(center=ExpScaleFamily(scale=estimate(E1.est0)),
neighbor=ContNeighborhood(radius=0.5))
E1.IC9 <- optIC(model=E1.Rob3, risk=asMSE())
-(E1.est1 <- oneStepEstimator(E1.x, IC=E1.IC9, start=E1.est0$estimate))
+(E1.est1 <- oneStepEstimator(E1.x, IC=E1.IC9, start=E1.est0))
## 4. one-step estimation: radius interval
E1.IC10 <- radiusMinimaxIC(L2Fam=ExpScaleFamily(rate=1/E1.est0$estimate),
Modified: branches/robast-0.6/pkg/ROptEst/inst/scripts/NormalLocationScaleModel.R
===================================================================
--- branches/robast-0.6/pkg/ROptEst/inst/scripts/NormalLocationScaleModel.R 2008-08-08 17:53:31 UTC (rev 157)
+++ branches/robast-0.6/pkg/ROptEst/inst/scripts/NormalLocationScaleModel.R 2008-08-10 16:39:02 UTC (rev 158)
@@ -3,10 +3,10 @@
###############################################################################
require(ROptEst)
-## generates normal location and scale family with mean = 0 and sd = 1
-N0 <- NormLocationScaleFamily(mean=0, sd=1)
+## generates normal location and scale family with mean = -2 and sd = 3
+N0 <- NormLocationScaleFamily(mean=-2, sd=3)
N0 # show G0
-plot(N0) # plot of Norm(mean = 0, sd = 1) and L_2 derivative
+plot(N0) # plot of Norm(mean = -2, sd = 3) and L_2 derivative
checkL2deriv(N0)
## classical optimal IC
@@ -64,7 +64,7 @@
## radius minimax IC
## (may take quite some time!)
-(N0.IC4 <- radiusMinimaxIC(L2Fam=N0, neighbor=ContNeighborhood(),
+system.time(N0.IC4 <- radiusMinimaxIC(L2Fam=N0, neighbor=ContNeighborhood(),
risk=asMSE(), loRad=0, upRad=Inf))
checkIC(N0.IC4)
Risks(N0.IC4)
More information about the Robast-commits
mailing list