RJMCMC Algorithm Research Articles

Bayesian analysis of the COVID-19 pandemic using a Poisson process with change-points

Abstract Analyzing COVID-19 data presents a challenge in Bayesian computations of the Poisson process because the experimental conditions are not under control. This lack of homogeneity can lead to inconsistent model parameters, which violates the assumptions of Bayesian inference. In this paper, we study the multiple change-point detection problem from this viewpoint for a non-homogeneous sample path of the Poisson process as the response variable. The rate parameters are linked to some explanatory using a generalized linear model. The number of change-points is considered to be unknown as well as their locations. We introduce a Bayesian paradigm to estimate the number and location of change-points. We also present an adaptive RJMCMC algorithm to generate pseudo-random samples from the posterior distributions. We apply the proposed model to analyze the COVID-19 infection curves from different countries and identify patterns of cases. We also assess the efficacy of interventions, such as vaccination and public health emergency responses, implemented by different countries. The results of the analysis provide valuable insights into the spread of COVID-19 and the effectiveness of interventions. The proposed model can be used to inform public health decision-making and help to improve the management of the pandemic.

Bayesian multiple change-points detection in autocorrelated binary process with application to COVID-19 infection pattern

Maintaining consistent model parameters in long-run experiments is challenging. We study the multiple change-points detection problem from this viewpoint for an autocorrelated sample path of a binary response variable. The probabilities of success and correlation structure are linked to explanatory or control variables using a generalized linear model. A Bayesian paradigm is introduced here to estimate the unknown number and locations of change-points and an adaptive RJMCMC algorithm is presented to generate pseudo-random samples from the posterior distributions. This approach allows researchers to incorporate the experts' ideas into the analysis of such data. The features of the proposed method are examined via some simulation studies. We also apply the presented method to COVID-19 data from four countries to detect the change-points in the infection procedures in the presence of some explanatory variables.

Trans-Dimensional Bayesian Inversion of Rayleigh Wave Dispersion Curves for Shallow Surface Two-Dimensional Imaging

In this study, we achieve 2-D spatial structure directly from dispersion measurements using the reversible jump Markov chain Monte Carlo (rj-McMC) algorithm with a fully 2-D model parameterization. New parsimonious parameterization with low dimensional representation is introduced, which enable us to accurately reproduce the wave velocity field and the associated uncertainties. Field data test demonstrates that the proposed method estimates directly the shear-wave velocity profile with better lateral resolution. The associated uncertainties provide directly interpretable uncertainty on structures of interest.

Selecting and averaging relaxed clock models in Bayesian tip dating of Mesozoic birds

AbstractRelaxed clock models are fundamental in Bayesian clock dating, but a single distribution characterizing the clock variation is typically selected. Hence, I developed a new reversible-jump Markov chain Monte Carlo (rjMCMC) algorithm for drawing posterior samples between the independent lognormal (ILN) and independent gamma rates (IGR) clock models. The ability of the rjMCMC algorithm to infer the true model was verified through simulations. I then applied the algorithm to the Mesozoic bird data previously analyzed under the white noise (WN) clock model. In comparison, averaging over the ILN and IGR models provided more reliable estimates of the divergence times and evolutionary rates. The ILN model showed slightly better fit than the IGR model and much better fit than the autocorrelated lognormal (ALN) clock model. When the data were partitioned, different partitions showed heterogeneous model fit for ILN and IGR clocks. The implementation provides a general framework for selecting and averaging relaxed clock models in Bayesian dating analyses.

Gravity and magnetic joint inversion for basement and salt structures with the reversible-jump algorithm

SUMMARYGravity and magnetic data resolve the Earth with variable spatial resolution, and Earth structure exhibits both discontinuous and gradual features. Therefore, model parametrization complexity should be able to address such variability by locally adapting to the resolving power of the data. The reversible-jump Markov chain Monte Carlo (rjMcMC) algorithm provides variable spatial resolution that is consistent with data information. To address the prevalent non-uniqueness in joint inversion of potential field data, we use a novel spatial partitioning with nested Voronoi cells that is explored by rjMcMC sampling. The nested Voronoi parametrization partitions the subsurface in terms of rock types, such as sedimentary, salt and basement rocks. Therefore, meaningful prior information can be specified for each type which reduces non-uniqueness. We apply nonoverlapping prior distributions for density contrast and susceptibility between rock types. In addition, the choice of noise parametrization can lead to significant trade-offs with model resolution and complexity. We adopt an empirical estimation of full data covariance matrices that include theory and observational errors to account for spatially correlated noise. The method is applied to 2-D gravity and magnetic data to study salt and basement structures. We demonstrate that meaningful partitioning of the subsurface into sediment, salt, and basement structures is achieved by these advances without requiring regularization. Multiple simulated- and field-data examples are presented. Simulation results show clear delineation of salt and basement structures while resolving variable length scales. The field data show results that are consistent with observations made in the simulations. In particular, we resolve geologically plausible structures with varying length scales and clearly differentiate salt structure and basement topography.

AN ENERGY SEGMENTATION METHOD OF HIGH-RESOLUTION SAR IMAGE BASED ON MULTIPLE FEATURES

Abstract. To achieve the optimal image segmentation, an energy segmentation method based on multiple features and blocks for high resolution Synthetic Aperture Radar (SAR) image is proposed in this paper. First of all, a feature vector of pixel is formed with the texture feature extracted by curvelet transform and means function, the boundary feature extracted by curvelet transform and Canny Operator, and the original spectral feature; a feature set is formed by all feature vectors of pixels in the image. The feature vector is considered as segmentation basis, and its domain is partitioned by regular tessellation. On the partitioned image domain, a label variable is assigned to a regular block; each homogeneous region is fitted by one or more regular blocks; Obviously, a label field is formed by all the label variables of regular blocks. The model of label field is built by using energy function of neighborhood relationship. The feature set is considered as a realization of a random filed of multiple features (multiple features field for short). A heterogeneous energy function is used to establish the model of multiple features field. Then the established models of the label field and multiple features field are combined to define global energy function of image segmentation, and non-constrained Gibbs probability distribution is used to describe the global energy function to build the energy segmentation model based on multiple features. Further, a RJMCMC algorithm is designed to simulate from the model to segment SAR image. To verify the feasibility and superiority of the proposed approach, real SAR images are tested.

Application of Reversible Jump Markov Chain Monte Carlo Algorithms to Elastic and Petrophysical Amplitude-Versus-Angle Inversions

We infer the elastic and petrophysical properties from pre-stack seismic data through a transdimensional Bayesian inversion. In this approach the number of model parameters (i.e. the number of layers) is treated as an unknown, and a reversible jump Markov Chain Monte Carlo (rjMCMC) algorithm is used to sample the variable-dimension model space. This inversion scheme provides a parsimonious solution, and reliably quantifies the uncertainties affecting the estimated model parameters. Parallel tempering, which employs a sequence of interacting Markov chains in which the likelihood function is successively relaxed, is used to improve the efficiency of the probabilistic sampling. In addition, the delayed rejection updating scheme is employed to speed up the convergence of the rjMCMC algorithm to the stationary regime. Both elastic and petrophysical inversions invert the amplitude versus angle responses and employ a convolutional forward modelling based on the exact Zoeppritz equations. First, synthetic tests are used to assess the reliability of the implemented rjMCMC algorithms, then their applicability is demonstrated by inverting field seismic data acquired onshore. In this case the inversion was aimed at inferring the elastic and petrophysical properties around a gas-saturated reservoir hosted in a shale-sand sequence. In this case, the final outcomes provided by the rjMCMC algorithms are also compared with the predictions of linear Bayesian elastic and petrophysical inversions. The synthetic and field data examples demonstrate that the implemented algorithms can successfully estimate model uncertainty, model dimensionality and subsurface parameters with an affordable computational cost.

Convergence Tests for Transdimensional Markov Chains in Geoscience Imaging

Classic inversion methods adjust a model with a predefined number of parameters to the observed data. With transdimensional inversion algorithms such as the reversible-jump Markov Chain Monte Carlo (rjMCMC), it is possible to vary this number during the inversion and to interpret the observations in a more flexible way. Geoscience imaging applications use this behaviour to automatically adjust model resolution to the inhomogeneities of the investigated system, while keeping the model parameters on an optimal level. The rjMCMC algorithm produces an ensemble as result, a set of model realizations which together represent the posterior probability distribution of the investigated problem. The realizations are evolved via sequential updates from a randomly chosen initial solution, and converge toward the target posterior distribution of the inverse problem. Up to a point in the chain, the realizations may be strongly biased by the initial model, and have to be discarded from the final ensemble. With convergence assessment techniques, this point in the chain can be identified. Transdimensional MCMC methods produce ensembles which are not suitable for classic convergence assessment techniques because of the changes in parameter numbers. To overcome this hurdle, three solutions are introduced to convert model realizations to a common dimensionality while maintaining the statistical characteristics of the ensemble. A scalar, a vector and a matrix representation for models is presented, inferred from tomographic subsurface investigations, and three classic convergence assessment techniques are applied on them. It is shown that appropriately chosen scalar conversions of the models could retain similar statistical ensemble properties as geologic projections created by rasterization.

Prior specification for binary Markov mesh models

We propose prior distributions for all parts of the specification of a Markov mesh model. In the formulation we define priors for the sequential neighborhood, for the parametric form of the conditional distributions and for the parameter values. By simulating from the resulting posterior distribution when conditioning on an observed scene, we thereby obtain an automatic model selection procedure for Markov mesh models. To sample from such a posterior distribution, we construct a reversible jump Markov chain Monte Carlo algorithm (RJMCMC). We demonstrate the usefulness of our prior formulation and the limitations of our RJMCMC algorithm in two examples.

A reversible-jump Markov chain Monte Carlo algorithm for 1D inversion of magnetotelluric data

This paper presents a new computer code developed to solve the 1D magnetotelluric (MT) inverse problem using a Bayesian trans-dimensional Markov chain Monte Carlo algorithm. MT data are sensitive to the depth-distribution of rock electric conductivity (or its reciprocal, resistivity). The solution provided is a probability distribution - the so-called posterior probability distribution (PPD) for the conductivity at depth, together with the PPD of the interface depths. The PPD is sampled via a reversible-jump Markov Chain Monte Carlo (rjMcMC) algorithm, using a modified Metropolis-Hastings (MH) rule to accept or discard candidate models along the chains. As the optimal parameterization for the inversion process is generally unknown a trans-dimensional approach is used to allow the dataset itself to indicate the most probable number of parameters needed to sample the PPD. The algorithm is tested against two simulated datasets and a set of MT data acquired in the Clare Basin (County Clare, Ireland). For the simulated datasets the correct number of conductive layers at depth and the associated electrical conductivity values is retrieved, together with reasonable estimates of the uncertainties on the investigated parameters. Results from the inversion of field measurements are compared with results obtained using a deterministic method and with well-log data from a nearby borehole. The PPD is in good agreement with the well-log data, showing as a main structure a high conductive layer associated with the Clare Shale formation.In this study, we demonstrate that our new code go beyond algorithms developend using a linear inversion scheme, as it can be used: (1) to by-pass the subjective choices in the 1D parameterizations, i.e. the number of horizontal layers in the 1D parameterization, and (2) to estimate realistic uncertainties on the retrieved parameters. The algorithm is implemented using a simple MPI approach, where independent chains run on isolated CPU, to take full advantage of parallel computer architectures. In case of a large number of data, a master/slave appoach can be used, where the master CPU samples the parameter space and the slave CPUs compute forward solutions.

Bayesian Single Changepoint Estimation in a Parameter‐driven Model

AbstractIn this paper, we consider the problem of estimating a single changepoint in a parameter‐driven model. The model – an extension of the Poisson regression model – accounts for serial correlation through a latent process incorporated in its mean function. Emphasis is placed on the changepoint characterization with changes in the parameters of the model. The model is fully implemented within the Bayesian framework. We develop a RJMCMC algorithm for parameter estimation and model determination. The algorithm embeds well‐devised Metropolis–Hastings procedures for estimating the missing values of the latent process through data augmentation and the changepoint. The methodology is illustrated using data on monthly counts of claimants collecting wage loss benefit for injuries in the workplace and an analysis of presidential uses of force in the USA.

Transdimensional Love-wave tomography of the British Isles and shear-velocity structure of the East Irish Sea Basin from ambient-noise interferometry

We present the first Love-wave group velocity and shear velocity maps of the British Isles obtained from ambient noise interferometry and fully non-linear inversion. We computed interferometric inter-station Green's functions by cross-correlating the transverse component of ambient noise records retrieved by 61 seismic stations across the UK and Ireland. Group velocity measurements along each possible inter-station path were obtained using frequency-time analysis and converted into a series of inter-station traveltime datasets between 4 and 15 seconds period. Traveltime uncertainties estimated from the standard deviation of dispersion curves constructed by stacking randomly-selected subsets of daily cross-correlations, were observed to be too low to allow reasonable data fits to be obtained during tomography. Data uncertainties were therefore estimated again during the inversion as distance-dependent functionals. We produced Love-wave group velocity maps within 8 different period bands using a fully non-linear tomography method which combines the transdimensional reversible-jump Markov chain Monte Carlo (rj-McMC) algorithm with an eikonal raytracer. By modelling exact raypaths at each step of the Markov chain we ensured that the non-linear character of the inverse problem was fully and correctly accounted for. Between 4 and 10 seconds period, the group velocity maps show remarkable agreement with the known geology of the British Isles and correctly identify a number of low-velocity sedimentary basins and high-velocity features. Longer period maps, in which most sedimentary basins are not visible, are instead mainly representative of basement rocks. In a second stage of our study we used the results of tomography to produce a series of Love-wave group velocity dispersion curves across a grid of geographical points focussed around the East Irish Sea sedimentary basin. We then independently inverted each curve using a similar rj-McMC algorithm to obtain a series of one-dimensional shear velocity profiles. By merging all 1D profiles, we created a fully three-dimensional model of the crust beneath the East Irish Sea. The depth to basement in this model compares well with that averaged from seismic reflection profiles. This result is the first 3-dimensional model in the UK with fully quantified uncertainties: it shows basin depths and basement structures, and their concomitant uncertainties.

COUPLING REGULAR TESSELLATION WITH RJMCMC ALGORITHM TO SEGMENT SAR IMAGE WITH UNKNOWN NUMBER OF CLASSES

Abstract. This paper presents a Synthetic Aperture Radar (SAR) image segmentation approach with unknown number of classes, which is based on regular tessellation and Reversible Jump Markov Chain Monte Carlo (RJMCMC') algorithm. First of all, an image domain is portioned into a set of blocks by regular tessellation. The image is modeled on the assumption that intensities of its pixels in each homogeneous region satisfy an identical and independent Gamma distribution. By Bayesian paradigm, the posterior distribution is obtained to build the region-based image segmentation model. Then, a RJMCMC algorithm is designed to simulate from the segmentation model to determine the number of homogeneous regions and segment the image. In order to further improve the segmentation accuracy, a refined operation is performed. To illustrate the feasibility and effectiveness of the proposed approach, two real SAR image is tested.

COUPLING REGULAR TESSELLATION WITH RJMCMC ALGORITHM TO SEGMENT SAR IMAGE WITH UNKNOWN NUMBER OF CLASSES

This paper presents a Synthetic Aperture Radar (SAR) image segmentation approach with unknown number of classes, which is based on regular tessellation and Reversible Jump Markov Chain Monte Carlo (RJMCMC') algorithm. First of all, an image domain is portioned into a set of blocks by regular tessellation. The image is modeled on the assumption that intensities of its pixels in each homogeneous region satisfy an identical and independent Gamma distribution. By Bayesian paradigm, the posterior distribution is obtained to build the region-based image segmentation model. Then, a RJMCMC algorithm is designed to simulate from the segmentation model to determine the number of homogeneous regions and segment the image. In order to further improve the segmentation accuracy, a refined operation is performed. To illustrate the feasibility and effectiveness of the proposed approach, two real SAR image is tested.

Using hierarchical centering to facilitate a reversible jump MCMC algorithm for random effects models

Hierarchical centering has been described as a reparameterization method applicable to random effects models. It has been shown to improve mixing of models in the context of Markov chain Monte Carlo (MCMC) methods. A hierarchical centering approach is proposed for reversible jump MCMC (RJMCMC) chains which builds upon the hierarchical centering methods for MCMC chains and uses them to reparameterize models in an RJMCMC algorithm. Although these methods may be applicable to models with other error distributions, the case is described for a log-linear Poisson model where the expected value λ includes fixed effect covariates and a random effect for which normality is assumed with a zero-mean and unknown standard deviation. For the proposed RJMCMC algorithm including hierarchical centering, the models are reparameterized by modeling the mean of the random effect coefficients as a function of the intercept of the λ model and one or more of the available fixed effect covariates depending on the model. The method is appropriate when fixed-effect covariates are constant within random effect groups. This has an effect on the dynamics of the RJMCMC algorithm and improves model mixing. The methods are applied to a case study of point transects of indigo buntings where, without hierarchical centering, the RJMCMC algorithm had poor mixing and the estimated posterior distribution depended on the starting model. With hierarchical centering on the other hand, the chain moved freely over model and parameter space. These results are confirmed with a simulation study. Hence, the proposed methods should be considered as a regular strategy for implementing models with random effects in RJMCMC algorithms; they facilitate convergence of these algorithms and help avoid false inference on model parameters.

Segmentation of high-resolution SAR image with unknown number of classes based on regular tessellation and RJMCMC algorithm

This article presents a statistics- and region-based approach to segmentation of synthetic aperture radar (SAR) images. The proposed approach can automatically determine the number of classes and segment the image simultaneously. First of all, an image domain is partitioned into a set of blocks by regular tessellation and the image is modelled on the assumption that intensities of its pixels in each homogeneous region satisfy an identical and independent gamma distribution. The Bayesian paradigm is followed to build an image segmentation model. Then, a Reversible Jump Markov Chain Monte Carlo scheme is designed to simulate the segmentation model, which determines the number of classes and segments the image roughly. Furthermore, in order to improve the accuracy of the segmentation results, refined operation is performed. The results obtained from both real and simulated SAR images show that the proposed approach works well and efficient.

Efficient trans-dimensional Bayesian inversion for geoacoustic profile estimation

This paper considers sampling efficiency of trans-dimensional (trans-D) Bayesian inversion based on the reversible-jump Markov-chain Monte Carlo (rjMCMC) algorithm, with application to seabed acoustic reflectivity inversion. Trans-D inversion is applied to sample the posterior probability density over geoacoustic parameters for an unknown number of seabed layers, providing profile estimates with uncertainties that include the uncertainty in the model parameterization. However, the approach is computationally intensive. The efficiency of rjMCMC sampling is largely determined by the proposal schemes applied to perturb existing parameters and to assign values for parameters added to the model. Several proposal schemes are examined, some of which appear new for trans-D geoacoustic inversion. Perturbations of existing parameters are considered in a principal-component space based on an eigen-decomposition of the unit-lag parameter covariance matrix (computed from successive models along the Markov chain, a diminishing adaptation). The relative efficiency of proposing new parameters from the prior versus a Gaussian distribution focused near existing values is considered. Parallel tempering, which employs a sequence of interacting Markov chains with successively relaxed likelihoods, is also considered to increase the acceptance rate of new layers. The relative efficiency of various proposal schemes is compared through repeated inversions with a pragmatic convergence criterion.

Reliable and Fast Estimation of Recombination Rates by Convergence Diagnosis and Parallel Markov Chain Monte Carlo.

Genetic recombination is an essential event during the process of meiosis resulting in an exchange of segments between paired chromosomes. Estimating recombination rate is crucial for understanding the process of recombination. Experimental methods are normally difficult and limited to small scale estimations. Thus statistical methods using population genetics data are important for large-scale analysis. LDhat is an extensively used statistical method using rjMCMC algorithm to predict recombination rates. Due to the complexity of rjMCMC scheme, LDhat may take a long time for large SNP data sets. In addition, rjMCMC parameters should be manually defined in the original program which directly impact results. To address these issues, we designed an improved algorithm based on LDhat implementing MCMC convergence diagnostic algorithms to automatically predict values of parameters and monitor the mixing process. Then parallel computation methods were employed to further accelerate the new program. The new algorithms have been tested on ten samples from HapMap phase 2 data set. The results were compared with previous code and showed nearly identical output. However, our new methods achieved significant acceleration proving that they are more efficient and reliable for the estimation of recombination rates. The stand-alone package is freely available for download http://www.ntu.edu.sg/home/zhengjie/software/CPLDhat.

Hidden diversity in the Andes: Comparison of species delimitation methods in montane marsupials

Cryptic genetic diversity is a significant challenge for systematists faced with ever-increasing amounts of DNA sequence data. Computationally intensive coalescent-based analyses involving multiple unlinked loci are the only currently viable methods by which to assess the extent to which phenotypically similar populations (or metapopulations) are genetically distinct lineages. Although coalescent-based approaches have been tested extensively via simulations, few empirical studies have examined the impact of prior assumptions and dataset size on the ability to assess genetic isolation (evolutionary independence) using molecular data alone. Here, we consider the efficacy of two coalescent-based approaches (BPP and SpeDeSTEM) for testing the evolutionary independence of cryptic mtDNA haplogroups within three morphologically diagnosable species of Andean mouse opossums (Thylamys pallidior, T. sponsorius, and T. venustus). Fourteen anonymous nuclear loci, one X-linked nuclear intron, and one mitochondrial gene were analyzed for multiple individuals within each haplogroup of interest. We inferred individual gene trees for each locus and considered all of the nuclear loci jointly in a species-tree analysis. Using only the nuclear loci, we performed “species validation” tests for the cryptic mitochondrial lineages in SpeDeSTEM and BPP. For BPP, we also tested a wide range of prior assumptions, assessed performance of the rjMCMC algorithm, and examined how many loci were necessary to confidently delimit lineages. Results from BPP provided strong support for two independent evolutionary lineages each within T. pallidior, T. sponsorius, and T. venustus, whereas SpeDeSTEM results did not support splitting out mtDNA haplogroups as distinct evolutionary units. For most tests, BPP was robust to prior assumptions, although priors were shown to have an effect on both the strength of lineage recognition among T. venustus haplotypes and on the efficiency of the rjMCMC algorithm. Comparisons of results from datasets with different numbers of loci revealed that some cryptic lineages could be confidently delimited with as few as two loci.

Image segmentation using voronoi polygons and MCMC, with application to muscle fibre images

We investigate a Bayesian method for the segmentation of muscle fibre images. The images are reasonably well approximated by a Dirichlet tessellation, and so we use a deformable template model based on Voronoi polygons to represent the segmented image. We consider various prior distributions for the parameters and suggest an appropriate likelihood. Following the Bayesian paradigm, the mathematical form for the posterior distribution is obtained (up to an integrating constant). We introduce a Metropolis–Hastings algorithm and a reversible jump Markov chain Monte Carlo algorithm (RJMCMC) for simulation from the posterior when the number of polygons is fixed or unknown. The particular moves in the RJMCMC algorithm are birth, death and position/colour changes of the point process which determines the location of the polygons. Segmentation of the true image was carried out using the estimated posterior mode and posterior mean. A simulation study is presented which is helpful for tuning the hyperparameters and to assess the accuracy. The algorithms work well on a real image of a muscle fibre cross-section image, and an additional parameter, which models the boundaries of the muscle fibres, is included in the final model.