Stable Distributions for Animations

In many time-controlled creative coding projects like animations, I need pseudo-random parameters to move around smoothly without straying too far from the origin. For a simple approach, I assign each value a position $x_t$ velocity $v_t$ . For a discrete time variable $t$ , the position updates as usual according to the Semi-implicit Euler method $x_{t+1}=x_t+v_{t+1}$ , this brings numerical stability and is more straightforward to implemen. For the velocity, I choose

v_{t+1}=\gamma v_t−kx_t+\epsilon_t.

where $0<\gamma<1$ is a friction term, $0<k$ is a spring constant which keeps parameters close to $0$ , and $\epsilon$ is noise from a Gaussian distribution with variance $\sigma^2$ .

Often we want to start the animation with a random constellation drawn from a distribution that remains stable over time, so the parameters on average shouldn’t come closer, spread further or change speed during the process.

To describe the spread of the system, we use the variances of $x$ and $v$ together with their covariance. Remember, this is what our system looks like:

\begin{align*} x_{t+1}&=x_t+v_{t+1}=(1-k)x_t+\gamma v_t+\epsilon_t\\ v_{t+1}&=\gamma v_t-k x_t+\epsilon_t\\ \end{align*}

By following the simple computation rule for covariance of linear combinations, we arrive at the following covariance matrix:

\begin{align*} \var[x_{t+1}]&=(1-k)^2\var[x_t]+\gamma^2\var[v_t]+\gamma(2-2k)\cov[x_t,v_t]+\sigma^2\\ \var[v_{t+1}]&=k^2\var[x_t]+\gamma^2\var[v_t]-2\gamma k\cov[x_t,v_t]+\sigma^2\\ \cov[x_{t+1},v_{t+1}]&=(k^2-k)\var[x_t]+\gamma^2\var[v_t]+\gamma(1-2k)\cov[x_t,v_t]+\sigma^2 \end{align*}

Assuming that variance and covariance do not depend $t$ (stable distribution), we can ask WolframAlpha to solve this system for $V:=\var[x]$ , $W:=\var[v]$ and $C:=\cov[x,v]$ .

Writing down the solution and simplify a bit, we get a commond denominator $D:=(1-\gamma)(2(1+\gamma)-k)$ and then:

V=\frac{(\gamma+1)\sigma^2}{kD},\quad W=\frac{2\sigma^2}{D},\quad C=\frac{\sigma^2}{D}

This means that if we start the system in a distribution satisfying this covariance matrix, the matrix won’t change during the run. However, a distribution is more than just its covariance matrix, so the distribution itself could still change.

However, we don’t need to worry when choosing a bivariate normal distribution, which you can think of (assuming all means are zero) as two variables constructed as linear combinations from two independent standard normal distributions. This family of distributions has the useful property of being closed under linear combinations and also fully specified by its covariance matrix. This means, if we start with a bivariate distribution satisfying our stable matrix, we get indeed a truly stable distribution.

To implement this distribution, it’s better to work with the standard deviations $\sigma_x:=\sqrt{\var[x]}$ , $\sigma_v:=\sqrt{\var[v]}$ and correlation $\rho_{xv}:=\cov[x,v]/(\sigma_x\sigma_v)$ . The correlation always lives in $[-1,1]$ , and vice versa, every number in this interval is a valid correlation, as long as the standard deviations are positive, while the covariance matrix must be positive definite, which is more difficult to check.

We can express a bivariate distribution as linear combination of two independent standard normal distributions $z_0$ and $z_1$ as follows (this is not the unique decomposition, but a canonical form):

x=\sigma_x z_0,\quad v=\sigma_v(\rho_{xv} z_0+\sqrt{1-\rho_{xv}^2} z_1)

Before implementing this stable distribution, one more consideration is that in some environments, like the web browser’s requestAnimationFrame() method, time intervals between frames can vary. To handle this, let $\Delta$ denote the elapsed time between frame $t$ and $t+1$ . In this case, we may update position according to $x_{t+1}=x_t+\Delta v_{t+1}$ , and velocity according to $v_{t+1}=\gamma^\Delta v_t−\Delta kx_t+\sqrt{\Delta}\epsilon$ . This is a good heuristic choice, if we set two of the three parameters to zero, we can check that different subdivisions of the same time interval lead to the exact same outcome (in the case of $\epsilon$ , same outcome in a stochastic sense).

Here is a simple implementation in JavaScript, click to reset the simulation:

function boxMuller() {
  const phi = 2 * Math.PI * Math.random();
  const r = Math.sqrt(-2 * Math.log(Math.random()));
  return [r * Math.cos(phi), r * Math.sin(phi)];
}

function bivariate({ sigma0, sigma1, rho }) {
  const [z0, z1] = boxMuller();
  return [sigma0 * z0, sigma1 * (rho * z0 + Math.sqrt(1 - rho * rho) * z1)];
}

function stableParams({ gamma, k, sigma }) {
  const S = sigma * sigma;
  const D = (1 - gamma) * (2 * (1 + gamma) - k);
  const varX = ((1 + gamma) * S) / (k * D);
  const varV = (2 * S) / D;
  const covXV = S / D;
  const sigma0 = Math.sqrt(varX);
  const sigma1 = Math.sqrt(varV);
  return { sigma0, sigma1, rho: covXV / (sigma0 * sigma1) };
}

function resetParticles(particles, params) {
  for (let p of particles) {
    const [x0, v0] = bivariate(params);
    const [x1, v1] = bivariate(params);
    p.x = [x0, x1];
    p.v = [v0, v1];
  }
}

function updateParticles(particles, vars, dt) {
  const gammaDt = Math.pow(vars.gamma, dt);
  const kDt = dt * vars.k;
  const sigmaDt = Math.sqrt(dt) * vars.sigma;

  for (let p of particles) {
    const [z0, z1] = boxMuller();
    p.v[0] = gammaDt * p.v[0] - kDt * p.x[0] + sigmaDt * z0;
    p.v[1] = gammaDt * p.v[1] - kDt * p.x[1] + sigmaDt * z1;
    p.x[0] += dt * p.v[0];
    p.x[1] += dt * p.v[1];
    p.el.style.transform = `translate(${[
      `calc(${p.x[0]} * var(--dim))`,
      `calc(${p.x[1]} * var(--dim))`,
    ].join(",")}`;
  }
}