指数函数中单个变量的 Scipy 曲线答案

【问题标题】：Scipy Curvefit for a Single Variable in an Exponential Function指数函数中单个变量的 Scipy 曲线
【发布时间】：2020-04-21 14:06:42
【问题描述】：

所以我有一组 X 和 Y 数据，我从输出文件中读取它们然后截断。我将在下面给出一个示例，以防万一有人希望它作为可能测试我的问题的参考（抱歉，它太多了）

此数据显示时间和记录，并代表指数衰减类型的趋势。记录的所有数据都是相似的，但是每条记录都有一个变化的系数，所以我正在尝试开发一个代码来找出那个单一的系数是什么。我使用的等式是：

Y*((exp(-TMA*(log(2.)/HL110))) + (Xexp(-TMA(log(2.)/HL108))))) + b

这里改变的变量是 Y。其他的都是已知的。这是我想要拟合的变量（Y）。我在 Excel 中做了一些工作，可以说它处于 9000 年代的高位（这只是内存不足）。其他情况在 4000 和 7000 中。所以它范围很广，这就是为什么我需要一个代码来完成它，否则我每次都必须手动完成，而且我们有成千上万的记录需要分析。我写了一个代码，但它是扁平的，并没有真正提供合适的。我会在下面提供它。它还包含上面提到的所有常量，这些常量不会更改。

### Section 1 ###
from scipy import *
from matplotlib import pyplot
from scipy.optimize import minimize_scalar
from scipy.optimize import curve_fit
import numpy as np


### Section 2 ###
data = np.loadtxt('Ag - Near_7_2026.txt') ### LOAD FILE DATA HERE ###
data_trunc = data[25:len(data)] ### TRUNCATED DATA UP TO 104 SEC ###
TM = data_trunc[:,0] ### TIME MARK ###
TMA = TM + 4 ### CORRECTED TIME ARRAY, ELAPSED TIME ###
Counts = data_trunc[:,1]
Sigma = sqrt([Counts])


### DEFINE PROBLEM CONSTANTS ###
HL110 = 24.6 ### ENDF ACCEPTED HL ###
HL108 = 142.92 ### ENDF ACCEPTED HL ###
b = 1.3333
X = 0.02955 ### FROM MCNP MODEL ###


### Function Handel ###
def func(TMA, Y, X, HL110, HL108, b):
    return Y*((exp(-TMA*(log(2.)/HL110))) + (X*exp(-TMA*(log(2.)/HL108)))) + b ### MODEL FUNCTION ###

f = func(TMA, 5000, X, HL110, HL108, b) ### CALLABLE NEEDED FOR CURVE_FIT ###

# Data plotting ###
pyplot.plot(TMA, f, '.b', label = 'data')
pyplot.legend(fontsize = 'large')


### Curve Fitting and plotting ###
popt, pcov = curve_fit(func, TMA, f)
pyplot.plot(TMA, func(TMA, *popt), 'r-', label = 'fit')
pyplot.tick_params(labelsize='large')
pyplot.legend(fontsize='large')
pyplot.xlabel('Adjusted Time')
pyplot.ylabel('Counts')
pyplot.show()

我已尽力希望注释掉这里的大部分代码，以帮助任何帮助我理解什么是什么的人。当我这样做时，我使用了https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html 作为我的参考。（我与他们或任何事情无关，我只是想展示我的思维过程，以防万一它提供任何关于我哪里出错的额外信息）。

非常感谢任何帮助，我可以提供任何要求的澄清信息！

【问题讨论】：

标签： python python-3.x scipy curve-fitting scipy-optimize

【解决方案1】：

这似乎只需要很少的改动。我使用您的数据创建了一个文本文件，对于代码本身，我没有将常量传递给函数。

### Section 1 ###
from scipy import *
from matplotlib import pyplot
from scipy.optimize import minimize_scalar
from scipy.optimize import curve_fit
import numpy as np


### Section 2 ###
#data = np.loadtxt('Ag - Near_7_2026.txt') ### LOAD FILE DATA HERE ###
data = np.loadtxt('temp.dat')
data_trunc = data[25:len(data)] ### TRUNCATED DATA UP TO 104 SEC ###
TM = data_trunc[:,0] ### TIME MARK ###
TMA = TM + 4 ### CORRECTED TIME ARRAY, ELAPSED TIME ###
Counts = data_trunc[:,1]
Sigma = sqrt([Counts])


### DEFINE PROBLEM CONSTANTS ###
HL110 = 24.6 ### ENDF ACCEPTED HL ###
HL108 = 142.92 ### ENDF ACCEPTED HL ###
b = 1.3333
X = 0.02955 ### FROM MCNP MODEL ###


### Function Handel ###
def func(TMA, Y): # no need to pass constants
    return Y*((exp(-TMA*(log(2.)/HL110))) + (X*exp(-TMA*(log(2.)/HL108)))) + b ### MODEL FUNCTION ###

# # no need to pass constants
f = func(TMA, 5000) ### CALLABLE NEEDED FOR CURVE_FIT ###

# Data plotting ###
pyplot.plot(TMA, f, '.b', label = 'data')
pyplot.legend(fontsize = 'large')


### Curve Fitting and plotting ###
popt, pcov = curve_fit(func, TMA, f)
print('Fitted parameters:', popt)

pyplot.plot(TMA, func(TMA, *popt), 'r-', label = 'fit')
pyplot.tick_params(labelsize='large')
pyplot.legend(fontsize='large')
pyplot.xlabel('Adjusted Time')
pyplot.ylabel('Counts')
pyplot.show()

EDIT -- 求解 Y 的代码

### Section 1 ###
from scipy import *
from matplotlib import pyplot
from scipy.optimize import minimize_scalar
from scipy.optimize import curve_fit
import numpy as np


### Section 2 ###
#data = np.loadtxt('Ag - Near_7_2026.txt') ### LOAD FILE DATA HERE ###
data = np.loadtxt('temp.dat')
data_trunc = data[25:len(data)] ### TRUNCATED DATA UP TO 104 SEC ###
TM = data_trunc[:,0] ### TIME MARK ###
TMA = TM + 4 ### CORRECTED TIME ARRAY, ELAPSED TIME ###
Counts = data_trunc[:,1]
Sigma = sqrt(Counts)


### DEFINE PROBLEM CONSTANTS ###
HL110 = 24.6 ### ENDF ACCEPTED HL ###
HL108 = 142.92 ### ENDF ACCEPTED HL ###
b = 1.3333
X = 0.02955 ### FROM MCNP MODEL ###


### Function Handel ###
def func(TMA, Y): # no need to pass constants
    return Y*((exp(-TMA*(log(2.)/HL110))) + (X*exp(-TMA*(log(2.)/HL108)))) + b ### MODEL FUNCTION ###

# # no need to pass constants
#f = func(TMA, 5000) ### CALLABLE NEEDED FOR CURVE_FIT ###

# Data plotting ###
pyplot.plot(TMA, Counts, '.b', label = 'data')
pyplot.legend(fontsize = 'large')


### Curve Fitting and plotting ###
popt, pcov = curve_fit(func, TMA, Counts)
print('Fitted parameters:', popt)

pyplot.plot(TMA, func(TMA, *popt), 'r-', label = 'fit')
pyplot.tick_params(labelsize='large')
pyplot.legend(fontsize='large')
pyplot.xlabel('Adjusted Time')
pyplot.ylabel('Counts')
pyplot.show()

【讨论】：

我的问题是我们如何确定确切的 Y 值是什么？在代码中，我将“5000”作为占位符，但是如何找出最适合的 Y 呢？
我已经在标记为“编辑 - 解决 Y 的代码”的注释下添加了代码，用于执行此操作。我使用“计数”而不是“f”。顺便说一句，我的背景是核工程。
哦，太棒了！在什么方面？那么这里的问题是它给出的协方差相当高。如果我错了，请原谅我，但这需要最小化，对吗？当通过 Excel 工作表完成时，给定数据集的 Y 值约为 5890。 solution 给出的 Y 远大于 20,000，协方差远大于 6000。
编辑：可以看出，通过截断，建议的 Y 值甚至更低，大约一个数量级。
在我发布的代码中，如果将 curve_fit() 中的拟合值替换为“popt=[5890.0]”，那么（在技术术语中）情节很糟糕。数据、Python 方程或值 5890 不正确。