load("/home/boyazhang/repos/unifdist/code/beta_para_map_64_4.RData")
out <- as.data.frame(out)
colnames(out) <- c("n", "dim", "shape1", "shape2","ise")
mse <- aggregate(out$ise, by = list(shape1 = out$shape1, shape2 = out$shape2), FUN = mean)
sd <- aggregate(out$ise, by = list(shape1 = out$shape1, shape2 = out$shape2), FUN = sd)
mse <- cbind(mse, sd = sd$x)
colnames(mse) <- c("shape1", "shape2", "mse_mean", "mse_sd")
# output table sorted my mse
head(mse[order(mse$mse_mean),], 20)
## shape1 shape2 mse_mean mse_sd
## 75 5.854651 5.676395 0.1849375 0.3608567
## 67 5.094195 4.882296 0.2024808 0.3217116
## 62 6.576493 4.755063 0.2728293 0.3796693
## 80 6.698160 6.147944 0.3102960 0.4794132
## 89 7.634800 6.622671 0.3510862 0.5071781
## 50 5.097106 3.639907 0.4115033 0.5373890
## 76 7.736977 5.749178 0.4689290 0.6081755
## 92 8.845609 7.022838 0.5515619 0.7286760
## 90 5.632049 6.798564 0.6228796 0.6658565
## 100 7.787993 7.865598 0.6799258 0.7772241
## 104 8.739144 8.026981 0.6804976 0.7216081
## 74 4.880530 5.674156 0.8185374 0.5395742
## 111 9.379400 8.545423 0.9034170 0.8502848
## 114 8.355985 8.808863 0.9464322 0.8653307
## 116 7.293426 8.852920 0.9542185 0.7837166
## 108 6.228584 8.193279 1.1441909 0.7530847
## 96 9.968921 7.343646 1.2081063 1.1204892
## 123 9.016193 9.665028 1.2798873 0.9491312
## 107 10.564148 8.190308 1.5733895 1.1088938
## 122 7.602228 9.616542 1.5920360 0.9993876
library(hetGP)
X <- as.matrix(out[,3:4])
Z <- out$ise
nvar = 2
## Model fitting
# settings <- list(return.hom = TRUE) # To keep homoskedastic model used for training
# model <- mleHetGP(X = X, Z = Z, lower = rep(10^(-4), nvar), upper = rep(50, nvar),
# covtype = "Gaussian", settings = settings)
model <- mleHetGP(X = X, Z = log(Z), lower = rep(10^(-4), 2), upper = rep(50, 2),
covtype = "Matern5_2", maxit=1999)
## A quick view of the fit
summary(model)
## N = 20000 n = 200 d = 2
## Matern5_2 covariance lengthscale values of the main process: 1.376277 1.587044
## Variance/scale hyperparameter: 0.5115415
## Summary of Lambda values:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.07144 0.10532 0.16298 0.75100 0.44318 6.87966
## Estimated constant trend value: 1.163563
## Matern5_2 covariance lengthscale values of the log-noise process: 1.977682 2.28055
## Nugget of the log-noise process: 1e-06
## Estimated constant trend value of the log-noise process: -1.528281
## MLE optimization:
## Log-likelihood = -9172.687 ; Nb of evaluations (obj, gradient) by L-BFGS-B: 94 94 ; message: CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH
xx <- seq(0, 15, length.out = 501)
xgrid <- as.matrix(expand.grid(xx,xx))
predictions <- predict(x = xgrid, object = model)
par(mfrow=c(1,1)); cols <- heat.colors(10)
## Display mean predictive surface
image(x = xx, y = xx, z = matrix(predictions$mean, ncol = length(xx)), xlab="shape1",ylab="shape2", col=cols,
main = "n=64,dim=4: mean surface with raw data ranks and lowest 20 prediction locs")
text(mse$shape1, mse$shape2, labels = rank(mse$mse_mean))
dat_pred <- data.frame(shape1 = xgrid[,1], shape2 = xgrid[,2], pred_mean = predictions$mean, pred_sd2 = predictions$sd2)
dat_pred_top <- head(dat_pred[order(dat_pred$pred_mean),], 20)
points(dat_pred_top[,1], dat_pred_top[,2])
head(dat_pred_top,1)
## shape1 shape2 pred_mean pred_sd2
## 83357 5.7 4.98 -3.254058 0.03481939
image(x = xx, y = xx, z = matrix(predictions$sd2, ncol = length(xx)), xlab="shape1",ylab="shape2", col=cols, main = "n=64,dim=4: sd2 surface")
text(mse$shape1, mse$shape2, labels = rank(mse$mse_sd))
#------------------------------------linear regression-------------------------------
mse <- transform(mse,beta_m = shape1/(shape1+shape2), beta_var = shape1*shape2/(shape1+shape2+1)/(shape1+shape2)^2, y = log(mse$mse_mean))
## plot beta mean v.s log(mse_mean)
# with raw data
par(mfcol = c(1,2))
plot(mse$beta_m, mse$y,xlim = c(0,1), xlab = "shape1/(shape1 + shape2)", ylab = "log(mse_mean)", main = "n=64, dim=4: log_mse v.s. mean of target beta distr")
# with gp predictions
plot(xgrid[,1]/(xgrid[,1]+xgrid[,2]+.Machine$double.eps),predictions$mean,
xlab = "shape1/(shape1 + shape2)", ylab = "log(mse_mean)")
plot(mse$y, mse$mse_sd, xlab = "log(mse_mean)", ylab = "mse_sd")
