> pdf("muscle1.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)))
      [,1]     [,2]
 [1,]  177 177.4163
 [2,]  279 269.5436
 [3,]  346 312.9811
 [4,]  160 173.8721
 [5,]  193 196.0639
 [6,]  280 291.0358
 [7,]  335 342.5129
 [8,]  169 174.6943
 [9,]  212 226.0206
[10,]  244 259.4179
[11,]  285 292.6506
[12,]  181 179.0599
[13,]  298 297.4362
[14,]  212 205.2187
[15,]  317 299.2128
[16,]  347 352.2495
[17,]  186 184.3223
[18,]  216 207.3476
[19,]  265 269.5157
[20,]  306 302.4635
[21,]  348 354.9175
[22,]  209 188.7710
[23,]  324 308.7242
[24,]  352 364.7397
> 
> PRESS <- sum((Y-Yhat_i)^2)
> MSEP <- PRESS/length(Y)
> 
> print(c(PRESS,MSEP))
[1] 3442.4490  143.4354
> 
> dev.off()