{'fraction_non_alphanumeric': 0.03913385826771654, 'fraction_numerical': 0.0165748031496063, 'mean_word_length': 4.327390939597316, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Annealed Important Sampling for Models with Latent Variables

This paper is concerned with Bayesian inference when the likelihood is analytically intractable but can be unbiasedly estimated. We propose an annealed importance sampling procedure for estimating expectations with respect to the posterior. The proposed algorithm is useful in cases where finding a good proposal density is challenging, and when estimates of the marginal likelihood are required. The effect of likelihood estimation is investigated, and the results provide guidelines on how to set up the precision of the likelihood estimation in order to optimally implement the procedure. The methodological results are empirically demonstrated in several simulated and real data examples.


Introduction
Models incorporating latent variables are very popular in many statistical applications. For example, generalized linear mixed models (Fitzmaurice et al., 2011), that use latent variables to account for dependence between observations, appear in the genetics, social and medical sciences literatures as well as many other areas of statistics. State space models (Durbin and Koopman, 2001), whose latent variables follow a Markov process, are used in economics, finance and engineering. Gaussian process classifiers (see, e.g., Filippone, 2013;Rasmussen and Williams, 2006), that use a set of latent variables distributed as a Gaussian process to account for uncertainty in predictions, are used in computer science. Inference about the model parameters θ in latent variable models can be challenging because the likelihood is expressed as an integral over the latent variables. This integral is analytically intractable in general. It can also be computationally challenging when the dimension of the latent variables is high. Recent work by Beaumont (2003); Andrieu and Roberts (2009); Andrieu et al. (2010) shows that it is possible to carry out Bayesian inference in latent variable models by using an unbiased estimate of the likelihood within a Markov chain Monte Carlo (MCMC) sampling scheme. This method is known as particle MCMC. However, the resulting Markov chain can often be trapped in local modes and it is difficult to asses if it has converged. We find that even for simple models it is difficult for the Markov chain to mix adequately and the chain can take a long time to converge if the log of the estimated likelihood is too variable. The marginal likelihood is often used to choose between models.
Another drawback of MCMC in general, and particle MCMC in particular, is that it is often difficult to use it to estimate the marginal likelihood.
Another approach to Bayesian inference for models with latent variables is importance sampling squared (IS 2 ) proposed in Tran et al. (2013). They show that importance sampling (IS) with the likelihood replaced by its unbiased estimate is still valid for estimating expectations with respect to the posterior. IS 2 offers several advantages over particle MCMC: (1) It is easy to estimate the standard errors of the estimators; (2) It is straightforward to parallelize the computation; (3) It is straightforward to estimate the marginal likelihood. However, as is typical of importance sampling algorithms, a potential drawback with IS 2 is that its performance may depend heavily on the proposal density for θ. A good proposal density may be difficult to obtain in complex models.
When it is possible to evaluate the likelihood, annealed importance sampling (AIS) (Neal, 2001) is a useful method for estimating expectations with respect to the posterior for θ. AIS is an importance sampling method in which samples are first drawn from an easily-generated distribution and then moved towards the distribution of interest through Markov kernels. AIS explores the parameter space efficiently and is useful in cases where the target distribution is multimodal and/or when choosing an appropriate proposal density is challenging.
This article proposes an AIS algorithm for Bayesian inference when working with an estimated likelihood, which we denote as AISEL (annealed importance sampling with an estimated likelihood). Our first contribution is to show that the algorithm is valid for estimating expectations with respect to the exact posterior when the likelihood is estimated unbiasedly.
As with particle MCMC and IS 2 , it is important to understand the effect of estimating the likelihood on the resulting inference. The second contribution of this article is to answer this question by comparing the efficiency of AISEL with the efficiency of the corresponding AIS procedure that assumes that the likelihood is available. We show that the ratio of the efficiency of AISEL to that of AIS is smaller than or equal to 1, and the ratio is equal to 1 if and only if the estimate of the likelihood is exact. This ratio decreases exponentially with the product of the variance of the log of the estimated likelihood and a term that depends on the the annealing schedule in the AIS algorithm. The term based on the annealing schedule is small if the annealing schedule evolves slowly. This result allows us to understand how much accuracy is lost when working with an estimated likelihood. An attractive feature of AISEL is that it is more robust than IS 2 and particle MCMC to the variability of the log likelihood estimate. This is important when only highly variable estimates of the likelihood are available, which often occurs if it is expensive to obtain accurate estimates of the likelihood.
The third contribution of the article is to provide theory and practical guidelines for optimally choosing the number of particles to estimate the likelihood so as to minimize the overall computational cost for a given precision. The fourth contribution is to describe an efficient yet simple method to compute the marginal likelihood, which is important for model choice.
The SMC 2 algorithm of Chopin et al. (2013) sequentially updates the posterior of the model parameters as new observations arrive. The validity of the method is justified because the likelihood estimated by the particle filter is unbiased. In contrast, our AIS algorithm uses all the data and is static. As discussed in their paper (Chopin et al., 2013, Section 5.2), the reasoning used to justify SMC 2 does not apply to tempered sampling in the spirit of the AIS, because a tempered likelihood estimator is not an unbiased estimator of the corresponding tempered likelihood. However, Section 3.1 uses variable augmentation to justify the validity of the AIS method when working with an unbiased likelihood estimate.
We illustrate the proposed methodology through a simulated example, as well as the analysis of a Pound/Dollar exchange rate dataset using a stochastic volatility model. We show in these examples that the AIS method leads to efficient inference when optimally implemented.
The article is organized as follows. Section 2 reviews the original AIS of Neal (2001).
Section 3 presents the main results. Section 4 presents the examples and Section 5 concludes.
The technical proofs are in the Appendix.

Annealed importance sampling
Let p(y|θ) be the density of the data y, where θ is a parameter vector belonging to a space Θ ⊂ R d . Let p(θ) be the prior for θ and π(θ) ∝ p(θ)p(y|θ) its posterior. We are interested in the case where the likelihood p(y|θ) is analytically intractable but can be unbiasedly estimated.
• Set θ i = θ (T ) and compute the unnormalized weight Note that the algorithm is parallelizable. Neal (2001) shows that the above algorithm is an IS procedure operating on the extended space Θ T with the artificial target density proportional to f (θ (1) , ..., θ (T ) ) = η a T (θ (T ) )L ξa T −1 (θ (T ) , θ (T −1) )...L ξa 1 (θ (2) , θ (1) ), where is a backward kernel density, and the proposal density The original target density π = ξ a T is the last marginal of f (θ (1) ,...,θ (T ) ) because for all t. This shows that AIS is a valid IS method with the weight w ∝ f /g. Hence, the The AIS procedure explores the parameter space efficiently, and is useful when the target distribution is multimodal and when choosing an appropriate proposal density is challenging.
3 Annealed importance sampling with an estimated likelihood This section presents an AISEL algorithm for estimating the integral (1) when the likelihood is analytically intractable but can be estimated unbiasedly. Let p N (y|θ) denote an estimator of p(y|θ) with N the number of particles used to estimate the likelihood. We define a sequence of functions η at (θ), t = 0,...,T , by We propose the following algorithm for generating M weighted samples approximate the posterior π(θ); see Section 3.1 Algorithm 1 (AISEL). Let π 0 (θ) be some easily-generated density.
2. For t = 1,...,T (i) Weighting: compute the unnormalized weights and set the new normalized weights (iii) Markov move: for each i = 1,...,M, move the sample θ i according to a Metropolis-Hastings step as follows. Let q t (θ|θ c ) be a proposal with θ c = θ i the current state.
Generate θ p ∼ q t (θ|θ c ) and set θ i = θ p with probability Otherwise set θ i = θ c .
If p N (y|θ)=p(y|θ), the above algorithm is a special case of the SMC sampler in Del Moral et al. The algorithm is also closely related to the resample-move algorithm of Gilks and Berzuini (2001), except that they perform the resampling step in every iteration t.

Formal justification
The output of Algorithm 1 is weighted samples approximate π(θ), we make the following assumption.
Let us write p N (y|θ) as p(y|θ)e z where z = log p N (y|θ)−log p(y|θ) is a random variable whose distribution is governed by the randomness occurring when estimating the likelihood p(y|θ). Let g N (z|θ) be the density of z. Assumption 1 implies that We define as the joint density of θ and z on the extended space Θ= Θ⊗R. Then its first marginal is the We define the following sequence of interpolation densities on Θ = Θ⊗R with Note that ξ a T (θ,z) = π N (θ,z), which is our new target density defined on Θ. Algorithm 1 is entirely equivalent to the following procedure.
Remark 1. It is advisable to perform the Markov move step over a few burn-in iterations so that the samples move closer to the equilibrium distribution.
Remark 2. If ϕ satisfies the conditions in Theorem 1 of Chopin (2004), then with the asymptotic variance σ 2 AISEL (ϕ) defined recursively as in Chopin (2004). Except for the special case in which no resampling steps in Algorithm 1 are performed, the asymptotic variance σ 2 AISEL (ϕ), and therefore the variance of ϕ AISEL , does not admit a closed form. A natural and potential technique to estimate Var( ϕ AISEL ) is to run Algorithm 1 in batches independently and in parallel. Then, we have several independent batches of weighted ..,R with r M r = M, and the corresponding R independent estimates ϕ (r) AISEL of E π (ϕ). The variance of the estimator ϕ AISEL can be estimated by If no resampling steps are performed, then the particles θ i are independent and we have a closed form expression for estimating the asymptotic variance of ϕ AISEL It is straightforward to show that this estimate is consistent. However, it is important to perform resampling if necessary.

Estimating the marginal likelihood
Marginal likelihood (2) is important for model comparison purposes. Except for some trivial cases, computing the marginal likelihood is challenging because of its integral form. Friel and Pettitt (2008) propose a very efficient method, called power posterior method, for estimating the marginal likelihood, that exploits the temparing sampling framework as in AIS. This section extends the power posterior method to the case with latent variables.
We consider for now a continuous sequence of interpolation densities (6) as follows We have the following result. The proof is in the Appendix.
Let f (s) denote the integrand in the right side of (9). The scalar integral in (9) can be deterministically approximated by with {a t , t = 0,...,T } the annealing schedule as in Section 2. The function f (a t ) can be estimated by the output of Algorithm 1 after iteration t.
This approach of estimating the marginal likelihood fits naturally to the AISEL procedure and is straightforward to implement. Note that the values p N (y|θ (t) i ) in (11) can be used for the weighting step in iteration t+1 of Algorithm 1, hence no extra computation is needed except calculation in (10) and (11).

The effect of estimating the likelihood
The efficiency of an IS procedure with proposal density g and weights w i is often measured by the effective sample size defined by (see, e.g. Neal (2001) and Liu (2001, Chapter 2 is often called the coefficient of variation of the unnormalized weights w i . The bigger the ESS the more efficient the IS procedure. This section investigates how much the ESS is reduced when working with an estimated likelihood. We consider the case of the original AIS procedure, i.e. Algorithm 1 without the resampling step, and work with the notation in Algorithm 1'. Write x=(θ,z) and x (t) =(θ (t) ,z (t) ). Without the resampling step, Algorithm 1' can be written as • Generate x (1) = (θ (1) ,z (1) ) ∼ ξ a 0 (θ,z) = π 0 (θ)g N (z|θ).
• Set (θ i ,z i ) = (θ (T ) ,z (T ) ) and compute the corresponding unnormalized weight Denote by ESS AIS and ESS AISEL the effective sample sizes of the AIS procedures when the likelihood is given and when it is estimated, respectively.
We make the following assumption which is satisfied in almost cases.
Following Pitt et al. (2012), we make the following further assumptions.

Assumption 3(i) is justified by Lemma 1. Assumption 3(ii) keeps the variance Var(z)
constant across different values of θ, thus making it easy to associate the ESS with σ. Under Assumption 3, the density g N (z|θ) depends only on σ and is denoted by g(z|σ).
As in Neal (2001), this assumption separates out the effect related to Markov chain convergence and allows us to study the effect of estimating the likelihood on the sequential sampling scheme.
The theorem, whose proof is in the Appendix, shows that the efficiency is reduced by the factor exp(τ σ 2 ) when working with an estimated likelihood. If a t = 1/T , then τ = 1/T for all t, and the theorem shows that increasing T and thus making the ξ at closer to each other helps improve efficiency.
In the IS 2 approach, Tran et al. (2013) show that ESS IS 2 /ESS IS = exp(−σ 2 ), where ESS IS 2 and ESS IS are the effective sample sizes of IS when the likelihood is estimated and given, respectively. Similarly in particle MCMC, Pitt et al. (2012) show that efficiency is reduced by a factor of approximately exp(−σ 2 ) when working with an estimated likelihood. We can see that the accuracy of likelihood estimation is less important in AISEL than in IS 2 and particle MCMC, because the factor exp(−τ σ 2 ) can be made small when τ is decreased. This means that AISEL can be more robust than IS 2 and particle MCMC in cases where we only have a rough estimate of the likelihood, or it is expensive to obtain an accurate estimate of the likelihood.

Practical guidelines on selecting the number of particles
This section studies how to select the number of particles N optimally. A large number of particles N results in a precise likelihood estimate, and therefore an accurate estimate of E π (ϕ), but at a greater computational cost. A small N leads to a large variance of the likelihood estimator, so we need a larger number of importance samples M in order to obtain the desired accuracy of the AISEL estimator. In either case, the computation is expensive. It is important to select an optimal value of N that minimizes the computational cost.
The time to compute the likelihood estimate p N (y|θ) can be written as τ 0 +N(θ)τ 1 where τ 0 ≥0 and τ 1 > 0 (Tran et al., 2013). For example, if p N (y|θ) is estimated by IS, then τ 0 is the overhead cost spent on estimating the proposal density and τ 1 is the computing time used to generate each sample and compute the weight. Note that under Assumption 3, N depends The variance of the AISEL estimator is approximated as with Var π (ϕ) = E π (ϕ−E π (ϕ)) 2 . From (12), Let P * be a prespecified precision. Then we need approximately particles in order to have that precision. The required computing time to run AISEL is in characterizes the computing as a function of σ 2 , which is minimized at The optimal number of particles N is such that Var N,θ (z) = Var(log p N (y|θ)) = σ 2 opt .
Let Var N,θ (z) be an estimate of Var N,θ (z), which can be obtained by using, e.g., the delta method or the jackknife. See Tran et al. (2013) for more details. We suggest the following practical guidelines for tuning the optimal number of particles N. Note that N generally depends on θ but this dependence is suppressed for notational simplicity.
The case τ 0 >0. First, we need to estimateγ 2 . Let {θ 1 ,...,θ J } be a few initial draws from the initial density ξ 0 (θ). Then, start with some large N 0 ,γ 2 can be initially estimated by as Var N 0 ,θ j (z)= γ 2 (θ j )/N 0 . By substituting this estimate ofγ 2 into (16) we obtain an estimate σ 2 opt of σ 2 opt . We now can start Algorithm 1 and update γ 2 (and therefore σ 2 opt ) as we go. For each draw of θ, we start with some small N and increase N if Var N,θ (z) > σ 2 opt .

Time normalized variance.
In the examples in Section 4 we use the time normalized variance (TNV) as a measure of efficiency (Tran et al., 2013) of a sampling procedure. The TNV of the AISEL estimator ϕ AISEL is defined as where τ (M,N) is the total CPU time used to run the AISEL procedure with M importance samples and N particles. From (13) and (14), and is propotional to CT * (σ 2 ).
Remark 3. Letting the optimal N opt = N opt (θ) depend on θ is theoretically interesting but might in some cases be ultimately inefficient. The reason is that extra computing time is needed to tune N for each θ. A simple strategy is to make approximation that γ 2 (θ) ≈γ 2 .
Then, the optimal number of particles, which is constant across θ, is determined as We follow this strategy in the examples below.
We first run Algorithm 1 to generate M =100 importance samples with N =10, and obtain τ 0 = 7.2×10 −3 ,τ 1 = 5.9×10 −4 andγ 2 = 17.7. This gives σ 2 opt = 2.6 and the optimal number of particles (as constant across θ) N opt = 7. We then run Algorithm 1 to generate M = 5000 importance samples for five values of N, N = 1,7,10,20 and N = 50. In order to be able to estimate the variance of the estimator, Algorithm 1 is run in parallel as described in Remark 2 with R = 20 batches. Figure 1 plots the time normalized variance in (18) versus N. The TNV is averaged over the posterior mean estimates of the four parameters β 0 , β, σ 2 1 and σ 2 2 . The TNV appears to be minimized at N =10, close to the theoretical optimal value N =7. The results suggests that the TNV is weakly sensitive around the optimal value of N. The efficiency decreases linearly when N is higher than the optimal value, whereas the efficiency can deteriorate exponentially when N is below the optimal. This phenomenon is also observed in the IS 2 method (Tran et al., 2013). In practice, it is therefore advisable to use for N a value which is slightly bigger than N opt . Table 1 reports the estimate of the posterior mean when using N = 10 particles.

Real Data Example: The SV model
We analyze the Pound/Dollar data set (Kim et al. (1998)) using both the standard SV model and the SV model with leverage effect. For the standard SV model the measurement equation, for the t th observation y t , is given by where ε t follows the standard normal distribution, and h t is a latent variable defined for t = 1,2,...,n−1, as where µ is the unconditional mean, φ is the level of persistence, σ η is the scaling parameter, for the latent process and η t follows the standard normal distribution. The SV model given in (20) and (21) is completed by the initial state, given as For the analysis we assume, a priori, that µ is normally distributed, such that µ∼N(0,100), we assume that the log of variance, h t , is generated by a stationary process with positive autocorrelation, and where the persistence parameter, φ, follows a Beta distribution where φ ∼ Be(15,1.5) and we assume σ η has an inverted gamma prior IG(10,0.1).
To implement Algorithm 1, we use the bootstrap particle filter of Gordon et al. (1993) to obtain an unbiased estimate of the log-likelihood. For this specific implementation τ 0 = 0, for which (16) implies that the variance of the estimated log-likelihood should be 1 τ for the annealed IS scheme that we implement. In the analysis of the SV model we set T = 15. To estimate the optimal number particles, we compute the log-likelihood one thousand times, evaluated at parameter values that are typical for the SV model. Specifically, we set µ=−0.6, φ = 0.98 and we set σ η = 0.16. This leads us to the conclusion that 24 is the optimal number of particles for this model and data set. Algorithm 1 also requires us to set M, which we set to M = 1000.
We also need to specify step (iii) of Algorithm 1, which is the Markov move step. Specifically, we employ 5 random walk Metropolis Hastings (RWMH) steps for each Markov move; see Robert and Casella (1999) for further details on the RWMH algorithm. Similarly, to Chopin (2002) (in the sequential Monte Carlo context), we take the covariance of the RWMH algorithm as the covariance of the current set of (annealed IS) particles. The covariance is scaled by the parameter α, where α is adjusted at each move step, based on the acceptance rate of the previous move step. Specifically, at each move step we update the scale parameter, where MF denotes a multiplication factor that scales α from the previous period. The multiplication factor is determined by the acceptance rate (AR) from the previous Markov move step.  We report in Table 2 the multiplication factor used in (23), given the value of the acceptance rate from the previous period. We find adapting in this fashion works well in all the examples we have considered so far, including the ones we consider in this paper.
The output for Figure 2 is produced by running Algorithm 1, using the method described in Remark 3 with R = 200 batches. For this part of the analysis here we set a t = t/n. As in the simulated example, the theoretically optimal number of simulated particles, is close the what is empirically observed as optimal. We also observe that the penalty for using too few particles, with respect to the time corrected measure of accuracy, can be much greater than using too many particles.
We use the annealed IS algorithm for estimation of the SV model on the Pound/Dollar data set. Here we set a t = (t/n) 3 . This ensures that we move away from the prior very  slowly initially, and ensures that the effective sample size is high as we move over the entire extended state space. This is important when estimating the marginal likelihood as we require accurate estimates of the expected value of the log-likelihood across the entire temperature range, and not just at the end of the estimation process, which is the case for parameter estimates. The results are reported in Table 3. We report both the posterior mean and standard deviation for each of the parameters. The time taken for estimation is 176 seconds, using the Julia programming language (Bezanson et al. (2012)), on a Core i7 Linux box, with 16 Gigabytes of RAM. Note the code has not been parallelized, so even better performance could be achieved, with a more highly optimized implementation.
Unlike the standard SV model, in (20), (21) and (22), in which the measurement and state disturbances are independent, the SV model with leverage effect allows for correlation between ε t and η t ; see Omori et al. (2007), for further details. Specifically, it is assumed that where |ρ| < 1. As for the standard SV model we use the bootstrap particle filter to obtain an estimate of the log-likelihood. We set the number of particles to 20, which corresponds to our estimate of the optimal number of particles. The other algorithmic parameters remain the same as for the standard SV model.  Output from the analysis of the Pound/Dollar data set using the SV model with leverage is reported in Table 4. The analysis took 172 seconds on using the Julia programming language (Bezanson et al. (2012)), on a Core i7 Linux box, with 16 Gigabytes of RAM. There is little evidence of leverage effect, for the Pound/Dollar data set, based on this analysis. In particular, the log marginal likelihood strongly favours the standard SV model, which isn't surprising given that the estimate of ρ is close to zero.

Conclusions
We have presented the annealed IS algorithm for Bayesian inference in models with latent variables. The proposed AISEL method can be considered as a supplement to existing Monte Carlo methods for latent variable models, including particle MCMC (Beaumont, 2003;Andrieu and Roberts, 2009;Andrieu et al., 2010), SMC 2 (Chopin et al., 2013) and IS 2 (Tran et al., 2013). The theory and methodology presented in this paper are useful for Bayesian inference in latent variable models where the posterior distribution is multimodal and choosing an appropriate proposal density is challenging. An estimate of the log marginal likelihood is obtained as a byproduct of the estimation procedure.