## R-code (September 13, 2007)


## make 4 plots for pairs of data with various correlations

cor.plot <- function(rho, n){
library(MASS)

par(mfrow=c(2,2))

Sigma <- matrix(c(1,rho,rho,1),2,2)
Mu <- matrix(c(0,0),2,1)

for(i in 1:4){
yx <- mvrnorm(n, Mu, Sigma)
plot(yx[,2], yx[,1])
ell <- ellipse(Sigma, level=0.95)
lines(ell[,1],ell[,2],lwd=2)
}

}


cor.plot(0, 1000)


## quadratic example
x <- seq(-5,5, by=0.1)
y <-  -x^2 + rnorm(length(x))
plot(x, y)
cor(x,y)
lm.a <- lm(y~x)
summary(lm.a)
abline(lm.a)

lm.b <- lm(y ~ x + I(x^2))
summary(lm.b)


## read in the mathmarks data again!
X <- read.table("mathmarks.txt")

## test whether X is a matrix
is.matrix(X)

## test whether X is a data frame
is.data.frame(X)

## set X to be a matrix
X <- as.matrix(X)
is.matrix(X)

n <- dim(X)[1]
p <- dim(X)[2]

x.bar <- apply(X, 2, mean)  ## take the mean for each column of X.
S <- cov(X)
R <- cor(X)

## calculate the sample mean by hand using R
one <- matrix(rep(1, n), nrow=n, ncol=1)
x.bar.2 <- (1/n)*t(X)%*%one

## calculate the sample covariance matrix by hand using R
Ident <-diag(1,n)
S.2 <- (1/(n-1))*t(X)%*%(Ident-(1/n)*one%*%t(one))%*%X

## calculate the sample correlation matrix by hand using R
D <- diag( diag(S.2), p)
D.sqrt <- sqrt(D)
R.2 <- solve(D.sqrt)%*%S.2%*%solve(D.sqrt)