4
$\begingroup$

I'm using the tgp package in R for fully Bayesian Gaussian Process Regression, and it's great! I'm currently performing regression for experimental data coming from turbomachinery testing, and I'm using the bgp function. This function uses a GP prior with either a linear mean or a constant mean (respectively, option meanfn="linear" or meanfn="constant", which is the default). Note that tgp allows the use of treed Gaussian priors, but for now I'm staying simple, so I'm using the bgp function which doesn't use regression trees, just ordinary Gaussian Processes.

I would like my posterior predictive mean to go to zero away from the training set data, for physical reasons. How can I impose that? I was thinking to set the prior over $\beta_0$ to a Normal distribution centered at 0 and with an extremely small variance, but I'm not sure how to do that. From help(btgp)

bprior Linear (beta) prior, default is "bflat"; alternates include "b0" hierarchical Normal
prior, "bmle" empirical Bayes Normal prior, "b0not" Bayesian treed LMstyle
prior from Chipman et al. (same as "b0" but without tau2), "bmzt" a independent
Normal prior (mean zero) with inverse-gamma variance (tau2), and
"bmznot" is the same as "bmznot" without tau2. The default "bflat" gives
an “improper” prior which can perform badly when the signal-to-noise ratio is
low. In these cases the “proper” hierarchical specification "b0" or independent
"bmzt" or "bmznot" priors may perform better

Default is the improper prior "bflat", which is not what I want. If I use the "b0" hierarchical Normal prior, I guess I cannot set the mean and the variance because they should become additional hyperparameters to be determined in the Bayesian paradigm. Thus, I may go for "bmzt", the independent Normal prior with zero mean. However, with this prior I cannot set the variance, which is again an hyperparameter. Basically, I want my prior mean function to be zero, so that away from the data, also the posterior predictive mean will be zero. Is there a way to achieve that?

EDIT: nobody wants to have a try? :) As my actual case is quite complicated, I wrote a small test case which illustrates the main problem, with the help of the tgp package author. NOTE: unless you have an optimized version of R, you may want to set BTE = c(1000,10000,2) in the call to bgp, or you may have to wait for a very long time to get an answer.

# clear the workspace
rm(list=ls())
gc()
graphics.off()

# set seed for reproducibility
set.seed(825)

# load required packages
library(tgp)
library(ggplot2)

# simulated data
x <- seq(-1,1,len=100)
eps <- rnorm(n=100,mean=0,sd=0.5)
y <- -5*x^2+eps
ymean <- mean(y)

# prediction points
xpred <- seq(-20,20,len=100)

# fit GP
GPModel <- bgp(X=x,Z=y,XX=xpred,meanfn = "constant", bprior="bmzt", 
BTE = c(2000,52000,2), tau2.p=c(1,10000), tau2.lam="fixed")    
ypred <- GPModel$ZZ.mean 

# plots
ymean_vector <- rep(ymean,100)
df <- data.frame(x,y,xpred,ypred,ymean_vector)
p <- ggplot(data=df)
p <- p + geom_point(aes(x=x,y=y)) + 
    geom_line(aes(x=xpred,y=ypred),col="blue") +
    geom_line(aes(x=xpred,y=ymean_vector),col="red") +
    geom_line(aes(x=xpred,y=GPModel$ZZ.q1), col="green") + 
    geom_line(aes(x=xpred,y=GPModel$ZZ.q2), col="green")
p

The resulting plot is

enter image description here

The mean response is the red line: the blue line is the GP posterior predictive mean, and the green lines give the 90% credible interval.Thus, outside the training data range, the data mean is indeed included in the 90% credible interval, but I would like the predictive mean to converge to it...I think that if I could find a way to set the standard deviation of the prior for $\beta_0$ to some extremely small value, I would achieve what I want, but I don't know how to do it.

EDIT2: I can use either a multiplicative (separable) squared exponential kernel or an additive squared exponential kernel.

sep_Gaussian_Kernel <- function(x,y,sigma,l) {
    prod(sigma*exp(-0.5*(abs(x-y)/l)^2))
}    

add_Gaussian_Kernel <- function(x,y,sigma,l) {
    sum(sigma*exp(-0.5*(abs(x-y)/l)^2))/length(x)
} 
$\endgroup$
  • $\begingroup$ A fully Bayesian model over which you want fine-grained control to tinker with the model specification? It sounds like you'd be interested in using stan. $\endgroup$ – Sycorax Mar 29 '16 at 15:53
  • $\begingroup$ Probably! What's stan? install.packages("stan") tells me there's no stan package. I'd really like an R solution, though. It has been already quite a pain to wrangle my data and import them in R, I'd rather not export them and import them in some other environment, unless there's really no other way... $\endgroup$ – DeltaIV Mar 29 '16 at 16:06
  • $\begingroup$ Stan is multi-platform. Try install.packages("rstan") ;-) or visit www.mc-stan.org. $\endgroup$ – Sycorax Mar 29 '16 at 16:07
  • $\begingroup$ will do, thanks! Hope it doesn't have a steep learning curve ;) any chance you could redo my example in rstan? That would be a great way for me to learn, and it would help me tackle my actual problem, which is of course more complex than the one shown here. $\endgroup$ – DeltaIV Mar 29 '16 at 16:17
  • $\begingroup$ thanks! That would be fantastic :) either multiplicative squared exponential or additive squared exponential would do. In my example they make no difference because there's just one predictor, but in the real case I have more. I'm adding their definition in my question. $\endgroup$ – DeltaIV Mar 29 '16 at 16:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.