数据的指数拟合（Python）

Question

数据的指数拟合（Python）

4

你好，我正在尝试使用多项式或指数函数拟合我的数据，但两者都失败了。我使用的代码如下：

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/journal/ajp/83/8/10.1119/1.4923022 @Oliver回答后，其余部分已经进行了修改:

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

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

我得到的结果是（红色：铜，虚线蓝色：氩，蓝色：乘法）。

这是我希望得到的结果，但通过使用曲线函数，这将是一个平滑的曲线，而我想最终得到这个结果（在@oliver的答案下已经发表了评论，指出了什么是错误或误解）。

- jackaraz

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Oliver W. · Accepted Answer

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。因此，您应该将第二列作为自变量（在函数文档字符串中标记为xdata），然后将f1作为因变量传递给curve_fit。像这样：

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

这将非常完美地发挥作用。

现在，当您想绘制两个图形乘积的平滑版本时，您应该从独立变量的一个公共区间开始。对于您来说，这是光子能量。两个数据文件中的第二列都取决于此：有一个函数（分别为氩和铜），将μ/ρ与光子能量相关联。因此，如果您拥有大量能量的数据点，并且已经获得了这些函数，您将拥有许多μ/ρ的数据点。尽管这些函数是未知的，但我所能做的最好的事情只是简单地插值。然而，数据是对数的，因此需要对数插值，而不是默认的线性插值。

因此，现在通过获取许多光子能量数据点来继续。在数据集中，能量点呈指数增长，因此您可以使用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)')

总结一下：

生成足够数量的自变量数据点以获得平滑结果
插值得到关系式 光子能量 -> μ/ρ
将函数映射到插值后的 μ/ρ 上