R Code for Climate Mathematics:
Theory and Applications
A Cambridge University Press book By
SSP Shen and RCJ Somerville
Version 1.0 released in July 2019 and coded by Dr. Samuel Shen, Distinguished Professor
San Diego State University, California, USA
Email: sshen@sdsu.edu
Version 2.0 compiled and updated by Momtaza Sayd
San Diego State University May 2021
Version 3.0 compiled and updated by Joaquin Stawsky
San Diego State University June 2022
Chapter 4: Climate Data Matrices and Linear Algebra
Matrix algebra
#Matrix subtraction
matrix(c(1,1,1,-1), nrow = 2) - matrix(c(1,3,2,0), nrow = 2)## [,1] [,2]
## [1,] 0 -1
## [2,] -2 -1#Matrix multiplication
matrix(c(1,1,1,-1), nrow = 2) %*% matrix(c(1,2,3,4), nrow = 2)## [,1] [,2]
## [1,] 3 7
## [2,] -1 -1matrix(c(1,2,3,4), nrow = 2) %*% matrix(c(1,1,1,-1), nrow = 2) ## [,1] [,2]
## [1,] 4 -2
## [2,] 6 -2#Find the inverse of a matrix
solve(matrix(c(1,1,1,-1), nrow = 2))## [,1] [,2]
## [1,] 0.5 0.5
## [2,] 0.5 -0.5#Verify the inverse matrix
matrix(c(0.5,0.5,0.5,-0.5), nrow = 2) %*% matrix(c(1,1,1,-1), nrow = 2) ## [,1] [,2]
## [1,] 1 0
## [2,] 0 1#Solve a group of linear equations
solve(matrix(c(1,1,1,-1),nrow = 2),c(20,4))## [1] 12 8#This is the result x1 = 12, and x2 = 8.
#Compute a covariance matrix
A <- matrix(c(1,-1,2,0,3,1),nrow = 2)
A## [,1] [,2] [,3]
## [1,] 1 2 3
## [2,] -1 0 1covm <- (1/(dim(A)[2]))*A%*%t(A)
covm #is the covariance matrix.## [,1] [,2]
## [1,] 4.6666667 0.6666667
## [2,] 0.6666667 0.6666667#Matrix multiplication of a vector changes the vector direction
u <- c(1,-1)
v <- covm%*%u
v## [,1]
## [1,] 4
## [2,] 0#u and v are in different directions.
#Matrix mutiplication of an eigenvector does change direction
ew <- eigen(covm)
eigen(covm)$values## [1] 4.7748518 0.5584816eigen(covm)$vectors## [,1] [,2]
## [1,] -0.9870875 0.1601822
## [2,] -0.1601822 -0.9870875#Verify the eigenvectors and eigenvalues
covm%*%ew$vectors[,1]/ew$values[1]## [,1]
## [1,] -0.9870875
## [2,] -0.1601822#This is the first eigenvector
#SVD: Decomposition and recovery of a space-time data matrix
#Develop a 2-by-3 space-time data matrix for SVD
A <- matrix(c(1,-1,2,0,3,1),nrow = 2)
A## [,1] [,2] [,3]
## [1,] 1 2 3
## [2,] -1 0 1#Perform SVD calculation
msvd <- svd(A)
msvd## $d
## [1] 3.784779 1.294390
##
## $u
## [,1] [,2]
## [1,] -0.9870875 -0.1601822
## [2,] -0.1601822 0.9870875
##
## $v
## [,1] [,2]
## [1,] -0.2184817 -0.8863403
## [2,] -0.5216090 -0.2475023
## [3,] -0.8247362 0.3913356msvd$d## [1] 3.784779 1.294390msvd$u## [,1] [,2]
## [1,] -0.9870875 -0.1601822
## [2,] -0.1601822 0.9870875msvd$v## [,1] [,2]
## [1,] -0.2184817 -0.8863403
## [2,] -0.5216090 -0.2475023
## [3,] -0.8247362 0.3913356#One can verify that A = UDV', where V' is transpose of V.
verim <- msvd$u%*%diag(msvd$d)%*%t(msvd$v)
verim## [,1] [,2] [,3]
## [1,] 1 2.000000e+00 3
## [2,] -1 1.665335e-16 1round(verim)## [,1] [,2] [,3]
## [1,] 1 2 3
## [2,] -1 0 1#This is the original data matrix A
#Covariance matrix and its eigen-problem
covm <- (1/(dim(A)[2]))*A%*%t(A)
eigcov <- eigen(covm)
eigcov$values## [1] 4.7748518 0.5584816eigcov$vectors## [,1] [,2]
## [1,] -0.9870875 0.1601822
## [2,] -0.1601822 -0.9870875#Eigenvalues from SVD and from a covariance matrix
((msvd$d)^2)/(dim(A)[2])## [1] 4.7748518 0.5584816eigcov$values## [1] 4.7748518 0.5584816
Plot Fig. 4.2
#SOI and the Standardized SLP data at Darwin and Tahiti
setwd("~/sshen/climmath")
#setEPS()
#postscript("fig0402a.eps", height = 4.9, width = 7)
#par(mar = c(4.0,4.2,1.5,0.5))
Pta <- read.table("data/PSTANDtahiti.txt", header = F)
# Remove the first column that is the year
ptamon <- Pta[,seq(2,13)]
#Convert the matrix into a vector according to mon: Jan 1951, Feb 1951, ..., Dec 2015
ptamonv <- c(t(ptamon))
#Generate time ticks from Jan 1951 to Dec 2015
xtime <- seq(1951, 2016-1/12, 1/12)
# Plot the Tahiti standardized SLP anomalies
plot(xtime, ptamonv,type = "l",xlab = "Year",
ylab = "Pressure Anomalies", cex.lab = 1.2, cex.axis = 1.2,
main = "Standardized Tahiti SLP Anomalies", col = "red",
xlim = range(xtime), ylim = c(-4,4))
text(1952,3.5, "(a)", cex = 1.2)
# dev.off()
Do the same for Darwin SLP
setwd("~/sshen/climmath/")
#setEPS()
#postscript("fig0402b.eps", height = 4.9, width = 7)
#par(mar = c(4.0,4.2,1.5,0.5))
Pda <- read.table("data/PSTANDdarwin.txt", header = F)
pdamon <- Pda[,seq(2,13)]
pdamonv <- c(t(pdamon))
plot(xtime, pdamonv,type = "l",cex.lab = 1.2, cex.axis = 1.2,
xlab = "Year",ylab = "Pressure",
main = "Standardized Darwin SLP Anomalies", col = "blue",
xlim = range(xtime), ylim = c(-4,4))
text(1952,3.5, "(b)", cex = 1.2)
# dev.off()
Plot the SOI index
#setEPS()
#postscript("fig0402c.eps", height = 4.9, width = 7)
#par(mar = c(4.0,4.2,2.0,0.5))
SOI <- ptamonv-pdamonv
plot(xtime, SOI ,type = "l", cex.lab = 1.2, cex.axis = 1.2,
xlab = "Year",ylab = "SOI index",
col = "black",xlim = range(xtime), ylim = c(-6,6), lwd = 1)
#Add ticks on top edge of the plot box
axis (3, at = seq(1951,2015,4), labels = seq(1951,2015,4))
#Add ticks on the right edge of the plot box
axis (4, at = seq(-4,4,2), labels = seq(-4,4,2))
lines(xtime, rep(0,length(xtime)))
text(1985, 5, "SOI = Tahiti Pressure - Darwin Pressure", cex = 1.2)
text(1952,5.6, "(c)", cex = 1.2)
#abline(lm(SOI ~ xtime), col = "red", lwd = 2)
# dev.off()
CSOI and AMO time series comparison
setwd("~/sshen/climmath/")
# setEPS()
#postscript("fig0402d.eps", height = 4.9, width = 7)
par(mar = c(4,4.2,1.6,4))
cnegsoi <- -cumsum(ptamonv-pdamonv)
#par(mgp = c(2,2,4,0))
plot(xtime, cnegsoi,type = "l", cex.lab = 1.2, cex.axis = 1.2,
xlab = "Year",ylab = "Negative CSOI index",
main = "CSOI and AMO Index Comparison",
col = "purple",xlim = range(xtime), ylim = range(cnegsoi), lwd = 1.5)
legend(1960,15, col = c("purple"),lty = 1,lwd = 2.0,
legend = c("CSOI"),bty = "n",text.font = 2,cex = 1.2)
text(1951,12, "(d)", cex = 1.2)
#AMO data and plot
amodat <- read.table("data/AMO1951-2015.txt", header = FALSE)
amots <- as.vector(t(amodat[,2:13]))
par(new = TRUE)
plot(xtime, amots,type = "l",col = "darkgreen",
cex.lab = 1.2, cex.axis = 1.2,
lwd = 1.5,axes = FALSE,xlab = "",ylab = "")
legend(1960,0.45, col = c("darkgreen"),lty = 1,lwd = 2.0,
legend = c("AMO index"),bty = "n",text.font = 2,cex = 1.2)
#Suppress the axes and assign the y-axis to side 4
axis(4, cex.axis = 1.3)
mtext("AMO index",side = 4,line = 3, cex = 1.2)
# dev.off()
Compute the weighted SOI by SVD
#Time-space data format
ptada <- cbind(ptamonv,pdamonv)
#Space-time data format
ptada <- t(ptada)
#Implement SVD for the sapce-time data
svdptd <- svd(ptada)
recontd <- svdptd$u%*%diag(svdptd$d[1:2])%*%t(svdptd$v)
recontd <- ptada
#Two EOFs from the space-time data
U <- svdptd$u
U## [,1] [,2]
## [1,] -0.6146784 0.7887779
## [2,] 0.7887779 0.6146784#These are the weights for the WSOI
svdptd$d## [1] 31.34582 22.25421D <- diag(svdptd$d)
D## [,1] [,2]
## [1,] 31.34582 0.00000
## [2,] 0.00000 22.25421#PCs: The principal components
V <- svdptd$v
V## [,1] [,2]
## [1,] -5.820531e-02 1.017018e-02
## [2,] -4.026198e-02 -4.419324e-02
## [3,] -2.743069e-03 -8.276652e-02
## [4,] 7.583129e-03 -4.294790e-02
## [5,] 2.379220e-02 -1.376248e-02
## [6,] -1.476952e-02 -9.027043e-02
## [7,] 3.278087e-02 -6.656126e-02
## [8,] 4.772043e-03 -2.324615e-02
## [9,] 3.444712e-02 -4.764183e-02
## [10,] 3.470782e-02 -1.887153e-02
## [11,] 2.630857e-02 -1.100041e-02
## [12,] 3.326827e-02 4.221348e-02
## [13,] 2.941437e-02 -5.316599e-02
## [14,] 2.072041e-02 -2.283104e-02
## [15,] -1.568766e-02 2.835519e-02
## [16,] 1.006549e-02 1.104831e-02
## [17,] -3.160202e-02 -2.329405e-02
## [18,] -2.575316e-02 1.730688e-02
## [19,] -1.790933e-02 3.129287e-03
## [20,] -5.848866e-03 -4.060094e-02
## [21,] 7.288417e-03 -2.048407e-02
## [22,] -1.817003e-02 -2.564102e-02
## [23,] -6.109569e-03 -6.937124e-02
## [24,] 4.421790e-02 -1.412969e-02
## [25,] -8.693959e-03 3.033495e-02
## [26,] 1.790933e-02 -3.129287e-03
## [27,] 4.216628e-03 -2.955263e-02
## [28,] -1.088161e-02 -4.612509e-02
## [29,] 6.323806e-02 -4.646029e-03
## [30,] -6.143578e-03 -1.813711e-02
## [31,] 2.221662e-03 2.522591e-02
## [32,] 3.755291e-02 -8.980741e-02
## [33,] 4.617886e-02 -1.767409e-02
## [34,] -1.960958e-03 3.544399e-03
## [35,] 1.258187e-02 1.381039e-02
## [36,] 1.705920e-02 1.302806e-02
## [37,] -3.045719e-02 -6.191523e-02
## [38,] 1.006549e-02 1.104831e-02
## [39,] -1.313728e-02 -2.011686e-02
## [40,] -3.267885e-02 -8.714114e-02
## [41,] -2.290806e-02 -5.362900e-02
## [42,] -7.254409e-03 -3.075006e-02
## [43,] -1.676449e-02 -3.549189e-02
## [44,] -4.307306e-02 -2.449149e-02
## [45,] -7.549120e-03 -8.286232e-03
## [46,] -1.036021e-02 1.141552e-02
## [47,] -2.811085e-03 1.970175e-02
## [48,] -4.954536e-02 3.106937e-02
## [49,] 2.597985e-02 6.269755e-02
## [50,] -7.634133e-02 -6.670497e-02
## [51,] -2.990177e-02 -5.560876e-02
## [52,] -1.076823e-03 -6.384709e-02
## [53,] -4.388918e-02 -5.956827e-02
## [54,] -5.117760e-02 -3.908419e-02
## [55,] -5.987156e-02 -8.749246e-03
## [56,] -5.706047e-02 -2.845100e-02
## [57,] -4.784510e-02 -1.245336e-03
## [58,] -5.624435e-02 6.625785e-03
## [59,] -4.784510e-02 -1.245336e-03
## [60,] -3.026450e-02 6.932334e-02
## [61,] -5.258314e-02 -2.923332e-02
## [62,] -6.153780e-02 -2.766868e-02
## [63,] -5.036148e-02 -4.007413e-03
## [64,] -3.719019e-02 -3.512468e-02
## [65,] -5.928213e-02 -5.367690e-02
## [66,] -4.362848e-02 -3.079796e-02
## [67,] -4.640555e-02 -6.233034e-02
## [68,] -4.866122e-02 -3.632212e-02
## [69,] -6.404281e-03 -4.690741e-02
## [70,] -6.438290e-02 4.326721e-02
## [71,] -6.993705e-03 -1.979755e-03
## [72,] -4.222293e-02 -4.064884e-02
## [73,] -2.072041e-02 2.283104e-02
## [74,] -2.607036e-04 -2.877031e-02
## [75,] -1.062091e-02 -1.735479e-02
## [76,] -1.202645e-02 -7.503910e-03
## [77,] 2.545845e-02 5.156945e-03
## [78,] -8.104536e-03 -1.459271e-02
## [79,] -1.062091e-02 -1.735479e-02
## [80,] 1.513225e-02 -3.466167e-02
## [81,] 3.385769e-02 -2.714175e-03
## [82,] 3.887908e-03 4.414534e-02
## [83,] 3.581865e-02 -6.258575e-03
## [84,] 1.147104e-02 1.197434e-03
## [85,] 7.219272e-02 -6.210673e-03
## [86,] 1.679850e-02 -1.574224e-02
## [87,] -8.954663e-03 1.564644e-03
## [88,] -1.676449e-02 -3.549189e-02
## [89,] 1.291059e-02 -5.988758e-02
## [90,] -1.450882e-02 -6.150012e-02
## [91,] -1.676449e-02 -3.549189e-02
## [92,] -3.830102e-02 -4.773764e-02
## [93,] 5.916882e-03 -6.186733e-02
## [94,] -1.003149e-02 -6.228244e-02
## [95,] 1.679850e-02 -1.574224e-02
## [96,] 2.405290e-02 1.500782e-02
## [97,] 2.663729e-02 -8.469837e-02
## [98,] 5.454410e-02 2.568892e-02
## [99,] -4.840052e-02 -7.551812e-03
## [100,] -1.872545e-02 -3.194749e-02
## [101,] -2.124182e-02 -3.470957e-02
## [102,] 3.366501e-03 -1.339527e-02
## [103,] 1.036021e-02 -1.141552e-02
## [104,] 2.947118e-04 -2.246383e-02
## [105,] -1.666246e-03 -1.891943e-02
## [106,] -1.679850e-02 1.574224e-02
## [107,] -3.552394e-02 -1.620525e-02
## [108,] -3.607936e-02 -2.251173e-02
## [109,] -6.993705e-03 -1.979755e-03
## [110,] -2.516373e-03 -2.762077e-03
## [111,] -3.944586e-02 -9.116456e-03
## [112,] -3.467381e-02 -3.236261e-02
## [113,] -2.542444e-02 -5.639108e-02
## [114,] -6.143578e-03 -1.813711e-02
## [115,] -2.013099e-02 -2.209662e-02
## [116,] -2.993578e-02 -4.374623e-03
## [117,] -3.045719e-02 -6.191523e-02
## [118,] -4.477332e-03 7.823218e-04
## [119,] -1.875945e-02 1.928664e-02
## [120,] -2.660329e-02 3.346424e-02
## [121,] 1.062091e-02 1.735479e-02
## [122,] -3.163603e-02 2.794008e-02
## [123,] 6.297735e-02 -3.341633e-02
## [124,] -2.882495e-02 8.238330e-03
## [125,] -1.398741e-02 -3.959511e-03
## [126,] -2.516373e-03 -2.762077e-03
## [127,] -9.510078e-03 -4.741833e-03
## [128,] -6.993705e-03 -1.979755e-03
## [129,] -2.516373e-03 -2.762077e-03
## [130,] 1.009950e-02 -4.018582e-02
## [131,] -1.402142e-02 4.727462e-02
## [132,] -5.039549e-02 4.722672e-02
## [133,] -7.052647e-02 2.513010e-02
## [134,] 1.258187e-02 1.381039e-02
## [135,] -3.366501e-03 1.339527e-02
## [136,] -6.438290e-03 4.326721e-03
## [137,] -4.251765e-02 -1.818501e-02
## [138,] -3.075190e-02 -3.945140e-02
## [139,] -7.515112e-03 -5.952037e-02
## [140,] -2.212595e-02 3.268191e-02
## [141,] -1.928086e-02 -3.825397e-02
## [142,] -3.889044e-02 -2.809979e-03
## [143,] -1.369270e-02 -2.642334e-02
## [144,] -6.438290e-03 4.326721e-03
## [145,] -3.722419e-02 1.610945e-02
## [146,] -1.153905e-02 1.012708e-01
## [147,] -4.558943e-02 -2.725356e-02
## [148,] -3.467381e-02 -3.236261e-02
## [149,] -1.787532e-02 -4.810485e-02
## [150,] 1.735391e-02 -9.435764e-03
## [151,] 4.477332e-03 -7.823218e-04
## [152,] -4.998739e-03 -5.675829e-02
## [153,] 2.183124e-02 -1.021809e-02
## [154,] 4.755039e-02 2.370916e-02
## [155,] 3.467381e-02 3.236261e-02
## [156,] 4.951135e-02 2.016476e-02
## [157,] 1.509824e-02 1.657246e-02
## [158,] -8.501272e-04 1.615735e-02
## [159,] -4.392319e-02 -8.334134e-03
## [160,] -3.918515e-02 1.965385e-02
## [161,] -1.284257e-02 -4.258069e-02
## [162,] -2.908565e-02 -2.053197e-02
## [163,] -2.486902e-02 -5.008460e-02
## [164,] -5.568894e-02 1.293226e-02
## [165,] -5.062218e-02 -3.277772e-02
## [166,] -4.728969e-02 5.061141e-03
## [167,] -1.143703e-02 -5.243157e-02
## [168,] 1.091562e-02 -5.109043e-03
## [169,] 1.206046e-02 -4.373022e-02
## [170,] -1.928086e-02 -3.825397e-02
## [171,] -3.075190e-02 -3.945140e-02
## [172,] 2.431361e-02 4.377813e-02
## [173,] -7.549120e-03 -8.286232e-03
## [174,] 2.516373e-02 2.762077e-02
## [175,] 6.938163e-02 1.349108e-02
## [176,] 2.323678e-02 -2.006896e-02
## [177,] 4.869523e-02 -1.491202e-02
## [178,] 3.774560e-02 4.143116e-02
## [179,] 5.679977e-02 -3.193082e-04
## [180,] -9.770782e-03 -3.351214e-02
## [181,] 4.702899e-02 -3.383145e-02
## [182,] 7.288417e-03 -2.048407e-02
## [183,] 3.052520e-02 -4.055303e-02
## [184,] 5.882874e-03 -1.063320e-02
## [185,] 1.565366e-02 2.287894e-02
## [186,] -8.138544e-03 3.664143e-02
## [187,] -4.738035e-03 -2.798798e-02
## [188,] -2.545845e-02 -5.156945e-03
## [189,] 7.288417e-03 -2.048407e-02
## [190,] 6.698993e-03 2.444358e-02
## [191,] 0.000000e+00 0.000000e+00
## [192,] 1.287658e-02 -8.653442e-03
## [193,] -6.323806e-02 4.646029e-03
## [194,] -5.905544e-02 2.632754e-02
## [195,] -4.477332e-02 7.823218e-03
## [196,] 3.887908e-03 4.414534e-02
## [197,] -2.607036e-04 -2.877031e-02
## [198,] -2.379220e-02 1.376248e-02
## [199,] -9.215367e-03 -2.720566e-02
## [200,] -2.630857e-02 1.100041e-02
## [201,] -1.346600e-02 5.358110e-02
## [202,] -1.960958e-03 3.544399e-03
## [203,] 1.372671e-02 -2.481079e-02
## [204,] 2.183124e-02 -1.021809e-02
## [205,] -2.846222e-02 -1.166938e-01
## [206,] -5.202772e-02 -2.292684e-02
## [207,] -2.516373e-03 -2.762077e-03
## [208,] -3.071789e-03 -9.068554e-03
## [209,] -5.117760e-02 -3.908419e-02
## [210,] -3.918515e-02 1.965385e-02
## [211,] -2.545845e-02 -5.156945e-03
## [212,] -1.117632e-02 -2.366126e-02
## [213,] 9.804790e-03 -1.772200e-02
## [214,] 3.071789e-03 9.068554e-03
## [215,] 1.620907e-02 2.918542e-02
## [216,] -6.438290e-03 4.326721e-03
## [217,] 4.261967e-02 -1.355174e-01
## [218,] 1.487155e-02 -6.343197e-02
## [219,] -1.372671e-02 2.481079e-02
## [220,] 1.676449e-02 3.549189e-02
## [221,] 5.622171e-03 -3.940350e-02
## [222,] -6.993705e-03 -1.979755e-03
## [223,] 1.624308e-02 -2.204872e-02
## [224,] 3.921916e-03 -7.088798e-03
## [225,] 2.915367e-02 -8.193629e-02
## [226,] 3.526324e-02 -1.256505e-02
## [227,] 3.071789e-03 9.068554e-03
## [228,] -1.398741e-02 -3.959511e-03
## [229,] 4.503402e-02 2.094709e-02
## [230,] 3.500253e-02 -4.133536e-02
## [231,] -2.320278e-02 -3.116517e-02
## [232,] 2.516373e-03 2.762077e-03
## [233,] -1.620907e-02 -2.918542e-02
## [234,] -3.692948e-02 -6.354378e-03
## [235,] 1.343199e-02 -2.346965e-03
## [236,] -2.682998e-02 -4.654020e-02
## [237,] -4.781110e-02 -5.247947e-02
## [238,] -4.585014e-02 -5.602387e-02
## [239,] -5.820531e-02 1.017018e-02
## [240,] -7.385897e-02 -1.270876e-02
## [241,] -1.620907e-02 -2.918542e-02
## [242,] -7.637534e-02 -1.547083e-02
## [243,] -8.618013e-02 2.251162e-03
## [244,] -6.434889e-02 -7.966924e-03
## [245,] -3.300757e-02 -1.344318e-02
## [246,] -1.121033e-02 2.757287e-02
## [247,] -1.032620e-02 -3.981861e-02
## [248,] -5.735518e-02 -5.987168e-03
## [249,] -5.232244e-02 -4.630137e-04
## [250,] -5.905544e-02 2.632754e-02
## [251,] -1.846474e-02 -3.177189e-03
## [252,] -1.561965e-02 -7.411307e-02
## [253,] -2.013099e-02 -2.209662e-02
## [254,] -4.418389e-02 -3.710444e-02
## [255,] -2.990177e-02 -5.560876e-02
## [256,] 7.254409e-03 3.075006e-02
## [257,] 6.209322e-02 3.397515e-02
## [258,] 1.709321e-02 -3.820607e-02
## [259,] 4.872924e-02 -6.614615e-02
## [260,] 1.790933e-02 -3.129287e-03
## [261,] 5.036148e-02 4.007413e-03
## [262,] 3.885643e-02 5.404411e-02
## [263,] 1.620907e-02 2.918542e-02
## [264,] 5.228843e-02 5.169715e-02
## [265,] 1.313728e-02 2.011686e-02
## [266,] 5.036148e-02 4.007413e-03
## [267,] -2.823552e-02 -3.668933e-02
## [268,] -6.404281e-03 -4.690741e-02
## [269,] -1.705920e-02 -1.302806e-02
## [270,] -4.529472e-02 -4.971739e-02
## [271,] -2.460832e-02 -2.131430e-02
## [272,] -5.062218e-02 -3.277772e-02
## [273,] -4.366248e-02 2.043617e-02
## [274,] -2.938036e-02 1.931854e-03
## [275,] -1.004282e-01 -3.047865e-02
## [276,] -6.271665e-02 6.218664e-02
## [277,] -9.176829e-02 -9.579469e-03
## [278,] -7.889171e-02 -1.823291e-02
## [279,] -9.287912e-02 -2.219242e-02
## [280,] -3.692948e-02 -6.354378e-03
## [281,] -3.163603e-02 2.794008e-02
## [282,] -1.650378e-02 -6.721588e-03
## [283,] -4.170153e-02 1.689177e-02
## [284,] -3.411840e-02 -2.605613e-02
## [285,] -4.140681e-02 -5.572056e-03
## [286,] -3.748490e-02 -1.266085e-02
## [287,] -2.187653e-03 -7.646004e-02
## [288,] -1.450882e-02 -6.150012e-02
## [289,] 1.402142e-02 -4.727462e-02
## [290,] -2.493704e-02 5.238367e-02
## [291,] -5.987156e-02 -8.749246e-03
## [292,] -4.977205e-02 -4.893507e-02
## [293,] -2.712469e-02 -2.407637e-02
## [294,] -4.895593e-02 -1.385829e-02
## [295,] -7.163730e-02 1.251715e-02
## [296,] -7.526451e-02 -2.857881e-03
## [297,] -7.611464e-02 1.329947e-02
## [298,] -6.434889e-02 -7.966924e-03
## [299,] -4.170153e-02 1.689177e-02
## [300,] -7.611464e-02 1.329947e-02
## [301,] -5.091689e-02 -1.031389e-02
## [302,] -6.212723e-02 1.725898e-02
## [303,] -7.075317e-02 -5.487434e-02
## [304,] -9.804790e-03 1.772200e-02
## [305,] -8.988671e-03 5.279878e-02
## [306,] -9.510078e-03 -4.741833e-03
## [307,] 3.023049e-02 -1.808921e-02
## [308,] 2.938036e-02 -1.931854e-03
## [309,] 4.059070e-02 -2.950473e-02
## [310,] -1.817003e-02 -2.564102e-02
## [311,] -2.908565e-02 -2.053197e-02
## [312,] 7.322425e-03 -7.171821e-02
## [313,] 9.838799e-03 -6.895613e-02
## [314,] -4.532873e-02 1.516742e-03
## [315,] 1.653779e-02 -4.451255e-02
## [316,] 1.094963e-02 -5.634318e-02
## [317,] 1.457683e-02 -4.096815e-02
## [318,] 3.219145e-02 -2.163360e-02
## [319,] 3.836903e-02 -5.473063e-02
## [320,] 2.660329e-02 -3.346424e-02
## [321,] 2.882495e-02 -8.238330e-03
## [322,] 4.111210e-02 2.803589e-02
## [323,] 4.673427e-02 -1.136762e-02
## [324,] 3.666878e-02 -2.241593e-02
## [325,] 1.954156e-02 6.702428e-02
## [326,] 9.683505e-02 -3.613051e-02
## [327,] 1.395340e-02 5.519364e-02
## [328,] 1.091562e-02 -5.109043e-03
## [329,] -5.621034e-02 -4.460835e-02
## [330,] -3.101260e-02 -6.822171e-02
## [331,] -2.820152e-02 -8.792346e-02
## [332,] -2.094711e-02 -5.717340e-02
## [333,] -5.293450e-03 -3.429446e-02
## [334,] 1.121033e-02 -2.757287e-02
## [335,] 1.110831e-03 1.261295e-02
## [336,] 1.960958e-03 -3.544399e-03
## [337,] 1.039421e-02 -6.264965e-02
## [338,] -2.885896e-02 5.947246e-02
## [339,] -6.143578e-03 -1.813711e-02
## [340,] 2.811085e-03 -1.970175e-02
## [341,] -1.735391e-02 9.435764e-03
## [342,] -2.127583e-02 1.652456e-02
## [343,] -5.480480e-02 -5.445922e-02
## [344,] 5.327459e-03 -1.693967e-02
## [345,] -6.698993e-03 -2.444358e-02
## [346,] 5.032747e-03 5.524155e-03
## [347,] 1.679850e-02 -1.574224e-02
## [348,] 2.908565e-02 2.053197e-02
## [349,] -1.761461e-02 -1.933454e-02
## [350,] -1.258187e-02 -1.381039e-02
## [351,] 8.172552e-03 -8.787556e-02
## [352,] 1.849875e-02 -4.805694e-02
## [353,] -1.110831e-03 -1.261295e-02
## [354,] -2.221662e-03 -2.522591e-02
## [355,] 8.501272e-04 -1.615735e-02
## [356,] -1.176575e-02 2.126639e-02
## [357,] 1.790933e-02 -3.129287e-03
## [358,] -5.554154e-04 -6.306476e-03
## [359,] 1.287658e-02 -8.653442e-03
## [360,] 1.405543e-03 -9.850876e-03
## [361,] -1.954156e-02 -6.702428e-02
## [362,] 5.894237e-04 -4.492766e-02
## [363,] 4.895593e-02 1.385829e-02
## [364,] 3.921916e-03 -7.088798e-03
## [365,] -3.049119e-02 -1.068110e-02
## [366,] -4.643956e-02 -1.109621e-02
## [367,] -3.578464e-02 -4.497556e-02
## [368,] -2.853024e-02 -1.422550e-02
## [369,] -1.565366e-02 -2.287894e-02
## [370,] 1.343199e-02 -2.346965e-03
## [371,] -8.659951e-03 -2.089918e-02
## [372,] -1.513225e-02 3.466167e-02
## [373,] -3.892445e-02 4.842415e-02
## [374,] -6.733001e-03 2.679055e-02
## [375,] -2.405290e-02 -1.500782e-02
## [376,] -8.501272e-04 1.615735e-02
## [377,] 1.398741e-02 3.959511e-03
## [378,] 3.833503e-02 -3.496497e-03
## [379,] 5.457811e-02 -2.554521e-02
## [380,] 6.353277e-02 -2.710986e-02
## [381,] 6.575443e-02 -1.883952e-03
## [382,] 6.379347e-02 1.660447e-03
## [383,] 9.542950e-02 -2.627963e-02
## [384,] 8.199751e-02 -2.393267e-02
## [385,] 1.292872e-01 -2.899381e-02
## [386,] 1.348754e-01 -1.716318e-02
## [387,] 8.732497e-02 -4.087234e-02
## [388,] 3.052520e-02 -4.055303e-02
## [389,] -1.653779e-02 4.451255e-02
## [390,] -2.947118e-04 2.246383e-02
## [391,] 2.460832e-02 2.131430e-02
## [392,] -3.400828e-05 5.123413e-02
## [393,] -2.800883e-02 4.331511e-02
## [394,] -1.483754e-02 1.219784e-02
## [395,] 7.809824e-03 3.705654e-02
## [396,] 5.554154e-03 6.306476e-02
## [397,] -3.661213e-03 3.585910e-02
## [398,] -3.522923e-02 -3.866908e-02
## [399,] 6.733001e-03 -2.679055e-02
## [400,] -1.065492e-02 3.387935e-02
## [401,] -1.006549e-02 -1.104831e-02
## [402,] 1.091562e-02 -5.109043e-03
## [403,] -6.993705e-03 -1.979755e-03
## [404,] -1.983628e-02 -4.456045e-02
## [405,] -6.143578e-03 -1.813711e-02
## [406,] 1.454283e-02 1.026599e-02
## [407,] -5.622171e-03 3.940350e-02
## [408,] -9.702765e-03 -1.359804e-01
## [409,] 1.369270e-02 2.642334e-02
## [410,] -4.643956e-02 -1.109621e-02
## [411,] -3.463980e-02 -8.359674e-02
## [412,] -5.143830e-02 -6.785450e-02
## [413,] -2.427960e-02 -9.501226e-02
## [414,] 1.009950e-02 -4.018582e-02
## [415,] 5.032747e-03 5.524155e-03
## [416,] -3.637407e-02 -4.790184e-05
## [417,] -3.887908e-03 -4.414534e-02
## [418,] 1.121033e-02 -2.757287e-02
## [419,] 1.144839e-03 -3.862118e-02
## [420,] -6.733001e-03 2.679055e-02
## [421,] -4.137281e-02 -5.680619e-02
## [422,] 3.389170e-02 -5.394831e-02
## [423,] -1.317129e-02 3.111727e-02
## [424,] -1.258187e-02 -1.381039e-02
## [425,] 8.399248e-03 -7.871120e-03
## [426,] -4.333376e-02 -5.326179e-02
## [427,] -1.313728e-02 -2.011686e-02
## [428,] 4.545348e-03 -1.032506e-01
## [429,] 1.568766e-02 -2.835519e-02
## [430,] -2.101512e-02 4.529487e-02
## [431,] 4.144082e-02 -4.566208e-02
## [432,] 5.872672e-02 4.737043e-02
## [433,] 2.379220e-02 -1.376248e-02
## [434,] 4.395720e-02 -4.290000e-02
## [435,] 5.284384e-02 5.800362e-02
## [436,] 5.343327e-02 1.307597e-02
## [437,] 4.673427e-02 -1.136762e-02
## [438,] 3.640808e-02 -5.118623e-02
## [439,] 4.925065e-02 -8.605540e-03
## [440,] 3.385769e-02 -2.714175e-03
## [441,] 3.944586e-02 9.116456e-03
## [442,] 1.398741e-02 3.959511e-03
## [443,] 0.000000e+00 0.000000e+00
## [444,] 1.790933e-02 -3.129287e-03
## [445,] 9.476070e-03 5.597597e-02
## [446,] 1.317129e-02 -3.111727e-02
## [447,] -2.934635e-02 -4.930228e-02
## [448,] -2.811085e-03 1.970175e-02
## [449,] -3.274686e-02 1.532713e-02
## [450,] 3.037781e-03 6.030269e-02
## [451,] -3.666878e-02 2.241593e-02
## [452,] -5.346727e-02 3.815817e-02
## [453,] -6.715997e-02 1.173483e-02
## [454,] -5.343327e-02 -1.307597e-02
## [455,] -6.268264e-02 1.095251e-02
## [456,] -4.366248e-02 2.043617e-02
## [457,] -5.764989e-02 1.647666e-02
## [458,] -3.921916e-02 7.088798e-02
## [459,] -3.833503e-02 3.496497e-03
## [460,] -6.460959e-02 -3.673723e-02
## [461,] -4.840052e-02 -7.551812e-03
## [462,] -2.408691e-02 3.622631e-02
## [463,] -3.411840e-02 -2.605613e-02
## [464,] 7.288417e-03 -2.048407e-02
## [465,] -1.846474e-02 -3.177189e-03
## [466,] -2.993578e-02 -4.374623e-03
## [467,] 6.733001e-03 -2.679055e-02
## [468,] 2.542444e-02 5.639108e-02
## [469,] 3.071789e-03 9.068554e-03
## [470,] 7.493579e-02 7.655584e-02
## [471,] 1.568766e-02 -2.835519e-02
## [472,] -6.993705e-03 -1.979755e-03
## [473,] -4.532873e-02 1.516742e-03
## [474,] -1.258187e-02 -1.381039e-02
## [475,] -2.098111e-02 -5.939266e-03
## [476,] 6.438290e-03 -4.326721e-03
## [477,] 2.908565e-02 2.053197e-02
## [478,] -9.249375e-03 2.402847e-02
## [479,] 1.679850e-02 -1.574224e-02
## [480,] 9.510078e-03 4.741833e-03
## [481,] -2.127583e-02 1.652456e-02
## [482,] -9.510078e-03 -4.741833e-03
## [483,] 2.964107e-02 2.683845e-02
## [484,] 1.960958e-02 -3.544399e-02
## [485,] 3.889044e-02 2.809979e-03
## [486,] 5.588162e-03 1.183063e-02
## [487,] 0.000000e+00 0.000000e+00
## [488,] 1.705920e-02 1.302806e-02
## [489,] 5.457811e-02 -2.554521e-02
## [490,] 3.500253e-02 -4.133536e-02
## [491,] 2.768011e-02 3.038285e-02
## [492,] 6.242194e-02 -3.972281e-02
## [493,] 1.093829e-01 2.891401e-02
## [494,] 3.385769e-02 -2.714175e-03
## [495,] 7.778088e-02 5.619958e-03
## [496,] 3.889044e-02 2.809979e-03
## [497,] -1.202645e-02 -7.503910e-03
## [498,] 1.738792e-02 -6.066990e-02
## [499,] 2.682998e-02 4.654020e-02
## [500,] -1.454283e-02 -1.026599e-02
## [501,] -6.959697e-03 -5.321389e-02
## [502,] 4.732370e-02 -5.629527e-02
## [503,] 2.434762e-02 -7.456008e-03
## [504,] 1.709321e-02 -3.820607e-02
## [505,] 2.996979e-02 -4.685951e-02
## [506,] 2.130984e-02 -6.775870e-02
## [507,] 2.457431e-02 7.254843e-02
## [508,] 4.699498e-02 1.740269e-02
## [509,] 1.424811e-02 3.272982e-02
## [510,] 3.245215e-02 7.136700e-03
## [511,] 2.800883e-02 -4.331511e-02
## [512,] 4.026198e-02 4.419324e-02
## [513,] 2.964107e-02 2.683845e-02
## [514,] 4.140681e-02 5.572056e-03
## [515,] 2.516373e-03 2.762077e-03
## [516,] -1.006549e-02 -1.104831e-02
## [517,] 3.921916e-03 -7.088798e-03
## [518,] -9.510078e-03 -4.741833e-03
## [519,] 2.712469e-02 2.407637e-02
## [520,] 4.562344e-02 -2.398057e-02
## [521,] 2.434762e-02 -7.456008e-03
## [522,] 1.568766e-02 -2.835519e-02
## [523,] 4.647357e-02 -4.013792e-02
## [524,] 4.699498e-02 1.740269e-02
## [525,] 5.709448e-02 -2.278314e-02
## [526,] 4.281236e-02 -4.278819e-03
## [527,] 2.209195e-02 1.855222e-02
## [528,] 4.866122e-02 3.632212e-02
## [529,] 1.539295e-02 -5.891365e-03
## [530,] 8.920655e-03 4.966949e-02
## [531,] -2.908565e-02 -2.053197e-02
## [532,] 1.990429e-02 -5.790782e-02
## [533,] 1.232116e-02 -1.495992e-02
## [534,] -3.366501e-03 1.339527e-02
## [535,] -1.957557e-02 -1.579014e-02
## [536,] -1.287658e-02 8.653442e-03
## [537,] -8.954663e-03 1.564644e-03
## [538,] -2.777077e-03 -3.153238e-02
## [539,] -2.516373e-03 -2.762077e-03
## [540,] 1.624308e-02 -2.204872e-02
## [541,] -3.833503e-02 3.496497e-03
## [542,] -1.313728e-02 -2.011686e-02
## [543,] -3.974057e-02 1.334737e-02
## [544,] -2.993578e-02 -4.374623e-03
## [545,] -6.472298e-03 5.556085e-02
## [546,] -4.199624e-02 3.935560e-02
## [547,] -2.405290e-02 -1.500782e-02
## [548,] -3.049119e-02 -1.068110e-02
## [549,] -3.156802e-02 -7.452819e-02
## [550,] -2.516373e-02 -2.762077e-02
## [551,] -8.161190e-04 -3.507678e-02
## [552,] -3.830102e-02 -4.773764e-02
## [553,] -1.235517e-02 6.619405e-02
## [554,] -6.490430e-02 -1.427340e-02
## [555,] 1.483754e-02 -1.219784e-02
## [556,] 2.960706e-02 7.807259e-02
## [557,] 4.536274e-02 -5.275088e-02
## [558,] 5.398868e-02 1.938244e-02
## [559,] 2.993578e-02 4.374623e-03
## [560,] 5.594964e-02 1.583804e-02
## [561,] 5.732117e-02 5.722130e-02
## [562,] 6.235392e-02 6.274546e-02
## [563,] 5.395467e-02 7.061658e-02
## [564,] 3.748490e-02 1.266085e-02
## [565,] 9.183631e-02 -9.288880e-02
## [566,] 8.389045e-02 7.499120e-02
## [567,] 8.676955e-02 -4.717882e-02
## [568,] 5.317256e-02 -1.569434e-02
## [569,] -7.843832e-03 1.417760e-02
## [570,] -3.581865e-02 6.258575e-03
## [571,] -4.366248e-02 2.043617e-02
## [572,] -4.059070e-02 2.950473e-02
## [573,] -3.666878e-02 2.241593e-02
## [574,] -4.140681e-02 -5.572056e-03
## [575,] -4.166752e-02 -3.434236e-02
## [576,] -5.287785e-02 -6.769490e-03
## [577,] -6.771539e-02 5.428351e-03
## [578,] -3.470782e-02 1.887153e-02
## [579,] -5.117760e-02 -3.908419e-02
## [580,] -5.121160e-02 1.214994e-02
## [581,] -4.511340e-03 5.201646e-02
## [582,] -6.733001e-03 2.679055e-02
## [583,] -1.875945e-02 1.928664e-02
## [584,] -6.472298e-03 5.556085e-02
## [585,] 5.554154e-04 6.306476e-03
## [586,] -3.359699e-02 3.148448e-02
## [587,] -3.777961e-02 9.802974e-03
## [588,] -4.898994e-02 3.737585e-02
## [589,] -2.157054e-02 3.898839e-02
## [590,] -5.405670e-02 8.308582e-02
## [591,] -4.758440e-02 2.752497e-02
## [592,] -4.588415e-02 -4.789735e-03
## [593,] -1.065492e-02 3.387935e-02
## [594,] 8.659951e-03 2.089918e-02
## [595,] 4.216628e-03 -2.955263e-02
## [596,] -2.964107e-02 -2.683845e-02
## [597,] -2.856424e-02 3.700864e-02
## [598,] -4.362848e-02 -3.079796e-02
## [599,] -7.300884e-02 -2.886611e-02
## [600,] -4.022797e-02 -9.542737e-02
## [601,] -2.859825e-02 8.824277e-02
## [602,] -7.156929e-02 -8.995112e-02
## [603,] -3.163603e-02 2.794008e-02
## [604,] -8.104536e-03 -1.459271e-02
## [605,] 1.902016e-02 9.483666e-03
## [606,] -1.565366e-02 -2.287894e-02
## [607,] 7.843832e-03 -1.417760e-02
## [608,] 1.317129e-02 -3.111727e-02
## [609,] -8.104536e-03 -1.459271e-02
## [610,] -2.482365e-03 -5.399621e-02
## [611,] -2.630857e-02 1.100041e-02
## [612,] 2.911966e-02 -3.070216e-02
## [613,] -1.094963e-02 5.634318e-02
## [614,] -3.444712e-02 4.764183e-02
## [615,] 1.173174e-02 2.996774e-02
## [616,] 2.516373e-03 2.762077e-03
## [617,] 3.078591e-02 -1.178273e-02
## [618,] 8.954663e-03 -1.564644e-03
## [619,] 2.375819e-02 3.747165e-02
## [620,] 3.304158e-02 -3.779096e-02
## [621,] 2.349749e-02 8.701344e-03
## [622,] 1.261588e-02 -3.742375e-02
## [623,] 1.817003e-02 2.564102e-02
## [624,] 4.362848e-02 3.079796e-02
## [625,] 6.993705e-03 1.979755e-03
## [626,] 2.045971e-02 -5.160134e-02
## [627,] 8.399248e-03 -7.871120e-03
## [628,] 9.215367e-03 2.720566e-02
## [629,] 9.804790e-03 -1.772200e-02
## [630,] 2.016500e-02 -2.913752e-02
## [631,] -9.804790e-03 1.772200e-02
## [632,] -3.627204e-03 -1.537503e-02
## [633,] 3.366501e-03 -1.339527e-02
## [634,] -2.221662e-03 -2.522591e-02
## [635,] 1.202645e-02 7.503910e-03
## [636,] -4.336777e-02 -2.027657e-03
## [637,] 4.758440e-02 -2.752497e-02
## [638,] -4.451261e-02 3.659352e-02
## [639,] -2.150252e-02 -6.347988e-02
## [640,] 3.078591e-02 -1.178273e-02
## [641,] -3.777961e-02 9.802974e-03
## [642,] 3.522923e-02 3.866908e-02
## [643,] 1.790933e-02 -3.129287e-03
## [644,] 1.036021e-02 -1.141552e-02
## [645,] 1.620907e-02 2.918542e-02
## [646,] 3.627204e-03 1.537503e-02
## [647,] 2.323678e-02 -2.006896e-02
## [648,] 2.493704e-02 -5.238367e-02
## [649,] -1.147104e-02 -1.197434e-03
## [650,] 1.130781e-01 -5.817923e-02
## [651,] -1.121033e-02 2.757287e-02
## [652,] 2.571915e-02 3.392725e-02
## [653,] 2.741940e-02 1.612546e-03
## [654,] -2.013099e-02 -2.209662e-02
## [655,] -6.993705e-03 -1.979755e-03
## [656,] 1.036021e-02 -1.141552e-02
## [657,] -1.343199e-02 2.346965e-03
## [658,] -4.366248e-02 2.043617e-02
## [659,] 3.105797e-03 -4.216558e-02
## [660,] 3.627204e-03 1.537503e-02
## [661,] -6.294334e-02 -1.781780e-02
## [662,] -1.144839e-03 3.862118e-02
## [663,] -6.742068e-02 -1.703548e-02
## [664,] -4.640555e-02 -6.233034e-02
## [665,] 2.098111e-02 5.939266e-03
## [666,] 1.062091e-02 1.735479e-02
## [667,] 2.460832e-02 2.131430e-02
## [668,] 3.944586e-02 9.116456e-03
## [669,] 2.516373e-02 2.762077e-02
## [670,] 4.977205e-02 4.893507e-02
## [671,] 5.214071e-04 5.754061e-02
## [672,] 1.646978e-02 5.795572e-02
## [673,] 2.549245e-02 -4.607719e-02
## [674,] 3.071789e-03 9.068554e-03
## [675,] -1.117632e-02 -2.366126e-02
## [676,] 7.254409e-03 3.075006e-02
## [677,] 5.588162e-03 1.183063e-02
## [678,] -2.486902e-02 -5.008460e-02
## [679,] 1.369270e-02 2.642334e-02
## [680,] -1.372671e-02 2.481079e-02
## [681,] -9.770782e-03 -3.351214e-02
## [682,] -2.741940e-02 -1.612546e-03
## [683,] -3.411840e-02 -2.605613e-02
## [684,] -6.212723e-02 1.725898e-02
## [685,] -6.575443e-02 1.883952e-03
## [686,] -9.487409e-02 3.258611e-02
## [687,] -4.454662e-02 8.782766e-02
## [688,] -2.157054e-02 3.898839e-02
## [689,] 5.293450e-03 3.429446e-02
## [690,] -1.653779e-02 4.451255e-02
## [691,] -7.843832e-03 1.417760e-02
## [692,] -3.278087e-02 6.656126e-02
## [693,] -4.395720e-02 4.290000e-02
## [694,] -4.006929e-02 8.704534e-02
## [695,] -4.702899e-02 3.383145e-02
## [696,] -5.150632e-02 3.461377e-02
## [697,] -3.862974e-02 2.596033e-02
## [698,] -6.912093e-02 1.527923e-02
## [699,] -1.317129e-02 3.111727e-02
## [700,] -2.826953e-02 1.454481e-02
## [701,] -5.214071e-04 -5.754061e-02
## [702,] -5.032747e-03 -5.524155e-03
## [703,] -9.510078e-03 -4.741833e-03
## [704,] 3.661213e-03 -3.585910e-02
## [705,] -1.424811e-02 -3.272982e-02
## [706,] 4.255165e-02 -3.304912e-02
## [707,] 1.960958e-02 -3.544399e-02
## [708,] 2.630857e-02 -1.100041e-02
## [709,] 3.140934e-02 -1.079445e-01
## [710,] 5.402269e-02 -3.185169e-02
## [711,] 2.434762e-02 -7.456008e-03
## [712,] -4.088541e-02 5.196855e-02
## [713,] -3.607936e-02 -2.251173e-02
## [714,] -1.065492e-02 3.387935e-02
## [715,] -6.186652e-02 4.602929e-02
## [716,] -6.493831e-02 3.696073e-02
## [717,] -7.977585e-02 4.915857e-02
## [718,] -6.212723e-02 1.725898e-02
## [719,] -4.536274e-02 5.275088e-02
## [720,] -1.069005e-01 2.508220e-02
## [721,] -8.229222e-02 4.639650e-02
## [722,] -9.964613e-02 5.583226e-02
## [723,] -9.039676e-02 3.180379e-02
## [724,] -6.578844e-02 5.311809e-02
## [725,] -8.988671e-03 5.279878e-02
## [726,] -5.066755e-03 4.570998e-02
## [727,] -3.248616e-02 4.409743e-02
## [728,] -1.261588e-02 3.742375e-02
## [729,] -3.167004e-02 7.917422e-02
## [730,] -2.575316e-02 1.730688e-02
## [731,] -3.696349e-02 4.487976e-02
## [732,] -9.261842e-02 6.577883e-03
## [733,] -3.862974e-02 2.596033e-02
## [734,] -1.402142e-02 4.727462e-02
## [735,] -3.156802e-02 -7.452819e-02
## [736,] 1.173174e-02 2.996774e-02
## [737,] -3.071789e-03 -9.068554e-03
## [738,] 1.539295e-02 -5.891365e-03
## [739,] -3.332492e-03 -3.783886e-02
## [740,] 9.215367e-03 2.720566e-02
## [741,] -5.882874e-03 1.063320e-02
## [742,] -1.176575e-02 2.126639e-02
## [743,] -5.622171e-03 3.940350e-02
## [744,] 2.127583e-02 -1.652456e-02
## [745,] -2.482365e-03 -5.399621e-02
## [746,] 7.549120e-03 8.286232e-03
## [747,] -5.261715e-02 2.200081e-02
## [748,] -8.954663e-03 1.564644e-03
## [749,] -2.882495e-02 8.238330e-03
## [750,] -4.810581e-02 -3.001564e-02
## [751,] -2.797482e-02 -7.919022e-03
## [752,] -9.510078e-03 -4.741833e-03
## [753,] -1.509824e-02 -1.657246e-02
## [754,] -5.554154e-04 -6.306476e-03
## [755,] -2.964107e-02 -2.683845e-02
## [756,] -2.947118e-04 2.246383e-02
## [757,] -5.343327e-02 -1.307597e-02
## [758,] -6.404281e-03 -4.690741e-02
## [759,] 3.526324e-02 -1.256505e-02
## [760,] -3.049119e-02 -1.068110e-02
## [761,] -1.653779e-02 4.451255e-02
## [762,] -8.104536e-03 -1.459271e-02
## [763,] 6.143578e-03 1.813711e-02
## [764,] 2.908565e-02 2.053197e-02
## [765,] 2.797482e-02 7.919022e-03
## [766,] 2.405290e-02 1.500782e-02
## [767,] 3.441311e-02 3.592301e-03
## [768,] 1.820404e-02 -2.559312e-02
## [769,] 2.438162e-02 -5.869014e-02
## [770,] -2.845094e-03 7.093588e-02
## [771,] 2.630857e-02 -1.100041e-02
## [772,] 5.293450e-03 3.429446e-02
## [773,] 3.186273e-02 5.206436e-02
## [774,] 1.931487e-02 -1.298016e-02
## [775,] 4.503402e-02 2.094709e-02
## [776,] 5.428339e-02 -3.081385e-03
## [777,] 6.127710e-02 -1.101630e-03
## [778,] 7.241942e-02 7.379377e-02
## [779,] 1.846474e-02 3.177189e-03
## [780,] 1.542696e-02 -5.712550e-02
Plot Fig. 4.3
U <- svdptd$u #The matrix of the EOF spatial patterns
V <- svdptd$v #The matrix of the PC temporal patterns
D <- diag(svdptd$d) #The diagonal matrix of standard deviations
recontada <- U%*%D%*%t(V) #Verify the SVD decomposition
svd_error <- ptada - recontada
max(svd_error)## [1] 3.996803e-15min(svd_error)## [1] -4.218847e-15
Plot WSOI1
#setEPS()
#postscript("fig0403a.eps", height = 4.9, width = 7)
par(mar = c(4,4.2,2.0,0.5))
xtime <- seq(1951, 2016-1/12, 1/12)
wsoi1 <- D[1,1]*t(V)[1,]
plot(xtime, wsoi1,type = "l", cex.lab = 1.2, cex.axis = 1.2,
xlab = "Year",ylab = "WSOI1",
col = "black",xlim = range(xtime), ylim = c(-4,4))
axis (3, at = seq(1951,2015,4), labels = seq(1951,2015,4))
text(1985,-3.7,"Weighted SOI1 Index", cex = 1.2)
text(1952,3.5, "(a)", cex = 1.2)
# dev.off()
Plot WSOI2
#setEPS()
#postscript("fig0403b.eps", height = 4.9, width = 7)
par(mar = c(4,4.2,2.0,0.5))
xtime <- seq(1951, 2016-1/12, 1/12)
wsoi2 <- D[2,2]*t(V)[2,]
plot(xtime, wsoi2,type = "l", cex.lab = 1.3, cex.axis = 1.3,
xlab = "Year",ylab = "WSOI2",
col = "blue",xlim = range(xtime), ylim = c(-4,4))
axis (3, at = seq(1951,2015,4), labels = seq(1951,2015,4))
text(1985,-3.7,"Weighted SOI2 Index", cex = 1.5)
text(1952,3.5, "(b)", cex = 1.3)
# dev.off()
Cumulative WSOIs
setwd("~/sshen/climmath/")
#setEPS()
#postscript("fig0403c.eps", height = 4.9, width = 7)
par(mar = c(4,4.2,2.0,4.2))
cwsoi1 <- cumsum(wsoi1)
plot(xtime, cwsoi1,type = "l",cex.lab = 1.2, cex.axis = 1.2,
xlab = "Year",ylab = "Weighted SOI1",
col = "black",lwd = 3, ylim = c(-200,50),
main = "Comparison between CWSOI1 and the Smoothed AMO Index")
text(1951,40, "(c)", cex = 1.3)
#axis (3, at = seq(1951,2015,4), labels = seq(1951,2015,4))
legend(1970,20, col = c("black"),lty = 1,lwd = 3.0,
legend = c("CWSOI1"),bty = "n",text.font = 2,cex = 1.2)
#Superimpose CSOI time series on this CWSOI1
cnegsoi <- -cumsum(ptamonv-pdamonv)
lines(xtime, cnegsoi,type = "l",col = "purple", lwd = 3.0)
legend(1970,2, col = c("purple"),lty = 1,lwd = 3.0,
legend = c("CSOI"),bty = "n",text.font = 2,cex = 1.2)
#24-month ahead moving average of the monthly AMO index
amodat <- read.table("data/AMO1951-2015.txt", header = FALSE)
amots <- as.vector(t(amodat[,2:13]))
#install.packages("TTR")
library("TTR")## Warning: package 'TTR' was built under R version 4.1.3amomv <- SMA(amots,n = 24, fill = NA)
#Average of the previous n points
par(new = TRUE)
xtime <- seq(1951,2015,len = 780)
plot(xtime, amomv[37:816],type = "l",col = "darkgreen",
lwd = 3,axes = FALSE,xlab = "",ylab = "",
cex.lab = 1.2, cex.axis = 1.2,
ylim = c(-0.5, 0.5))
legend(1970,0.23, col = c("darkgreen"),lty = 1,lwd = 3.0,
legend = c("AMO index"),bty = "n",text.font = 2,cex = 1.2)
#Suppress the axes and assign the y-axis to side 4
axis(4, cex.axis = 1.2, col.axis = "darkgreen")
mtext("Smoothed AMO index",
side = 4,line = 3, cex = 1.2, col = "darkgreen")
# dev.off()
Plot cumulative WSOI2: CWSOI2
#setEPS()
#postscript("fig0403d.eps", height = 4.9, width = 7)
par(mar = c(4,4.2,2.0,0.5))
cwsoi1 <- cumsum(wsoi1)
wsoi2 <- D[2,2]*t(V)[2,]
cwsoi2 <- cumsum(wsoi2)
plot.new()
# par(mar = c(4,4,3,4))
plot(xtime, cwsoi2,type = "l",xlab = "Year",ylab = "CWSOI2",
cex.lab = 1.2, cex.axis = 1.2,
col = "blue",lwd = 3, ylim = c(-200,50),
main = "CWSOI2 Index: The Cumulative PC2")
text(1951,40, "(d)", cex = 1.3)
# dev.off()
Plot Fig. 4.4
#Display the two ENSO modes on a world map
library(maps)## Warning: package 'maps' was built under R version 4.1.3library(mapdata)## Warning: package 'mapdata' was built under R version 4.1.3plot.new()
par(mfrow = c(2,1))
par(mar = c(0,0,0,0)) #Zero space between (a) and (b)
map(database = "world2Hires",ylim = c(-70,70), mar = c(0,0,0,0))
grid(nx = 12,ny = 6)
points(231, -18,pch = 16,cex = 2, col = "red")
text(231, -30, "Tahiti 0.61", col = "red")
points(131, -12,pch = 16,cex = 2.6, col = "blue")
text(131, -24, "Darwin -0.79", col = "blue")
axis(2, at = seq(-70,70,20),
col.axis = "black", tck = -0.05, las = 2, line = -0.9,lwd = 0)
axis(1, at = seq(0,360,60),
col.axis = "black",tck = -0.05, las = 1, line = -0.9,lwd = 0)
text(180,30, "El Nino Southern Oscillation Mode 1",col = "purple",cex = 1.3)
text(10,-60,"(a)", cex = 1.4)
box()
# par(mar = c(0,0,0,0)) #Plot mode 2
map(database = "world2Hires", ylim = c(-70,70), mar = c(0,0,0,0))
grid(nx = 12,ny = 6)
points(231, -18,pch = 16,cex = 2.6, col = "red")
text(231, -30, "Tahiti 0.79", col = "red")
points(131, -12,pch = 16,cex = 2, col = "red")
text(131, -24, "Darwin 0.61", col = "red")
text(180,30, "El Nino Southern Oscillation Mode 2",col = "purple",cex = 1.3)
axis(2, at = seq(-70,70,20),
col.axis = "black", tck = -0.05, las = 2, line = -0.9,lwd = 0)
axis(1, at = seq(0,360,60),
col.axis = "black",tck = -0.05, las = 1, line = -0.9,lwd = 0)
text(10,-60,"(b)", cex = 1.4)
box()
An example of multivariate regression
x1 <- c(1,2,3) #Given the coordinates of the 3 points
x2 <- c(2,1,3)
y <- c(-1,2,1)
df <- data.frame(x1,x2,y) #Put data into the data.frame format
fit <- lm(y ~ x1 + x2, data = df)
fit #Show the regression results##
## Call:
## lm(formula = y ~ x1 + x2, data = df)
##
## Coefficients:
## (Intercept) x1 x2
## -5.128e-16 1.667e+00 -1.333e+001.667*x1-1.333*x2 #Verify that 3 points determining a plane## [1] -0.999 2.001 1.002
Multivariate regression with four data points
u <- c(1,2,3,1)
v <- c(2,4,3,-1)
w <- c(1,-2,3,4)
mydata <- data.frame(u,v,w)
myfit <- lm(w ~ u + v, data = mydata)
summary(myfit) #Show the result##
## Call:
## lm(formula = w ~ u + v, data = mydata)
##
## Residuals:
## 1 2 3 4
## 1.0 -1.0 0.5 -0.5
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.0000 1.8708 0.535 0.687
## u 2.0000 1.2472 1.604 0.355
## v -1.5000 0.5528 -2.714 0.225
##
## Residual standard error: 1.581 on 1 degrees of freedom
## Multiple R-squared: 0.881, Adjusted R-squared: 0.6429
## F-statistic: 3.7 on 2 and 1 DF, p-value: 0.345
Another example of multivariate regregression
dat <- matrix(rnorm(40),nrow = 10, dimnames = list(c(letters[1:10]), c(LETTERS[23:26])))
fdat <- data.frame(dat)
fit <- lm(Z ~ W + X + Y, data = fdat)
summary(fit)##
## Call:
## lm(formula = Z ~ W + X + Y, data = fdat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.68252 -0.30759 -0.06481 0.22976 1.03628
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.5234 0.2147 2.438 0.0506 .
## W -0.2324 0.1949 -1.192 0.2782
## X -0.3785 0.1705 -2.220 0.0682 .
## Y 0.3289 0.2896 1.136 0.2995
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6385 on 6 degrees of freedom
## Multiple R-squared: 0.4694, Adjusted R-squared: 0.2042
## F-statistic: 1.77 on 3 and 6 DF, p-value: 0.2526