In this paper, we discuss an extension of the Split Hamiltonian Monte Carlo
(Split HMC) method for Gaussian process model (GPM). This method is based on
splitting the Hamiltonian in a way that allows much of the movement around the
state space to be done at low computational cost. To this end, we approximate
the negative log density (i.e., the energy function) of the distribution of
interest by a quadratic function U0 for which Hamiltonian dynamics can be
solved analytically. The overall energy function U is then written as U0 + U1,
where U1 is the approximation error.

Regularization is widely used in statistics and machine learning to prevent
overfitting and gear solution towards prior information. In general, a
regularized estimation problem minimizes the sum of a loss function and a
penalty term. The penalty term is usually weighted by a tuning parameter and
encourages certain constraints on the parameters to be estimated. Particular
choices of constraints lead to the popular lasso, fused-lasso, and other
generalized $l_1$ penalized regression methods.

In this paper we describe the basic features of the Bergm package for the
open-source R software which provides a comprehensive framework for Bayesian
analysis for exponential random graph models: tools for parameter estimation,
model selection and goodness-of-fit diagnostics. We illustrate the capabilities
of this package describing the algorithms that drive the package through a
tutorial analysis of two well-known network datasets.

Exponential random graph models are a class of widely used exponential family
models for social networks. The topological structure of an observed network is
modeled by the relative prevalence of a set of local sub-graph configurations
termed network statistics. One of the key tasks in the application of these
models is which network statistics to include in the model. This can be thought
of as statistical model selection problem.

This paper addresses the problem of summarizing the posterior distributions
that typically arise, in a Bayesian framework, when dealing with signal
decomposition problems with unknown number of components. Such posterior
distributions are defined over union of subspaces of differing dimensionality
and can be sampled from using modern Monte Carlo techniques, for instance the
increasingly popular RJ-MCMC method. No generic approach is available, however,
to summarize the resulting variable-dimensional samples and extract from them
component-specific parameters.

Sequential Monte Carlo (SMC) methods, also known as particle filters, are
simulation-based recursive algorithms for the approximation of the a posteriori
probability measures generated by state-space dynamical models. At any given
(discrete) time t, a particle filter produces a set of samples over the state
space of the system of interest. These samples are referred to as "particles"
and can be used to build a discrete approximation of the a posteriori
probability distribution of the state, conditional on a sequence of available
observations.

Hamiltonian Monte Carlo (HMC) is a Markov chain Monte Carlo (MCMC) algorithm
that avoids the random walk behavior and sensitivity to correlated parameters
that plague many MCMC methods by taking a series of steps informed by
first-order gradient information. These features allow it to converge to
high-dimensional target distributions much more quickly than simpler methods
such as random walk Metropolis or Gibbs sampling. However, HMC's performance is
highly sensitive to two user-specified parameters: a step size {\epsilon} and a
desired number of steps L.

Bundle adjustment (BA) is the problem of refining a visual reconstruction to
produce better structure and viewing parameter estimates. This problem is often
formulated as a nonlinear least squares problem, where data arises from
interest point matching. Mismatched interest points cause serious problems in
this approach, as a single mismatch will affect the entire reconstruction. In
this paper, we propose a novel robust Student's t BA algorithm (RST-BA).

We extend Tropp's analysis of Orthogonal Matching Pursuit (OMP) using the
Exact Recovery Condition (ERC) to a first exact recovery analysis of Orthogonal
Least Squares (OLS). We show that when ERC is met, OLS is guaranteed to exactly
recover the unknown support. Moreover, we provide a closer look at the analysis
of both OMP and OLS when ERC is not fulfilled. We show that there exist
dictionaries for which some subsets are never recovered with OMP. This
phenomenon, which also appears with $\ell_1$ minimization, does not occur for
OLS.

Methods for choosing a fixed set of knot locations in additive spline models
are fairly well established in the statistical literature. While most of these
methods are in principle directly extendable to non-additive surface models,
they are likely to be less successful in that setting because of the curse of
dimensionality, especially when there are more than a couple of covariates.

There are two main algorithms for fitting constrained marginal models to
discrete data available in the literature. We show that they are equivalent,
each being advantageous in different circumstances, and give some results on
their convergence properties. Extensions of the method are provided to allow
the inclusion of individual covariates, and for maximization under
$L_1$-penalties.

Let $K$ be a closed convex polyhedron defined by a finite number of linear
inequalities. In this paper we refine the theory of abstract tubes (Naiman and
Wynn, 1997) associated with $K$ when $K$ is perturbed. In particular, we focus
on the perturbation that is lexicographic and in an outer direction. An
algorithm for constructing the abstract tube by means of linear programming and
its implementation are discussed.

The problem of rank-one approximation of symmetric tensors is important in
independent components analysis, also known as blind source separation, as well
as polynomial optimization. The rank-one approximation problem reduces to
finding the principal Z eigenvector, and an effective algorithm was recently
discovered for computing symmetric tensor eigenvectors.

We present an iterative sampling method which delivers upper and lower
bounding processes for the Brownian path. We develop such processes with
particular emphasis on being able to unbiasedly simulate them on a personal
computer. The dominating processes converge almost surely in the supremum and
$L_1$ norms.

Statistical analysis of max-stable processes used to model spatial extremes
has been limited by the difficulty in calculating the joint likelihood
function. This precludes all standard likelihood-based approaches, including
Bayesian approaches. In this paper we present a Bayesian approach through the
use of approximate Bayesian computing. This circumvents the need for a joint
likelihood function by instead relying on simulations from the (unavailable)
likelihood. This method is compared with an alternative approach based on the
composite likelihood.

While statisticians are well-accustomed to performing exploratory analysis in
the modeling stage of an analysis, the notion of conducting preliminary
general-purpose exploratory analysis in the Monte Carlo stage (or more
generally, the model-fitting stage) of an analysis is an area which we feel
deserves much further attention. Towards this aim, this paper proposes a
general-purpose algorithm for automatic density exploration.

We present a new approach to Bayesian inference that entirely avoids Markov
chain simulation, by constructing a map that pushes forward the prior measure
to the posterior measure. Existence and uniqueness of a suitable
measure-preserving map is established by formulating the problem in the context
of optimal transport theory. We discuss various means of explicitly
parameterizing the map and computing it efficiently through solution of an
optimization problem, exploiting gradient information from the forward model
when possible.

Accuracy and Availability computations for a GNSS System - or combination of
Systems - through Service Volume Simulations take considerable time. Therefore,
the computation of the accuracy in 2D and 3D are often simplified by an
approximate solution. The drawback is that such simplifications can lead to
accuracy results that are too conservative (up to 25% in the 2D case and up to
43% in the 3D case, for a 95% confidence level), which in turn translates into
pessimistic System Availability.

In this paper, we develop a new method to sample from posterior distributions
in hierarchical models without using Markov chain Monte Carlo. This method is
generally applicable to high-dimensional models involving large data sets.
Illustrative analysis exemplifies the ease with which one could implement our
method, which results in independent samples from the posterior distributions
of interest.

We address the issue of knots selection for Gaussian predictive process
methodology. Predictive process approximation provides an effective solution to
the cubic order computational complexity of Gaussian process models. This
approximation crucially depends on a set of points, called knots, at which the
original process is retained, while the rest is approximated via a
deterministic extrapolation. Knots should be few in number to keep the
computational complexity low, but provide a good coverage of the process domain
to limit approximation error.

The correlation of the result lists provided by search engines is fundamental
and it has deep and multidisciplinary ramifications. Here, we present automatic
and unsupervised methods to assess whether or not search engines provide
results that are comparable or correlated. We have two main contributions:
First, we provide evidence that for more than 80% of the input queries -
independently of their frequency - the two major search engines share only
three or fewer URLs in their search results, leading to an increasing
divergence.

The paper builds upon a recent approach to find the approximate bounds of a
real function using Polynomial Chaos expansions. Given a function of random
variables with compact support probability distributions, the intuition is to
quantify the uncertainty in the response using Polynomial Chaos expansion and
discard all the information provided about the randomness of the output and
extract only the bounds of its compact support.

computation (ABC). Here we discuss how to construct the perturbation kernels
that are required in ABC-SMC approaches, in order to construct a set of
distributions that start out from a suitably defined prior and converge towards
the unknown posterior. We derive optimality criteria for different kernels,
which are based on the Kullback-Leibler divergence between intermediate
distributions and the posterior distribution. We will show that for many
complicated posterior distributions locally adapted kernels tend to show the
best performance.

In this article we propose a novel MCMC method based on deterministic
transformations T : X x D --> X where X is the state-space and D is some set
which may or may not be a subset of X. We refer to our new methodology as
Transformation-based Markov chain Monte Carlo (TMCMC).

The entropy is one of the most applicable uncertainty measures in many
statistical and en- gineering problems. In statistical literature, the entropy
is used in calculation of the Kullback- Leibler (KL) information which is a
powerful mean for performing goodness of fit tests. Ranked Set Sampling (RSS)
seems to provide improved estimators of many parameters of the popu- lation in
the huge studied problems in the literature. It is developed for situations
where the variable of interest is difficult or expensive to measure, but where
ranking in small sub-samples is easy.

Reversible jump Markov chain Monte Carlo (RJMCMC) extends ordinary MCMC
methods for use in Bayesian multimodel inference. We show that RJMCMC can be
implemented as Gibbs sampling with alternating updates of a model indicator and
a vector-valued "palette" of parameters denoted $\bm \psi$. Like an artist uses
the palette to mix dabs of color for specific needs, we create model-specific
parameters from the set available in $\bm \psi$.

A new methodology for model determination in decomposable graphical Gaussian
models is developed. The Bayesian paradigm is used and, for each given graph, a
hyper inverse Wishart prior distribution on the covariance matrix is
considered. This prior distribution depends on hyper-parameters. It is
well-known that the models's posterior distribution is sensitive to the
specification of these hyper-parameters and no completely satisfactory method
is registered.

We recast importance sampling in terms of minimizing a Monte Carlo estimate
of the Kullback-Leibler divergence over a well-chosen exponential family. We
further examine two alternative Kullback-Leibler approximations that do not
require the target distribution to be normalized. One of these makes the
variational problem equivalent to normalized importance sampling, while the
other leads to a new moment estimation technique that works remarkably well if
the target distribution is reasonably close to the approximating family even
though the proposal distribution is mismatched.

Approximate Bayesian computation (ABC) is a class of algorithmic methods in
Bayesian inference using statistical summaries and computer simulations. ABC
has become popular in evolutionary genetics and in other branches of biology.
However model selection under ABC algorithms has been a subject of intense
debate during the recent years. Here we propose novel approaches to model
selection based on posterior predictive distributions and approximations of the
deviance.

Importance sampling is a common technique for Monte Carlo approximation,
including Monte Carlo approximation of p-values. Here it is shown that a simple
correction of the usual importance sampling p-values creates valid p-values,
meaning that a hypothesis test created by rejecting the null when the p-value
is <= alpha will also have a type I error rate <= alpha. This correction uses
the importance weight of the original observation, which gives valuable
diagnostic information under the null hypothesis.

We present a computational approach for generating Markov bases for multi-way
contin- gency tables whose cells counts might be constrained by fixed marginals
and by lower and upper bounds. Our framework includes tables with structural
zeros as a particular case. In- stead of computing the entire Markov basis in
an initial step, our framework finds sets of local moves that connect each
table in the reference set with a set of neighbor tables. We construct a Markov
chain on the reference set of tables that requires only a set of local moves at
each iteration.

This paper addresses finite sample stability properties of sequential Monte
Carlo methods for approximating sequences of probability distributions. The
results presented herein are applicable in the scenario where the start and end
distributions in the sequence are fixed and the number of intermediate steps is
a parameter of the algorithm. Under assumptions which hold on non-compact
spaces, it is shown that the effect of the initial distribution decays
exponentially fast in the number of intermediate steps and the corresponding
stochastic error is stable in \mathbb{L}_{p} norm.

Many least squares problems involve affine equality and inequality
constraints. Although there are variety of methods for solving such problems,
most statisticians find constrained estimation challenging. The current paper
proposes a new path following algorithm for quadratic programming based on
exact penalization. Similar penalties arise in $l_1$ regularization in model
selection. Classical penalty methods solve a sequence of unconstrained problems
that put greater and greater stress on meeting the constraints.

An off-line handwritten alphabetical character recognition system using
multilayer feed forward neural network is described in the paper. A new method,
called, diagonal based feature extraction is introduced for extracting the
features of the handwritten alphabets. Fifty data sets, each containing 26
alphabets written by various people, are used for training the neural network
and 570 different handwritten alphabetical characters are used for testing.

Global sensitivity analysis is often impracticable for complex and time
demanding numerical models, as it requires a large number of runs. The
reduced-basis approach provides a way to replace the original model by a much
faster to run code. In this paper, we are interested in the information loss
induced by the approximation on the estimation of sensitivity indices. We
present a method to provide a robust error assessment, hence enabling
significant time savings without sacrifice on precision and rigourousness.

Online (also called "recursive" or "adaptive") estimation of fixed model
parameters in hidden Markov models is a topic of much interest in times series
modelling. In this work, we propose an online parameter estimation algorithm
that combines two key ideas. The first one, which is deeply rooted in the
Expectation-Maximization (EM) methodology consists in reparameterizing the
problem using complete-data sufficient statistics. The second ingredient
consists in exploiting a purely recursive form of smoothing in HMMs based on an
auxiliary recursion.

We consider various versions of adaptive Gibbs and Metropolis- within-Gibbs
samplers, which update their selection probabilities (and per- haps also their
proposal distributions) on the y during a run, by learning as they go in an
attempt to optimise the algorithm.We present a cautionary example of how even a
simple-seeming adaptive Gibbs sampler may fail to converge.We then present
various positive results guaranteeing convergence of adaptive Gibbs samplers
under certain conditions.

Approximate Bayesian computation (ABC), also known as likelihood-free
methods, have become a favourite tool for the analysis of complex stochastic
models, primarily in population genetics but also in financial analyses. We
advocated in Grelaud et al.

Clustering with fast algorithms large samples of high dimensional data is an
important challenge in computational statistics. Borrowing ideas from MacQueen
(1967) who introduced a sequential version of the $k$-means algorithm, we
propose in this paper a new class of recursive stochastic gradient algorithms
designed for the $k$-medians loss criterion. By their recursive nature, these
algorithms are very fast and are well adapted to deal with large samples of
data that are allowed to arrive sequentially.

This paper describes a method for estimating the marginal likelihood or Bayes
factors of Bayesian models using non-parametric importance sampling ("arrogance
sampling"). This method can also be used to compute the normalizing constant of
probability distributions. Because the required inputs are samples from the
distribution to be normalized and the scaled density at those samples, this
method may be a convenient replacement for the harmonic mean estimator. The
method has been implemented in the open source R package margLikArrogance.

The Adaptive Multiple Importance Sampling (AMIS) algorithm is aimed at an
optimal recycling of past simulations in an iterated importance sampling
scheme. The difference with earlier adaptive importance sampling
implementations like Population Monte Carlo is that the importance weights of
all simulated values, past as well as present, are recomputed at each
iteration, following the technique of the deterministic multiple mixture
estimator of Owen and Zhou (2000).

A general purpose variance reduction technique for Markov chain Monte Carlo
estimators based on the zero-variance principle introduced in the physics
literature by Assaraf and Caffarel (1999, 2003), is proposed. Conditions for
unbiasedness of the zero-variance estimator are derived. A central limit
theorem is also proved under regularity conditions. The potential of the new
idea is illustrated with real applications to Bayesian inference for probit,
logit and GARCH models.

L1 -penalized regression methods such as the Lasso (Tibshirani 1996) that
achieve both variable selection and shrinkage have been very popular. An
extension of this method is the Fused Lasso (Tibshirani and Wang 2007), which
allows for the incorporation of external information into the model. In this
article, we develop new and fast algorithms for solving the Fused Lasso which
are based on coordinate-wise optimization. This class of algorithms has
recently been applied very successfully to solve L1 -penalized problems very
quickly (Friedman et al. 2007).

This paper proposes approaches for the analysis of multiple changepoint
models when dependency in the data is modelled through a hierarchical Gaussian
Markov random field model. Integrated nested Laplace approximations are used to
approximate data quantities, and an approximate filtering recursions approach
is proposed for savings in compuational cost when detecting changepoints. All
of these methods are simulation free. Analysis of real data demonstrates the
usefulness of the approach in general.

In this paper, we study the alternating direction method for finding the
Dantzig selectors, which are first introduced in [8]. In particular, at each
iteration we apply the nonmonotone gradient method proposed in [17] to
approximately solve one subproblem of this method. We compare our approach with
a first-order method proposed in [3]. The computational results show that our
approach usually outperforms that method in terms of CPU time while producing
solutions of comparable quality.

This paper presents a Markov chain Monte Carlo method to generate approximate
posterior samples in retrospective multiple changepoint problems where the
number of changes is not known in advance. The method uses conjugate models
whereby the marginal likelihood for the data between consecutive changepoints
is tractable. Inclusion of hyperpriors gives a near automatic algorithm
providing a robust alternative to popular filtering recursions approaches in
cases which may be sensitive to prior information. Three real examples are used
to demonstrate the proposed approach.

The EM algorithm is a widely used methodology for penalized likelihood
estimation. Provable monotonicity and convergence are the hallmarks of the EM
algorithm and these properties are well established for smooth likelihood and
smooth penalty functions. However, many relaxed versions of variable selection
penalties are not smooth. The goal of this paper is to introduce a new class of
Space Alternating Penalized Kullback Proximal extensions of the EM algorithm
for nonsmooth likelihood inference.

We present a very fast algorithm for general matrix factorization of a data
matrix for use in the statistical analysis of high-dimensional data via latent
factors. Such data are prevalent across many application areas and generate an
ever-increasing demand for methods of dimension reduction in order to undertake
the statistical analysis of interest. Our algorithm uses a gradient-based
approach which can be used with an arbitrary loss function provided the latter
is differentiable.

This technical report is the union of two contributions to the discussion of
the Read Paper "Riemann manifold Langevin and Hamiltonian Monte Carlo methods"
by B. Calderhead and M. Girolami, presented in front of the Royal Statistical
Society on October 13th 2010 and to appear in the Journal of the Royal
Statistical Society Series B. The first comment establishes a parallel and
possible interactions with Adaptive Monte Carlo methods.

In Fern\'andez-Dur\'an (2004), a new family of circular distributions based
on nonnegative trigonometric sums (NNTS models) is developed. Because the
parameter space of this family is the surface of the hypersphere, an efficient
Newton-like algorithm on manifolds is generated in order to obtain the maximum
likelihood estimates of the parameters.

We study quasi-Newton methods from the viewpoint of information geometry
induced associated with Bregman divergences. Fletcher has studied a variational
problem which derives the approximate Hessian update formula of the
quasi-Newton methods. We point out that the variational problem is identical to
optimization of the Kullback-Leibler divergence, which is a discrepancy measure
between two probability distributions. The Kullback-Leibler divergence for the
multinomial normal distribution corresponds to the objective function Fletcher
has considered.

In this paper, we consider the implications of the fact that parallel
raw-power can be exploited by a generic Metropolis--Hastings algorithm if the
proposed values are independent. In particular, we present improvements to the
independent Metropolis--Hastings algorithm that significantly decrease the
variance of any estimator derived from the MCMC output, for a null computing
cost since those improvements are based on a fixed number of target density
evaluations.

Stochastic differential equations (SDEs) are established tools to model
physical phenomena whose dynamics are affected by random noise. By estimating
parameters of an SDE intrinsic randomness of a system around its drift can be
identified and separated from the drift itself. When it is of interest to model
dynamics within a given population, i.e. to model simultaneously the
performance of several experiments or subjects, mixed-effects modelling allows
for the distinction of between and within experiment variability.

In the following paper we consider a simulation technique for stochastic
trees. One of the most important areas in computational genetics is the
calculation and subsequent maximization of the likelihood function associated
to such models. This typically consists of using importance sampling (IS) and
sequential Monte Carlo (SMC) techniques. The approach proceeds by simulating
the tree, backward in time from observed data, to a most recent common ancestor
(MRCA). However, in many cases, the computational time and variance of
estimators are often too high to make standard approaches useful.

The problem of matching unlabelled point sets using Bayesian inference is
considered. Two recently proposed models for the likelihood are compared, based
on the Procrustes size-and-shape and the full configuration. Bayesian inference
is carried out for matching point sets using Markov chain Monte Carlo
simulation. An improvement to the existing Procrustes algorithm is proposed
which improves convergence rates, using occasional large jumps in the burn-in
period.

In this paper we study interior point (IP) methods for solving optimal design
problems. In particular, we propose a primal IP method for solving the problems
with general convex optimality criteria and establish its global convergence.
In addition, we reformulate the problems with A-, D- and E-criterion into
linear or log-determinant semidefinite programs (SDPs) and apply standard
primal-dual IP solvers such as SDPT3 [21,25] to solve the resulting SDPs. We
also compare the IP methods with the widely used multiplicative algorithm
introduced by Silvey et al. [18].

The parameters of a discrete stationary Markov model are transition
probabilities between states. Traditionally, data consists in sequences of
observed states for a given number of individuals over the whole observation
period. In such a case, the estimation of transition probabilities is
straightforwardly made by counting one-step moves from a given state to
another. In many real-life problems, however, the inference is much more
difficult as state sequences are not fully observed, namely the state of each
individual is known only for some given values of the time variable $t$.

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.

An overview is provided for the problem of calculating Bayes factors (BF) for
model selection using methods based on variants of importance sampling and
variants of Laplace Approximation, specifically path sampling (PS). These
methods are compared on several simulated and two real life data set. It is
noted that there are unexpected instabilities in PS and we provide both a
theoretical explanation and suggestions for improvement. Estimates of the same
BF may differ widely across methods.

Most surrogate models for computer experiments are interpolators, and the
most common interpolator is a Gaussian process (GP) that deliberately omits a
small-scale (measurement) error term called the nugget. The explanation is that
computer experiments are, by definition, "deterministic", and so there is no
measurement error. We think this is too narrow a focus for a computer
experiment and a statistically inefficient way to model them.

rdering of regression or classification coefficients occurs in many
real-world applications. Fused Lasso exploits this ordering by explicitly
regularizing the differences between neighboring coefficients through an
$\ell_1$ norm regularizer. However, due to nonseparability and nonsmoothness of
the regularization term, solving the fused Lasso problem is computationally
demanding. Existing solvers can only deal with problems of small or medium
size, or a special case of the fused Lasso problem in which the predictor
matrix is identity matrix.

We introduce a new class of Sequential Monte Carlo (SMC) methods, which we
call free energy SMC. This class is inspired by free energy methods, which
originate from Physics, and where one samples from a biased distribution such
that a given function $\xi(\theta)$ of the state $\theta$ is forced to be
uniformly distributed over a given interval.

The Gaussian process (GP) is a popular way to specify dependencies between
random variables in a probabilistic model. In the Bayesian framework the
covariance structure can be specified using unknown hyperparameters.
Integrating over these hyperparameters considers different possible
explanations for the data when making predictions. This integration is often
performed using Markov chain Monte Carlo (MCMC) sampling. However, with
non-Gaussian observations standard hyperparameter sampling approaches require
careful tuning and may converge slowly.

The classic chi-squared statistic for testing goodness-of-fit has long been a
cornerstone of modern statistical practice. The statistic consists of a sum in
which each summand involves multiplying by the inverse of (i.e., dividing by)
the probability associated with the corresponding bin in the distribution being
tested for goodness-of-fit. This inversion typically precipitates rebinning to
uniformize the probabilities associated with the bins, in order to make the
test reasonably powerful. With the now widespread availability of computers,
there is no longer any need for this.

We consider the augmented Lagrangian method (ALM) as a solver for the fused
lasso signal approximator (FLSA) problem. The ALM is a dual method in which
squares of the constraint functions are added as penalties to the Lagrangian.
In order to apply this method to FLSA, two types of auxiliary variables are
introduced to transform the original unconstrained minimization problem into a
linearly constrained minimization problem. Each updating in this iterative
algorithm consists of just a simple one-dimensional convex programming problem,
with closed form solution in many cases.

We consider the problem of cardinality estimation in data stream
applications, focusing on two techniques that use pseudo-random variates to
form low-dimensional data sketches. We apply conventional statistical methods
to compare algorithms based on storing either selected order statistics or
random projections. We derive estimators of the cardinality in both cases and
show that the maximal-term estimator is recursively computable and has
exponentially decreasing error bounds.

In this paper, we describe two computational methods for calculating the
cumulative distribution function and the upper quantiles of the maximal
difference between a Brownian bridge and its concave majorant. The first method
has two different variants that are both based on a Monte Carlo approach,
whereas the second uses the Gaver-Stehfest (GS) algorithm for numerical
inversion of Laplace transform.

A general methodology is presented for the construction and effective use of
control variates for reversible MCMC samplers. The values of the coefficients
of the optimal linear combination of the control variates are computed, and
adaptive, consistent MCMC estimators are derived for these optimal
coefficients. All methodological and asymptotic arguments are rigorously
justified. Numerous MCMC simulation examples from Bayesian inference
applications demonstrate that the resulting variance reduction can be quite
dramatic.

This paper tackles the problem of image deconvolution with joint estimation
of PSF parameters and hyperparameters. Within a Bayesian framework, the
solution is inferred via a global a posteriori law for unknown parameters and
object. The estimate is chosen as the posterior mean, numerically calculated by
means of a Monte-Carlo Markov chain algorithm. The estimates are efficiently
computed in the Fourier domain and the effectiveness of the method is shown on
simulated examples.

We study the group test for DNA library screening based on probabilistic
approach. Group test is a method of detecting a few positive items from among a
large number of items, and has wide range of applications. In DNA library
screening, positive item corresponds to the clone having a specified DNA
segment, and it is necessary to identify and isolate the positive clones for
compiling the libraries. In the group test, a group of items, called pool, is
assayed in a lump in order to save the cost of testing, and positive items are
detected based on the observation from each pool.

In the paradigm of computer experiments, the choice of an experimental design
is an important issue. When no information is available about the black-box
function to be approximated, an exploratory design have to be used. In this
context, two dispersion criteria are usually considered: the minimax and the
maximin ones. In the case of a hypercube domain, a standard strategy consists
of taking the maximin design within the class of latin hypercube designs.
However, in a non hypercube context, it does not make sense to use the latin
hypercube strategy.

This paper focuses on utilizing two different Bayesian methods to deal with a
variety of toy problems which occur in data analysis. In particular we
implement the Variational Bayesian and Nested Sampling methods to tackle the
problems of polynomial selection and Gaussian Mixture Models, comparing the
algorithms in terms of processing speed and accuracy. In the problems tackled
here it is the Variational Bayesian algorithms which are the faster though both
results give similar results.

We describe two slice sampling methods for taking multivariate steps using
the crumb framework. These methods use the gradients at rejected proposals to
adapt to the local curvature of the log-density surface, a technique that can
produce much better proposals when parameters are highly correlated. We
evaluate our methods on four distributions and compare their performance to
that of a non-adaptive slice sampling method and a Metropolis method. The
adaptive methods perform favorably on low-dimensional target distributions with
highly-correlated parameters.

Because of their multimodality, mixture posterior densities 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 density is multimodal. In a second step, the marginal
log-density of the reaction coordinate is estimated; this quantity is called
"free energy" in the computational statistical physics literature.

Implementing Bayesian variable selection for linear Gaussian regression
models for analysing high dimensional data sets is of current interest in many
fields. In order to make such analysis operational, we propose a new sampling
algorithm based upon Evolutionary Monte Carlo and designed to work under the
"large p, small n" paradigm, thus making fully Bayesian multivariate analysis
feasible, for example, in genetics/genomics experiments. Two real data examples
in genomics are presented, demonstrating the performance of the algorithm in a
space of up to 10,000 covariates.

While Robert and Rousseau (2010) addressed the foundational aspects of
Bayesian analysis, the current chapter details its practical aspects through a
review of the computational methods available for approximating Bayesian
procedures. Recent innovations like Monte Carlo Markov chain, sequential Monte
Carlo methods and more recently Approximate Bayesian Computation techniques
have considerably increased the potential for Bayesian applications and they
have also opened new avenues for Bayesian inference, first and foremost
Bayesian model choice.

In many applications, such as physiology and finance, large time series data
bases are to be analyzed requiring the computation of linear, nonlinear and
other measures. Such measures have been developed and implemented in commercial
and freeware softwares rather selectively and independently. The Measures of
Analysis of Time Series ({\tt MATS}) {\tt MATLAB} toolkit is designed to handle
an arbitrary large set of scalar time series and compute a large variety of
measures on them, allowing for the specification of varying measure parameters
as well.

In this paper, we propose a general class of algorithms for optimizing an
extensive variety of nonsmoothly penalized objective functions that satisfy
certain regularity conditions. The proposed framework utilizes the
majorization-minimization (MM) algorithm as its core optimization engine. The
resulting algorithms rely on iterated soft-thresholding, implemented
componentwise, allowing for fast, stable updating that avoids the need for any
high-dimensional matrix inversion.

We consider various versions of adaptive Gibbs and Metropolis within-Gibbs
samplers, which update their selection probabilities (and perhaps also their
proposal distributions) on the fly during a run, by learning as they go in an
attempt to optimise the algorithm. We present a cautionary example of how even
a simple-seeming adaptive Gibbs sampler may fail to converge. We then present
various positive results guaranteeing convergence of adaptive Gibbs samplers
under certain conditions.

Bayesian phylogenetic methods are generating noticeable enthusiasm in the
field of molecular systematics. Several phylogenetic models are often at stake
and different approaches are used to compare them within a Bayesian framework.
The Bayes factor, defined as the ratio of the marginal likelihoods of two
competing models, plays a key role in Bayesian model selection. However, its
computation is still a challenging problem.

Many probabilistic models introduce strong dependencies between variables
using a latent multivariate Gaussian distribution or a Gaussian process. We
present a new Markov chain Monte Carlo algorithm for performing inference in
models with multivariate Gaussian priors. Its key properties are: 1) it has
simple, generic code applicable to many models, 2) it has no free parameters,
3) it works well for a variety of Gaussian process based models.

$\alpha$-stable distributions are utilised as models for heavy-tailed noise
in many areas of statistics, finance and signal processing engineering.

However, in general, neither univariate nor multivariate $\alpha$-stable
models admit closed form densities which can be evaluated pointwise. This
complicates the inferential procedure.

As a result, $\alpha$-stable models are practically limited to the univariate
setting under the Bayesian paradigm, and to bivariate models under the
classical framework.

For the most popular sequential change detection rules such as CUSUM, EWMA,
and the Shiryaev-Roberts test, we develop integral equations and a concise
numerical method to compute a number of performance metrics, including average
detection delay and average time to false alarm. We pay special attention to
the Shiryaev-Roberts procedure and evaluate its performance for various
initialization strategies.

The computational complexity of different steps of the basic SSA is
discussed. It is shown that the use of the general-purpose "blackbox" routines
(e.g. found in packages like LAPACK) leads to huge waste of time resources
since the special Hankel structure of the trajectory matrix is not taken into
account. We outline several state-of-the-art algorithms (for example,
Lanczos-based truncated SVD) which can be modified to exploit the structure of
the trajectory matrix. The key components here are hankel matrix-vector
multiplication and hankelization operator.

We present a sequential Monte Carlo sampler variant of the partial rejection
control algorithm, and show that this variant can be considered as a sequential
Monte Carlo sampler with a modified mutation kernel. We prove that the new
sampler can reduce the variance of the incremental importance weights when
compared with standard sequential Monte Carlo samplers. We provide a study of
theoretical properties of the new algorithm, and make connections with some
existing algorithms.

This is the compilation of our comments submitted to the Journal of the Royal
Statistical Society, Series B, to be published within the discussion of the
Read Paper of Andrieu, Doucet and Hollenstein.

There is a growing interest in the literature for adaptive Markov Chain Monte
Carlo methods based on sequences of random transition kernels $\{P_n\}$ where
the kernel $P_n$ is allowed to have an invariant distribution $\pi_n$ not
necessarily equal to the distribution of interest $\pi$ (target distribution).
These algorithms are designed such that as $n\to\infty$, $P_n$ converges to
$P$, a kernel that has the correct invariant distribution $\pi$. Typically, $P$
is a kernel with good convergence properties, but one that cannot be directly
implemented.

A "cocktail algorithm" is proposed for numerical computation of (approximate)
D-optimal designs. This new algorithm combines and extends the multiplicative
algorithm of Silvey et al. (1978) and the vertex exchange method (VEM) of
Bohning (1986), and shares their simplicity and monotonic convergence
properties. Numerical examples show that the cocktail algorithm can lead to
dramatically improved speed, sometimes by orders of magnitude, relative to
either VEM or the multiplicative algorithm.

In this tutorial we schematically illustrate four algorithms:

(1) ABC rejection for parameter estimation

(2) ABC SMC for parameter estimation

(3) ABC rejection for model selection on the joint space

(4) ABC SMC for model selection on the joint space.

We consider Bayesian analysis of a class of multiple changepoint models.
While there are a variety of efficient ways to analyse these models if the
parameters associated with each segment are independent, there are few general
approaches for models where the parameters are dependent. Under the assumption
that the dependence is Markov, we propose an efficient online algorithm for
sampling from an approximation to the posterior distribution of the number and
position of the changepoints. In a simulation study, we show that the
approximation introduced is negligible.

The Lasso is a very well known penalized regression model, which adds an
$L_{1}$ penalty with parameter $\lambda_{1}$ on the coefficients to the squared
error loss function. The Fused Lasso extends this model by also putting an
$L_{1}$ penalty with parameter $\lambda_{2}$ on the difference of neighboring
coefficients, assuming there is a natural ordering. In this paper, we develop a
fast path algorithm for solving the Fused Lasso Signal Approximator that
computes the solutions for all values of $\lambda_1$ and $\lambda_2$.

The bivariate distribution with exponential conditionals (BEC) is introduced
by Arnold and Strauss [Bivariate distributions with exponential conditionals,
J. Amer. Statist. Assoc. 83 (1988) 522--527]. This work presents a simple and
fast algorithm for simulating random variates from this density.

Scaling of proposals for Metropolis algorithms is an important practical
problem in MCMC implementation. Criteria for scaling based on empirical
acceptance rates of algorithms have been found to work consistently well across
a broad range of problems. Essentially, proposal jump sizes are increased when
acceptance rates are high and decreased when rates are low. In recent years,
considerable theoretical support has been given for rules of this type which
work on the basis that acceptance rates around 0.234 should be preferred.

Recently developed adaptive Markov chain Monte Carlo (MCMC) methods have been
applied successfully to many problems in Bayesian statistics. Grapham is a new
open source implementation covering several such methods, with emphasis on
graphical models for directed acyclic graphs.

In some studies requiring predictive and CPU-time consuming numerical models,
the sampling design of the model input variables has to be chosen with caution.
For this purpose, Latin hypercube sampling has a long history and has shown its
robustness capabilities. In this paper we propose and discuss a new algorithm
to build a Latin hypercube sample (LHS) taking into account inequality
constraints between the sampled variables. This technique, called constrained
Latin hypercube sampling (cLHS), consists in doing permutations on an initial
LHS to honor the desired monotonic constraints.

The CUSUM procedure is known to be optimal for detecting a change in
distribution under a minimax scenario, whereas the Shiryaev-Roberts procedure
is optimal for detecting a change that occurs at a distant time horizon. As a
simpler alternative to the conventional Monte Carlo approach, we propose a
numerical method for the systematic comparison of the two detection schemes in
both settings, i.e., minimax and for detecting changes that occur in the
distant future.