This experiment compares intergrated squared error (ISE) of d, d \(\in (0,1)^{dim}\), on beta distributions with different parameters, dim = 2,3, …, 10
par(mfcol = c(2,3))
xx <- seq(0,1,length = 100)
plot(xx, dbeta(xx, 5, 2), type = "l", main = "beta(5,2)")
plot(xx, dbeta(xx, 2, 2), type = "l", main = "beta(2,2)")
plot(xx, dbeta(xx, 2, 5), type = "l", main = "beta(2,5)")
plot(xx, dbeta(xx, 3, 6), type = "l", main = "beta(3,6)")
plot(xx, dbeta(xx, 6, 3), type = "l", main = "beta(6,3)")
plot(xx, dbeta(xx, 1/2, 1/2), type = "l", main = "beta(1/2,1/2)")
### This experiment compares intergrated squared error (ISE) of d, d $\in (0,1)^{dim}$, on beta distributions with different parameters, dim = 2,3,4
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 <- bd(n=n, m=dim, shape1 = 5,shape2 = 2)$X
}
else if(method == 2){
x <- bd(n=n, m=dim, shape1 = 2,shape2 = 2)$X
}
else if(method == 3){
x <- bd(n=n, m=dim, shape1 = 2,shape2 = 5)$X
}
else if(method == 4){
x <- bd(n=n, m=dim, shape1 = 3,shape2 = 6)$X
}
else if(method == 5){
x <- bd(n=n, m=dim, shape1 = 6,shape2 = 3)$X
}
else if(method == 6){
x <- bd(n=n, m=dim, shape1 = .5,shape2 = .5)$X
}
else if(method == 7){
x <- matrix(runif(n * dim), ncol = dim)
}
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 = 1:7, n = c(8,16), dim = 2:4)
set <- matrix(rep(t(A),100), 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('beta_para.RData')
#-------------------------------------plot-----------------------------------------------------------------------------------
load('/home/boyazhang/repos/unifdist/code/beta_para.RData')
out <- as.data.frame(out)
colnames(out) <- c("n", "dim", "method", "ise")
out$method <- replace(out$method, out$method == 1, "beta(5,2)" )
out$method <- replace(out$method, out$method == 2, "beta(2,2)" )
out$method <- replace(out$method, out$method == 3, "beta(2,5)" )
out$method <- replace(out$method, out$method == 4, "beta(3,6)" )
out$method <- replace(out$method, out$method == 5, "beta(6,3)" )
out$method <- replace(out$method, out$method == 6, "beta(0.5,0.5)" )
out$method <- replace(out$method, out$method == 7, "random" )
par(mfcol = c(1,1))
boxplot(ise~method, data = out[out$dim == 2 & out$n == 8,], main = "dim = 2, n=8, ise comparison", ylab = "ISE")
boxplot(ise~method, data = out[out$dim == 3 & out$n == 8,], main = "dim = 3, n=8, ise comparison", ylab = "ISE")
boxplot(ise~method, data = out[out$dim == 4 & out$n == 8,], main = "dim = 4, n=8, ise comparison", ylab = "ISE")
boxplot(ise~method, data = out[out$dim == 2 & out$n == 16,], main = "dim = 2, n=16, ise comparison", ylab = "ISE")
boxplot(ise~method, data = out[out$dim == 3 & out$n == 16,], main = "dim = 3, n=16, ise comparison", ylab = "ISE")
boxplot(ise~method, data = out[out$dim == 4 & out$n == 16,], main = "dim = 4, n=16, ise comparison", ylab = "ISE")
par(mfcol = c(2,4))
boxplot(ise~n, data = out[out$dim == 2 & out$method == "beta(0.5,0.5)",], main = "dim = 2, beta(0.5,0.5), ise comparison", ylab = "ISE")
boxplot(ise~n, data = out[out$dim == 2 & out$method == "beta(5,2)",], main = "dim = 2, beta(5,2), ise comparison", ylab = "ISE")
boxplot(ise~n, data = out[out$dim == 2 & out$method == "beta(2,5)",], main = "dim = 2, beta(2,5), ise comparison", ylab = "ISE")
boxplot(ise~n, data = out[out$dim == 2 & out$method == "beta(6,3)",], main = "dim = 2, beta(6,3), ise comparison", ylab = "ISE")
boxplot(ise~n, data = out[out$dim == 2 & out$method == "beta(3,6)",], main = "dim = 2, beta(3,6), ise comparison", ylab = "ISE")
boxplot(ise~n, data = out[out$dim == 2 & out$method == "beta(2,2)",], main = "dim = 2, beta(2,2), ise comparison", ylab = "ISE")
boxplot(ise~n, data = out[out$dim == 2 & out$method == "random",], main = "dim = 2, random, ise comparison", ylab = "ISE")
par(mfcol = c(2,4))
boxplot(ise~n, data = out[out$dim == 3 & out$method == "beta(0.5,0.5)",], main = "dim = 3, beta(0.5,0.5), ise comparison", ylab = "ISE")
boxplot(ise~n, data = out[out$dim == 3 & out$method == "beta(5,2)",], main = "dim = 3, beta(5,2), ise comparison", ylab = "ISE")
boxplot(ise~n, data = out[out$dim == 3 & out$method == "beta(2,5)",], main = "dim = 3, beta(2,5), ise comparison", ylab = "ISE")
boxplot(ise~n, data = out[out$dim == 3 & out$method == "beta(6,3)",], main = "dim = 3, beta(6,3), ise comparison", ylab = "ISE")
boxplot(ise~n, data = out[out$dim == 3 & out$method == "beta(3,6)",], main = "dim = 3, beta(3,6), ise comparison", ylab = "ISE")
boxplot(ise~n, data = out[out$dim == 3 & out$method == "beta(2,2)",], main = "dim = 3, beta(2,2), ise comparison", ylab = "ISE")
boxplot(ise~n, data = out[out$dim == 3 & out$method == "random",], main = "dim = 3, random, ise comparison", ylab = "ISE")
par(mfcol = c(2,4))
boxplot(ise~n, data = out[out$dim == 4 & out$method == "beta(0.5,0.5)",], main = "dim = 4, beta(0.5,0.5), ise comparison", ylab = "ISE")
boxplot(ise~n, data = out[out$dim == 4 & out$method == "beta(5,2)",], main = "dim = 4, beta(5,2), ise comparison", ylab = "ISE")
boxplot(ise~n, data = out[out$dim == 4 & out$method == "beta(2,5)",], main = "dim = 4, beta(2,5), ise comparison", ylab = "ISE")
boxplot(ise~n, data = out[out$dim == 4 & out$method == "beta(6,3)",], main = "dim = 4, beta(6,3), ise comparison", ylab = "ISE")
boxplot(ise~n, data = out[out$dim == 4 & out$method == "beta(3,6)",], main = "dim = 4, beta(3,6), ise comparison", ylab = "ISE")
boxplot(ise~n, data = out[out$dim == 4 & out$method == "beta(2,2)",], main = "dim = 4, beta(2,2), ise comparison", ylab = "ISE")
boxplot(ise~n, data = out[out$dim == 4 & out$method == "random",], main = "dim = 4, random, ise comparison", ylab = "ISE")
