Notes for Using R

© 2008 Krishna Myneni


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.

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")

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
	>

Moving Data Between Vectors and Matrices

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
	>

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 odesolve 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 odesolve 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)                      # initial conditions
times <- c(0, 0.1*(1:100))         # time steps at which to compute the solution

require("odesolve")

# 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 full R script for the above example is given in odesolve-example.R.