python cross section curve fitting

Question

I have a set of points describing a cross section of a simple surface cavity/bulge. Polynomial approximation would be enough, but numpy.polyfit needs a certain degree. I thought of several iterations with different degrees to chose one with the smallest mean residual. Is there any existing function for this or a better way to get a good curve? Calculation time is of great importance: data sets are small (about 20 points), but there are thousands of them.

The initial task was to find this cavities geodesics on a very discrete surface - do any simplier methods exist?

Answer 1

Unless I'm mistaken, you'll always get the smallest residual with the largest degree allowed. If you allow the degree d to be chosen by the optimizer, it will choose infinitely large d, though in practice it will stop when d = dof where dof is the number of degrees of freedom of your fit (basically, the number of data points) at which point you'll have zero residual.

If you're not really interested in the true functional form of the fit and just want curves, you could look into interpolation, using the scipy.interpolate module, which is more flexible for the type of things you might be doing.

If speed matters, and format isn't that important, just try np.polynomial.polynomial.polyfit with some reasonable degree, and scipy.interpolate.UnivariateSpline, and see which is faster.

Looks to me like the spline is much faster, and for my example they give the same result (remember, a spline is basically just polynomials strung together)

import numpy as np
from numpy.polynomial import polynomial as poly
from scipy import interpolate as interp

import matplotlib.pyplot as plt

n = 20
x = np.linspace(0, 2*np.pi, n)
a = np.sin(x) + np.random.uniform(-.2, .2, n)

s = interp.UnivariateSpline(x, a)
p = poly.polyfit(x, a, 3)
p = poly.Polynomial(p)

plt.figure()
plt.plot(x, np.sin(x), '-', x, a, 'o', x, s(x), '--', x, p(x), '.-')

In [32]: timeit s = interp.UnivariateSpline(x, a)
10000 loops, best of 3: 22.1 µs per loop

In [33]: timeit p = poly.Polynomial(poly.polyfit(x, a, 3))
1000 loops, best of 3: 392 µs per loop

In [34]: timeit p = poly.polyfit(x, a, 3)
1000 loops, best of 3: 311 µs per loop

Including evaluation:

In [35]: timeit interp.UnivariateSpline(x, a)(x)
10000 loops, best of 3: 44.9 µs per loop

In [37]: timeit poly.Polynomial(poly.polyfit(x, a, 3))(x)
1000 loops, best of 3: 470 µs per loop

For fun, and to illustrate the concept of overfitting, this is a polynomial with d >= dof:

p = poly.Polynomial(poly.polyfit(x, a, x.size-1))
plot(x, np.sin(x), '-', x, a, 'o', x, p(x), ':')