rand_design_pattern_facts.R

## random design patterns with different d range
gpi.mse <- function(x, drange){
  D <- distance(x)
  eps <- sqrt(.Machine$double.eps)
  dtrue <- runif(I, min = drange[1], max = drange[2])
  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)
}

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

n <- 8
dim <- 2
I <- 1000

drange <- c(0,0.1)
rand_design_1 <- list()
mse1 <- rep(NA, 10000)
for(i in 1:10000){
  x <- matrix(runif(n*dim), ncol = dim)
  rand_design_1[[i]] <- x
  mse1[i] <- gpi.mse(x, drange)
}

drange <- c(1.3, sqrt(2))
rand_design_2 <- list()
mse2 <- rep(NA, 10000)
for(i in 1:10000){
  x <- matrix(runif(n*dim), ncol = dim)
  rand_design_2[[i]] <- x
  mse2[i] <- gpi.mse(x, drange)
}

drange <- c(0, sqrt(2))
rand_design_3 <- list()
mse3 <- rep(NA, 10000)
for(i in 1:10000){
  x <- matrix(runif(n*dim), ncol = dim)
  rand_design_3[[i]] <- x
  mse3[i] <- gpi.mse(x, drange)
}

save.image("rand_design_pattern_facts_1.RData")

##---------------------------------------------------------------------------plot--------------------------------------------------------------
setwd("/home/boyazhang/repos/unifdist/code")
load("rand_design_pattern_facts_1.RData")

rand1top50 <- c(1:10000)[rank(mse1)<=50]
dst1 <- dist(rand_design_1[[rand1top50[1]]])
for(i in rand1top50){
  dst1 <- c(dst1, dist(rand_design_1[[i]]))
}

rand1end50 <- c(1:10000)[rank(mse1)>=9951]
dse1 <- dist(rand_design_1[[rand1end50[1]]])
for(i in rand1end50){
  dse1 <- c(dse1, dist(rand_design_1[[i]]))
}

rand2top50 <- c(1:10000)[rank(mse2)<=50]
dst2 <- dist(rand_design_2[[rand2top50[1]]])
for(i in rand2top50){
  dst2 <- c(dst2, dist(rand_design_2[[i]]))
}

rand2end50 <- c(1:10000)[rank(mse2)>=9951]
dse2 <- dist(rand_design_2[[rand2end50[1]]])
for(i in rand2end50){
  dse2 <- c(dse2, dist(rand_design_2[[i]]))
}

rand3top50 <- c(1:10000)[rank(mse3)<=50]
dst3 <- dist(rand_design_3[[rand3top50[1]]])
for(i in rand3top50){
  dst3 <- c(dst3, dist(rand_design_3[[i]]))
}

rand3end50 <- c(1:10000)[rank(mse3)>=9951]
dse3 <- dist(rand_design_3[[rand3end50[1]]])
for(i in rand3end50){
  dse3 <- c(dse3, dist(rand_design_3[[i]]))
}

plot(density(dst1, from = 0, to =sqrt(2)),  col = 1)
lines(density(dse1, from = 0, to = sqrt(2)), col = 2)
lines(density(dst2, from = 0, to =sqrt(2)),  col = 3)
lines(density(dse2, from = 0, to =sqrt(2)),  col = 4)
lines(density(dst3, from = 0, to =sqrt(2)),  col = 5)
lines(density(dse3, from = 0, to =sqrt(2)),  col = 6)
legend("topright", c("d in (0, 0.1) top50", "d in (0, 0.1) end50", "d in (1.3, sqrt(2)) top50","d in (1.3, sqrt(2)) end50",
                     "d in (0, sqrt(2)) top50","d in (0, sqrt(2)) end50"), col = 1:6, lty = rep(1,6))

rand3top50 <- c(1:10000)[rank(mse3)<=50]
d1 <- dist(rand_design_3[[rand3top50[1]]])
for(i in rand3top50){
  d1 <- c(d1, dist(rand_design_3[[i]]))
}

rand3top50 <- c(1:10000)[rank(mse3)>50 & rank(mse3)<=100]
d2 <- dist(rand_design_3[[rand3top50[1]]])
for(i in rand3top50){
  d2 <- c(d2, dist(rand_design_3[[i]]))
}

rand3top50 <- c(1:10000)[rank(mse3)>100 & rank(mse3)<=150]
d3 <- dist(rand_design_3[[rand3top50[1]]])
for(i in rand3top50){
  d3 <- c(d3, dist(rand_design_3[[i]]))
}

rand3top50 <- c(1:10000)[rank(mse3)>5000 & rank(mse3)<=5050]
d4 <- dist(rand_design_3[[rand3top50[1]]])
for(i in rand3top50){
  d4 <- c(d4, dist(rand_design_3[[i]]))
}

rand3end50 <- c(1:10000)[rank(mse3)>=9951]
d5 <- dist(rand_design_3[[rand3end50[1]]])
for(i in rand3end50){
  d5 <- c(d5, dist(rand_design_3[[i]]))
}

plot(density(d1, from =0, to = sqrt(2)), col = 1, main = "distance density one different mse ranks")
lines(density(d2, from =0, to = sqrt(2)), col = 2)
lines(density(d3, from =0, to = sqrt(2)), col = 3)
lines(density(d4, from =0, to = sqrt(2)), col = 4)
lines(density(d5, from =0, to = sqrt(2)), col = 5)
legend("topright", c("1-50", "51-100", "101-150", "5001-5050", "9951-10000"), col = 1:5, lty = rep(1,5))