1require(robustbase) 2source(system.file("test-tools-1.R", package="Matrix", mustWork=TRUE)) 3showProc.time() 4 5data(hbk); hbk.x <- data.matrix(hbk[, 1:3]) 6 7covComed(hbk.x) 8covComed(hbk.x, n.iter=4) 9showProc.time() 10 11data(radarImage) 12covComed(radarImage) 13covComed(radarImage[,3:5], n.iter = 5) 14showProc.time() 15 16data(bushfire) ; covComed(bushfire) 17data(heart); covComed(heart[, 1:2]) 18data(starsCYG); covComed(starsCYG) 19data(stackloss); covComed(stack.x) 20showProc.time() 21 22if(!robustbase:::doExtras()) quit() 23 24## if ( doExtras ) ----------------------------------------------------------------- 25## ============== 26 27 28 29i.rr <- c("raw.cov", "raw.center", "cov", "center") 30n <- 1024 ; p <- 7 31set.seed(47) 32showSys.time( 33 rX <- replicate(100, covComed(matrix(rnorm(n*p), n,p))[i.rr], 34 simplify=FALSE)) 35## Computing simulation-average (cov / center) <==> looking at Bias 36## _FIXME_ Really look at "MSE = Var + Bias^2" -- or something like 37## "simulation-average Squared Error or other Loss" 38C0 <- Reduce("+", lapply(rX, `[[`, "raw.cov")) / length(rX) 39C. <- Reduce("+", lapply(rX, `[[`, "cov")) / length(rX) 40round(1000 * C0) 41round(1000 * C.) 42assert.EQ(C0, diag(p), tol= 0.04, giveRE=TRUE) #-> 0.02805 43assert.EQ(C., diag(p), tol= 0.09, giveRE=TRUE) #-> 0.06475 44## Hmm.. raw.cov is better than cov ?? 45c00 <- Reduce("+", lapply(rX, `[[`, "raw.center")) / length(rX) 46c0 <- Reduce("+", lapply(rX, `[[`, "center")) / length(rX) 47stopifnot(print(sqrt(mean( (c00 - rep(0, p))^2 ))) < 0.005)# 0.004188 48stopifnot(print(sqrt(mean( (c0 - rep(0, p))^2 ))) < 0.005)# 0.003434 49 50n <- 4096 ; p <- 11 51set.seed(17) 52showSys.time( 53 r4 <- replicate(64, covComed(matrix(10+rnorm(n*p), n,p))[i.rr], 54 simplify=FALSE)) 55C0 <- Reduce("+", lapply(r4, `[[`, "raw.cov")) / length(r4) 56C. <- Reduce("+", lapply(r4, `[[`, "cov")) / length(r4) 57round(1000 * C0) 58round(1000 * C.) 59assert.EQ(C0, diag(p), tol = 0.025, giveRE=TRUE) # 0.0162 60assert.EQ(C., diag(p), tol = 0.06 , giveRE=TRUE) # 0.0486 61## Again... raw.cov better than cov ?? 62c00 <- Reduce("+", lapply(r4, `[[`, "raw.center")) / length(r4) 63c0 <- Reduce("+", lapply(r4, `[[`, "center")) / length(r4) 64assert.EQ(c00, rep(10, p), tol = 2e-4, giveRE=TRUE)# 7.97267e-05 = "raw" is better ? 65assert.EQ(c0 , rep(10, p), tol = 2e-4, giveRE=TRUE)# 0.0001036 66 67