############# # A real data example ############# library("SAVE") ############# # load data ############# data(spotweldfield, package = "SAVE") data(spotweldmodel, package = "SAVE") ############## # create the SAVE object which describes the problem and # compute the corresponding mle estimates ############## sw <- SAVE(response.name = "diameter", controllable.names = c("current", "load", "thickness"), calibration.names = "tuning", field.data = spotweldfield, model.data = spotweldmodel, mean.formula = as.formula(~1), bestguess = list(tuning = 4)) # summary of the results summary(sw) ############## # obtain the posterior distribution of the unknown parameters ############## set.seed(0) swbayes <- bayesfit(object = sw, prior = c(uniform("tuning", upper = 8, lower = 0.8)), n.iter = 20000, n.burnin = 100, n.thin = 2) # summary of the results summary(swbayes) # traceplots plot(swbayes, option = "trace", col = 2) # posterior of the calibration parameter (Fig 3) plot(swbayes, option = "calibration", col = 4, lty = 3) # posterior of the measurement error precision and bias function precision (Fig # 4) plot(swbayes, option = "precision", col = 2, lty = 4) ############## # validate the computer model at chosen set of controllable # inputs ############### load <- c(4, 5.3) curr <- seq(from = 20, to = 30, length = 20) g <- c(1, 2) xnew <- expand.grid(current = curr, load = load, thickness = g) set.seed(0) valsw <- validate(object = swbayes, newdesign = xnew, n.burnin = 100) # summary of results summary(valsw) # plot results (Fig 5) plot(valsw) # customized plot (Fig 6) av <- (valsw@validate)[, "pure.model"] tau <- (valsw@validate)[, "tau.pm"] par(oma = c(1, 1, 2, 0), mfrow = c(2, 2)) for (i in 1:2) { for (j in 1:2) { v <- ((i - 1) * 40 + (j - 1) * 20 + 1):((i - 1) * 40 + j * 20) plot(curr, av[v], type = "l", ylim = c(3, 9), xlab = "current", ylab = "weld diameter") lines(curr, av[v] + tau[v], lty = 3) lines(curr, av[v] - tau[v], lty = 3) text(22, 9, paste("thickness=", g[i], ", load=", load[j], sep = ""), cex = 0.8, pos = 1) # field data that correspond to this situation v <- ((i - 1) * 60 + (j - 1) * 30 + 1):((i - 1) * 60 + j * 30) data <- spotweldfield$diameter[v] inputs <- spotweldfield$current[v] points(inputs, data) } } mtext("Pure-model predictions", side = 3, outer = T, cex = 1.2) par(mfrow = c(1, 1)) ##### # another application - derivative wrt current ##### # design load <- 4 curr <- seq(from = 20, to = 30, length = 80) g <- 1 xnew <- expand.grid(current = curr, load = load, thickness = g) # Obtain samples set.seed(0) prdersw <- predictreality(object = swbayes, newdesign = xnew, tol = 1e-12) delta <- diff(curr)[1] model <- prdersw@modelpred dmodel <- diff(t(model))/delta bias <- prdersw@biaspred dbias <- diff(t(bias))/delta dreal <- dmodel + dbias # bias-corrected prediction dav <- apply(dreal, 1, mean) # tolerance bounds tmpdata <- matrix(dav, ncol = dim(dreal)[2], nrow = dim(dreal)[1], byrow = F) tmpdata <- dreal - tmpdata tmpdata <- abs(tmpdata) tau.real <- apply(tmpdata, 1, quantile, 0.9) # pure-model prediction construct design at which to emulate the model u <- 3.2 load <- 4 curr <- seq(from = 20, to = 30, length = 80) g <- 1 xnewpure <- expand.grid(current = curr, load = load, thickness = g, tuning = u) set.seed(0) pmsw <- predictcode(object = swbayes, newdesign = xnewpure, n.iter = 20000) # samples and estimate puremodel <- pmsw@samples dpuremodel <- diff(t(puremodel))/delta dpure <- apply(dpuremodel, 1, mean) # plot (Fig 7) plot(curr[-1], dav, ty = "l", ylim = c(min(dav - tau.real), max(dav + tau.real)), xlab = "current", ylab = "derivative") lines(curr[-1], dav + tau.real, lty = 2) lines(curr[-1], dav - tau.real, lty = 2) lines(curr[-1], dpure, lty = 4, lwd = 2) ############# # Simulated example ############# library("SAVE") ####### # load data ####### data(synthfield, package = "SAVE") data(synthmodel, package = "SAVE") ############## # create the SAVE object which describes the problem and # compute the corresponding mle estimates ############## synth <- SAVE(response.name = "y", controllable.names = "x", calibration.names = "v", field.data = synthfield, model.data = synthmodel, mean.formula = ~1 + x, bestguess = list(v = 1.5)) ############## # Bayesian fit ############## set.seed(0) synth <- bayesfit(object = synth, prior = uniform(var.name = "v", lower = 0, upper = 3), n.iter = 20000) summary(synth) # Fig 9 Histogram of the posterior for v plot(synth, option = "calibration") # Fig 10 Histogram of the precisions plot(synth, option = "precision") ############# # validate at a grid of x points ############## xnew <- data.frame(x = seq(from = 0.05, to = 3.05, length = 25)) valsynth <- validate(object = synth, newdesign = xnew, n.burnin = 100) # summary and plot of the validation process summary(valsynth) plot(valsynth) # Fig 11 A plot to compare the bias corrected prediction, the actual reality # curve the computer model evaluated at the posterior mean and at the least # square estimate par(mfrow = c(1, 1)) a <- 0 b <- 5 av.real <- (valsynth@validate)[, "bias.corrected"] delta <- (valsynth@validate)[, "tau.bc"] # Bias corrected prediction and tolerance bars plot(xnew$x, av.real, ty = "n", ylim = c(a, b), main = "Predictions", xlab = "x", ylab = "y") lines(xnew$x, av.real, lty = 1) lines(xnew$x, av.real + delta, lty = 2) lines(xnew$x, av.real - delta, lty = 2) points(synthfield$x, synthfield$y, pch = "*") # Representation of reality yR <- function(x) { 3.5 * exp(-1.7 * x) + 1.5 } lines(xnew$x, yR(xnew$x), col = 2) # Representation of Computer model at the least square estimate yM <- function(x) { 5 * exp(-x[1] * x[2]) } vls <- 0.63 tmp <- cbind(xnew$x, vls) lines(xnew$x, apply(tmp, 1, FUN = yM), col = 3) # Representation of Computer model at the posterior mean vmean <- mean(synth@mcmcsample[-(1:100), "v"]) tmp <- cbind(xnew$x, vmean) lines(xnew$x, apply(tmp, 1, yM), col = 4)