The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

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 #-- Jul 10 18:59:40 --# nbasis: 0 
#> MCMC iteration 1000 #-- Jul 10 18:59:40 --# nbasis: 1 
#> MCMC iteration 2000 #-- Jul 10 18:59:40 --# nbasis: 1 
#> MCMC iteration 3000 #-- Jul 10 18:59:40 --# nbasis: 1 
#> MCMC iteration 4000 #-- Jul 10 18:59:40 --# nbasis: 2 
#> MCMC iteration 5000 #-- Jul 10 18:59:40 --# nbasis: 1 
#> MCMC iteration 6000 #-- Jul 10 18:59:40 --# nbasis: 1 
#> MCMC iteration 7000 #-- Jul 10 18:59:40 --# nbasis: 1 
#> MCMC iteration 8000 #-- Jul 10 18:59:40 --# nbasis: 2 
#> MCMC iteration 9000 #-- Jul 10 18:59:40 --# nbasis: 1 
#> MCMC iteration 10000 #-- Jul 10 18:59:40 --# nbasis: 1 
#> MCMC Start #-- Jul 10 18:59:40 --# nbasis: 0 
#> MCMC iteration 1000 #-- Jul 10 18:59:40 --# nbasis: 5 
#> MCMC iteration 2000 #-- Jul 10 18:59:41 --# nbasis: 1 
#> MCMC iteration 3000 #-- Jul 10 18:59:41 --# nbasis: 0 
#> MCMC iteration 4000 #-- Jul 10 18:59:41 --# nbasis: 0 
#> MCMC iteration 5000 #-- Jul 10 18:59:41 --# nbasis: 2 
#> MCMC iteration 6000 #-- Jul 10 18:59:41 --# nbasis: 2 
#> MCMC iteration 7000 #-- Jul 10 18:59:41 --# nbasis: 2 
#> MCMC iteration 8000 #-- Jul 10 18:59:41 --# nbasis: 2 
#> MCMC iteration 9000 #-- Jul 10 18:59:41 --# nbasis: 1 
#> MCMC iteration 10000 #-- Jul 10 18:59:41 --# nbasis: 1
#> 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 #-- Jul 10 18:59:42 --# nbasis: 0 
#> MCMC iteration 1000 #-- Jul 10 18:59:42 --# nbasis: 6 
#> MCMC iteration 2000 #-- Jul 10 18:59:42 --# nbasis: 6 
#> MCMC iteration 3000 #-- Jul 10 18:59:42 --# nbasis: 6 
#> MCMC iteration 4000 #-- Jul 10 18:59:42 --# nbasis: 6 
#> MCMC iteration 5000 #-- Jul 10 18:59:42 --# nbasis: 6 
#> MCMC iteration 6000 #-- Jul 10 18:59:42 --# nbasis: 6 
#> MCMC iteration 7000 #-- Jul 10 18:59:42 --# nbasis: 6 
#> MCMC iteration 8000 #-- Jul 10 18:59:42 --# nbasis: 6 
#> MCMC iteration 9000 #-- Jul 10 18:59:43 --# nbasis: 6 
#> MCMC iteration 10000 #-- Jul 10 18:59:43 --# nbasis: 6 
#> MCMC Start #-- Jul 10 18:59:43 --# nbasis: 0 
#> MCMC iteration 1000 #-- Jul 10 18:59:43 --# nbasis: 2 
#> MCMC iteration 2000 #-- Jul 10 18:59:43 --# nbasis: 3 
#> MCMC iteration 3000 #-- Jul 10 18:59:43 --# nbasis: 3 
#> MCMC iteration 4000 #-- Jul 10 18:59:43 --# nbasis: 2 
#> MCMC iteration 5000 #-- Jul 10 18:59:43 --# nbasis: 8 
#> MCMC iteration 6000 #-- Jul 10 18:59:43 --# nbasis: 5 
#> MCMC iteration 7000 #-- Jul 10 18:59:43 --# nbasis: 6 
#> MCMC iteration 8000 #-- Jul 10 18:59:43 --# nbasis: 6 
#> MCMC iteration 9000 #-- Jul 10 18:59:43 --# nbasis: 5 
#> MCMC iteration 10000 #-- Jul 10 18:59:43 --# nbasis: 6
#> 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
  )

}

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.