[Robast-commits] r954 - in branches/robast-1.1/pkg/RobLox: . R inst man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Jul 17 10:01:07 CEST 2018
Author: ruckdeschel
Date: 2018-07-17 10:01:06 +0200 (Tue, 17 Jul 2018)
New Revision: 954
Modified:
branches/robast-1.1/pkg/RobLox/DESCRIPTION
branches/robast-1.1/pkg/RobLox/R/showdown.R
branches/robast-1.1/pkg/RobLox/inst/NEWS
branches/robast-1.1/pkg/RobLox/man/0RobLox-package.Rd
Log:
[RobLox] in branch 1.1 + wherever possible also use q.l internally instead of q to
provide functionality in IRKernel
Modified: branches/robast-1.1/pkg/RobLox/DESCRIPTION
===================================================================
--- branches/robast-1.1/pkg/RobLox/DESCRIPTION 2018-07-17 07:55:52 UTC (rev 953)
+++ branches/robast-1.1/pkg/RobLox/DESCRIPTION 2018-07-17 08:01:06 UTC (rev 954)
@@ -1,18 +1,19 @@
Package: RobLox
-Version: 1.0
-Date: 2015-05-03
+Version: 1.1.0
+Date: 2018-07-17
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.14.0), stats, distrMod(>= 2.5.2), RobAStBase(>= 0.9)
Imports: lattice, RColorBrewer, Biobase, RandVar(>= 0.9.2), distr(>= 2.5.2)
Suggests: MASS
-Authors at R: c(person("Matthias", "Kohl", role=c("cre", "cph"), email="Matthias.Kohl at stamats.de"),
- person("Peter", "Ruckdeschel", role=c("aut", "cph")))
+Authors at R: c(person("Matthias", "Kohl", role=c("cre", "cph"),
+ email="Matthias.Kohl at stamats.de"), person("Peter", "Ruckdeschel", role=c("aut",
+ "cph")))
ByteCompile: yes
License: LGPL-3
Encoding: latin1
URL: http://robast.r-forge.r-project.org/
LastChangedDate: {$LastChangedDate$}
LastChangedRevision: {$LastChangedRevision$}
-SVNRevision: -Inf
+SVNRevision: 940
Modified: branches/robast-1.1/pkg/RobLox/R/showdown.R
===================================================================
--- branches/robast-1.1/pkg/RobLox/R/showdown.R 2018-07-17 07:55:52 UTC (rev 953)
+++ branches/robast-1.1/pkg/RobLox/R/showdown.R 2018-07-17 08:01:06 UTC (rev 954)
@@ -1,122 +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)
- rad <- rbinom(n*M, prob = eps, size = 1)
- Mid <- rnorm(n*M)
- Mcont <- r(contD)(n*M)
- Mre <- matrix((1-rad)*Mid + rad*Mcont, ncol = n)
- ind <- rowSums(matrix(rad, ncol = n)) >= n/2
- while(any(ind)){
- M1 <- sum(ind)
- cat("Samples to re-simulate:\t", M1, "\n")
- rad <- rbinom(n*M1, prob = eps, size = 1)
- Mid <- rnorm(n*M1)
- Mcont <- r(contD)(n*M1)
- Mre[ind,] <- (1-rad)*Mid + rad*Mcont
- ind[ind] <- rowSums(matrix(rad, ncol = n)) >= n/2
- }
- rm(Mid, Mcont, rad, 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
-}
-
+###############################################################################
+## 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.l(contD)(1e-15))
+ to <- max(6, q.l(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)
+ rad <- rbinom(n*M, prob = eps, size = 1)
+ Mid <- rnorm(n*M)
+ Mcont <- r(contD)(n*M)
+ Mre <- matrix((1-rad)*Mid + rad*Mcont, ncol = n)
+ ind <- rowSums(matrix(rad, ncol = n)) >= n/2
+ while(any(ind)){
+ M1 <- sum(ind)
+ cat("Samples to re-simulate:\t", M1, "\n")
+ rad <- rbinom(n*M1, prob = eps, size = 1)
+ Mid <- rnorm(n*M1)
+ Mcont <- r(contD)(n*M1)
+ Mre[ind,] <- (1-rad)*Mid + rad*Mcont
+ ind[ind] <- rowSums(matrix(rad, ncol = n)) >= n/2
+ }
+ rm(Mid, Mcont, rad, 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-1.1/pkg/RobLox/inst/NEWS
===================================================================
--- branches/robast-1.1/pkg/RobLox/inst/NEWS 2018-07-17 07:55:52 UTC (rev 953)
+++ branches/robast-1.1/pkg/RobLox/inst/NEWS 2018-07-17 08:01:06 UTC (rev 954)
@@ -8,6 +8,14 @@
information)
#######################################
+version 1.1
+#######################################
+
+under the hood:
++ wherever possible also use q.l internally instead of q to
+ provide functionality in IRKernel
+
+#######################################
version 1.0
#######################################
@@ -97,4 +105,4 @@
+ bug in computation of as. covariance in roblox corrected.
+ bug in rsOptIC corrected.
+ introduction of NEWS-file
-+ update of CITATION-file (based on code provided by A. Zeileis on R help)
\ No newline at end of file
++ update of CITATION-file (based on code provided by A. Zeileis on R help)
Modified: branches/robast-1.1/pkg/RobLox/man/0RobLox-package.Rd
===================================================================
--- branches/robast-1.1/pkg/RobLox/man/0RobLox-package.Rd 2018-07-17 07:55:52 UTC (rev 953)
+++ branches/robast-1.1/pkg/RobLox/man/0RobLox-package.Rd 2018-07-17 08:01:06 UTC (rev 954)
@@ -12,15 +12,15 @@
\details{
\tabular{ll}{
Package: \tab RobLox \cr
-Version: \tab 1.0 \cr
-Date: \tab 2015-05-03 \cr
+Version: \tab 1.1.0 \cr
+Date: \tab 2018-07-08 \cr
Depends: \tab R(>= 2.14.0), stats, distrMod(>= 2.5.2), RobAStBase(>= 0.9) \cr
Imports: \tab lattice, RColorBrewer, Biobase, RandVar(>= 0.9.2), distr(>= 2.5.2) \cr
Suggests: \tab MASS\cr
ByteCompile: \tab yes \cr
License: \tab LGPL-3 \cr
URL: \tab http://robast.r-forge.r-project.org/\cr
-SVNRevision: \tab -Inf \cr
+SVNRevision: \tab 940 \cr
}
}
\author{Matthias Kohl \email{matthias.kohl at stamats.de}}
More information about the Robast-commits
mailing list