0
# Convergence Analysis
1
2
emcee provides statistical tools for assessing MCMC chain convergence through autocorrelation analysis. These functions help determine when chains have converged and how many effective samples have been obtained, which is crucial for reliable Bayesian inference.
3
4
## Capabilities
5
6
### Autocorrelation Time Estimation
7
8
Estimate the integrated autocorrelation time, which measures how correlated successive samples are in the chain.
9
10
```python { .api }
11
def integrated_time(x, c: int = 5, tol: int = 50, quiet: bool = False,
12
has_walkers: bool = True):
13
"""
14
Estimate integrated autocorrelation time of time series.
15
16
Args:
17
x: Input time series [steps, ...] or [steps, nwalkers, ndim]
18
c: Window factor for automatic windowing (default: 5)
19
tol: Minimum chain length as multiple of autocorr time (default: 50)
20
quiet: Suppress convergence warnings (default: False)
21
has_walkers: Whether input has walker dimension (default: True)
22
23
Returns:
24
ndarray: Autocorrelation time for each parameter
25
26
Raises:
27
AutocorrError: If autocorrelation time cannot be reliably estimated
28
"""
29
```
30
31
### Autocorrelation Function
32
33
Compute the normalized autocorrelation function for 1D time series.
34
35
```python { .api }
36
def function_1d(x):
37
"""
38
Estimate normalized autocorrelation function of 1D series.
39
40
Args:
41
x: 1D time series as numpy array
42
43
Returns:
44
ndarray: Autocorrelation function of the time series
45
46
Raises:
47
ValueError: If input is not 1D
48
"""
49
```
50
51
### Autocorrelation Exceptions
52
53
Exception raised when autocorrelation analysis fails or is unreliable.
54
55
```python { .api }
56
class AutocorrError(Exception):
57
"""
58
Raised when autocorrelation time estimation fails.
59
60
This typically occurs when:
61
- Chain is too short relative to autocorrelation time
62
- Chain has not converged
63
- Numerical issues in autocorrelation computation
64
"""
65
```
66
67
## Usage Examples
68
69
### Basic Autocorrelation Analysis
70
71
```python
72
import emcee
73
import numpy as np
74
75
def log_prob(theta):
76
return -0.5 * np.sum(theta**2)
77
78
# Run MCMC
79
sampler = emcee.EnsembleSampler(32, 2, log_prob)
80
pos = np.random.randn(32, 2)
81
sampler.run_mcmc(pos, 1000)
82
83
# Compute autocorrelation time
84
try:
85
tau = sampler.get_autocorr_time()
86
print(f"Autocorrelation time: {tau}")
87
print(f"Mean tau: {np.mean(tau):.1f}")
88
except emcee.autocorr.AutocorrError as e:
89
print(f"Autocorr analysis failed: {e}")
90
```
91
92
### Chain Length Assessment
93
94
```python
95
# Check if chain is long enough
96
chain = sampler.get_chain()
97
nsteps = chain.shape[0]
98
99
try:
100
tau = sampler.get_autocorr_time(tol=20) # Less strict tolerance
101
102
# Rule of thumb: chain should be >> 50 * tau
103
min_length = 50 * np.max(tau)
104
if nsteps > min_length:
105
print(f"Chain length ({nsteps}) is sufficient (need > {min_length:.0f})")
106
else:
107
print(f"Chain too short! Need > {min_length:.0f} steps, have {nsteps}")
108
109
except emcee.autocorr.AutocorrError:
110
print("Chain too short for reliable autocorr estimation")
111
```
112
113
### Progressive Autocorrelation Monitoring
114
115
```python
116
# Monitor autocorrelation during sampling
117
def monitor_autocorr(sampler, nsteps=100, check_interval=50):
118
pos = np.random.randn(sampler.nwalkers, sampler.ndim)
119
120
for i, state in enumerate(sampler.sample(pos, iterations=nsteps)):
121
if i % check_interval == 0 and i > 0:
122
try:
123
tau = sampler.get_autocorr_time(tol=10, quiet=True)
124
print(f"Step {i}: tau = {np.mean(tau):.1f}")
125
126
# Check convergence criterion
127
if i > 50 * np.max(tau):
128
print(f"Converged at step {i}")
129
break
130
131
except emcee.autocorr.AutocorrError:
132
print(f"Step {i}: tau estimation failed")
133
134
sampler = emcee.EnsembleSampler(32, 2, log_prob)
135
monitor_autocorr(sampler, nsteps=2000)
136
```
137
138
### Effective Sample Size
139
140
```python
141
# Calculate effective number of independent samples
142
chain = sampler.get_chain(flat=True)
143
nsamples = chain.shape[0]
144
145
try:
146
tau = sampler.get_autocorr_time()
147
148
# Effective sample size
149
n_eff = nsamples / (2 * tau)
150
print(f"Total samples: {nsamples}")
151
print(f"Effective samples per parameter: {n_eff}")
152
print(f"Minimum effective samples: {np.min(n_eff):.0f}")
153
154
except emcee.autocorr.AutocorrError:
155
print("Cannot compute effective sample size - autocorr failed")
156
```
157
158
### 1D Autocorrelation Function
159
160
```python
161
from emcee.autocorr import function_1d
162
163
# Get single parameter chain
164
chain = sampler.get_chain(flat=True)
165
param_chain = chain[:, 0] # First parameter
166
167
# Compute autocorrelation function
168
acf = function_1d(param_chain)
169
170
# Plot autocorrelation function
171
import matplotlib.pyplot as plt
172
plt.figure()
173
plt.plot(acf[:100]) # First 100 lags
174
plt.axhline(y=0, color='k', linestyle='--', alpha=0.3)
175
plt.xlabel('Lag')
176
plt.ylabel('Autocorrelation')
177
plt.title('Autocorrelation Function')
178
plt.show()
179
180
# Find where ACF drops to 1/e
181
import numpy as np
182
cutoff = 1/np.e
183
lag_1e = np.argmax(acf < cutoff)
184
print(f"Autocorr drops to 1/e at lag: {lag_1e}")
185
```
186
187
### Burn-in Determination
188
189
```python
190
# Use autocorrelation to determine burn-in
191
def estimate_burnin(sampler, min_length=100):
192
chain = sampler.get_chain()
193
nsteps = chain.shape[0]
194
195
# Try different burn-in lengths
196
burnin_candidates = np.arange(min_length, nsteps//2, 50)
197
tau_estimates = []
198
199
for burnin in burnin_candidates:
200
try:
201
# Compute tau for post-burn-in chain
202
post_burnin_chain = chain[burnin:]
203
tau = emcee.autocorr.integrated_time(post_burnin_chain, quiet=True)
204
tau_estimates.append(np.mean(tau))
205
except emcee.autocorr.AutocorrError:
206
tau_estimates.append(np.inf)
207
208
# Find where tau stabilizes
209
tau_estimates = np.array(tau_estimates)
210
stable_idx = np.argmin(np.abs(tau_estimates - np.median(tau_estimates[:-10])))
211
212
recommended_burnin = burnin_candidates[stable_idx]
213
print(f"Recommended burn-in: {recommended_burnin} steps")
214
215
return recommended_burnin
216
217
burnin = estimate_burnin(sampler)
218
```
219
220
### Thinning Based on Autocorrelation
221
222
```python
223
# Thin chain to get approximately independent samples
224
try:
225
tau = sampler.get_autocorr_time()
226
thin_factor = int(2 * np.max(tau)) # Thin by ~2x autocorr time
227
228
print(f"Autocorrelation time: {np.max(tau):.1f}")
229
print(f"Recommended thinning factor: {thin_factor}")
230
231
# Get thinned chain
232
thinned_chain = sampler.get_chain(discard=burnin, thin=thin_factor, flat=True)
233
print(f"Thinned chain shape: {thinned_chain.shape}")
234
235
# Effective sample size after thinning
236
n_eff_thinned = thinned_chain.shape[0]
237
print(f"Effective independent samples: {n_eff_thinned}")
238
239
except emcee.autocorr.AutocorrError:
240
print("Using default thinning factor of 10")
241
thinned_chain = sampler.get_chain(discard=burnin, thin=10, flat=True)
242
```
243
244
### Convergence Diagnostics
245
246
```python
247
def convergence_diagnostics(sampler, target_n_eff=1000):
248
"""Comprehensive convergence assessment."""
249
250
chain = sampler.get_chain()
251
nsteps, nwalkers, ndim = chain.shape
252
253
print(f"Chain shape: {nsteps} steps × {nwalkers} walkers × {ndim} dims")
254
print(f"Total samples: {nsteps * nwalkers}")
255
256
try:
257
# Autocorrelation analysis
258
tau = sampler.get_autocorr_time(tol=20)
259
max_tau = np.max(tau)
260
261
print(f"\nAutocorrelation times: {tau}")
262
print(f"Maximum tau: {max_tau:.1f}")
263
264
# Chain length assessment
265
min_required = 50 * max_tau
266
if nsteps > min_required:
267
print(f"✓ Chain length sufficient ({nsteps} > {min_required:.0f})")
268
else:
269
print(f"✗ Chain too short (need {min_required:.0f}, have {nsteps})")
270
271
# Effective sample size
272
n_eff = nsteps * nwalkers / (2 * tau)
273
min_n_eff = np.min(n_eff)
274
275
print(f"Effective samples per param: {n_eff}")
276
print(f"Minimum effective samples: {min_n_eff:.0f}")
277
278
if min_n_eff > target_n_eff:
279
print(f"✓ Sufficient effective samples ({min_n_eff:.0f} > {target_n_eff})")
280
else:
281
print(f"✗ Need more samples (want {target_n_eff}, have {min_n_eff:.0f})")
282
283
# Acceptance fraction
284
acc_frac = np.mean(sampler.acceptance_fraction)
285
print(f"\nAcceptance fraction: {acc_frac:.3f}")
286
287
if 0.2 <= acc_frac <= 0.5:
288
print("✓ Acceptance fraction in good range")
289
else:
290
print("⚠ Acceptance fraction outside optimal range (0.2-0.5)")
291
292
except emcee.autocorr.AutocorrError as e:
293
print(f"✗ Autocorrelation analysis failed: {e}")
294
print("Chain likely too short or not converged")
295
296
# Run diagnostics
297
convergence_diagnostics(sampler)
298
```