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