###############  Handout #1 for ST434/534  ###############

###introduction to R####################

##R package can be downloaded from website for free http://www.r-project.org##

##use "help" to know the functions##

help(sqrt) #r will tell you the usage of sqrt function and give you exmaples how to use the function

## or use help.search("***")

###### Matrix Algebra #########


## Using a '#' comments out any text afterwards on a line

######### Vectors ###########

## Concatenate function 'c' to create a vector
u <- c(1,2,3,4)

## print u
print(u)

##just typing u also prints
u

u <- 1:4  ## makes the same vector as above
u
v <- seq(0,1.5,by=.5)
v
v <- seq(0,1.5,length=4) ## same thing
w <- rep(u, times=2)
w

## Scalar mult
2*u

## Vector Addition
u+v

## square every element of u
u^2

## multiply each element of u by the corresponding elements of v
u*v

## sum up every element of u
sum(u)

## subsetting
v[1]        ## get the first element of v
v[1:3]      ## get the first 3 elements of v
v[c(2,4)]   ## get the 2nd and 4th elements of v
v[-4]       ## get all but the 4th element of v

## Random numbers
x <- rnorm(n=10, mean=10, sd=1)
y <- rnorm(n=10)                ## default generates N(0,1)

## plotting.  
plot(x,y)



######### Matrices ##########

A <- rbind(u,v)
A
B <- cbind(u,v)
B

## subsetting
A[1,1]         ## get the (1,1) element of A
A[1:2,1:3]     ## get the first 3 rows of A with the first three columns
A[,1:2]        ## get the 1st and 2nd columns of A (with all rows)
B[-3,]         ## get all but the 3rd row of B (with all columns)
               ## (e.g. eliminate he 3rd observation)

## transpose
t(A)

## Matrix multiplication
C <- A%*%B
C
A%*%C  ## doesn't work b.c. A is (2x4) and C is (2x2)
C%*%A  ## does work

## Element by element multiplication
C*C   ## is not the same as C%*%C
C%*%C
A*C   ## doesn't work, A and C must have the same dimensions

## Matrix Inverse
C.inv <- solve(C)
C.inv
C%*%C.inv

## Note: C^(-1) does NOT work
C%*%C^(-1)    ## does not equal the identity

## Eigenvalues/vectors
decomp <- eigen(C)
decomp

## Spectral Representation
Lambda <- diag(decomp$values)
Lambda
P <- decomp$vectors
C.spec <- P%*%Lambda%*%t(P)

## Square Root of C
C.sqrt <- P%*%sqrt(Lambda)%*%t(P)
C.sqrt%*%C.sqrt

## Note: Not the same as element by element square root
sqrt(C)
sqrt(C)%*%sqrt(C)   ## does not equal C

## Generalized inverse
D <- B%*%A
D
det(D)  ## no inverse exists; columns are linearly dependent.

## Create Ginv for D
decomp <- eigen(D)
P <- decomp$vectors
lambda <- decomp$values
lambda
lambda.recip <- rep(0, length(lambda))
lambda.recip
logic <- (abs(lambda) > 1E-10)
logic
lambda[logic]
lambda.recip[logic] <- lambda[logic]^(-1)
lambda.recip
Lam.ginv <- diag(lambda.recip)
D.ginv <- P%*%Lam.ginv%*%t(P)
D%*%D.ginv%*%D                ## this equals D
D


####### Data Frames #########

## Data frames are basically the same as matrices

## Turn a matrix into a data frame
B.data <- as.data.frame(B)

## read data into a data frame from a webpage
ex.data <- read.table(file="http://www.stat.unm.edu/**/**.txt", header=T)
ex.data

## Read data from a local file
ex.data <- read.table(file="W:/teaching/stat445545anova/CH16PR07.txt")
ex.data


######## Functions ##########

## You won't use these much in this class but 
## they are very useful in general.

## Creating a function that adds two vectors together

my.add <- function(a,b){
  return(a + b)
}

my.add(u,v)


my.add2 <- function(a,b){  ## performs the same function as above
  na <- length(a)
  c <- rep(0, na)
  for(i in 1:na){
    c[i] <- a[i] + b[i]
  }
  return(c)
}

my.add2(u,v)