Lorentzian scipy.optimize.leastsq fit to data fails

Question

Since I took a lecture on Python I wanted to use it to fit my data. Although I have been trying for a while now, I still have no idea why this is not working.

What I would like to do

Take one data-file after another from a subfolder (here called: 'Test'), transform the data a little bit and fit it with a Lorentzian function.

Problem description

When I run the code posted below, it does not fit anything and just returns my initial parameters after 4 function calls. I tried scaling the data, playing around with ftol and maxfev after checking the python documentation over and over again, but nothing improved. I also tried changing the lists to numpy.arrays explicitely, as well as the solution given to the question scipy.optimize.leastsq returns best guess parameters not new best fit, x = x.astype(np.float64). No improvement. Strangely enough, for few selected data-files this same code worked at some point, but for the majority it never did. It can definitely be fitted, since a Levenberg-Marquard fitting routine gives reasonably good results in Origin.

Can someone tell me what is going wrong or point out alternatives...?

import numpy,math,scipy,pylab
from scipy.optimize import leastsq
import glob,os
for files in glob.glob("*.txt"):
    x=[]
    y=[]
    z=[]
    f = open(files, 'r')
    raw=f.readlines()
    f.close()
    del raw[0:8]       #delete Header
    for columns in ( raw2.strip().split() for raw2 in raw ):  #data columns
        x.append(float(columns[0]))
        y.append(float(columns[1]))
        z.append(10**(float(columns[1])*0.1)) #transform data for the fit
    def lorentz(p,x):
        return (1/(1+(x/p[0] - 1)**4*p[1]**2))*p[2]
    def errorfunc(p,x,z):
        return lorentz(p,x)-z

    p0=[3.,10000.,0.001]

    Params,cov_x,infodict,mesg,ier = leastsq(errorfunc,p0,args=(x,z),full_output=True)
    print Params
    print ier

Answer 1

Without seeing your data it is hard to tell what is going wrong. I generated some random noise and used your code to perform a fit to it. Everything works okay. This algorithm does not allow for parameter boundaries so you may run into problems if your p0 is close to zero. I did the following:

import numpy as np
from scipy.optimize import leastsq
import matplotlib.pyplot as plt

def lorentz(p,x):
    return p[2] / (1.0 + (x / p[0] - 1.0)**4 * p[1]**2)

def errorfunc(p,x,z):
        return lorentz(p,x)-z

p = np.array([0.5, 0.25, 1.0], dtype=np.double)
x = np.linspace(-1.5, 2.5, num=30, endpoint=True)
noise = np.random.randn(30) * 0.05
z = lorentz(p,x)
noisyz = z + noise

p0 = np.array([-2.0, -4.0, 6.8], dtype=np.double) #Initial guess
solp, ier = leastsq(errorfunc, 
                    p0, 
                    args=(x,noisyz),
                    Dfun=None,
                    full_output=False,
                    ftol=1e-9,
                    xtol=1e-9,
                    maxfev=100000,
                    epsfcn=1e-10,
                    factor=0.1)

plt.plot(x, z, 'k-', linewidth=1.5, alpha=0.6, label='Theoretical')
plt.scatter(x, noisyz, c='r', marker='+', color='r', label='Measured Data')
plt.plot(x, lorentz(solp,x), 'g--', linewidth=2, label='leastsq fit')
plt.xlim((-1.5, 2.5))
plt.ylim((0.0, 1.2))
plt.grid(which='major')
plt.legend(loc=8)
plt.show()

This yielded a solution of:
solp = array([ 0.51779002, 0.26727697, 1.02946179])
Which is close to the theoretical value:
np.array([0.5, 0.25, 1.0])