Scipy.Curve 适应指数函数

Scipy.Curve fit struggling with exponential function

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

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 项的值太大,代码没有产生有用的输出,但我不知道如何诊断溢出的位置来自或如何解决问题。任何帮助将不胜感激,谢谢。

溢出恰好发生在错误消息告诉您的位置:在 fnct_to_optreturn 表达式中。我要求您在错误点之前打印有问题的值;这会告诉你问题所在。

在错误点,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

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

例如:

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)

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

这更接近预期。