So I don't know how to do this, and I'm not sure if there's some math that I'm completely forgetting right now but I'm at a loss. In short, I have some fairly simple initial code using scipy and numpy, and I want to fit an exponential curve to it:
from scipy.optimize import curve_fit
import numpy as np
# sample data
x = np.array([7620., 7730., 7901., 8139., 8370., 8448., 8737., 8824., 9089., 9233., 9321., 9509., 9568., 9642., 9756., 9915., 10601., 10942.])
y = np.array([0.01228478, 0.01280466, 0.01363498, 0.01493918, 0.01530108, 0.01569484, 0.01628133, 0.01547824, 0.0171548, 0.01743745, 0.01776848, 0.01773898, 0.01839569, 0.01823377, 0.01843686, 0.01875542, 0.01881426, 0.01977975])
# define type of function to search
def model_func(x, a, k, b):
return a * np.exp(-k*x) + b
# curve fit
p0 = (2.e-6,300,1)
opt, pcov = curve_fit(model_func, x, y, p0)
a, k, b = opt
# test result
x2 = np.linspace(7228, 11000, 3000)
y2 = model_func(x2, a, k, b)
fig, ax = plt.subplots()
ax.plot(x2, y2, color='r', label='Fit. func: $f(x) = %.3f e^{%.3f x} %+.3f$' % (a,k,b))
ax.plot(x, y, 'bo', label='data with noise')
ax.legend(loc='best')
plt.show()
My issue is try as I may I can't figure out the initial parameters for p0- I've tried a range of values, but frankly I have no idea what I'm doing so I'm not getting a solution here. Can someone suggest how to do it? Thank you!
To me this data is more close to a lign than an exponential curve, are you sure your model is right?
About the initial guesses, I assume you don't have any further knowledge about the function, if so use it:
For x->\inf the fucntion approaches b. So I would use a guess of about 0.025 for b.
For the other two variables, take a subsample of x and y and solve the equation explicit:
a * np.exp(-k*x[0]) + 0.025 = y[0]
a * np.exp(-k*x[-1]) + 0.025 = y[-1]
Solving this gives:
a = (y[0]-0.025)/np.exp(-k*x[0])
e^(-k*(x[-1]-x[0])=(y[-1]-0.025)/(y[0]-0.025) # and then take logarithm
k = -np.log((y[-1]-0.025)/(y[0]-0.025))/(x[-1]-x[0])
a = (y[0]-0.025)/np.exp(-k*x[0])
gives k=0.00026798747972760543 and a=-0.80114087848462689
This method can be used generaly, throw away as many points as necessarry to solve the equation exact and then use these values for starting optimal values