accepted version

Author: ChihLi

2d_jet_blade_true.Rdata

@@ -3,36 +3,35 @@ eps <- 1e-8
 
 # fit a Gaussian process model
 GP.fit <- function(X, y, k, g=eps, 
-                  init=rep(10, ncol(X)),
-                  lower=0.1, upper=1000){
+                   init=rep(10, ncol(X)),
+                   lower=0.1, upper=1000){
   # X is the input
   # y is output
   # k is the smoothness parameter in matern kernel
   # g is the nugget effect
   # init is the initial value of lengthscale parameter when performing optimization
   # lower is the lower bound of lengthscale parameter when performing optimization
   # upper is the upper bound of lengthscale parameter when performing optimization
-  
   if(is.null(dim(X))) X <- matrix(X, ncol = 1)
 
   y <- scale(y, scale=FALSE)
 
   n <- length(y)
   if(!is.null(k)){ # if the smoothness is given
-    nlsep <- function(par, X, Y) 
+    obj <- function(par, X, Y) 
     {
       theta <- par
       R <- sqrt(distance(t(t(X)*theta)))
       K <- matern.kernel(R, k=k)
       Ki <- solve(K+diag(g,n))
-      ldetK <- determinant(K, logarithm=TRUE)$modulus
-      ll <- - (n/2)*log(t(Y) %*% Ki %*% Y) - (1/2)*ldetK
-      return(-ll)
+
+      loocv <- mean((diag(1/diag(Ki)) %*% Ki %*% Y)^2)
+      return(loocv)
     }
 
     tic <- proc.time()[3]
 
-    outg <- optim(init, nlsep, 
+    outg <- optim(init, obj, 
                   method="L-BFGS-B", lower=lower, upper=upper, X=X, Y=y)
     toc <- proc.time()[3]
 
@@ -43,30 +42,25 @@ GP.fit <- function(X, y, k, g=eps,
     loocv.save <- rep(0,length(k.candidate))
 
     for(i in 1:length(k.candidate)){
-      nlsep <- function(par, X, Y) 
+      obj <- function(par, X, Y) 
       {
         theta <- par
         R <- sqrt(distance(t(t(X)*theta)))
         K <- matern.kernel(R, k=k.candidate[i])
         Ki <- solve(K+diag(g,n))
-        ldetK <- determinant(K, logarithm=TRUE)$modulus
-        ll <- - (n/2)*log(t(Y) %*% Ki %*% Y) - (1/2)*ldetK
-        return(-ll)
+        loocv <- mean((diag(1/diag(Ki)) %*% Ki %*% Y)^2)
+        return(loocv)
       }
 
       tic <- proc.time()[3]
 
-      outg <- optim(init, nlsep,
+      outg <- optim(init, obj,
                     method="L-BFGS-B", lower=lower, upper=upper, X=X, Y=y)
       toc <- proc.time()[3]
 
       theta.save[i,] <- outg$par
 
-      R <- sqrt(distance(t(t(X)*theta.save[i,])))
-      K <- matern.kernel(R, k=k.candidate[i])
-      Ki <- solve(K+diag(g,n))
-      
-      loocv.save[i] <- mean((diag(1/diag(Ki)) %*% Ki %*% y)^2)
+      loocv.save[i] <- outg$value
     }
     k <- k.candidate[which.min(loocv.save)]
     theta <- theta.save[which.min(loocv.save),]

GP.R

@@ -1,6 +1,8 @@
-"X2...X...1","rep.MM.tt.l..n."
-0,0.0125
-0.5,0.0125
--0.5,0.0125
--0.25,0.0125
-0.75,0.0125
+"X1","X2","rep.MM.tt.l..n."
+0.453125,0.703125,0.05
+0.390625,0.390625,0.05
+0.640625,0.640625,0.05
+0.515625,0.265625,0.05
+0.265625,0.515625,0.05
+0.2734375,0.3828125,0.05
+0.5234375,0.6328125,0.05

Rmatlab_files/.DS_Store

@@ -1,5 +1,7 @@
-0,0.012500000000000000693889390390723,0.39971446424469275759960851246433,1.311069999999999957651652948698
-0.5,0.012500000000000000693889390390723,0.51639448598742210272405372961657,0.56689658300000000945573219723883
--0.5,0.012500000000000000693889390390723,0.31309325845750018979529727403133,0.49697162499999997242738913882931
--0.25,0.012500000000000000693889390390723,0.35323817107758859190624889379251,0.40993395900000001441654262634984
-0.75,0.012500000000000000693889390390723,0.58955524643800005790694740426261,0.56776395899999998473361983997165
+0.453125,0.703125,0.050000000000000002775557561562891,21.12135924143782617079523333814,1.6989625420000000755038627175963
+0.390625,0.390625,0.050000000000000002775557561562891,4.4355911525197075917503752862103,0.84547937500000003296918293926865
+0.640625,0.640625,0.050000000000000002775557561562891,7.2743694901323232571144217217807,0.81774108300000003524132807797287
+0.515625,0.265625,0.050000000000000002775557561562891,18.263094761612002514539199182764,0.70163470800000005134222647029674
+0.265625,0.515625,0.050000000000000002775557561562891,19.733075703317034310657618334517,0.9584084579999999631283458256803
+0.2734375,0.3828125,0.050000000000000002775557561562891,9.8450206901014389870852028252557,0.7251950000000000340349970429088
+0.5234375,0.6328125,0.050000000000000002775557561562891,12.053113630961497904081625165418,0.67546404100000001502479562986991