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