Data Fitting with Polynomials and Splines

The first two examples are from the tutorial for the SciPy interpolation module scipy.interpolate.

Interpolation of 1D data

In [2]:
%matplotlib inline
import numpy as np
from scipy.interpolate import interp1d
In [37]:
x = np.linspace(0, 10, 10)
y = np.cos(-x**2/8.0)
f = interp1d(x, y)
f2 = interp1d(x, y, kind='cubic')
In [38]:
xnew = np.linspace(0, 10, 40)
import matplotlib.pyplot as plt
plt.plot(x,y,'o',xnew,f(xnew),'-', xnew, f2(xnew),'--')
plt.legend(['data', 'linear', 'cubic'], loc='best')
plt.show()

Interpolation of multivariate data

In [3]:
def func(x, y):
    return x*(1-x)*np.cos(4*np.pi*x) * np.sin(4*np.pi*y**2)**2
In [21]:
grid_y, grid_x = np.meshgrid(np.linspace(0,1,200), np.linspace(0,1,100))
In [22]:
points = np.random.rand(1000, 2)
values = func(points[:,0], points[:,1])
In [23]:
from scipy.interpolate import griddata
grid_z0 = griddata(points, values, (grid_x, grid_y), method='nearest')
grid_z1 = griddata(points, values, (grid_x, grid_y), method='linear')
grid_z2 = griddata(points, values, (grid_x, grid_y), method='cubic')
In [26]:
import matplotlib.pyplot as plt
plt.subplot(221)
plt.imshow(func(grid_x, grid_y).T, extent=(0,1,0,1), origin='lower')
plt.plot(points[:,0], points[:,1], 'k.', ms=1)
plt.title('Original')
plt.subplot(222)
plt.imshow(grid_z0.T, extent=(0,1,0,1), origin='lower')
plt.title('Nearest')
plt.subplot(223)
plt.imshow(grid_z1.T, extent=(0,1,0,1), origin='lower')
plt.title('Linear')
plt.subplot(224)
plt.imshow(grid_z2.T, extent=(0,1,0,1), origin='lower')
plt.title('Cubic')
plt.gcf().set_size_inches(6, 6)
plt.show()

Splines

In [28]:
from math import pi
import numpy as np
import matplotlib.pyplot as plt
from scipy import interpolate
In [29]:
x = np.arange(0, 2*pi+pi/4, 2*pi/8)
y = np.sin(x)
s = interpolate.InterpolatedUnivariateSpline(x, y)
xnew = np.arange(0, 2*pi, pi/50)
ynew = s(xnew)
In [31]:
plt.figure()
plt.plot(x, y, 'x', xnew, ynew, xnew, np.sin(xnew), x, y, 'b')
plt.legend(['Linear', 'spline', 'True'])
plt.axis([-0.05, 6.33, -1.05, 1.05])
plt.title('InterpolatedUnivariateSpline')
plt.show()
In [32]:
ynew = s.derivative()(xnew)
plt.figure()
plt.plot(xnew, ynew)
plt.show()
In [49]:
print s.integral(0, pi)
2.00182213189
In [50]:
print s.roots()
[ 3.14159265]

Lets compare NumPy's polynomial interpolation with spline interpolation.

In [51]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import InterpolatedUnivariateSpline as Spline

x = np.array([1., 2., 3., 4., 5., 6.])
y = np.array([2., 6., 4., 5., 0., 5.])
p = np.poly1d(np.polyfit(x, y, 5))  # polyfit gives you the points 
                                    # and poly1d turns it into a 
                                    # convenient object for evaluation.
spl = Spline(x, y)

t = np.linspace(1, 6, 50)

plt.figure()
plt.plot(t, p(t), 'g-')
plt.plot(t, spl(t), 'r--')
plt.scatter(x, y)
plt.show()