Properly scale random movement in simulation with variable time delta

Viewed 153

I want to write a simulation of a particle moving in 1D. The simulation has the following parameters:

  • max_speed: the max speed the particle can obtain (in meters/second)
  • delta: the time interval between two simulation steps (in seconds)
  • iters: the number of steps to run the simulation

Below is an example of a working simulation. It returns two lists of values: one for the particle's position at each time step, and one for the corresponding seconds at each time step.

import random

def simulate_particle(max_speed, delta, iters):
    current_position = 0
    positions = [current_position]
    seconds = [0]
    for frame in range(1, iters+1):
        max_distance_possible = max_speed * delta
        current_position += max_distance_possible * random.uniform(-1, 1)
        positions.append(current_position)
        seconds.append(frame * delta)
    return positions, seconds

As you can see, the particle's position at the next frame is determined by:

  • calculating the max distance possible that the particle is able to travel between two frames (taking into account the max_speed and delta parameters)
  • multiplying this value by a value sampled from an uniform distribution over [-1, 1] (in order to add some randomness)

The issue here is that the particle movement becomes tied to the delta parameter. Since larger deltas allow for a larger max distance, it becomes easier for particles in simulations with large deltas to move further away from the starting position.

Here's a graphical representation of this issue:

import matplotlib.pyplot as plt

for delta, iters, color in zip(
    [600, 60, 6],  # decreasing deltas
    [10, 100, 1000],  # larger number of iterations to achieve same end time
    ['blue', 'orange', 'green'],
):
    for repeat in range(50):  # run simulations 50 times for each delta
        ys, xs = simulate_particle(max_speed=5, delta=delta, iters=iters)
        label = f'{delta=}, {iters=}' if repeat == 0 else None
        plt.plot(xs, ys, color=color, label=label)

plt.xlabel('Time (s)')
plt.ylabel('Position')
plt.legend()
plt.show()

output: enter image description here

I understand why this happens: It's much harder to reach extreme values when sampling many times from a narrow distribution, compared to sampling a few times from a broader one.

My question is how do I fix this? Is there some kind of normalization I can do? I tried changing the uniform distribution to other distributions (e.g. gaussian), but eventually the same effect happens.

2 Answers

As you already know you are simulating different processes, which is why you obtain different behaviours.

A way of "rescaling" your curves is to plot the results in reduced units: divide the time (x) array by delta (which sets the timescale) and the position (y) array by delta * max_speed (which sets the lengthscale) to obtain something like this (note that I zoomed in to highlight the region of the plot where the overlap between the three sets of traces is a bit more visible):

new plot

The main thing that you are running into is that composing multiple random runs headlong into the Central Limit Theorem. Namely, as you have seen, when you add up multiple uniform draws, you approach a normal distribution in the final particle locations rather than a single larger uniform distribution. Drawing from the uniform distribution in order to skip timesteps means that you are going to be effectively sampling from a distribution with much larger shoulders than you should. Hence the blue lines blowing up.

To fix this, you should figure out what your true fundamental timestep delta0 is and run with that. If you want to speed up by some factor, then you should draw from the Irwin–Hall distribution, and scale & shift that draw by the same transformation you would make to a [0, 1] uniform distribution to get it in the [-1, 1] range. And whether you speed up or not, you should change your max_distance_possible equation to rely on the fundamental timestep.

So for example, here's a bit of pseudocode to replace the main loop:

delta0 = 1
speedupfactor = 10
delta = speedupfactor * delta0

# IrwinHall scaled and shifted as if the underlying uniform draw was over [-1,1]
randdraw = 2*IrwinHall(n=speedupfactor) - speedupfactor  
max_distance_possible = max_speed * delta0
current_position += max_distance_possible * randdraw

positions.append(current_position)
seconds.append(frame * delta)

If your process is truly continuous, then the limiting case of delta0 = 0 means this will be equivalent to drawing from a purely normal distribution, and you can do that instead. This means you will have to define movement in terms of a standard deviation rather than a max speed.

Unfortunately scipy.stats does not have Irwin-Hall as a probability distribution, but it is implemented in sympy as sympy.stats.UniformSum()

Related