tests/testthat/test-loo.R

In DiceKriging: Kriging Methods for Computer Experiments

x <- c(0, 0.4, 0.6, 0.8, 1)
y <- c(-0.3, 0, -0.8, 0.5, 0.9)
n <- length(x)

theta <- 0.01; sigma <- 3;

list.type <- list("SK", "SK", "UK", "UK")
list.trend.reestim <- list(FALSE, TRUE, FALSE, TRUE)

list.case <- list("", "with known nugget") #, "with noisy observations")
list.nugget <- list(NULL, sigma/10) #, NULL)
list.noise.var <- list(NULL, NULL) #, 1:n)

precision <- 1e-6


for (b in 1:length(list.type)) {
  if (list.trend.reestim[[b]]) {
    trend <- NULL
  } else {
    trend <- c(-1,2)
  }
    
for (a in 1:length(list.case)) {
  
  case <- paste(list.type[[b]], list.case[[a]], "; trend.reestim=", list.trend.reestim[[b]])
  print(case)
  
  m <- km(~x, design=data.frame(x=x), response=data.frame(y=y), 
          covtype="matern5_2", coef.trend=trend, coef.cov=theta, coef.var=sigma^2,
          nugget=list.nugget[[a]], noise.var=list.noise.var[[a]])

  loo.mean <- loo.sd <- rep(NA, n)
  for (i in 1:n){
    mloo <- km(~x, design=data.frame(x=x[-i]), response=data.frame(y=y[-i]), 
               covtype="matern5_2", coef.trend=trend, coef.cov=theta, coef.var=sigma^2,
               nugget=list.nugget[[a]], noise.var=list.noise.var[[a]][-i])
    p <- predict(mloo, newdata=x[i], type=list.type[[b]], checkNames=FALSE)
    loo.mean[i] <- p$mean
    loo.sd[i] <- p$sd
  }

  loo <- leaveOneOut.km(m, type=list.type[[b]], trend.reestim=list.trend.reestim[[b]])
  
  test_that(desc=paste("Check LOO mean for", list.type[[b]], list.case[[a]]), expect_true(max(abs(loo$mean-loo.mean)/loo.mean) < precision))
  test_that(desc=paste("Check LOO sd for", list.type[[b]], list.case[[a]]), expect_true(max(abs(loo$sd-loo.sd)/loo.sd) < precision))
  sum((loo$mean - loo.mean)^2)
  sum((loo$sd - loo.sd)^2)
}
  
}