python - Scipy: how to curve_fit an implicit scalar function?

Question

Is it possible to concatenate scipy.optimize.curve_fit with scipy.optimize.bisect (or fsolve, or whatever) for implicit scalar functions?

In practice, have a look at this Python code where I try to define an implicit function and pass it to curve_fit to obtain the best fit for a parameter:

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]

Now, I am aware of all the thing that may go wrong when trying to automatically find roots (although bisection should be stable, provided there's a root between a and b); however, the error I get seems to concern the dimentionality of the output of 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

My guess is that curve_fit basically treats the output of the input function as a vector having the same dimensionality of the input data; I though I could easily work around this by 'vectorizing' the implicit function, or func_egg, although it does not seem as trivial as I thought.

Am I missing something? Is there a simple workaround?

question from:https://stackoverflow.com/questions/66051108/scipy-how-to-curve-fit-an-implicit-scalar-function

Answer 1

Waitting for answers