2d_regression.R

# Y = data matrix. 
#     yearly u anomalies for the northern hemisphere between 1948 and 2010. 
#     63 rows and 5328 columns (2.5 degree grid)
# 
# X = timeseries matrix. 
#     column vector with values 1948 to 2010
#     63 rows 1 column
N <- nrow(X)

# Singular value decomposition
eof <- svd(Y)
U <- eof$u
S <- diag(eof$d, nrow = length(eof$d))
V <- eof$v

# Normalization
F <- sqrt(N)*(U)
E <- V%*%S/sqrt(N)

# Eq 24:
B_tilde <- c(solve(t(X) %*% X) %*% t(X) %*% F)

# Eq 25:
B <- B_tilde %*% t(E)

# Linear regression point by point
regr <- data[, FitLm(var, year), by = .(lon, lat)][term == "year"]

# Merge everything
X_data$coldims[, B := c(B)]
regr <- X_data$coldims[regr, on = c("lon", "lat")]

# Comparing poinwise regression with eof regression
with(regr, plot(estimate, B, xlab = "Pointwise regression", ylab = "EOF regression"))  

with(regr, mean(estimate/B))
# > 11845.08