# Calibrate a bassPCA or bassBasis Model to Data 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. ```r 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. ## Returns An object ## 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. ## Examples ```r ## 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) ``` ## See Also predict.bassBasis for prediction and sobolBasis for sensitivity analysis.

calibrate.bassBasis function