# R FUNCTIONS FOR WEEK 2 DEMONSTRATION. # GAUSSIAN BASIS FUNCTION COMPUTATION. Takes a matrix of inputs for # training cases, the width of the Gaussian basis functions, and the # interval the basis functions are spaced at (default same as width) # as arguments. Returns the matrix of Gaussian basis function values # for these training cases. Phi <- function (x.train, s, b=s) { mu <- seq (-2*s, 1+2*s, by=b) Phi <- matrix (NA, length(x.train), length(mu)) for (j in seq_along(mu)) { Phi[,j] <- exp(-0.5*(x.train-mu[j])^2/s^2) } Phi } # FUNCTION FOR PENALIZED LEAST SQUARES ESTIMATION. Takes a matrix of # basis function values for training cases, a vector of responses # values for training cases, and a penalty magnitude as arguments. # Returns the penalized least squares estimate of the parameters, with # quadratic penalty. The basis functions passed should not include # the all-one basis function, which is added automatically here, with # its coefficient not penalized. # # There are faster and more accurate ways of finding the estimate than # than by explicitly finding the inverse, as here, but this way is OK # in most circumstances. penalized.least.squares <- function (Phi, y.train, lambda) { Phi <- cbind(1,Phi) Istar <- diag(ncol(Phi)) Istar[1,1] <- 0 S <- lambda*Istar + t(Phi)%*%Phi solve(S) %*% t(Phi) %*% y.train } # FIND THE VALIDATION SQUARED ERROR. The arguments are the matrix of # basis function values in the training set, the vector of responses # in the training set, the penalty magnitude, and a vector of indexes # of validation cases. val.err <- function (Phi, y.train, lambda, val.ix) { beta <- penalized.least.squares (Phi[-val.ix,], y.train[-val.ix], lambda) predictions <- cbind(1,Phi[val.ix,]) %*% beta mean ((y.train[val.ix]-predictions)^2) } # FUNCTION TO FIND ARRAY OF VALIDATION ERRORS FOR VARIOUS LAMBDA AND S. val.array <- function (x.train, y.train, val.ix, try.lambda, try.s) { V <- matrix (NA, length(try.s), length(try.lambda)) rownames(V) <- paste("s=",try.s,sep="") colnames(V) <- paste("lambda=",try.lambda,sep="") for (i in 1:length(try.s)) { Phi.s <- Phi(x.train,try.s[i]) for (j in 1:length(try.lambda)) { V[i,j] <- val.err (Phi.s, y.train, try.lambda[j], val.ix) } } V }