Evaluate model with lmfit

Question

I want to plot calculate and plot the predictions of my model out of the defined range of the input variables, and I am trying to do this with lmfit. Here is the code that I am working on. The model is defined as a function of just one independent variable t, but the function itself uses also another independent set of observations. I can calculate the new_times to which evaluate the model, but can not do the same for the other set of observations. Besides the terrible fit (not an issue per se here), I have highlighted the errors I get in case I use lmfit as I thought it worked:

import numpy as np
import matplotlib.pyplot as plt

from lmfit import Model
import scipy.integrate as it
import scipy.constants as scc

times=np.asarray([58418, 58422, 58424, 58426, 58428, 58430, 58540, 58542, 58650, 58654, 58656, 58658, 58660, 58662, 58664, 58666, 58668, 58670, 58672, 58674, 58768, 58770, 58772, 58774, 58776, 58778, 58780, 58782, 58784, 58786, 58788, 58790, 58792, 58794, 58883, 58884, 58886, 58888, 58890, 58890, 58892, 58894, 58896, 58898, 58904])

y_obs =np.asarray([ 0.00393986,0.00522288,0.00820794,0.01102782,0.00411525,0.00297762, 0.00463183,0.00602662,0.0114886, 0.00176694,0.01241464,0.01316199, 0.01108201, 0.01056611, 0.0107585, 0.00723887,0.0082614, 0.01239229, 0.00148118,0.00407329,0.00626722,0.01026926,0.01408419,0.02638901, 0.02284189, 0.02142943, 0.02274845, 0.01315814, 0.01155898, 0.00985705, 0.00476936,0.00130343,0.00350376,0.00463576, 0.00610933, 0.00286234, 0.00845177,0.00849791,0.0151215, 0.0151215, 0.00967625,0.00802465, 0.00291534, 0.00819779,0.00366089])

y_obs_err = np.asarray([6.12189334e-05, 6.07487598e-05, 4.66365096e-05, 4.48781264e-05, 5.55250430e-05, 6.18699105e-05, 6.35339947e-05, 6.21108524e-05, 5.55636135e-05, 7.66087180e-05, 4.34256323e-05, 3.61131000e-05, 3.30783270e-05, 2.41312040e-05, 2.85080015e-05, 2.96644612e-05, 4.58662869e-05, 5.19419065e-05, 6.00479888e-05, 6.62586953e-05, 3.64830945e-05, 2.58120956e-05, 1.83249104e-05, 1.59433858e-05, 1.33375408e-05, 1.29714326e-05, 1.26025166e-05, 1.47293107e-05, 2.17933175e-05, 2.21611713e-05, 2.42946630e-05, 3.61296843e-05, 4.23009806e-05, 7.23405476e-05, 5.59390368e-05, 4.68144974e-05, 3.44773949e-05, 2.32907036e-05, 2.23491451e-05, 2.23491451e-05, 2.92956472e-05, 3.28665479e-05, 4.41214301e-05, 4.88142073e-05, 7.19116984e-05])

p= np.asarray([ 2.82890497,3.75014266,5.89347542,7.91821558,2.95484056,2.13799544, 3.32575733,4.32724456,8.2490644, 1.26870083,8.91397925,9.45059128, 7.95712563, 7.58669608, 7.72483557,5.19766853,5.93186433,8.89793105, 1.06351782,2.92471065,4.49999613,7.37354766, 10.11275281, 18.94787684, 16.40097363, 15.38679306, 16.33387783, 9.44782842, 8.29959664,7.07757293, 3.42450524,0.93588962,2.515773,3.32857547,7.180216, 2.05522399, 6.06855409,6.1016838,10.8575614,10.8575614, 6.94775991,5.76187014, 2.09327787, 5.88619335,2.62859611])

def new_f_function(t, sum, f0, a, b, c, T0):

  obs_f = f0 + it.cumtrapz(-a * p**c + b, t-T0, initial=0)
  new_f = obs_f*(1+sum/scc.c)

  return new_f

# Create Model
model = Model(new_f_function, independent_vars=['t'])

# Initialize Parameter
params = model.make_params()

params['sum'].value = 1.483 
params['sum'].min = 1.47
params['sum'].max = 1.50

params['f0'].value = 1.483 
params['f0'].min = 1.47
params['f0'].max = 1.50

params['a'].value = 1.483 
params['a'].min = 1.47
params['a'].max = 1.50

params['b'].value = 1.483 

params['c'].value = 1.483 

params['T0'].value = 1.483 

result = model.fit(y_obs, params, weights=(1./y_obs_err), t=times, scale_covar=False)
print result.fit_report()

# New x-values to evaluate the model 
x_fit = np.linspace(min(times)-10., max(times)+10, 1e4)  

fig, (ax1, ax2) = plt.subplots(nrows=2, ncols=1, figsize=(3, 6), sharex='all')
ax1.scatter(x=times, y=p, marker='+', c='k')
ax2.errorbar(x=times, y=y_obs, yerr=y_obs_err, marker='.', ls='', label='DATA')
ax2.plot(times, result.best_fit, label='best fit')

new_predictions = result.eval(**result.best_values)
#ax2.plot(x_fit, new_predictions, label='extended fit') # This gives error: `ValueError: x and y must have same first dimension, but have shapes (10000,) and (45,)`

predicted = result.eval(t=x_fit) # This gives error: `raise ValueError("If given, length of x along axis must be the "ValueError: If given, length of x along axis must be the same as y.`

plt.legend()
plt.show()

What am I missing here?

Answer 1

Your error message is coming from scipy.integrate.cumtrapz(). Reading the full traceback would show you that. The message is saying that the two arguments of scipy.integrate.cumtrapz() must be the same size.

Indeed, at it.cumtrapz(-a * p**c + b, t-T0, initial=0) the first argument will have a size set by your variable p, which is a fixed-length array and the second will be set by t, the independent variable of your model function.

So, when you do result.eval(t=NEW_ARRAY) you would definitely get that error message if your new array was not the same length as your array p. That is the actual error. How would you solve it? Well, you would have to somehow pass in an array for p that was the same length as the new t array.

It's generally a bad idea to mix global arrays and local arrays. This illustrates one of the problems in doing so. There are several other weird things going on that will make the fit not work very well. Like a) why are all initial values the same, and b) why are you placing extremely tight (and the same) bounds on several of the parameters? I guess those are not really the question here, but they are probably going to make the fit not work very well.