[SciPy-User] Interpolating discrete values with splprep and splev

[josef.pktd at gmail.com]

On Thu, Aug 19, 2010 at 5:08 PM, Matthew Cragun <mcragun at gmail.com> wrote:
> I'm trying to fit a curve to some data points using splprep and splev.  Then
> I would like to take a few discrete values and get the inerpolated values
> back using the spline.
> I am having a problem getting the correct y values out based on the spline
> fit.  In my code below I would like a plot of the saex values against their
> interpolated values, but can't get them out properly.
> saex2 gives: [  37.27079882   55.12775627  106.64397968  142.10925103
>  158.67651813 176.27838508]
> saey  gives:  [  1.76685635e-01  -3.89842080e+00  -5.48621169e+01
>  -1.42948834e+02  -2.03504645e+02  -2.83162159e+02]
> Any thoughts?  Here is my code:
> from numpy import arange, cos, linspace, pi, sin, random, r_, array
> from scipy.interpolate import splprep, splev
> # Original points
> x = array([0,3,6,9,12,15])
> y = array([0.530,0.546,0.586,0.658,0.727,0.801])
> #Interpolated X locations
> saex = array([2.484,3.671,7.063,9.354,10.408,11.515])
> # spline parameters
> s=3.0 # smoothness parameter
> k=3 # spline order
> nest=-1 # estimate of number of knots needed (-1 = maximal)
> # find the knot points
> tckp,u = splprep([x,y],s=s,k=k,nest=nest)
> # evaluate spline
> splinex,spliney= splev(linspace(0,1,400),tckp)
> #evaluate at discrete points
> saex2,saey= splev(saex,tckp,der=0)
> print saex2
> #plot
> import pylab
> #plot the data
> data,=pylab.plot(x,y,'o',label='data')
> #plot the fitted curve
> fit,=pylab.plot(splinex,spliney,'r-',label='fit')
> #plot the discrete locations
> #sae,=pylab.plot(saex2,saey,'x',label='sae points')
> pylab.legend()
> pylab.xlabel('x')
> pylab.ylabel('y')
>
> pylab.savefig('splprep_demo.png')


I'm not quite sure what you are asking for, nor why you are parametric
splines, but if you want the interpolated values for the original x,
then this seems to work up to 4, 5 decimals

>>> res = splev(u,tckp)
>>> res
[array([ -1.40838043e-05,   3.00005731e+00,   5.99991159e+00,
         9.00006221e+00,   1.19999820e+01,   1.50000010e+01]), array([
0.53034055,  0.54429613,  0.58940979,  0.65458787,  0.7287074 ,
        0.80065826])]
>>> x-res[0]
array([  1.40838043e-05,  -5.73105777e-05,   8.84070598e-05,
        -6.22059350e-05,   1.80088155e-05,  -9.83166856e-07])
>>> y-res[1]
array([-0.00034055,  0.00170387, -0.00340979,  0.00341213, -0.0017074 ,
        0.00034174])


>>> u
array([ 0.        ,  0.19996446,  0.39994386,  0.59996306,
0.79997756,  1.        ])
>>> x
array([ 0,  3,  6,  9, 12, 15])
>>> 1.*x/x.max()
array([ 0. ,  0.2,  0.4,  0.6,  0.8,  1. ])



>>> res = splev((x-x.min())*1.0/x.max(),tckp)
>>> res
[array([ -1.40838043e-05,   3.00059043e+00,   6.00075376e+00,
         9.00061633e+00,   1.20003186e+01,   1.50000010e+01]), array([
0.53034055,  0.54430171,  0.5894258 ,  0.65460107,  0.72871581,
        0.80065826])]
>>> x-res[0]
array([  1.40838043e-05,  -5.90429918e-04,  -7.53755099e-04,
        -6.16330116e-04,  -3.18593347e-04,  -9.83166856e-07])
>>> y-res[1]
array([-0.00034055,  0.00169829, -0.0034258 ,  0.00339893, -0.00171581,
        0.00034174])

non-parametric might work easier in this case, but I don't know about
the details.
Are the parametric splines using a symmetric least-squares distance
and the non-parametric a one-sided ?

>>> tckp2 = interpolate.splrep(x,y,s=s,k=k)
>>> res2 = splev(x,tckp2)
>>> y
array([ 0.53 ,  0.546,  0.586,  0.658,  0.727,  0.801])
>>> res - y
array([ 0.00034127, -0.00170635,  0.0034127 , -0.0034127 ,  0.00170635,
       -0.00034127])
>>> res = splev(u,tckp)
>>> res[1] - res2
array([ -7.21620543e-07,   2.47969098e-06,  -2.90360301e-06,
         5.65085112e-07,   1.05423256e-06,  -4.73785097e-07])

Josef


> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
>