# R/tgp.default.params.R In tgp: Bayesian Treed Gaussian Process Models

```#*******************************************************************************
#
# Bayesian Regression and Adaptive Sampling with Gaussian Process Trees
# Copyright (C) 2005, University of California
#
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
# 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 (rbgramacy@ams.ucsc.edu)
#
#*******************************************************************************

## tgp.default.params:
##
## create a default parameter vector for tgp-class models with
## the specified dimension, mean function, correllation model
## and other augmentations specified in ...
##
## Note that the choice of bprior can negate the usefulness of
## (or override) some of the parameters, particularly the hierarchical
## parameters tau2.p and tau2.lam

"tgp.default.params" <-
function(d, meanfn=c("linear", "constant") ,
corr=c("expsep", "exp", "mrexpsep", "matern", "sim", "twovar"),
splitmin=1, basemax=d, ...)
{
## check the d argument, other check in tgp.check.params
if(length(d) != 1)
stop("d should be an integer scalar >= 1")

## check the splitmin argument, other check in tgp.check.params
if(length(splitmin) != 1)
stop("splitmin should be an integer scalar >= 1")

## check the basemax argument, other check in tgp.check.params
if(length(basemax) != 1)
stop("basemax should be an integer scalar >= 1")

## setting of col, the dim of (1,X) based on the mean function
meanfn <- match.arg(meanfn)
if(meanfn == "linear") { col <- d+1 }
else if(meanfn == "constant"){
col <- 1
## not sure why I ever had this code here
## if(basemax != d) {
##    warning("must have basemax = d for constant mean function")
##    basemax <- d
## }
}

## adjust the starting beta and Wi values on basemax
beta <- rep(0, min(col, basemax+1))
Wi <- diag(1, length(beta))

## check the corr argument, and augment splitmin
## if fitting a multi-resolution model
corr <- match.arg(corr)
## PERHAPS THIS SHOULD BE DONE IN THE C CODE SO THAT WHEN WE PRINT WITHIN C IT MAKES MORE SENSE
if(corr=="mrexpsep")  splitmin <- splitmin + 1

## parameters shared by all models
params <-
list(
tree=c(0.5,2,max(c(10,basemax+2)),   # tree prior params <alpha>,<beta>,<minpart>
splitmin, basemax),          # (continued)
col=col,                       # defined above, based on meanfn
meanfn=meanfn,                 # one of "linear" or "constant"
bprior="bflat",		# linear prior (b0, bmle, bflat, b0not, or bmzt)
beta=beta, 		        # start/prior vals for beta
Wi=Wi,                         # start/prior vals for Wi
s2tau2=c(1,1),                 # start vals for s2, and tau2
s2.p=c(5,10),			# s2 prior params (initial values) <a0> and <g0>
s2.lam=c(0.2,10),		# s2 hierarc inv-gamma prior params (or "fixed")
tau2.p=c(5,10),	       	# tau2 prior params (initial values) <a0> and <g0>
tau2.lam=c(0.2,0.1),		# tau2 hierarch inv-gamma prior params (or "fixed")
corr=corr,			# correlation model (exp, expsep, matern, sim)
gd=c(0.1, 0.5),                # start vals for nug and d
nug.p=c(1,1,1,1),		# nug gamma-mix prior params (initial values)
nug.lam="fixed",		# nug hierarch gamma-mix prior params (or "fixed")
gamma=c(10,0.2,0.7),		# gamma linear pdf parameter
d.p=c(1.0,20.0,10.0,10.0),	# d gamma-mix prior params (initial values)
delta.p=c(),                   # delta parameter for high fidelity variance
nugf.p=c(),                    # residual process nugget gamma-mix prior params
d.lam="fixed",			# d lambda hierarch gamma-mix prior params (or "fixed")
dp.sim=diag(0.2,basemax),      # d-proposal covariance for sim correlaton
nu=c()                         # matern correlation smoothing parameter
)

## parameters specific to multi-resolution corr model
if(corr == "mrexpsep"){
mrd.p <- c(1,10,1,10)               # gamma-mix params for the discr process (this for 'wigl')
params\$d.p <- c(params\$d.p, mrd.p)
params\$delta.p <- c(1,1,1,1)
params\$nugf.p <- c(1,1,1,1)
}

## Replace the parameters with ellipsis arguments,
## these should match the entries of params, or be "minpart"
plist <- list( ... )
args <- names(plist)
if(length(plist)>0) {
pmatch <- match(args, c(names(params), "minpart"))
for(i in 1:length(plist)){
if(args[i] == "minpart") params\$tree[3] <- plist[[i]]
else if(!is.na(pmatch[[i]])) params[[pmatch[i]]]<- plist[[i]]
else stop(paste("your argument \"", args[i], "\" is not recognized", sep=""))
}
}

return(params)
}

## tgp.check.params:
##
## check that the parameter list describes a proper hierarchical parameter
## vector (of dimension d) -- and simultabiously convert the list into
## a double-vector to be passed to the C-side of tgp via .C

"tgp.check.params" <-
function(params, d)
{
## check the number of parameters
if(is.null(params)) return(matrix(-1));
if(length(params) != 22) {
stop(paste("Number of params should be 22, you have", length(params), "\n"));
}

## tree prior parameters
if(length(params\$tree) != 5) {
stop(paste("length of params\$tree should be 5, you have",
length(params\$tree), "\n"));
}

## check tree minpart is bigger than input dimension
basemax <- params\$tree[5]
if(params\$tree[3] <= basemax) {
stop(paste("tree minpart", params\$tree[3], "should be > basemax =", basemax, "\n"));
}

## check tree splitmin is <= than input dimension
if(params\$tree[4] < 1 || params\$tree[4] > d) {
stop(paste("tree splitmin", params\$tree[4], "should be >= 1 and <= d =", d, "\n"));
}

## check tree basemax is > splitmin and <= than input dimension
if(basemax < 1 || params\$tree[5] > d) {
stop(paste("tree basemax", basemax, "should be >= 1 and <= d =", d, "\n"));
}

## tree alpha and beta parameters
p <- c(as.numeric(params\$tree))

## assign the mean function
if(params\$meanfn == "linear") {
meanfn <- 0;
if(params\$col != d+1)
stop(paste("col=", params\$col, " should be d+1=", d+1,
"with linear mean function",  sep=""))
} else if(params\$meanfn == "constant"){
meanfn <- 1;
if(params\$col != 1)
stop(paste("col=", params\$col,
" should be 1 with constant mean function",  sep=""))
} else { cat(paste("params\$meanfn =", params\$meanfn, "not valid\n")); meanfn <- 0; }
p <- c(p, meanfn)

## beta linear prior model
## check the type of beta prior, and possibly augment by p0
if(params\$bprior == "b0") { p <- c(p,0);  }
else if(params\$bprior == "bmle") { p <- c(p, 1); }
else if(params\$bprior == "bflat") { p <- c(p, 2); }
else if(params\$bprior == "b0not") { p <- c(p, 3); }
else if(params\$bprior == "bmzt") { p <- c(p, 4); }
else if(params\$bprior == "bmznot") { p <- c(p, 5); }
else { stop(paste("params\$bprior =", params\$bprior, "not valid\n")); }

## initial settings of beta linear prior mean parameters
if(length(params\$beta) != min(params\$col, basemax+1)) {
stop(paste("length of params\$beta should be", min(params\$col, basemax+1),
"you have", length(params\$beta), "\n"));
}

## finally, set the params\$beta
p <- c(p, as.numeric(params\$beta))

## initial settings of the beta linear prior correlation parameters
if(nrow(params\$Wi) != length(params\$beta) && ncol(params\$Wi) != nrow(params\$Wi)) {
stop(paste("params\$Wi should be", length(params\$beta), "x", length(params\$beta),
"you have", nrow(params\$Wi), "x", ncol(params\$Wi), "\n"));
}

## finally, set the params\$Wi
p <- c(p, as.numeric(params\$Wi))

## initial settings of variance parameters
if(length(params\$s2tau2) != 2) {
stop(paste("length of params\$s2tau2 should be 2 you have",
length(params\$s2tau2), "\n"));
}
p <- c(p, as.numeric(params\$s2tau2))

## sigma^2 prior parameters
if(length(params\$s2.p) != 2) {
stop(paste("length of params\$s2.p should be 2 you have",
length(params\$s2.p), "\n"));
}
p <- c(p, as.numeric(params\$s2.p))

## hierarchical prior parameters for sigma^2 (exponentials) or "fixed"
if(length(params\$s2.lam) != 2 && params\$s2.lam[1] != "fixed") {
stop(paste("length of params\$s2.lam should be 2 or fixed, you have",
params\$s2.lam, "\n"));
}
if(params\$s2.lam[1] == "fixed") p <- c(p, rep(-1, 2))
else p <- c(p, as.numeric(params\$s2.lam))

## tau^2 prior parameters
if(length(params\$tau2.p) != 2) {
stop(paste("length of params\$tau2.p should be 2 you have",
length(params\$tau2.p),"\n"));
}
p <- c(p, as.numeric(params\$tau2.p))

## hierarchical prior parameters for tau^2 (exponentials) or "fixed"
if(length(params\$tau2.lam) != 2 && params\$tau2.lam[1] != "fixed") {
stop(paste("length of params\$s2.lam should be 2 or fixed, you have",
params\$tau2.lam, "\n"));
}
if(params\$tau2.lam[1] == "fixed") p <- c(p, rep(-1, 2))
else p <- c(p, as.numeric(params\$tau2.lam))

## correllation model
if(params\$corr == "exp") { p <- c(p, 0); }
else if(params\$corr == "expsep") { p <- c(p, 1); }
else if(params\$corr == "matern") { p <- c(p, 2); }
else if(params\$corr == "mrexpsep") { p <- c(p,3) }
else if(params\$corr == "sim") { p <- c(p,4) }
else if(params\$corr == "twovar") { p <- c(p,5) }
else { stop(paste("params\$corr =", params\$corr, "not valid\n")); }

## initial settings of variance parameters
if(length(params\$gd) != 2) {
stop(paste("length of params\$gd should be 2 you have",
length(params\$gd), "\n"));
}
p <- c(p, as.numeric(params\$gd))

## mixture of gamma (initial) prior parameters for nug
if(length(params\$nug.p) == 1 && params\$nug.p[1] == 0) params\$nug.p <- rep(0,4)
if(length(params\$nug.p) != 4) {
stop(paste("length of params\$nug.p should be 4 you have",
length(params\$nug.p),"\n"));
}
if(params\$nug.p[1] == 0) params\$nug.p[2] <- params\$gd[1]
p <- c(p, as.numeric(params\$nug.p))

## hierarchical prior params for nugget g (exponentials) or "fixed"
if(length(params\$nug.lam) != 4 && params\$nug.lam[1] != "fixed") {
stop(paste("length of params\$nug.lam should be 4 or fixed, you have",
params\$nug.lam, "\n"));
}
if(params\$nug.lam[1] == "fixed") p <- c(p, rep(-1, 4))
else p <- c(p, as.numeric(params\$nug.lam))

## gamma theta1 theta2 LLM prior params
if(length(params\$gamma) != 3) {
stop(paste("length of params\$gamma should be 3, you have",
length(params\$gamma),"\n"));
}
if(params\$gamma[1] > 0 && params\$corr == "sim") stop("cannot have sim corr with LLM")
if(!prod(params\$gamma[2:3] > 0)) stop("all params\$gamma[2:3] must be positive\n")
if(sum(params\$gamma[2:3]) >= 1.0) stop("sum(gamma[2:3]) > 1 not allowed\n")
p <- c(p, as.numeric(params\$gamma))

## mixture of gamma (initial) prior parameters for range parameter d

## if(length(params\$d.p) == 1 && params\$d.p[1] == 0) params\$d.p <- rep(0,4)
if(length(params\$d.p) != 8 && params\$corr == "mrexpsep") {
stop(paste("length of params\$d.p should be 8 you have",
length(params\$d.p),"\n"));
}
else if( length(params\$d.p) != 4 && params\$corr != "mrexpsep" ) {
stop(paste("length of params\$d.p should be 4 you have",
length(params\$d.p),"\n"));
}
if(params\$d.p[1] == 0) params\$d.p[2] <- params\$gd[2]
if(length(params\$d.p) == 8 && params\$d.p[5] == 0) params\$d.p[6] <- params\$gd[2]

## finally, set the params\$d.p
p <- c(p, as.numeric(params\$d.p))

## delta.p -- only do this if we are using mrexpsep
if(length(params\$delta.p) != 4 && params\$corr == "mrexpsep") {
stop(paste("length of params\$delta.p should be 4 you have",
length(params\$delta.p),"\n"));
}
if(params\$corr == "mrexpsep") p<- c(p, as.numeric(params\$delta.p))

## nugf.p -- only do this if we are using mrexpsep
if(length(params\$nugf.p) != 4 && params\$corr == "mrexpsep") {
stop(paste("length of params\$delta.p should be 4 you have",
length(params\$nug.p),"\n"));
}
if(params\$corr == "mrexpsep") p<- c(p, as.numeric(params\$nugf.p))

## hierarchical prior params for range d (exponentials) or "fixed"
if(length(params\$d.lam) != 4 && params\$d.lam[1] != "fixed") {
stop(paste("length of params\$d.lam should be 4 or fixed, you have",
length(params\$d.lam),"\n"));
}
if(params\$d.lam[1] == "fixed") p <- c(p, rep(-1, 4))
else p <- c(p, as.numeric(params\$d.lam))

## add sd for proposals for sim-d parameters
if(params\$corr == "sim") {
if(nrow(params\$dp.sim) != basemax || ncol(params\$dp.sim) != basemax)
stop("dp.sim should be ", basemax, "x", basemax, "\n");
p <- c(p,as.numeric(params\$dp.sim))
}

## nu smoothness parameter for Matern correlation function
if(params\$corr == "matern") {
if(params\$nu < 0) stop(paste("nu should be greater than zero, you have",
params\$nu, "\n"))
}
p <- c(p, as.numeric(params\$nu))

## return the constructed double-vector of parameters for C
return(p)
}
```

## Try the tgp package in your browser

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

tgp documentation built on Sept. 7, 2020, 5:10 p.m.