使用lmfit库拟合双缝实验数据时遇报错求助
双缝实验数据拟合问题(基于lmfit库)
我有双缝实验的实验数据,希望用lmfit库拟合数值曲线,因为相比Python其他拟合工具,我更偏好它的拟合结果输出形式。以下是我的代码和遇到的问题:
初始代码与Key Error: N报错
from PIL import Image from scipy.optimize import curve_fit import matplotlib.pyplot as plt import numpy as np import pandas as pd from scipy.signal import argrelextrema from lmfit.models import Model #Determine the observed intensity distribution img_data = pd.read_csv('Image31 intensity profile.csv') pixels = np.array(img_data['Distance_(_)']) counts = np.array(img_data['Gray_Value']) norm_intensity = counts / np.linalg.norm(counts) # Double slit interferometry coherent beam def keV2m(keV): """Converts Photon Energy [keV] to wavelength [m]. .. note:: Calculation in Vacuum. """ wl = 1./(keV*1000)*4.1356*(10**(-7))*2.9998 return wl ener = 10#kev I_0 = 1 a = 2.0e-6 # slit width in m d = 10e-6 # slit seperation in m L = 6.0 # slit to detec dist in m detec_pixel_size = 3.1e-6 #pixel size x = np.arange(-500,500)*detec_pixel_size + 1e-9 def coh_inter_patt(x): lam = keV2m(ener) sinq_part = ((np.sin((np.pi*x*a)/(lam*L))/((np.pi*x*a)/(lam*L)))**2) sinq_part[sinq_part == np.nan] = 1 inter_pat = I_0*(np.cos((np.pi*x*d)/(lam*L))**2)*sinq_part return inter_pat def coh_inter_patt_sft(x,detec_pixel_size,sft=10): lam = keV2m(ener) sinq_part = ((np.sin((np.pi*(x+(sft*detec_pixel_size))*a)/(lam*L))/((np.pi*(x+(sft*detec_pixel_size))*a)/(lam*L)))**2) sinq_part[sinq_part == np.nan] = 1 inter_pat = I_0*(np.cos((np.pi*(x+(sft*detec_pixel_size))*d)/(lam*L))**2)*sinq_part return inter_pat inter_pat = coh_inter_patt(x) inter_pat_sft = coh_inter_patt_sft(x,detec_pixel_size,sft=10) max_vals = np.sort(argrelextrema(inter_pat, np.greater)[0]) min_vals = np.sort(argrelextrema(inter_pat, np.less)[0]) visibility = (inter_pat[max_vals] - inter_pat[min_vals[1:]])/(inter_pat[max_vals] + inter_pat[min_vals[1:]]) # 'Blurred' patterns extension = 40 # 10*detec_pixel_size long extended source blur_patt = np.zeros_like(coh_inter_patt(x)) for i in range(extension): blur_patt += coh_inter_patt_sft(x,detec_pixel_size,sft=i) blur_patt += coh_inter_patt_sft(x,detec_pixel_size,sft=-i) blur_patt = blur_patt/np.linalg.norm(blur_patt) max_vals = np.sort(argrelextrema(inter_pat, np.greater)[0]) min_vals = np.sort(argrelextrema(inter_pat, np.less)[0]) visibility = (inter_pat[max_vals] - inter_pat[min_vals[1:]])/(inter_pat[max_vals] + inter_pat[min_vals[1:]]) blur_max_vals = np.sort(argrelextrema(blur_patt, np.greater)[0]) blur_min_vals = np.sort(argrelextrema(blur_patt, np.less)[0]) blur_visibility = (blur_patt[max_vals] - blur_patt[min_vals[1:]])/(blur_patt[max_vals] + blur_patt[min_vals[1:]]) #Look at plot output fig = plt.figure(figsize=(10,6)) ax1 = fig.add_subplot(111) ax1.plot(pixels, norm_intensity, color='red', label='Observed') ax1.plot(blur_patt,label='Analytical, vis %f'%blur_visibility.mean(), color='k') ax1.legend() #iterative curve fitting def Blurred_fringes(N): extension = N # 10*detec_pixel_size long extended source blur_patt = np.zeros_like(coh_inter_patt(x)) for i in range(extension): blur_patt += coh_inter_patt_sft(x,detec_pixel_size,sft=i) blur_patt += coh_inter_patt_sft(x,detec_pixel_size,sft=-i) #inter_pat = inter_pat/np.linalg.norm(inter_pat) blur_patt = blur_patt/np.linalg.norm(blur_patt) return blur_patt model = Model(Blurred_fringes) params1 = model.make_params(N=35) results1 = model.fit(counts, params1, x=pixels)
运行上述代码时出现Key Error: N报错。我清楚lmfit的fit()函数需要make_params()生成的字典对象,且模型函数的第一个参数会被默认视为自变量。
修改后的代码与NaN值报错
为解决上述问题,我修改了曲线拟合部分的代码:
X = np.arange(0,1000,1) def Blurred_fringes(x, N): extension = N # 10*detec_pixel_size long extended source blur_patt = np.zeros_like(coh_inter_patt(x)) for i in range(int(extension)): blur_patt += coh_inter_patt_sft(x,detec_pixel_size,sft=i) blur_patt += coh_inter_patt_sft(x,detec_pixel_size,sft=-i) #inter_pat = inter_pat/np.linalg.norm(inter_pat) blur_patt = blur_patt/np.linalg.norm(blur_patt) return blur_patt test_y = Blurred_fringes(x, 20) plt.plot(X,test_y) plt.show() model = Model(Blurred_fringes) params1 = model.make_params(N=35) results1 = model.fit(counts, params1, x=pixels)
修改后运行出现**ValueError: The model function generated NaN values and the fit aborted! Please check your model function and/or set boundaries on parameters where applicable.**报错,同时伴随以下运行警告:
:30: RuntimeWarning: invalid value encountered in divide sinq_part = ((np.sin((np.pi*x*a)/(lam*L))/((np.pi*x*a)/(lam*L)))**2) :41: RuntimeWarning: invalid value encountered in divide sinq_part = ((np.sin((np.pi*(x+(sft*detec_pixel_size))*a)/(lam*L))/((np.pi*(x+(sft*detec_pixel_size))*a)/(lam*L)))**2)
我知道sinc函数容易出现除零错误,尝试用sinq_part[sinq_part == np.nan] = 1处理,但问题仍未解决。
最终目标
先让函数完成拟合,再添加迭代步骤,通过查看一系列N值的置信水平找到最优拟合结果。N是唯一的拟合参数,其他提前定义的变量均为固定实验参数(可能因误差略有变动)。恳请提供解决思路。
内容的提问来源于stack exchange,提问作者Baba Booey
相关产品推荐
相关产品推荐

