distance_dists

100 realizations of 8-point 2-d random design; I merge all distances together and plot the distance density, in order to look at the mean curve of all densities

setwd("/home/boyazhang/repos/unifdist/code")
source('ud.R')

## Loading required package: mvtnorm

## Loading required package: tgp

## 
## Attaching package: 'EnvStats'

## The following objects are masked from 'package:stats':
## 
##     predict, predict.lm

## The following object is masked from 'package:base':
## 
##     print.default

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

## 
## Attaching package: 'laGP'

## The following object is masked from 'package:plgp':
## 
##     distance

ds1 <- NULL
n <- 8
dim <- 2
for(i in 1:100){
  x <- matrix(runif(n*dim), ncol = dim)
  ds1 <- c(ds1, dist(x)) 
}
plot(density(ds1, from = min(ds1), to = max(ds1)))

100 realizations of 8-point 2-d maximin design; I also merge all distances together and plot the distance density, in order to look at the mean curve of all densities

ds2 <- NULL
n <- 8
dim <- 2
for(i in 1:100){
  x <- maximin(n, dim)
  ds2 <- c(ds2, dist(x)) 
}
plot(density(ds2, from = min(ds2), to = max(ds2)))

20 realizations of 8-point 2-d unifdist design with random initial

ds3 <- NULL
n <- 8
dim <- 2
for(i in 1:10){
  x <- ud(n,dim,1,10^5,1)$X
  ds3 <- c(ds3, dist(x)) 
}
plot(density(ds3, from = min(ds3), to = max(ds3)))

20 realizations of 8-point 2-d unifdist design with maximin initial

ds4 <- NULL
n <- 8
dim <- 2
for(i in 1:10){
  x <- ud.1(n,dim,1,10^5,1)$X
  ds4 <- c(ds4, dist(x)) 
}
plot(density(ds4, from = min(ds4), to = max(ds4)))

50 better realizations of 8-point 2-d random design with mse < 0.15

n <- 8
dim <- 2
I <- 1000
gpi.mse <- function(x){
  D <- distance(x)
  eps <- sqrt(.Machine$double.eps)
  dtrue <- runif(I)
  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(mse)
}

rand_design <- list()
mser <- rep(NA, 50)
ks <- rep(NA, 50)
i <- 1
tot <- 0
du <- seq(from = .Machine$double.eps , to = sqrt(dim), length.out = 1000)  ## true uniform distances
while(i <= 50){
  tot <- tot+1
  x <- matrix(runif(n*dim), ncol = dim)
  mse1 <- gpi.mse(x)
  if(mse1 < 0.15){
    mser[i] <- mse1
    rand_design[[i]] <- x
    dx <- dist(x)
    ks[i] <- ks.test(dx, du)$statistic
    i <- i+1
  }
}

ds5 <- NULL
for(i in 1:50){
  x <- rand_design[[i]]
  ds5 <- c(ds5, dist(x)) 
}
plot(density(ds5, from = min(ds5), to = max(ds5)))

plot(density(ds1, from = min(ds1), to = max(ds1)), ylim = c(0,3))
lines(density(ds2, from = min(ds2), to = max(ds2)), col = 2)
lines(density(ds3, from = min(ds3), to = max(ds3)), col = 3)
lines(density(ds4, from = min(ds4), to = max(ds4)), col = 4)
lines(density(ds5, from = min(ds5), to = max(ds5)), col = 5)

legend("topright", c("random", "maximin", "unifdist_rand", "unifdist_maximin", "better_rand_design"), col = 1:5, lty = rep(1,5))