Appendix A: A Tutorial of R and R Studio

Create your working directory named LinAlg and go there

Run the following, replacing “alyss/OneDrive/Documents” with your directory, to set the working directory: setwd(“/Users/alyss/OneDrive/Documents/LinAlg”).

Run getwd() to verify that you are in the directory/folder.

To save a figure

Before plotting the figure, run “setEPS(),” then run “postscript(”figname.eps”, width = w, height = h),” replacing figname, w, and h with the desired parameters. The figure will not show up as an output, and will instead be a file on your computer.

Basic Commands and Calculations

1+4
## [1] 5
#[1] 5  This is the result
2+pi/4-0.8
## [1] 1.985398
x<-1
y<-2
z<-4
t<-2*x^y-z
t
## [1] -2
u=2        # "=" sign and "<-" are almost equivalent
v=3        # The text behind the "#" sign is comments
u+v
## [1] 5
sin(u*v)    # u*v = 6 in the sine function is considered radian by R 
## [1] -0.2794155
#Enter temperature data in c() 
tmax <- c(77, 72, 75, 73, 66, 64, 59)
#Show the data
tmax
## [1] 77 72 75 73 66 64 59
#Generate a sequence using different methods
seq(1,8)
## [1] 1 2 3 4 5 6 7 8
seq(8)
## [1] 1 2 3 4 5 6 7 8
seq(1,8, by=1)
## [1] 1 2 3 4 5 6 7 8
seq(1,8, length=8)
## [1] 1 2 3 4 5 6 7 8
seq(1,8, length.out =8)
## [1] 1 2 3 4 5 6 7 8
#Define a function
samfctn <- function(x) x*x
samfctn(4)
## [1] 16
fctn2 <- function(x,y,z) x+y-z/2
fctn2(1,2,3)
## [1] 1.5

Basic Plots

#Plot temperature data
plot(1:7, c(77, 72, 75, 73, 66, 64, 59))

#More plot examples
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
x <- seq(-1, 1, length=100)
y <- seq(-1, 1, length=100)
z <- outer(x, y, function(x, y)(1-x^2-y^2)) #outer (x,y, function) renders z function on the x, y grid
persp(x,y,z, theta=330)  # yields a 3D surface with perspective angle 330 deg

#Contour plot
contour(x,y,z) #lined contours

filled.contour(x,y,z) #color map of contours

Symbolic calculations

D(expression(x^2,'x'), 'x') 
## 2 * x
# Take derivative of x^2 w.r.t. x 
2 * x #The answer is 2x
##   [1] -2.00000000 -1.95959596 -1.91919192 -1.87878788 -1.83838384 -1.79797980
##   [7] -1.75757576 -1.71717172 -1.67676768 -1.63636364 -1.59595960 -1.55555556
##  [13] -1.51515152 -1.47474747 -1.43434343 -1.39393939 -1.35353535 -1.31313131
##  [19] -1.27272727 -1.23232323 -1.19191919 -1.15151515 -1.11111111 -1.07070707
##  [25] -1.03030303 -0.98989899 -0.94949495 -0.90909091 -0.86868687 -0.82828283
##  [31] -0.78787879 -0.74747475 -0.70707071 -0.66666667 -0.62626263 -0.58585859
##  [37] -0.54545455 -0.50505051 -0.46464646 -0.42424242 -0.38383838 -0.34343434
##  [43] -0.30303030 -0.26262626 -0.22222222 -0.18181818 -0.14141414 -0.10101010
##  [49] -0.06060606 -0.02020202  0.02020202  0.06060606  0.10101010  0.14141414
##  [55]  0.18181818  0.22222222  0.26262626  0.30303030  0.34343434  0.38383838
##  [61]  0.42424242  0.46464646  0.50505051  0.54545455  0.58585859  0.62626263
##  [67]  0.66666667  0.70707071  0.74747475  0.78787879  0.82828283  0.86868687
##  [73]  0.90909091  0.94949495  0.98989899  1.03030303  1.07070707  1.11111111
##  [79]  1.15151515  1.19191919  1.23232323  1.27272727  1.31313131  1.35353535
##  [85]  1.39393939  1.43434343  1.47474747  1.51515152  1.55555556  1.59595960
##  [91]  1.63636364  1.67676768  1.71717172  1.75757576  1.79797980  1.83838384
##  [97]  1.87878788  1.91919192  1.95959596  2.00000000
fx= expression(x^2,'x')  #assign a function  
D(fx,'x') #differentiate the function w.r.t. x
## 2 * x
2 * x  #The answer is 2x
##   [1] -2.00000000 -1.95959596 -1.91919192 -1.87878788 -1.83838384 -1.79797980
##   [7] -1.75757576 -1.71717172 -1.67676768 -1.63636364 -1.59595960 -1.55555556
##  [13] -1.51515152 -1.47474747 -1.43434343 -1.39393939 -1.35353535 -1.31313131
##  [19] -1.27272727 -1.23232323 -1.19191919 -1.15151515 -1.11111111 -1.07070707
##  [25] -1.03030303 -0.98989899 -0.94949495 -0.90909091 -0.86868687 -0.82828283
##  [31] -0.78787879 -0.74747475 -0.70707071 -0.66666667 -0.62626263 -0.58585859
##  [37] -0.54545455 -0.50505051 -0.46464646 -0.42424242 -0.38383838 -0.34343434
##  [43] -0.30303030 -0.26262626 -0.22222222 -0.18181818 -0.14141414 -0.10101010
##  [49] -0.06060606 -0.02020202  0.02020202  0.06060606  0.10101010  0.14141414
##  [55]  0.18181818  0.22222222  0.26262626  0.30303030  0.34343434  0.38383838
##  [61]  0.42424242  0.46464646  0.50505051  0.54545455  0.58585859  0.62626263
##  [67]  0.66666667  0.70707071  0.74747475  0.78787879  0.82828283  0.86868687
##  [73]  0.90909091  0.94949495  0.98989899  1.03030303  1.07070707  1.11111111
##  [79]  1.15151515  1.19191919  1.23232323  1.27272727  1.31313131  1.35353535
##  [85]  1.39393939  1.43434343  1.47474747  1.51515152  1.55555556  1.59595960
##  [91]  1.63636364  1.67676768  1.71717172  1.75757576  1.79797980  1.83838384
##  [97]  1.87878788  1.91919192  1.95959596  2.00000000
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)
2 * x * sin(x) + x^2 * cos(x)
##   [1] 2.2232442755 2.1621237197 2.1001569223 2.0374552454 1.9741295366
##   [6] 1.9102900272 1.8460462307 1.7815068435 1.7167796456 1.6519714028
##  [11] 1.5871877703 1.5225331975 1.4581108337 1.3940224359 1.3303682779
##  [16] 1.2672470605 1.2047558237 1.1429898609 1.0820426339 1.0220056908
##  [21] 0.9629685847 0.9050187953 0.8482416517 0.7927202577 0.7385354186
##  [26] 0.6857655712 0.6344867148 0.5847723450 0.5366933900 0.4903181485
##  [31] 0.4457122306 0.4029385012 0.3620570247 0.3231250141 0.2861967807
##  [36] 0.2513236874 0.2185541047 0.1879333686 0.1595037417 0.1333043769
##  [41] 0.1093712835 0.0877372965 0.0684320482 0.0514819428 0.0369101336
##  [46] 0.0247365036 0.0149776477 0.0076468594 0.0027541183 0.0003060825
##  [51] 0.0003060825 0.0027541183 0.0076468594 0.0149776477 0.0247365036
##  [56] 0.0369101336 0.0514819428 0.0684320482 0.0877372965 0.1093712835
##  [61] 0.1333043769 0.1595037417 0.1879333686 0.2185541047 0.2513236874
##  [66] 0.2861967807 0.3231250141 0.3620570247 0.4029385012 0.4457122306
##  [71] 0.4903181485 0.5366933900 0.5847723450 0.6344867148 0.6857655712
##  [76] 0.7385354186 0.7927202577 0.8482416517 0.9050187953 0.9629685847
##  [81] 1.0220056908 1.0820426339 1.1429898609 1.2047558237 1.2672470605
##  [86] 1.3303682779 1.3940224359 1.4581108337 1.5225331975 1.5871877703
##  [91] 1.6519714028 1.7167796456 1.7815068435 1.8460462307 1.9102900272
##  [96] 1.9741295366 2.0374552454 2.1001569223 2.1621237197 2.2232442755
fxy = expression(x^2+y^2, 'x','y') 
#One can define a function of 2 or more variables
fxy #renders an expression of the function in terms of x and y
## expression(x^2 + y^2, "x", "y")
D(fxy,'x') #yields the partial derivative with respect to x: 2 * x
## 2 * x
D(fxy,'y') #yields the partial derivative with respect to y: 2 * y
## 2 * y
square = function(x) x^2
integrate (square, 0,1) #Integrate x^2 from 0 to 1 equals to 1/3 with details: 0.3333333 with absolute error < 3.7e-15
## 0.3333333 with absolute error < 3.7e-15
integrate(cos,0,pi/2) 
## 1 with absolute error < 1.1e-14
#Integrate cos(x) from 0 to pi/2 equals to 1 with details: 1 with absolute error < 1.1e-14

Vectors and matrices

c(1,6,3,pi,-3) #c() gives a vector, considered a 4X1 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 2 increment 
## [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  # This multiplication * multiples each pair of elements 
## [1]  1 -2  3 -4
x%*%y #This is the dot product of two vectors and yields 
##      [,1]
## [1,]   -2
t(x)  # Transforms x into a row 1X4 vector
##      [,1] [,2] [,3] [,4]
## [1,]    1   -1    1   -1
t(x)%*%y #This is equivalent to dot product and forms 1X1 matrix
##      [,1]
## [1,]   -2
x%*%t(y) #This column times row 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
my=matrix(y,ncol=2) 
#Convert a vector into a matrix of the same number of elements
#The matrix elements go by column, first column, second, etc
#Commands matrix(y,ncol=2, nrow=2)  or matrix(y,2)  
#or matrix(y,2,2) does the same job
my
##      [,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, again via columns
## [1] 1 2 3 4
mx <- matrix(c(1,1,-1,-1), byrow=TRUE,nrow=2)
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 #This is the real matrix multiplication in matrix theory
##      [,1] [,2]
## [1,]    3    7
## [2,]   -3   -7
det(my) #determinant
## [1] -2
myinv = solve(my) #yields the inverse of a matrix
myinv
##      [,1] [,2]
## [1,]   -2  1.5
## [2,]    1 -0.5
myinv%*%my #verifies the inverse of a matrix
##      [,1] [,2]
## [1,]    1    0
## [2,]    0    1
diag(my) #yields the diagonal vector of a matrix
## [1] 1 4
myeig=eigen(my) #yields eigenvalues and unit eigenvectors
myeig
## eigen() decomposition
## $values
## [1]  5.3722813 -0.3722813
## 
## $vectors
##            [,1]       [,2]
## [1,] -0.5657675 -0.9093767
## [2,] -0.8245648  0.4159736
myeig$values
## [1]  5.3722813 -0.3722813
myeig$vectors
##            [,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 a rectangular matrix of mXn
mysvd$d
## [1] 5.4649857 0.3659662
mysvd$u
##            [,1]       [,2]
## [1,] -0.5760484 -0.8174156
## [2,] -0.8174156  0.5760484
mysvd$v
##            [,1]       [,2]
## [1,] -0.4045536  0.9145143
## [2,] -0.9145143 -0.4045536
ysol=solve(my,c(1,3)) 
#solve linear equations matrix %*% x = b
ysol  #solve(matrix, b)
## [1]  2.5 -0.5
my%*%ysol #verifies the solution
##      [,1]
## [1,]    1
## [2,]    3

Simple statistics

x=rnorm(10) #generate 10 normally distributed numbers
x
##  [1] -0.2171938 -0.2911395  1.1124764 -0.5185530 -1.3022137 -0.9662170
##  [7] -0.5040798  0.7256663  1.9473589  0.1538707
mean(x)
## [1] 0.01399754
var(x)
## [1] 0.9847093
sd(x)
## [1] 0.9923252
median(x)
## [1] -0.2541666
quantile(x)
##         0%        25%        50%        75%       100% 
## -1.3022137 -0.5149347 -0.2541666  0.5827174  1.9473589
range(x) #yields the min and max of x
## [1] -1.302214  1.947359
max(x)
## [1] 1.947359
boxplot(x) #yields the box plot of x

w=rnorm(1000)

summary(rnorm(12)) #statistical summary of the data sequence
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -1.94974 -0.36142  0.12761 -0.01368  0.59512  1.00767
hist(w) #yields the histogram of 1000 random numbers with a normal distribution

#Linear regression and linear trend line

#2007-2016 data of the global temperature anomalies
#Source: NOAAGlobalTemp data
t=2007:2016
T=c(.36,.30, .39, .46, .33, .38, .42, .50, .66, .70)
lm(T ~ t) #Linear regression model of temp vs time
## 
## Call:
## lm(formula = T ~ t)
## 
## Coefficients:
## (Intercept)            t  
##   -73.42691      0.03673
#Temperature change rate is 0.03673  deg C/yr or 0.37 deg C/decade
plot(t,T, type="o",xlab="Year",ylab="Temperature [deg C]",
     main="2007-2016 Global Temperature Anomalies 
     and Their Linear Trend [0.37 deg C/decade] ")
abline(lm(T ~ t), lwd=2, col="red") #Regression line

#The R packages and the datasets used in this book are listed below and can be downloaded and installed first before proceeding to the R codes in the rest of the book. 
#The R packages: animation, chron, e1071, fields, ggplot2, lattice, latticeExtra, maps, mapdata, mapproj, matrixStats, ncdf, NLRoot, RColorBrewer, rgdal, rasterVis, raster, sp, TTR

#To load the package "animation", you can do library(animation)

#You can also load all these packages in one shot using pacman
#install.packages("pacman")
library(pacman)
## Warning: package 'pacman' was built under R version 4.4.2
pacman::p_load(animation, chron, e1071, fields, ggplot2, lattice, 
               latticeExtra, maps, mapdata, mapproj, matrixStats, ncdf4, 
               NLRoot, RColorBrewer, rgdal, rasterVis, raster, sp, TTR)
## Installing package into 'C:/Users/alyss/AppData/Local/R/win-library/4.4'
## (as 'lib' is unspecified)
## Warning: package 'NLRoot' is not available for this version of R
## 
## A version of this package for your version of R might be available elsewhere,
## see the ideas at
## https://cran.r-project.org/doc/manuals/r-patched/R-admin.html#Installing-packages
## Warning: unable to access index for repository http://www.stats.ox.ac.uk/pub/RWin/bin/windows/contrib/4.4:
##   cannot open URL 'http://www.stats.ox.ac.uk/pub/RWin/bin/windows/contrib/4.4/PACKAGES'
## Warning: 'BiocManager' not available.  Could not check Bioconductor.
## 
## Please use `install.packages('BiocManager')` and then retry.
## Warning in p_install(package, character.only = TRUE, ...):
## Warning in library(package, lib.loc = lib.loc, character.only = TRUE,
## logical.return = TRUE, : there is no package called 'NLRoot'
## Installing package into 'C:/Users/alyss/AppData/Local/R/win-library/4.4'
## (as 'lib' is unspecified)
## Warning: package 'rgdal' is not available for this version of R
## 
## A version of this package for your version of R might be available elsewhere,
## see the ideas at
## https://cran.r-project.org/doc/manuals/r-patched/R-admin.html#Installing-packages
## Warning: unable to access index for repository http://www.stats.ox.ac.uk/pub/RWin/bin/windows/contrib/4.4:
##   cannot open URL 'http://www.stats.ox.ac.uk/pub/RWin/bin/windows/contrib/4.4/PACKAGES'
## Warning: 'BiocManager' not available.  Could not check Bioconductor.
## 
## Please use `install.packages('BiocManager')` and then retry.
## Warning in p_install(package, character.only = TRUE, ...):
## Warning in library(package, lib.loc = lib.loc, character.only = TRUE,
## logical.return = TRUE, : there is no package called 'rgdal'
## Warning in pacman::p_load(animation, chron, e1071, fields, ggplot2, lattice, : Failed to install/load:
## NLRoot, rgdal

#The zipped data file: https://cambridge.org/climatemathematics/data.zip

#On your computer, you can create a directory called LinAlg under your user name. 
#The one used in the book is Users/sshen/LinAlg  
#You unzip the data and move the data folder under the Users/(your user)/LinAlg directory. 
#A data folder will be created: Users/(your user)/LinAlg/data. 
#The data folder contains about 400 MB of data. 
#Place all the R codes in the directory  Users/(your user)/climmath. 
#Then, you can run all the codes in this book after replacing (your user) by your user name on your own computer.