0

我正在使用此处评分最高的代码

并对其进行修改以适合 3 个光谱的 4 个高斯给出。当我一次只处理一个光谱但我想自动化代码以便我对多个光谱进行曲线拟合时，该代码有效。这是我到目前为止仅用于 3 个光谱的代码，但我计划做更多。请注意，我的 y_real 是 3 个不同通量的列表，而我的 xfit 记录了所有 3 个光谱范围相同的波长。我只是在这个例子中使用了一些数据点。我的问题出现在 res 语句（ValueError：操作数不能与形状一起广播），在开始绘图之前我不知道如何解决。 

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)
4

0 回答 0