# Calculate the smoother matrix of a spline
  # The fitted values of a spline (or any linear smoother) are given by
  # m = w y
  # where y is n*1 matrix of observed responses, m is the n*1 matrix of
  # fitted values, and w is the n*n smoother matrix
    # Other names for w = smoothing matrix, influence matrix, hat matrix
  # In a smoothing spline, w will depend on the x values and on the curvature
  # penalty lambda
  # This function calculates/recovers the w matrix for a smoothing spline
  # How it works is suppressed to not spoil a part of the homework
# Inputs: a fitted smoothing spline model (a.spline)
  # the vector of x values used for fitting, in the original order (x)
    # a.spline$x will be re-organized into increasing order for plotting
# Output: a square matrix with one row per entry in x
# Presumes: a.spline was estimated by the smooth.spline function
  # no non-default settings chosen when estimating a.spline
  # no missing values in x
smoother.matrix <- function(a.spline, x) {
    n <- length(x)
    w <- matrix(0, nrow=n, ncol=n)
    for (i in 1:n) {
        y <- rep_len(0, n)  # Equivalent to rep(0, length.out=n) but faster
        y[i] <- 1
        w[,i] <- fitted(smooth.spline(x, y, lambda=a.spline$lambda))
    }
    return(w)
}