polyfit refining: setting polynomial to be always possitive

Question

I am trying to fit a polynomial to my data, e.g.

import scipy as sp

x = [1,6,9,17,23,28]

y = [6.1, 7.52324, 5.71, 5.86105, 6.3, 5.2]

and say I know the degree of polynomial (e.g.: 3), then I just use scipy.polyfit method to get the polynomial of a given degree:

+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

fittedModelFunction = sp.polyfit(x, y, 3)

func = sp.poly1d(fittedModelFunction) 

++++++++++++++++++++++++++++++ QUESTIONS: ++++++++++++++++++++++++++++++

1) How can I tell in addition that the resulting function func must be always positive (i.e. f(x) >= 0 for any x)?

2) How can I further define a constraint (e.g. number of (local) min and max points, etc.) in order to get a better fitting?

Is there smth like this: http://mail.scipy.org/pipermail/scipy-user/2007-July/013138.html but more accurate?

Answer 1

Always Positve

I haven't been able to find a scipy reference that determines if a function is positive-definite, but an indirect way would be to find the all the roots - Scipy Roots - of the function and inspect the limits near those roots. There are a few cases to consider:

No roots at all

Pick any x and evaluate the function. Since the function does not cross the x-axis because of a lack of roots, any positive result will indicate the function is positive!

Finite number of roots

This is probably the most likely case. You would have to inspect the limits before and after each root - Scipy Limits. You would have to specify your own minimum acceptable delta for the limit however. I haven't seen a 2-sided limit method provided by Scipy, but it looks simple enough to make your own.

from sympy import limit

// f: function, v: variable to limit, p: point, d: delta
// returns two limit values
def twoSidedLimit(f, v, p, d):
    return limit(f, v, p-d), limit(f, v, p+d)

Infinite roots

I don't think that polyfit would generate an oscillating function, but this is something to consider. I don't know how to handle this with the method I have already offered... Um, hope it does not happen?

Constraints

The only built-in form of constraints seems to be limited to the optimize library of SciPy. A crude way to enforce constraints for polyfit would be to get the function from polyfit, generate a vector of values for various x, and try to select values from the vector that violate the constraint. If you try to use filter, map, or lambda it may be slow with large vectors since python's filter makes a copy of the list/vector being filtered. I can't really help in this regard.