From: Spencer Graves <spencer.graves_at_pdf.com>

Date: Sun, 20 Jul 2008 08:51:31 -0700

*>
*

> Duncan Murdoch

*>
*

*>>
*

*>> However, computations using natural splines are numerically
*

*>> unstable. The standard solution to this problem is to use B-splines,
*

*>> which are 0 outside a finite interval.
*

*>> Let's look at your example:
*

*>> n <- 9
*

*>> x <- 1:n
*

*>> y <- rnorm(n)
*

*>> plot(x, y, main = paste("spline[fun](.) through", n, "points"))
*

*>> spl <- smooth.spline(x,y)
*

*>> lines(spl)
*

*>>
*

*>> The 'smooth.spline' function uses B-splines. To see what they
*

*>> look like, let's do the following:
*

*>> library(fda)
*

*>> Bspl.basis <- create.bspline.basis(unique(spl$fit$knot))
*

*>>
*

*>> # Check to make sure: all.equal(knots(Bspl.basis, interior=FALSE),
*

*>> spl$fit$knot)
*

*>> # TRUE
*

*>>
*

*>> # What do B-splines look like? plot(Bspl.basis)
*

*>> abline(v=knots(Bspl.basis), lty='dotted', col='red')
*

*>> # 7 interior knots, 2 end knots replicated 4 times each, for a
*

*>> spline of order 4, degree 3 (cubic splines) # total of 15 knots
*

*>> # Each spline uses 5 consecutive knots, which means there will be 11
*

*>> basis functions. # NOTE: 'smooth.spline' rescaled the interval
*

*>> [1, 9] to [0, 1]. # Evaluate the 11 B-splines at 'x'
*

*>> Bspl.basis.x <- eval.basis((x-1)/8, Bspl.basis)
*

*>>
*

*>> round(Bspl.basis.x, 4)
*

*>>
*

*>> # Now the manual computation: y.spl <- Bspl.basis.x %*% spl$fit$coef
*

*>>
*

*>> # Plot to confirm: plot(x, y, main = paste("spline[fun](.) through",
*

*>> n, "points"))
*

*>> spl.xy <- spline(x, y)
*

*>> lines(spl.xy)
*

*>> points(x, y.spl, pch=2, col='red')
*

*>>
*

*>> Hope this helps. Spencer
*

*>>
*

*>> rkevinburton_at_charter.net wrote:
*

*>>> Fair enough. FOr a spline interpolation I can do the following:
*

*>>>
*

*>>>
*

*>>>> n <- 9
*

*>>>> x <- 1:n
*

*>>>> y <- rnorm(n)
*

*>>>> plot(x, y, main = paste("spline[fun](.) through", n, "points"))
*

*>>>> lines(spline(x, y))
*

*>>>>
*

*>>> Then look at the coefficients generated as:
*

*>>>
*

*>>>
*

*>>>> f <- splinefun(x, y)
*

*>>>> ls(envir = environment(f))
*

*>>>>
*

*>>> [1] "ties" "ux" "z"
*

*>>>> splinecoef <- get("z", envir = environment(f))
*

*>>>> slinecoef
*

*>>>>
*

*>>> $method
*

*>>> [1] 3
*

*>>>
*

*>>> $n
*

*>>> [1] 9
*

*>>>
*

*>>> $x
*

*>>> [1] 1 2 3 4 5 6 7 8 9
*

*>>>
*

*>>> $y
*

*>>> [1] 0.93571604 0.44240485 0.45451903 -0.96207396 -1.13246522
*

*>>> -0.60032698
*

*>>> [7] -1.77506105 -0.09171419 -0.23262573
*

*>>>
*

*>>> $b
*

*>>> [1] -1.53673409 0.22775629 -0.81788209 -1.16966436 0.73558677
*

*>>> -0.68744178
*

*>>> [7] 0.08639287 1.86770869 -2.92992167
*

*>>>
*

*>>> $c
*

*>>> [1] 1.3657783 0.3987121 -1.4443504 1.0925682 0.8126830
*

*>>> -2.2357115 3.0095462
*

*>>> [8] -1.2282303 -3.5694000
*

*>>>
*

*>>> $d
*

*>>> [1] -0.32235542 -0.61435416 0.84563953 -0.09329507 -1.01613149
*

*>>> 1.74841922
*

*>>> [7] -1.41259217 -0.78038989 -0.78038989
*

*>>>
*

*>>> WHen I look at ?spline there is even an example of "manually" using
*

*>>> these coefficeients:
*

*>>>
*

*>>> ## Manual spline evaluation --- demo the coefficients :
*

*>>> .x <- get("ux", envir = environment(f))
*

*>>> u <- seq(3,6, by = 0.25)
*

*>>> (ii <- findInterval(u, .x))
*

*>>> dx <- u - .x[ii]
*

*>>> f.u <- with(splinecoef,
*

*>>> y[ii] + dx*(b[ii] + dx*(c[ii] + dx* d[ii])))
*

*>>> stopifnot(all.equal(f(u), f.u))
*

*>>>
*

*>>>
*

*>>> For the smooth.spline as
*

*>>>
*

*>>> spl <- smooth.spline(x,y)
*

*>>>
*

*>>> I can also look at the coefficients:
*

*>>>
*

*>>> spl$fit
*

*>>> $knot
*

*>>> [1] 0.000 0.000 0.000 0.000 0.125 0.250 0.375 0.500 0.625 0.750
*

*>>> 0.875 1.000
*

*>>> [13] 1.000 1.000 1.000
*

*>>>
*

*>>> $nk
*

*>>> [1] 11
*

*>>>
*

*>>> $min
*

*>>> [1] 1
*

*>>>
*

*>>> $range
*

*>>> [1] 8
*

*>>>
*

*>>> $coef
*

*>>> [1] 0.90345898 0.73823276 0.40777431 -0.08046715 -0.54625461
*

*>>> -0.85205147
*

*>>> [7] -0.96233408 -0.91373830 -0.66529714 -0.47674774 -0.38246971
*

*>>>
*

*>>> attr(,"class")
*

*>>> [1] "smooth.spline.fit"
*

*>>>
*

*>>> But there isn't an example on how to "manual" use these
*

*>>> coefficients. This is what I was asking about. Once I hae the
*

*>>> coefficients how do I "manually" interpolate using the coefficients
*

*>>> given and x.
*

*>>>
*

*>>> Thank you.
*

*>>>
*

*>>> Kevin
*

*>>>
*

*>>>
*

*>>> ---- Spencer Graves <spencer.graves_at_pdf.com> wrote:
*

*>>>> PLEASE do read the posting guide
*

*>>>> http://www.R-project.org/posting-guide.html and provide commented,
*

*>>>> minimal, self-contained, reproducible code.
*

*>>>>
*

*>>>> I do NOT know how to do what you want, but with a
*

*>>>> self-contained example, I suspect many people on this list --
*

*>>>> probably including me -- could easily solve the problem. Without
*

*>>>> such an example, there is a high probability that any answer might
*

*>>>> (a) not respond to your need, and (b) take more time to develop,
*

*>>>> just because we don't know enough of what you are asking.
*

*>>>> Spencer
*

*>>>>
*

*>>>> rkevinburton_at_charter.net wrote:
*

*>>>>
*

*>>>>> Like I indicated. I understand the coefficients in a B-spline
*

*>>>>> context. If I use the the 'spline' or 'splinefun' I can get the
*

*>>>>> coefficients and they are grouped as 'a', 'b', 'c', and 'd'
*

*>>>>> coefficients. But the coefficients for smooth.spline is just an
*

*>>>>> array. I basically want to take these coefficients and outside of
*

*>>>>> 'R' use them to form an interpolation. In other words I want 'R'
*

*>>>>> to do the hard work and then export the results so they can be
*

*>>>>> used else where.
*

*>>>>>
*

*>>>>> Thank you.
*

*>>>>>
*

*>>>>> Kevin
*

*>>>>>
*

*>>>> Spencer Graves wrote:
*

*>>>>
*

*>>>>> I believe that a short answer to your question is that the
*

*>>>>> "smooth" is a linear combination of B-spline basis functions, and
*

*>>>>> the coefficients are the weights assigned to the different
*

*>>>>> B-splines in that basis.
*

*>>>>> Before offering a much longer answer, I would want to know
*

*>>>>> what problem you are trying to solve and why you want to know.
*

*>>>>> For a brief description of B-splines, see
*

*>>>>> "http://en.wikipedia.org/wiki/B-spline". For a slightly longer
*

*>>>>> commentary on them I suggest the "scripts\ch01.R" in the
*

*>>>>> DierckxSpline package: That script computes and displays some
*

*>>>>> B-splines using "splineDesign", "spline.des" in the 'splines'
*

*>>>>> package plus comparable functions in the 'fda' package. For more
*

*>>>>> info on this, I found the first chapter of Paul Dierckx (1993)
*

*>>>>> Curve and Surface Fitting with Splines (Oxford U. Pr.). Beyond
*

*>>>>> that, I've learned a lot from the 'fda' package and the two
*

*>>>>> companion volumes by Ramsay and Silverman (2006) Functional Data
*

*>>>>> Analysis, 2nd ed. and (2002) Applied Functional Data Analysis
*

*>>>>> (both Springer).
*

*>>>>> If you'd like more help from this listserve, PLEASE do read
*

*>>>>> the posting guide http://www.R-project.org/posting-guide.html and
*

*>>>>> provide commented, minimal, self-contained, reproducible code.
*

*>>>>> Hope this helps. Spencer Graves
*

*>>>>>
*

*>>>>> rkevinburton_at_charter.net wrote:
*

*>>>>>
*

*>>>>>> I like what smooth.spline does but I am unclear on the output. I
*

*>>>>>> can see from the documentation that there are fit.coef but I am
*

*>>>>>> unclear what those coeficients are applied to.With spline I
*

*>>>>>> understand the "noraml" coefficients applied to a cubic
*

*>>>>>> polynomial. But these coefficients I am not sure how to
*

*>>>>>> interpret. If I had a description of the algorithm maybe I could
*

*>>>>>> figure it out but as it is I have this question. Any help?
*

*>>>>>>
*

*>>>>>> Kevin
*

*>>>>>>
*

*>>>>>> ______________________________________________
*

*>>>>>> R-help_at_r-project.org mailing list
*

*>>>>>> https://stat.ethz.ch/mailman/listinfo/r-help
*

*>>>>>> PLEASE do read the posting guide
*

*>>>>>> http://www.R-project.org/posting-guide.html
*

*>>>>>> and provide commented, minimal, self-contained, reproducible code.
*

*>>>>>>
*

*>>> ______________________________________________
*

*>>> R-help_at_r-project.org mailing list
*

*>>> https://stat.ethz.ch/mailman/listinfo/r-help
*

*>>> PLEASE do read the posting guide
*

*>>> http://www.R-project.org/posting-guide.html
*

*>>> and provide commented, minimal, self-contained, reproducible code.
*

*>>>
*

*>>
*

*>> ______________________________________________
*

*>> R-help_at_r-project.org mailing list
*

*>> https://stat.ethz.ch/mailman/listinfo/r-help
*

*>> PLEASE do read the posting guide
*

*>> http://www.R-project.org/posting-guide.html
*

*>> and provide commented, minimal, self-contained, reproducible code.
*

*>
*

R-help_at_r-project.org mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Sun 20 Jul 2008 - 15:55:04 GMT

Date: Sun, 20 Jul 2008 08:51:31 -0700

<in line>

Duncan Murdoch wrote:

> On 20/07/2008 11:11 AM, Spencer Graves wrote:

*>> Are you aware that there are many different kinds of splines?
**>> With "spline" and "splinefun", you can use method = "fmm" (Forsyth,
**>> Malcolm and Moler), "natural", or "periodic". I'm not familiar with
**>> "fmm", but it seems to be adequately explained by the "Manual spline
**>> evaluation" you quoted from the documentation.
**>> Natural splines are perhaps the simplest: I(x-x0)*(x-x0)^j,
**>> where x0 is a knot, and I(z) = 1 if z>0 and 0 otherwise.
**>
**> That's not what R means by "natural spline" in this context. Here it
**> means that the function becomes linear outside the range of the knots.
**>
**> I would call the I(x-x0)*(x-x0)^j splines the "truncated power basis"
**> for polynomial splines; B-splines are a different basis for the same
**> set of splines (assuming the knots and degrees match). Natural
**> splines are a subspace of these (since linear functions are a subspace
**> of polynomials). I don't know of a simple basis for them.
*

Thanks for the correction. I erred by writing this from memory. Dierckx (1993, p. 4) says, "A natural spline function is a spline of odd degree k = 2*m-1 (m>=2) which satisfies the additional constraints

(D^(m+j))s(a) = (D^(m+j))s(b) = 0, j = 0, 1, ..., m-2.

He further (p. 5) defines the "truncated power functions", which is what I mistakenly called "natural splines".

Thanks again for the correction. Spencer

> Duncan Murdoch

R-help_at_r-project.org mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Sun 20 Jul 2008 - 15:55:04 GMT

Archive maintained by Robert King, hosted by
the discipline of
statistics at the
University of Newcastle,
Australia.

Archive generated by hypermail 2.2.0, at Sun 20 Jul 2008 - 16:31:53 GMT.

*
Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help.
Please read the posting
guide before posting to the list.
*