#' Principal component analysis
#'
#' Calculates the principal components accorindg to the algorithm by Wold
#'
#' @param X matrix or numeric dataframe.
#' @param a integer, number of principal components to calculate.
#' @param it integer, maximum number of iterations if error tolerance is not
#' reached.
#' @param tolerance float, precision tolerance.
#'
#' @return list, T: PCA scores. P: PCA loadings. pcvar: proportion of varation
#' explained by PCA components
#'
#' @export
#' @examples
#' data(iris)
#' nipals(iris[, 1:4], a = 2)
nipals <- function(X, a, it = 10, tolerance = 1e-4){
if (a > ncol(X)) {
stop("cannot calculate more components than variables", call. = FALSE)
}
Xh <- scale(X, center = TRUE, scale = FALSE) #mean-centering of data matrix X
nr <- 0
T <- NULL
P <- NULL
pcvar <- NULL
varTotal <- sum(diag(var(Xh)))
currVar <- varTotal
precision <- c()
for (h in 1:a) {
th <- Xh[,1] # starting value for th is 1st column of Xh
ende <- FALSE
# 3 inner steps of NIPALS algorithm
while (!ende){
nr <- nr + 1
# the result of matrix multiplication operation (%*%) is a matrix of a single
# valule. A matrix cannot multiply another using scalar multiplication (*).
# as.vector convert a value of class matrix to a value of class double.
# (A'*B)' = B'*A
ph <- t((t(th) %*% Xh) * as.vector(1/(t(th) %*% th))) #LS regression for ph
ph <- ph * as.vector(1/sqrt(t(ph) %*% ph)) #normalization of ph
thnew <- t(t(ph) %*% t(Xh) * as.vector(1 / (t(ph) %*% ph))) #LS regression for th
prec <- t(th - thnew) %*% (th - thnew) # calculate precision
precision <- append(precision, prec) # record precision
th <- thnew #refresh th in any case
#check convergence of th
if (prec <= (tolerance ^ 2)) {
ende <- TRUE
} else if (it <= nr) { #too many iterations
ende <- TRUE
cat("\nWARNING! Iteration stop in h=",h," without convergence!\n\n")
}
}
Xh <- Xh - (th %*% t(ph)) #calculate new Xh
T <- cbind(T,th) #build matrix T
P <- cbind(P,ph) #build matrix P
oldVar <- currVar
currVar <- sum(diag(var(Xh)))
pcvar <- c(pcvar,(oldVar - currVar) / varTotal)
nr <- 0
}
list(T = T, P = P, pcvar = pcvar, precision = precision)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.