# Functions for lab2 # by Tom Minka par(mfrow=c(2,1)) dev.copy = function() { # create a new device which is configured the same way as the old opar = par(no.readonly=T) genv = globalenv() eval(call(.Device),genv) par(opar) } ################### shuffle <- function(x) x[order(runif(length(x)))] # companions to xinch,yinch inchx = function(x,warn.log=T) { if (warn.log && par("xlog")) warning("x log scale: inchx() is non-sense") x / diff(par("usr")[1:2]) * par("pin")[1] } inchy = function (y, warn.log = TRUE) { if (warn.log && par("ylog")) warning("y log scale: inchy() is non-sense") y / diff(par("usr")[3:4]) * par("pin")[2] } find.gap <- function(h0,y,h,symmetric=T) { # returns the smallest position of a box of size h0 which will not collide # with boxes of size h centered on y. # sort y ord = order(y) y = y[ord] h = h[ord] n = length(y) # compute all gaps dy = diff(y)-(h[1:(n-1)]+h[2:n])/2 # dy must be positive since there are no collisions among y's i = which(dy > h0) if(length(i) == 0) { # only position is at ends y.pos = y[n]+(h[n]+h0)/2 y.neg = y[1]-(h[1]+h0)/2 if(symmetric && abs(y.neg) < abs(y.pos)) y.neg else y.pos } else { # this isn't quite optimal y.new = y[i]+(h[i]+h0)/2 y.new[which.min(abs(y.new))] } } move.collisions.vertical <- function(x,w,h,symmetric=T) { # find the minimal vertical shifts so that boxes of width w and height h # centered on x do not overlap. n = length(x) if(length(w) == 1) w = rep(w,n) if(length(h) == 1) h = rep(h,n) y = numeric(n) step = median(h)/4 # loop in reverse order y[n] = 0 for(i in (n-1):1) { # find nearby boxes (those for which overlap is possible) is = (i+1):n j = is[abs(x[i]-x[is]) < (w[i]+w[is])/2] # now only need to avoid vertical overlap if(length(j) == 0) y[i] = 0 else y[i] = find.gap(h[i],y[j],h[j],symmetric) } y } ################# jitterplot <- function(x,scale=1,ylim=c(1-0.1/scale,1+0.1/scale),pch=1, xlab,ylab="jitter",...) { # xlab, scale, cex # jitter range is 1 +/- 0.1 if(missing(xlab)) xlab = deparse(substitute(x)) stripchart(x,method="jitter",ylim=ylim,pch=pch,xlab=xlab,ylab=ylab,...) } histodot <- function(x,symmetric=F,xlim=range(x), ylim=if(symmetric) c(-1,1) else c(0,1), cex=par("cex"),pch=par("pch"), main="",xlab,ylab="",...) { # uses: shuffle, yinch, move.collisions if(missing(xlab)) xlab = deparse(substitute(x)) x = shuffle(x) plot.new() plot.window(xlim,ylim) # this is for pch="H" symbol.inch = par("csi")*cex/1.5 if(pch == 1) symbol.inch = symbol.inch/2 y.inch = move.collisions.vertical(inchx(x),symbol.inch,symbol.inch,symmetric) y = yinch(y.inch) if(min(y) < ylim[1]) ylim = ylim/ylim[1]*min(y) if(max(y) > ylim[2]) ylim = ylim/ylim[2]*max(y) plot.window(xlim,ylim); box(); axis(1) points(x,y,cex=cex,pch=pch,...) title(main=main,xlab=xlab,ylab=ylab) } density.curve <- function(x,adjust=1,bw="bcv",kernel="g",cut=0, xlab,main="",add=F,...) { if(missing(xlab)) xlab = deparse(substitute(x)) d = density(x,bw=bw,adjust=adjust,kernel=kernel,cut=cut,na.rm=T) if(add) { lines(d,...) } else { plot(d,xlab=xlab,main=main,type="l",...) } invisible(d) }