This experiment compares Intergrated squared error of d on unifdist.rand, random design, maximin and betadist(3,6) designs

For each design realization, I (= 1000) d’s are generated from unif(0, \(\sqrt dim\))

library(plgp)
library(mvtnorm)
library(laGP)

library(foreach)
library(doParallel)
source("ud.R")
gpi.mse <- function(method, n, dim, I = 1000, seed = 0){
  set.seed(seed)
  if(method ==1){
    x <- ud.m(n, dim, k = 1, T = 10^5, starter_num = 5)$X
  }
  else if(method == 2){
    x <- matrix(runif(n*dim), ncol = dim)
  }
  else if(method == 3){
    x <- maximin(n, dim)
  }
  else if(method == 4){
    x <- bd(n, dim)$X
  }
  D <- distance(x)
  eps <- sqrt(.Machine$double.eps)
  dtrue <- runif(I, min = eps, max = sqrt(dim))
  dhat <- rep(NA, I)
  for(i in 1:I){
    sigma <- exp(-D/dtrue[i] + diag(eps, n))
    y <- rmvnorm(1, sigma = sigma)
    gpi <- newGP(x, y, d = 0.1, g = eps, dK = T)
    dhat[i] <- mleGP(gpi, param = "d", tmax = 10)$d
    deleteGP(gpi)
  }
  mse <- mean((dtrue - dhat)^2)
  return(c(n, dim, method, mse, seed))
}

A <- expand.grid(method = c(1,2,3,4), n = 8, dim = 2:10)
set <- matrix(rep(t(A),40), ncol=ncol(A), byrow=TRUE)
set <- cbind(set, 1:nrow(set))
colnames(set) <- c("method", "n", "dim", "seed")
cl <- makeCluster(8)
registerDoParallel(cl)
out <- foreach(i = 1:nrow(set), .errorhandling = "remove", .combine='rbind', .packages=c('plgp', 'mvtnorm', 'laGP' )) %dopar% { 
              gpi.mse(method =set[i,1], n = set[i,2], dim = set[i,3], seed = set[i,4]) }

stopCluster(cl)
save.image('GP_mse_rand_unifdist_cp_1.RData') 
#-------------------------------------plot-----------------------------------------------------------------------------------                                                                                                                                                                                                              
load('/home/boyazhang/repos/unifdist/code/GP_mse_rand_unifdist_cp_1.RData')
out <- as.data.frame(out)
colnames(out) <- c("n", "dim", "method", "ise", "seed")
out$method <- replace(out$method, out$method == 1, "5-start unifdist" )
out$method <- replace(out$method, out$method == 2, "random" )
out$method <- replace(out$method, out$method == 3, "maximin" )
out$method <- replace(out$method, out$method == 4, "beta(3,6)" )

par(mfcol = c(1,3))
boxplot(ise~method, data = out[out$dim == 2,], main = "dim = 2, ise comparison", ylab = "ISE")
boxplot(ise~method, data = out[out$dim == 3,], main = "dim = 3, ise comparison", ylab = "ISE")
boxplot(ise~method, data = out[out$dim == 4,], main = "dim = 4, ise comparison", ylab = "ISE")

boxplot(ise~method, data = out[out$dim == 5,], main = "dim = 5, ise comparison", ylab = "ISE")
boxplot(ise~method, data = out[out$dim == 6,], main = "dim = 6, ise comparison", ylab = "ISE")
boxplot(ise~method, data = out[out$dim == 7,], main = "dim = 7, ise comparison", ylab = "ISE")

boxplot(ise~method, data = out[out$dim == 8,], main = "dim = 8, ise comparison", ylab = "ISE")
boxplot(ise~method, data = out[out$dim == 9,], main = "dim = 9, ise comparison", ylab = "ISE")
boxplot(ise~method, data = out[out$dim == 10,], main = "dim = 10, ise comparison", ylab = "ISE")

### low-dim performance is quite different from high-dim performance, maybe 10000 iterations are not enough to get a satisfying “unifdist” or “betadist” design #### dim = 10

 udtest <- ud(8,10,1)
 bdtest <- bd(8,10)
 par(mfcol = c(2,1))
 edfPlot(udtest)

 plot(density(udtest$X))
 xx <- seq(0,1,length=100)
 plot(xx, dbeta(xx, 3,6))

#### dim = 5

 udtest <- ud(8,5,1)
 bdtest <- bd(8,5)
 par(mfcol = c(2,1))
 edfPlot(udtest)

 plot(density(udtest$X))
 xx <- seq(0,1,length=100)
 plot(xx, dbeta(xx, 3,6))

#### dim = 3

 udtest <- ud(8,3,1)
 bdtest <- bd(8,3)
 par(mfcol = c(2,1))
 edfPlot(udtest)

 plot(density(udtest$X))
 xx <- seq(0,1,length=100)
 plot(xx, dbeta(xx, 3,6))