#
# Problem on Page 191
#
source('exercise_5_1.R') # Load data and does computations from the example in Section 5.1.2.
y_mean = 0.0
y_std = 1.0 # assumed
c_names =c ('X1', 'X2')

#
# PC regression:
#

# Eigen analysis:
#
ES = eigen(XtX)
print('Eigenvalues')
print(round(ES$values, 3))
P = ES$vectors
print('Eigenvectors')
print(round(P, 3))

beta_hat = solve( XtX, Xty )
beta_c = beta_hat
PT_times_beta_s = t(P) %*% beta_c
PT_times_beta_s[2] = 0. # [ I_k, 0; 0, 0 ] times this is only taking the first row as we are dropping the eigenvector with the smallest eigenvalue
beta_pc_c = ( P %*% PT_times_beta_s ) * y_std
print(beta_pc_c)

pt_1 = sprintf('PC model: Y_PC = %+.3f', y_mean)
pt_2 = paste(sprintf('%+.3f %s', beta_pc_c, c_names), sep='', collapse=' ')
print(sprintf('%s %s', pt_1, pt_2))
print(sprintf('Model PC: y_hat= y(0.2, 0.7)= %f', sum( beta_pc_c * c( 0.2, 0.7 ) )))
print(sprintf('Model PC: y_hat= y(0.5, 0.5)= %f', sum( beta_pc_c * c( 0.5, 0.5 ) )))

#
# Ridge regression:
#
cv = 0.05

beta_ridge_z = solve(XtX + cv*diag(m), Xty) # the coefficients of the standardized predictors

pt_1 = sprintf('Ridge model (c=%.2f): Y= %.3f', cv, y_mean)
pt_2 = paste(sprintf('%+.3f %s', beta_ridge_z, c_names), sep='', collapse=' ')
print(sprintf('%s %s', pt_1, pt_2))

# The predictions with this model:
#
print(sprintf('Ridge model (c=%.2f): y_hat= y(0.2, 0.7)= %f', cv, sum( beta_ridge_z * c( 0.2, 0.7 ) )))
print(sprintf('Ridge model (c=%.2f): y_hat= y(0.5, 0.5)= %f', cv, sum( beta_ridge_z * c( 0.5, 0.5 ) )))