Error in LMP/linear trend about libkriging HOT 1 CLOSED

# 1D function, small design
X <- as.matrix(c(0.0, 0.25, 0.5, 0.75, 1.0))
f <- function(x) 1 - 1 / 2 * (sin(12 * x) / (1 + x) + 2 * cos(7 * x) * x^5 + 0.7)
y <- f(X) # + 10*rnorm(nrow(X))

library(rlibkriging)
k <- Kriging(y, X, kernel="matern5_2", objective = "LMP", optim = "BFGS", regmodel = "linear")
plot( function(t) k$logMargPostFun(t)$logMargPost ,xlab=expression(theta))
abline(v=k$theta(), col="red")

# LL
obj_fun = function(theta,...) logLikelihoodFun(k,theta)$logLikelihood
obj_grad = function(theta,...) logLikelihoodFun(k,theta,grad = T)$logLikelihoodGrad
obj_grad_fd = function(theta,...) (obj_fun(theta+0.01)-obj_fun(theta-0.01))/0.02

plot(Vectorize(obj_fun),col='blue')
for (x in seq(0.01,1,,11)){
  fx = obj_fun(x)
  gx = obj_grad(x)
  gx_ref = obj_grad_fd(x)
  arrows(x,fx,x+.1,fx+.1*gx,col='red')
  arrows(x,fx,x+.11,fx+.11*gx_ref,col='blue')
}

# LMP
obj_fun = function(theta,...) logMargPostFun(k,theta)$logMargPost
obj_grad = function(theta,...) logMargPostFun(k,theta,grad = T)$logMargPostGrad
obj_grad_fd = function(theta,...) (obj_fun(theta+0.01)-obj_fun(theta-0.01))/0.02

plot(Vectorize(obj_fun),col='blue')
for (x in seq(0.01,1,,11)){
  fx = obj_fun(x)
  gx = obj_grad(x)
  gx_ref = obj_grad_fd(x)
  arrows(x,fx,x+.1,fx+.1*gx,col='red', lty=2)
  arrows(x,fx,x+.11,fx+.11*gx_ref,col='blue')
}

library(RobustGaSP)
r = rgasp(design=X,response = y,kernel_type="matern_5_2",optimization = "lbfgs",trend=cbind(matrix(1,length(y),1),X))

obj_fun_rgasp = function(theta,...) 
  -RobustGaSP:::neg_log_marginal_post_approx_ref(param=log(1/theta),nugget=0, [email protected], 
                                              R0=r@R0,X=r@X, zero_mean=r@zero_mean,output=r@output, 
                                              CL=r@CL, 
                                              a=0.2,
                                              b=1/(length(r@output))^{1/dim(as.matrix(r@input))[2]}*(0.2+dim(as.matrix(r@input))[2]),
                                              kernel_type=rep(as.integer(3),ncol(X)),alpha=r@alpha)

obj_grad_rgasp = function(theta,...) 
  -RobustGaSP:::neg_log_marginal_post_approx_ref_deriv(param=log(1/theta),nugget=0, [email protected], 
                                                 R0=r@R0,X=r@X, zero_mean=r@zero_mean,output=r@output, 
                                                 CL=r@CL, 
                                                 a=0.2,
                                                 b=1/(length(r@output))^{1/dim(as.matrix(r@input))[2]}*(0.2+dim(as.matrix(r@input))[2]),
                                                 kernel_type=rep(as.integer(3),ncol(X)),alpha=r@alpha) / -theta

lines(seq(0,1,,101),(Vectorize(obj_fun_rgasp)(seq(0,1,,101))),col='red')
abline(v=1/r@beta_hat,col='red')
for (x in seq(0.01,1,,11)){
  fx = obj_fun_rgasp(x)
  gx = obj_grad_rgasp(x)
  arrows(x,fx,x+.1,fx+.1*gx,col='red', lty=1)
}