Markov Chain Monte Carlo Implementation Research Articles

Robust Constraint of Luminosity Function Evolution through MCMC Sampling

AbstractWe present a new galaxy survey simulation package, which combines the power of Markov Chain Monte Carlo (MCMC) sampling with a robust and adaptable model of galaxy evolution. The aim of this code is to aid in the characterization and study of new and existing galaxy surveys. In this paper we briefly describe the MCMC implementation and the survey simulation methodology and associated tools. A test case of this full suite was to constrain the evolution of the IR Luminosity Function (LF) based on the HerMES (Herschel SPIRE) survey of the Spitzer First Look Survey field. The initial results are consistent with previous studies, but our more general approach should be of wider benefit to the community.

A BAYESIAN METHOD FOR THE INTERCALIBRATION OF SPECTRA IN REVERBERATION MAPPING

Flux calibration of spectra in reverberation mapping (RM) is most often performed by assuming the flux constancy of some specified narrow emission lines, which stem from an extended region that is sometimes partially spatially resolved, in contrast to the point-like broad-line region and the central continuum source. The inhomogeneous aperture geometries used among different observation sets in a joint monitoring campaign introduce systematic deviations to the fluxes of broad lines and central continuum, and intercalibration over these data sets is required. As an improvement to the previous empirical correction performed by comparing the (nearly) contemporaneous observation points, we describe a feasible Bayesian method that obviates the need for (nearly) contemporaneous observations, naturally incorporates physical models of flux variations, and fully takes into account the measurement errors. In particular, it fits all the data sets simultaneously regardless of samplings and makes use of all of the information in the data sets. A Markov Chain Monte Carlo implementation is employed to recover the parameters and uncertainties for intercalibration. Application to the RM data sets of NGC 5548 with joint monitoring shows the high fidelity of our method.

A Bayesian Algorithm for Anxiety Index Prediction Based on Cerebral Blood Oxygenation in the Prefrontal Cortex Measured by Near Infrared Spectroscopy.

Stress-induced psychological and somatic diseases are virtually endemic nowadays. Written self-report anxiety measures are available; however, these indices tend to be time consuming to acquire. For medical patients, completing written reports can be burdensome if they are weak, in pain, or in acute anxiety states. Consequently, simple and fast non-invasive methods for assessing stress response from neurophysiological data are essential. In this paper, we report on a study that makes predictions of the state-trait anxiety inventory (STAI) index from oxyhemoglobin and deoxyhemoglobin concentration changes of the prefrontal cortex using a two-channel portable near-infrared spectroscopy device. Predictions are achieved by constructing machine learning algorithms within a Bayesian framework with nonlinear basis function together with Markov Chain Monte Carlo implementation. In this paper, prediction experiments were performed against four different data sets, i.e., two comprising young subjects, and the remaining two comprising elderly subjects. The number of subjects in each data set varied between 17 and 20 and each subject participated only once. They were not asked to perform any task; instead, they were at rest. The root mean square errors for the four groups were 6.20, 6.62, 4.50, and 6.38, respectively. There appeared to be no significant distinctions of prediction accuracies between age groups and since the STAI are defined between 20 and 80, the predictions appeared reasonably accurate. The results indicate potential applications to practical situations such as stress management and medical practice.

A BAYESIAN APPROACH TO ESTIMATE THE SIZE AND STRUCTURE OF THE BROAD-LINE REGION IN ACTIVE GALACTIC NUCLEI USING REVERBERATION MAPPING DATA

This is the first paper in a series devoted to systematic study of the size and structure of the broad-line region (BLR) in active galactic nuclei (AGNs) using reverberation mapping (RM) data. We employ a recently developed Bayesian approach that statistically describes the variabibility as a damped random walk process and delineates the BLR structure using a flexible disk geometry that can account for a variety of shapes, including disks, rings, shells, and spheres. We allow for the possibility that the line emission may respond non-linearly to the continuum, and we detrend the light curves when there is clear evidence for secular variation. We use a Markov Chain Monte Carlo implementation based on Bayesian statistics to recover the parameters and uncertainties for the BLR model. The corresponding transfer function is obtained self-consistently. We tentatively constrain the virial factor used to estimate black hole masses; more accurate determinations will have to await velocity-resolved RM data. Application of our method to RM data with Hbeta monitoring for about 40 objects shows that the assumed BLR geometry can reproduce quite well the observed emission-line fluxes from the continuum light curves. We find that the Hbeta BLR sizes obtained from our method are on average ~20% larger than those derived from the traditional cross-correlation method. Nevertheless, we still find a tight BLR size-luminosity relation with a slope of alpha=0.55\pm0.03 and an intrinsic scatter of ~0.18 dex. In particular, we demonstrate that our approach yields appropriate BLR sizes for some objects (such as Mrk 142 and PG 2130+099) where traditional methods previously encountered difficulties.

Bayesian random threshold estimation in a Cox proportional hazards cure model

In this paper, we develop a Bayesian approach to estimate a Cox proportional hazards model that allows a threshold in the regression coefficient, when some fraction of subjects are not susceptible to the event of interest. A data augmentation scheme with latent binary cure indicators is adopted to simplify the Markov chain Monte Carlo implementation. Given the binary cure indicators, the Cox cure model reduces to a standard Cox model and a logistic regression model. Furthermore, the threshold detection problem reverts to a threshold problem in a regular Cox model. The baseline cumulative hazard for the Cox model is formulated non-parametrically using counting processes with a gamma process prior. Simulation studies demonstrate that the method provides accurate point and interval estimates. Application to a data set of oropharynx cancer patients suggests a significant threshold in age at diagnosis such that the effect of gender on disease-specific survival changes after the threshold.

Markov Chain Monte Carlo: an introduction for epidemiologists

Markov Chain Monte Carlo (MCMC) methods are increasingly popular among epidemiologists. The reason for this may in part be that MCMC offers an appealing approach to handling some difficult types of analyses. Additionally, MCMC methods are those most commonly used for Bayesian analysis. However, epidemiologists are still largely unfamiliar with MCMC. They may lack familiarity either with he implementation of MCMC or with interpretation of the resultant output. As with tutorials outlining the calculus behind maximum likelihood in previous decades, a simple description of the machinery of MCMC is needed. We provide an introduction to conducting analyses with MCMC, and show that, given the same data and under certain model specifications, the results of an MCMC simulation match those of methods based on standard maximum-likelihood estimation (MLE). In addition, we highlight examples of instances in which MCMC approaches to data analysis provide a clear advantage over MLE. We hope that this brief tutorial will encourage epidemiologists to consider MCMC approaches as part of their analytic tool-kit.

Optimal Scaling of the Random Walk Metropolis: General Criteria for the 0.234 Acceptance Rule

Scaling of proposals for Metropolis algorithms is an important practical problem in Markov chain Monte Carlo implementation. Analyses of the random walk Metropolis for high-dimensional targets with specific functional forms have shown that in many cases the optimal scaling is achieved when the acceptance rate is approximately 0.234, but that there are exceptions. We present a general set of sufficient conditions which are invariant to orthonormal transformation of the coordinate axes and which ensure that the limiting optimal acceptance rate is 0.234. The criteria are shown to hold for the joint distribution of successive elements of a stationary pth-order multivariate Markov process.

Optimal Scaling of the Random Walk Metropolis: General Criteria for the 0.234 Acceptance Rule

Scaling of proposals for Metropolis algorithms is an important practical problem in Markov chain Monte Carlo implementation. Analyses of the random walk Metropolis for high-dimensional targets with specific functional forms have shown that in many cases the optimal scaling is achieved when the acceptance rate is approximately 0.234, but that there are exceptions. We present a general set of sufficient conditions which are invariant to orthonormal transformation of the coordinate axes and which ensure that the limiting optimal acceptance rate is 0.234. The criteria are shown to hold for the joint distribution of successive elements of a stationary pth-order multivariate Markov process.

Population Forecasts for Bangladesh, Using a Bayesian Methodology

Population projection for many developing countries could be quite a challenging task for the demographers mostly due to lack of availability of enough reliable data. The objective of this paper is to present an overview of the existing methods for population forecasting and to propose an alternative based on the Bayesian statistics, combining the formality of inference. The analysis has been made using Markov Chain Monte Carlo (MCMC) technique for Bayesian methodology available with the software WinBUGS. Convergence diagnostic techniques available with the WinBUGS software have been applied to ensure the convergence of the chains necessary for the implementation of MCMC. The Bayesian approach allows for the use of observed data and expert judgements by means of appropriate priors, and a more realistic population forecasts, along with associated uncertainty, has been possible.

Improving Efficiency of the Bayesian Approach to Water Distribution Contaminant Source Characterization with Support Vector Regression

There are multiple sources of uncertainties in urban water-distribution systems, e.g., nodal water demand and sensor measurement error. All of these uncertainties increase the complexity of contaminant source identification in a sparse sensor network. The large number of attributes (e.g., contaminant source location, magnitude, injection starting time, and duration) of a contaminant event profile cannot be identified given limited sensor data. Instead, the uncertainties in the contaminant event profile need to be characterized. Markov chain Monte Carlo (MCMC) methods for Bayesian analyses allow for the characterization of the uncertainty in the contamination event profile. To account for stochastic water demands, which has been shown in some circumstances to be necessary if the contaminant event is to be properly characterized, the evaluation of the likelihood function is the most computationally expensive part of the MCMC implementation. Previous work applied Monte Carlo methods for error propagation (MCEP). The research reported in this paper investigates the application of support vector regression (SVR) to speed the likelihood evaluation. This coupled MCMC-SVR approach enables probabilistic inference of the contaminant event. An SVR model, which maps from the contaminant event space to the likelihood space, is built for each node in the network to evaluate the likelihood function during the MCMC chain evolution. For the case study investigated, MCMC-SVR is computationally feasible and robust in inferring the contaminant event. A comparison between MCMC-SVR and MCMC-MCEP reveals that there is no substantial difference between the inferences provided by the two models, whereas the former is more computationally efficient.

Assessing Convergence and Mixing of MCMC Implementations via Stratification

Some posterior distributions lead to Markov chain Monte Carlo (MCMC) chains that are naturally viewed as collections of subchains. Examples include mixture models, regime-switching models, and hidden Markov models. We obtain MCMC-based estimators of posterior expectations by combining different subgroup (subchain) estimators using stratification and poststratification methods. Variance estimates of the limiting distributions of such estimators are developed. Based on these variance estimates, we propose a test statistic to aid in the assessment of convergence and mixing of chains. We compare our diagnostic with other commonly used methods. The approach is illustrated in two examples: a latent variable model for arsenic concentration in public water systems in Arizona and a Bayesian hierarchical model for Pacific sea surface temperatures. Supplementary materials, which include MATLAB codes for the proposed method, are available online.

Forecasting size-fractionated particle number concentrations in the urban atmosphere

Abstract Airborne particulate matter affects human health, especially in urban areas where air pollutant concentrations are high. In order to reduce exposure to particulates and other pollutants it is essential to forecast concentrations of these. In this work we introduce a statistical model to forecast size-fractionated particle number concentrations. Our forecasting approach is based on a parametric regression model that utilises traffic intensity and meteorological parameters as covariates and has an autoregressive dependence structure for the error terms. We use a Bayesian framework with a Markov Chain Monte Carlo (MCMC) implementation to derive the forecasts numerically from a given set of learning data. The particle concentration forecast is provided as probability distributions for a few days ahead. For the development and test of the model we used weather and particle size distribution data from an urban background station in Helsinki, and traffic data. Comparison of forecast distributions and measurements shows that the adopted probabilistic characterisation of particle number concentrations is adequate. In particular, forecasts of the log number concentration of ultra-fine particles (diameter R 2 equal to 0.67 (at 3 h time resolution). Our model is flexible and it may be implemented for other urban locations, provided that measurements of number concentrations, as well as measurements and forecasts of covariates are available.

An adaptive sequential Monte Carlo method for approximate Bayesian computation

Approximate Bayesian computation (ABC) is a popular approach to address inference problems where the likelihood function is intractable, or expensive to calculate. To improve over Markov chain Monte Carlo (MCMC) implementations of ABC, the use of sequential Monte Carlo (SMC) methods has recently been suggested. Most effective SMC algorithms that are currently available for ABC have a computational complexity that is quadratic in the number of Monte Carlo samples (Beaumont et al., Biometrika 86:983---990, 2009; Peters et al., Technical report, 2008; Toni et al., J. Roy. Soc. Interface 6:187---202, 2009) and require the careful choice of simulation parameters. In this article an adaptive SMC algorithm is proposed which admits a computational complexity that is linear in the number of samples and adaptively determines the simulation parameters. We demonstrate our algorithm on a toy example and on a birth-death-mutation model arising in epidemiology.

On thinning of chains in MCMC

Summary1. Markov chain Monte Carlo (MCMC) is a simulation technique that has revolutionised the analysis of ecological data, allowing the fitting of complex models in a Bayesian framework. Since 2001, there have been nearly 200 papers using MCMC in publications of the Ecological Society of America and the British Ecological Society, including more than 75 in the journal Ecology and 35 in the Journal of Applied Ecology.2. We have noted that many authors routinely ‘thin’ their simulations, discarding all but every kth sampled value; of the studies we surveyed with details on MCMC implementation, 40% reported thinning.3. Thinning is often unnecessary and always inefficient, reducing the precision with which features of the Markov chain are summarised. The inefficiency of thinning MCMC output has been known since the early 1990’s, long before MCMC appeared in ecological publications.4. We discuss the background and prevalence of thinning, illustrate its consequences, discuss circumstances when it might be regarded as a reasonable option and recommend against routine thinning of chains unless necessitated by computer memory limitations.

Analyzing spatial point patterns subject to measurement error

We address the issue of inference for a noisy point pattern. The unobserved true point process is modelled as a nonhomogeneous Poisson process. For modeling the underlying intensity surface we use a scaled Gaussian mixture distribution. The noise that creeps in during the measurement procedure causes random displacement of the true locations. We consider two settings. With a bounded region of interest, (i) this displacement may cause a true location within the boundary to be associated with an ‘observed’ location outside of the region and thus missed and (ii) we have the possibility in (i) but also vice versa; the displacement may bring in an observed location whose true location lies outside the region. Under (i), we can only lose points and, depending on the variability in the measurement error as well as the number of true locations close to boundary, this can cause a significant number of locations to be lost from our recorded set of data. Estimation of the intensity surface from the observed data can be misleading especially near the boundary of our domain of interest. Under (ii), the modeling problem is more difficult; points can be both lost and gained and it is challenging to characterize how we may gain points with no data on the underlying intensity outside the domain of interest. In both cases, we work within a hierarchical Bayes framework, modeling the latent point pattern using a Cox process and, given the process realization, introducing a suitable measurement error model. Hence, the specification includes the true number of points as an unknown. We discuss choice of measurement error model as well as identifiability problems which arise. Models are fitted using an markov chain Monte Carlo implementation. After validating our method against several synthetic datasets we illustrate its application for two ecological datasets.

Bayesian estimation of weak material dispersion: theory and experiment

This work considers the estimation of dispersion in materials via an interferometric technique. At its core, the problem involves extracting the quadratic variation in phase over a range of wavelengths based on measured optical intensity. The estimation problem becomes extremely difficult for weakly dispersive materials where the quadratic nonlinearity is very small relative to the uncertainty inherent in experiment. This work provides a means of estimating dispersion in the face of such uncertainty. Specifically, we use a Markov Chain Monte Carlo implementation of Bayesian analysis to provide both the dispersion estimate and the associated confidence interval. The interplay between various system parameters and the size of the resulting confidence interval is discussed. The approach is then applied to several different experimental samples.

A robust multivariate measurement error model with skew-normal/independent distributions and Bayesian MCMC implementation

Skew-normal/independent distributions are a class of asymmetric thick-tailed distributions that include the skew-normal distribution as a special case. In this paper, we explore the use of Markov Chain Monte Carlo (MCMC) methods to develop a Bayesian analysis in multivariate measurement errors models. We propose the use of skew-normal/independent distributions to model the unobserved value of the covariates (latent variable) and symmetric normal/independent distributions for the random errors term, providing an appealing robust alternative to the usual symmetric process in multivariate measurement errors models. Among the distributions that belong to this class of distributions, we examine univariate and multivariate versions of the skew-normal, skew- t , skew-slash and skew-contaminated normal distributions. The results and methods are applied to a real data set.

Monte Carlo Inference for State–Space Models of Wild Animal Populations

We compare two Monte Carlo (MC) procedures, sequential importance sampling (SIS) and Markov chain Monte Carlo (MCMC), for making Bayesian inferences about the unknown states and parameters of state-space models for animal populations. The procedures were applied to both simulated and real pup count data for the British grey seal metapopulation, as well as to simulated data for a Chinook salmon population. The MCMC implementation was based on tailor-made proposal distributions combined with analytical integration of some of the states and parameters. SIS was implemented in a more generic fashion. For the same computing time MCMC tended to yield posterior distributions with less MC variation across different runs of the algorithm than the SIS implementation with the exception in the seal model of some states and one of the parameters that mixed quite slowly. The efficiency of the SIS sampler greatly increased by analytically integrating out unknown parameters in the observation model. We consider that a careful implementation of MCMC for cases where data are informative relative to the priors sets the gold standard, but that SIS samplers are a viable alternative that can be programmed more quickly. Our SIS implementation is particularly competitive in situations where the data are relatively uninformative; in other cases, SIS may require substantially more computer power than an efficient implementation of MCMC to achieve the same level of MC error.

Identifiability of parameters and behaviour of MCMC chains: a case study using the reaction norm model

Markov chain Monte Carlo (MCMC) enables fitting complex hierarchical models that may adequately reflect the process of data generation. Some of these models may contain more parameters than can be uniquely inferred from the distribution of the data, causing non-identifiability. The reaction norm model with unknown covariates (RNUC) is a model in which unknown environmental effects can be inferred jointly with the remaining parameters. The problem of identifiability of parameters at the level of the likelihood and the associated behaviour of MCMC chains were discussed using the RNUC as an example. It was shown theoretically that when environmental effects (covariates) are considered as random effects, estimable functions of the fixed effects, (co)variance components and genetic effects are identifiable as well as the environmental effects. When the environmental effects are treated as fixed and there are other fixed factors in the model, the contrasts involving environmental effects, the variance of environmental sensitivities (genetic slopes) and the residual variance are the only identifiable parameters. These different identifiability scenarios were generated by changing the formulation of the model and the structure of the data and the models were then implemented via MCMC. The output of MCMC sampling schemes was interpreted in the light of the theoretical findings. The erratic behaviour of the MCMC chains was shown to be associated with identifiability problems in the likelihood, despite propriety of posterior distributions, achieved by arbitrarily chosen uniform (bounded) priors. In some cases, very long chains were needed before the pattern of behaviour of the chain may signal the existence of problems. The paper serves as a warning concerning the implementation of complex models where identifiability problems can be difficult to detect a priori. We conclude that it would be good practice to experiment with a proposed model and to understand its features before embarking on a full MCMC implementation.

IntCal09 and Marine09 Radiocarbon Age Calibration Curves, 0–50,000 Years cal BP

The IntCal04 and Marine04 radiocarbon calibration curves have been updated from 12 cal kBP (cal kBP is here defined as thousands of calibrated years before AD 1950), and extended to 50 cal kBP, utilizing newly available data sets that meet the IntCal Working Group criteria for pristine corals and other carbonates and for quantification of uncertainty in both the 14C and calendar timescales as established in 2002. No change was made to the curves from 0–12 cal kBP. The curves were constructed using a Markov chain Monte Carlo (MCMC) implementation of the random walk model used for IntCal04 and Marine04. The new curves were ratified at the 20th International Radiocarbon Conference in June 2009 and are available in the Supplemental Material at www.radiocarbon.org.