sobol_grnn.R

data(Boston, package = "MASS")

grnn.fit <- function(x, y, sigma) {
  return(grnn::smooth(grnn::learn(data.frame(y, x)), sigma))
}

grnn.predict <- function(nn, x) {
  c <- parallel::detectCores() - 1
  return(do.call(rbind, 
                 parallel::mcMap(function(i) grnn::guess(nn, as.matrix(x[i, ])), 
                                 1:nrow(x), mc.cores = c))[,1])
}  

r2 <- function(act, pre) {
  rss <- sum((pre - act) ^ 2)
  tss <- sum((act - mean(act)) ^ 2)
  return(1 - rss / tss)
}

grnn.cv <- function(nn, sigmas, nfolds, seed) {
  dt <- nn$set
  set.seed(seed)
  folds <- caret::createFolds(1:nrow(dt), k = nfolds, list = FALSE)
  cv <- function(s) {
    r <- do.call(rbind,
                 lapply(1:nfolds,
                        function(i) data.frame(Ya = nn$Ya[folds == i], 
                                               Yp = grnn.predict(grnn.fit(nn$Xa[folds != i, ], nn$Ya[folds != i], s), 
                                                                 data.frame(nn$Xa[folds == i,])))))
    return(data.frame(sigma = s, R2 = r2(r$Ya, r$Yp)))
  }
  r2_lst <- Reduce(rbind, Map(cv, sigmas))
  return(r2_lst[r2_lst$R2 == max(r2_lst$R2), ])
}

gen_sobol <- function(min, max, n, seed) {
  return(round(min + (max - min) * randtoolbox::sobol(n, dim = 1, scrambling = 1, seed = seed), 4))
}

gen_unifm <- function(min, max, n, seed) {
  set.seed(seed)
  return(round(min + (max - min) * runif(n), 4))
}

net <- grnn.fit(Boston[, -14], Boston[, 14], sigma = 2)
                        
sobol_out <- Reduce(rbind, Map(function(x) grnn.cv(net, gen_sobol(5, 10, 10, x), 4, 2019), seq(1, 10)))

unifm_out <- Reduce(rbind, Map(function(x) grnn.cv(net, gen_unifm(5, 10, 10, x), 4, 2019), seq(1, 10)))

out <- rbind(cbind(type = rep("sobol", 10), sobol_out), 
             cbind(type = rep("unifm", 10), unifm_out)) 

boxplot(R2 ~ type, data = out, main = "Sobol Sequence vs. Uniform Random", 
        ylab = "CV RSquare", xlab = "Sequence Type")