Part 4: Computing/Programming140.776 Statistical ComputingIngo RuczinskiThanks to Thomas Lumley and Robert Gentleman of the R-core group (http://www.r-project.org/) for providingsome tex files that appear in part in the following slides.optionsoptions package:base R DocumentationOptions SettingsDescription:Allow the user to set and examine a variety of global "options"which affect the way in which R computes and displays its results.Usage:options(...)getOption(x).OptionsArguments:...: any options can be defined, using ’name = value’.However, only the ones below are used in "base R".Further, ’options(’name’) == options()[’name’]’, see theexample.x: a character string holding an option name.Customizing your R environmentCreating a .Renviron file in your root directory is a good idea.For example, when you use the postscript command to createa figure, the default for the paper format is ’A4’. You can change itto format ’letter’ by typingpostscript("test.ps",paper="letter")but you had to do this every time you want to create a postscriptfile. However, if you put the lineR_PAPERSIZE=letterinto your .Renviron file, R will use the ’letter’ format as default.And it works in all subdirectories!Customizing your R environmentCreating a .Rprofile file in your root directory is an equally goodidea. For example, the commandsload("/home/iruczins/code/R/functions/.RData")options(width=200,defaultPackages=c(getOption("defaultPackages"),"nlme","rpart","survival","tree"))ensure that my own R functions are available when I start an Rsession, that the width of the screen output is 200 characters wide,and that the default libraries plus the libraries nlme, rpart, andsurvival are loaded. Works in all subdirectories!Type ?Startup to see what exactly R is doing.The BATCH modeThe command source reads in R code, but it is generally not agood idea to use it to run large simulations. This should be doneusing the R BATCH mode:R BATCH inputfile outputfile &In R, you have to use the & to run the job in the background. InSplus, this was not necessary.If you run jobs on a cluster and share CPU time with other users, don’t forget tonice them!systemThe function system allows you execute Unix/Linux commandsinside an R session. For example,> system("ls")> system("pwd")> system("ping www.google.com")lists the files in your current directory, shows you the path, and’pings’ Google.The function system.time (see also proc.time) lets you deter-mines how much time the currently running R process has alreadyconsumed.DebuggingR has some built-in program debugging tools. Check out the helpfiles for:browserdebugtracetracebackrecoverDefining functionsA function definition looks likemedian <- function(x, na.rm = FALSE){... lots of code...## a return valuesort(x, partial = half)[half]}This function has two arguments, x and na.rm. The second ar-gument has a default value of FALSE, so it is optional. The firstargument is required.The last line of the function computes a return value, which is notassigned to anything.Functionsstd.dev = function(x) sqrt(var(x))t.test.p = function(x,mu=0) {n=length(x)tv=sqrt(n)*(mean(x)-mu)/std.dev(x)return(2*(1-pt(abs(tv),n-1)))}t.stat = function(x,mu=0) {n=length(x)tv=sqrt(n)*(mean(x)-mu)/std.dev(x)list(t=tv,p=2*(1-pt(abs(tv),n-1)))}> z=rnorm(300,1,2)> t.stat(z)$t[1] 8.019781$p[1] 2.420286e-14Function argumentsSuppose the function is called asmedian(height*width)and height and width are c(1,2,3) and 10 respectively.The function argument x is now the value of height*width,c(10,20,30). It still has this value even if there are new vari-ables height and width defined inside the function.The argument na.rm has its default value FALSE, since itsvalue wasn’t specified.ExampleWe could complete the median function asmedian <- function (x, na.rm = FALSE){if (na.rm) {remove missing values}x <- x[!is.na(x)]else if (any(is.na(x)))return(NA)n <- length(x)half <- (n + 1)/2if (n%%2 == 1) { ## odd nsort(x)[half]}else { ## even nsum(sort(x)[c(half, half+ 1)])/2}}It sorts x and then returns the middle observation or the averageof the middle two.ScopeThere may be variables with the same name (eg x, n) inside andoutside a function. Rules for working out which one to use arecalledscoping rules. R’s are:First look for a variable inside the function.Then look for a variable inside the function that the function wasdefined in (if any), and so on up.Finally look in the global environment, the variables visible atthe command line.This is different from S, where the second step doesn’t happen, butthe difference only matters in some specialized cases. The thirdstep is usually accidental for looking up variables (it’s important forlooking up functions).Scope exampleThe value of this function is itself a functionpower <- function(lambda){function(x) {xˆlambda}}square <- power(2)cube <- power(3)> square(1:2)[1] 1 4> cube(1:2)[1] 1 8Inside the square function, what is lambda? It isn’t a local vari-able, so R looks at the function where square was defined. Herelambda exists. Its value was 2.Unevaluated argumentsEarlier we said that inmedian(height*width)the function argument x just stored the value of height*width.This isn’t quite true. Until you look at x it stores the whole expres-sion height*width. Graphics commands use this to get plotlabels, since the substitute function lets you copy the expres-sion without looking at it.> label <- function(x) {list(value=x,actual=substitute(x))}> label(1+1)$value[1] 2$actual1 + 1Thislazy evaluation is also used in handling model formulas.Avoiding iterationThe canonical bad R program looks like this:## multiply two vectorsfor(i in 1:n){d[i]<-a[i]*b[i]}## compute the inner products<-0for(i in 1:n){s<-s+d[i]}The right way to do this is:s<-sum(a*b)Multiplication, like many operations and functions, is vectorized: itworks elementwise on whole objects.Why avoid iteration?There are two reasons to replace loops with vectorized calcula-tions:Speed: the for() loop is much slower since the expressionmust be evaluated by the interpreter every time.Clarity: it is much easier to see what sum(a*b) does.Vectorized functionsThese includeThe operators &, |, !, +, -, *, /, ˆ, %%.Mathematical functions such as log, sin, pnorm, choose,gamma, and many more.Random number generators such as rnorm, rpois, ...ifelse for vectorized conditionals.Vector recyclingRecall the recycling rules:If a and b are vectors of the same