Best fit algebraic curve

Intro

The example is inspired by ​this presentation

The get_best_fit(pts, max_deg) function defined below expects two arguments:

• a list of points pts
• an integer max_deg

It computes and returns the algebraic curve of degree at most max_deg which fits best the list of points in the sense that if f(x, y) = 0 is the equation of the curve then the euclidian norm of the vector [f(x_1, y_1), ... f(x_n, y_n)] is minimised, (x_1, y_1) ... (x_n, y_n) being the coordinates of the points pts. For the minimising to be possible, we need to fix the norm of the vector made of the coefficients of the curve: in the method described below, we require the maximum norm of this vector to be 1.

The Test class allows playing with a parametrically defined curve by creating a sample list of points and allowing to add best fit curves of varying degrees (see Try it section below).

Code

Copy and paste this into the Python input field

from linalg import Vector, Matrix, SVDecomposition
import math

def get_coords(pts):
    return [(p.x.value, p.y.value) for p in pts]

def powers(n):
    return ((i, d-i) for d in xrange(n, -1, -1) for i in xrange(d, -1, -1))

def monomials(x, y, n):
    return [x**i * y**j for i, j in powers(n)]

def get_polyfunc(coeffs, n):
    def f(x, y):
        return sum(c*m for c, m in zip(coeffs, monomials(x, y, n)) if c)
    return Function(f)

def find_coeffs(coords, n):
    matrix_data = [monomials(x, y, n) for x, y in coords]
    last_row = [0]*len(matrix_data[0])
    matrix_data.append(last_row)
    B = Vector([0]*len(coords) + [1])
    norm_min = 1e10
    for i, p in enumerate(powers(n)):
        last_row[i] = 1
        last_row[i-1] = 0
        A = Matrix(matrix_data)
        X = SVDecomposition.solve(A, B)
        norm = (A*X - B).norm
        if norm < norm_min:
            coeffs = X
            norm_min = norm
    return coeffs

def get_best_fit(pts, max_deg):
    coords = get_coords(pts)
    coeffs = find_coeffs(coords, max_deg)
    poly = get_polyfunc(coeffs, max_deg)
    return poly.implicitcurve

class Test(object):
    def __init__(self, f, n=10, start=0, end=2*math.pi):
        self.pts = [Point(*f(start + (end-start)*i/n)) for i in range(n)]
    def best_fit(self, n):
        return get_best_fit(self.pts, n)

def circle_eqn(r):
    return lambda t: (a*math.cos(t), b*math.sin(t))

def limacon(a, b):
    return lambda t: (a*(1+math.cos(2*t))+b*math.cos(t),
                      b*math.sin(t)+a*math.sin(2*t)/2)

Try it

Try 20 points on a limacon r=1+2cost. Type this

lim12 = Test(limacon(1, 2), 20)

Now you can see the points in the euclidian view. Now find the best fit algebraic curves of degree 1, 2 and 3

lim12.best_fit(1)
lim12.best_fit(2)
lim12.best_fit(3)

We're getting closer!

lim12.best_fit(4)

Now we've got it! In the screenshot below, I've shown

• the linear best fit in yellow
• the quadratic one in blue
• the cubic in green
• the quartic in red.

Make your own

Find the parametric equation of your favourite curve and try it!