[Robast-commits] r412 - in branches/robast-0.8/pkg/ROptEst: R inst/scripts
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Aug 29 12:27:58 CEST 2010
Author: ruckdeschel
Date: 2010-08-29 12:27:58 +0200 (Sun, 29 Aug 2010)
New Revision: 412
Added:
branches/robast-0.8/pkg/ROptEst/inst/scripts/MBRE.R
Modified:
branches/robast-0.8/pkg/ROptEst/R/LowerCaseMultivariate.R
branches/robast-0.8/pkg/ROptEst/R/getInfLM.R
branches/robast-0.8/pkg/ROptEst/R/getInfRobIC_asAnscombe.R
branches/robast-0.8/pkg/ROptEst/R/getInfRobIC_asBias.R
Log:
ROptEst: fixed an issue with lower case solution in information/self-standardization;
now Anscombe-Risk also works with information standardization ( but is not suited for self-standardization!)
Modified: branches/robast-0.8/pkg/ROptEst/R/LowerCaseMultivariate.R
===================================================================
--- branches/robast-0.8/pkg/ROptEst/R/LowerCaseMultivariate.R 2010-08-27 11:13:15 UTC (rev 411)
+++ branches/robast-0.8/pkg/ROptEst/R/LowerCaseMultivariate.R 2010-08-29 10:27:58 UTC (rev 412)
@@ -15,12 +15,13 @@
if(is.null(z.comp))
z.comp <- rep(TRUE, nrow(trafo))
+ force(normtype)
lA.comp <- sum(A.comp)
- abs.fct <- function(x, L2, stand, cent, normtype){
+ abs.fct <- function(x, L2, stand, cent, normtype.0){
X <- evalRandVar(L2, as.matrix(x))[,,1] - cent
Y <- stand %*% X
- return(fct(normtype)(Y))
+ return(fct(normtype.0)(Y))
}
itermin <- 0
@@ -29,7 +30,10 @@
p <- nrow(trafo)
k <- ncol(trafo)
A <- matrix(0, ncol = k, nrow = p)
+
A[A.comp] <- param[1:lA.comp]
+ A.max <- max(abs(A.comp))
+ A <- A/A.max
z <- numeric(k)
z[z.comp] <- param[(lA.comp+1):length(param)]
@@ -48,32 +52,46 @@
neighbor = neighbor, biastype = biastype,
Distr = Distr, V.comp = A.comp,
cent = z, stand = A, w = w0)
+ weight(w0) <- minbiasweight(w0, neighbor = neighbor,
+ biastype = biastype,
+ normW = normtype)
+
+ w <<- w0
}
E1 <- E(object = Distr, fun = abs.fct, L2 = L2deriv, stand = A,
- cent = z, normtype = normtype, useApply = FALSE)
- stA <- if (is(normtype,"QFnorm"))
+ cent = z, normtype.0 = normtype, useApply = FALSE)
+ stA <- if (is(normtype,"QFNorm"))
QuadForm(normtype)%*%A else A
+# erg <- E1/sum(diag(stA %*% t(trafo)))
+# print(list(A,stA,E1))
erg <- E1/sum(diag(stA %*% t(trafo)))
clip(w0) <- 1/erg
w <<- w0
if(verbose && itermin %% 15 == 1){
- cat("trying to find lower case solution;\n")
+# if(verbose && itermin %% 2 == 1){
+ cat("trying to find lower case solution;\n")
cat("current Lagrange Multiplier value:\n")
print(list(A=A, z=z,erg=erg))
}
-
+
return(erg)
}
A.vec <- as.vector(A.start[A.comp])
- p.vec <- c(A.vec, z.start[z.comp])
- force(normtype)
+ A.max <- max(abs(A.vec))
+ A.vec <- A.vec/A.max
+ p.vec <- c(A.vec, z.start[z.comp]/A.max)
+
erg <- optim(p.vec, bmin.fct, method = "Nelder-Mead",
control = list(reltol = tol, maxit = 100*maxiter),
L2deriv = L2deriv, Distr = Distr, trafo = trafo)
+ A.max <- max(abs(stand(w)))
+ stand(w) <- stand(w)/A.max
+ weight(w) <- minbiasweight(w, neighbor = neighbor,
+ biastype = biastype,
+ normW = normtype)
-
return(list(erg=erg, w=w, normtype = normtype, z.comp = z.comp, itermin = itermin))
}
Modified: branches/robast-0.8/pkg/ROptEst/R/getInfLM.R
===================================================================
--- branches/robast-0.8/pkg/ROptEst/R/getInfLM.R 2010-08-27 11:13:15 UTC (rev 411)
+++ branches/robast-0.8/pkg/ROptEst/R/getInfLM.R 2010-08-29 10:27:58 UTC (rev 412)
@@ -75,7 +75,7 @@
if(verbose && iter>1 && iter < maxiter && iter%%5 == 1){
cat("current precision in IC algo:\t", prec, "\n")
- print(round(c(A=A,a=a),3))
+ print(round(c(A=prec[1],a=prec[2]),3))
}
if(prec < tol) break
if(iter > maxiter){
Modified: branches/robast-0.8/pkg/ROptEst/R/getInfRobIC_asAnscombe.R
===================================================================
--- branches/robast-0.8/pkg/ROptEst/R/getInfRobIC_asAnscombe.R 2010-08-27 11:13:15 UTC (rev 411)
+++ branches/robast-0.8/pkg/ROptEst/R/getInfRobIC_asAnscombe.R 2010-08-29 10:27:58 UTC (rev 412)
@@ -134,7 +134,6 @@
}
std <- if(is(normtype,"QFNorm"))
QuadForm(normtype) else diag(p)
- print(std)
trV.ML <- sum(diag(std%*%FI0))
if(is.null(upper))
@@ -147,14 +146,14 @@
L2derivSymm = L2derivSymm,
L2derivDistrSymm = L2derivDistrSymm,
z.start = z.start, A.start = A.start,
- trafo = trafo, maxiter = maxi,
- tol = toli,
+ trafo = trafo, maxiter = maxiter,
+ tol = tol,
warn = FALSE, Finfo = Finfo,
QuadForm = std, verbose = verbose)
if(is.null(lower)||(lower< lowBerg$b))
{lower <- lowBerg$b
- print(lowBerg$risk$asAnscombe)
+# print(lowBerg$risk$asAnscombe)
f.low <- lowBerg$risk$asAnscombe - eff
} else {
risk.b <- asHampel(bound = lower, biastype = biastype,
Modified: branches/robast-0.8/pkg/ROptEst/R/getInfRobIC_asBias.R
===================================================================
--- branches/robast-0.8/pkg/ROptEst/R/getInfRobIC_asBias.R 2010-08-27 11:13:15 UTC (rev 411)
+++ branches/robast-0.8/pkg/ROptEst/R/getInfRobIC_asBias.R 2010-08-29 10:27:58 UTC (rev 412)
@@ -40,16 +40,18 @@
stop("Not yet implemented.")
normtype <- normtype(risk)
+
+# if(FALSE){
if(is(normtype,"SelfNorm")){
warntxt <- paste(gettext(
"Using self-standardization, there are problems with the existence\n"
),gettext(
"of a minmax Bias IC. Instead we compute the optimal MSE-solution\n"
),gettext(
- "to a large radius (r = 10)\n"
+ "to a large radius (r = 15)\n"
))
if(warn) cat(warntxt)
- neighbor at radius <- 10
+ neighbor at radius <- 15
res <- getInfRobIC(L2deriv = L2deriv,
risk = asMSE(normtype = normtype),
neighbor = neighbor, Distr = Distr,
@@ -58,6 +60,19 @@
trafo = trafo, onesetLM = FALSE, z.start = z.start,
A.start = A.start, upper = 1e4, maxiter = maxiter,
tol = tol, warn = warn, verbose = verbose)
+
+ A.max <- max(abs(res$A))
+ res$A <- res$A/A.max
+ res$a <- res$a/A.max
+ w <- res$w
+ A.max.w <- max(abs(stand(w)))
+ stand(w) <- stand(w)/A.max.w
+ normtype <- res$normtype
+ weight(w) <- minbiasweight(w, neighbor = neighbor,
+ biastype = biastype(risk),
+ normW = normtype)
+ res$w <- w
+
res$risk$asBias <- list(value = sqrt(nrow(trafo)),
biastype = symmetricBias(),
normtype = normtype,
@@ -65,9 +80,10 @@
remark = gettext("value is only a bound"))
return(res)
}
+# }
FI <- solve(trafo%*%solve(Finfo)%*%t(trafo))
- if(is(normtype,"InfoNorm"))
+ if(is(normtype,"QFNorm"))
{QuadForm(normtype) <- PosSemDefSymmMatrix(FI);
normtype(risk) <- normtype}
@@ -191,6 +207,8 @@
k <- ncol(trafo)
A <- matrix(0, ncol=k, nrow=p)
A[DA.comp] <- matrix(param[1:lA.comp], ncol=k, nrow=p)
+ A.max <- max(abs(A))
+ A <- A/A.max
z <- numeric(k)
z[z.comp] <- param[(lA.comp+1):length(param)]
a <- as.vector(A %*% z)
Added: branches/robast-0.8/pkg/ROptEst/inst/scripts/MBRE.R
===================================================================
--- branches/robast-0.8/pkg/ROptEst/inst/scripts/MBRE.R (rev 0)
+++ branches/robast-0.8/pkg/ROptEst/inst/scripts/MBRE.R 2010-08-29 10:27:58 UTC (rev 412)
@@ -0,0 +1,21 @@
+if(FALSE){
+require(ROptEst)
+options("newDevice"=TRUE)
+
+## generates normal location and scale family with mean = -2 and sd = 3
+ ### 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)
+}
\ No newline at end of file
More information about the Robast-commits
mailing list