# 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
}