Can someone explain to me what a cubic spline is, and how we could use it to interpolate a function?
I have searched on the internet but I would like a simple explanation.
Solved – Explanation of cubic spline interpolation
interpolationnumericssplines
Related Solutions
This paper presented at UseR! 2009 seems to address a similar problem
http://www.r-project.org/conferences/useR-2009/slides/Roustant+Ginsbourger+Deville.pdf
It suggests the DiceKriging package http://cran.r-project.org/web/packages/DiceKriging/index.html
In particular, check the functions km and predict.
Here is a an example of three dimensional interpolation. It looks to be straightforward to generalise.
x <- c(0, 0.4, 0.6, 0.8, 1)
y <- c(0, 0.2, 0.3, 0.4, 0.5)
z <- c(0, 0.3, 0.4, 0.6, 0.8)
model <- function(param){
2*param[1] + 3*param[2] +4*param[3]
}
model.in <- expand.grid(x,y,z)
names(model.in) <- c('x','y','z')
model.out <- apply(model.in, 1, model)
# fit a kriging model
m.1 <- km(design=model.in, response=model.out, covtype="matern5_2")
# estimate a response
interp <- predict(m.1, newdata=data.frame(x=0.5, y=0.5, z=0.5), type="UK", se.compute=FALSE)
# check against model output
interp$mean
# [1] 4.498902
model(c(0.5,0.5,0.5))
# [1] 4.5
# check we get back what we put in
interp <- predict(m.1, newdata=model.in, type="UK", se.compute=FALSE)
all.equal(model.out, interp$mean)
# TRUE
I have made progress toward answers to my questions. Here I post what I have learned so far and I can modify my post as I learn more.
Firstly, the best reference I have found so far is Numerical Computing with MATLAB by Cleve Moler, as I mentioned in an edit to the question and as Mark L. Stone mentioned in a comment.
I have not yet derived the algorithm shown in Wikipedia. Nevertheless, I assumed my R code was correct and expanded that code to handle the eight points Wood (2006:124) used in Figure 3.3 of Generalized Additive Models: An Introduction with R.
Those eight points, to the best of my eyesight, are:
my.data <- read.table(text = '
x y
0.04 0.60
0.18 0.01
0.32 0.30
0.415 0.28
0.58 0.79
0.68 0.27
0.81 0.70
0.88 0.83
', header = TRUE, stringsAsFactors = FALSE, na.strings = 'NA')
In the process I discovered an error in my R code, described below. After correcting that error I discovered that the example in Wikipedia, and my R code, does indeed reproduce Wood's (2006) Figure 3.3. I was also able to reproduce that same plot using this built-in R command:
lines(spline(x, y, method = c("natural"), n = 100), lwd = 2, col = 'blue')
More specifically:
spline(x, y, method = c("natural"), n = 100)
In other words, the Wikipedia example is using the natural cubic spline as near as I can tell. Not the FMM spline. The FMM spline seems to stand for the Forsythe, Malcolm, Moler spline described in Forsythe, G. E., Malcolm, M. A. and Moler, C. B. (1977) Computer Methods for Mathematical Computations. Wiley.
Wolfram-Alpha does not reproduce that plot. At least I have not been able to do it. When I use this command in Wolfram-Alpha:
interpolating polynomial {0.04,0.60},{0.18,0.01},{0.32,0.30},{0.415,0.28},{0.58,0.79},{0.68,0.27},{0.81,0.70},{0.88,0.83}
Wolfram-Alpha returns a single polynomial:
6.29713 - 221.199 x + 2402.99 x^2 - 12051.2 x^3 + 32012.2 x^4 - 46374.4 x^5 + 34557.4 x^6 - 10355.5 x^7
Note that this is a 7th-order polynomial, not a cubic polynomial.
I did try to create a 3rd-order polynomial in Wolfram-Alpha using:
BSplineCurve[{{0.04,0.60},{0.18,0.01},{0.32,0.30},{0.415,0.28},{0.58,0.79},{0.68,0.27},{0.81,0.70},{0.88,0.83}}, SplineDegree -> 3]
However, this returns no results. This might be Mathematica code. Regardless, the Wolfram-Alpha code I have tried does not return a natural cubic spline.
Once I derive the example equations in Wikipedia I will post my interpretation of them, in particular regarding a11 and a22. My limited understanding of the difference between the splines R offers is that it has to do with how they treat the second derivative of the first and last cubic spline.
Lastly, for now, regarding the error in my R code, the line:
a33 = 2 * ((1 / (x[3] - x[2])) + (1 / (x[3] - x[2])))
should have been.
a33 = 2 * ((1 / (x[3] - x[2])) + (1 / (x[4] - x[3])))
Best Answer
If you have a function $f(x)$ on some interval $[a,b]$, which is divided on $[x_{i-1}, x_i]$ such as $a=x_0< x_1< ... <x_N=b$ then you can interpolate this function by a cubic spline $S(x)$. $S(x)$ is a piecewise function: on each $h_i = x_i - x_{i-1}$ it's a cubic polynomial, which can be written for simplicity as $S_i(x) = a_i + b_i(x - x_i) + {c_i\over2}(x-x_i)^2 + {d_i\over6}(x - x_i)^3 \,\!$. It has to satisfy the next constraints:
1) passing through the knots : $S_i\left(x_{i}\right) = f(x_{i})$
2) be continuous up to the 2nd derivative:
$S_i\left(x_{i-1}\right) = S_{i-1}(x_{i-1}) \\ S'_i\left(x_{i-1}\right) = S'_{i-1}(x_{i-1}) \\ S''_i\left(x_{i-1}\right) = S''_{i-1}(x_{i-1})$
3) for natural splines:
$S''(a) = S''(b) = 0.$
These equations will uniquely define spline coefficients.
A good way to understand this is to take e.g. 3 points and manually solve systems for coefficients of $S_1(x)$ and $S_2(x)$
Finally you should get the next system:
$a_i = f\left(x_{i}\right) \,\!$
$h_ic_{i-1} + 2(h_i + h_{i+1})c_i + h_{i+1}c_{i+1} = 6\left({{f_{i+1} - f_i}\over{h_{i+1}}} - {{f_{i} - f_{i-1}}\over{h_{i}}}\right) \,\!$
$d_i = {{c_i - c_{i-1}}\over{h_i}} \,\!$
$b_i = {1\over2}h_ic_i - {1\over6}h_i^2d_i + {{f_i - f_{i-1}}\over{h_i}}= {{f_i - f_{i-1}}\over{h_i}} + {{h_i(2c_i + c_{i-1})}\over6} \,\!$