calibrate.bassBasis: Calibrate a bassPCA or bassBasis Model to Data

View source: R/calib.R

calibrate.bassBasisR 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)

BASS documentation built on July 9, 2023, 6:57 p.m.