Оценить модель с помощью lmfit

Question

Я хочу построить расчет и построить прогноз моей модели из заданного диапазона входных переменных, и я пытаюсь сделать это с lmfit. Вот код, над которым я работаю. Модель определяется как функция только от одной независимой переменной t, но сама функция использует также другой независимый набор наблюдений. Я могу рассчитать new_times для оценки модели, но не могу сделать то же самое для другого набора наблюдений. Помимо ужасной подгонки (здесь не проблема сама по себе), я выделил ошибки, которые я получаю, если использую lmfit, как мне показалось, что это работает:

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()

Что мне здесь не хватает?

Answer 1

Ваше сообщение об ошибке приходит от scipy.integrate.cumtrapz(). Чтение полной трассировки покажет вам это. В сообщении говорится, что два аргумента scipy.integrate.cumtrapz() должны быть одинакового размера.

Действительно, при it.cumtrapz(-a * p**c + b, t-T0, initial=0) первый аргумент будет иметь размер, установленный вашей переменной p, которая является фиксированной Длина массива и второй будет установлена ​​t, независимой переменной вашей модельной функции.

Итак, когда вы делаете result.eval(t=NEW_ARRAY), вы обязательно получите это сообщение об ошибке, если ваш новый массив не будет иметь такую ​​же длину, как ваш массив p. Это настоящая ошибка. Как бы вы решили это? Ну, вам нужно было бы как-то передать массив для p, который был бы такой же длины, как новый массив t.

Как правило, смешивать глобальные и локальные массивы - плохая идея. Это иллюстрирует одну из проблем в этом. Происходит несколько других странных вещей, которые делают подгонку не очень хорошей. Например: а) почему все начальные значения одинаковы, и б) почему вы устанавливаете чрезвычайно жесткие (и одинаковые) границы для нескольких параметров? Я думаю, что это не совсем вопрос, но они, вероятно, приведут к тому, что подгонка не будет работать очень хорошо.