Hello, I would need some help or advice on the following point : I'm trying to fit some X-ray scattering data with either gaussian or lorentzian functions using the least square algorithm (optimize.leastsq). I started with something quite simple, i.e. defining the residuals for the different functions in a similar way as in the tutorial and then calling the leastsq routine. In the end, I find the results a bit puzzling : if I define the data to be fitted as a gaussian + random noise (as in the tutorial), the fit with gaussian functions is perfect. The fit of the same data with lorentzian is not very good, I have the impression that the fit position is quite OK and the integrated intensity as well, but the maximum intensity is way off... If I use real data where I have two peaks and linear background, as in the example below, then I simply cannot get correct fit, no matter how good the initial estimates are... Can anyone give me a bit of light on this one ? Many thanks Aure def multifuncresiduals(multifunc,y,x,gaunum=2,lornum=0,linenum=0): err = y for i in range(gaunum): #gaussian params come first in list print multifunc[2*i],multifunc[2*i+1] err = err-1/(sqrt(2*pi)*multifunc[2*i+1])*\ exp(-(x-multifunc[2*i])**2/(2*multifunc[2*i+1]**2)) for i in range(lornum): #lorentzian params come first in list i += gaunum #need to correct for index print multifunc[2*i],multifunc[2*i+1] err = err-1/pi*0.5*multifunc[2*i+1]/((x-multifunc[2*i])**2+(0.5*multifunc[2*i+1])**2) for i in range(linenum): i += gaunum + lornum #need to correct for index print multifunc[2*i],multifunc[2*i+1] err = err-(multifunc[2*i+1]+multifunc[2*i]*x) print '----------------' return err def multifunc(y,x,multifuncinit): result = leastsq(multifuncresiduals,multifuncinit,args=(y,x)) return result if __name__ == '__main__': ######################################## #test for multiple gaussian fit based on input file # #read chiplot file content datafilename = 'D:/Boulot/Software/Python/azimuth.chi' x,y =readChi(datafilename) #define init values multigauinit = [90,100,270,100] #initial estimates result= multifunc(y,x,multigauinit)