polyfit

A follow-up to the previous post.

Polynomial fitting is also very easy with the numpy packages polyfit and poly1d.

In [196]: x = range(100)
In [197]: y = randn(100)
In [198]: plot(x,y)
Out[198]: []

Here I am asking polyfit to fit me a 2nd degree polynomial.

In [199]: polyfit(x,y,2)
Out[199]: array([-0.00018313,  0.01669275, -0.09621319])

The polyfit function returns the polynomial coefficients in a list.
If I want to use them directly as a fit function, just embed them in a new polynomial object:

In [200]: fitfunc = poly1d(polyfit(x,y,2))
In [201]: plot(x,fitfunc(x))
Out[201]: []

Saving the plot like this

In [202]: savefig('/Users/maye/Desktop/blog_polyfit.png')

and looks like this:

Polynomials with Python

Seriously, can it be any easier? ;)

If you are not in a pylab session, import the module like this:

In [148]: from numpy import poly1d

 Otherwise, just "import poly1d" should work.
Now let's get a polynomial for the coefficients of [3,2,1] (always in decreasing order!):

In [149]: p = poly1d([3,2,1])

Printing it provides a semi-analytical printout:

In [150]: print p
   2
3 x + 2 x + 1

Applying new x values to it is easy, because the poly1d object is a function:

In [152]: newx = linspace(0,10,10)
In [153]: p(newx)
Out[153]:
array([   1.        ,    6.92592593,   20.25925926,   41.        ,
         69.14814815,  104.7037037 ,  147.66666667,  198.03703704,
        255.81481481,  321.        ])

Lots of other things are possible with this object. IPython's object inspection makes it easy to discover them:

In [154]: p.
p.coeffs    p.deriv     p.integ     p.order     p.variable

In [155]: p.deriv()
Out[155]: poly1d([6, 2])
In [156]: pderiv = p.deriv()
In [157]: print pderiv
6 x + 2

Roots for this polynomial can be either determined by the roots function that is imported in a pylab session (or importable like from numpy import roots)

In [158]: roots(p)
Out[158]: array([-0.33333333+0.47140452j, -0.33333333-0.47140452j])
In [159]: p.r
Out[159]: array([-0.33333333+0.47140452j, -0.33333333-0.47140452j])

PS: One of these days I really have to find out how to do code high-lighting in Blogger, or, preferably, go all the way and do IPython notebook posts.

Intuitive division in Python

When one is hot-coding some idea, it is all about writing down an idea fast, not pretty (even so it's quite hard to write ugly code in Python due to its indental paradigm), not clean or precise, just trying out an idea.

It is then that it could happen to overlook that a division of 2 integer numbers is forcing an integer solution of this calculation (the underlying C-paradigm for divisions in Python is enforcing this), so 3/4 is 0, not 0.75, because one forgot to write something like 3/4.0 to enforce the floating point division.

Many errors in un-counted Python (and C) classes were caused by this intricacy, maybe even some good ideas have been discredited because the code produced weird results. And because of the caused frustrations, the Python inventors have decided to get rid of this, beginning with the 3.x versions of Python, so 3/4 is actually 0.75.

And because this was such a great idea, it was decided to make it available as an option in Python 2.x (before you judge this as a minor change, admit to yourself how often you have been caught by this bug!?).

So, to enable this 'safe' division in Python 2.x code, you have to write the following line before all other imports:

from __future__ import division

(thats 2 underscores in front AND behind the future)

Enjoy your hot-coding now! ;)