##-----------------------------------------------------------------##
##            Script for Writing R Packages                        ##
##                        John Fox                                 ##
##                     IQS Barcelona                               ##
##                     January  2015                               ##
##-----------------------------------------------------------------##

  # Functions for a simple matrix-algebra package

source("http://socserv.socsci.mcmaster.ca/jfox/Courses/R/IQSBarcelona/matrixDemos.R")

objects()

# some examples

A <- matrix(c(1,2,3,4,5,6,7,8,10), 3, 3) # a nonsingular matrix
A
RREF(A) # the reduced row-echelon form of A

b <- 1:3
RREF(A, b)  # solving the matrix equation Ax = b

RREF(A, diag(3)) # inverting A

B <- matrix(1:9, 3, 3) # a singular matrix
B
RREF(B)

RREF(B, b)

RREF(B, diag(3))

library(MASS) # for fractions

Ginv(A, fractions=TRUE)  # a generalized inverse of A = inverse of A
round(Ginv(A) %*% A, 6)  # check

Ginv(B, fractions=TRUE)  # a generalized inverse of B
B %*% Ginv(B) %*% B   # check

C <- matrix(c(1,2,3,2,5,6,3,6,10), 3, 3) # nonsingular, symmetric
C
cholesky(C)
cholesky(C) %*% t(cholesky(C))  # check

EC <- eig(C) # eigenanalysis of C
EC$vectors %*% diag(EC$values) %*% t(EC$vectors) # check

(D <- A[, -3])

(SD <- SVD(D))  # singular value decomposition of D
SD$U %*% diag(SD$d) %*% t(SD$V) # check

(QRD <- QR(D))  # QR decomposition of D
QRD$Q %*% QRD$R  # check

      # LS by the SVD

data("Duncan", package="car")
X <- cbind(1, as.matrix(Duncan[, c("income", "education")]))
head(X)
y <- Duncan$prestige
head(y)
(svd <- SVD(X))
(b <- svd$V %*% diag(1/svd$d) %*%t(svd$U) %*% y)
b - coef(lm(prestige ~ income + education, data=Duncan))

    # principal-components analysis

data("DavisThin", package="car")
(R <- cor(DavisThin))
eig(R)

    # Create the package "skeleton"

objects()
remove(A, b, B, C, D, EC, SD, svd, QRD, X, y)
objects()

setwd("c:/temp") # make a mess in the temp directory

package.skeleton("matrixDemos")

# Check and install package using the tools in the RStudio Build tab
#    or:
# to check package (at OS prompt): R CMD check matrixDemos
# to install package (at OS prompt): R CMD INSTALL matrixDemos
# Then, in a fresh session:

library(matrixDemos)
help(package="matrixDemos")
?MatrixDecompositions
example("MatrixDecompositions")