EFMBC Exmaple

library(impala)
library(BASS)
library(mvBayes)
library(fdasrvf)
library(bayesplot)
#> This is bayesplot version 1.15.0
#> - Online documentation and vignettes at mc-stan.org/bayesplot
#> - bayesplot theme set to bayesplot::theme_default()
#>    * Does _not_ affect other ggplot2 plots
#>    * See ?bayesplot_theme_set for details on theme setting

Generate Example Data

f <- function(x) {
  dnorm(seq(0, 1, length.out = 99),
        sin(2 * pi * x[1] ^ 2) / 4 - x[1] / 10 + 0.5,
        0.05) * x[2]
}

n = 100
nt = 99
p = 3
x_train = matrix(runif(n * p), n)
x_test = matrix(runif(1000 * p), 1000)
e = rnorm(n * 99)
y_train = matrix(0, n, nt)
for (i in 1:n) {
  y_train[i, ] = f(x_train[i, ])
}
y_test = matrix(0, 1000, nt)
for (i in 1:1000) {
  y_test[i, ] = f(x_test[i, ])
}

Generate observation and align using EFDA framework

x_true = c(0.1028, 0.5930)
ftilde_obs = f(x_true)
gam_obs = seq(0, 1, length.out = nt)
vv_obs = gam_to_v(gam_obs)

tt = seq(0, 1, length.out = nt)
out = multiple_align_functions(t(y_train), tt, ftilde_obs, 0.01)
#> ℹ Using lambda = 0.01

#> Aligning 100 functions in SRSF space...

gam_train = out$warping_functions
vv_train = gam_to_v(gam_train)
ftilde_train = out$fn
qtilde_train = out$qn

matplot(tt, t(y_train), type = "l", col = "gray")
lines(tt, ftilde_obs, col = "black")
legend(
  x = 0,
  y = 8,
  legend = c("simulation", "experiment"),
  col = c("gray", "black"),
  lty = 1
)


matplot(tt, ftilde_train, type = "l", col = "gray")
lines(tt, ftilde_obs, col = "black")
legend(
  x = 0,
  y = 8,
  legend = c("simulation", "experiment"),
  col = c("gray", "black"),
  lty = 1
)


matplot(tt, gam_train, type = "l", col = "gray")
lines(tt, gam_obs, col = "black")
legend(
  x = 0,
  y = 8,
  legend = c("simulation", "experiment"),
  col = c("gray", "black"),
  lty = 1
)


matplot(tt, vv_train, type = "l", col = "gray")
lines(tt, vv_obs, col = "black")
legend(
  x = 0,
  y = 8,
  legend = c("simulation", "experiment"),
  col = c("gray", "black"),
  lty = 1
)

Train Emulators

emu_ftilde = mvBayes(bass, x_train, t(ftilde_train), nBasis=2)
#> Starting mvBayes with 2 components, using 1 cores.
#> MCMC Start #-- Jun 12 10:31:11 --# nbasis: 0 
#> MCMC iteration 1000 #-- Jun 12 10:31:11 --# nbasis: 2 
#> MCMC iteration 2000 #-- Jun 12 10:31:11 --# nbasis: 2 
#> MCMC iteration 3000 #-- Jun 12 10:31:11 --# nbasis: 2 
#> MCMC iteration 4000 #-- Jun 12 10:31:12 --# nbasis: 2 
#> MCMC iteration 5000 #-- Jun 12 10:31:12 --# nbasis: 2 
#> MCMC iteration 6000 #-- Jun 12 10:31:12 --# nbasis: 2 
#> MCMC iteration 7000 #-- Jun 12 10:31:12 --# nbasis: 2 
#> MCMC iteration 8000 #-- Jun 12 10:31:12 --# nbasis: 3 
#> MCMC iteration 9000 #-- Jun 12 10:31:12 --# nbasis: 2 
#> MCMC iteration 10000 #-- Jun 12 10:31:12 --# nbasis: 2 
#> MCMC Start #-- Jun 12 10:31:12 --# nbasis: 0 
#> MCMC iteration 1000 #-- Jun 12 10:31:12 --# nbasis: 2 
#> MCMC iteration 2000 #-- Jun 12 10:31:12 --# nbasis: 2 
#> MCMC iteration 3000 #-- Jun 12 10:31:12 --# nbasis: 2 
#> MCMC iteration 4000 #-- Jun 12 10:31:12 --# nbasis: 1 
#> MCMC iteration 5000 #-- Jun 12 10:31:12 --# nbasis: 2 
#> MCMC iteration 6000 #-- Jun 12 10:31:13 --# nbasis: 2 
#> MCMC iteration 7000 #-- Jun 12 10:31:13 --# nbasis: 3 
#> MCMC iteration 8000 #-- Jun 12 10:31:13 --# nbasis: 1 
#> MCMC iteration 9000 #-- Jun 12 10:31:13 --# nbasis: 3 
#> MCMC iteration 10000 #-- Jun 12 10:31:13 --# nbasis: 2
#> Generating 'samples' attribute, since it was absent in 'bmList[[1]]'
#> Approximating 'residSD', since 'residSDExtract' is None
plot(emu_ftilde)
#> Warning in predict.mvBayes(x, Xtest, returnPostCoefs = TRUE, idxSamples =
#> idxSamples): 'idxSamples' is not an argument of the bayesModel predict
#> function...setting idxSamples='default'


emu_vv = mvBayes(bass, x_train, t(vv_train), nBasis=2)
#> Starting mvBayes with 2 components, using 1 cores.
#> MCMC Start #-- Jun 12 10:31:13 --# nbasis: 0 
#> MCMC iteration 1000 #-- Jun 12 10:31:14 --# nbasis: 8 
#> MCMC iteration 2000 #-- Jun 12 10:31:14 --# nbasis: 10 
#> MCMC iteration 3000 #-- Jun 12 10:31:14 --# nbasis: 10 
#> MCMC iteration 4000 #-- Jun 12 10:31:14 --# nbasis: 9 
#> MCMC iteration 5000 #-- Jun 12 10:31:14 --# nbasis: 9 
#> MCMC iteration 6000 #-- Jun 12 10:31:14 --# nbasis: 9 
#> MCMC iteration 7000 #-- Jun 12 10:31:14 --# nbasis: 9 
#> MCMC iteration 8000 #-- Jun 12 10:31:14 --# nbasis: 9 
#> MCMC iteration 9000 #-- Jun 12 10:31:14 --# nbasis: 9 
#> MCMC iteration 10000 #-- Jun 12 10:31:14 --# nbasis: 9 
#> MCMC Start #-- Jun 12 10:31:14 --# nbasis: 0 
#> MCMC iteration 1000 #-- Jun 12 10:31:14 --# nbasis: 2 
#> MCMC iteration 2000 #-- Jun 12 10:31:14 --# nbasis: 2 
#> MCMC iteration 3000 #-- Jun 12 10:31:14 --# nbasis: 3 
#> MCMC iteration 4000 #-- Jun 12 10:31:14 --# nbasis: 2 
#> MCMC iteration 5000 #-- Jun 12 10:31:15 --# nbasis: 2 
#> MCMC iteration 6000 #-- Jun 12 10:31:15 --# nbasis: 3 
#> MCMC iteration 7000 #-- Jun 12 10:31:15 --# nbasis: 2 
#> MCMC iteration 8000 #-- Jun 12 10:31:15 --# nbasis: 2 
#> MCMC iteration 9000 #-- Jun 12 10:31:15 --# nbasis: 2 
#> MCMC iteration 10000 #-- Jun 12 10:31:15 --# nbasis: 2
#> Generating 'samples' attribute, since it was absent in 'bmList[[1]]'
#> Approximating 'residSD', since 'residSDExtract' is None
plot(emu_vv)
#> Warning in predict.mvBayes(x, Xtest, returnPostCoefs = TRUE, idxSamples =
#> idxSamples): 'idxSamples' is not an argument of the bayesModel predict
#> function...setting idxSamples='default'

Run Calibration using Impala

input_names = c("theta0", "theta1", "theta2")
bounds = list()
bounds[['theta0']] = c(0, 1)
bounds[['theta1']] = c(0, 1)
bounds[['theta2']] = c(0, 1)

setup = CalibSetup(bounds, cf_bounds)

model_ftilde = ModelmvBayes(emu_ftilde, input_names)
model_vv = ModelmvBayes(emu_vv, input_names)

setup = addVecExperiments(setup, t(ftilde_obs), model_ftilde, 0.01, 20, rep(1, nt))
setup = addVecExperiments(setup, t(vv_obs), model_vv, 0.01, 20, rep(1, nt))
setup = setTemperatureLadder(setup, 1.05 ^ (0:3))
setup = setMCMC(setup, 4000, 2000, 1, 10)
out_cal = calibPool(setup)

Look at posteriors and chain diagnostics

uu = seq(2500, 4000, 2)
theta = tran_unif(out_cal$theta[uu, 1, ], setup$bounds_mat, names(setup$bounds))
expnums = 1

cnt = 1
ftilde_pred_obs = vector(mode = "list", length = setup$nexp)
gam_pred_obs = vector(mode = "list", length = setup$nexp)
for (idx in expnums) {
  time_new = seq(0, 1, length.out = nt)

  fn = setup$models[[1]]$input_names
  parmat_array = matrix(0, length(theta[[fn[1]]]), length(fn))
  for (i in 1:length(fn)) {
    parmat_array[, i] = theta[[fn[i]]]
  }

  ftilde_pred_obs[[idx]] = evalm(setup$models[[cnt]], theta)
  vv_pred_obs = evalm(setup$models[[cnt + 1]], theta)
  cnt = cnt + 2
  gam_pred_obs[[idx]] = v_to_gam(t(vv_pred_obs))

  matplot(
    time_new,
    ftilde_train,
    type = "l",
    col = "gray",
    lty = 1
  )
  matplot(
    time_new,
    t(ftilde_pred_obs[[idx]]),
    type = "l",
    col = "lightblue",
    lty = 1,
    add = TRUE
  )
  lines(time_new, ftilde_obs, col = "black")
  legend(
    x = 0,
    y = 8,
    legend = c("simulation", "prediction", "experiment"),
    col = c("gray", "lightblue", "black"),
    lty = 1
  )

  matplot(
    time_new,
    vv_train,
    type = "l",
    col = "gray",
    lty = 1
  )
  matplot(
    time_new,
    t(vv_pred_obs),
    type = "l",
    col = "lightblue",
    lty = 1,
    add = TRUE
  )
  lines(time_new, vv_obs, col = "black")
  legend(
    x = 0,
    y = 8,
    legend = c("simulation", "prediction", "experiment"),
    col = c("gray", "lightblue", "black"),
    lty = 1
  )

  matplot(
    time_new,
    gam_train,
    type = "l",
    col = "gray",
    lty = 1
  )
  matplot(
    time_new,
    gam_pred_obs[[idx]],
    type = "l",
    col = "lightblue",
    lty = 1,
    add = TRUE
  )
  lines(time_new, gam_obs, col = "black")
  legend(
    x = 0,
    y = 8,
    legend = c("simulation", "prediction", "experiment"),
    col = c("gray", "lightblue", "black"),
    lty = 1
  )
}


samps = data.frame(theta)
mcmc_pairs(samps, off_diag_fun="hex", diag_fun = "dens")
#> Warning: Only one chain in 'x'. This plot is more useful with multiple chains.


mcmc_hist(samps)
#> `stat_bin()` using `bins = 30`. Pick better value `binwidth`.


mcmc_trace(samps)


plot(out_cal$llik, type="l", main="Log-Likelihood")


# misaligned prediction
for (idx in expnums){
  obspred = matrix(0, nt, length(uu))
  for (j in 1:length(uu)){
    obspred[,j] = warp_f_gamma(ftilde_pred_obs[[idx]][j,], time_new, invertGamma(gam_pred_obs[[idx]][,j]))
  }

  matplot(
    time_new,
    t(y_train),
    type = "l",
    col = "gray",
    lty = 1
  )
  matplot(
    time_new,
    obspred,
    type = "l",
    col = "lightblue",
    lty = 1,
    add = TRUE
  )
  lines(time_new, ftilde_obs, col = "black")
  legend(
    x = 0,
    y = 8,
    legend = c("simulation", "prediction", "experiment"),
    col = c("gray", "lightblue", "black"),
    lty = 1
  )

}