Scipy.Curve 拟合指数函数答案

【问题标题】：Scipy.Curve fit struggling with exponential functionScipy.Curve 拟合指数函数
【发布时间】：2021-03-20 09:09:28
【问题描述】：

我正在尝试拟合方程的曲线：

y = ( (np.exp(-k2*(t+A))) - ((k1/v)*Co) )/ -k2

where  A = (-np.log((k1/v)*Co))/k2

主管给了我一个数据集，该数据集看起来像一个粗略的指数，在其顶部变平为一条水平直线。当我拟合方程时，我只收到曲线拟合的直线和相应的警告：

<ipython-input-24-7e57039f2862>:36: RuntimeWarning: overflow encountered in exp
 return ( (np.exp(-k2*(t+A))) - ((k1/v)*Co) )/ -k2

我使用的代码如下：

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy.optimize import differential_evolution

xData_forfit = [1.07683e+13, 1.16162e+13, 1.24611e+13, 1.31921e+13, 1.40400e+13, 2.65830e+13,
 2.79396e+13, 2.86676e+13, 2.95155e+13, 3.03605e+13, 3.12055e+13, 3.20534e+13,
 3.27814e+13, 3.36293e+13, 3.44772e+13, 3.53251e+13, 3.61730e+13, 3.77459e+13,
 3.85909e+13, 3.94388e+13, 4.02838e+13, 4.11317e+13, 4.19767e+13, 4.27076e+13,
 5.52477e+13, 5.64143e+13, 5.72622e+13, 5.81071e+13, 5.89550e+13, 5.98000e+13,
 6.05280e+13, 6.13759e+13, 6.22209e+13, 6.30658e+13, 6.39137e+13, 6.46418e+13,
 6.55101e+13, 6.63551e+13, 6.72030e+13, 6.80480e+13, 6.88929e+13, 6.97408e+13,
 7.04688e+13, 7.13167e+13, 7.21617e+13, 8.50497e+13, 8.58947e+13, 8.67426e+13,
 8.75876e+13, 8.83185e+13, 9.00114e+13, 9.08563e+13, 9.17013e+13]
yData_forfit = [1375.409524, 1378.095238, 1412.552381, 1382.904762, 1495.2, 1352.4,
 1907.971429, 1953.52381,  1857.352381, 1873.990476, 1925.114286, 1957.085714,
 2030.52381,  1989.8,      2042.733333, 2060.095238, 2134.361905, 2200.742857,
 2342.72381,  2456.047619, 2604.542857, 2707.971429 ,2759.87619,  2880.52381,
 3009.590476, 3118.771429, 3051.52381,  3019.771429, 3003.561905, 3083.0,
 3082.885714, 2799.866667, 3012.419048, 3013.266667, 3106.714286, 3090.47619,
 3216.638095, 3108.447619, 3199.304762, 3154.257143, 3112.419048, 3284.066667,
 3185.942857, 3157.380952, 3158.47619,  3464.257143, 3434.67619,  3291.457143,
 2851.371429, 3251.904762, 3056.152381, 3455.07619,  3386.942857]

def fnct_to_opt(t, k2, k1):
    #EXPERIMENTAL CONSTANTS
    v = 105
    Co = 1500

    A = (-np.log((k1/v)*Co))/k2
    return ( (np.exp(-k2*(t+A))) - ((k1/v)*Co) )/ -k2


initial_k2k1 = [100, 1*10**-3]
constants = curve_fit(fnct_to_opt, xData_forfit, yData_forfit, p0=initial_k2k1)
k2_fit = constants[0][0]
k1_fit = constants[0][1]


fit = []
for i in xData_forfit:
    fit.append(fnct_to_opt(i,k2_fit,k1_fit))


plt.plot(xData_forfit, yData_forfit, 'or', ms='2')
plt.plot(xData_forfit, fit)

这给了我这个情节：

据我所知，由于 np.exp 项的值太大，代码没有产生有用的输出，但我不知道如何诊断此溢出的来源或如何解决问题。任何帮助将不胜感激，谢谢。

【问题讨论】：

请提供预期的minimal, reproducible example (MRE)。我们应该能够复制和粘贴您的代码的连续块，执行该文件，并重现您的问题以及跟踪问题点的输出。这让我们可以根据您的测试数据和所需的输出来测试我们的建议。我们希望您执行基本诊断以包含在您的帖子中。至少，在错误点打印可疑值并追溯它们的来源。
由于几个未定义的符号，您发布的代码无法运行。
我的错，已添加到数据集中，因此现在可以运行。我不知道如何去诊断溢出到底发生在哪里，因为它被埋在不透明的curve_fit的某个地方。
您有特定的上下文吗？喜欢行业或领域

标签： python numpy scipy curve-fitting scipy-optimize

【解决方案1】：

溢出发生在错误消息告诉你的地方：在fnct_to_opt 的return 表达式中。我要求您在错误点之前打印有问题的值；这会告诉你问题。

在错误点，A 中的值在 e+13 到 e+14 的范围内。 t 微不足道； k2 略低于-10000.0

因此，np.exp 的参数中的值远远超出了函数可以处理的域。只需在您的函数中添加一行并观察结果：

def fnct_to_opt(t, k2, k1):
    #EXPERIMENTAL CONSTANTS
    v = 105
    Co = 1500

    A = (-np.log((k1/v)*Co))/k2
    print("TRACE", "\nk2", k2, "\nt", t, "\nA", A, "\nother", k1, v, Co)
    return ( (np.exp(-k2*(t+A))) - ((k1/v)*Co) )/ -k2

【讨论】：

我明白你现在关于打印的意思了。我认为我将不得不研究方程式本身 - 我假设有了这个溢出，我不会得到任何结果。也许这是一个理论错误。感谢您的帮助
对。我怀疑某处的方程式有错误。您使用 e^(10^15) 生成的数量对于任何实际用途来说都有点大。

【解决方案2】：

我认为问题可能出在优化函数上，从某种意义上说可能是一个错误。

例如：

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy.optimize import differential_evolution

xData_forfit = [1.07683e+13, 1.16162e+13, 1.24611e+13, 1.31921e+13, 1.40400e+13, 2.65830e+13,
 2.79396e+13, 2.86676e+13, 2.95155e+13, 3.03605e+13, 3.12055e+13, 3.20534e+13,
 3.27814e+13, 3.36293e+13, 3.44772e+13, 3.53251e+13, 3.61730e+13, 3.77459e+13,
 3.85909e+13, 3.94388e+13, 4.02838e+13, 4.11317e+13, 4.19767e+13, 4.27076e+13,
 5.52477e+13, 5.64143e+13, 5.72622e+13, 5.81071e+13, 5.89550e+13, 5.98000e+13,
 6.05280e+13, 6.13759e+13, 6.22209e+13, 6.30658e+13, 6.39137e+13, 6.46418e+13,
 6.55101e+13, 6.63551e+13, 6.72030e+13, 6.80480e+13, 6.88929e+13, 6.97408e+13,
 7.04688e+13, 7.13167e+13, 7.21617e+13, 8.50497e+13, 8.58947e+13, 8.67426e+13,
 8.75876e+13, 8.83185e+13, 9.00114e+13, 9.08563e+13, 9.17013e+13]
yData_forfit = [1375.409524, 1378.095238, 1412.552381, 1382.904762, 1495.2, 1352.4,
 1907.971429, 1953.52381,  1857.352381, 1873.990476, 1925.114286, 1957.085714,
 2030.52381,  1989.8,      2042.733333, 2060.095238, 2134.361905, 2200.742857,
 2342.72381,  2456.047619, 2604.542857, 2707.971429 ,2759.87619,  2880.52381,
 3009.590476, 3118.771429, 3051.52381,  3019.771429, 3003.561905, 3083.0,
 3082.885714, 2799.866667, 3012.419048, 3013.266667, 3106.714286, 3090.47619,
 3216.638095, 3108.447619, 3199.304762, 3154.257143, 3112.419048, 3284.066667,
 3185.942857, 3157.380952, 3158.47619,  3464.257143, 3434.67619,  3291.457143,
 2851.371429, 3251.904762, 3056.152381, 3455.07619,  3386.942857]

def fnct_to_opt(t, k2, k1):
    #EXPERIMENTAL CONSTANTS
    v = 105
    Co = 1500

    #A = (-np.log((k1/v)*Co))/k2
    #return ( (np.exp(-k2*(t+A))) - ((k1/v)*Co) )/ -k2
    #A = (np.log((k1/v)*Co))/k2
    return k2/np.log(t) + k1


initial_k2k1 = [10, 1]
constants = curve_fit(fnct_to_opt, xData_forfit, yData_forfit, p0=initial_k2k1)
k2_fit = constants[0][0]
k1_fit = constants[0][1]
#v_fit = constants[0][2]
#Co_fit = constants[0][3]

fit = []
for i in xData_forfit:
    fit.append(fnct_to_opt(i,k2_fit,k1_fit))


plt.plot(xData_forfit, yData_forfit, 'or', ms='2')
plt.plot(xData_forfit, fit)

所以我放置了一个更简单的函数，但背后有一些更清晰的直觉。例如，在原版中，我认为这些符号和指数根本不会实现形状。但是在我看来，指数放错了位置，所以我将其更改为对数。添加一个常数和一个比例参数。我建议仔细检查原始功能。推导可能存在问题。我不认为这是一个计算问题。

这更接近预期。

【讨论】：

谢谢，我认为这一定是我试图拟合的方程的误解。