[Robast-commits] r211 - in pkg/ROptEst: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Nov 15 13:55:22 CET 2008
Author: stamats
Date: 2008-11-15 13:55:22 +0100 (Sat, 15 Nov 2008)
New Revision: 211
Added:
pkg/ROptEst/R/cniperCont.R
pkg/ROptEst/man/cniperCont.Rd
Modified:
pkg/ROptEst/DESCRIPTION
pkg/ROptEst/NAMESPACE
pkg/ROptEst/R/AllGeneric.R
Log:
added functions to compute and plot cniper contamination and cniper points, respectively.
Modified: pkg/ROptEst/DESCRIPTION
===================================================================
--- pkg/ROptEst/DESCRIPTION 2008-11-12 21:33:45 UTC (rev 210)
+++ pkg/ROptEst/DESCRIPTION 2008-11-15 12:55:22 UTC (rev 211)
@@ -1,6 +1,6 @@
Package: ROptEst
Version: 0.6.1
-Date: 2008-10-31
+Date: 2008-11-13
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
Modified: pkg/ROptEst/NAMESPACE
===================================================================
--- pkg/ROptEst/NAMESPACE 2008-11-12 21:33:45 UTC (rev 210)
+++ pkg/ROptEst/NAMESPACE 2008-11-15 12:55:22 UTC (rev 211)
@@ -23,7 +23,8 @@
"lowerCaseRadius",
"minmaxBias", "getBiasIC",
"getL1normL2deriv",
- "getModifyIC")
+ "getModifyIC",
+ "cniperCont", "cniperPoint", "cniperPointPlot")
exportMethods("updateNorm")
export("getL2normL2deriv")
export("roptest")
Modified: pkg/ROptEst/R/AllGeneric.R
===================================================================
--- pkg/ROptEst/R/AllGeneric.R 2008-11-12 21:33:45 UTC (rev 210)
+++ pkg/ROptEst/R/AllGeneric.R 2008-11-15 12:55:22 UTC (rev 211)
@@ -74,3 +74,12 @@
if(!isGeneric("getModifyIC")){
setGeneric("getModifyIC", function(L2FamIC, neighbor, risk) standardGeneric("getModifyIC"))
}
+if(!isGeneric("cniperCont")){
+ setGeneric("cniperCont", function(IC1, IC2, L2Fam, neighbor, risk, ...) standardGeneric("cniperCont"))
+}
+if(!isGeneric("cniperPoint")){
+ setGeneric("cniperPoint", function(L2Fam, neighbor, risk, ...) standardGeneric("cniperPoint"))
+}
+if(!isGeneric("cniperPointPlot")){
+ setGeneric("cniperPointPlot", function(L2Fam, neighbor, risk, ...) standardGeneric("cniperPointPlot"))
+}
Added: pkg/ROptEst/R/cniperCont.R
===================================================================
--- pkg/ROptEst/R/cniperCont.R (rev 0)
+++ pkg/ROptEst/R/cniperCont.R 2008-11-15 12:55:22 UTC (rev 211)
@@ -0,0 +1,70 @@
+setMethod("cniperCont", signature(IC1 = "IC",
+ IC2 = "IC",
+ L2Fam = "L2ParamFamily",
+ neighbor = "ContNeighborhood",
+ risk = "asMSE"),
+ function(IC1, IC2, L2Fam, neighbor, risk, lower, upper, n = 101){
+ R1 <- Risks(IC1)[["trAsCov"]]
+ if(is.null(R1)) R1 <- getRiskIC(IC1, risk = trAsCov(), L2Fam = L2Fam)
+ if(length(R1) > 1) R1 <- R1$value
+
+ R2 <- Risks(IC2)[["trAsCov"]]
+ if(is.null(R2)) R2 <- getRiskIC(IC2, risk = trAsCov(), L2Fam = L2Fam)
+ if(length(R2) > 1) R2 <- R2$value
+
+ r <- neighbor at radius
+
+ fun <- function(x){
+ y1 <- evalIC(IC1, x)
+ y2 <- evalIC(IC2, x)
+ R1 - R2 + r^2*(as.vector(y1 %*% y1) - as.vector(y2 %*% y2))
+ }
+ x <- seq(from = lower, to = upper, length = n)
+ y <- sapply(x, fun)
+ plot(x, y, type = "l", main = "Cniper region plot",
+ xlab = "Dirac point", ylab = "Asymptotic MSE difference (IC1 - IC2)")
+# text(min(x), max(y)/2, "IC2", pos = 4)
+# text(min(x), min(y)/2, "IC1", pos = 4)
+ abline(h = 0)
+ invisible()
+ })
+setMethod("cniperPoint", signature(L2Fam = "L2ParamFamily",
+ neighbor = "ContNeighborhood",
+ risk = "asMSE"),
+ function(L2Fam, neighbor, risk, lower, upper){
+ tr.invF <- sum(diag(solve(FisherInfo(L2Fam))))
+ psi <- optIC(model = L2Fam, risk = asCov())
+ robMod <- InfRobModel(center = L2Fam, neighbor = neighbor)
+ eta <- optIC(model = robMod, risk = asMSE())
+ maxMSE <- Risks(eta)$asMSE
+ Delta <- sqrt(maxMSE - tr.invF)/neighbor at radius
+ fun <- function(x){
+ y <- evalIC(psi, x)
+ sqrt(as.vector(y %*% y)) - Delta
+ }
+ res <- uniroot(fun, lower = lower, upper = upper)$root
+ names(res) <- "cniper point"
+ res
+ })
+setMethod("cniperPointPlot", signature(L2Fam = "L2ParamFamily",
+ neighbor = "ContNeighborhood",
+ risk = "asMSE"),
+ function(L2Fam, neighbor, risk, lower, upper, n = 101){
+ tr.invF <- sum(diag(solve(FisherInfo(L2Fam))))
+ psi <- optIC(model = L2Fam, risk = asCov())
+ robMod <- InfRobModel(center = L2Fam, neighbor = neighbor)
+ eta <- optIC(model = robMod, risk = asMSE())
+ maxMSE <- Risks(eta)$asMSE
+ fun <- function(x){
+ 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)")
+# text(min(x), max(y)/2, "Robust", pos = 4)
+# text(min(x), min(y)/2, "Classic", pos = 4)
+ abline(h = 0)
+ invisible()
+ })
Added: pkg/ROptEst/man/cniperCont.Rd
===================================================================
--- pkg/ROptEst/man/cniperCont.Rd (rev 0)
+++ pkg/ROptEst/man/cniperCont.Rd 2008-11-15 12:55:22 UTC (rev 211)
@@ -0,0 +1,99 @@
+\name{cniperCont}
+\alias{cniperCont}
+\alias{cniperCont-methods}
+\alias{cniperCont,IC,IC,L2ParamFamily,ContNeighborhood,asMSE-method}
+\alias{cniperPoint}
+\alias{cniperPoint-methods}
+\alias{cniperPoint,L2ParamFamily,ContNeighborhood,asMSE-method}
+\alias{cniperPointPlot}
+\alias{cniperPointPlot-methods}
+\alias{cniperPointPlot,L2ParamFamily,ContNeighborhood,asMSE-method}
+\title{ Generic Functions for Computation and Plot of Cniper Contamination
+ and Cniper Points. }
+\description{
+ These generic functions and their methods can be used to determine cniper
+ contamination as well as cniper points. That is, under which (Dirac)
+ contamination is the risk of one procedure larger than the risk of some
+ other procedure.
+}
+\usage{
+cniperCont(IC1, IC2, L2Fam, neighbor, risk, ...)
+\S4method{cniperCont}{IC,IC,L2ParamFamily,ContNeighborhood,asMSE}(IC1,
+ IC2, L2Fam, neighbor, risk, lower, upper, n = 101)
+
+cniperPoint(L2Fam, neighbor, risk, ...)
+\S4method{cniperPoint}{L2ParamFamily,ContNeighborhood,asMSE}(L2Fam,
+ neighbor, risk, lower, upper)
+
+cniperPointPlot(L2Fam, neighbor, risk, ...)
+\S4method{cniperPointPlot}{L2ParamFamily,ContNeighborhood,asMSE}(L2Fam,
+ neighbor, risk, lower, upper, n = 101)
+}
+\arguments{
+ \item{IC1}{ object of class \code{IC} }
+ \item{IC2}{ object of class \code{IC} }
+ \item{L2Fam}{ object of class \code{L2ParamFamily} }
+ \item{neighbor}{ object of class \code{Neighborhood} }
+ \item{risk}{ object of class \code{RiskType} }
+ \item{\dots}{ additional parameters. }
+ \item{lower, upper}{ the lower and upper end points of the
+ contamination interval. }
+ \item{n}{ number of points between \code{lower} and \code{upper}}
+}
+\details{
+ In case of \code{cniperCont} the difference between the risks of two ICs
+ is plotted.
+
+ The function \code{cniperPoint} can be used to determine cniper
+ points. That is, points such that the optimally robust estimator
+ has smaller minimax risk than the classical optimal estimator under
+ contamination with Dirac measures at the cniper points.
+
+ As such points might be difficult to find, we provide the
+ function \code{cniperPointPlot} which can be used to obtain a plot
+ of the risk difference.
+
+ For more details about cniper contamination and cniper points we refer
+ to Section~3.5 of Kohl et al. (2008) as well as Ruckdeschel (2004) and
+ the Introduction of Kohl (2005).
+}
+\value{invisible() resp. cniper point is returned.}
+\references{
+ Kohl, M. and Ruckdeschel, H. and Rieder, H. (2008). Infinitesimally
+ Robust Estimation in General Smoothly Parametrized Models. Unpublished
+ Manuscript.
+
+ Kohl, M. (2005) \emph{Numerical Contributions to the Asymptotic Theory of Robustness}.
+ Bayreuth: Dissertation.
+
+ Ruckdeschel, P. (2004). Higher Order Asymptotics for the MSE of M-Estimators
+ on Shrinking Neighborhoods. Unpublished Manuscript.
+}
+\author{Matthias Kohl \email{Matthias.Kohl at stamats.de}}
+%\note{}
+%\seealso{ ~~objects to See Also as \code{\link{help}}, ~~~ }
+\examples{
+## cniper contamination
+P <- PoisFamily(lambda = 4)
+RobP1 <- InfRobModel(center = P, neighbor = ContNeighborhood(radius = 0.1))
+IC1 <- optIC(model=RobP1, risk=asMSE())
+RobP2 <- InfRobModel(center = P, neighbor = ContNeighborhood(radius = 1))
+IC2 <- optIC(model=RobP2, risk=asMSE())
+cniperCont(IC1 = IC1, IC2 = IC2, L2Fam = P,
+ neighbor = ContNeighborhood(radius = 0.5),
+ risk = asMSE(),
+ lower = 0, upper = 8, n = 101)
+
+## cniper point plot
+cniperPointPlot(P, neighbor = ContNeighborhood(radius = 0.5),
+ risk = asMSE(), lower = 0, upper = 10)
+
+## cniper point
+cniperPoint(P, neighbor = ContNeighborhood(radius = 0.5),
+ risk = asMSE(), lower = 0, upper = 4)
+cniperPoint(P, neighbor = ContNeighborhood(radius = 0.5),
+ risk = asMSE(), lower = 4, upper = 8)
+}
+\concept{cniper contamination}
+\concept{cniper point}
+\keyword{robust}
More information about the Robast-commits
mailing list