可以使用scipy minimize拟合曲线，但无法使用scipy curve_fit拟合

Question

我正在尝试使用scipy curve_fit将函数y= 1-a(1-bx)**n拟合到一些实验数据中。该模型仅适用于y>0，因此我剪切计算出的值以强制执行此条件。下面是代码：

import numpy as np
import scipy.optimize
import matplotlib.pyplot as plt

# Driver function for scipy.minimize

def driver_func(x, xobs, yobs):

    # Evaluate the fit function with the current parameter estimates

    ynew = myfunc(xobs, *x)
    yerr = np.sum((ynew - yobs) ** 2)

    return yerr

# Define function

def myfunc(x, a, b, n):

    y = 1.0 - a * np.power(1.0 - b * x, n) 
    y = np.clip(y, 0.00, None )

    return y

if __name__ == "__main__":

    # Initialise data

    yobs = np.array([0.005, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.004, 
                    0.048, 0.119, 0.199, 0.277, 0.346, 0.395, 0.444, 0.469, 
                    0.502, 0.527, 0.553, 0.582, 0.595, 0.603, 0.612, 0.599])
    xobs = np.array([0.013, 0.088, 0.159, 0.230, 0.292, 0.362, 0.419, 0.471,
                    0.528, 0.585, 0.639, 0.687, 0.726, 0.772, 0.814, 0.854,
                    0.889, 0.924, 0.958, 0.989, 1.015, 1.045, 1.076, 1.078])

    # Initial guess

    p0 = [2.0, 0.5, 2.0]

    # Check fit pre-regression

    yold = myfunc(xobs, *p0)
    plt.plot(xobs, yobs, 'ko', label='data', fillstyle='none')
    plt.plot(xobs, yold, 'g-', label='pre-fit: a=%4.2f, b=%4.2f, n=%4.2f' % tuple(p0))

    # Fit curve using SCIPY CURVE_FIT

    try:
        popt, pcov = scipy.optimize.curve_fit(myfunc, xobs, yobs, p0=p0)
    except:
        print("Could not fit data using SCIPY curve_fit")
    else:
        ynew = myfunc(xobs, *popt)
        plt.plot(xobs, ynew, 'r-', label='post-curve_fit: a=%4.2f, b=%4.2f, n=%4.2f' % tuple(popt))

    # Fit curve using SCIPY MINIMIZE

    res = scipy.optimize.minimize(driver_func, p0, args=(xobs, yobs), method='Nelder-Mead')
    ynw2 = myfunc(xobs, *res.x)
    plt.plot(xobs, ynw2, 'y-', label='post-minimize: a=%4.2f, b=%4.2f, n=%4.2f' % tuple(res.x))

    plt.legend()
    plt.show()

我还使用了 SCIPY MINIMIZE 来实现同样的功能。如下图所示，MINIMIZE 能够正常工作，但 CURVE_FIT 基本上用完了所有的评估次数并放弃了，尽管起始猜测与 MINIMIZE 解决方案（至少从视觉上）并不相差太远。希望能得到任何关于为什么 curve_fit 在这里似乎无法工作的想法。

谢谢！

更新： 根据mikuszefski的评论，我做了以下调整： 1. 从fit函数中删除了剪辑。

def myfunc_noclip(x, a, b, n):
    y = 1.0 - a * np.power(1.0 - b * x, n) 
    return y

2. introduced clipped arrays by removing data below a threshold

   ymin = 0.01
   xclp = xobs[np.where(yobs >= ymin)]
   yclp = yobs[np.where(yobs >= ymin)]
   

3. improved the initial guess (again visually)

   p0 = [1.75, 0.5, 2.0]
   

4. updated the call to curve_fit

   popt, pcov = scipy.optimize.curve_fit(myfunc_noclip, xclp, yclp, p0=p0)
   

但这似乎没有帮助，如下图所示：

其他Stack Overflow上的帖子似乎表明，scipy curve_fit在拟合参数中有指数时存在困难，例如 SciPy curve_fit not working when one of the parameters to fit is a power 所以我猜我也遇到了同样的问题。不确定如何解决...

- PetGriffin

Answer 1

这个问题是由函数定义中的剪裁引起的。两种最小化方法的工作原理根本不同，因此对这种剪裁的反应也非常不同。在这里，minimize 使用了 Nelder-Mead，这是一种无梯度方法。因此，该算法不计算数值梯度，也不估计任何雅可比矩阵。相比之下，least-squares（最小二乘法）确实这样做。然而，如果函数不连续，近似梯度和从此得出的雅可比矩阵就有些可疑了。如前所述，这种不连续性是由 np.clip 引入的。当它被移除后，可以很容易地看到，与剪辑一起使用时，P0 的猜测并不像它看起来那么好。尽管如此，curve_fit 在增加 maxfev=5000 后仍能收敛，而将方法更改为 method='CG' 时，minimize 立即失败。为了了解算法的困难之处，可以尝试手动提供 jac。

一些注意事项： 1）关于裁剪，删除被裁剪的数据可能是一个好主意，这样可以避免相关问题。 2）查看协方差矩阵，n的误差和与其他值的相关性非常高。

所以这是我得到的内容。

import numpy as np
import scipy.optimize
import matplotlib.pyplot as plt

# Driver function for scipy.minimize
def driver_func( x, xobs, yobs ):
    # Evaluate the fit function with the current parameter estimates
    ynew = myfunc( xobs, *x)
    yerr = np.sum( ( ynew - yobs ) ** 2 )
    return yerr

# Define functions
def myfunc( x, a, b, n ):
    y = 1.0 - a * np.power( 1.0 - b * x, n ) 
    y = np.clip( y, 0.00, None )
    return y

def myfunc_noclip( x, a, b, n ):
    y = 1.0 - a * np.power( 1.0 - b * x, n ) 
    return y

if __name__ == "__main__":

    # Initialise data
    yobs = np.array([
        0.005, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.004, 
        0.048, 0.119, 0.199, 0.277, 0.346, 0.395, 0.444, 0.469, 
        0.502, 0.527, 0.553, 0.582, 0.595, 0.603, 0.612, 0.599
    ])
    xobs = np.array([
        0.013, 0.088, 0.159, 0.230, 0.292, 0.362, 0.419, 0.471,
        0.528, 0.585, 0.639, 0.687, 0.726, 0.772, 0.814, 0.854,
        0.889, 0.924, 0.958, 0.989, 1.015, 1.045, 1.076, 1.078
    ])

    # Clipped data
    ymin = 0.01
    xclp = xobs[ np.where( yobs >= ymin ) ]
    yclp = yobs[ np.where( yobs >= ymin ) ]

    # Initial guess
    p0 = [ 2.0, 0.5, 2.0 ]

    # Check fit pre-regression
    yold = myfunc( xobs, *p0 )
    plt.plot( xobs, yobs, 'ko', label='data', fillstyle='none' )
    plt.plot( xobs, yold, 'g-', label='pre-fit: a=%4.2f, b=%4.2f, n=%4.2f' % tuple( p0 ) )

    # Fit curve using SCIPY CURVE_FIT
    try:
        popt, pcov = scipy.optimize.curve_fit( myfunc, xobs, yobs, p0=p0, maxfev=5000 )
    except:
        print("Could not fit data using SCIPY curve_fit")
    else:
        ynew = myfunc( xobs, *popt )
        plt.plot( xobs, ynew, 'r-', label="curve-fit: a=%4.2f, b=%4.2e, n=%4.2f" % tuple( popt ) )

    # Fit curve using SCIPY CURVE_FIT on clipped data
    p0 = [ 1.75, 1e-4, 1e3 ]
    try:
        popt, pcov = scipy.optimize.curve_fit( myfunc_noclip, xclp, yclp, p0=p0 )
    except:
        print("Could not fit data using SCIPY curve_fit")
    else:
        ynew = myfunc_noclip( xobs, *popt )
        plt.plot( xobs, ynew, 'k-', label="curve-fit clipped data: a=%4.2f, b=%4.2e, n=%4.2f" % tuple( popt ) )

    # Fit curve using SCIPY MINIMIZE
    p0 = [ 2.0, 0.5, 2.0 ]
    res = scipy.optimize.minimize( driver_func, p0, args=( xobs, yobs ), method='Nelder-Mead' )
    # ~res = scipy.optimize.minimize(driver_func, p0, args=(xobs, yobs), method='CG')
    ynw2 = myfunc( xobs, *res.x )
    plt.plot( xobs, ynw2, 'y--', label='Nelder-Mead 1: a=%4.2f, b=%4.2f, n=%4.2f' % tuple( res.x ) )
    p0 = [ 2.4, 3.6e-4, 5.6e3 ]
    res = scipy.optimize.minimize( driver_func, p0, args=( xobs, yobs ), method='Nelder-Mead' )
    # ~res = scipy.optimize.minimize(driver_func, p0, args=(xobs, yobs), method='CG')
    ynw2 = myfunc( xobs, *res.x )
    plt.plot( xobs, ynw2, 'b:', label='Nelder-Mead 2: a=%4.2f, b=%4.2e, n=%4.2e' % tuple( res.x ) )

    plt.legend( loc=2 )
    plt.ylim( -0.05, 0.7 )
    plt.grid()
    plt.show()

所以我会说它运行得还可以。不过，我收到了一个溢出警告。

- mikuszefski

---

嗨mikuszefski，谢谢你的评论。我已经更新了帖子，包括你在（1）中建议的更改结果，但不幸的是它们并没有起到作用。不确定如何解决（2）-这个方程是通用的，在其他数据集上似乎工作正常，所以我猜我只是在示例数据中运气不好。我在其他一些（x，y）数据集上也遇到了同样的问题。 - PetGriffin

@PetGriffin 已经制作并根据您的编辑进行了调整。在这里运行良好。 - mikuszefski

嗨@mikuszefski - 感谢你花时间进一步调查。我运行了你的代码版本，它工作得很好，但是如果我注释掉语句p0 = [1.75, 1e-4, 1e3]，那么它就不工作了。可以使用数据的其他属性来估计初始参数，并且这些起始值（不幸地）是非物理的。所以问题似乎是我的数据和模型，而不是曲线拟合。 - PetGriffin

@PetGriffin 我认为问题在于该模型中参数之间的极高相关性。卡方最小值可能在参数空间中以某种对角线方式被拉伸得非常长。因此，数据的微小变化可能导致拟合结果发生重大变化。从这个意义上说，这是一个模型问题。你可以看到初始斜率与 b * n 相关，所以你可以轻松地在它们之间进行转换。根据问题背后的物理学，你可以尝试不同的饱和模型。 - mikuszefski