HybridGibbs#

class cuqi.experimental.mcmc.HybridGibbs(target, sampling_strategy, num_sampling_steps=None)#

Hybrid Gibbs sampler for sampling a joint distribution.

Gibbs sampling samples the variables of the distribution sequentially, one variable at a time. When a variable represents a random vector, the whole vector is sampled simultaneously.

The sampling of each variable is done by sampling from the conditional distribution of that variable given the values of the other variables. This is often a very efficient way of sampling from a joint distribution if the conditional distributions are easy to sample from.

Hybrid Gibbs sampler is a generalization of the Gibbs sampler where the conditional distributions are sampled using different MCMC samplers.

When the conditionals are sampled exactly, the samples from the Gibbs sampler converge to the joint distribution. See e.g. Gelman et al. “Bayesian Data Analysis” (2014), Third Edition for more details.

In each Gibbs step, the corresponding sampler has the initial_point and initial_scale (if applicable) set to the value of the previous step and the sampler is reinitialized. This means that the sampling is not fully stateful at this point. This means samplers like NUTS will lose their internal state between Gibbs steps.

Parameters:

  • target (cuqi.distribution.JointDistribution) – Target distribution to sample from.

  • sampling_strategy (dict) – Dictionary of sampling strategies for each variable. Keys are variable names. Values are sampler objects.

  • num_sampling_steps (dict, optional) – Dictionary of number of sampling steps for each variable. The sampling steps are defined as the number of times the sampler will call its step method in each Gibbs step. Default is 1 for all variables.

Example

import cuqi
import numpy as np

# Model and data
A, y_obs, probinfo = cuqi.testproblem.Deconvolution1D(phantom='sinc').get_components()
n = A.domain_dim

# Define distributions
d = cuqi.distribution.Gamma(1, 1e-4)
l = cuqi.distribution.Gamma(1, 1e-4)
x = cuqi.distribution.GMRF(np.zeros(n), lambda d: d)
y = cuqi.distribution.Gaussian(A, lambda l: 1/l)

# Combine into a joint distribution and create posterior
joint = cuqi.distribution.JointDistribution(d, l, x, y)
posterior = joint(y=y_obs)

# Define sampling strategy
sampling_strategy = {
    'x': cuqi.experimental.mcmc.LinearRTO(maxit=15),
    'd': cuqi.experimental.mcmc.Conjugate(),
    'l': cuqi.experimental.mcmc.Conjugate(),
}

# Define Gibbs sampler
sampler = cuqi.experimental.mcmc.HybridGibbs(posterior, sampling_strategy)

# Run sampler
sampler.warmup(200)
sampler.sample(1000)

# Get samples removing burn-in
samples = sampler.get_samples().burnthin(200)

# Plot results
samples['x'].plot_ci(exact=probinfo.exactSolution)
samples['d'].plot_trace(figsize=(8,2))
samples['l'].plot_trace(figsize=(8,2))
__init__(target, sampling_strategy, num_sampling_steps=None)#

Methods

__init__(target, sampling_strategy[, ...])

get_samples()

sample(Ns)

Sample from the joint distribution using Gibbs sampling

step()

Sequentially go through all parameters and sample them conditionally on each other

tune(skip_len, update_count)

Run a single tuning step on each of the samplers in the Gibbs sampling scheme

validate_targets()

Validate each of the conditional targets used in the Gibbs steps

warmup(Nb[, tune_freq])

Warmup (tune) the samplers in the Gibbs sampling scheme