Re: smoothing matrix

[Martin Maechler <maechler@stat.math.ethz.ch>]

>>>>> "Lin" == L Y Hin <lyhin@netvigator.com> writes:

    Lin> Dear S-plus users,
    Lin> I have a problem with smoothing spline that I'd be grateful for your
    Lin> advice.
    Lin> I've been advised that I can obtain the smoother matrix, ie., S=(I +
    Lin> lambda*K)^-1, from the smooth.spline() function as below:

(I'm using " <- " instead of "_" which is already more readable)

    Lin> x <- 1:20
    Lin> hat <- diag(x)
    Lin> spar <- smooth.spline(x)$spar
    Lin> for (i in 1:length(x))
    Lin> hat[,i] <- smooth.spline(smooth.spline,hat[,i],spar=spar)$y

which probably should rather be

         hat[,i] <- smooth.spline(x, hat[,i],spar=spar)$y
         

    Lin> and hat should be the smoother matrix S.


    Lin> However, this does not work for cases with ties, eg.,

    Lin> x <- rep(1:5,(2:6))

    Lin> Kindly please advise me whether there is any way for me to obtain the
    Lin> smoother matrix S as shown in the Reinsch form above?

Yes, the following is even more generally useful; and I have
only recently patched it to work for non-unique x such that
dim(S) |-> c(n,n)  where n <- length(x)

It was written for R but should work in all implementations of
S, incl. S-PLUS:
-------------------------------------------------------------------------

## for R :
require(modreg)

hatMat <- function(x, trace = FALSE,
                   pred.sm = function(x,y,...)
                   predict(smooth.spline(x,y, ...), x = x)$y,
                   ...)
{
    ## Purpose: Return Hat matrix of a smoother -- very general (but slow)
    ## -------------------------------------------------------------------------
    ## Arguments: x: the design points
    ##          trace: if TRUE, only the trace of the matrix is returned, i.e.,
    ##                  sum(diag( <matrix> ))
    ##          pred.sm: a function(x,y,..), predict(smoother(x,y,...))
    ##          ....   : smoothing parameters, for the pred.sm() function.
    ## -------------------------------------------------------------------------
    ## Author: Martin Maechler, Date:  7 Mar 2001, 11:12
    n <- length(x)
    y <- pred.sm(x, numeric(n), ...)
    if(!is.numeric(y) || length(y) !=n)
        stop("`pred.sm' does not return a numeric length n vector")
    H <- if(trace) 0 else matrix(as.numeric(NA), n,n)
    for (i in 1:n) {
        y <- numeric(n) ; y[i] <- 1
        y <- pred.sm(x, y, ...)
        if(trace) H <- H + y[i] else H[,i] <- y
    }
    H
}

## Example `pred.sm' arguments for hatMat() :
pspl <- function(x,y,...) predict(smooth.spline(x,y, ...), x = x)$y
ploes <-function(x,y,...) predict(loess(y ~ x, ...))
pksm <- function(x,y,...) ksmooth(x,y,"normal", x.points=x, ...)$y
pRmean <- function(x,y,...) run.mean(y, ...)
pRline <- function(x,y,...) run.line(x,y, ...)$y

all.equal(sum(diag((hatMat(c(1:4, 7:12), df = 4)))),
                    hatMat(c(1:4, 7:12), df = 4, trace = TRUE), tol = 1e-12)
## TRUE

----

I hope this helps; it's really not slow anymore on current
computers and recent versions of R, as long as n <- length(x) is
below 100, say.
E.g with the newly recent R 1.4 and a Linux Pentium III, 900MHz:

> x <- sort(runif(100))
> system.time(H <- hatMat(x, spar = .5))
[1] 1.70 0.04 2.25 0.00 0.00

{where system.time is I think unix.time() or dos.time() for S+}