尝试推断Lotka-Volterra模型的模型参数

Question

def derivative(X, t, A, B, C, D):
    x, y = X
    dotx = x * (A - B * y)
    doty = y * (-D + C * x)
    return np.array([dotx, doty])

def integration(t,A,B,C,D,X0):
    res = odeint(derivative, X0, t, args = (A,B,C,D))
    return res

X0 = [30, 4]

X = array([[30. ,  4. ],
       [47.2,  6.1],
       [70.2,  9.8],
       [77.4, 35.2],
       [36.3, 59.4],
       [20.6, 41.7],
       [18.1, 19. ],
       [21.4, 13. ],
       [22. ,  8.3],
       [25.4,  9.1],
       [27.1,  7.4],
       [40.3,  8. ],
       [57. , 12.3],
       [76.6, 19.5],
       [52.3, 45.7],
       [19.5, 51.1],
       [11.2, 29.7],
       [ 7.6, 15.8],
       [14.6,  9.7],
       [16.2, 10.1],
       [24.7,  8.6]])

t = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0, 20.0]

XData = t
YData = X

curve_fit(integration,XData,YData)

所以 X 是我的数据，第一列是物种 x，第二列是物种 y。 我尝试使用颂歌和曲线拟合来推断该 Lotka-Volterra 模型的参数。 该错误表明没有足够的值来解压（预期为 2，实际为 1） 我实际上什至不确定是否应该以这种方式推断参数。

任何人都可以帮我解决这个问题吗？有没有更好的推断参数的方法。 预先感谢！

Answer 1

请注意，在

scipy.optimize.curve_fit

中，

xdata

和

ydata

都需要是平面数组，您使用多维结构。虽然强烈建议

xdata

的每个元素包含一个输入值或向量，但没有任何要求。

ydata

是一个常数，也可以通过其他方式传递。它的存在只是为了方便标准回归任务，因此只需更改这两个数组即可在一个位置实现输入数据的完全更改。因此，

xdata

的长度是

ydata

的两倍也没有问题。只需将

xdata

应用于二维数组即可。接下来，参数列表必须是标量列表，因此添加

.flatten()

并传递初始向量

Y0

。对代码的这些更正和更改共同产生了“结果”。这不太有说服力。

我得到了更好的结果，但仍然不完美，使用多重拍摄方法，获取

[X0,Y0]

中的点并积分时间步长

X[:-1]

，将收集的端点列表与

1

进行比较。这可以更好地找到与幅度和频率相匹配的参数，但会产生轻微的速度差异，通过对系数进行 3% 校正后效果会更好。人们可能需要结合使用两种方法来尊重本地和全球特征。

事实上，这种组合方法是有效的，给出了参数

X[1:]

---

组合拟合程序代码：

对于残差计算，首先封装求解器，以避免求解器参数重复。然后使用它首先在整个区间上与可变初始点积分，然后在时间步 1 段上积分。

A,B = 0.5215206964006734,  0.02567364947581818
C,D = 0.02493663631623848, 0.8476224408838039

X0,Y0 =  34.53872014350661, 4.653177640949391

这显然需要以类似的方式准备参考数组

def solver(XY,t,para): 
    return odeint(derivative, XY, t, args = para, atol=1e-8, rtol=1e-11)

def integration(XY_arr,*para):
    XY0 = para[4:]
    para = para[:4]
    T = np.arange(len(XY_arr))
    res0 = solver(XY0,T, para)
    res1 = [ solver(XY,[t,t+1],para)[-1] 
             for t,XY in enumerate(XY_arr[:-1]) ]
    return np.concatenate([res0,res1]).flatten()

之后曲线拟合过程调用保持不变，所有变化都发生在之前

XData = X
YData = np.concatenate([ X,X[1:]]).flatten()

p0 =[ 0.5215,  0.02567,  
      0.02493,  0.8476,
      34.53, 4.653]