数据的指数拟合 (python)
Exponential fit of the data (python)
您好，我正在尝试用多项式或指数函数拟合我的数据，但我都失败了。我使用的代码如下：
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 的回答下关于错误或误解的评论）
之所以 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 中数据的第 2 列。它与第一列无关，尽管没有什么可以阻止您绘制第二列与第一列的修改版本当然，这就是您在绘制时所做的：
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 的参数未被使用。传递给拟合算法总是一个坏主意。摆脱它。
技巧来了：您根据第 2 列使用指数模型制作了 f1。因此，您应该传递 curve_fit 第二列作为自变量（在 function's doc-string 中标记为 xdata），然后将 f1 作为因变量。像这样：
popt1, pcov = curve_fit(model, argon[:,1], f1)
popt2, pcov = curve_fit(model, cupper[:,1], f2)
这将非常有效。
现在，当您想要绘制 2 个图形乘积的平滑版本时，您应该从自变量中的 common 区间开始。对你来说，这就是光子能量。两个数据文件中的第二列取决于：有一个函数（一个用于氩，另一个用于铜）将 μ/ρ 与光子能量相关联。因此，如果您有很多能量数据点，并且您设法获得了这些功能，那么您将有很多 μ/ρ 的数据点。由于这些函数是未知的，所以我能做的最好的事情就是简单地插值。但是数据是对数的，所以需要对数插值，而不是默认的线性。
所以现在，继续获取光子能量的大量数据点。在数据集中，能量点呈指数增长，因此您可以使用 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)')
好了。总结一下：
- 生成足够数量的自变量数据点以获得平滑的结果
- 插入关系
photon energy -> μ/ρ
- 将您的函数映射到插值
μ/ρ
您好，我正在尝试用多项式或指数函数拟合我的数据，但我都失败了。我使用的代码如下：
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()
这给了我
#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
我最终得到的是（红色：铜，蓝色虚线：氩，蓝色：乘法）
之所以 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 中数据的第 2 列。它与第一列无关，尽管没有什么可以阻止您绘制第二列与第一列的修改版本当然，这就是您在绘制时所做的：
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 的参数未被使用。传递给拟合算法总是一个坏主意。摆脱它。
技巧来了：您根据第 2 列使用指数模型制作了 f1。因此，您应该传递 curve_fit 第二列作为自变量（在 function's doc-string 中标记为 xdata），然后将 f1 作为因变量。像这样：
popt1, pcov = curve_fit(model, argon[:,1], f1)
popt2, pcov = curve_fit(model, cupper[:,1], f2)
这将非常有效。
现在，当您想要绘制 2 个图形乘积的平滑版本时，您应该从自变量中的 common 区间开始。对你来说，这就是光子能量。两个数据文件中的第二列取决于：有一个函数（一个用于氩，另一个用于铜）将 μ/ρ 与光子能量相关联。因此，如果您有很多能量数据点，并且您设法获得了这些功能，那么您将有很多 μ/ρ 的数据点。由于这些函数是未知的，所以我能做的最好的事情就是简单地插值。但是数据是对数的，所以需要对数插值，而不是默认的线性。
所以现在，继续获取光子能量的大量数据点。在数据集中，能量点呈指数增长，因此您可以使用 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)')
好了。总结一下：
- 生成足够数量的自变量数据点以获得平滑的结果
- 插入关系
photon energy -> μ/ρ - 将您的函数映射到插值
μ/ρ