[Robast-commits] r82 - in pkg: . ROptEst/R ROptEst/chm ROptEst/inst/scripts ROptEst/man ROptEst.Rcheck RobAStBase RobAStBase/R RobAStBase/chm RobAStBase/man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Mar 29 05:54:45 CET 2008
Author: ruckdeschel
Date: 2008-03-29 05:54:44 +0100 (Sat, 29 Mar 2008)
New Revision: 82
Added:
pkg/ROptEst.Rcheck/
pkg/ROptEst.Rcheck/ROptEst/
pkg/ROptEst/R/getComp.R
Modified:
pkg/ROptEst/R/AllGeneric.R
pkg/ROptEst/R/LowerCaseMultivariate.R
pkg/ROptEst/R/getAsRisk.R
pkg/ROptEst/R/getIneffDiff.R
pkg/ROptEst/R/getInfRobIC_asBias.R
pkg/ROptEst/R/getInfRobIC_asGRisk.R
pkg/ROptEst/R/getInfRobIC_asHampel.R
pkg/ROptEst/R/getInfRobIC_asUnOvShoot.R
pkg/ROptEst/R/getInfV.R
pkg/ROptEst/R/getRiskIC.R
pkg/ROptEst/R/leastFavorableRadius.R
pkg/ROptEst/R/lowerCaseRadius.R
pkg/ROptEst/R/radiusMinimaxIC.R
pkg/ROptEst/R/updateNorm.R
pkg/ROptEst/chm/00Index.html
pkg/ROptEst/chm/ROptEst.chm
pkg/ROptEst/chm/ROptEst.toc
pkg/ROptEst/chm/getAsRisk.html
pkg/ROptEst/chm/getInfRobIC.html
pkg/ROptEst/chm/minmaxBias.html
pkg/ROptEst/chm/updateNorm-methods.html
pkg/ROptEst/inst/scripts/NormalLocationScaleModel.R
pkg/ROptEst/inst/scripts/PoissonModel.R
pkg/ROptEst/man/getAsRisk.Rd
pkg/ROptEst/man/getInfRobIC.Rd
pkg/ROptEst/man/minmaxBias.Rd
pkg/ROptEst/man/updateNorm-methods.Rd
pkg/RobAStBase/NAMESPACE
pkg/RobAStBase/R/TotalVarIC.R
pkg/RobAStBase/R/Weights.R
pkg/RobAStBase/R/generateICfct.R
pkg/RobAStBase/R/getRiskIC.R
pkg/RobAStBase/R/infoPlot.R
pkg/RobAStBase/chm/RobAStBase.chm
pkg/RobAStBase/chm/getweight.html
pkg/RobAStBase/man/getweight.Rd
Log:
******************************************
********** Major release step ************
******************************************
major steps:
* self- and info-standardization;
* onesided, asymmetric ICs
* everything with weight functions
... everything compiles now!
except for minor warnings:
+ in distrMod we use a class union with "function", which
seems to be a problem
+ in RobAStBase, there is some warning about the replacement function "weight<-" which I do not understand
+ in ROptEst still some consistency checks with prior versions are due
+ it seems that lower Case solutions and minimax radius solutions need not exist for self-standardized ICs ...
NEXT Step[s]: control objects for easier interfaces...
Modified: pkg/ROptEst/R/AllGeneric.R
===================================================================
--- pkg/ROptEst/R/AllGeneric.R 2008-03-28 11:36:49 UTC (rev 81)
+++ pkg/ROptEst/R/AllGeneric.R 2008-03-29 04:54:44 UTC (rev 82)
@@ -71,3 +71,7 @@
if(!isGeneric("updateNorm")){
setGeneric("updateNorm", function(normtype, ...) standardGeneric("updateNorm"))
}
+if(!isGeneric("getBiasIC")){
+ setGeneric("getBiasIC",
+ function(L2deriv, neighbor, biastype, ...) standardGeneric("minmaxBias"))
+}
Modified: pkg/ROptEst/R/LowerCaseMultivariate.R
===================================================================
--- pkg/ROptEst/R/LowerCaseMultivariate.R 2008-03-28 11:36:49 UTC (rev 81)
+++ pkg/ROptEst/R/LowerCaseMultivariate.R 2008-03-29 04:54:44 UTC (rev 82)
@@ -1,65 +1,67 @@
.LowerCaseMultivariate <- function(L2deriv, neighbor, biastype,
- normtype, Distr, L2derivDistrSymm, trafo, z.start,
- A.start, maxiter, tol){
+ normtype, Distr, trafo, z.start,
+ A.start, z.comp, A.comp, maxiter, tol){
w <- new("HampelWeight")
-
if(is.null(z.start)) z.start <- numeric(ncol(trafo))
if(is.null(A.start)) A.start <- trafo
+ if(is.null(A.comp))
+ A.comp <- matrix(TRUE, nrow = nrow(trafo), ncol = ncol(trafo))
+ if(is.null(z.comp))
+ z.comp <- rep(TRUE, nrow(trafo))
+ lA.comp <- sum(A.comp)
+
abs.fct <- function(x, L2, stand, cent, normtype){
X <- evalRandVar(L2, as.matrix(x))[,,1] - cent
Y <- stand %*% X
return(fct(normtype)(Y))
}
- bmin.fct <- function(param, L2deriv, Distr, trafo, z.comp){
+ bmin.fct <- function(param, L2deriv, Distr, trafo){
p <- nrow(trafo)
k <- ncol(trafo)
- A <- matrix(param[1:(p*k)], ncol=k, nrow=p)
+ A <- matrix(0, ncol = k, nrow = p)
+ A[A.comp] <- param[1:lA.comp]
z <- numeric(k)
- z[z.comp] <- param[(p*k+1):length(param)]
-
+ z[z.comp] <- param[(lA.comp+1):length(param)]
+
+# if(is(normtype,"SelfNorm"))
+# A <- A/max(A)
+
+ w0 <- w
+ cent(w0) <- z
+ stand(w0) <- A
+ weight(w0) <- minbiasweight(w0, neighbor = neighbor,
+ biastype = biastype,
+ normW = normtype)
+ w <<- w0
if (is(normtype,"SelfNorm")){
- w0 <- w
- cent(w0) <- z
- stand(w0) <- A
- weight(w0) <- minbiasweight(w0, neighbor = neighbor,
- biastype = biastype,
- normtype = normtype)
- w <<- w0
normtype <<- updateNorm(normtype = normtype, L2 = L2deriv,
neighbor = neighbor, biastype = biastype,
- Distr = Distr, V.comp = matrix(TRUE, p,p),
- cent = z, stand = A, w = w, ...)
-
+ Distr = Distr, V.comp = A.comp,
+ cent = z, stand = A, w = w0)
}
E1 <- E(object = Distr, fun = abs.fct, L2 = L2deriv, stand = A,
cent = z, normtype = normtype, useApply = FALSE)
stA <- if (is(normtype,"QFnorm"))
QuadForm(normtype)%*%A else A
-
- return(E1/sum(diag(stA %*% t(trafo))))
+ erg <- E1/sum(diag(stA %*% t(trafo)))
+ clip(w0) <- 1/erg
+ w <<- w0
+ return(erg)
}
- nrvalues <- length(L2deriv)
- z.comp <- rep(TRUE, nrvalues)
- for(i in 1:nrvalues)
- if(is(L2derivDistrSymm[[i]], "SphericalSymmetry"))
- if(L2derivDistrSymm[[i]]@SymmCenter == 0)
- z.comp[i] <- FALSE
-
- A.vec <- as.vector(A.start)
+ A.vec <- as.vector(A.start[A.comp])
+ p.vec <- c(A.vec, z.start[z.comp])
force(normtype)
-
- erg <- optim(c(A.vec, z.start[z.comp]), bmin.fct, method = "Nelder-Mead",
+ erg <- optim(p.vec, bmin.fct, method = "Nelder-Mead",
control = list(reltol = tol, maxit = 100*maxiter),
- L2deriv = L2deriv, Distr = Distr, trafo = trafo, z.comp = z.comp)
+ L2deriv = L2deriv, Distr = Distr, trafo = trafo)
- return(list(erg=erg, w=w, normtype = normtype))
+ return(list(erg=erg, w=w, normtype = normtype, z.comp = z.comp))
}
-
Modified: pkg/ROptEst/R/getAsRisk.R
===================================================================
--- pkg/ROptEst/R/getAsRisk.R 2008-03-28 11:36:49 UTC (rev 81)
+++ pkg/ROptEst/R/getAsRisk.R 2008-03-29 04:54:44 UTC (rev 82)
@@ -57,20 +57,26 @@
L2deriv = "RealRandVariable",
neighbor = "ContNeighborhood",
biastype = "ANY"),
- function(risk, L2deriv, neighbor, biastype, Distr,
+ function(risk, L2deriv, neighbor, biastype, Distr, DistrSymm, L2derivSymm,
L2derivDistrSymm, trafo, z.start, A.start, maxiter, tol){
normtype <- normtype(risk)
biastype <- biastype(risk)
+
+ comp <- .getComp(L2deriv, DistrSymm, L2derivSymm, L2derivDistrSymm)
+ z.comp <- comp$"z.comp"
+ A.comp <- comp$"A.comp"
+
eerg <- .LowerCaseMultivariate(L2deriv, neighbor, biastype,
- normtype, Distr, L2derivDistrSymm, trafo, z.start,
- A.start, maxiter, tol)
+ normtype, Distr, trafo, z.start,
+ A.start, z.comp = z.comp, A.comp = A.comp, maxiter, tol)
erg <- eerg$erg
bias <- 1/erg$value
return(list(asBias = bias))
})
+
###############################################################################
## asymptotic covariance
###############################################################################
Added: pkg/ROptEst/R/getComp.R
===================================================================
--- pkg/ROptEst/R/getComp.R (rev 0)
+++ pkg/ROptEst/R/getComp.R 2008-03-29 04:54:44 UTC (rev 82)
@@ -0,0 +1,24 @@
+.getComp <- function(L2deriv, DistrSymm, L2derivSymm,
+ L2derivDistrSymm){
+ nrvalues <- length(L2deriv)
+ z.comp <- rep(TRUE, nrvalues)
+ A.comp <- matrix(TRUE, ncol = nrvalues, nrow = nrvalues)
+ for(i in 1:nrvalues){
+ if(is(L2derivDistrSymm[[i]], "SphericalSymmetry"))
+ if(L2derivDistrSymm[[i]]@SymmCenter == 0)
+ z.comp[i] <- FALSE
+ }
+ for(i in 1:(nrvalues-1))
+ for(j in (i+1):nrvalues){
+ if(is(DistrSymm, "SphericalSymmetry")){
+ if((is(L2derivSymm[[i]], "OddSymmetric") & is(L2derivSymm[[j]], "EvenSymmetric"))
+ | (is(L2derivSymm[[j]], "OddSymmetric") & is(L2derivSymm[[i]], "EvenSymmetric")))
+ if((L2derivSymm[[i]]@SymmCenter == L2derivSymm[[j]]@SymmCenter)
+ & (L2derivSymm[[i]]@SymmCenter == DistrSymm at SymmCenter))
+ A.comp[i,j] <- FALSE
+ }
+ }
+ A.comp[col(A.comp) < row(A.comp)] <- A.comp[col(A.comp) > row(A.comp)]
+ return(list(A.comp = A.comp, z.comp = z.comp))
+
+}
Modified: pkg/ROptEst/R/getIneffDiff.R
===================================================================
--- pkg/ROptEst/R/getIneffDiff.R 2008-03-28 11:36:49 UTC (rev 81)
+++ pkg/ROptEst/R/getIneffDiff.R 2008-03-29 04:54:44 UTC (rev 82)
@@ -21,7 +21,6 @@
else
ineffUp <- (as.vector(res$A)*trafo - res$b^2*(radius^2-upRad^2))/upRisk
assign("ineff", ineffUp, envir = sys.frame(which = -4))
-
return(ineffUp - ineffLo)
}else{
if(is(L2Fam at distribution, "UnivariateDistribution")){
@@ -45,6 +44,7 @@
}
}
trafo <- L2Fam at param@trafo
+ p <- nrow(trafo)
neighbor at radius <- radius
res <- getInfRobIC(L2deriv = L2deriv, neighbor = neighbor, risk = risk,
Distr = L2Fam at distribution, DistrSymm = L2Fam at distrSymm,
@@ -52,11 +52,17 @@
Finfo = L2Fam at FisherInfo, trafo = trafo, z.start = z.start,
A.start = A.start, upper = upper.b, maxiter = MaxIter,
tol = eps, warn = warn)
- ineffLo <- (sum(diag(res$A%*%t(trafo))) - res$b^2*(radius^2-loRad^2))/loRisk
+ normtype(risk) <- res$normtype
+ std <- if(is(normtype(risk),"QFNorm")) QuadForm(normtype(risk)) else diag(p)
+
+
+ ineffLo <- (sum(diag(std%*%res$A%*%t(trafo))) -
+ res$b^2*(radius^2-loRad^2))/loRisk
if(upRad == Inf)
ineffUp <- res$b^2/upRisk
else
- ineffUp <- (sum(diag(res$A%*%t(trafo))) - res$b^2*(radius^2-upRad^2))/upRisk
+ ineffUp <- (sum(diag(std%*%res$A%*%t(trafo))) -
+ res$b^2*(radius^2-upRad^2))/upRisk
assign("ineff", ineffUp, envir = sys.frame(which = -4))
cat("current radius:\t", radius, "\tMSE-inefficiency difference:\t", ineffUp - ineffLo, "\n")
Modified: pkg/ROptEst/R/getInfRobIC_asBias.R
===================================================================
--- pkg/ROptEst/R/getInfRobIC_asBias.R 2008-03-28 11:36:49 UTC (rev 81)
+++ pkg/ROptEst/R/getInfRobIC_asBias.R 2008-03-29 04:54:44 UTC (rev 82)
@@ -5,19 +5,34 @@
risk = "asBias",
neighbor = "UncondNeighborhood"),
function(L2deriv, risk, neighbor, symm, trafo, maxiter,
- tol){
+ tol, ...){
minmaxBias(L2deriv, neighbor, biastype(risk), symm,
trafo, maxiter, tol)
})
setMethod("getInfRobIC", signature(L2deriv = "RealRandVariable",
risk = "asBias",
neighbor = "ContNeighborhood"),
- function(L2deriv, risk, neighbor, Distr, L2derivDistrSymm, z.start,
- A.start, trafo, maxiter, tol){
+ function(L2deriv, risk, neighbor, Distr, DistrSymm, L2derivSymm,
+ L2derivDistrSymm, z.start,
+ A.start, Finfo, trafo, maxiter, tol, ...){
+
+ normtype <- normtype(risk)
+
+ FI <- solve(trafo%*%solve(Finfo)%*%t(trafo))
+ if(is(normtype,"InfoNorm") || is(normtype,"SelfNorm") )
+ {QuadForm(normtype) <- PosSemDefSymmMatrix(FI);
+ normtype(risk) <- normtype}
+
+ comp <- .getComp(L2deriv, DistrSymm, L2derivSymm,
+ L2derivDistrSymm)
+
+ z.comp <- comp$"z.comp"
+ A.comp <- comp$"A.comp"
+
minmaxBias(L2deriv = L2deriv, neighbor = neighbor,
biastype = biastype(risk), normtype = normtype(risk),
- Distr = Distr, L2derivDistrSymm = L2derivDistrSymm,
- z.start = z.start, A.start = A.start, trafo = trafo,
+ Distr = Distr, z.start = z.start, A.start = A.start,
+ z.comp = z.comp, A.comp = A.comp, trafo = trafo,
maxiter = maxiter, tol = tol)
})
@@ -25,7 +40,7 @@
setMethod("minmaxBias", signature(L2deriv = "UnivariateDistribution",
neighbor = "ContNeighborhood",
biastype = "BiasType"),
- function(L2deriv, neighbor, biastype = symmetricBias(), symm,
+ function(L2deriv, neighbor, biastype, symm,
trafo, maxiter, tol){
zi <- sign(as.vector(trafo))
A <- as.matrix(zi)
@@ -49,8 +64,8 @@
cent(w) <- z
stand(w) <- A
clip(w) <- b
- weight(w) <- minbiasweight(w, neighbor = neighbor, biastype = biastype,
- normtype = NormType())
+ weight(w) <- minbiasweight(w, neighbor = neighbor, biastype = biastype,
+ normW= NormType())
return(list(A = A, a = zi*z, b = b, d = d, risk = Risk, info = info,
w = w, biastype = biastype, normtype = NormType()))
@@ -59,7 +74,7 @@
setMethod("minmaxBias", signature(L2deriv = "UnivariateDistribution",
neighbor = "TotalVarNeighborhood",
biastype = "BiasType"),
- function(L2deriv, neighbor, biastype = symmetricBias(),
+ function(L2deriv, neighbor, biastype,
symm, trafo,
maxiter, tol){
zi <- sign(as.vector(trafo))
@@ -74,7 +89,7 @@
if(zi == 1)
a <- -b*(1-p0)/(1-ws0)
else
- a <- b*(p0-ws0)/(1-ws0)
+ a <- -b*(p0-ws0)/(1-ws0)
info <- c("minimum asymptotic bias (lower case) solution")
Risk <- list(asCov = a^2*(p0-ws0) + (zi*a+b)^2*(1-p0), asBias = b)
@@ -82,30 +97,33 @@
w <- new("BdStWeight")
stand(w) <- A
clip(w) <- c(a, a+b)
- weight(w) <- minbiasweight(w, neighbor = neighbor, biastype = biastype,
- normtype = NormType())
+ weight(w) <- minbiasweight(w, neighbor = neighbor, biastype = biastype)
- return(list(A = A, a = a, b = b, d = 1, risk = Risk, info = info,
+ return(list(A = A, a = a, b = b, d = 0, risk = Risk, info = info,
w = w, biastype = biastype, normtype = NormType()))
})
setMethod("minmaxBias", signature(L2deriv = "RealRandVariable",
neighbor = "ContNeighborhood",
biastype = "BiasType"),
- function(L2deriv, neighbor, biastype, normtype, Distr, L2derivDistrSymm,
- z.start, A.start, trafo, maxiter, tol){
+ function(L2deriv, neighbor, biastype, normtype, Distr,
+ z.start, A.start, z.comp, A.comp, trafo, maxiter, tol){
+
eerg <- .LowerCaseMultivariate(L2deriv, neighbor, biastype,
- normtype, Distr, L2derivDistrSymm, trafo, z.start,
- A.start, maxiter, tol)
+ normtype, Distr, trafo, z.start,
+ A.start, z.comp = z.comp, A.comp = A.comp, maxiter, tol)
erg <- eerg$erg
b <- 1/erg$value
param <- erg$par
+ lA.comp <- sum(A.comp)
+
p <- nrow(trafo)
k <- ncol(trafo)
- A <- matrix(param[1:(p*k)], ncol=k, nrow=p)
+ A <- matrix(0, ncol=k, nrow=p)
+ A[A.comp] <- matrix(param[1:lA.comp], ncol=k, nrow=p)
z <- numeric(k)
- z[erg$z.comp] <- param[(p*k+1):length(param)]
+ z[z.comp] <- param[(lA.comp+1):length(param)]
a <- as.vector(A %*% z)
d <- numeric(p)
# computation of 'd', in case 'L2derivDistr' not abs. cont.
@@ -166,16 +184,15 @@
biastype = "onesidedBias"),
function(L2deriv, neighbor, biastype, symm,
trafo, maxiter, tol){
-
infotxt <- c("minimum asymptotic bias (lower case) solution")
noIC <- function(){
warntxt <- gettext("There exists no IC attaining the infimal maxBias.")
warning(warntxt)
- return(list(A = 1, a = 1, d = 0, b = 0,
+ return(list(A = matrix(1), a = 1, d = 0, b = 0,
Risk = list(asBias = 0, asCov = 0, warning = warntxt),
info = infotxt, w = NULL, biastype = biastype,
normtype = NormType()))}
- if(!is(L2deriv, "DiscreteDistribution"))
+ if(is(L2deriv, "DiscreteDistribution"))
{ if(is.finite(lowerCaseRadius(L2deriv, neighbor, risk = asMSE(), biastype)))
{
sign <- sign(biastype)
@@ -201,7 +218,8 @@
d0 <- 0
asCov <- v1^2*(1-p0)+v2^2*p0
Risk0 <- list(asBias = b0, asCov = asCov)
-
+ A0 <- matrix(A0,1,1)
+
w <- new("HampelWeight")
cent(w) <- z0
stand(w) <- A0
Modified: pkg/ROptEst/R/getInfRobIC_asGRisk.R
===================================================================
--- pkg/ROptEst/R/getInfRobIC_asGRisk.R 2008-03-28 11:36:49 UTC (rev 81)
+++ pkg/ROptEst/R/getInfRobIC_asGRisk.R 2008-03-29 04:54:44 UTC (rev 82)
@@ -33,7 +33,7 @@
## new
lower0 <- getL1normL2deriv(L2deriv = L2deriv, cent = z) /
(1 + neighbor at radius^2)
- upper0 <- sqrt( ( Finfo + z^2 )/(( 1 + neighbor at radius^2)^2 - 1) )
+ upper0 <- sqrt( as.numeric( Finfo + z^2 )/(( 1 + neighbor at radius^2)^2 - 1) )
if (!is.null(upper)|(iter == 1))
{lower <- .Machine$double.eps^0.75
}else{ lower <- lower0; upper <- upper0}
@@ -80,17 +80,23 @@
biastype = biastype, clip = b, cent = a, stand = A,
trafo = trafo)
Cov <- getInfV(L2deriv = L2deriv, neighbor = neighbor,
- biastype = biastype, clip = b/A, cent = z, stand = A)
+ biastype = biastype, clip = c0, cent = z, stand = A)
Risk <- c(Risk, list(asBias = b, asCov = Cov))
- w <- new("HampelWeight")
- cent(w) <- z
- stand(w) <- A
- clip(w) <- b
+ if(is(neighbor,"ContNeighborhood"))
+ {w <- new("HampelWeight")
+ clip(w) <- b
+ cent(w) <- z
+ stand(w) <- A
+ }else{
+ w <- new("BdStWeight")
+ clip(w) <- c(0,b)+a
+ stand(w) <- A
+ }
weight(w) <- getweight(w, neighbor = neighbor, biastype = biastype,
- normtype = normtype(risk))
+ normW = NormType())
return(list(A = A, a = a, b = b, d = NULL, risk = Risk, info = info, w = w,
biastype = biastype, normtype = normtype(risk)))
@@ -113,7 +119,8 @@
FI <- solve(trafo%*%solve(Finfo)%*%t(trafo))
if(is(normtype,"InfoNorm") || is(normtype,"SelfNorm") )
- {QuadForm(normtype) <- PosSemDefSymmMatrix(FI); normtype(risk) <- normtype}
+ {QuadForm(normtype) <- PosSemDefSymmMatrix(FI);
+ normtype(risk) <- normtype}
if(is.null(z.start)) z.start <- numeric(ncol(trafo))
if(is.null(A.start)) A.start <- trafo %*% solve(Finfo)
@@ -130,26 +137,13 @@
res$risk <- c(Risk, res$risk)
return(res)
}
- nrvalues <- length(L2deriv)
- z.comp <- rep(TRUE, nrvalues)
- A.comp <- matrix(TRUE, ncol = nrvalues, nrow = nrvalues)
- for(i in 1:nrvalues){
- if(is(L2derivDistrSymm[[i]], "SphericalSymmetry"))
- if(L2derivDistrSymm[[i]]@SymmCenter == 0)
- z.comp[i] <- FALSE
- }
- for(i in 1:(nrvalues-1))
- for(j in (i+1):nrvalues){
- if(is(DistrSymm, "SphericalSymmetry")){
- if((is(L2derivSymm[[i]], "OddSymmetric") & is(L2derivSymm[[j]], "EvenSymmetric"))
- | (is(L2derivSymm[[j]], "OddSymmetric") & is(L2derivSymm[[i]], "EvenSymmetric")))
- if((L2derivSymm[[i]]@SymmCenter == L2derivSymm[[j]]@SymmCenter)
- & (L2derivSymm[[i]]@SymmCenter == DistrSymm at SymmCenter))
- A.comp[i,j] <- FALSE
- }
- }
- A.comp[col(A.comp) < row(A.comp)] <- A.comp[col(A.comp) > row(A.comp)]
+ comp <- .getComp(L2deriv, DistrSymm, L2derivSymm,
+ L2derivDistrSymm)
+
+ z.comp <- comp$"z.comp"
+ A.comp <- comp$"A.comp"
+
w <- new("HampelWeight")
z <- z.start
A <- A.start
@@ -163,14 +157,7 @@
##
cent(w) <- z
stand(w) <- A
-
- if ((iter == 1)||is(normtype,"SelfNorm"))
- {normtype(risk) <- normtype <- updateNorm(normtype = normtype,
- FI = FI, L2 = L2deriv, neighbor = neighbor, biastype = biastype,
- Distr = Distr, V.comp = A.comp, cent = z, stand = A, w = w)}
-
- weight(w) <- getweight(w, neighbor = neighbor, biastype = biastype,
- normtype = normtype)
+
## new
lower0 <- getL1normL2deriv(L2deriv = L2deriv, cent = z, stand = A,
Distr = Distr, normtype = normtype)/(1+neighbor at radius^2)
@@ -181,7 +168,8 @@
{lower <- .Machine$double.eps^0.75;
if(is.null(upper)) upper <- 10*upper0
}else{ lower <- lower0; upper <- upper0}
- print(c(iter, lower,upper, lower0, upper0))
+
+
##
b <- try(uniroot(getInfClip,
## new
@@ -202,21 +190,20 @@
neighbor = neighbor, Distr = Distr, L2derivDistrSymm = L2derivDistrSymm,
z.start = z.start, A.start = A.start, trafo = trafo,
maxiter = maxiter, tol = tol)
- Risk <- getAsRisk(risk = risk, L2deriv = L2deriv, neighbor = neighbor,
+ normtype(risk) <- res$normtype
+ Risk <- getAsRisk(risk = risk, L2deriv = L2deriv, neighbor = neighbor,
biastype = biastype, clip = NULL,
cent = res$a, stand = res$A, trafo = trafo)
res$risk <- c(Risk, res$risk)
return(res)
}
clip(w) <- b
-
- if (is(normtype,"SelfNorm"))
- {normtype(risk) <- normtype <- updateNorm(normtype = normtype,
- FI = FI, L2 = L2deriv, neighbor = neighbor, biastype = biastype,
- Distr = Distr, V.comp = A.comp, cent = z, stand = A, w = w)}
+
weight(w) <- getweight(w, neighbor = neighbor, biastype = biastype,
- normtype = normtype)
+ normW = normtype)
+
+
z <- getInfCent(L2deriv = L2deriv, neighbor = neighbor,
biastype = biastype, Distr = Distr, z.comp = z.comp,
w = w)
@@ -224,6 +211,13 @@
biastype = biastype, Distr = Distr, A.comp = A.comp,
cent = z, trafo = trafo, w = w)
+ normtype.old <- normtype
+
+ if (is(normtype,"SelfNorm"))
+ {normtype(risk) <- normtype <- updateNorm(normtype = normtype,
+ L2 = L2deriv, neighbor = neighbor, biastype = biastype,
+ Distr = Distr, V.comp = A.comp, cent = z, stand = A, w = w)}
+
prec <- max(abs(b-b.old), max(abs(A-A.old)), max(abs(z-z.old)))
cat("current precision in IC algo:\t", prec, "\n")
if(prec < tol) break
@@ -235,14 +229,10 @@
if (onesetLM){
cent(w) <- z
stand(w) <- A
- if (is(normtype,"SelfNorm"))
- {normtype(risk) <- normtype <- updateNorm(normtype = normtype,
- FI = FI, L2 = L2deriv, neighbor = neighbor, biastype = biastype,
- Distr = Distr, V.comp = A.comp, cent = z, stand = A, w = w)}
-
weight(w) <- getweight(w, neighbor = neighbor, biastype = biastype,
- normtype = normtype)
- }
+ normW = normtype)
+ } else normtype <- normtype.old
+
a <- as.vector(A %*% z)
info <- paste("optimally robust IC for", sQuote(class(risk)[1]))
Risk <- getAsRisk(risk = risk, L2deriv = L2deriv, neighbor = neighbor,
Modified: pkg/ROptEst/R/getInfRobIC_asHampel.R
===================================================================
--- pkg/ROptEst/R/getInfRobIC_asHampel.R 2008-03-28 11:36:49 UTC (rev 81)
+++ pkg/ROptEst/R/getInfRobIC_asHampel.R 2008-03-29 04:54:44 UTC (rev 82)
@@ -70,13 +70,18 @@
# biastype = biastype, clip = b, cent = a, stand = A)$asCov
Risk <- list(asCov = Cov, asBias = b, asMSE = Cov + neighbor at radius^2*b^2)
- w <- new("HampelWeight")
- clip(w) <- b
- cent(w) <- z
- stand(w) <- A
-
+ if(is(neighbor,"ContNeighborhood"))
+ {w <- new("HampelWeight")
+ clip(w) <- b
+ cent(w) <- z
+ stand(w) <- A
+ }else{
+ w <- new("BdStWeight")
+ clip(w) <- c(0,b)+a
+ stand(w) <- A
+ }
weight(w) <- getweight(w, neighbor = neighbor, biastype = biastype,
- normtype = NormType())
+ normW = NormType())
return(list(A = A, a = a, b = b, d = NULL, risk = Risk, info = info,
w = w, biastype = biastype, normtype = NormType()))
@@ -91,11 +96,14 @@
biastype <- biastype(risk)
normtype <- normtype(risk)
+ p <- nrow(trafo)
FI <- solve(trafo%*%solve(Finfo)%*%t(trafo))
if(is(normtype,"InfoNorm") || is(normtype,"SelfNorm") )
{QuadForm(normtype) <- PosSemDefSymmMatrix(FI); normtype(risk) <- normtype}
+ std <- if(is(normtype,"QFNorm")) QuadForm(normtype) else diag(p)
+
if(is.null(z.start)) z.start <- numeric(ncol(trafo))
if(is.null(A.start)) A.start <- trafo
@@ -127,32 +135,15 @@
"=> the minimum asymptotic bias (lower case) solution is returned\n")
return(res)
}
- nrvalues <- length(L2deriv)
- z.comp <- rep(TRUE, nrvalues)
- A.comp <- matrix(1, ncol = nrvalues, nrow = nrvalues)
- for(i in 1:nrvalues){
- if(is(L2derivDistrSymm[[i]], "SphericalSymmetry"))
- if(L2derivDistrSymm[[i]]@SymmCenter == 0){
- z.comp[i] <- FALSE
- A.comp[i,i] <- 2
- }
- }
- for(i in 1:(nrvalues-1))
- for(j in (i+1):nrvalues){
- if(z.comp[i] | z.comp[j]){
- if(is(DistrSymm, "SphericalSymmetry")){
- if((is(L2derivSymm[[i]], "OddSymmetric") & is(L2derivSymm[[j]], "EvenSymmetric"))
- | (is(L2derivSymm[[j]], "OddSymmetric") & is(L2derivSymm[[i]], "EvenSymmetric")))
- if((L2derivSymm[[i]]@SymmCenter == L2derivSymm[[j]]@SymmCenter)
- & (L2derivSymm[[i]]@SymmCenter == DistrSymm at SymmCenter))
- A.comp[i,j] <- 0
- }else{
- A.comp[i,j] <- 2
- }
- }
- }
- A.comp[col(A.comp) < row(A.comp)] <- A.comp[col(A.comp) > row(A.comp)]
+ comp <- .getComp(L2deriv, DistrSymm, L2derivSymm,
+ L2derivDistrSymm)
+
+ z.comp <- comp$"z.comp"
+ A.comp <- comp$"A.comp"
+
+ w <- new("HampelWeight")
+ clip(w) <- b
z <- z.start
A <- A.start
iter <- 0
@@ -162,19 +153,26 @@
A.old <- A
cent(w) <- z
stand(w) <- A
- normtype(risk) <- normtype <- updateNorm(normtype = normtype, FI = FI,
- L2 = L2deriv, neighbor = neighbor, biastype = biastype,
- Distr = Distr, V.comp = A.comp, cent = z, stand = A, w = w)
weight(w) <- getweight(w, neighbor = neighbor, biastype = biastype,
- normtype = normtype)
+ normW = normtype)
+
z <- getInfCent(L2deriv = L2deriv, neighbor = neighbor,
biastype = biastype, Distr = Distr, z.comp = z.comp,
w = w)
A <- getInfStand(L2deriv = L2deriv, neighbor = neighbor,
biastype = biastype, Distr = Distr, A.comp = A.comp,
cent = z, trafo = trafo, w = w)
+
+ normtype.old <- normtype
+ if(is(normtype,"SelfNorm")){
+ normtype(risk) <- normtype <- updateNorm(normtype = normtype,
+ L2 = L2deriv, neighbor = neighbor, biastype = biastype,
+ Distr = Distr, V.comp = A.comp, cent = z, stand = A, w = w)
+ std <- QuadForm(normtype)
+ }
+
prec <- max(max(abs(A-A.old)), max(abs(z-z.old)))
cat("current precision in IC algo:\t", prec, "\n")
if(prec < tol) break
@@ -186,13 +184,12 @@
if (onesetLM){
cent(w) <- z
stand(w) <- A
- normtype(risk) <- normtype <- updateNorm(normtype = normtype, FI = FI,
- L2 = L2deriv, neighbor = neighbor, biastype = biastype,
- Distr = Distr, V.comp = A.comp, cent = z, stand = A, w = w)
weight(w) <- getweight(w, neighbor = neighbor, biastype = biastype,
- normtype = normtype)
+ normW = normtype)
}
+ else normtype <- normtype.old
+
info <- paste("optimally robust IC for 'asHampel' with bound =", round(b,3))
a <- as.vector(A %*% z)
Cov <- getInfV(L2deriv = L2deriv, neighbor = neighbor,
@@ -202,8 +199,8 @@
#getAsRisk(risk = asCov(), L2deriv = L2deriv, neighbor = neighbor,
# biastype = biastype, Distr = Distr, clip = b, cent = a,
# stand = A)$asCov
- Risk <- list(trAsCov = sum(diag(Cov)), asCov = Cov, asBias = b,
- asMSE = sum(diag(Cov)) + neighbor at radius^2*b^2)
+ Risk <- list(trAsCov = sum(diag(std%*%Cov)), asCov = Cov, asBias = b,
+ asMSE = sum(diag(std%*%Cov)) + neighbor at radius^2*b^2)
return(list(A = A, a = a, b = b, d = NULL, risk = Risk, info = info,
w = w, biastype = biastype, normtype = normtype))
Modified: pkg/ROptEst/R/getInfRobIC_asUnOvShoot.R
===================================================================
--- pkg/ROptEst/R/getInfRobIC_asUnOvShoot.R 2008-03-28 11:36:49 UTC (rev 81)
+++ pkg/ROptEst/R/getInfRobIC_asUnOvShoot.R 2008-03-29 04:54:44 UTC (rev 82)
@@ -92,7 +92,7 @@
if(!is.numeric(c0)){
if(warn) cat("The IC algorithm did not converge!\n",
"=> the minimum asymptotic bias (lower case) solution is returned\n")
- res <- getInfRobIC(L2deriv = L2deriv, risk = asBias(),
+ res <- getInfRobIC(L2deriv = L2deriv, risk = asBias(biastype = biastype),
neighbor = neighbor, Finfo = Finfo,
symm = symm, trafo = trafo, upper = upper,
maxiter = maxiter, tol = tol, warn = warn)
Modified: pkg/ROptEst/R/getInfV.R
===================================================================
--- pkg/ROptEst/R/getInfV.R 2008-03-28 11:36:49 UTC (rev 81)
+++ pkg/ROptEst/R/getInfV.R 2008-03-29 04:54:44 UTC (rev 82)
@@ -92,9 +92,10 @@
nu2 <- nu(biastype)[2]
c1 <- cent - clip/nu1
c2 <- cent + clip/nu2
- return(stand^2*(m2df(L2deriv, c2) - m2df(L2deriv, c1)
- + 2 * cent *(m1df(L2deriv, c1) - m1df(L2deriv, c2))
- + cent^2 * (p(L2deriv)(c2) -p(L2deriv)(c1))
- + clip^2 * (1-p(L2deriv)(c2)/nu2^2 +p(L2deriv)(c1)/nu1^2)
- ))
+ V0 <- m2df(L2deriv, c2) - m2df(L2deriv, c1)
+ V1 <- m1df(L2deriv, c2) - m1df(L2deriv, c1)
+ V2 <- p(L2deriv)(c2) -p(L2deriv)(c1)
+ V3 <- (1-p(L2deriv)(c2))/nu2^2 +p(L2deriv)(c1)/nu1^2
+ V <- stand^2*( V0 - 2 * cent * V1 + cent^2 * V2 + clip^2 * V3)
+ return(V)
})
Modified: pkg/ROptEst/R/getRiskIC.R
===================================================================
--- pkg/ROptEst/R/getRiskIC.R 2008-03-28 11:36:49 UTC (rev 81)
+++ pkg/ROptEst/R/getRiskIC.R 2008-03-29 04:54:44 UTC (rev 82)
@@ -7,7 +7,7 @@
L2Fam = "missing"),
function(IC, risk){
L2Fam <- eval(IC at CallL2Fam)
- getRiskIC(IC, risk, L2Fam)
+ getRiskIC(IC = IC, risk = risk, L2Fam = L2Fam)
})
setMethod("getRiskIC", signature(IC = "HampIC",
@@ -16,7 +16,7 @@
L2Fam = "missing"),
function(IC, risk, L2Fam){
Cov <- IC at Risks[["asCov"]]
- return(list(asCov = list(distribution = .getDistr(IC at L2Fam), value = Cov)))
+ return(list(asCov = list(distribution = .getDistr(L2Fam), value = Cov)))
})
@@ -25,10 +25,11 @@
###############################################################################
setMethod("getBiasIC", signature(IC = "HampIC",
neighbor = "UncondNeighborhood"),
- function(IC, neighbor, L2Fam){
+ function(IC, neighbor, L2Fam,...){
if(missing(L2Fam))
- {misF <- TRUE; L2Fam <- eval(IC at CallL2Fam)}
-
+ {misF <- TRUE;
+ L2Fam <- eval(IC at CallL2Fam)}
+ print(L2Fam)
return(list(asBias = list(distribution = .getDistr(L2Fam),
neighborhood = neighbor at type, value = IC at Risks[["asBias"]])))
})
Modified: pkg/ROptEst/R/leastFavorableRadius.R
===================================================================
--- pkg/ROptEst/R/leastFavorableRadius.R 2008-03-28 11:36:49 UTC (rev 81)
+++ pkg/ROptEst/R/leastFavorableRadius.R 2008-03-29 04:54:44 UTC (rev 82)
@@ -14,6 +14,14 @@
stop("'rho' not in (0,1)")
biastype <- biastype(risk)
+ normtype <- normtype(risk)
+
+ FI0 <- L2Fam at param@trafo%*%solve(L2Fam at FisherInfo)%*%t(L2Fam at param@trafo)
+ FI <- solve(FI0)
+ if(is(normtype,"InfoNorm") || is(normtype,"SelfNorm") )
+ {QuadForm(normtype) <- PosSemDefSymmMatrix(FI);
+ normtype(risk) <- normtype}
+
L2derivDim <- numberOfMaps(L2Fam at L2deriv)
if(L2derivDim == 1){
leastFavFct <- function(r, L2Fam, neighbor, risk, rho,
@@ -40,7 +48,8 @@
}
if(upRad == Inf){
- bmin <- getAsRisk(risk = asBias(), L2deriv = L2Fam at L2derivDistr[[1]],
+ bmin <- getAsRisk(risk = asBias(biastype = biastype),
+ L2deriv = L2Fam at L2derivDistr[[1]],
neighbor = neighbor, biastype = biastype,
trafo = L2Fam at param@trafo, symm = L2Fam at L2derivSymm[[1]])
upRisk <- bmin^2
@@ -95,6 +104,10 @@
L2derivDistrSymm <- new("DistrSymmList", L2)
}
}
+
+ std <- if(is(normtype,"QFNorm"))
+ QuadForm(normtype) else diag(nrow(L2Fam at param@trafo))
+
leastFavFct <- function(r, L2Fam, neighbor, risk, rho,
z.start, A.start, upper.b, MaxIter, eps, warn){
loRad <- r*rho
@@ -106,7 +119,7 @@
trafo <- L2Fam at param@trafo
if(identical(all.equal(loRad, 0), TRUE)){
loRad <- 0
- loRisk <- sum(diag(solve(L2Fam at FisherInfo)))
+ loRisk <- sum(diag(std%*%FI0))
}else{
neighbor at radius <- loRad
resLo <- getInfRobIC(L2deriv = L2deriv, neighbor = neighbor, risk = risk,
@@ -115,17 +128,24 @@
Finfo = L2Fam at FisherInfo, trafo = trafo, z.start = z.start,
A.start = A.start, upper = upper.b, maxiter = MaxIter,
tol = eps, warn = warn)
- loRisk <- getAsRisk(risk = risk, L2deriv = L2deriv, neighbor = neighbor,
+ riskLo <- risk
+ normtype(riskLo) <- resLo$normtype
+ loRisk <- getAsRisk(risk = riskLo, L2deriv = L2deriv, neighbor = neighbor,
biastype = biastype, clip = resLo$b, cent = resLo$a,
stand = resLo$A, trafo = trafo)[[1]]
}
if(upRad == Inf){
- bmin <- getAsRisk(risk = asBias(), L2deriv = L2deriv, neighbor = neighbor,
- biastype = biastype,
- Distr = L2Fam at distribution, L2derivDistrSymm = L2Fam at L2derivDistrSymm,
- trafo = trafo, z.start = z.start, A.start = A.start,
- maxiter = maxiter, tol = tol)$asBias
+ bmin <- getAsRisk(risk = asBias(biastype = biastype(risk),
+ normtype = normtype), L2deriv = L2deriv,
+ neighbor = neighbor, biastype = biastype,
+ Distr = L2Fam at distribution,
+ DistrSymm = L2Fam at distrSymm,
+ L2derivSymm = L2derivSymm,
+ L2derivDistrSymm= L2derivDistrSymm,
+ trafo = trafo, z.start = z.start,
+ A.start = A.start,
+ maxiter = maxiter, tol = tol)$asBias
upRisk <- bmin^2
}else{
neighbor at radius <- upRad
@@ -135,7 +155,9 @@
Finfo = L2Fam at FisherInfo, trafo = trafo, z.start = z.start,
A.start = A.start, upper = upper.b, maxiter = maxiter,
tol = tol, warn = warn)
- upRisk <- getAsRisk(risk = risk, L2deriv = L2deriv, neighbor = neighbor,
+ riskUp <- risk
+ normtype(riskUp) <- resUp$normtype
+ upRisk <- getAsRisk(risk = riskUp, L2deriv = L2deriv, neighbor = neighbor,
biastype = biastype,
clip = resUp$b, cent = resUp$a, stand = resUp$A, trafo = trafo)[[1]]
}
Modified: pkg/ROptEst/R/lowerCaseRadius.R
===================================================================
--- pkg/ROptEst/R/lowerCaseRadius.R 2008-03-28 11:36:49 UTC (rev 81)
+++ pkg/ROptEst/R/lowerCaseRadius.R 2008-03-29 04:54:44 UTC (rev 82)
@@ -5,7 +5,7 @@
neighbor = "ContNeighborhood",
risk = "asMSE",
biastype = "ANY"),
- function(L2Fam, neighbor, risk, biastype = symmetricBias()){
+ function(L2Fam, neighbor, risk, biastype){
if(length(L2Fam at param) != 1) stop("not yet implemented")
D1 <- L2Fam at distribution
@@ -55,7 +55,7 @@
neighbor = "TotalVarNeighborhood",
risk = "asMSE",
biastype = "ANY"),
- function(L2Fam, neighbor, risk, biastype = symmetricBias()){
+ function(L2Fam, neighbor, risk, biastype){
if(length(L2Fam at param) != 1) stop("not yet implemented")
D1 <- L2Fam at distribution
Modified: pkg/ROptEst/R/radiusMinimaxIC.R
===================================================================
--- pkg/ROptEst/R/radiusMinimaxIC.R 2008-03-28 11:36:49 UTC (rev 81)
+++ pkg/ROptEst/R/radiusMinimaxIC.R 2008-03-29 04:54:44 UTC (rev 82)
@@ -15,6 +15,10 @@
stop("'upRad < loRad' is not fulfilled")
biastype <- biastype(risk)
L2derivDim <- numberOfMaps(L2Fam at L2deriv)
+
+ if(is(normtype(risk),"SelfNorm")||is(normtype(risk),"InfoNorm"))
+ upRad <- min(upRad,10)
+
if(L2derivDim == 1){
ow <- options("warn")
options(warn = -1)
@@ -38,8 +42,10 @@
}
if(upRad == Inf){
- bmin <- getAsRisk(risk = asBias(), L2deriv = L2Fam at L2derivDistr[[1]],
- neighbor = neighbor, biastype = biastype, trafo = L2Fam at param@trafo)$asBias
+ bmin <- getAsRisk(risk = asBias(biastype = biastype),
+ L2deriv = L2Fam at L2derivDistr[[1]],
+ neighbor = neighbor, biastype = biastype,
+ trafo = L2Fam at param@trafo)$asBias
upRisk <- bmin^2
}else{
neighbor at radius <- upRad
@@ -93,7 +99,21 @@
L2derivDistrSymm <- new("DistrSymmList", L2)
}
}
+ normtype <- normtype(risk)
+
+ Finfo <- L2Fam at FisherInfo
trafo <- L2Fam at param@trafo
+
+ p <- nrow(trafo)
+ FI0 <- trafo%*%solve(Finfo)%*%t(trafo)
+ FI <- solve(FI0)
+
+ if(is(normtype,"InfoNorm") || is(normtype,"SelfNorm") )
+ {QuadForm(normtype) <- PosSemDefSymmMatrix(FI);
+ normtype(risk) <- normtype}
+ std <- if(is(normtype,"QFNorm")) QuadForm(normtype) else diag(p)
+
+
ow <- options("warn")
options(warn = -1)
upper.b <- upper
@@ -102,7 +122,7 @@
if(identical(all.equal(loRad, 0), TRUE)){
loRad <- 0
- loRisk <- sum(diag(solve(L2Fam at FisherInfo)))
+ loRisk <- sum(diag(std%*%FI0))
}else{
neighbor at radius <- loRad
resLo <- getInfRobIC(L2deriv = L2deriv, neighbor = neighbor, risk = risk,
@@ -111,14 +131,24 @@
Finfo = L2Fam at FisherInfo, trafo = trafo, z.start = z.start,
A.start = A.start, upper = upper.b, maxiter = maxiter,
tol = tol, warn = warn)
- loRisk <- getAsRisk(risk = risk, L2deriv = L2deriv, neighbor = neighbor,
- biastype = biastype, clip = resLo$b, cent = resLo$a, stand = resLo$A, trafo = trafo)[[1]]
+ riskLo <- risk
+ normtype(riskLo) <- resLo$normtype
+ loRisk <- getAsRisk(risk = riskLo, L2deriv = L2deriv,
+ neighbor = neighbor, biastype = biastype,
+ clip = resLo$b, cent = resLo$a,
+ stand = resLo$A, trafo = trafo)[[1]]
}
if(upRad == Inf){
- bmin <- getAsRisk(risk = asBias(), L2deriv = L2deriv, neighbor = neighbor,
- biastype = biastype, Distr = L2Fam at distribution, L2derivDistrSymm = L2Fam at L2derivDistrSymm,
- trafo = trafo, z.start = z.start, A.start = A.start,
+ bmin <- getAsRisk(risk = asBias(biastype = biastype(risk),
+ normtype = normtype), L2deriv = L2deriv,
+ neighbor = neighbor, biastype = biastype,
+ Distr = L2Fam at distribution,
+ DistrSymm = L2Fam at distrSymm,
+ L2derivSymm = L2derivSymm,
+ L2derivDistrSymm= L2derivDistrSymm,
+ trafo = trafo, z.start = z.start,
+ A.start = A.start,
maxiter = maxiter, tol = tol)$asBias
upRisk <- bmin^2
}else{
@@ -129,10 +159,11 @@
Finfo = L2Fam at FisherInfo, trafo = trafo, z.start = z.start,
A.start = A.start, upper = upper.b, maxiter = maxiter,
tol = tol, warn = warn)
- upRisk <- getAsRisk(risk = risk, L2deriv = L2deriv, neighbor = neighbor,
+ riskUp <- risk
+ normtype(riskUp) <- resUp$normtype
+ upRisk <- getAsRisk(risk = riskUp, L2deriv = L2deriv, neighbor = neighbor,
biastype = biastype, clip = resUp$b, cent = resUp$a, stand = resUp$A, trafo = trafo)[[1]]
}
-
leastFavR <- uniroot(getIneffDiff, lower = lower, upper = upper,
tol = .Machine$double.eps^0.25, L2Fam = L2Fam, neighbor = neighbor,
z.start = z.start, A.start = A.start, upper.b = upper.b, risk = risk,
Modified: pkg/ROptEst/R/updateNorm.R
===================================================================
--- pkg/ROptEst/R/updateNorm.R 2008-03-28 11:36:49 UTC (rev 81)
+++ pkg/ROptEst/R/updateNorm.R 2008-03-29 04:54:44 UTC (rev 82)
@@ -1,11 +1,12 @@
-setMethod("updateNorm", "NormType", function(normtype, ...) normtype)
-setMethod("updateNorm", "InfoNorm", function(normtype, FI, ...)
- {QuadForm(normtype) <- PosSemDefSymmMatrix(FI); normtype})
+#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, ...)
+ 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(solve(Cv)); normtype})
+ QuadForm(normtype) <- PosSemDefSymmMatrix(solve(Cv))
+ normtype})
\ No newline at end of file
Modified: pkg/ROptEst/chm/00Index.html
===================================================================
--- pkg/ROptEst/chm/00Index.html 2008-03-28 11:36:49 UTC (rev 81)
+++ pkg/ROptEst/chm/00Index.html 2008-03-29 04:54:44 UTC (rev 82)
@@ -293,10 +293,6 @@
<table width="100%">
<tr><td width="25%"><a href="updateNorm-methods.html">updateNorm</a></td>
<td>Methods for Function updateNorm in Package ‘ROptEst’ </td></tr>
-<tr><td width="25%"><a href="updateNorm-methods.html">updateNorm,InfoNorm-method</a></td>
-<td>Methods for Function updateNorm in Package ‘ROptEst’ </td></tr>
-<tr><td width="25%"><a href="updateNorm-methods.html">updateNorm,NormType-method</a></td>
-<td>Methods for Function updateNorm in Package ‘ROptEst’ </td></tr>
<tr><td width="25%"><a href="updateNorm-methods.html">updateNorm,SelfNorm-method</a></td>
<td>Methods for Function updateNorm in Package ‘ROptEst’ </td></tr>
<tr><td width="25%"><a href="updateNorm-methods.html">updateNorm-methods</a></td>
Modified: pkg/ROptEst/chm/ROptEst.chm
===================================================================
(Binary files differ)
Modified: pkg/ROptEst/chm/ROptEst.toc
===================================================================
--- pkg/ROptEst/chm/ROptEst.toc 2008-03-28 11:36:49 UTC (rev 81)
+++ pkg/ROptEst/chm/ROptEst.toc 2008-03-29 04:54:44 UTC (rev 82)
@@ -494,14 +494,6 @@
<param name="Local" value="updateNorm-methods.html">
</OBJECT>
<LI> <OBJECT type="text/sitemap">
-<param name="Name" value="updateNorm,InfoNorm-method">
-<param name="Local" value="updateNorm-methods.html">
-</OBJECT>
-<LI> <OBJECT type="text/sitemap">
-<param name="Name" value="updateNorm,NormType-method">
-<param name="Local" value="updateNorm-methods.html">
-</OBJECT>
-<LI> <OBJECT type="text/sitemap">
<param name="Name" value="updateNorm,SelfNorm-method">
<param name="Local" value="updateNorm-methods.html">
</OBJECT>
Modified: pkg/ROptEst/chm/getAsRisk.html
===================================================================
--- pkg/ROptEst/chm/getAsRisk.html 2008-03-28 11:36:49 UTC (rev 81)
+++ pkg/ROptEst/chm/getAsRisk.html 2008-03-29 04:54:44 UTC (rev 82)
@@ -68,8 +68,8 @@
## S4 method for signature 'asBias, RealRandVariable,
## ContNeighborhood, ANY':
-getAsRisk(risk, L2deriv, neighbor, biastype, Distr, L2derivDistrSymm, trafo,
- z.start, A.start, maxiter, tol, norm = EuclideanNorm)
+getAsRisk(risk, L2deriv, neighbor, biastype, Distr, DistrSymm, L2derivSymm,
+ L2derivDistrSymm, trafo, z.start, A.start, maxiter, tol, norm = EuclideanNorm)
## S4 method for signature 'asCov, UnivariateDistribution,
## ContNeighborhood, ANY':
@@ -136,6 +136,12 @@
<tr valign="top"><td><code>Distr</code></td>
<td>
object of class <code>"Distribution"</code>. </td></tr>
+<tr valign="top"><td><code>DistrSymm</code></td>
+<td>
+object of class <code>"DistributionSymmetry"</code>. </td></tr>
+<tr valign="top"><td><code>L2derivSymm</code></td>
+<td>
+object of class <code>"FunSymmList"</code>. </td></tr>
<tr valign="top"><td><code>L2derivDistrSymm</code></td>
<td>
object of class <code>"DistrSymmList"</code>. </td></tr>
Modified: pkg/ROptEst/chm/getInfRobIC.html
===================================================================
--- pkg/ROptEst/chm/getInfRobIC.html 2008-03-28 11:36:49 UTC (rev 81)
+++ pkg/ROptEst/chm/getInfRobIC.html 2008-03-29 04:54:44 UTC (rev 82)
@@ -57,8 +57,8 @@
## S4 method for signature 'RealRandVariable, asBias,
## ContNeighborhood':
-getInfRobIC(L2deriv, risk, neighbor, Distr,
- L2derivDistrSymm, z.start, A.start, trafo, maxiter, tol)
+getInfRobIC(L2deriv, risk, neighbor, Distr, DistrSymm, L2derivSymm,
+ L2derivDistrSymm, z.start, A.start, Finfo, trafo, maxiter, tol)
## S4 method for signature 'UnivariateDistribution,
## asHampel, UncondNeighborhood':
Modified: pkg/ROptEst/chm/minmaxBias.html
===================================================================
--- pkg/ROptEst/chm/minmaxBias.html 2008-03-28 11:36:49 UTC (rev 81)
+++ pkg/ROptEst/chm/minmaxBias.html 2008-03-29 04:54:44 UTC (rev 82)
@@ -56,7 +56,7 @@
## S4 method for signature 'RealRandVariable,
## ContNeighborhood, BiasType':
minmaxBias(L2deriv, neighbor, biastype, Distr,
- L2derivDistrSymm, z.start, A.start, trafo, maxiter, tol)
+ L2derivDistrSymm, z.start, A.start, z.comp, A.comp, trafo, maxiter, tol)
</pre>
@@ -92,6 +92,12 @@
<tr valign="top"><td><code>A.start</code></td>
<td>
initial value for the standardizing matrix. </td></tr>
+<tr valign="top"><td><code>z.comp</code></td>
+<td>
+<code>logical</code> indicator which indices need to be computed and which are 0 due to symmetry. </td></tr>
+<tr valign="top"><td><code>A.comp</code></td>
+<td>
+<code>matrix</code> of <code>logical</code> indicator which indices need to be computed and which are 0 due to symmetry.</td></tr>
<tr valign="top"><td><code>trafo</code></td>
<td>
matrix: transformation of the parameter. </td></tr>
Modified: pkg/ROptEst/chm/updateNorm-methods.html
===================================================================
--- pkg/ROptEst/chm/updateNorm-methods.html 2008-03-28 11:36:49 UTC (rev 81)
+++ pkg/ROptEst/chm/updateNorm-methods.html 2008-03-29 04:54:44 UTC (rev 82)
@@ -7,8 +7,6 @@
<table width="100%"><tr><td>updateNorm-methods(ROptEst)</td><td align="right">R Documentation</td></tr></table><object type="application/x-oleobject" classid="clsid:1e2a7bd0-dab9-11d0-b93a-00c04fc99f9e">
<param name="keyword" value="R: updateNorm-methods">
<param name="keyword" value="R: updateNorm">
-<param name="keyword" value="R: updateNorm,NormType-method">
-<param name="keyword" value="R: updateNorm,InfoNorm-method">
<param name="keyword" value="R: updateNorm,SelfNorm-method">
<param name="keyword" value=" Methods for Function updateNorm in Package ‘ROptEst’">
</object>
@@ -27,13 +25,9 @@
<h3>Usage</h3>
<pre>updateNorm(normtype, ...)
-## S4 method for signature 'NormType':
-updateNorm(normtype, ...)
-## S4 method for signature 'InfoNorm':
-updateNorm(normtype, FI, ...)
## S4 method for signature 'SelfNorm':
updateNorm(normtype, L2, neighbor, biastype, Distr, V.comp,
- cent, stand, w, ...)
+ cent, stand, w)
</pre>
@@ -46,9 +40,6 @@
<tr valign="top"><td><code>...</code></td>
<td>
further arguments to be passed to specific methods.</td></tr>
-<tr valign="top"><td><code>FI</code></td>
-<td>
-matrix: Fisher Information</td></tr>
<tr valign="top"><td><code>L2</code></td>
<td>
L2derivative</td></tr>
@@ -80,8 +71,8 @@
<p>
<code>updateNorm</code> is used internally in the opt-IC-algorithm to be
-able to work with a norm that depends on the Fisher information at a certain
-parameter (<code>InfoType</code>) or on the current covariance (<code>SelfNorm</code>)
+able to work with a norm that depends on the current covariance
+(<code>SelfNorm</code>)
</p>
@@ -98,10 +89,6 @@
<h3>Methods</h3>
<dl>
-<dt>updateNorm</dt><dd><code>signature(normtype = "NormType")</code>: leaves the norm unchanged;</dd>
-<dt>updateNorm</dt><dd><code>signature(normtype = "InfoNorm")</code>:
-udates the norm in the information-standardized case; just used
-internally in the opt-IC-Algorithm. </dd>
<dt>updateNorm</dt><dd><code>signature(normtype = "SelfNorm")</code>:
udates the norm in the self-standardized case; just used
internally in the opt-IC-Algorithm. </dd>
Modified: pkg/ROptEst/inst/scripts/NormalLocationScaleModel.R
===================================================================
--- pkg/ROptEst/inst/scripts/NormalLocationScaleModel.R 2008-03-28 11:36:49 UTC (rev 81)
+++ pkg/ROptEst/inst/scripts/NormalLocationScaleModel.R 2008-03-29 04:54:44 UTC (rev 82)
@@ -27,6 +27,18 @@
plot(N0.IC1)
infoPlot(N0.IC1)
+system.time(N0.IC1.i <- optIC(model = N0.Rob1, risk = asMSE(normtype=InfoNorm())), gcFirst = TRUE)
+checkIC(N0.IC1.i)
+Risks(N0.IC1.i)
+plot(N0.IC1.i)
+infoPlot(N0.IC1.i)
+
+system.time(N0.IC1.s <- optIC(model = N0.Rob1, risk = asMSE(normtype=SelfNorm())), gcFirst = TRUE)
+checkIC(N0.IC1.s)
+Risks(N0.IC1.s)
+plot(N0.IC1.s)
+infoPlot(N0.IC1.s)
+
# lower case solutions
(N0.IC2 <- optIC(model = N0.Rob1, risk = asBias(), tol = 1e-10))
checkIC(N0.IC2)
@@ -34,6 +46,12 @@
plot(N0.IC2)
infoPlot(N0.IC2)
+(N0.IC2.i <- optIC(model = N0.Rob1, risk = asBias(normtype=InfoNorm()), tol = 1e-10))
+checkIC(N0.IC2.i)
+Risks(N0.IC2.i)
+plot(N0.IC2.i)
+infoPlot(N0.IC2.i)
+
# Hampel solution
(N0.IC3 <- optIC(model = N0.Rob1, risk = asHampel(bound = clip(N0.IC1))))
checkIC(N0.IC3)
@@ -50,6 +68,13 @@
plot(N0.IC4)
infoPlot(N0.IC4)
+(N0.IC4.i <- radiusMinimaxIC(L2Fam=N0, neighbor=ContNeighborhood(),
+ risk=asMSE(normtype=InfoNorm()), loRad=0, upRad=Inf))
+checkIC(N0.IC4.i)
+Risks(N0.IC4.i)
+plot(N0.IC4.i)
+infoPlot(N0.IC4.i)
+
# least favorable radius
# (may take quite some time!)
#N0.r.rho1 <- leastFavorableRadius(L2Fam=N0, neighbor=ContNeighborhood(),
Modified: pkg/ROptEst/inst/scripts/PoissonModel.R
===================================================================
--- pkg/ROptEst/inst/scripts/PoissonModel.R 2008-03-28 11:36:49 UTC (rev 81)
+++ pkg/ROptEst/inst/scripts/PoissonModel.R 2008-03-29 04:54:44 UTC (rev 82)
@@ -34,6 +34,16 @@
Risks(IC1)
plot(IC1)
+(IC1.p <- optIC(model=RobP1, risk=asMSE(biastype=positiveBias())))
+checkIC(IC1.p)
+Risks(IC1.p)
+plot(IC1.p)
+
+(IC1.a <- optIC(model=RobP1, risk=asMSE(biastype=asymmetricBias(nu = c(1,0.2)))))
+checkIC(IC1.a)
+Risks(IC1.a)
+plot(IC1.a)
+
(IC2 <- optIC(model=RobP2, risk=asMSE()))
checkIC(IC2)
Risks(IC2)
@@ -46,6 +56,17 @@
Risks(IC3)
plot(IC3)
+(IC3.p <- optIC(model=RobP1, risk=asBias(biastype=positiveBias())))
+checkIC(IC3.p)
+Risks(IC3.p)
+plot(IC3.p)
+
+(IC3.a <- optIC(model=RobP1, risk=asBias(biastype=asymmetricBias(nu = c(1,0.2)))))
+checkIC(IC3.a)
+Risks(IC3.a)
+plot(IC3.a)
+
+
(IC4 <- optIC(model=RobP2, risk=asBias()))
checkIC(IC4)
Risks(IC4)
@@ -70,6 +91,18 @@
Risks(IC7)
plot(IC7)
+(IC7.p <- radiusMinimaxIC(L2Fam=P, neighbor=ContNeighborhood(),
+ risk=asMSE(biastype=positiveBias()), loRad=0, upRad=0.5))
+checkIC(IC7.p)
+Risks(IC7.p)
+plot(IC7.p)
+
+(IC7.a <- radiusMinimaxIC(L2Fam=P, neighbor=ContNeighborhood(),
+ risk=asMSE(biastype=asymmetricBias(nu = c(1,0.2))), loRad=0, upRad=0.5))
+checkIC(IC7.a)
+Risks(IC7.a)
+plot(IC7.a)
+
(IC8 <- radiusMinimaxIC(L2Fam=P, neighbor=TotalVarNeighborhood(),
risk=asMSE(), loRad=0, upRad=Inf))
checkIC(IC8)
Modified: pkg/ROptEst/man/getAsRisk.Rd
===================================================================
--- pkg/ROptEst/man/getAsRisk.Rd 2008-03-28 11:36:49 UTC (rev 81)
+++ pkg/ROptEst/man/getAsRisk.Rd 2008-03-29 04:54:44 UTC (rev 82)
@@ -37,8 +37,8 @@
\S4method{getAsRisk}{asBias,UnivariateDistribution,TotalVarNeighborhood,ANY}(risk, L2deriv, neighbor, biastype, trafo)
-\S4method{getAsRisk}{asBias,RealRandVariable,ContNeighborhood,ANY}(risk, L2deriv, neighbor, biastype, Distr, L2derivDistrSymm, trafo,
- z.start, A.start, maxiter, tol, norm = EuclideanNorm)
+\S4method{getAsRisk}{asBias,RealRandVariable,ContNeighborhood,ANY}(risk, L2deriv, neighbor, biastype, Distr, DistrSymm, L2derivSymm,
+ L2derivDistrSymm, trafo, z.start, A.start, maxiter, tol, norm = EuclideanNorm)
\S4method{getAsRisk}{asCov,UnivariateDistribution,ContNeighborhood,ANY}(risk, L2deriv, neighbor, biastype, clip, cent, stand)
@@ -67,6 +67,8 @@
\item{stand}{ standardizing matrix. }
\item{trafo}{ matrix: transformation of the parameter. }
\item{Distr}{ object of class \code{"Distribution"}. }
+ \item{DistrSymm}{ object of class \code{"DistributionSymmetry"}. }
+ \item{L2derivSymm}{ object of class \code{"FunSymmList"}. }
\item{L2derivDistrSymm}{ object of class \code{"DistrSymmList"}. }
\item{z.start}{ initial value for the centering constant. }
\item{A.start}{ initial value for the standardizing matrix. }
Modified: pkg/ROptEst/man/getInfRobIC.Rd
===================================================================
--- pkg/ROptEst/man/getInfRobIC.Rd 2008-03-28 11:36:49 UTC (rev 81)
+++ pkg/ROptEst/man/getInfRobIC.Rd 2008-03-29 04:54:44 UTC (rev 82)
@@ -30,8 +30,8 @@
\S4method{getInfRobIC}{UnivariateDistribution,asBias,UncondNeighborhood}(L2deriv, risk, neighbor, symm, trafo,
maxiter, tol)
-\S4method{getInfRobIC}{RealRandVariable,asBias,ContNeighborhood}(L2deriv, risk, neighbor, Distr,
- L2derivDistrSymm, z.start, A.start, trafo, maxiter, tol)
+\S4method{getInfRobIC}{RealRandVariable,asBias,ContNeighborhood}(L2deriv, risk, neighbor, Distr, DistrSymm, L2derivSymm,
+ L2derivDistrSymm, z.start, A.start, Finfo, trafo, maxiter, tol)
\S4method{getInfRobIC}{UnivariateDistribution,asHampel,UncondNeighborhood}(L2deriv, risk, neighbor, symm, Finfo, trafo,
upper, maxiter, tol, warn)
Modified: pkg/ROptEst/man/minmaxBias.Rd
===================================================================
--- pkg/ROptEst/man/minmaxBias.Rd 2008-03-28 11:36:49 UTC (rev 81)
+++ pkg/ROptEst/man/minmaxBias.Rd 2008-03-29 04:54:44 UTC (rev 82)
@@ -29,7 +29,7 @@
maxiter, tol)
\S4method{minmaxBias}{RealRandVariable,ContNeighborhood,BiasType}(L2deriv, neighbor, biastype, Distr,
- L2derivDistrSymm, z.start, A.start, trafo, maxiter, tol)
+ L2derivDistrSymm, z.start, A.start, z.comp, A.comp, trafo, maxiter, tol)
}
\arguments{
@@ -43,6 +43,8 @@
\item{L2derivDistrSymm}{ object of class \code{"DistrSymmList"}. }
\item{z.start}{ initial value for the centering constant. }
\item{A.start}{ initial value for the standardizing matrix. }
+ \item{z.comp}{ \code{logical} indicator which indices need to be computed and which are 0 due to symmetry. }
+ \item{A.comp}{ \code{matrix} of \code{logical} indicator which indices need to be computed and which are 0 due to symmetry.}
\item{trafo}{ matrix: transformation of the parameter. }
\item{maxiter}{ the maximum number of iterations. }
\item{tol}{ the desired accuracy (convergence tolerance).}
Modified: pkg/ROptEst/man/updateNorm-methods.Rd
===================================================================
--- pkg/ROptEst/man/updateNorm-methods.Rd 2008-03-28 11:36:49 UTC (rev 81)
+++ pkg/ROptEst/man/updateNorm-methods.Rd 2008-03-29 04:54:44 UTC (rev 82)
@@ -2,24 +2,19 @@
\docType{methods}
\alias{updateNorm-methods}
\alias{updateNorm}
-\alias{updateNorm,NormType-method}
-\alias{updateNorm,InfoNorm-method}
\alias{updateNorm,SelfNorm-method}
\title{ Methods for Function updateNorm in Package `ROptEst' }
\description{updateNorm-methods to update norm in IC-Algo}
\usage{updateNorm(normtype, ...)
-\S4method{updateNorm}{NormType}(normtype, ...)
-\S4method{updateNorm}{InfoNorm}(normtype, FI, ...)
\S4method{updateNorm}{SelfNorm}(normtype, L2, neighbor, biastype, Distr, V.comp,
- cent, stand, w, ...)
+ cent, stand, w)
}
\arguments{
\item{normtype}{normtype of class \code{NormType}}
\item{\dots}{ further arguments to be passed to specific methods.}
- \item{FI}{matrix: Fisher Information}
\item{L2}{L2derivative}
\item{neighbor}{ object of class \code{"Neighborhood"}. }
\item{biastype}{ object of class \code{"BiasType"} }
@@ -31,10 +26,6 @@
\item{w}{object of class \code{RobWeight}; current weight}
}
\section{Methods}{\describe{
-\item{updateNorm}{\code{signature(normtype = "NormType")}: leaves the norm unchanged;}
-\item{updateNorm}{\code{signature(normtype = "InfoNorm")}:
- udates the norm in the information-standardized case; just used
- internally in the opt-IC-Algorithm. }
\item{updateNorm}{\code{signature(normtype = "SelfNorm")}:
udates the norm in the self-standardized case; just used
internally in the opt-IC-Algorithm. }
@@ -44,8 +35,8 @@
}
\details{\code{updateNorm} is used internally in the opt-IC-algorithm to be
- able to work with a norm that depends on the Fisher information at a certain
- parameter (\code{InfoType}) or on the current covariance (\code{SelfNorm})}
+ able to work with a norm that depends on the current covariance
+ (\code{SelfNorm})}
\author{Peter Ruckdeschel \email{Peter.Ruckdeschel at uni-bayreuth.de}}
\seealso{\code{\link[distrMod]{NormType-class}}}
%\examples{}
Modified: pkg/RobAStBase/NAMESPACE
===================================================================
--- pkg/RobAStBase/NAMESPACE 2008-03-28 11:36:49 UTC (rev 81)
+++ pkg/RobAStBase/NAMESPACE 2008-03-29 04:54:44 UTC (rev 82)
@@ -41,12 +41,10 @@
exportMethods("oneStepEstimator",
"locMEstimator")
exportMethods("weight", "weight<-",
- "clip", "clip<-",
- "stand", "stand<-",
- "cent", "cent<-",
- "getweight", "minbiasweight",
+ "getweight",
+ "minbiasweight",
"generateIC.fct",
- "makeIC")
+ "makeIC", "normtype", "biastype")
exportMethods("getRiskIC")
exportMethods("getBiasIC")
export("ContNeighborhood", "TotalVarNeighborhood")
Modified: pkg/RobAStBase/R/TotalVarIC.R
===================================================================
--- pkg/RobAStBase/R/TotalVarIC.R 2008-03-28 11:36:49 UTC (rev 81)
+++ pkg/RobAStBase/R/TotalVarIC.R 2008-03-29 04:54:44 UTC (rev 82)
@@ -3,6 +3,7 @@
Curve = EuclRandVarList(RealRandVariable(Map = c(function(x){x}), Domain = Reals())),
Risks, Infos, clipLo = -Inf, clipUp = Inf, stand = as.matrix(1),
lowerCase = NULL, neighborRadius = 0, w = new("BdStWeight")){
+
if(missing(name))
name <- "IC of total variation type"
if(missing(Risks))
@@ -41,15 +42,16 @@
L2Fam = "L2ParamFamily"),
function(neighbor, L2Fam, res){
A <- res$A
- a <- sign(as.vector(A))*res$a
+ clipLo <- sign(as.vector(A))*res$a
b <- res$b
w <- res$w
ICfct <- vector(mode = "list", length = 1)
Y <- as(A %*% L2Fam at L2deriv, "EuclRandVariable")
- if((a == -Inf) & (b == Inf))
+ if((clipLo == -Inf) & (b == Inf))
clipUp <- Inf
else
- clipUp <- a + b
+ clipUp <- clipLo + b
+
return(TotalVarIC(
name = "IC of total variation type",
CallL2Fam = call("L2ParamFamily",
@@ -68,10 +70,10 @@
FisherInfo.fct = L2Fam at FisherInfo.fct),
Curve = generateIC.fct(neighbor, L2Fam, res),
clipUp = clipUp,
- clipLo = a,
+ clipLo = clipLo,
stand = A,
lowerCase = res$d,
- weight = w,
+ w = w,
neighborRadius = neighbor at radius,
Risks = res$risk,
Infos = matrix(res$info, ncol = 2,
Modified: pkg/RobAStBase/R/Weights.R
===================================================================
--- pkg/RobAStBase/R/Weights.R 2008-03-28 11:36:49 UTC (rev 81)
+++ pkg/RobAStBase/R/Weights.R 2008-03-29 04:54:44 UTC (rev 82)
@@ -27,13 +27,13 @@
setMethod("getweight",
signature(Weight = "HampelWeight", neighbor = "ContNeighborhood",
biastype = "BiasType"),# normtype = "NormType"),
- function(Weight, neighbor, biastype, normtype)
+ function(Weight, neighbor, biastype, normW)
{A <- stand(Weight)
b <- clip(Weight)
z <- cent(Weight)
function(x){
y <- A%*%(x-z)
- norm0 <- fct(normtype)(y)
+ norm0 <- fct(normW)(y)
ind2 <- (norm0 < b/2)
norm1 <- ind2*b/2 + (1-ind2)*norm0
ind1 <- (norm0 < b)
@@ -83,9 +83,8 @@
function(Weight, neighbor, biastype, ...)
{A <- stand(Weight)
b <- clip(Weight)
- a <- A * cent(Weight)
- b1 <- a
- b2 <- b + a
+ b1 <- -b[1]
+ b2 <- b[2]
function(x){
y <- A*x
norm1 <- pmax(-y,b1/2)
@@ -98,13 +97,13 @@
setMethod(minbiasweight,
signature(Weight = "HampelWeight", neighbor = "ContNeighborhood",
biastype = "BiasType"),# norm = "NormType"),
- function(Weight, neighbor, biastype, normtype)
+ function(Weight, neighbor, biastype, normW)
{A <- stand(Weight)
b <- clip(Weight)
z <- cent(Weight)
function(x){
y <- A%*%(x-z)
- norm0 <- fct(normtype)(y)
+ norm0 <- fct(normW)(y)
ind <- 1-.eq(norm0)
ind*b/(norm0+1-ind)
}
@@ -118,16 +117,15 @@
function(Weight, neighbor, biastype, ...)
{A <- stand(Weight)
b <- clip(Weight)
- b1 <- b/nu(biastype)[1]
- b2 <- b/nu(biastype)[2]
+ b1 <- -b[1]
+ b2 <- b[2]
z <- cent(Weight)
function(x){
y <- A*(x-z)
- norm1 <- abs(y[y>0])
- norm2 <- abs(y[y<0])
- ind1 <- 1-.eq(norm1)
- ind2 <- 1-.eq(norm2)
- ind1*b1/(norm1+1-ind1) + ind2*b2/(norm2+1-ind2)
+ indp <- (y>0)
+ ind0 <- .eq(y)
+ indm <- (y<0)
+ indm*b1/(y+ind0) + indp*b2/(y+ind0)
}
}
)
@@ -141,9 +139,9 @@
z <- cent(Weight)
function(x){
y <- A*(x-z)
- norm1 <- abs(y[y * sign(biastype)>0])
- ind1 <- 1-.eq(norm1)
- ind1*b1/(norm1+1-ind1)
+ ind <- (y*sign(biastype) >0)
+ ind0 <- .eq(y)
+ ind*b/(y+ind0)+(1-ind)
}
}
)
@@ -151,20 +149,18 @@
setMethod(minbiasweight,
signature(Weight = "BdStWeight", neighbor = "TotalVarNeighborhood",
- biastype = "BiasType"),# norm = "missing"),
+ biastype = "BiasType"),
function(Weight, neighbor, biastype, ...)
{A <- stand(Weight)
b <- clip(Weight)
- a <- A * cent(Weight)
- b1 <- a
- b2 <- b + a
+ b1 <- b[1]
+ b2 <- b[2]
function(x){
y <- A*x
- norm1 <- abs(y[y>0])
- norm2 <- abs(y[y<0])
- ind1 <- 1-.eq(norm1)
- ind2 <- 1-.eq(norm2)
- ind1*b1/(norm1+1-ind1) + ind2*b2/(norm2+1-ind2)
+ indp <- (y>0)
+ ind0 <- .eq(y)
+ indm <- (y<0)
+ indm*b1/(y+ind0) + indp*b2/(y+ind0)
}
}
)
Modified: pkg/RobAStBase/R/generateICfct.R
===================================================================
--- pkg/RobAStBase/R/generateICfct.R 2008-03-28 11:36:49 UTC (rev 81)
+++ pkg/RobAStBase/R/generateICfct.R 2008-03-29 04:54:44 UTC (rev 82)
@@ -2,8 +2,8 @@
## for internal use only!
setMethod("generateIC.fct", signature(neighbor = "UncondNeighborhood", L2Fam = "L2ParamFamily"),
function(neighbor, L2Fam, res){
- A <- res$A
- a <- res$a
+ A <- as.matrix(res$A)
+ a <- if(is(neighbor,"TotalVarNeighborhood")) 0 else res$a
b <- res$b
d <- res$d
w <- weight(res$w)
@@ -12,34 +12,35 @@
ICfct <- vector(mode = "list", length = nrvalues)
Y <- as(A %*% L2Fam at L2deriv - a, "EuclRandVariable")
L <- as(diag(dim)%*%L2Fam at L2deriv, "EuclRandVariable")
+ L.fct <- function(x) evalRandVar(L,x)
if(nrvalues == 1){
if(!is.null(d)){
ICfct[[1]] <- function(x){}
body(ICfct[[1]]) <- substitute(
{ ind <- 1-.eq(Y(x))
- b*(Y(x)*w(L(x)) + zi*(1-ind)*d) },
- list(Y = Y at Map[[1]], L = L at Map[[1]], w = w, b = b, d = d,
+ Y(x)*w(L(x)) + zi*(1-ind)*d*b },
+ list(Y = Y at Map[[1]], L = L.fct, w = w, b = b, d = d,
zi = sign(L2Fam at param@trafo), .eq = .eq))
}else{
ICfct[[1]] <- function(x){}
body(ICfct[[1]]) <- substitute({ Y(x)*w(L(x)) },
- list(Y = Y at Map[[1]], L = L at Map[[1]], w = w))
+ list(Y = Y at Map[[1]], L = L.fct, w = w))
}
}else{
if(!is.null(d))
for(i in 1:nrvalues){
ICfct[[i]] <- function(x){}
- body(ICfct[[i]]) <- substitute({ind <- 1-.eq(Y(x))
- ind*b*Yi(x)*w(L(x)) + (1-ind)*d
+ body(ICfct[[i]]) <- substitute({ind <- 1-.eq(Yi(x))
+ ind*Yi(x)*w(L(x)) + (1-ind)*d
},
- list(Yi = Y at Map[[i]], L = L at Map[[1]], w = w,
- b = b, d = d[i], .eq = .eq))
+ list(Yi = Y at Map[[i]], L = L.fct, w = w,
+ b = b, d = d[i]))#, .eq = .eq))
}
else
for(i in 1:nrvalues){
ICfct[[i]] <- function(x){}
body(ICfct[[i]]) <- substitute({ Yi(x)*w(L(x)) },
- list(Yi = Y at Map[[i]], L = L at Map[[1]], w = w))
+ list(Yi = Y at Map[[i]], L = L.fct, w = w))
}
}
return(EuclRandVarList(EuclRandVariable(Map = ICfct, Domain = Y at Domain,
Modified: pkg/RobAStBase/R/getRiskIC.R
===================================================================
--- pkg/RobAStBase/R/getRiskIC.R 2008-03-28 11:36:49 UTC (rev 81)
+++ pkg/RobAStBase/R/getRiskIC.R 2008-03-29 04:54:44 UTC (rev 82)
@@ -93,14 +93,17 @@
neighbor = "UncondNeighborhood",
L2Fam = "missing"),
function(IC, risk, neighbor, tol = .Machine$double.eps^0.25){
- getBiasIC(IC, neighbor, biastype(risk), normtype(risk), tol)
+ getBiasIC(IC = IC, neighbor = neighbor,
+ biastype = biastype(risk), normtype = normtype(risk), tol = tol)
})
setMethod("getRiskIC", signature(IC = "IC",
risk = "asBias",
neighbor = "UncondNeighborhood",
L2Fam = "L2ParamFamily"),
function(IC, risk, neighbor, L2Fam, tol = .Machine$double.eps^0.25){
- getBiasIC(IC, neighbor, L2Fam, biastype(risk), normtype(risk), tol)
+ getBiasIC(IC = IC, neighbor = neighbor, L2Fam = L2Fam,
+ biastype = biastype(risk), normtype = normtype(risk),
+ tol = tol)
})
###############################################################################
## asymptotic MSE
@@ -110,10 +113,11 @@
neighbor = "UncondNeighborhood",
L2Fam = "missing"),
function(IC, risk, neighbor, tol = .Machine$double.eps^0.25){
- L2fam <- eval(IC at CallL2Fam)
+ L2Fam <- eval(IC at CallL2Fam)
getRiskIC(IC = IC, risk = risk, neighbor = neighbor,
L2Fam = L2Fam, tol = tol)
})
+
setMethod("getRiskIC", signature(IC = "IC",
risk = "asMSE",
neighbor = "UncondNeighborhood",
@@ -126,8 +130,7 @@
if(rad == Inf) return(Inf)
trCov <- getRiskIC(IC = IC, risk = trAsCov(), L2Fam = L2Fam)
- Bias <- getRiskIC(IC = IC, risk = asBias(), neighbor = neighbor, L2Fam = L2Fam,
- biastype = biastype(risk))
+ Bias <- getRiskIC(IC = IC, risk = asBias(), neighbor = neighbor, L2Fam = L2Fam)
prec <- checkIC(IC, L2Fam, out = FALSE)
if(prec > tol)
Modified: pkg/RobAStBase/R/infoPlot.R
===================================================================
--- pkg/RobAStBase/R/infoPlot.R 2008-03-28 11:36:49 UTC (rev 81)
+++ pkg/RobAStBase/R/infoPlot.R 2008-03-29 04:54:44 UTC (rev 82)
@@ -22,12 +22,33 @@
trafo <- L2Fam at param@trafo
dims <- nrow(trafo)
+
+ QFc <- diag(dims)
+ if(is(object,"ContIC") & dims>1 )
+ {if (is(normtype(object),"QFNorm")) QFc <- QuadForm(normtype(object))
+ QFc0 <- solve( trafo %*% solve(L2Fam at FisherInfo) %*% t(trafo ))
+ if (is(normtype(object),"SelfNorm")|is(normtype(object),"InfoNorm"))
+ QFc <- QFc0
+ }
+ print(QFc)
+ QFc.5 <- sqrt(PosSemDefSymmMatrix(QFc))
+ print(QFc.5)
+
classIC <- as(trafo %*% solve(L2Fam at FisherInfo) %*% L2Fam at L2deriv, "EuclRandVariable")
- absInfoClass <- classIC %*% classIC
+ absInfoClass <- t(classIC) %*% QFc %*% classIC
absInfoClass <- sapply(x.vec, absInfoClass at Map[[1]])
+
+ QF <- diag(dims)
+ if(is(object,"ContIC") & dims>1 )
+ {if (is(normtype(object),"QFNorm")) QF <- QuadForm(normtype(object))}
+ print(QF)
+ QF.5 <- sqrt(PosSemDefSymmMatrix(QF))
+ print(QF.5)
+
IC1 <- as(diag(dims) %*% object at Curve, "EuclRandVariable")
- absInfo <- IC1 %*% IC1
+ absInfo <- t(IC1) %*% QF %*% IC1
absInfo <- sapply(x.vec, absInfo at Map[[1]])
+
plot(x.vec, absInfoClass, type = plty, lty = "dashed",
ylim = c(0, 2*max(absInfo, na.rm = TRUE)), xlab = "x",
ylab = "absolute information", col = grey(0.5))
@@ -52,11 +73,15 @@
opar <- par()
get(getOption("device"))()
par(mfrow = c(nrows, ncols))
+ IC1.i.5 <- QF.5%*%IC1
+ classIC.i.5 <- QFc.5%*%classIC
for(i in 1:dims){
- y.vec <- sapply(x.vec, IC1 at Map[[i]])^2/absInfo
+ y.vec <- sapply(x.vec, IC1.i.5 at Map[[i]])^2/absInfo
plot(x.vec, y.vec, type = plty, lty = lty, lwd = 2,
xlab = "x", ylab = "relative information", ylim = c(0, 1.1))
- lines(x.vec, sapply(x.vec, classIC at Map[[i]])^2/absInfoClass, type = plty,
+
+ yc.vec <- sapply(x.vec, classIC.i.5 at Map[[i]])^2/absInfoClass
+ lines(x.vec, yc.vec, type = plty,
lty = "dashed", col = grey(0.5))
legend(max(x.vec), 1.1, xjust = 1, cex = 0.6,
legend = c("class. opt. IC"), lty = "dashed", col = c(grey(0.5)))
Modified: pkg/RobAStBase/chm/RobAStBase.chm
===================================================================
(Binary files differ)
Modified: pkg/RobAStBase/chm/getweight.html
===================================================================
--- pkg/RobAStBase/chm/getweight.html 2008-03-28 11:36:49 UTC (rev 81)
+++ pkg/RobAStBase/chm/getweight.html 2008-03-29 04:54:44 UTC (rev 82)
@@ -38,31 +38,34 @@
minbiasweight(Weight, neighbor, biastype, ...)
## S4 method for signature 'HampelWeight, ContNeighborhood,
## BiasType':
-getweight(Weight, neighbor, biastype, normtype)
+getweight(Weight, neighbor, biastype, normW)
## S4 method for signature 'HampelWeight, ContNeighborhood,
## BiasType':
-minbiasweight(Weight, neighbor, biastype, normtype)
+minbiasweight(Weight, neighbor, biastype, normW)
## S4 method for signature 'HampelWeight, ContNeighborhood,
## onesidedBias':
getweight(Weight, neighbor, biastype, ...)
## S4 method for signature 'HampelWeight, ContNeighborhood,
## onesidedBias':
-minbiasweight(Weight, neighbor, biastype,...)
+minbiasweight(Weight, neighbor, biastype, ...)
## S4 method for signature 'HampelWeight, ContNeighborhood,
## asymmetricBias':
getweight(Weight, neighbor, biastype, ...)
## S4 method for signature 'HampelWeight, ContNeighborhood,
## asymmetricBias':
-minbiasweight(Weight, neighbor, biastype,...)
+minbiasweight(Weight, neighbor, biastype, ...)
+## S4 method for signature 'BdStWeight,
+## TotalVarNeighborhood, BiasType':
+getweight(Weight, neighbor, biastype, ...)
+## S4 method for signature 'BdStWeight,
+## TotalVarNeighborhood, BiasType':
+minbiasweight(Weight, neighbor, biastype, ...)
</pre>
<h3>Arguments</h3>
<table summary="R argblock">
-<tr valign="top"><td><code>...</code></td>
-<td>
-additional arguments </td></tr>
<tr valign="top"><td><code>Weight</code></td>
<td>
Object of class <code>"RobWeight"</code>. </td></tr>
@@ -72,7 +75,7 @@
<tr valign="top"><td><code>biastype</code></td>
<td>
Object of class <code>"BiasType"</code>. </td></tr>
-<tr valign="top"><td><code>normtype</code></td>
+<tr valign="top"><td><code>normW</code></td>
<td>
Object of class <code>"NormType"</code> — only for signature <code>HampelWeight,ContNeighborhood,BiasType</code>. </td></tr>
<tr valign="top"><td><code>...</code></td>
Modified: pkg/RobAStBase/man/getweight.Rd
===================================================================
--- pkg/RobAStBase/man/getweight.Rd 2008-03-28 11:36:49 UTC (rev 81)
+++ pkg/RobAStBase/man/getweight.Rd 2008-03-29 04:54:44 UTC (rev 82)
@@ -20,19 +20,20 @@
\usage{
getweight(Weight, neighbor, biastype, ...)
minbiasweight(Weight, neighbor, biastype, ...)
-\S4method{getweight}{HampelWeight,ContNeighborhood,BiasType}(Weight, neighbor, biastype, normtype)
-\S4method{minbiasweight}{HampelWeight,ContNeighborhood,BiasType}(Weight, neighbor, biastype, normtype)
+\S4method{getweight}{HampelWeight,ContNeighborhood,BiasType}(Weight, neighbor, biastype, normW)
+\S4method{minbiasweight}{HampelWeight,ContNeighborhood,BiasType}(Weight, neighbor, biastype, normW)
\S4method{getweight}{HampelWeight,ContNeighborhood,onesidedBias}(Weight, neighbor, biastype, ...)
-\S4method{minbiasweight}{HampelWeight,ContNeighborhood,onesidedBias}(Weight, neighbor, biastype,...)
+\S4method{minbiasweight}{HampelWeight,ContNeighborhood,onesidedBias}(Weight, neighbor, biastype, ...)
\S4method{getweight}{HampelWeight,ContNeighborhood,asymmetricBias}(Weight, neighbor, biastype, ...)
-\S4method{minbiasweight}{HampelWeight,ContNeighborhood,asymmetricBias}(Weight, neighbor, biastype,...)
+\S4method{minbiasweight}{HampelWeight,ContNeighborhood,asymmetricBias}(Weight, neighbor, biastype, ...)
+\S4method{getweight}{BdStWeight,TotalVarNeighborhood,BiasType}(Weight, neighbor, biastype, ...)
+\S4method{minbiasweight}{BdStWeight,TotalVarNeighborhood,BiasType}(Weight, neighbor, biastype, ...)
}
\arguments{
- \item{\dots}{ additional arguments }
\item{Weight}{ Object of class \code{"RobWeight"}. }
\item{neighbor}{ Object of class \code{"Neighborhood"}. }
\item{biastype}{ Object of class \code{"BiasType"}. }
- \item{normtype}{ Object of class \code{"NormType"} --- only for signature \code{HampelWeight,ContNeighborhood,BiasType}. }
+ \item{normW}{ Object of class \code{"NormType"} --- only for signature \code{HampelWeight,ContNeighborhood,BiasType}. }
\item{\dots}{possibly additional (unused) arguments --- like in a call to the less specific methods.}
}
%\details{}
More information about the Robast-commits
mailing list