# How can I determine if a system is equilibrated?

Cross-posted in CrossValidated.SE and MMSE

I am experimenting with a new MCMC protocol and new research.

In the context of Monte Carlo simulation, a "state of equilibrium" refers to a condition where the system being simulated has reached a stable configuration or distribution of properties that does not significantly change over time. This state represents a balance or equilibrium between the various forces, interactions, and constraints acting on the system.

I wrote the following function to test if a system is equilibrated:

import numpy as np

def is_equilibrated(data, shave_size=50000):
returns = False
mid_index = len(data) // 2

first_half = data[:mid_index]
second_half = data[mid_index:]

decimals = 6

while True:
first_auto_corr = np.correlate(first_half, first_half, mode='full')
second_auto_corr = np.correlate(second_half, second_half, mode='full')

corr1 = np.corrcoef(first_auto_corr, second_auto_corr)[0, 1]

print(f"Size {len(first_half)} vs. {len(second_half)}")
print(f"corr  = {corr1} ")

first_half = first_half[shave_size:]
second_half = second_half[shave_size:]

if corr1 == 1.0:
break

return returns


The function iteratively compares the autocorrelation of two halves of the input data, shortening the dataset at each iteration, until perfect correlation (equilibrium) is achieved between the two halves or until the loop is manually terminated.

I wrote this function to check if my simulation data is "matured" enough to analyze, i.e., good enough to computationally extract further information. In other words, it lets me know if I ran the simulation long enough to obtain a useful dataset.

The drawback of this algorithm is that it uses autocorrelation computation, which takes a lot of time, irrespective of the programming language used.

Can you propose any alternative technique?

• I think it would be useful if you outlined in formulas what it is you are trying to do. Commented Mar 5 at 14:35
• Are you familiar with Gelman-Rubin statistics? Commented Mar 5 at 16:35
• @whpowell96, Are you familiar with Gelman-Rubin statistics? --- No. Commented Mar 5 at 18:10

Let $$x_i^j$$ be the $$i$$th sample from the $$j$$th chain, with $$1\leq i \leq M$$ and $$1\leq j \leq N$$, so we have $$N$$ chains with $$M$$ samples per chain. We then define the following statistics: \begin{aligned} \bar{x}^j &= \frac{1}{M} \sum_{i=1}^M x_i^j\\ \bar{x} &= \frac{1}{N}\sum_{j=1}^N \bar{x}^j \\ B &= \frac{M}{N-1}\sum_{j=1}^N (\bar{x}^j - \bar{x})^2 \\ W &= \frac{1}{N}\sum_{j=1}^M\left(\frac{1}{M-1}\sum_{i=1}^N (x_i^j - \bar{x}^j)\right). \end{aligned} These statistics represent the mean of chain $$j$$, the mean of means of all chains, the variance between chains, and the total variance of all chains, respectively. The Gelma-Rubin statistic is then given by $$R = \frac{\frac{M-1}{M}W + \frac{1}{L}B}{W}.$$ This quantity, or sometimes its square root, termed the potential scale reduction factor, can be used to determine convergence of multiple MCMC runs by noticing that when each chain has converged to the same stationary distribution, $$B=0$$ and $$R\to1$$ as $$M\to\infty$$. It is common to compute this statistic for multiple chains after burn-in periods have been discarded and sometimes you split each chain in half and consider them as $$2N$$ chains. Usually $$R<1.1$$ is considered good. There are problems with Gelman-Rubin statistics and modifications, improvements, and alternative posterior analysis tools are used in tandem to determine convergence. For a good starting place, see this section of the Stan reference manual.