[Robast-commits] r309 - in branches/robast-0.7/pkg/ROptEst: . R inst inst/scripts man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Jun 30 14:45:01 CEST 2009
Author: stamats
Date: 2009-06-30 14:45:00 +0200 (Tue, 30 Jun 2009)
New Revision: 309
Added:
branches/robast-0.7/pkg/ROptEst/inst/CITATION
branches/robast-0.7/pkg/ROptEst/inst/NEWS
Modified:
branches/robast-0.7/pkg/ROptEst/DESCRIPTION
branches/robast-0.7/pkg/ROptEst/R/roptest.R
branches/robast-0.7/pkg/ROptEst/inst/scripts/GammaModel.R
branches/robast-0.7/pkg/ROptEst/man/0ROptEst-package.Rd
branches/robast-0.7/pkg/ROptEst/man/roptest.Rd
Log:
merged branch and trunk
Modified: branches/robast-0.7/pkg/ROptEst/DESCRIPTION
===================================================================
--- branches/robast-0.7/pkg/ROptEst/DESCRIPTION 2009-06-30 04:32:35 UTC (rev 308)
+++ branches/robast-0.7/pkg/ROptEst/DESCRIPTION 2009-06-30 12:45:00 UTC (rev 309)
@@ -1,6 +1,6 @@
Package: ROptEst
Version: 0.7
-Date: 2008-11-27
+Date: 2009-04-22
Title: Optimally robust estimation
Description: Optimally robust estimation in general smoothly
parameterized models using S4 classes and methods.
Modified: branches/robast-0.7/pkg/ROptEst/R/roptest.R
===================================================================
--- branches/robast-0.7/pkg/ROptEst/R/roptest.R 2009-06-30 04:32:35 UTC (rev 308)
+++ branches/robast-0.7/pkg/ROptEst/R/roptest.R 2009-06-30 12:45:00 UTC (rev 309)
@@ -1,8 +1,8 @@
###############################################################################
## Optimally robust estimation
###############################################################################
-roptest <- function(x, L2Fam, eps, eps.lower, eps.upper, initial.est,
- neighbor = ContNeighborhood(), risk = asMSE(), steps = 1,
+roptest <- function(x, L2Fam, eps, eps.lower, eps.upper, fsCor = 1, initial.est,
+ neighbor = ContNeighborhood(), risk = asMSE(), steps = 1L,
distance = CvMDist, startPar = NULL, verbose = FALSE,
useLast = getRobAStBaseOption("kStepUseLast"), ...){
es.call <- match.call()
@@ -41,6 +41,11 @@
if((eps < 0) || (eps > 0.5))
stop("'eps' has to be in (0, 0.5]")
}
+ if(fsCor <= 0)
+ stop("'fsCor' has to be positive")
+ if(length(fsCor) != 1){
+ stop("'fsCor' has to be of length 1")
+ }
if(!is.integer(steps))
steps <- as.integer(steps)
if(steps < 1){
@@ -66,9 +71,13 @@
r.upper <- sqrtn*eps.upper
ICstart <- radiusMinimaxIC(L2Fam = L2FamStart, neighbor = neighbor, risk = risk,
loRad = r.lower, upRad = r.upper, verbose = verbose)
+ if(!isTRUE(all.equal(fsCor, 1, tol = 1e-3))){
+ neighbor at radius <- neighborRadius(ICstart)*fsCor
+ infMod <- InfRobModel(center = L2FamStart, neighbor = neighbor)
+ ICstart <- optIC(model = infMod, risk = risk, verbose = verbose)
+ }
}else{
- r <- sqrtn*eps
- neighbor at radius <- r
+ neighbor at radius <- sqrtn*eps*fsCor
infMod <- InfRobModel(center = L2FamStart, neighbor = neighbor)
ICstart <- optIC(model = infMod, risk = risk, verbose = verbose)
}
Copied: branches/robast-0.7/pkg/ROptEst/inst/CITATION (from rev 308, pkg/ROptEst/inst/CITATION)
===================================================================
--- branches/robast-0.7/pkg/ROptEst/inst/CITATION (rev 0)
+++ branches/robast-0.7/pkg/ROptEst/inst/CITATION 2009-06-30 12:45:00 UTC (rev 309)
@@ -0,0 +1,20 @@
+if(!exists("meta") || is.null(meta)) meta <- packageDescription("ROptEst")
+year <- sub("-.*", "", meta$Date)
+note <- sprintf("R package version %s", meta$Version)
+
+citHeader("To cite package ROptEst in publications use:")
+
+citEntry(entry="Manual",
+ title = "ROptEst: Optimally robust estimation",
+ author = personList(as.person("M. Kohl"),
+ as.person("P. Ruckdeschel")),
+ language = "English",
+ year = year,
+ note = note,
+ type = "R package",
+ url = "http://robast.r-forge.r-project.org/",
+ textVersion = paste("Kohl, M., and Ruckdeschel, P.",
+ sprintf("(%s).", year),
+ "ROptEst: Optimally robust estimation.",
+ paste(note, ".", sep = ""),
+ "URL http://robast.r-forge.r-project.org/"))
Copied: branches/robast-0.7/pkg/ROptEst/inst/NEWS (from rev 308, pkg/ROptEst/inst/NEWS)
===================================================================
--- branches/robast-0.7/pkg/ROptEst/inst/NEWS (rev 0)
+++ branches/robast-0.7/pkg/ROptEst/inst/NEWS 2009-06-30 12:45:00 UTC (rev 309)
@@ -0,0 +1,19 @@
+###############################################################################
+## News: to package ROptEst
+###############################################################################
+
+#######################################
+version 0.7
+#######################################
++ added finite-sample correction factor to function roptest
+
+#######################################
+version 0.6.2
+#######################################
++ added functions to compute and plot cniper contamination and cniper points,
+ respectively.
++ introduced option("newDevice") to control new opening of graphic devices
++ use of on.exit() to restore old settings for options() and par() at the end
+ of functions
++ introduction of NEWS-file
++ update of CITATION-file (based on code provided by A. Zeileis on R help)
\ No newline at end of file
Modified: branches/robast-0.7/pkg/ROptEst/inst/scripts/GammaModel.R
===================================================================
--- branches/robast-0.7/pkg/ROptEst/inst/scripts/GammaModel.R 2009-06-30 04:32:35 UTC (rev 308)
+++ branches/robast-0.7/pkg/ROptEst/inst/scripts/GammaModel.R 2009-06-30 12:45:00 UTC (rev 309)
@@ -4,14 +4,64 @@
require(ROptEst)
options("newDevice"=TRUE)
+
## generates Gamma Family with
-## scale = 2 and shape = 0.5
-G <- GammaFamily(scale = 2, shape = 0.5)
+## scale = 2 and shape = 0.1
+G <- GammaFamily(scale = 2, shape = 0.1)
G # show G
-plot(G) # plot of Gammad(scale = 2, shape = 0.5) and L_2 derivative
+plot(G) # plot of Gammad(scale = 2, shape = 0.1) and L_2 derivative
distrExOptions(ErelativeTolerance = 1e-8) # increase precision for E
checkL2deriv(G)
+## 30.06.09: new method for "E" with
+## signature(object = "Gammad", fun = "function", cond = "missing")
+## in package distrEx introduced which slightly reduces the problem
+## documented below.
+
+## more precisely:
+## numerical integration gives
+E(Gammad(scale = 2, shape = 0.1), function(x) (log(x/2)-digamma(0.1))^2)
+
+## whereas
+trigamma(0.1)
+
+## Problem is more or less caused by integration of log(x/2)^2
+## occurs for shape parameter small (<< 1)
+E(Gammad(scale = 2, shape = 0.1), function(x) log(x/2)^2)
+distrExOptions(ErelativeTolerance = 1e-10)
+E(Gammad(scale = 2, shape = 0.1), function(x) log(x/2)^2)
+distrExOptions(ErelativeTolerance = 1e-7)
+E(Gammad(scale = 2, shape = 0.1), function(x) log(x/2)^2)
+
+## result should be
+res1 <- 2*digamma(0.1)*E(Gammad(scale = 2, shape = 0.1), function(x) log(x/2))
+res2 <- digamma(0.1)^2
+trigamma(0.1) + res1 - res2
+
+fun <- function(x) log(x/2)^2*dgamma(x, scale = 2, shape = 0.1)
+integrate(fun, lower = 0, upper = Inf, rel.tol = 1e-6)
+integrate(fun, lower = 0, upper = Inf, rel.tol = 1e-3)
+integrate(fun, lower = 1e-9, upper = Inf, rel.tol = 1e-6)
+integrate(fun, lower = 1e-13, upper = Inf, rel.tol = 1e-6)
+## problem at zero!
+fun(0)
+curve(fun, from = 0, to = 1, n = 501)
+
+fun <- function(x) (log(x/2)-digamma(0.1))^2*dgamma(x, scale = 2, shape = 0.1)
+curve(fun, from = 0, to = 1, n = 501)
+
+GLIntegrate(fun, lower = 1e-9, upper = qgamma(1-1e-9, scale = 2, shape = 0.1), order = 500)
+
+
+## generates Gamma Family with
+## scale = 2 and shape = 1
+G <- GammaFamily(scale = 2, shape = 1)
+G # show G
+plot(G) # plot of Gammad(scale = 2, shape = 1) and L_2 derivative
+distrExOptions(ErelativeTolerance = 1e-8) # increase precision for E
+checkL2deriv(G)
+
+
## classical optimal IC
IC0 <- optIC(model = G, risk = asCov())
IC0 # show IC
Modified: branches/robast-0.7/pkg/ROptEst/man/0ROptEst-package.Rd
===================================================================
--- branches/robast-0.7/pkg/ROptEst/man/0ROptEst-package.Rd 2009-06-30 04:32:35 UTC (rev 308)
+++ branches/robast-0.7/pkg/ROptEst/man/0ROptEst-package.Rd 2009-06-30 12:45:00 UTC (rev 309)
@@ -13,9 +13,8 @@
\tabular{ll}{
Package: \tab ROptEst\cr
Version: \tab 0.7 \cr
-Date: \tab 2008-11-27 \cr
-Depends: \tab R(>= 2.7.0), methods, distr(>= 2.0), distrEx(>= 2.0),
-distrMod(>= 2.0), RandVar(>= 0.6.4), RobAStBase\cr
+Date: \tab 2009-04-22 \cr
+Depends: \tab R(>= 2.7.0), methods, distr(>= 2.0), distrEx(>= 2.0),distrMod(>= 2.0), RandVar(>= 0.6.4), RobAStBase\cr
LazyLoad: \tab yes\cr
License: \tab LGPL-3\cr
URL: \tab http://robast.r-forge.r-project.org/\cr
Modified: branches/robast-0.7/pkg/ROptEst/man/roptest.Rd
===================================================================
--- branches/robast-0.7/pkg/ROptEst/man/roptest.Rd 2009-06-30 04:32:35 UTC (rev 308)
+++ branches/robast-0.7/pkg/ROptEst/man/roptest.Rd 2009-06-30 12:45:00 UTC (rev 309)
@@ -6,8 +6,8 @@
parametric families via k-step construction.
}
\usage{
-roptest(x, L2Fam, eps, eps.lower, eps.upper, initial.est,
- neighbor = ContNeighborhood(), risk = asMSE(), steps = 1,
+roptest(x, L2Fam, eps, eps.lower, eps.upper, fsCor = 1, initial.est,
+ neighbor = ContNeighborhood(), risk = asMSE(), steps = 1L,
distance = CvMDist, startPar = NULL, verbose = FALSE,
useLast = getRobAStBaseOption("kStepUseLast"), ...)
}
@@ -20,6 +20,8 @@
lower bound for the amount of gross errors. See details below. }
\item{eps.upper}{ positive real (\code{eps.lower} <= \code{eps.upper} <= 0.5):
upper bound for the amount of gross errors. See details below. }
+ \item{fsCor}{ positive real: factor used to correct the neighborhood radius;
+ see details. }
\item{initial.est}{ initial estimate for unknown parameter. If missing
minimum distance estimator is computed. }
\item{neighbor}{ object of class \code{"UncondNeighborhood"} }
@@ -65,6 +67,14 @@
If \code{eps} is missing, the radius-minimax estimator in sense of
Rieder et al. (2001, 2008), respectively Section 2.2 of Kohl (2005) is returned.
+ Finite-sample and higher order results suggest that the asymptotically
+ optimal procedure is to liberal. Using \code{fsCor} the radius can be
+ modified - as a rule enlarged - to obtain a more conservative estimate.
+ In case of normal location and scale there is function
+ \code{\link[RobLox]{finiteSampleCorrection}} which returns a finite-sample
+ corrected (enlarged) radius based on the results of large Monte-Carlo
+ studies.
+
The default value of argument \code{useLast} is set by the
global option \code{kStepUseLast} which by default is set to
\code{FALSE}. In case of general models \code{useLast}
@@ -194,6 +204,17 @@
## plot of relative and absolute information; cf. Kohl (2005)
infoPlot(pIC(robest))
+## finite-sample correction
+if(require(RobLox)){
+ n <- length(chem)
+ r <- 0.05*sqrt(n)
+ r.fi <- finiteSampleCorrection(n = n, r = r)
+ fsCor <- r.fi/r
+ robest <- roptest(chem, NormLocationScaleFamily(), eps = 0.05,
+ fsCor = fsCor, steps = 3)
+ estimate(robest)
+}
+
## compute optimally robust estimator (unknown contamination)
## takes some time -> use package RobLox!
robest1 <- roptest(chem, NormLocationScaleFamily(), eps.lower = 0.05,
More information about the Robast-commits
mailing list