load("/home/boyazhang/repos/unifdist/code/beta_para_map_32_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
## 108 4.857314 7.830210 0.06925382 0.006400629
## 99 4.563656 7.094364 0.06983989 0.006013910
## 89 4.728287 6.129442 0.07285811 0.006789681
## 141 6.168468 10.295455 0.07316357 0.008413966
## 171 6.804890 12.729438 0.07370263 0.007156088
## 152 6.802167 11.219320 0.07476046 0.007692728
## 175 7.747490 13.204800 0.07561378 0.008057637
## 194 7.964042 14.738094 0.07583452 0.007836511
## 129 6.541631 9.418658 0.07638154 0.006996699
## 98 5.419220 7.072942 0.07670117 0.007230928
## 73 4.458450 4.878517 0.07880667 0.007755214
## 163 8.040653 12.076008 0.07883051 0.008290077
## 173 8.682310 12.746455 0.07905545 0.008400883
## 185 8.530622 13.839621 0.07912558 0.006971934
## 148 7.452969 10.819510 0.07934247 0.006991802
## 183 9.507787 13.567227 0.07952101 0.006234426
## 186 10.250136 14.091506 0.08036477 0.007570373
## 170 9.806924 12.688982 0.08049283 0.008944372
## 106 6.308692 7.542981 0.08053234 0.008990862
## 131 7.446944 9.616523 0.08102313 0.008457208
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.254951 2.054654
## Variance/scale hyperparameter: 1.218393
## Summary of Lambda values:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.004378 0.008209 0.072504 0.367553 0.578538 1.864967
## Estimated constant trend value: -0.7521984
## Matern5_2 covariance lengthscale values of the log-noise process: 1.254951 2.054654
## Nugget of the log-noise process: 1e-06
## Estimated constant trend value of the log-noise process: -2.436364
## MLE optimization:
## Log-likelihood = -5701.927 ; Nb of evaluations (obj, gradient) by L-BFGS-B: 58 58 ; 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=32,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
## 53711 3.09 3.21 -3.777628 0.0998953
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")
