[Robast-commits] r427 - in branches/robast-0.8/pkg/ROptEst: inst/scripts man

[noreply at r-forge.r-project.org]

Author: ruckdeschel
Date: 2010-09-18 18:17:46 +0200 (Sat, 18 Sep 2010)
New Revision: 427

Added:
   branches/robast-0.8/pkg/ROptEst/inst/scripts/AnscombeOrNot.R
Modified:
   branches/robast-0.8/pkg/ROptEst/man/getReq.Rd
Log:
ROptEst: Integrated new example to getReq.Rd and new script AnscombeOrNot.R 

Added: branches/robast-0.8/pkg/ROptEst/inst/scripts/AnscombeOrNot.R
===================================================================
--- branches/robast-0.8/pkg/ROptEst/inst/scripts/AnscombeOrNot.R	                        (rev 0)
+++ branches/robast-0.8/pkg/ROptEst/inst/scripts/AnscombeOrNot.R	2010-09-18 16:17:46 UTC (rev 427)
@@ -0,0 +1,60 @@
+###############################################################
+#
+# check out how much bias protection is obtainable for
+#     95% efficiency in ideal model in particular models
+#
+###############################################################
+
+require(ROptEst)
+
+# some check function
+checkOut <- function(L2M, nbd, extraICs = NULL, biastype = symmetricBias(),
+                     normtype= NormType()){
+  ## L2M : L2 model
+  ## nbd : neighborhood
+  ## extraICs : list of further ICs)
+  ## biastype and normtype (...)
+  
+  force(normtype)
+  force(biastype)
+  RobM <- InfRobModel(center = L2M, neighbor = nbd)
+  ## rmx - IC
+  IC.rmx <- radiusMinimaxIC(L2Fam=L2M, neighbor=nbd,
+                            risk=asMSE(biastype=biastype, normtype=normtype),
+                            upRad=15,loRad=0.005)
+  ## Anscombe - IC
+  IC.ans <- optIC(model = RobM, risk = asAnscombe(biastype=biastype,
+                                                  normtype=normtype))
+  ## Lower case / Most Bias robust Estimator
+  IC.mbe <- optIC(model = RobM, risk = asBias(biastype=biastype,
+                                              normtype=normtype), tol = 1e-10)
+
+  todo <- list(rmx=IC.rmx,ans=IC.ans,mbe=IC.mbe)
+  if(!is.null(extraICs)){
+     namICs <- names(extraICs)
+     extraICs <- lapply(extraICs, function(x) makeIC(x,L2M))
+     names(extraICs) <- namICs
+     todo <- c(todo,extraICs)
+  }
+  ie <- 1/c(unlist(lapply(todo,getMaxIneff, neighbor=nbd)))
+}
+
+contnb <- ContNeighborhood(radius = 0.5)
+totvnb <- TotalVarNeighborhood(radius = 0.5)
+medianmad <- list(function(x)sign(x),function(x)sign(abs(x)-qnorm(.75)))
+### Normal location and scale --- takes ~2min:
+system.time({print(round(mineff.ls <- checkOut(L2M = NormLocationScaleFamily(), nbd = contnb,
+         extraICs = list(medmad=medianmad)),3))})
+### Normal location
+system.time(print(round(mineff.l <- checkOut(L2M = NormLocationFamily(), nbd = contnb),3)))
+### Normal scale convex contamination:
+system.time(print(round(mineff.sc <- checkOut(L2M = NormScaleFamily(), nbd = contnb),3)))
+### Normal scale total variation
+system.time(print(round(mineff.sv <- checkOut(L2M = NormScaleFamily(), nbd = totvnb),3)))
+### Poisson(lambda=1) convex contamination:
+system.time(print(round(mineff.pc <- checkOut(L2M = PoisFamily(lambda = 1), nbd = contnb),3)))
+### Poisson(lambda=1) scale convex contamination:
+system.time(print(round(mineff.pv <- checkOut(L2M = PoisFamily(lambda = 1), nbd = totvnb),3)))
+
+
+

Modified: branches/robast-0.8/pkg/ROptEst/man/getReq.Rd
===================================================================
--- branches/robast-0.8/pkg/ROptEst/man/getReq.Rd	2010-09-18 14:30:04 UTC (rev 426)
+++ branches/robast-0.8/pkg/ROptEst/man/getReq.Rd	2010-09-18 16:17:46 UTC (rev 427)
@@ -48,6 +48,16 @@
 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)
+
+### when to use MAD and when Qn 
+##  for Qn, see C. Croux, P. Rousseeuw (1993). Alternatives to the Median 
+##      Absolute Deviation, JASA 88(424):1273-1283
+L2M <- NormScaleFamily()
+IC.mad <- makeIC(function(x)sign(abs(x)-qnorm(.75)),L2M)
+d.qn <- (2^.5*qnorm(5/8))^-1
+IC.qn <- makeIC(function(x) d.qn*(1/4 - pnorm(x+1/d.qn) + pnorm(x-1/d.qn)), L2M)
+getReq(asMSE(), neighbor, IC.mad, IC.qn)
+# => MAD is better once r > 0.5144 (i.e. for more than 2 outliers for n = 30)
 }
 \concept{Hampel risk}
 \concept{risk}