5

csaps() in matlab does a cubic spline according to a particular definition of the smoothing parameter p. Here is some matlab code and its result:

     % x variable
    age = 75:99  

    % y variable
    diffs = [-39   -2 -167  -21  -13   32  -37 -132 -143  -91  -93  -88  -62 -112  -95  -28  -90  -40  -27  -23  -28  -11   -8   -6    1]

    % 0.0005 is the parameter p, and the later specification of 
    % age are the desired x for prediction
    csaps(age,diffs,0.0005,age)
    % result (column headers removed):
     -63.4604  -64.0474  -64.6171  -65.1397  -65.6111  -66.0165  -66.3114  
     -66.4123  -66.2229  -65.6726  -64.7244  -63.3582  -61.5676  -59.3568  
     -56.7364  -53.7382  -50.4086  -46.7922  -42.9439  -38.9183  -34.7629  
     -30.5180  -26.2186  -21.8912  -17.5532

I'd like to get the same result in R. I've tried base::smooth.spline(), but the smoothing parameter spar is specified in a different way that I can't seem to relate to matlab's p (can you?). The closest result I've been able to get has been with the smooth.Pspline() function of the pspline package. Here is some code to get things rolling in R:

age <- 75:99
diffs <- c(-39L, -2L, -167L, -21L, -13L, 32L, -37L, -132L, -143L, -91L, 
-93L, -88L, -62L, -112L, -95L, -28L, -90L, -40L, -27L, -23L, 
-28L, -11L, -8L, -6L, 1L)
predict(pspline::smooth.Pspline(
                           x = age,
                           y = diffs, 
                           norder = 2, 
                           method = 1,
                           spar = 1 / 0.0005     # p given in MP and matlab as 0.0005
                         ),age)
# which gives something close, but not exactly the same:
 [1] -63.46487 -64.05103 -64.61978 -65.14158 -65.61214 -66.01662 -66.31079
 [8] -66.41092 -66.22081 -65.67009 -64.72153 -63.35514 -61.56447 -59.35372
[15] -56.73367 -53.73584 -50.40680 -46.79098 -42.94333 -38.91850 -34.76393
[22] -30.51985 -26.22131 -21.89474 -17.55757

The csaps() help page is here

smooth.spline() help can be found here (code not given because I think maybe the relation between spar and p is pretty hairy, so maybe not worth going down this path)

pspline::smooth.Pspline() help is here

This other person's quest from 2008 appears to have gone unanswered, making me feel like this guy.

R is chock full of spline doings, so if the saavy amongst ye can point me to the one that does the same thing as matlab's csaps() (or a trick along those lines) I'd be most appreciative.

[EDIT 19-8-2013] spar needs to be specified as (1-p)/p (rather than 1/p) and then results will agree to as far as numerical precision can take you. See answer below.

tim riffe
  • 5,651
  • 1
  • 26
  • 40
  • You already have a method that is in agreement to 5 decimal places (on the ratios). I fail to see the value in trying to get exact agreement to an approximation to an undefined "something". – IRTFM Aug 17 '13 at 00:51
  • @DWin This is a small sprocket in a big matlab -> R conversion project of some legacy code. Lots of stuff downstream depends on this, and the function is used more than once in the original code. It's important to us because downstream estimates need to be exact, as in *exact* – tim riffe Aug 17 '13 at 00:57
  • Exact? You just drew two splines through that set of scattered points. How can you assert that either one of them is "exact"? – IRTFM Aug 17 '13 at 01:20
  • @DWin yeah, it's like fingerpainting, I get it. I don't claim that either exactly represents the data. What I need is to reproduce a result *exactly*, not find the truth exactly. Task not whim. The question is valid. If it's not implemented in R then I'll make do, but I'll still fish for an answer if I may. – tim riffe Aug 17 '13 at 02:05
  • You should examine how sensitive the R result is to variations in "spar". – IRTFM Aug 17 '13 at 02:17
  • If you really must produce identical values then `edit csaps` and look at the code. – horchler Aug 17 '13 at 15:48
  • thank you DWin and horchler. Got the answer-- I needed to specify spar differently – tim riffe Aug 19 '13 at 15:24

2 Answers2

6

My colleague found the answer: One converts matlab's p to pspline::smooth.Pspline()'s spar not as 1/p, but as (1-p)/p, and then results will agree out to whatever the degree of numerical precision is:

c(predict(pspline::smooth.Pspline(
                          x = age,
                          y = diffs, 
                          norder = 2, 
                          method = 1,
                          spar = (1-0.0005) / 0.0005     # p given in MP and matlab as 
                         ),age))
 [1] -63.46035 -64.04741 -64.61705 -65.13972 -65.61114 -66.01646 -66.31144
 [8] -66.41232 -66.22285 -65.67263 -64.72443 -63.35823 -61.56761 -59.35675
[15] -56.73643 -53.73821 -50.40864 -46.79221 -42.94387 -38.91828 -34.76291
[22] -30.51801 -26.21863 -21.89122 -17.55320
tim riffe
  • 5,651
  • 1
  • 26
  • 40
  • thanks, really helpfull 10Y later (there are indeed two main cubic smoothing spline algo and it was difficuly to identify the correct one in R : https://stackoverflow.com/questions/75221215/is-there-an-r-spline-that-matches-the-schoenberg-algorithm/75235153#75235153 – Nielsou Akbrg Jan 25 '23 at 14:07
-1

Here is what I found in p. 16 of MATLAB/R Reference by David Hiebeler. [I don't use Matlab,however].

Fit natural cubic spline(S′′(x) = 0 at both endpoints) to points (xi, yi)whose coordinates are in vectors x and y; evaluate at points whose x coordinates are in vector xx, storing corresponding y’s in yy

Matlab: 

pp=csape(x,y,’variational’);
yy=ppval(pp,xx) but note that
csape is in Matlab’s Spline
Toolbox

R

tmp=spline(x,y,method=’natural’,
xout=xx); yy=tmp$y
Metrics
  • 15,172
  • 7
  • 54
  • 83
  • Thanks for trying, but no dice. Looking to reproduce `csaps()`. I adapted your (Hiebeler's) R code to my problem and it doesn't do the trick. – tim riffe Aug 17 '13 at 00:33