API

Implementation of Markov chain Monte Carlo (MCMC) approaches to sample from multivariate Gaussian distributions.

class pygauss.mcmc_sampling.sampler_MS(mu, Q, ini, b, band=True, seed=None, size=1)[source]

Algorithm dedicated to sample from a multivariate real-valued Gaussian distribution \(\mathcal{N}(\boldsymbol{\mu},\mathbf{Q}^{-1})\) where \(\mathbf{Q}\) is a symmetric and positive definite precision matrix. We assume here that the matrix splitting scheme \(\mathbf{Q} = \mathbf{M} - \mathbf{N}\) holds.

__init__(mu, Q, ini, b, band=True, seed=None, size=1)[source]
Parameters:
  • mu (1-D array_like, of length d) –
  • Q (2-D array_like, of shape (d,d)) – Precision matrix.
  • ini (1-D array_like, of length d. Initialization of the Markov chain.) –
  • b (int) – Bandwidth of the precision matrix Q.
  • band (boolean, optional) – Indicates if the precision matrix is banded with bandwidth b.
  • seed (int, optional) – Random seed to reproduce experimental results.
  • size (int, optional) – Given a size of for instance T, T independent and identically distributed (i.i.d.) samples are returned.
exact_MS(method='Gauss-Seidel')[source]

The samplers considered here are exact.

Parameters:method (string, optional) – Matrix splitting approach to choose within [‘Gauss-Seidel’,’Richardson’,’Jacobi’,’SOR’,’SSOR’,’Cheby-SSOR’].
Returns:theta – The drawn samples, of shape (d,size), if that was provided. If not, the shape is (d,1).
Return type:ndarray, of shape (d,size)

Examples

>>> import mcmc_sampling as mcmc
>>> d = 10
>>> mu = np.zeros(d)
>>> ini = np.zeros(d)
>>> Q = np.eye(d)
>>> b = 1
>>> band = True
>>> S = mcmc.sampler_MS(mu,Q,ini=ini,b=b,band=True,seed=2022,size=1)
>>> theta = S.exact_MS(method="Gauss-Seidel")
approx_MS(method='Clone-MCMC', omega=1)[source]

The samplers considered here are approximate.

Parameters:
  • method (string, optional) – Matrix splitting approach to choose within [‘Clone-MCMC’,’Hogwild’].
  • omega (float, optional) – Tuning parameter appearing in some approximate matrix splitting Gibbs samplers.
Returns:

theta – The drawn samples, of shape (d,size), if that was provided. If not, the shape is (d,1).

Return type:

ndarray, of shape (d,size)

Examples

>>> import mcmc_sampling as mcmc
>>> d = 10
>>> mu = np.zeros(d)
>>> ini = np.zeros(d)
>>> Q = np.eye(d)
>>> b = 1
>>> band = True
>>> S = mcmc.sampler_MS(mu,Q,ini=ini,b=b,band=True,seed=2022,size=1)
>>> theta = S.approx_MS(method="Gauss-Seidel",omega=1)
class pygauss.mcmc_sampling.sampler_DA(mu, seed=None, size=1)[source]

Algorithm dedicated to sample from a multivariate real-valued Gaussian distribution \(\mathcal{N}(\boldsymbol{\mu},\mathbf{Q}^{-1})\) where \(\mathbf{Q}\) is a symmetric and positive definite precision matrix. We assume here that \(\mathbf{Q} = \mathbf{G}_1^T\mathbf{\Lambda}_1^{-1}\mathbf{G}_1 + \mathbf{G}_2^T\mathbf{\Lambda}_2^{-1}\mathbf{G}_2\). Sampling from the corresponding multivariate Gaussian distribution is done with an MCMC algorithm based on a data augmentation scheme.

__init__(mu, seed=None, size=1)[source]
Parameters:
  • mu (1-D array_like, of length d) –
  • seed (int, optional) – Random seed to reproduce experimental results.
  • size (int, optional) – Given a size of for instance T, T independent and identically distributed (i.i.d.) samples are returned.
exact_DA_circu_diag_band(Lamb1, g, M, N, Q2, b2, method='GEDA')[source]

The samplers considered here are exact. We further assume here that \(\mathbf{G}_1\) is a circulant matrix, \(\mathbf{\Lambda}_1\) is diagonal, \(\mathbf{G}_2\) is the identity matrix and \(\mathbf{Q}_2 = \mathbf{\Lambda}_2^{-1}\) is a band matrix.

Parameters:
  • Lamb1 (1-D array_like, of length d) – Diagonal elements of \(\mathbf{\Lambda}_1\).
  • g (2-D array_like, of shape (N, M)) – Vector built by stacking the first columns associated to the \(M\) blocks of size \(N\) of the matrix \(\mathbf{G}_1\).
  • M (int) – Number of different blocks in \(\mathbf{G}_1\).
  • N (int) – Dimension of each block in \(\mathbf{G}_1\).
  • Q2 (2-D array_like, of shape (d, d)) – Precision matrix \(\mathbf{Q}_2\).
  • b2 (int) – Bandwidth of \(\mathbf{Q}_2\).
  • method (string, optional) – Data augmentation approach to choose within [‘EDA’,’GEDA’].
Returns:

theta – The drawn samples, of shape (d,size), if that was provided. If not, the shape is (d,1).

Return type:

ndarray, of shape (d,size)

Examples

>>> import mcmc_sampling as mcmc
>>> d = 15
>>> Lamb1 = np.random.normal(2,0.1,d)
>>> g =  np.reshape(np.random.normal(2,0.1,d),(d,1))
>>> M = 1
>>> N = d
>>> Q2 = np.diag(np.random.normal(2,0.1,d))
>>> b2 = 0
>>> S = mcmc.sampler_DA(mu,seed=2022,size=1)
>>> theta = S.exact_DA_circu_diag_band(Lamb1,g,M,N,
                                       Q2,b2,method="EDA")
exact_DA_circu_diag_circu(Lamb1, LambG1, LambG2, A, method='GEDA')[source]

The samplers considered here are exact. We further assume here that \(\mathbf{G}_1\) is a circulant matrix, \(\mathbf{\Lambda}_1\) is diagonal, \(\mathbf{\Lambda}_2\) is the identity matrix and \(\mathbf{G}_2\) is a circulant matrix.

Parameters:
  • Lamb1 (1-D array_like, of length d) – Diagonal elements of \(\mathbf{\Lambda}_1\).
  • LambG1 (1-D array_like, of length d) – Diagonal elements of the Fourier counterpart matrix associated to the matrix \(\mathbf{G}_1\).
  • LambG2 (1-D array_like, of length d) – Diagonal elements of the Fourier counterpart matrix associated to the matrix \(\mathbf{G}_2\).
  • A (function) – Linear operator returning the matrix-vector product \(\mathbf{Qx}\) where \(\mathbf{x} \in \mathbb{R}^d\).
  • method (string, optional) – Data augmentation approach to choose within [‘EDA’,’GEDA’].
Returns:

theta – The drawn samples, of shape (d,size), if that was provided. If not, the shape is (d,1).

Return type:

ndarray, of shape (d,size)

Examples

>>> import mcmc_sampling as mcmc
>>> d = 15
>>> Lamb1 = np.random.normal(2,0.1,d)
>>> g =  np.reshape(np.random.normal(2,0.1,d),(d,1))
>>> M = 1
>>> N = d
>>> Q2 = np.diag(np.random.normal(2,0.1,d))
>>> b2 = 0
>>> S = mcmc.sampler_DA(mu,seed=2022,size=1)
>>> theta = S.exact_DA_circu_diag_band(Lamb1,g,M,N,
                                       Q2,b2,method="EDA")