SciPy - Interpolate

SciPy - Interpolation (scipy.interpolate)

Interpolation is a fundamental technique used to estimate unknown values that fall between known data points. It’s like ‘drawing a line’ or ‘fitting a curve’ through your existing data to predict values at intermediate points. This chapter explores how SciPy’s scipy.interpolate module facilitates various interpolation methods.

What is Interpolation?

Imagine you have a set of measurements, for example, temperature readings at specific times. Interpolation helps you estimate the temperature at a time for which you don’t have a direct measurement, assuming this time falls within the range of your observed times. It’s a valuable tool in statistics, science, engineering, and finance for tasks like data smoothing, resampling, and filling in missing data.

Let’s generate some sample data to demonstrate interpolation:

import numpy as np
from scipy import interpolate
import matplotlib.pyplot as plt # For visualization

# Generate some sample data points (x, y)
# np.linspace(start, stop, num_points)
x_known = np.linspace(0, 10, 11) # 11 points from 0 to 10
y_known = np.cos(-x_known**2 / 9.0) + 0.5 * np.sin(x_known / 2.0)

print("Known x values:", x_known)
print("Known y values:", y_known)

# Visualize the known data points
# plt.plot(x_known, y_known, 'o', label='Known Data Points')
# plt.xlabel("x")
# plt.ylabel("y")
# plt.title("Original Data Points")
# plt.legend()
# plt.grid(True)
# plt.show()

This code creates x_known and y_known arrays. If plotted, these points would appear scattered, following a complex curve. For example, the output values would be:

Known x values: [ 0.  1.  2.  3.  4.  5.  6.  7.  8.  9. 10.]
Known y values: [ 1.          1.22690959  1.12199856  0.62092957 -0.13149013 -0.77131049
 -0.9996715  -0.75164093 -0.11188687  0.66798572  1.1840709 ]

If you were to plot these x_known and y_known points, you would see a scatter plot representing the function $y = \cos(-x^2/9) + 0.5 \sin(x/2)$ at discrete x values. The task of interpolation is to estimate y values for x’s not in x_known.

1-D Interpolation with interp1d

The interp1d class in scipy.interpolate is a convenient way to create an interpolation function from fixed data points. This function can then be evaluated at any point within the domain of the original data.

Let’s create interpolation functions using linear and cubic methods:

# Create interpolation functions
# 'linear' interpolation: connects known points with straight lines
interp_linear = interpolate.interp1d(x_known, y_known, kind='linear')

# 'cubic' interpolation: fits a cubic spline through the known points for a smoother curve
interp_cubic = interpolate.interp1d(x_known, y_known, kind='cubic')

# Generate new x values where we want to interpolate y values
x_new = np.linspace(0, 10, 101) # 101 points for a smoother curve

# Interpolate y values at x_new using the created functions
y_linear_interp = interp_linear(x_new)
y_cubic_interp = interp_cubic(x_new)

# Visualize the results (requires matplotlib)
# plt.figure(figsize=(10, 6))
# plt.plot(x_known, y_known, 'o', label='Known Data Points', markersize=8)
# plt.plot(x_new, y_linear_interp, '-', label='Linear Interpolation')
# plt.plot(x_new, y_cubic_interp, '--', label='Cubic Interpolation')
# plt.xlabel("x")
# plt.ylabel("y")
# plt.title("1-D Interpolation Examples")
# plt.legend()
# plt.grid(True)
# plt.show()

interp1d takes the known x and y values and a kind argument specifying the interpolation method. Common kinds include: ‘linear’, ‘nearest’, ‘zero’, ‘slinear’ (spline linear), ‘quadratic’, ‘cubic’. The returned objects (interp_linear, interp_cubic) are callable functions.

If you were to plot this, you’d see the original data points as circles. The ‘Linear Interpolation’ would be a series of straight lines connecting these points. The ‘Cubic Interpolation’ would be a much smoother curve passing through all the original points, generally providing a more accurate representation of an underlying smooth function.

The textual output for the first few interpolated y-values would look like:

# First 5 y-values from linear interpolation for x_new[0] to x_new[4]
# print("Linear interp y_new[:5]:", y_linear_interp[:5])
# First 5 y-values from cubic interpolation for x_new[0] to x_new[4]
# print("Cubic interp y_new[:5]:", y_cubic_interp[:5])

# Example Output Snippet:
# Linear interp y_new[:5]: [1.         1.02269096 1.04538192 1.06807287 1.09076383]
# Cubic interp y_new[:5]:  [1.         1.03396658 1.06506318 1.09264928 1.1160844 ]

Spline Interpolation

Splines are special functions defined piecewise by polynomials. They are widely used for fitting smooth curves through data points. Historically, drafters used flexible strips of wood (mechanical splines), bent to pass through specified points (knots), to draw smooth curves. Digital splines achieve a similar effect mathematically.

Univariate Spline with UnivariateSpline

The UnivariateSpline class in scipy.interpolate fits a smoothing spline to a given set of 1-D data points. Unlike interp1d which (for kinds like ‘cubic’) passes exactly through all points, UnivariateSpline can create a smoothed approximation, controlled by a smoothing factor s.

Key parameters for UnivariateSpline(x, y, w=None, bbox=[None, None], k=3, s=None, ext=0, check_finite=False):

• x, y: The data points to fit.
• w: Weights for spline fitting. If None, all points have equal weight.
• k: Degree of the smoothing spline (e.g., k=1 for linear, k=3 for cubic). Default is 3.
• s: Smoothing factor. If s=0, the spline will interpolate through all data points (if possible). If s is larger, the spline will be smoother but may not pass through all points. If None, s is chosen automatically based on the data.
• ext: Controls extrapolation behavior for points outside the original data range.

Let’s generate noisy data and fit a spline:

import numpy as np
from scipy.interpolate import UnivariateSpline
import matplotlib.pyplot as plt # For visualization

# Generate noisy data
np.random.seed(42)
x_noisy = np.linspace(-3, 3, 60)
y_noisy = np.exp(-x_noisy**2) + 0.15 * np.random.randn(len(x_noisy))

# Fit a UnivariateSpline (default k=3, s might be auto-chosen or 0 if data allows)
# s=0 forces interpolation through all points if no weights are given and task is possible
spline_interp = UnivariateSpline(x_noisy, y_noisy, s=0) # s=0 attempts to interpolate

# Fit a spline with some smoothing (s > 0)
# A larger 's' means more smoothing.
spline_smooth = UnivariateSpline(x_noisy, y_noisy, s=0.5) # s=0.5 for smoothing

# Generate points for plotting the spline curves
x_plot = np.linspace(-3, 3, 200)
y_spline_interp_plot = spline_interp(x_plot)
y_spline_smooth_plot = spline_smooth(x_plot)

# Visualize (requires matplotlib)
# plt.figure(figsize=(10, 6))
# plt.plot(x_noisy, y_noisy, 'o', label='Noisy Data Points', alpha=0.6)
# plt.plot(x_plot, y_spline_interp_plot, '-', label='Spline (s=0, interpolating)')
# plt.plot(x_plot, y_spline_smooth_plot, '--', label='Spline (s=0.5, smoothing)')
# plt.xlabel("x")
# plt.ylabel("y")
# plt.title("Univariate Spline Interpolation and Smoothing")
# plt.legend()
# plt.grid(True)
# plt.show()

If you plotted this: The ‘Noisy Data Points’ would be scattered. The ‘Spline (s=0, interpolating)’ curve would try to pass through all these noisy points, potentially resulting in a wobbly line. The ‘Spline (s=0.5, smoothing)’ curve would be smoother, capturing the general trend of the data rather than fitting every noise fluctuation.

You can adjust the smoothing factor s using spline_object.set_smoothing_factor(new_s_value) after creating the spline object if you want to experiment with different levels of smoothing.

Interpolation is critical in many domains: signal processing (resampling audio or images), financial modeling (estimating asset prices between trades), weather forecasting (estimating conditions at unmonitored locations), and computer graphics (creating smooth paths or surfaces).

scipy.interpolate also offers more advanced interpolation techniques, including multivariate interpolation (for data in 2D, 3D, or higher dimensions) using functions like griddata, LinearNDInterpolator, CloughTocher2DInterpolator, and Radial Basis Function (Rbf) interpolation.