[Robast-commits] r295 - pkg/RobLoxBioC/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Apr 22 13:08:04 CEST 2009
Author: stamats
Date: 2009-04-22 13:08:04 +0200 (Wed, 22 Apr 2009)
New Revision: 295
Added:
pkg/RobLoxBioC/R/0AllGeneric.R
pkg/RobLoxBioC/R/IlluminaSimStudyFunction.R
Removed:
pkg/RobLoxBioC/R/AllGeneric.R
Log:
added IlluminaSimStudy function, renamed AllGeneric.R to 0AllGeneric.R to be on first position ...
Copied: pkg/RobLoxBioC/R/0AllGeneric.R (from rev 294, pkg/RobLoxBioC/R/AllGeneric.R)
===================================================================
--- pkg/RobLoxBioC/R/0AllGeneric.R (rev 0)
+++ pkg/RobLoxBioC/R/0AllGeneric.R 2009-04-22 11:08:04 UTC (rev 295)
@@ -0,0 +1,14 @@
+############# preparations ################
+.onLoad <- function(lib, pkg) {
+ require("methods", character = TRUE, quietly = TRUE)
+}
+
+if(!isGeneric("robloxbioc")){
+ setGeneric("robloxbioc",
+ function(x, ...) standardGeneric("robloxbioc"))
+}
+
+if(!isGeneric("KolmogorovMinDist")){
+ setGeneric("KolmogorovMinDist",
+ function(x, D, ...) standardGeneric("KolmogorovMinDist"))
+}
Property changes on: pkg/RobLoxBioC/R/0AllGeneric.R
___________________________________________________________________
Name: svn:mergeinfo
+
Deleted: pkg/RobLoxBioC/R/AllGeneric.R
===================================================================
--- pkg/RobLoxBioC/R/AllGeneric.R 2009-04-22 09:23:58 UTC (rev 294)
+++ pkg/RobLoxBioC/R/AllGeneric.R 2009-04-22 11:08:04 UTC (rev 295)
@@ -1,14 +0,0 @@
-############# preparations ################
-.onLoad <- function(lib, pkg) {
- require("methods", character = TRUE, quietly = TRUE)
-}
-
-if(!isGeneric("robloxbioc")){
- setGeneric("robloxbioc",
- function(x, ...) standardGeneric("robloxbioc"))
-}
-
-if(!isGeneric("KolmogorovMinDist")){
- setGeneric("KolmogorovMinDist",
- function(x, D, ...) standardGeneric("KolmogorovMinDist"))
-}
Added: pkg/RobLoxBioC/R/IlluminaSimStudyFunction.R
===================================================================
--- pkg/RobLoxBioC/R/IlluminaSimStudyFunction.R (rev 0)
+++ pkg/RobLoxBioC/R/IlluminaSimStudyFunction.R 2009-04-22 11:08:04 UTC (rev 295)
@@ -0,0 +1,126 @@
+###############################################################################
+## Function to perform simulation study comparing Illumina's default method
+## with rmx estimators
+###############################################################################
+IlluminaSimStudy <- function(n, M, eps, seed = 123,
+ eps.lower = 0, eps.upper = 0.05,
+ steps = 3L, fsCor = TRUE, contD,
+ 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 <- sample(1:M, min(M, 20))
+ if(plot1) dev.new()
+ print(
+ stripplot(rep(1:20, each = 20) ~ as.vector(Mre[ind,]),
+ ylab = "samples", xlab = "x", pch = 20,
+ main = "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)
+ ## Illumina method
+ ind <- (Mre < (Median - 3*Mad)) | (Mre > (Median + 3*Mad))
+ x.ill <- Mre
+ x.ill[ind] <- NA
+ Illum.mean <- rowMeans(x.ill, na.rm = TRUE)
+ n.ill <- rowSums(!is.na(x.ill))
+ n.ill[n.ill < 1] <- NA
+ Illum.sd <- sqrt(rowSums((x.ill - Illum.mean)^2, na.rm = TRUE)/(n.ill-1))
+ Illum <- cbind(Illum.mean, Illum.sd)
+
+ ## 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, Illum[,1], RadMinmax[,1])
+ Ergebnis2 <- list(Sd, Mad, Illum[,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", "Illumina", "rmx"),
+ fill = myCol, ncol = 5, cex = 1.5)
+ on.exit(par(op))
+ }
+
+ ## ML-estimator
+ MSE1.1 <- n*mean(Mean^2)
+ ## Median + MAD
+ MSE2.1 <- n*mean(Median^2)
+ ## Illumina's default method
+ MSE3.1 <- n*mean(Illum[,1]^2)
+ ## Radius-minimax
+ MSE4.1 <- n*mean(RadMinmax[,1]^2)
+ empMSE <- data.frame(ML = MSE1.1, Med = MSE2.1, Illumina = 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)
+ ## Illumina's default method
+ MSE3.2 <- n*mean((Illum[,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
+}
+
+
More information about the Robast-commits
mailing list