0

I'm using the top rated code from here

and modified it to fit 4 gaussian gives for 3 spectra. The code works when I'm just doing working one spectra at a time but I want to automate the code so that I do curve fitting on multiple spectra. Here is the code I have so far just for 3 spectra but I plan to do more. Note that my y_real is a list for 3 different fluxes while my xfit notes the wavelength which is the same range for all 3 spectra. I'm just using a few data points just for this example. My problem arises for res statement(ValueError: operands could not be broadcast together with shapes) which I do not know how to fix before I start plotting. 

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


xfit = np.array([6520.0,  6545.0,  6570.0,  6595.0,  6620.0])
y_real = np.array([[0.0, 8.9013813844724581, 33.114149624958607, 13.297472921668451,3.8540928722192049], [0.0, 11.491419225489725, 18.866478352891686, 9.0782365151426738, 1.8757625220616632], [0.0, 10.680394752562883, 21.394692831684814, 11.661265465293802, 0.9362924799953376]])


Lamda = np.array([ 6564.61389433,  6564.61389433,  6564.61389433])
delLamda1 = np.array([ 14.75496508,  14.75496508,  14.75496508])
delLamda2 = np.array([ 20.65455091,  20.65455091,  20.65455091])


def norm(xfit, mean, sd, a):
  norm = []
  for i in range(xfit.size):
    norm += [a*np.exp(-(xfit[i] - mean)**2/(2*sd**2))]
  return np.array(norm)

mean1, mean2 = 0, -2
std1, std2 = 0.5, 1

m, dm1, dm2, sd1, sd2, sd3, sd4, a, a2, a3, ab= [Lamda, -delLamda1, delLamda2, 1.0,1.0, 1.0, 2.0, 20.0, 20.0, 20.0, 10.0]
p = [m, dm1, dm2, sd1, sd2, sd3, sd4, a, a2, a3, ab] # Initial guesses for leastsq                                                                                       
y_init = norm(xfit, m, sd1, a) + norm(xfit, m + dm1, sd2, a2) + norm(xfit, m+dm2, sd3, a3) +  norm(xfit, m, sd4, ab) # For final comparison plot                        \


def res(p, y_real, xfit):
  m, dm1, dm2, sd1, sd2, sd3, sd4, a, a2, a3, ab= p
  m1 = m
  m2 = m1 + dm1
  m3 = m1 +dm2
  y_fit = norm(xfit, m1, sd1, a) + norm(xfit, m2, sd2, a2) + norm(xfit, m3, sd3, a3) +   norm(xfit, m1, sd4, ab)
  err = y_real - y_fit
  return err

plsq = leastsq(res, p, args = (y_real, xfit), ftol=1.0e-09, gtol=1.0e-09, xtol=1.0e-09, maxfev=2000)
Community
  • 1
  • 1
user1821176
  • 1,141
  • 2
  • 18
  • 29
  • I would suggest to try first to fit the spectra consecutively in a loop. The error comes from the fact tat y_real and y_fit in res() do not have the same shape. You can fix this by transposing one of them, e.g. y_fit.T but then you run into other problems which seems to be related to the fact that you use less data points than fit parameters, but to be honest, I do not really understand your code. – Christian K. Oct 25 '13 at 18:11
  • You might also want to have a look at peak-o-mat (http://lorentz.sf.net) which is a peak fitting application (on top of scipy) and which can be scripted entirely. It would take just a few lines to solve your fitting task. I can tell you how, if you like. – Christian K. Oct 25 '13 at 18:14
  • I am wondering how do you correctly tranpose y_fit in my code though? – user1821176 Oct 25 '13 at 20:23
  • A ndarray has a transpose attribute which returns the array with the two first axes swapped: y_fit.T – Christian K. Oct 26 '13 at 13:00

0 Answers0