Scipy：如何曲线拟合隐式标量函数？

Question

对于隐式标量函数，是否可以将 scipy.optimize.curve_fit 与 scipy.optimize.bisect（或 fsolve 或其他）连接起来？

实际上，看看这段 Python 代码，我在其中尝试定义一个隐式函数并将其传递给curve_fit 以获得参数的最佳拟合：

import numpy as np
import scipy.optimize as opt
import scipy.special as spc

# Estimate of initial parameter (not really important for this example)
fact, _, _, _ = spc.airy(-1.0188)
par0 = -np.log(2.0*fact*(18**(1.0/3.0))*np.pi*1e-6)

# Definition of an implicit parametric function f(c,t;b)=0
def func_impl(c, t, p) :
    return ( c - ((t**3)/9.0) / ( np.log(t*(c**(1.0/3.0))) + p ) )

# definition of the function I believe should be passed to curve_fit
def func_egg(t, p) :
    x_st, _ = opt.bisect( lambda x : func_impl(x, t, p), a=0.01, b=0.3 )
    return x_st

# Some data points
t_data = np.deg2rad(np.array([95.0, 69.1, 38.8, 14.7]))
c_data = np.array([0.25, 0.10, 0.05, 0.01])

# Call to curve_fit
popt, pcov = opt.curve_fit(func_egg, t_data, c_data, p0=par0)
b = popt[0]

现在，我知道在尝试自动查找根时可能会出错的所有事情（尽管二分法应该是稳定的，前提是 a 和 b);但是，我得到的错误似乎与func_impl 的输出维度有关：

Traceback (most recent call last):
  File "example_fit.py", line 23, in <module>
    popt, pcov = opt.curve_fit(func_egg, t_data, c_data, p0=par0)
  File "/usr/local/lib/python3.7/site-packages/scipy/optimize/minpack.py", line 752, in curve_fit
    res = leastsq(func, p0, Dfun=jac, full_output=1, **kwargs)
  File "/usr/local/lib/python3.7/site-packages/scipy/optimize/minpack.py", line 383, in leastsq
    shape, dtype = _check_func('leastsq', 'func', func, x0, args, n)
  File "/usr/local/lib/python3.7/site-packages/scipy/optimize/minpack.py", line 26, in _check_func
    res = atleast_1d(thefunc(*((x0[:numinputs],) + args)))
  File "/usr/local/lib/python3.7/site-packages/scipy/optimize/minpack.py", line 458, in func_wrapped
    return func(xdata, *params) - ydata
  File "example_fit.py", line 15, in func_egg
    x_st, _ = opt.bisect( lambda x : func_impl(x, t, p), a=0.01, b=0.3 )
  File "/usr/local/lib/python3.7/site-packages/scipy/optimize/zeros.py", line 550, in bisect
    r = _zeros._bisect(f, a, b, xtol, rtol, maxiter, args, full_output, disp)
  File "example_fit.py", line 15, in <lambda>
    x_st, _ = opt.bisect( lambda x : func_impl(x, t, p), a=0.01, b=0.3 )
  File "example_fit.py", line 11, in func_impl
    return ( c - ((t**3)/9.0) / ( np.log(t*(c**(1.0/3.0))) + p ) )
TypeError: only size-1 arrays can be converted to Python scalars

我的猜测是curve_fit基本上将输入函数的输出视为与输入数据具有相同维度的向量；虽然我可以通过“矢量化”隐式函数或func_egg 轻松解决这个问题，尽管它看起来不像我想象的那么微不足道。

我错过了什么吗？有没有简单的解决方法？

Answer 1

我想我最终会回答我自己的问题...我希望这对其他人有用。

我们先选择一个更简单的隐函数，在本例中为 f(c,t;b)=c-b*t^3（原因后面会说清楚）：

import numpy as np
import scipy.optimize as opt
import scipy.special as spc
import matplotlib.pyplot as plt

# Definition of an implicit parametric function f(c,t;b)=0
def func_impl(c, t, p) :
    return (c-p*t**3)

让我们矢量化它：

v_func_impl = np.vectorize(func_impl)

与问题中的脚本相同，但现在 (1) func_egg 已矢量化并且 (2) 我使用 newton 而不是 bisect（我发现提供 x0 而不是 @987654329 更容易@):

# Definition of the function I believe should be passed to curve_fit
def func_egg(t, p) :
    x_st = opt.newton( lambda x : func_impl(x, t, p), x0=0.05 )
    return x_st

v_func_egg = np.vectorize(func_egg)

# Some data points
t_data = np.deg2rad(np.array([127.0, 95.0, 69.1, 38.8]))
c_data = np.array([0.6, 0.25, 0.10, 0.05])

# Call to curve_fit
popt, pcov = opt.curve_fit(v_func_egg, t_data, c_data, p0=par0)
b = popt[0]

现在可以了！

plt.plot(t_data, c_data)
plt.plot(np.linspace(0.5, 2.5), b*np.linspace(0.5, 2.5)**3)

所以，本质上：

• 为了连接 scipy 曲线拟合和求根，需要确保每个函数都是矢量化的（或者可以处理 numpy 数组作为输入和输出）
• 确保您的函数不是“太难看”，否则即使串联工作，寻根过程本身也可能无法找到结果（这涉及到数值数学；我应该检查一下我原来的规律性函数）