如何20次运行粒子模拟函数并计算各步平均值与误差？

需求说明 #

需要将已实现的粒子模拟函数run运行20次，对每一步的种群数量计算平均值，同时得到各步的误差（标准差或标准误）。

原模拟函数run #

def run(A, B, C, t1, t2, t3, time, steps):
    dt = time/steps
    
    p_A = 1/t1 * dt
    p_B = 1/t2 * dt
    p_C = 1/t3 * dt

    testrules_N = [
        ('A', 'B', p_A),
        ('B', 'C', p_B),
        ('C', 'A', p_C)
    ]

    testrules = [
        ('A', 'B', p_A),
        ('B', 'C', p_B)
    ]
    A1,B1,C1 = evolve_system(A, B, C, testrules_N, steps)
    A2,B2,C2 = evolve_system(A1[-1], B1[-1], C1[-1], testrules, steps)
    return A1, B1, C1, A2, B2, C2

原有代码问题 #

你提供的代码存在以下问题，无法实现需求：

• 仅调用了一次run函数，未完成20次循环运行
• 存储数组populations维度定义错误，未匹配run返回的6组时间序列数据
• 未正确指定numpy.average和numpy.std的计算轴，且未保存统计结果
• 赋值逻辑混乱，无法将每次运行的结果正确存入数组

修正后的代码 #

1. 导入依赖与定义参数 #

import numpy as np

# 替换为你的实际参数值
A0, B0, C0 = 1000, 0, 0
t_half_A, t_half_B, t_half_C = 10, 20, 30
t_total = 50
nsteps = 100
n_runs = 20  # 运行次数设置为20

2. 初始化结果存储数组 #

run返回6组时间序列（A1/B1/C1/A2/B2/C2），每组长度为nsteps + 1（包含初始状态），因此创建对应维度的数组存储所有运行结果：

# 初始化存储数组，形状为(运行次数, 步数+1)
A1_all = np.zeros((n_runs, nsteps + 1))
B1_all = np.zeros((n_runs, nsteps + 1))
C1_all = np.zeros((n_runs, nsteps + 1))
A2_all = np.zeros((n_runs, nsteps + 1))
B2_all = np.zeros((n_runs, nsteps + 1))
C2_all = np.zeros((n_runs, nsteps + 1))

3. 循环运行模拟 #

for i in range(n_runs):
    # 单次运行模拟
    A1, B1, C1, A2, B2, C2 = run(A0, B0, C0, t_half_A, t_half_B, t_half_C, t_total, nsteps)
    # 保存本次运行结果
    A1_all[i] = A1
    B1_all[i] = B1
    C1_all[i] = C1
    A2_all[i] = A2
    B2_all[i] = B2
    C2_all[i] = C2

4. 计算平均值与误差 #

# 计算每一步的平均值（axis=0表示对所有运行次数取平均）
A1_mean = np.mean(A1_all, axis=0)
B1_mean = np.mean(B1_all, axis=0)
C1_mean = np.mean(C1_all, axis=0)
A2_mean = np.mean(A2_all, axis=0)
B2_mean = np.mean(B2_all, axis=0)
C2_mean = np.mean(C2_all, axis=0)

# 计算每一步的标准差（衡量种群数量的离散程度）
A1_std = np.std(A1_all, axis=0)
B1_std = np.std(B1_all, axis=0)
C1_std = np.std(C1_all, axis=0)
A2_std = np.std(A2_all, axis=0)
B2_std = np.std(B2_all, axis=0)
C2_std = np.std(C2_all, axis=0)

# 若需要平均值的标准误（误差棒常用），可计算如下
A1_se = A1_std / np.sqrt(n_runs)
B1_se = B1_std / np.sqrt(n_runs)
C1_se = C1_std / np.sqrt(n_runs)
A2_se = A2_std / np.sqrt(n_runs)
B2_se = B2_std / np.sqrt(n_runs)
C2_se = C2_std / np.sqrt(n_runs)

关键说明 #

• 确保evolve_system返回的序列长度为nsteps + 1（包含初始状态），若实际长度不同，需调整存储数组的第二维度
• 使用axis=0指定对运行次数维度进行统计计算，得到每一步的平均和误差值
• 标准差反映单次运行结果的离散程度，标准误反映平均值的可靠性，可根据需求选择使用

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