load("/home/boyazhang/repos/unifdist/code/beta_para_map_16_2.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
## 37 2.216328 2.675833 0.04196201 0.02827014
## 77 4.894074 5.506707 0.04288450 0.02021158
## 71 4.136959 5.149799 0.04339780 0.02543473
## 45 3.339937 3.160866 0.04519418 0.02263051
## 63 3.252284 4.683300 0.04578504 0.03790452
## 67 5.580657 4.931695 0.04640385 0.02567290
## 55 4.412013 4.057383 0.04664979 0.03158139
## 56 2.796066 4.058591 0.04796092 0.03922280
## 83 6.469290 5.963655 0.04866442 0.03960064
## 97 6.071751 7.043815 0.05121782 0.03746309
## 86 4.525203 6.237224 0.05160917 0.03746184
## 32 3.294018 2.383757 0.05452840 0.04464137
## 114 8.213256 8.384534 0.05578339 0.04471785
## 109 7.235359 8.058196 0.05786790 0.04206288
## 24 1.779214 1.627886 0.05930125 0.05200242
## 91 5.387115 6.720263 0.05952749 0.04847739
## 102 6.573748 7.681288 0.05987341 0.04873732
## 46 4.518677 3.216644 0.06088164 0.04020663
## 139 10.833878 10.426468 0.06123205 0.03973242
## 135 14.619250 10.010862 0.06185755 0.04314956
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.867616 2.201087
## Variance/scale hyperparameter: 1.116138
## Summary of Lambda values:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.05346 0.13229 0.23473 0.24236 0.32518 0.62202
## Estimated constant trend value: -1.995137
## Matern5_2 covariance lengthscale values of the log-noise process: 1.867616 2.201087
## Nugget of the log-noise process: 1e-06
## Estimated constant trend value of the log-noise process: -1.640833
## MLE optimization:
## Log-likelihood = -14202.01 ; Nb of evaluations (obj, gradient) by L-BFGS-B: 86 86 ; 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 = "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
## 38160 2.49 2.28 -3.67841 0.008453076
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")
