by Fedor Iskhakov, ANU
Description: How to approximate functions which are only defined on grid of points. Spline and polynomial interpolation.
between the data points
Extrapolation is computing the approximated function outside of the original data interval
Should be avoided in general
Spline = curve composed of independent pieces
Definition A function $ s(x) $ on $ [a,b] $ is a spline of order $ n $ ( = degree $ n-1 $) iff
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
np.random.seed(2008) # fix random number sequences
x = np.sort(np.random.uniform(-5,10,12)) # sorted random numbers on [-5,10]
xr = np.linspace(-5,10,12) # regular grid on [-5,10]
func=lambda x: np.exp(-x/4)*np.sin(x) + 1/(1+x**2) # function to interpolate
def plot1(ifunc,fdata=(x,func(x)),f=func,color='b',label='',extrapolation=False):
'''helper function to make plots'''
xd = np.linspace(-5,10,1000) # for making continuous lines
plt.figure(num=1, figsize=(10,8))
plt.scatter(fdata[0],fdata[1],color='r') # interpolation data
plt.plot(xd,f(xd),color='grey') # true function
if extrapolation:
xdi = xd
else:
# restriction for interpolation only
xdi=xd[np.logical_and(xd>=fdata[0][0],xd<=fdata[0][-1])]
if ifunc:
plt.plot(xdi,ifunc(xdi),color=color,label=label)
if label:
plt.legend()
elif label:
plt.title(label)
plot1(None,label='True function')
from scipy import interpolate # Interpolation routines
fi = interpolate.interp1d(x,func(x)) # returns the interpolation function
plot1(fi,label='interp1d')
help(interpolate.interp1d)
Help on class interp1d in module scipy.interpolate.interpolate: class interp1d(scipy.interpolate.polyint._Interpolator1D) | interp1d(x, y, kind='linear', axis=-1, copy=True, bounds_error=None, fill_value=nan, assume_sorted=False) | | Interpolate a 1-D function. | | `x` and `y` are arrays of values used to approximate some function f: | ``y = f(x)``. This class returns a function whose call method uses | interpolation to find the value of new points. | | Note that calling `interp1d` with NaNs present in input values results in | undefined behaviour. | | Parameters | ---------- | x : (N,) array_like | A 1-D array of real values. | y : (...,N,...) array_like | A N-D array of real values. The length of `y` along the interpolation | axis must be equal to the length of `x`. | kind : str or int, optional | Specifies the kind of interpolation as a string | ('linear', 'nearest', 'zero', 'slinear', 'quadratic', 'cubic', | 'previous', 'next', where 'zero', 'slinear', 'quadratic' and 'cubic' | refer to a spline interpolation of zeroth, first, second or third | order; 'previous' and 'next' simply return the previous or next value | of the point) or as an integer specifying the order of the spline | interpolator to use. | Default is 'linear'. | axis : int, optional | Specifies the axis of `y` along which to interpolate. | Interpolation defaults to the last axis of `y`. | copy : bool, optional | If True, the class makes internal copies of x and y. | If False, references to `x` and `y` are used. The default is to copy. | bounds_error : bool, optional | If True, a ValueError is raised any time interpolation is attempted on | a value outside of the range of x (where extrapolation is | necessary). If False, out of bounds values are assigned `fill_value`. | By default, an error is raised unless ``fill_value="extrapolate"``. | fill_value : array-like or (array-like, array_like) or "extrapolate", optional | - if a ndarray (or float), this value will be used to fill in for | requested points outside of the data range. If not provided, then | the default is NaN. The array-like must broadcast properly to the | dimensions of the non-interpolation axes. | - If a two-element tuple, then the first element is used as a | fill value for ``x_new < x[0]`` and the second element is used for | ``x_new > x[-1]``. Anything that is not a 2-element tuple (e.g., | list or ndarray, regardless of shape) is taken to be a single | array-like argument meant to be used for both bounds as | ``below, above = fill_value, fill_value``. | | .. versionadded:: 0.17.0 | - If "extrapolate", then points outside the data range will be | extrapolated. | | .. versionadded:: 0.17.0 | assume_sorted : bool, optional | If False, values of `x` can be in any order and they are sorted first. | If True, `x` has to be an array of monotonically increasing values. | | Attributes | ---------- | fill_value | | Methods | ------- | __call__ | | See Also | -------- | splrep, splev | Spline interpolation/smoothing based on FITPACK. | UnivariateSpline : An object-oriented wrapper of the FITPACK routines. | interp2d : 2-D interpolation | | Examples | -------- | >>> import matplotlib.pyplot as plt | >>> from scipy import interpolate | >>> x = np.arange(0, 10) | >>> y = np.exp(-x/3.0) | >>> f = interpolate.interp1d(x, y) | | >>> xnew = np.arange(0, 9, 0.1) | >>> ynew = f(xnew) # use interpolation function returned by `interp1d` | >>> plt.plot(x, y, 'o', xnew, ynew, '-') | >>> plt.show() | | Method resolution order: | interp1d | scipy.interpolate.polyint._Interpolator1D | builtins.object | | Methods defined here: | | __init__(self, x, y, kind='linear', axis=-1, copy=True, bounds_error=None, fill_value=nan, assume_sorted=False) | Initialize a 1D linear interpolation class. | | ---------------------------------------------------------------------- | Data descriptors defined here: | | __dict__ | dictionary for instance variables (if defined) | | __weakref__ | list of weak references to the object (if defined) | | fill_value | The fill value. | | ---------------------------------------------------------------------- | Methods inherited from scipy.interpolate.polyint._Interpolator1D: | | __call__(self, x) | Evaluate the interpolant | | Parameters | ---------- | x : array_like | Points to evaluate the interpolant at. | | Returns | ------- | y : array_like | Interpolated values. Shape is determined by replacing | the interpolation axis in the original array with the shape of x. | | ---------------------------------------------------------------------- | Data descriptors inherited from scipy.interpolate.polyint._Interpolator1D: | | dtype
fi = interpolate.interp1d(x,func(x),kind='linear')
plot1(fi,label='Linear')
for knd, clr in ('previous','m'),('next','b'),('nearest','g'):
fi = interpolate.interp1d(x,func(x),kind=knd)
plot1(fi,label=knd,color=clr)
plt.show()
for knd, clr in ('slinear','m'),('quadratic','b'),('cubic','g'):
fi = interpolate.interp1d(x,func(x),kind=knd)
plot1(fi,color=clr,label=knd)
# Approximation errors
# x = np.sort(np.random.uniform(-5,10,11)) # generate new data
for knd, clr in ('slinear','m'),('quadratic','b'),('cubic','g'):
fi = interpolate.interp1d(x,func(x),kind=knd,bounds_error=False)
xd = np.linspace(-5,10,1000)
erd=np.abs(func(xd)-fi(xd))
plt.plot(xd,erd,color=clr)
print('Max error with %s splines is %1.5e'%(knd,np.nanmax(erd)))
Max error with slinear splines is 1.05142e+00 Max error with quadratic splines is 3.89974e-01 Max error with cubic splines is 4.35822e-01
# Approximation errors for regular grid
for knd, clr in ('slinear','m'),('quadratic','b'),('cubic','g'):
fi = interpolate.interp1d(xr,func(xr),kind=knd,bounds_error=False)
xd = np.linspace(-5,10,1000)
erd=np.abs(func(xd)-fi(xd))
plt.plot(xd,erd,color=clr)
print('Max error with %s splines is %1.5e'%(knd,np.nanmax(erd)))
Max error with slinear splines is 4.63043e-01 Max error with quadratic splines is 3.48546e-01 Max error with cubic splines is 1.89578e-01
How to reduce approximation errors?
In economic models we usually can control all of these
Back to the beginning to explore the idea of replacing original $ f(x) $ with simpler $ g(x) $
Does polynomial $ g(x) $ converge to $ f(x) $ when there are more points?
from numpy.polynomial import polynomial
degree = len(x)-1 # passing through all dots
p = polynomial.polyfit(x,func(x),degree)
fi = lambda x: polynomial.polyval(x,p)
plot1(fi,label='Polynomial of degree %d'%degree,extrapolation=True)
# now with regular grid
degree = len(x)-1 # passing through all dots
p = polynomial.polyfit(xr,func(xr),degree)
fi = lambda x: polynomial.polyval(x,p)
plot1(fi,fdata=(xr,func(xr)),label='Polynomial of degree %d'%degree,extrapolation=True)
# how number of points affect the approximation (with degree=n-1)
for n, clr in (5,'m'),(10,'b'),(15,'g'),(25,'r'):
x2 = np.linspace(-5,10,n)
p = polynomial.polyfit(x2,func(x2),n-1)
fi = lambda x: polynomial.polyval(x,p)
plot1(fi,fdata=(x2,func(x2)),label='%d points'%n,color=clr,extrapolation=True)
plt.show()
# how locations of points affect the approximation (with degree=n-1)
np.random.seed(2025)
n=8
for clr in 'b','g','c':
x2 = np.linspace(-4,9,n) + np.random.uniform(-1,1,n) # perturb points a little
p = polynomial.polyfit(x2,func(x2),n-1)
fi = lambda x: polynomial.polyval(x,p)
plot1(fi,fdata=(x2,func(x2)),label='%d points'%n,color=clr,extrapolation=True)
plt.show()
# how degree of the polynomial affects the approximation
for degree, clr in (7,'b'),(9,'g'),(11,'m'):
p=polynomial.polyfit(xr,func(xr),degree)
fi=lambda x: polynomial.polyval(x,p)
plot1(fi,fdata=(xr,func(xr)),label='Polynomial of degree %d'%degree,color=clr,extrapolation=True)
We could also go back to function approximation and fit polynomials of lower degree
Inner product
$$ \langle f,g \rangle = \int_D f(x)g(x)w(x)dx $$$ \{\phi_i\} $ is a family of orthogonal polynomials w.r.t. $ w(x) $ iff
$$ \langle \phi_i,\phi_j \rangle = 0, i\ne j $$Let $ \mathcal{P}_n $ denote the space of all polynomials of degree $ n $ over $ D $
$$ \lVert f - p \rVert_2 = \inf_{q \in \mathcal{P}_n} \lVert f - q \rVert_2 = \inf_{q \in \mathcal{P}_n} \left[ \int_D ( f(x)-g(x) )^2 dx \right]^{\tfrac{1}{2}} $$if and only if
$$ \langle f-p,q \rangle = 0, \text{ for all } q \in \mathcal{P}_n $$Orthogonal projection is the best approximating polynomial in L2-norm
Measures the absolute difference over the whole domain $ D $
What is the best polynomial approximation in the uniform (infinity, sup) norm?
$$ \lVert f - p \rVert_{\infty} = \inf_{q \in \mathcal{P}_n} \lVert f - q \rVert_{\infty} = \inf_{q \in \mathcal{P}_n} \sup_{x \in D} | f(x) - g(x) | $$Chebyshev proved existence and uniqueness of the best approximating polynomial in uniform norm.
Suppose $ f: [-1,1]\rightarrow R $ is $ C^k $ function for some $ k\ge 1 $, and let $ I_n $ be the degree $ n $ polynomial interpolation of $ f $ with nodes at zeros of $ T_{n+1}(x) $. Then
$$ \lVert f - I_n \rVert_{\infty} \le \left( \frac{2}{\pi} \log(n+1) +1 \right) \frac{(n-k)!}{n!}\left(\frac{\pi}{2}\right)^k \lVert f^{(k)}\rVert_{\infty} $$📖 Judd (1988) Numerical Methods in Economics
import numpy.polynomial.chebyshev as cheb
for degree, clr in (7,'b'),(9,'g'),(11,'m'):
fi=cheb.Chebyshev.interpolate(func,degree,[-5,10])
plot1(fi,fdata=(None,None),color=clr,label='Chebyshev with n=%d'%degree,extrapolation=True)
Generally much harder!