title: Fitting multiple Lorentzians to Brillouin Spectrum using Scipy in Python 3
body: <p>I am trying to fit Brillouin Spectra (with several peaks) using scipy.optimize.curve_fit. 
  I have have multiple spectra with several peaks and I am trying to fit them with lorentzian functions (one Lorentzian per peak). I am trying to automate the process for bulk analysis (i.e., using the peak finding algorithm of scipy to get peak positions, peak widths and peaks heights and use them as initial guesses for the fit). I am now working on one spectrum to see if the general idea works, then I will extend it to be automatic and work with all the spectra I have. So far, I have done this:</p>
  
  <pre><code>import numpy as np
  import matplotlib.pyplot as plt
  from scipy.signal import find_peaks
  from scipy.optimize import curve_fit
  
  #import y data from linked google sheet 
  y_data = np.loadtxt( 'y_peo.tsv' )
  #define x data 
  x_data = np.arange( len( y_data ) )
  
  #find peak properties (peak position, amplitude, full width half maximum ) to use as 
  #initial guesses for the curve_fit function 
  pk, properties = find_peaks(y_data, height = 3, width = 3, prominence=0.1 ) #pk returns peaks position, properties returns 
  #other properties associated with the peaks
  I = properties ['peak_heights'] #amplitude
  fwhm = (properties['widths']) #full width half maximum 
  
  #define sum of lorentzians
  def func(x, *params): 
      y = np.zeros_like(x)
      for i in range(0, len(params), 3):
          x_0 = params[i]
          I = params[i+1]
          y_0 = params[i+2]
          y = y + (I*y_0**2)/(y_0**2 +(x-x_0)**2) 
      return y
  
  #initial guess list, to be populated with parameters found above (pk, I, fwhm)
  guess = [] 
  
  for i in range(len(pk)): 
      guess.append(pk[i])
      guess.append(I[i])
      guess.append(fwhm[i]) 
  
  #convert list to array
  guess=np.array(guess)
  
  #fit
  popt, pcov = curve_fit(func, x_data, y_data, p0=guess, method = 'lm',  maxfev=1000000)
  print(popt)
  fit = func(x_data, *popt)
  
  #plotting
  fig= plt.figure(figsize=(10,5))
  ax= fig.add_subplot(1,1,1)
  ax.plot(pk, y_data[pk], 'o', ms=5)
  ax.plot(x_data, y_data, 'ok', ms=1)
  ax.plot(x_data, fit , 'r--', lw=0.5)
  </code></pre>
  
  <p>Where <em>y_peo</em> is the dependent variable (here are the values as a google sheet file: <a href="https://docs.google.com/spreadsheets/d/1UB2lqs0Jmhthed7B0U9rznqcbpEMRmRX99c-iOBHLeE/edit?usp=sharing" rel="nofollow noreferrer">https://docs.google.com/spreadsheets/d/1UB2lqs0Jmhthed7B0U9rznqcbpEMRmRX99c-iOBHLeE/edit?usp=sharing</a>) and <em>pixel</em> is the independent variable (an arbitrary array of values created in the code). This is what I get: <a href="https://i.stack.imgur.com/bgl0s.png" rel="nofollow noreferrer">Result of spectrum fit</a>. Any idea why the last triplet is not fitted as expected? I checked and all the peaks are found correctly by the scipy.signal.find_peaks function - hence I think the problem is with the scipy.optimize.curve_fit, since I had to increase the number of maxfev for the fit to 'work'. Any idea on how to tackle this problem in a smarter way? </p>
  
  <p>Thanks in advance, </p>
  
  <p>Giuseppe. </p>