Chapter 2: Basics of R Programming

R as a Smart Calculator

1 + 4 #The text behind the symbol # is a comment in R

## [1] 5

2 + pi/4-0.8 #pi is circumference ratio which is approximately 3.1415926

## [1] 1.985398

x <- 1 
# <- is the assignment symbol: it assigns 1 to x.
y <- 2
z <- 4
t <- 2*x^y-z #This statement is equivalent to 2*1^2-4 = -2
t

## [1] -2

u = 2 
v = 3
# Symbols "=" and "<-" mean the same in most cases
u + v

## [1] 5

sin(u*v) # u*v = 6 is considered a radian value

## [1] -0.2794155

Write a function in R

square <- function(x) x*x 
# This statement creates function where it is a block of code
# that runs when it is called.
# It may pass in data which are called parameters 
square(4)

## [1] 16

fctn <- function(x,y,z) {x+y-z/2} 
fctn(1,2,3)

## [1] 1.5

Plot with R

plot(sin, -pi, 2*pi) #Plot the curve of y = sin(x) from -pi to 2*pi

square <- function(x) x*x #Define a function
plot(square, -3, 2)# Plot the defined function

# Plot a 3D surface using the following code
x <- y <- seq(-1, 1, length=100)
Z <- outer(x, y, function(x, y){1-x^2-y^2})
# outer(x,y, function) gives the z values defined on x, y grid
persp(x = x, y = y, z = Z, theta = 310)

# yields a 3D surface with perspective angle 310 deg

Symbolic Calculation for Calculus by R

D(expression(x^2,'x'), 'x') # Take derivative of x^2 which is 2x

## 2 * x

fx <- expression(x^2,'x') 
#define a symbolic function with variable x 
D(fx,'x') #differentiate the function with respect to x

## 2 * x

fx <- expression(x^2*sin(x),'x')
#Change the expression and use the same derivative command
D(fx,'x')

## 2 * x * sin(x) + x^2 * cos(x)

fxyz <- expression(x^2+y^2+z^2, 'x','y','z')
#define a function of two or more variables
fxyz #This gives the expression of the function in x, y and z

## expression(x^2 + y^2 + z^2, "x", "y", "z")

D(fxyz,'x') #gives the partial derivative with respect to x: 2 * x

## 2 * x

D(fxyz,'y') #gives the partial derivative with respect to y: 2 * y

## 2 * y

square <- function(x) x^2
integrate(square, 0, 1)

## 0.3333333 with absolute error < 3.7e-15

#Numerically integrating x^2 from 0 to 1 
# gives 0.3333333 with absolute error < 3.7e-15

integrate(cos, 0, pi/2)

## 1 with absolute error < 1.1e-14

#Integrates cos(x) from 0 to pi/2 
# gives 1 with absolute error < 1.1e-14

Vectors and Matrices

c(1, 6, 3, pi, -3) #Enter data inside c() to create a 5x1 column vector

## [1]  1.000000  6.000000  3.000000  3.141593 -3.000000

seq(2, 6) #Generate a sequence from 2 to 6

## [1] 2 3 4 5 6

seq(1, 10, 2) # Generate a sequence from 1 to 10 with an increment of 2

## [1] 1 3 5 7 9

x <- c(1, -1, 1, -1)
x + 1 #1 is added to each element of x

## [1] 2 0 2 0

2*x #2 multiplies each element of x

## [1]  2 -2  2 -2

x/2 #Each element of x is divided by 2

## [1]  0.5 -0.5  0.5 -0.5

y <- seq(1, 4)
x*y # * multiples each pair of elements

## [1]  1 -2  3 -4

t(x) # Transpose of x matrix

##      [,1] [,2] [,3] [,4]
## [1,]    1   -1    1   -1

t(x)%*%y #Matrix multiplication: 1x4 matrix times a 4x1 matrix

##      [,1]
## [1,]   -2

#equivalent to a dot product

x%*%t(y) #4x1 matrix times a 1x4 matrix yields a 4x4 matrix

##      [,1] [,2] [,3] [,4]
## [1,]    1    2    3    4
## [2,]   -1   -2   -3   -4
## [3,]    1    2    3    4
## [4,]   -1   -2   -3   -4

mx <- matrix(x, nrow=2)
#Convert a vector y into a matrix of 2 rows. 
mx

##      [,1] [,2]
## [1,]    1    1
## [2,]   -1   -1

#or simply
my <- matrix(y,2)
my #The matrix elements go by column, first column, second, etc

##      [,1] [,2]
## [1,]    1    3
## [2,]    2    4

dim(my) #Find dimensions of a matrix

## [1] 2 2

as.vector(my) #Convert a matrix to a vector, also via columns

## [1] 1 2 3 4

mx*my #Multiplication between each pair of elements

##      [,1] [,2]
## [1,]    1    3
## [2,]   -2   -4

mx/my #Division between each pair of elements

##      [,1]       [,2]
## [1,]  1.0  0.3333333
## [2,] -0.5 -0.2500000

mx-2*my

##      [,1] [,2]
## [1,]   -1   -5
## [2,]   -5   -9

mx%*%my

##      [,1] [,2]
## [1,]    3    7
## [2,]   -3   -7

det(my) #Determinant of a square matrix

## [1] -2

myinv <- solve(my) #Find inverse of a matrix 
myinv

##      [,1] [,2]
## [1,]   -2  1.5
## [2,]    1 -0.5

myinv%*%my #verifies the inverse of a matrix and results in identity matrix

##      [,1] [,2]
## [1,]    1    0
## [2,]    0    1

diag(my) #outputs the diagonal vector of a matrix

## [1] 1 4

myeig <- eigen(my)
myeig #yields the two eigenvalues and eigenvectors

## eigen() decomposition
## $values
## [1]  5.3722813 -0.3722813
## 
## $vectors
##            [,1]       [,2]
## [1,] -0.5657675 -0.9093767
## [2,] -0.8245648  0.4159736

myeig$values #outputs only eigenvalues

## [1]  5.3722813 -0.3722813

myeig$vectors #outputs only eigenvectors

##            [,1]       [,2]
## [1,] -0.5657675 -0.9093767
## [2,] -0.8245648  0.4159736

mysvd <- svd(my) #SVD decomposition of a matrix M = UDV'
#SVD can be done for any m-by-n rectangular matrix 
mysvd

## $d
## [1] 5.4649857 0.3659662
## 
## $u
##            [,1]       [,2]
## [1,] -0.5760484 -0.8174156
## [2,] -0.8174156  0.5760484
## 
## $v
##            [,1]       [,2]
## [1,] -0.4045536  0.9145143
## [2,] -0.9145143 -0.4045536

#outputs d, U, and V where d is the SVD eigenvalues
# U is the spatial eigenvectors, and V is the temporal eigenvectors

mysvd$d #outputs d only, as a vector

## [1] 5.4649857 0.3659662

U <- mysvd$u 
U #outputs the U matrix

##            [,1]       [,2]
## [1,] -0.5760484 -0.8174156
## [2,] -0.8174156  0.5760484

D <- diag(mysvd$d) 
#generates the diagonal D matrix

V <- mysvd$v 
V #outputs the V matrix

##            [,1]       [,2]
## [1,] -0.4045536  0.9145143
## [2,] -0.9145143 -0.4045536

U%*%D%*%t(V) #recovers the original matrix 'my'

##      [,1] [,2]
## [1,]    1    3
## [2,]    2    4

b <- c(1, 3)
ysol <- solve(my, b) 
#solve linear equations my %*% x = b 
ysol

## [1]  2.5 -0.5

my%*%ysol #verifies the solution, equal to b

##      [,1]
## [1,]    1
## [2,]    3

Statistics

x <- rnorm(8) 
#generates 8 random numbers in (-infty, infty)
#with normal distribution with default mean = 0, and default sd = 1
#If you want to set a specific mean and sd, you need to place in the
#function call
x

## [1]  1.6561469  0.4224978 -1.7913717  0.2452649 -0.6139012 -1.9181532 -0.6218588
## [8] -0.3628891

mean(x) #mean

## [1] -0.373033

var(x) #variance

## [1] 1.378835

sd(x) #standard deviation

## [1] 1.174238

median(x) #median

## [1] -0.4883951

quantile(x) #0%         25%         50%         75%        100%

##         0%        25%        50%        75%       100% 
## -1.9181532 -0.9142370 -0.4883951  0.2895731  1.6561469

range(x) #yields the min and max of x

## [1] -1.918153  1.656147

max(x) #maximum value of x

## [1] 1.656147

boxplot(x) #yields the box plot of x

w <- rnorm(1000) 
#generate 1000 random numbers in N(0,1) distribution

z <- rnorm(10000, mean = 10, sd = 5) 
#generate 100 random numbers following N(10,5^2)
#mean = 10, standard deviation = 5

summary(rnorm(12)) #statistical summary of the data sequence

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -1.1107 -0.5656 -0.1532  0.0207  0.4940  1.5796

hist(w) #plot the histogram of 1000 random numbers

#Find the current directory
getwd()

## [1] "/Users/momtaza/Desktop/RMathModel"

##Change directory 
#setwd("/Users/sshen/Desktop/MyDocs/teach") #for Mac 
#setwd("D:/My Documents/Desktop/teach") #for PC

R Tutorials

setwd("~/Desktop/RMathModel/data")

##Commands to input data files:

mydata <- read.csv("mydata.csv") 
# read csv file named "mydata.csv"

mydata <- read.table("mydata.txt") 
# read text file named "mydata.txt"

library(gdata)

## 
## Attaching package: 'gdata'

## The following object is masked from 'package:stats':
## 
##     nobs

## The following object is masked from 'package:utils':
## 
##     object.size

## The following object is masked from 'package:base':
## 
##     startsWith

# load the gdata package 

# Install packages by using install.packages
# install.packages("readxl")
library("readxl")

mydata <- read_excel("mydata.xlsx") 
# read an excel file

library(foreign) 
# load the foreign package 

#mydata <- read.mtp("mydata.mtp") 
# read from .mtp file


#mydata <- read.spss("myfile", to.data.frame=TRUE)
#read a spss file

ff <- tempfile()
cat(file = ff, "123456", "987654", sep = "\n") 
read.fortran(ff, c("F2.1","F2.0","I2"))

##    V1 V2 V3
## 1 1.2 34 56
## 2 9.8 76 54

# read a fotran file

library(ncdf4)
# load the ncdf4 package 

ncin <- ncdf4::nc_open("air.mon.mean.nc") 
# open a NetCDF file

lon <- ncvar_get(ncin, "lon") 
#read data "lon" from a netCDF file into R

library("rjson")
# load the rjson package 

jd <- fromJSON(file = "argo4903285_59.json")
# read data from a JSON file argo4903285_59.json

library(ncdf4) 
#Suppose the following error message pops up when you load ncdf4 package
#Error in library(ncdf4) : there is no package called ncdf4 

#Then you can install the R package by
#install.packages("ncdf4")

#R packages: animation, chron, e1071, fields, ggplot2, lattice,
#latticeExtra, maps, mapdata, mapproj, matrixStats, ncdf4,
#NLRoot, RColorBrewer, rgdal, rasterVis, raster, rjson, sp, TTR

#The zipped data for this book can be downloaded from: 
#climatemathematics.sdsu.edu/data.zip

#To load a single package, such as "animation", you can do
library(animation)

#You can also load all these packages at once using pacman package
#install.packages("pacman")
library(pacman)
pacman::p_load(animation, chron, e1071, fields, ggplot2, lattice,
               latticeExtra, maps, mapdata, mapproj, matrixStats, ncdf4,
               NLRoot, RColorBrewer, rasterVis, raster, rjson, sp, TTR)

Introduction to Modern Mathematical Modeling with R:

A User’s Manual to Train Mathematical Consultants

A Cambridge University Press book by

SSP Shen

Chapter 2: Basics of R Programming

R as a Smart Calculator

Write a function in R

Plot with R

Symbolic Calculation for Calculus by R

Vectors and Matrices

Statistics

R Tutorials