Description Usage Arguments Details Value Note Author(s) References See Also Examples
Facilitates localized Gaussian process inference and prediction at a large
set of predictive locations, by (essentially) calling laGP
at each location, and returning the moments of the predictive
equations, and indices into the design, thus obtained
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  aGP(X, Z, XX, start = 6, end = 50, d = NULL, g = 1/10000,
method = c("alc", "alcray", "mspe", "nn", "fish"), Xi.ret = TRUE,
close = min((1000+end)*if(method[1] == "alcray") 10 else 1, nrow(X)),
numrays = ncol(X), num.gpus = 0, gpu.threads = num.gpus,
omp.threads = if (num.gpus > 0) 0 else 1,
nn.gpu = if (num.gpus > 0) nrow(XX) else 0, verb = 1)
aGP.parallel(cls, XX, chunks = length(cls), X, Z, start = 6, end = 50,
d = NULL, g = 1/10000, method = c("alc", "alcray", "mspe", "nn", "fish"),
Xi.ret = TRUE,
close = min((1000+end)*if(method[1] == "alcray") 10 else 1, nrow(X)),
numrays = ncol(X), num.gpus = 0, gpu.threads = num.gpus,
omp.threads = if (num.gpus > 0) 0 else 1,
nn.gpu = if (num.gpus > 0) nrow(XX) else 0, verb = 1)
aGP.R(X, Z, XX, start = 6, end = 50, d = NULL, g = 1/10000,
method = c("alc", "alcray", "mspe", "nn", "fish"), Xi.ret = TRUE,
close = min((1000+end) *if(method[1] == "alcray") 10 else 1, nrow(X)),
numrays = ncol(X), laGP=laGP.R, verb = 1)
aGPsep(X, Z, XX, start = 6, end = 50, d = NULL, g = 1/10000,
method = c("alc", "alcray", "nn"), Xi.ret = TRUE,
close = min((1000+end)*if(method[1] == "alcray") 10 else 1, nrow(X)),
numrays = ncol(X), omp.threads = 1, verb = 1)
aGPsep.R(X, Z, XX, start = 6, end = 50, d = NULL, g = 1/10000,
method = c("alc", "alcray", "nn"), Xi.ret = TRUE,
close = min((1000+end)*if(method[1] == "alcray") 10 else 1, nrow(X)),
numrays = ncol(X), laGPsep=laGPsep.R, verb = 1)
aGP.seq(X, Z, XX, d, methods=rep("alc", 2), M=NULL, ncalib=0, ...)

X 
a 
Z 
a vector of responses/dependent values with 
XX 
a 
start 
the number of Nearest Neighbor (NN) locations to start each
independent call to 
end 
the total size of the local designs; 
d 
a prior or initial setting for the lengthscale
parameter in a Gaussian correlation function; a (default)

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

method 
specifies the method by which 
methods 
for 
Xi.ret 
a scalar logical indicating whether or not a 
close 
a nonnegative integer 
numrays 
a scalar integer indicating the number of rays for each
greedy search; only relevant when 
laGP 
applicable only to the Rversion 
laGPsep 
applicable only to the Rversion 
num.gpus 
applicable only to the Cversion 
gpu.threads 
applicable only to the Cversion 
omp.threads 
applicable only to the Cversion 
nn.gpu 
a scalar nonnegative integer between 
verb 
a nonnegative integer specifying the verbosity level; 
cls 
a cluster object created by 
chunks 
a scalar integer indicating the number of chunks to break 
M 
an optional function taking two matrix inputs, of 
ncalib 
an integer between 1 and 
... 
other arguments passed from 
This function invokes laGP
with argument Xref
= XX[i,]
for each i=1:nrow(XX)
, building up local designs,
inferring local correlation parameters, and
obtaining predictive locations independently for each location. For
more details see laGP
.
The function aGP.R
is a prototype Ronly version for
debugging and transparency purposes. It is slower than
aGP
, which is primarily in C. However it may be
useful for developing new programs that involve similar subroutines.
Note that aGP.R
may provide different output than aGP
due to differing library subroutines deployed in R and C.
The function aGP.parallel
allows aGP
to be called on segments
of the XX
matrix distributed to a cluster created by parallel.
It breaks XX
into chunks
which are sent to aGP
workers pointed to by the entries of cls
. The aGP.parallel
function
collects the outputs from each chunk before returning an object
almost identical to what would have been returned from a single aGP
call. On a single (SMP) node, this represents is a poorman's version of
the OpenMP version described below. On multiple nodes both can be used.
If compiled with OpenMP flags, the independent calls to
laGP
will be
farmed out to threads allowing them to proceed in parallel  obtaining
nearly linear speedups. At this time aGP.R
does not
facilitate parallel computation, although a future version may exploit
the parallel functionality for clustered parallel execution.
If num.gpus > 0
then the ALC part of the independent
calculations performed by each thread will be offloaded to a GPU.
If both gpu.threads >= 1
and omp.threads >= 1
,
some of the ALC calculations will be done on the GPUs, and some
on the CPUs. In our own experimentation we have not found this
to lead to large speedups relative to omp.threads = 0
when
using GPUs. For more details, see Gramacy, Niemi, & Weiss (2014).
The aGP.sep
function is provided primarily for use in calibration
exercises, see Gramacy, et al. (2015). It automates a sequence of
aGP
calls, each with a potentially different method,
successively feeding the previous estimate of local lengthscale (d
)
in as an initial set of values for the next call. It also allows the
use of aGP
to be bypassed, feeding the inputs into a usersupplied
M
function instead. This feature is enabled when
methods = FALSE
. The M
function takes two matrices
(same number of rows) as inputs, where the first ncol(X)  ncalib
columns represent
“field data” inputs shared by the physical and computer model
(in the calibration context), and the remaining ncalib
are
the extra tuning or calibration parameters required to evalue the computer
model. For examples illustrating aGP.seq
please see the
documentation file for discrep.est
and demo("calib")
The output is a list
with the following components.
mean 
a vector of predictive means of length 
var 
a vector of predictive variances of length

llik 
a vector indicating the log likelihood/posterior
probability of the data/parameter(s) under the chosen subdesign for
each predictive location in 
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 
The aGP.seq
function only returns the output from the final aGP
call.
When methods = FALSE
and M
is supplied, the returned object
is a data frame with a mean
column indicating the output of the computer
model run, and a var
column, which at this time is zero
aGPsep
provides the same functionality as aGP
but deploys
a separable covariance function. Criteria (method
s) EFI and
MSPE are not supported by aGPsep
at this time.
Note that using method="NN"
gives the same result as specifying
start=end
, however at some extra computational expense.
At this time, this function provides no facility to find local designs
for the subset of predictive locations XX
jointly, i.e.,
providing a matrix Xref
to laGP
. See laGP
for more details/support for this alternative.
The use of OpenMP threads with aGPsep
is not as efficient as with
aGP
when calculating MLEs with respect to the lengthscale (i.e.,
d=list(mle=TRUE, ...)
). The reason is that the lbfgsb
C
entry point uses static variables, and is therefore not thread safe.
To circumvent this problem, an OpenMP critical
pragma is used,
which can create a small bottle neck
Robert B. Gramacy rbg@vt.edu
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 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
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
vignette("laGP")
,
laGP
, alcGP
, mspeGP
, alcrayGP
,
makeCluster
, clusterApply
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  ## first, a "computer experiment"
## 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 < expand.grid(x, x)
Z < f2d(X)
## predictive grid with NN=400 locations,
## change NN to 10K (length=100) to mimic setup in Gramacy & Apley (2014)
## the low NN set here is for fast CRAN checks
xx < seq(1.975, 1.975, length=10)
XX < expand.grid(xx, xx)
ZZ < f2d(XX)
## get the predictive equations, first based on Nearest Neighbor
out < aGP(X, Z, XX, method="nn", verb=0)
## RMSE
sqrt(mean((out$mean  ZZ)^2))
## Not run:
## refine with ALC
out2 < aGP(X, Z, XX, method="alc", d=out$mle$d)
## RMSE
sqrt(mean((out2$mean  ZZ)^2))
## visualize the results
par(mfrow=c(1,3))
image(xx, xx, matrix(out2$mean, nrow=length(xx)), col=heat.colors(128),
xlab="x1", ylab="x2", main="predictive mean")
image(xx, xx, matrix(out2$meanZZ, nrow=length(xx)), col=heat.colors(128),
xlab="x1", ylab="x2", main="bias")
image(xx, xx, matrix(sqrt(out2$var), nrow=length(xx)), col=heat.colors(128),
xlab="x1", ylab="x2", main="sd")
## refine with MSPE
out3 < aGP(X, Z, XX, method="mspe", d=out2$mle$d)
## RMSE
sqrt(mean((out3$mean  ZZ)^2))
## End(Not run)
## version with ALCray which is much faster than the ones not
## run above
out2r < aGP(X, Z, XX, method="alcray", d=out$mle$d, verb=0)
sqrt(mean((out2r$mean  ZZ)^2))
## a simple example with estimated nugget
library(MASS)
## motorcycle data and predictive locations
X < matrix(mcycle[,1], ncol=1)
Z < mcycle[,2]
XX < matrix(seq(min(X), max(X), length=100), ncol=1)
## first stage
out < aGP(X=X, Z=Z, XX=XX, end=30, g=list(mle=TRUE), verb=0)
## plot smoothed versions of the estimated parameters
par(mfrow=c(2,1))
df < data.frame(y=log(out$mle$d), XX)
lo < loess(y~., data=df, span=0.25)
plot(XX, log(out$mle$d), type="l")
lines(XX, lo$fitted, col=2)
dfnug < data.frame(y=log(out$mle$g), XX)
lonug < loess(y~., data=dfnug, span=0.25)
plot(XX, log(out$mle$g), type="l")
lines(XX, lonug$fitted, col=2)
## second stage design
out2 < aGP(X=X, Z=Z, XX=XX, end=30, verb=0,
d=list(start=exp(lo$fitted), mle=FALSE),
g=list(start=exp(lonug$fitted)))
## plot the estimated surface
par(mfrow=c(1,1))
plot(X,Z)
df < 20
s2 < out2$var*(df2)/df
q1 < qt(0.05, df)*sqrt(s2) + out2$mean
q2 < qt(0.95, df)*sqrt(s2) + out2$mean
lines(XX, out2$mean)
lines(XX, q1, col=1, lty=2)
lines(XX, q2, col=1, lty=2)
## compare to the singleGP result provided in the mleGP documentation

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