load("/home/boyazhang/repos/unifdist/code/beta_para_map_16_3.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
## 56 2.420988 4.346548 0.1274603 0.01993396
## 151 4.772067 11.499783 0.1275764 0.01507356
## 79 3.046271 6.163298 0.1284363 0.02037836
## 137 4.400685 10.453686 0.1287662 0.01838880
## 171 4.980460 13.258985 0.1313422 0.02138326
## 157 4.956067 12.263731 0.1316991 0.02016727
## 84 3.631188 6.698245 0.1340903 0.01851228
## 167 4.280799 12.921115 0.1342958 0.07069001
## 98 4.257004 7.657890 0.1352799 0.02351583
## 27 1.509478 1.859418 0.1356197 0.02152429
## 178 5.718563 13.545401 0.1370476 0.02114637
## 111 4.360240 8.541571 0.1373413 0.02409659
## 193 5.652910 14.774445 0.1377309 0.06161640
## 130 5.194884 9.978201 0.1406753 0.01914846
## 153 5.714754 11.684453 0.1413412 0.02173378
## 117 5.096528 9.121518 0.1419350 0.03137810
## 189 6.609499 14.340459 0.1422058 0.02175157
## 108 5.109193 8.286095 0.1426439 0.02634031
## 71 3.578203 5.308953 0.1457548 0.02630696
## 89 4.508022 6.858183 0.1469942 0.02346975
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.155082 2.249499
## Variance/scale hyperparameter: 0.3219894
## Summary of Lambda values:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.04570 0.06394 0.07787 0.42702 0.66055 2.71455
## Estimated constant trend value: -1.20324
## Matern5_2 covariance lengthscale values of the log-noise process: 1.155082 2.249499
## Nugget of the log-noise process: 1e-06
## Estimated constant trend value of the log-noise process: -1.418743
## MLE optimization:
## Log-likelihood = -795.7699 ; Nb of evaluations (obj, gradient) by L-BFGS-B: 51 51 ; 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=16,dim=3: 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
## 83258 2.73 4.98 -2.237657 0.01113246
image(x = xx, y = xx, z = matrix(predictions$sd2, ncol = length(xx)), xlab="shape1",ylab="shape2", col=cols, main = "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
plot(mse$beta_m, mse$y, xlab = "shape1/(shape1 + shape2)", ylab = "log(mse_mean)")
# 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")
