Likelihood-free parallel tempering

Abstract

Approximate Bayesian Computational (ABC) methods, or likelihood-free methods, have appeared in the past fifteen years as useful methods to perform Bayesian analysis when the likelihood is analytically or computationally intractable. Several ABC methods have been proposed: MCMC methods have been developed by Marjoram et al. (2003) and by Bortot et al. (2007) for instance, and sequential methods have been proposed among others by Sisson et al. (2007), Beaumont et al. (2009) and Del Moral et al. (2012). Recently, sequential ABC methods have appeared as an alternative to ABC-PMC methods (see for instance McKinley et al., 2009; Sisson et al., 2007). In this paper a new algorithm combining population-based MCMC methods with ABC requirements is proposed, using an analogy with the parallel tempering algorithm (Geyer 1991). Performance is compared with existing ABC algorithms on simulations and on a real example.

Appendices

Appendix A: Matrices of exchange rates for ABC-PT with and without rings

Tables 7 and 8 compare accepted exchange rates between chains based on the modified toy example, using ABC-PT algorithms with or without rings. The rate between chains i and j is equal to the number of accepted exchange moves between these chains divided by the total number of iterations. It is observed that the use of rings allows more accepted exchange moves, especially between high tempered chains.

Table 7 ABC-PT with 3 rings for the modified toy example: the diagonal gives local acceptance rates, and outside the diagonal are given the proportions of accepted exchange moves

Table 8 ABC-PT for the modified toy example: the diagonal gives local acceptance rates, and outside the diagonal are given the proportions of accepted exchange moves

Appendix B: Formula for the tuberculosis example

We used the same notations than Tanaka et al. (2006). The number of cases of genotype i at time t is denoted by X _i (t), G(t) is the number of distinct genotypes that have existed in the population up to and including time t, and N(t) is the total number of cases at time t.

$$N(t)=\sum_{i=1}^{G(t)} X_i(t).$$

The genotypes are labeled 1,2,3,… for convenience, although the ordering has no meaning, except that i=1 represents the parental type from which others are descended (directly or indirectly). The three rates of the system are the birth rate per case per year α, the death rate per case per year δ, and the mutation rate per case per year θ. Tanaka et al. (2006) define the following probabilities:

$$P_{i,x}(t) = P\bigl(X_i(t)=x\bigr),$$

$$\bar{P}_n(t) = P\bigl(N(t)=n\bigr)\quad \hbox{and}\quad \tilde{P}_g(t)= P\bigl(G(t)=g\bigr).$$

The time evolution of P _i,x (t) is described by the following differential equations:

Initially there is only one copy of the ancestral genotype, hence the initial conditions are: P _i,x (0)=0 for all (i,x), except P _1,1(0)=1, and for i=2,3,4,…,P _i,0(0)=1. To take into account the creation of new genotypes, the probability $\tilde{P}_{g}(t)$ is described by the following differential equations (only a mutation can create a new genotype):

The initial condition is G(0)=1. Let t _g be the time when a new genotype g is created, we have P _g,1(t _g )=P(X _g (t _g )=1)=1, and P _g,x (t _g )=P(X _g (t _g )=x)=0 for x≠1.

The total number of cases N(t) is described by the following differential equations (only a birth or a death influence changes in this number):

The initial conditions are $\bar{P}_{1}(0)=1$ and $\bar{P}_{n}(0)=0$ for n≠1.