Interpolation using Scipy

Scipy is a very versatile and extremely efficient python library written on top of Numpy. SciPy works great for all kinds of scientific programming projects (science, mathematics, and engineering). It offers efficient numerical routines such as numerical optimization, integration, and others in submodules. The extensive documentation makes working with this library really easy.

Interpolation & why do we need this

Interpolation is trying to figure out what the result would be for a new x value for a given set of x & y values. Here is an example of what an interpolation is

#importing matplotlib and setting some variables
import matplotlib.pyplot as plt
%matplotlib inline
plt.rcParams['figure.figsize']=[12,8]
plt.style.use('dark_background')
#importing np, scipy and interpolate module from scipy
import scipy
import numpy as np
from scipy import interpolate
#generating some x values and using exp function to get y values
x = np.arange(5,20)
y = np.exp(x/3.0)
#using interpolate.interp1d to get some "interpolated" values
f = interpolate.interp1d(x,y)
x1 = np.arange(6,18)
y1 = f(x1)
#plotting these together to show the interpolation
plt.plot(x,y,'o',x1,y1,'r--')
plt.show()

Obviously the scipy library is very good and correctly maps to all the points and seems to getting the interpolation right too.

Extrapolation: Out of sample interpolation
The bigger question is how would scipy interpolation (rather extrapolation) would perform on data that it has not seen and is not within the given range of x values.

#generating x2,y2. 
#To see how extrapolation works we need x2 should be more than x
x2 = np.arange(5,30)
y2 = np.exp(x2/3.0)
#Here we use UnivariateSpline function to get extrapolated values
f2 = interpolate.UnivariateSpline(x,y)

x1 = np.arange(6,30)
y1 = f2(x1)
#plotting all x,y, x1,y1 & x2,y2 values
plt.plot(x1,y1,'y--',x2,y2,'wo',x,y,'ro')
plt.show()

As you can see from above, we see that all values within the x range have correct y values ie x1,y1 are right on target. But once we go out of range ie x>20 we can see that we are extrapolating here and the values slowly diverge. so need to be careful with the values we get when we are extrapolating.

Interpolation for higher order functions
Here we will look at how we can interpolate better for higher order functions. Here is an example graph of interpolation using a ‘linear’ model.

We can get better results using cubic spline interpolation. This is done in two steps, first we define the splines using the scipy interpolation modules “splrep” function to create the spline representations. We then use the “splev” function to evaluate the new y’s based on x’s & the tuple returned by splrep.

from scipy import interpolate
x = np.arange(0,360,5)
y = np.sin(np.deg2rad(x))**2-np.cos(np.deg2rad(x))**2
x2 = np.arange(0,360,10)
#Here we create the splines tuple using splrep
splines = interpolate.splrep(x,y)
#Here we get y2 from x2 & the splrep tuple y2=interpolate.splev(x2,splines)
plt.plot(x,y,'o',x2,y2,'--x')

Finally lets see how we can do this for multivariate functions.

We create a simple surface using the function
z = sin(x) + cos(y)

#basic matplotlib imports
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D

#helper function to draw 3d-surfaces
def plot_3dsurface(x, y, z):
    fig = plt.figure(figsize=(10,8))
    ax = Axes3D(fig)
    surface_plot = ax.plot_surface(x, y, z, 
                   cmap=cm.Spectral,linewidth=0.1)
    plt.show()
x = np.arange(-3, 3, 1) 
y = np.arange(-3, 3, 1) 
xm, ym = np.meshgrid(x, y)
zm = np.sin(xm) + np.cos(ym)
plot_3dsurface(xm,ym,zm)

Now the same surface but with ‘linear’ interpolation.

5#Using 'linear' interpolation
f = interpolate.interp2d(xm, ym, zm, kind='linear')
x2 = np.arange(-3, 3, 0.1)
y2 = np.arange(-3, 3, 0.1)
x2m, y2m = np.meshgrid(x2, y2)
z2 = f(x2,y2)
plot_3dsurface(x2m,y2m,z2)

Here is a much smoother interpolation using cubic spline interpolation

#Using 'cubic' interpolation
f2 = interpolate.interp2d(xm, ym, zm, kind='cubic')
x3 = np.arange(-3, 3, 0.1)
y3 = np.arange(-3, 3, 0.1)
x3m, y3m = np.meshgrid(x3, y3)
z3 = f2(x3,y3)
plot_3dsurface(x3m,y3m,z3)

Conclusions
1. You can use scipy for both interpolation and extrapolation.
2. Extrapolation needs to be with in reasonable bounds.
3. Finally when using higher order functions use cubic spline interpolation.