[Robast-commits] r437 - in pkg: ROptEst ROptEst/R ROptEst/inst ROptEst/inst/scripts ROptEst/man ROptEstOld ROptEstOld/inst ROptRegTS ROptRegTS/inst RandVar RandVar/R RandVar/inst RandVar/man RobAStBase RobAStBase/R RobAStBase/inst RobAStBase/man RobLox RobLox/R RobLox/inst RobLox/man RobLoxBioC RobLoxBioC/R RobLoxBioC/inst RobLoxBioC/man RobRex RobRex/inst
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Dec 2 15:24:43 CET 2010
Author: ruckdeschel
Date: 2010-12-02 15:24:42 +0100 (Thu, 02 Dec 2010)
New Revision: 437
Added:
pkg/ROptEst/R/AllShow.R
pkg/ROptEst/R/asAnscombe.R
pkg/ROptEst/R/asL1L4.R
pkg/ROptEst/R/getAsGRiskfct.R
pkg/ROptEst/R/getInfRobIC_asAnscombe.R
pkg/ROptEst/R/getMaxIneff.R
pkg/ROptEst/R/getReq.R
pkg/ROptEst/inst/scripts/AnscombeOrNot.R
pkg/ROptEst/inst/scripts/MBRE.R
pkg/ROptEst/inst/scripts/NbinomModel.R
pkg/ROptEst/inst/scripts/NormalLocationModel.R
pkg/ROptEst/man/asAnscombe-class.Rd
pkg/ROptEst/man/asAnscombe.Rd
pkg/ROptEst/man/asL1-class.Rd
pkg/ROptEst/man/asL1.Rd
pkg/ROptEst/man/asL4-class.Rd
pkg/ROptEst/man/asL4.Rd
pkg/ROptEst/man/getAsGRiskFct-methods.Rd
pkg/ROptEst/man/getMaxIneff.Rd
pkg/ROptEst/man/getReq.Rd
Modified:
pkg/ROptEst/DESCRIPTION
pkg/ROptEst/NAMESPACE
pkg/ROptEst/R/AllClass.R
pkg/ROptEst/R/AllGeneric.R
pkg/ROptEst/R/LowerCaseMultivariate.R
pkg/ROptEst/R/cniperCont.R
pkg/ROptEst/R/getAsRisk.R
pkg/ROptEst/R/getInfClip.R
pkg/ROptEst/R/getInfGamma.R
pkg/ROptEst/R/getInfLM.R
pkg/ROptEst/R/getInfRobIC_asBias.R
pkg/ROptEst/R/getInfRobIC_asGRisk.R
pkg/ROptEst/R/getInfV.R
pkg/ROptEst/R/leastFavorableRadius.R
pkg/ROptEst/R/radiusMinimaxIC.R
pkg/ROptEst/inst/NEWS
pkg/ROptEst/inst/scripts/BinomialModel.R
pkg/ROptEst/inst/scripts/GammaModel.R
pkg/ROptEst/inst/scripts/NormalLocationScaleModel.R
pkg/ROptEst/inst/scripts/NormalScaleModel.R
pkg/ROptEst/inst/scripts/PoissonModel.R
pkg/ROptEst/man/0ROptEst-package.Rd
pkg/ROptEst/man/cniperCont.Rd
pkg/ROptEst/man/getAsRisk.Rd
pkg/ROptEst/man/getInfClip.Rd
pkg/ROptEst/man/getInfGamma.Rd
pkg/ROptEst/man/getInfRobIC.Rd
pkg/ROptEst/man/internals.Rd
pkg/ROptEst/man/optIC.Rd
pkg/ROptEst/man/roptest.Rd
pkg/ROptEstOld/DESCRIPTION
pkg/ROptEstOld/inst/NEWS
pkg/ROptRegTS/DESCRIPTION
pkg/ROptRegTS/inst/NEWS
pkg/RandVar/DESCRIPTION
pkg/RandVar/R/util.R
pkg/RandVar/inst/NEWS
pkg/RandVar/man/0RandVar-package.Rd
pkg/RobAStBase/DESCRIPTION
pkg/RobAStBase/R/AllPlot.R
pkg/RobAStBase/R/IC.R
pkg/RobAStBase/R/comparePlot.R
pkg/RobAStBase/R/ddPlot.R
pkg/RobAStBase/R/ddPlot_utils.R
pkg/RobAStBase/R/getRiskIC.R
pkg/RobAStBase/R/infoPlot.R
pkg/RobAStBase/R/qqplot.R
pkg/RobAStBase/inst/NEWS
pkg/RobAStBase/man/0RobAStBase-package.Rd
pkg/RobAStBase/man/IC.Rd
pkg/RobAStBase/man/comparePlot.Rd
pkg/RobAStBase/man/cutoff-class.Rd
pkg/RobAStBase/man/cutoff.Rd
pkg/RobAStBase/man/ddPlot-methods.Rd
pkg/RobAStBase/man/infoPlot.Rd
pkg/RobAStBase/man/internals_ddPlot.Rd
pkg/RobAStBase/man/makeIC-methods.Rd
pkg/RobAStBase/man/qqplot.Rd
pkg/RobLox/DESCRIPTION
pkg/RobLox/R/roblox.R
pkg/RobLox/inst/NEWS
pkg/RobLox/man/0RobLox-package.Rd
pkg/RobLoxBioC/DESCRIPTION
pkg/RobLoxBioC/R/AffySimStudyFunction.R
pkg/RobLoxBioC/R/IlluminaSimStudyFunction.R
pkg/RobLoxBioC/inst/NEWS
pkg/RobLoxBioC/man/0RobLoxBioC-package.Rd
pkg/RobRex/DESCRIPTION
pkg/RobRex/inst/NEWS
Log:
merged branch robast-2.3 into trunk
Modified: pkg/ROptEst/DESCRIPTION
===================================================================
--- pkg/ROptEst/DESCRIPTION 2010-12-02 13:59:11 UTC (rev 436)
+++ pkg/ROptEst/DESCRIPTION 2010-12-02 14:24:42 UTC (rev 437)
@@ -1,14 +1,18 @@
Package: ROptEst
-Version: 0.7
-Date: 2009-10-16
+Version: 0.8
+Date: 2010-12-02
Title: Optimally robust estimation
-Description: Optimally robust estimation in general smoothly parameterized models using S4 classes and methods.
-Depends: R(>= 2.7.0), methods, distr(>= 2.0), distrEx(>= 2.0), distrMod(>= 2.0), RandVar(>= 0.6.4), RobAStBase
+Description: Optimally robust estimation in general smoothly
+ parameterized models using S4 classes and methods.
+Depends: R(>= 2.7.0), methods, distr(>= 2.0), distrEx(>= 2.0),
+ distrMod(>= 2.0), RandVar(>= 0.6.4), RobAStBase
Author: Matthias Kohl, Peter Ruckdeschel
Maintainer: Matthias Kohl <Matthias.Kohl at stamats.de>
LazyLoad: yes
License: LGPL-3
URL: http://robast.r-forge.r-project.org/
Encoding: latin1
-LastChangedDate: {$LastChangedDate$}
+LastChangedDate: {$LastChangedDate: 2009-11-01 11:47:36 +0100 (So, 01
+ Nov 2009) $}
LastChangedRevision: {$LastChangedRevision$}
+SVNRevision: 436
Modified: pkg/ROptEst/NAMESPACE
===================================================================
--- pkg/ROptEst/NAMESPACE 2010-12-02 13:59:11 UTC (rev 436)
+++ pkg/ROptEst/NAMESPACE 2010-12-02 14:24:42 UTC (rev 437)
@@ -4,6 +4,7 @@
import("distrMod")
import("RobAStBase")
+exportClasses("asAnscombe", "asL1", "asL4")
exportMethods("optIC",
"getInfRobIC",
"getFixRobIC",
@@ -25,6 +26,9 @@
"getL1normL2deriv",
"getModifyIC",
"cniperCont", "cniperPoint", "cniperPointPlot")
-exportMethods("updateNorm", "scaleUpdateIC")
-export("getL2normL2deriv")
+exportMethods("updateNorm", "scaleUpdateIC", "eff",
+ "get.asGRisk.fct")
+export("getL2normL2deriv",
+ "asAnscombe", "asL1", "asL4",
+ "getReq", "getMaxIneff")
export("roptest","getLagrangeMultByOptim","getLagrangeMultByIter")
Modified: pkg/ROptEst/R/AllClass.R
===================================================================
--- pkg/ROptEst/R/AllClass.R 2010-12-02 13:59:11 UTC (rev 436)
+++ pkg/ROptEst/R/AllClass.R 2010-12-02 14:24:42 UTC (rev 437)
@@ -8,3 +8,21 @@
}
+## asymptotic Anscombe risk
+setClass("asAnscombe", representation(eff = "numeric"),
+ prototype = prototype(eff = .95,
+ type = "optimal bias robust IC for given ARE in the ideal model"),
+ contains = "asRiskwithBias",
+ validity = function(object){
+ if(any(object at eff <= 0|object at eff > 1))
+ stop("'eff' has to be in (0,1]")
+ else TRUE
+ })
+
+
+## asymptotic L4 error
+setClass("asL4", contains = "asGRisk",
+ prototype = prototype(type = "asymptotic mean power 4 error"))
+## asymptotic L1 error
+setClass("asL1", contains = "asGRisk",
+ prototype = prototype(type = "asymptotic mean absolute error"))
Modified: pkg/ROptEst/R/AllGeneric.R
===================================================================
--- pkg/ROptEst/R/AllGeneric.R 2010-12-02 13:59:11 UTC (rev 436)
+++ pkg/ROptEst/R/AllGeneric.R 2010-12-02 14:24:42 UTC (rev 437)
@@ -86,3 +86,9 @@
if(!isGeneric("cniperPointPlot")){
setGeneric("cniperPointPlot", function(L2Fam, neighbor, risk, ...) standardGeneric("cniperPointPlot"))
}
+if(!isGeneric("eff")){
+ setGeneric("eff", function(object) standardGeneric("eff"))
+}
+if(!isGeneric("get.asGRisk.fct")){
+ setGeneric("get.asGRisk.fct", function(Risk) standardGeneric("get.asGRisk.fct"))
+}
Copied: pkg/ROptEst/R/AllShow.R (from rev 436, branches/robast-0.8/pkg/ROptEst/R/AllShow.R)
===================================================================
--- pkg/ROptEst/R/AllShow.R (rev 0)
+++ pkg/ROptEst/R/AllShow.R 2010-12-02 14:24:42 UTC (rev 437)
@@ -0,0 +1,6 @@
+setMethod("show", "asAnscombe",
+ function(object){
+ cat(paste("An object of class", dQuote(class(object)), "\n"))
+ cat("risk type:\t", object at type, "\n")
+ cat("ARE in the ideal model:\t", object at eff, "\n")
+ })
Modified: pkg/ROptEst/R/LowerCaseMultivariate.R
===================================================================
--- pkg/ROptEst/R/LowerCaseMultivariate.R 2010-12-02 13:59:11 UTC (rev 436)
+++ pkg/ROptEst/R/LowerCaseMultivariate.R 2010-12-02 14:24:42 UTC (rev 437)
@@ -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){
+# 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))
}
Copied: pkg/ROptEst/R/asAnscombe.R (from rev 436, branches/robast-0.8/pkg/ROptEst/R/asAnscombe.R)
===================================================================
--- pkg/ROptEst/R/asAnscombe.R (rev 0)
+++ pkg/ROptEst/R/asAnscombe.R 2010-12-02 14:24:42 UTC (rev 437)
@@ -0,0 +1,8 @@
+## generating function
+asAnscombe <- function(eff = .95, biastype = symmetricBias(),
+ normtype = NormType()){
+ new("asAnscombe", eff = eff, biastype = biastype,
+ normtype = normtype) }
+
+## access method
+setMethod("eff", "asAnscombe", function(object) object at eff)
Copied: pkg/ROptEst/R/asL1L4.R (from rev 436, branches/robast-0.8/pkg/ROptEst/R/asL1L4.R)
===================================================================
--- pkg/ROptEst/R/asL1L4.R (rev 0)
+++ pkg/ROptEst/R/asL1L4.R 2010-12-02 14:24:42 UTC (rev 437)
@@ -0,0 +1,11 @@
+## generating function
+
+asL4 <-
+## The function is currently defined as
+function(biastype = symmetricBias(), normtype = NormType()){
+ new("asL4", biastype = biastype, normtype = normtype) }
+
+asL1 <-
+## The function is currently defined as
+function(biastype = symmetricBias(), normtype = NormType()){
+ new("asL1", biastype = biastype, normtype = normtype) }
Modified: pkg/ROptEst/R/cniperCont.R
===================================================================
--- pkg/ROptEst/R/cniperCont.R 2010-12-02 13:59:11 UTC (rev 436)
+++ pkg/ROptEst/R/cniperCont.R 2010-12-02 14:24:42 UTC (rev 437)
@@ -50,7 +50,9 @@
setMethod("cniperPointPlot", signature(L2Fam = "L2ParamFamily",
neighbor = "ContNeighborhood",
risk = "asMSE"),
- function(L2Fam, neighbor, risk, lower, upper, n = 101){
+ function(L2Fam, neighbor, risk, lower, upper, n = 101, ...){
+ dots <- as.list(match.call(call = sys.call(sys.parent(1)),
+ expand.dots = FALSE)$"...")
D <- trafo(L2Fam at param)
tr.invF <- sum(diag(D %*% solve(FisherInfo(L2Fam)) %*% t(D)))
psi <- optIC(model = L2Fam, risk = asCov())
@@ -61,12 +63,44 @@
y <- evalIC(psi, x)
tr.invF + as.vector(y %*% y)*neighbor at radius^2 - maxMSE
}
- x <- seq(from = lower, to = upper, length = n)
- y <- sapply(x, fun)
- plot(x, y, type = "l", main = "Cniper point plot",
- xlab = "Dirac point", ylab = "Asymptotic MSE difference (classic - robust)")
+ dots$x <- x <- seq(from = lower, to = upper, length = n)
+ dots$y <- sapply(x, fun)
+ colSet <- ltySet <- lwdSet <- FALSE
+ if(!is.null(dots$col)) {colSet <- TRUE; colo <- eval(dots$col)}
+ if(colSet) {
+ colo <- rep(colo,length.out=2)
+ dots$col <- colo[1]
+ }
+ if(!is.null(dots$lwd)) {lwdSet <- TRUE; lwdo <- eval(dots$lwd)}
+ if(lwdSet) {
+ lwdo <- rep(lwdo,length.out=2)
+ dots$lwd <- lwdo[1]
+ }
+ if(!is.null(dots$lty)) {ltySet <- TRUE; ltyo <- eval(dots$lty)}
+ if(ltySet && ((!is.numeric(ltyo) && length(ltyo)==1)||
+ is.numeric(ltyo))){
+ ltyo <- list(ltyo,ltyo)
+ dots$lty <- ltyo[[1]]
+ }else{ if (ltySet && !is.numeric(ltyo) && length(ltyo)==2){
+ dots$lty <- ltyo[[1]]
+ }
+ }
+ if(is.null(dots$main)) dots$main <- gettext("Cniper point plot")
+ if(is.null(dots$xlab)) dots$xlab <- gettext("Dirac point")
+ if(is.null(dots$ylab))
+ dots$ylab <- gettext("Asymptotic MSE difference (classic - robust)")
+ dots$type <- "l"
+ do.call(plot, dots)
# text(min(x), max(y)/2, "Robust", pos = 4)
# text(min(x), min(y)/2, "Classic", pos = 4)
- abline(h = 0)
+ dots$x <- dots$y <- dots$xlab <- dots$ylab <- dots$main <- dots$type <- NULL
+ dots$h <- 0
+ if(colSet) dots$col <- colo[2]
+ if(lwdSet) dots$lwd <- lwdo[2]
+ if(ltySet) dots$lty <- ltyo[[2]]
+ do.call(abline, dots)
invisible()
})
+#
+#cniperPointPlot(L2Fam=N0, neighbor=ContNeighborhood(radius = 0.5), risk=asMSE(),lower=-12, n =30, upper=8, lwd=c(2,4),lty=list(c(5,1),3),col=c(2,4))
+
\ No newline at end of file
Copied: pkg/ROptEst/R/getAsGRiskfct.R (from rev 436, branches/robast-0.8/pkg/ROptEst/R/getAsGRiskfct.R)
===================================================================
--- pkg/ROptEst/R/getAsGRiskfct.R (rev 0)
+++ pkg/ROptEst/R/getAsGRiskfct.R 2010-12-02 14:24:42 UTC (rev 437)
@@ -0,0 +1,14 @@
+
+setMethod("get.asGRisk.fct", signature(Risk = "asL1"),
+ function(Risk){return( function(r,s,b){
+ rb <- r*b; w <- rb/s;
+ rb*(2*pnorm(w)-1)+2*s*dnorm(w)})})
+
+setMethod("get.asGRisk.fct", signature(Risk = "asMSE"),
+ function(Risk){ return(function(r,s,b){
+ rb <- r*b; rb^2+s^2})})
+
+setMethod("get.asGRisk.fct", signature(Risk = "asL4"),
+ function(Risk){ return(function(r,s,b){
+ rb <- r*b; rb^4+6*s^2*rb^2+3*s^4})})
+
Modified: pkg/ROptEst/R/getAsRisk.R
===================================================================
--- pkg/ROptEst/R/getAsRisk.R 2010-12-02 13:59:11 UTC (rev 436)
+++ pkg/ROptEst/R/getAsRisk.R 2010-12-02 14:24:42 UTC (rev 437)
@@ -5,7 +5,7 @@
L2deriv = "UnivariateDistribution",
neighbor = "Neighborhood",
biastype = "ANY"),
- function(risk, L2deriv, neighbor, biastype, clip = NULL, cent = NULL, stand, trafo){
+ function(risk, L2deriv, neighbor, biastype, normtype = NULL, clip = NULL, cent = NULL, stand, trafo, ...){
if(!is.finite(neighbor at radius))
mse <- Inf
else
@@ -17,7 +17,7 @@
L2deriv = "EuclRandVariable",
neighbor = "Neighborhood",
biastype = "ANY"),
- function(risk, L2deriv, neighbor, biastype, clip = NULL, cent = NULL, stand, trafo){
+ function(risk, L2deriv, neighbor, biastype, normtype = NULL, clip = NULL, cent = NULL, stand, trafo, ...){
if(!is.finite(neighbor at radius))
mse <- Inf
else{
@@ -37,7 +37,7 @@
L2deriv = "UnivariateDistribution",
neighbor = "ContNeighborhood",
biastype = "ANY"),
- function(risk, L2deriv, neighbor, biastype, trafo){
+ function(risk, L2deriv, neighbor, biastype, normtype = NULL, clip = NULL, cent = NULL, stand = NULL, trafo, ...){
z <- q(L2deriv)(0.5)
bias <- abs(as.vector(trafo))/E(L2deriv, function(x, z){abs(x - z)},
useApply = FALSE, z = z)
@@ -48,7 +48,7 @@
L2deriv = "UnivariateDistribution",
neighbor = "TotalVarNeighborhood",
biastype = "ANY"),
- function(risk, L2deriv, neighbor, biastype, trafo){
+ function(risk, L2deriv, neighbor, biastype, normtype = NULL, clip = NULL, cent = NULL, stand = NULL, trafo, ...){
bias <- abs(as.vector(trafo))/(-m1df(L2deriv, 0))
return(list(asBias = bias))
@@ -57,9 +57,10 @@
L2deriv = "RealRandVariable",
neighbor = "ContNeighborhood",
biastype = "ANY"),
- function(risk, L2deriv, neighbor, biastype, Distr, DistrSymm, L2derivSymm,
+ function(risk, L2deriv, neighbor, biastype, normtype = NULL,
+ clip = NULL, cent = NULL, stand = NULL, Distr, DistrSymm, L2derivSymm,
L2derivDistrSymm, Finfo, trafo, z.start, A.start, maxiter, tol, warn,
- verbose = NULL){
+ verbose = NULL, ...){
if(missing(verbose)|| is.null(verbose))
verbose <- getRobAStBaseOption("all.verbose")
@@ -95,9 +96,10 @@
L2deriv = "RealRandVariable",
neighbor = "TotalVarNeighborhood",
biastype = "ANY"),
- function(risk, L2deriv, neighbor, biastype, Distr, DistrSymm, L2derivSymm,
+ function(risk, L2deriv, neighbor, biastype, normtype = NULL,
+ clip = NULL, cent = NULL, stand = NULL, Distr, DistrSymm, L2derivSymm,
L2derivDistrSymm, Finfo, trafo, z.start, A.start, maxiter, tol, warn,
- verbose = NULL){
+ verbose = NULL, ...){
if(missing(verbose)|| is.null(verbose))
verbose <- getRobAStBaseOption("all.verbose")
@@ -124,44 +126,49 @@
L2deriv = "UnivariateDistribution",
neighbor = "ContNeighborhood",
biastype = "ANY"),
- function(risk, L2deriv, neighbor, biastype, clip, cent, stand){
+ function(risk, L2deriv, neighbor, biastype, normtype = NULL,clip, cent, stand, trafo = NULL, ...){
# c0 <- clip/abs(as.vector(stand))
# D1 <- L2deriv - cent/as.vector(stand)
# Cov <- (clip^2*(p(D1)(-c0) + 1 - p(D1)(c0))
# + as.vector(stand)^2*(m2df(D1, c0) - m2df(D1, -c0)))
- return(list(asCov =
- getInfV(L2deriv, neighbor, biastype(risk), clip/abs(as.vector(stand)),
+ Cov <- getInfV(L2deriv, neighbor, biastype(risk), clip/abs(as.vector(stand)),
cent/abs(as.vector(stand)), abs(as.vector(stand)))
- ))
+
+ if(!is.null(trafo)) Cov <- trafo%*%Cov%*%t(trafo)
+ return(list(asCov = Cov ))
})
setMethod("getAsRisk", signature(risk = "asCov",
L2deriv = "UnivariateDistribution",
neighbor = "TotalVarNeighborhood",
biastype = "ANY"),
- function(risk, L2deriv, neighbor, biastype, clip, cent, stand){
+ function(risk, L2deriv, neighbor, biastype, normtype = NULL, clip, cent, stand, trafo = NULL, ...){
+ if(missing(biastype)||is.null(biastype)) biastype <- biastype(risk)
# g0 <- cent/abs(as.vector(stand))
# c0 <- clip/abs(as.vector(stand))
# Cov <- (abs(as.vector(stand))^2*(g0^2*p(L2deriv)(g0)
# + (g0+c0)^2*(1 - p(L2deriv)(g0+c0))
# + m2df(L2deriv, g0+c0) - m2df(L2deriv, g0)))
# return(list(asCov = Cov))
- return(list(asCov =
- getInfV(L2deriv, neighbor, biastype(risk), clip/abs(as.vector(stand)),
+ Cov <- getInfV(L2deriv, neighbor, biastype, clip/abs(as.vector(stand)),
cent/abs(as.vector(stand)), abs(as.vector(stand)))
- ))
+ if(!is.null(trafo)) Cov <- trafo%*%Cov%*%t(trafo)
+ return(list(asCov = Cov))
})
setMethod("getAsRisk", signature(risk = "asCov",
L2deriv = "RealRandVariable",
neighbor = "ContNeighborhood",
biastype = "ANY"),
- function(risk, L2deriv, neighbor, biastype, Distr, cent,
- stand, V.comp = matrix(TRUE, ncol = nrow(stand), nrow = nrow(stand)),
- w){
- return(getInfV(L2deriv = L2deriv, neighbor = neighbor,
- biastype = biastype(risk), Distr = Distr,
+ function(risk, L2deriv, neighbor, biastype, normtype = NULL, clip=NULL, cent,
+ stand, Distr, trafo = NULL, V.comp = matrix(TRUE, ncol = nrow(stand), nrow = nrow(stand)),
+ w, ...){
+ if(missing(biastype)||is.null(biastype)) biastype <- biastype(risk)
+ Cov <- getInfV(L2deriv = L2deriv, neighbor = neighbor,
+ biastype = biastype, Distr = Distr,
V.comp = V.comp, cent = cent,
- stand = stand, w = w))
+ stand = stand, w = w)
+ if(!is.null(trafo)) Cov <- trafo%*%Cov%*%t(trafo)
+ return(list(asCov = Cov))
})
# Y <- as(stand %*% L2deriv - cent, "EuclRandVariable")
# absY <- sqrt(Y %*% Y)
@@ -181,6 +188,7 @@
+
###############################################################################
## trace of asymptotic covariance
###############################################################################
@@ -188,20 +196,29 @@
L2deriv = "UnivariateDistribution",
neighbor = "UncondNeighborhood",
biastype = "ANY"),
- function(risk, L2deriv, neighbor, biastype, clip, cent, stand){
+ function(risk, L2deriv, neighbor, biastype, normtype = NULL, clip, cent, stand, trafo=NULL, ...){
+ if(missing(biastype)||is.null(biastype)) biastype <- biastype(risk)
+ if(missing(normtype)||is.null(normtype)) normtype <- normtype(risk)
Cov <- getAsRisk(risk = asCov(), L2deriv = L2deriv, neighbor = neighbor,
- biastype = biastype(risk), clip = clip, cent = cent, stand = stand)$asCov
+ biastype = biastype, clip = clip, cent = cent, stand = stand)$asCov
+ std <- if(is(normtype,"QFNorm"))
+ QuadForm(normtype) else diag(nrow(as.matrix(Cov)))
- return(list(trAsCov = as.vector(Cov)))
+ return(list(trAsCov = sum(diag(std%*%as.matrix(Cov)))))
})
+
setMethod("getAsRisk", signature(risk = "trAsCov",
L2deriv = "RealRandVariable",
neighbor = "ContNeighborhood",
biastype = "ANY"),
- function(risk, L2deriv, neighbor, biastype, Distr, clip, cent, stand, normtype){
+ function(risk, L2deriv, neighbor, biastype, normtype, clip, cent, stand, Distr, trafo = NULL,
+ V.comp = matrix(TRUE, ncol = nrow(stand), nrow = nrow(stand)), w,...){
+ if(missing(biastype)||is.null(biastype)) biastype <- biastype(risk)
+ if(missing(normtype)||is.null(normtype)) normtype <- normtype(risk)
Cov <- getAsRisk(risk = asCov(), L2deriv = L2deriv, neighbor = neighbor,
- biastype = biastype(risk), Distr = Distr, clip = clip,
- cent = cent, stand = stand)$asCov
+ biastype = biastype, Distr = Distr, clip = clip,
+ cent = cent, stand = stand, trafo = trafo,
+ V.comp = V.comp, w = w)$asCov
p <- nrow(stand)
std <- if(is(normtype,"QFNorm")) QuadForm(normtype) else diag(p)
@@ -210,13 +227,45 @@
})
###############################################################################
+## Anscombe criterion
+###############################################################################
+setMethod("getAsRisk", signature(risk = "asAnscombe",
+ L2deriv = "UnivariateDistribution",
+ neighbor = "UncondNeighborhood",
+ biastype = "ANY"),
+ function(risk, L2deriv, neighbor, biastype, normtype = NULL, clip, cent, stand, trafo = NULL, FI, ... ){
+ if(missing(biastype)||is.null(biastype)) biastype <- biastype(risk)
+ if(missing(normtype)||is.null(normtype)) normtype <- normtype(risk)
+ trAsCov.0 <- getAsRisk(risk = trAsCov(), L2deriv = L2deriv, neighbor = neighbor,
+ biastype = biastype, normtype = normtype,
+ clip = clip, cent = cent, stand = stand)$trAsCov
+ return(list(asAnscombe = FI/trAsCov.0))
+ })
+setMethod("getAsRisk", signature(risk = "asAnscombe",
+ L2deriv = "RealRandVariable",
+ neighbor = "ContNeighborhood",
+ biastype = "ANY"),
+ function(risk, L2deriv, neighbor, biastype, normtype, clip, cent, stand, Distr, trafo = NULL,
+ V.comp = matrix(TRUE, ncol = nrow(stand), nrow = nrow(stand)), FI,
+ w, ...){
+ if(missing(biastype)||is.null(biastype)) biastype <- biastype(risk)
+ if(missing(normtype)||is.null(normtype)) normtype <- normtype(risk)
+ trAsCov.0 <- getAsRisk(risk = trAsCov(), L2deriv = L2deriv, neighbor = neighbor,
+ biastype = biastype, normtype = normtype,
+ Distr = Distr, clip = clip,
+ cent = cent, stand = stand, V.comp = V.comp,
+ w = w)$trAsCov
+ return(list(asAnscombe = FI/trAsCov.0))
+ })
+
+###############################################################################
## asymptotic under-/overshoot risk
###############################################################################
setMethod("getAsRisk", signature(risk = "asUnOvShoot",
L2deriv = "UnivariateDistribution",
neighbor = "UncondNeighborhood",
biastype = "ANY"),
- function(risk, L2deriv, neighbor, biastype, clip, cent, stand, trafo){
+ function(risk, L2deriv, neighbor, biastype, normtype = NULL,clip, cent, stand, trafo, ...){
if(identical(all.equal(neighbor at radius, 0), TRUE))
return(list(asUnOvShoot = pnorm(-risk at width/sqrt(as.vector(stand)))))
@@ -236,7 +285,8 @@
L2deriv = "UnivariateDistribution",
neighbor = "ContNeighborhood",
biastype = "onesidedBias"),
- function(risk, L2deriv, neighbor, biastype, trafo){
+ function(risk, L2deriv, neighbor, biastype, normtype = NULL,
+ clip = NULL, cent = NULL, stand = NULL, trafo, ...){
D1 <- L2deriv
if(!is(D1, "DiscreteDistribution"))
@@ -265,7 +315,8 @@
L2deriv = "UnivariateDistribution",
neighbor = "ContNeighborhood",
biastype = "asymmetricBias"),
- function(risk, L2deriv, neighbor, biastype, trafo){
+ function(risk, L2deriv, neighbor, biastype, normtype = NULL,
+ clip = NULL, cent = NULL, stand = NULL, trafo, ...){
nu1 <- nu(biastype)[1]
nu2 <- nu(biastype)[2]
num <- nu2/(nu1+nu2)
@@ -284,8 +335,8 @@
L2deriv = "UnivariateDistribution",
neighbor = "Neighborhood",
biastype = "onesidedBias"),
- function(risk, L2deriv, neighbor, biastype,
- clip, cent, stand, trafo){
+ function(risk, L2deriv, neighbor, biastype, normtype = NULL,
+ clip, cent, stand, trafo, ...){
A <- as.vector(stand)*as.vector(trafo)
r <- neighbor at radius
b <- clip*A
@@ -300,3 +351,39 @@
semvar <- (v+r^2*b^2)*pnorm(sv)+ r*b*sqrt(v)*dnorm(sv)
return(list(asSemivar = semvar))
})
+
+###############################################################################
+## asymptotic L1 risk
+###############################################################################
+setMethod("getAsRisk", signature(risk = "asL1",
+ L2deriv = "UnivariateDistribution",
+ neighbor = "Neighborhood",
+ biastype = "ANY"),
+ function(risk, L2deriv, neighbor, biastype, normtype = NULL, clip = NULL, cent = NULL, stand, trafo, ...){
+ if(!is.finite(neighbor at radius))
+ L1 <- Inf
+ else{
+ s <- getInfV(L2deriv, neighbor, biastype, clip, cent, stand)
+ r <- neighbor at radius
+ L1 <- get.asGRisk.fct(risk)(r,s=s^.5,b=clip)
+ }
+ return(list(asL1 = L1))
+ })
+
+###############################################################################
+## asymptotic L4 risk
+###############################################################################
+setMethod("getAsRisk", signature(risk = "asL4",
+ L2deriv = "UnivariateDistribution",
+ neighbor = "Neighborhood",
+ biastype = "ANY"),
+ function(risk, L2deriv, neighbor, biastype, normtype = NULL, clip = NULL, cent = NULL, stand, trafo, ...){
+ if(!is.finite(neighbor at radius))
+ L4 <- Inf
+ else{
+ s <- getInfV(L2deriv, neighbor, biastype, clip, cent, stand)
+ r <- neighbor at radius
+ L4 <- get.asGRisk.fct(risk)(r,s=s^.5,b=clip)
+ }
+ return(list(asL4 = L4))
+ })
Modified: pkg/ROptEst/R/getInfClip.R
===================================================================
--- pkg/ROptEst/R/getInfClip.R 2010-12-02 13:59:11 UTC (rev 436)
+++ pkg/ROptEst/R/getInfClip.R 2010-12-02 14:24:42 UTC (rev 437)
@@ -1,6 +1,7 @@
###############################################################################
## optimal clipping bound for asymptotic MSE
###############################################################################
+
setMethod("getInfClip", signature(clip = "numeric",
L2deriv = "UnivariateDistribution",
risk = "asMSE",
@@ -11,6 +12,7 @@
getInfGamma(L2deriv = L2deriv, risk = risk,
neighbor = neighbor, biastype = biastype, cent = cent, clip = clip))
})
+
setMethod("getInfClip", signature(clip = "numeric",
L2deriv = "UnivariateDistribution",
risk = "asMSE",
@@ -27,6 +29,7 @@
neighbor = neighbor, biastype = biastype, cent = cent, clip = clip))
}
})
+
setMethod("getInfClip", signature(clip = "numeric",
L2deriv = "EuclRandVariable",
risk = "asMSE",
@@ -40,6 +43,92 @@
})
###############################################################################
+## optimal clipping bound for asymptotic L1risk
+###############################################################################
+
+setMethod("getInfClip", signature(clip = "numeric",
+ L2deriv = "UnivariateDistribution",
+ risk = "asL1",
+ neighbor = "ContNeighborhood"),
+ function(clip, L2deriv, risk, neighbor, biastype,
+ cent, symm, trafo){
+ s <- getInfV(L2deriv, neighbor, biastype, clip, cent, stand=1)
+ r <- neighbor at radius
+ w <- r * clip / s^.5
+ dp <- 2*dnorm(w)
+ pp <- 2*pnorm(w)-1
+ return(s^.5*r*pp/dp +
+ getInfGamma(L2deriv = L2deriv, risk = risk,
+ neighbor = neighbor, biastype = biastype, cent = cent, clip = clip))
+ })
+
+setMethod("getInfClip", signature(clip = "numeric",
+ L2deriv = "UnivariateDistribution",
+ risk = "asL1",
+ neighbor = "TotalVarNeighborhood"),
+ function(clip, L2deriv, risk, neighbor, biastype,
+ cent, symm, trafo){
+ s <- getInfV(L2deriv, neighbor, biastype, clip, cent, stand=1)
+ r <- neighbor at radius
+ w <- r * clip / s^.5
+ dp <- 2*dnorm(w)
+ pp <- 2*pnorm(w)-1
+ lhs <- s^.5*r*pp/dp
+ if(symm){
+ return(lhs +
+ getInfGamma(L2deriv = sign(as.vector(trafo))*L2deriv, risk = risk,
+ neighbor = neighbor, biastype = biastype, cent = -clip/2, clip = clip))
+ }else{
+ return(lhs +
+ getInfGamma(L2deriv = sign(as.vector(trafo))*L2deriv, risk = risk,
+ neighbor = neighbor, biastype = biastype, cent = cent, clip = clip))
+ }
+ })
+
+
+###############################################################################
+## optimal clipping bound for asymptotic L4 risk
+###############################################################################
+
+setMethod("getInfClip", signature(clip = "numeric",
+ L2deriv = "UnivariateDistribution",
+ risk = "asL4",
+ neighbor = "ContNeighborhood"),
+ function(clip, L2deriv, risk, neighbor, biastype,
+ cent, symm, trafo){
+ s <- getInfV(L2deriv, neighbor, biastype, clip, cent, stand=1)
+ r <- neighbor at radius
+ mse <- r^2 *clip^2 + s
+ mse4 <- (r^2 *clip^2/3 + s)/mse
+ return(r^2*clip*mse4 +
+ getInfGamma(L2deriv = L2deriv, risk = risk,
+ neighbor = neighbor, biastype = biastype, cent = cent, clip = clip))
+ })
+
+setMethod("getInfClip", signature(clip = "numeric",
+ L2deriv = "UnivariateDistribution",
+ risk = "asL4",
+ neighbor = "TotalVarNeighborhood"),
+ function(clip, L2deriv, risk, neighbor, biastype,
+ cent, symm, trafo){
+ s <- getInfV(L2deriv, neighbor, biastype, clip, cent, stand=1)
+ r <- neighbor at radius
+ mse <- r^2 *clip^2 + s
+ mse4 <- (r^2 *clip^2/3 + s)/mse
+ if(symm){
+ return(r^2*clip*mse4 +
+ getInfGamma(L2deriv = sign(as.vector(trafo))*L2deriv, risk = risk,
+ neighbor = neighbor, biastype = biastype, cent = -clip/2, clip = clip))
+ }else{
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/robast -r 437
More information about the Robast-commits
mailing list