## averaged acceptance ratios

Posted in Statistics with tags , , , , , , , , , , , , , on January 15, 2021 by xi'an

In another recent arXival, Christophe Andrieu, Sinan Yıldırım, Arnaud Doucet, and Nicolas Chopin study the impact of averaging estimators of acceptance ratios in Metropolis-Hastings algorithms. (It is connected with the earlier arXival rephrasing Metropolis-Hastings in terms of involutions discussed here.)

“… it is possible to improve performance of this algorithm by using a modification where the acceptance ratio r(ξ) is integrated with respect to a subset of the proposed variables.”

This interpretation of the current proposal makes it a form of Rao-Blackwellisation, explicitly mentioned on p.18, where, using a mixture proposal, with an adapted acceptance probability, it depends on the integrated acceptance ratio only. Somewhat magically using this ratio and its inverse with probability ½. And it increases the average Metropolis-Hastings acceptance probability (albeit with a larger number of simulations). Since the ideal averaging is rarely available, the authors implement a Monte Carlo averaging version. With applications to the exchange algorithm and to reversible jump MCMC. The major application is to pseudo-marginal settings with a high complexity (in the number T of terms) and where the authors’ approach does scale efficiently with T. There is even an ABC side to the story as one illustration is made of the ABC approximation to the posterior of an α-stable sample. As an encompassing proposal for handling Metropolis-Hastings environments with latent variables and several versions of the acceptance ratios, this is quite an interesting paper that I think we will study in further detail with our students.

## pseudo-marginal MCMC with minimal replicas

Posted in Books, Statistics, University life with tags , , , , , , on November 16, 2016 by xi'an

A week ago, Chris Sherlock, Alexandre Thiery and Anthony Lee (Warwick) arXived a short paper where they show that it is most efficient to use only one random substitute to the likelihood when considering a pseudo-marginal algorithm. Thus generalising an earlier paper of Luke Bornn and co-authors I commented earlier. A neat side result in the paper is that the acceptance probability for m replicas [in the pseudo-marginal approximation] is at most m/s the acceptance probability for s replicas when s<m. The main result is as follows:

There is a (superficial?) connection with the fact that when running Metropolis-within-Gibbs there is no advantage in doing more than one single Metropolis iteration, as the goal is to converge to the joint posterior, rather than approximating better the full conditional…

## a pseudo-marginal perspective on the ABC algorithm

Posted in Mountains, pictures, Statistics, University life with tags , , , , , , , , on May 5, 2014 by xi'an

My friends Luke Bornn, Natesh Pillai and Dawn Woodard just arXived along with Aaron Smith a short note on the convergence properties of ABC. When compared with acceptance-rejection or regular MCMC. Unsurprisingly, ABC does worse in both cases. What is central to this note is that ABC can be (re)interpreted as a pseudo-marginal method where the data comparison step acts like an unbiased estimator of the true ABC target (not of the original ABC target, mind!). From there, it is mostly an application of Christophe Andrieu’s and Matti Vihola’s results in this setup. The authors also argue that using a single pseudo-data simulation per parameter value is the optimal strategy (as compared with using several), when considering asymptotic variance. This makes sense in terms of simulating in a larger dimensional space but what of the cost of producing those pseudo-datasets against the cost of producing a new parameter? There are a few (rare) cases where the datasets are much cheaper to produce.

## MCMC using an approximation

Posted in R, Statistics with tags , , , , on April 17, 2014 by xi'an

This paper by Andrès Christen and Colin Fox is an reexamination of an earlier paper by Fox and Nicholls (1997) where the idea of testing for acceptance before computing the exact likelihood was first suggested. In the event of a low acceptance rate, and of a good enough approximation, the overall cost can be significantly reduced obviously. The issue is then to recover the exact target. In practice, the original proposal q is used to generate a value θ’ that is tested against the approximate target f⁰. If accepted, θ’ is then tested against the true target f, using a pseudo proposal q⁰ which is simply translating the earlier step. A sort of delayed acceptance that only efficiently reduces the computing cost when the approximation is good enough. Since, overall, the probability of acceptance of a given θ’ is smaller with this scheme. (Or, as stated within the paper, the original kernel dominates the new one in Peskun’s sense.) The most interesting question about this paper is how to compute the approximation. Subsampling the original data and hence using a small enough part of the likelihood function reduces the computing cost but induces variability.

#prior p~Be(7.5,.5)
#data  x~B(N,p)
#posterior p|x~Be(x+7.5,N+.5-x)
T=10^4;N=100;y=32;a=7.5;b=.5

p=rep(.2,T)
for (t in 2:T){
pp=p[t-1]+runif(1,-.1,.1) #random walk
p[t]=p[t-1] #by default rejection
if ((pp>0)&(pp<1)){
ratio=dbinom(y,N,pp)/dbinom(y,N,p[t-1])
if ((runif(1)<ratio)&(runif(1)<
dbeta(pp,a,b)/dbeta(p[t-1],a,b)))
p[t]=pp
}}


One application of the method I checked during my stay in Leuven was to separate the likelihood from the prior and compare only one term to a uniform at a time. For instance, with a costly prior, the approximate target f⁰ could be the likelihood, in which case the first acceptance probability would be based on the ratio of the likelihoods and the second acceptance probability on the ratio of the priors. As shown by the above graph, on a standard normal-normal example, the true posterior and the density estimate of the simulated sample almost coincide (for 10⁵ MCMC iterations)…

Actually, the principle also applies to a product of individual terms as in a likelihood so each term can be evaluated separately! This increases the variability and the potential for rejection but if each likelihood term is costly to compute this could bring some improvement. The graphs and the R code below illustrate an implementation in the Beta-binomial case when the binomial observation is replaced with a sequence of Bernoulli observations. The fit is adequate on 10⁵ iterations, but the autocorrelation in the sequence is high (the ACF is for the 100 times thinned sequence) and the acceptance rate down to 5%.

#even better, use the likelihood decomposition
T=10^5
yful=rep(0,N)
yful[sample(1:N,y)]=1 #pseudo Bernoulli sequence

p=rep(.2,T)
for (t in 2:T){
pp=p[t-1]+runif(1,-.1,.1) #random walk
p[t]=p[t-1] #by default rejection

if ((pp>0)&(pp<1)){
n=1
ratio=dbinom(yful[n],1,pp)/dbinom(yful[n],1,p[t-1])
while ((n<N)&(runif(1)<ratio)){
n=n+1
ratio=dbinom(yful[n],1,pp)/dbinom(yful[n],1,p[t-1])}
if ((n==N)&(runif(1)<ratio)){
if (runif(1)<dbeta(pp,a,b)/dbeta(p[t-1],a,b))
p[t]=pp
}
}}


$\prod_{i=1}^n\frac{L(\theta^\prime|x_i)}{L(\theta|x_i)}.$