# # Written by: # -- # John L. Weatherwax 2009-04-21 # # email: wax@alum.mit.edu # # Please send comments and especially bug reports to the # above email address. # # EPage 314 # #----- if( ! require("ISLR") ){ install.packages("ISLR") } if( ! require("boot") ){ install.packages("boot") } vt100ClearScreen <- function(...) cat("\033[2J") vt100ClearScreen() set.seed(0) p = 100 n = 1000 # Generate some regression coefficients beta_0, beta_1, ..., beta_p beta_truth = rnorm( p+1 ) # Generate some data (append a column of ones): Xs = c( rep(1,n), rnorm( n*p ) ) X = matrix( data=Xs, nrow=n, ncol=(p+1), byrow=FALSE ) # Produce the response: Y = X %*% beta_truth + 0.1 * rnorm(n) # Get the true estimated coefficient estimates using lm: # m = lm( Y ~ X - 1 ) beta_lm = m\$coeff # Estimate beta_i using backfitting: # beta_hat = rnorm( p+1 ) # initial estimate of beta's is taken to be random n_iters = 10 beta_estimates = matrix( data=rep(NA,n_iters*(p+1)), nrow=n_iters, ncol=(p+1) ) beta_differences_with_truth = rep(NA,n_iters) beta_differences_with_LS = rep(NA,n_iters) for( ii in 1:n_iters ){ for( pi in 0:p ){ # for beta_0, beta_1, ... beta_pi ... beta_p # Perform simple linear regression on the variable X_pi (assuming we know all other values of beta_pi): # a = Y - X[,-(pi+1)] %*% beta_hat[-(pi+1)] # remove all predictors except beta_0 if( pi==0 ){ m = lm( a ~ 1 ) # estimate a constant beta_hat[pi+1] = m\$coef[1] }else{ m = lm( a ~ X[,pi+1] ) # estimate the slope on X_pi beta_hat[pi+1] = m\$coef[2] } } beta_estimates[ii,] = beta_hat beta_differences_with_truth[ii] = sqrt( sum( ( beta_hat - beta_truth )^2 ) ) beta_differences_with_LS[ii] = sqrt( sum( ( beta_hat - beta_lm )^2 ) ) } #postscript("../../WriteUp/Graphics/Chapter7/prob_12_convergence_plot.eps", onefile=FALSE, horizontal=FALSE) plot( 1:n_iters, beta_differences_with_LS, main='||beta_truth-beta_LS||_2', pch=19, type='b' ) grid() #dev.off()