指数拟合数据(python)

Question

嗨,我正在尝试使用多项式或指数函数来拟合我的数据,这两种方法都失败了.我正在使用的代码如下:

with open('argon.dat','r') as f:
    argon=f.readlines()

eng1 = np.array([float(argon[argon.index(i)].split('\n')[0].split('  ')[0])*1000 for i in argon])
II01 = np.array([1-math.exp(-float(argon[argon.index(i)].split('\n')[0].split('  ')[1])*(1.784e-3*6.35)) for i in argon])

with open('copper.dat','r') as f:
    copper=f.readlines()

eng2 = [float(copper[copper.index(i)].split('\n')[0].split('  ')[0])*1000 for i in copper]
II02 = [math.exp(-float(copper[copper.index(i)].split('\n')[0].split('  ')[1])*(8.128e-2*8.96)) for i in copper]

fig, ax1 = plt.subplots(figsize=(12,10))
ax2 = ax1.twinx()
ax1.set_yscale('log')
ax2.set_yscale('log')

arg = ax2.plot(eng1, II01, 'b--', label='Argon gas absorption at STP (6.35 cm)')
cop = ax1.plot(eng2, II02, 'r', label='Copper wall transp. (0.81 mm)')
plot = arg+cop

labs = [l.get_label() for l in plot]
ax1.legend(plot,labs,loc='lower right', fontsize=14)
ax1.set_ylim(1e-6,1)
ax2.set_ylim(1e-6,1)
ax1.set_xlim(0,160)
ax1.set_ylabel(r'$\displaystyle I/I_0$', fontsize=18)
ax2.set_ylabel(r'$\displaystyle 1-I/I_0$', fontsize=18)
ax1.set_xlabel('Photon Energy [keV]', fontsize=18)
plt.show()

哪能给我 我想要做的就是这样的他们适合的指数曲线和繁殖这些曲线与检测效率,最终绘制数据,而不是(我试图通过元素乘元素,但我没有足够的数据点有一个平滑的曲线)我试图使用polyfit并尝试定义一个指数函数来看它的工作情况但是我在两种情况下都得到了一条线

#def func(x, a, c, d):
#    return a*np.exp(-c*x)+d
#    
#popt, pcov = curve_fit(func, eng1, II01)
#plt.plot(eng1, func(eng1, *popt), label="Fitted Curve")

和

model = np.polyfit(eng1, II01 ,5) 
y = np.poly1d(model)
#splineYs = np.exp(np.polyval(model,eng1)) # also tried this but didnt work
ax2.plot(eng1,y)

如有需要,可从http://www.nist.gov/pml/data/xraycoef/index.cfm获取数据 .图3中也可以找到类似的工作:http://scitation.aip.org/content/ AAPT /日记/ AJP/83/8/10.1119/1.4923022

在@Oliver的回答之后,休息是偏见的:

我使用现有数据进行乘法运算:

i = 0
eff1 = []
while i < len(arg):
    eff1.append(arg[i]*cop[i])
    i += 1 

我最终得到的是(红色:铜色,蓝色虚线:氩气,蓝色:乘法) 这是我想拿到而是通过使用功能曲线,这将是我想结束了平滑曲线(评论已根据@奥利弗的回答提出关于什么是错的或者误解)

Answer 1

原因curvefit是给你一个常数(一条平线),是因为你使用你定义的模型传递了一个不相关的数据集!

让我先重新创建你的设置:

argon = np.genfromtxt('argon.dat')
copper = np.genfromtxt('copper.dat')

f1 = 1 - np.exp(-argon[:,1] * 1.784e-3 * 6.35)
f2 = np.exp(-copper[:,1] * 8.128e-2 * 8.96)

现在请注意,f1它基于文件中数据的第二列argon.dat.它与第一列无关,虽然没有什么可以阻止你绘制第二列的修改版本与第一列的比较,这就是你绘制时所做的:

import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

plt.semilogy(copper[:,0]*1000, f2, 'r-')  # <- f2 was not based on the first column of that file, but on the 2nd. Nothing stops you from plotting those together though...
plt.semilogy(argon[:,0]*1000, f1, 'b--')
plt.ylim(1e-6,1)
plt.xlim(0, 160)

def model(x, a, b, offset):
    return a*np.exp(-b*x) + offset

备注:在您的模型中,您有一个名为b未使用的参数.传递给拟合算法总是一个坏主意.摆脱它.

现在,这就是诀窍:你是f1基于第二列,使用指数模型制作的.所以你应该传递curve_fit第二列作为自变量(xdata在函数的doc-string中标记),然后f1作为因变量.像这样:

popt1, pcov = curve_fit(model, argon[:,1], f1)
popt2, pcov = curve_fit(model, cupper[:,1], f2)

这将非常有效.

现在,当您想要绘制2个图形的乘积的平滑版本时,您应该从自变量中的公共区间开始.对你而言,这是光子能量.两个数据文件中的第二列取决于:有一个函数(一个用于氩,另一个用于铜)?/?与光子能量相关.所以,如果你有很多能量数据点,并且你设法获得了这些功能,那么你将拥有许多数据点?/?.虽然这些功能是未知的,但我能做的最好的事情就是简单地进行插值.但是,数据是对数的,因此需要对数插值,而不是默认的线性.

所以现在,继续获得大量的光子能量数据点.在数据集中,能量点呈指数增长,因此您可以使用以下方法创建一组不错的新点np.logspace:

indep_var = argon[:,0]*1000
energy = np.logspace(np.log10(indep_var.min()),
                     np.log10(indep_var.max()),
                     512)  # both argon and cupper have the same min and max listed in the "energy" column.

它的优势在于两个数据集中的能量具有相同的最小值和最大值.否则,您将不得不减少此日志空间的范围.

接下来,我们(对数)插值关系energy -> ?/?:

interpolated_mu_rho_argon = np.power(10, np.interp(np.log10(energy), np.log10(indep_var), np.log10(argon[:,1]))) # perform logarithmic interpolation
interpolated_mu_rho_copper = np.power(10, np.interp(np.log10(energy), np.log10(copper[:,0]*1000), np.log10(copper[:,1])))

这是对刚刚完成的工作的直观表示:

f, ax = plt.subplots(1,2, sharex=True, sharey=True)
ax[0].semilogy(energy, interpolated_mu_rho_argon, 'gs-', lw=1)
ax[0].semilogy(indep_var, argon[:,1], 'bo--', lw=1, ms=10)
ax[1].semilogy(energy, interpolated_mu_rho_copper, 'gs-', lw=1)
ax[1].semilogy(copper[:,0]*1000, copper[:,1], 'bo--', lw=1, ms=10)
ax[0].set_title('argon')
ax[1].set_title('copper')
ax[0].set_xlabel('energy (keV)')
ax[0].set_ylabel(r'$\mu/\rho$ (cm²/g)')

标有蓝点的原始数据集已经过精细插值.

现在,最后的步骤变得简单了.因为已经找到了映射?/?到某个指数变量(我已重命名为f1和的函数)的模型参数f2,所以它们可用于生成存在的数据的平滑版本,以及两者的乘积这些功能:

plt.figure()
plt.semilogy(energy, model(interpolated_mu_rho_argon, *popt1), 'b-', lw=1)
plt.semilogy(argon[:,0]*1000, f1, 'bo ')

plt.semilogy(copper[:,0]*1000, f2, 'ro ',)
plt.semilogy(energy, model(interpolated_mu_rho_copper, *popt2), 'r-', lw=1) # same remark here!

argon_copper_prod = model(interpolated_mu_rho_argon, *popt1)*model(interpolated_mu_rho_copper, *popt2)
plt.semilogy(energy, argon_copper_prod, 'g-')

plt.ylim(1e-6,1)
plt.xlim(0, 160)
plt.xlabel('energy (keV)')
plt.ylabel(r'$\mu/\rho$ (cm²/g)')

你去吧 总结一下:

1. 生成足够数量的自变量数据点以获得平滑的结果
2. 插入关系 photon energy -> ?/?
3. 将您的函数映射到插值 ?/?