TypeError from scipy.optimize.curve_fit

Question

I am trying to fit a data set to an exponential model using scipy. However, the covariance matrix that is returned is 'inf' and I receive the following error:

Traceback (most recent call last): File "C:\Users\Christopher\Desktop\advancedLab\03_ML_Mogni\analyzeTHIS.py", line 21, in print 'A = ', ans[0], '+/-', cov[0][0] TypeError: 'float' object has no attribute '__getitem__'

Now, I don't understand why my covariance matrix comes back as infinite... I know that is the source of this error message. My data set has times, which are between 400 and 20000 units, probabilities, which are between 0.005 and 0.3 units, and their respective standard deviations. Here is my code:

import numpy
import scipy.optimize, scipy.stats.stats
import pylab

use = raw_input('File: ')

time = numpy.loadtxt(use + 'time.txt', unpack = 'True')
prob = numpy.loadtxt(use + 'prob.txt', unpack = 'True')
std = numpy.loadtxt(use + 'std.txt', unpack = 'True')
stdProb = numpy.loadtxt(use + 'stdProb.txt', unpack = 'True')

guess = prob
guess.tolist()

def muon_model(x,A):
    return A*numpy.exp(-A*x)

ans, cov = scipy.optimize.curve_fit(muon_model, time, prob, p0 = 2200, sigma = stdProb, maxfev = 10000)

print cov
print 'A = ', ans[0], '+/-', cov[0][0]

numpy.savetxt(use + '_counts_best_fit.txt', ans)
numpy.savetxt(use + '_counts_cov.txt', cov)

x = numpy.linspace(0,22000,500)
y = ans[0]*numpy.exp(-ans[0]*x)

model = muon_model(time,ans[0])
pylab.plot(x,y,label='Model')
pylab.errorbar(time, prob, xerr = std, yerr = stdProb, label = 'Data', fmt = 'o')
pylab.legend()
pylab.xlabel('Time (nanoseconds)')
pylab.ylabel('Probability')
pylab.savefig(use + '_counts_with_fit.png')
pylab.clf()

mod = []
ind = 0
while ind < len(time):
    mod.append(math.exp(-1000*ans[0]*ind) - math.exp(-1000*ans[0]*(ind + 1)))
residual = prob - mod

chisq = (residual**2/std**2).sum()
dof = len(counts)-len(ans)
print 'Chi-Sq = ', chisq
print 'dof = ', dof
print 'Reduced Chi-Sq = ', chisq/dof
print 'PTE = ', scipy.stats.stats.chisqprob(chisq,dof)

Any help would be greatly appreciated.

EDIT: I changed my guess parameter to 0.5, because I put the wrong guess in. I still ended up with the same answer. Unfortunately, my data points are not equally spaced. What I did was bin my times in 1000 unit intervals, and took my data point in that bin to be the average value of all the values in that bin. I am still getting the same error no matter which guess I use.

Answer 1

welcome to SciComp. First of all, if you need purely exponential fitting and have equidistant data, Prony's method (see expfit for an implementation in Octave which can easily be ported to Python) is much more appropriate (and numerically stable).

That said, if I look at your initial guess for the parameter (2200) and you say that you have time spanning from 0 to 22000, your exponential model will give zeros a whole lot of the time (in double precision). And this will lead to numerical trouble. Are you sure about this parameter (no scaling error)? If this is not the issue, can you make (a smaller version of) your data set available somewhere?

Disclaimer: I am the original author of Octave's expfit function (but it has been improved by others later on).