[R] smooth.spline
Duncan Murdoch
murdoch at stats.uwo.ca
Sun Jul 20 17:34:54 CEST 2008
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.
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.
More information about the R-help
mailing list