bassPCA: Bayesian Adaptive Spline Surfaces (BASS) with PCA...

In lanl/BASS: Bayesian Adaptive Spline Surfaces

Bayesian Adaptive Spline Surfaces (BASS) with PCA decomposition of response

Description

Decomposes a multivariate or functional response onto a principal component basis and fits a BASS model to each basis coefficient.

Usage

bassPCA(
  xx = NULL,
  y = NULL,
  dat = NULL,
  n.pc = NULL,
  perc.var = 99,
  n.cores = 1,
  parType = "fork",
  center = T,
  scale = F,
  ...
)

Arguments

xx	a data frame or matrix of predictors with dimension n x p. Categorical predictors should be included as factors.
y	a response matrix (functional response) with dimension n x m.
dat	optional (for more control) list with elements xx (same as above), y (same as above), n.pc (number of principal components used), basis (principal components with dimension m x n.pc), newy (reduced dimension y with dimension n.pc x n), trunc.error (optional truncation error with dimension n x m), y.m (vector mean removed before PCA with dimension m), y.s (vector sd scaled before PCA with dimension m). If dat is specified, xx, y and n.pc do not need to be specified.
n.pc	number of principal components to use
perc.var	optionally specify percent of variance to explain instead of n.pc
n.cores	integer number of cores (threads) to use
parType	either "fork" or "socket". Forking is typically faster, but not compatible with Windows. If n.cores==1, parType is ignored.
center	logical whether to subtract the mean before getting the principal components, or else a numeric vector of dimension m for the center to be used
scale	logical whether to divide by the standard deviation before getting the principal components, or else a numeric vector of dimension m for the scale to be used
...	arguements to be passed to bass function calls.

Details

Gets the PCA decomposition of the response y, and fits a bass model to each PCA basis coefficient, bass(dat$xx,dat$newy[i,],...) for i in 1 to n.pc, possibly in parallel.

Value

An object of class 'bassBasis' with two elements:

mod.list	list (of length n.pc) of individual bass models
dat	same as dat above

See Also

predict.bassBasis for prediction and sobolBasis for sensitivity analysis.

Examples

## Not run: 
  ## simulate data (Friedman function)
f<-function(x){
  10*sin(pi*x[,1]*x[,2])+20*(x[,3]-.5)^2+10*x[,4]+5*x[,5]
}
## simulate data (Friedman function with first variable as functional)
sigma<-.1 # noise sd
n<-500 # number of observations
nfunc<-50 # size of functional variable grid
xfunc<-seq(0,1,length.out=nfunc) # functional grid
x<-matrix(runif(n*9),n,9) # 9 non-functional variables, only first 4 matter
X<-cbind(rep(xfunc,each=n),kronecker(rep(1,nfunc),x)) # to get y
y<-matrix(f(X),nrow=n)+rnorm(n*nfunc,0,sigma)

## fit BASS
library(parallel)
mod<-bassPCA(x,y,n.pc=5,n.cores=min(5,parallel::detectCores()))
plot(mod$mod.list[[1]])
plot(mod$mod.list[[2]])
plot(mod$mod.list[[3]])
plot(mod$mod.list[[4]])
plot(mod$mod.list[[5]])

hist(mod$dat$trunc.error)

## prediction
npred<-100
xpred<-matrix(runif(npred*9),npred,9)
Xpred<-cbind(rep(xfunc,each=npred),kronecker(rep(1,nfunc),xpred))
ypred<-matrix(f(Xpred),nrow=npred)
pred<-predict(mod,xpred,mcmc.use=1:1000) # posterior predictive samples (each is a curve)
matplot(ypred,apply(pred,2:3,mean),type='l',xlab='observed',ylab='mean prediction')
abline(a=0,b=1,col=2)
matplot(t(ypred),type='l') # actual
matplot(t(apply(pred,2:3,mean)),type='l') # mean prediction

## sensitivity
sens<-sobolBasis(mod,int.order = 2,ncores = max(parallel::detectCores()-2,1),
                 mcmc.use=1000) # for speed, only use a few samples
plot(sens)


## calibration
x.true<-runif(9,0,1) # what we are trying to learn
yobs<-f(cbind(xfunc,kronecker(rep(1,nfunc),t(x.true)))) +
  rnorm(nfunc,0,.1) # calibration data (with measurement error)
plot(yobs)

cal<-calibrate.bassBasis(y=yobs,mod=mod,
                         discrep.mean=rep(0,nfunc),
                         discrep.mat=diag(nfunc)[,1:2]*.0000001,
                         sd.est=.1,
                         s2.df=50,
                         s2.ind=rep(1,nfunc),
                         meas.error.cor=diag(nfunc),
                         bounds=rbind(rep(0,9),rep(1,9)),
                         nmcmc=10000,
                         temperature.ladder=1.05^(0:30),type=1)

nburn<-5000
uu<-seq(nburn,10000,5)

pairs(rbind(cal$theta[uu,1,],x.true),col=c(rep(1,length(uu)),2),ylim=c(0,1),xlim=c(0,1))

pred<-apply(predict(mod,cal$theta[uu,1,],nugget = T,trunc.error = T,
      mcmc.use = cal$ii[uu]),3,function(x) diag(x)+rnorm(length(uu),0,sqrt(cal$s2[uu,1,1])))
qq<-apply(pred,2,quantile,probs=c(.025,.975))
matplot(t(qq),col='lightgrey',type='l')
lines(yobs,lwd=3)



## End(Not run)
## minimal example for CRAN testing
mod<-bassPCA(1:10,matrix(1:20,10),n.pc=2,nmcmc=2,nburn=1)

lanl/BASS documentation built on May 15, 2024, 6:40 p.m.