ODE系统求解中For循环结果存入列表的问题排查

存在的错误点 #

• 未初始化存储容器：acc和rejected列表没有在循环前定义，导致rejected.append()会触发未定义错误；且acc在每次满足条件时被重新赋值，无法累积存储所有符合条件的结果。
• 错误的存储逻辑：
  • acc被赋值为包含重复布尔判断结果的列表，没有存储需要的beta、gamma、最终患病率和发病率数值。
  • rejected.append(beta_samples)会把整个beta_samples数组添加进去，而非当前被拒绝的beta和gamma参数。
• 冗余的条件判断：acc.append(320 < I[-1]*100000 < 480)重复了相同的条件判断，没有实际意义。

修正后的代码 #

import numpy as np
from scipy.integrate import odeint

# 初始化存储列表
accepted = []
rejected = []

beta_samples = np.random.uniform(0, 30, 50)
gamma_samples = np.random.uniform(0, 2, 50)

for beta, gamma in zip(beta_samples, gamma_samples):
    # Total population, N.
    N = 1
    # Initial number of infected and recovered individuals, I0 and R0.
    I0, R0 = 0.001, 0
    # Everyone else, S0, is susceptible to infection initially.
    U0 = N - I0 - R0
    J0 = I0
    Lf0, Ls0 = 0, 0
    mu, muTB, sigma, rho = 1/80, 1/6, 1/6, 0.03
    u, v, w = 0.88, 0.083, 0.0006
    t = np.linspace(0, 500, 500+1)
    
    # The SIR model differential equations.
    def deriv(y, t, N, beta, gamma, mu, muTB, sigma, rho, u, v, w):
        U, Lf, Ls, I, R, cInc = y
        b = (mu * (U + Lf + Ls + R)) + (muTB * I)
        lamda = beta * I
        clamda = 0.2 * lamda
        dU = b - ((lamda + mu) * U)
        dLf = (lamda*U) + ((clamda)*(Ls + R)) - ((u + v + mu) * Lf)
        dLs = (u * Lf) - ((w + clamda + mu) * Ls)
        dI = w*Ls + v*Lf - ((gamma + muTB + sigma) * I) + (rho * R)
        dR = ((gamma + sigma) * I) - ((rho + clamda + mu) * R)
        cI = w*Ls + v*Lf + (rho * R)
        return dU, dLf, dLs, dI, dR, cI
    
    # Integrate the SIR equations over the time grid, t.
    solve = odeint(deriv, (U0, Lf0, Ls0, I0, R0, J0), t, args=(N, beta, gamma, mu, muTB, sigma, rho, u, v, w))
    U, Lf, Ls, I, R, cInc = solve.T

    # 计算最终患病率和发病率（转换为每10万人口）
    final_prevalence = I[-1] * 100000
    final_incidence = (cInc[1:] - cInc[:-1])[-1] * 100000

    if 320 < final_prevalence < 480 and 240 < final_incidence < 360:
        # 存储符合条件的参数和计算值
        accepted.append({
            'beta': beta,
            'gamma': gamma,
            'final_prevalence': final_prevalence,
            'final_incidence': final_incidence
        })
        print(f"接受 beta={beta:.2f}, gamma={gamma:.2f} | 患病率: {final_prevalence:.2f}, 发病率: {final_incidence:.2f}")
    else:
        # 存储被拒绝的参数
        rejected.append({'beta': beta, 'gamma': gamma})
        print(f"拒绝 beta={beta:.2f}, gamma={gamma:.2f}")

# 查看结果
print("\n接受的参数组数量:", len(accepted))
print("拒绝的参数组数量:", len(rejected))

修正说明 #

1. 初始化存储列表：在循环前定义accepted和rejected列表，用于累积所有符合/不符合条件的结果。
2. 优化存储结构：使用字典存储每组参数和对应的计算值，方便后续查看和处理。
3. 避免重复计算：将最终患病率和发病率提前计算并赋值给变量，简化条件判断和代码可读性。
4. 修正拒绝逻辑：将当前的beta和gamma添加到rejected列表，而非整个样本数组。

内容的提问来源于stack exchange，提问作者Landon