Behavior of scipy's splrep

Question

I have a set of data points and would like to approximate them with a spline function. I used two different functions:

1. splrep from scipy
2. and a cubic spline function that I found here.

The results look like this.

The code is as follows:

from matplotlib.pyplot import *
from numpy import *
from scipy import interpolate
#----------------------------------------------
s = arange(257)/256.0
z = s[::-1]
b = transpose(array((z*z*z,
                 3*z*z*s, 
                 3*z*s*s,
                 s*s*s)))
def cubicspline(c,t): 
return dot(b[t],c)
#----------------------------------------------

A = array([
   [ -126.041   ,  246.867004],
   [ -113.745003,   92.083   ],
   [  208.518997, -183.796997],
   [  278.859009, -190.552994]])

a1 = A[:,0]
a2 = A[:,1] 
cs = reshape(A, (-1, 4, 2))
X = []
Y = []
#spline with cubicspline()
for (x,y) in [cubicspline(c,16*t) for c in cs for t in arange(17)]:
X.append(x)
Y.append(y)

# spline with splrep
tck = interpolate.splrep( a1, a2)

xnew = np.arange( min(a1), max(a1), 5)
ynew = interpolate.splev(xnew, tck)
plot(a1, a2, "--ob", ms = 9,  label = "points")
plot(X, Y, "r", lw=2, label = "cubicspline")
plot(xnew, ynew, "g", lw=2, label = "splrep")
legend(); savefig("image.png"); show()

As you may see the results of splrep are far from being satisfying. Can someone please explain this behavior and how to get reasonable approximation from splrep?

Answer 1

You need to define what you mean by "satisfying". Clearly, your cubic spline is not interpolating through the points, whereas the splrep result does (and is perfectly satisfactory in that sense). Note also that your 'cubicspline' is actually just a single polynomial rather than a spline (which are polynomials with breakpoints).

You need to explicitly tell splrep that the spline doesn't need to go through the points --- pass in a nonzero s smoothing parameter. How to choose this properly, see this question: scipy.interpolate.UnivariateSpline not smoothing regardless of parameters