Bootstrap confidence bands for the CDF using ranked-set sampling

In ranked-set sampling (RSS), a stratification by ranks is used to obtain a sample that tends to be more informative than a simple random sample of the same size. Previous work has shown that if the rankings are perfect, then one can use RSS to obtain Kolmogorov–Smirnov type confidence bands for the CDF that are narrower than those obtained under simple random sampling. Here we develop Kolmogorov–Smirnov type confidence bands that work well whether the rankings are perfect or not. These confidence bands are obtained by using a smoothed bootstrap procedure that takes advantage of special features of RSS. We show through a simulation study that the coverage probabilities are close to nominal even for samples with just two or three observations. A new algorithm allows us to avoid the bootstrap simulation step when sample sizes are relatively small.     

Nonparametric particle filtering and smoothing with quasi-Monte Carlo sampling

Sequential Monte Carlo methods (also known as particle filters and smoothers) are used for filtering and smoothing in general state-space models. These methods are based on importance sampling. In practice, it is often difficult to find a suitable proposal which allows effective importance sampling. This article develops an original particle filter and an original particle smoother which employ nonparametric importance sampling. The basic idea is to use a nonparametric estimate of the marginally optimal proposal. The proposed algorithms provide a better approximation of the filtering and smoothing distributions than standard methods. The methods’ advantage is most distinct in severely nonlinear situations. In contrast to most existing methods, they allow the use of quasi-Monte Carlo (QMC) sampling. In addition, they do not suffer from weight degeneration rendering a resampling step unnecessary. For the estimation of model parameters, an efficient on-line maximum-likelihood (ML) estimation technique is proposed which is also based on nonparametric approximations. All suggested algorithms have almost linear complexity for low-dimensional state-spaces. This is an advantage over standard smoothing and ML procedures. Particularly, all existing sequential Monte Carlo methods that incorporate QMC sampling have quadratic complexity. As an application, stochastic volatility estimation for high-frequency financial data is considered, which is of great importance in practice. The computer code is partly available as supplemental material.     

Perfect simulation of positive Gaussian distributions

We provide an exact simulation algorithm that produces variables from truncated Gaussian distributions on ( ₊) ^p via a perfect sampling scheme, based on stochastic ordering and slice sampling, since accept-reject algorithms like the one of Geweke (1991) and Robert (1995) are difficult to extend to higher dimensions.     

Free energy methods for Bayesian inference: efficient exploration of univariate Gaussian mixture posteriors

Because of their multimodality, mixture posterior distributions are difficult to sample with standard Markov chain Monte Carlo (MCMC) methods. We propose a strategy to enhance the sampling of MCMC in this context, using a biasing procedure which originates from computational Statistical Physics. The principle is first to choose a “reaction coordinate”, that is, a “direction” in which the target distribution is multimodal. In a second step, the marginal log-density of the reaction coordinate with respect to the posterior distribution is estimated; minus this quantity is called “free energy” in the computational Statistical Physics literature. To this end, we use adaptive biasing Markov chain algorithms which adapt their targeted invariant distribution on the fly, in order to overcome sampling barriers along the chosen reaction coordinate. Finally, we perform an importance sampling step in order to remove the bias and recover the true posterior. The efficiency factor of the importance sampling step can easily be estimated a priori once the bias is known, and appears to be rather large for the test cases we considered.     

Reversible jump and the label switching problem in hidden Markov models

Reversible jump Markov chain Monte Carlo (RJMCMC) algorithms can be efficiently applied in Bayesian inference for hidden Markov models (HMMs), when the number of latent regimes is unknown. As for finite mixture models, when priors are invariant to the relabelling of the regimes, HMMs are unidentifiable in data fitting, because multiple ways to label the regimes can alternate during the MCMC iterations; this is the so-called label switching problem. HMMs with an unknown number of regimes are considered here and the goal of this paper is the comparison, both applied and theoretical, of five methods used for tackling label switching within a RJMCMC algorithm; they are: post-processing, partial reordering, permutation sampling, sampling from a Markov prior and rejection sampling. The five strategies we compare have been proposed mostly in the literature of finite mixture models and only two of them, i.e. rejection sampling and partial reordering, have been presented in RJMCMC algorithms for HMMs. We consider RJMCMC algorithms in which the parameters are updated by Gibbs sampling and the dimension of the model changes in split-and-merge and birth-and-death moves. Finally, an example illustrates and compares the five different methodologies.     

A Jonckheere-Terpstra-type test for perfect ranking in balanced ranked set sampling

Many methods based on ranked set sampling (RSS) assume perfect ranking of the samples. Here, by using the data measured by a balanced RSS scheme, we propose a nonparametric test for the assumption of perfect ranking. The test statistic that we use formally corresponds to the Jonckheere-Terpstra-type test statistic. We show formal relations of the proposed test for perfect ranking to other methods proposed recently in the literature. Through an empirical power study, we demonstrate that the proposed method performs favorably compared to many of its competitors.     

Rate estimation in partially observed Markov jump processes with measurement errors

We present a simulation methodology for Bayesian estimation of rate parameters in Markov jump processes arising for example in stochastic kinetic models. To handle the problem of missing components and measurement errors in observed data, we embed the Markov jump process into the framework of a general state space model. We do not use diffusion approximations. Markov chain Monte Carlo and particle filter type algorithms are introduced which allow sampling from the posterior distribution of the rate parameters and the Markov jump process also in data-poor scenarios. The algorithms are illustrated by applying them to rate estimation in a model for prokaryotic auto-regulation and the stochastic Oregonator, respectively.     

Inference in molecular population genetics

Full likelihood-based inference for modern population genetics data presents methodological and computational challenges. The problem is of considerable practical importance and has attracted recent attention, with the development of algorithms based on importance sampling (IS) and Markov chain Monte Carlo (MCMC) sampling. Here we introduce a new IS algorithm. The optimal proposal distribution for these problems can be characterized, and we exploit a detailed analysis of genealogical processes to develop a practicable approximation to it. We compare the new method with existing algorithms on a variety of genetic examples. Our approach substantially outperforms existing IS algorithms, with efficiency typically improved by several orders of magnitude. The new method also compares favourably with existing MCMC methods in some problems, and less favourably in others, suggesting that both IS and MCMC methods have a continuing role to play in this area. We offer insights into the relative advantages of each approach, and we discuss diagnostics in the IS framework.     

Efficient heterogeneous sampling for stochastic simulation with an illustration in health care applications

In modeling disease transmission, contacts are assumed to have different infection rates. A proper simulation must model the heterogeneity in the transmission rates. In this article, we present a computationally efficient algorithm that can be applied to a population with heterogeneous transmission rates. We conducted a simulation study to show that the algorithm is more efficient than other algorithms for sampling the disease transmission in a subset of the heterogeneous population. We use a valid stochastic model of pandemic influenza to illustrate the algorithm and to estimate the overall infection attack rates of influenza A (H1N1) in a Canadian city.     

Maximizing generalized linear mixed model likelihoods with an automated Monte Carlo EM algorithm

Two new implementations of the EM algorithm are proposed for maximum likelihood fitting of generalized linear mixed models. Both methods use random (independent and identically distributed) sampling to construct Monte Carlo approximations at the E-step. One approach involves generating random samples from the exact conditional distribution of the random effects (given the data) by rejection sampling, using the marginal distribution as a candidate. The second method uses a multivariate t importance sampling approximation. In many applications the two methods are complementary. Rejection sampling is more efficient when sample sizes are small, whereas importance sampling is better with larger sample sizes. Monte Carlo approximation using random samples allows the Monte Carlo error at each iteration to be assessed by using standard central limit theory combined with Taylor series methods. Specifically, we construct a sandwich variance estimate for the maximizer at each approximate E-step. This suggests a rule for automatically increasing the Monte Carlo sample size after iterations in which the true EM step is swamped by Monte Carlo error. In contrast, techniques for assessing Monte Carlo error have not been developed for use with alternative implementations of Monte Carlo EM algorithms utilizing Markov chain Monte Carlo E-step approximations. Three different data sets, including the infamous salamander data of McCullagh and Nelder, are used to illustrate the techniques and to compare them with the alternatives. The results show that the methods proposed can be considerably more efficient than those based on Markov chain Monte Carlo algorithms. However, the methods proposed may break down when the intractable integrals in the likelihood function are of high dimension.     

Recursive computation of inclusion probabilities in ranked-set sampling

We derive recursive algorithms for computing first-order and second-order inclusion probabilities for ranked-set sampling from a finite population. These algorithms make it practical to compute inclusion probabilities even for relatively large sample and population sizes. As an application, we use the inclusion probabilities to examine the performance of Horvitz-Thompson estimators under different varieties of balanced ranked-set sampling. We find that it is only for balanced Level 2 sampling that the Horvitz-Thompson estimator can be relied upon to outperform the simple random sampling mean estimator.     

Paired double-ranked set sampling

Abstract

In environmental monitoring and assessment, the main focus is to achieve observational economy and to collect data with unbiased, efficient and cost-effective sampling methods. Ranked set sampling (RSS) is one traditional method that is mostly used for accomplishing observational economy. In this article, we propose an unbiased sampling scheme, named paired double RSS (PDRSS) for estimating the population mean. We study the performance of the mean estimators under PDRSS based on perfect and imperfect rankings. It is shown that, for perfect ranking, the variance of the mean estimator under PDRSS is always less than the variance of mean estimator based on simple random sampling, paired RSS and RSS. The mean estimators under RSS, median RSS, PDRSS, and double RSS are also compared with the regression estimator of population mean based on SRS. The procedure is also illustrated with a case study using a real data set.     

Best linear unbiased and invariant estimation in location-scale families based on double-ranked set sampling

Abstract

In this article, we propose the best linear unbiased estimators (BLUEs) and best linear invariant estimators (BLIEs) for the unknown parameters of location-scale family of distributions based on double-ranked set sampling (DRSS) using perfect and imperfect rankings. These estimators are then compared with the BLUEs and BLIEs based on ranked set sampling (RSS). It is shown that under perfect ranking, the proposed estimators are uniformly better than the BLUEs and BLIEs obtained via RSS. We also propose the best linear unbiased quantile (BLUQ) and the best linear invariant quantile (BLIQ) estimators for normal distribution under DRSS. It is observed that the proposed quantile estimators are more efficient than the BLUQ and BLIQ estimators based on RSS for both perfect and imperfect orderings.     

Convergence in the Wasserstein Metric for Markov Chain Monte Carlo Algorithms with Applications to Image Restoration

Abstract

In this paper, we show how the time for convergence to stationarity of a Markov chain can be assessed using the Wasserstein metric, rather than the usual choice of total variation distance. The Wasserstein metric may be more easily applied in some applications, particularly those on continuous state spaces. Bounds on convergence time are established by considering the number of iterations required to approximately couple two realizations of the Markov chain to within ε tolerance. The particular application considered is the use of the Gibbs sampler in the Bayesian restoration of a degraded image, with pixels that are a continuous grey-scale and with pixels that can only take two colours. On finite state spaces, a bound in the Wasserstein metric can be used to find a bound in total variation distance. We use this relationship to get a precise O(N log N) bound on the convergence time of the stochastic Ising model that holds for appropriate values of its parameter as well as other binary image models. Our method employing convergence in the Wasserstein metric can also be applied to perfect sampling algorithms involving coupling from the past to obtain estimates of their running times.     

Model-based adaptive spatial sampling for occurrence map construction

In many environmental management problems, the construction of occurrence maps of species of interest is a prerequisite to their effective management. However, the construction of occurrence maps is a challenging problem because observations are often costly to obtain (thus incomplete) and noisy (thus imperfect). It is therefore critical to develop tools for designing efficient spatial sampling strategies and for addressing data uncertainty. Adaptive sampling strategies are known to be more efficient than non-adaptive strategies. Here, we develop a model-based adaptive spatial sampling method for the construction of occurrence maps. We apply the method to estimate the occurrence of one of the world’s worst invasive species, the red imported fire ant, in and around the city of Brisbane, Australia. Our contribution is threefold: (i) a model of uncertainty about invasion maps using the classical image analysis probabilistic framework of Hidden Markov Random Fields (HMRF), (ii) an original exact method for optimal spatial sampling with HMRF and approximate solution algorithms for this problem, both in the static and adaptive sampling cases, (iii) an empirical evaluation of these methods on simulated problems inspired by the fire ants case study. Our analysis demonstrates that the adaptive strategy can lead to substantial improvement in occurrence mapping.     

An omnibus two-sample test for ranked-set sampling data

We develop an omnibus two-sample test for ranked-set sampling (RSS) data. The test statistic is the conditional probability of seeing the observed sequence of ranks in the combined sample, given the observed sequences within the separate samples. We compare the test to existing tests under perfect rankings, finding that it can outperform existing tests in terms of power, particularly when the set size is large. The test does not maintain its level under imperfect rankings. However, one can create a permutation version of the test that is comparable in power to the basic test under perfect rankings and also maintains its level under imperfect rankings. Both tests extend naturally to judgment post-stratification, unbalanced RSS, and even RSS with multiple set sizes. Interestingly, the tests have no simple random sampling analog.     

Some results concerning off-training-set and IID error for the Gibbs and the Bayes optimal generalizers

In this paper we analyse the average behaviour of the Bayes-optimal and Gibbs learning algorithms. We do this both for off-training-set error and conventional IID (independent identically distributed) error (for which test sets overlap with training sets). For the IID case we provide a major extension to one of the better known results. We also show that expected IID test set error is a non-increasing function of training set size for either algorithm. On the other hand, as we show, the expected off-training-set error for both learning algorithms can increase with training set size, for non-uniform sampling distributions. We characterize the relationship the sampling distribution must have with the prior for such an increase. We show in particular that for uniform sampling distributions and either algorithm, the expected off-training-set error is a non-increasing function of training set size. For uniform sampling distributions, we also characterize the priors for which the expected error of the Bayes-optimal algorithm stays constant. In addition we show that for the Bayes-optimal algorithm, expected off-training-set error can increase with training set size when the target function is fixed, but if and only if the expected error averaged over all targets decreases with training set size. Our results hold for arbitrary noise and arbitrary loss functions.     

Parametric models for response-biased sampling

Suppose that subjects in a population follow the model f   ( y ^* x ^*; ) where y ^* denotes a response, x ^* denotes a vector of covariates and is the parameter to be estimated. We consider response-biased sampling, in which a subject is observed with a probability which is a function of its response. Such response-biased sampling frequently occurs in econometrics, epidemiology and survey sampling. The semiparametric maximum likelihood estimate of is derived, along with its asymptotic normality, efficiency and variance estimates. The estimate proposed can be used as a maximum partial likelihood estimate in stratified response-selective sampling. Some computation algorithms are also provided.     

Sample-based Maximum Likelihood Estimation of the Autologistic Model

New recursive algorithms for fast computation of the normalizing constant for the autologistic model on the lattice make feasible a sample-based maximum likelihood estimation (MLE) of the autologistic parameters. We demonstrate by sampling from 12 simulated 420×420 binary lattices with square lattice plots of size 4×4, …, 7×7 and sample sizes between 20 and 600. Sample-based results are compared with ‘benchmark’ MCMC estimates derived from all binary observations on a lattice. Sample-based estimates are, on average, biased systematically by 3%–7%, a bias that can be reduced by more than half by a set of calibrating equations. MLE estimates of sampling variances are large and usually conservative. The variance of the parameter of spatial association is about 2–10 times higher than the variance of the parameter of abundance. Sample distributions of estimates were mostly non-normal. We conclude that sample-based MLE estimation of the autologistic parameters with an appropriate sample size and post-estimation calibration will furnish fully acceptable estimates. Equations for predicting the expected sampling variance are given.     

On MCMC sampling in hierarchical longitudinal models

Markov chain Monte Carlo (MCMC) algorithms have revolutionized Bayesian practice. In their simplest form (i.e., when parameters are updated one at a time) they are, however, often slow to converge when applied to high-dimensional statistical models. A remedy for this problem is to block the parameters into groups, which are then updated simultaneously using either a Gibbs or Metropolis-Hastings step. In this paper we construct several (partially and fully blocked) MCMC algorithms for minimizing the autocorrelation in MCMC samples arising from important classes of longitudinal data models. We exploit an identity used by Chib (1995) in the context of Bayes factor computation to show how the parameters in a general linear mixed model may be updated in a single block, improving convergence and producing essentially independent draws from the posterior of the parameters of interest. We also investigate the value of blocking in non-Gaussian mixed models, as well as in a class of binary response data longitudinal models. We illustrate the approaches in detail with three real-data examples.