可以用 scipy 最小化拟合曲线，但不能用 scipy curve_fit
Can fit curve with scipy minimize but not with scipy curve_fit
我正在尝试使用 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函数的clipping如下：
def myfunc_noclip(x, a, b, n):
y = 1.0 - a * np.power(1.0 - b * x, n)
return y
通过删除低于阈值的数据引入了裁剪数组
ymin = 0.01
xclp = xobs[np.where(yobs >= ymin)]
yclp = yobs[np.where(yobs >= ymin)]
改进了初始猜测（再次在视觉上）
p0 = [1.75, 0.5, 2.0]
将调用更新为 curve_fit
popt, pcov = scipy.optimize.curve_fit(myfunc_noclip, xclp, yclp, p0=p0)
但这似乎没有帮助，如下图所示：
Whosebug 上的其他帖子似乎表明 scipy curve_fit 在拟合曲线时遇到问题，其中一个拟合参数是指数，例如
SciPy curve_fit not working when one of the parameters to fit is a power
所以我猜我有同样的问题。虽然不确定如何解决它...
这个问题是由函数定义中的裁剪引起的。这两种最小化方法的工作原理根本不同，因此对这种裁剪的反应也大不相同。这里minimize和Nelder-Mead一起使用，是一种无梯度的方法。因此，该算法不计算数值梯度，也不估计任何雅可比矩阵。相比之下，最终由 curve_fit 调用的 least-squares 正是这样做的。但是，如果函数不是连续的，则逼近梯度并由此得出任何雅可比行列式都有些问题。如前所述，这种不连续性是由 np.clip 引入的。删除后，可以很容易地看到，P0 猜测并不像包含裁剪的情况下看起来那么好。 curve_fit 确实随着 maxfev=5000 的增加而收敛，而 minimize 在将方法更改为 method='CG' 时立即失败。要查看算法的困难，可以尝试手动提供 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()
所以我会说它工作得很好。不过，我收到溢出警告。
我正在尝试使用 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函数的clipping如下：
def myfunc_noclip(x, a, b, n):
y = 1.0 - a * np.power(1.0 - b * x, n)
return y
通过删除低于阈值的数据引入了裁剪数组
ymin = 0.01 xclp = xobs[np.where(yobs >= ymin)] yclp = yobs[np.where(yobs >= ymin)]改进了初始猜测（再次在视觉上）
p0 = [1.75, 0.5, 2.0]将调用更新为 curve_fit
popt, pcov = scipy.optimize.curve_fit(myfunc_noclip, xclp, yclp, p0=p0)
但这似乎没有帮助，如下图所示：
Whosebug 上的其他帖子似乎表明 scipy curve_fit 在拟合曲线时遇到问题，其中一个拟合参数是指数，例如
SciPy curve_fit not working when one of the parameters to fit is a power
所以我猜我有同样的问题。虽然不确定如何解决它...
这个问题是由函数定义中的裁剪引起的。这两种最小化方法的工作原理根本不同，因此对这种裁剪的反应也大不相同。这里minimize和Nelder-Mead一起使用，是一种无梯度的方法。因此，该算法不计算数值梯度，也不估计任何雅可比矩阵。相比之下，最终由 curve_fit 调用的 least-squares 正是这样做的。但是，如果函数不是连续的，则逼近梯度并由此得出任何雅可比行列式都有些问题。如前所述，这种不连续性是由 np.clip 引入的。删除后，可以很容易地看到，P0 猜测并不像包含裁剪的情况下看起来那么好。 curve_fit 确实随着 maxfev=5000 的增加而收敛，而 minimize 在将方法更改为 method='CG' 时立即失败。要查看算法的困难，可以尝试手动提供 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()
所以我会说它工作得很好。不过，我收到溢出警告。