关键词：
【中文标题】Python问题：ODE系统的参数估计，使得系数是函数【英文标题】：Python question: parameter estimation of ODE systems such that coefficients are functions 【发布时间】：2021-08-06 13:09:39 【问题描述】：我的模型是 ODE 系统，系数是具有参数的函数。
我的目标是系数中的参数估计。
dC1dt = c_11(t)*C1 + c_12(t)*C2
dC2dt = c_22(t)*C2 + c_22(t)*C2
我知道如何用 R 解决它，但我必须用 Python 解决它。
首先，我尝试在固定条件下使用 odeint。
完整代码：
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
# fixed values
V_room = 40 #room volume[m^3]
v = 0.01 #velocity of material [m/s]
h = 2.5 # height of room [m/s]
Q = 1 # average flux [m^3/min]
V_0 =1 # initial value of nf [m^3]
Q_0 = 0.1 #initial value of Q_ff
#initial values
C_0 = [0.0002, 2e-05]
#parameter
k= 0.06
beta = 0.001
alpha = 0.005
# measured time
ts = np.linspace(0,60,60)
# fixed functions
def V_nf(t):
return V_0 + (V_room - V_0)*(1-np.exp(-alpha*t))
def V_ff(t):
return V_room - V_nf(t)
def dVnf(t):
return alpha*(V_nf(t)-V_0)
def dVff(t):
return -dVnf(t)
def S_nf(t):
return V_nf(t)**(2/3)
def Q_ff(t):
return Q_0*np.exp(-beta*t)
=====
coefficients = 'coef_11' : coef11, 'coef_12' : coef12, 'coef_21' : coef21, 'coef_22' : coef22
def coef11(t):
return (-k*S_nf(t)-dVnf(t))/V_nf(t) - Q/V_room - v/h
def coef12(t):
return (k*S_nf(t)+Q_ff(t))/V_nf(t) - Q*V_ff(t)/(V_room * V_nf(t))
def coef21(t):
return k*S_nf(t)/V_ff(t)
def coef22(t):
return (-dVff(t)-Q_ff(t)-k*S_nf(t))/V_ff(t) - v/h
def dCdt(C,t):
C_nf = C[0]
C_ff = C[1]
c_11 = coefficients['coef_11'](t)
c_12 = coefficients['coef_12'](t)
c_21 = coefficients['coef_21'](t)
c_22 = coefficients['coef_22'](t)
dCnf = c_11* C_nf + c_12 * C_ff
dCff = c_21* C_nf + c_22 * C_ff
return [dCnf, dCff]
C = odeint(dCdt, C_0, ts, args= (coefficients,))
但是，我得到了
TypeError Traceback (most recent call last)
<ipython-input-114-7cece3e6a166> in <module>
13 return [dCnf, dCff]
14
---> 15 C = odeint(dCdt, C_0, ts, args= (coefficients,))
~\Anaconda3\lib\site-packages\scipy\integrate\odepack.py in odeint(func, y0, t, args, Dfun, col_deriv, full_output, ml, mu, rtol, atol, tcrit, h0, hmax, hmin, ixpr, mxstep, mxhnil, mxordn, mxords, printmessg, tfirst)
242 full_output, rtol, atol, tcrit, h0, hmax, hmin,
243 ixpr, mxstep, mxhnil, mxordn, mxords,
--> 244 int(bool(tfirst)))
245 if output[-1] < 0:
246 warning_msg = _msgs[output[-1]] + " Run with full_output = 1 to get quantitative information."
TypeError: dCdt() takes 2 positional arguments but 3 were given
不知道这个方法是否正确，为什么会出现这个问题。
正如我所说，我的目标是估计参数（k，alpha，beta）。
因此，如果您知道如何做到这一点，请使用数据。
数据：
C1 = array([2.00000e-04, 3.52006e+00, 1.01378e+00, 8.47760e-01, 6.19000e-01,
6.09940e-01, 4.35010e-01, 3.49150e-01, 3.87830e-01, 3.24830e-01,
1.97040e-01, 1.84630e-01, 1.72520e-01, 1.35980e-01, 1.38430e-01,
1.21520e-01, 1.46680e-01, 9.08900e-02, 1.09650e-01, 8.71000e-02,
9.24400e-02, 1.11200e-01, 7.96600e-02, 9.26300e-02, 8.08700e-02,
8.04000e-02, 7.69200e-02, 9.06400e-02, 7.30600e-02, 7.22200e-02,
6.57900e-02, 7.87200e-02, 7.67700e-02, 7.28100e-02, 6.65600e-02,
6.87300e-02, 7.65600e-02, 6.87200e-02, 7.52700e-02, 8.66000e-02,
6.90700e-02, 6.51900e-02, 5.39200e-02, 6.46600e-02, 5.96400e-02,
6.90100e-02, 6.23700e-02, 1.00230e-01, 9.05900e-02, 4.96300e-02,
6.75400e-02, 5.50200e-02, 4.75700e-02, 5.82800e-02, 5.34400e-02,
5.30200e-02, 4.10700e-02, 4.10800e-02, 4.91300e-02, 4.16300e-02])
C2 = array([5.000e-05, 2.000e-05, 1.000e-04, 1.250e-03, 1.500e-03, 2.630e-03,
2.600e-04, 5.540e-03, 2.160e-03, 9.420e-03, 8.030e-03, 1.369e-02,
9.620e-03, 1.527e-02, 1.209e-02, 9.600e-03, 1.081e-02, 1.414e-02,
1.522e-02, 1.244e-02, 2.223e-02, 2.312e-02, 2.683e-02, 2.398e-02,
2.658e-02, 2.841e-02, 2.123e-02, 3.052e-02, 3.349e-02, 4.027e-02,
3.110e-02, 3.134e-02, 3.185e-02, 2.964e-02, 4.078e-02, 2.530e-02,
4.443e-02, 2.856e-02, 2.459e-02, 2.568e-02, 2.104e-02, 2.477e-02,
2.297e-02, 2.512e-02, 2.133e-02, 2.002e-02, 1.427e-02, 3.033e-02,
1.946e-02, 2.173e-02, 1.800e-02, 1.346e-02, 2.039e-02, 2.132e-02,
1.416e-02, 1.376e-02, 1.079e-02, 6.640e-03, 1.324e-02, 1.056e-02])
【问题讨论】：
您的意思是linspace(0,60,61) 而不是linspace(0,60,60)？使用 (0,60,61)，您得到的时间步长正好为 1.0，但使用 (0,60,60)，您得到的时间步长为 1.0169。另外，你确定你的 C1 数据是好的，还是有错字？在我看来，第一个数据点不属于。
【参考方案1】：
您收到 TypeError 是因为您忘记了 dCdt 函数的定义中的参数。应该是def dCdt(C,t,coefficients):。
至于参数优化，有很多不同的方法和途径。你在R中使用了什么方法？在 python 中找到良好参数的最简单方法之一是使用scipy.optimize.minimize 最小化误差平方和。
def SumSquaredErr(x):
k,alpha,beta = x
# Insert all your function definitions here #
C = odeint(dCdt, C_0, ts, args= (coefficients,))
return ((C[1:,0] - C1)**2).sum() + ((C[1:,1] - C2)**2).sum()
x0 = [k,alpha,beta] # your best guess for the parameter values
bounds = [(None,None),(0,None),(None,None)] # alpha>0 so V_nf!=0
res = minimize(SumSquaredErr, x0, bounds=bounds)
print('best [k,alpha,beta]: ', res.x)
我在这里假设 ts = np.linspace(0,60,61) 和 C1 和 C2 是时间 ts[1:] 的数据，即 1,2,3,...,59,60。
【讨论】：
非常感谢。我在 R 中使用 library(deSolve) # 求解微分方程的库 library(minpack.lm) # 使用 levenberg-marquardt 算法进行最小二乘拟合的库 @Galois 。在这种情况下，您可以将scipy.optimize.least_squares 与method='lm' 一起使用，即“在MINPACK 中实现的Levenberg-Marquardt 算法”，因此实际上与您在R 中使用的算法相同。偏最小二乘法的岭回归分析
...术A岭回归分析是一种修正的最小二乘估计法，当自变量系统中存在多重相关性时，它可以提供一个比最小二乘法更为稳定的估计，并且回归系数的标准差也比最小二乘估计的要小。根据高斯——马尔科夫定理，多重相关性并不... 查看详情
广义最小二乘估计是啥？
...计能获得较精确的结果。 假设所讨论的单输入单输出系统的差分方程模型是式中uk和yk分别是输入和输出序列：和是算子多项式,它们的系数是需要通过估计来求出的未知数;z-1是单位延迟算子；ek是误差序列，它是零均值平稳... 查看详情
python中的ODE求解系统；在本次通话中完成的多余工作
】python中的ODE求解系统；在本次通话中完成的多余工作【英文标题】：SolvingsystemofODEsinpython;excessworkdoneonthiscall【发布时间】：2021-06-1806:30:42【问题描述】：我正在尝试在python中解决一个耦合ODE系统以获得不同的潜力。它适用于... 查看详情
通过scipy和numpy使用python将数据拟合到ode系统(代码片段)
我在通过Scipy＆Numpy将我的MATLAB代码翻译成Python时遇到了一些麻烦。我坚持如何为我的ODE系统找到最佳参数值（k0和k1）以适应我的十个观察数据点。我目前对k0和k1有一个初步猜测。在MATLAB中，我可以使用一种叫做“fminsearch”的... 查看详情
用matlab估计armax模型参数的问题
1、使用ARMAX函数估计出来的结果是下面这样的形式：A(q)=1-1.284q^-1+0.7177q^-2C(q)=1-1.191q^-1+0.6625q^-2请问为什么没有常数项的估计结果？2、使用ARMAX估计ARMA（2,2）和使用garchfit函数估计ARMA(2,2)/GARCH（0,0）效果是一样的吗？3、如果想要... 查看详情
在python中生成和解决同时的ode(代码片段)
我对Python比较陌生，在编写一段生成然后解决微分方程系统的代码时遇到了一些问题。我这样做的方法是在函数var（）的列表中重新创建一组变量和系数，（x0，x1，...，xn）和（c0，c1，...，cn）。然后在EOM1（）中构造方程。初... 查看详情
使用 ODE 系统进行参数优化
...meterOptimisationwithsystemofODEs【发布时间】：2022-01-0711:11:34【问题描述】：我有一对ODE，我目前正试图将它们拟合到我拥有的一个小数据集，但是我在优化两个参数（a和c）时遇到了一些问题。ODE采用稍微改变的Lotka-Volterra形式，由... 查看详情
利用python进行回归分析(代码片段)
通常大家会认为曲线拟合和回归分析类似，但其实回归分析中是包含曲线拟合的。拟合是研究因变量和自变量的函数关系的。而回归是研究随机变量间的相关关系的。拟合侧重于调整参数，使得与给出的数据相符合。而... 查看详情
ls信道估计,mmse信道估计以及cs信道估计算法的误码率对比仿真(代码片段)
...的过程。如果信道是线性的话，那么信道估计就是对系统冲激响应进行估计。需强调的是信道估计是信道对输入信号影响的一种数学表示，而“好”的信道估计则是使得某种估计误差最小化的估计算法。无线通信系统的... 查看详情
如何在 R 中求解具有时间相关参数的 ODE 系统？
...DEwithtimedependentparametersinR?【发布时间】：2021-12-1821:47:09【问题描述】：我正在尝试通过deSolve求解这个ODE系统，dX/dt=-X*a+(YX)b+c和dY/dt=-Ya+(XY)*b对于时间[0,200]，a=0.30，b=0.2但对于时间[50,70]，c为1，否则为0。 查看详情
在python（odeint）中求解具有时间相关系数的常微分方程
】在python（odeint）中求解具有时间相关系数的常微分方程【英文标题】：solveordinarydifferentialequationwithtimedependentcoefficientsinpython(odeint)【发布时间】：2017-11-1813:03:13【问题描述】：我想使用scipy的odeint函数求解具有15个时间相关系... 查看详情
使用最小平方的非线性 ODE 优化
...linearODEsoptimizationwithleastsq【发布时间】：2021-11-3006:28:55【问题描述】：[更新]我正在研究非线性ODE系统优化并将其拟合到实验数据。我有一个由5个模型ODE组成的系统，必须通过17个参数进行优化。我的方法是计算求解的ODE和实验... 查看详情
频率学派贝叶斯学派估计的区别
这里的频率学派，认为参数θ是一个常量 ，只有属于置信区间，或者∉置信区间，没有属于这个某个置信区间的概率是0.9的说法。第一个意思是整体分布的一个参数θ，取θ的某一个先验分布，计算在该先验分... 查看详情
python中具有无限初始条件的ODE
】python中具有无限初始条件的ODE【英文标题】：ODEswithinfiniteinitlalconditioninpython【发布时间】：2013-01-0318:56:45【问题描述】：我有一个二阶微分方程，我想在python中求解它。问题是对于其中一个变量，我在0中没有初始条件，而只... 查看详情
在回归分析中,部分系数没有通过显著性检验,该如何处理?
...是否序列自相关，如果严重偏离2，则认为存在序列相关问题。F统计值是衡量回归方程整体显著性的假设检验，越大越显著。追问您好！如果大于0.05的变量该咋办？ 查看详情
将最大似然估计的系数放入观星表中
...imumlikelihoodintoastargazertable【发布时间】：2014-02-1519:05:00【问题描述】：Stargazer为lm（和其他）对象生成了非常好的乳胶表。假设我已经通过最大似然拟合了一个模型。我希望stargazer为我的估计生成一个类似lm的表。我该怎么做？... 查看详情
如何使用 Python 估计幂律分布的指数？
】如何使用Python估计幂律分布的指数？【英文标题】：HowtoestimatetheexponentofapowerlawdistributionusingPython?【发布时间】：2012-01-2818:24:25【问题描述】：我目前正在做一些网络分析，我想估计平均聚类系数与节点度的幂律分布的指数。... 查看详情
ode仿真引擎使用
使用ODE，必须记住参数，ERP(误差减少参数)和CFM(约束力混合)。1.ERP(误差减少参数) ERP是一个修正联合误差的参数。随着仿真的不断重复，由于各种内部逼近，使得关节误差不断增大。ERP可以修复错误。该值必须... 查看详情