# RPD -- Examples 3.1+ -- X = Ozone, y = Yield

X <- cbind(1, c(0.02, 0.07, 0.11, 0.15))

y <- rbind(242, 237, 231, 201)

XtX <- t(X) %*% X

Xty <- t(X) %*% y

XtXinv <- solve(XtX)

betahat <- XtXinv %*% Xty

P <- X %*% XtXinv %*% t(X)

yhat <- P %*% y

e <- (diag(1,4) - P) %*% y

print(X)
print(y)
print(XtX)
print(Xty)
print(XtXinv)
print(betahat)
print(P)
print(yhat)
print(e)


# Much more computationally efficient ways to obtain betahat, yhat, and e:

betahat <- solve(XtX,Xty)  # X = solve(A,B) solves A %*% X = B

yhat <- X %*% betahat

e <- y - yhat

print(betahat)
print(yhat)
print(e)