--- title: "Homework 2 Solutions" subtitle: "Advanced Statistical Computing (STAT 6984)" author: "Robert B. Gramacy ( : )
Department of Statistics, Virginia Tech" output: html_document --- {r, echo=FALSE} options(width=80) library(knitr) knit_hooks$set(no.main = function(before, options, envir) { if (before) par(mar = c(4.1, 4.1, 1.1, 1.1)) # smaller margin on top })  ### Problem 2: Immutable objects (30 pts) *Consider the following function which swaps elements i and j of a vector v.* {r} swap <- function(v, i, j) { tmp <- v[i] v[i] <- v[j] v[j] <- tmp }  *A disadvantage of this implementation is that it copies the entire vector, v, in order to work with just two of its elements. Consider the following example.* {r} v <- 1:1000000000 system.time(swap(v, i=1, j=2))  - *Nearly a second to "touch" three numbers isn't super speedy.* *Your tasks on this problem are the following.* a. *Report on how much time it takes to swap two elements (i=1; j=2) directly on the command line, i.e., without wrapping in a function.* Here is the setup of such an expression, which is quoted below so it can be evaluated for timing purposes later. {r} i <- 1 j <- 2 expr <- quote({tmp <- v[1]; v[1] <- v[2]; v[2] <- tmp})  Here are the times from five consecutive evaluations. {r} rbind(system.time(eval(expr)), system.time(eval(expr)), system.time(eval(expr)))  - Notice how the first evaluation takes a substantial amount of time, but subsequent evaluations are essentially instantaneous. - The same thing does not happen with consecutive evaluations of the original swap function. {r} rbind(system.time(swap(v, i=1, j=2)), system.time(swap(v, i=1, j=2)), system.time(swap(v, i=1, j=2)))  - Each evaluation consistently takes about 8/10 of a second. b. *Write a new version of the swap function, called swap.eval, which uses quote and eval to perform the calculation just like in part a. but within the function environment and without copying v by working on v in the parent.frame. Although this is a toy example, a similar code might be useful if, say, indicies i and j required substantial pre-calculation within the function before the swap occurred. Demonstrate your swap.eval with i=1; j=2 and report on the time.* The function below takes the same idea behind the quoted function above, but puts it inside a function. The eval call puts the scope of evaluation back in the parent frame so that v need not be copied in at great expense. Note however that all variables must come from the parent frame, including i and j. {r} swap.eval <- function() { e <- quote({tmp <- v[i]; v[i] <- v[j]; v[j] <- tmp}) eval(e, envir=parent.frame()) }  Now we can try the swap. {r} rbind(system.time(swap.eval()), system.time(swap.eval()), system.time(swap.eval()))  - Very similar behavior to part a.: slow the first time, but essentially instantaneous afterwards. c. *Write a similar function named swap.do which can be called via do.call that similarly accesses v in the parent frame. Add a print statement at the end of swap.do to show the first five elements of v after the swap occurs. Demonstrate swap.do with i=1; j=2 and report on the time. Are there any disadvantages to swap.do compared to swap.eval?* The function below suites our specifications. However, I have omitted the print statement for now. {r} swap.do <- function(i, j) { tmp <- v[i] v[i] <- v[j] v[j] <- tmp }  Lets call it three times and collect timings, just as above. {r} expr <- quote(do.call("swap.do", args=list(i=1, j=2), quote=TRUE)) rbind(system.time(eval(expr)), system.time(eval(expr)), system.time(eval(expr)))  - So the behavior is quite similar to the original swap function, perhaps ever-so-slightly faster. - That's one disadvantage: not nearly as fast as swap.eval. Now lets add the print statement to explore other behavior. {r} swap.do <- function(i, j) { tmp <- v[i] v[i] <- v[j] v[j] <- tmp print(v[1:5]) }  Consider the state of v before and after the swap.do call ... Before: {r} print(v[1:5])  During: {r} do.call("swap.do", args=list(i=1, j=2), quote=TRUE)  - Notice that v has changed *within* the function. After: {r} v[1:5]  - Back the way it was at the start. - So the a change has not taken effect in the parent frame. That's another disadvantage. Consider the swap.eval alternative: {r} swap.eval() v[1:5]  - v is indeed changed! That seems more useful. ### Problem 3: Bisection broadening (30 pts) *Recall the bisection algorithm and S3 object-oriented implementation from class. The bisection method can be generalized to deal with the case$f(x_l) f(x_r) > 0$, by broadening the bracket. That is,* - *reduce$x_l$and/or increase$x_r, and try again.* - *A reasonable choice is to double the width of the interval, i.e.,* \begin{aligned} m &\leftarrow (x_l + x_r)/2, & w &\leftarrow x_r - x_l \\ x_l&\leftarrow m - w, & x_r &\leftarrow m + w. \end{aligned} *Your tasks are the following.* a. *Incorporate bracketing into the \R{bisection()} function we coded. Note that broadening is not guaranteed to findx_l$and$x_r$such that$f(x_l) f(x_r) \leq 0$, so you should include a limit on the number of times broadening is successively tried with a sensible default.* See the revised [bisection2.R](http://bobby.gramacy.com/teaching/asc/bisection2.R) file, which is a version of the original [bisection.R](http://bobby.gramacy.com/teaching/asc/bisection.R), augmented to include broadening. Some highlights of the modifications include - An if statement at the top of the bisection fucntion checking if broadening is required, and if so implementing at most 100 (by default) iterations of the algorithm described above via a new function called broaden. That function has a similar verbosity feature allowing progress to be printed to the screen in real time. - Changes to the S3 print and summary methods allow that progress to be visualized in text form after the fact. Actually, the summary method is unchanged, but since it calls print it benefits from changes therein to report broadening and bracketing iterations separately. - Changes to the plot S3 method similarly augment visualizations of progress. Broadening iterations are noted by a "+" mark above iteration numbers. To load those functions in we have to source the new R file. {r} source("bisection2.R")  b. *Use your modified function to find a root of the (original) function$f(x)$we used in class, but with a different starting interval of$x_l = 2$and$x_r = 3], i.e., not containing the root we found in class.* First define the function, as in class. {r} f <- function(x) log(x) - exp(-x)  Then call the function and demonstrate the print output. Here the default verb=0 is used so as not to clutter up the document, but feel free to try verb=1. {r} fr <- bisection(f, 2, 3) fr  - Observe how it separately reports on broadening and bracketing iterations. - Notice that this is the same answer we got with earlier with a less pathological initialization. The summary output contains more information. {r} summary(fr)  - Besides the updated print, the output is pretty similar to our earlier version. The three plot commands below, arranged horizontally, give a nice visualization of progress. {r, dev.args = list(bg = 'transparent'), fig.width=10, fig.height=4, fig.align="center"} par(mfrow=c(1,3)) plot(fr, main="bisection progress", ylab="f(x)", xlab="x") plot(fr, main="bisection progress [zoom]", after=20, ylab="f(x)", xlab="x") plot(fr, main="bisection progress [zoom]", after=25, ylab="f(x)", xlab="x")  - Observe how the first plot has indicies with a "+" above them, indicating broadening iterations. c. *Use your modified function find the root of* $$h(x) = (x - 1)^3 - 2x^2 + 10 - \sin(x),$$ starting withx_l = 1$and$x_r = 2$. First define the new h function. {r} h <- function(x) (x-1)^3 - 2*x^2 - sin(x)  Here are the results, skipping summary. {r} hr <- bisection(h, 1, 2) hr  - Looks like we found a root of r signif(hr$ans, 3) after r hr$broad and r nrow(hr$prog)-hrbroad bracketing iterations. The three plot commands below, arranged horizontally, give a nice visualization of progress. {r, dev.args = list(bg = 'transparent'), fig.width=10, fig.height=4, fig.align="center"} par(mfrow=c(1,3)) plot(hr, main="bisection progress", ylab="f(x)", xlab="x") plot(hr, main="bisection progress [zoom]", after=20, ylab="f(x)", xlab="x") plot(hr, main="bisection progress [zoom]", after=25, ylab="f(x)", xlab="x")  - We can see that substantial broadening was required. - And also notice how "close" we got to finding a root nearx=0.5\$, but the algorithm wasn't fooled!. ### Problem 4: R scripts from the Unix prompt (25 pts) The two interfaces I was thinking of are Rscript and R CMD BATCH. although R -e offers a more limited yet related capability. The main difference is that BATCH is somewhat more customizable in terms of how it handles input, writes output, etc. To keep it simple here, I will demonstrate Rscript. - I hope your own solution was somewhat more verbose than the above. But I didn't see any need to regurgitate the relevent help files here. I made a file called [rendpurl.R](http://bobby.gramacy.com/teaching/asc/rendpurl.R) which has the following two commands. {r} lines <- readLines("rendpurl.R") cat(paste(lines, collapse="\n"))  - Ignore the ## above, which comes from Rmarkdown, not from the R script. - render creates the .html file and purl creates the .R version. My build.sh script has the following lines. {r} lines <- readLines("build.sh") cat(paste(lines, collapse="\n"))  Since the rendpurl.R` script works on this very Rmarkdown document, running it here would create an infinite loop, so lets not do that. Simply download the files if you want to try them.