0
# MCMC Sampling Methods
1
2
PyMC3 provides state-of-the-art Markov Chain Monte Carlo (MCMC) algorithms for sampling from posterior distributions. The library includes advanced samplers like the No-U-Turn Sampler (NUTS), various Metropolis methods, and specialized algorithms, with automatic step method assignment and adaptive tuning.
3
4
## Capabilities
5
6
### Main Sampling Functions
7
8
Core functions for MCMC sampling and posterior analysis.
9
10
```python { .api }
11
def sample(draws=1000, step=None, init='auto', n_init=200000,
12
initvals=None, trace=None, chain_idx=0, chains=None, cores=None,
13
tune=1000, progressbar=True, model=None, random_seed=None,
14
discard_tuned_samples=True, compute_convergence_checks=True,
15
callback=None, jitter_max_retries=10, start=None,
16
return_inferencedata=None, idata_kwargs=None, mp_ctx=None,
17
pickle_backend='pickle', **kwargs):
18
"""
19
Draw samples from the posterior distribution using MCMC.
20
21
The main interface for MCMC sampling in PyMC3. Automatically selects
22
appropriate step methods, performs adaptive tuning, and returns results
23
in ArviZ InferenceData format for analysis.
24
25
Parameters:
26
- draws: int, number of samples to draw from posterior
27
- step: step method or list of step methods (auto-assigned if None)
28
- init: str or dict, initialization method ('auto', 'adapt_diag', 'jitter+adapt_diag', 'advi', 'advi_map', 'map')
29
- n_init: int, number of initialization attempts (default 200000)
30
- initvals: PointType or Sequence, initial values for variables
31
- trace: Backend, trace backend for storing samples
32
- chain_idx: int, chain index for multiprocessing
33
- chains: int, number of parallel chains (defaults to number of cores)
34
- cores: int, number of CPU cores to use (all available if None)
35
- tune: int, number of tuning steps before sampling
36
- progressbar: bool, display progress bar
37
- model: Model, model to sample from (current context if None)
38
- random_seed: int or list, random seeds for chains
39
- discard_tuned_samples: bool, exclude tuning samples from results
40
- compute_convergence_checks: bool, compute R-hat and effective sample size
41
- callback: function, callback function called after each sample
42
- jitter_max_retries: int, maximum retries for jittered initialization
43
- start: dict, starting values for variables (deprecated, use initvals)
44
- return_inferencedata: bool, return ArviZ InferenceData object
45
- idata_kwargs: dict, arguments for InferenceData construction
46
- mp_ctx: multiprocessing context
47
- pickle_backend: str, backend for pickling ('pickle' or 'dill')
48
49
Returns:
50
- InferenceData or MultiTrace: posterior samples and diagnostics
51
"""
52
53
def sample_posterior_predictive(trace, samples=None, model=None,
54
var_names=None, size=None,
55
keep_size=False, random_seed=None,
56
progressbar=True, return_inferencedata=True,
57
extend_inferencedata=True,
58
predictions=False, **kwargs):
59
"""
60
Generate posterior predictive samples from model.
61
62
Uses posterior parameter samples to generate predictions from the
63
likelihood, enabling model checking and out-of-sample prediction.
64
65
Parameters:
66
- trace: InferenceData or MultiTrace, posterior samples
67
- samples: int, number of posterior predictive samples (all posterior samples if None)
68
- model: Model, model for prediction (current context if None)
69
- var_names: list, variables to sample (all observed RVs if None)
70
- size: int, size of each posterior predictive sample
71
- keep_size: bool, preserve original variable sizes
72
- random_seed: int, random seed for reproducibility
73
- progressbar: bool, display progress bar
74
- return_inferencedata: bool, return as InferenceData
75
- extend_inferencedata: bool, add to existing InferenceData
76
- predictions: bool, sample from out-of-sample predictions
77
78
Returns:
79
- dict or InferenceData: posterior predictive samples
80
"""
81
82
def sample_posterior_predictive_w(traces, samples=None, models=None,
83
weights=None, **kwargs):
84
"""
85
Generate weighted posterior predictive samples from multiple models.
86
87
Combines posterior predictive samples from multiple models using
88
specified weights, useful for Bayesian model averaging.
89
90
Parameters:
91
- traces: list, posterior traces from multiple models
92
- samples: int, number of samples per model
93
- models: list, models corresponding to traces
94
- weights: array, model weights (uniform if None)
95
- kwargs: additional arguments for sample_posterior_predictive
96
97
Returns:
98
- dict: weighted posterior predictive samples
99
"""
100
101
def sample_prior_predictive(samples=500, model=None, var_names=None,
102
random_seed=None):
103
"""
104
Generate samples from the prior predictive distribution.
105
106
Samples from the prior distributions of variables to simulate data
107
before observing any evidence, useful for prior predictive checks
108
and model validation.
109
110
Parameters:
111
- samples: int, number of samples from the prior predictive (default 500)
112
- model: Model, model to sample from (current context if None)
113
- var_names: list, variables to sample (all observed and unobserved RVs if None)
114
- random_seed: int, random seed for reproducibility
115
116
Returns:
117
- dict: prior predictive samples for each variable
118
"""
119
120
def iter_sample(draws, step, start=None, trace=None, chain=0,
121
tune=None, model=None, random_seed=None, callback=None):
122
"""
123
Generator for iterative MCMC sampling.
124
125
Provides fine-grained control over sampling process, yielding
126
individual samples for custom processing or online analysis.
127
128
Parameters:
129
- draws: int, number of samples to draw
130
- step: step method for sampling
131
- start: dict, starting parameter values
132
- trace: BaseTrace, trace object to store samples
133
- chain: int, chain identifier
134
- tune: int, number of tuning steps (0 if None)
135
- model: Model, model to sample from
136
- random_seed: int, random seed
137
- callback: function, called after each sample
138
139
Yields:
140
- MultiTrace: updated trace after each sample
141
"""
142
143
def stop_tuning(step):
144
"""
145
Stop tuning phase for step methods.
146
147
Parameters:
148
- step: step method or CompoundStep
149
"""
150
```
151
152
### Step Method Assignment
153
154
Functions for automatic and manual step method assignment.
155
156
```python { .api }
157
def assign_step_methods(model, step=None, methods=None,
158
step_kwargs=None, **kwargs):
159
"""
160
Assign appropriate step methods to model variables.
161
162
Automatically selects optimal step methods based on variable properties,
163
or uses user-specified methods for custom sampling strategies.
164
165
Parameters:
166
- model: Model, model to assign step methods for
167
- step: step method or list, user-specified step methods
168
- methods: list, available step method classes
169
- step_kwargs: dict, keyword arguments for step methods
170
171
Returns:
172
- list: assigned step methods for each variable
173
"""
174
175
def instantiate_steppers(model, steps, selected, step_kwargs=None):
176
"""
177
Create step method instances for assigned variables.
178
179
Parameters:
180
- model: Model, model context
181
- steps: list, step method classes
182
- selected: list, variables assigned to each step method
183
- step_kwargs: dict, keyword arguments for step methods
184
185
Returns:
186
- list: instantiated step method objects
187
"""
188
```
189
190
### Hamiltonian Monte Carlo Methods
191
192
Gradient-based MCMC methods that use Hamiltonian dynamics for efficient exploration of parameter space.
193
194
```python { .api }
195
class NUTS:
196
"""
197
No-U-Turn Sampler (NUTS) - adaptive Hamiltonian Monte Carlo.
198
199
NUTS automatically tunes step size and path length to efficiently
200
sample from complex posterior distributions. Default sampler for
201
continuous variables in PyMC3.
202
"""
203
204
def __init__(self, vars=None, Emax=1000, target_accept=0.8,
205
gamma=0.05, k=0.75, t0=10, model=None,
206
scaling=None, is_cov=False, potential=None,
207
integrator='leapfrog', **kwargs):
208
"""
209
Initialize NUTS sampler.
210
211
Parameters:
212
- vars: list, variables to sample (all continuous if None)
213
- Emax: float, maximum energy change for U-turn detection
214
- target_accept: float, target acceptance probability for adaptation
215
- gamma: float, adaptation parameter for step size
216
- k: float, adaptation parameter for step size
217
- t0: float, adaptation parameter for step size
218
- model: Model, model context
219
- scaling: array, scaling matrix for mass matrix
220
- is_cov: bool, whether scaling is covariance (precision if False)
221
- potential: Potential, custom potential function
222
- integrator: str, integration method ('leapfrog')
223
"""
224
225
class HamiltonianMC:
226
"""
227
Hamiltonian Monte Carlo sampler with fixed path length.
228
229
Classical HMC implementation with user-specified step size
230
and path length parameters.
231
"""
232
233
def __init__(self, vars=None, path_length=2, step_rand=None,
234
step_scale=0.25, model=None, blocked=True,
235
potential=None, integrator='leapfrog', **kwargs):
236
"""
237
Initialize HMC sampler.
238
239
Parameters:
240
- vars: list, variables to sample
241
- path_length: float, length of Hamiltonian trajectory
242
- step_rand: function, random step size generator
243
- step_scale: float, step size scaling factor
244
- model: Model, model context
245
- blocked: bool, block variables together
246
- potential: Potential, custom potential function
247
- integrator: str, integration method
248
"""
249
```
250
251
### Metropolis Methods
252
253
Random walk and proposal-based MCMC methods for discrete and constrained variables.
254
255
```python { .api }
256
class Metropolis:
257
"""
258
Metropolis-Hastings sampler with customizable proposals.
259
260
General-purpose MCMC method using random walk or custom
261
proposal distributions for parameter updates.
262
"""
263
264
def __init__(self, vars=None, S=None, proposal_dist=None,
265
lambda_=None, model=None, blocked=True, **kwargs):
266
"""
267
Initialize Metropolis sampler.
268
269
Parameters:
270
- vars: list, variables to sample
271
- S: array, scaling matrix for proposals
272
- proposal_dist: distribution, proposal distribution class
273
- lambda_: float, scaling parameter for proposals
274
- model: Model, model context
275
- blocked: bool, update variables jointly
276
"""
277
278
class BinaryMetropolis:
279
"""
280
Metropolis sampler for binary variables.
281
282
Specialized for Bernoulli and binary discrete variables
283
using bit-flip proposals.
284
"""
285
286
class BinaryGibbsMetropolis:
287
"""
288
Gibbs-Metropolis sampler for binary variables.
289
290
Combines Gibbs sampling with Metropolis updates for
291
efficient sampling of binary variable vectors.
292
"""
293
294
class CategoricalGibbsMetropolis:
295
"""
296
Gibbs-Metropolis sampler for categorical variables.
297
298
Efficient sampling for discrete categorical variables
299
using category-wise Gibbs updates.
300
"""
301
302
class DEMetropolis:
303
"""
304
Differential Evolution Metropolis sampler.
305
306
Population-based MCMC method that uses differential
307
evolution for proposal generation.
308
"""
309
310
def __init__(self, vars=None, lamb=None, gamma=1, epsilon=1e-6,
311
save_history=False, model=None, **kwargs):
312
"""
313
Initialize DE Metropolis sampler.
314
315
Parameters:
316
- vars: list, variables to sample
317
- lamb: float, lambda parameter (auto-tuned if None)
318
- gamma: float, differential weight parameter
319
- epsilon: float, noise parameter
320
- save_history: bool, save population history
321
- model: Model, model context
322
"""
323
324
class DEMetropolisZ:
325
"""
326
Differential Evolution Metropolis with Z adaptation.
327
328
Enhanced DE-MCMC with adaptive Z parameter tuning
329
for improved mixing and convergence.
330
"""
331
```
332
333
### Metropolis Proposal Distributions
334
335
Proposal distributions for Metropolis-Hastings updates.
336
337
```python { .api }
338
class NormalProposal:
339
"""
340
Multivariate normal proposal distribution.
341
342
Standard random walk proposal using multivariate normal
343
distribution centered at current state.
344
"""
345
346
def __init__(self, s):
347
"""
348
Parameters:
349
- s: array, scaling matrix or standard deviations
350
"""
351
352
class CauchyProposal:
353
"""
354
Multivariate Cauchy proposal distribution.
355
356
Heavy-tailed proposals for robust exploration of
357
posterior distributions with long tails.
358
"""
359
360
class LaplaceProposal:
361
"""
362
Multivariate Laplace proposal distribution.
363
364
Symmetric double exponential proposals for sparse
365
or regularized parameter estimation.
366
"""
367
368
class PoissonProposal:
369
"""
370
Poisson proposal distribution for count variables.
371
372
Specialized proposals for discrete count parameters
373
using Poisson distribution.
374
"""
375
376
class MultivariateNormalProposal:
377
"""
378
Multivariate normal proposal with full covariance.
379
380
Correlated proposals that adapt to posterior geometry
381
for efficient exploration of complex parameter spaces.
382
"""
383
384
class UniformProposal:
385
"""
386
Uniform proposal distribution.
387
388
Uniform random proposals within specified bounds
389
for constrained parameter spaces.
390
"""
391
```
392
393
### Specialized Samplers
394
395
Advanced and specialized sampling methods for specific model types.
396
397
```python { .api }
398
class Slice:
399
"""
400
Slice sampler for univariate distributions.
401
402
Adaptive sampler that automatically adjusts to local
403
posterior structure without tuning parameters.
404
"""
405
406
def __init__(self, vars=None, w=1.0, tune=True, model=None, **kwargs):
407
"""
408
Initialize slice sampler.
409
410
Parameters:
411
- vars: list, variables to sample (one variable only)
412
- w: float, initial width parameter
413
- tune: bool, adapt width parameter during tuning
414
- model: Model, model context
415
"""
416
417
class EllipticalSlice:
418
"""
419
Elliptical slice sampler for Gaussian process priors.
420
421
Specialized sampler for models with Gaussian process
422
priors or multivariate normal latent variables.
423
"""
424
425
def __init__(self, vars=None, prior_cov=None, model=None, **kwargs):
426
"""
427
Initialize elliptical slice sampler.
428
429
Parameters:
430
- vars: list, variables with Gaussian prior
431
- prior_cov: array, prior covariance matrix
432
- model: Model, model context
433
"""
434
435
class ElemwiseCategorical:
436
"""
437
Element-wise categorical sampler.
438
439
Efficient Gibbs sampling for categorical variables
440
by updating each element independently.
441
"""
442
443
class PGBART:
444
"""
445
Particle Gibbs sampler for BART models.
446
447
Specialized sampler for Bayesian Additive Regression
448
Trees using particle Gibbs methods.
449
"""
450
```
451
452
### Compound Step Methods
453
454
Methods for combining multiple step methods and handling complex model structures.
455
456
```python { .api }
457
class CompoundStep:
458
"""
459
Compound step method combining multiple samplers.
460
461
Coordinates multiple step methods for models with
462
different variable types or sampling requirements.
463
"""
464
465
def __init__(self, methods):
466
"""
467
Initialize compound step method.
468
469
Parameters:
470
- methods: list, individual step method objects
471
"""
472
473
@property
474
def vars(self):
475
"""Variables handled by all component methods."""
476
477
@property
478
def generates_stats(self):
479
"""Whether any component generates sampling statistics."""
480
```
481
482
### Multi-Level Delayed Acceptance
483
484
Advanced methods for expensive likelihood evaluations and multi-level modeling.
485
486
```python { .api }
487
class MLDA:
488
"""
489
Multi-Level Delayed Acceptance sampler.
490
491
Hierarchical sampling method for models with expensive
492
likelihood evaluations using coarse approximations.
493
"""
494
495
def __init__(self, vars=None, coarse_models=None, base_sampler=None,
496
base_kwargs=None, model=None, **kwargs):
497
"""
498
Initialize MLDA sampler.
499
500
Parameters:
501
- vars: list, variables to sample
502
- coarse_models: list, coarse model approximations
503
- base_sampler: step method class for base level
504
- base_kwargs: dict, arguments for base sampler
505
- model: Model, fine model context
506
"""
507
508
class DEMetropolisZMLDA:
509
"""
510
DE Metropolis Z with Multi-Level Delayed Acceptance.
511
512
Combines differential evolution with multi-level
513
approximation for expensive population-based sampling.
514
"""
515
516
class MetropolisMLDA:
517
"""
518
Metropolis with Multi-Level Delayed Acceptance.
519
520
Multi-level Metropolis sampling using hierarchical
521
model approximations for computational efficiency.
522
"""
523
524
class RecursiveDAProposal:
525
"""
526
Recursive delayed acceptance proposal mechanism.
527
528
Implements recursive evaluation of proposals through
529
hierarchy of model approximations.
530
"""
531
```
532
533
## Usage Examples
534
535
### Basic MCMC Sampling
536
537
```python
538
import pymc3 as pm
539
import numpy as np
540
541
# Basic model sampling
542
with pm.Model() as model:
543
mu = pm.Normal('mu', mu=0, sigma=10)
544
sigma = pm.HalfNormal('sigma', sigma=5)
545
y_obs = pm.Normal('y_obs', mu=mu, sigma=sigma, observed=data)
546
547
# Sample with default settings (NUTS)
548
trace = pm.sample(draws=1000, tune=1000, chains=4)
549
550
# Sample with specific parameters
551
trace = pm.sample(draws=2000, tune=1500, target_accept=0.9,
552
chains=4, cores=2, random_seed=42)
553
```
554
555
### Custom Step Method Assignment
556
557
```python
558
# Manual step method assignment
559
with pm.Model() as mixed_model:
560
# Continuous parameters
561
theta = pm.Beta('theta', alpha=1, beta=1)
562
563
# Discrete parameters
564
k = pm.Poisson('k', mu=5)
565
566
# Binary parameter
567
switch = pm.Bernoulli('switch', p=0.5)
568
569
# Likelihood
570
y_obs = pm.NegativeBinomial('y_obs', mu=theta*k*switch,
571
alpha=1, observed=count_data)
572
573
# Assign specific step methods
574
step1 = pm.NUTS([theta]) # NUTS for continuous
575
step2 = pm.Metropolis([k]) # Metropolis for discrete
576
step3 = pm.BinaryMetropolis([switch]) # Binary Metropolis
577
578
trace = pm.sample(1000, step=[step1, step2, step3])
579
```
580
581
### Advanced NUTS Configuration
582
583
```python
584
# Customized NUTS sampler
585
with pm.Model() as nuts_model:
586
x = pm.Normal('x', mu=0, sigma=1, shape=10)
587
y_obs = pm.Normal('y_obs', mu=x.sum(), sigma=0.1, observed=5.0)
588
589
# NUTS with custom settings
590
nuts_step = pm.NUTS(vars=[x],
591
target_accept=0.95, # Higher acceptance rate
592
max_treedepth=12) # Deeper trees
593
594
trace = pm.sample(1000, step=nuts_step,
595
init='advi+adapt_diag_grad')
596
```
597
598
### Differential Evolution Metropolis
599
600
```python
601
# Population-based sampling with DE-MCMC
602
with pm.Model() as de_model:
603
# Multi-modal posterior
604
mu1 = pm.Normal('mu1', mu=-2, sigma=1)
605
mu2 = pm.Normal('mu2', mu=2, sigma=1)
606
607
# Mixture likelihood creating multi-modality
608
components = [pm.Normal.dist(mu=mu1, sigma=0.5),
609
pm.Normal.dist(mu=mu2, sigma=0.5)]
610
611
y_obs = pm.Mixture('y_obs', w=[0.3, 0.7],
612
comp_dists=components, observed=data)
613
614
# DE-MCMC for multi-modal exploration
615
step = pm.DEMetropolis(vars=[mu1, mu2], lamb=2.38)
616
617
# Multiple chains for population
618
trace = pm.sample(1000, step=step, chains=10, cores=1)
619
```
620
621
### Slice Sampling
622
623
```python
624
# Slice sampler for univariate parameters
625
with pm.Model() as slice_model:
626
# Parameter with complex posterior shape
627
tau = pm.Gamma('tau', alpha=1, beta=1)
628
629
# Hierarchical structure
630
sigmas = pm.InverseGamma('sigmas', alpha=2, beta=tau, shape=5)
631
632
y_obs = pm.Normal('y_obs', mu=0, sigma=sigmas, observed=group_data)
633
634
# Slice sampling for tau (univariate)
635
step_tau = pm.Slice([tau])
636
637
# NUTS for sigmas (multivariate)
638
step_sigmas = pm.NUTS([sigmas])
639
640
trace = pm.sample(1000, step=[step_tau, step_sigmas])
641
```
642
643
### Posterior Predictive Sampling
644
645
```python
646
# Comprehensive posterior predictive workflow
647
with pm.Model() as pred_model:
648
alpha = pm.Normal('alpha', mu=0, sigma=1)
649
beta = pm.Normal('beta', mu=0, sigma=1)
650
sigma = pm.HalfNormal('sigma', sigma=1)
651
652
mu = alpha + beta * x_data
653
y_obs = pm.Normal('y_obs', mu=mu, sigma=sigma, observed=y_data)
654
655
# Sample posterior
656
trace = pm.sample(1000, tune=1000)
657
658
# In-sample posterior predictive checks
659
ppc = pm.sample_posterior_predictive(trace, samples=500)
660
661
# Out-of-sample prediction with new data
662
pm.set_data({'x_data': x_new})
663
ppc_new = pm.sample_posterior_predictive(trace, samples=500)
664
665
# Custom variables for prediction
666
with pred_model:
667
# Add prediction variable
668
y_pred = pm.Normal('y_pred', mu=mu, sigma=sigma)
669
670
ppc_custom = pm.sample_posterior_predictive(trace,
671
var_names=['y_pred'],
672
samples=1000)
673
```
674
675
### Sequential and Online Sampling
676
677
```python
678
# Online/sequential sampling with iter_sample
679
with pm.Model() as online_model:
680
theta = pm.Beta('theta', alpha=1, beta=1)
681
682
# Initial sampling
683
step = pm.Metropolis([theta])
684
trace = pm.backends.NDArray(model=online_model)
685
686
# Sequential sampling with processing
687
for i, current_trace in enumerate(pm.iter_sample(100, step,
688
trace=trace,
689
tune=50)):
690
if i % 10 == 0:
691
# Process intermediate results
692
current_mean = current_trace['theta'].mean()
693
print(f"Step {i}: current mean = {current_mean:.3f}")
694
695
# Optional: early stopping based on convergence
696
if i > 50 and abs(current_mean - 0.5) < 0.01:
697
break
698
```
699
700
### Custom Initialization
701
702
```python
703
# Custom initialization strategies
704
with pm.Model() as init_model:
705
mu = pm.Normal('mu', mu=0, sigma=10)
706
sigma = pm.HalfNormal('sigma', sigma=5)
707
y_obs = pm.Normal('y_obs', mu=mu, sigma=sigma, observed=data)
708
709
# MAP initialization
710
map_estimate = pm.find_MAP()
711
trace_map = pm.sample(1000, init=map_estimate)
712
713
# ADVI initialization
714
trace_advi = pm.sample(1000, init='advi', n_init=50000)
715
716
# Custom starting values
717
start = {'mu': 2.0, 'sigma_log__': np.log(1.5)}
718
trace_custom = pm.sample(1000, init=start)
719
720
# Multiple chain initialization
721
n_chains = 4
722
starts = []
723
for i in range(n_chains):
724
starts.append({
725
'mu': np.random.normal(0, 2),
726
'sigma_log__': np.random.normal(0, 1)
727
})
728
729
trace_multi = pm.sample(1000, init=starts, chains=n_chains)
730
```