python - 拟合后如何获取lmfit参数呢？

Question

我编写了一个程序来拟合一些拉曼光谱峰。我需要返回拟合参数(位置、振幅、HWHM)。

我使用 modul lmfit 创建了一个带有约束的洛伦兹峰。

根据我的图形图，我在拟合峰和原始数据之间有很好的一致性。但是当涉及到 after 拟合参数时，我遇到了一个问题，程序只返回初始值。

我绑定(bind)了“report_fit 模块”并更改了初始参数但没有成功。参数值不会演变。

令我困扰的是，该程序在我同事的 PC 上运行，但在我的 PC 上运行不正常。所以问题可能出在我的python版本上。

我使用的是 spyder 2.3.9，python 3.4 在 windows 10 下安装了 anaconda。lmfit 模块 0.9.3 似乎部分工作，因为我可以获得良好的拟合协议(protocol)(来自图 plt.plot)。但是我无法在拟合后返回参数值。

---

这是我的代码:

#!/usr/bin/python3
# -*- coding:utf-8 -*-

import os
import numpy as np
import matplotlib.pyplot as plt
from math import factorial
from scipy.interpolate import interp1d
from lmfit import minimize, Parameters  #,report_fit

##############################################################################
# Fit of Raman peaks

def fmin150(pars,x,y):  
    amp= pars['Amp_150'].value
    cent=pars['Cent_150'].value
    hwhm=pars['Wid_150'].value
    a=pars['a_150'].value
    b=pars['b_150'].value
    peak = (amp*hwhm)/(((x-cent)**2)+(hwhm**2)) + ((a*x)+b)
    return peak - y    

def fmin220(pars,x,y):  
    amp= pars['Amp_220'].value
    cent=pars['Cent_220'].value
    hwhm=pars['Wid_220'].value
    a=pars['a_220'].value
    b=pars['b_220'].value
    peak = (amp*hwhm)/(((x-cent)**2)+(hwhm**2)) + ((a*x)+b)
    return peak - y   

def fmin2d(pars,x,y):  
    amp= pars['Amp_2D'].value
    cent=pars['Cent_2D'].value
    hwhm=pars['Wid_2D'].value
    a=pars['a_2D'].value
    b=pars['b_2D'].value
    peak = (amp*hwhm)/(((x-cent)**2)+(hwhm**2))  + ((a*x)+b)
    return peak - y 

##############################################################################
def fit(filename):
    """
    Fit Raman spectrum file
    Return filename, for each peak  (*f*whm, position, height) 
    PL (Position, Intensity, Area)
    """
    print ("----------------------------")
    print("Treating file : ")
    print(filename)
    try:
        data = np.loadtxt(filename)
    except:
        print("Unable to load file")

    xx=data[:,0]
    yy=data[:,1]
 #### Define fitting interval  ##### 
    # Cu oxides (unités en cm-1)
    xminim150 = 120
    xmaxim150 = 170
    xminim220 = 175
    xmaxim220 = 275
    xminim300 = 280
    xmaxim300 = 340
    xminim640 = 345
    xmaxim640 = 800
    # Graphene
    xminimG = 1540
    xmaximG = 1680
    xminim2D = 2600
    xmaxim2D = 2820

    # Definition Baground = place without the fitted peaks
    zone1 = (xx > min(xx)) & (xx < xminim150)
    zone2 = (xx > xmaxim150) & (xx < xminim220)
    zone3 = (xx > xmaxim220) & (xx < xminim300)
    zone4 = (xx > xmaxim300) & (xx < xminim640)
    zone5 = (xx > xmaxim640) & (xx < xminimG)
    zone6 = (xx > xmaximG) & (xx < xminim2D)
    zone7 = (xx > xmaxim2D) & (xx < max(xx))

    x_BG = np.concatenate((xx[zone1],xx[zone2],xx[zone3],xx[zone4],xx[zone5],xx[zone6],xx[zone7]))
    y_BG = np.concatenate((yy[zone1],yy[zone2],yy[zone3],yy[zone4],yy[zone5],yy[zone6],yy[zone7]))


    #Creation de l'interpolation lineaire
    f1 = interp1d(x_BG, y_BG, kind='linear')
    xinterpmin= x_BG[0]   # valeur de x_BG min
    xinterpmax= x_BG[len(x_BG)-1]  # valeur de x_BG max
    nbxinterp = len(xx) * 4 #(nb de point en x)* 4 pour une interpolation correcte

    xnew = np.linspace(xinterpmin, xinterpmax, num=nbxinterp, endpoint=True)
    ynew= f1(xnew)
##########################  Fit 2D peaks  ############################### 
    # create a set of Parameters
    pars150 = Parameters()
    pars220 = Parameters()
    pars300 = Parameters()
    pars640 = Parameters()
    parsg = Parameters()  
    pars2d = Parameters()

            #####   Cu2O pic 150 cm-1     #####
    pars150.add('Amp_150', value=10, min=0.0001, max=100000)  # Amplitude ~intensity
    pars150.add('Cent_150', value=150, min=140, max=160)      # Center  ~position
    pars150.add('Wid_150', value=5, min=4, max=15 )           # Width is the HWHM
    pars150.add('a_150', value=1, min=-100000, max=100000) 
    pars150.add('b_150', value=10, min=-100000, max=100000)
            #####   Cu2O pic 220 cm-1     #####
    pars220.add('Amp_220', value=10, min=0.0001, max=100000)
    pars220.add('Cent_220', value=220, min=200, max=230)
    pars220.add('Wid_220', value=5, min=4, max=15 )
    pars220.add('a_220', value=1, min=-100000, max=100000) 
    pars220.add('b_220', value=10, min=-100000, max=100000)
        #####    Graphene 2D     #####
    pars2d.add('Amp_2D', value=15, min=0.0001, max=100000)
    pars2d.add('Cent_2D', value=2711, min=2670, max=2730)
    pars2d.add('Wid_2D', value=15, min=4, max=25 )
    pars2d.add('a_2D', value=1, min=-100000, max=100000) 
    pars2d.add('b_2D', value=10, min=-100000, max=100000)

    # define x for each peaks
    interv_150 = (xx > xminim150) & (xx < xmaxim150)
    x_150 = xx[interv_150]
    y_150 = yy[interv_150]
    interv_220 = (xx > xminim220) & (xx < xmaxim220)
    x_220 = xx[interv_220]
    y_220 = yy[interv_220]

    interv_2D = (xx > xminim2D) & (xx < xmaxim2D)
    x_2D = xx[interv_2D]
    y_2D = yy[interv_2D]
    ###########################################################
    # Performe FIT with leastsq model  ###########
    result_150 = minimize(fmin150, pars150, args=(x_150, y_150))
    result_220 = minimize(fmin220, pars220, args=(x_220, y_220))
    result_2D = minimize(fmin2d, pars2d, args=(x_2D, y_2D))

    # calculate final result
    final_150 = y_150 + result_150.residual
    final_220 = y_220 + result_220.residual
    final_2D = y_2D + result_2D.residual

    ###########################################################
    # Parameter after fit #
    amp_150= pars150['Amp_150'].value
    cent_150=pars150['Cent_150'].value
    fwhm_150=2*pars150['Wid_150'].value

    amp_220= pars220['Amp_220'].value
    cent_220=pars220['Cent_220'].value
    fwhm_220=2*pars220['Wid_220'].value

    amp_2D= pars2d['Amp_2D'].value
    cent_2D=pars2d['Cent_2D'].value
    fwhm_2D=2*pars2d['Wid_2D'].value

    ###########
    #Plot data#
    plt.plot(xx, yy, 'k+' ,x_150, final_150, 'r', x_220, final_220,'r', x_2D, final_2D,'b')
    plt.xlabel(r'Raman shift (cm$^{-1}$)', fontsize=14)
    plt.ylabel('Intensity (a.u.)', fontsize=14)
    plt.xlim(0,3000)
    plt.title(filename, fontsize=16)
    savename=filename[:-4]+".png"
    print(savename)
    plt.savefig(savename)
    plt.clf()

    return filename, amp_150, cent_150, fwhm_150, amp_220, cent_220, fwhm_220,  amp_2D, cent_2D, fwhm_2D


def main():
    """main program loop"""
    print("----------------------------")
    liste = []
    for filename in os.listdir(os.getcwd()):
      if filename.endswith(".txt"):
        liste.append(filename)

    f1 = open("TestFit_all.dat","w")
    header = "Filename\tI_150\tCentre_150\tFWHM_150\tI_220\tCentre_220\tFWHM_220\tI_300\tCentre_300\tFWHM_300\tI_640\tCentre_640\tFWHM_640\tI_G\tCentre_G\tFWHM_G\tI_2D\tCentre_2D\tFWHM_2D\tI_PL\tI_PL_err\tCent_PL\tCent_PL_err\tArea1000_PL\n"
    f1.write(header)
    for i in liste:
        txt=fit(i)
        print(txt)
        #text = str(txt[0])+"\t"+str(txt[1])+"\t"+str(txt[2])+"\t"+str(txt[3])+"\t"+str(txt[4])+"\t"+"\n"
        text = str(txt[0])+"\t"+str(txt[1])+"\t"+str(txt[2])+"\t"+str(txt[3])+"\t"+str(txt[4])+"\t"+str(txt[5])+"\t"+str(txt[6])+"\t"+str(txt[7])+"\t"+str(txt[8])+"\t"+str(txt[9])+"\t"+"\t"+"\n"
        f1.write(text)

    f1.close()

    print("----------------------------")
    print("Done")
###################################################################
# start

if __name__ == "__main__":
  main()  

谢谢你的帮助:)

A working fit exemple

请将您的回复发送至 deniz.cakir@etu.umontpellier.fr

Answer 1

您的脚本可以更短以说明问题。但是，最终参数保存在 result_150.params 中，依此类推。传递给 minimize() 的参数在拟合中没有改变。

好吧，对于 lmfit 版本 0.9.0 和更高版本来说都是如此。也许您的同事在 lmfit 上有旧版本。