numpy exp 旧代码中的新错误
new error in old code in numpy exp
最近我正在处理一些数据,在保存绘图后我能够使用 curve_fit 获得曲线,并且我 return 编辑了相同的代码之后才发现获得的值它不起作用。
#! python 3.5.2
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats
from scipy.optimize import curve_fit
data= np.array([
[24, 0.176644513],
[27, 0.146382841],
[30, 0.129891534],
[33, 0.105370908],
[38, 0.077820511],
[50, 0.047407538]])
x, y = np.array([]), np.array([])
for val in data:
x = np.append(x, val[0])
y = np.append(y, (val[1]/(1-val[1])))
def f(x, a, b):
return (np.exp(-a*x)**b)
# The original a and b values obtained
a = -0.2 # after rounding
b = -0.32 # after rounding
plt.scatter(x, y)
Xcurve = np.linspace(x[0], x[-1], 500)
plt.plot(Xcurve, f(Xcurve,a,b), ls='--', color='k', lw=1)
plt.show()
# the original code to get the values
a = b = 1
popt, pcov = curve_fit(f, x, y, (a, b))
而之前 curve_fit return 编辑了值 a、b = -0.2、-0.32 现在 returns:
警告(来自警告模块):
文件“C:/Users ... 第 22 行
return (np.exp(-a*x)**b)
RuntimeWarning:在 exp
中遇到溢出
据我所知代码没有改变。谢谢
在不知道代码发生了什么变化的情况下,很难说出 "working" 和 "not working" 之间发生了什么变化。可能是您使用的 scipy 版本的变化给出了不同的结果:在过去几年中 curve_fit()
中的底层实现发生了变化。
而且:curve_fit()
(以及底层 python 和它使用的 Fortran 代码)需要对参数进行相当好的初始猜测才能解决许多问题。如果对参数的猜测不正确,许多问题都会失败。
指数衰减问题似乎对 Levenberg-Marquardt 算法(以及 curve_fit()
使用的实现)特别具有挑战性,并且确实需要合理的起点。也很容易进入参数的一部分 space 其中函数的计算结果为零,参数值的更改无效。
如果可能的话,如果您的问题涉及指数衰减,那么使用对数 space 会很有帮助。也就是说,模型 log(f)
,而不是 f
本身。特别是对于您的问题,您的模型函数是 exp(-a*x)**b
。你真的是这个意思吗? a
和 b
将完全相关。
此外,您可能会发现 lmfit
有帮助。它有一个 Model
class 用于曲线拟合,使用类似的底层代码,但允许修复或设置任何参数的界限。你的问题的一个例子是(大约):
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats
from scipy.optimize import curve_fit
import lmfit
data= np.array([
[24, 0.176644513],
[27, 0.146382841],
[30, 0.129891534],
[33, 0.105370908],
[38, 0.077820511],
[50, 0.047407538]])
x, y = np.array([]), np.array([])
for val in data:
x = np.append(x, val[0])
y = np.append(y, (val[1]/(1-val[1])))
def f(x, a, b):
print("In f: a, b = " , a, b)
return (np.exp(-a*x)**b)
fmod = lmfit.Model(f)
params = fmod.make_params(a=-0.2, b=-0.4)
# set bounds on parameters
params['a'].min = -2
params['a'].max = 0
params['b'].vary = False
out = fmod.fit(y, params, x=x)
print(out.fit_report())
plt.plot(x, y)
plt.plot(x, out.best_fit, '--')
plt.show()
最近我正在处理一些数据,在保存绘图后我能够使用 curve_fit 获得曲线,并且我 return 编辑了相同的代码之后才发现获得的值它不起作用。
#! python 3.5.2
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats
from scipy.optimize import curve_fit
data= np.array([
[24, 0.176644513],
[27, 0.146382841],
[30, 0.129891534],
[33, 0.105370908],
[38, 0.077820511],
[50, 0.047407538]])
x, y = np.array([]), np.array([])
for val in data:
x = np.append(x, val[0])
y = np.append(y, (val[1]/(1-val[1])))
def f(x, a, b):
return (np.exp(-a*x)**b)
# The original a and b values obtained
a = -0.2 # after rounding
b = -0.32 # after rounding
plt.scatter(x, y)
Xcurve = np.linspace(x[0], x[-1], 500)
plt.plot(Xcurve, f(Xcurve,a,b), ls='--', color='k', lw=1)
plt.show()
# the original code to get the values
a = b = 1
popt, pcov = curve_fit(f, x, y, (a, b))
而之前 curve_fit return 编辑了值 a、b = -0.2、-0.32 现在 returns: 警告(来自警告模块): 文件“C:/Users ... 第 22 行 return (np.exp(-a*x)**b) RuntimeWarning:在 exp
中遇到溢出据我所知代码没有改变。谢谢
在不知道代码发生了什么变化的情况下,很难说出 "working" 和 "not working" 之间发生了什么变化。可能是您使用的 scipy 版本的变化给出了不同的结果:在过去几年中 curve_fit()
中的底层实现发生了变化。
而且:curve_fit()
(以及底层 python 和它使用的 Fortran 代码)需要对参数进行相当好的初始猜测才能解决许多问题。如果对参数的猜测不正确,许多问题都会失败。
指数衰减问题似乎对 Levenberg-Marquardt 算法(以及 curve_fit()
使用的实现)特别具有挑战性,并且确实需要合理的起点。也很容易进入参数的一部分 space 其中函数的计算结果为零,参数值的更改无效。
如果可能的话,如果您的问题涉及指数衰减,那么使用对数 space 会很有帮助。也就是说,模型 log(f)
,而不是 f
本身。特别是对于您的问题,您的模型函数是 exp(-a*x)**b
。你真的是这个意思吗? a
和 b
将完全相关。
此外,您可能会发现 lmfit
有帮助。它有一个 Model
class 用于曲线拟合,使用类似的底层代码,但允许修复或设置任何参数的界限。你的问题的一个例子是(大约):
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats
from scipy.optimize import curve_fit
import lmfit
data= np.array([
[24, 0.176644513],
[27, 0.146382841],
[30, 0.129891534],
[33, 0.105370908],
[38, 0.077820511],
[50, 0.047407538]])
x, y = np.array([]), np.array([])
for val in data:
x = np.append(x, val[0])
y = np.append(y, (val[1]/(1-val[1])))
def f(x, a, b):
print("In f: a, b = " , a, b)
return (np.exp(-a*x)**b)
fmod = lmfit.Model(f)
params = fmod.make_params(a=-0.2, b=-0.4)
# set bounds on parameters
params['a'].min = -2
params['a'].max = 0
params['b'].vary = False
out = fmod.fit(y, params, x=x)
print(out.fit_report())
plt.plot(x, y)
plt.plot(x, out.best_fit, '--')
plt.show()