Maximum likelihood from spatial random effects models via the stochastic approximation expectation maximization algorithm

We introduce a class of spatial random effects models that have Markov random fields (MRF) as latent processes. Calculating the maximum likelihood estimates of unknown parameters in SREs is extremely difficult, because the normalizing factors of MRFs and additional integrations from unobserved random effects are computationally prohibitive. We propose a stochastic approximation expectation-maximization (SAEM) algorithm to maximize the likelihood functions of spatial random effects models. The SAEM algorithm integrates recent improvements in stochastic approximation algorithms; it also includes components of the Newton-Raphson algorithm and the expectation-maximization (EM) gradient algorithm. The convergence of the SAEM algorithm is guaranteed under some mild conditions. We apply the SAEM algorithm to three examples that are representative of real-world applications: a state space model, a noisy Ising model, and segmenting magnetic resonance images (MRI) of the human brain. The SAEM algorithm gives satisfactory results in finding the maximum likelihood estimate of spatial random effects models in each of these instances.     

Bayesian texture segmentation of weed and crop images using reversible jump Markov chain Monte Carlo methods

Summary. A Bayesian method for segmenting weed and crop textures is described and implemented. The work forms part of a project to identify weeds and crops in images so that selective crop spraying can be carried out. An image is subdivided into blocks and each block is modelled as a single texture. The number of different textures in the image is assumed unknown. A hierarchical Bayesian procedure is used where the texture labels have a Potts model (colour Ising Markov random field) prior and the pixels within a block are distributed according to a Gaussian Markov random field, with the parameters dependent on the type of texture. We simulate from the posterior distribution by using a reversible jump Metropolis–Hastings algorithm, where the number of different texture components is allowed to vary. The methodology is applied to a simulated image and then we carry out texture segmentation on the weed and crop images that motivated the work.     

A cluster‐sample approach for Monte Carlo integration using multiple samplers

A computational problem in many fields is to estimate simultaneously multiple integrals and expectations, assuming that the data are generated by some Monte Carlo algorithm. Consider two scenarios in which draws are simulated from multiple distributions but the normalizing constants of those distributions may be known or unknown. For each scenario, existing estimators can be classified as using individual samples separately or using all the samples jointly. The latter pooled‐sample estimators are statistically more efficient but computationally more costly to evaluate than the separate‐sample estimators. We develop a cluster‐sample approach to obtain computationally effective estimators, after draws are generated for each scenario. We divide all the samples into mutually exclusive clusters and combine samples from each cluster separately. Furthermore, we exploit a relationship between estimators based on samples from different clusters to achieve variance reduction. The resulting estimators, compared with the pooled‐sample estimators, typically yield similar statistical efficiency but have reduced computational cost. We illustrate the value of the new approach by two examples for an Ising model and a censored Gaussian random field. The Canadian Journal of Statistics 41: 151–173; 2013 © 2012 Statistical Society of Canada     

Fast sampling of Gaussian Markov random fields

This paper demonstrates how Gaussian Markov random fields (conditional autoregressions) can be sampled quickly by using numerical techniques for sparse matrices. The algorithm is general and efficient, and expands easily to various forms for conditional simulation and evaluation of normalization constants. We demonstrate its use by constructing efficient block updates in Markov chain Monte Carlo algorithms for disease mapping.     

Bayesian local influence analysis of skew-normal spatial dynamic panel data models

The existing studies on spatial dynamic panel data model (SDPDM) mainly focus on the normality assumption of response variables and random effects. This assumption may be inappropriate in some applications. This paper proposes a new SDPDM by assuming that response variables and random effects follow the multivariate skew-normal distribution. A Markov chain Monte Carlo algorithm is developed to evaluate Bayesian estimates of unknown parameters and random effects in skew-normal SDPDM by combining the Gibbs sampler and the Metropolis–Hastings algorithm. A Bayesian local influence analysis method is developed to simultaneously assess the effect of minor perturbations to the data, priors and sampling distributions. Simulation studies are conducted to investigate the finite-sample performance of the proposed methodologies. An example is illustrated by the proposed methodologies.     

Gaussian Markov random field spatial models in GAMLSS

This paper describes the modelling and fitting of Gaussian Markov random field spatial components within a Generalized AdditiveModel for Location, Scale and Shape (GAMLSS) model. This allows modelling of any or all the parameters of the distribution for the response variable using explanatory variables and spatial effects. The response variable distribution is allowed to be a non-exponential family distribution. A new package developed in R to achieve this is presented. We use Gaussian Markov random fields to model the spatial effect in Munich rent data and explore some features and characteristics of the data. The potential of using spatial analysis within GAMLSS is discussed. We argue that the flexibility of parametric distributions, ability to model all the parameters of the distribution and diagnostic tools of GAMLSS provide an ideal environment for modelling spatial features of data.     

Adaptive Gaussian Markov random fields with applications in human brain mapping

Summary.  Functional magnetic resonance imaging has become a standard technology in human brain mapping. Analyses of the massive spatiotemporal functional magnetic resonance imaging data sets often focus on parametric or non-parametric modelling of the temporal component, whereas spatial smoothing is based on Gaussian kernels or random fields. A weakness of Gaussian spatial smoothing is underestimation of activation peaks or blurring of high curvature transitions between activated and non-activated regions of the brain. To improve spatial adaptivity, we introduce a class of inhomogeneous Markov random fields with stochastic interaction weights in a space-varying coefficient model. For given weights, the random field is conditionally Gaussian, but marginally it is non-Gaussian. Fully Bayesian inference, including estimation of weights and variance parameters, can be carried out through efficient Markov chain Monte Carlo simulation. Although motivated by the analysis of functional magnetic resonance imaging data, the methodological development is general and can also be used for spatial smoothing and regression analysis of areal data on irregular lattices. An application to stylized artificial data and to real functional magnetic resonance imaging data from a visual stimulation experiment demonstrates the performance of our approach in comparison with Gaussian and robustified non-Gaussian Markov random-field models.     

Approximate computations for binary Markov random fields and their use in Bayesian models

Discrete Markov random fields form a natural class of models to represent images and spatial datasets. The use of such models is, however, hampered by a computationally intractable normalising constant. This makes parameter estimation and a fully Bayesian treatment of discrete Markov random fields difficult. We apply approximation theory for pseudo-Boolean functions to binary Markov random fields and construct approximations and upper and lower bounds for the associated computationally intractable normalising constant. As a by-product of this process we also get a partially ordered Markov model approximation of the binary Markov random field. We present numerical examples with both the pairwise interaction Ising model and with higher-order interaction models, showing the quality of our approximations and bounds. We also present simulation examples and one real data example demonstrating how the approximations and bounds can be applied for parameter estimation and to handle a fully Bayesian model computationally.     

The definition of tree-indexed Markov chains in random environment and their existence

This paper gives the definition of tree-indexed Markov chains in random environment with discrete state space, and then studies some equivalent theorems of tree-indexed Markov chains in random environment. Finally, we give the equivalence on tree-indexed Markov chains in Markov environment and double Markov chains indexed by a tree.     

On Exact Simulation of Markov Random Fields Using Coupling from the Past

A general framework for exact simulation of Markov random fields using the Propp–Wilson coupling from the past approach is proposed. Our emphasis is on situations lacking the monotonicity properties that have been exploited in previous studies. A critical aspect is the convergence time of the algorithm; this we study both theoretically and experimentically. Our main theoretical result in this direction says, roughly, that if interactions are sufficiently weak, then the expected running time of a carefully designed implementation is O ( N log N ), where N is the number of interacting components of the system. Computer experiments are carried out for random q -colourings and for the Widom–Rowlinson lattice gas model.     

Multiple-site updates in maximum a posteriori and marginal posterior modes image estimation

We describe standard single-site Monte Carlo Markov chain methods, the Hastings and Metropolis algorithms, the Gibbs sampler and simulated annealing, for maximum a posteriori and marginal posterior modes image estimation. These methods can experience great difficulty in traversing the whole image space in a finite time when the target distribution is multi-modal. We present a survey of multiple-site update methods, including Swendsen and Wang's algorithm, coupled Markov chains and cascade algorithms designed to tackle the problem of moving between modes of the posterior image distribution. We compare the performance of some of these algorithms for sampling from degraded and non-degraded Ising models     

Disagreement Loop and Path Creation/Annihilation Algorithms for Binary Planar Markov Fields with Applications to Image Segmentation

Abstract. We introduce a class of Gibbs–Markov random fields built on regular tessellations that can be understood as discrete counterparts of Arak–Surgailis polygonal fields. We focus first on consistent polygonal fields, for which we show consistency, Markovianity and solvability by means of dynamic representations. Next, we develop disagreement loop as well as path creation and annihilation dynamics for their general Gibbsian modifications, which cover most lattice‐based Gibbs–Markov random fields subject to certain mild conditions. Applications to foreground–background image segmentation problems are discussed.     

On the Value of derivative evaluations and random walk suppression in Markov Chain Monte Carlo algorithms

Two strategies that can potentially improve Markov Chain Monte Carlo algorithms are to use derivative evaluations of the target density, and to suppress random walk behaviour in the chain. The use of one or both of these strategies has been investigated in a few specific applications, but neither is used routinely. We undertake a broader evaluation of these techniques, with a view to assessing their utility for routine use. In addition to comparing different algorithms, we also compare two different ways in which the algorithms can be applied to a multivariate target distribution. Specifically, the univariate version of an algorithm can be applied repeatedly to one-dimensional conditional distributions, or the multivariate version can be applied directly to the target distribution.     

Prior specification of neighbourhood and interaction structure in binary Markov random fields

We formulate a prior distribution for the energy function of stationary binary Markov random fields (MRFs) defined on a rectangular lattice. In the prior we assign distributions to all parts of the energy function. In particular we define priors for the neighbourhood structure of the MRF, what interactions to include in the model, and for potential values. We define a reversible jump Markov chain Monte Carlo (RJMCMC) procedure to simulate from the corresponding posterior distribution when conditioned to an observed scene. Thereby we are able to learn both the neighbourhood structure and the parametric form of the MRF from the observed scene. We circumvent evaluations of the intractable normalising constant of the MRF when running the RJMCMC algorithm by adopting a previously defined approximate auxiliary variable algorithm. We demonstrate the usefulness of our prior in two simulation examples and one real data example.     

On the use of quasi-likelihood for estimation in hidden-Markov random fields

This paper illustrates the use of quasi-likelihood methods of inference for hidden Markov random fields. These are simple to use and can be employed under circumstances where only the model form and its covariance structure are specified. In particular they can be used to derive the same estimating equations as the E-M algorithm or change of measure methods, which make full distributional assumptions.     

Multiple window discrete scan statistic for higher-order Markovian sequences

Accurate and efficient methods to detect unusual clusters of abnormal activity are needed in many fields such as medicine and business. Often the size of clusters is unknown; hence, multiple (variable) window scan statistics are used to identify clusters using a set of different potential cluster sizes. We give an efficient method to compute the exact distribution of multiple window discrete scan statistics for higher-order, multi-state Markovian sequences. We define a Markov chain to efficiently keep track of probabilities needed to compute p-values for the statistic. The state space of the Markov chain is set up by a criterion developed to identify strings that are associated with observing the specified values of the statistic. Using our algorithm, we identify cases where the available approximations do not perform well. We demonstrate our methods by detecting unusual clusters of made free throw shots by National Basketball Association players during the 2009–2010 regular season.     

A sequential procedure for testing the existence of a random walk model in finite samples

Given the random walk model, we show, for the traditional unrestricted regression used in testing stationarity, that no matter what the initial value of the random walk is or its drift or its error standard deviation, the sampling distributions of certain statistics remain unchanged. Using Monte Carlo simulations, we estimate, for different finite samples, the sampling distributions of these statistics. After smoothing the percentiles of the empirical sampling distributions, we come up with a new set of critical values for testing the existence of a random walk, if each statistic is being used on an individual base. Combining the new sets of critical values, we finally suggest a general methodology for testing for a random walk model.     

Complete convergence for weighted sums of END random variables and its application to nonparametric regression models

In this article, the complete convergence for weighted sums of extended negatively dependent (END, for short) random variables is investigated. Some sufficient conditions for the complete convergence are provided. In addition, the Marcinkiewicz–Zygmund type strong law of large numbers for weighted sums of END random variables is obtained. The results obtained in the article generalise and improve the corresponding one of Wang et al. [(2014b), ‘On Complete Convergence for an Extended Negatively Dependent Sequence’, Communications in Statistics-Theory and Methods, 43, 2923–2937]. As an application, the complete consistency for the estimator of nonparametric regression model is established.     

Variational Bayes for estimating the parameters of a hidden Potts model

Hidden Markov random field models provide an appealing representation of images and other spatial problems. The drawback is that inference is not straightforward for these models as the normalisation constant for the likelihood is generally intractable except for very small observation sets. Variational methods are an emerging tool for Bayesian inference and they have already been successfully applied in other contexts. Focusing on the particular case of a hidden Potts model with Gaussian noise, we show how variational Bayesian methods can be applied to hidden Markov random field inference. To tackle the obstacle of the intractable normalising constant for the likelihood, we explore alternative estimation approaches for incorporation into the variational Bayes algorithm. We consider a pseudo-likelihood approach as well as the more recent reduced dependence approximation of the normalisation constant. To illustrate the effectiveness of these approaches we present empirical results from the analysis of simulated datasets. We also analyse a real dataset and compare results with those of previous analyses as well as those obtained from the recently developed auxiliary variable MCMC method and the recursive MCMC method. Our results show that the variational Bayesian analyses can be carried out much faster than the MCMC analyses and produce good estimates of model parameters. We also found that the reduced dependence approximation of the normalisation constant outperformed the pseudo-likelihood approximation in our analysis of real and synthetic datasets.     

A note on a universal random variate generator for integer-valued random variables

A universal generator for integer-valued square-integrable random variables is introduced. The generator relies on a rejection technique based on a generalization of the inversion formula for integer-valued random variables. This approach allows to create a dominating probability function, whose evaluation solely involves two integrals depending on the characteristic function of the random variable to be generated. The proposal gives rise to a simple algorithm which may be implemented in a few code lines and which may show good performance when the classical families of distributions—such as the Poisson and the Binomial—are considered. In addition, applications to the Poisson-Tweedie and the Luria-Delbrück distributions are provided.