使用 lmfit 评估模型
Evaluate model with lmfit
我想在输入变量的定义范围之外绘制计算和绘制我的模型的预测，我正在尝试使用 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()
我在这里错过了什么？
您的错误消息来自 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 数组长度相同的数组。
混合全局数组和局部数组通常不是一个好主意。这说明了这样做的问题之一。还有其他一些奇怪的事情会导致合身效果不佳。比如 a) 为什么所有初始值都相同，以及 b) 为什么要对几个参数设置非常严格（且相同）的界限？我想这些并不是真正的问题，但它们可能会使合身效果不佳。
我想在输入变量的定义范围之外绘制计算和绘制我的模型的预测，我正在尝试使用 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()
我在这里错过了什么？
您的错误消息来自 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 数组长度相同的数组。
混合全局数组和局部数组通常不是一个好主意。这说明了这样做的问题之一。还有其他一些奇怪的事情会导致合身效果不佳。比如 a) 为什么所有初始值都相同，以及 b) 为什么要对几个参数设置非常严格（且相同）的界限？我想这些并不是真正的问题，但它们可能会使合身效果不佳。