R version 2.9.0 (2009-04-17)
Copyright (C) 2009 The R Foundation for Statistical Computing
ISBN 3-900051-07-0

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

>  # 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
 [1]  1  2  3  4  5  6  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27
[26] 28 29 30 31 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56
[51] 57 58 59 60 61
>    y <- all[,2] 
>    y
 [1]  74.590  78.720  75.931  82.470  89.840  84.000  88.050  88.620  89.620
[10]  94.480  90.950  98.230 101.130  95.904  97.545  99.482  97.585 100.448
[19]  96.629 104.114 104.162 104.863 106.240 109.069 113.580 117.200 115.035
[28] 114.277 115.117 114.820 116.338 119.814 121.327 124.002 126.244 128.922
[37] 128.740 130.840 128.209 128.490 135.601 133.791 138.857 138.767 139.130
[46] 140.293 143.137 147.350 151.388 144.317 151.207 152.882 156.867 155.749
[55] 157.735
>  # construct design matrix W
>  #   first some tools
>    h <- diff(knots)              # length k
>    lamsm <- h[-1]/(h[-1]+h[-k])  # length k-1
>    lamsm
[1] 0.8000000 0.5140196
>    lambig <- diag(2,kp1)                  # fill diagonal
>    diag(lambig[2:k,1:(k-1)]) <- 1-lamsm           # subdiagonal
>    diag(lambig[2:k,3:kp1]) <- lamsm               # superdiagonal
>    lambig[1,2] <- -4                    # end conditions
>    lambig[kp1,k] <- -4                  # for Poirier
>    lambig
     [,1]       [,2] [,3]      [,4]
[1,]  2.0 -4.0000000  0.0 0.0000000
[2,]  0.2  2.0000000  0.8 0.0000000
[3,]  0.0  0.4859804  2.0 0.5140196
[4,]  0.0  0.0000000 -4.0 2.0000000
>    theta <- matrix(0,kp1,kp1)
>    diag(theta[2:k,2:k]) <- -6/(h[-k]*h[-1])       # diagonal
>    diag(theta[2:k,1:(k-1)]) <- 6/(h[-k]*(h[-1]+h[-k]))  # subdiagonal
>    diag(theta[2:k,3:kp1]) <- 6/(h[-1]*(h[-1]+h[-k]))  # superdiagonal
>    theta   
           [,1]         [,2]         [,3]        [,4]
[1,] 0.00000000  0.000000000  0.000000000 0.000000000
[2,] 0.02840237 -0.035502959  0.007100592 0.000000000
[3,] 0.00000000  0.004313436 -0.008391578 0.004078142
[4,] 0.00000000  0.000000000  0.000000000 0.000000000
>  # 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
          [,1]
[1,]  76.03371
[2,]  86.70855
[3,] 116.11677
[4,] 158.64801
>    qrstuff <- qr(W)
>    backsolve(qrstuff$qr[1:kp1,1:kp1],qr.qty(qrstuff,y)[1:kp1])
[1]  76.03371  86.70855 116.11677 158.64801
>    sum(qr.qty(qrstuff,y)[(kp1+1):nobs]^2)
[1] 397.674
>    qr.qty(qrstuff,cbind(rep(1,nobs),x))[(kp1+1):nobs,]
                                x
 [1,]  1.387779e-17  3.552714e-15
 [2,]  1.942890e-16  5.329071e-15
 [3,]  5.551115e-16  3.552714e-15
 [4,]  2.081668e-16  2.664535e-15
 [5,]  3.191891e-16 -8.881784e-16
 [6,]  5.412337e-16  0.000000e+00
 [7,]  4.718448e-16 -8.881784e-16
 [8,]  3.469447e-16 -8.881784e-16
 [9,]  3.469447e-16 -5.329071e-15
[10,]  4.163336e-16 -7.105427e-15
[11,]  4.718448e-16  8.881784e-16
[12,]  2.983724e-16 -8.881784e-15
[13,]  4.163336e-16 -9.769963e-15
[14,]  1.908196e-16 -8.659740e-15
[15,]  2.792905e-16 -1.398881e-14
[16,]  1.179612e-16 -2.031708e-14
[17,]  4.302114e-16 -1.110223e-14
[18,]  2.844947e-16 -1.465494e-14
[19,]  2.220446e-16 -2.309264e-14
[20,]  2.636780e-16 -2.042810e-14
[21,] -1.387779e-17 -2.575717e-14
[22,]  1.387779e-16 -2.664535e-14
[23,]  1.387779e-16 -2.309264e-14
[24,]  1.110223e-16 -2.486900e-14
[25,] -3.053113e-16 -3.730349e-14
[26,] -3.885781e-16 -2.842171e-14
[27,] -1.387779e-16 -2.842171e-14
[28,] -2.775558e-16 -3.019807e-14
[29,] -2.498002e-16 -3.375078e-14
[30,] -4.996004e-16 -4.796163e-14
[31,] -5.551115e-16 -3.907985e-14
[32,] -4.163336e-16 -3.197442e-14
[33,] -5.828671e-16 -4.085621e-14
[34,] -4.024558e-16 -3.907985e-14
[35,] -6.383782e-16 -3.641532e-14
[36,] -5.412337e-16 -4.440892e-14
[37,] -4.705437e-16 -3.250178e-14
[38,] -2.359224e-16 -3.286260e-14
[39,] -3.053113e-16 -3.641532e-14
[40,] -3.053113e-16 -3.019807e-14
[41,] -1.387779e-16 -2.842171e-14
[42,] -2.220446e-16 -2.309264e-14
[43,] -1.110223e-16 -2.486900e-14
[44,] -1.110223e-16 -2.131628e-14
[45,]  5.551115e-17 -1.421085e-14
[46,]  0.000000e+00 -1.421085e-14
[47,]  1.110223e-16 -7.105427e-15
[48,]  1.110223e-16 -7.105427e-15
[49,]  1.110223e-16  0.000000e+00
[50,]  2.220446e-16  7.105427e-15
[51,]  2.220446e-16  7.105427e-15
>  # clean up & go
>   rm(list=ls())
>   q()