Re: [R] smooth.spline

Sun, 20 Jul 2008 08:54:01 -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




      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

[EMAIL PROTECTED] 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 <[EMAIL PROTECTED]> 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

[EMAIL PROTECTED] 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

[EMAIL PROTECTED] 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@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@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@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@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.






















					Reply via email to