R/LOO.R

In maatouk/constrKriging: Kriging with inequality constraints

coef.cov_LOO =function(object) {

  fctGamma = object$fctGamma
  input=object$call$design
  output=object$call$response
  n=length(output)
  
  fctinvGamma <- function(theta){
    invGamma1 = chol(fctGamma(theta))
    invGamma = chol2inv(invGamma1)
    return(invGamma)
  }
  
  Ai <- function(i){
    object$A[-i,]
  }
  
  Amat1 <- diag(nrow(object$Gamma))
  Amat1[1,1] <- 0
  Amati <- function(i){
    rbind(Ai(i),Amat1)
  }
  
  fctzetoil <- function(theta,i){
    solve.QP(fctinvGamma(theta),dvec=rep(0,nrow(object$Gamma)),Amat=t(Amati(i)),bvec=c(output[-i],rep(0,nrow(object$Gamma))),meq=n-1)$solution
  }
  
  fctoptim <- function(theta){
    sum = 0
    for(i in 1 : length(input)){
      sum = sum + (output[i]-(object$A%*%fctzetoil(theta,i))[i])^2
    }
    return(sum/length(input))
  }
  
  fctoptim.try =  function(theta) {
    f=NaN
    try(f <- fctoptim(theta))
    if(is.nan(f)) cat("(!) ",theta," -> ",f)
    return(f)
  }
  
  range_design=range(input)
  lower = min(dist(input))*1E-5
  upper = (range_design[2]-range_design[1])*10
  
  #plot(Vectorize(fctoptim.try), xlim=c(lower,upper))
  
  ## 1-D optim, slow to converge, and often leads to high theta
  #tol=(upper-lower)/1000
  #min = optimize(fctoptim.try,lower=lower,upper = upper,tol = tol)$minimum
  #min = optim(par = object$call$coef.cov,fctoptim.try,method = "L-BFGS-B",lower=lower,upper =upper,abstol=tol)$par
  
  ## search for a not(too-hight theta, which gives a reasonable LOO. Faster and more stable.
  thetas=seq(lower,upper,,20)
  mins = Vectorize(fctoptim.try)(thetas)
  if(all(is.nan(mins))) stop("Could not evaluate LOO at any theta !")
  #plot(thetas,mins)
  lower_bound_mins = min(mins)+(max(mins)-min(mins))*.01
  lower_theta = min(thetas[which(mins<lower_bound_mins)])
  
  return(lower_theta)
}