[Robast-commits] r419 - in branches/robast-0.8/pkg/ROptEst: . R inst/scripts man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Sep 3 13:55:08 CEST 2010
Author: ruckdeschel
Date: 2010-09-03 13:55:08 +0200 (Fri, 03 Sep 2010)
New Revision: 419
Added:
branches/robast-0.8/pkg/ROptEst/R/getMaxIneff.R
branches/robast-0.8/pkg/ROptEst/man/getMaxIneff.Rd
Modified:
branches/robast-0.8/pkg/ROptEst/NAMESPACE
branches/robast-0.8/pkg/ROptEst/inst/scripts/NormalLocationModel.R
branches/robast-0.8/pkg/ROptEst/man/getReq.Rd
Log:
ROptEst: + new function getMaxIneff() to compute for any IC of class 'IC' the maximal inefficiency on r in [0,Inf)
+ added calls to getMaxIneff() and getReq() in script to Normal Location
Modified: branches/robast-0.8/pkg/ROptEst/NAMESPACE
===================================================================
--- branches/robast-0.8/pkg/ROptEst/NAMESPACE 2010-09-03 01:46:48 UTC (rev 418)
+++ branches/robast-0.8/pkg/ROptEst/NAMESPACE 2010-09-03 11:55:08 UTC (rev 419)
@@ -26,6 +26,9 @@
"getL1normL2deriv",
"getModifyIC",
"cniperCont", "cniperPoint", "cniperPointPlot")
-exportMethods("updateNorm", "scaleUpdateIC", "eff", "get.asGRisk.fct")
-export("getL2normL2deriv","asAnscombe", "asL1", "asL4", "getReq")
+exportMethods("updateNorm", "scaleUpdateIC", "eff",
+ "get.asGRisk.fct")
+export("getL2normL2deriv",
+ "asAnscombe", "asL1", "asL4",
+ "getReq", "getMaxIneff")
export("roptest","getLagrangeMultByOptim","getLagrangeMultByIter")
Added: branches/robast-0.8/pkg/ROptEst/R/getMaxIneff.R
===================================================================
--- branches/robast-0.8/pkg/ROptEst/R/getMaxIneff.R (rev 0)
+++ branches/robast-0.8/pkg/ROptEst/R/getMaxIneff.R 2010-09-03 11:55:08 UTC (rev 419)
@@ -0,0 +1,76 @@
+getMaxIneff <- function(IC, neighbor, biastype = symmetricBias(),
+ normtype = NormType(), z.start = NULL,
+ A.start = NULL, maxiter = 50,
+ tol = .Machine$double.eps^0.4,
+ warn = TRUE, verbose = NULL){
+ if(!is(IC,"IC"))
+ stop("Argument IC must be of class 'IC'.")
+
+ ow <- options("warn")
+ on.exit(options(ow))
+
+ sb <- .getSB(IC,neighbor)
+ si <- sb$s^2
+ bi <- sb$b^2
+
+ risk <- asBias(biastype=biastype, normtype=normtype)
+
+ L2Fam <- eval(IC at CallL2Fam)
+ trafo <- trafo(L2Fam at param)
+ symm <- L2Fam at L2derivDistrSymm[[1]]
+ Finfo <- L2Fam at FisherInfo
+ L2derivDim <- numberOfMaps(L2Fam at L2deriv)
+
+ FI0 <- trafo%*%solve(Finfo)%*%t(trafo)
+ std <- if(is(normtype,"QFNorm"))
+ QuadForm(normtype) else diag(nrow(trafo))
+ s0 <- sum(diag(std%*%FI0))
+
+ if(L2derivDim==1){
+ L2deriv <- L2Fam at L2derivDistr[[1]]
+ b0 <- getInfRobIC(L2deriv = L2deriv, risk = risk,
+ neighbor = neighbor, symm = symm,
+ trafo = trafo, maxiter = maxiter,
+ tol = tol, warn = warn, Finfo = Finfo,
+ verbose = verbose)$risk$asBias$value^2
+ }else{
+ if(is(L2Fam at distribution, "UnivariateDistribution")){
+ if((length(L2Fam at L2deriv) == 1) & is(L2Fam at L2deriv[[1]], "RealRandVariable")){
+ L2deriv <- L2deriv[[1]]
+ L2derivSymm <- L2Fam at L2derivSymm
+ L2derivDistrSymm <- L2Fam at L2derivDistrSymm
+ }else{
+ L2deriv <- diag(dimension(L2deriv)) %*% L2deriv
+ L2deriv <- RealRandVariable(Map = L2deriv at Map,
+ Domain = L2deriv at Domain)
+ nrvalues <- numberOfMaps(L2deriv)
+ if(numberOfMaps(L2deriv) != nrvalues){
+ L1 <- vector("list", nrvalues)
+ L2 <- vector("list", nrvalues)
+ for(i in 1:nrvalues){
+ L1[[i]] <- NonSymmetric()
+ L2[[i]] <- NoSymmetry()
+ }
+ L2derivSymm <- new("FunSymmList", L1)
+ L2derivDistrSymm <- new("DistrSymmList", L2)
+ }
+ }
+ if(!warn) options(warn = -1)
+
+ b0 <- getInfRobIC(L2deriv = L2deriv, neighbor = neighbor,
+ risk = risk, 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, warn = warn,
+ verbose = verbose)$risk$asBias$value^2
+ }else{
+ stop("not yet implemented")
+ }
+
+ }
+ return(max(si/s0,bi/b0))
+}
+
+
+
\ No newline at end of file
Modified: branches/robast-0.8/pkg/ROptEst/inst/scripts/NormalLocationModel.R
===================================================================
--- branches/robast-0.8/pkg/ROptEst/inst/scripts/NormalLocationModel.R 2010-09-03 01:46:48 UTC (rev 418)
+++ branches/robast-0.8/pkg/ROptEst/inst/scripts/NormalLocationModel.R 2010-09-03 11:55:08 UTC (rev 419)
@@ -114,6 +114,22 @@
Risks(N0.IC8)
plot(N0.IC8)
+getReq(asMSE(),neighbor,N0.ICA,N0.IC1,n=1)
+getReq(asMSE(),neighbor,N0.ICA,N0.IC1,n=30)
+getReq(asL1(),neighbor,N0.ICA,N0.IC1,n=30)
+getReq(asL4(),neighbor,N0.ICA,N0.IC1,n=30)
+getReq(asMSE(),neighbor,N0.ICA,N0.IC7,n=30)
+getReq(asL1(),neighbor,N0.ICA,N0.IC7,n=30)
+getReq(asL4(),neighbor,N0.ICA,N0.IC7,n=30)
+getReq(asMSE(),neighbor,N0.IC1,N0.IC7,n=30)
+
+
+getMaxIneff(N0.ICA,neighbor=ContNeighborhood())
+getMaxIneff(N0.IC1,neighbor=ContNeighborhood())
+getMaxIneff(N0.IC7,neighbor=ContNeighborhood())
+
+
+
## least favorable radius
## (may take quite some time!)
(N0.r.rho1 <- leastFavorableRadius(L2Fam=N0, neighbor=ContNeighborhood(),
Added: branches/robast-0.8/pkg/ROptEst/man/getMaxIneff.Rd
===================================================================
--- branches/robast-0.8/pkg/ROptEst/man/getMaxIneff.Rd (rev 0)
+++ branches/robast-0.8/pkg/ROptEst/man/getMaxIneff.Rd 2010-09-03 11:55:08 UTC (rev 419)
@@ -0,0 +1,65 @@
+\name{getMaxIneff}
+\alias{getMaxIneff}
+
+\title{getMaxIneff -- computation of the maximal inefficiency of an IC}
+\description{
+ computes the maximal inefficiency of an IC for the radius range [0,Inf).
+}
+\usage{getMaxIneff(IC, neighbor, biastype = symmetricBias(),
+ normtype = NormType(), z.start = NULL,
+ A.start = NULL, maxiter = 50,
+ tol = .Machine$double.eps^0.4,
+ warn = TRUE, verbose = NULL)}
+\arguments{
+ \item{IC}{some IC of class \code{IC}}
+ \item{neighbor}{ object of class \code{Neighborhood};
+ the neighborhood at which to compute the bias. }
+ \item{biastype}{a bias type of class \code{BiasType}}
+ \item{normtype}{ a norm type of class \code{NormType}}
+ \item{z.start}{ initial value for the centering constant. }
+ \item{A.start}{ initial value for the standardizing matrix. }
+ \item{maxiter}{ the maximum number of iterations. }
+ \item{tol}{ the desired accuracy (convergence tolerance).}
+ \item{warn}{ logical: print warnings. }
+ \item{verbose}{ logical: if \code{TRUE}, some messages are printed }
+}
+%\details{}
+\value{The maximal inefficiency, i.e.; a number in [1,Inf).}
+\references{
+ Hampel et al. (1986) \emph{Robust Statistics}.
+ The Approach Based on Influence Functions. New York: Wiley.
+
+ Rieder, H. (1994) \emph{Robust Asymptotic Statistics}. New York: Springer.
+
+ Rieder, H., Kohl, M. and Ruckdeschel, P. (2008) The Costs of not Knowing
+ the Radius. Statistical Methods and Applications \emph{17}(1) 13-40.
+
+ Rieder, H., Kohl, M. and Ruckdeschel, P. (2001) The Costs of not Knowing
+ the Radius. Submitted. Appeared as discussion paper Nr. 81.
+ SFB 373 (Quantification and Simulation of Economic Processes),
+ Humboldt University, Berlin; also available under
+ \url{www.uni-bayreuth.de/departments/math/org/mathe7/RIEDER/pubs/RR.pdf}
+
+}
+\author{Peter Ruckdeschel \email{peter.ruckdeschel at fraunhofer.itwm.de}}
+%\note{}
+\examples{
+N0 <- NormLocationFamily(mean=2, sd=3)
+## L_2 family + infinitesimal neighborhood
+neighbor <- ContNeighborhood(radius = 0.5)
+N0.Rob1 <- InfRobModel(center = N0, neighbor = neighbor)
+## OBRE solution (ARE 95%)
+N0.ICA <- optIC(model = N0.Rob1, risk = asAnscombe(.95))
+## OMSE solution radius 0.5
+N0.ICM <- optIC(model=N0.Rob1, risk=asMSE())
+## RMX solution
+N0.ICR <- radiusMinimaxIC(L2Fam=N0, neighbor=neighbor,risk=asMSE())
+
+getMaxIneff(N0.ICA,neighbor)
+getMaxIneff(N0.ICM,neighbor)
+getMaxIneff(N0.ICR,neighbor)
+
+}
+\concept{Inefficiency}
+\concept{risk}
+\keyword{robust}
Modified: branches/robast-0.8/pkg/ROptEst/man/getReq.Rd
===================================================================
--- branches/robast-0.8/pkg/ROptEst/man/getReq.Rd 2010-09-03 01:46:48 UTC (rev 418)
+++ branches/robast-0.8/pkg/ROptEst/man/getReq.Rd 2010-09-03 11:55:08 UTC (rev 419)
@@ -37,11 +37,17 @@
N0.ICA <- optIC(model = N0.Rob1, risk = asAnscombe(.95))
## MSE solution
N0.ICM <- optIC(model=N0.Rob1, risk=asMSE())
+## RMX solution
+N0.ICR <- radiusMinimaxIC(L2Fam=N0, neighbor=neighbor,risk=asMSE())
getReq(asMSE(),neighbor,N0.ICA,N0.ICM,n=1)
getReq(asMSE(),neighbor,N0.ICA,N0.ICM,n=30)
getReq(asL1(),neighbor,N0.ICA,N0.ICM,n=30)
getReq(asL4(),neighbor,N0.ICA,N0.ICM,n=30)
+getReq(asMSE(),neighbor,N0.ICA,N0.ICR,n=30)
+getReq(asL1(),neighbor,N0.ICA,N0.ICR,n=30)
+getReq(asL4(),neighbor,N0.ICA,N0.ICR,n=30)
+getReq(asMSE(),neighbor,N0.ICM,N0.ICR,n=30)
}
\concept{Hampel risk}
\concept{risk}
More information about the Robast-commits
mailing list