# Edps/Soc 584 and Psych 594 example on how to 
#
#  Read in data set into a matiex
#  Write functions to compute
#   (i) n, xbar
#   (ii) S
#   (iii) R
#   (iv) n, xbar, S, and R
#

#
#
# First get a sample data set
#
scores <- read.table('scores_data.txt', header=TRUE)
attach(scores)
names(scores)                           # look at the variable names
head(scores)                            # see the first 6 lines of data
tail(scores)                            # see the last 6 lines of data

X <- as.matrix(scores[,2:5])            # put the variables for the test scores 
                                        # into matrix X

p <- ncol(X)                            # Check to make sure things worked OK
p
n <- nrow(X)
head(X)
tail(X)

#
# Define function to compute mean vector that prints out the
#  n = sample size 
#  xbar = column vector with means in it
#
meanvector <- function(X){
    n <- nrow(X)
    one <- array(1, dim=c(n))
    xbar <- t(X) %*% one/n
    list(n=n, xbar=xbar)
    }

#
# Examples of running the function:
#
# (a) Print n and xbar to console

meanvector(X)

# (b) Print only xbar to console
meanvector(X)$xbar

# (c) Print only n to console
meanvector(X)$n

# (d) Set xbartest equal to the mean vector 
xbartest <-meanvector(X)$xbar


#
# Function that computes  S = sample covariance matrix
#
Smatrix  <- function(X) {
    n <- nrow(X)
    one <- array(1, dim=c(n))
    xbar <- t(X) %*% one/n
    W <- t(X - one %*% t(xbar)) %*% (X - one %*% t(xbar))
    S <- W / (n-1) 
    list(S=S)
     }

#
# Function that computes R = sample covariance matrix
#
#  Note that this one calls the function Smatrix 
#
Rmatrix <- function(X){
    S <- Smatrix(X)$S
    Dsqrt <- diag(diag(sqrt(S)) )
    invD  <- solve(Dsqrt)
    R <- invD %*% S %*% invD
    list (R=R)
    }    

#
# Function that computes 
#  n   = sample size
#  xbar= column vector of means
#  S   = sample covariance matrix
#  R   = sampe correlation matirx
#

samplestats <- function(X) {
    n <- nrow(X)
    one <- array(1, dim=c(n))
    xbar <- t(X) %*% one/n
    W <- t(X - one %*% t(xbar)) %*% (X - one %*% t(xbar))
    S <- W / (n-1) 
    Dsqrt <- diag(diag(sqrt(S)) )
    invD  <- solve(Dsqrt)
    R <- invD %*% S %*% invD
    list (n = n, xbar=xbar, S = S, R=R )
    }

#
# Read male and female data into different matrices
#  and compute mean vectors and the difference between them
#
male   <- as.matrix(scores[1:32,2:5])
female <- as.matrix(scores[33:64,2:5])

maleall <- samplestats(male)
mbar <- maleall$xbar


femaleall <- samplestats(female)
fbar <- femaleall$xbar

difference <- mbar - fbar