2

This is my code (it's an example of a larger piece of code):

from scipy.optimize import curve_fit

def func(x, a, b):
   return a + b*x

xlist = [10, 30, 50, 70, 90, 110, 130, 150, 170, 190, 210, 230]
ylist = [0.0074999999999999997, 0.011875, 0.0057812499999999999, 0.0036458333333333334, 0.0020312500000000001, 0.0013125000000000001, 0.00098958333333333342, 0.00089285714285714283, 0.00074218750000000001, 0.00093749999999999997, 0.00071874999999999999, 0.00088068181818181821]

popt, pcov = curve_fit(func, xlist, ylist)

print(popt[0], popt[1])

As you can see I'm attempting a very simple fit with a very simple a + b*x function. The issue is that this returns the values:

(-119.99689110581872, 1.0)

for a=popt[0]and b=popt[1] but the same fit with zunzun.com gives much more reasonable values:

a = 7.8372289537296004E-03
b = -3.9402329475524466E-05

for the exact same function.

What am I doing wrong here?

Edit

Warren below stated this might be a bug. Should I report it as such or is this the expected behaviour?

Gabriel
  • 40,504
  • 73
  • 230
  • 404

2 Answers2

3

It works if you change the values in xlist to a numpy array:

In [38]: popt, pcov = curve_fit(func, array(xlist, dtype=float), ylist)

In [39]: popt
Out[39]: array([  7.83722896e-03,  -3.94023294e-05])

At first this looks like a bug, but what's happening is that the underlying code takes the argument xdata and passes it unchanged to your function. In your example, this means in the expression a + b*x, x is a Python list. That means that b*x is not doing the calculation that you want.

So to make your definition of func work, the argument xdata must be a numpy array. Or, you could define func this way:

def func(x, a, b):
    return a + b*np.asarray(x)
Warren Weckesser
  • 110,654
  • 19
  • 194
  • 214
  • Wait, I can also change the values in xlist to floating point simply assigning them as floating point (10., 30., etc) but that doesn't work. Why your method works but using floating points from the start does not? – Gabriel Jun 05 '13 at 00:05
  • I'm not sure--looking into it now. The problem exists in the function `leastsq`, for which `curve_fit` is a wrapper. – Warren Weckesser Jun 05 '13 at 00:15
2

The default initial parameters for scipy's curve fitting routine are all 1.0, if I recall correctly. While this works well for simple cases such as straight lines - since error space is simpler for straight line curve fitting -if you try to fit where error space is "bumpy and uneven" you will likely again see differences between your results and those of zunzun.com.

If you have some idea of the correct values (or range of values) where the optimum solution should exist, starting the scipy routine with such values will in almost all cases give better results than starting with all ones.

While zunzun.com allows entry of initial parameter "estimates" - or guesses - these are not required by the web site as a genetic algorithm is used to determine initial parameter estimates. The BSD-licensed Python fitting source code for the zunzun.com web site is at

http://code.google.com/p/pyeq2/

and it comes with a very wide range of examples. Look at the bottom of the zunzun.com web pages for links to the code as well. Let me know if you have any questions and I will be glad to help.

James Phillips zunzun@zunzun.com

James Phillips
  • 4,526
  • 3
  • 13
  • 11