[Robast-commits] r67 - in pkg/RobLox: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Feb 22 11:00:04 CET 2008
Author: stamats
Date: 2008-02-22 11:00:03 +0100 (Fri, 22 Feb 2008)
New Revision: 67
Modified:
pkg/RobLox/DESCRIPTION
pkg/RobLox/R/roblox.R
pkg/RobLox/man/roblox.Rd
Log:
removed redundant arguments, slightly modified return value
Modified: pkg/RobLox/DESCRIPTION
===================================================================
--- pkg/RobLox/DESCRIPTION 2008-02-21 10:46:14 UTC (rev 66)
+++ pkg/RobLox/DESCRIPTION 2008-02-22 10:00:03 UTC (rev 67)
@@ -1,6 +1,6 @@
Package: RobLox
Version: 0.6.0
-Date: 2008-02-19
+Date: 2008-02-22
Title: Optimally robust influence curves for location and scale
Description: functions for the determination of optimally
robust influence curves in case of normal
Modified: pkg/RobLox/R/roblox.R
===================================================================
--- pkg/RobLox/R/roblox.R 2008-02-21 10:46:14 UTC (rev 66)
+++ pkg/RobLox/R/roblox.R 2008-02-22 10:00:03 UTC (rev 67)
@@ -2,8 +2,7 @@
## computation of radius-minimax IC
## using predefined functions included in "sysdata.rda"
###############################################################################
-.getlsInterval <- function(r, rlo, rup, delta, A.loc.start, a.sc.start,
- A.sc.start, bUp, itmax){
+.getlsInterval <- function(r, rlo, rup){
if(r > 10){
b <- 1.618128043
const <- 1.263094656
@@ -37,7 +36,7 @@
return(effup-efflo)
}
-.getlInterval <- function(r, rlo, rup, bUp){
+.getlInterval <- function(r, rlo, rup){
if(r > 10){
b <- sqrt(pi/2)
A <- b^2*(1+r^2)
@@ -63,7 +62,7 @@
return(effup-efflo)
}
-.getsInterval <- function(r, rlo, rup, delta, bUp, itmax){
+.getsInterval <- function(r, rlo, rup){
if(r > 10){
b <- 1/(4*qnorm(0.75)*dnorm(qnorm(0.75)))
A <- b^2*(1+r^2)
@@ -90,6 +89,11 @@
return(effup-efflo)
}
+
+
+###############################################################################
+## computation of k-step construction
+###############################################################################
.onestep.loc <- function(x, initial.est, A, b, sd){
u <- A*(x-initial.est)/sd^2
IC <- mean(u*pmin(1, b/abs(u)), na.rm = TRUE)
@@ -147,13 +151,13 @@
return(est)
}
}
+
+
###############################################################################
## optimally robust estimator for normal location and/or scale
###############################################################################
-roblox <- function(x, mean, sd, eps, eps.lower, eps.upper, initial.est, k = 1,
- tol = 1e-6, A.loc.start = 1, a.sc.start = 0, A.sc.start = 0.5,
- bUp = 1000, itmax = 100, returnIC = FALSE){
-
+roblox <- function(x, mean, sd, eps, eps.lower, eps.upper, initial.est, k = 1,
+ returnIC = FALSE){
if(missing(x))
stop("'x' is missing with no default")
if(missing(eps) && missing(eps.lower) && missing(eps.upper)){
@@ -165,11 +169,15 @@
eps.upper <- 0.5
if(missing(eps.lower) && !missing(eps.upper))
eps.lower <- 0
+ if(length(eps.lower) != 1 || length(eps.upper) != 1)
+ stop("'eps.lower' and 'eps.upper' have to be of length 1")
if(!is.numeric(eps.lower) || !is.numeric(eps.upper) || eps.lower >= eps.upper)
stop("'eps.lower' < 'eps.upper' is not fulfilled")
if((eps.lower < 0) || (eps.upper > 0.5))
stop("'eps.lower' and 'eps.upper' have to be in [0, 0.5]")
}else{
+ if(length(eps) != 1)
+ stop("'eps' has to be of length 1")
if(eps == 0)
stop("'eps = 0'! => use functions 'mean' and 'sd' for estimation")
if((eps < 0) || (eps > 0.5))
@@ -184,9 +192,6 @@
k <- as.integer(k)
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(missing(initial.est)){
mean <- median(x, na.rm = TRUE)
sd <- mad(x, na.rm = TRUE)
@@ -231,9 +236,7 @@
r <- (rlo + rup)/2
}else{
r <- uniroot(.getlsInterval, lower = rlo+1e-8, upper = rup,
- tol = .Machine$double.eps^0.25, rlo = rlo, rup = rup,
- delta = tol, A.loc.start = A.loc.start, a.sc.start = a.sc.start,
- A.sc.start = A.sc.start, bUp = bUp, itmax = itmax)$root
+ tol = .Machine$double.eps^0.25, rlo = rlo, rup = rup)$root
}
if(r > 10){
b <- sd*1.618128043
@@ -275,7 +278,10 @@
ncol = 2, dimnames = list(NULL, c("method", "message")))
return(list(optIC = IC1, mean = robEst[1], sd = robEst[2]))
}else
- return(list(mean = robEst[1], sd = robEst[2]))
+ return(list(mean = robEst[1], sd = robEst[2],
+ "contamination interval" = round(c(eps.lower, eps.upper), 3),
+ "least favorable contamination" = round(r/sqrtn, 3),
+ "maximum MSE-inefficiency" = round(ineff, 3)))
}
}else{
if(missing(mean)){
@@ -314,7 +320,7 @@
r <- (rlo+rup)/2
}else{
r <- uniroot(.getlInterval, lower = rlo+1e-8, upper = rup,
- tol = .Machine$double.eps^0.25, rlo = rlo, rup = rup, bUp = bUp)$root
+ tol = .Machine$double.eps^0.25, rlo = rlo, rup = rup)$root
}
if(r > 10){
b <- sd*sqrt(pi/2)
@@ -348,7 +354,10 @@
ncol = 2, dimnames = list(NULL, c("method", "message")))
return(list(optIC = IC1, mean = robEst, sd = sd))
}else
- return(list(mean = robEst, sd = sd))
+ return(list(mean = robEst, sd = sd,
+ "contamination interval" = round(c(eps.lower, eps.upper), 3),
+ "least favorable contamination" = round(r/sqrtn, 3),
+ "maximum MSE-inefficiency" = round(ineff, 3)))
}
}
if(missing(sd)){
@@ -389,8 +398,7 @@
r <- (rlo+rup)/2
}else{
r <- uniroot(.getsInterval, lower = rlo+1e-8, upper = rup,
- tol = .Machine$double.eps^0.25, rlo = rlo, rup = rup,
- delta = tol, bUp = bUp, itmax = itmax)$root
+ tol = .Machine$double.eps^0.25, rlo = rlo, rup = rup)$root
}
if(r > 10){
b <- sd/(4*qnorm(0.75)*dnorm(qnorm(0.75)))
@@ -426,7 +434,10 @@
ncol = 2, dimnames = list(NULL, c("method", "message")))
return(list(optIC = IC1, mean = mean, sd = robEst))
}else
- return(list(mean = mean, sd = robEst))
+ return(list(mean = mean, sd = robEst,
+ "contamination interval" = round(c(eps.lower, eps.upper), 3),
+ "least favorable contamination" = round(r/sqrtn, 3),
+ "maximum MSE-inefficiency" = round(ineff, 3)))
}
}
}
Modified: pkg/RobLox/man/roblox.Rd
===================================================================
--- pkg/RobLox/man/roblox.Rd 2008-02-21 10:46:14 UTC (rev 66)
+++ pkg/RobLox/man/roblox.Rd 2008-02-22 10:00:03 UTC (rev 67)
@@ -9,9 +9,7 @@
respectively.
}
\usage{
-roblox(x, mean, sd, eps, eps.lower, eps.upper, initial.est, k = 1,
- tol = 1e-06, A.loc.start = 1, a.sc.start = 0, A.sc.start = 0.5,
- bUp = 1000, itmax = 100, returnIC = FALSE)
+roblox(x, mean, sd, eps, eps.lower, eps.upper, initial.est, k = 1, returnIC = FALSE)
}
\arguments{
\item{x}{ vector \code{x} of data values }
@@ -26,16 +24,6 @@
\item{initial.est}{ initial estimate for \code{mean} and/or \code{sd}. If missing
median and/or MAD are used. }
\item{k}{ positive integer. k-step is used to compute the optimally robust estimator.}
- \item{tol}{ the desired accuracy (convergence tolerance). }
- \item{A.loc.start}{ positive real: starting value for
- the standardizing constant of the location part. }
- \item{a.sc.start}{ real: starting value for centering
- constant of the scale part. }
- \item{A.sc.start}{ positive real: starting value for
- the standardizing constant of the scale part. }
- \item{bUp}{ positive real: the upper end point of the
- interval to be searched for the clipping bound b. }
- \item{itmax}{ the maximum number of iterations. }
\item{returnIC}{ logical: should IC be returned. See details below. }
}
\details{
@@ -70,8 +58,14 @@
list of location and scale estimates
\item{mean }{Description of 'comp1'}
\item{sd }{Description of 'comp2'}
+
if 'returnIC' is 'TRUE' the list also contains
\item{optIC}{ object of class \code{"ContIC"}; optimally robust IC }
+
+ if 'returnIC' is 'FALSE' an 'eps' is missing the list also contains
+ \item{contamination interval}{ interval for the amount of gross errors }
+ \item{least favorable contamination}{ amount of gross errors used for the computations }
+ \item{maximum MSE-inefficiency}{ maximum (asymptotic) MSE-inefficiency }
}
\references{
Kohl, M. (2005) \emph{Numerical Contributions to the Asymptotic Theory of Robustness}.
More information about the Robast-commits
mailing list