[Robast-commits] r194 - in pkg/ROptEst: R man

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

Author: stamats
Date: 2008-10-31 11:06:52 +0100 (Fri, 31 Oct 2008)
New Revision: 194

Modified:
   pkg/ROptEst/R/getModifyIC.R
   pkg/ROptEst/R/getRiskIC.R
   pkg/ROptEst/man/roptest.Rd
Log:
incorporated the computation of confidence intervals, some minor modifications of code such that asympt. covariance gets computed.

Modified: pkg/ROptEst/R/getModifyIC.R
===================================================================
--- pkg/ROptEst/R/getModifyIC.R	2008-10-31 10:05:36 UTC (rev 193)
+++ pkg/ROptEst/R/getModifyIC.R	2008-10-31 10:06:52 UTC (rev 194)
@@ -123,8 +123,11 @@
                 b <- sdneu*clip(IC)/sdalt
                 a <- sdneu*cent(IC)/sdalt
                 mse <- sum(diag(A))
+                Cov <- sdneu^2*Risks(IC)$asCov/sdalt^2
+
                 res <- list(A = A, a = sdneu*cent(IC)/sdalt, b = b, d = NULL,
-                            risk = list(asMSE = mse, asBias = b, 
+                            risk = list(asCov = Cov,
+                                        asMSE = mse, asBias = b, 
                                         trAsCov = mse - r^2*b^2), 
                             info = Infos(IC), w = w,
                             normtype = normtype(IC), biastype = biastype(IC),

Modified: pkg/ROptEst/R/getRiskIC.R
===================================================================
--- pkg/ROptEst/R/getRiskIC.R	2008-10-31 10:05:36 UTC (rev 193)
+++ pkg/ROptEst/R/getRiskIC.R	2008-10-31 10:06:52 UTC (rev 194)
@@ -16,7 +16,10 @@
                                  L2Fam = "L2ParamFamily"),
     function(IC, risk, L2Fam){
         Cov <- IC at Risks[["asCov"]]
-        return(list(asCov = list(distribution = .getDistr(L2Fam), value = Cov)))
+        if(is.null(Cov))
+            return(getRiskIC(as(IC, "IC"), risk = risk, L2Fam = L2Fam))
+        else
+            return(list(asCov = list(distribution = .getDistr(L2Fam), value = Cov)))
     })
 
 setMethod("getRiskIC", signature(IC = "TotalVarIC", 

Modified: pkg/ROptEst/man/roptest.Rd
===================================================================
--- pkg/ROptEst/man/roptest.Rd	2008-10-31 10:05:36 UTC (rev 193)
+++ pkg/ROptEst/man/roptest.Rd	2008-10-31 10:06:52 UTC (rev 194)
@@ -111,15 +111,39 @@
 x <- rbinom(100, size=25, prob=(1-ind)*0.25 + ind*0.9)
 
 ## ML-estimate
-MLEstimator(x, BinomFamily(size = 25))
+MLest <- MLEstimator(x, BinomFamily(size = 25))
+estimate(MLest)
+confint(MLest)
 
 ## compute optimally robust estimator (known contamination)
-roptest(x, BinomFamily(size = 25), eps = 0.05, steps = 3)
+robest1 <- roptest(x, BinomFamily(size = 25), eps = 0.05, steps = 3)
+estimate(robest1)
+confint(robest1, method = symmetricBias())
+## neglecting bias
+confint(robest1)
+plot(pIC(robest1))
 
 ## compute optimally robust estimator (unknown contamination)
-roptest(x, BinomFamily(size = 25), eps.lower = 0, eps.upper = 0.2, steps = 3)
+robest2 <- roptest(x, BinomFamily(size = 25), eps.lower = 0, eps.upper = 0.2, steps = 3)
+estimate(robest2)
+confint(robest2, method = symmetricBias())
+plot(pIC(robest2))
 
+## total variation neighborhoods (known deviation)
+robest3 <- roptest(x, BinomFamily(size = 25), eps = 0.025, 
+                   neighbor = TotalVarNeighborhood(), steps = 3)
+estimate(robest3)
+confint(robest3, method = symmetricBias())
+plot(pIC(robest3))
 
+## total variation neighborhoods (unknown deviation)
+robest4 <- roptest(x, BinomFamily(size = 25), eps.lower = 0, eps.upper = 0.1, 
+                   neighbor = TotalVarNeighborhood(), steps = 3)
+estimate(robest4)
+confint(robest4, method = symmetricBias())
+plot(pIC(robest4))
+
+
 #############################
 ## 2. Poisson data
 #############################
@@ -129,28 +153,50 @@
        rep(10, 10), rep(11, 4), rep(12, 0), rep(13, 1), rep(14, 1))
 
 ## ML-estimate
-MLEstimator(x, PoisFamily())
+MLest <- MLEstimator(x, PoisFamily())
+estimate(MLest)
+confint(MLest)
 
 ## compute optimally robust estimator (unknown contamination)
-roptest(x, PoisFamily(), eps.upper = 0.1, steps = 3)
+robest <- roptest(x, PoisFamily(), eps.upper = 0.1, steps = 3)
+estimate(robest)
+confint(robest, symmetricBias())
+plot(pIC(robest))
 
+## total variation neighborhoods (unknown deviation)
+robest1 <- roptest(x, PoisFamily(), eps.upper = 0.05, 
+                  neighbor = TotalVarNeighborhood(), steps = 3)
+estimate(robest1)
+confint(robest1, symmetricBias())
+plot(pIC(robest1))
+
+
 #############################
 ## 3. Normal (Gaussian) location and scale
 #############################
-## Generate a contaminated sample
-ind <- rbinom(100, size=1, prob=0.05) 
-x <- rnorm(100, mean=0, sd=(1-ind) + ind*9)
+## 24 determinations of copper in wholemeal flour
+library(MASS)
+data(chem)
+plot(chem, main = "copper in wholemeal flour", pch = 20)
 
 ## ML-estimate
-MLEstimator(x, NormLocationScaleFamily())
+MLest <- MLEstimator(chem, NormLocationScaleFamily())
+estimate(MLest)
+confint(MLest)
 
 ## compute optimally robust estimator (known contamination)
-## takes some time
-roptest(x, NormLocationScaleFamily(), eps = 0.05, steps = 3)
+## takes some time -> you can use package RobLox for normal 
+## location and scale which is optimized for speed
+robest <- roptest(chem, NormLocationScaleFamily(), eps = 0.05, steps = 3)
+estimate(robest)
+confint(robest, symmetricBias())
 
 ## compute optimally robust estimator (unknown contamination)
-## takes some time
-roptest(x, NormLocationScaleFamily(), eps.upper = 0.2, steps = 3)
+## takes some time -> use package RobLox!
+robest1 <- roptest(chem, NormLocationScaleFamily(), eps.lower = 0.05, 
+                   eps.upper = 0.1, steps = 3)
+estimate(robest1)
+confint(robest1, symmetricBias())
 }
 \concept{k-step construction}
 \concept{optimally robust estimation}