Rで3次スプライン補間の実装

以下のように実装しました。原理は省略（後日書き足すかも）。

spline3.R

spline3 <- function(X, Y){
  N <- length(X) -1 
# N-1*N-1の行列Hを定義
  H <- matrix(0, nrow=(N-1), ncol=(N-1))
# ujを求める
  h_0 <- X[2]-X[1]
  h_1 <- X[3]-X[2]
  h_Nm2 <- X[N] - X[N-1]
  h_Nm1 <- X[N+1] - X[N]
  H[1,1:2] <- c(2*(h_0+h_1), h_1)
  H[N-1, c(N-2, N-1)] <- c(h_Nm2, 2*(h_Nm2 + h_Nm1))
  vi <- rep(0, N-1)
  vi[1]   <- 6 *((Y[3] - Y[2]) / h_1 - (Y[2] - Y[1]) / h_0)
  vi[N-1]  <- 6 *((Y[N+1] - Y[N]) / h_Nm1 - (Y[N] - Y[N-1]) / h_Nm2)
 
  for(i in 2:(N-2)){
    h_im1 <- X[i+1] - X[i]
    h_i   <- X[i+2] - X[i+1]
    vi[i]   <- 6*((Y[i+2] - Y[i+1]) / hi - (Y[i+1] - Y[i]) / hi_1)
    H[i, c(i-1, i, i+1)] <- c(h_im1, 2*(h_im1+h_i), h_i)
  }
  H_1 <- solve(H)
  u_ <- H_1 %*% vi
  u  <- c(0, u_, 0)
# 求まったujからaj, bj, cj, djを計算
  a <- NULL; b <- NULL; c <- NULL; d <- NULL
  for(j in 1:N){
    a[j] <- (u[j+1] - u[j]) / (X[j+1] - X[j]) / 6
    b[j] <- u[j]/2
    c[j] <- (Y[j+1] - Y[j]) / (X[j+1] - X[j]) - (X[j+1] - X[j]) * (u[j+1] + 2*u[j]) / 6
    d[j] <- Y[j]
  }
  output <- data.frame(a, b, c, d)
  return(output)
  }

実は、spline()関数を使えば自分で実装しなくてもよいです。
rnorm()で乱数を20個生成してプロットし、19区間のそれぞれに対して係数a,b,c,dを上のspline3()関数から計算して描画し、最後にRのspline()関数で計算した分（実線）も重ね描きしました。区間以外にも点線で描画されてしまってごちゃごちゃしていますが、実線と点線が重なっていることがわかります。

参考にしたページ
http://www.yamamo10.jp/yamamoto/lecture/2006/5E/interpolation/interpolation_html/node3.html