I am trying to use leastsq from the scipy.optimize module to find a best fit line, where there are 3 unknown parameters. I have written out the code however the program runs and returns the initial guess as the optimization parameters (essentially the leastsq function does nothing in my program).
Here is the simplified code, with the data I am using.
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize #Leastsq Levenberg-Marquadt Algorithm
a = [6.011737374832778931e+10, 1.253141174941418152e+11, 1.297179270983954620e+11, 1.577611269699349976e+11, 2.238721138568337708e+11, 4.315190768277650146e+11, 5.407543229455815430e+11, 5.382697825162881470e+11, 5.308844442309879150e+11, 4.528975799036213379e+11, 2.890679882365477905e+11, 2.798981319634357300e+11, 2.798981319634357300e+11]
b = [1.228900000000000006e+02, 1.465500000000000114e+02, 1.761399999999999864e+02, 2.057199999999999989e+02, 2.353100000000000023e+02,2.648999999999999773e+02, 2.945000000000000000e+02, 3.315000000000000000e+02, 3.758999999999999773e+02, 4.203199999999999932e+02, 4.647400000000000091e+02, 5.091700000000000159e+02, 5.980399999999999636e+02]
a_arr = np.asarray(a) #electron density
b_arr = np.asarray(b) #altitude
#function where p is the list of variables that go into the fitting function and x is the independent variable
def fitfunc (p,x):
return p[0] * np.exp((0.5 * (1- (x - p[1])/p[2])) - np.exp(-(x - p[1])/p[2]))
def errfunc (p,x,y):
err = y - fitfunc(p,x)
return err
#fitfunc = lambda p, x: p[0] * np.exp((0.5 * (1- (x - p[1])/p[2])) - np.exp(-(x - p[1])/p[2]))
#errfunc = lambda p, x, y: fitfunc(p,x)-y
p0 = [10**11.7,300.,50.]
#p = np.asarray(p0)
p1, success = optimize.leastsq(errfunc, p0[:], args=(a_arr, b_arr))
b_arr2 = np.linspace(b_arr.min(), b_arr.max(), 80)
print("estimated parameters: ", p1)
print("observed parameters: ", p0)
#with OPTIMIZATION
plt.plot(a_arr, b_arr,"o--", fitfunc(p1, b_arr2 ), b_arr2, '-')
plt.legend(['data', 'optimization'], loc='best')
plt.title('plot')
plt.xlabel('x')
plt.ylabel('y')
plt.show()
Everytime I get ('estimated parameters: ', array([ 5.01187234e+11, 3.00000000e+02, 5.00000000e+01])) ('observed parameters: ', [501187233627.2715, 300.0, 50.0]) which is the same thing.
I have tried use lambda functions, changing around the data types.
Is it that my errfunc (cost/objective function) is incorrect?
Thanks in advance!
If I understand the code correctly, I think you're passing the arguments in the wrong order to errfunc.
Instead of:
p1, success = optimize.leastsq(errfunc, p0[:], args=(a_arr, b_arr))
it should be
p1, success = optimize.leastsq(errfunc, p0[:], args=(b_arr, a_arr))
Your errfunc wants the b array first and a array second.
When I make the above change, I get the observed parameters as:
[7.43278755e+11, 2.77592278e+02, 8.88750029e+01]
And the plot looks much better as well.