基于Milstein格式的随机SIR类ODE系统实现问题排查

我需要为Ramos等人2021年提出的类SIR模型ODE系统添加噪声，采用Milstein格式实现了相关方程，但绘制易感者、感染者等变量结果时，图形异常，与确定性模型完全不符，且无噪声波动。代码如下：

# Create Brownian Motion
np.random.seed(1)
dS = np.sqrt(dt) * np.random.randn(tmax)
dE=np.sqrt(dt) * np.random.randn(tmax)
dI=np.sqrt(dt) * np.random.randn(tmax)
dIu=np.sqrt(dt) * np.random.randn(tmax)
dDu=np.sqrt(dt) * np.random.randn(tmax)
dHR=np.sqrt(dt) * np.random.randn(tmax)
dHD=np.sqrt(dt) * np.random.randn(tmax)
dB=[dS,dE,dI,dIu,dDu,dHR,dHD]
sigma=[0.5,0,0,0,0,0,0]
#Brief definition of systemf function evaluating the second term at time t for each variant i
for i in range(nvariants):
    newe = S[0]*(mbetae[i]*E[i] + mbetai[i]*I[i] + mbetaiu[i]*Iu[i] + mbetahr[i]*HR[i] + mbetahd[i]*HD[i])/totalpop
    newi = gammae * E[i]
    newhid = gammai  * I[i]
    newhiu = gammaiu * Iu[i]
    newr= gammahr * HR[i]
    newd = gammahd*HD[i]
    newq = gammaq  * Q[i]
    neweS = neweS + newe
    fE[i] = newe-newi
    fI[i] = newi - newhid
    fIu[i]= (1-theta[i]-omegau)*newhid - newhiu
    fHR[i]= p[i]*(theta[i]-fatrate[i])*newhid - newr
    fHD[i]= fatrate[i] *newhid - newd
    fQ[i] = (1-p[i])*(theta[i]-fatrate[i])*newhid + newr - newq
fS[0] =  - neweS -(vjRK[int(mt.floor(t))]) #fS.insert(0,-neweS -(vjRK[mt.floor(t)]))
return [fS, fE, fI, fIu, fHR, fHD, fQ]
for t in range(delayini,tmax-1):
      fsyseval=systemf(t,states[t],beta[2*t],gamma[2*t],frate[2*t],theta[2*t],p[2*t],omegau[2*t],vjsum)
#Running the scheme
for s in range(numstates):
        for i in range(nvariants):
            states[t+1][s][i] =states[t][s][i]+fsyseval[s][i]*dt+sigma[s]*dB[s][t]*states[t][s][i] + 0.5*sigma[s]**2 * states[t][s][i] * (dB[s][t] ** 2 - dt)

问题排查 #

1. 核心执行逻辑错误 #

代码中的return [fS, fE, fI, fIu, fHR, fHD, fQ]语句直接终止了函数，导致后续的Milstein迭代循环完全没有执行——相当于你只定义了确定性项的计算逻辑，但根本没跑随机迭代，自然不会有噪声波动，甚至确定性模型的迭代也可能未正确完成。

2. 噪声配置失效 #

你的sigma数组仅给S变量设置了0.5的噪声强度，其余变量全为0，这意味着只有S变量能被添加噪声。如果连S变量都没出现波动，就是上述执行逻辑错误导致的。

3. 其他潜在问题 #

• neweS未初始化就进行累加，会导致fS[0]计算错误，直接让确定性部分结果异常。
• 布朗运动数组dB的维度与states三维结构（时间步×状态变量×变异株）的匹配性未明确，可能出现索引错误。
• 未明确噪声类型：当前实现的是乘性噪声（噪声项与当前状态值成正比），若模型要求加性噪声，代码逻辑完全错误。

修正方案 #

1. 修复代码执行结构 #

将systemf定义为独立函数，把迭代逻辑放在函数外部执行，确保Milstein格式的代码能被正确运行：

# 定义确定性项计算函数
def systemf(t, states_t, beta_t, gamma_t, frate_t, theta_t, p_t, omegau_t, vjsum, totalpop, gammae, gammai, gammaiu, gammahr, gammahd, gammaq, mbetae, mbetai, mbetaiu, mbetahr, mbetahd):
    S = states_t[0]
    E = states_t[1]
    I = states_t[2]
    Iu = states_t[3]
    HR = states_t[4]
    HD = states_t[5]
    Q = states_t[6]
    nvariants = len(E)
    
    neweS = 0  # 必须初始化累加变量
    fE = np.zeros(nvariants)
    fI = np.zeros(nvariants)
    fIu = np.zeros(nvariants)
    fHR = np.zeros(nvariants)
    fHD = np.zeros(nvariants)
    fQ = np.zeros(nvariants)
    fS = np.zeros(1)
    
    for i in range(nvariants):
        newe = S[0]*(mbetae[i]*E[i] + mbetai[i]*I[i] + mbetaiu[i]*Iu[i] + mbetahr[i]*HR[i] + mbetahd[i]*HD[i])/totalpop
        newi = gammae * E[i]
        newhid = gammai  * I[i]
        newhiu = gammaiu * Iu[i]
        newr= gammahr * HR[i]
        newd = gammahd*HD[i]
        newq = gammaq  * Q[i]
        
        neweS += newe
        fE[i] = newe - newi
        fI[i] = newi - newhid
        fIu[i] = (1 - theta_t[i] - omegau_t[i])*newhid - newhiu
        fHR[i] = p_t[i]*(theta_t[i] - frate_t[i])*newhid - newr
        fHD[i] = frate_t[i] * newhid - newd
        fQ[i] = (1 - p_t[i])*(theta_t[i] - frate_t[i])*newhid + newr - newq
    
    fS[0] = - neweS - vjsum[int(np.floor(t))]
    return [fS, fE, fI, fIu, fHR, fHD, fQ]

# 初始化参数（需根据你的模型补充完整）
np.random.seed(1)
dt = 0.01  # 明确时间步长
tmax = 1000  # 明确总时间步数
totalpop = 1e6  # 总人口数
nvariants = 2  # 变异株数量
numstates = 7  # 状态变量数：S,E,I,Iu,HR,HD,Q

# 初始化布朗运动
dS = np.sqrt(dt) * np.random.randn(tmax)
dE = np.sqrt(dt) * np.random.randn(tmax)
dI = np.sqrt(dt) * np.random.randn(tmax)
dIu = np.sqrt(dt) * np.random.randn(tmax)
dDu = np.sqrt(dt) * np.random.randn(tmax)
dHR = np.sqrt(dt) * np.random.randn(tmax)
dHD = np.sqrt(dt) * np.random.randn(tmax)
dB = [dS, dE, dI, dIu, dDu, dHR, dHD]

# 配置噪声强度（根据模型需求调整，不要全为0）
sigma = [0.5, 0.1, 0.1, 0.1, 0.05, 0.05, 0.1]

# 初始化状态数组（补充你的初始条件）
states = np.zeros((tmax, numstates, nvariants))
states[0, 0, :] = totalpop * 0.9  # 初始易感者比例
states[0, 2, 0] = totalpop * 0.01  # 初始感染者数量

# 运行Milstein迭代
delayini = 0
for t in range(delayini, tmax-1):
    # 获取当前时间步的模型参数
    beta_t = beta[2*t]
    gamma_t = gamma[2*t]
    frate_t = frate[2*t]
    theta_t = theta[2*t]
    p_t = p[2*t]
    omegau_t = omegau[2*t]
    vjsum_t = vjRK[int(np.floor(t))]
    
    # 计算确定性项
    fsyseval = systemf(t, states[t], beta_t, gamma_t, frate_t, theta_t, p_t, omegau_t, vjsum_t, totalpop, gammae, gammai, gammaiu, gammahr, gammahd, gammaq, mbetae, mbetai, mbetaiu, mbetahr, mbetahd)
    
    # 更新状态变量
    for s in range(numstates):
        for i in range(nvariants):
            # Milstein格式：确定性项 + 一阶随机项 + 二阶随机项
            states[t+1][s][i] = states[t][s][i] + fsyseval[s][i]*dt \
                                + sigma[s] * dB[s][t] * states[t][s][i] \
                                + 0.5 * sigma[s]**2 * states[t][s][i] * (dB[s][t]**2 - dt)

2. 验证与调整 #

• 噪声有效性测试：先将sigma仅给S变量设非零值，运行后查看S曲线是否有波动，确认随机项生效。
• 噪声类型修正：若模型要求加性噪声，将随机项改为sigma[s]*dB[s][t]，去掉与states[t][s][i]的乘积。
• 数值稳定性检查：人口模型变量不能为负，若出现负数，需减小sigma或dt，或考虑将负数截断为0（需谨慎评估合理性）。

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