3
$\begingroup$

I am attempting to write R code for cubic splines to connect points on a graph. Specifically, I am attempting to reproduce Figure 3.3 of Wood (2006) Generalized Additive Models: An Introduction with R (page 124) where he uses seven cubic splines to connect eight points. I am starting with only three or four data points and have modified their values, so they may not closely match his data.

I have written R code based on equations provided at Wikipedia.

https://en.wikipedia.org/wiki/Spline_interpolation

I believe I have successfully programmed the equations provided in Wikipedia for the example there with three data points. However, I am unsure whether I have correctly modified that code for four data points.

What is the name of the algorithm shown in the above Wikipedia article and used in the example there? Is there a readily accessible citation for that algorithm?

I cannot see where the equations in the Wikipedia example follow from the derivation of the algorithm. Specifically, where do the equations for a11 and a22 come from and what do a11 and a22 represent?

The plot created by my R code for three data points appears to match the plot provided in the Wikipedia example. However, I cannot reproduce that plot using the spline function in R.

Wolphram-Alpha creates a plot that appears to match the plot returned by the FMM method in the spline function in R. What does FMM stand for and what would be a good citation for that method?

I suspect the reason the plot I obtain and the example plot in Wikipedia does not match the plot returned by Wolphram-Alpha or by the FMM function in R is that the example in Wikipedia uses two different cubic splines joined together at the middle data point, but that Wolphram-Alpha / FMM in R create one polynomial equation using all data points simultaneously. Is that correct?

In other words, I suspect Wikipedia and my R code use spline interpolation, but that Wolfram-Alpha and the FMM method in R use polynomial interpolation. Is that correct?

Originally I was hoping for verification that I had generated the correct equations for four data points based on the equations in Wikipedia for three data points. However, I have now shifted my focus to the questions in the preceding paragraphs and removed my R code.

If I can learn the answers to the questions above perhaps I can make progress following the algorithm shown in the above Wikipedia article and then perhaps post a programming question on Stack Overflow. Thank you for any answers to those questions.

EDIT

So far the best reference I have found is Numerical Computing with MATLAB by Cleve Moler, available here:

http://uk.mathworks.com/moler/chapters.html

If I make substantial progress I will post what I learn as an answer.

Improve this question
$\endgroup$
10
  • 1
    $\begingroup$ You should probably reconsider the whole idea of fitting a cubic spline to 4 points, which will almost certainly result in overfitting. $\endgroup$
    – dsaxton
    Jul 11, 2015 at 22:07
  • 2
    $\begingroup$ This is not a site for identifying problems in code. If you want to present details of your approach as algebra, it may be on topic, but if your question is purely one of programming and debugging, it will likely be off topic here. More broadly, you may be better off describing the underlying problem that you want to use splines to solve. $\endgroup$
    – Glen_b
    Jul 12, 2015 at 0:41
  • $\begingroup$ Gaussian processes are a flexible class of models to interpolate between points. More specifically, this seems like a case where full Bayes would be appropriate: small sample size, uncertainty about underlying function, etc. $\endgroup$
    – Sycorax ♦
    Jul 12, 2015 at 2:16
  • $\begingroup$ @Glen_b On page 124 of Wood (2006) Generalized Additive Models in Figure 3.3 he uses cubic splines to connect eight points. I was attempting to recreate that figure starting with four points. In posting this I was hoping someone could verify that I was using the correct equations, which I think are obtained from first and second derivatives. If this is viewed as a programming issue perhaps the post could be moved to Stack Overflow? $\endgroup$
    – Mark Miller
    Jul 12, 2015 at 11:00
  • 1
    $\begingroup$ In addition to this by Moler mathworks.com/moler/interp.pdf , Chapter by de Boor ftp.cs.wisc.edu/Approx/bsplbasic.pdf Book by de Boor amazon.com/Practical-Splines-Applied-Mathematical-Sciences/dp/… After you get done with that, you should know a fair bit about splines. $\endgroup$
    – Mark L. Stone
    Jul 12, 2015 at 21:02

1 Answer 1

Reset to default
1
$\begingroup$

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])))
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.