## code to accompany Lecture 7 of ## Statistical computing with R by ## Robert B. Gramacy ## so the code will wrap appropriately in the slides options(width=50) ## load good bisection code source("bisection.R") ## bisection method function WITHOUT SANITY CHECKS; ## output is a root of f in the interval (xl, xr) bisection2 <- function(f, xl, xr, tol=sqrt(.Machine\$double.eps), verb=0) { ## create an output object out <- list(f=f, tol=tol) class(out) <- "bisection" ## setup and check outputs fl <- f(xl) fr <- f(xr) out\$prog <- data.frame(xl=xl, xr=xr, fl=fl, fr=fr) if(fl == 0) { out\$ans <- xl; return(out) } else if(fr == 0) { out\$ans <- xr; return(out) } ## successively refine xl and xr n <- 1 while((xr - xl) > tol) { xm <- (xl + xr)/2 fm <- f(xm) if(fm == 0) { out\$ans <- xm; return(out) } else if (fl * fm < 0) { xr <- xm; fr <- fm } else { xl <- xm; fl <- fm } ## next iteration n <- n + 1 out\$prog[n,] <- c(xl, xr, fl, fr) if(verb > 0) cat("n=", n, ", (xl, xr)=(", xl, ", ", xr, ")\n", sep="") } ## return (approximate) root out\$ans <- (xl + xr)/2 return(out) } ## call the good one then the bad one f <- function(x) log(x) - exp(-x) ## results in an error fr <- bisection(f, 2, 1) ## results in nonsense output fr2 <- bisection2(f, 2, 1) fr2\$ans ## correct result fr3 <- bisection(f, 1, 2) fr3\$ans ## simple sanity check with stop doWork <- function(filename) { ## checking the type of file if(class(filename) != "character" || length(filename) != 1) stop("filenames should be a single character string") ## checking that the file exists if(! file.exists(filename)) stop("Could not open file: ", filename) ## otherwise do something with the file, e.g., read.delim(filename) } ## testing the function doWork("this") doWork(1) read.delim("this") read.delim(1) ## stopifnot x <- -2 stopifnot(x > 0) ## warnings if(c(TRUE, FALSE)) TRUE else FALSE x <- 1; attach(data.frame(x=2)) x ## suppressWarnings from regress.pls regress.pls ## illustrating try res <- try({ 1/2 }, silent=TRUE) res res <- try({ 1/"2" }, silent=TRUE) res ## a second chance res <- try({ 1/"2" }, silent=TRUE) if(class(res) == "try-error") { res <- 1/as.numeric("2") } res ## glmn.hr: ## ## estimate a glmnet model using X and y, then ## predict at XX and compare to yy. The penalty ## parameter is estimated first by cv.glmnet, ## then the 1-SE lambda is used to re-train the ## model. For safety, try is used in the re-estimation ## since the full data version may throw an error ## due to insufficient penalization glmn.hr <- function(X, y, XX, yy, nfolds=10) { ## estimate lambda via CV fit.cv <- cv.glmnet(X, y, family="binomial", nfolds=nfolds) lambda <- fit.cv\$lambda.1se ## be safe about re-estimating on full data while(1) { fit <- try(glmnet(X, y, family="binomial", lambda=lambda), silent=TRUE) if(class(fit)[1] == "try-error") { cat("glmnet error with lambda =", lambda, "\n") lambda <- 2*lambda } else break; } ## predict & calculate hits p <- predict(fit, newx=XX, s=lambda, type="response") ## calculate hits hit <- (round(p) == yy) return(list(hit=hit, lambda=lambda)) } ## illustrating try-catch tryCatch(1/"2", error=cat("wups\n"), warning=cat("oops\n"), finally=cat("done\n")) tryCatch(if(c(0,1)) 1, error=cat("wups\n"), warn=cat("oops\n"), finally=cat("done\n")) ## post-mortem debugging source("mind.R") ## finding minimum distance pair m <- rbind(c(0, 12, 5), c(12, 0, 8), c(5, 8, 0)) m ## try/not in slides mind(m) ## traceback can be helpful here traceback() ## adding print statmements imin <- function(x) { lx <- nrow(x) i <- x[lx] print(lx); print(i) j <- which.min(x[(i+1):(lx-1)]) return(c(j,x[j])) } ## re-run mind(m) ## working in the debugger instead source("mind.R") options(error=recover) mind(m) ## from here on out we are interactive in a way that ## can't be scripted ## fixing one error: imin <- function(x) { lx <- length(x) i <- x[lx] j <- which.min(x[(i+1):(lx-1)]) return(c(j,x[j])) } ## re-run mind(m) ## interactive again, return when we've decided to ## debug imin debug(imin) mind(m) ## another correction imin <- function(x) { lx <- length(x) i <- x[lx] j <- which.min(x[(i+1):(lx-1)]) k <- j + i return(c(k,x[k])) } mind(m) ## a fix to mind mind <- function(d) { n <- nrow(d) dd <- cbind(d,1:n) wmins <- apply(dd[-n,],1,imin) i <- which.min(wmins[2,]) j <- wmins[1,i] return(c(d[i,j],i,j)) } mind(m) ## slow code for computing powers of a matrix ## ## powers1: forms a matrix of powers of the ## vector x up to the degree dg powers1 <- function(x, dg) { pw <- matrix(x, nrow=length(x)) prod <- x ## current product for(i in 2:dg) { prod <- prod * x pw <- cbind(pw, prod) } return(pw) } ## running the code x <- runif(10000000) system.time(p1 <- powers1(x, 16)) ## now profiling Rprof() p1 <- powers1(x, 16) Rprof(NULL) summaryRprof() ## a faster powers function ## ## powers2: same as above but with pre-allocation powers2 <- function(x, dg) { pw <- matrix(nrow=length(x), ncol=dg) prod <- x ## current product pw[,1] <- prod for(i in 2:dg) { prod <- prod * x pw[,i] <- prod } } ## running the code system.time(p2 <- powers2(x, 16)) ## profiling to compare: try/not in slides Rprof() p2 <- powers2(x, 16) Rprof(NULL) summaryRprof() ## setting up for ar.gibbs source("ar.R") x <- ts(scan("ar3.txt", quiet=TRUE)) lx <- length(x) XKs <- make.XKs(x, 4) ## now for the profiling Rprof("ar3.Rprof") theta.3 <- ar.gibbs(100000, x[5:lx], Xk=XKs[[3]]) Rprof(NULL) ar3.Rprof <- summaryRprof("ar3.Rprof") ## try/not in slides ar3.Rprof\$by.total[1:10,] ar3.Rprof\$by.self[1:10,] ## lets try rmvnorm <- function(n,mu,sigma) { p <- length(mu) z <- matrix(rnorm(n * p),nrow=n) ch <- chol(sigma,pivot=T) piv <- attr(ch,"pivot") zz <- (z%*%ch) zzz <- 0*zz zzz[,piv] <- zz zzz + matrix(mu,nrow=n,ncol=p,byrow=T) } ## new profiling session Rprof("ar3.Rprof2") theta.3 <- ar.gibbs(100000, x[5:lx], Xk=XKs[[3]]) Rprof(NULL) ar3.Rprof2 <- summaryRprof("ar3.Rprof2") ## compare ar3.Rprof\$sampling.time ar3.Rprof2\$sampling.time ## Rprofmem is better in Rstudio Rprofmem()