Use of curve_fit to fit data

Question

I'm new to scipy and matplotlib, and I've been trying to fit functions to data. The first example in the Scipy Cookbook works fantastically, but when I am trying it with points read from a file, the initial coefficients I give (p0 below) never seem to actually change, and the covariance matrix is always INF.

I've tried to fit even data following a line, to no avail. Is it a problem with the way I am importing the data? If so, is there a better way to do it?

import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import scipy as sy

with open('data.dat') as f:
    noms = f.readline().split('\t')

    dtipus = [('x', sy.float32)] + [('y', sy.float32)]

    data = sy.loadtxt(f,delimiter='\t',dtype=dtipus)

    x = data['x']
    y = data['y']

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

    p0 = sy.array([1,1,1])

    coeffs, matcov = curve_fit(func, x, y, p0)

    yaj = func(x, coeffs[0], coeffs[1], coeffs[2])

    print(coeffs)
    print(matcov)

    plt.plot(x,y,'x',x,yaj,'r-')
    plt.show()

Thanks!

Solution 1

It seems to me that the problem is indeed in how you import your data. Faking this datafile:

$:~/temp$ cat data.dat
1.0  2.0
2.0  4.2
3.0  8.4
4.0  16.1

and using the pylab's loadtxt function for reading:

import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import scipy as sy
import pylab as plb  

data = plb.loadtxt('data.dat')  
x = data[:,0]
y= data[:,1]

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

p0 = sy.array([1,1,1])
coeffs, matcov = curve_fit(func, x, y, p0)

yaj = func(x, coeffs[0], coeffs[1], coeffs[2])
print(coeffs)
print(matcov)

plt.plot(x,y,'x',x,yaj,'r-')
plt.show()

works for me. By the way, you can use dtypes to name the columns.

Solution 2

The underlying problem with your load data is that you cast it to float32, but in scipy 0.10.1, curve_fit works with float64 but not float32 (it's a bug, not a feature). Your example works with float64.

Comments

[Ironil]

Yes, thank you! Loading the data with just loadtxt made it. It seems I was trying to do it the hard way, but I'll keep investigating what was going wrong.