# Numerical stability of taking the mean of outputs from the simulation of a discrete stochastic dynamical system

I am writing a simulation for a discrete stochastic dynamical system. Since the simulation is stochastic, I need to run the simulation multiple times and then average the values of each timestep. I have some simple julia code below to demonstrate the process.

My question was really about numerical stability when averaging data from say 100, 1000, 10000 simulations. The simulation is discrete and not continuous, so I don't need to worry as much about the fine grained discretization, where error accumulation really matters. But, I was not sure if there are issues with round-off error accumulating over time even for something as simple as a mean?

In the code below, I keep a running mean of some simulation output v. So for each run through the simulation, I just take the average of the running mean and the new value--in this case a random number rand(). Is this a numerically stable way to do this, or do I need to allocate and save each simulation to some large array or file, and then average them at the end?

n2 = 5
v = 0.0
for i in 1:n2
v = mean([v, rand()])
end


Any suggestions would really be appreciated. I am trying to find the right balance between numerical stability versus limiting unnecessary allocations of memory. Thanks.

• The formula for your running mean is not correct. With each new sample, you are averaging the previous average and the new value, but you need to weigh the previous average with the number of previously seen samples from which it was computed. The correct formula is $\bar x_n = \frac{n-1}{n} \bar x_{n-1} + \frac{1}{n} x_n$. Dec 14, 2020 at 18:12
• @WolfgangBangerth Oh thanks so much for correcting my formula. es, I have seen this formula mentioned before, but I did not connect it with my own code. Thanks for setting me straight. Dec 14, 2020 at 19:37