如何适配SERIVHD模型

Question

使用 Covid Act Now 数据 (states.timeseries.csv)，我想开发一个预测性 SEIRVHD 模型。我有代码的主要代码架构（见下文），但我需要帮助。

1. 在 COVID-19 数据上改进健身/训练我的 SEIRVHD 模型
2. 根据 COVID 数据优化参数
3. 绩效评估

我也已经删除了 nan 值——用零填充它们，并执行填充转发。

此代码定义并求解 SEIRVHD 模型以预测 COVID 19 的传播。SEIRVHD 模型是 SEIR（易感-暴露-感染-恢复）模型的变体，增加了用于接种疫苗的个体 (V)、住院治疗 ( H)、ICU 入院 (ICU) 和死亡 (D)。

seirvhd_model 函数定义了控制疾病传播的微分方程，其输入包括模型参数，例如传播率 (beta)、潜伏期 (alpha)、恢复率 (gamma)、疫苗接种率 (delta) 和其他。

模型的初始条件和时间范围是根据 DataFrame df 中的数据定义的，solve_seirvhd 函数使用 scipy.integrate 中的 odeint 来求解微分方程并为模型的不同部分生成预测值。

model_residuals 函数计算实际数据与模型预测之间的差异，然后在 scipy.optimize 的 curve_fit 函数中使用它来优化模型参数以拟合数据。最后，代码生成一个图来比较实际案例和模型案例。

# Define SEIRVHD model
def seirvhd_model(y, t, alpha, beta, gamma, delta, eta, epsilon, theta, kappa, lambda_H, lambda_ICU, mu_H, mu_ICU):
    S, E, I, V, H, ICU, D, R = y
    dSdt = -(beta * I * S) / N
    dEdt = (beta * I * S) / N - alpha * E
    dIdt = alpha * E - gamma * I - delta * V
    dVdt = delta * V - eta * V
    dHdt = theta * I - (kappa + lambda_H + mu_H) * H
    dICUdt = kappa * H - (lambda_ICU + mu_ICU) * ICU
    dDdt = lambda_H * H + lambda_ICU * ICU
    dRdt = gamma * I + eta * V + mu_H * H + mu_ICU * ICU
    return dSdt, dEdt, dIdt, dVdt, dHdt, dICUdt, dDdt, dRdt

# Define initial conditions and time range
#Define model parameters
N = 32820000  # Total population of US
N, I0, V0, R0, E0, D0, H0, ICU0, S0 = (32820000, df.loc[df.index[0], "actuals.cases"],
df.loc[df.index[0], "actuals.vaccinationsInitiated"],
df.loc[df.index[0], "actuals.vaccinationsCompleted"], 0,
df.loc[df.index[0], "actuals.deaths"],
df.loc[df.index[0], "actuals.hospitalBeds.capacity"],
df.loc[df.index[0], "actuals.icuBeds.capacity"],
N - df.loc[df.index[0], "actuals.cases"] -
df.loc[df.index[0], "actuals.vaccinationsInitiated"] -
df.loc[df.index[0], "actuals.vaccinationsCompleted"],)

y0 = [S0, E0, I0, V0, H0, ICU0, D0, R0]

#Define time range to solve over (start and end dates of data)
start_date = pd.Timestamp('2020-03-15')
df = df[df['date'] >= start_date]
end_date = df["date"].iloc[-1]
t_range = [0, (end_date - start_date).days]

# Define function to solve the SEIRVHD model and compute the residuals
def solve_seirvhd(params):
    alpha, beta, gamma, delta, eta, epsilon, theta, kappa, lambda_H, lambda_ICU, mu_H, mu_ICU = params
    sol = odeint(seirvhd_model, y0, t, args=(alpha, beta, gamma, delta, eta, epsilon, theta, kappa, lambda_H, lambda_ICU, mu_H, mu_ICU))
    infected = sol


from scipy.optimize import curve_fit

# Define a function to calculate the difference between the actual data and model predictions
def model_residuals(t, *params):
    alpha, beta, gamma, delta, eta, epsilon, theta, kappa, lambda_H, lambda_ICU, mu_H, mu_ICU = params
    sol = odeint(seirvhd_model, y0, t, args=(alpha, beta, gamma, delta, eta, epsilon, theta, kappa, lambda_H, lambda_ICU, mu_H, mu_ICU))
    return (sol[:, 2] + sol[:, 6] - df['actuals.cases'] - df['actuals.deaths']).ravel()

# Define initial parameter values
p0 = (0.2, 1.75, 0.5, 0.03, 0.2, 0.01, 0.1, 0.05, 0.1, 0.05, 0.05, 0.1)

# Fit the model to the data using curve_fit
popt, pcov = curve_fit(model_residuals, df.index.values, np.zeros(len(df)), p0=p0)

# Print the optimized parameter values
print('Optimized parameter values:')
print('alpha:', popt[0])
print('beta:', popt[1])
print('gamma:', popt[2])
print('delta:', popt[3])
print('eta:', popt[4])
print('epsilon:', popt[5])
print('theta:', popt[6])
print('kappa:', popt[7])
print('lambda_H:', popt[8])
print('lambda_ICU:', popt[9])
print('mu_H:', popt[10])
print('mu_ICU:', popt[11])

import matplotlib.pyplot as plt

# Get model predictions for cases
model_cases = sol[:, 2]

# Plot actual cases and model predictions
plt.figure(figsize=(10, 6))
plt.plot(df['date'], df['actuals.cases'], label='Actual cases')
plt.plot(df['date'], model_cases[:len(df)], label='Modeled cases')
plt.legend()
plt.xlabel('Date')
plt.ylabel('Number of cases')
plt.title('Comparison of actual and modeled cases')
plt.show()