Notes for Using R

© 2008--2013 Krishna Myneni



Index

  1. Basic Info
    1. Startup
    2. Exiting R
    3. Loading an R Program
    4. Help for Built-In Functions
  2. Functions
    1. Defining Functions
    2. Calling Functions
    3. Plotting Functions
    4. Roots of a Function
  3. Complex Numbers
    1. Setting a Complex Variable
    2. Parts of a Complex Variable or Number
    3. Functions for Complex Numbers
  4. Vectors, Matrices, and Higher Dimensional Arrays
    1. Creating Vectors and Matrices
    2. Properties of Vectors and Matrices
    3. Manipulating Vectors and Matrices
    4. Reading and Writing Vectors and Matrices From/To a File
  5. Plotting Data
    1. 2-D Plotting From Vectors
    2. 2-D Plotting From a Matrix
    3. Setting Plot Attributes
    4. Overlaying Multiple Plots of Data
  6. Matrix Computations
    1. Matrix Multiplication
    2. Matrix Inversion
    3. Trace of a Matrix
    4. Eigenvalues and Eigenvectors of a Matrix
  7. R Packages
    1. Description of Additional Packages
    2. Installing a Package
    3. Loading an Installed Package
    4. Example: Numerical Solution of Ordinary Differential Equations

I. Basic Info

R is a free and powerful computing environment, designed for numerical computations and visualizing data. It is available for a wide variety of systems, and may be obtained from http://www.r-project.org. R is also a language for expressing, and programming numerical calculations. The language of R is object-oriented and elegant, with a simple and intuitive syntax. While many of the built-in functions of R are useful for performing calculations in statistics, the large variety of functions, as well as package extensions, allow R to be applied to a wide range of problems in the physical sciences and applied mathematics. These notes illustrate some common computational tasks encountered by students, teachers, and researchers in the physical sciences.

Startup

On a Unix/Linux workstation, open a terminal window and type "R" at the prompt. A startup message such as the following will be displayed:

	R version 2.6.2 (2008-02-08)
	Copyright (C) 2008 The R Foundation for Statistical Computing
	ISBN 3-900051-07-0

	R is free software and comes with ABSOLUTELY NO WARRANTY.
	You are welcome to redistribute it under certain conditions.
	Type 'license()' or 'licence()' for distribution details.

	  Natural language support but running in an English locale

	R is a collaborative project with many contributors.
	Type 'contributors()' for more information and
	'citation()' on how to cite R or R packages in publications.

	Type 'demo()' for some demos, 'help()' for on-line help, or
	'help.start()' for an HTML browser interface to help.
	Type 'q()' to quit R.

	>


The > is R's command prompt.

Exiting R

As shown in the startup message, you may enter the command, q(), to quit and return to the system prompt.

Loading an R program (script)


	> source("filename")

where filename is the name of a file, typically with extension .R, containing the R commands to be executed.

Help for Built-In Functions

R provides a built in help system, similar to man pages on a Unix system. To obtain help on a built-in function, such as the seq function, which returns a sequence of values as a vector, type

        > help("seq")

If you are unsure of the exact function name, you may also search the documentation for a particular topic, using help.search. For example,
 
        > help.search("sequence")

will show relevant functions for sequence generation.

II. Functions

Defining Functions

Simple Function of One Variable


	> myfunc <- function(x) {
      		y <- cos(x) + sin(x)
      		return (y)
  	}
	>

Calling Functions

Call with real argument

 
	> myfunc(0.5)
	[1] 1.357008
	>

Call with complex argument (pure imaginary)

 
	> myfunc(0.5*1i)
	[1] 1.127626+0.521095i
	>

Call with complex argument (with real and imaginary parts)

 
	> myfunc(-2 + 5*1i)
	[1] -98.36115+36.59336i
	> 

Call with a 1-D Array (vector) of Values


	> x <- seq(-4, 4, by=0.01)
	> y <- myfunc(x)

x is a vector, generated by the seq function. It may be passed as an argument to the function, and a vector, y, is returned. Assignment of the variables x and y may also be done with the more conventional notation,

	> x = seq(-4, 4, by=0.01)
	> y = myfunc(x)

Plotting Functions


	> plot(myfunc, 0, 2*pi)

Overlaying another function plot


	> plot(myfunc2, 0, 2*pi, add=TRUE)

Setting the plot color


	> plot(myfunc2, 0, 2*pi, add=TRUE, col="red")

Roots of a Function

We will take as an example, the problem of finding the roots of the function,


          f(x) = f_g(x) - f_l(x)


where

          f_g(x) = exp(-a*x^2)

and

          f_l(x) = b/(x^2 + b)

For the example, we will take a=1, b=0.5.

The problem of finding the roots of f(x) is the same as the problem of finding the values of x at which the two functions, f_g(x) and f_l(x), intersect. The problem can be solved approximately by plotting the two functions, and visually determining the values of x at which they intersect --- the graphical method. A more precise method is to use an iterative search of $x$, to find the values at which f(x) = 0. R provides such a search method with the function uniroot. However, we must first obtain the approximate roots with the graphical method, so that the search interval can be bounded for uniroot.

Graphical method

Set the values of a and b:

      > a = 1
      > b = 0.5
Define the functions, f_g(x), and f_l(x), and plot them together on a graph:

      > fg <- function(x) { return( exp(-a*x^2) ) }
      > fl <- function(x) { return( b/(x^2 + b) ) }
      > plot(fg, -5, 5, col='blue')
      > plot(fl, -5, 5, col='red', add=TRUE)

The following graph will appear.



From the graph we see that there appear to be three points of intersection, one near x = -1, one near x = 0, and one near x = +1. We could read these values of x somewhat more precisely from the graph; however, now that we know the vicinty of the intersection points, we can find them precisely through an interative search, or "root finding".

Numerical method

First, define the function for which we want to find the root. For the example above,

	> f <- function(x) { return( fg(x) - fl(x) ) }

Now, find the precise value of a given root (value of x at which fg(x) and fl(x) intersect), using its approximate vicinity. From the graphical method, we know there is a root near x = -1. Using uniroot, a precise value may be found by specifying the interval in x over which to search for the root,

	> uniroot(f, c(-1.2, -0.8))
	$root
	[1] -1.120907

	$f.root
	[1] -4.398283e-08

	$iter
	[1] 3

	$estim.prec
	[1] 0.0001094043

uniroot searches for the root in the interval, x = (-1.2, -0.8), and outputs several values: the root itself, the value of the function f at the root position, the number of iterations required to find the root, and the estimated precision of the root. From the output above, we see that this root occurs at

x = -1.12091 +/- 0.0001


Additional arguments can be given to the uniroot function, if a more precise value of the root is desired. Use the built-in help function to find out how to set the tolerance for the search. The object returned by uniroot is a list, and the different values in the output above are the different elements of the list. If only the value of the root is desired, for example, to store in a variable for later use, one can retrieve only the root value from the list as follows,

	> root1 = uniroot(f, c(-1.2, -0.8))$root
	> root1
	[1] -1.120907



III. Complex Numbers

Setting a complex variable


	> z <- 3.2+4*1i

or

	> z = 3.2+4*1i

Parts of a complex variable or number

Example: Illustrate the functions Re, Im, Mod, abs, and Arg.

	> Re(z)
	[1] 3.2
	> Im(z)
	[1] 4
	> Mod(z)
	[1] 5.122499
	> abs(z)
	[1] 5.122499
	> Arg(z)
	[1] 0.8960554

The Mod and abs functions both return the length (modulus) of the number, and Arg returns the argument (phase).

Functions for Complex Numbers

All of the standard transcendental functions may be applied to complex numbers as well as to real numbers, such as sin, cos, exp, log, etc. A function specific to complex numbers is Conj, which returns the complex conjugate of the number:

	> x = 3 + 1i*2
	> Conj(x)
	[1] 3-2i
	>


IV. Vectors, Matrices, and Higher Dimensional Arrays

Creating Vectors and Matrices

Creating a vector

Use the combine function, c( ... ), to combine a series of values into a vector.

	> x <- c(1, 2, 3, 11, 12, 13)
	> x
	[1]  1  2  3 11 12 13
	> 

Creating a real matrix

Example: Use the matrix function to create a real 2x3 matrix.

	> mat1 <- matrix(c(1,2,3,11,12,13), nrow=2, byrow=TRUE)
	> mat1
	     [,1] [,2] [,3]
	[1,]    1    2    3
	[2,]   11   12   13
	>

Creating a complex matrix

Example: Use the matrix funtion to create a complex 2x2 matrix.

	> mat2 <- matrix(c(1+1i, 0, 0, 1-1i), nrow=2, byrow=TRUE)
	> mat2
	     [,1] [,2]
	[1,] 1+1i 0+0i
	[2,] 0+0i 1-1i
	> 

Properties of Vectors and Matrices

Length of a vector

Use the length function to determine the number of elements in a vector.
Example:
       > x <- seq(-4, 4, by=0.01)
       > length(x)
       [1] 801
shows that we have created an 801 element vector using the seq function.

Dimensions of a matrix

Use the dim function to determine the dimensions of a matrix, or higher dimensional array.
Example:
       > mat1 <- matrix(c(1,2,3,11,12,13), nrow=2, byrow=TRUE)
       > dim(mat1)
       [1] 2 3
dim returns a vector containing the dimensions of the matrix.

Manipulating Vectors and Matrices

Concatenating two vectors into a single vector

Example:
       > x <- 1:10
       > x
       [1]  1  2  3  4  5  6  7  8  9 10
       > y <- 11:20
       > y
       [1] 11 12 13 14 15 16 17 18 19 20
       > z <- c(x, y)
       > z
       [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20
In the above example, the vectors x and y are 10 element vectors. The combine function, c(), is used to combine (concatenate) the elements together to form the vector, z.

Making a vector from a subset of an existing vector

Example:
       > x <- seq(-2, 2, by=0.1)
       > length(x)
       [1] 41
       > y <- x[18:24]
       > y
       [1] -0.3 -0.2 -0.1  0.0  0.1  0.2  0.3
In this example, we create a 41 element vector, x, and use the range 18:24, to specify the index range for elements in x, which are returned to form the new vector y.

Building a Matrix From Two Vectors

Example: Use the array function to build a two column matrix from two vectors.

	> x <- seq(-4, 4, by=0.01)
	> y <- myfunc(x)
	> d <- array(c(x,y), c(length(x), 2))

x and y are vectors (1-d arrays), generated from a sequence, and a function, respectively.

d is a 2-column matrix, with column 1 containing x, and column 2 containing y.

Extracting Column of a Matrix Into a Vector


	> x <- d[,1]    
	> y <- d[,2] 

The first column of matrix d is placed in vector x, and the second column of d is placed in vector y.

Reading and Writing Vectors and Matrices from/to a File

Reading a vector from an ascii file


	> d <- scan("filename")

Reading a two column matrix from an ascii file


	> d <- matrix(scan("filename"), ncol=2, byrow=TRUE)

Writing a vector into a single column of a file


	> write(x, ncolumns=1, file="array.dat")

Writing a matrix into a file


	> write(t(mat1), ncolumns=3, file="matrix.dat")

Notice the use of the transpose function, t(). The matrix must be transposed, using t(), to make the columns in the file correspond to the columns of the matrix.

V. Plotting Data

2-D Plotting From Vectors

Example: Use the plot function to plot a vector of y-values versus a vector of x-values.

	> x <- seq(-4, 4, by=0.01)
	> y <- myfunc(x)
	> plot(x,y)

x and y are vectors of the same length.

2-D Plotting From a Matrix

Example: Use the plot function to plot the data stored in a two-column matrix, with the first column containing the x-values, and the second column containing the y-values.

	> d <- array(c(x,y), c(length(x), 2))
 	> plot(d)

Setting Plot Attributes

Example: Use the plot function to plot data from the matrix d, using a line plotting symbol, with color red, and axes labels for the x-axis and y-axis.

	> plot(d, type="l", col="red", xlab="Time (s)", ylab="Vout (V)")

Other plot attributes such as title, xlim, and ylim may also be specified.

Overlaying Multiple Plots of Data

Example: Use the plot and lines functions to overlay three plots on a single graph.

	> plot(d, xlim=c(-4, 4), ylim=c(-1,1))
	> lines(x, y, col="blue")
	> lines(d2, col="green")

The first plot command sets up the graph and plots the matrix d. Subsequent lines commands overlays plots of the respective data onto the graph.

VI. Matrix Computations

Matrix Multiplication

Example: Use the outer product operator, %*%, to multiply together two matrices.

	> a <- matrix(c(1, 0, -2, 0, 3, -1), nrow=2, byrow=TRUE)
	> b <- matrix(c(0, 3, -2, -1, 0, 4), nrow=3, byrow=TRUE)
	> c <- a%*%b
	> c
     		[,1] [,2]
	[1,]    0   -5
	[2,]   -6   -7
	> d <- b%*%a
	> d
     		[,1] [,2] [,3]
	[1,]    0    9   -3
	[2,]   -2   -3    5
	[3,]    0   12   -4
	> 

Matrix Inversion

Example: Use the function solve to invert a 3x3 matrix.
	> a <- matrix(c(1, 0, -2, 4, 1, 0, 1, 1, 7), nrow=3, byrow=TRUE)
	> a
     		[,1] [,2] [,3]
	[1,]    1    0   -2
	[2,]    4    1    0
	[3,]    1    1    7
	> b <- solve(a)
	> b
     		[,1] [,2] [,3]
	[1,]    7   -2    2
	[2,]  -28    9   -8
	[3,]    3   -1    1
	> 

Trace of a Matrix

Compute the trace of a matrix by retrieving the diagonal elements of the matrix into a vector, and computing the sum of the elements in the vector. This can be done in a single command, as follows.
	> a <- matrix(c(1, 0, -2, 4, 1, 0, 1, 1, 7), nrow=3, byrow=TRUE)
	> sum(diag(a))
        [1] 9
The diag function retrieves the diagonal elements of the matrix.

Eigenvalues and Eigenvectors of a Matrix

Example: Use the function eigen to compute the eigenvalues and eigenvectors a real 2x2 matrix, and a complex 3x3 matrix.

Real Matrix


	> sigma_x <- matrix(c(0,1,1,0), nrow=2, byrow=TRUE)
	> sigma_x
	     [,1] [,2]
	[1,]    0    1
	[2,]    1    0
	> eigen(sigma_x)
	$values
	[1]  1 -1

	$vectors
	          [,1]       [,2]
	[1,] 0.7071068 -0.7071068
	[2,] 0.7071068  0.7071068

	> 

Complex Matrix


	> Jy <- matrix(c(0,1,0,-1,0,1,0,-1,0), nr=3, byrow=TRUE)
	> Jy
	     [,1] [,2] [,3]
	[1,]    0    1    0
	[2,]   -1    0    1
	[3,]    0   -1    0
	> Jy <- Jy*(0+1i)/sqrt(2)
	> Jy
	             [,1]         [,2]         [,3]
	[1,] 0+0.0000000i 0+0.7071068i 0+0.0000000i
	[2,] 0-0.7071068i 0+0.0000000i 0+0.7071068i
	[3,] 0+0.0000000i 0-0.7071068i 0+0.0000000i
	> eigen(Jy)
	$values
	[1]  1.000000e+00 -2.376571e-16 -1.000000e+00

	$vectors
	                [,1]                       [,2]            [,3]
	[1,] -0.5+0.0000000i 7.071068e-01+0.000000e+00i -0.5+0.0000000i
	[2,]  0.0+0.7071068i 0.000000e+00+2.498002e-16i  0.0-0.7071068i
	[3,]  0.5+0.0000000i 7.071068e-01+0.000000e+00i  0.5+0.0000000i

	> 


VII. R Packages

Description of Additional Packages

Visit http://cran.r-project.org and click on Packages. An alphabetical list of available packages for R will be displayed. Click on a particular package name to obtain additional information about the package, including a pdf document which describes its use.

It is not necessary to download the package source for installation if the machine on which you installed R has internet access.

Installing a Package


	> install.packages()

A list of mirror sites for R packages will appear. Select an appropriate mirror site. Then a list of available packages will appear. Select the desired package. The package will be downloaded, built, and installed.

Loading an Installed Package

From within the R environment, or within an R script, you may load a package using,


	> require("packagename")

Example: Numerical Solution of Ordinary Differential Equations

We will use the extra R package deSolve to illustrate how to solve a system of ordinary differential equations. This package provides the robust ODE solver named lsoda taken from the well-known collection of solvers, ODEPACK, written by Alan Hindmarsh and others. Assuming that you have installed the deSolve package using the procedure given above, you may use it to find a solution to the following system of first order differential equations:


	du/dt = -u + alpha*v
	
	dv/dt = beta*u - v

where alpha and beta are two specified parameters. Although the above set of coupled, linear, first-order differential equations may be solved easily by linear methods, they are used to illustrate the general numerical procedure, which may be applied to a more complex system of linear or nonlinear differential equations.

First, define a function which computes the derivatives, given the current values of the variables, the independent variable (time, for example), and the values of the parameters:


# Function to evaluate the ODEs

derivs <- function(t, x, p) {
	udot = -x[1] + p[1]*x[2]
	vdot = p[2]*x[1] - x[2]	
	return (list(c(udot, vdot)))
}

The function named derivs takes three arguments: the value of the independent variable (time), a vector x containing the state of the system at the current time, and a vector containing the parameter values. Thus, at the current value of t,

	u = x[1]
	v = x[2]

	alpha = p[1]
	beta = p[2]

and derivs returns a list of the values of the derivatives. In order to compute a numerical solution of the equations using lsoda, we must supply it with a vector of time values for which we want the solution, and the initial conditions, i.e. the values of the variables at the starting time. The vector of parameter values must also be specified. Thus, we must do the following:

params <- c( 2, 0.1 )              # alpha = 2, beta = 0.1
x0 <- c(0, 0.05)                   # initial conditions
times <- c(0, 0.1*(1:100))         # time steps at which to compute the solution

require("deSolve")

# Solve the differential equations for this example
out <- lsoda( x0, times, derivs, params)

The solutions, u(t) and v(t), are stored in the array out. They may be plotted using the following commands.

> plot(out[,1], out[,2], type='l', col="blue", xlab="Time (s)", ylab="u(t), v(t)", ylim=c(0, .05))
> lines(out[,1], out[,3], col="red")

The following graph will be shown.



The full R script for the above example is given in odesolve-example.R.