<jupyter_start><jupyter_text>Autocorrelation Diagnostic (Alpha) Autocorrelation plots show an estimate of the independence of a time series signal with itself.  As a diagnostic for Markov Monte Carlo Chains they provides an estimate the independence of each draw from previous draws in the same chain. High auto-correlation does not necessarily imply lack of convergence, only loss of efficiency. The higher the autocorrelation the larger the sample size you will need to get the correct answer, but you can have convergence with high autocorrelation.In fact Markov chain theory guarantees that under certain (mild) conditions you will get to the correct stationary distribution. Unfortunately, this guarantee requires infinite samples. So the trick is to guess when the chain has "practically reached infinity" for a given problem. Autocorrelation plots help the practioner make that assessment.  Arviz Autocorrelation plotsBy default ArviZ will generate an autocorrelation plot for each variable and for each chain. The autocorrelation should decrease as the lag increases, which indicates a low level of correlation.<jupyter_code>import arviz as az
import numpy as np
import matplotlib.pyplot as plt

az.style.use('arviz-darkgrid')
data = az.load_arviz_data('centered_eight')
ax = az.plot_autocorr(data, var_names="mu")<jupyter_output><empty_output><jupyter_text>Autocorrelation deriviationAutocorrelation at each lag is defined as $$\large \rho_{lag} = \frac{\Sigma_{i}^{N-lag} (\theta_{i} - \overline{\theta})(\theta_{i + lag} - \overline{\theta})}{\Sigma_{i}^{N}(\theta_{i} - \overline{\theta})}$$ This formula may be recognized as a special form of Pearson's correlation coefficient, where instead of measuring the correlation between two arbitrary sets, we are instead measuring the correlation of different time slices of the same chain. We can manually calculate the autocorrelation using numpy<jupyter_code>chain_0 = data.posterior.sel(chain=0)
mu_0 = chain_0["mu"]
lag = np.arange(0,100)
auto_corr = []

for shift in lag:
    # Or None is to handle the case of indexing at -0
    corr = np.corrcoef(mu_0[:-shift or None],mu_0[shift:])[0,1]
    auto_corr.append(corr)
    
auto_corr[:5]
fig, ax = plt.subplots(figsize=(12,6))
ax.vlines(x=np.arange(1, 100),
                  ymin=0, ymax=auto_corr,
                  lw=1)
fig.suptitle("Autocorrelation of Mu Chain 0");
# Filter to just one chain
axes = az.plot_autocorr(chain_0, var_names="mu", figsize=(12,6))
axes[0][0].set_title("Autocorrelation of Mu Chain 0");<jupyter_output><empty_output><jupyter_text>Note that this plot is identical to chain 0 for `mu` in the posterior data  Additional examples MCMC chains with high autocorrelation can indicate that a lack of convergence to a stationary distribution. Take the following example where we try and use PyMC3 estimate the parameters of a generated distribution.<jupyter_code>import pymc3 as pm

# Generate observed distribution with fixed parameters
SD = 2
MU = -5
obs = np.random.normal(loc=MU, scale=SD, size=10000)
# Attempt to use pymc3 to estimate mean of the distribution
with pm.Model() as model:
    mu = pm.Normal("mu", mu=-5000, sd=1)
    y = pm.Normal("y", mu=mu, sd=SD, observed=obs)
    step = pm.Metropolis()
    trace = pm.sample(100, step, chains=2)<jupyter_output>Only 100 samples in chain.
Multiprocess sampling (2 chains in 2 jobs)
Metropolis: [mu]
Sampling 2 chains: 100%|██████████| 1200/1200 [00:00<00:00, 2729.68draws/s]
The gelman-rubin statistic is larger than 1.2 for some parameters.
The number of effective samples is smaller than 10% for some parameters.<jupyter_text>In this example`Mu`is the unknown parameter that is being estimated by PyMC3. The model is specified with a prior Normal prior with a mean of -5000, and only 100 steps are taken with a Metropolis Hasting sampler.<jupyter_code>axes = az.plot_autocorr(trace, max_lag=100)<jupyter_output><empty_output><jupyter_text>When plotting the autocorrelation of the chains however a high degree of correlation is present, indicating that the sampling did not converge to a stationary distribution.<jupyter_code># Attempt to use pymc3 to estimate mean of the distribution
with pm.Model() as model:
    mu = pm.Uniform("mu")
    y = pm.Normal("y", mu=mu, sd=SD, observed=obs)
    step = pm.Metropolis()
    trace = pm.sample(5000, step, chains=2)
    
axes = az.plot_autocorr(trace)<jupyter_output>Multiprocess sampling (2 chains in 2 jobs)
Metropolis: [mu]
Sampling 2 chains: 100%|██████████| 11000/11000 [00:02<00:00, 5021.88draws/s]
The number of effective samples is smaller than 25% for some parameters.
{'fraction_non_alphanumeric': 0.08935155063977446, 'fraction_numerical': 0.029277813923227064, 'mean_word_length': 5.076416337285902, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 19, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '196962', 'license': 'MIT', 'n_tokens_mistral': 1487, 'n_tokens_neox': 1342, 'n_words': 548, 'path': 'archive/plots/AutoCorrelation.ipynb', 'repo_name': 'rpgoldman/Exploratory-Analysis-of-Bayesian-Models', 'size': '4611'}