Variational Bayesian Blind Image Deconvolution: A review <s> D. Sampling methods <s> We make an analogy between images and statistical mechanics systems. Pixel gray levels and the presence and orientation of edges are viewed as states of atoms or molecules in a lattice-like physical system. The assignment of an energy function in the physical system determines its Gibbs distribution. Because of the Gibbs distribution, Markov random field (MRF) equivalence, this assignment also determines an MRF image model. The energy function is a more convenient and natural mechanism for embodying picture attributes than are the local characteristics of the MRF. For a range of degradation mechanisms, including blurring, non-linear deformations, and multiplicative or additive noise, the posterior distribution is an MRF with a structure akin to the image model. By the analogy, the posterior distribution defines another (imaginary) physical system. Gradual temperature reduction in the physical system isolates low-energy states (‘annealing’), or what is the same thing, the most probable states under the Gib... <s> BIB001 </s> Variational Bayesian Blind Image Deconvolution: A review <s> D. Sampling methods <s> We have sold 4300 copies worldwide of the first edition (1999). This new edition contains five completely new chapters covering new developments. <s> BIB002 </s> Variational Bayesian Blind Image Deconvolution: A review <s> D. Sampling methods <s> This paper presents a new Metropolis-adjusted Langevin algorithm (MALA) that uses convex analysis to simulate efficiently from high-dimensional densities that are log-concave, a class of probability distributions that is widely used in modern high-dimensional statistics and data analysis. The method is based on a new first-order approximation for Langevin diffusions that exploits log-concavity to construct Markov chains with favourable convergence properties. This approximation is closely related to Moreau-Yoshida regularisations for convex functions and uses proximity mappings instead of gradient mappings to approximate the continuous-time process. The proposed method complements existing MALA methods in two ways. First, the method is shown to have very robust stability properties and to converge geometrically for many target densities for which other MALA are not geometric, or only if the step size is sufficiently small. Second, the method can be applied to high-dimensional target densities that are not continuously differentiable, a class of distributions that is increasingly used in image processing and machine learning and that is beyond the scope of existing MALA and HMC algorithms. To use this method it is necessary to compute or to approximate efficiently the proximity mappings of the logarithm of the target density. For several popular models, including many Bayesian models used in modern signal and image processing and machine learning, this can be achieved with convex optimisation algorithms and with approximations based on proximal splitting techniques, which can be implemented in parallel. The proposed method is demonstrated on two challenging high-dimensional and non-differentiable models related to image resolution enhancement and low-rank matrix estimation that are not well addressed by existing MCMC methodology. <s> BIB003
Since the posterior distribution is not analytically available, sampling methods can be used to draw a large number of samples from it, and Monte Carlo integration techniques provide tools which allow performing inference on this dataset. To simulate the posterior distribution Markov chains are used to develop different sampling methods. Perhaps the most widely used sampling method is Gibbs sampling described by the Geman and Geman in BIB001 . More recent methods are the Metropolis adjusted Langevin algorithms BIB002 and Hamiltonian Monte Carlo  . To better understand the sampling methods let us see an example of Gibbs sampling. If we can write down analytic expressions for the conditional distributions of all the unknowns we wish to estimate, given the others, we simply draw samples from each of the distributions in turn, conditioned on the most recently generated samples values for the other parameters. In our case we want to simulate p(x, h, Ω|y); the iterations would proceed as follows:  , Ω (0) , y) Second iteration: x (2) ← p(x|h  , Ω  , y) h  ← p(h|x  , Ω  , y) where the symbol ← means that the values are drawn form the distribution on the right. Once enough samples have been collected, point estimates and other statistics of the distribution may be found using Monte Carlo integration, for example the Minimum Squared Error estimator of the mean can be obtained , where J is the number of drawn samples. Due to the expensive computational cost (which is even worse when it is applied to high-dimensional data, such as images), the use of sampling methods in BID is not very extended. The works in this field are focused on developing more efficient algorithms. Ge et al.  or Kail et al. [58] propose modified versions of the Gibbs sampling, and Pereyra BIB003 uses the Langevin algorithm which uses convex analysis to simulate efficiently the distributions.