# R/util.R In laGP: Local Approximate Gaussian Process Regression

#### Documented in distance

```#*******************************************************************************
#
# Local Approximate Gaussian Process Regression
# Copyright (C) 2013, The University of Chicago
#
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
# License as published by the Free Software Foundation; either
# version 2.1 of the License, or (at your option) any later version.
#
# This library is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
# Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public
# License along with this library; if not, write to the Free Software
# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
#
# Questions? Contact Robert B. Gramacy (rbg@vt.edu)
#
#*******************************************************************************

## rbetter_R:
##
## sample uniformly
## from the (diagonal) half of a rectable where the sum of
## the rows of X (the sample) are less than ystar

rbetter <- function(n, rect, ybest)
{
rr <- .C("rbetter_R",
n = as.integer(n),
m = as.integer(nrow(rect)),
rect = as.double(rect),
ystar = as.double(ybest),
X = double(n*nrow(rect)),
PACKAGE = "laGP")

return(matrix(rr\$X, ncol=nrow(rect), byrow=TRUE))
}

## distance:
##
## calculate the distance matrix between the rows of X1
## and X2, or between x1 and itself when X2=NULL

distance <- function(X1, X2=NULL)
{
## coerse arguments and extract dimensions
X1 <- as.matrix(X1)
n1 <- nrow(X1)
m <- ncol(X1)

if(is.null(X2)) {

## calculate D
outD <- .C("distance_symm_R",
X = as.double(t(X1)),
n = as.integer(n1),
m = as.integer(m),
D = double(n1 * n1),
PACKAGE = "laGP")

## return the distance matrix
return(matrix(outD\$D, ncol=n1, byrow=TRUE))

} else {

## coerse arguments and extract dimensions
X2 <- as.matrix(X2)
n2 <- nrow(X2)

## check inputs
if(ncol(X1) != ncol(X2)) stop("col dim mismatch for X1 & X2")

## calculate D
outD <- .C("distance_R",
X1 = as.double(t(X1)),
n1 = as.integer(n1),
X2 = as.double(t(X2)),
n2 = as.integer(n2),
m = as.integer(m),
D = double(n1 * n2),
PACKAGE = "laGP")

## return the distance matrix
return(matrix(outD\$D, ncol=n2, byrow=TRUE))
}
}

## llikGP.d:
##
## build a GP with X,Z,g,d and return the
## log likelihood, and the derivatives

llikGP.d <- function(dg, X=NULL, y=NULL, gpi.in=NULL,
ab=c(0,0), param=c("d", "g"))
{
## return error for negative values
## if(any(dg <= 1e-6)) return(-Inf)

## possibly init GP
if(is.null(gpi.in)) gpi <- newGP(X,y,dg[1,1],dg[1,2], dK=TRUE)
else {
gpi <- gpi.in
newparamsGP(gpi, dg[1,1], dg[1,2])
if(!(is.null(X) && is.null(y))) stop("must have gpi or X & y")
}

## check params and reconcile dab and gab
param <- match.arg(param)
dab <- gab <- c(0,0)
if(param == "d") dab <- ab
else gab <- ab

## make llik and dllik calculations
llik <- rep(NA, nrow(dg))
dllik <- data.frame(d=llik, d2=llik)
llik[1] <- llikGP(gpi, dab, gab)
dllik[1,] <- dllikGP(gpi, ab, param=param)
if(nrow(dg) > 1) {
for(i in 2:nrow(dg)) {
newparamsGP(gpi, dg[i,1], dg[i,2])
llik[i] <- llikGP(gpi, dab, gab)
dllik[i,] <- dllikGP(gpi, ab, param=param)
}
}

## clean up
if(is.null(gpi.in)) deleteGP(gpi)

## return
return(data.frame(llik, dllik))
}
```

## Try the laGP package in your browser

Any scripts or data that you put into this service are public.

laGP documentation built on March 31, 2023, 9:46 p.m.