Type: | Package |
Title: | Tools for Gaussian Process Modeling in Uncertainty Quantification |
Version: | 0.1.0-6 |
Date: | 2024-04-24 |
Maintainer: | Pulong Ma <mpulong@gmail.com> |
Description: | Gaussian processes ('GPs') have been widely used to model spatial data, 'spatio'-temporal data, and computer experiments in diverse areas of statistics including spatial statistics, 'spatio'-temporal statistics, uncertainty quantification, and machine learning. This package creates basic tools for fitting and prediction based on 'GPs' with spatial data, 'spatio'-temporal data, and computer experiments. Key characteristics for this GP tool include: (1) the comprehensive implementation of various covariance functions including the 'Matérn' family and the Confluent 'Hypergeometric' family with isotropic form, tensor form, and automatic relevance determination form, where the isotropic form is widely used in spatial statistics, the tensor form is widely used in design and analysis of computer experiments and uncertainty quantification, and the automatic relevance determination form is widely used in machine learning; (2) implementations via Markov chain Monte Carlo ('MCMC') algorithms and optimization algorithms for GP models with all the implemented covariance functions. The methods for fitting and prediction are mainly implemented in a Bayesian framework; (3) model evaluation via Fisher information and predictive metrics such as predictive scores; (4) built-in functionality for simulating 'GPs' with all the implemented covariance functions; (5) unified implementation to allow easy specification of various 'GPs'. |
License: | GPL-2 | GPL-3 [expanded from: GPL (≥ 2)] |
Encoding: | UTF-8 |
BugReports: | https://github.com/pulongma/GPBayes/issues |
Imports: | Rcpp (≥ 1.0.1), stats, methods |
LinkingTo: | Rcpp, RcppEigen, RcppProgress |
SystemRequirements: | GNU Scientific Library version >= 2.5 |
NeedsCompilation: | yes |
Packaged: | 2024-04-25 00:31:01 UTC; mapulong |
Author: | Pulong Ma [aut, cre] |
RoxygenNote: | 7.3.1 |
Repository: | CRAN |
Date/Publication: | 2024-04-25 02:40:03 UTC |
Tools for Gaussian Stochastic Process Modeling in Uncertainty Quantification
Description
Gaussian processes (GPs) have been widely used to model spatial data, spatio-temporal data, and computer experiments in diverse areas of statistics including spatial statistics, spatio-temporal statistics, uncertainty quantification, and machine learning. This package creates basic tools for fitting and prediction based on GPs with spatial data, spatio-temporal data, and computer experiments. Key characteristics for this GP tool include: (1) the comprehensive implementation of various covariance functions including the Matérn family and the Confluent Hypergeometric family with isotropic form, tensor form, and automatic relevance determination form, where the isotropic form is widely used in spatial statistics, the tensor form is widely used in design and analysis of computer experiments and uncertainty quantification, and the automatic relevance determination form is widely used in machine learning; (2) implementations via Markov chain Monte Carlo (MCMC) algorithms and optimization algorithms for GP models with all the implemented covariance functions. The methods for fitting and prediction are mainly implemented in a Bayesian framework; (3) model evaluation via Fisher information and predictive metrics such as predictive scores; (4) built-in functionality for simulating GPs with all the implemented covariance functions; (5) unified implementation to allow easy specification of various GPs.
Details
Data types: For many scientific applications, spatial data, spatio-temporal data, and computer experiments arise naturally. This package provides a comprehensive set of basic tools to fit GaSP models for univariate and multivariate spatial data, spatio-temporal data, computer experiments. Various covariance functions have been implemented including the Confluent Hypergeometric covariance functions, the Matérn covariance functions, the Gaussian covariance function, the generalized Cauchy covariance function. These covariance families can be in isotropic form, in tensor form, or in automatic relevance determination form. The routines
kernel
andikernel
contain the details of implementation.Model simulation: This package can simulate realizations from GaSP for different types of data including spatial data, spatio-temporal data, and computer experiments. This feature is quite useful in part because benchmarks are used to evaluate the performance of GaSP models. This functionality is implemented in the routine
gp.sim
for unconditional simulation andgp.condsim
for conditional simulation.Model fitting: Both maximum likelihood methods (or its variants) and Bayes estimation methods such as maximum a posterior (MAP) and Markov chain Monte Carlo (MCMC) methods are implemented. In this package, the nugget parameter is included in the model by default for the sake of better prediction performance and stable computation in practice. In addition, the smoothness parameter in covariance functions such as the Matérn class and the Confluent Hypergeometric class can be estimated. The routine
gp.optim
provides optimization based estimation approaches and the routinegp.mcmc
provides MCMC algorithms based estimation approaches.Model prediction: Prediction is made based on the parameter estimation procedure. If maximum likelihood estimation (MLE) methods are used for parameter estimation, the plug-in approach is used for prediction in the sense that MLEs of parameters are plugged into posterior predictive distributions. If partial Bayes methods (e.g., maximum a posterior) are used, the plug-in approach is used for prediction as well. If fully Bayes methods via MCMC algorithms are used, posterior samples are drawn from posterior predictive distributions. The routine
gp.mcmc
allows prediction to be made within the MCMC algorithms, while the routinegp.predict
generates prediction with estimated parameters.Model assessment: Tools for assessing model adequacy are included in a Bayesian context. Deviance information criteria (DIC), log pointwise predictive density, and log joint predictive density can be computed via the routine
gp.model.adequacy
.
Author(s)
Pulong Ma mpulong@gmail.com
References
Cressie, N. (1993). “Statistics for Spatial Data.” John Wiley & Sons, New York, revised edition.
Ma and Bhadra (2023). “Beyond Matérn: On a Class of Interpretable Confluent Hypergeometric Covariance Functions.” Journal of the American Statistical Association 118(543), 2045-2058.
Sacks, Jerome, William J Welch, Toby J Mitchell, and Henry P Wynn. (1989). “Design and Analysis of Computer Experiments.” Statistical Science 4(4). Institute of Mathematical Statistics: 409–435.
Santner, Thomas J., Brian J. Williams, and William I. Notz. (2018). “The Design and Analysis of Computer Experiments”; 2nd Ed. New York: Springer.
Stein, Michael L. (1999). “Interpolation of Spatial Data.” Springer Science & Business Media, New York.
See Also
Examples
#####################################################################
#####################################################################
############## Examples for fitting univariate GP models ############
## Set up the Sine example from the tgp package
code = function(x){
y = (sin(pi*x/5) + 0.2*cos(4*pi*x/5))*(x<=9.6) + (x/10-1)*(x>9.6)
}
n=100
input = seq(0, 20, length=n)
XX = seq(0, 20, length=99)
Ztrue = code(input)
set.seed(1234)
output = Ztrue + rnorm(length(Ztrue), sd=0.1)
df.data = data.frame(x=c(input), y=output, y.true=Ztrue)
## fitting a GaSP model with the Cauchy prior
fit = GaSP(formula=~1, output, input,
param=list(range=3, nugget=0.1, nu=2.5),
smooth.est=FALSE, input.new=XX,
cov.model=list(family="matern", form="isotropic"),
proposal=list(range=.35, nugget=.8, nu=0.8),
dtype="Euclidean", model.fit="Cauchy_prior", nsample=3000,
burnin=500, verbose=TRUE)
## fitting a GaSP model with the beta prior
fit = GaSP(formula=~1, output, input,
param=list(range=3, nugget=0.1, nu=2.5),
smooth.est=FALSE, input.new=XX,
cov.model=list(family="matern", form="isotropic"),
prior=list(range=list(a=1,b=1,lb=0,ub=20),
nugget=list(a=1,b=1,lb=0,ub=var(output)),
proposal=list(range=.35, nugget=.8, nu=0.8),
dtype="Euclidean", model.fit="Beta_prior", nsample=3000,
burnin=500, verbose=TRUE))
## fitting a GaSP model with the marginal maximum likelihood approach
fit = GaSP(formula=~1, output, input,
param=list(range=3, nugget=0.1, nu=2.5),
smooth.est=FALSE, input.new=XX,
cov.model=list(family="matern", form="isotropic"),
dtype="Euclidean", model.fit="MMLE", verbose=TRUE)
## fitting a GaSP model with the profile maximum likelihood approach
fit = GaSP(formula=~1, output, input,
param=list(range=3, nugget=0.1, nu=2.5),
smooth.est=FALSE, input.new=XX,
cov.model=list(family="matern", form="isotropic"),
dtype="Euclidean", model.fit="MPLE", verbose=TRUE)
Modified Bessel function of the second kind
Description
This function calls the GSL scientific library to evaluate the modified Bessel function of the second kind.
Usage
BesselK(nu, z)
Arguments
nu |
a real positive value |
z |
a real positive value |
Value
a numerical value
Author(s)
Pulong Ma mpulong@gmail.com
See Also
The Confluent Hypergeometric correlation function proposed by Ma and Bhadra (2023)
Description
This function computes the Confluent Hypergeometric correlation function given a distance matrix. The Confluent Hypergeometric correlation function is given by
C(h) = \frac{\Gamma(\nu+\alpha)}{\Gamma(\nu)}
\mathcal{U}\left(\alpha, 1-\nu, \biggr(\frac{h}{\beta}\biggr)^2 \right),
where \alpha
is the tail decay parameter. \beta
is the range parameter.
\nu
is the smoothness parameter. \mathcal{U}(\cdot)
is the confluent hypergeometric
function of the second kind. Note that this parameterization of the CH covariance
is different from the one in Ma and Bhadra (2023). For details about this covariance,
see Ma and Bhadra (2023; doi:10.1080/01621459.2022.2027775).
Usage
CH(d, range, tail, nu)
Arguments
d |
a matrix of distances |
range |
a numerical value containing the range parameter |
tail |
a numerical value containing the tail decay parameter |
nu |
a numerical value containing the smoothness parameter |
Value
a numerical matrix
Author(s)
Pulong Ma mpulong@gmail.com
See Also
GPBayes-package, GaSP
, gp, matern
, kernel
, ikernel
Building, fitting, predicting for a GaSP model
Description
This function serves as a wrapper to build, fit, and make prediction
for a Gaussian process model. It calls on functions gp
, gp.mcmc
,
gp.optim
, gp.predict
.
Usage
GaSP(
formula = ~1,
output,
input,
param,
smooth.est = FALSE,
input.new = NULL,
cov.model = list(family = "CH", form = "isotropic"),
model.fit = "Cauchy_prior",
prior = list(),
proposal = list(range = 0.35, tail = 2, nugget = 0.8, nu = 0.8),
nsample = 5000,
burnin = 1000,
opt = NULL,
bound = NULL,
dtype = "Euclidean",
verbose = TRUE
)
Arguments
formula |
an object of |
output |
a numerical vector including observations or outputs in a GaSP |
input |
a matrix including inputs in a GaSP |
param |
a list including values for regression parameters, covariance parameters, and nugget variance parameter. The specification of param should depend on the covariance model.
|
smooth.est |
a logical value indicating whether smoothness parameter will be estimated. |
input.new |
a matrix of new input locations |
cov.model |
a list of two strings: family, form, where family indicates the family of covariance functions including the Confluent Hypergeometric class, the Matérn class, the Cauchy class, the powered-exponential class. form indicates the specific form of covariance structures including the isotropic form, tensor form, automatic relevance determination form.
|
model.fit |
a string indicating the choice of priors on correlation parameters:
|
prior |
a list containing tuning parameters in prior distribution. This is used only if a subjective Bayes estimation method with informative priors is used. |
proposal |
a list containing tuning parameters in proposal distribution. This is used only if a Bayes estimation method is used. |
nsample |
an integer indicating the number of MCMC samples. |
burnin |
an integer indicating the burn-in period. |
opt |
a list of arguments to setup the
|
bound |
Default value is
for the Confluent Hypergeometric covariance and the Cauchy covariance.
|
dtype |
a string indicating the type of distance:
|
verbose |
a logical value. If it is |
Value
a list containing the S4
object gp and prediction results
Author(s)
Pulong Ma mpulong@gmail.com
See Also
GPBayes-package
, gp
, gp.mcmc
, gp.optim
, gp.predict
Examples
code = function(x){
y = (sin(pi*x/5) + 0.2*cos(4*pi*x/5))*(x<=9.6) + (x/10-1)*(x>9.6)
return(y)
}
n=100
input = seq(0, 20, length=n)
XX = seq(0, 20, length=99)
Ztrue = code(input)
set.seed(1234)
output = Ztrue + rnorm(length(Ztrue), sd=0.1)
# fitting a GaSP model with the objective Bayes approach
fit = GaSP(formula=~1, output, input,
param=list(range=3, nugget=0.1, nu=2.5),
smooth.est=FALSE, input.new=XX,
cov.model=list(family="matern", form="isotropic"),
proposal=list(range=.35, nugget=.8, nu=0.8),
dtype="Euclidean", model.fit="Cauchy_prior", nsample=50,
burnin=10, verbose=TRUE)
Confluent hypergeometric function of the second kind
Description
This function calls the GSL scientific library to evaluate the confluent hypergeometric function of the second kind; see Abramowitz and Stegun 1972, p.505.
Usage
HypergU(a, b, x)
Arguments
a |
a real value |
b |
a real value |
x |
a real value |
Value
a numerical value
Author(s)
Pulong Ma mpulong@gmail.com
See Also
The generalized Cauchy correlation function
Description
This function computes the generalized Cauchy correlation function given a distance matrix. The generalized Cauchy covariance is given by
C(h) = \left\{ 1 + \left( \frac{h}{\phi} \right)^{\nu}
\right\}^{-\alpha/\nu},
where \phi
is the range parameter. \alpha
is the tail decay parameter.
\nu
is the smoothness parameter.
The case where \nu=2
corresponds to the Cauchy covariance model, which is infinitely differentiable.
Usage
cauchy(d, range, tail, nu)
Arguments
d |
a matrix of distances |
range |
a numerical value containing the range parameter |
tail |
a numerical value containing the tail decay parameter |
nu |
a numerical value containing the smoothness parameter |
Value
a numerical matrix
Author(s)
Pulong Ma mpulong@gmail.com
See Also
Find the correlation parameter given effective range
Description
This function finds the correlation parameter given effective range
Usage
cor.to.par(
d,
param,
family = "CH",
cor.target = 0.05,
lower = NULL,
upper = NULL,
tol = .Machine$double.eps
)
Arguments
d |
a numerical value containing the effective range |
param |
a list containing correlation parameters. The specification of
param should depend on the covariance model. If the parameter value is
|
family |
a string indicating the type of covariance structure. The following correlation functions are implemented:
|
cor.target |
a numerical value. The default value is 0.05, which means that correlation parameters are searched such that the correlation is approximately 0.05. |
lower |
a numerical value. This sets the lower bound to find the
correlation parameter via the |
upper |
a numerical value. This sets the upper bound to find the
correlation parameter via the |
tol |
a numerical value. This sets the precision of the solution with default value
specified as the machine precision |
Value
a numerical value of correlation parameters
Author(s)
Pulong Ma mpulong@gmail.com
See Also
GPBayes-package, GaSP
, kernel
, ikernel
Examples
range = cor.to.par(1,param=list(tail=0.5,nu=2.5), family="CH")
tail = cor.to.par(1,param=list(range=0.5,nu=2.5), family="CH")
range = cor.to.par(1,param=list(nu=2.5),family="matern")
A wraper to construct the derivative of correlation matrix with respect to correlation parameters
Description
This function wraps existing built-in routines to construct the derivative of correlation matrix with respect to correlation parameters.
Usage
deriv_kernel(d, range, tail, nu, covmodel)
Arguments
d |
a matrix or a list of distances returned from |
range |
a vector of range parameters |
tail |
a vector of tail decay parameters |
nu |
a vector of smoothness parameters |
covmodel |
a list of two strings: family, form, where family indicates the family of covariance functions including the Confluent Hypergeometric class, the Matérn class, the Cauchy class, the powered-exponential class. form indicates the specific form of covariance structures including the isotropic form, tensor form, automatic relevance determination form.
|
Value
a list of matrices
Author(s)
Pulong Ma mpulong@gmail.com
See Also
CH
, matern
, kernel
, GPBayes-package, GaSP
Examples
input = seq(0,1,length=10)
d = distance(input,input,type="isotropic",dtype="Euclidean")
dR = deriv_kernel(d,range=0.5,tail=0.2,nu=2.5,
covmodel=list(family="CH",form="isotropic"))
Compute distances for two sets of inputs
Description
This function computes distances for two sets of inputs and returns
a R
object.
Usage
distance(input1, input2, type = "isotropic", dtype = "Euclidean")
Arguments
input1 |
a matrix of inputs |
input2 |
a matrix of inputs |
type |
a string indicating the form of distances with three froms supported currently: isotropic, tensor, ARD. |
dtype |
a string indicating distance type: Euclidean, GCD, where the latter indicates great circle distance. |
Value
a R object holding distances for two sets of inputs. If type is isotropic, a matrix of distances is returned; if type is tensor or ARD, a list of distance matrices along each input dimension is returned.
a numeric vector or matrix of distances
Author(s)
Pulong Ma mpulong@gmail.com
Examples
input = seq(0,1,length=20)
d = distance(input, input, type="isotropic", dtype="Euclidean")
Construct the S4
object gp
Description
This function constructs the S4
object gp that is used for Gaussian process
model fitting and prediction.
Usage
gp(
formula = ~1,
output,
input,
param,
smooth.est = FALSE,
cov.model = list(family = "CH", form = "isotropic"),
dtype = "Euclidean"
)
Arguments
formula |
an object of |
output |
a numerical vector including observations or outputs in a GaSP |
input |
a matrix including inputs in a GaSP |
param |
a list including values for regression parameters, covariance parameters, and nugget variance parameter. The specification of param should depend on the covariance model.
|
smooth.est |
a logical value indicating whether smoothness parameter will be estimated. |
cov.model |
a list of two strings: family, form, where family indicates the family of covariance functions including the Confluent Hypergeometric class, the Matérn class, the Cauchy class, the powered-exponential class. form indicates the specific form of covariance structures including the isotropic form, tensor form, automatic relevance determination form.
|
dtype |
a string indicating the type of distance:
|
Value
an S4
object of gp class
Author(s)
Pulong Ma mpulong@gmail.com
See Also
Examples
code = function(x){
y = (sin(pi*x/5) + 0.2*cos(4*pi*x/5))*(x<=9.6) + (x/10-1)*(x>9.6)
return(y)
}
n=100
input = seq(0, 20, length=n)
XX = seq(0, 20, length=99)
Ztrue = code(input)
set.seed(1234)
output = Ztrue + rnorm(length(Ztrue), sd=0.1)
obj = gp(formula=~1, output, input,
param=list(range=4, nugget=0.1,nu=2.5),
smooth.est=FALSE,
cov.model=list(family="matern", form="isotropic"))
The gp
class
Description
This is an S4 class definition for gp
in the GaSP
package.
Slots
formula
an object of
formula
class that specifies regressors; seeformula
for details.output
a numerical vector including observations or outputs in a GaSP
input
a matrix including inputs in a GaSP
param
a list including values for regression parameters, correlation parameters, and nugget variance parameter. The specification of param should depend on the covariance model.
The regression parameters are denoted by coeff. Default value is
\mathbf{0}
.The marginal variance or partial sill is denoted by sig2. Default value is 1.
The nugget variance parameter is denoted by nugget for all covariance models. Default value is 0.
For the Confluent Hypergeometric class, range is used to denote the range parameter
\beta
. tail is used to denote the tail decay parameter\alpha
. nu is used to denote the smoothness parameter\nu
.For the generalized Cauchy class, range is used to denote the range parameter
\phi
. tail is used to denote the tail decay parameter\alpha
. nu is used to denote the smoothness parameter\nu
.For the Matérn class, range is used to denote the range parameter
\phi
. nu is used to denote the smoothness parameter\nu
. When\nu=0.5
, the Matérn class corresponds to the exponential covariance.For the powered-exponential class, range is used to denote the range parameter
\phi
. nu is used to denote the smoothness parameter. When\nu=2
, the powered-exponential class corresponds to the Gaussian covariance.
cov.model
a list of two strings: family, form, where family indicates the family of covariance functions including the Confluent Hypergeometric class, the Matérn class, the Cauchy class, the powered-exponential class. form indicates the specific form of covariance structures including the isotropic form, tensor form, automatic relevance determination form.
- family
-
- CH
The Confluent Hypergeometric correlation function is given by
C(h) = \frac{\Gamma(\nu+\alpha)}{\Gamma(\nu)} \mathcal{U}\left(\alpha, 1-\nu, \left(\frac{h}{\beta}\right)^2\right),
where
\alpha
is the tail decay parameter.\beta
is the range parameter.\nu
is the smoothness parameter.\mathcal{U}(\cdot)
is the confluent hypergeometric function of the second kind. Note that this parameterization of the CH covariance is different from the one in Ma and Bhadra (2023). For details about this covariance, see Ma and Bhadra (2023; doi:10.1080/01621459.2022.2027775).- cauchy
The generalized Cauchy covariance is given by
C(h) = \left\{ 1 + \left( \frac{h}{\phi} \right)^{\nu} \right\}^{-\alpha/\nu},
where
\phi
is the range parameter.\alpha
is the tail decay parameter.\nu
is the smoothness parameter with default value at 2.- matern
The Matérn correlation function is given by
C(h)=\frac{2^{1-\nu}}{\Gamma(\nu)} \left( \frac{h}{\phi} \right)^{\nu} \mathcal{K}_{\nu}\left( \frac{h}{\phi} \right),
where
\phi
is the range parameter.\nu
is the smoothness parameter.\mathcal{K}_{\nu}(\cdot)
is the modified Bessel function of the second kind of order\nu
.- exp
The exponential correlation function is given by
C(h)=\exp(-h/\phi),
where
\phi
is the range parameter. This is the Matérn correlation with\nu=0.5
.- matern_3_2
The Matérn correlation with
\nu=1.5
.- matern_5_2
The Matérn correlation with
\nu=2.5
.- powexp
The powered-exponential correlation function is given by
C(h)=\exp\left\{-\left(\frac{h}{\phi}\right)^{\nu}\right\},
where
\phi
is the range parameter.\nu
is the smoothness parameter.- gauss
The Gaussian correlation function is given by
C(h)=\exp\left(-\frac{h^2}{\phi^2}\right),
where
\phi
is the range parameter.
- form
-
- isotropic
This indicates the isotropic form of covariance functions. That is,
C(\mathbf{h}) = C^0(\|\mathbf{h}\|; \boldsymbol \theta),
where
\| \mathbf{h}\|
denotes the Euclidean distance or the great circle distance for data on sphere.C^0(\cdot)
denotes any isotropic covariance family specified in family.- tensor
This indicates the tensor product of correlation functions. That is,
C(\mathbf{h}) = \prod_{i=1}^d C^0(|h_i|; \boldsymbol \theta_i),
where
d
is the dimension of input space.h_i
is the distance along thei
th input dimension. This type of covariance structure has been often used in Gaussian process emulation for computer experiments.- ARD
This indicates the automatic relevance determination form. That is,
C(\mathbf{h}) = C^0\left(\sqrt{\sum_{i=1}^d\frac{h_i^2}{\phi^2_i}}; \boldsymbol \theta \right),
where
\phi_i
denotes the range parameter along thei
th input dimension.
smooth.est
a logical value. If it is
TRUE
, the smoothness parameter will be estimated; otherwise the smoothness is not estimated.dtype
a string indicating the type of distance:
- Euclidean
Euclidean distance is used. This is the default choice.
- GCD
Great circle distance is used for data on sphere.
loglik
a numerical value containing the log-likelihood with current
gp
object.mcmc
a list containing MCMC samples if available.
prior
a list containing tuning parameters in prior distribution. This is used only if a Bayes estimation method with informative priors is used.
proposal
a list containing tuning parameters in proposal distribution. This is used only if a Bayes estimation method is used.
info
a list containing the maximum distance in the input space. It should be a vector if isotropic covariance is used, otherwise it is vector of maximum distances along each input dimension
Author(s)
Pulong Ma mpulong@gmail.com
See Also
Perform conditional simulation from a Gaussian process
Description
This function simulate from the Gaussian process model conditional on existing data (i.e., locations, observations). This is known as conditional simulation in geostatistics, which simulates realizations from the (posterior) predictive distribution of the process given the data.
Usage
gp.condsim(obj, XX, nsample = 1, seed = NULL)
Arguments
obj |
an |
XX |
a matrix of new locations where conditional simulation is performed. |
nsample |
number of conditional simulations default at 1 |
seed |
random generation seed default to be |
Value
an array (vector or matrix) of conditional simulations
Fisher information matrix
Description
This function computes the Fisher information matrix I(\sigma^2, \boldsymbol \theta)
for a
Gaussian process model y(\cdot) \sim \mathcal{GP}(h^\top(\mathbf{x})\mathbf{b}, \sigma^2 c(\cdot, \cdot) )
, where
c(\mathbf{x}_1, \mathbf{x}_2) = r(\mathbf{x}_1, \mathbf{x}_2) + \tau^2\mathbf{1}(\mathbf{x}_1=\mathbf{x}_2)
with correlation function
r(\cdot, \cdot)
and nugget parameter \tau^2
; \mathbf{b}
is a vector of regression coefficients,
\sigma^2
is the variance parameter (or partial sill).
Given n
data points that are assumed to be realizations from the GP model,
the standard likelihood is defined as
L(\mathbf{b}, \sigma^2, \boldsymbol \theta; \mathbf{y}) = \mathcal{N}_n(\mathbf{H}\mathbf{b}, \sigma^2 (\mathbf{R} + \tau^2\mathbf{I}) ),
where \mathbf{y}:=(y(\mathbf{x}_1), \ldots, y(\mathbf{x}_n))^\top
is a vector of n
observations.
\mathbf{H}
is a matrix of covariates, \boldsymbol \theta
contains correlation
parameters and nugget parameter, \mathbf{R}
denotes the n
-by-n
correlation matrix.
The integrated likelihood is defined as
L^{I}(\sigma^2, \boldsymbol \theta; \mathbf{y}) = \int L(\mathbf{b}, \sigma^2, \boldsymbol \theta; \mathbf{y}) \pi^{R}(\mathbf{b} \mid \sigma^2, \boldsymbol \theta) d \mathbf{b},
where \pi^{R}(\mathbf{b} \mid \sigma^2, \boldsymbol \theta)=1
is the conditional Jeffreys-rule (or reference prior)
in the model with the above standard likelihood when (\sigma^2, \boldsymbol \theta)
is assumed to be known.
For the Matérn class, current implementation only computes Fisher information matrix for variance parameter
\sigma^2
, range parameter\phi
, and nugget variance parameter\tau^2
. That is,I(\sigma^2, \boldsymbol \theta) = I(\sigma^2, \phi, \tau^2)
.For the Confluent Hypergeometric class, current implementation computes Fisher information matrix for variance parameter
\sigma^2
, range parameter\beta
, tail decay parameter\alpha
, smoothness parameter\nu
and nugget variance parameter\tau^2
. That is,I(\sigma^2, \boldsymbol \theta) = I(\sigma^2, \beta, \alpha, \nu, \tau^2)
.
Usage
gp.fisher(
obj = NULL,
intloglik = FALSE,
formula = ~1,
input = NULL,
param = NULL,
cov.model = NULL,
dtype = "Euclidean"
)
Arguments
obj |
a |
intloglik |
a logical value with default value |
formula |
an object of |
input |
a matrix including inputs in a GaSP |
param |
a list including values for regression parameters, covariance parameters, and nugget variance parameter. The specification of param should depend on the covariance model.
|
cov.model |
a list of two strings: family, form, where family indicates the family of covariance functions including the Confluent Hypergeometric class, the Matérn class, the Cauchy class, the powered-exponential class. form indicates the specific form of covariance structures including the isotropic form, tensor form, automatic relevance determination form.
|
dtype |
a string indicating the type of distance:
|
Value
a numerical matrix of Fisher information
Author(s)
Pulong Ma mpulong@gmail.com
See Also
GPBayes-package, GaSP
, gp
, kernel
, ikernel
,
Examples
n=100
input = seq(0, 20, length=n)
range = 1
tail = .5
nu = 1.5
sig2 = 1
nugget = 0.01
coeff = 0
par = list(range=range, tail=tail, nu=nu, sig2=sig2, nugget=nugget, coeff=coeff)
I = gp.fisher(formula=~1, input=input,
param=list(range=4, nugget=0.1,nu=2.5),
cov.model=list(family="CH", form="isotropic"))
get posterior summary for MCMC samples
Description
This function processes posterior samples in the gp
object.
Usage
gp.get.mcmc(obj, burnin = 500)
Arguments
obj |
a |
burnin |
a numerical value specifying the burn-in period for calculating posterior summaries. |
Value
a list of posterior summaries
See Also
GPBayes-package, GaSP
, gp
, gp.mcmc
A wraper to fit a Gaussian stochastic process model with MCMC algorithms
Description
This function is a wraper to estimate parameters via MCMC algorithms in the GaSP model with different choices of priors.
Usage
gp.mcmc(
obj,
input.new = NULL,
method = "Cauchy_prior",
prior = list(),
proposal = list(),
nsample = 10000,
verbose = TRUE
)
Arguments
obj |
an |
input.new |
a matrix of prediction locations. Default value is |
method |
a string indicating the Bayes estimation approaches with different choices of priors on correlation parameters:
|
prior |
a list containing tuning parameters in prior distributions. This is used only if a Bayes estimation method with subjective priors is used. |
proposal |
a list containing tuning parameters in proposal distributions. This is used only if a Bayes estimation method is used. |
nsample |
an integer indicating the number of MCMC samples. |
verbose |
a logical value. If it is |
Value
a gp
object with prior, proposal, MCMC samples included.
Author(s)
Pulong Ma mpulong@gmail.com
See Also
GPBayes-package, GaSP
, gp, gp.optim
Examples
code = function(x){
y = (sin(pi*x/5) + 0.2*cos(4*pi*x/5))*(x<=9.6) + (x/10-1)*(x>9.6)
return(y)
}
n=100
input = seq(0, 20, length=n)
XX = seq(0, 20, length=99)
Ztrue = code(input)
set.seed(1234)
output = Ztrue + rnorm(length(Ztrue), sd=0.1)
obj = gp(formula=~1, output, input,
param=list(range=4, nugget=0.1,nu=2.5),
smooth.est=FALSE,
cov.model=list(family="matern", form="isotropic"))
fit.mcmc = gp.mcmc(obj, method="Cauchy_prior",
proposal=list(range=0.3, nugget=0.8),
nsample=100, verbose=TRUE)
Model assessment based on Deviance information criterion (DIC), logarithmic pointwise predictive density (lppd), and logarithmic joint predictive density (ljpd).
Description
This function computes effective number of parameters (pD), deviance information criterion (DIC), logarithmic pointwise predictive density (lppd), and logarithmic joint predictive density (ljpd). For detailed introduction of these metrics, see Chapter 7 of Gelman et al. (2013).
The deviance function for a model with a vector of parameters
\boldsymbol \theta
is defined as
D(\boldsymbol \theta) = -2\log p(\mathbf{y} \mid \boldsymbol \theta),
where \mathbf{y}:=(y(\mathbf{x}_1), \ldots, y(\mathbf{x}_n))^\top
is a vector of n
observations.
The effective number of parameters (see p.172 of Gelman et al. 2013) is defined as
pD = E_{\boldsymbol \theta| \mathbf{y}}[D(\boldsymbol \theta)] - D(\hat{ \boldsymbol \theta }),
where
\hat{\boldsymbol \theta} = E_{\boldsymbol \theta | \mathbf{y}}[\boldsymbol \theta].
The interpretation is that the effective number of parameters is the “expected" deviance minus the “fitted" deviance. HigherpD
implies more over-fitting with estimate\hat{\boldsymbol \theta}
.The Deviance information criteria (DIC) (see pp. 172-173 of Gelman et al. 2013) is
DIC = E_{\boldsymbol \theta | \mathbf{y}}[D(\boldsymbol \theta)] + pD.
DIC approximates Akaike information criterion (AIC) and is more appropriate for hierarchical models than AIC and BIC.
The log predictive density (lpd) is defined as
p(y(\mathbf{x}_0) \mid \mathbf{y}) = \int p(y(\mathbf{x}_0) \mid \boldsymbol \theta, \mathbf{y}) p(\boldsymbol \theta \mid \mathbf{y}) d \boldsymbol \theta,
where
\mathbf{y}:=(y(\mathbf{x}_1), \ldots, y(\mathbf{x}_n))^\top
is a vector ofn
observations.\boldsymbol \theta
contains correlation parameters and nugget parameter. This predictive density should be understood as an update of the likelihood since data is treated as prior information now. With a set of prediction locations\mathcal{X}:=\{x_0^i: i=1, \ldots, m\}
. The log pointwise predictive density (lppd) is defined aslppd = \sum_{i=1}^m \log p(y(\mathbf{x}_0^i) \mid \mathbf{y}).
The log joint predictive density (ljpd) is defined as
ljpd = \log p(y(\mathcal{X})).
The
lppd
is connected to cross-validation, while theljpd
measures joint uncertainty across prediction locations.
Usage
gp.model.adequacy(
obj,
testing.input,
testing.output,
pointwise = TRUE,
joint = TRUE
)
Arguments
obj |
a |
testing.input |
a matrix of testing inputs |
testing.output |
a vector of testing outputs |
pointwise |
a logical value with default value |
joint |
a logical value with default value |
Value
a list containing pD, DIC, lppd, ljpd.
Author(s)
Pulong Ma mpulong@gmail.com
References
Gelman, Andrew, John B. Carlin, Hal S. Stern, David B. Dunson, Aki Vehtari, and Donald B. Rubin (2013). Bayesian Data Analysis, Third Edition. CRC Press.
See Also
A wraper to fit a Gaussian stochastic process model with optimization methods
Description
This function is a wraper to estimate parameters in the GaSP model with different choices of estimation methods using numerical optimization methods.
Usage
gp.optim(obj, method = "MMLE", opt = NULL, bound = NULL)
Arguments
obj |
an |
method |
a string indicating the parameter estimation method:
|
opt |
a list of arguments to setup the
|
bound |
Default value is
for the Confluent Hypergeometric covariance and the Cauchy covariance.
|
Value
a list of updated gp
object obj and
fitted information fit
Author(s)
Pulong Ma mpulong@gmail.com
See Also
GPBayes-package, GaSP
, gp, gp.mcmc
Examples
code = function(x){
y = (sin(pi*x/5) + 0.2*cos(4*pi*x/5))*(x<=9.6) + (x/10-1)*(x>9.6)
return(y)
}
n=100
input = seq(0, 20, length=n)
XX = seq(0, 20, length=99)
Ztrue = code(input)
set.seed(1234)
output = Ztrue + rnorm(length(Ztrue), sd=0.1)
obj = gp(formula=~1, output, input,
param=list(range=4, nugget=0.1,nu=2.5),
smooth.est=FALSE,
cov.model=list(family="matern", form="isotropic"))
fit.optim = gp.optim(obj, method="MPLE")
Prediction at new inputs based on a Gaussian stochastic process model
Description
This function provides the capability to make prediction based on a GaSP when different estimation methods are employed.
Usage
gp.predict(obj, input.new, method = "Bayes")
Arguments
obj |
an |
input.new |
a matrix of new input lomessageions |
method |
a string indicating the parameter estimation method:
|
Value
a list of predictive mean, predictive standard deviation, 95% predictive intervals
Author(s)
Pulong Ma mpulong@gmail.com
See Also
GPBayes-package, GaSP
, gp, gp.mcmc
, gp.optim
Examples
code = function(x){
y = (sin(pi*x/5) + 0.2*cos(4*pi*x/5))*(x<=9.6) + (x/10-1)*(x>9.6)
return(y)
}
n=100
input = seq(0, 20, length=n)
XX = seq(0, 20, length=99)
Ztrue = code(input)
set.seed(1234)
output = Ztrue + rnorm(length(Ztrue), sd=0.1)
obj = gp(formula=~1, output, input,
param=list(range=4, nugget=0.1,nu=2.5),
smooth.est=FALSE,
cov.model=list(family="matern", form="isotropic"))
fit.optim = gp.optim(obj, method="MMLE")
obj = fit.optim$obj
pred = gp.predict(obj, input.new=XX, method="MMLE")
Simulate from a Gaussian stochastic process model
Description
This function simulates realizations from Gaussian processes.
Usage
gp.sim(
formula = ~1,
input,
param,
cov.model = list(family = "CH", form = "isotropic"),
dtype = "Euclidean",
nsample = 1,
seed = NULL
)
Arguments
formula |
an object of |
input |
a matrix including inputs in a GaSP |
param |
a list including values for regression parameters, covariance parameters, and nugget variance parameter. The specification of param should depend on the covariance model.
|
cov.model |
a list of two strings: family, form, where family indicates the family of covariance functions including the Confluent Hypergeometric class, the Matérn class, the Cauchy class, the powered-exponential class. form indicates the specific form of covariance structures including the isotropic form, tensor form, automatic relevance determination form.
|
dtype |
a string indicating the type of distance:
|
nsample |
an integer indicating the number of realizations from a Gaussian process |
seed |
a number specifying random number seed |
Value
a numerical vector or a matrix
Author(s)
Pulong Ma mpulong@gmail.com
See Also
Examples
n=50
y.sim = gp.sim(input=seq(0,1,length=n),
param=list(range=0.5,nugget=0.1,nu=2.5),
cov.model=list(family="matern",form="isotropic"),
seed=123)
A wraper to build different kinds of correlation matrices between two sets of inputs
Description
This function wraps existing built-in routines to construct a covariance
matrix for two input matrices based on data type, covariance type, and distance type. The constructed
covariance matrix can be directly used for GaSP fitting and and prediction for spatial
data, spatio-temporal data, and computer experiments. This function explicitly takes inputs as arguments.
The prefix “i" in ikernel
standards for “input".
Usage
ikernel(input1, input2, range, tail, nu, covmodel, dtype = "Euclidean")
Arguments
input1 |
a matrix of input locations |
input2 |
a matrix of input locations |
range |
a vector of range parameters, which could be a scalar. |
tail |
a vector of tail decay parameters, which could be a scalar. |
nu |
a vector of smoothness parameters, which could be a scalar. |
covmodel |
a list of two strings: family, form, where family indicates the family of covariance functions including the Confluent Hypergeometric class, the Matérn class, the Cauchy class, the powered-exponential class. form indicates the specific form of covariance structures including the isotropic form, tensor form, automatic relevance determination form.
|
dtype |
a string indicating distance type: Euclidean, GCD, where the latter indicates great circle distance. |
Value
a correlation matrix
Author(s)
Pulong Ma mpulong@gmail.com
See Also
CH
, matern
, kernel
, GPBayes-package, GaSP
Examples
input = seq(0,1,length=10)
cormat = ikernel(input,input,range=0.5,tail=0.2,nu=2.5,
covmodel=list(family="CH",form="isotropic"))
A wraper to build different kinds of correlation matrices with distance as arguments
Description
This function wraps existing built-in routines to construct a covariance matrix based on data type, covariance type, and distance type with distances as inputs. The constructed covariance matrix can be directly used for GaSP fitting and and prediction for spatial data, spatio-temporal data, and computer experiments.
Usage
kernel(d, range, tail, nu, covmodel)
Arguments
d |
a matrix or a list of distances |
range |
a vector of range parameters, which could be a scalar. |
tail |
a vector of tail decay parameters, which could be a scalar. |
nu |
a vector of smoothness parameters, which could be a scalar. |
covmodel |
a list of two strings: family, form, where family indicates the family of covariance functions including the Confluent Hypergeometric class, the Matérn class, the Cauchy class, the powered-exponential class. form indicates the specific form of covariance structures including the isotropic form, tensor form, automatic relevance determination form.
|
Value
a correlation matrix
Author(s)
Pulong Ma mpulong@gmail.com
See Also
CH
, matern
, ikernel
, GPBayes-package, GaSP
Examples
input = seq(0,1,length=10)
d = distance(input,input,type="isotropic",dtype="Euclidean")
cormat = kernel(d,range=0.5,tail=0.2,nu=2.5,
covmodel=list(family="CH",form="isotropic"))
A wraper to compute the natural logarithm of the integrated likelihood function
Description
This function wraps existing built-in routines to construct the natural logarithm of the integrated likelihood function. The constructed loglikelihood can be directly used for numerical optimization
Usage
loglik(par, output, H, d, covmodel, smooth, smoothness_est)
Arguments
par |
a numerical vector, with which numerical optimization routine such as |
output |
a matrix of outputs |
H |
a matrix of regressors in the mean function of a GaSP model. |
d |
an R object holding the distances. It should be a distance matrix for constructing isotropic correlation matrix, or a list of distance matrices along each input dimension for constructing tensor or ARD types of correlation matrix. |
covmodel |
a list of two strings: family, form, where family indicates the family of covariance functions including the Confluent Hypergeometric class, the Matérn class, the Cauchy class, the powered-exponential class. form indicates the specific form of covariance structures including the isotropic form, tensor form, automatic relevance determination form.
|
smooth |
The smoothness parameter |
smoothness_est |
a logical value indicating whether the smoothness parameter is estimated. |
Value
The natural logarithm of marginal or integrated likelihood
Author(s)
Pulong Ma mpulong@gmail.com
See Also
CH
, matern
, gp.optim
, GPBayes-package, GaSP
The Matérn correlation function proposed by Matérn (1960)
Description
This function computes the Matérn correlation function given a distance matrix. The Matérn correlation function is given by
C(h)=\frac{2^{1-\nu}}{\Gamma(\nu)} \left(\frac{h}{\phi} \right)^{\nu}
\mathcal{K}_{\nu}\left( \frac{h}{\phi} \right),
where \phi
is the range parameter. \nu
is the smoothness parameter.
\mathcal{K}_{\nu}(\cdot)
is the modified Bessel function of the second kind of order \nu
.
The form of covariance includes the following special cases by specifying \nu
to be 0.5, 1.5, 2.5.
\nu=0.5
corresponds to the exponential correlation function (exp) of the formC(h) = \exp\left\{ - \frac{h}{\phi} \right\}
\nu=1.5
corresponds to the Matérn correlation function with smoothness parameter 1.5 (matern_3_2) of the formC(h) = \left( 1 + \frac{h}{\phi} \right) \exp\left\{ - \frac{h}{\phi} \right\}
\nu=2.5
corresponds to the Matérn correlation function with smoothness parameter 2.5 (matern_5_2) of the formC(h) = \left\{ 1 + \frac{h}{\phi} + \frac{1}{3}\left(\frac{h}{\phi}\right)^2 \right\} \exp\left\{ - \frac{h}{\phi} \right\}
Usage
matern(d, range, nu)
Arguments
d |
a matrix of distances |
range |
a numerical value containing the range parameter |
nu |
a numerical value containing the smoothness parameter |
Value
a numerical matrix
Author(s)
Pulong Ma mpulong@gmail.com
See Also
GPBayes-package, GaSP
, gp, CH
, kernel
, ikernel
The powered-exponential correlation function
Description
This function computes the powered-exponential correlation function given a distance matrix. The powered-exponential correlation function is given by
C(h)=\exp\left\{-\left(\frac{h}{\phi}\right)^{\nu}\right\},
where \phi
is the range parameter. \nu
is the smoothness parameter.
The case \nu=2
corresponds to the well-known Gaussian correlation.
Usage
powexp(d, range, nu)
Arguments
d |
a matrix of distances |
range |
a numerical value containing the range parameter |
nu |
a numerical value containing the smoothness parameter |
Value
a numerical matrix
Author(s)
Pulong Ma mpulong@gmail.com
See Also
Print the information an object of the gp
class
Description
Print the information an object of the gp
class
Usage
## S4 method for signature 'gp'
show(object)
Arguments
object |
an object of |