[Robast-commits] r52 - in pkg/RobLox: R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Feb 20 06:13:44 CET 2008
Author: stamats
Date: 2008-02-20 06:13:44 +0100 (Wed, 20 Feb 2008)
New Revision: 52
Modified:
pkg/RobLox/R/roblox.R
pkg/RobLox/man/roblox.Rd
Log:
RobLox checks and installs without any errors or warnings on R version 2.6.2 Patched (2008-02-10 r44423)
Modified: pkg/RobLox/R/roblox.R
===================================================================
--- pkg/RobLox/R/roblox.R 2008-02-19 14:14:10 UTC (rev 51)
+++ pkg/RobLox/R/roblox.R 2008-02-20 05:13:44 UTC (rev 52)
@@ -40,7 +40,7 @@
}else{
A1up <- .getA1.locsc(rup)
A2up <- .getA2.locsc(rup)
- effup <- (A1 + A2 - b^2*(r^2 - rlo^2))/(A1up + A2up)
+ effup <- (A1 + A2 - b^2*(r^2 - rup^2))/(A1up + A2up)
}
return(effup-efflo)
@@ -115,20 +115,20 @@
###############################################################################
## optimally robust estimator for normal location and/or scale
###############################################################################
-roblox <- function(x, mean, sd, eps, eps.lower, eps.upper, initial.est = "med",
+roblox <- function(x, mean, sd, eps, eps.lower, eps.upper, initial.est,
tol = 1e-6, A.loc.start = 1, a.sc.start = 0, A.sc.start = 0.5,
bUp = 1000, itmax = 100, returnIC = FALSE){
if(missing(x))
stop("'x' is missing with no default")
- if(missing(eps) & missing(eps.lower) & missing(eps.upper)){
+ if(missing(eps) && missing(eps.lower) && missing(eps.upper)){
eps.lower <- 0
eps.upper <- 0.5
}
if(missing(eps)){
- if(!missing(eps.lower) & missing(eps.upper))
+ if(!missing(eps.lower) && missing(eps.upper))
eps.upper <- 0.5
- if(missing(eps.lower) & !missing(eps.upper))
+ if(missing(eps.lower) && !missing(eps.upper))
eps.lower <- 0
if(!is.numeric(eps.lower) || !is.numeric(eps.upper) || eps.lower >= eps.upper)
stop("'eps.lower' < 'eps.upper' is not fulfilled")
@@ -140,26 +140,19 @@
if((eps < 0) || (eps > 0.5))
stop("'eps' has to be in (0, 0.5]")
}
- if((initial.est != "ksMD") && (initial.est != "med"))
- stop("invalid 'initial.est'")
- if(missing(mean) & missing(sd)){
+ if(missing(mean) && missing(sd)){
if(!is.numeric(A.loc.start) || !is.numeric(a.sc.start) || !is.numeric(A.sc.start))
stop("Starting values 'A.loc.start', 'a.sc.start' and 'A.sc.start' have to be numeric")
- if(initial.est == "ksMD"){
- KSdist <- function(param, x){
- if(param[2] <= 0) return(Inf)
- return(ks.test(x, "pnorm", mean = param[1], sd = param[2])$statistic)
- }
- res <- optim(c(0, 1), f = KSdist, method = "Nelder-Mead",
- control=list(reltol = tol), x = x)$par
- mean <- res[1]
- sd <- res[2]
- }
- if(initial.est == "med"){
+ if(missing(initial.est)){
mean <- median(x, na.rm = TRUE)
sd <- mad(x, na.rm = TRUE)
+ }else{
+ if(!is.numeric(initial.est) || length(initial.est) != 2)
+ stop("'initial.est' needs to be a numeric vector of length 2 or missing")
+ mean <- initial.est[1]
+ sd <- initial.est[2]
}
if(!missing(eps)){
@@ -176,7 +169,7 @@
A <- sd^2*diag(c(.getA1.locsc(r), .getA2.locsc(r)))
a <- sd*c(0, .geta.locsc(r))
b <- sd*.getb.locsc(r)
- mse <- sd^2*(a1 + a3)
+ mse <- sum(diag(A))
IC1 <- generateIC(neighbor = ContNeighborhood(radius = r),
L2Fam = NormLocationScaleFamily(mean = mean, sd = sd),
res = list(A = A, a = a, b = b, d = NULL,
@@ -192,7 +185,7 @@
sqrtn <- sqrt(length(x))
rlo <- sqrtn*eps.lower
rup <- sqrtn*eps.upper
- r <- uniroot(.getlsInterval, lower = rlo+1e-5, upper = rup,
+ r <- uniroot(.getlsInterval, lower = rlo+1e-8, upper = rup,
tol = .Machine$double.eps^0.25, rlo = rlo, rup = rup,
mean = mean, sd = sd, delta = tol, A.loc.start = A.loc.start,
a.sc.start = a.sc.start, A.sc.start = A.sc.start,
@@ -206,11 +199,12 @@
A <- sd^2*diag(c(.getA1.locsc(r), .getA2.locsc(r)))
a <- sd*c(0, .geta.locsc(r))
b <- sd*.getb.locsc(r)
- mse <- sd^2*(a1 + a3)
+ mse <- sum(diag(A))
IC1 <- generateIC(neighbor = ContNeighborhood(radius = r),
L2Fam = NormLocationScaleFamily(mean = mean, sd = sd),
res = list(A = A, a = a, b = b, d = NULL,
- risk = list(asMSE = mse, asBias = b, asCov = mse - r^2*b^2)))
+ risk = list(asMSE = mse, asBias = b, asCov = mse - r^2*b^2),
+ info = c("roblox", "optimally robust IC for AL estimators and 'asMSE'")))
}
if(rlo == 0){
ineff <- (sum(diag(stand(IC1))) - clip(IC1)^2*r^2)/(1.5*sd^2)
@@ -240,14 +234,13 @@
}
}else{
if(missing(mean)){
- if(initial.est == "ksMD"){
- KSdist.mean <- function(mean, x, sd){
- return(ks.test(x, "pnorm", mean = mean, sd = sd)$statistic)
- }
- mean <- optimize(f = KSdist.mean, interval = c(min(x), max(x)),
- tol = tol, x = x, sd = sd)$minimum
+ if(missing(initial.est)){
+ mean <- median(x, na.rm = TRUE)
+ }else{
+ if(!is.numeric(initial.est) || length(initial.est) != 1)
+ stop("'initial.est' needs to be a numeric vector of length 1 or missing")
+ mean <- initial.est
}
- if(initial.est == "med") mean <- median(x, na.rm = TRUE)
if(!missing(eps)){
r <- sqrt(length(x))*eps
@@ -274,7 +267,7 @@
sqrtn <- sqrt(length(x))
rlo <- sqrtn*eps.lower
rup <- sqrtn*eps.upper
- r <- uniroot(.getlInterval, lower = rlo+1e-5, upper = rup,
+ r <- uniroot(.getlInterval, lower = rlo+1e-8, upper = rup,
tol = .Machine$double.eps^0.25, rlo = rlo, rup = rup,
mean = mean, sd = sd, bUp = bUp)$root
if(r > 80){
@@ -285,7 +278,8 @@
IC1 <- generateIC(neighbor = ContNeighborhood(radius = r),
L2Fam = NormLocationFamily(mean = mean, sd = sd),
res = list(A = as.matrix(A), a = 0, b = b, d = NULL,
- risk = list(asMSE = A, asBias = b, asCov = A - r^2*b^2)))
+ risk = list(asMSE = A, asBias = b, asCov = A - r^2*b^2),
+ info = c("roblox", "optimally robust IC for AL estimators and 'asMSE'")))
}
if(rlo == 0){
ineff <- (as.vector(stand(IC1)) - clip(IC1)^2*r^2)/sd^2
@@ -312,15 +306,13 @@
}
}
if(missing(sd)){
- if(initial.est == "ksMD"){
- KSdist.sd <- function(sd, x, mean){
- return(ks.test(x, "pnorm", mean = mean, sd = sd)$statistic)
- }
- sd <- optimize(f = KSdist.sd,
- interval = c(.Machine$double.eps^0.5, max(x)-min(x)),
- tol = tol, x = x, mean = mean)$minimum
+ if(missing(initial.est)){
+ sd <- mad(x, na.rm = TRUE)
+ }else{
+ if(!is.numeric(initial.est) || length(initial.est) != 1)
+ stop("'initial.est' needs to be a numeric vector of length 1 or missing")
+ sd <- initial.est
}
- if(initial.est == "med") sd <- mad(x, na.rm = TRUE)
if(!missing(eps)){
r <- sqrt(length(x))*eps
@@ -349,7 +341,7 @@
sqrtn <- sqrt(length(x))
rlo <- sqrtn*eps.lower
rup <- sqrtn*eps.upper
- r <- uniroot(.getsInterval, lower = rlo+1e-5, upper = rup,
+ r <- uniroot(.getsInterval, lower = rlo+1e-8, upper = rup,
tol = .Machine$double.eps^0.25, rlo = rlo, rup = rup,
mean = mean, sd = sd, delta = tol, bUp = bUp,
itmax = itmax)$root
@@ -363,7 +355,8 @@
IC1 <- generateIC(neighbor = ContNeighborhood(radius = r),
L2Fam = NormScaleFamily(mean = mean, sd = sd),
res = list(A = as.matrix(A), a = a, b = b, d = NULL,
- risk = list(asMSE = A, asBias = b, asCov = A - r^2*b^2)))
+ risk = list(asMSE = A, asBias = b, asCov = A - r^2*b^2),
+ info = c("roblox", "optimally robust IC for AL estimators and 'asMSE'")))
}
if(rlo == 0){
ineff <- (as.vector(stand(IC1)) - clip(IC1)^2*r^2)/(0.5*sd^2)
Modified: pkg/RobLox/man/roblox.Rd
===================================================================
--- pkg/RobLox/man/roblox.Rd 2008-02-19 14:14:10 UTC (rev 51)
+++ pkg/RobLox/man/roblox.Rd 2008-02-20 05:13:44 UTC (rev 52)
@@ -9,7 +9,7 @@
respectively.
}
\usage{
-roblox(x, mean, sd, eps, eps.lower, eps.upper, initial.est = "med",
+roblox(x, mean, sd, eps, eps.lower, eps.upper, initial.est,
tol = 1e-06, A.loc.start = 1, a.sc.start = 0, A.sc.start = 0.5,
bUp = 1000, itmax = 100, returnIC = FALSE)
}
@@ -23,10 +23,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{initial.est}{ indicates which initial estimate should be used. Must be
- one of \code{"ksMD"} or \code{"med"}. In case \code{"ksMD"} the
- Kolmogorov(-Smirnov) minimum distance estimator and in case \code{"med"}
- median and/or mad are computed, respectively. }
+ \item{initial.est}{ initial estimate for \code{mean} and/or \code{sd}. If missing
+ median and/or MAD are used. }
\item{tol}{ the desired accuracy (convergence tolerance). }
\item{A.loc.start}{ positive real: starting value for
the standardizing constant of the location part. }
@@ -41,25 +39,29 @@
}
\details{
Computes the optimally robust estimator for location with scale specified,
- scale with location specified, or both if neither is specified.
-
+ scale with location specified, or both if neither is specified. The computation
+ uses a one-step construction with an appropriate initial estimate for location
+ or scale or location and scale, respectively. Valid candidates are e.g.
+ median and/or MAD (default) as well as Kolmogorov(-Smirnov) or von Mises minimum
+ distance estimators; cf. Rieder (1994) and Kohl (2005).
+
If the amount of gross errors (contamination) is known, it can be
specified by \code{eps}. The radius of the corresponding infinitesimal
contamination neighborhood is obtained by multiplying \code{eps}
by the square root of the sample size.
-
+
If the amount of gross errors (contamination) is unknown, try to find a
rough estimate for the amount of gross errors, such that it lies
between \code{eps.lower} and \code{eps.upper}.
-
+
In case \code{eps.lower} is specified and \code{eps.upper} is missing,
\code{eps.upper} is set to 0.5. In case \code{eps.upper} is specified and
\code{eps.lower} is missing, \code{eps.lower} is set to 0.
-
+
If neither \code{eps} nor \code{eps.lower} and/or \code{eps.upper} is
specified, \code{eps.lower} and \code{eps.upper} are set to 0 and 0.5,
respectively.
-
+
If \code{eps} is missing, the radius-minimax estimator in sense of
Rieder et al. (2001), respectively Section 2.2 of Kohl (2005) is returned.
}
@@ -122,15 +124,23 @@
# most robust
c(median(x), mad(x))
-# Kolmogorov(-Smirnov) minimum distance estimator (robust)
-MDEstimator(x, ParamFamily = NormLocationScaleFamily(), distance = KolmogorovDist)
-
# optimally robust (amount of gross errors known)
c(res1$mean, res1$sd)
# optimally robust (amount of gross errors unknown)
c(res2$mean, res2$sd)
+# Kolmogorov(-Smirnov) minimum distance estimator (robust)
+(ks.est <- MDEstimator(x, ParamFamily = NormLocationScaleFamily(), distance = KolmogorovDist))
+
+# optimally robust (amount of gross errors known)
+roblox(x, eps = 0.05, initial.est = ks.est$estimate)
+
+# von Mises minimum distance estimator (robust)
+(vM.est <- MDEstimator(x, ParamFamily = NormLocationScaleFamily(), distance = vonMisesDist))
+
+# optimally robust (amount of gross errors known)
+roblox(x, eps = 0.05, initial.est = vM.est$estimate)
}
\concept{normal location}
\concept{normal scale}
More information about the Robast-commits
mailing list