Chapter 3 31A.

###
Natural.CSpline = function(x, 
                           y, 
                           x.new){
  n = length(x)
  m = length(x.new)
  h = xvec[-1] - xvec[-n]
  A = matrix(rep(0,(n)^2), ncol=n)
  Y = rep(0,n)
  y.new = NULL
  Y[c(1,n)] = 0
  A[1,1] = 1
  A[n,n] = 1
  for (i in 2:(n-1)){
     A[i,(i-1):(i+1)] = c(h[i-1], 2*(h[i-1]+h[i]), h[i])
     Y[i] =3*(y[i+1]-y[i])/h[i]-3*(y[i]-y[i-1])/h[i-1]
  }
  cc = as.vector(solve(A)%*%Y)
  aa = y                # input y values
  ID = 1:(n-1)
  bb = dd = NULL
  for (i in 1:(n-1)){
    dd[i] = (cc[i+1] - cc[i])/(3*h[i])
    bb[i] = (aa[i+1] - aa[i])/h[i] - h[i]*(2*cc[i] + cc[i+1])/3
  }
  cc = cc[-n]
  aa = aa[-n]
  coef = cbind(ID = ID, aa = aa, bb = bb, cc = cc, dd = dd)
  #predicted value
  dS = NULL                       # derivative of the C-spline.
  for (j in 1:length(x.new)){
  for (i in 1:n){
    if (x.new[j]<= x[i]){
       next
      }
     y.new[j] = aa[i] + bb[i]*(x.new[j] - x[i]) + cc[i]*(x.new[j] - x[i])^2 + dd[i]*(x.new[j] - x[i])^3
     dS[j] = bb[i] + 2*cc[i]*(x.new[j] - x[i]) +3* dd[i]*(x.new[j] - x[i])^2
   }
  }
 list(coef = coef, y.new = y.new, d.S = dS)
}
xvec=c(0,6,10,13,17,20,28)
yvec=c(6.67,17.33,42.67,37.33,30.10,29.31,28.74)
###
Natural.CSpline(x = xvec, y = yvec, x.new = .25)
## $coef
##      ID    aa          bb            cc            dd
## [1,]  1  6.67 -0.44686998 -2.773082e-16  0.0617649070
## [2,]  2 17.33  6.22373997  1.111768e+00 -0.2709883292
## [3,]  3 42.67  2.11044677 -2.140092e+00  0.2810920124
## [4,]  4 37.33 -3.14061865  3.897365e-01 -0.0141142057
## [5,]  5 30.10 -0.70020864  2.203660e-01 -0.0249136386
## [6,]  6 29.31 -0.05068078 -3.856730e-03  0.0001606971
## 
## $y.new
## [1] 6.559248
## 
## $d.S
## [1] -0.4352891
xvec=c(0,6,10,13,17,20,28)
yvec=c(6.67,16.11,18.89,15,10.56,9.44,8.89)
###
Natural.CSpline(x = xvec, y = yvec, x.new = 0.25)
## $coef
##      ID    aa         bb            cc           dd
## [1,]  1  6.67  1.6628737 -8.964075e-17 -0.002487232
## [2,]  2 16.11  1.3942526 -4.477018e-02 -0.032510741
## [3,]  3 18.89 -0.5244245 -4.348991e-01  0.059161669
## [4,]  4 15.00 -1.5364539  9.755594e-02  0.002264382
## [5,]  5 10.56 -0.6473160  1.247285e-01 -0.011133651
## [6,]  6  9.44 -0.1995535  2.452566e-02 -0.001021902
## 
## $y.new
## [1] 7.08568
## 
## $d.S
## [1] 1.662407

Sample 1 polynomial 6.67 - 0.44686998x - 2.773082e-16x^2 + 0.0617649070x^3 for x [0,6] 17.33 + 6.22373997(x-6) + 1.111768e+00(x-6)^2 -0.2709883292(x-6)^3 for x [6,10] 42.67 + 2.11044677(x-10) -2.140092e+00(x-10)^2 + 0.2810920124(x-10)^3 for x [10,13] 37.33 -3.14061865(x-13) + 3.897365e-01(x-13)^2 -0.0141142057(x-13)^3 for x [13,17] 30.10 -0.70020864(x-17) + 2.203660e-01(x-17)^2 -0.0249136386(x-17)^3 for x [17,20] 29.31 -0.05068078(x-20) -3.856730e-03(x-20)^2 0.0001606971(x-20)^3 for x [20,28]

Sample 2 polynomial 6.67 + 1.6628737x -8.964075e-17x^2 -0.002487232x^3 for x [0,6] 16.11 +1.3942526(x-6) -4.477018e-02(x-6)^2 -0.032510741(x-6)^3 for x [6,10] 18.89 -0.5244245(x-10) -4.348991e-01(x-10)^2 + 0.059161669(x-10)^3 for x [10,13] 15.00 -1.5364539(x-13) + 9.755594e-02(x-13)^2 + 0.002264382(x-13)^3 for x [13,17] 10.56 -0.6473160(x-17) + 1.247285e-01(x-17)^2 -0.011133651(x-17)^3 for x [17,20] 9.44 -0.1995535(x-20) + 2.452566e-02(x-20)^2 -0.001021902(x-20)^3 for x [20,28]

31B. Sample 1 max average value

[0,6] is 17.33 on day 6 [6,10] is 31.832 on day 8.867 [10,13] is 43.23 on day 10.553 [13,17] is 37.33 on day 13 [17,20] is 20.1 on day 17 [20,28] is 29.31 on day 20

Sample 2 max average value

[0,6] is 16.11 on day 6 [6,10] is 19.023 on day 9.323 [10,13] is 18.89 on day 10 [13.17] is 15 on day 13 [17.20] is 10.56 on day 17 [20,28] is 9.44 on day 20

  1. Take the derivative of S_0(x) which is S’_0(x) = 2-3x^2 and take 2nd derivative which is -6x Then take the derivative of S_1(x) which is s’_1(x) = b+2c(x-1)+3d(x-1)^2 and take 2nd derivative which is 2c+6d(x-1) Now let 2-3x^2 = b+2c(x-1)+3d(x-1)^2 and let x = 1 which makes b = -1 Now let -6x = 2c+6d(x-1) and let x = 1 which makes c = -3 Now we find for d by making -6x = 2(-3)+6d(x-1) while x = 0 which gives d = 1 b= -1 c = -3 d = 1

3C.

xvec=c(-.5,-.25,0)
yvec=c(-.0247500,.3349375,1.101)
###
Natural.CSpline(x = xvec, y = yvec, x.new = 0.25)
## $coef
##      ID         aa       bb     cc     dd
## [1,]  1 -0.0247500 1.032375 0.0000  6.502
## [2,]  2  0.3349375 2.251500 4.8765 -6.502
## 
## $y.new
## [1] NA
## 
## $d.S
## [1] NA

The polynomial is -0.0247500 + 1.032375(x+.5) + 6.502(x+.5)^3 for x [-.5,-.25] 0.3349375 + 2.251500(x+.25) + 4.8765(x+.25)^2 -6.502(x+.25)^3 [-.25,0] 3D.

xvec=c(.1,.2,.3,.4)
yvec=c(-.62049958,-0.28398668,0.00660095,0.24842440)
###
Natural.CSpline(x = xvec, y = yvec, x.new = 0.25)
## $coef
##      ID          aa       bb        cc         dd
## [1,]  1 -0.62049958 3.455087  0.000000 -8.9957933
## [2,]  2 -0.28398668 3.185213 -2.698738 -0.9463033
## [3,]  3  0.00660095 2.617076 -2.982629  9.9420967
## 
## $y.new
## [1] -0.1315912
## 
## $d.S
## [1] 2.908242

-0.62049958 + 3.455087(x-.1) -8.9957933(x-.1)^3 for x [.1,.2] -0.28398668 + 3.185213(x-.2) -2.698738(x-.2)^2 -0.9463033(x-.2)^3 for x [.2,.3] 0.00660095 + 2.617076(x-.3) -2.982629(x-.3)^2 + 9.9420967(x-.3)^3 for x [.3,.4]