CODEWITHSUNDEEP
python-2.7numpy
AstroInt
AstroInt

Reputation: 63

How do I fit n data points with an (n-1)-degree polynomial

I'm trying to fit, for example, 5 data points to a polynomial of degree 4. This is my code:

linex = array([1,2,3,4,5])
liney = array([23.0, 28.521, 36.31542, 49.2312312, 78.231312])
p,sp = polyfit(linex, liney, deg=4, cov=True)
print p
print sp

When I try to run something as simple as that, I get

ValueError: operands could not be broadcast together with shapes (5,5) (0,)

Supposedly it is an error with line 608 in polynomial.py in polyfit, where it says 

return c, Vbase*fac

The fitting works if I change deg = 3, but that is not what I am looking for. Both arrays linex and liney have 5 points each as well. A brief search did not find a solution that addresses my problem.

Upvotes: 1

Views: 969

Answers (1)

featuredpeow
featuredpeow

Reputation: 2201

Internally polyfit finds solution Ax = y equation using least squares where A is Vandermonde matrix. For degree 4 and 5 points this equation has exact solution (A has shape 5x5, x - 5x1, y - 5x1). numpy.linalg.lstsq function returns numpy array with shape () in this case as estimate of error instead of returning vector with zeroes and code in polyfit fails to recognize that. However cov matrix (as it is computed there now) would be zero matrix even if it would recognized that. So to workaround it you can assume that cov matrix is zero matrix in your code for all degrees deg >= number_of_points - 1 or just call polyfit(linex, liney, deg=4) if you don't need covariance matrix.

Excerpts from docs:

polyfit

The solution minimizes the squared error

.. math ::
    E = \\sum_{j=0}^k |p(x_j) - y_j|^2

in the equations::

    x[0]**n * p[0] + ... + x[0] * p[n-1] + p[n] = y[0]
    x[1]**n * p[0] + ... + x[1] * p[n-1] + p[n] = y[1]
    ...
    x[k]**n * p[0] + ... + x[k] * p[n-1] + p[n] = y[k]

The coefficient matrix of the coefficients `p` is a Vandermonde matrix.

and numpy.linalg.lstsq:

Return the least-squares solution to a linear matrix equation.

Solves the equation `a x = b` by computing a vector `x` that
minimizes the Euclidean 2-norm `|| b - a x ||^2`.  The equation may
be under-, well-, or over- determined (i.e., the number of
linearly independent rows of `a` can be less than, equal to, or
greater than its number of linearly independent columns).  If `a`
is square and of full rank, then `x` (but for round-off error) is
the "exact" solution of the equation.

Parameters
----------
a : (M, N) array_like
    "Coefficient" matrix.

... returns: ... residuals : {(), (1,), (K,)} ndarray Sums of residuals; squared Euclidean 2-norm for each column in b - a*x. If the rank of a is < N or M <= N, this is an empty array. If b is 1-dimensional, this is a (1,) shape array. Otherwise the shape is (K,).

Upvotes: 2

Related Questions