python - 如何使用 python 的 lmfit 检索常微分方程 (ODE) 的最佳约束系数率？
问题描述
我正在尝试使用 python 最小化常微分方程（ODE）问题的损失函数。
此损失函数旨在检索用户定义的 ODE 的最佳系数率。
尽管我在互联网上找到了我的努力和示例，但所有结果都以 ODE 结尾，其解决方案返回的系数值超出了各自的界限。
我知道这可能有几个原因：最小化器的选择、错误的 ODE 定义、集成步长等等。
尽管如此，我的主要问题是找到一个示例，其中最小化器可以支持 ODE 集成中的集成“中断”。
例如，如果我对简单的 ODE 模型感兴趣，例如 Prey-Predator ODE (Lotka-Volterra)，则必须设置的一条规则是没有种群可以达到负值。
因此，如果 ODE 最小化器由于某种给定的原因为所建模的两个种群（Prey 和 Predator）中的任何一个返回负值，则最小化器必须停止其积分。
然而，看起来，这种规则不受常见最小化的支持，至少对于我测试的 lmfit 的 Minimizer 对象而言。
即使为每个 Parameter 对象设置了边界规则，Minimizer 仍会继续插入超出给定阈值的数据（即：负总体）。
例如：
假设我想根据一些经验数据找到 Prey-Predator ODE 的最佳系数率。
在这种情况下，如果我使用 python 的 lmfit 库来解决我的问题，我可以执行以下操作：
from lmfit import Parameters, report_fit, Minimizer
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import ode
end = '\n'*2 + '-'*50 + '\n'*2
def ode_f(t, xs, ps):
"""Lotka-Volterra predator-prey model."""
if isinstance(ps, Parameters):
Prey_Natality = ps[r'Prey_Natality'].value
Predation_rate = ps['Predation_rate'].value
Predator_Natality = ps['Predator_Natality'].value
Predator_Mortality = ps['Predator_Mortality'].value
Prey_natural_mortality = ps['Prey_natural_mortality'].value
else:
(Prey_Natality, Predation_rate, Predator_Natality,
Predator_Mortality, Prey_natural_mortality) = ps
Prey, Predator = xs
dPrey_dt = ( + Prey_Natality*Prey
- Predation_rate*(Prey)*Predator
- Prey_natural_mortality*Prey)
dPred_dt = ( + Predation_rate* (Prey ) *Predator
+ Predator_Natality*Predator
- Predator_Mortality*Predator)
return [dPrey_dt, dPred_dt]
def solout(y):
if np.any( y <=0):
return -1 # stop integration
else:
return 0
def ode_solver(t,
x0,
ps,
mode='vode',
nsteps=500,
method='bdf'):
"""
Solution to the ODE x'(t) = f(t,x,k) with initial condition x(0) = x0
"""
r = ode(ode_f).set_integrator(mode,
nsteps=nsteps,
method=method)
t0 = t.min()
tmax = t.max()
dt = np.diff(t)[0]
r.set_initial_value(x0, t0).set_f_params(ps)
y = []
times = []
integration_time = r.t
while r.successful() and integration_time < tmax:
r.integrate(integration_time + dt)
integration_time = r.t
yi = r.y
y.append(yi)
times.append(integration_time)
if solout(yi) == -1:
print('Time stoped at: {0}'.format(integration_time))
break
return (np.array(y)).astype(float), times
def ODE_solver_residual_evaluator(ps, ts, data):
x0 = [ps['Prey_Pop'].value, ps['Predator_Pop'].value]
model, times = ode_solver(ts, x0, ps)
# if data.shape[0] <= model.shape[0]:
# data.resize((model.shape[0],data.shape[1]),refcheck=False)
# else:
# model.resize((data.shape[0],data.shape[1]),refcheck=False)
return ( model[:len(data )] - data[:len(model)] )
def residual_dim_reducer(residual_array):
return np.square(residual_array).sum()
if '__main__' == __name__:
########
dt = 10**(-4)
t = np.arange(0, 100, dt)
Prey_initial_pop = 1200
Predator_initial_pop = 50
x0 = np.array([Prey_initial_pop,
Predator_initial_pop])
Prey_Natality = 2.6
Predation_rate = 0.12
Predator_Natality = 0.401
Predator_Mortality = 0.0025
Prey_natural_mortality = 0.001
true_params = np.array((Prey_Natality,
Predation_rate,
Predator_Natality,
Predator_Mortality,
Prey_natural_mortality))
data, times = ode_solver(t, x0, true_params)
data += np.random.lognormal(size=data.shape)*0.5
Param_Names = ['Prey_Natality',
'Predation',
'Predator_Natality',
'Predator_Mortality',
'Prey_natural_mortality']
Populations = ['Prey population',
'Predator population']
for i in range(data.shape[1]):
plt.plot(times, np.real(data[:,i]), 'o', label=r'original {0}'.format(Populations[i]))
plt.legend()
plt.show()
import os
to_save = os.path.join(os.getcwd(), 'original data.png')
plt.savefig(to_save)
print('Creating the minizer object', end=end)
################3
# set parameters incluing bounds
params = Parameters()
params.add('Prey_Pop', value= 100, min=1, max=600)
params.add('Predator_Pop', value=10, min=1, max=400)
params.add('Prey_Natality', value=0.03, min=0.000001, max=3.5)
params.add('Predation_rate', value=0.02, min=0.00003, max=3.5)
params.add('Predator_Natality', value=0.0004, min=0.00001, max=3.2)
params.add('Predator_Mortality', value=0.003, min=0.000001, max=3.2)
params.add('Prey_natural_mortality', value=0.001, min=0.0000001, max=0.2)
fitter = Minimizer(ODE_solver_residual_evaluator,
params,
fcn_args=(t, data),
iter_cb=None,
scale_covar=True,
nan_policy='propagate',
reduce_fcn=residual_dim_reducer,
calc_covar=True,
)
fitted_ODE = fitter.minimize(method='Nelder-Mead')
Optimum_params = fitted_ODE.params.valuesdict()
x0_optimum = np.array([Optimum_params.pop('Prey_Pop'),
Optimum_params.pop('Predator_Pop')])
Y_fitted, times_fitted = ode_solver(t, x0_optimum, Optimum_params.values())
Param_Names = list(params.keys())
print(end, 'Param_Names: {0}'.format(Param_Names))
##################
data = data[:len(Y_fitted)]
Y_fitted = Y_fitted[:len(data)]
from sklearn.metrics import (explained_variance_score,
r2_score,
mean_squared_error)
explained_variance_score = explained_variance_score(Y_fitted, data)
R2 = r2_score(Y_fitted, data)
RMSE = np.sqrt(mean_squared_error(Y_fitted, data))
print('Explained variance: {0} \n R2: {1} \n RMSE: {2}'.format(explained_variance_score,
R2, RMSE))
print(end, 'Fitting by Minizer is complete', end=end)
# display fitted statistics
report_fit(fitted_ODE)
print(end)
######################3
fig2 = plt.figure()
n_ref_markers = 12
markers_on = np.linspace(0, data[:,i].size-1, n_ref_markers).astype(int).tolist()
for i in range(Y_fitted.shape[1]):
plt.plot(times_fitted[:len(times)],
Y_fitted[:len(times),i],
'-',
linewidth=1.1,
label=r'fitted {0} '.format(Param_Names[i]))
plt.plot(np.arange(len(data)), data[:,i], 'o',
markevery=markers_on,
label=r'original {0}'.format(Param_Names[i]))
fig2.legend()
fig2.show()
返回值为：
- 两个数字：A) 原始数据 B) 原始数据加上 ODE 的建模数据，根据 lmfit 返回的最佳系数比率
- 包含拟合测量和参数统计描述的回报表。
返回输出：
--------------------------------------------------
Fitting by Minizer is complete
--------------------------------------------------
Models fitting error:
Explained variance: -90.1682809468072
R2: -3125.4358694840084
RMSE: 785.9581933129715
--------------------------------------------------
[[Fit Statistics]]
# fitting method = Nelder-Mead
# function evals = 66437
# data points = 364
# variables = 7
chi-square = 2.2485e+08
reduced chi-square = 629842.649
Akaike info crit = 4867.50583
Bayesian info crit = 4894.78590
## Warning: uncertainties could not be estimated:
[[Variables]]
Prey_Pop: 117.479436 +/- 16.3998313 (13.96%) (init = 100)
Predator_Pop: 397.552948 +/- nan (nan%) (init = 10)
Prey_Natality: 1.05567443 +/- nan (nan%) (init = 0.03)
Predation_rate: 3.46543190 +/- 0.05124925 (1.48%) (init = 0.02)
Predator_Natality: 0.48528830 +/- nan (nan%) (init = 0.0004)
Predator_Mortality: 2.76733581 +/- 0.06777831 (2.45%) (init = 0.003)
Prey_natural_mortality: 0.03928761 +/- 0.00503378 (12.81%) (init = 0.001)
[[Correlations]] (unreported correlations are < 0.100)
C(Predation_rate, Prey_natural_mortality) = -0.596
C(Prey_Pop, Predator_Mortality) = 0.179
C(Predation_rate, Predator_Mortality) = 0.141
C(Prey_Pop, Prey_natural_mortality) = 0.127
C(Predator_Mortality, Prey_natural_mortality) = -0.112
这个问题来自从这里提取的一个教学示例：
任何有关该主题的帮助将不胜感激。
真挚地，
解决方案
推荐阅读
- javascript - 无法通过 Rosbridge 订阅里程计 ROS 消息
- php - 按类别分类的 WordPress 相关内容为自定义帖子类型中的每个帖子产生相同的结果
- visual-studio-code - 如何在离线网络中安装 VSCode 帮助文件
- azure-devops - 在 Azure Devops 上附加到 PATH？
- javascript - Push javascript array objects from mysql query
- gnuplot - 如何根据星期几在轴上标记抽动
- r - R：在 XY 网格中查找缺失的组合
- angular - ExpressionChangedAfterItHasBeenCheckedError Angular *ngIf 指令导致错误
- android - GStreamer 管道正在播放，但没有填充回调（仅在 Android 上发生，在 Windows 和 Linux 上有效）
- sql - 无法在 Oracle 查询中按日期分组