1
2MCV <- function(lambdas, formula, data, tau = 0.5, k = 10){
3    F <- Munge(formula, lambdas = lambdas)
4    f <- rqss(F, data, tau = tau)
5    n <- f$n
6    m <- length(f$qss)
7    y <- f$y[1:n]
8    folds = sample(rep(1:k, length = n))
9    U = NULL
10    for(i in 1:k){
11	s = which(folds != i)
12	M = rqss(F, data = data[s,], tau = tau)
13	nd = data[-s,]
14	G = matrix(0,nrow(nd),m)
15	for(j in 1:m){ #remove extrapolates, if any
16	    g = f$qss[[j]]$xyz[,1]
17	    G[,j] = (min(g[s]) < g[-s]) & (g[-s] < max(g[s]))
18	}
19	h = as.logical(apply(G,1,prod))
20	u = predict(M, newdata = nd[h,]) - (y[-s])[h]
21	U = c(U,(u * (tau - (u < 0))))
22    }
23    mean(U)
24}
25set.seed(1729)
26n <- 200
27x <- sort(runif(n, 0, 20))
28g0 <- function(x, tau)
29    log(x) + 0.2*(log(x))^3 + log(x) * qnorm(tau)/4
30y <- g0(x, runif(n))
31D <- data.frame(y = y, x = x)
32lams <- mcvs <- seq(.02, 5, by = 0.2)
33for(i in 1:length(mcvs))
34    mcvs[i] <- MCV(lams[i], y ~ qss(x, lambda = lambdas[1]), D)
35par(mfrow = c(1,2))
36plot(lams, mcvs, cex = .5, lwd = 2, type = 'l',
37    xlab = expression(lambda), ylab = expression(MCV( lambda )))
38lambdastar <- lams[which.min(mcvs)]
39
40plot(x, y, cex = .5, col = "grey")
41f <- rqss(y ~ qss(x, lambda = lambdastar), data = D)
42plot(f, add = TRUE, lwd = 2)
43lines(x,g0(x, 0.5), col = "red", lwd = 2)
44text(10, 1,bquote(lambda == ~  .(lambdastar)))
45
46