# October 22, 2002

# The brute force cubic spline smoother

cubic.sp <- function(lam,x,y,tk,k,br=F,cex=1) {

    if(missing(tk)) { tk<-seq(min(x)-.001,max(x)+.001,l=k+1)} else {k<-length(tk)-1}
    plot(x,y,pch=16,cex=cex,xlim=range(x,tk),main=paste("lambda = ",as.character(lam)))
        abline(v=tk,lty=2); n<-length(x)
    hk <- diff(tk); if(length(y)!=n) stop("y wrong length"); k1<-k-1

			# the big A matrix
    A <- matrix(0,n,4*k); sk <- seq(0,by=k,length=4)
    for(i in seq(n)) { j <- max(seq(k)[x[i]-tk[1:k]>=0]); z<-x[i]-tk[j]
	A[i,j+sk] <- c(1,z,z^2/2,z^3/6) }

    Z <- matrix(0,2*k,2*k)
    W12 <- diag(hk^2)/2; W11 <- diag(hk); W22 <- diag(hk^3)/3
    W <- rbind( cbind(Z,Z), cbind(Z, rbind( cbind(W11,W12), cbind(W12,W22) ) ) )

B11 <- cbind(diag(k-1),0); B11[col(B11)-row(B11)==1] <- -1
B12 <- cbind(diag(hk[-k],nrow=k1),0); B13 <- cbind(diag(hk[-k]^2,nrow=k1)/2,0)
B14 <- cbind(diag(hk[-k]^3,nrow=k1)/6,0); Z <- matrix(0,k-1,k)

B <-rbind( cbind(B11,B12,B13,B14), cbind(Z,B11,B12,B13), cbind(  Z,  Z,B11,B12) )

browser()
Ti <- solve(t(A)%*%A+lam*W)
mu <- 2*solve(B%*%Ti%*%t(B), B%*%Ti%*%t(A)%*%y)

th <- Ti %*% ( t(A)%*%y - .5*t(B)%*%mu )

  for(j in 1:k) { xx <- seq(0,hk[j],l=25)
   yy <- th[j]+th[j+k]*xx+th[j+2*k]*xx^2/2+th[j+3*k]*xx^3/6
   lines(xx+tk[j],yy,col=2,lwd=2) }
 if(br)browser()
}

par(err=-1)

ll = c(10,1,.1,.01,.001,.0001,.00001,.000001,.0000001)

for(lam in ll) { cubic.sp(lam,c(-3,-2,-1,0,1,2,3),c(-27,-8,-1,0,1,8,27),c(-4,0,4),br=T) }

for(lam in ll) { cubic.sp(lam,seq(.05,.95,.05), sqrt(seq(.05,.95,.05)), c(0,.2,.5,1),br=T) }

set.seed(425); xx<-runif(100); yy <- 25*(xx*(1-xx))^3 + rnorm(100,0,.05)

for(lam in ll) { cubic.sp(lam,xx,yy,c(0,.2,.5,1),br=T) }
for(lam in ll) { cubic.sp(lam,xx,yy,seq(0,1,l=11),br=T) }
for(lam in ll) { cubic.sp(lam,xx,yy,sort(xx)[seq(1,length(xx),length=15)],br=T) }

set.seed(427); xx<-runif(1000,-1,1); yy <- 1/(1+exp(-10*xx)) + rnorm(1000,0,.05)
for(lam in ll) { cubic.sp(lam,xx,yy,seq(-1,1,.5),br=T,cex=.4) }
for(lam in ll) { cubic.sp(lam,xx,yy,seq(-1,1,.1),br=T,cex=.4) }

   mesh = c( seq(-1,-.3,,3), seq(-.2,.2,.1), seq(.3,1,,3) )
for(lam in ll) { cubic.sp(lam,xx,yy,mesh,br=T,cex=.4) }