Проблемы с подбором лоренцевского пика в Python 2 с использованием lmfit

import numpy as np
import os
import matplotlib.pyplot as plt
from lmfit.models import GaussianModel, ExponentialModel, LorentzianModel, VoigtModel
import scipy
from scipy.optimize import leastsq
###############################################################################
###############################################################################

os.chdir('D:/reference')

def readfio(filename):
"""
READ *.fio file into dictionary of arrays
"""
    cols = []
# Open file
    f = open(filename +'.fio')
    counter = 0
    for n, line in enumerate(f):
        if line.strip().startswith('Col'):
            linestoskip = n +  1
            cols.append(line.split()[2]) #To have only the column headers in a list
            counter = counter + 1 #gives the number of columns, in principle it is not necessary
    f.close() # close the file 
    # Read numerical data without header
    #print (cols)
    data = np.genfromtxt(filename+'.fio',skip_header=linestoskip, skip_footer =1)
    return data, cols
###############################################################################
###############################################################################
def fitting():
    data, cols = readfio(filename)
    x = (2*np.pi/wavelength(energy))*(np.sin(np.radians(data[:,cols.index('om')]))+
    np.sin(np.radians(data[:,cols.index('tt')])-np.radians(data[:,cols.index('om')])))  
    y = data[:,cols.index('signalcounter_atten')]/data[:,cols.index('petra_beamcurrent')]

    exp_mod = ExponentialModel(prefix='exp_')
    pars = exp_mod.guess(y, x=x)

    gauss1  = LorentzianModel(prefix='g1_')
    pars.update( gauss1.make_params())

    pars['g1_center'].set(3.33, min=3.2, max=3.45)
    pars['g1_sigma'].set(0.02, min=0.01)
    pars['g1_amplitude'].set(2000, min=1)

    gauss2  = LorentzianModel(prefix='g2_')

    pars.update(gauss2.make_params())

    pars['g2_center'].set(3.59, min=3.45, max=3.7)
    pars['g2_sigma'].set(0.01, min=0.001)
    pars['g2_amplitude'].set(10000, min=1)

    mod = gauss1 + gauss2 + exp_mod


    init = mod.eval(pars, x=x)
    plt.semilogy(x, y)
    plt.semilogy(x, init, 'k--')

    out = mod.fit(y, pars, x=x)

    print(out.fit_report(min_correl=0.01))

    plt.semilogy(x, out.best_fit, 'r-')
    plt.show()
###############################################################################    ###############################################################################

def wavelength(energy):
    h = 6.626e-34  #joules
    c = 2.998e8    #m/sec
    eV = 1.602e-19 #Joules
    wavelength = h*c/(energy*1000*eV)*1e10       
    return wavelength

##############################################################################    
energy = 10.000 # in KeV    
filename = 'alignment_00241'
fitting()

У меня есть эта программа, чтобы подогнать два пика с помощью функции Лоренца. Однако я не подхожу наилучшим образом. Кто-нибудь может мне с этим помочь. Это тот файл, который я использовал для примерки. Файл данных можно загрузить по этой ссылке


person Karthick    schedule 19.05.2016    source источник
comment
Предполагая, что вы используете правильную модель и оптимальное значение находится в диапазоне параметров, но соответствие возвращает неоптимальное значение, затем либо установите другое начальное значение, либо используйте нелокальный оптимизатор, установив method = diff_evolution в качестве опции. в мод.фит.   -  person Neapolitan    schedule 19.05.2016


Ответы (1)


Я считаю, что в ваших y-данных есть NaN, что мешает подгонке. Вы, вероятно, увидите NaN для значений параметров. Если да, то сделайте что-нибудь вроде ::

for ix in np.where(np.isnan(y)):
    y[ix]  = (y[ix-1] + y[ix+1]) / 2.0

перед выполнением подгонки может помочь.

person M Newville    schedule 20.05.2016