0
$\begingroup$

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?

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

1 Answer 1

2
$\begingroup$

Gelman-Rubin statistics are a common heuristic for determining the convergence of sampling schemes such as Markov Chain Monte Carlo (MCMC) methods. The key idea is to simulate multiple independent sampling sets, termed chains, then compare the variance within a chain and the variance between chains. If there is little variance between chains, then we consider the chains to have converged.

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.