# generalized approximate cross-validation criterion (Yuan 2006) (GACV): cv.qenv <- function (X, Y, dims=NULL,tau=0.5, fold = 5, maxiter=1e2,ftol=1e-4,seed=NULL,method=2,parallel=FALSE,core=5) { time1=proc.time() X <- as.matrix(X) Y <- as.matrix(Y) n <- nrow(X) p <- ncol(X) r <- ncol(Y) if(is.null(dims)) dims=1:p set.seed(seed) if(!is.null(seed)){ set.seed(seed) tmp=sample(1:n) } else{ tmp=1:n } foldi <- split(tmp, rep(1:fold, length = n)) if(parallel){ #------ library(doParallel) registerDoParallel(cores=core) tmp<-foreach(i = 1:fold,.combine = rbind,.verbose=TRUE,.export=c('sechol','qenv_gmm','qr_gmm','n','foldi','tau')) %dopar% { omit <- foldi[[i]] tempX <- X[-omit,,drop=FALSE] tempY <- Y[-omit,,drop=FALSE] testX <- X[omit,,drop=FALSE] testY <- Y[omit,,drop=FALSE] out = rep(0,2*length(dims)) for(j in 1:length(dims)){ cat('Fold =',i, ', u=',dims[j] ,'.\n') m0 <- qenv_gmm(tempX,tempY,dims[j],tau=tau,method=method) pred <- matrix(1, nrow(testX),1) %*% t(m0$alpha) + testX %*% m0$beta resi <- testY - pred out[j] <-mean(resi^2) out[j+length(dims)]= mean(resi*(tau-(resi<0))) } return(out) } # end fold #------ mspemati=tmp[,1:length(dims)] RCVi=tmp[,-(1:length(dims))] } else { mspemati = matrix(0,fold,length(dims),dimnames=list(fold=1:fold,u=dims)) RCVi = matrix(0,fold,length(dims),dimnames=list(fold=1:fold,u=dims)) for (i in 1:fold) { omit <- foldi[[i]] tempX <- as.matrix(X[-omit,]) tempY <- as.matrix(Y[-omit,]) testX <- as.matrix(X[omit,]) testY <- as.matrix(Y[omit,]) for(j in 1 : length(dims)){ cat('----------fold=',i,',dims=',dims[j],'\n') #init=grams((pls::plsr(tempY~tempX,dims[j]))$projection) m0 <- qenv_gmm(tempX,tempY,dims[j],tau=tau,method=method) pred <- matrix(1, nrow(testX),1) %*% t(m0$alpha) + testX %*% m0$beta resi <- testY - pred RCVi[i,j] <-mean(resi*(tau-(resi<0))) mspemati[i,j]= mean(resi^2) } # end dims }#end fold } mspemat<- apply(mspemati, 2, mean) mspe_u = dims[which.min(mspemat)] RCV<- apply(RCVi, 2, mean) RCV_u = dims[which.min(RCV)] cv <- list(mspemat = mspemat,mspe_u=mspe_u,RCV=RCV,RCV_u=RCV_u,fit.time=(proc.time()-time1)[3]) cv }