[Robast-commits] r324 - in branches/robast-0.7/pkg/ROptEst: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Jul 15 13:45:40 CEST 2009
Author: stamats
Date: 2009-07-15 13:45:40 +0200 (Wed, 15 Jul 2009)
New Revision: 324
Modified:
branches/robast-0.7/pkg/ROptEst/DESCRIPTION
branches/robast-0.7/pkg/ROptEst/R/getInfRobIC_asHampel.R
branches/robast-0.7/pkg/ROptEst/R/optIC.R
branches/robast-0.7/pkg/ROptEst/man/0ROptEst-package.Rd
branches/robast-0.7/pkg/ROptEst/man/getInfRobIC.Rd
branches/robast-0.7/pkg/ROptEst/man/optIC.Rd
Log:
added argument checkBounds, setting this to FALSE clearly speeds up computations for asHampel Risk
Modified: branches/robast-0.7/pkg/ROptEst/DESCRIPTION
===================================================================
--- branches/robast-0.7/pkg/ROptEst/DESCRIPTION 2009-07-15 10:06:41 UTC (rev 323)
+++ branches/robast-0.7/pkg/ROptEst/DESCRIPTION 2009-07-15 11:45:40 UTC (rev 324)
@@ -1,6 +1,6 @@
Package: ROptEst
Version: 0.7
-Date: 2009-04-22
+Date: 2009-07-15
Title: Optimally robust estimation
Description: Optimally robust estimation in general smoothly
parameterized models using S4 classes and methods.
Modified: branches/robast-0.7/pkg/ROptEst/R/getInfRobIC_asHampel.R
===================================================================
--- branches/robast-0.7/pkg/ROptEst/R/getInfRobIC_asHampel.R 2009-07-15 10:06:41 UTC (rev 323)
+++ branches/robast-0.7/pkg/ROptEst/R/getInfRobIC_asHampel.R 2009-07-15 11:45:40 UTC (rev 324)
@@ -5,63 +5,68 @@
risk = "asHampel",
neighbor = "UncondNeighborhood"),
function(L2deriv, risk, neighbor, symm, Finfo, trafo,
- upper, maxiter, tol, warn, noLow = FALSE, verbose = FALSE){
+ upper, maxiter, tol, warn, noLow = FALSE, verbose = FALSE,
+ checkBounds = TRUE){
biastype <- biastype(risk)
normtype <- normtype(risk)
A <- trafo / E(L2deriv, function(x){x^2})
b <- risk at bound
- bmax <- abs(as.vector(A))*max(abs(q(L2deriv)(0)), q(L2deriv)(1))
- if(b >= bmax){
- if(warn) cat("'b >= maximum asymptotic bias' => (classical) optimal IC\n",
- "in sense of Cramer-Rao bound is returned\n")
- res <- getInfRobIC(L2deriv = L2deriv, risk = asCov(),
- neighbor = neighbor, Finfo = Finfo, trafo = trafo,
- verbose = verbose)
- res <- c(res, list(biastype = biastype, normtype = NormType()))
- Cov <- res$risk$asCov
- r <- neighbor at radius
- res$risk$asBias <- list(value = b, biastype = biastype,
- normtype = normtype,
- neighbortype = class(neighbor))
- res$risk$asMSE <- list(value = Cov + r^2*b^2,
- r = r,
- at = neighbor)
- return(res)
- }
+ if(checkBounds){
+ bmax <- abs(as.vector(A))*max(abs(q(L2deriv)(0)), q(L2deriv)(1))
+ if(b >= bmax){
+ if(warn) cat("'b >= maximum asymptotic bias' => (classical) optimal IC\n",
+ "in sense of Cramer-Rao bound is returned\n")
+ res <- getInfRobIC(L2deriv = L2deriv, risk = asCov(),
+ neighbor = neighbor, Finfo = Finfo, trafo = trafo,
+ verbose = verbose)
+ res <- c(res, list(biastype = biastype, normtype = NormType()))
+ Cov <- res$risk$asCov
+ r <- neighbor at radius
+ res$risk$asBias <- list(value = b, biastype = biastype,
+ normtype = normtype,
+ neighbortype = class(neighbor))
+ res$risk$asMSE <- list(value = Cov + r^2*b^2,
+ r = r,
+ at = neighbor)
+ return(res)
+ }
- if(!noLow){
- res <- getInfRobIC(L2deriv = L2deriv, risk = asBias(biastype = biastype),
- neighbor = neighbor, symm = symm,
- trafo = trafo, maxiter = maxiter, tol = tol, Finfo = Finfo,
- warn = warn, verbose = verbose)
- bmin <- res$b
- cat("minimal bound:\t", bmin, "\n")
- }else bmin <- b/2
+ if(!noLow){
+ res <- getInfRobIC(L2deriv = L2deriv, risk = asBias(biastype = biastype),
+ neighbor = neighbor, symm = symm,
+ trafo = trafo, maxiter = maxiter, tol = tol, Finfo = Finfo,
+ warn = warn, verbose = verbose)
+ bmin <- res$b
+ cat("minimal bound:\t", bmin, "\n")
+ }else{
+ bmin <- b/2
+ }
- if(b <= bmin){
- if(warn) cat("'b <= minimum asymptotic bias'\n",
- "=> the minimum asymptotic bias (lower case) solution is returned\n")
- Risk <- list(asMSE = res$risk$asCov + neighbor at radius^2*bmin^2)
- res$risk <- c(Risk, res$risk)
- return(res)
+ if(b <= bmin){
+ if(warn) cat("'b <= minimum asymptotic bias'\n",
+ "=> the minimum asymptotic bias (lower case) solution is returned\n")
+ Risk <- list(asMSE = res$risk$asCov + neighbor at radius^2*bmin^2)
+ res$risk <- c(Risk, res$risk)
+ return(res)
+ }
+ # bmin <- getAsRisk(risk = asBias(biastype = biastype, normtype = normtype),
+ # L2deriv = L2deriv, neighbor = neighbor,
+ # biastype = biastype, trafo = trafo, Finfo = Finfo,
+ # warn = warn)$asBias
+ # if(b <= bmin){
+ # if(warn) cat("'b <= minimum asymptotic bias'\n",
+ # "=> the minimum asymptotic bias (lower case) solution is returned\n")
+ # res <- getInfRobIC(L2deriv = L2deriv, risk = asBias(biastype = biastype),
+ # neighbor = neighbor, symm = symm,
+ # trafo = trafo, maxiter = maxiter, tol = tol, Finfo = Finfo,
+ # warn = warn)
+ # Risk <- list(asMSE = res$risk$asCov + neighbor at radius^2*bmin^2)
+ # res$risk <- c(Risk, res$risk)
+ # return(res)
+ # }
}
-# bmin <- getAsRisk(risk = asBias(biastype = biastype, normtype = normtype),
-# L2deriv = L2deriv, neighbor = neighbor,
-# biastype = biastype, trafo = trafo, Finfo = Finfo,
-# warn = warn)$asBias
-# if(b <= bmin){
-# if(warn) cat("'b <= minimum asymptotic bias'\n",
-# "=> the minimum asymptotic bias (lower case) solution is returned\n")
-# res <- getInfRobIC(L2deriv = L2deriv, risk = asBias(biastype = biastype),
-# neighbor = neighbor, symm = symm,
-# trafo = trafo, maxiter = maxiter, tol = tol, Finfo = Finfo,
-# warn = warn)
-# Risk <- list(asMSE = res$risk$asCov + neighbor at radius^2*bmin^2)
-# res$risk <- c(Risk, res$risk)
-# return(res)
-# }
c0 <- b/as.vector(A)
if(is(symm, "SphericalSymmetry"))
S <- symm at SymmCenter == 0
@@ -128,7 +133,8 @@
neighbor = "ContNeighborhood"),
function(L2deriv, risk, neighbor, Distr, DistrSymm, L2derivSymm,
L2derivDistrSymm, Finfo, trafo, onesetLM = FALSE,
- z.start, A.start, upper, maxiter, tol, warn, verbose = FALSE){
+ z.start, A.start, upper, maxiter, tol, warn, verbose = FALSE,
+ checkBounds = TRUE){
biastype <- biastype(risk)
normtype <- normtype(risk)
@@ -142,57 +148,58 @@
if(is.null(z.start)) z.start <- numeric(ncol(trafo))
if(is.null(A.start)) A.start <- trafo
-
- ClassIC <- trafo %*% solve(Finfo) %*% L2deriv
- lower <- q(Distr)(getdistrOption("TruncQuantile"))
- upper <- q(Distr)(1-getdistrOption("TruncQuantile"))
- x <- seq(from = lower, to = upper, by = 0.01)
- bmax <- evalRandVar(ClassIC, as.matrix(x))^2
- bmax <- sqrt(max(colSums(bmax)))
b <- risk at bound
- cat("numerical approximation of maximal bound:\t", bmax, "\n")
- if(b >= bmax){
- if(warn) cat("'b >= maximum asymptotic bias' => (classical) optimal IC\n",
- "in sense of Cramer-Rao bound is returned\n")
- res <- getInfRobIC(L2deriv = L2deriv, risk = asCov(), neighbor = neighbor,
- Distr = Distr, Finfo = Finfo, trafo = trafo,
- QuadForm = std, verbose = verbose)
- res <- c(res, list(biastype = biastype, normtype = normtype))
- trAsCov <- sum(diag(std%*%res$risk$asCov));
- r <- neighbor at radius
- res$risk$trAsCov <- list(value = trAsCov, normtype = normtype)
- res$risk$asBias <- list(value = b, biastype = biastype,
- normtype = normtype,
- neighbortype = class(neighbor))
- res$risk$asMSE <- list(value = trAsCov + r^2*b^2,
- r = r,
- at = neighbor)
- return(res)
- }
- res <- getInfRobIC(L2deriv = L2deriv,
- risk = asBias(biastype = biastype, normtype = normtype),
- neighbor = neighbor, Distr = Distr, DistrSymm = DistrSymm,
- L2derivSymm = L2derivSymm, L2derivDistrSymm = L2derivDistrSymm,
- z.start = z.start, A.start = A.start, trafo = trafo,
- maxiter = maxiter, tol = tol, warn = warn, Finfo = Finfo,
- verbose = verbose)
- bmin <- res$b
+ if(checkBounds){
+ ClassIC <- trafo %*% solve(Finfo) %*% L2deriv
+ lower <- q(Distr)(getdistrOption("TruncQuantile"))
+ upper <- q(Distr)(1-getdistrOption("TruncQuantile"))
+ x <- seq(from = lower, to = upper, by = 0.01)
+ bmax <- evalRandVar(ClassIC, as.matrix(x))^2
+ bmax <- sqrt(max(colSums(bmax)))
+ cat("numerical approximation of maximal bound:\t", bmax, "\n")
+ if(b >= bmax){
+ if(warn) cat("'b >= maximum asymptotic bias' => (classical) optimal IC\n",
+ "in sense of Cramer-Rao bound is returned\n")
+ res <- getInfRobIC(L2deriv = L2deriv, risk = asCov(), neighbor = neighbor,
+ Distr = Distr, Finfo = Finfo, trafo = trafo,
+ QuadForm = std, verbose = verbose)
+ res <- c(res, list(biastype = biastype, normtype = normtype))
+ trAsCov <- sum(diag(std%*%res$risk$asCov));
+ r <- neighbor at radius
+ res$risk$trAsCov <- list(value = trAsCov, normtype = normtype)
+ res$risk$asBias <- list(value = b, biastype = biastype,
+ normtype = normtype,
+ neighbortype = class(neighbor))
+ res$risk$asMSE <- list(value = trAsCov + r^2*b^2,
+ r = r,
+ at = neighbor)
+ return(res)
+ }
- cat("minimal bound:\t", bmin, "\n")
- if(b <= bmin){
- if(warn) cat("'b <= minimum asymptotic bias'\n",
- "=> the minimum asymptotic bias (lower case) solution is returned\n")
+ res <- getInfRobIC(L2deriv = L2deriv,
+ risk = asBias(biastype = biastype, normtype = normtype),
+ neighbor = neighbor, Distr = Distr, DistrSymm = DistrSymm,
+ L2derivSymm = L2derivSymm, L2derivDistrSymm = L2derivDistrSymm,
+ z.start = z.start, A.start = A.start, trafo = trafo,
+ maxiter = maxiter, tol = tol, warn = warn, Finfo = Finfo,
+ verbose = verbose)
+ bmin <- res$b
- asMSE <- sum(diag(std%*%res$risk$asCov)) + neighbor at radius^2*bmin^2
- if(!is.null(res$risk$asMSE)) res$risk$asMSE <- asMSE
- else res$risk <- c(list(asMSE = asMSE), res$risk)
+ cat("minimal bound:\t", bmin, "\n")
+ if(b <= bmin){
+ if(warn) cat("'b <= minimum asymptotic bias'\n",
+ "=> the minimum asymptotic bias (lower case) solution is returned\n")
- return(res)
+ asMSE <- sum(diag(std%*%res$risk$asCov)) + neighbor at radius^2*bmin^2
+ if(!is.null(res$risk$asMSE)) res$risk$asMSE <- asMSE
+ else res$risk <- c(list(asMSE = asMSE), res$risk)
+
+ return(res)
+ }
}
- comp <- .getComp(L2deriv, DistrSymm, L2derivSymm,
- L2derivDistrSymm)
+ comp <- .getComp(L2deriv, DistrSymm, L2derivSymm, L2derivDistrSymm)
z.comp <- comp$"z.comp"
A.comp <- comp$"A.comp"
Modified: branches/robast-0.7/pkg/ROptEst/R/optIC.R
===================================================================
--- branches/robast-0.7/pkg/ROptEst/R/optIC.R 2009-07-15 10:06:41 UTC (rev 323)
+++ branches/robast-0.7/pkg/ROptEst/R/optIC.R 2009-07-15 11:45:40 UTC (rev 324)
@@ -4,7 +4,7 @@
setMethod("optIC", signature(model = "InfRobModel", risk = "asRisk"),
function(model, risk, z.start = NULL, A.start = NULL, upper = 1e4,
maxiter = 50, tol = .Machine$double.eps^0.4, warn = TRUE,
- noLow = FALSE, verbose = FALSE){
+ noLow = FALSE, verbose = FALSE, ...){
L2derivDim <- numberOfMaps(model at center@L2deriv)
ow <- options("warn")
on.exit(options(ow))
@@ -15,7 +15,7 @@
symm = model at center@L2derivDistrSymm[[1]],
Finfo = model at center@FisherInfo, trafo = trafo(model at center@param),
upper = upper, maxiter = maxiter, tol = tol, warn = warn,
- noLow = noLow, verbose = verbose)
+ noLow = noLow, verbose = verbose, ...)
res$info <- c("optIC", res$info)
res <- c(res, modifyIC = getModifyIC(L2FamIC = model at center,
neighbor = model at neighbor,
@@ -49,7 +49,7 @@
L2derivDistrSymm = L2derivDistrSymm, Finfo = model at center@FisherInfo,
trafo = trafo(model at center@param), z.start = z.start, A.start = A.start,
upper = upper, maxiter = maxiter, tol = tol, warn = warn,
- verbose = verbose)
+ verbose = verbose, ...)
options(ow)
res$info <- c("optIC", res$info)
res <- c(res, modifyIC = getModifyIC(L2FamIC = model at center,
Modified: branches/robast-0.7/pkg/ROptEst/man/0ROptEst-package.Rd
===================================================================
--- branches/robast-0.7/pkg/ROptEst/man/0ROptEst-package.Rd 2009-07-15 10:06:41 UTC (rev 323)
+++ branches/robast-0.7/pkg/ROptEst/man/0ROptEst-package.Rd 2009-07-15 11:45:40 UTC (rev 324)
@@ -13,7 +13,7 @@
\tabular{ll}{
Package: \tab ROptEst\cr
Version: \tab 0.7 \cr
-Date: \tab 2009-04-22 \cr
+Date: \tab 2009-07-15 \cr
Depends: \tab R(>= 2.7.0), methods, distr(>= 2.0), distrEx(>= 2.0),distrMod(>= 2.0), RandVar(>= 0.6.4), RobAStBase\cr
LazyLoad: \tab yes\cr
License: \tab LGPL-3\cr
Modified: branches/robast-0.7/pkg/ROptEst/man/getInfRobIC.Rd
===================================================================
--- branches/robast-0.7/pkg/ROptEst/man/getInfRobIC.Rd 2009-07-15 10:06:41 UTC (rev 323)
+++ branches/robast-0.7/pkg/ROptEst/man/getInfRobIC.Rd 2009-07-15 11:45:40 UTC (rev 324)
@@ -35,10 +35,11 @@
L2derivDistrSymm, z.start, A.start, Finfo, trafo, maxiter, tol, warn, verbose = FALSE, ...)
\S4method{getInfRobIC}{UnivariateDistribution,asHampel,UncondNeighborhood}(L2deriv, risk, neighbor, symm, Finfo, trafo,
- upper, maxiter, tol, warn, noLow = FALSE, verbose = FALSE)
+ upper, maxiter, tol, warn, noLow = FALSE, verbose = FALSE, checkBounds = TRUE)
\S4method{getInfRobIC}{RealRandVariable,asHampel,ContNeighborhood}(L2deriv, risk, neighbor, Distr, DistrSymm, L2derivSymm,
- L2derivDistrSymm, Finfo, trafo, onesetLM = FALSE, z.start, A.start, upper, maxiter, tol, warn, verbose = FALSE)
+ L2derivDistrSymm, Finfo, trafo, onesetLM = FALSE, z.start, A.start, upper, maxiter, tol, warn,
+ verbose = FALSE, checkBounds = TRUE)
\S4method{getInfRobIC}{UnivariateDistribution,asGRisk,UncondNeighborhood}(L2deriv, risk, neighbor, symm, Finfo, trafo,
upper, maxiter, tol, warn, noLow = FALSE, verbose = FALSE)
@@ -74,6 +75,8 @@
\code{PosSemDefSymmMatrix} for use of different
(standardizing) norm }
\item{verbose}{ logical: if \code{TRUE}, some messages are printed }
+ \item{checkBounds}{ logical: if \code{TRUE}, minimal and maximal clipping bound are
+ computed to check if a valid bound was specified. }
}
%\details{}
\value{The optimally robust IC is computed.}
Modified: branches/robast-0.7/pkg/ROptEst/man/optIC.Rd
===================================================================
--- branches/robast-0.7/pkg/ROptEst/man/optIC.Rd 2009-07-15 10:06:41 UTC (rev 323)
+++ branches/robast-0.7/pkg/ROptEst/man/optIC.Rd 2009-07-15 11:45:40 UTC (rev 324)
@@ -15,7 +15,7 @@
\S4method{optIC}{InfRobModel,asRisk}(model, risk, z.start = NULL, A.start = NULL,
upper = 1e4, maxiter = 50,
tol = .Machine$double.eps^0.4, warn = TRUE,
- noLow = FALSE, verbose = FALSE)
+ noLow = FALSE, verbose = FALSE, ...)
\S4method{optIC}{InfRobModel,asUnOvShoot}(model, risk, upper = 1e4, maxiter = 50,
tol = .Machine$double.eps^0.4, warn = TRUE)
More information about the Robast-commits
mailing list