#R Notes for the Book "Climate Mathematics With R"
#by SSP Shen and RCJ Somerville
#Cambridge University Press, 2018
#
#R Programming for Climate Data Analysis and Visualization
#
#
#Samuel S.P. Shen used these R codes to teach
#Short Courses on R Programming for Climate Science
#at NOAA-NCEI, Asheville, USA in May 2017 and
#Chinese Academy of Sciences, Beijing, China in June 2017
###################################################
#Chapter 2 of the Book "Climate Mathematics With R"
###################################################
#
#Define a sequence
#The same sequence 1,2, ..., 8 can be generated by different commands
1:8
## [1] 1 2 3 4 5 6 7 8
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
#Evaluate the function value
samfctn(4)
## [1] 16
fctn2 <- function(x,y,z) x+y-z/2
fctn2(1,2,3)
## [1] 1.5
#Simple plots
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 calculation
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
## [6] -1.79797980 -1.75757576 -1.71717172 -1.67676768 -1.63636364
## [11] -1.59595960 -1.55555556 -1.51515152 -1.47474747 -1.43434343
## [16] -1.39393939 -1.35353535 -1.31313131 -1.27272727 -1.23232323
## [21] -1.19191919 -1.15151515 -1.11111111 -1.07070707 -1.03030303
## [26] -0.98989899 -0.94949495 -0.90909091 -0.86868687 -0.82828283
## [31] -0.78787879 -0.74747475 -0.70707071 -0.66666667 -0.62626263
## [36] -0.58585859 -0.54545455 -0.50505051 -0.46464646 -0.42424242
## [41] -0.38383838 -0.34343434 -0.30303030 -0.26262626 -0.22222222
## [46] -0.18181818 -0.14141414 -0.10101010 -0.06060606 -0.02020202
## [51] 0.02020202 0.06060606 0.10101010 0.14141414 0.18181818
## [56] 0.22222222 0.26262626 0.30303030 0.34343434 0.38383838
## [61] 0.42424242 0.46464646 0.50505051 0.54545455 0.58585859
## [66] 0.62626263 0.66666667 0.70707071 0.74747475 0.78787879
## [71] 0.82828283 0.86868687 0.90909091 0.94949495 0.98989899
## [76] 1.03030303 1.07070707 1.11111111 1.15151515 1.19191919
## [81] 1.23232323 1.27272727 1.31313131 1.35353535 1.39393939
## [86] 1.43434343 1.47474747 1.51515152 1.55555556 1.59595960
## [91] 1.63636364 1.67676768 1.71717172 1.75757576 1.79797980
## [96] 1.83838384 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
## [6] -1.79797980 -1.75757576 -1.71717172 -1.67676768 -1.63636364
## [11] -1.59595960 -1.55555556 -1.51515152 -1.47474747 -1.43434343
## [16] -1.39393939 -1.35353535 -1.31313131 -1.27272727 -1.23232323
## [21] -1.19191919 -1.15151515 -1.11111111 -1.07070707 -1.03030303
## [26] -0.98989899 -0.94949495 -0.90909091 -0.86868687 -0.82828283
## [31] -0.78787879 -0.74747475 -0.70707071 -0.66666667 -0.62626263
## [36] -0.58585859 -0.54545455 -0.50505051 -0.46464646 -0.42424242
## [41] -0.38383838 -0.34343434 -0.30303030 -0.26262626 -0.22222222
## [46] -0.18181818 -0.14141414 -0.10101010 -0.06060606 -0.02020202
## [51] 0.02020202 0.06060606 0.10101010 0.14141414 0.18181818
## [56] 0.22222222 0.26262626 0.30303030 0.34343434 0.38383838
## [61] 0.42424242 0.46464646 0.50505051 0.54545455 0.58585859
## [66] 0.62626263 0.66666667 0.70707071 0.74747475 0.78787879
## [71] 0.82828283 0.86868687 0.90909091 0.94949495 0.98989899
## [76] 1.03030303 1.07070707 1.11111111 1.15151515 1.19191919
## [81] 1.23232323 1.27272727 1.31313131 1.35353535 1.39393939
## [86] 1.43434343 1.47474747 1.51515152 1.55555556 1.59595960
## [91] 1.63636364 1.67676768 1.71717172 1.75757576 1.79797980
## [96] 1.83838384 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")
#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)
## 0.3333333 with absolute error < 3.7e-15
#Integrate x^2 from 0 to 1 equals to 1/3 with details below
#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 below
#1 with absolute error < 1.1e-14
#Vectors, matrices and solve linear equations
c(1,6,3,pi,-3) #c() gives a vector and is 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
#[1] 2 2
as.vector(my) #Convert a matrix to a vector, again via columns
## [1] 1 2 3 4
#[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
#Basic statistics of a sequence
x=rnorm(10) #generate 10 normally distributed numbers
x
## [1] 0.2164979 1.2602845 -0.1657241 -0.9490608 1.1999920 0.2255151
## [7] -0.4123322 -1.3635690 -0.9687768 0.2872702
mean(x)
## [1] -0.06699032
var(x)
## [1] 0.7844135
sd(x)
## [1] 0.8856712
median(x)
## [1] 0.02538694
quantile(x)
## 0% 25% 50% 75% 100%
## -1.36356905 -0.81487866 0.02538694 0.27183141 1.26028454
range(x) #yields the min and max of x
## [1] -1.363569 1.260285
max(x)
## [1] 1.260285
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.
## -2.02527 -0.55336 0.21350 -0.07699 0.50271 0.95319
hist(w)
#yields the histogram of 1000 random numbers with a normal distribution