| calibrate.bassBasis {BASS} | R Documentation |
Calibrate a bassPCA or bassBasis Model to Data
Description
Robust modular calibration of a bassPCA or bassBasis emulator using adaptive Metropolis, tempering, and decorrelation steps in an effort to be free of any user-required tuning.
Usage
calibrate.bassBasis(
y,
mod,
type,
sd.est,
s2.df,
s2.ind,
meas.error.cor,
bounds,
discrep.mean,
discrep.mat,
nmcmc = 10000,
temperature.ladder = 1.05^(0:30),
decor.step.every = 100,
verbose = T
)
Arguments
y |
vector of calibration data |
mod |
a emulator of class bassBasis, whose predictions should match y (i.e., predictions from mod should be the same length as y) |
type |
one of c(1,2). 1 indicates a model that uses independent truncation error variance, no measurement error correlation, and discrepancy on a basis while type 2 indicates a model that uses a full truncation error covariance matrix, a full measurement error correlation matrix, a fixed full discrepancy covariance matrix, and a fixed discrepancy mean. 1 is for situations where computational efficiency is important (because y is dense), while 2 is only for cases where y is a short vector. |
sd.est |
vector of prior estimates of measurement error standard deviation |
s2.df |
vector of degrees of freedom for measurement error sd prior estimates |
s2.ind |
index vector, same length as y, indicating which sd.est goes with which y |
meas.error.cor |
a fixed correlation matrix for the measurement errors |
bounds |
a 2xp matrix of bounds for each input parameter, where p is the number of input parameters. |
discrep.mean |
discrepancy mean (fixed), only used if type=2 |
discrep.mat |
discrepancy covariance (fixed, for type 2) or basis (if not square, for type 1) |
nmcmc |
number of MCMC iterations. |
temperature.ladder |
an increasing vector, all greater than 1, for tempering. Geometric spacing is recommended, so that you have (1+delta)^(0:ntemps), where delta is small (typically between 0.05 and 0.2) and ntemps is the number of elements in the vector. |
decor.step.every |
integer number of MCMC iterations between decorrelation steps. |
verbose |
logical, whether to print progress. |
Details
Fits a modular Bayesian calibration model, with
y = Kw(\theta) + Dv + \epsilon, ~~\epsilon \sim N(0,\sigma^2 R)
f(x) = a_0 + \sum_{m=1}^M a_m B_m(x)
and B_m(x) is a BASS basis function (tensor product of spline basis functions). We use priors
a \sim N(0,\sigma^2/\tau (B'B)^{-1})
M \sim Poisson(\lambda)
as well as the priors mentioned in the arguments above.
Value
An object
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)