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

发布于 2025-01-18 16:01:13 字数 1068 浏览 4 评论 0原文

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）实际上我什至不确定是否应该以这种方式推断参数。

谁能帮我解决这个问题，有没有更好的推断参数的方法。提前致谢！

原文

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)

So X is my data, the first column is species x, and second column is species y.
I tried to infer parameters for this Lotka-Volterra model using ode and curve fit.
The error says not enough values to unpack (expected 2, got 1)
I am actually not even sure whether I should infer parameter this way.

Can anyone help me with this, are there any better methods of infering parameters.
Thanks in advance!

分享到QQ

分享到微博

如果你对这篇内容有疑问，欢迎到本站社区发帖提问参与讨论，获取更多帮助，或者扫码二维码加入 Web 技术交流群。

发布评论

需要登录才能够评论，你可以免费注册一个本站的账号。

乱世争霸 2025-01-25 16:01:14

请注意，在中，scipy.optimize.curve_fit actuemnts xdata和ydata ydata 是平坦的阵列，您使用多维结构。虽然强烈建议XDATA包含一个输入值或每个元素的ydata的向量，但不需要它。 XDATA是一个常数，也可以通过其他方式传递。它仅是为了方便标准回归任务，因此只需更改这两个阵列就可以在一个位置实现输入数据的完整更改。

因此，拥有ydata是xdata的两倍也没问题。只需将.flatten（）应用于二维数组即可。

接下来，参数列表必须是标量列表，因此添加y0并传递初始矢量[x0，y0]。

这些更正和对代码的更改一起导致“结果”。这不是很有说服力。

我得到了更好的结果，但仍然不完美，在使用多个拍摄方法中，在x [： - 1]中获得积分，并集成为时间步长1，将收集的终点列表与x [1：]进行比较。这可以更好地查找与幅度和频率相匹配的参数，但会产生轻微的速度差异，通过对系数进行3％的校正，看起来更好。

人们可能需要两种方法的混合物来获得尊重本地和全球特征。

确实，这种组合的方法起作用，给出参数

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

X0,Y0 =  34.53872014350661, 4.653177640949391

代码为该组合的拟合程序：用于剩余计算，首先封装求解器以避免重复求解器参数。然后使用它首先在整个间隔上与可变初始点进行集成，然后在时间步骤1段中进行集成。

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]

显然，

params, info = curve_fit(integration,XData,YData,p0=p0)
XY0, para = params[4:], params[:4]
print(XY0,tuple(para))

t_plot = np.linspace(0,len(X),500)
x_plot = solver(XY0, t_plot, tuple(para))

Note that in scipy.optimize.curve_fitboth arguemnts xdata and ydata are required to be flat arrays, you use a multi-dimensional structure. While it is strongly suggested that xdata contains one input value or vector per element of ydata, there is no requirement for it. xdata is a constant that could also have been passed some other way. It is there just for the convenience in standard regression tasks, so a complete change in the input data can be effected in one location by just changing these two arrays.

Thus it is also no problem to have ydata twice as long as xdata. Just apply .flatten() to the 2-dimensional arrays.

Next, the parameter list has to be a list of scalars, so add Y0 and pass the initial vector [X0,Y0].

Together these corrections and changes to your code lead to a "result". Which is not very convincing.

I got a better result, but still not perfect, in using a multiple shooting approach, taking the points in X[:-1] and integrating for the time step 1, comparing the collected list of end-points to X[1:]. This works better in finding parameters that match amplitude and frequency, but produces a slight speed difference that looks better with a 3% correction of the coefficients.

One would probably need a mix of both approaches to get the local as well as global characteristics respected.

And indeed, this combined approach works, giving parameters

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

X0,Y0 =  34.53872014350661, 4.653177640949391

Code for that combined fitting program: For the residual computation, first encapsulate the solver to avoid repetition of solver parameters. Then use that to first integrate over the full interval with the variable initial point, and then over the time step 1 segments.

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()

This obviously needs the reference array prepared in a similar fashion

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

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

After that the curve fitting procedure call remains the same, all changes happened before

params, info = curve_fit(integration,XData,YData,p0=p0)
XY0, para = params[4:], params[:4]
print(XY0,tuple(para))

t_plot = np.linspace(0,len(X),500)
x_plot = solver(XY0, t_plot, tuple(para))

回复收藏 0 原文

~没有更多了~