Description Usage Arguments Details Value Note Author(s) References See Also Examples
Build a subdesign of X
of size end
, and infer parameters,
for approximate Gaussian process prediction at reference location(s)
Xref
. Return the moments of those predictive equations, and indices
into the local design
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24  laGP(Xref, start, end, X, Z, d = NULL, g = 1/10000,
method = c("alc", "alcopt", "alcray", "mspe", "nn", "fish"), Xi.ret = TRUE,
close = min((1000+end)*if(method[1] %in% c("alcray", "alcopt")) 10 else 1, nrow(X)),
alc.gpu = FALSE, numstart = if(method[1] == "alcray") ncol(X) else 1,
rect = NULL, lite = TRUE, verb = 0)
laGP.R(Xref, start, end, X, Z, d = NULL, g = 1/10000,
method = c("alc", "alcopt", "alcray", "mspe", "nn", "fish"),
Xi.ret = TRUE, pall = FALSE,
close = min((1000+end)*if(method[1] %in% c("alcray", "alcopt")) 10 else 1, nrow(X)),
parallel = c("none", "omp", "gpu"),
numstart = if(method[1] == "alcray") ncol(X) else 1,
rect = NULL, lite = TRUE, verb = 0)
laGPsep(Xref, start, end, X, Z, d = NULL, g = 1/10000,
method = c("alc", "alcopt", "alcray", "nn"), Xi.ret = TRUE,
close = min((1000+end)*if(method[1] %in% c("alcray", "alcopt")) 10 else 1, nrow(X)),
alc.gpu = FALSE, numstart = if(method[1] == "alcray") ncol(X) else 1,
rect = NULL, lite = TRUE, verb=0)
laGPsep.R(Xref, start, end, X, Z, d = NULL, g = 1/10000,
method = c("alc", "alcopt", "alcray", "nn"),
Xi.ret = TRUE, pall = FALSE,
close = min((1000+end)*if(method[1] %in% c("alcray", "alcopt")) 10 else 1, nrow(X)),
parallel = c("none", "omp", "gpu"),
numstart = if(method[1] == "alcray") ncol(X) else 1,
rect = NULL, lite = TRUE, verb = 0)

Xref 
a vector of length 
start 
the number of Nearest Neighbor (NN) locations for initialization; must
specify 
end 
the total size of the local designs; must have 
X 
a 
Z 
a vector of responses/dependent values with 
d 
a prior or initial setting for the lengthscale
parameter for a Gaussian correlation function; a (default)

g 
a prior or initial setting for the nugget parameter; a

method 
Specifies the method by which 
Xi.ret 
A scalar logical indicating whether or not a vector of indices
into 
pall 
a scalar logical (for 
close 
a nonnegative integer 
alc.gpu 
a scalar 
parallel 
a switch indicating if any parallel calculation of
the criteria is desired. Currently parallelization at this level is only
provided for option 
numstart 
a scalar integer indicating the number of rays for each
greedy search when 
rect 
an optional 
lite 
Similar to the 
verb 
a nonnegative integer specifying the verbosity level; 
A subdesign of X
of size end
is builtup according to
the criteria prescribed by the method
and then used to predict at
Xref
. The first start
locations are NNs in order to
initialize the first GP, via newGP
or newGPsep
,
and thereby initialize the
search. Each subsequent addition is found via the chosen criterion
(method
), and the GP fit is updated via updateGP
or updateGPsep
The runtime is cubic in end
, although
the multiplicative “constant” out front depends on the
method
chosen, the size of the design X
, and
close
. The "alcray"
method has a smaller constant
since it does not search over all candidates exhaustively.
After building the subdesign, local MLE/MAP lengthscale (and/or
nugget) parameters are estimated, depending on the d
and
g
arguments supplied. This is facilitated by calls to
mleGP
or jmleGP
.
Finally predGP
is called on the resulting local GP
model, and the parameters of the resulting Studentt distribution(s)
are returned. Unless Xi.ret = FALSE
, the indices of the
local design are also returned.
laGP.R
and laGPsep.R
are a prototype Ronly version for
debugging and transparency purposes. They are slower than
laGP
and laGPsep
, which are primarily in C, and may not
provide identical output in all cases due to differing library implementations
used as subroutines; see note below for an example. laGP.R
and other
.R
functions in the package may be useful for developing new programs
that involve similar subroutines. The current version of laGP.R
allows OpenMP and/or GPU parallelization of the criteria (method
) if
the package is compiled with the appropriate flags. See README/INSTALL in
the package source for more information. For algorithmic details, see
Gramacy, Niemi, & Weiss (2014)
The output is a list
with the following components.
mean 
a vector of predictive means of length 
s2 
a vector of Studentt scale parameters of length

df 
a Studentt degrees of freedom scalar (applies to all

llik 
a scalar indicating the maximized log likelihood or log posterior probability of the data/parameter(s) under the chosen subdesign; provided up to an additive constant 
time 
a scalar giving the passage of wallclock time elapsed for (substantive parts of) the calculation 
method 
a copy of the 
d 
a fulllist version of the 
g 
a fulllist version of the 
mle 
if 
Xi 
when 
close 
a copy of the input argument 
laGPsep
provides the same functionality as laGP
but deploys
a separable covariance function. However criteria (method
s) EFI and
MSPE are not supported. This is considered “beta” functionality
at this time.
Note that using method="NN"
gives the same result as specifying
start=end
, however at some extra computational expense.
Handling multiple reference locations
(nrow(Xref) > 1
) is “beta” functionality. In this case
the initial start
locations are chosen by applying NN to the
average distances to all Xref
locations. Using
method="alcopt"
causes exhaustive search to be approximated by
a continuous analog via closed form derivatives.
See alcoptGP
for more details. Although the approximation
provided has a spirit similar to method="alcray"
, in that
both methods are intended to offer a thrifty alternative,
method="alcray"
is not applicable when nrow(Xref) > 1
.
Differences between the C qsort
function and R's
order
function may cause chosen designs returned from
laGP
and laGP.R
(code and laGPsep
and laGPsep.R
)
to differ when multiple X
values are equidistant to Xref
Robert B. Gramacy [email protected] and Furong Sun [email protected]
F. Sun, R.B. Gramacy, B. Haaland, E. Lawrence, and A. Walker (2017) Emulating satellite drag from large simulation experiments; preprint on arXiv:1712.00182. http://arxiv.org/abs/1712.00182
R.B. Gramacy (2016). laGP: LargeScale Spatial Modeling via
Local Approximate Gaussian Processes in R., Journal of Statistical
Software, 72(1), 146; or see vignette("laGP")
R.B. Gramacy and B. Haaland (2016). Speeding up neighborhood search in local Gaussian process prediction. Technometrics, 58(3), pp. 294303; preprint on arXiv:1409.0074 http://arxiv.org/abs/1409.0074
R.B. Gramacy and D.W. Apley (2015). Local Gaussian process approximation for large computer experiments. Journal of Computational and Graphical Statistics, 24(2), pp. 561678; preprint on arXiv:1303.0383; http://arxiv.org/abs/1303.0383
R.B. Gramacy, J. Niemi, R.M. Weiss (2014). Massively parallel approximate Gaussian process regression. SIAM/ASA Journal on Uncertainty Quantification, 2(1), pp. 568584; preprint on arXiv:1310.5182; http://arxiv.org/abs/1310.5182
vignette("laGP")
,
aGP
, newGP
, updateGP
,
predGP
, mleGP
, jmleGP
,
alcGP
, mspeGP
, alcrayGP
,
randLine
## pathbased local prediction via laGP
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100  ## examining a particular laGP call from the larger analysis provided
## in the aGP documentation
## A simple 2d test function used in Gramacy & Apley (2014);
## thanks to Lee, Gramacy, Taddy, and others who have used it before
f2d < function(x, y=NULL)
{
if(is.null(y)) {
if(!is.matrix(x) && !is.data.frame(x)) x < matrix(x, ncol=2)
y < x[,2]; x < x[,1]
}
g < function(z)
return(exp((z1)^2) + exp(0.8*(z+1)^2)  0.05*sin(8*(z+0.1)))
z < g(x)*g(y)
}
## build up a design with N=~40K locations
x < seq(2, 2, by=0.02)
X < as.matrix(expand.grid(x, x))
Z < f2d(X)
## optional first pass of nearest neighbor
Xref < matrix(c(1.725, 1.725), nrow=TRUE)
out < laGP(Xref, 6, 50, X, Z, method="nn")
## second pass via ALC, ALCOPT, MSPE, and ALCray respectively,
## conditioned on MLE dvalues found above
out2 < laGP(Xref, 6, 50, X, Z, d=out$mle$d)
out2.alcopt < laGP(Xref, 6, 50, X, Z, d=out2$mle$d, method="alcopt")
out2.mspe < laGP(Xref, 6, 50, X, Z, d=out2$mle$d, method="mspe")
out2.alcray < laGP(Xref, 6, 50, X, Z, d=out2$mle$d, method="alcray")
## look at the different designs
plot(rbind(X[out2$Xi,], X[out2.mspe$Xi,]), type="n",
xlab="x1", ylab="x2", main="comparing local designs")
points(Xref[1], Xref[2], col=2, cex=0.5)
text(X[out2$Xi,], labels=1:50, cex=0.6)
text(X[out2.alcopt$Xi,], labels=1:50, cex=0.6, col="forestgreen")
text(X[out2.mspe$Xi,], labels=1:50, cex=0.6, col="blue")
text(X[out2.alcray$Xi,], labels=1:50, cex=0.6, col="red")
legend("right", c("ALC", "ALCOPT", "MSPE", "ALCRAY"),
text.col=c("black", "forestgreen", "blue", "red"), bty="n")
## compare computational time
c(nn=out$time, alc=out2$time, alcopt=out2.alcopt$time,
mspe=out2.mspe$time, alcray=out2.alcray$time)
## Not run:
## A comparison between ALCex, ALCopt and NN
## defining a predictive path
wx < seq(0.85, 0.45, length=100)
W < cbind(wx0.75, wx^3+0.51)
## three comparators from Sun, et al. (2017)
## largerthandefault "close" argument to capture locations nearby path
p.alc < laGPsep(W, 6, 100, X, Z, method="alc", close=10000, lite=FALSE)
p.alcopt < laGPsep(W, 6, 100, X, Z, method="alcopt", lite=FALSE)
## note that close=10*(1000+end) would be the default for method = "alcopt"
p.nn < laGPsep(W, 6, 100, X, Z, method="nn", close=10000, lite=FALSE)
## Mahalanobis distance to combine RMSE and covariance estimates
mahdist < function(Y, mu, Sigma) {
Ymmu < Y  mu
Sigmai < solve(Sigma)
return(sqrt(t(Ymmu) %*% Sigmai %*% Ymmu))
}
## comparison by Mahalanobis distance
WY < f2d(W)
c(alc=mahdist(WY, p.alc$mean, p.alc$Sigma),
alcopt=mahdist(WY, p.alcopt$mean, p.alcopt$Sigma),
nn=mahdist(WY, p.nn$mean, p.nn$Sigma))
## comparison via time
c(alc=p.alc$time, alcopt=p.alcopt$time, nn=p.nn$time)
## visualization
## first ALCex
plot(W, type="l", bty="n", lwd=2, xlab="x1", ylab="x2",
xlim=c(2.25,0), ylim=c(0.75,1.25), main="")
points(W[length(wx)/2,1], W[length(wx)/2,2], col=2, pch=19)
points(X[p.alc$Xi,], col="blue", cex=0.6)
legend("bottomright", c("x0", "ALCopt", "ALCex", "NN"),
pch=c(19, 23, 21, 22), bty="n", cex=1.25,
col=c("red", "purple", "blue", "forestgreen"))
legend("topleft", "W(x0)", lty=1, col=1, lwd=2, bty="n",
cex=1.25)
## then shifted NN
lines(W[,1]+0.25, W[,2]0.25)
points(W[length(wx)/2,1]+0.25, W[length(wx)/2,2]0.25, col=2, pch=19)
points(X[p.nn$Xi,1]+0.25, X[p.nn$Xi,2]0.25, pch=22,
col="forestgreen", cex=0.6)
## finally shifted (in the other direction) ALCopt
lines(W[,1]0.25, W[,2]+0.25)
points(W[length(wx)/2,1]0.25, W[length(wx)/2,2]+0.25, col=2, pch=19)
points(X[p.alcopt$Xi,1]0.25, X[p.alcopt$Xi,2]+0.25, pch=23,
col="purple", cex=0.6)
## End(Not run)

Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.