# splrgm.r          cubic spline interpolation
 #
   kp1 <- 4                        # knots of spline including endpts
   k <- kp1 - 1
   nobs <- 55
   knots <- c(1,7.5,33.5,61.0001)
   all <-read.table("indy500.dat",nrows=55)
   x <- all[,1] - 1910               # subtract, start at 1911
   x
   y <- all[,2] 
   y
 # construct design matrix W
 #   first some tools
   h <- diff(knots)              # length k
   lamsm <- h[-1]/(h[-1]+h[-k])  # length k-1
   lamsm
   lambig <- matrix(0,kp1,kp1)        # zero matrix
   goofy <- matrix(0,kp1+1,kp1)       # zero matrix
   goofy[1,(1:(k-1))] <- 1-lamsm 
   goofy[2,(1:(k-1))] <- 2 
   goofy[3,(1:(k-1))] <- lamsm 
   lambig <- matrix(c(rep(0,kp1),goofy),kp1,kp1,byrow=TRUE)
   diag(lambig) <- rep(2,kp1)
   lambig[1,2] <- -4                    # end conditions
   lambig[kp1,k] <- -4                  # for Poirier
   lambig
   goofy <- matrix(0,kp1+1,kp1)       # zero matrix
   goofy[1,(1:(k-1))] <- 6/(h[-k]*(h[-1]+h[-k])) 
   goofy[2,(1:(k-1))] <- -6/(h[-k]*h[-1])
   goofy[3,(1:(k-1))] <- 6/(h[-1]*(h[-1]+h[-k])) 
   theta <- matrix(c(rep(0,kp1),goofy,rep(0,kp1)),kp1,kp1,byrow=TRUE)
   theta   
 # now P and Q matrices
   P <- matrix(0,nobs,kp1)
   Q <- matrix(0,nobs,kp1)
   for (i in(1:nobs)) {
   xi <- x[i]
   j <- findInterval(xi,knots)
   P[i,j] <- (knots[j+1]-xi)*((knots[j+1]-xi)^2-h[j]^2)/(6*h[j])
   P[i,j+1] <- (xi-knots[j])*((xi-knots[j])^2-h[j]^2)/(6*h[j])
   Q[i,j] <- (knots[j+1]-xi)/h[j]
   Q[i,j+1] <- (xi-knots[j])/h[j]
                    }
   W <- P%*%solve(lambig,theta) + Q
 # regression computations
   solve(t(W)%*%W)%*%t(W)%*%y
   qrstuff <- qr(W)
   backsolve(qrstuff$qr[1:kp1,1:kp1],qr.qty(qrstuff,y)[1:kp1])
   sum(qr.qty(qrstuff,y)[(kp1+1):nobs]^2)
   qr.qty(qrstuff,cbind(rep(1,nobs),x))[(kp1+1):nobs,]
 # clean up & go
  rm(list=ls())
  q()