In OakleyJ/MUCM: Gaussian process emulator methods based on the MUCM toolkit

Training data

Inputs and ranges:

xr <- t(apply(fit$training.inputs, 2, range))
colnames(xr) <- c("minimum", "maximum")
xr <- signif(xr, 3)
knitr::kable(xr, padding = 0)
#stargazer::stargazer(xr, type = "html", header=FALSE)

Number of training runs: r fit$n.train

Prior mean function

Formula for the prior mean function:

r format(fit$formula)

Posterior mean of the coefficients $\hat{\beta}$:

stargazer::stargazer(fit$betahat, type = "html", header=FALSE)

Comparison of fitted values $h(x)^T\hat{\beta}$ with training outputs

f <- fit$H.training %*% fit$betahat
y <- fit$training.outputs
plot(y, f, pch = 16,
     xlab = "training output",
     ylab = "prior mean fitted value",
     xlim = range(f, y),
     ylim = range(f, y))
abline(0, 1, lty = 2)

Correlation function

# define correlation function used
if (isTRUE(all.equal(corMatern2.5, fit$cor.function))) {
    cor.funct <- "Matern 2.5 and the equation is " 
    equation <- paste0("\\[c(\\mathbf{x,x'}) = \\sigma^2 \\frac{1}{\\Gamma(\\nu) 2^{\\nu - 1}} \\left(\\sqrt{2\\nu}r \\right)^\\nu K_\\nu \\left(\\sqrt{2\\nu}r \\right) \\] \\n  where $\\Gamma$ is the gamma function, $K_\\nu$ is the modified Bessel function, where $\\nu = 2.5$ and $r$ is distance between inputs $x$ and $x'$, scaled by $\\mathbf{\\delta}$ as:
\\[ \\Big[ \\frac{(x_{1}-x'_{1})^2}{\\delta_1^2} +...+ \\frac{(x_{p}-x'_{p})^2}{\\delta_p^2}  \\Big]^{0.5}.\\]")

} else if (isTRUE(all.equal(corGaussian, fit$cor.function))) {
    cor.funct <- "Gaussian and the equation is " 
    equation <- paste0("\\[c(\\mathbf{x,x'}) = \\prod_{i=1}^{p} \\exp \\left[ -\\left( \\dfrac{\\mathbf{x}_i-\\mathbf{x}_i'}{\\mathbf{\\delta}_i} \\right) ^2\\right].\\]") 

} else if (isTRUE(all.equal(corGaussianPeriodic, fit$cor.function))) {
    cor.funct <- "Gaussian Periodic"
    equation <- NULL
} else { 
    cor.funct <- "user defined."
    equation <- NULL
}

• Correlation function:
  r paste(cor.funct, equation)

w <- apply(fit$training.inputs, 2, range)
delta <- exp(fit$phi.hat)
delta.scaled <- delta / (w[2, ] - w[1, ])
stargazer::stargazer(data.frame(delta = delta, 
                                delta.scaled = delta.scaled), type = "html", summary = FALSE)

knitr::kable(data.frame(data.frame(delta = delta, 
                                delta.scaled = delta.scaled)))

• $4\hat{\sigma} = r signif(4 * sqrt(fit$sigmasq.hat), 3)$
• output range: r signif(max(fit$training.outputs) - min(fit$training.outputs), 3)
• Posterior mean $\hat{\sigma}^2$

stargazer::stargazer(fit$sigmasq.hat,
                     type = "html", nobs = FALSE, header=FALSE, iqr = TRUE)

if (!is.null(fit$opt.convergence)) {
    if (fit$opt.convergence == 0) {
        conv <- "optimisation was successful"
    } else if (fit$opt.convergence == 1) {
        conv <- "optimisation reached iteration limit"
    } else 
        conv <- paste("optimisation may not have succeeded - Convergence code:", fit$opt.convergence, "(See ?optim for more details)")
    # conv <- cat("Convergence: ", conv, "\n")
    # if (!is.null(fit$opt.message)) {
    #     cat("Message:", fit$opt.message, "\n")
}

• The r fit$optim.method method was used for optimisation. The r conv.

Scatter plots

• Plots of output against each input, emulator mean and point-wise 95% intervals

emPlot1D(fit, input = 1:ncol(fit$training.inputs))

Cross validation analysis

• Leave one out cross validation (LOOCV) was carried out. The following is a plot of the predictions (calculated using LOOCV) against the true output.

CV.fit <- crossVal(fit)

+ The root mean square error was calculated as r format(CV.fit$mean.RMSE) and the normalised root mean square error was r format(CV.fit$mean.RMSE/(max(fit$training.outputs) - min(fit$training.outputs))).

Session Info

Please note that the MUCM package version used is r packageVersion("MUCM"). The package depends on the following packages. The version number indicates the version number of the package being used.

installedPackages <- installed.packages()
Package_Name <- unlist(strsplit(installedPackages[which(rownames(installedPackages)=="MUCM"), "Imports", drop = F],", ")[[1]])
version <- vector(mode = "character")
for (i in Package_Name){
    version <- c(version, paste(packageVersion(i)))
}
dependencies <- data.frame("Package Name" = Package_Name, "Package Version" = version)

stargazer::stargazer(dependencies, summary = FALSE,
                     type = "html", title = "Dependencies on the MUCM package", header=FALSE)

OakleyJ/MUCM documentation built on May 7, 2019, 9:01 p.m.