Introduction to Linear Algebra in R	STP 526, Fall 2009, S. Lohr

The statistical software package R may be downloaded from www.r-project.org; 
it has many functions that will perform linear algebra calculations.
Try the following commands in R to see how they work. (Don't just cut and
paste; type them in yourself so you start to learn them!)
Note that '#' is the comment character in R; it ignores everything
behind the #'s.
The symbol '<-' is the assignment operator.

v1 <- c(3,4,5)  	#Assigns the vector (3,4,5)' to v1
v1		#Prints the contents of v1 (Note:  this is a column
		#vector, but R displays the elements in a row.)
v2 <- c(1,2,-1)	#Assigns the vector (1,2,-1)' to v2
v1[3]		#Selects third element of v1

mat1 <- cbind(v1,v2)   #Matrix mat1 is [v1 v2] (3 x 2 matrix)
mat2 <- rbind(v1,v2)   #Matrix mat2 is mat1'  (2 x 3 matrix)
mat3 <- t(mat1)	#Transpose of mat1
mat2[1,2]	#Selects (1,2) element of mat2
mat3[,1:2]	#Selects first two columns of mat3

4*mat1		#Scalar multiplication
v1+v2		#Vector (or matrix) addition
mat3 <- mat1 %*% mat2	#Matrix multiplication
mat4 <- mat2 %*% mat1
t(v1) %*% v2	#Calculates v1' v2
crossprod(v1,v2)#Same as above, but faster computationally
dim(mat1)	#Gives number of rows and columns in mat1
diag(mat3)	#Extracts diagonal elements of mat3, puts in vector
id5 <- diag(5)    #Creates 5 x 5 identity matrix
diag		#Prints the function "diag"
v3 <- rep(1,10)	#Creates 10 x 1 vector of ones
v3 %*% t(v3)	#Creates 10 x 10 matrix of ones (= J)

det(mat3)	#Calculates the determinant of mat3
sum(diag(mat3)) #Calculates trace of a matrix
kronecker(v1,mat1)    #Finds Kronecker product of v1 and mat1

solve(mat4)	#Finds inverse of mat4
solve(mat4,c(1,2))      #Finds solution x to mat4 %*% x = (1,2)'
help("solve")	#Gives the help information for function "solve"
help.start()	#Use when you don't know the name of the command you want

# Note that solve(mat4) is in decimal form. Sometimes you may want to
# express terms as fractions. To do this, load the package "MASS"
# from the packages menu in R.  Then try

fractions(solve(mat4))

# Note: "fractions" is not always exact, so check your work when using it.

apply(mat2,2,sum)  #Sums each column of mat2 (column is dimension 2)
apply(mat2,1,sum)  #Sums each row of mat2 (row is first dimension)
		#Note:  can substitute any other function of vectors
		#for sum---try it with mean or var instead of sum.

scale(mat3)	#Centers and scales the columns of mat3, so each
		#column in the result has mean 0 and std dev 1

eigen(mat3)	#Finds eigenvalues and eigenvectors of mat3
svd(mat3)	#Creates singular value decomposition of mat3;
		#mat3 = u %*% diag(d) %*% t(v)

mat5 <- chol(mat4)#Returns Cholesky decomposition of mat4, if mat 4
	#is symmetric and non-negative definite.
	#Gives the upper triangular matrix: t(mat5)%*%mat5=mat4

qr(mat3)	#Returns qr decomposition of mat3; mat3 = Q%*%R, with
	#Q orthogonal and R upper triangular; also gives rank of matrix.
	#Does not give Q and R explicitly unless you give options
qr(mat3)$rank # Gives rank of matrix
qr.solve(mat3[,1:2],mat3[,3])  #Solve a system of equations

var(mat1)	#Finds variance matrix for columns of mat1.
cor(mat1)	#Finds correlation matrix for columns of mat1.

y <- c(5,4,5)
lm(y ~ v1+v2)	#Fits linear regression model:  y =b0+ b1 v1 + b2 v2	
temp <- lm(y~v1+v2)
summary(temp) 	# to see statistics
        #Note this is a saturated model so all SE's = 0

# Generating random variates; doing simple plots and statistics

u1 <- rnorm(200,5,3) # generates 100 normals with mean 5, s.d. 3
hist(u1) 	# draws a histogram
u2 <- runif(200,0,2) # generates 50 Unif(0,2) variates
hist(u2)
u3 <- rchisq(200,4) # generates 50 chi-square(4) variates
hist(u3)
plot(u2,u3)	#Note the x-axis variable is first:  plot(x,y)
boxplot(u1,u2,u3) # side by side boxplots
mean(u1)
median(u1)
quantile(u1)
var(u1)
cor(u1,u2)
cov(u1,u2)	#Covariance of u1 and u2
cov(cbind(u1,u2)) #Covariance matrix of u1 and u2
summary(u1)




# Writing R functions
# You'll often want to write your own function to calculate something not 
# already in R, or to do several steps in succession. 

# Here's an example of a function to permute two rows of a matrix, indexed
# by the quantities row1 and row2.
# The last line, amat, indicates that the value amat is returned at
# the end of the function

permuterow <- function(amat, row1, row2) {
    temp <- amat[row1,]
    amat[row1,] <- amat[row2,]
    amat[row2,] <- temp
    amat
}

permuterow(mat1,2,3)