[Robast-commits] r495 - in branches/robast-0.9/pkg/RobLox: . R man tests/Examples
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Jul 1 15:46:36 CEST 2012
Author: stamats
Date: 2012-07-01 15:46:36 +0200 (Sun, 01 Jul 2012)
New Revision: 495
Added:
branches/robast-0.9/pkg/RobLox/R/showdown.R
branches/robast-0.9/pkg/RobLox/man/showdown.Rd
Modified:
branches/robast-0.9/pkg/RobLox/DESCRIPTION
branches/robast-0.9/pkg/RobLox/NAMESPACE
branches/robast-0.9/pkg/RobLox/man/0RobLox-package.Rd
branches/robast-0.9/pkg/RobLox/tests/Examples/RobLox-Ex.Rout.save
Log:
merge of trunk and robast 0-9
Modified: branches/robast-0.9/pkg/RobLox/DESCRIPTION
===================================================================
--- branches/robast-0.9/pkg/RobLox/DESCRIPTION 2012-07-01 04:49:40 UTC (rev 494)
+++ branches/robast-0.9/pkg/RobLox/DESCRIPTION 2012-07-01 13:46:36 UTC (rev 495)
@@ -1,18 +1,18 @@
Package: RobLox
Version: 0.9
-Date: 2010-12-03
+Date: 2012-07-01
Title: Optimally robust influence curves and estimators for location and scale
Description: Functions for the determination of optimally robust influence curves and
estimators in case of normal location and/or scale
-Depends: R(>= 2.7.0), stats, distrMod(>= 2.0.1), RobAStBase(>= 0.1.1)
-Suggests: Biobase
+Depends: R(>= 2.7.0), stats, lattice, RColorBrewer, Biobase, distr, distrMod, RobAStBase
+Suggests: MASS
Author: Matthias Kohl <Matthias.Kohl at stamats.de>
Maintainer: Matthias Kohl <Matthias.Kohl at stamats.de>
+LazyLoad: yes
ByteCompile: yes
-LazyLoad: yes
License: LGPL-3
Encoding: latin1
URL: http://robast.r-forge.r-project.org/
LastChangedDate: {$LastChangedDate$}
LastChangedRevision: {$LastChangedRevision$}
-SVNRevision: 439
+SVNRevision: 454
Modified: branches/robast-0.9/pkg/RobLox/NAMESPACE
===================================================================
--- branches/robast-0.9/pkg/RobLox/NAMESPACE 2012-07-01 04:49:40 UTC (rev 494)
+++ branches/robast-0.9/pkg/RobLox/NAMESPACE 2012-07-01 13:46:36 UTC (rev 495)
@@ -1,3 +1,5 @@
+importFrom(distr, q)
+
export(finiteSampleCorrection,
rlsOptIC.AL,
rlsOptIC.M,
@@ -21,4 +23,5 @@
rsOptIC,
roblox,
rowRoblox,
- colRoblox)
+ colRoblox,
+ showdown)
Copied: branches/robast-0.9/pkg/RobLox/R/showdown.R (from rev 494, pkg/RobLox/R/showdown.R)
===================================================================
--- branches/robast-0.9/pkg/RobLox/R/showdown.R (rev 0)
+++ branches/robast-0.9/pkg/RobLox/R/showdown.R 2012-07-01 13:46:36 UTC (rev 495)
@@ -0,0 +1,122 @@
+###############################################################################
+## Function to perform simulation study comparing some estimator with
+## rmx estimators
+###############################################################################
+showdown <- function(n, M, eps, contD, seed = 123, estfun, estMean, estSd,
+ eps.lower = 0, eps.upper = 0.05, steps = 3L, fsCor = TRUE,
+ plot1 = FALSE, plot2 = FALSE, plot3 = FALSE){
+ stopifnot(n >= 3)
+ stopifnot(eps >= 0, eps <= 0.5)
+ if(plot1){
+ from <- min(-6, q(contD)(1e-15))
+ to <- max(6, q(contD)(1-1e-15))
+ curve(pnorm, from = from, to = to, lwd = 2, n = 201,
+ main = "Comparison: ideal vs. real", ylab = "cdf")
+ fun <- function(x) (1-eps)*pnorm(x) + eps*p(contD)(x)
+ curve(fun, from = from, to = to, add = TRUE, col = "orange",
+ lwd = 2, n = 201, ylab = "cdf")
+ legend("topleft", legend = c("ideal", "real"),
+ fill = c("black", "orange"))
+ }
+
+ set.seed(seed)
+ r <- rbinom(n*M, prob = eps, size = 1)
+ Mid <- rnorm(n*M)
+ Mcont <- r(contD)(n*M)
+ Mre <- matrix((1-r)*Mid + r*Mcont, ncol = n)
+ ind <- rowSums(matrix(r, ncol = n)) >= n/2
+ while(any(ind)){
+ M1 <- sum(ind)
+ cat("Samples to re-simulate:\t", M1, "\n")
+ r <- rbinom(n*M1, prob = eps, size = 1)
+ Mid <- rnorm(n*M1)
+ Mcont <- r(contD)(n*M1)
+ Mre[ind,] <- (1-r)*Mid + r*Mcont
+ ind[ind] <- rowSums(matrix(r, ncol = n)) >= n/2
+ }
+ rm(Mid, Mcont, r, ind)
+
+
+ if(plot2){
+ ind <- if(M <= 20) 1:M else sample(1:M, 20)
+ if(plot1) dev.new()
+ M1 <- min(M, 20)
+ print(
+ stripplot(rep(1:M1, each = n) ~ as.vector(Mre[ind,]),
+ ylab = "samples", xlab = "x", pch = 20,
+ main = ifelse(M <= 20, "Samples", "20 randomly chosen samples"))
+ )
+ }
+
+ ## ML-estimator: mean and sd
+ Mean <- rowMeans(Mre)
+ Sd <- sqrt(rowMeans((Mre-Mean)^2))
+ ## Median and MAD
+ Median <- rowMedians(Mre)
+ Mad <- rowMedians(abs(Mre - Median))/qnorm(0.75)
+ ## Competitor
+ if(missing(estfun)){
+ Comp <- apply(Mre, 1, estMean)
+ Comp <- cbind(Comp, apply(Mre, 1, estSd))
+ }else
+ Comp <- t(apply(Mre, 1, estfun))
+
+ ## Radius-minimax estimator
+ RadMinmax <- estimate(rowRoblox(Mre, eps.lower = eps.lower,
+ eps.upper = eps.upper, k = steps,
+ fsCor = fsCor))
+
+ if(plot3){
+ Ergebnis1 <- list(Mean, Median, Comp[,1], RadMinmax[,1])
+ Ergebnis2 <- list(Sd, Mad, Comp[,2], RadMinmax[,2])
+ myCol <- brewer.pal(4, "Dark2")
+ if(plot1 || plot2) dev.new()
+ layout(matrix(c(1, 1, 1, 1, 3, 2, 2, 2, 2, 3), ncol = 2))
+ boxplot(Ergebnis1, col = myCol, pch = 20, main = "Location")
+ abline(h = 0)
+ boxplot(Ergebnis2, col = myCol, pch = 20, main = "Scale")
+ abline(h = 1)
+ op <- par(mar = rep(2, 4))
+ plot(c(0,1), c(1, 0), type = "n", axes = FALSE)
+ legend("center", c("ML", "Med/MAD", "competitor", "rmx"),
+ fill = myCol, ncol = 4, cex = 1.5)
+ on.exit(par(op))
+ }
+
+ ## ML-estimator
+ MSE1.1 <- n*mean(Mean^2)
+ ## Median + MAD
+ MSE2.1 <- n*mean(Median^2)
+ ## Tukey
+ MSE3.1 <- n*mean(Comp[,1]^2)
+ ## Radius-minimax
+ MSE4.1 <- n*mean(RadMinmax[,1]^2)
+ empMSE <- data.frame(ML = MSE1.1, Med = MSE2.1, Competitor = MSE3.1, "rmx" = MSE4.1)
+ rownames(empMSE) <- "n x empMSE (loc)"
+ relMSE <- empMSE[1,]/empMSE[1,4]
+ empMSE <- rbind(empMSE, relMSE)
+ rownames(empMSE)[2] <- "relMSE (loc)"
+
+ ## ML-estimator
+ MSE1.2 <- n*mean((Sd-1)^2)
+ ## Median + MAD
+ MSE2.2 <- n*mean((Mad-1)^2)
+ ## Tukey
+ MSE3.2 <- n*mean((Comp[,2]-1)^2)
+ ## Radius-minimax
+ MSE4.2 <- n*mean((RadMinmax[,2]-1)^2)
+ empMSE <- rbind(empMSE, c(MSE1.2, MSE2.2, MSE3.2, MSE4.2))
+ rownames(empMSE)[3] <- "n x empMSE (scale)"
+ relMSE <- empMSE[3,]/empMSE[3,4]
+ empMSE <- rbind(empMSE, relMSE)
+ rownames(empMSE)[4] <- "relMSE (scale)"
+ empMSE <- rbind(empMSE, c(MSE1.1 + MSE1.2, MSE2.1 + MSE2.2, MSE3.1 + MSE3.2,
+ MSE4.1 + MSE4.2))
+ rownames(empMSE)[5] <- "n x empMSE (loc + scale)"
+ relMSE <- empMSE[5,]/empMSE[5,4]
+ empMSE <- rbind(empMSE, relMSE)
+ rownames(empMSE)[6] <- "relMSE (loc + scale)"
+
+ empMSE
+}
+
Modified: branches/robast-0.9/pkg/RobLox/man/0RobLox-package.Rd
===================================================================
--- branches/robast-0.9/pkg/RobLox/man/0RobLox-package.Rd 2012-07-01 04:49:40 UTC (rev 494)
+++ branches/robast-0.9/pkg/RobLox/man/0RobLox-package.Rd 2012-07-01 13:46:36 UTC (rev 495)
@@ -13,7 +13,7 @@
\tabular{ll}{
Package: \tab RobLox \cr
Version: \tab 0.9 \cr
-Date: \tab 2010-12-03 \cr
+Date: \tab 2012-07-1 \cr
Depends: \tab R(>= 2.7.0), stats, distrMod(>= 2.0.1), RobAStBase(>= 0.1.1) \cr
LazyLoad: \tab yes \cr
License: \tab LGPL-3 \cr
Copied: branches/robast-0.9/pkg/RobLox/man/showdown.Rd (from rev 494, pkg/RobLox/man/showdown.Rd)
===================================================================
--- branches/robast-0.9/pkg/RobLox/man/showdown.Rd (rev 0)
+++ branches/robast-0.9/pkg/RobLox/man/showdown.Rd 2012-07-01 13:46:36 UTC (rev 495)
@@ -0,0 +1,78 @@
+\name{showdown}
+\Rdversion{1.1}
+\alias{showdown}
+\title{Estimator Showdown by Monte-Carlo Study.}
+\description{
+ The function \code{showdown} can be used to perform Monte-Carlo studies
+ comparing a competitor with rmx estimators in case of normal location and scale.
+ In addition, maximum likelihood (ML) estimators (mean and sd) and median and
+ MAD are computed. The comparison is based on the empirical MSE.
+}
+\usage{
+showdown(n, M, eps, contD, seed = 123, estfun, estMean, estSd,
+ eps.lower = 0, eps.upper = 0.05, steps = 3L, fsCor = TRUE,
+ plot1 = FALSE, plot2 = FALSE, plot3 = FALSE)
+}
+%- maybe also 'usage' for other objects documented here.
+\arguments{
+ \item{n}{integer; sample size, should be at least 3.}
+ \item{M}{integer; Monte-Carlo replications.}
+ \item{eps}{amount of contamination in [0, 0.5].}
+ \item{contD}{object of class \code{"UnivariateDistribution"}; contaminating distribution.}
+ \item{seed}{random seed.}
+ \item{estfun}{function to compute location and scale estimator; see details below.}
+ \item{estMean}{function to compute location estimator; see details below.}
+ \item{estSd}{function to compute scale estimator; see details below.}
+ \item{eps.lower}{used by rmx estimator.}
+ \item{eps.upper}{used by rmx estimator.}
+ \item{steps}{integer; steps used for estimator construction.}
+ \item{fsCor}{logical; use finite-sample correction.}
+ \item{plot1}{logical; plot cdf of ideal and real distribution.}
+ \item{plot2}{logical; plot 20 (or M if M < 20) randomly selected samples.}
+ \item{plot3}{logical; generate boxplots of the results.}
+}
+\details{
+Normal location and scale with mean = 0 and sd = 1 is used as ideal model (without
+restriction due to equivariance).
+
+Since there is no estimator which yields reliable results if 50 percent or more of the
+observations are contaminated, we use a modification where we re-simulate all samples
+including at least 50 percent contaminated data.
+
+If \code{estfun} is specified it has to compute and return a location and scale estimate
+(vector of length 2). One can also specify the location and scale estimator separately
+by using \code{estMean} and \code{estSd} where \code{estMean} computes and returns
+the location estimate and \code{estSd} the scale estimate.
+
+We use funtion \code{\link{rowRoblox}} for the computation of the rmx estimator.
+}
+\value{Data.frame including empirical MSE (standardized by sample size n) and
+relMSE with respect to the rmx estimator.
+}
+\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.
+ Extended version: \url{http://www.stamats.de/RRlong.pdf}
+}
+\author{Matthias Kohl \email{Matthias.Kohl at stamats.de}}
+%\note{}
+\seealso{\code{\link{rowRoblox}}}
+\examples{
+library(MASS)
+## compare with Huber's Proposal 2
+showdown(n = 20, M = 100, eps = 0.02, contD = Norm(mean = 3, sd = 3),
+ estfun = function(x){ unlist(hubers(x)) },
+ plot1 = TRUE, plot2 = TRUE, plot3 = TRUE)
+
+## compare with Huber M estimator with MAD scale
+showdown(n = 20, M = 100, eps = 0.02, contD = Norm(mean = 3, sd = 3),
+ estfun = function(x){ unlist(huber(x)) },
+ plot1 = TRUE, plot2 = TRUE, plot3 = TRUE)
+}
+\concept{Monte-Carlo study}
+\keyword{robust}
Modified: branches/robast-0.9/pkg/RobLox/tests/Examples/RobLox-Ex.Rout.save
===================================================================
--- branches/robast-0.9/pkg/RobLox/tests/Examples/RobLox-Ex.Rout.save 2012-07-01 04:49:40 UTC (rev 494)
+++ branches/robast-0.9/pkg/RobLox/tests/Examples/RobLox-Ex.Rout.save 2012-07-01 13:46:36 UTC (rev 495)
@@ -22,14 +22,37 @@
> source(file.path(R.home("share"), "R", "examples-header.R"))
> options(warn = 1)
> library('RobLox')
-Loading required package: distrMod
+Loading required package: lattice
+Loading required package: RColorBrewer
+Loading required package: Biobase
+Loading required package: BiocGenerics
+
+Attaching package: ‘BiocGenerics’
+
+The following object(s) are masked from ‘package:stats’:
+
+ xtabs
+
+The following object(s) are masked from ‘package:base’:
+
+ Filter, Find, Map, Position, Reduce, anyDuplicated, cbind,
+ colnames, duplicated, eval, get, intersect, lapply, mapply, mget,
+ order, paste, pmax, pmax.int, pmin, pmin.int, rbind, rep.int,
+ rownames, sapply, setdiff, table, tapply, union, unique
+
+Welcome to Bioconductor
+
+ Vignettes contain introductory material; view with
+ 'browseVignettes()'. To cite Bioconductor, see
+ 'citation("Biobase")', and for packages 'citation("pkgname")'.
+
+Loading required package: distr
Loading required package: startupmsg
:startupmsg> Utilities for start-up messages (version 0.8)
:startupmsg>
:startupmsg> For more information see ?"startupmsg",
:startupmsg> NEWS("startupmsg")
-Loading required package: distr
Loading required package: sfsmisc
Loading required package: SweaveListingUtils
:SweaveListingUtils> Utilities for Sweave together with
@@ -84,6 +107,7 @@
df, qqplot, sd
+Loading required package: distrMod
Loading required package: distrEx
:distrEx> Extensions of package distr (version 2.4)
:distrEx>
@@ -316,32 +340,6 @@
> ind <- rbinom(200, size=1, prob=0.05)
> X <- matrix(rnorm(200, mean=ind*3, sd=(1-ind) + ind*9), nrow = 2)
> rowRoblox(X)
-Loading required package: Biobase
-Loading required package: BiocGenerics
-
-Attaching package: ‘BiocGenerics’
-
-The following object(s) are masked from ‘package:RandVar’:
-
- Map
-
-The following object(s) are masked from ‘package:stats’:
-
- xtabs
-
-The following object(s) are masked from ‘package:base’:
-
- Filter, Find, Map, Position, Reduce, anyDuplicated, cbind,
- colnames, duplicated, eval, get, intersect, lapply, mapply, mget,
- order, paste, pmax, pmax.int, pmin, pmin.int, rbind, rep.int,
- rownames, sapply, setdiff, table, tapply, union, unique
-
-Welcome to Bioconductor
-
- Vignettes contain introductory material; view with
- 'browseVignettes()'. To cite Bioconductor, see
- 'citation("Biobase")', and for packages 'citation("pkgname")'.
-
Evaluations of Optimally robust estimate:
-----------------------------------------
An object of class “Estimate”
@@ -365,9 +363,6 @@
>
>
> cleanEx()
-
-detaching ‘package:Biobase’, ‘package:BiocGenerics’
-
> nameEx("finiteSampleCorrection")
> ### * finiteSampleCorrection
>
@@ -1783,32 +1778,6 @@
> ind <- rbinom(200, size=1, prob=0.05)
> X <- matrix(rnorm(200, mean=ind*3, sd=(1-ind) + ind*9), nrow = 2)
> rowRoblox(X)
-Loading required package: Biobase
-Loading required package: BiocGenerics
-
-Attaching package: ‘BiocGenerics’
-
-The following object(s) are masked from ‘package:RandVar’:
-
- Map
-
-The following object(s) are masked from ‘package:stats’:
-
- xtabs
-
-The following object(s) are masked from ‘package:base’:
-
- Filter, Find, Map, Position, Reduce, anyDuplicated, cbind,
- colnames, duplicated, eval, get, intersect, lapply, mapply, mget,
- order, paste, pmax, pmax.int, pmin, pmin.int, rbind, rep.int,
- rownames, sapply, setdiff, table, tapply, union, unique
-
-Welcome to Bioconductor
-
- Vignettes contain introductory material; view with
- 'browseVignettes()'. To cite Bioconductor, see
- 'citation("Biobase")', and for packages 'citation("pkgname")'.
-
Evaluations of Optimally robust estimate:
-----------------------------------------
An object of class “Estimate”
@@ -2138,9 +2107,6 @@
>
>
> cleanEx()
-
-detaching ‘package:Biobase’, ‘package:BiocGenerics’
-
> nameEx("rsOptIC")
> ### * rsOptIC
>
@@ -2184,13 +2150,56 @@
>
>
>
+> cleanEx()
+> nameEx("showdown")
+> ### * showdown
+>
+> flush(stderr()); flush(stdout())
+>
+> ### Name: showdown
+> ### Title: Estimator Showdown by Monte-Carlo Study.
+> ### Aliases: showdown
+> ### Keywords: robust
+>
+> ### ** Examples
+>
+> library(MASS)
+> ## compare with Huber's Proposal 2
+> showdown(n = 20, M = 100, eps = 0.02, contD = Norm(mean = 3, sd = 3),
++ estfun = function(x){ unlist(hubers(x)) },
++ plot1 = TRUE, plot2 = TRUE, plot3 = TRUE)
+dev.new(): using pdf(file="Rplots1.pdf")
+ ML Med Competitor rmx
+n x empMSE (loc) 1.151729 1.571687 1.0217407 1.0412422
+relMSE (loc) 1.106111 1.509434 0.9812710 1.0000000
+n x empMSE (scale) 1.086260 1.493500 0.7204599 0.6113725
+relMSE (scale) 1.776757 2.442864 1.1784303 1.0000000
+n x empMSE (loc + scale) 2.237989 3.065187 1.7422006 1.6526147
+relMSE (loc + scale) 1.354211 1.854750 1.0542086 1.0000000
+>
+> ## compare with Huber M estimator with MAD scale
+> showdown(n = 20, M = 100, eps = 0.02, contD = Norm(mean = 3, sd = 3),
++ estfun = function(x){ unlist(huber(x)) },
++ plot1 = TRUE, plot2 = TRUE, plot3 = TRUE)
+dev.new(): using pdf(file="Rplots2.pdf")
+dev.new(): using pdf(file="Rplots3.pdf")
+ ML Med Competitor rmx
+n x empMSE (loc) 1.151729 1.571687 1.091427 1.0412422
+relMSE (loc) 1.106111 1.509434 1.048197 1.0000000
+n x empMSE (scale) 1.086260 1.493500 1.493499 0.6113725
+relMSE (scale) 1.776757 2.442864 2.442862 1.0000000
+n x empMSE (loc + scale) 2.237989 3.065187 2.584925 1.6526147
+relMSE (loc + scale) 1.354211 1.854750 1.564143 1.0000000
+>
+>
+>
> ### * <FOOTER>
> ###
> cat("Time elapsed: ", proc.time() - get("ptime", pos = 'CheckExEnv'),"\n")
-Time elapsed: 56.499 0.168 56.895 0 0
+Time elapsed: 54.099 0.72 55.019 0 0
> grDevices::dev.off()
-null device
- 1
+pdf
+ 2
> ###
> ### Local variables: ***
> ### mode: outline-minor ***
More information about the Robast-commits
mailing list