Bayesian Inference in the Presence of Intractable Normalizing Functions

TL;DR

This paper reviews and compares Monte Carlo algorithms for Bayesian inference in models with intractable normalizing functions, highlighting their efficiency and theoretical foundations through simulations and real data examples.

Contribution

It provides a unifying framework for understanding existing algorithms and offers practical recommendations for their application in complex statistical models.

Findings

Different algorithms vary in computational efficiency

Connections among algorithms are clarified

Guidelines for practitioners are provided

Abstract

Models with intractable normalizing functions arise frequently in statistics. Common examples of such models include exponential random graph models for social networks and Markov point processes for ecology and disease modeling. Inference for these models is complicated because the normalizing functions of their probability distributions include the parameters of interest. In Bayesian analysis they result in so-called doubly intractable posterior distributions which pose significant computational challenges. Several Monte Carlo methods have emerged in recent years to address Bayesian inference for such models. We provide a framework for understanding the algorithms and elucidate connections among them. Through multiple simulated and real data examples, we compare and contrast the computational and statistical efficiency of these algorithms and discuss their theoretical bases. Our study provides practical recommendations for practitioners along with directions for future research for MCMC methodologists.

Tables

Table 1. Table 1: Inference results for an Ising model on a 10 × \times 10 lattice. 20,000 MCMC samples are generated from each algorithm. The highest posterior density (HPD) is calculated by using coda package in R . The calculation of Effective Sample Size (ESS) follows Kass et al., ( 1998 ), Robert and Casella, ( 2013 ) . “Acc” represents acceptance probability.

$\theta =0.2$	Mean	95%HPD	ESS	Acc	Time(second)	ESS/Time
AVM	0.20	(0.08,0.32)	1527.88	0.37	12.46	122.62
Exchange	0.20	(0.08,0.33)	2061.23	0.49	15.29	134.81
DMH	0.21	(0.08,0.34)	2068.29	0.48	1.78	1161.96
AEX	0.21	(0.07,0.35)	1778.15	0.47	23.72	74.96
ALR	0.20	(0.08,0.34)	3647.34	0.59	9.33	390.93
RussianR	0.20	(0.06,0.33)	2650.83	0.50	13849.10	0.19
NoisyDMH	0.20	(0.06,0.33)	3347.76	0.59	71.02	47.14
Gold	0.20	(0.08,0.33)	9845.89

Table 2. Table 2: Inference results for 2-star in ERGM for Florentine business dataset. 30,000 MCMC samples are generated from each algorithm. The highest posterior density (HPD) is calculated by using coda package in R . The calculation of Effective Sample Size (ESS) follows Kass et al., ( 1998 ), Robert and Casella, ( 2013 ) . “Acc” represents acceptance probability.

${\theta }_2$	Mean	95%HPD	ESS	Acc	Time(second)	ESS/Time
DMH	1.24	(0.02,2.57)	1026.67	0.24	10.45	98.25
AEX	1.24	(0.17,2.58)	991.87	0.20	126.46	7.84
ALR	1.25	(0.16,2.52)	1456.57	0.33	2500.37	0.58
RussianR	1.27	(0.03,2.68)	1433.60	0.31	33534.96	0.04
NoisyDMH	1.28	(0.03,2.59)	1600.90	0.32	297.35	5.38
Gold	1.27	(0.08,2.50)	9655.90

Table 3. Table 3: Inference results for graphlet orbit factor for “NL” in EMON dataset. 40,000 MCMC samples are generated from each algorithm. The highest posterior density (HPD) is calculated by using the coda package in R . The calculation of Effective Sample Size (ESS) follows Kass et al., ( 1998 ), Robert and Casella, ( 2013 ) . “Acc” represents acceptance probability.

${\theta }_9\times {10}^2$	Mean	95%HPD	ESS	Acc	Time(minute)	ESS/Time
DMH	1.93	(0.61,2.98)	583.08	0.05	16.10	36.21
NoisyDMH	1.93	(0.08,3.45)	796.01	0.08	1026.15	0.78
Gold	1.89	(0.58,2.98)	5494.31

Table 4. Table 4: Inference results for outputs for θ 1 subscript 𝜃 1 \theta_{1} in attraction repulsion point process model ( n ≈ 200 ) n\approx 200) . 40,000 MCMC samples are generated from each algorithm. The highest posterior density (HPD) is calculated by using coda package in R . The calculation of Effective Sample Size (ESS) follows Kass et al., ( 1998 ), Robert and Casella, ( 2013 ) . “Acc” represents acceptance probability.

${\theta }_1$	Mean	95%HPD	ESS	Acc	Time(minute)	ESS/Time
DMH	1.20	(1.03,1.33)	1343.11	0.27	3.93	341.35
RussianR	1.19	(1.04,1.35)	1867.54	0.32	3299.76	0.57
NoisyDMH	1.20	(1.04,1.34)	2676.49	0.42	77.73	34.44
Gold	1.19	(1.03,1.34)	4397.39

Table 5. Table 5: Comparison between algorithms

		Exact or
Method	Class	inexact	Adaptation	Requirement	Ease of use
AVM	A	Exact	No	Perfect sampling	3
Exchange	A	Exact	No	Perfect sampling	2
DMH	A	Inexact	No	Nothing	1
AEX	A/L	Exact	Yes	Sufficient statistics	4
ALR	L	Exact	Yes	Sufficient statistics	4
Russian Roulette	L	Exact	No	Nothing	4
Noisy DMH	L	Inexact	No	Nothing	2

Table 6. Table 6: Summary of notations used for complexity calculations.

Symbol	Definition
$n$	the number of data points
$m$	current iteration index in MCMC
$N$	the number of samples used to construct importance sampling estimate
$N_0$	the number of preliminary runs to attain particles (AEX, ALR)
$N_1$	the number of preliminary runs for the auxiliary chain in AEX
$b$	the thinning interval for MCMC
$c$	the number of possible cores for parallel computing
$\tau$	random stopping time (ALR, Russian Roulette)
${\tilde{\tau }}^R$	expected stopping time for geometric Russian Roulette algorithm
$p$	dimension of the parameter
$d$	the number of particles
$S(\cdot )$	sufficient statistics
$L(\cdot )$	complexity function of evaluating $h(\cdot |\theta )$
$G(\cdot )$	complexity function of inner sampler

Table 7. Table 7: Computational complexity and memory size of the algorithms in the simulation examples. If costs are negligible, we didn’t include in the table. C, M, FC, FM denote complexity, memory, fixed complexity, fixed memory respectively. Asterisk symbols indicate calculated complexity is expected number.

DMH	Exponential family	Point process
C	$G\left(n\right)+L\left(n\right)$	$G\left(n^2\right)+3L\left(n^2\right)$
FM	$p$	$n^2$
NoisyDMH	Exponential family	Point process
C	$N\left[G\left(n\right)+L\left(n\right)\right]/c$	$NG\left(n^2\right)/c+\left[2N/c+1\right]L\left(n^2\right)$
FM	$p$	$n^2$
$\ast$AEX	Exponential family	Point process
C	$0.5\left[G\left(n\right)+L\left(n\right)\right]$	$0.5G\left(n^2\right)+\left[N_1/b+m+5\right]L\left(n^2\right)$
M	$p+2$	$n^2+3$
FC	$\left[N_0+0.5N_1\right]\left[G\left(n\right)+L\left(n\right)\right]$	$\left[N_0+0.5N_1\right]G\left(n^2\right)+\left[3N_0+N_1\right]L\left(n^2\right)$
FM	$p+dp+d^2+\left(p+2\right)N_1/b$	$n^2+dp+d^2+\left(n^2+3\right)N_1/b$
$\ast$ALR	Exponential family	Point process
C	$G\left(n\right)+L\left(n\right)$	$G\left(n^2\right)+\left[\left(d+2m+2\right)/c+1\right]L\left(n^2\right)$
M	$p+1$	$n^2+d+1$
FC	$\left(N_0+\tau \right)\left[G\left(n\right)+L\left(n\right)\right]$	$\left[N_0+\tau \right]G\left(n^2\right)+\left[3N_0+\tau d/c\right]L\left(n^2\right)$
FM	$p+dp+d$	$n^2+dp+d$
$\ast$Roulette	Exponential family	Point process
C	$N\left({\tilde{\tau }}^R+1\right)\left[G\left(n\right)+L\left(n\right)\right]/c$	$N\left({\tilde{\tau }}^R+1\right)G\left(n^2\right)/c+\left[2N\left({\tilde{\tau }}^R+1\right)/c+1\right]L\left(n^2\right)$
FM	$p$	$n^2$

Equations

$$\pi (\theta |\mathbf{x})\propto p(\theta )\frac{h(\mathbf{x}|\theta )}{Z(\theta )}.$$

$$\alpha =\min \left\{1,\frac{f({\mathbf{y}}^{\prime }|{\theta }^{\prime },\mathbf{x})p({\theta }^{\prime })h(\mathbf{x}|{\theta }^{\prime })\bcancel{Z(\theta )}h(\mathbf{y}|\theta )\bcancel{Z({\theta }^{\prime })}q(\theta |{\theta }^{\prime })}{f(\mathbf{y}|\theta ,\mathbf{x})p(\theta )h(\mathbf{x}|\theta )\bcancel{Z({\theta }^{\prime })}h({\mathbf{y}}^{\prime }|{\theta }^{\prime })\bcancel{Z(\theta )}q({\theta }^{\prime }|\theta )}\right\}.$$

$$\pi (\theta ,{\theta }^{\prime },\mathbf{y}|\mathbf{x})\propto p(\theta )L(\theta |\mathbf{x})q({\theta }^{\prime }|\theta )L({\theta }^{\prime }|\mathbf{y})=p(\theta )\frac{h(\mathbf{x}|\theta )}{Z(\theta )}q({\theta }^{\prime }|\theta )\frac{h(\mathbf{y}|{\theta }^{\prime })}{Z({\theta }^{\prime })}.$$

$$\min \left\{1,\frac{\bcancel{s(\{\theta ,{\theta }^{\prime }\}|\{{\theta }^{\ast },{\theta }^{\prime \ast }\})}\pi ({\theta }^{\prime },\theta ,\mathbf{y}|\mathbf{x})}{\bcancel{s(\{{\theta }^{\ast },{\theta }^{\prime \ast }\}|\{\theta ,{\theta }^{\prime }\})}\pi (\theta ,{\theta }^{\prime },\mathbf{y}|\mathbf{x})}\right\}=\min \left\{1,\frac{p({\theta }^{\prime })h(\mathbf{x}|{\theta }^{\prime })\bcancel{Z(\theta )}h(\mathbf{y}|\theta )\bcancel{Z({\theta }^{\prime })}q(\theta |{\theta }^{\prime })}{p(\theta )h(\mathbf{x}|\theta )\bcancel{Z({\theta }^{\prime })}h(\mathbf{y}|{\theta }^{\prime })\bcancel{Z(\theta )}q({\theta }^{\prime }|\theta )}\right\},$$

$$f_k(\mathbf{y}|\theta ,{\theta }^{\prime })=h(\mathbf{y}|{\theta }^{\prime }{)}^{1-{\beta }_k}h(\mathbf{y}|\theta {)}^{{\beta }_k},{\beta }_k=\frac{k}{K+1},\text{leavevmode}k=1,\dots ,K.$$

$$h(\mathbf{x}|{\theta }^{\prime })T_{{\theta }^{\prime }}^m(\mathbf{y}|\mathbf{x})=h(\mathbf{y}|{\theta }^{\prime })T_{{\theta }^{\prime }}^m(\mathbf{x}|\mathbf{y}).$$

$$\alpha =\min \left\{1,\frac{p({\theta }^{\prime })h(\mathbf{y}|\theta )T_{{\theta }^{\prime }}^m(\mathbf{x}|\mathbf{y})q(\theta |{\theta }^{\prime })}{p(\theta )h(\mathbf{x}|\theta )T_{{\theta }^{\prime }}^m(\mathbf{y}|\mathbf{x})q({\theta }^{\prime }|\theta )}\right\},$$

$$P(\mathbf{y}={\mathbf{x}}_l|{\theta }^{\prime })\propto \sum _{j=1}^{n+1}{w}_j^{(I_j)}\frac{h({\mathbf{x}}_j|{\theta }^{\prime })}{h({\mathbf{x}}_j|{\theta }^{(I_j)})}1\{{\mathbf{x}}_j={\mathbf{x}}_l\},\text{leavevmode}\text{leavevmode}l=1,...,n+1.$$

$$\begin{array}{rl}Z(\theta )& ={\int }_{\mathcal{X}}h(\mathbf{x}|\theta )d\mathbf{x}\\ & =\sum _{i=1}^{d}k(\theta ,{\theta }^{(i)}){\int }_{\mathcal{X}}\frac{h(\mathbf{x}|\theta )}{h(\mathbf{x}|{\theta }^{(i)})}h(\mathbf{x}|{\theta }^{(i)})d\mathbf{x}\\ & =\sum _{i=1}^{d}k(\theta ,{\theta }^{(i)})Z({\theta }^{(i)}){\int }_{\mathcal{X}}\frac{h(\mathbf{x}|\theta )}{h(\mathbf{x}|{\theta }^{(i)})}\frac{h(\mathbf{x}|{\theta }^{(i)})}{Z({\theta }^{(i)})}d\mathbf{x}.\end{array}$$

$${\hat{Z}}_{n+1}(\theta )=\sum _{i=1}^{d}k(\theta ,{\theta }^{(i)})\exp (c_{n+1}^{(i)})\frac{{\sum }_{k=1}^{n+1}\frac{h({\mathbf{x}}_k|\theta )}{h({\mathbf{x}}_k|{\theta }^{(i)})}1\{I_k=i\}}{{\sum }_{k=1}^{n+1}1\{I_k=i\}},$$

$$\alpha =\min \left\{1,\frac{p({\theta }^{\prime })h(\mathbf{x}|{\theta }^{\prime })q(\theta |{\theta }^{\prime })}{p(\theta )h(\mathbf{x}|\theta )q({\theta }^{\prime }|\theta )}\frac{1}{N}\sum _{i=1}^N\frac{h({\mathbf{y}}_i|\theta )}{h({\mathbf{y}}_i|{\theta }^{\prime })}\right\}.$$

$$L(\theta |\mathbf{x})=\frac{1}{Z(\theta )}\exp \left\{\theta S(\mathbf{x})\right\},\text{leavevmode}\text{leavevmode}S(\mathbf{x})=\sum _{i=1}^m\sum _{j=1}^{n-1}{x}_{i,j}{x}_{i,j+1}+\sum _{i=1}^{m-1}\sum _{j=1}^n{x}_{i,j}{x}_{i+1,j},$$

$$L(\theta |\mathbf{x})=\frac{1}{Z(\theta )}\exp \left\{{\theta }_1{S}_1(\mathbf{x})+{\theta }_2{S}_2(\mathbf{x})+{\theta }_3{S}_3(\mathbf{x})+{\theta }_4{S}_4(\mathbf{x})\right\},$$

$$S_l(\mathbf{x})=\sum _{i=1}^n\left(\genfrac{}{}{0ex}{}{x_{i+}}{l}\right),l=1,2,3;\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}{S}_4(\mathbf{x})=\sum _{i<j<k}{x}_{i,j}{x}_{j,k}{x}_{k,i}.$$

$$L(\theta |\mathbf{x})=\frac{1}{Z(\theta )}\exp \left\{\sum _{l=1}^9{\theta }_l{S}_l(\mathbf{x})\right\},$$

$$S_1(\mathbf{x})=\sum _{i=1}^n\left(\genfrac{}{}{0ex}{}{x_{i+}}{1}\right),\text{leavevmode}\text{leavevmode}{S}_2(\mathbf{x})=\sum _{i<j<k}{x}_{i,j}{x}_{j,k}{x}_{k,i}$$

$$S_3(\mathbf{x})=\sum _{i<j}{x}_{i,j}(1\{\text{mbox}{sponsor}_i=\text{mbox}county\}+1\{\text{mbox}{sponsor}_j=\text{mbox}county\})$$

$$S_4(\mathbf{x})=\sum _{i<j}{x}_{i,j}(1\{\text{mbox}{sponsor}_i=\text{mbox}federal\}+1\{\text{mbox}{sponsor}_j=\text{mbox}federal\})$$

$$S_5(\mathbf{x})=\sum _{i<j}{x}_{i,j}(1\{\text{mbox}{sponsor}_i=\text{mbox}private\}+1\{\text{mbox}{sponsor}_j=\text{mbox}private\})$$

$$S_6(\mathbf{x})=\sum _{i<j}{x}_{i,j}(1\{\text{mbox}{sponsor}_i=\text{mbox}state\}+1\{\text{mbox}{sponsor}_j=\text{mbox}state\})$$

$$S_7(\mathbf{x})=\sum _{i<j}{x}_{i,j}(\text{mbox}{command}_i+\text{mbox}{command}_j)$$

$$S_8(\mathbf{x})=\sum _{x_{i,j}\in L}G{D}_4,\text{leavevmode}\text{leavevmode}{S}_9(\mathbf{x})=\sum _{x_{i,j}\in NL}G{D}_4$$

$$\varphi (D)=\{\begin{array}{ll}0& 0\le D\le R\\ {\theta }_1-{\left(\frac{\sqrt{{\theta }_1}}{{\theta }_2-R}(D-{\theta }_2)\right)}^2& R<D\le D_1\\ 1+\frac{1}{({\theta }_3(D-D_2){)}^2}& D>D_1\end{array}$$

$$L(\theta |\mathbf{x})=\frac{{\lambda }^n\left[{\prod }_{i=1}^n\exp \left\{\min \left({\sum }_{i\ne j}\log (\varphi (D_{i,j})),1.2\right)\right\}\right]}{Z(\theta )},\text{leavevmode}\text{leavevmode}\theta =\{\lambda ,{\theta }_1,{\theta }_2,{\theta }_3\}.$$

$$\frac{f({\mathbf{y}}^{\prime }|{\theta }^{\prime },\mathbf{x})h(\mathbf{y}|\theta )}{h({\mathbf{y}}^{\prime }|{\theta }^{\prime })f(\mathbf{y}|\theta ,\mathbf{x})}=\frac{h({\mathbf{y}}^{\prime }|\hat{\theta })h(\mathbf{y}|\theta )}{h({\mathbf{y}}^{\prime }|{\theta }^{\prime })h(\mathbf{y}|\hat{\theta })}.$$

$$\alpha =\min \left\{1,\frac{p({\theta }^{\prime })h(\mathbf{x}|{\theta }^{\prime })h(\mathbf{y}|{\theta }_n)q({\theta }_n|{\theta }^{\prime })}{p({\theta }_n)h(\mathbf{x}|{\theta }_n)h(\mathbf{y}|{\theta }^{\prime })q({\theta }^{\prime }|{\theta }_n)}\right\},$$

$$L(\theta |\mathbf{x})=\frac{h(\mathbf{x}|\theta )}{\tilde{Z}(\theta )}c(\theta )\left[1+\sum _{n=1}^{\infty }k(\theta {)}^n\right],\text{leavevmode}\text{leavevmode}k(\theta )=1-c(\theta )\frac{Z(\theta )}{\tilde{Z}(\theta )},$$

$$\hat{L}(\theta |\mathbf{x})=\frac{h(\mathbf{x}|\theta )}{\tilde{Z}(\theta )}c(\theta )\left[1+\sum _{n=1}^{\infty }\prod _{i=1}^n{\hat{k}}_i(\theta )\right],\text{leavevmode}\text{leavevmode}{\hat{k}}_i(\theta )=1-c(\theta )\frac{{\hat{Z}}_i(\theta )}{\tilde{Z}(\theta )}.$$

$$E[S_{\tau }]=\sum _{j=1}^{\infty }{S}_{j}P(\tau =j)=1+\sum _{n=1}^{\infty }\prod _{i=1}^n{\hat{k}}_i(\theta ).$$

$${\hat{L}}_{\tau }(\theta |\mathbf{x})=\frac{h(\mathbf{x}|\theta )}{\tilde{Z}(\theta )}c(\theta )S_{\tau },\text{leavevmode}\text{leavevmode}{S}_{\tau }=\left[1+\sum _{n=1}^{\tau -1}\frac{{\prod }_{i=1}^n{\hat{k}}_i(\theta )}{{\prod }_{j=1}^{n-1}{q}_j}\right],$$

$$E[\tau ]=\sum _{x=1}^{\infty }xp(\tau =x)=\sum _{x=1}^{\infty }x{k}^{x(x-1)/2}(1-k^x)=\sum _{x=1}^{\infty }{k}^{x(x-1)/2}={\tilde{\tau }}^R$$

Full text

Keywords: Markov chain Monte Carlo, doubly intractable distributions, exponential random graph models, Markov point processes, importance sampling

1 Introduction

Markov chain Monte Carlo (MCMC) has been used to routinely carry out Bayesian inference for an enormous variety of complicated models (cf. Brooks et al.,, 2011). However, inference for Bayesian models with intractable normalizing functions, where the normalizing constant of the model is itself a function of parameters of interest, is still far from routine. There are many well known models that have intractable functions, for instance the class of exponential family random graph models and its variants (cf. Robins et al.,, 2007, Hunter and Handcock,, 2012) which are popular models for social networks. Non-Gaussian Markov random field models in spatial statistics (cf. Besag,, 1974, Hughes et al.,, 2011, for a review) also have intractable normalizing functions. In fact, it is worth noting that the Ising model (Lenz,, 1920, Ising,, 1925), which is a non-Gaussian Markov random field, appears in the landmark paper on the Metropolis algorithm (Metropolis et al.,, 1953). While the Metropolis algorithm provides an elegant way to simulate from this model for a given parameter value, the paper did not consider the more difficult problem of performing inference for this model.

Consider $h(\mathbf{x}|\theta)$, an unnormalized probability model for a random variable $\mathbf{x}\in\mathcal{X}$ given a parameter vector $\theta\in\Theta$, with a normalizing function $Z(\theta)=\int_{\mathcal{X}}h(\mathbf{x}|\theta)d\mathbf{x}$. Let $p(\theta)$ be the prior density for $\theta$. The likelihood function, $L(\theta|\mathbf{x})$ is $h(\mathbf{x}|\theta)/Z(\theta)$ and the posterior density of $\theta$ is

[TABLE]

The problem in constructing an MCMC algorithm stems from the fact that $Z(\theta)$ cannot be easily evaluated and the acceptance probability at each step of the Metropolis-Hastings (MH) algorithm (Metropolis et al.,, 1953, Hastings,, 1970) involves evaluating $Z(\theta)$ both at the current and proposed value of $\theta$.

There have been several proposals to address the intractable normalizing function problem in a maximum likelihood context. Besag, (1974) proposed the maximum pseudolikelihood estimate (MPLE) which maximizes the pseudolikelihood, a particular likelihood approximation that does not require $Z(\theta)$. MPLE does not, in general, define a likelihood function. Furthermore, MPLE may be a reasonable estimator when there is weak dependence among data points relative to the size of the data set but in other cases its performance is unsatisfactory. Younes, (1988) proposed a stochastic gradient algorithm to solve normal equations to find the maximum likelihood estimate (MLE) for models with intractable normalizing functions. However the step size and starting point must be selected carefully, otherwise the algorithm becomes slow and does not converge (Ibáñez and Simó,, 2003). Geyer and Thompson, (1992) propose MCMC-MLE which is based on maximizing a Monte Carlo approximation to the likelihood; this approximation is based on an importance sampling approximation of $Z(\theta)$. This is an elegant and theoretically justified algorithm. In practice, MCMC-MLE suffers from some of the usual challenges faced by importance sampling approaches, namely that for an accurate approximation to MLE, the initial value for the algorithm should be reasonably close to the MLE. Some of these issues may be partially addressed by umbrella sampling (Torrie and Valleau,, 1977, Geyer,, 2011), which involves sampling from mixtures of importance sampling distributions. However approximating standard errors can also be difficult in situations where analytical gradients for the unnormalized likelihood are unavailable and using bootstrap techniques may be computationally infeasible (cf. Goldstein et al.,, 2015). In such cases and also in situations where there is an interest in incorporating prior information about the parameters or avoiding model degeneracies (cf. Handcock,, 2003), Bayesian alternatives may be preferable. In this manuscript, we focus on Bayesian inference; we refer readers to Geyer, (2011) for a general overview and Hunter et al., (2012) for a review of recent MCMC-MLE methods for exponential random graph models.

Several MCMC algorithms have recently been proposed for Bayesian inference in the presence of intractable normalizing functions. These algorithms may be broadly classified into two general if somewhat overlapping categories: (1) algorithms where the introduction of a well chosen auxiliary variable results in the normalizing function (or a ratio of normalizing functions) canceling out in the Metropolis-Hastings acceptance probability, and (2) directly approximating the normalizing function (or a ratio of normalizing functions), and substituting the approximation into the acceptance probability. Here we will refer to the first as an auxiliary variable approach and the second as a likelihood approximation approach (see also a discussion in Liang et al.,, 2016). In addition, MCMC algorithms may also be classified as “asymptotically exact” or “asymptotically inexact”. For asymptotically exact algorithms the Markov chain’s stationary distribution is exactly equal to the desired posterior distribution. On the other hand, asymptotically inexact or “noisy” algorithms generate Markov chains without this property; even asymptotically the samples generated only follow the target distribution approximately.

In what follows we discuss several MCMC algorithms. We provide an explanation of the ideas underpinning each algorithm along with figures that summarize them and make it easier to see how they are related. We also discuss theoretical justifications and practical implementation issues. We carry out a comparative study by using three different examples: an Ising model, a social network model, and a spatial point process. We provide some guidance about potential advantages and disadvantages of each algorithm along with connections among them, providing some future avenues for research. The remainder of this paper is organized as follows. In Section 2 we discuss several auxiliary variable algorithms. In Section 3 we cover several likelihood approximation algorithms. In Section 4 we describe the application of the algorithms in the context of three different case studies and provide some insights based on our results. In particular, we discuss in detail the computational complexity of the algorithms in Section 5. We point out connections between the algorithms in Section 6 along with general guidelines and recommendations based on our study. We conclude with a summary in Section 7.

2 Auxiliary Variable Approaches

In this section, we review several auxiliary variable approaches. Here, the target distribution includes both the parameter of interest as well as an auxiliary variable. By a clever choice of the auxiliary variable proposal the intractable functions get cancelled in the acceptance probability of the Metropolis-Hastings algorithm. What distinguishes the different auxiliary variable algorithms from each other is how the auxiliary variable is sampled.

2.1 Auxiliary Variable MCMC

Møller et al., (2006) introduce an auxiliary variable $\mathbf{y}$ with the conditional density $f(\mathbf{y}|\theta,\mathbf{x})$ so that the intractable terms are cancelled in the acceptance probability of the Metropolis-Hastings algorithm. Suppose the original target density is $\pi(\theta|\mathbf{x})\propto p(\theta)h(\mathbf{x}|\theta)/Z(\theta)$. Then the augmented target density is $\pi(\theta,\mathbf{y}|\mathbf{x})\propto f(\mathbf{y}|\theta,\mathbf{x})p(\theta)h(\mathbf{x}|\theta)/Z(\theta)$ whose marginal density becomes $\int_{\mathcal{X}}\pi(\theta,\mathbf{y}|\mathbf{x})d\mathbf{y}=\pi(\theta|\mathbf{x})$, the original target density. Now consider a joint proposal density for $\{\theta,\mathbf{y}\}$ updates which can be factorized as $q(\theta^{\prime},\mathbf{y}^{\prime}|\theta,\mathbf{y})=q(\mathbf{y}^{\prime}|\theta^{\prime})q(\theta^{\prime}|\theta)$. Møller et al., (2006) take the proposal for the auxiliary variable $q(\mathbf{y}^{\prime}|\theta^{\prime})$ as $h(\mathbf{y}^{\prime}|\theta^{\prime})/Z(\theta^{\prime})$ which is the same as the model for the data $\mathbf{x}$ given the parameter value $\theta^{\prime}$ proposed from $q(\theta^{\prime}|\theta)$. The resulting algorithm for the above joint density with this proposal density has Metropolis-Hastings acceptance probability

[TABLE]

The resulting Metropolis-Hastings acceptance probability (2) does not contain the normalizing functions because they get cancelled out. Since $\int_{\mathcal{X}}\pi(\theta,\mathbf{y}|\mathbf{x})d\mathbf{y}=\pi(\theta|\mathbf{x})$, the marginal $\theta$ samples follow the original target distribution. We will henceforth use AVM to refer to this auxiliary variable MCMC algorithm. We summarize the algorithm in Figure 1. We should note that updating $\mathbf{y}$ requires drawing a sample with exact distribution $h(\cdot|\theta^{\prime})/Z(\theta^{\prime})$. For some models this can be achieved via perfect sampling (Propp and Wilson,, 1996), a clever method that uses bounding Markov chains to construct a sampler where the draws are exactly (not just asymptotically) from the target distribution. Although perfect sampling is possible for some models, for instance certain Markov random field (MRF) models, it is not trivial to construct a perfect sampler in general. This is a major practical limitation of the AVM.

Components to be tuned: Møller et al., (2006) report that the choice of $f(\mathbf{y}|\theta,\mathbf{x})$ impacts the mixing of the chain. Suppose $Z(\theta)$ is known so that we can choose $f(\mathbf{y}|\theta,\mathbf{x})=h(\mathbf{y}|\theta)/Z(\theta)$. Then it is easily seen that (2) is identical to the acceptance probability of the Metropolis-Hastings algorithm with the stationary density $\pi(\theta|\mathbf{x})$, which implies that this chain has the same convergence properties as the Markov chain where the normalizing function is known and the same proposal $q(\theta^{\prime}|\theta)$ is used for $\theta$. Therefore, a reasonable choice of conditional density for the auxiliary variable is $f(\mathbf{y}|\theta,\mathbf{x})$ that approximates $h(\mathbf{y}|\theta)/Z(\theta)$. A simple choice is $f(y|\theta,\mathbf{x})=h(\mathbf{y}|\widehat{\theta})/Z(\widehat{\theta})$, where $\widehat{\theta}$ may be an approximation to the MLE. $\widehat{\theta}$ should be predetermined before implementing the AVM algorithm. For example, the maximum pseudolikelihood estimate (MPLE) proposed by Besag, (1974) may be an option though for some problems the MPLE may be a poor approximation to the MLE. Then $f(\mathbf{y}^{\prime}|\theta^{\prime},\mathbf{x})/f(\mathbf{y}|\theta,\mathbf{x})=h(\mathbf{y}^{\prime}|\widehat{\theta})/h(\mathbf{y}|\widehat{\theta})$, which makes it possible to calculate (2) because there are no intractable terms. Another possible choice of $f(\mathbf{y}|\theta,\mathbf{x})$ is using a normalizable density without $Z(\theta)$. AVM mixes better as $f(\mathbf{y}|\theta,\mathbf{x})$ more closely resembles $h(\mathbf{y}|\theta)/Z(\theta)$.

Theoretical justification: The Markov chain satisfies the detailed balance condition if $\pi(\theta|\mathbf{x})T(\theta^{\prime}|\theta)=\pi(\theta^{\prime}|\mathbf{x})T(\theta|\theta^{\prime})$, where $T(\theta^{\prime}|\theta)$ is the transition kernel of the Markov chain with the target density $\pi(\theta|\mathbf{x})$. Since the Markov Chain satisfies detailed balance with respect to the augmented target distribution, the marginal distribution of which is the posterior distribution of $\theta$, AVM is an asymptotically exact MCMC algorithm. However, to achieve the detailed balance condition, we need to sample $\mathbf{y}$ exactly from the likelihood function $h(\mathbf{y}|\theta^{\prime})/Z(\theta^{\prime})$.

2.2 The Exchange Algorithm

Appearing almost simultaneously with Møller et al., (2006), the exchange algorithm (Murray et al.,, 2006) also constructs an augmented target distribution and updates the augmented state via the Metropolis-Hastings algorithm. Consider the auxiliary variable $\mathbf{y}$ which follows $h(\mathbf{y}|\theta^{\prime})/Z(\theta^{\prime})$ and conditional density of $\theta^{\prime}$ for given $\theta$ as $q(\theta^{\prime}|\theta)$ where $\theta$ is the parameter setting of data $\mathbf{x}$. Then the augmented joint density is

[TABLE]

For this augmented density, $\{\theta^{\prime},\mathbf{y}\}$ is updated through block-Gibbs samplers; $\theta^{\prime}$ is generated from the proposal $q(\cdot|\theta)$ and the auxiliary variable $\mathbf{y}$ is generated from $h(\cdot|\theta^{\prime})/Z(\theta^{\prime})$. The first two arrows in Figure 2 correspond to the update of $\{\theta^{\prime},\mathbf{y}\}$. Then $\theta$ is updated through exchanging parameter settings. Let $\{\theta,\theta^{\prime}\}$ be the current parameter settings for $\{\mathbf{x},\mathbf{y}\}$. Consider a swapping proposal $s(\{\theta^{*},\theta^{\prime*}\}|\{\theta,\theta^{\prime}\})=\delta(\theta^{*}-\theta^{\prime})\delta(\theta^{\prime*}-\theta)$, where $\delta$ denotes the Dirac delta function. After swapping is proposed, data $\mathbf{x}$ follows $\theta^{\prime}$ instead of $\theta$ and the auxiliary variable $\mathbf{y}$ follows $\theta$ instead of $\theta^{\prime}$. The symmetric swapping proposal results in the acceptance probability ($\alpha$),

[TABLE]

with intractable terms cancelled. Since $\int_{\mathcal{X}}\int_{\Theta}\pi(\theta,\theta^{\prime},\mathbf{y}|\mathbf{x})d\theta^{\prime}d\mathbf{y}=\pi(\theta|\mathbf{x})$, the marginal $\theta$ samples follow the original target distribution, as shown in Figure 2. This idea of swapping parameter settings is connected to the parallel tempering algorithm (Geyer,, 1991). Parallel tempering swaps the states while preserving joint distributions in the product space. The exchange algorithm also swaps parameter settings between $\mathbf{x}$ and $\mathbf{y}$, while preserving the augmented distribution. Since we do not need to keep $\{\theta^{\prime},\mathbf{y}\}$ in the actual implementation of the algorithm, the exchange algorithm may be simply written as Algorithm 2.

The exchange algorithm can also be extended through “bridging” (Murray et al.,, 2006), Algorithm 3. The idea may be summarized as follows. In the acceptance probability of Algorithm 2, one may think of the ratio $h(\mathbf{x}|\theta^{\prime})/h(\mathbf{x}|\theta)$ as measuring whether the proposed $\theta^{\prime}$ explains the data $\mathbf{x}$ better than current $\theta$ or not. Also, $h(\mathbf{y}|\theta)/h(\mathbf{y}|\theta^{\prime})$ represents how well the auxiliary variable $\mathbf{y}$ is supported under $\theta$ versus $\theta^{\prime}$. Therefore, even if $\theta^{\prime}$ is a much better candidate than $\theta$ under the data $\mathbf{x}$, swapping can be rejected if the auxiliary variable $\mathbf{y}$ is improbable under $\theta$, which can lead to slow mixing of the chain. To improve this, annealed importance sampling (AIS) (Neal,, 1996, 2001) can be combined with the exchange algorithm. Instead of sampling a single auxiliary variable, a series of auxiliary variables $\{\mathbf{y}_{0},\mathbf{y}_{1},...,\mathbf{y}_{K}\}$ are sampled from the intermediate densities between $h(\mathbf{y}|\theta^{\prime})$ and $h(\mathbf{y}|\theta)$. Intermediate densities can be constructed as

[TABLE]

After proposing $\theta^{\prime}$, generate $\mathbf{y}_{0}$ from $h(\mathbf{y}|\theta^{\prime})/Z(\theta^{\prime})$. Then generate each $\mathbf{y}_{k}$ from $T_{k}(\mathbf{y}_{k}|\mathbf{y}_{k-1})$, the Metropolis-Hastings (MH) transition kernel whose stationary density is $f_{k}(\mathbf{y}|\theta,\theta^{\prime})$. This bridging step helps to generate $\mathbf{y}_{K}$ that is more probable under $\theta$.

Components to be tuned: The basic exchange algorithm does not require any tuning besides the usual tuning of the proposal for $\theta$. For the extended version, there are many options for bridging but even here (5) provides an automated schedule. Larger $K$ in (5) can lead to better mixing of the chain at the expense of computing time. Effective sample size (ESS) (Kass et al.,, 1998) provides a rough diagnostic for how well the chain is mixing. Therefore, effective sample size per unit time (ESS/T) can be used to determine $K$ that provides a compromise between mixing and computational cost (Murray,, 2007).

Theoretical justification: Both the exchange algorithm and the extended exchange algorithm with bridging are asymptotically exact MCMC in that constructed Markov chains satisfy detailed balance condition with respect to the target distribution. However just like AVM, the exchange algorithm is of limited applicability because it requires perfect sampling to generate the auxiliary variable.

2.3 The Double Metropolis-Hastings Sampler

Both AVM and the exchange algorithm require perfect sampling which can be very expensive or impossible for complicated probability models. Liang, (2010) replaces perfect sampling for the auxiliary variable with a standard Metropolis-Hastings algorithm; the last state of the resulting Markov chain is then treated like a draw from the perfect sampler. This is called Double Metropolis-Hastings (DMH) because a Metropolis-Hastings algorithm is used within another Metropolis-Hastings algorithm. As in Figure 3, one is the “outer sampler” to generate $\theta$ draws while the other is the “inner sampler” used to generate the auxiliary variable $\mathbf{y}$.

Define $T^{m}_{\theta^{\prime}}(\mathbf{y}|\mathbf{x})$ as $m$-MH updates from $\mathbf{x}$ to $\mathbf{y}$ under $\theta^{\prime}$ whose stationary density is $h(\mathbf{y}|\theta^{\prime})/Z(\theta^{\prime})$; this satisfies the detailed balanced condition

[TABLE]

Plugging this result into (4) yields

[TABLE]

which is the acceptance probability of the DMH algorithm. DMH approximates perfect sampling through $m$ steps of an (inner) Metropolis-Hastings algorithm. A similar approximate approach which replaces perfect sampling with MH updates for the social network models is discussed in Caimo and Friel, (2011). As we show later, DMH is simpler to implement and computationally more efficient relative to all other algorithms. However, its computational efficiency directly depends on the efficiency of the inner sampler. If the inner sampler update $T^{m}_{\theta^{\prime}}(\mathbf{y}|\mathbf{x})$ is expensive or needs large $m$ to get probable auxiliary variable $\mathbf{y}$, then DMH becomes computationally expensive. Since $m$ is usually proportional to the dimension of the data $\mathbf{x}$, the inner sampler is the main computational bottleneck for large data.

Components to be tuned: Users need to decide the number of MH updates $m$ to generate $\mathbf{y}$. There are numerous stopping rules for MCMC if computing expectations is the goal (cf. Flegal et al.,, 2008, Gong and Flegal,, 2015). However, here determining convergence to the stationary distribution is of interest (Jones and Hobert,, 2001, Rosenthal,, 1995), which is a very challenging problem in general. A simple heuristic for determining $m$ is to set $m$ to be proportional to the size of the data, for instance $m=5n$, where $n$ is the size of the data. We recommend that DMH should be run for larger $m$ values, for instance $m=10n$ to confirm that the results do not change much as $m$ is increased.

Theoretical justification: DMH is an asymptotically inexact algorithm because the detailed balance condition does not hold for the outer sampler unless the inner sampler length ($m$) approaches infinity; in practice the inner sampler length is of course finite.

2.4 Adaptive Exchange Algorithm

In order to address the asymptotic inexactness of DMH, Liang et al., (2016) propose the adaptive exchange algorithm (AEX). AEX runs two chains simultaneously: the auxiliary chain and the target chain (see Figure 4). The auxiliary chain generates a sample from a family of distributions, $\{h(\mathbf{x}|\theta^{(1)})/Z(\theta^{(1)}),...,h(\mathbf{x}|\theta^{(d)})/Z(\theta^{(d)})\}$, where $\{\theta^{(1)},...,\theta^{(d)}\}$ are pre-specified “particles” covering a parameter space. The generated sample is stored at each iteration. Then the target chain generates a posterior sample from the $\pi(\theta|\mathbf{x})$ via the exchange algorithm. The difference is that the auxiliary variable $\mathbf{y}$ is sampled from the cumulative samples in the auxiliary chain until current iteration (resampling). With increasing iterations, cumulative samples from the auxiliary chain grow, so that the resampling of $\mathbf{y}$ in the target chain becomes identical to exact sampling of $\mathbf{y}$. Then similar to the exchange algorithm, the marginal density of $\theta$ converges to the target $\pi(\theta|\mathbf{x})$. AEX is therefore an asymptotically exact algorithm that does not require perfect sampling. We begin with some notation for the $n$th iteration of AEX.

• •

Particles: $\{\theta^{(1)},...,\theta^{(d)}\}$, each $\theta^{(i)}\in\Theta$. These particles are fixed through the algorithm.

• •

Particle index: $I_{n}\in\{1,...,d\}$ returns index of a chosen particle at $n$th iteration.

• •

Auxiliary data: $\mathbf{x}_{n}\in\mathcal{X}$ is an approximate sample from the probability model at selected particle $\theta^{(I_{n})}$. That is, $\mathbf{x}_{n}\sim h(\cdot|\theta^{(I_{n})})/Z(\theta^{(I_{n})})$ .

• •

Normalizing function approximation at each particle: $\mathbf{w}_{n}=\{w^{(1)}_{n},...,w^{(d)}_{n}\}\in W$. For $i=1,...,d$, as $n$ gets large $w^{(i)}_{n}$ converges to $Z(\theta^{(i)})$ via the stochastic approximation algorithm (Liang et al.,, 2007).

• •

Cumulative information: $H_{n}=\cup_{j=1}^{n}\{I_{j},\mathbf{x}_{j},\mathbf{w}_{j}\}$ is necessary for constructing the algorithm. Therefore, $\{I_{n},\mathbf{x}_{n},\mathbf{w}_{n}\}$ should be stored at each iteration.

• •

Posterior sample: $\theta_{n}\in\Theta$

AEX updates a non-Markovian stochastic process $\{I_{n+1},\mathbf{x}_{n+1},\mathbf{w}_{n+1},\theta_{n+1}\}\in\{1,...,d\}\times\mathcal{X}\times R^{d}\times\Theta$ in the $(n+1)$st iteration. For each iteration, the auxiliary chain generates a sample $\mathbf{x}_{n+1}$ from the mixture density, $(1/d)\sum_{i=1}^{d}h(\mathbf{x}|\theta^{(i)})/Z(\theta^{(i)})$ through stochastic approximation Monte Carlo (SAMC) (Liang et al.,, 2007). Particles are $d$-number of points in the grid over a parameter space $\Theta$. Particles are chosen so that they can cover the important region of the parameter space and the sample spaces of neighboring distributions $\{h(\mathbf{x}|\theta^{(i)})/Z(\theta^{(i)})\}$ overlap each other. A clever choice of particles can improve mixing of the auxiliary chain and is important for AEX to be a practical algorithm. We refer the reader to the supplementary material for strategies for choosing particles. For the $(n+1)$st update of AEX, cumulative information $H_{n+1}$ is constructed through the auxiliary chain. For $(n+1)$st iteration of the algorithm, $\{I_{n+1},\mathbf{x}_{n+1},\mathbf{w}_{n+1}\}$ is updated for given $\{I_{n},\mathbf{x}_{n},\mathbf{w}_{n}\}$ using SAMC and then added to $H_{n}$ to construct $H_{n+1}$. Hence, the accumulated information $(H_{n+1})$ becomes larger with each iteration of AEX. See Algorithm 5 for details.

Then the target chain generates $\theta_{n+1}$ from the posterior using the exchange algorithm. For proposed $\theta^{\prime}$, the auxiliary variable $\mathbf{y}$ is sampled from $\{\mathbf{x}_{1},...,\mathbf{x}_{n+1}\}$ through a dynamic importance sampling procedure. The resampling probability is

[TABLE]

In (8), $h(\mathbf{x}_{j}|\theta^{(I_{j})})/w^{(I_{j})}_{j}$ is the importance function at $j$th iteration. Since the importance function changes as $I_{j}$ varies, it is called a dynamic importance function. Here, $\mathbf{x}_{l}$’s more likely under $\theta^{\prime}$ are more likely to be sampled. Liang et al., (2016) show that the distribution of resampled $\mathbf{y}$ converges to $h(\mathbf{y}|\theta^{\prime})/Z(\theta)$ as $n$ increases. Since the proposal for $\mathbf{y}$ in (8) changes with each iteration, AEX appears at first to be an adaptive MCMC algorithm. However, unlike basic adaptive MCMC algorithms the target distribution also varies as $H_{n+1}$ grows.

In Step 2(a) of Algorithm 5, $\{a_{n}\}$ is called the gain factor, which is the step size of update for the abundance factor $\mathbf{w}_{n}$. A large gain factor can lead to faster movement around all $\{h(\mathbf{x}|\theta^{(i)})/Z(\theta^{(i)})\}$. Since $a_{n}=n_{0}/\max(n_{0},n)$ becomes smaller as $n$ increases, movement between particles becomes slower with increasing iterations. Liang et al., (2007) suggest larger $n_{0}$ for complicated problems so that makes it fast to move around all state space. Step 2(b) illustrates varying truncation (Jin and Liang,, 2013), which can help the convergence of $\mathbf{w}_{n}$ regardless of starting point. For $\mathbf{w}_{n}\in W$, consider the sequence of set $\{K_{s}\}\in W$, which satisfies $\cup_{s\geq 0}K_{s}=W$ and $K_{s}\subset K_{s+1}^{\mathrm{o}}$, where $K_{s+1}^{\mathrm{o}}$ denotes the interior of set $K_{s+1}$. A truncation function $\mathrm{T}$ maps a point in $\mathcal{X}\times W$ to a random point in $\mathcal{X}_{0}\times K_{0}$ for $\mathcal{X}_{0}\subset\mathcal{X}$, and $\sigma_{n}$ is the number of truncations that occurred by the $n$th iteration ($\sigma_{0}=0$). From such truncation, when an updated $\mathbf{w}_{n+1/2}$ is outside the interesting target region $(K_{s})$, it is reinitialized with a smaller gain factor and expanded target region $(K_{s+1}\supset K_{s})$. This varying truncation can help SAMC to find an appropriate step size $\{a_{n}\}$ and starting point of abundance factor automatically (Jin and Liang,, 2013). The reader may understand such varying truncation method as a safeguard for when the abundance factor becomes degenerate.

In practice, Liang et al., (2016) suggest that only the auxiliary chain (Algorithm 5 Part 1) is implemented $N_{1}$ iterations at first to construct the database, and then the auxiliary and target chain can be run simultaneously. In detail, through $N_{1}$ preliminary iterations of the auxiliary chain, information: $\cup_{j=1}^{N_{1}}\{I_{j},\mathbf{x}_{j},\mathbf{w}_{j}\}$ is constructed. Using this information as initial $H_{0}=\cup_{j=1}^{N_{1}}\{I_{j},\mathbf{x}_{j},\mathbf{w}_{j}\}$, we can implement the entire Algorithm 5. In this case, the size of cumulative information become $|H_{0}|=N_{1}$, and $|H_{n+1}|=N_{1}+n+1$ (without preliminary step, $|H_{0}|=0$, and $|H_{n+1}|=n+1$). The reason of conducting preliminary step is because the early performance of the target chain can be affected by initial starting values of the auxiliary chain. This preliminary step can construct an initial database $H_{0}$ and improves the mixing of the target chain.

Although we classify AEX as an auxiliary variable approach, AEX also shares characteristics with likelihood approximation approaches. The abundance factor $\{w^{(i)}_{j}\}$ converges to $\{Z(\theta^{(i)})\}$ for each particle, which is guaranteed by SAMC. This implies that likelihood approximation approaches are applied to estimate $Z(\theta^{(i)})$ directly. Then the auxiliary variable is resampled to cancel out $Z(\theta)$ in the target chain. A strength of AEX is its broad applicability, while still remaining asymptotically exact. It is reasonably efficient computationally. However because of its adaptive structure, there can be serious memory issues. When there are low-dimensional sufficient statistics, as is typically the case with exponential family models, only the sufficient statistics of $\mathbf{x}_{n+1}$ need to be stored in Step 3 of the Algorithm 5. However, in the absence of such sufficient statistics, $\mathbf{x}_{n+1}$ itself need to be stored at each step. Furthermore, without sufficient statistics, resampling probability calculations in Step 5 becomes expensive because $h(\mathbf{x}_{j}|\theta^{\prime})$ should be recalculated; when there are sufficient statistics, one can simply take the product of $\theta^{\prime}$ and the sufficient statistic.

Components to be tuned: We provide details about choosing particles (number of particles, strategies for choosing particles) in the supplementary material. In the auxiliary chain, the number of neighbors for updating $\theta^{(I_{n})}$, the number of preliminary runs for the auxiliary chain ($N_{1}$), the number of MH updates ($m$), gain factor component ($n_{0}$), and the target region ($K_{s}$) of abundance factor should be tuned. For particles, neighbor is defined using the distance from each other. We can define the closest 20 points as neighboring particles, and this choice appears to work well both in simulated examples in this manuscript as well as the examples in Liang et al., (2016). For updating $\mathbf{x}_{n}$, $m=1$ cycle of MH updates are enough, because we don’t need to generate independent samples. For complicated problems Liang et al., (2016) recommend larger values for $N_{1}$ and $n_{0}$. For the simulation examples, Liang et al., (2016) use around $N_{1}=10,000d$, $n_{0}=25,000$ where $d$ is the number of particles. We used similar values in our simulated examples. If the preliminary auxiliary chain does not appear to converge, larger $n_{0}$ and $N_{1}$ should be used. The convergence of the auxiliary chain can be checked by comparing sampling frequencies $v(i)/N_{1}$ with the target probabilities $1/d$, where $v(i)$ is the number of visitations at each particle. If $v(i)/N_{1}\approx 1/d$ for each particle, then the preliminary run of auxiliary chain can be diagnosed as having converged. For the choice of compact subset of $W$, Liang et al., (2016) set $K_{s}=[0,10^{100+10s}]^{d}$ and $\mathcal{X}_{0}=\mathcal{X}$ in their examples. Detailed simulation settings about the social network model also can be checked in Jin et al., (2013).

Theoretical justification: The AEX is an asymptotically exact algorithm and the ergodic average satisfies the Weak Law of Large Numbers. Since AEX is adaptive in both proposal and target, conventional theory for adaptive MCMC does not directly apply. Liang et al., (2016) extend a result in Roberts and Rosenthal, (2007)[Theorem 1], which shows the ergodicity of the adaptive chain with the same target distribution. Although the original proof (Liang et al.,, 2016) assumes compactness of the parameter space $\Theta$, Jin et al., (2013) extend the results without a compactness assumption for $\Theta$. The details of the assumptions for the proof(s) are provided across several results in Jin et al., (2013)[Lemma 3.1, 3.2 and Theorem 3.1]. In order to make it easy for readers to understand the assumptions from a practical point of view, we have attempted to distill the key assumptions as follows: (1) both $\mathcal{X}$ and $W$ are compact, (2) $h(\mathbf{x}|\theta)$ is bounded away from [math] and $\infty$, (3) $\lim_{n\rightarrow\infty}a_{n}=0,\leavevmode\nobreak\ \sum_{n=1}^{\infty}a_{n}=\infty,\leavevmode\nobreak\ \sum_{n=1}^{\infty}a^{\eta}_{n}<\infty$ for some $\eta\in(1,2]$, (4) Doeblin condition holds for kernel of auxiliary chain. Some of these conditions are easily verified in practice. The sample space $\mathcal{X}$ for data is compact for finite lattice problem and point process with finite number of points on bounded domain. Sample space $W$ for abundance factor is also compact, if we set $K_{s}$ as suggested in the components to be tuned above. The second assumption is also satisfied for realistic parameter settings. The third is a technical assumption for SAMC, and is satisfied if we construct $a_{n}=n_{0}/\max(n_{0},n)$. The sufficient condition for the last assumption is local positiveness of $m$-MH updates in Step 1(b) of Algorithm 5, which means $\exists\leavevmode\nobreak\ \delta>0,\epsilon>0$ such that, $\forall\mathbf{x}\in\mathcal{X},|\mathbf{x}-\mathbf{y}|\leq\delta$ implies $T(\mathbf{y}|\mathbf{x})\geq\epsilon$. This is satisfied for a simple Gibbs sampler, birth-death sampler in point process, and tie-no-tie sampler in social network models. Therefore, broadly speaking, the assumptions hold in most problems. We point the readers to the manuscript in order to find the complete set of assumptions for these results.

3 Likelihood Approximation Methods

While the auxiliary variable approaches we have discussed so far avoid approximating $Z(\theta)$, the approaches based on likelihood approximations we discuss in this section directly approximate $Z(\theta)$ through Monte Carlo and substitute the approximation into the Metropolis-Hastings acceptance probability.

3.1 Atchade, Lartillot and Robert’s Adaptive MCMC

Atchade et al., (2008) provide an adaptive MCMC algorithm, henceforth the ALR algorithm, which approximates $Z(\theta)$ through importance sampling in the acceptance probability. ALR constructs a non-Markovian stochastic process and approximates $Z(\theta)$ at every iteration using the entire sample path of the process. ALR is an asymptotically exact algorithm that does not require perfect sampling.

The ALR algorithm utilizes importance sampling ideas developed in MCMC-MLE (Geyer and Thompson,, 1992). Consider $\{\mathbf{x}_{1},...,\mathbf{x}_{n}\}$, samples from $h(\mathbf{x}|\theta^{(0)})/Z(\theta^{(0)})$ for some $\theta^{(0)}$. Then $Z(\theta)/Z(\theta^{(0)})$ is an expectation that can be approximated by importance sampling. However the approximation might be poor if $\theta$ is far from $\theta^{(0)}$. Atchade et al., (2008) introduces multiple particles, $\{\theta^{(1)},...,\theta^{(d)}\}$, and a linear combination of importance sampling estimates for $\{Z(\theta)/Z(\theta^{(1)}),...,Z(\theta)/Z(\theta^{(d)})\}$ approximates $Z(\theta)$. Larger weights are assigned to the estimate of $Z(\theta)/Z(\theta^{(i)})$ when $\theta$ and $\theta^{(i)}$ are closer to each other. This idea of sampling from a mixture is related to umbrella sampling (Torrie and Valleau,, 1977) which provides a more robust approximation than a single point importance sampling estimate, as long as one of the $\theta^{(i)}$s is close to $\theta$. We begin with some notation for the $n$th iteration of ALR.

• •

Particles: $\{\theta^{(1)},...,\theta^{(d)}\}$, each $\theta^{(i)}\in\Theta$. These particles are fixed through the algorithm.

• •

Particle index: $I_{n}\in\{1,...,d\}$ returns the index of the selected particle at $n$th iteration.

• •

Auxiliary data: $\mathbf{x}_{n}\in\mathcal{X}$ is an approximate sample from the probability model at a selected particle $\theta^{(I_{n})}$. i.e. $\mathbf{x}_{n}\sim h(\cdot|\theta^{(I_{n})})/Z(\theta^{(I_{n})})$ .

• •

Normalizing function approximation at each particle: $\mathbf{c}_{n}=\{c^{(1)}_{n},...,c^{(d)}_{n}\}\in R^{d}$. For $i=1,...,d$, as $n$ gets large $c^{(i)}_{n}$ converges to $\log{Z(\theta^{(i)})}$ by stochastic approximation (Wang and Landau,, 2001, Atchade and Liu,, 2004).

• •

Cumulative information: $H_{n}=\cup_{j=1}^{n}\{I_{j},\mathbf{x}_{j}\}$ is necessary for constructing the algorithm. Therefore, $\{I_{n},\mathbf{x}_{n}\}$ should be stored at each iteration.

• •

Posterior sample: $\theta_{n}\in\Theta$.

By extending a result in Wang and Landau, (2001), Atchade and Liu, (2004), the ALR algorithm constructs a non-Markovian stochastic process $\{\mathbf{x}_{n+1},I_{n+1},\mathbf{c}_{n+1},\theta_{n+1}\}\in\mathcal{X}\times\{1,...,d\}\times R^{d}\times\Theta$. $Z(\theta)$ is approximated through multiple importance sampling and plugged into the acceptance probability to generate samples from $\pi(\theta|\mathbf{x})$. A stochastic approximation for estimating marginal densities, similar to the approximation used in ALR, is also described in Liang, (2007). The decomposition of $Z(\theta)$ is

[TABLE]

This suggests the following Monte Carlo approximation of $Z(\theta)$,

[TABLE]

a weighted sum of the importance sampling estimate $Z(\theta^{(i)})$ at each particle. A “similarity kernel” $k(\theta,\theta^{(i)})$ measures the distance between $\theta$ and $\theta^{(i)}$ and assigns larger weights to particles closer to $\theta$. Similar to AEX, particles are pre-specified before starting the algorithm and are chosen so that they can cover the parameter space. In the examples in Atchade et al., (2008), particles are collected from the stochastic approximation (Younes,, 1988), though particles could also be chosen as in AEX (see supplementary material for details). Algorithm 7 provides ALR details.

We split the ALR algorithms into two parts (see Figure 5. The purpose of Part 1 of the Algorithm 7 is to obtain a reasonable starting value for the normalizing function approximation without the memory and computational costs associated with updating the entire process $\{\mathbf{x}_{n},I_{n},\mathbf{c}_{n},\theta_{n}\}$. The step size $\{\gamma_{n}\}$ plays an important role in the convergence of $\mathbf{c}_{n}$ (analogous to how the gain factor $a_{n}$ is important for the abundance factor $\mathbf{w}_{n}$ in AEX). If $\gamma_{n}$ is too small, it takes a long time for convergence. If $\gamma_{n}$ is too large it can lead to large variance of $\mathbf{c}_{n}$. Therefore in the beginning of the algorithm, $\gamma_{n}$ is set to be a large value so that it can move around the space fast; it slowly becomes smaller after moving around the state space. In Algorithm 7 (Part 1), the number of visits to each particle is recorded through $\mathbf{v}$. If $\mathbf{v}$ is close to uniform distribution ($\{I_{n}\}$ has visited $\{1,...,d\}$ about equally), $\gamma_{n}$ becomes smaller. If step size $\gamma_{n}$ is smaller than some convergence criteria $\epsilon_{1}$, we can assume $\mathbf{c}_{n}$ has converged to a reasonable value.

The strengths and weakness of ALR are similar to those of AEX: ALR is asymptotically exact without requiring perfect sampling. However, the computational and memory costs can be large. Without low-dimensional sufficient statistics, the entire $\mathbf{x}_{n}$ chain needs to be stored with each iteration to calculate $\widehat{Z}_{n+1}(\theta)$ in Step 4 of the Algorithm 8. Furthermore, without sufficient statistics, calculations in Step 3 and Step 4 become expensive, because $h(\mathbf{x}_{k}|\theta),h(\mathbf{x}_{k}|\theta^{(i)})$ need to be recalculated; with sufficient statistics, one can simply take the product of $\theta$ or $\theta^{(i)}$ and the sufficient statistic of $\mathbf{x}_{k}$.

Components to be tuned: In this algorithm, step size ($\gamma_{n}$), the convergence checking components ($\epsilon_{1},\epsilon_{2}$), number of MH updates ($m$), number of particles ($d$) and kernel $k(\theta,\theta^{(i)})$ need to be tuned. Atchade et al., (2008) set initial $\gamma_{0}$ as 1, $\epsilon_{1}=0.001,\epsilon_{2}=0.2$. Under these settings, consider the sequence of bounded stopping times $0=\tau_{0}<\tau_{1}<...<\tau_{10}=\tau$, where $\tau$ is the total number of iteration in Algorithm 7 (Part 1). Until $\tau_{1}$, initial $\gamma_{0}$ is used as step size. For this step size, $\tau_{1}$ is the stopping time until each particle has been visited equally (i.e. $\max_{i=1,...,d}|v(i)-\bar{v}|\leq\epsilon_{2}\bar{v}$). Then $\gamma_{\tau_{1}+1}$ becomes $\gamma_{0}/2=1/2$ and is kept until $\tau_{2}$. This is repeated until $\tau_{10}$ where $\gamma_{\tau_{10}+1}=1/2^{10}<0.001=\epsilon_{1}$. Likewise, step size is controlled to hasten convergence of $\mathbf{c}_{n}$. Once $\mathbf{c}_{n}$ appears to have converged, Algorithm 8 (Part 2) is implemented. $\gamma_{n+1}$ is updated as the deterministic sequence $0.001/(n+1)^{0.7}$ in Algorithm 8 (Part 2). Similar to AEX, $m=1$ cycle of MH updates are enough in practice. For the number of particles, $d=100p$ appears to work well in practice to cover the $p$-dimensional parameter space. Although various choice of kernel $k(\theta,\theta^{(i)})$ is possible, uniform kernel with bandwidth $h$ is used in this manuscript. This kernel gives $1/h$ weights for the $h$ closest particles and gives 0 weights for others. Bandwidth should be determined by trials and errors and we used $h=20$.

Theoretical justification: The ALR is an asymptotically exact algorithm and the ergodic average from the chain satisfies the Strong Law of Large Numbers. Atchade et al., (2008) provide a proof that $\{\theta_{n}\}$ from the generated process converges to the target distribution exactly. Three assumptions are required: (1) $h(\mathbf{x}|\theta)$ is bounded away from 0 and $\infty$, (2) $q^{n_{0}}(\theta^{\prime}|\theta)\geq\epsilon$ for all $\theta,\theta^{\prime}\in\Theta$ where $q$ is proposal density and $n_{0}\geq 1$ is an integer, (3) $\{\gamma_{n}\}$ is random and adaptively updated based on the previous generated process, which satisfies $\sum\gamma_{n}=\infty$ and $\sum\gamma^{2}_{n}<\infty$ almost surely. The first assumption holds for finite $\mathcal{X}$ and realistic parameter settings. The second assumption also holds for most symmetric kernels. The third assumption is technical, and taken from Wang and Landau, (2001). Because we can typically control the sequence $\{\gamma_{n}\}$, this assumption is also generally easy to satisfy. Therefore, the assumptions appear to hold in most cases in practice.

3.2 Pseudo-marginal MCMC

Suppose $\widehat{L}(\theta|\mathbf{x})$ is a positive and unbiased Monte Carlo approximation of $L(\theta|\mathbf{x})$. An MCMC algorithm that uses $\widehat{L}(\theta|\mathbf{x})$ in place of $L(\theta|\mathbf{x})$ is called pseudo-marginal MCMC (Beaumont,, 2003, Andrieu and Roberts,, 2009) and its stationary distribution is equal to the desired target posterior distribution.

The pseudo-marginal framework in our context requires an unbiased approximation for $1/Z(\theta)$. An unbiased approximation of $Z(\theta)$ can be obtained using importance sampling estimate $\widehat{Z}(\theta)$ via MCMC samples from the likelihood. However, $1/\widehat{Z}(\theta)$ is a consistent but biased approximation for $1/Z(\theta)$. The Russian Roullette algorithm (Lyne et al.,, 2015) addresses this bias through a clever geometric series correction (details are in the supplement). This is an asymptotically exact algorithm and the theoretical assumptions are generally satisfied for general forms of the probability model $h(\mathbf{x}|\theta)$. However, it requires multiple $\widehat{Z}_{i}(\theta)$s at each iteration and each $\widehat{Z}_{i}(\theta)$ approximation itself requires multiple MCMC samples from $h(\mathbf{x}|\theta)$, making it computationally very expensive.

3.3 Noisy MCMC and Hybrids

If a Markov chain with transition kernel $P$ satisfies detailed balance with respect to the target $\pi$, it is asymptotically exact. However, when $P$ is approximated by $\widehat{P}$, the samples generated may only approximately follow the target $\widehat{\pi}$. Such algorithms may be generically referred to as “noisy MCMC”. Alquier et al., (2016) describes a broad class of noisy MCMC algorithms and uses total variational distance to quantify the distance between the asymptotically exact and inexact chain. Noisy MCMC in its broader sence is a large class of algorithms, including, for instance, pseudo-marginal MCMC and ALR. Here we discuss a noisy MCMC algorithm that builds upon the exchange algorithm; hence this algorithm may be thought of as a hybrid between auxiliary and likelihood approximation-based methods. In the exchange algorithm, a single auxiliary variable $\mathbf{y}$ is generated from $h(\mathbf{y}|\theta^{\prime})/Z(\theta^{\prime})$. Therefore the $h(\mathbf{y}|\theta)/h(\mathbf{y}|\theta^{\prime})$ term in the (4) may be thought of as a one-sample unbiased importance sampling estimate of $Z(\theta)/Z(\theta^{\prime})$. Instead of using a single $\mathbf{y}$, if multiple auxiliary variables $\{\mathbf{y}_{1},...,\mathbf{y}_{N}\}$ are generated from the $h(\mathbf{y}|\theta^{\prime})/Z(\theta^{\prime})$ the resulting importance sampling approximation will have smaller variance. This is called the noisy exchange algorithm with acceptance probability

[TABLE]

For $1<N<\infty$, the algorithm is asymptotically inexact because the detailed balance condition does not hold. Liang and Jin, (2013) also proposes a similar approach called Monte Carlo MH (MCMH). Alquier et al., (2016) reports that the noisy exchange algorithm shows better mixing than the exchange algorithm, and proposes a noisy Metropolis adjusted Langevin (MALA) exchange algorithm, where the gradient of the intractable distribution is approximated to construct proposals, that further improves upon mixing. Also, the estimate for the ergodic average has smaller bias than the exchange algorithm in empirical studies. Of course, there may be a tradeoff between improving mixing and increasing computational costs per iteration. We skip details but mention the general noisy MCMC approach as it may lead to other useful algorithms. For instance, DMH with multiple auxiliary variables, say “noisy DMH”, would be a simple extension.

Components to be tuned: The number of auxiliary variables, $N$, is an additional component to be tuned compared to the exchange algorithm or DMH. The choice of $N$ depends on some tradeoffs: as $N$ becomes large, the constructed chain can have better mixing and estimates from the chain can have lower variance at the expense of computing time. Parallel computing may be helpful for sampling $N$ $\mathbf{y}_{i}$s, and evaluating $h(\mathbf{y}_{i}|\theta)/h(\mathbf{y}_{i}|\theta^{\prime})$ independently.

Theoretical justification: Alquier et al., (2016) provides an upper bound for total variation norm distance between $P$, transition kernel for the exchange algorithm, and $\widehat{P}$, transition kernel for the noisy exchange algorithm, based on a result in Mitrophanov, (2005)[Corollary 3.1], which requires uniform ergodicity of $P$. Two assumptions are used in this derivation: (1) the prior $p(\theta)$ is bounded away from [math] and $\infty$, (2) $q(\theta^{\prime}|\theta)$ is bounded away from [math] and $\infty$; both assumptions are typically satisfied.

4 Simulated and Real Data Examples

We now study the algorithms in the context of three models that are of general interest: (1) the Ising model, (2) a social network model, and (3) an attraction-repulsion point process model. These models also present different computational challenges. The code for this is implemented in R (Ihaka and Gentleman,, 1996) and C++, using the Rcpp and RcppArmadillo packages (Eddelbuettel et al.,, 2011). We use the examples to compare the efficiency as well as practical implementation challenges of the algorithms. Simulation settings for the Russian roulette algorithm are provided in the supplementary material.

4.1 The Ising Model

Consider an Ising model (Lenz,, 1920, Ising,, 1925) on an $m$ by $n$ lattice with parameter $\theta$. The observed lattice $\mathbf{x}=\{x_{ij}\}$ has binary values $x_{i,j}=\{-1,1\}$, where $i,j$ denotes row and column location in the lattice. The model is

[TABLE]

where spatial dependence is imposed via $S(\mathbf{x})$. The larger $\theta$ becomes, the stronger the interactions between the data on the lattice. Summation over all $2^{mn}$ possible configurations of this model is required for the calculation of the normalizing function, which is computationally expensive even for lattices of moderate size. The simulations are conducted using perfect sampling (Propp and Wilson,, 1996) on a 10$\times$10 lattice with uniform prior [0,1] with parameter settings representing moderate dependence, with $\theta=0.2$, and strong dependence, with $\theta=0.43$.

All algorithms are tuned according to the methods described in the previous section. $h(\mathbf{y}|\widehat{\theta})/Z(\widehat{\theta})$ is chosen as the conditional density of the auxiliary variable y for the AVM algorithm, where $\widehat{\theta}$ is MPLE. The auxiliary variable is generated by 10 cycles of Gibbs updates in DMH. For AEX, 100 particles are selected among 5,000 samples from fractional DMH (descriptions are in the supplementary material). Then the preliminary run of the auxiliary chain is implemented for 420,000 iterations with the first 20,000 iterations discarded for burn-in. To expedite computing and reduce memory management issues, the resulting samples are thinned (at equally spaced intervals of 20) to obtain 20,000 samples. In the auxiliary chain, we set $n_{0}=20,000$, $K_{s}=[0,100^{100+10s}]^{100}$, and $\mathbf{x}_{n}$ is updated by a single cycle of Gibbs updates. For ALR, 100 particles are drawn from the uniform prior and $\mathbf{x}_{n}$ is updated by a single cycle of Gibbs updates. The kernel giving equal weights for $h=10$ nearest particles and 0 for others is used. In noisy DMH, 100 samples are used to produce importance sampling estimate and same inner sampler is constructed as DMH. We note that in order to make the computations feasible, we used parallel computing to obtain importance sampling estimates for both noisy DMH and the Russian Roulette algorithms. The parallel computing was implemented through OpenMP with the samples generated in parallel across 8 processors. We treated a run from the exchange algorithm as our gold standard; it was run for 1,010,000 iterations with 10,000 discarded for burn-in and 10,000 thinned samples are obtained from the remaining 1,000,000. Same simulation settings are used for $\theta=0.43$ case. Since perfect sampling takes very long for $\theta=0.43$ case, we use the ALR algorithm as the gold standard. All algorithms were run until the Monte Carlo standard errors calculated by batch means (Jones et al.,, 2006, Flegal et al.,, 2008) are below 0.01.

Table 1 shows that the estimates from different algorithms are well matched to the gold standard when $\theta=0.2$. Posterior densities in Figure 6 also indicate agreement. Likelihood approximation approaches show larger ESS than auxiliary variable approaches. On the other hand, auxiliary variable approaches have smaller computational costs than likelihood approximation approaches. In particular, DMH shows the shortest computing time. For this parameter setting, inference results from both asymptotically exact and inexact algorithms are accurate.

However when $\theta=0.43$, inference results are biased for asymptotically inexact algorithms. In Table 1 and Figure 6, it is observed that the inference results from DMH and Noisy DMH do not match the gold standard. Due to the strong dependence at this parameter setting, mixing of the inner sampler to generate samples from the $h(\mathbf{x}|\theta)/Z(\theta)$ is slower than for the $\theta=0.2$ case. Therefore a large number of Gibbs updates are necessary for accurate inference. Although Noisy DMH provides a closer estimate to the true value than DMH because it uses multiple samples for estimating $Z(\theta)/Z(\theta^{\prime})$, it is still biased. The number of iterations for the inner sampler in the Russian Roulette algorithm also has to be increased due to the slow mixing (see supplementary material for details). We can also observe that AVM and the exchange algorithm take several hours, because perfect sampling takes longer to achieve coalescence. On the other hand, there are no extra costs for AEX and ALR in the presence of strong dependence. For both algorithms, still a single cycle of Gibbs update is enough for updating $\mathbf{x}_{n}$. This is because both algorithms require updating not sampling $\mathbf{x}_{n}$ from the probability model.

Summary: In the case where the Ising model has only moderate dependence, all the algorithms work reasonably well. Here, DMH may be preferable because it is the easiest to construct and the fastest to run. On the other hand, in the context of Ising models that exhibit strong dependence, the computational problem becomes more challenging. For the asymptotically inexact algorithms like DMH, in order to obtain accurate results in the presence of strong dependence in the Ising model, the number of iterations for the inner sampler needs to be very large. This results in a much larger computational burden. AVM and the exchange algorithms also become impractical because perfect sampling takes too long. AEX and ALR can be useful for spatial autologistic models because both algorithms can attain accurate estimates with moderate computing speed even in the presence of strong dependence.

4.2 Social Network Models

Exponential random graph models (ERGM) (Robins et al.,, 2007, Hunter et al.,, 2008) provide an approach for modeling relationships among nodes of a network. Consider the undirected ERGM with $n$ nodes. For all $i\neq j$, $x_{i,j}=1$ if the $i$th node and $j$th node are connected, otherwise $x_{i,j}=0$ and $x_{i,i}$ is defined as 0. This forms an $n$ by $n$ adjacency matrix, $\mathbf{x}$. Calculation of the normalizing function requires summation over all $2^{n(n-1)/2}$ configuration, which is infeasible. Here we investigate two ERGM examples, the first with 4 parameters and the second with 9 parameters; the latter example also involves a more complicated summary statistic which poses an additional computational challenge. We have not implemented AVM and the exchange algorithms because perfect samplers are possible only for a few special cases (Butts,, 2012).

4.2.1 A Florentine Business Network

In this example, we study the Florentine business dataset (Padgett,, 1994), which describes the business networks among 16 Florentine families. Consider the undirected ERGM, where the probability model is

[TABLE]

[TABLE]

Sufficient statistics $S(\mathbf{x})=\left\{S_{1}(\mathbf{x}),S_{2}(\mathbf{x}),S_{3}(\mathbf{x}),S_{4}(\mathbf{x})\right\}$ represent the number of edges, two-stars, three-stars and triangles respectively, where $x_{i+}$ indicates row sum of adjacency matrix. Triangle represents the number of cyclic relationships and $k$-star indicates the number of nodes which have exactly $k$ relationships.

Each of the algorithms is tuned according to the guidelines in the previous sections. The auxiliary variable is generated by 10 cycles of Gibbs updates in DMH. In AEX, 200 particles are selected among 5,000 samples from fractional DMH, as described in the supplementary material. Then the preliminary run of the auxiliary chain is implemented for 630,000 iterations with the first 30,000 iterations discarded for burn-in. The resulting samples are thinned (at equally spaced intervals of 20) to obtain 30,000 samples. In the auxiliary chain, we set $n_{0}=20,000$, $K_{s}=[0,100^{100+10s}]^{200}$, and $\mathbf{x}_{n}$ is updated by a single cycle of Gibbs updates. For ALR, 400 particles are chosen from the short run of DMH, and a single cycle of Gibbs updates are used to update $\mathbf{x}_{n}$. The same kernel with $h=20$ is used as previous example. 100 samples are used for the importance sampling estimate in noisy DMH. Importance sampling estimates in Russian Roulette and noisy DMH algorithms are obtained in parallel as in the previous example. To obtain a gold standard for comparisons, we run the AEX algorithm 10 times independently. Each run consists of 101,000 iterations with 1,000 iterations discarded as burn-in. Then 10,000 samples are obtained by thinning from the remaining 10 sets of 100,000 samples. All algorithms were run until the Monte Carlo standard error is at or below 0.02.

Here we only provide the inference results regarding $\theta_{2}$ because similar results are observed for the other parameters. Table 2 indicates that the estimates from the different algorithms are similar to the those of the gold standard. This is because 10 cycles of Gibbs updates are enough to generate auxiliary samples from the likelihood in this example. However, for fewer iterations, say 1 to 2 cycles, the asymptotically inexact algorithms are biased. As in the Ising model example, less correlated samples can be generated from likelihood approximation approaches, at additional computational expense. Both ALR and AEX are computationally efficient compare to the Russian Roulette algorithm. This is because we can effectively store the previous samples using sufficient statistics in the likelihood. However, ALR is relatively expensive here (unlike in the Ising model), because it takes longer to visit the entire state space equally with increasing dimensions. The performance of AEX is relatively robust in the multidimensional case and is the fastest among asymptotically exact algorithms if the particles are carefully chosen.

Summary: AEX is the most reliable approach because it is asymptotically exact while at the same time retaining some computational efficiency because this model has low-dimensional statistics. On the other hand, DMH is simple and computationally efficient and as long as the length of the inner sampler is reasonable, it also provides accurate inference.

4.2.2 An Emergent Multi-organizational

Network (EMON)

In this example, we study the Mount St. Helens emergent multi-organizational network (EMON) dataset (Drabek et al.,, 1981), which describes communication networks for the search and rescue operations among 27 organizations. There are three node attributes in the model: (1) sponsorship level (city, county, federal, private, and state), (2) command rank score (higher score indicates higher rank), (3) whether headquarters are sited locally to the impact area or not (L/NL). Consider the undirected ERGM, where the probability model is

[TABLE]

[TABLE]

The sufficient statistics are $S_{1}(\mathbf{x})$ (the number of edges), $S_{2}(\mathbf{x})$ (triangles), $S_{3}(\mathbf{x})-S_{6}(\mathbf{x})$ (node factor for sponsorship level), $S_{7}(\mathbf{x})$ (node covariance for command rank score), and $S_{8}(\mathbf{x})-S_{9}(\mathbf{x})$ (graphlet orbit factor for location of headquarters). Graphlet statistics are small, connected subgraphs which represent the certain topological structure of a network (Pržulj et al.,, 2004, Pržulj,, 2007). Here we used Graphlet 4 with automorphism orbit 7 for describing brokerage roles between the 27 organizations (Yaveroglu et al.,, 2014).

In this model there are practical implementation issues for both particle-based algorithms (ALR and AEX): (1) Both algorithms require increasing number of particles to cover the high-dimensional parameter space. (2) The algorithms used for generating particles can exhibit slow mixing. For example the default implementation of DMH mixes slowly, which implies that it takes longer chains to generate the necessary particles for the ALR or AEX algorithms. (3) Even after the particles are generated, the stochastic approximation algorithm can still be slow when the number of parameters is large. Even with large gain factor (AEX) or step size (ALR), visiting each state space equally is infeasible when the parameter dimension is high (as in the 9-dimensional case above). This issue is perhaps even more problematic than issues (1) and (2) above.

The Russian Roulette algorithm is also computationally expensive because of the complexity of summary statistics in this example (wall time of roughly 1 month). Therefore, we study only the DMH and noisy DMH algorithms for this example. The details of our implementation are as follows. We generate auxiliary variables via 10 cycles of Gibbs updates, and 300 samples are used to construct an importance sampling estimate in noisy DMH. Importance sampling approximations in noisy DMH are evaluated through parallel methods, as in previous examples. We use the DMH algorithm with 20 cycles as the gold standard; it was run for 101,000 iterations with 1,000 samples discarded for burn-in and 10,000 thinned samples are obtained from the remaining samples. All algorithms are run until the Monte Carlo standard error is no larger than 0.03.

Here we provide inference results for $\theta_{9}$ because similar results are observed for the other parameters. Table 3 shows that posterior means from the different algorithms are well matched to the gold standard. It is observed that the acceptance rate is much lower than in the simpler Florentine business data example. The highest posterior density (HPD) interval obtained from Noisy DMH is slightly wider than that of the gold standard. As in the previous examples, a Noisy DMH algorithm can generate less correlated samples with higher acceptance rates compared to the DMH algorithm.

Summary: Particle-based algorithms (AEX and ALR) are infeasible for the 9-dimensional parameter case even though there are summary statistics to be stored. Because of computationally expensive graphlet statistics, Russian Roulette is infeasible. With the choice of an appropriate length for the inner sampler, DMH can provide accurate inference, with the highest ESS/T. This fact suggests that for these challenging cases, DMH may still be a practical approach.

4.3 Spatial Interaction Point Process

A spatial point process $\mathbf{x}=\{x_{1},...,x_{n}\}$ is a realisation of random points in a bounded plane $S\subset R^{2}$. By introducing an interaction function $\phi(D_{ij})$ which is the function of distance between the coordinates of $x_{i}$ and $x_{j}$, a probability model may be used to describe spatial patterns among the point. Extending the Strauss process (Strauss,, 1975) which explains repulsion patterns among point, Goldstein et al., (2015) develop a point process model to explain both attraction and repulsion patterns of the cells infected with human respiratory syncytial virus (RSV). The interaction function is

[TABLE]

and the probability model is

[TABLE]

There are four parameters $\{\lambda,\theta_{1},\theta_{2},\theta_{3}\}$ in the model: $\lambda$ controls intensity of the point process and $\{\theta_{1},\theta_{2},\theta_{3}\}$ are the parameters controlling the interaction function. $\theta_{1}$ is the peak value of $\phi$, $\theta_{2}$ is value of $D$ at the peak of $\phi$ and $\theta_{3}$ represents descent rate after the peak. When the points are too close to each other, $\phi$ value in (15) is less than $1$ which means points have a tendency to remain apart. However as the distance between points is larger, $\phi$ value become increased which means that points clump together. Attraction patterns become smaller as the distance between the points become larger. Likewise, this model can capture attraction repulsion spatial association among infected cells. Calculation of the normalizing function requires integration over the continuous domain $S$, which is intractable.

Goldstein et al., (2015) implement DMH for three independent replicates of 3,200 points in a well with radius 1,350 pixels, which is about 10,000 points. Even with code written in C, inference took roughly 12 hours. This is because the number of points $n$ is large, and for each iteration of the DMH, thousands of birth-death MCMC steps are required to generate the auxiliary variable. To allow for a comparison with other algorithms, which are computationally more expensive than DMH, we work with a smaller point pattern; however, this pattern is still computationally challenging enough to serve as a good testbed for the various algorithms we consider. Simulations are conducted on a well with 337.5 radius pixels and number of points $n\approx 200$ without replicates. Point process $\mathbf{x}$ is generated through a long run of birth-death MCMC. We follow RSV-B simulation settings in Goldstein et al., (2015), where the true parameter is $\{\lambda\times 10^{4},\theta_{1},\theta_{2},\theta_{3}\}=\{4,1.2,15,0.3\}$. We use uniform priors on $\theta$. Since the number of data points is small, descent rate parameter $\theta_{3}$ is not recovered well. Therefore $\theta_{3}$ is fixed at $0.3$ and we infer the other three parameters.

Implementing AVM or the exchange algorithm is infeasible because perfect sampling is impossible for this example. Although both AEX and ALR can be implemented theoretically, it is not practical because there are no summary statistics. Without summary statistics, we need to store the distance matrix of cumulative point process samples with varying dimension around 200 by 200, which is a burden on memory. Therefore, we study DMH, noisy DMH and the Russian Roulette algorithms for this example. Each algorithm is tuned according to the previous sections. For all the algorithms, samples are generated from the likelihood through 2,000 iterations of birth-death MCMC. For noisy DMH, 100 samples are used to construct an importance sampling approximation with each iteration. Importance sampling approximations in noisy DMH and Russian Roulette are evaluated through parallel methods, as in the past. We use the Russian Roulette algorithm as the gold standard; it was run for 101,000 iterations with 1,000 samples discarded for burn-in and 10,000 thinned samples are obtained from the remaining samples. All algorithms were run until the Monte Carlo standard error is at or below 0.01.

Here we only provide results for $\theta_{1}$ because similar results are observed for the other parameters. Table 4 shows that inference results from the different algorithms are well matched to the gold standard and the highest posterior density (HPD) intervals cover true values. As in the previous examples, likelihood approximation approaches can generate less correlated samples with higher acceptance rates.

Summary: DMH has significant advantages in terms of computational efficiency. Compared to the Russian Roulette algorithm which takes about 55 hours, DMH only takes several minutes to be implemented. This is because both calculation of the $h(\mathbf{x}|\theta)$ and the inner sampler to generate auxiliary variables are expensive compared to previous examples. Furthermore, considering ESS/T, DMH can generate effective sample within the shortest time, while simulated samples are close to the gold standard. This fact demonstrates that for computationally expensive problem, especially the model without low-dimensional sufficient statistics, we recommend DMH as a practical approach, though asymptotically exact inference is not guaranteed theoretically. Considering that the example is 16 times smaller than the original, it is infeasible to implement other algorithms for the original problem.

5 Computational Complexity

Here we investigate the computational complexity of the algorithms, that is, how the algorithms scale as we increase the number of data points. We believe this is particularly relevant in modern statistics since a common and important question to ask is how well an algorithm works as data sizes get large. Algorithms requiring perfect sampling will not be considered, because perfect sampling can only be constructed for limited cases. Let $L(\cdot)$ be the complexity of evaluating $h(\mathbf{x}|\theta)$, $G(\cdot)$ be the complexity of inner sampler, and $n$ be the number of points in data $\mathbf{x}$.

We begin with a few caveats. It is challenging to calculate complexity because some algorithms contain random quantities (stopping times), or have many components that need to be tuned specifically for different problems. Also, the mixing of the Markov chains and the quality of asymptotically inexact algorithms are influenced by the data itself (for instance, dependence) and not just the size of the data. We also note that to simplify calculations, we assume the dimensions of the data $\mathbf{x}$ and the auxiliary variable $\mathbf{y}$ are assumed to be same ($n$). This is not strictly accurate for the point process example where the auxiliary variable is generated through birth-death MCMC with varying dimensions. However, the dimensions are approximately similar to original data’s dimension $n$. In the point process example we can store data $\mathbf{x}$ itself or distance matrix of $\mathbf{x}$ in the adaptive algorithms (AEX, ALR) to evaluate $h(\mathbf{x}|\theta)$. The latter can avoid re-computations at the expense of memory. Here we assume that distance matrix of $\mathbf{x}$ will be stored. We assume that the length of the inner sampler is proportional to $n$. It is true when the inner sampler is used for update (AEX, ALR) which requires a single cycle of update (proportional to $n$). However, it may not be correct if the purpose of the inner sampler is sampling (DMH, Russian Roulette, Noisy DMH) from the probability model. The mixing of the inner sampler can be different depending on the parameter settings. However, several cycles of update (proportional to $n$) appears to work in practice.

All the algorithms require sampling (DMH, Noisy DMH, Russian Roultte) or updating (AEX, ALR) from the probability model and evaluating $h(\mathbf{x}|\theta)$. In exponential family models such as Ising model or ERGM, complexity of $L(\cdot)$ is of order $n$, because the number of calculations for $S(\mathbf{x})$ is proportional to $n$. For the inner sampler, when a single point is proposed to be changed, only neighboring points (fixed number regardless of $n$) are affected due to the Markovian assumption. Because the length of the inner sampler is assumed to increase with $n$, $G(\cdot)$ is also of order $n$. On the other hand, complexity of $L(\cdot)$ for point process model is of order $n^{2}$, because evaluation of the $h(\mathbf{x}|\theta)$ requires calculating an $n$ by $n$ distance matrix for $n$ data points, and evaluating interaction function $\phi(\cdot)$ on the corresponding distance matrix. For the inner sampler, when a single point is proposed to be added (birth) or deleted (death), the distance of a proposed point from other $n$ points should be calculated, and then $\phi(\cdot)$ should be evaluated at each point. Since the length of the inner sampler is assumed to be proportional to $n$, $G(\cdot)$ is of order $n^{2}$.

There is a big difference in calculating $h(\mathbf{x}|\theta)$ for exponential family and point process model. For the exponential family model, once we evaluate $S(\mathbf{x})$, we can simply take the product of $\theta$ and $S(\mathbf{x})$ for evaluation of $h(\mathbf{x}|\theta)$ in different $\theta$. On the other hand for point process model, $h(\mathbf{x}|\theta)$ should be recalculated; $\phi(\cdot)$ should be evaluated at the distance matrix of $\mathbf{x}$ with different parameters. Here we provide our main observations (see supplement for details): (1) Complexity of exponential family is $\mathcal{O}(n)$ and $\mathcal{O}(n^{2})$ in point process for all the algorithms. (2) Although algorithms have the same complexity, there are major differences in calculations per iteration. The amount of calculations per iteration of AEX, ALR and DMH are similar for exponential family models. Considering complexity and memory costs, DMH is the cheapest algorithm for both models. (3) Adaptive algorithms (AEX and ALR) require more memory. With increasing iterations, both algorithms become slower because algorithms use cumulative samples $H_{n+1}$ in calculations per iteration (AEX: Step 5 of the Algorithm 6, ALR: Step 4 of the Algorithm 8). Memory may not be an issue for the exponential family model. However for models without low-dimensional summary statistics, memory costs can become prohibitively expensive.

Figure 7 is the observed computing time for several algorithms with different scales in both models. We only include parts of results for Noisy DMH and the Russian Roulette algorithm. This is because both algorithms are too expensive to compare with other algorithms for large scales. For different scales, mixing of each algorithm is determined to be similar based on effective sample sizes, once we use appropriate step size (covariance) for proposal. Step size can be tuned to achieve similar acceptance rate for different scales; the larger data size becomes, the smaller step is used. We can also adaptively update step size (Atchadé,, 2006, Atchade et al.,, 2008). Figure 7 supports our calculations about $\mathcal{O}(n)$ complexity for exponential family and $\mathcal{O}(n^{2})$ complexity for point process model. Also in the exponential family model, slopes of AEX, ALR and DMH are similar to each other. However DMH is the fastest because of memory requirements of both adaptive algorithms. These facts are consistent with our calculations.

6 Connections and Summary of Results

Here, we point out connections between the algorithms. Also, based on our study, we provide some conclusions about advantages and disadvantages of each algorithm.

6.1 Connections and Observations

All the algorithms require sampling from the probability model either approximately (MCMC or particle methods) or exactly (perfect sampling). These samples, in turn, are used in some form of an importance sampling estimate. This is clear in the likelihood approximation approaches; for the auxiliary variable approaches one may think of the acceptance probability of the Metropolis-Hastings algorithm as containing a single sample importance sampling approximation of $Z(\theta)/Z(\theta^{\prime})$. Let the conditional density of the auxiliary variable in AVM be $f(\mathbf{y}|\theta,\mathbf{x})=h(\mathbf{y}|\widehat{\theta})/Z(\widehat{\theta})$, where $\widehat{\theta}$ is the MPLE. The part of the acceptance probability only related to the auxiliary variable is

[TABLE]

Murray et al., (2006) point out that since $\mathbf{y}\sim h(\cdot|\theta)/Z(\theta)$ and $\mathbf{y}^{\prime}\sim h(\cdot|\theta^{\prime})/Z(\theta^{\prime})$, $h(\mathbf{y}^{\prime}|\widehat{\theta})/h(\mathbf{y}^{\prime}|\theta^{\prime})$ and $h(\mathbf{y}|\widehat{\theta})/h(\mathbf{y}|\theta)$ may be thought of as one-sample unbiased importance sampling approximations of $Z(\widehat{\theta})/Z(\theta^{\prime})$ and $Z(\widehat{\theta})/Z(\theta)$ respectively. Therefore, (17) is the ratio of two unbiased estimates which is a biased estimate of $Z(\theta)/Z(\theta^{\prime})$. Murray et al., (2006) explain that compared to AVM, the exchange algorithm is more direct because $h(\mathbf{y}|\theta)/h(\mathbf{y}|\theta^{\prime})$ in the (4) is one-sample unbiased importance sampling estimate for $Z(\theta)/Z(\theta^{\prime})$ where $\mathbf{y}\sim h(\cdot|\theta^{\prime})/Z(\theta^{\prime})$. Though sampling schemes for the auxiliary variable are different, the same logic is applied to DMH and AEX. The only difference is that DMH generates $\mathbf{y}$ via MCMC, and AEX generates $\mathbf{y}$ via resampling (dynamic importance sampling from the mixture distribution). Both classes of algorithms are clearly connected through importance sampling. We summarize these connections in Figure 8.

AEX lies at the intersection of both classes of algorithms. As in the likelihood approximation approach, in the auxiliary chain $w^{(i)}_{n}$ approximates (up to a constant) $Z(\theta^{(i)})$ for each $\theta^{(i)}$. As in the auxiliary variable methods, the intractable functions are cancelled in the target chain at each iteration. Furthermore, AEX is closely connected to ALR. Both algorithms approximate $Z(\theta^{(i)})$ at each particle – $w_{n}^{(i)}$ in AEX and $c_{n}^{(i)}$ in ALR. Using their respective approximations, both collect a sample from a family of distributions $\{h(\mathbf{x}|\theta^{(1)})/Z(\theta^{(1)}),...,h(\mathbf{x}|\theta^{(d)})/Z(\theta^{(d)})\}$ via (different) stochastic approximations and save a sample with each iteration. As the number of iterations grows, the accumulated dataset $H_{n+1}$ grows, which makes both algorithm asymptotically exact. In AEX the resampling distribution of $\mathbf{y}$ becomes close to $h(\mathbf{y}|\theta^{\prime})/Z(\theta^{\prime})$ and in ALR the approximation $\widehat{Z}_{n+1}(\theta)$ converges to the truth. The difference is that AEX cancels out $Z(\theta)$ while ALR plugs-in $\widehat{Z}_{n+1}(\theta)$ into the acceptance probability.

In general, likelihood approximation approaches have better mixing than auxiliary variable approaches. However auxiliary variable approaches are less expensive per iteration resulting in higher effective sample size per second. Hence, auxiliary variable approaches tend to often be faster than likelihood approximation approaches. However the efficiency of algorithms can change depending on the model, parameter settings, and the data, as well as decisions about tuning components in each algorithm. In particular, strong dependence in observations tends to, slow sample generation, both with approximately and exact samplers. Perfect sampling takes longer to achieve coalescence, which can slow the AVM and exchange algorithm. Also, MCMC sampling requires longer chains, increasing computing time for the DMH and the Russian Roulette algorithms as well. In this case, AEX and ALR show relatively robust performance because both algorithms require updating $\mathbf{x}_{n}$, not sampling $\mathbf{x}_{n}$ from the probability model. For both algorithms, the choice of particles is crucial in determining both statistical and computational efficiency. If the parameter has higher dimensions, fewer particles are likely to be in the higher density regions of the posterior distribution of $\theta$. Then the approximation of $Z(\theta)$ might be inaccurate (ALR) or resampled auxiliary variable might be improbable (AEX), which can lead to poor mixing. Requiring an increase in the number of particles also slows computing. Therefore, careful choice of particles is important (see the supplementary material), especially for multidimensional parameters.

6.2 Guidelines and Recommendations

Comparing asymptotically exact and inexact algorithms is challenging in general because we have to simultaneously account for two different aspects of the algorithms: (1) mixing of the Markov chain, which affects the rate at which sample-based approximations converge to the true values, and (2) the quality of the target approximation in the case of asymptotically inexact algorithms. In order to provide reasonable comparisons, we always check the final results in each case using a variety of diagnostics, for example plots of marginal distributions as the Monte Carlo sample size increases, in order to convince ourselves that the approximation we obtain finally is close enough to the truth so that the above issues are avoided. Once this is the case, we can compare the efficiency of the algorithms via effective sample size (ESS) and effective samples per unit time (ESS/T). We have used ESS as a practical criteria of comparing mixing of the different algorithms; ESS/T can measure the computational efficiency as well.

Table 5 summarizes algorithms based on the following criteria:

• •

Class: which class the each algorithm falls in. A,L indicates auxiliary variable approaches and likelihood approximation approaches respectively and A/L indicates both.

• •

Exact or inexact: whether the stationary distribution of the Markov chain is exactly equal to the target posterior in the long run or not.

• •

Adaptation: whether the algorithm is adaptive (stationary distribution changes) or not.

• •

Requirement: key requirement(s) for the algorithm to be practical.

• •

Ease of use (rank): how much users need to make decision to run the algorithm. The smallest number indicates the simplest one.

AVM and the exchange algorithm propose the novel idea of cancelling out $Z(\theta)$ and thus have inspired many other algorithms. However, it is difficult to implement both algorithms in practice because the algorithms require perfect sampling. Even if we can implement perfect sampling it still has limited applicability because perfect sampling is very slow, either in the presence of strong dependence or for large scale data.

AEX, ALR, and Russian Roulette are asymptotically exact without perfect sampling. AEX is the fastest, especially when there are low-dimensional sufficient statistics for the model. AEX is faster than ALR for multidimensional $\theta$ cases, because AEX doesn’t have random stopping time ($\tau$). For ALR it is observed that $\tau$ in Algorithm 7 can be large for multidimensional $\theta$ problems, because it becomes difficult to move around the state space equally. This makes ALR slower than AEX for the complex parameter space. Therefore we recommend AEX for the model with low-dimensional summary statistics such as spatial autologistic model and ERGM. However as we describe in Section 4.2.2, both algorithms become less practical as the parameter space increases.

Without low-dimensional summary statistics, AEX poses memory burdens. In principle, the Russian Roulette algorithm can be applicable for the general form of the likelihood without sufficient statistics. The Russian Roulette algorithm does not require summary statistics or the assumption that $h(\mathbf{x}|\theta)$ should be bounded. However, it may not be practical because constructing unbiased likelihood estimate requires $\tau$ number of unbiased estimates $\widehat{Z}_{i}(\theta)$ per iteration. If $\widehat{Z}_{i}(\theta)$ requires $N$ samples, $\tau N$ samples need to be generated from the likelihood using inner sampler per iteration. For this reason, the Russian Roulette algorithm is computationally expensive compared to other methods.

We recommend DMH for computationally expensive likelihoods without low-dimensional summary statistics. Since only a single sample from the likelihood is required per iteration, DMH is computationally cheaper than any other algorithms. Furthermore, the algorithm is widely applicable without requiring any assumptions and relatively robust for high-dimensional parameter space. Once, we choose the appropriate number of the inner sampler updates, we can obtain plausible inference results in practice as well as the effective samples within the shortest time (ESS/T). However exactness is not guaranteed in DMH. Finally, if we can afford some additional computing expenses, noisy DMH may be useful because it can result in fast mixing chains.

7 Discussion

We have discussed MCMC algorithms for likelihoods with intractable normalizing functions, explaining some theoretical underpinnings as well as examining practical implementation issues and efficiencies. We find that auxiliary variable approaches tend to be more efficient than likelihood approximation approaches, though efficiencies vary quite a bit depending on the dimension of the problem, and on the availability of low-dimensional sufficient statistics for models (especially for AEX and ALR algorithms). The specifics of the dataset also matter, for instance data sets with strong dependence can slow down algorithms for inferring parameters of spatial point processes. There is some overlap among the central ideas for most algorithms and all have a clear connection to importance sampling.

For users, practical implementation challenges are as important as theoretical justifications. chopin2015leave also point out ease of use (fewer components to be tuned) and generalizability of code as important criteria for comparing Bayesian computation algorithms. From this perspective, AEX and ALR may be difficult choices for practitioners though both algorithms are asymptotically exact and moderately fast. This is because users have to manually tune lots of components, and do so differently for different problems. AEX seems to be a good choice for social network examples, at least in exponential random graph models where models have low-dimensional sufficient statistics. However, even here, Jin et al., (2013) discuss the need to carefully tune AEX in order to address some well known degeneracy problems (Handcock,, 2003, Schweinberger,, 2011).

Overall, we suggest starting with the DMH algorithm for its ease of implementation and the fact that it is typically computationally efficient. If instabilities in results are observed, for instance when the inner sampler length is increased, then AEX may be a more robust approach, even though it is much harder to implement. A preliminary run of DMH can be helpful for both tuning (e.g. selecting particles) and debugging of the AEX algorithm. Sophisticated users of MCMC and particle-based algorithms could also consider ALR as an alternative, and noisy MCMC methods may also be useful if there are severe mixing issues.

Finally, we hope that we have convinced readers that there are a number of very ingenious approaches and ideas to consider for tackling inference with intractable normalizing function problems. There are opportunities for extending or even combining some of these algorithms – it is not hard to believe that with some adjustments to the algorithms and clever uses of, say, parallel computing, any of these algorithms, for example the Russian roulette algorithm, is likely to become more efficient and hence more practical. Many of the ideas here relate to the even broader modern research problem of inference with models where likelihood functions are entirely unavailable or intractable. In the intractable normalizing function case it is apparent that there is still a dearth of generic, automated approaches for constructing efficient MCMC algorithms. For many models inference is infeasible even for moderately high-dimensional problems (cf. Goldstein et al.,, 2015). Furthermore, careful approaches for comparing the efficiency and diagnosing convergence of asymptotically inexact algorithms is still an open problem. We hope that this manuscript will encourage further research on all of these challenging topics.

**Supplementary Material for Bayesian Inference in the Presence of Intractable Normalizing Functions

**

Jaewoo Park and Murali Haran

The supplementary material provides details about how to select particles in the adaptive exchange (AEX) algorithm and a description of the Russian roulette algorithm. It also provides details about how computational complexity calculations for various algorithm was calculated.

Appendix A Choice of particles in AEX

Liang et al., (2016) suggests an approach for choosing particles. Here, we provide the algorithm in detail , including information on how various algorithm components are to be tuned for choosing particles. In Step 1 of Algorithm 9, $N_{0}$ number of candidate particles are generated through fractional double Metropolis-Hastings (DMH). Fractional DMH is DMH with larger acceptance probability. By using an acceptance probability $\alpha^{\zeta}$ with $0<\zeta<1$ where

[TABLE]

which is the acceptance probability of the double Metropolis-Hastings algorithm, we can draw samples from a wider region than the actual parameter space; samples can cover the support of the target posterior. Then a max-min procedure (see Steps 2 and 3 of Algorithm 9) is conducted to choose $d$ particles among $N_{0}$ candidates. Through this procedure, the $d$ selected particles have almost the same coverage in parameter space as $N_{0}$ candidates. We may use other strategies for choosing particles. For a social network model example, Jin et al., (2013) selects particles through Approximate Bayesian Computation (ABC) (Beaumont et al.,, 2002, Marin et al.,, 2012), an approximate sampling scheme that does not involve evaluating the likelihood.

Components to be tuned: To select particles, we have to decide the number of candidate particles ($N_{0}$), and the number of particles ($d$) in Algorithm 9. As the dimension of $\theta$ increases, we recommend using a larger $d$ to move around parameter space smoothly in the auxiliary chain. In this manuscript, we use $d=100$ for a one-dimension case, and $d=200$ for the four-dimension case. In both cases, a few thousand, $N_{0}=5,000$, appear to be adequate.

Appendix B The Russian Roulette Algorithm

Let $\widetilde{Z}(\theta)$ be the upper bound of $Z(\theta)$ or any estimate of $Z(\theta)$. Then the likelihood is represented through multiplication with infinite geometric series as

[TABLE]

where $c(\theta)$ controls $|k(\theta)|<1$. Therefore with an infinite supply of independent unbiased estimates $\widehat{Z}_{i}(\theta)$ for $Z(\theta)$, an unbiased approximation of the likelihood can be constructed as follows.

[TABLE]

Another possible choice is using multiplication correction with a Maclaurin series expansion of the exponential function (Lyne et al.,, 2015). However both the geometric and Maclaurin series corrections require an infinite number of unbiased estimates $\widehat{Z}_{i}(\theta)$, which makes them impractical. To address this problem, Lyne et al., (2015) propose stochastic truncation of the infinite series called the Russian Roulette, originally used in physics (Carter and Cashwell,, 1975). Using simulated finite random stopping time $\tau$, only a $\tau$ number of $\widehat{Z}_{i}(\theta)$ are required to construct an unbiased likelihood approximation. This leads to an asymptotically exact algorithm.

Let the random stopping time be $\tau\geq 1$ with $p_{n}=P(\tau\geq n)$ for all $n\geq 1$, and $p_{1}=1$. The stopping time can be defined as $\tau=\inf\left\{m\geq 1:U_{j}\geq q_{j}\right\}$ where $U_{j}\sim$ i.i.d. $u(0,1),q_{j}\in(0,1]$, which makes $p_{n}=\prod_{j=1}^{n-1}q_{j}$. Then the partial sum $S_{\tau}=1+\sum_{n=1}^{\tau-1}\prod_{i=1}^{n}\widehat{k}_{i}(\theta)/p_{n}$ with $S_{1}=1$ is called the Russian Roulette estimate. Lyne et al., (2015) show that $S_{\tau}$ is the unbiased estimate of infinite series as

[TABLE]

This indicates that the infinite series can be represented as the expectation of a finite series via the Russian Roulette truncation. In this way, an unbiased approximation of the likelihood is achieved by

[TABLE]

which makes the expectation of $\widehat{L}_{\tau}(\theta|\mathbf{x})$ equal (20). $\widehat{L}_{\tau}(\theta|\mathbf{x})$ will be plugged into the acceptance probability of the Metropolis-Hastings algorithm as in Figure 9. This finite series truncation makes the algorithm practical, since the finite $\tau$ number of $\widehat{Z}_{i}(\theta)$ is used to construct unbiased approximation. The choice of $q_{j}$ directly has an impact on both variance of the likelihood approximation and the computing speed of the algorithm. Smaller $q_{j}$ can lead to smaller $\tau$ which is the required number of the $\widehat{Z}_{i}(\theta)$ for constructing $\widehat{L}_{\tau}(\theta|\mathbf{x})$. However Lyne et al., (2015) provides a proof that if $q_{j}$ is close to zero, variance of likelihood approximation become large or even infinite. Lyne et al., (2015) suggests that $q_{j}=\prod_{i=1}^{j}\widehat{k}_{i}(\theta)$ appears to work well in practice.

Algorithm 10 summarizes the Russian Roulette algorithm. Although the Russian Roulette algorithm is an asymptotically exact MCMC and applicable for general forms of $h(\mathbf{x}|\theta)$, the algorithm is computationally expensive. It requires multiple ($\tau$) $\widehat{Z}_{i}(\theta)$s with each iteration (Step 3(b)). Each $\widehat{Z}_{i}(\theta)$ approximation in turn requires multiple MCMC samples from the likelihood, making it computationally expensive.

Components to be tuned: When a geometric correction is used for the Russian Roulette algorithm, five components of the algorithm need to be tuned: raw estimate $\widetilde{Z}(\theta)$, $c(\theta)$ which controls $|k(\theta)|<1$, ($N$) the number of samples for constructing the importance sampling estimate at each iteration, ($m$) length of inner sampler to generate MCMC samples, and the importance function $h(\mathbf{x}|\theta^{(0)})/Z(\theta^{(0)})$. If a tight upper bound of $Z(\theta)$ can be derived, it is a computationally efficient choice for $\widetilde{Z}(\theta)$, because it can lead to small $\widehat{k}_{i}(\theta)$ and small $q_{j}$, which implies small $\tau$. However, in most situation it is difficult to get a tight upper bound of $Z(\theta)$. Therefore, importance sampling approximations can be used for $\widetilde{Z}(\theta)$. $c(\theta)$ and $N$ should be determined through trial and error to make $|\widehat{k}_{i}(\theta)|<1$. In our simulation examples, we set $N$ between 1,000 and 2,000. $c(\theta)$=0.4 is used for every example. The length of the inner sampler $m$ should be larger as there is strong dependence among the data to generate more reliable samples. A simple choice for importance function is using maximum pseudolikelihood estimates as $\theta^{(0)}$. Other types of importance sampling estimate (e.g. annealed importance sampling) can also be used to construct $\widehat{Z}_{i}(\theta)$. To reduce computing time, $\widehat{Z}_{i}(\theta)$ may be obtained in parallel.

Theoretical justification: The Russian Roulette algorithm is based on the pseudo-marginal framework. If $h(\mathbf{x}|\theta)$ is bounded, we can attain a positive unbiased approximation of the likelihood and implement the algorithm without any theoretical difficulties. Even when $h(\mathbf{x}|\theta)$ is unbounded, $|\widehat{L}_{\tau}(\theta|\mathbf{x})|$ can be substituted instead of $\widehat{L}_{\tau}(\theta|\mathbf{x})$ in the Algorithm 10 and weighted average of $f(\theta_{i})$ with sign of $\widehat{L}_{\tau}(\theta|\mathbf{x})$ is used to estimate $E[f(\theta)]$ without losing exactness. See “sign problem” in Lyne et al., (2015) for details. Therefore, this algorithm is generally applicable.

B.1 Details of the simulation settings used in the examples

Ising model: Simulation settings in the Russian Roulette algorithm follow Lyne et al., (2015). When $\theta=0.2$, an annealed importance sampling (AIS) estimate of $Z(\theta)$ is constructed using 1,000 samples, and each sample is generated with 2,000 temperatures from an uniform lattice using the schedule. Geometric series correction with $c(\theta)=0.4$ is used to construct an unbiased estimate of likelihood. Same simulation settings are used for $\theta=0.43$ case except that 4,000 temperatures are used. Temperature schedules should be increased due to the slow mixing, which makes the algorithm more expensive.

Social network model: 2,000 samples are generated through 20 cycles of Gibbs updates and used for importance sampling estimate whose importance function is $h(\mathbf{x}|\theta^{(0)})/Z(\theta^{(0)})$ where $\theta^{(0)}$ is the MPLE. A geometric series correction with $c(\theta)=0.4$ is used to construct an unbiased estimate of the likelihood.

Spatial point process: 1,000 samples are generated through 2,000 iterations of birth-death MCMC and used for importance sampling approximations $\widehat{Z}_{i}(\theta)$ and $\widetilde{Z}(\theta)$. As an importance function, $h(\mathbf{x}|\theta^{(0)})/Z(\theta^{(0)})$ is used where $\theta^{(0)}$ is sample mean from short run of DMH. A geometric series correction with $c(\theta)=0.4$ is used for the Russian Roulette algorithm.

Appendix C Complexity calculations

Here we provide details about calculating complexity and memory costs for the algorithms. Table 7 summarizes these calculations. In Table 7, complexity indicates computational cost for a single run of the algorithms. Memory represents the magnitude of the quantities to be saved per each iteration to avoid re-computations. Fixed complexity measures all computational costs and fixed memory represents the size of variables to be saved before implementing the entire algorithm. We summarize our notation in Table 6. All other notations follow previous sections.

C.1 Complexity of DMH

In the double Metropolis-Hastings algorithm

Exponential family

• •

Complexity

Step 2: A single auxiliary variable $\mathbf{y}$ is generated.

Step 3: $S(\mathbf{y})$ is calculated once.

Total: $G(n)+L(n)$

• •

Fixed memory

Data: $p$-dimensional $S(\mathbf{x})$ should be stored.

Total: $p$ 2. 2.

Point process

• •

Complexity

Step 2: A single auxiliary variable $\mathbf{y}$ is generated.

Step 3: $h(\mathbf{x}|\theta^{\prime}),h(\mathbf{y}|\theta),h(\mathbf{y}|\theta^{\prime})$ are calculated.

Total: $G(n^{2})+3L(n^{2})$

• •

Fixed memory

Data: $n$ by $n$ distance matrix of $\mathbf{x}$ should be stored.

Total: $n^{2}$

C.2 Complexity of Noisy DMH

In the noisy double Metropolis-Hastings algorithm

Exponential family

• •

Complexity

Step 2: $N$ number of auxiliary variables $\mathbf{y}_{j}$ are generated.

Step 3: $N$ number of corresponding $S(\mathbf{y}_{j})$ are calculated.

Total: Steps can be parallelized, $N[G(n)+L(n)]/c$

• •

Fixed memory

Data: $p$-dimensional $S(\mathbf{x})$ should be stored.

Total: $p$ 2. 2.

Point process

• •

Complexity

Step 2: $N$ number of auxiliary variables $\mathbf{y}_{j}$ are generated.

Step 3: $N$ number of $h(\mathbf{y}_{j}|\theta),h(\mathbf{y}_{j}|\theta^{\prime})$ and a single $h(\mathbf{x}|\theta^{\prime})$ are calculated.

Total: Generating $N$ number of $\mathbf{y}_{j}$ and evaluating $N$ number of $h(\mathbf{y}_{j}|\theta),h(\mathbf{y}_{j}|\theta^{\prime})$ can be parallelized, $N[G(n^{2})+2L(n^{2})]/c+L(n^{2})$

• •

Fixed memory

Data: $n$ by $n$ distance matrix of $\mathbf{x}$ should be stored.

Total: $n^{2}$.

C.3 Complexity of AEX

In the adaptive exchange algorithm

Exponential family

• •

Complexity

Step 1(b): With $0.5$ probability, one MH update of $\mathbf{x}^{\prime}$ and $S(\mathbf{x}^{\prime})$ is calculated once.

Total: $0.5[G(n)+L(n)]$

• •

Memory

$p$-dimensional $S(\mathbf{x}_{m+1})$ and scalar $w_{m+1}^{(I_{m+1})}$, $I_{m+1}$ are stored.

Total: $p+2$

• •

Fixed complexity

Choosing particles: $N_{0}$ iterations of fractional DMH, $N_{0}[G(n)+L(n)]$

Preliminary run for auxiliary chain: $0.5N_{1}[G(n)+L(n)]$

Total: $[N_{0}+0.5N_{1}][G(n)+L(n)]$

• •

Fixed memory

Data: $p$-dimensional $S(\mathbf{x})$ should be stored.

Particles: $d$ number of $p$-dimensional particles, and $d$ by $d$ distance matrix of particles to define neighbors should be stored.

Preliminary run for auxiliary chain: $p$-dimensional $S(\mathbf{x}_{m+1})$ and scalar $w_{m+1}^{(I_{m+1})}$, $I_{m+1}$ are stored for $N_{1}$ iterations, and can be stored with $b$ thinning interval.

Total: $p+dp+d^{2}+(p+2)N_{1}/b$ 2. 2.

Point process

• •

Complexity

Step 1(a): With $0.5$ probability, $h(\mathbf{x}_{m}|\theta^{(I^{\prime})})$ is calculated once.

Step 1(b): With $0.5$ probability, one MH update of $\mathbf{x}^{\prime}$ and $h(\mathbf{x}^{\prime}|\theta^{(I_{m})})$ is calculated once.

Step 5: $|H_{m+1}|$ number of $h(\mathbf{x}_{j}|\theta^{\prime})$ should be calculated. With $N_{1}$ number of preliminary run and $b$ thinning interval for the auxiliary chain, $|H_{m+1}|$ is equal to $N_{1}/b+m+1$.

Step 6: $h(\mathbf{x}|\theta^{\prime}),h(\mathbf{y}|\theta),h(\mathbf{y}|\theta^{\prime})$ are calculated.

Total: $0.5G(n^{2})+[N_{1}/b+m+5]L(n^{2})$

• •

Memory

$n$ by $n$ distance matrix of $\mathbf{x}_{m+1}$, scalar $w_{m+1}^{(I_{m+1})}$, $I_{m+1}$ and evaluated $h(\mathbf{x}_{m+1}|\theta^{(I_{m+1})})$ are stored.

Total: $n^{2}+3$

• •

Fixed complexity

Choosing particles: $N_{0}$ iterations of fractional DMH, $N_{0}[G(n^{2})+3L(n^{2})]$

Preliminary run for auxiliary chain: $N_{1}[0.5L(n^{2})+0.5[G(n^{2})+L(n^{2})]]$

Total: $[N_{0}+0.5N_{1}]G(n^{2})+[3N_{0}+N_{1}]L(n^{2})$

• •

Fixed memory

Data: $n$ by $n$ distance matrix of $\mathbf{x}$ should be stored.

Particles: $d$ number of $p$-dimensional particles, and $d$ by $d$ distance matrix of particles to define neighbors should be stored.

Preliminary run for auxiliary chain: $n$ by $n$ distance matrix of $\mathbf{x}_{m+1}$, scalar $w_{m+1}^{(I_{m+1})}$, $I_{m+1}$ and evaluated $h(\mathbf{x}_{m+1}|\theta^{(I_{m+1})})$ are stored for $N_{1}$ iterations, and can be stored with $b$ thinning interval.

Total: $n^{2}+dp+d^{2}+(n^{2}+3)N_{1}/b$

In practice, AEX consists of three parts: (1) choosing particles, (2) preliminary run for the auxiliary chain, and (3) implementing the entire algorithm. Fixed complexity and fixed memory are evaluated from the first two parts.

C.4 Complexity of ALR

In the ALR algorithm

Exponential family

• •

Complexity

Step 1: $\mathbf{x}_{m+1}$ is generated once.

Step 2: $S(\mathbf{x}_{m+1})$ is calculated once.

Total: $G(n)+L(n)$

• •

Memory

$p$-dimensional $S(\mathbf{x}_{m+1})$ and scalar $I_{m+1}$ is stored.

Total: $p+1$

• •

Fixed complexity

Choosing particles: $N_{0}$ iterations of DMH, $N_{0}[G(n)+L(n)]$

Stopping time: Step 1-3 are repeated until $\tau$, $\tau[G(n)+L(n)]$

Total: $(N_{0}+\tau)[G(n)+L(n)]$

• •

Fixed memory

Data: $p$-dimensional $S(\mathbf{x})$ should be stored.

Particles: $d$ number of $p$-dimensional particles are stored.

Stopping time: $d$-dimensional $\mathbf{c}_{m}$ is stored.

Total: $p+dp+d$ 2. 2.

Point process

• •

Complexity

Step 1: $\mathbf{x}_{m+1}$ is generated once.

Step 2: $d$ number of $h(\mathbf{x}_{m+1}|\theta^{(i)})$ are calculated, which can be parallelized.

Step 4: $h(\mathbf{x}_{j}|\theta^{\prime}),h(\mathbf{x}_{j}|\theta)$ are evaluated for $m+1$ number of cumulative $\mathbf{x}_{j}$, which can be parallelized.

Step 5: $h(\mathbf{x}|\theta^{\prime})$ is evaluated once.

Total: $G(n^{2})+[(d+2m+2)/c+1]L(n^{2})$

• •

Memory

$n$ by $n$ distance matrix of $\mathbf{x}_{m+1}$, scalar $I_{m+1}$ and $d$ number of $h(\mathbf{x}_{m+1}|\theta^{(i)})$ should be stored.

Total: $n^{2}+d+1$

• •

Fixed complexity

Choosing particles: $N_{0}$ iterations of DMH, $N_{0}[G(n^{2})+3L(n^{2})]$

Stopping time: Step 1-3 are repeated until $\tau$, $\tau[G(n^{2})+dL(n^{2})/c]$

Total: $[N_{0}+\tau]G(n^{2})+[3N_{0}+\tau d/c]L(n^{2})$

• •

Fixed memory

Data: $n$ by $n$ distance matrix of $\mathbf{x}$ should be stored.

Particles: $d$ number of $p$-dimensional particles are stored.

Stopping time: $d$-dimensional $\mathbf{c}_{m}$ is stored.

Total: $n^{2}+dp+d$

Similar to AEX, the algorithm consists of three parts: (1) choosing particles, (2) updating $\{\mathbf{x}_{m},I_{m},\mathbf{c}_{m}\}$, until $\mathbf{c}_{m}$ converges, and (3) updating the entire process. Fixed complexity and fixed memory are evaluated from the first two parts.

It is observed that random stopping time $\tau$ is affected by the complexity of the parameter space and choice of particles rather than $n$. Intuitively, as the particles cover important regions of the parameter space, it is easier to move around the state space quickly. Therefore, proper choice of particles can lead to smaller $\tau$. We observe that if we choose particles covering interesting region of parameter space as in Algorithm 9, $\tau$ is almost fixed even with increasing $n$.

C.5 Complexity of the Russian Roulette algorithm

In the Russian roulette algorithm

Exponential family

• •

Complexity

Step 2: Importance sampling estimate, $\widetilde{Z}(\theta^{\prime})$ is constructed, which requires generating $N$ number of $\mathbf{y}_{j}$ and calculating $S(\mathbf{y}_{j})$.

Step 3: $\tau$ number of $\widehat{Z}_{i}(\theta^{\prime})$ are constructed, which require generating $N$ number of $\mathbf{y}_{j}$ and calculating $S(\mathbf{y}_{j})$ per each $\widehat{Z}_{i}(\theta^{\prime})$.

Total: Steps can be parallelized, $N(\tau+1)[G(n)+L(n)]/c$

• •

Fixed Memory

Data: $p$-dimensional $S(\mathbf{x})$ should be stored.

Total: $p$ 2. 2.

Point process

• •

Complexity

Step 2: Importance sampling estimate, $\widetilde{Z}(\theta^{\prime})$ is constructed, which requires generating $N$ number of $\mathbf{y}_{j}$ and calculating $h(\mathbf{y}_{j}|\theta^{\prime}),h(\mathbf{y}_{j}|\theta^{(0)})$.

Step 3: $\tau$ number of $\widehat{Z}_{i}(\theta^{\prime})$ are constructed, which require generating $N$ number of $\mathbf{y}_{j}$ and calculating $h(\mathbf{y}_{j}|\theta^{\prime}),h(\mathbf{y}_{j}|\theta^{(0)})$ per each $\widehat{Z}_{i}(\theta^{\prime})$. Then $h(\mathbf{x}|\theta^{\prime})$ is calculated.

Total: Step 2-3 can be parallelized, $N(\tau+1)[G(n^{2})+2L(n^{2})]/c+L(n^{2})$

• •

Fixed Memory

Data: $n$ by $n$ distance matrix of $\mathbf{x}$ should be stored.

Total: $n^{2}$.

[TABLE]

Because the stopping time $\tau$ is a random variable, we propose an approximate evaluation of the expected stopping time. Consider $q_{j}$, where $p(\tau\geq x)=\prod_{j=1}^{x-1}q_{j}$ as we defined in the previous section. We can use $q_{j}=\prod_{i=1}^{j}\widehat{k}_{i}(\theta)$, where $\widehat{k}_{i}(\theta)=1-c(\theta)\widehat{Z}_{i}(\theta)/\widetilde{Z}(\theta)$ as suggested in Lyne et al., (2015). For large enough importance sampling size $N$, $\widehat{Z}(\theta)_{i}/\widetilde{Z}(\theta)$ will be close to 1 when $\widetilde{Z}(\theta)$ is importance sampling estimate as well. In this case, $\widehat{k}_{i}(\theta)\approx k=0.6$ is almost constant because $c(\theta)=0.4$ is used for our examples. Therefore, $q_{j}\approx k^{j}$ and expected stopping time can be represented as (23). For $c(\theta)=0.4$, $\widetilde{\tau}^{R}\approx 1.86$.

Then remaining issue is how large $N$ should be depending on size of data $n$. $N$ should be tuned to make $|\widehat{k}_{i}(\theta)|<1$ for convergence of geometric series. As $N$ increases, it is oberved that the mixing of the chain become better (higher effective sample size) and ergodic average with lower Monte Carlo standard error can be attained at the expense of increased computational costs. How to tune $N$ efficiently is still debatable. In pseudo-marginal MCMC, Doucet et al., (2015) propose the rule of thumb to choose $N$ so that standard deviation of the log-likelihood estimate is around 1.2. Although this can be a one guideline to tune $N$, it is found that suggested $N$ is too large in our point process example which can lead unnecessary costs. Also, it is found that using Annealed importance sampling estimate can achieve lower variance in many problems than using simple importance sampling estimate. Considering all possibilities, it is challenging to generalize how $N$ should be tuned depending on $n$, because it can be problem specific.