How to use scipy least_squares to get the estimation of unknow variables-CodePudding

I am a newbie in using scipy.optimize. I have the following function calls func. I have x and y values given as a list and need to get the estimated value of a, b and c. I could use curve_fit to get the estimation of a, b and c. However, I want to explore the possibilities of using least_squares. When I run the following code, I get the below error. It'd be great if anyone could point me to the right direction.

import numpy as np
from scipy.optimize import curve_fit
from scipy.optimize import least_squares

np.random.seed(0)

x = np.random.randint(0, 100, 100) # sample dataset for independent variables
y = np.random.randint(0, 100, 100) # sample dataset for dependent variables  

def func(x,a,b,c):
    return a*x**2   b*x   c

def result(list_x, list_y):
      popt = curve_fit(func, list_x, list_y)
    
sol = least_squares(result,x, args=(y,),method='lm',jac='2-point',max_nfev=2000)   

TypeError: ufunc 'isfinite' not supported for the input types, and the inputs could not be 
safely coerced to any supported types according to the casting rule ''safe''

CodePudding user response：

To use least_squares you need a residual function and not the curve_fit. Also, least_squares requires a guess for the parameters that you are fitting (i.e. a,b,c). In your case, if you want to use least_squares you can write something similar (I just used random values for the guess)

import numpy as np
from scipy.optimize import least_squares

np.random.seed(0)

x = np.random.randint(0, 100, 100) # sample dataset for independent variables
y = np.random.randint(0, 100, 100) # sample dataset for dependent variables  

def func(x,a,b,c):
    return a*x**2   b*x   c

def residual(p, x, y):
    return y - func(x, *p)

guess = np.random.rand(3)

sol = least_squares(residual, guess, args=(x, y,),method='lm',jac='2-point',max_nfev=2000)

CodePudding user response：

The following code uses the least_squares() routine for optimization. The most important change in comparison to your code is ensuring that func() returns a vector of residuals. I also compared the solution to the linear algebra result to ensure correctness.

import numpy as np
from scipy.optimize import curve_fit
from scipy.optimize import least_squares

np.random.seed(0)

x = np.random.randint(0, 100, 100) # sample dataset for independent variables
y = np.random.randint(0, 100, 100) # sample dataset for dependent variables

def func(theta, x, y):
    # Return residual = fit-observed
    return (theta[0]*x**2   theta[1]*x   theta[2]) - y

# Initial parameter guess
theta0 = np.array([0.5, -0.1, 0.3])

# Compute solution providing initial guess theta0, x input, and y input
sol = least_squares(func, theta0, args=(x,y))
print(sol.x)

#------------------- OPTIONAL -------------------#
# Compare to linear algebra solution
temp = x.reshape((100,1))
X = np.hstack( (temp**2, temp, np.ones((100,1))) )
OLS = np.linalg.lstsq(X, y.reshape((100,1)), rcond=None)
print(OLS[0])