#### worksheet to accompany Setting Up Your Simulator ## Jonathan Rougier. The latest version of this worksheet should ## be available at http://www.maths.bris.ac.uk/~mazjcr/SUYSworksheet.R ## The SUYS document should be available at SUYSdocument.pdf, and there ## should be slides at SUYSslides.pdf. Go to ## http://www.maths.bris.ac.uk/~mazjcr/#SUYS for more information. ## set 'makefiles <- FALSE' if you don't want files to appear ## in your working directory makefiles <- TRUE ## this is just to allow 'source("SUYSworksheet.R")' to work as expected ## by resetting the 'par' values if (!exists("op")) op <- par(no.readonly = TRUE) par(op) ## intercept calls to dev.print to insert additional information pdf.print <- function(fnm, width = 5.5, ...) { fnm <- gsub("\\.pdf$", "", fnm) # in case passed in as 'fnm.pdf' dev.print(pdf, paste0("SUYS_", fnm, ".pdf"), width = width, ...) } ## so now all pdf files will start SUYS_ #### Some useful functions #### ## If you are familiar with R you should look through these to ## make sure you understand how they work. Otherwise, come back ## and have a look at these later, but -- for now -- just paste ## them in and pick up the thread at the definition of the ## simulator. ## When you do look at them, think up and implement a little test to ## make sure they work as expected -- the same applies for all of ## the functions below. #### Here are some basic functions to move between a box and unit cube ## includes a little bit of error checking on the ranges, ## 'stopifnot' is a /very/ useful function! Make sure that the ## exit condition is easy to understand, eg ## 'Error: lower - eps < X is not TRUE' ## Note the inclusion of a little bit of fuzz on the limits: this is ## to allow for numerical issues when going forward and backward. ## eps = 1e-6 is usual, being a bit larger than the sqrt of the machine ## epsilon ## these are written to preserve colnames on 'X', if they exist ## also useful here to preserve class 'axial' for 'lohi' function, ## see below, 'axial' and 'lohi' functions box2cube <- function(X, lower, upper, check = TRUE, eps = 1e-6) { stopifnot(is.matrix(X)) is.axial <- inherits(X, "axial") p <- ncol(X) stopifnot(length(lower) == p, length(upper) == p) X <- t(X) if (check) stopifnot(lower - eps < X, X < upper + eps) robj <- t((X - lower) / (upper - lower)) if (is.axial) class(robj) <- c("matrix", "axial") robj } cube2box <- function(X, lower, upper, check = TRUE, eps = 1e-6) { is.axial <- inherits(X, "axial") stopifnot(is.matrix(X)) p <- ncol(X) stopifnot(length(lower) == p, length(upper) == p) X <- t(X) if (check) stopifnot(0 - eps < X, X < 1 + eps) robj <- t(X * (upper - lower) + lower) if (is.axial) class(robj) <- c("matrix", "axial") robj } #### sanity check runs #### function to construct mid-point and axial runs on [0, 1]^p ## have this classed as 'axial' to help 'lohi' function axial <- function(p, midp = TRUE) { blank <- matrix(0.5, 2, p) robj <- do.call("rbind", lapply(1:p, function(j) { blank[, j] <- c(0, 1) blank })) if (midp) { robj <- rbind(blank[1, ], robj) class(robj) <- c("matrix", "axial") } robj } #### and one to construct the lohi runs for one output ## y is a vector of outputs corresponding to the rows of X, which ## is a set of axial runs (with midpoint) ## preserve the colnames on 'X' if they exist lohi <- function(X, y) { stopifnot(inherits(X, "axial"), nrow(X) == length(y)) p <- ncol(X) robj <- sapply(1:p, function(j) { i <- c(1, 2 * j + c(0, 1)) x <- X[i, j]; y <- y[i] x[c(which.min(y), which.max(y))] }) if (!is.null(nms <- colnames(X))) colnames(robj) <- nms robj } #### Pick up the thread from here! #### #### Here is a toy simulator with inputs on [0, 1]^5 and 1 output ## can call with a 5-vector or with a matrix with 5 columns simulator <- function(theta) { if (is.vector(theta)) { if (length(theta) == 5) dim(theta) <- c(1, 5) else stop("Odd shape for \'theta\'") } stopifnot(is.matrix(theta), ncol(theta) == 5) robj <- 3 * sqrt(theta[, 1]) - (1 + 4 * (theta[, 1] - 0.5)^2) * theta[, 2] + theta[, 3] * (1 - theta[, 4]) * theta[, 5] round(robj, 2) } ## although this simulator has inputs on [0, 1]^5, I'll use the ## more general code, for inputs on [lower, upper]. Here are the ## original limits orig <- list(lower = rep(0, 5), upper = rep(1, 5)) ## later on, I'll hold the current limits in 'curr' #### OK, now we have a chance to try everything out. Here's where #### the experiment starts [section 'Sanity check' in the document] ## each time you create a new object, type its name to see what ## it looks like. Sometimes I do this explicitly with a 'print' ## command ## create an axial design on 5 inputs, give it some colnames ## I will keep X in its original units throughout, although this will ## not matter for the choice of orig above Xaxial <- cube2box(axial(5), orig$lower, orig$upper) colnames(Xaxial) <- paste0("th_", 1:ncol(Xaxial)) ## definitely do the following (see above): get into the habit :) print(Xaxial) # just 'Xaxial' at the prompt print(head(Xaxial)) # just 'head(Xaxial)' at the prompt ## run the simulator on the axial design yaxial <- simulator(Xaxial) print(cbind(Xaxial, y = yaxial)) # the 'y =' is just for the col name ## now add the lohi runs, use 'lohi*1*' for the 1st output (might add more ## outputs later) Xlohi1 <- lohi(Xaxial, yaxial) ylohi1 <- simulator(Xlohi1) print(cbind(Xlohi1, y = ylohi1)) ## go through each column of Xlohi to make sure you understand ## what is going on ... eg why is the lo value of th1 equal to 0 and the ## hi value of th1 equal to 1? ## OK -- that's the end of the sanity check, which is the first set ## of runs for your simulator. #### Before starting the first wave, load 'randtoolbox' for #### the 'sobol' function ## If you have not installed 'randtoolbox', then do # install.packages("randtoolbox", dependencies = TRUE) ## (or equivalent, eg if you are using a GUI or Rstudio), *just once*. ## To load 'randtoolbox' (ie to make its namespace available) do library(randtoolbox) # for the 'sobol' function, see '?sobol' ## if you get an error here, you need to go back to the install.packages #### fooling around with sobol sequences for fun if (FALSE) { # OK to skip this ... ## generate a 100-point design on 5 inputs X <- sobol(100, 5) ## visualise as a pairs plot pairs(X[1:25, ], pch = 16) # pause, examine plot, pairs(X, pch = 16) # compare to this one ## ** advanced users only ## install & load the rgl package, and use it to visualise the ## first three columns of X library(rgl) open3d() plot3d(X[, 1:3], type = "s", radius = 0.02, col = rainbow(nrow(X)), xlab = "X_1", ylab = "X_2", zlab = "X_3") ## also try plotting the first 25 on the same colour scheme plot3d(X[1:25, 1:3], type = "s", radius = 0.02, col = rainbow(nrow(X)), xlab = "X_1", ylab = "X_2", zlab = "X_3") ## filling a high-dimensional space is not easy! Try some permutations ## on the above, with different numbers of rows and columns. Also ## experiment with the 'scramble' argument in 'sobol', and with other ## types of quasi-random sequence, see '?sobol' for more details. ## end ** advanced users only } # end if (FALSE) #### start of the first wave (sanity check is wav0) #### Xwav0 <- rbind(Xaxial, Xlohi1) ywav0 <- c(yaxial, ylohi1) print(cbind(Xwav0, y = ywav0)) ## write this out to a tex file in tabular format if (makefiles) { tmp <- cbind(formatC(Xwav0, digits = 1, width = 3, flag = "-"), y = formatC(ywav0, format = "f", digits = 2, width = 5)) tmp[] <- paste0("$", tmp, "$") write.table(tmp, file = "SUYS_sanity1.tex", quote = FALSE, sep = " & ", eol = " \\\\\n", row.names = FALSE, col.names = FALSE) } ## if makefiles = TRUE, then you now have a file called SUYS_sanity1.tex ## in your working directory; use getwd() to see what that is. #### compute the pi & beta vectors ['Proceeding in waves' in the document] ## this is a procedure for figuring out which observables to use, and how ## many runs to do. You can just paste this in, or, if you have some R ## experience, you can see if you can figure out what is going on. The ## SUYSdocument.pdf will explain everything. I'll insert a commentary too. tmp <- box2cube(Xwav0, orig$lower, orig$upper) # onto [0,1]^p data <- as.data.frame(cbind(tmp, y = ywav0)) ppwav0 <- lm(y ~ ., data = data) # linear regression, used in 'predict' below print(summary(ppwav0)) ## ppwav0 is a linear model to predict y using the parameter values. Now ## make 'actual' observations ... yobs <- 1.84 # observation and sigma <- 0.1 # tolerance for this output ## SECRET! I am using 'yobs <- simulator(rep(0.7, 5))' piwav0 <- coef(ppwav0)[-1] / sigma # not interested in the intercept bewav0 <- abs(piwav0) print(round(bewav0, 1)) ## this shows that th_1 is an active parameter, and maybe th_2 as well ## now to compute the size of the NROY region. I use a scrambled ## Sobol sequence to do the integration on [0,1]^p numerically grid <- sobol(1e4, 5, scrambling = 1) colnames(grid) <- colnames(Xwav0) pp <- predict(ppwav0, newdata = as.data.frame(grid)) # predict inCtildeJ <- abs(yobs - pp) / sigma <= 3 vJ <- mean(inCtildeJ) print(c(vJ = vJ)) nJ <- 10 * ceiling(max(1 / vJ, 1 / (1 - vJ))) print(c(nJ = nJ)) ## that's the number of runs we want to do in the first wave, still ## on the original ranges. So here's the first wave ... Xwav1 <- cube2box(sobol(nJ, 5), orig$lower, orig$upper) colnames(Xwav1) <- colnames(Xwav0) # pass on column names print(head(Xwav1)) ywav1 <- simulator(Xwav1) # and run the simulator #### Now analyse all of the runs so far (wav0 and wav1) ## NROY = 'Not Ruled Out Yet': a parameter value for which the ## implausibility is <= 3 X <- rbind(Xwav0, Xwav1) y <- c(ywav0, ywav1) impl <- abs(yobs - y) / sigma # implausibility print(sum(impl <= 3)) # number of NROY in wav0 + wav1 print(sum(impl[1:nrow(Xwav0)] <= 3)) # number of NROY in wav0 ## time for a parallel coordinates plot. Can use 'MASS:::parcoord' ## but I am going to use a modified version for more control parcoord <- function (X, xlab = "", col = 1, lty = 1, lwd = 1, var.lab = colnames(X), var.lim = FALSE) { X <- as.matrix(X) n <- nrow(X); p <- ncol(X); var.lab <- var.lab rx <- apply(X, 2, range, na.rm = TRUE) X <- (t(X) - rx[1, ]) / (rx[2, ] - rx[1, ]) # now transposed on (0, 1) ## set up the plotting frame matplot(1:p, X, type = "n", xlab = xlab, ylab = "", axes = FALSE) # set up the plot if (!is.null(var.lab)) axis(1, at = 1:p, labels = var.lab, tick = FALSE) for (i in 1:p) { lines(c(i, i), c(0, 1), col = "grey70") if (var.lim) text(c(i, i), c(0, 1), labels = format(rx[, i], digits = 3), xpd = NA, offset = 0.3, pos = c(1, 3), cex = 0.9) } ## add the lines and return matlines(1:p, X, type = "l", col = col, lty = lty, lwd = lwd) invisible() } ## sort out the colours. I define a palette of 3 colours, and then ## assign each run a colour according to the breaks 0, 3, 6, Inf; using ## the cut function palette <- c("blue", "green", "yellow") breaks <- c(0, 3, 6, Inf) col <- palette[cut(impl, breaks, include.lowest = TRUE)] ## this is the basic plot parcoord(X, col = col, var.lim = TRUE) ## here is my souped-up version: see if you can spot the improvements! lty <- ifelse(1:nrow(X) <= nrow(Xwav0), 2, 1) lwd <- ifelse(col == "blue", 2, 1) oo <- order(-impl) labels <- expression(theta[1], theta[2], theta[3], theta[4], theta[5]) par(mar = c(3, 1, 1, 1)) parcoord(X[oo, ], col = col[oo], lty = lty[oo], lwd = lwd[oo], var.lim = TRUE, var.lab = labels) if (makefiles) pdf.print("PCP1", height = 4.5) ## see below (last gasp design) for specifying lower and upper ## limits for each axis ## If you wanted to, you could extend wav1 with, say, 20 more runs, using # Xwav12 <- cube2box(sobol(20, 5, init = FALSE), orig$lower, orig$upper) # colnames(Xwav12) <- colnames(Xwav1) # ywav12 <- simulator(Xwav12) # X <- rbind(X, Xwav12) # y <- c(y, ywav12) ## and then carry on as above. NB 'init = FALSE' in 'sobol' to continue ## the current sequence #### aggressively reduce the ranges on the active variables th_1 and th_2 keep <- impl <= 3 tmp <- apply(X[keep, ], 2, range) curr <- orig print(curr) curr$lower[1:2] <- tmp[1, 1:2] curr$upper[1:2] <- tmp[2, 1:2] print(curr) # these are the new limits ## let's compute the new volume oldv <- prod(orig$upper - orig$lower) newv <- prod(curr$upper - curr$lower) print(newv / oldv) #### now would be a good time to pause for a while ... ## these techniques are all embellishments on the procedure described ## above #### last gasp design #### ## this uses X and impl defined for wav0 + wav1 ## this is a quick Euclidean distance function for two matrices mydist <- function(X, Y) { stopifnot(is.matrix(X), is.matrix(Y), ncol(X) == ncol(Y)) m <- nrow(X); n <- nrow(Y) i <- rep(1L:m, n) j <- rep(1L:n, rep(m, n)) robj <- sqrt(rowSums((X[i, ] - Y[j, ])^2)) dim(robj) <- c(m, n) robj } ## I'll test this one myself! if (FALSE) { foo <- matrix(rnorm(101 * 5), ncol = 5) dd <- unname(as.matrix(dist(foo))[1:10, -(1:10)]) dd1 <- mydist(foo[1:10, ], foo[-(1:10), ]) all.equal.numeric(dd, dd1) # TRUE } # end if (FALSE) ## map X into the cube according to its own limits. Note that these ## distance calculations have to be done on standardised parameters ## to avoid unit effects rx <- apply(X, 2, range) # will use this again later Xtmp <- box2cube(X, rx[1, ], rx[2, ]) NROY <- impl <= 3 ## Here's a big space-filling design on [0, 1]^5. I use ## 'scrambling = 1' just to be sure that there is no duplication ## of the points in X; no harm in setting my own random seed # grid <- sobol(1e3, 5, scrambling = 1, seed = 1001) ## ## ## ## currently different from SUYSdocument.pdf => ## another option, which actually I prefer, on balance, is to infill ## the existing points. This needs a little function to put 'nfill' ## points between each pair of points in 'X' infill <- function(X, nfill = 1, nodup = TRUE) { stopifnot(is.matrix(X), length(nfill) == 1L, nfill >= 1) dfill <- 1 / (nfill + 1) fillseq <- seq(from = dfill, by = dfill, length.out = nfill) robj <- matrix(numeric(0), 0L, ncol(X)) if (!is.null(nms <- colnames(X))) colnames(robj) <- nms n <- nrow(X) if (n < 2L) return(robj) for (i in 1L:(n-1)) for (j in (i+1L):n) { dxy <- X[j, ] - X[i, ] robj <- rbind(robj, sweep(outer(fillseq, dxy, "*"), 2L, X[i, ], "+")) } ## remove duplicates if (nodup) { hash <- apply(rbind(X, robj), 1L, paste, collapse = ":") dup <- duplicated(hash)[-(1L:n)] robj <- robj[!dup, , drop=FALSE] } robj } grid <- infill(Xtmp, nfill = 1) # nfill = 2 would give many more points print(str(grid)) # but 2516 is probably enough points to work with ## now find out who the nearest point belongs to, and then thin ## the space-filling design accordingly ddin <- mydist(grid, Xtmp[NROY, ]) ddout <- mydist(grid, Xtmp[!NROY, ]) keep <- apply(ddin, 1, min) < apply(ddout, 1, min) gasp <- grid[keep, ] # still on the cube at this point print(str(gasp)) ## still too many rows: need a little function to find a good subset, ## eg maximin ## this will grind to a halt at dist() with too many rows (>1e4) ## but I doubt that >1e4 rows will be necessary for setting up maximinsubset <- function(X, n, rand = TRUE, reps = 1e2, verbose = TRUE) { stopifnot(is.matrix(X)) dd <- unname(as.matrix(dist(X))) ddmax <- max(dd) # largest interpoint distance diag(dd) <- Inf # take this out of play m <- nrow(X) if (m <= n) { warning("m <= n, returning X") return(X) } ## maintain the current best n-subset in curr (sorted) all <- 1L:m curr <- if (rand) sort(sample.int(m, n)) else 1L:n ddcurr <- dd[curr, curr, drop=FALSE] sccurr <- min(ddcurr) for (i in 1L:reps) { ## point to drop has smallest distance and smallest ## second smallest distance small <- apply(ddcurr, 1, function(x) sort(x)[1L:2]) drop <- order(small[1, ], small[2, ])[1] keep <- curr[-drop] ## point to add has biggest min-interpoint distance cand <- setdiff(all, curr) ddcand <- dd[keep, cand, drop=FALSE] add <- which.max(apply(ddcand, 2, min)) new <- sort(c(keep, cand[add])) ## check for early exit ddnew <- dd[new, new, drop=FALSE] scnew <- min(ddnew) if (scnew <= sccurr) break if (verbose) cat(sprintf("** iteration %3i, ppn max dist = %4.1f%%\n", i, 100 * scnew / ddmax)) curr <- new; ddcurr <- ddnew; sccurr <- scnew } attr(curr, "mindd") <- sccurr curr } ## here we go: get a space-filling subset by maximin set.seed(1001) # control the random seed ii <- maximinsubset(gasp, 30) # these are the best 30 subset indices gasp <- gasp[ii, ] # and these are the parameter values we will use ## back to original units, let's call this wav2 Xwav2 <- cube2box(gasp, rx[1, ], rx[2, ]) # using rx again ywav2 <- simulator(Xwav2) ## this next bit of code is similar to the above, except just wav2 impl <- abs(yobs - ywav2) / sigma print(sum(impl <= 3)) ## sweet! Lots of NROY to take forwards ... draw the PCP of wav2 col <- palette[cut(impl, breaks, include.lowest = TRUE)] lwd <- ifelse(col == "blue", 2, 1) oo <- order(-impl) ## specify lower and upper values for the PCP by adding dummy rows Xtmp <- rbind(Xwav2[oo, ], rx) col <- c(col[oo], NA, NA) lwd <- c(lwd[oo], NA, NA) parcoord(Xtmp, col = col, lwd = lwd, var.lim = TRUE, var.lab = labels) if (makefiles) pdf.print("PCP2", height = 4.5) ## draw the PCP of the whole lot as well, to see if any further adjustment ## in the ranges is possible X <- rbind(X, Xwav2) y <- c(y, ywav2) impl <- abs(yobs - y) / sigma col <- palette[cut(impl, breaks, include.lowest = TRUE)] lwd <- ifelse(col == "blue", 2, 1) oo <- order(-impl) parcoord(X[oo, ], col = col[oo], lwd = lwd[oo], var.lim = TRUE, var.lab = labels) if (makefiles) pdf.print("PCP3", height = 4.5) ## SECRET! Let's put the 'true' value on as well abline(h = 0.7, col = "red") ## adjustment to be made to the ranges? print(curr) rx <- apply(X[impl <= 3, ], 2, range) curr$lower[1:2] <- rx[1, 1:2] curr$upper[1:2] <- rx[2, 1:2] print(curr) ## No. #### Nudging to increase variation in the mean ## a couple of shorthands for matrix multiplication "%t*%" <- function(x, y) crossprod(x, y) "%*t%" <- function(x, y) tcrossprod(x, y) ## this function nudges X so that A %*% t(X) = t(Z) condK <- function(U, A, Z) { stopifnot(is.matrix(U), is.matrix(A), is.matrix(Z)) n <- nrow(U) p <- ncol(U) k <- nrow(A) stopifnot(ncol(A) == p, nrow(Z) == n, ncol(Z) == k) stopifnot(0 <= U & U <= 1, 0 <= Z & Z <= 1) U <- qnorm(t(U)); Z <- qnorm(t(Z)) Z <- A %*% U - Z U <- U - A %t*% solve(A %*t% A, Z) # NOT the right way! t(pnorm(U)) } ## do a couple of pairs plots p <- 20 n <- 100 X <- sobol(n, p) colnames(X) <- paste0("th_", 1:ncol(X)) A <- matrix(1 / p, 1, p) # just consider the overall mean Z <- matrix(sobol(n, 1)) Xp <- condK(X, A, Z) pairs(cbind(X[, 17:20], Mean = rowMeans(X)), pch = 16, cex = 0.8, xlim = c(0, 1), ylim = c(0, 1)) pairs(cbind(Xp[, 17:20], Mean = rowMeans(Xp)), pch = 16, cex = 0.8, xlim = c(0, 1), ylim = c(0, 1)) ## by all means create a better set of pair plots (see Fig 4 and 5), ## here is my plotting function mypairs <- function(X, xlim, ylim, offset = 0) { X <- as.matrix(X) n <- ncol(X) op <- par(mfrow = c(n, n), mar = c(1, 1, 1, 0.1), cex = 0.75, mgp = c(2, 0.7, 0), xpd = NA, las = 1); on.exit(par(op)) xax <- c(0, 0.5, 1) for (i in 1:n) for (j in 1:n) { if (j < i) { plot.new() next } if (i == j) { plot.new() plot.window(xlim = c(0, 1), ylim = c(0, 1), asp = 1) if (i < n) nm <- substitute(theta[nn], list(nn = offset + i)) else nm <- "Mean" text(0.5, 0.5, nm, cex = 1.25) rect(0, 0, 1, 1) next } plot.new() plot.window(xlim = c(0, 1), ylim = c(0, 1), asp = 1) points(X[, c(j, i)], pch = 16) axis(1, xax, pos = 0); axis(2, xax, pos = 0); rect(0, 0, 1, 1) } invisible(NULL) } mypairs(cbind(X[, 17:20], Mean = rowMeans(X)), offset = 16) if (makefiles) pdf.print("trans1", width = 6) mypairs(cbind(Xp[, 17:20], mean = rowMeans(Xp)), offset = 16) if (makefiles) pdf.print("trans2", width = 6) #### If you get this far and would like a bigger challenge ... ## Suppose that another output has y2 <- theta_3 * sqrt(abs(y1)), ## for which yobs2 <- 0.95 and sigma2 <- 0.075. Let y1 be the first ## output (also yobs1, sigma1). ## 1. Using only the wav0 runs, find the pi tableau, a 2 by 5 matrix of ## pi values for (y1, y2) by (th_1, ..., th_5). This is used to ## identify good choices of outputs in complicated experiments where ## there are lots of outputs. ## 2. Starting with the new reduced ranges (based on wav0, wav1, and y1), ## do another experiment using y2 to reduce the ranges further. ## (a) find vJ and nJ (remember to remove legacy runs) ## (b) run the simulator nJ times, compute y2 ## (c) compute the implausibilities based on y2 and redraw the PCP ## (d) update the ranges and report the new volume