### R code from vignette source 'advanced.Rnw' ################################################### ### code chunk number 1: advanced.Rnw:10-11 ################################################### options(continue=" ") ################################################### ### code chunk number 2: advanced.Rnw:70-72 ################################################### options(digits=3) ds = data.frame(x = rnorm(8), id=1:8) ################################################### ### code chunk number 3: advanced.Rnw:96-99 ################################################### options(digits=3) ds = ds[order(ds$x),] ds ################################################### ### code chunk number 4: advanced.Rnw:161-165 ################################################### diffx = diff(ds$x) min(diffx) with(ds, id[which.min(diffx)]) # first val with(ds, id[which.min(diffx) + 1]) # second val ################################################### ### code chunk number 5: advanced.Rnw:240-246 ################################################### p = .78 + (3 * 1:7)/100 allprobs = matrix(nrow=length(p), ncol=11) for (i in 1:length(p)) { allprobs[i,] = round(dbinom(0:10, 10, p[i]),2) } table = cbind(p, allprobs) ################################################### ### code chunk number 6: advanced.Rnw:253-254 ################################################### table ################################################### ### code chunk number 7: advanced.Rnw:306-314 ################################################### runave = function(n, gendist, ...) { x = gendist(n, ...) avex = numeric(n) for (k in 1:n) { avex[k] = mean(x[1:k]) } return(data.frame(x, avex)) } ################################################### ### code chunk number 8: advanced.Rnw:363-367 ################################################### vals = 1000 set.seed(1984) cauchy = runave(vals, rcauchy) t4 = runave(vals, rt, 4) ################################################### ### code chunk number 9: runaver (eval = FALSE) ################################################### ## plot(c(cauchy$avex, t4$avex), xlim=c(1, vals), type="n") ## lines(1:vals, cauchy$avex, lty=1, lwd=2) ## lines(1:vals, t4$avex, lty=2, lwd=2) ## abline(0, 0) ## legend(vals*.6, -1, legend=c("Cauchy", "t with 4 df"), ## lwd=2, lty=c(1, 2)) ################################################### ### code chunk number 10: advanced.Rnw:443-444 ################################################### plot(c(cauchy$avex, t4$avex), xlim=c(1, vals), type="n") lines(1:vals, cauchy$avex, lty=1, lwd=2) lines(1:vals, t4$avex, lty=2, lwd=2) abline(0, 0) legend(vals*.6, -1, legend=c("Cauchy", "t with 4 df"), lwd=2, lty=c(1, 2)) ################################################### ### code chunk number 11: advanced.Rnw:485-493 ################################################### len = 10 fibvals = numeric(len) fibvals[1] = 1 fibvals[2] = 1 for (i in 3:len) { fibvals[i] = fibvals[i-1] + fibvals[i-2] } fibvals ################################################### ### code chunk number 12: advanced.Rnw:629-640 ################################################### # read in the data input = readLines("c:/book/co25_d00.dat", n=-1) # figure out how many counties, and how many entries num = length(grep("END", input)) allvals = length(input) numentries = allvals-num # create vectors to store data county = numeric(numentries); lat = numeric(numentries) long = numeric(numentries) ################################################### ### code chunk number 13: advanced.Rnw:657-678 ################################################### curval = 0 # number of counties seen so far # loop through each line for (i in 1:allvals) { if (input[i]=="END") { curval = curval + 1 } else { # remove extraneous spaces nospace = gsub("[ ]+", " ", input[i]) # remove space in first column nospace = gsub("^ ", "", nospace) splitstring = as.numeric(strsplit(nospace, " ")[[1]]) len = length(splitstring) if (len==3) { # new county curcounty = splitstring[1]; county[i-curval] = curcounty lat[i-curval] = splitstring[2]; long[i-curval] = splitstring[3] } else if (len==2) { # continue current county county[i-curval] = curcounty; lat[i-curval] = splitstring[1] long[i-curval] = splitstring[2] } } } ################################################### ### code chunk number 14: advanced.Rnw:702-707 ################################################### # read county names countynames = read.table("c:/book/co25_d00a.dat", header=FALSE) names(countynames) = c("county", "countyname") ################################################### ### code chunk number 15: advanced.Rnw:790-801 (eval = FALSE) ################################################### ## xvals = c(min(lat), max(lat)) ## yvals = c(range(long)) ## plot(xvals, yvals, pch=" ", xlab="", ylab="", xaxt="n", yaxt="n") ## counties = unique(county) ## for (i in 1:length(counties)) { ## # first element is an internal point ## polygon(lat[county==counties[i]][-1], long[county==counties[i]][-1]) ## # plot name of county using internal point ## text(lat[county==counties[i]][1], long[county==counties[i]][1], ## countynames$countyname[i]) ## } ################################################### ### code chunk number 16: advanced.Rnw:851-852 ################################################### options(digits=4) ################################################### ### code chunk number 17: advanced.Rnw:858-862 ################################################### library(ggmap) options(digits=4) amherst = c(lon=-72.52, lat=42.36) mymap = get_map(location=amherst, zoom=13, color="bw") ################################################### ### code chunk number 18: bikeplotmap ################################################### myride = read.csv("http://www.amherst.edu/~nhorton/sasr2/datasets/cycle.csv") head(myride, 2) ggmap(mymap) + geom_point(aes(x=Longitude, y=Latitude), data=myride) ################################################### ### code chunk number 19: advanced.Rnw:938-939 ################################################### options(digits=3) ################################################### ### code chunk number 20: advanced.Rnw:967-972 ################################################### library(foreign) USArrests.st = transform(USArrests, region = tolower(rownames(USArrests))) write.dta(USArrests.st, "c:\\book\\USArrests.dta", convert.factors="string") ################################################### ### code chunk number 21: choro ################################################### library(ggmap) USArrests.st = transform(USArrests.st, murder = cut_number(Murder, 5)) us_state_map = map_data('state') map_data = merge(USArrests.st, us_state_map, by="region") map_data = map_data[order(map_data$order),] p0 = ggplot(map_data, aes(x=long, y=lat, group=group)) + geom_polygon(aes(fill = murder)) + geom_path(colour='black') + theme(legend.position = "bottom", panel.background=element_rect(fill="transparent", color=NA)) + scale_fill_grey(start=1, end =.1) + coord_map(); plot(p0) ################################################### ### code chunk number 22: advanced.Rnw:1172-1188 (eval = FALSE) ################################################### ## # grab contents of web page ## urlcontents = readLines("http://tinyurl.com/cartoonguide") ## ## # find line with sales rank ## linenum = suppressWarnings(grep("See Top 100 in Books", urlcontents)) ## ## # split line into multiple elements ## linevals = strsplit(urlcontents[linenum], ' ')[[1]] ## ## # find element with sales rank number ## entry = grep("#", linevals) ## charrank = linevals[entry] # snag that entry ## charrank = substr(charrank, 2, nchar(charrank)) # kill '#' at start ## charrank = gsub(',' ,'', charrank) # remove commas ## salesrank = as.numeric(charrank) # make it numeric ## cat("salesrank=", salesrank, "\n") ################################################### ### code chunk number 23: advanced.Rnw:1316-1332 ################################################### library(RCurl) myurl = getURL("https://www3.amherst.edu/~nhorton/sasr2/datasets/cartoon.txt", ssl.verifypeer=FALSE) file = readLines(textConnection(myurl)) n = length(file)/2 rank = numeric(n) timeval = as.POSIXlt(rank, origin="1960-01-01") for (i in 1:n) { timeval[i] = as.POSIXlt(gsub('EST', '', gsub('EDT', '', file[(i-1)*2+1])), tz="EST5EDT", format="%a %b %d %H:%M:%S %Y") rank[i] = as.numeric(gsub('NA', '', gsub('salesrank= ','', file[i*2]))) } timerank = data.frame(timeval, rank) ################################################### ### code chunk number 24: advanced.Rnw:1343-1344 ################################################### head(timerank, 4) ################################################### ### code chunk number 25: advanced.Rnw:1398-1402 ################################################### library(lubridate) timeofday = hour(timeval) night = rep(0,length(timeofday)) # vector of zeroes night[timeofday < 8 | timeofday > 18] = 1 ################################################### ### code chunk number 26: cartoonplot (eval = FALSE) ################################################### ## plot(timeval, rank, type="l", xlab="", ylab="Amazon Sales Rank") ## points(timeval[night==1], rank[night==1], pch=20, col="black") ## points(timeval[night==0], rank[night==0], pch=20, col="red") ## legend(as.POSIXlt("2013-10-03 00:00:00 EDT"), 6000, ## legend=c("day","night"), col=c("red","black"), pch=c(20,20)) ## abline(v=as.numeric(as.POSIXlt("2013-10-01 00:00:00 EST")), lty=2) ################################################### ### code chunk number 27: advanced.Rnw:1442-1443 ################################################### plot(timeval, rank, type="l", xlab="", ylab="Amazon Sales Rank") points(timeval[night==1], rank[night==1], pch=20, col="black") points(timeval[night==0], rank[night==0], pch=20, col="red") legend(as.POSIXlt("2013-10-03 00:00:00 EDT"), 6000, legend=c("day","night"), col=c("red","black"), pch=c(20,20)) abline(v=as.numeric(as.POSIXlt("2013-10-01 00:00:00 EST")), lty=2) ################################################### ### code chunk number 28: advanced.Rnw:1555-1566 ################################################### truerand = function(numrand) { read.table(as.character(paste("http://www.random.org/integers/?num=", numrand, "&min=0&max=1000000000&col=1&base=10&format=plain&rnd=new", sep="")))/1000000000 } quotacheck = function() { line = as.numeric(readLines( "http://www.random.org/quota/?format=plain")) return(line) } ################################################### ### code chunk number 29: advanced.Rnw:1571-1573 ################################################### truerand(7) quotacheck() ################################################### ### code chunk number 30: advanced.Rnw:1806-1812 ################################################### # Define constants and useful functions weight = c(0.3, 0.2, 2.0) volume = c(2.5, 1.5, 0.2) value = c(3000, 1800, 2500) maxwt = 25 maxvol = 25 ################################################### ### code chunk number 31: advanced.Rnw:1827-1851 ################################################### # minimize the grid points we need to calculate max.items = floor(pmin(maxwt/weight, maxvol/volume)) # useful functions getvalue = function(n) sum(n*value) getweight = function(n) sum(n*weight) getvolume = function(n) sum(n*volume) # main function: return 0 if constraints not met, # otherwise return the value of the contents, and their weight findvalue = function(x) { thisweight = apply(x, 1, getweight) thisvolume = apply(x, 1, getvolume) fits = (thisweight <= maxwt) & (thisvolume <= maxvol) vals = apply(x, 1, getvalue) return(data.frame(I=x[,1], II=x[,2], III=x[,3], value=fits*vals, weight=thisweight, volume=thisvolume)) } # Find and evaluate all possible combinations combs = expand.grid(lapply(max.items, function(n) seq.int(0, n))) values = findvalue(combs) ################################################### ### code chunk number 32: advanced.Rnw:1859-1861 ################################################### max(values$value) values[values$value==max(values$value),]