################################################### ### chunk number 1: ################################################### options(width=50) noPrompt <- function() { options(prompt=" ") options(continue=" ") } doPrompt <- function() { options(prompt="> ") options(continue="+ ") } ################################################### ### chunk number 2: ################################################### noPrompt() ################################################### ### chunk number 3: ################################################### 2 + 2 ################################################### ### chunk number 4: ################################################### ## a possible learning curve y <- 1 x = c( "commands" = y, "vectorization"= (y <- y+1), "functions" = (y <- y+3), "Apply functions" = (y <- y+3), "default arguments" = (y <- y+1), "generic functions" = (y <- y+3), "S3 methods" = (y <- y+2), "S4 methods" = (y <- y+2), "scope" = (y <- y+4), "computing on the language" = (y <- y+10) ) out <- spline(x, 1:length(x)) plot(out, main="learning 'curve' -- function viewpoint", xlab="t", ylab = "useR <--> programmerR <--> wizaRd", type="l", col="gray70", lwd=5) y <- c(0,x) for(i in seq_along(x)) text(y[i-1] + 3,i,names(x)[i], pos=1) ################################################### ### chunk number 5: ################################################### x <- rnorm(100) prob <- ((1:100) - 0.5)/100 qp <- qnorm(prob) plot(sort(x), qp, pch=1, xlab="sorted data", ylab="Standard normal quantiles") ################################################### ### chunk number 6: ################################################### x <- runif(100) prob <- ((1:100) - 0.5)/100 qp <- qnorm(prob) plot(sort(x), qp, pch=1, xlab="sorted data", ylab="Standard normal quantiles") ################################################### ### chunk number 7: ################################################### ourqqplot <- function(x, xlab="sorted data", ylab = "standard normal quantiles", ...) { n <- length(x) prob <- seq(0,1,length=(n+2))[c(-1,-(n+2))] qp <- qnorm(prob) plot(sort(x), qp, pch=1, xlab=xlab, ylab=ylab,...) } ourqqplot(rnorm(100)) # call once ourqqplot(runif(100)) # call twice ################################################### ### chunk number 8: ################################################### f <- function(x) x^2 f ################################################### ### chunk number 9: ################################################### rm(f) ## clear out "f" source("f.R") ## read in f ################################################### ### chunk number 10: ################################################### f <- function(x) x^2 ################################################### ### chunk number 11: ################################################### f <- function(x) {x <- x[!is.na(x)]; mean(x)} f(c(1,1,NA,3,4)) ################################################### ### chunk number 12: ################################################### require(lattice) f <- function(x) print(bwplot(x)) # lattice issue f(rt(100, df=3)) ################################################### ### chunk number 13: ################################################### args(rnorm) ################################################### ### chunk number 14: ################################################### rnorm(n = 2, sd = 10) ################################################### ### chunk number 15: ################################################### rnorm(n = 2, me = 2) ################################################### ### chunk number 16: ################################################### rnorm(2,2) ################################################### ### chunk number 17: ################################################### args(stats:::t.test.default) ################################################### ### chunk number 18: eval=FALSE ################################################### ## alternative <- match.arg(alternative) ################################################### ### chunk number 19: ################################################### f <- function(x, y) { if(missing(y)) boxplot(x) else plot(x,y) } f(rt(100, df = 1)) ################################################### ### chunk number 20: ################################################### args(plot.default) ################################################### ### chunk number 21: ################################################### f <- function(x, label = deparse(x), doLabel = FALSE) { if(doLabel) label # evaluates label x <- x + 1 print(label) } f(1) # does not evaluate label first, so label is 2 f(1, doLabel = TRUE) ################################################### ### chunk number 22: ################################################### plotl <- function(x,y,...) plot(x,y,type="l",...) x <- seq(0,2*pi,length=100) plotl(x,sin(x), col="blue") # (issue with type="p", say) ################################################### ### chunk number 23: ################################################### args(t.test) args(stats:::t.test.default) ################################################### ### chunk number 24: ################################################### noCol <- function(x,y,...) { ## stop col argument a <- list(...) a$col <- NULL; a$xlab="x"; a$ylab="y" a$x <- x; a$y <- y do.call("plot",a) } x <- seq(0,2*pi,length=100); noCol(x,sin(x), col="blue") ################################################### ### chunk number 25: ################################################### m <- matrix(rnorm(10*15), nrow=10); ret <- 0 for(i in 1:10) { # nested for loops for(j in 1:15) { res <- ret + m[i,j] # err, sum(m) } } ################################################### ### chunk number 26: ################################################### findT <- function(x) mean(x)/sd(x)*sqrt(length(x)) ################################################### ### chunk number 27: ################################################### f1 <- function(m = 1000, n = 10) { res <- numeric(m) for(i in 1:m) { x <- rnorm(n) # data creation res[i] <- findT(x) # compute summary } res } system.time(f1()) # system.time() ################################################### ### chunk number 28: ################################################### f2 <- function(m = 1000, n = 10) { tmp <- matrix(rnorm(m*n), ncol=m) # data creation res <- apply(tmp, 2, FUN = findT) # summarize res } system.time(f2()) ################################################### ### chunk number 29: ################################################### m <- matrix(rnorm(10*100), ncol=100) res <- sapply(as.data.frame(m), findT) # apply to each col. ################################################### ### chunk number 30: ################################################### m <- 100; n <- 10 res <- sapply(1:m, function(i) findT(rnorm(n))) ################################################### ### chunk number 31: ################################################### require(MASS) with(Cars93, tapply(MPG.city, Cylinders, mean)) ################################################### ### chunk number 32: ################################################### mapply("rnorm", 1,2^(1:4),1) ################################################### ### chunk number 33: ################################################### require(doBy) summaryBy(MPG.highway ~ Cylinders, data=Cars93, FUN = function(x,...) { c(xbar=mean(x,...), n=length(x)) }, na.rm=TRUE) # not passed to length ################################################### ### chunk number 34: ################################################### d <- rbind("array"=paste("a",c("a","d","l","_"),"ply", sep=""), "data.frame"=paste("d",c("a","d","l","_"),"ply", sep=""), "list"=paste("l",c("a","d","l","_"),"ply", sep="")) dimnames(d) <- list("input" = c("array","data.frame","list"), "output" = c("array","data.frame","list","nothing")) ################################################### ### chunk number 35: ################################################### d ################################################### ### chunk number 36: ################################################### require(plyr) out <- ddply(baseball, variables. = c("year"), # as with tapply fun. = function(df) { # df -- extracted df ## mean if atleast 1 mean(subset(df, subset=X3b >= 1, select=X3b)) }) ################################################### ### chunk number 37: ################################################### require(ggplot2) print( with(out, qplot(year, X3b, main = "Whither the triple? Average triples per year for triples hitters.", ylab="Triples")) ) ################################################### ### chunk number 38: ################################################### require(plyr) ## note function returns data frame -- with doubles, triples, hr, hits out <- ddply(baseball, variables. = c("year"), # as with tapply fun. = function(df) { # df -- extracted df ## mean if atleast 1 hits <- subset(df, select=c("X2b","X3b","hr","h")) hits$singles <- with(hits, h - (X2b + X3b - hr)) mean(hits) # each column }) ################################################### ### chunk number 39: ################################################### d <- dim(out)[2] matplot(out$year, scale(out[,-1]), type="l", main = "average hits normalized") legend(1980,4,names(out)[2:d], col=1:(d-1), lty=1:(d-1)) ################################################### ### chunk number 40: ################################################### library(reshape) ## for melt tmp <- sapply(out,scale) # scale to compare, now a matrix tmp[,1] <- out[,"year"] # unscale year grph <- xyplot(value ~ year | variable, data = melt(as.data.frame(tmp), "year"), ## melt from reshape package type="l", main="Rise of HR?") print(grph) ################################################### ### chunk number 41: ################################################### 1/1+1 1/(1+1) ################################################### ### chunk number 42: ################################################### x <- 2 ################################################### ### chunk number 43: ################################################### unboundArg <- 1 f <- function(formalArg = 2, ...) { localArg <- 3 print(c(localArg, formalArg, unboundArg)) } f(0) ################################################### ### chunk number 44: ################################################### f1 <- function() { cat("env 1:");print(environment()) cat("parent 1:");print(parent.env(environment())) f2() } f2 <- function() { cat("env 2:");print(environment()) cat("parent 2:");print(parent.env(environment())) } f1() ## note both parents are global env. ################################################### ### chunk number 45: ################################################### x <- 5 sq <- function() { y <- x^2 y } myfunc <- function() { x <- 10 sq() } ################################################### ### chunk number 46: ################################################### myfunc() ################################################### ### chunk number 47: ################################################### x <- 5 myfunc <- function() { x <- 10 ## sq is now local to myfunc ## its parent environment is myfunc's -- x is bound to 10 sq <- function() { y <- x^2 y } sq() } myfunc() ## 100 = 10^2 ################################################### ### chunk number 48: ################################################### x <- "global" f <- function(x = "default", lazy.y = x) { print(x) x <- "local" print(lazy.y) ## changed binding in environment of function } f() ################################################### ### chunk number 49: ################################################### a <- 10 f <- function(a) { a <- 15 } c(f(a), a) # assignment did not change a ################################################### ### chunk number 50: ################################################### a <- 10 f <- function(a) assign("a", 15, envir=.GlobalEnv) c(f(a),a) ################################################### ### chunk number 51: ################################################### a <- 10 f <- function(a) a <<- 15 c(f(a),a) ################################################### ### chunk number 52: ################################################### make.negLL <- function(data, fixed = c(FALSE, FALSE)) { op <- fixed; force(data) # avoid change of data via lazy load function(p) { op[!fixed] <- p # adjust fixed parameters mu <- op[1]; sigma <- op[2] a <- -(1/2) * length(data) * log(2*pi*sigma^2) b <- -(1/2) * sum((data-mu)^2)/ sigma^2 -(a + b) }} # return a closure = fn + env. ################################################### ### chunk number 53: ################################################### set.seed(1) ################################################### ### chunk number 54: ################################################### normals <- rnorm(100, mean=1, sd=2) # artificial data c(xbar=mean(normals), s = sd(normals)) # summaries ## nLL <- make.negLL(normals) # nll is a function closure ls(environment(nLL)) # notice data is there ################################################### ### chunk number 55: ################################################### normals <- rnorm(100, mean = 100, sd=2) # modify normals in GlobalEnv ################################################### ### chunk number 56: ################################################### out <- optim(par=c(mu=0, sigma=1), fn = nLL, method="BFGS") out[['par']] ################################################### ### chunk number 57: ################################################### nLL1 <- make.negLL(normals, fixed=c(mu=100,FALSE)) out <- optim(c(sigma = 1), nLL1, method="BFGS") out[['par']] ################################################### ### chunk number 58: ################################################### data(CO2) # CO2 data class(CO2$Type) # factor class(CO2$conc) # numeric class(lm(uptake ~ conc, data=CO2)) # lm class(CO2$Plant) # length two! ################################################### ### chunk number 59: ################################################### par(mfrow=c(2,4)) plot(CO2$Type) # barplot plot(CO2$conc) # numeric plot(lm(uptake ~ conc, data=CO2)) # lm (4 plots) plot(CO2$Plant) # barplot ################################################### ### chunk number 60: ################################################### x <- "some string" class(x) attr(x,"class") # class is implicit class(x) <- c("String", class(x)) attr(x, "class") # explicit attribute ################################################### ### chunk number 61: ################################################### ## constructor of objects of class "String" - same name String <- function(x) { x <- paste(unlist(x), collapse="") class(x) <- c("String",class(x)) return(x) } String("bart") ################################################### ### chunk number 62: ################################################### print.String <- function(x) print(unclass(x)) String("Bart") ################################################### ### chunk number 63: ################################################### "+.String" <- function(e1,e2) { # Not Abelian! (x + y neq y + x) out <- paste(e1,e2, sep=" ") String(out) } x <- "Bart"; y <- "Simpson" (try(x+y)) class(x) <- "String" x + y ################################################### ### chunk number 64: extraction ################################################### "[.String" <- function(x,i,j,...,drop=TRUE) { tmp <- unlist(strsplit(x,"")) if(missing(i)) out <- tmp[] else out <- tmp[i] if(drop) String(out) else out } x[2:3] x[2:3, drop = FALSE] ################################################### ### chunk number 65: ################################################### squeeze <- function(x,size=20,...) UseMethod("squeeze") ################################################### ### chunk number 66: ################################################### squeeze.default <- function(x,size=20, ...) x # nothing ################################################### ### chunk number 67: ################################################### squeeze.String <- function(x, size=20, ...) { if(nchar(x) <= size) NextMethod(x) # example of next method ind <- which(" " == x[,,drop=FALSE]) if(length(ind) == 0) { if(nchar(x) <= size) return(x) else return(x[1:size] + "...") } else { i <- max(ind[which(ind <= size)]) return(x[1:i] + "...") }} ################################################### ### chunk number 68: ################################################### x <- String("The quick brown fox jumped over the lazy dog") squeeze(x) ################################################### ### chunk number 69: ################################################### setClass("Polygon", representation(x = "numeric", y = "numeric", gefgw = "numeric" # a tolerance ), prototype(x = numeric(0), y = numeric(0), gefgw = 0.0001) ) ################################################### ### chunk number 70: ################################################### QT <- setValidity("Polygon", function(object) { if(length(object@x) == length(object@y)) return(TRUE) else return("unequal length coordinates") }) ################################################### ### chunk number 71: ################################################### defPoly <- new("Polygon") ################################################### ### chunk number 72: ################################################### ## constructor Polygon <- function(x,y) { new("Polygon",x=x,y=y) } ##try it out triangle <- Polygon(x=c(0,2,1), y=c(0,0,1)) (try(Polygon(x = c(0,2,1), y = c(0)))) # different sizes ################################################### ### chunk number 73: ################################################### setAs("matrix","Polygon",function(from,to) { # coerce matrix if(nrow(from) > 1) new("Polygon",x=from[1,], y=from[2,]) else stop("Need 2 or more rows to define a polygon") }) setAs("data.frame","Polygon", function(from,to) { # data frame if(ncol(from) <= 1) stop("need two or more columns") if(! (is.numeric(from[,1]) && is.numeric(from[,2]))) stop("need first two columns to be numeric") as(t(as.matrix(from[,1:2])), "Polygon") }) ################################################### ### chunk number 74: ################################################### data(worldLL, package="PBSmapping") ## continent as polygons, 0 is -- you guess cont <- worldLL[worldLL[,1] == 0,3:4] class(cont) ## no as for PolySet. Could coerce, we promote: setOldClass("PolySet") # S3 class to S4 setIs("PolySet","data.frame") # force S3 inheritance p <- as(cont,"Polygon") ################################################### ### chunk number 75: ################################################### setMethod("show", signature = c(object="Polygon"), function(object) { cat(class(object),"object with\n") cat(" Num vertices:", length(object@x),"\n") invisible(object) }) ################################################### ### chunk number 76: ################################################### p # calls show(p) ################################################### ### chunk number 77: ################################################### QT <- setGeneric("print") # S3 generic, no def needed QT <- setMethod("print",signature(x="Polygon"), function(x, ...) { # match formal args print(rbind(x=x@x, y=x@y)) }) ################################################### ### chunk number 78: ################################################### print(triangle) ################################################### ### chunk number 79: ################################################### QT <- setGeneric("plot") QT <- setMethod("plot","Polygon",function(x,y,...) { ## use local functions to split up "..." localPlot <- function(...,col) plot(...) # pass col to polygon localPolygon <- function(...,main,sub,xlab,ylab,xlim,ylim) polygon(...) localPlot(NA,NA, xlim = range(x@x), ylim = range(x@y),...) localPolygon(x@x, x@y,...) }) ################################################### ### chunk number 80: ################################################### QT <- setGeneric("lines") QT <- setMethod("lines","Polygon", function(x,y,...) polygon(x@x,x@y, ...) ) ################################################### ### chunk number 81: ################################################### plot(p) # that's a continent? ################################################### ### chunk number 82: ################################################### ## new S4 generic -- use setGeneric with a function template QT <- setGeneric("rotate", def = function(object, theta) { standardGeneric("rotate") }) ## define a method to dispatch when generic is called QT <- setMethod("rotate", signature(object="Polygon",theta="numeric"), function(object,theta) { A <- rbind(c(cos(theta), -sin(theta)), c(sin(theta), cos(theta))) m0 = rbind(object@x, object@y) m1 = A %*% m0 return(new("Polygon",x = m1[1,], y = m1[2,])) }) ################################################### ### chunk number 83: ################################################### plot(rotate(p,pi), main="Australian view") ################################################### ### chunk number 84: ################################################### QT <- setGeneric("rotate<-", def = function(object, value) { standardGeneric("rotate<-") }) ## define a method to dispatch when generic is called QT <- setReplaceMethod("rotate", signature(object="Polygon",value="numeric"), function(object, value) { object <- rotate(object, value) return(object) }) ################################################### ### chunk number 85: ################################################### par(mfrow=c(1,2)) plot(triangle,main="triangle") rotate(triangle) <- pi/2 plot(triangle, main="triangle 2") ################################################### ### chunk number 86: ################################################### .findArea <- function(object) { x.mat <- cbind(object@x, object@y) if(nrow(x.mat) < 3) return(0); x.segmat <- cbind(x.mat, rbind(x.mat[2:nrow(x.mat), ], x.mat[1, ])); area <- abs(sum(x.segmat[,1] * x.segmat[,4] - x.segmat[,3] * x.segmat[,2])) / 2 return(area) } ################################################### ### chunk number 87: ################################################### ## find the area QT <- setGeneric("Area", function(object,...) standardGeneric("Area")) QT <- setMethod("Area","Polygon", function(object,...) { area <- .findArea(object) return(area) }) ################################################### ### chunk number 88: ################################################### library(diagram) m <- matrix(nrow=5, ncol=5,byrow=TRUE, data =c( ##a, b, c, d, e 0, 0, 0, 0, 0, #a 1, 0, 0, 0, 0, #b 1, 0, 0, 0, 0, #c 0, 0, 1, 0, 0, #d 0, 0, 1, 0, 0 #e )) nms <- letters[1:5] posm <- (cbind(c(2,1,3,2,4),c(1,2,2,3,3)) - 0.5)/4 plotmat(m, pos=posm, name=nms, box.type="square", box.prop=0.5, curve=0, main="map of a(b+c(d+e))") ################################################### ### chunk number 89: ################################################### library(diagram) m <- matrix(nrow=6, ncol=6,byrow=TRUE, data =c( ##sq re rh pa qu po 0, 1, 0, 0, 0, 0, #sq 0, 0, 0, 1, 0, 0, #re 0, 0, 0, 1, 0, 0, #rh 0, 0, 0, 0, 1, 0, #pa 0, 0, 0, 0, 0, 1, #qu 0, 0, 0, 0, 0, 0 #po )) nms <- c("square","rectangle","rhombus","parallellogram","quadrilateral","polygon") posm <- c(1,2,1,1,1) posm <- (cbind(c(4,4,0,2,2,2), c(4,3,3,2,1,0)) + .5)/5 plotmat(m, pos=posm,name=nms, box.type="square", box.prop=0.5, curve=0, main="4-sided polygons") ################################################### ### chunk number 90: ################################################### QT <- setClass("Quadrilateral", contains = "Polygon", ## subclass same slots, but stricter defn validity = function(object) { if(length(object@x) == 4) return(TRUE) else return("quadrilaterals have 4 sides, so need 4 vertices") }) ################################################### ### chunk number 91: ################################################### ##define some functions, but don't print them on screen .equalSideLengths <- function(object) { dxs <- diff(c(object@x, object@x[1])) dys <- diff(c(object@y, object@y[1])) if(abs(sum(c(dxs[1],dys[1])^2) - sum(c(dxs[2],dys[2])^2)) < object@gefgw) return(TRUE) else return("not all side lengths are equal") } .parallelSides <- function(object) { dxs <- diff(c(object@x, object@x[1])) dys <- diff(c(object@y, object@y[1])) if(det(rbind( c(dxs[1],dys[1]), c(dxs[3],dys[3]))) + det(rbind( c(dxs[2],dys[2]), c(dxs[4],dys[4]))) < object@gefgw) return(TRUE) else return("not all sides are parallel") } .equalAngles <- function(object) { dxs <- diff(c(object@x, object@x[1])) dys <- diff(c(object@y, object@y[1])) if(abs(sum(c(dxs[1],dys[1]) * c(dxs[2],dys[2])) + sum(c(dxs[3],dys[3]) * c(dxs[2],dys[2]))) < object@gefgw) return(TRUE) else return("not all angles are equal") } ################################################### ### chunk number 92: parallelogram ################################################### QT <- setClass("Parallelogram", contains = "Quadrilateral", validity = function(object) { .parallelSides(object) }) ################################################### ### chunk number 93: ################################################### pgram <- new("Parallelogram", x=c(0,1,2,1), y = c(0,0,1,1)) Area(pgram) plot(pgram, col="blue") ################################################### ### chunk number 94: rhombus ################################################### QT <- setClass("Rhombus", contains = "Parallelogram", validity = function(object) { .equalSideLengths(object) }) ################################################### ### chunk number 95: ################################################### rhombus <- new("Rhombus", x = c(-2,0,2,-0), y = c(0,1,0,-1)) plot(rhombus, col="red") ################################################### ### chunk number 96: rectangle ################################################### QT <- setClass("Rectangle",contains = "Parallelogram", validity = function(object) { .equalAngles(object) }) ################################################### ### chunk number 97: ################################################### rect <- new("Rectangle", x=c(0,2,2,0), y = c(0,0,1,1)) plot(rect, col="green") ################################################### ### chunk number 98: square ################################################### QT <- setClass("Square", contains = c("Rectangle"), validity = function(object) { .equalSideLengths(object) }) ################################################### ### chunk number 99: ################################################### sq <- new("Square", x = c(0,1,1,0), y = c(0,0,1,1)) plot(sq, col="magenta") ################################################### ### chunk number 100: ################################################### doPrompt()