scipy SIR模型耦合常微分方程的积分和拟合

cvxl0en2 于 2022-11-10 发布在其他

关注(0)|答案(2)|浏览(186)

在这种情况下，有3个常微分方程描述SIR模型。问题来了，我想从x_axis和y_axis值中计算哪些beta和gamma值最适合数据点。我目前使用的方法是使用scipy库中的odeint和同一库中的curve_fit方法对常微分方程进行积分。在这种情况下，如何计算beta和gamma的值以拟合数据点？
P.S.当前的错误是：ValueError: operands could not be broadcast together with shapes (3,) (14,)


# initial values

S_I_R = (0.762/763, 1/763, 0)

x_axis = [m for m in range(1,15)]
y_axis = [3,8,28,75,221,291,255,235,190,125,70,28,12,5]

# ODE's that describe the system

def equation(SIR_Values,t,beta,gamma):
    Array = np.zeros((3))
    SIR = SIR_Values
    Array[0] = -beta * SIR[0] * SIR[1]
    Array[1] = beta * SIR[0] * SIR[1] - gamma * SIR[1]
    Array[2] = gamma * SIR[1]
    return Array

# Results = spi.odeint(equation,S_I_R,time)

# fitting the values

beta_values,gamma_values = curve_fit(equation, x_axis,y_axis)

scipy

来源：https://stackoverflow.com/questions/73838531/integrating-and-fitting-coupled-odes-for-sir-modelling

2条答案

按热度按时间

xzabzqsa1#


# Starting values

S0 = 762/763
I0 = 1/763
R0 = 0

x_axis = np.array([m for m in range(0,15)])
y_axis = np.array([1,3,8,28,75,221,291,255,235,190,125,70,28,12,5])
y_axis = np.divide(y_axis,763)

def sir_model(y, x, beta, gamma):
    S = -beta * y[0] * y[1] 
    R = gamma * y[1]
    I = beta * y[0] * y[1] - gamma * y[1]
    return S, I, R

def fit_odeint(x, beta, gamma):
    return spi.odeint(sir_model, (S0, I0, R0), x, args=(beta, gamma))[:,1]

popt, pcov = curve_fit(fit_odeint, x_axis, y_axis)
beta,gamma = popt
fitted = fit_odeint(x_axis,*popt)
plt.plot(x_axis,y_axis, 'o', label = "infected per day")
plt.plot(x_axis, fitted, label = "fitted graph")
plt.xlabel("Time (in days)")
plt.ylabel("Fraction of infected")
plt.title("Fitted beta and gamma values")
plt.legend()
plt.show()

赞(0）回复(0）举报 2022-11-10

x33g5p2x2#

作为scipy文档中的this example，function必须输出与x_axis和y_axis大小相同的数组。

赞(0）回复(0）举报 2022-11-10