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。你真的是这个意思吗? ab 将完全相关。

此外,您可能会发现 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()