Math330 R Lab

This worksheet will re-cap some of the key R techniques you will need for this course.

Function writing

The general structure for writing a function in R is

myfunc<-function(<ARGUMENTS>){
 <CONTENTS that take the arguments and produce a result>
 return(result)
}

A simple example is a function that returns the square of a number:

myfunc<-function(x){
 result<-x^2
 return(result)
}

Alternatively and more concisely:

square<-function(x){
 x^2
}

since R returns, by default, the result of the last evaluation in the function.

Functions can accept vectors: e.g. try square(1:10). Functions can also take multiple arguments, for example the following function adds two numbers a and b:

add<-function(a,b){
 return(a+b)
}

Exercise: Try writing a function that takes two numbers as arguments and returns the square of the difference between them.

Function plotting

In order to plot a function in R, we need to create a vector of values at which to evaluate the function. We use the seq function to do this, which constructs a sequence starting from the specified from value, to the specified to value, of length length. See what is returned if you type seq(from=3,to=7,length=10).

The below code will produce a plot $y=x^2$ between $-5$ and $5$.

x<-seq(from=-5,to=5,length=1000)
plot(x,square(x),type=’l’)

Here the type=’l’ part tells R to join the points: see what happens if you omit this.

Once you have a plot it’s easy to add straight lines to it, using the abline function:

abline(h=2,col=’red’)
abline(v=1,col=’green’)
abline(a=10,b=-2,col=’blue’)

Be careful with multiple arguments that are vectors:

a<-seq(from=1,to=10,length=10)
b<-seq(from=1,to=5,length=5)
plot(a,add(a,b),type=’l’)

Try to ensure that your likelihood and deviance functions only have the unknown parameter $\theta$ as an argument (otherwise things get tricky!)

Root finding

It is often the case in statistics that we cannot find a root to an equation. One common example is finding the solutions to $D(\theta )=3.84$, i.e. finding the deviance-based confidence interval. Note that ‘roots’ are values where the function has value 0, so we need to re-arrange this to finding the value $\theta$ such that $f(\theta )=D(\theta )-3.84=0$.

Suppose we are interested in finding the roots of $f(x)=e^x-x-2$:

#specify the function
f<-function(x){
 return(exp(x) - x - 2)
}

#plot to find approx. root locations
x<-seq(from=-3,to=3,length=1000)
plot(x,f(x),type=’l’)
abline(h=0,lty=2)

We see there are two roots to the function: one is between $-1$ and $-2$, the other is between 1 and 2. R has a simple root finding function called uniroot, which takes a function f and looks for roots between lower and upper:

#first root
uniroot(f=f,lower=-2,upper=-1)

#second root
uniroot(f=f,lower=1,upper=2)

#error message
uniroot(f=f,lower=-2,upper=2)

The last line shows that uniroot can’t cope with two roots between lower and upper.