我希望对以下ODE拟合一条曲线:
dA/dt = k1*profit + k2
我有变量A
和profit
的观测时间序列,并希望使用Python中的曲线拟合技术获取k1
和k2
的最优值。我可以编写下面的代码进行操作,但是解决方案拟合效果不好,或者我的方法可能有误。
import numpy as np
from scipy.optimize import curve_fit
from scipy.integrate import odeint
def fitfunc(t, k1, k2):
'Function that returns A computed from an ODE for k1 and k2'
def myode(area, t):
profit = get_profit(t)
return k1*profit + k2
A0 = 200000 #initial value of A
out = odeint(myode, A0, t)
return out[:,0]
k_fit, kcov = curve_fit(fitfunc, time_span, obs_A) #time span from 1999-2019 and obs_A is the observed values of A
modeled_A = fitfunc(time_span, k_fit[0], k_fit[1])
这20年期间的利润和obs_A数据如下:
profit = [ 7.65976374e+06, -6.13172279e+06, 1.03946093e+07, 2.59937877e+06,
-7.88358386e+06, -1.38918115e+04, -3.13403157e+06, -4.74348806e+06,
1.87296164e+07, 4.13680709e+07, -1.77191198e+07, 2.39249499e+06,
1.38521564e+07, 6.52548348e+07, -5.78102494e+07, -5.72469988e+07,
-5.99056006e+06, -1.72424523e+07, 1.78509987e+07, 9.27860105e+06,
-9.96709853e+06]
obs_A = [200000., 165000., 150000., 180000., 190000., 195000., 200000.,
165000., 280000., 235000., 250000., 250000., 250000., 295000.,
295000., 285000., 245000., 315000., 235000., 245000., 305000.]
time_span = np.arange(1999,2020)
这里的get_profit
是一个函数,它会在给定的t
时刻输出利润值。该函数是通过以下方式对观测到的profit
数据进行插值创建的:
profit_fun = interp1d(t, profit.values, 1, fill_value="extrapolate")
def get_profit(t):
return profit_fin(t)
我不确定如何在这里使用profit
变量,因为它会在每个时间步骤发生变化。我的方法正确吗?
t
的插值(我不知道它有什么用处,但我不懂ODEs),结果可能会好一些,但看起来像是平滑曲线。从这里开始,我不知道该怎么做,抱歉。 - K.Cl