pdf("muscle4.pdf")


XY1 <- matrix(c(43.7,19,177,    43.7,43,279,    43.7,56,346,    54.6,13,160,
54.6,19,193,    54.6,43,280,    54.6,56,335,    55.7,13,169,    55.7,26,212,
55.7,34.5,244,  55.7,43,285,    58.8,13,181,    58.8,43,298,    60.5,19,212,
60.5,43,317,    60.5,56,347,    61.9,13,186,    61.9,19,216,    61.9,34.5,265,
61.9,43,306,    61.9,56,348,    66.7,13,209,    66.7,43,324,    66.7,56,352),
byrow=T,ncol=3)

# XY1 is a matrix containing n=24 rows and 3 columns (X1,X2,Y) we must form the X and Y matrices from it
# This is similar to when we input a raw data frame (dataset) and create a matrix

X0 <- matrix(rep(1,nrow(XY1)),ncol=1)   #  Create a Column of 1s for the intercept (same number of rows as XY1)

X <- cbind(X0,XY1[,1],XY1[,2])    # Bind the columns for Intercept (X0), X1 (XY[,1]), X2 (XY[,2]) into X

Y <- XY1[,3]

Yhat_i <- numeric(length(Y))                    

for (i in 1:length(Y)) {

X_i <- X[-i,]
Y_i <- Y[-i]

betahat_i <- solve(t(X_i) %*% X_i) %*% t(X_i) %*% Y_i
Yhat_i[i] <- X[i,] %*% betahat_i

}

print(cbind(as.matrix(Y),as.matrix(Yhat_i)))

PRESS <- sum((Y-Yhat_i)^2)
MSEP <- PRESS/length(Y)

print(c(PRESS,MSEP))

dev.off()