scipy.optimize 误差幂函数答案

【问题标题】：scipy.optimize error power functionscipy.optimize 误差幂函数
【发布时间】：2017-11-24 00:16:11
【问题描述】：

我很难用给定的函数来拟合一组值： f(x)= const*(1-(x/a)**b)**c

我正在使用带有以下代码的 python 3.6.3：

import numpy as np
import scipy.optimize as opt
from scipy.optimize import curve_fit

x=[0.,0.4,0.8,1.6,2.,2.4]
y=[0.09882902,0.07298427,0.05111438,0.01679405,0.00517385,0.00065633]

def func(x,a,b,c):
    return y[0] * ( 1 - (x/a)**b )**c

x0=np.array([2.0,0.9,1.5])
opt.curve_fit(func,x,y,p0=x0)

我收到以下错误消息：

RuntimeWarning: invalid value encountered in power
  return y[0] * ( 1 - (x/a)**b )**c
///: OptimizeWarning: Covariance of the parameters could not be estimated
  category=OptimizeWarning)

问题似乎是因为最后一个指数，因为以下函数可以正常工作：

def func(x,a,b,c):
    return y[0] * ( 1 - (x/a)**b )*c

【问题讨论】：

标签： python scipy data-fitting

【解决方案1】：

只需执行经典的调试步骤：打印您的组件以检查数学运算是否定义明确。

修改你的函数为：

def func(x,a,b,c):
    print(x/a)
    print((x/a)**b)
    print((1-(x/a)**b))

    result = y[0] * ( 1 - (x/a)**b )**c
    print(result)

    if not np.isfinite(result):
        assert False
    else:
        return result

你会看到：

[ 0.   0.2  0.4  0.8  1.   1.2]
[ 0.          0.23492379  0.43838329  0.81805215  1.          1.17831965]
[ 1.          0.76507621  0.56161671  0.18194785  0.         -0.17831965]
...-py:13: RuntimeWarning: invalid value encountered in power
  result = y[0] * ( 1 - (x/a)**b )**c
[ 0.09882902  0.06613655  0.04159532  0.00767017  0.                 nan]

对于该负值，求幂不起作用并且引入了nan（可能会在后面的步骤中引入更多）。

备注：初始点默认为全1（也便于打印）。

解决此问题取决于您实际尝试执行的操作（更改模型；使用边界，...）。

【讨论】：

谢谢。实际上，这个函数在 gnuplot 拟合中工作得很好，我应该使用它，虽然它对于 x -> a 变得无限小如果我没记错的话，“NaN”似乎是一个低于 1e-13 的数字。有没有办法克服这个限制？
NaN 是 不是数字 :-) 同样：这完全取决于您可以访问的许多模型决策。如果 gnuplot 有效，那么假设或数据会有所不同，因为数学问题不会消失（对于一般情况）。
可能是 this is relevant to you（可能与 gnuplot 有所不同）。
嗯，我知道nan是什么意思。如果我运行曲线拟合，从 [a=2.4,b=0.9,c=1.5] 开始，func(x) 的值在等于 1e-13 的步骤之后变为“不是数字”。那么也许位数变得太重要了？
成为 nan 与门槛无关。

【解决方案2】：

我找到了解决方案，谢谢 Sascha。当 x > a 时，函数确实是不确定的。为了克服这个问题，我将 a 的值限制为始终优于 x 的最大值：

fit=opt.curve_fit(func,x,y,p0=x0,bounds=([x[-1]+0.001,0.0,1.00],[5.0,1.0,2.0]))

但是，在 gnuplot 中，拟合可能是找到 x > a。我不知道为什么。也许它占用了 func(x>a) 的实际部分，但我不知道这样做是否真的很简单。

【讨论】：