So I am trying to fit some electrochemical data into a Michaelis Menten model (ODE) with python and the library lmfit. I put the code of the fitting below, which comes from an example at Text:
from lmfit import minimize, Parameters, Parameter, report_fit
from scipy.integrate import odeint
df = pd.read_csv("example3.csv", header=None, nrows=300)
data=df[2].values
t = time_x
def f(xs, t, ps):
R=1
try:
alpha = ps['alpha'].value
Vmax = ps['Vmax'].value
Km=ps['Km'].value
except:
alpha, Vmax, Km = ps
S = xs
dsdt=R-alpha*(Vmax*S/(Km+S))
return dsdt
def g(t, x0, ps):
solution = odeint(f, x0, t, args=(ps,))
return solution
def residual(ps, ts, data):
x0 = ps['x0'].value
model = g(ts, x0, ps)
return (model - data).ravel()
# set parameters incluing bounds
parameters = Parameters()
parameters.add('x0', value=float(data[0]), min=0, max=100)
parameters.add('Vmax',value=18,min=0, max=100)
parameters.add('alpha',value=1,min=0, max=1)
parameters.add('Km',value=5,min=0,max=100)
# fit model and find predicted values
result = minimize(residual, parameters, args=(t, data), method='leastsq')
final = data + result.residual.reshape(data.shape)
# plot data and fitted curves
plt.plot(t, data, 'o')
plt.plot(t, final, '--', linewidth=2, c='blue');
# display fitted statistics
report_fit(result)
And I get the following error:
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
<ipython-input-199-1bd65276a496> in <module>
38 # fit model and find predicted values
39 result = minimize(residual, parameters, args=(t, data), method='leastsq')
---> 40 final = data + result.residual.reshape(data.shape)
41
42 # plot data and fitted curves
ValueError: cannot reshape array of size 90000 into shape (300,)
I know what it means, residual array is of larger size that data and it cannot reshape it. But technically residuals should be the same size as data. If someone is familiar with the lmfit library it would be a massive help. I hope it is a dumb question, I just cannot see the error. Thank you in advance.
Most likely you are having a problem with the broadcast of arrays, the result of odeint is a column vector, the data treated as row, the difference is broadcast as matrix having 300*300=90000 elements.
Reduce the shape of the odeint result at some place, for instance in
return solution[:,0]