如何迭代含5%随机误差的微分方程图模型100次并收集数据？

不用每次创建新模型，只需将单次模拟+参数拟合的逻辑封装成可复用函数，再循环指定次数（比如100次），用数组统一存储每次的结果即可。以下是具体实现：

步骤说明 #

• 封装单次模拟流程：生成带5%随机误差的数据、运行曲线拟合，返回需要的结果（拟合参数、协方差矩阵、带误差的数据等）
• 初始化存储容器：用NumPy数组预先分配空间，高效存储多次迭代的结果
• 循环迭代：调用封装好的函数，批量完成100次模拟与数据收集

修改后的完整代码 #

import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate, optimize

def dose(t, y, b, s, c, p, d):
    target, infectious, virus = y
    return np.array([
        -b * target * virus,
        b * target * virus - s * infectious,
        (1. / (d + 1.)) * p * infectious - c * virus
    ])

def model(D, b, s, c, p):
    solutions = []
    for d in D:
        solution = integrate.solve_ivp(
            dose, [0, 5], y0=[1, 0, 0.01],
            t_eval=[2.8828828828828827],
            args=(b, s, c, p, d)
        )
        data = solution.y[2, 0] / 0.01950269536785707  
        solutions.append(data)
    return np.array(solutions)

# 封装单次模拟与拟合的函数
def run_single_simulation(true_params, D):
    b_true, s_true, c_true, p_true = true_params
    # 生成无误差的模型输出
    z = model(D, b_true, s_true, c_true, p_true)
    # 生成5%随机误差（修正：用np.random.normal生成数组，对应每个数据点的误差）
    sigma = z * 0.05
    noise = np.random.normal(0, sigma, size=len(z))
    zn = z + noise
    # 运行曲线拟合
    popt, pcov = optimize.curve_fit(
        model, D, zn, p0=[1e-5, 1, 1, 1e6],
        method="trf", bounds=(0, np.inf),
        sigma=sigma, absolute_sigma=True
    )
    return popt, pcov, zn

# 初始化参数与输入
true_params = (0.00001, 4, 4, 2000000)
D = np.logspace(-3, 3, 7)
n_iterations = 100

# 初始化存储结果的数组
all_popt = np.zeros((n_iterations, 4))  # 存储每次的拟合参数
all_pcov = np.zeros((n_iterations, 4, 4))  # 存储每次的协方差矩阵
all_zn = np.zeros((n_iterations, len(D)))  # 存储每次带误差的模拟数据

# 循环迭代100次
for i in range(n_iterations):
    popt, pcov, zn = run_single_simulation(true_params, D)
    all_popt[i] = popt
    all_pcov[i] = pcov
    all_zn[i] = zn

# 示例：查看拟合参数的统计结果
print("拟合参数的均值：")
print(np.mean(all_popt, axis=0))
print("\n拟合参数的标准差：")
print(np.std(all_popt, axis=0))

# 可选：可视化某次的拟合结果（比如第0次）
Dlog = np.logspace(-3, 3, 200)
fig, axe = plt.subplots()
axe.scatter(D, all_zn[0], label="带误差数据")
axe.semilogx(Dlog, model(Dlog, *all_popt[0]), label="拟合曲线")
axe.semilogx(Dlog, model(Dlog, *true_params), "--", label="真实模型")
axe.legend()
axe.grid()
plt.show()

关键细节说明 #

1. 修正误差生成逻辑：原代码中random.gauss只能生成单个值，改用np.random.normal生成与数据点数量匹配的误差数组，确保每个数据点都添加对应比例的随机噪声。
2. 变量名冲突修复：原代码中s同时表示模型参数和误差标准差，修改为sigma避免覆盖。
3. 高效数据存储：用NumPy数组预先分配空间，比列表追加更高效，尤其适合大规模迭代。
4. 可扩展性：如果需要收集其他数据（比如拟合的残差），只需在run_single_simulation函数中添加返回值，并扩展存储数组即可。

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