#Problem #5c
#Residuals

plot(X,Y, xlab="Year", ylab="log(Annual Rainfall)")
abline(lsfit(X,Y), lwd=2)

N=length(Y)

for(i in 1:5){		#Perform 5 itereations


Yhat = XX %*% Beta
epsilon = Y - Yhat 

#Weight the Residuals
C = 3
S = diff(quantile(epsilon, c(0.25, 0.75)))
S = 4*sd(epsilon)/3

U = epsilon/(C*S)
U[abs(U) > 1]=1
w=(1-(U*U))^2

Wt=matrix(0,ncol=N,nrow=N)
diag(Wt)=w
Beta = solve(t(XX) %*% Wt %*% XX) %*% t(XX) %*% Wt %*% Y
Beta
Yhat = XX %*% Beta

#Check with lsfit command
zz=lsfit(X, Y, wt=as.vector(w))
zz$coef

lines(X, Yhat, col="red")

}

#GCV of the original Model
zz=lsfit(X,Y)
epsilon = residuals(zz)
hat = XX %*% solve(t(XX) %*% XX) %*% t(XX)
GCV=sum(epsilon*epsilon)*N/((N - sum(diag(hat)))^2)
GCV
#[1] 0.02910873

#GCV of the last model
Yhat = XX %*% Beta
epsilon = Y - Yhat
hat = XX %*% solve(t(XX) %*% Wt %*% XX) %*% t(XX) %*% Wt
GCV=sum(epsilon*epsilon)*N/((N - sum(diag(hat)))^2)
GCV
#[1] 0.029204

#Very little difference between the two..