[R] smooth.spline
Spencer Graves
spencer.graves at pdf.com
Sun Jul 20 17:51:31 CEST 2008
<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
>
>>
>> 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