[Robast-commits] r50 - in pkg/RobLox: R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Feb 19 15:09:42 CET 2008
Author: stamats
Date: 2008-02-19 15:09:42 +0100 (Tue, 19 Feb 2008)
New Revision: 50
Added:
pkg/RobLox/R/sysdata.rda
Modified:
pkg/RobLox/R/roblox.R
pkg/RobLox/man/roblox.Rd
Log:
use saved values for Lagrange multipliers and approxfun to compute optimally robust ICs
Modified: pkg/RobLox/R/roblox.R
===================================================================
--- pkg/RobLox/R/roblox.R 2008-02-19 12:37:59 UTC (rev 49)
+++ pkg/RobLox/R/roblox.R 2008-02-19 14:09:42 UTC (rev 50)
@@ -1,65 +1,121 @@
###############################################################################
## computation of radius-minimax IC
+## using predefined functions included in "sysdata.rda"
###############################################################################
.getlsInterval <- function(r, rlo, rup, mean, sd, delta, A.loc.start,
a.sc.start, A.sc.start, bUp, itmax){
- Ab <- rlsOptIC.AL(r = r, mean = mean, sd = sd, A.loc.start = A.loc.start,
- a.sc.start = a.sc.start, A.sc.start = A.sc.start,
- bUp = bUp, delta = delta, itmax = itmax, computeIC = FALSE)
+ if(r > 80){
+ Ab <- rlsOptIC.AL(r = r, mean = mean, sd = sd, A.loc.start = A.loc.start,
+ a.sc.start = a.sc.start, A.sc.start = A.sc.start,
+ bUp = bUp, delta = delta, itmax = itmax, computeIC = FALSE)
+ A1 <- Ab$A[1,1]
+ A2 <- Ab$A[2,2]
+ b <- Ab$b
+ }else{
+ A1 <- .getA1.locsc(r)
+ A2 <- .getA2.locsc(r)
+ b <- .getb.locsc(r)
+ }
if(rlo == 0){
- efflo <- (sum(diag(Ab$A)) - Ab$b^2*r^2)/(1.5*sd^2)
+ efflo <- (A1 + A2 - b^2*r^2)/(1.5*sd^2)
}else{
- Ablo <- rlsOptIC.AL(r = rlo, mean = mean, sd = sd, A.loc.start = A.loc.start,
+ if(rlo > 80){
+ Ablo <- rlsOptIC.AL(r = rlo, mean = mean, sd = sd, A.loc.start = A.loc.start,
+ a.sc.start = a.sc.start, A.sc.start = A.sc.start,
+ bUp = bUp, delta = delta, itmax = itmax, computeIC = FALSE)
+ efflo <- (A1 + A2 - b^2*(r^2 - rlo^2))/sum(diag(Ablo$A))
+ }else{
+ A1lo <- .getA1.locsc(rlo)
+ A2lo <- .getA2.locsc(rlo)
+ efflo <- (A1 + A2 - b^2*(r^2 - rlo^2))/(A1lo + A2lo)
+ }
+ }
+
+ if(rup > 80){
+ Abup <- rlsOptIC.AL(r = rup, mean = mean, sd = sd, A.loc.start = A.loc.start,
a.sc.start = a.sc.start, A.sc.start = A.sc.start,
bUp = bUp, delta = delta, itmax = itmax, computeIC = FALSE)
- efflo <- (sum(diag(Ab$A)) - Ab$b^2*(r^2 - rlo^2))/sum(diag(Ablo$A))
+ effup <- (A1 + A2 - b^2*(r^2 - rup^2))/sum(diag(Abup$A))
+ }else{
+ A1up <- .getA1.locsc(rup)
+ A2up <- .getA2.locsc(rup)
+ effup <- (A1 + A2 - b^2*(r^2 - rlo^2))/(A1up + A2up)
}
- Abup <- rlsOptIC.AL(r = rup, mean = mean, sd = sd, A.loc.start = A.loc.start,
- a.sc.start = a.sc.start, A.sc.start = A.sc.start,
- bUp = bUp, delta = delta, itmax = itmax, computeIC = FALSE)
- effup <- (sum(diag(Ab$A)) - Ab$b^2*(r^2 - rup^2))/sum(diag(Abup$A))
-
return(effup-efflo)
}
.getlInterval <- function(r, rlo, rup, mean, sd, bUp){
- Ab <- rlOptIC(r = r, mean = mean, sd = sd, bUp = bUp, computeIC = FALSE)
+ if(r > 80){
+ Ab <- rlOptIC(r = r, mean = mean, sd = sd, bUp = bUp, computeIC = FALSE)
+ A <- Ab$A
+ b <- Ab$b
+ }else{
+ A <- .getA.loc(r)
+ b <- .getb.loc(r)
+ }
if(rlo == 0){
- efflo <- (Ab$A - Ab$b^2*r^2)/sd^2
+ efflo <- (A - b^2*r^2)/sd^2
}else{
- Ablo <- rlOptIC(r = rlo, mean = mean, sd = sd, bUp = bUp, computeIC = FALSE)
- efflo <- (Ab$A - Ab$b^2*(r^2 - rlo^2))/Ablo$A
+ if(rlo > 80){
+ Ablo <- rlOptIC(r = rlo, mean = mean, sd = sd, bUp = bUp, computeIC = FALSE)
+ efflo <- (Ab$A - Ab$b^2*(r^2 - rlo^2))/Ablo$A
+ }else{
+ Alo <- .getA.loc(rlo)
+ efflo <- (A - b^2*(r^2 - rlo^2))/Alo
+ }
}
- Abup <- rlOptIC(r = rup, mean = mean, sd = sd, bUp = bUp, computeIC = FALSE)
- effup <- (Ab$A - Ab$b^2*(r^2 - rup^2))/Abup$A
+ if(rup > 80){
+ Abup <- rlOptIC(r = rup, mean = mean, sd = sd, bUp = bUp, computeIC = FALSE)
+ effup <- (Ab$A - Ab$b^2*(r^2 - rup^2))/Abup$A
+ }else{
+ Aup <- .getA.loc(rup)
+ effup <- (A - b^2*(r^2 - rup^2))/Aup
+ }
return(effup-efflo)
}
.getsInterval <- function(r, rlo, rup, mean, sd, delta, bUp, itmax){
- Ab <- rsOptIC(r = r, mean = mean, sd = sd, bUp = bUp, delta = delta,
- itmax = itmax, computeIC = FALSE)
+ if(r > 80){
+ Ab <- rsOptIC(r = r, mean = mean, sd = sd, bUp = bUp, delta = delta,
+ itmax = itmax, computeIC = FALSE)
+ A <- Ab$A
+ b <- Ab$b
+ }else{
+ A <- .getA.sc(r)
+ b <- .getb.sc(r)
+ }
if(rlo == 0){
- efflo <- (Ab$A - Ab$b^2*r^2)/(0.5*sd^2)
+ efflo <- (A - b^2*r^2)/(0.5*sd^2)
}else{
- Ablo <- rsOptIC(r = rlo, mean = mean, sd = sd, bUp = bUp, delta = delta,
+ if(rlo > 80){
+ Ablo <- rsOptIC(r = rlo, mean = mean, sd = sd, bUp = bUp, delta = delta,
+ itmax = itmax, computeIC = FALSE)
+ efflo <- (A - b^2*(r^2 - rlo^2))/Ablo$A
+ }else{
+ Alo <- .getA.sc(rlo)
+ efflo <- (A - b^2*(r^2 - rlo^2))/Alo
+ }
+ }
+
+ if(rup > 80){
+ Abup <- rsOptIC(r = rup, mean = mean, sd = sd, bUp = bUp, delta = delta,
itmax = itmax, computeIC = FALSE)
- efflo <- (Ab$A - Ab$b^2*(r^2 - rlo^2))/Ablo$A
+ effup <- (A - b^2*(r^2 - rup^2))/Abup$A
+ }else{
+ Aup <- .getA.sc(rup)
+ effup <- (A - b^2*(r^2 - rup^2))/Aup
}
- Abup <- rsOptIC(r = rup, mean = mean, sd = sd, bUp = bUp, delta = delta,
- itmax = itmax, computeIC = FALSE)
- effup <- (Ab$A - Ab$b^2*(r^2 - rup^2))/Abup$A
-
return(effup-efflo)
}
###############################################################################
## optimally robust estimator for normal location and/or scale
###############################################################################
-roblox <- function(x, mean, sd, eps, eps.lower, eps.upper, initial.est = "ksMD",
+roblox <- function(x, mean, sd, eps, eps.lower, eps.upper, initial.est = "med",
tol = 1e-6, A.loc.start = 1, a.sc.start = 0, A.sc.start = 0.5,
bUp = 1000, itmax = 100, returnIC = FALSE){
@@ -107,13 +163,26 @@
}
if(!missing(eps)){
- IC1 <- rlsOptIC.AL(r = sqrt(length(x))*eps, mean = mean, sd = sd,
- A.loc.start = A.loc.start, a.sc.start = a.sc.start,
- A.sc.start = A.sc.start, bUp = bUp, delta = tol,
- itmax = itmax)
- Infos(IC1) <- matrix(c("roblox",
- "optimally robust IC for AL estimators and 'asMSE'"),
- ncol = 2, dimnames = list(NULL, c("method", "message")))
+ r <- sqrt(length(x))*eps
+ if(r > 80){
+ IC1 <- rlsOptIC.AL(r = r, mean = mean, sd = sd,
+ A.loc.start = A.loc.start, a.sc.start = a.sc.start,
+ A.sc.start = A.sc.start, bUp = bUp, delta = tol,
+ itmax = itmax)
+ Infos(IC1) <- matrix(c("roblox",
+ "optimally robust IC for AL estimators and 'asMSE'"),
+ ncol = 2, dimnames = list(NULL, c("method", "message")))
+ }else{
+ 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)
+ 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),
+ info = c("roblox", "optimally robust IC for AL estimators and 'asMSE'")))
+ }
robEst <- oneStepEstimator(x, IC1, c(mean, sd))
if(returnIC)
return(list(optIC = IC1, mean = robEst[1], sd = robEst[2]))
@@ -128,17 +197,34 @@
mean = mean, sd = sd, 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
- IC1 <- rlsOptIC.AL(r = r, mean = mean, sd = sd,
- A.loc.start = A.loc.start, a.sc.start = a.sc.start,
- A.sc.start = A.sc.start, bUp = bUp, delta = tol,
- itmax = itmax)
+ if(r > 80){
+ IC1 <- rlsOptIC.AL(r = r, mean = mean, sd = sd,
+ A.loc.start = A.loc.start, a.sc.start = a.sc.start,
+ A.sc.start = A.sc.start, bUp = bUp, delta = tol,
+ itmax = itmax)
+ }else{
+ 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)
+ 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)))
+ }
if(rlo == 0){
ineff <- (sum(diag(stand(IC1))) - clip(IC1)^2*r^2)/(1.5*sd^2)
}else{
- Ablo <- rlsOptIC.AL(r = rlo, mean = mean, sd = sd, A.loc.start = A.loc.start,
- a.sc.start = a.sc.start, A.sc.start = A.sc.start,
- bUp = bUp, delta = tol, itmax = itmax, computeIC = FALSE)
- ineff <- (sum(diag(stand(IC1))) - clip(IC1)^2*(r^2 - rlo^2))/sum(diag(Ablo$A))
+ if(rlo > 80){
+ Ablo <- rlsOptIC.AL(r = rlo, mean = mean, sd = sd, A.loc.start = A.loc.start,
+ a.sc.start = a.sc.start, A.sc.start = A.sc.start,
+ bUp = bUp, delta = tol, itmax = itmax, computeIC = FALSE)
+ ineff <- (sum(diag(stand(IC1))) - clip(IC1)^2*(r^2 - rlo^2))/sum(diag(Ablo$A))
+ }else{
+ A1lo <- .getA1.locsc(rlo)
+ A2lo <- .getA2.locsc(rlo)
+ ineff <- (sum(diag(stand(IC1))) - clip(IC1)^2*(r^2 - rlo^2))/(A1lo + A2lo)
+ }
}
Infos(IC1) <- matrix(c(rep("roblox", 3),
paste("radius-minimax IC for radius interval [",
@@ -164,11 +250,21 @@
if(initial.est == "med") mean <- median(x, na.rm = TRUE)
if(!missing(eps)){
- IC1 <- rlOptIC(r = sqrt(length(x))*eps, mean = mean, sd = sd,
- bUp = bUp)
- Infos(IC1) <- matrix(c("roblox",
- "optimally robust IC for AL estimators and 'asMSE'"),
- ncol = 2, dimnames = list(NULL, c("method", "message")))
+ r <- sqrt(length(x))*eps
+ if(r > 80){
+ IC1 <- rlOptIC(r = r, mean = mean, sd = sd, bUp = bUp)
+ Infos(IC1) <- matrix(c("roblox",
+ "optimally robust IC for AL estimators and 'asMSE'"),
+ ncol = 2, dimnames = list(NULL, c("method", "message")))
+ }else{
+ A <- .getA.loc(r)
+ b <- .getb.loc(r)
+ 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),
+ info = c("roblox", "optimally robust IC for AL estimators and 'asMSE'")))
+ }
robEst <- oneStepEstimator(x, IC1, mean)
if(returnIC)
return(list(optIC = IC1, mean = robEst, sd = sd))
@@ -181,12 +277,26 @@
r <- uniroot(.getlInterval, lower = rlo+1e-5, upper = rup,
tol = .Machine$double.eps^0.25, rlo = rlo, rup = rup,
mean = mean, sd = sd, bUp = bUp)$root
- IC1 <- rlOptIC(r = r, mean = mean, sd = sd, bUp = bUp)
+ if(r > 80){
+ IC1 <- rlOptIC(r = r, mean = mean, sd = sd, bUp = bUp)
+ }else{
+ A <- .getA.loc(r)
+ b <- .getb.loc(r)
+ 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)))
+ }
if(rlo == 0){
ineff <- (as.vector(stand(IC1)) - clip(IC1)^2*r^2)/sd^2
}else{
- Ablo <- rlOptIC(r = rlo, mean = mean, sd = sd, bUp = bUp, computeIC = FALSE)
- ineff <- (as.vector(stand(IC1)) - clip(IC1)^2*(r^2 - rlo^2))/Ablo$A
+ if(rlo > 80){
+ Ablo <- rlOptIC(r = rlo, mean = mean, sd = sd, bUp = bUp, computeIC = FALSE)
+ ineff <- (as.vector(stand(IC1)) - clip(IC1)^2*(r^2 - rlo^2))/Ablo$A
+ }else{
+ Alo <- .getA.loc(rlo)
+ ineff <- (as.vector(stand(IC1)) - clip(IC1)^2*(r^2 - rlo^2))/Alo
+ }
}
Infos(IC1) <- matrix(c(rep("roblox", 3),
paste("radius-minimax IC for radius interval [",
@@ -213,11 +323,23 @@
if(initial.est == "med") sd <- mad(x, na.rm = TRUE)
if(!missing(eps)){
- IC1 <- rsOptIC(r = sqrt(length(x))*eps, mean = mean, sd = sd,
- bUp = bUp, delta = tol, itmax = itmax)
- Infos(IC1) <- matrix(c("roblox",
- "optimally robust IC for AL estimators and 'asMSE'"),
- ncol = 2, dimnames = list(NULL, c("method", "message")))
+ r <- sqrt(length(x))*eps
+ if(r > 80){
+ IC1 <- rsOptIC(r = r, mean = mean, sd = sd,
+ bUp = bUp, delta = tol, itmax = itmax)
+ Infos(IC1) <- matrix(c("roblox",
+ "optimally robust IC for AL estimators and 'asMSE'"),
+ ncol = 2, dimnames = list(NULL, c("method", "message")))
+ }else{
+ A <- .getA.sc(r)
+ a <- .geta.sc(r)
+ b <- .getb.sc(r)
+ 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),
+ info = c("roblox", "optimally robust IC for AL estimators and 'asMSE'")))
+ }
robEst <- oneStepEstimator(x, IC1, sd)
if(returnIC)
return(list(optIC = IC1, mean = mean, sd = robEst))
@@ -231,14 +353,29 @@
tol = .Machine$double.eps^0.25, rlo = rlo, rup = rup,
mean = mean, sd = sd, delta = tol, bUp = bUp,
itmax = itmax)$root
- IC1 <- rsOptIC(r = r, mean = mean, sd = sd, bUp = bUp,
- delta = tol, itmax = itmax)
+ if(r > 80){
+ IC1 <- rsOptIC(r = r, mean = mean, sd = sd, bUp = bUp,
+ delta = tol, itmax = itmax)
+ }else{
+ A <- .getA.sc(r)
+ a <- .geta.sc(r)
+ b <- .getb.sc(r)
+ 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)))
+ }
if(rlo == 0){
ineff <- (as.vector(stand(IC1)) - clip(IC1)^2*r^2)/(0.5*sd^2)
}else{
- Ablo <- rsOptIC(r = rlo, mean = mean, sd = sd, bUp = bUp, delta = tol,
- itmax = itmax, computeIC = FALSE)
- ineff <- (as.vector(stand(IC1)) - clip(IC1)^2*(r^2 - rlo^2))/Ablo$A
+ if(rlo > 80){
+ Ablo <- rsOptIC(r = rlo, mean = mean, sd = sd, bUp = bUp, delta = tol,
+ itmax = itmax, computeIC = FALSE)
+ ineff <- (as.vector(stand(IC1)) - clip(IC1)^2*(r^2 - rlo^2))/Ablo$A
+ }else{
+ Alo <- .getA.sc(rlo)
+ ineff <- (as.vector(stand(IC1)) - clip(IC1)^2*(r^2 - rlo^2))/Alo
+ }
}
Infos(IC1) <- matrix(c(rep("roblox", 3),
paste("radius-minimax IC for radius interval [",
Added: pkg/RobLox/R/sysdata.rda
===================================================================
(Binary files differ)
Property changes on: pkg/RobLox/R/sysdata.rda
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Modified: pkg/RobLox/man/roblox.Rd
===================================================================
--- pkg/RobLox/man/roblox.Rd 2008-02-19 12:37:59 UTC (rev 49)
+++ pkg/RobLox/man/roblox.Rd 2008-02-19 14:09:42 UTC (rev 50)
@@ -9,7 +9,7 @@
respectively.
}
\usage{
-roblox(x, mean, sd, eps, eps.lower, eps.upper, initial.est = "ksMD",
+roblox(x, mean, sd, eps, eps.lower, eps.upper, initial.est = "med",
tol = 1e-06, A.loc.start = 1, a.sc.start = 0, A.sc.start = 0.5,
bUp = 1000, itmax = 100, returnIC = FALSE)
}
@@ -70,17 +70,20 @@
if 'returnIC' is 'TRUE' the list also contains
\item{optIC}{ object of class \code{"ContIC"}; optimally robust IC }
}
-\references{
+\references{
+ Kohl, M. (2005) \emph{Numerical Contributions to the Asymptotic Theory of Robustness}.
+ Bayreuth: Dissertation.
+
Rieder, H. (1994) \emph{Robust Asymptotic Statistics}. New York: Springer.
+ Rieder, H., Kohl, M. and Ruckdeschel, P. (2008) The Costs of not Knowing
+ the Radius. Statistical Methods and Applications \emph{17}(1) 13-40.
+
Rieder, H., Kohl, M. and Ruckdeschel, P. (2001) The Costs of not Knowing
the Radius. Submitted. Appeared as discussion paper Nr. 81.
SFB 373 (Quantification and Simulation of Economic Processes),
Humboldt University, Berlin; also available under
\url{www.uni-bayreuth.de/departments/math/org/mathe7/RIEDER/pubs/RR.pdf}
-
- Kohl, M. (2005) \emph{Numerical Contributions to the Asymptotic Theory of Robustness}.
- Bayreuth: Dissertation.
}
\author{Matthias Kohl \email{Matthias.Kohl at stamats.de}}
%\note{}
More information about the Robast-commits
mailing list