Least-square spline interpolation forcing interpolant to pass through specific points

0
votes

I am having issues in implementing some less-than-usual interpolation problem. I have some (x,y) data points scattered along some curve which a priori I don't know, and I want to reconstruct this curve at my best, interpolating my point with min square error. I thought of using scipy.interpolate.splrep for this purpose (but maybe there are better options you would advise to use). The additional difficulty in my case, is that I want to constrain the spline curve to pass through some specific points of my original data. I assume that playing with knots and weights could make the trick, but I don't know how to do so (I am procrastinating avoidance of spline interpolation theory besides basic fitting procedures). Also, for some undisclosed reasons, when I try to setup knots in my splrep I get the same error of this post, which keeps complicating things. The following is my sample code:

from __future__ import division
import numpy as np
import scipy.interpolate as spi
import matplotlib.pylab as plt

# Some surrogate sample data
f = lambda x : x**2 - x/2.
x = np.arange(0.,20.,0.1)
y = f(4*(x + np.random.normal(size=np.size(x))))

# I want to use spline interpolation with least-square fitting criterion, making sure though that the spline starts
# from the origin (or in general passes through a precise point of my dataset). 
# In my case for example I would like the spline to originate from the point in x=0. So I attempted to include as first knot x=0... 
# but it won't work, nor I am sure this is the right procedure...

fy = spi.splrep(x,y)
fy = spi.splrep(x,y,t=fy[0])
yy = spi.splev(x,fy)

plt.plot(x,y,'-',x,yy,'--')
plt.show()

which despite the fact I am even passing knots computed from a first call of splrep, it will give me:

  File "/usr/lib64/python2.7/site-packages/scipy/interpolate/fitpack.py", line 289, in splrep
    res = _impl.splrep(x, y, w, xb, xe, k, task, s, t, full_output, per, quiet)
  File "/usr/lib64/python2.7/site-packages/scipy/interpolate/_fitpack_impl.py", line 515, in splrep
    raise _iermess[ier][1](_iermess[ier][0])
ValueError: Error on input data
2
pythonscipyconstraintsinterpolationspline

2 Answers

1
votes

You use the weights argument of splrep: can give these points you need fixed very large weights. This is a workaround for sure, so keep an eye on the fit quality and stability.

0
votes

Setting high weights for specific points is indeed a working solution as suggested by @ev-br. In addition, because there is no direct way to match derivatives at the extrema of the curve, the same rationale can be applied in this case as well. Say you want the derivative in y[0] and y[-1] match the derivative of your data points, then you add large weights also for y[1] and y[-2], i.e.

weights = np.ones(len(x))
weights[[0,-1]] = 100  # Promote spline interpolant through first and last point
weights[[1,-2]] = 50   # Make spline interpolant derivative tend to derivatives at first/last point
fy = spi.splrep(x,y,w=weights,s=0.1)
yy = spi.splev(x,fy)