如何使用SciPy曲线拟合隐式标量函数？

对于隐式标量函数，是否可以将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进行“矢量化”来轻松地解决这个问题，尽管它看起来并不像我想象的那样微不足道。
我错过了什么吗？
是否有简单的解决方法？

我想我最终回答了我自己的问题。我希望这对其他人有用。
让我们首先选择一个更简单的隐函数，在本例中，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比提供[a,b]更容易）：


# 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

par0 = 0.05
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)
plt.show()

所以，从本质上说：

为了连接scipy曲线拟合和求根，需要确保每个函数都是矢量化的（或者可以处理numpy数组作为输入和输出）。
确保你的函数不是“太丑”，否则即使连接工作的根查找过程本身可能无法找到一个结果（这进入数值数学;我应该检查我的原始函数的规律性）。

如何使用SciPy曲线拟合隐式标量函数？

1条答案

相关问题

热门标签

最新问答