Parallel Tempering via Simulated Tempering Without Normalizing Constants

TL;DR

This paper introduces a novel Bayesian sampling method that enhances Simulated Tempering by eliminating the need for normalizing constants, allowing for continuous temperature schedules and efficient joint parameter and marginal likelihood estimation.

Contribution

It develops a new Bayesian methodology that modifies Simulated Tempering to avoid calculating normalizing constants, enabling continuous temperature schedules and improved efficiency.

Findings

Effective in mixture models of Gaussian distributions

Applicable to ordinary differential equation models

Provides accurate marginal likelihood estimates

Abstract

In this paper we develop a new general Bayesian methodology that simultaneously estimates parameters of interest and the marginal likelihood of the model. The proposed methodology builds on Simulated Tempering, which is a powerful algorithm that enables sampling from multi-modal distributions. However, Simulated Tempering comes with the practical limitation of needing to specify a prior for the temperature along a chosen discretization schedule that will allow calculation of normalizing constants at each temperature. Our proposed model defines the prior for the temperature so as to remove the need for calculating normalizing constants at each temperature and thereby enables a continuous temperature schedule, while preserving the sampling efficiency of the Simulated Tempering algorithm. The resulting algorithm simultaneously estimates parameters while estimating marginal likelihoods through thermodynamic integration. We illustrate the applicability of the new algorithm to different examples involving mixture models of Gaussian distributions and ordinary differential equation models.

Figures

Tables

Table 1. Table 1: Parameter estimates v.s. theoretical values

	Parameter estimates	Theoretical results
$\mu$	1.405(0.0019); -1.413(0.0019)	1.409 ; -1.409
${\sigma }^2$	1.026(0.015)	1.148

Table 2. Table 2: Marginal log-likelihood of the bimodal model

	Analytic solution	PT-STWNC	TI-PT-B	TI-PT-NB
	-38.946	-37.767(0.019)
T =30			-38.254(0.006)	-38.261(0.006)
T =60	-	-	-38.247(0.004)	-38.254(0.005)
T =100	-	-	-38.247(0.004)	-38.254(0.003)

Table 3. Table 3: Log marginal likelihood estimates of Galaxy data

Model fitted	PT-STWNC	TI-PT-B	TI-PT-NB
1. 2 components equal variances	-241.99(0.25)	-238.02(0.03)	-238.03(0.02)
2. 3 components unequal variances	-228.67(0.39)	-224.26(0.03)	-224.28(0.03)
3. 3 components equal variances	-236.10(0.47)	-224.20(0.06)	-224.23(0.05)
4. 4 components unequal variances	-222.72(0.28)	-223.88(0.02)	-223.89(0.02)
5. 5 components unequal variances	-221.62(0.50)	-223.85(0.02)	-223.85(0.02)

Table 4. Table 4: log Bayes factors

	PT-STWNC	TI-PT-NB	TI-PT-B
$\log BF21$	13.32	13.74	13.75
$\log BF31$	5.89	13.79	13.82
$\log BF41$	19.27	14.14	14.13
$\log BF51$	20.37	14.17	14.16
$\log BF32$	-7.44	0.05	0.06
$\log BF42$	5.85	0.39	0.38
$\log BF52$	7.05	0.43	0.41
$\log BF43$	13.381	0.34	0.31
$\log BF53$	14.49	0.37	0.34
$\log BF54$	1.1	0.03	0.02

Equations

$$P(\mathbf{\theta }\mid \tau ,\mathbf{Y})\propto P(\mathbf{Y}\mid \mathbf{\theta },\tau )P(\mathbf{\theta })=P(\mathbf{Y}\mid \mathbf{\theta }{)}^{\tau }P(\mathbf{\theta }).$$

$$P_t(\mathbf{\theta }\mid \mathbf{Y})=\frac{P(\mathbf{Y}\mid \mathbf{\theta }{)}^{{\tau }_t}P(\mathbf{\theta })}{z(\mathbf{Y}\mid {\tau }_t)},\text{mbox}t\in \{1,..,\text{T}\},$$

$$z(\mathbf{Y}\mid {\tau }_t)=\underset{\mathbf{\Theta }}{\int }P(\mathbf{Y}\mid \mathbf{\theta }{)}^{{\tau }_t}P(\mathbf{\theta })d\mathbf{\theta }$$

$$\min \left\{1,\frac{P_t({\mathbf{\theta }}_l\mid \mathbf{Y})}{P_t({\mathbf{\theta }}_t\mid \mathbf{Y})}\frac{P_l({\mathbf{\theta }}_t\mid \mathbf{Y})}{P_l({\mathbf{\theta }}_l\mid \mathbf{Y})}\right\}$$

$$P(\mathbf{\theta },\tau \mid \mathbf{Y})=\frac{P(\mathbf{Y}\mid \mathbf{\theta }{)}^{\tau }P(\mathbf{\theta })P(\tau )}{P(\mathbf{Y})},$$

$$P(\mathbf{Y})=\underset{0}{\overset{1}{\int }}\underset{\mathbf{\Theta }}{\int }P(\mathbf{Y}\mid \mathbf{\theta }{)}^{\tau }P(\mathbf{\theta })P(\tau )d\mathbf{\theta }d\tau .$$

$$\begin{array}[]{ll}P(\tau\mid\boldsymbol{Y})=\int\limits_{\boldsymbol{\Theta}}P(\tau,\boldsymbol{\theta}\mid\boldsymbol{Y})d\boldsymbol{\theta}&\propto P(\tau)\int\limits_{\boldsymbol{\Theta}}P(\boldsymbol{Y}\mid\boldsymbol{\theta})^{\tau}P(\boldsymbol{\theta})d\boldsymbol{\theta}\\ \\ &\propto P(\tau)z(\boldsymbol{Y}\mid\tau)\\ \end{array}$$

$$P(\tau )\propto \frac{1}{z(\mathbf{Y}\mid \tau )},$$

$$\min \left\{1,\frac{P({\mathbf{\theta }}^{(i+1)},{\tau }^{\ast }\mid \mathbf{Y})}{P({\mathbf{\theta }}^{(i+1)},{\tau }^{(i)}\mid \mathbf{Y})}\right\}$$

$${\mathbf{\theta }}_{\max }(\tau )=\arg \underset{\mathbf{\theta }}{\max }P(\mathbf{\theta }\mid \tau ,\mathbf{Y}),$$

$$P(\tau ,{\mathbf{\theta }}_{\max }(\tau )\mid \mathbf{Y})=U(0,1).$$

$$P({\tau }_i,{\mathbf{\theta }}_{\max }({\tau }_i)\mid \mathbf{Y})=P({\tau }_j,{\mathbf{\theta }}_{\max }({\tau }_j)\mid \mathbf{Y}).$$

$$P(\tau ,{\mathbf{\theta }}_{\max }(\tau )\mid \mathbf{Y})$$

$$\frac{P(\mathbf{Y}\mid {\mathbf{\theta }}_{\max }(\tau ){)}^{\tau }P({\mathbf{\theta }}_{\max }(\tau ))P(\tau )}{P_{\text{prof}}(\mathbf{Y})}$$

$$P(\tau )=\frac{P(\tau ,{\mathbf{\theta }}_{\max }(\tau )\mid \mathbf{Y})P_{\text{prof}}(\mathbf{Y})}{P(\mathbf{Y}\mid {\mathbf{\theta }}_{\max }(\tau ){)}^{\tau }P({\mathbf{\theta }}_{\max }(\tau ))}.$$

$$\min \left\{1,\frac{P(\mathbf{Y}\mid \mathbf{\theta }{)}^{{\tau }^{\ast }}P({\tau }^{\ast })}{P(\mathbf{Y}\mid \mathbf{\theta }{)}^{{\tau }^{(i)}}P({\tau }^{(i)})}\right\}$$

$${\alpha }_{\mathbf{\theta }}=\frac{P({\mathbf{\theta }}^{\ast },{\tau }^{(i)}\mid \mathbf{Y})}{P({\mathbf{\theta }}^{(i)},{\tau }^{(i)}\mid \mathbf{Y})}=\frac{P({\mathbf{\theta }}^{\ast }\mid {\tau }^{(i)},\mathbf{Y})}{P({\mathbf{\theta }}^{(i)}\mid {\tau }^{(i)},\mathbf{Y})};$$

$${\alpha }_{\tau }$$

$${\alpha }_{\varphi }=\frac{P(\mathbf{Y}\mid {\mathbf{\theta }}_2^{(i)}{)}^{{\tau }_1^{(i)}}P(\mathbf{Y}\mid {\mathbf{\theta }}_1^{(i)}{)}^{{\tau }_2=1}}{P(\mathbf{Y}\mid {\mathbf{\theta }}_1^{(i)}{)}^{{\tau }_1^{(i)}}P(\mathbf{Y}\mid {\mathbf{\theta }}_2^{(i)}{)}^{{\tau }_2=1}}$$

$$\frac{P({\text{M}}_1\mid \mathbf{Y})}{P({\text{M}}_2\mid \mathbf{Y})}=\frac{P(\mathbf{Y}\mid {\text{M}}_1)}{P(\mathbf{Y}\mid {\text{M}}_2)}\frac{P({\text{M}}_1)}{P({\text{M}}_2)}.$$

$${\text{B}}_{12}=\frac{P(\mathbf{Y}\mid {\text{M}}_1)}{P(\mathbf{Y}\mid {\text{M}}_2)}$$

$$P(\mathbf{Y}\mid {\text{M}}_j)=\underset{{\mathbf{\Theta }}_j}{\int }P(\mathbf{Y}\mid {\mathbf{\theta }}_j,{\text{M}}_j)P({\mathbf{\theta }}_j\mid {\text{M}}_j)d{\mathbf{\theta }}_j,$$

$$\log (P(\mathbf{Y}\mid {\text{M}}_j))=\underset{0}{\overset{1}{\int }}{E}_{\mathbf{\theta }\mid \mathbf{Y},\tau ,{\text{M}}_j}[\log (P(\mathbf{Y}\mid \mathbf{\theta },{\text{M}}_j))]d\tau ,$$

$$\log (P(\mathbf{Y}\mid {\text{M}}_j))=$$

$$\approx$$

$$\log (P(\mathbf{Y}\mid {\text{M}}_j))\approx \underset{\text{Approximation}}{\underbrace{\frac{1}{2}\sum _{t=2}^{\text{T}}\Delta {\tau }_t(E_{t,{\text{M}}_j}+E_{t-1,{\text{M}}_j})}}+$$

$$\underset{\text{Bias}}{\underbrace{\frac{1}{2}\sum _{t=2}^{\text{T}}[\text{KL}(p_{t-1,{\text{M}}_j}\Vert p_{t,{\text{M}}_j})-\text{KL}(p_{t,{\text{M}}_j}\Vert p_{t-1,{\text{M}}_j})]}},$$

$$\log P(\mathbf{Y}\mid {\text{M}}_j)=\underset{0}{\overset{1}{\int }}{E}_{\mathbf{\theta }\mid \mathbf{Y},\tau ,{\text{M}}_j}[\log (P(\mathbf{Y}\mid \mathbf{\theta },{\text{M}}_j))]d\tau$$

$$\approx \sum _{t=2}^{\text{T}}\Delta {\tau }_t{E}_{t,{\text{M}}_j},$$

$$\mathrm{Var}(\log P(\mathbf{Y}\mid {\text{M}}_j))\approx \sum _{t=2}^{\text{T}}\mathrm{Var}(\Delta {\tau }_t{E}_{t,{\text{M}}_j}).$$

$$\mathrm{Var}(\log P(\mathbf{Y}\mid {\text{M}}_j))\approx \sum _{t=2}^{\text{T}}\frac{1}{{\text{N}}^2}\mathrm{Var}(E_{t,{\text{M}}_j}).$$

$$\frac{1}{\text{N}}\mathrm{Var}(\mathbf{\theta }\mid \mathbf{Y},\tau =0,{\text{M}}_j)$$

Full text

methodology that simultaneously estimates parameters of interest and the marginal likelihood of the model. The proposed methodology builds on Simulated Tempering, which is a powerful algorithm that enables sampling from multi-modal distributions. However, Simulated Tempering comes with the practical limitation of needing to specify a prior for the temperature along a chosen discretization schedule that will allow calculation of normalizing constants at each temperature. Our proposed model defines the prior for the temperature so as to remove the need for calculating normalizing constants at each temperature and thereby enables a continuous temperature schedule, while preserving the sampling efficiency of the Simulated Tempering algorithm.

The resulting algorithm simultaneously estimates parameters while estimating marginal likelihoods through thermo-

dynamic integration. We illustrate the applicability of the new algorithm to different examples involving mixture models of Gaussian distributions and ordinary differential equation models.

1 Introduction

Basic random walk Markov Chain Monte Carlo (MCMC) methods are inefficient when faced with multi-modal posterior distributions, especially when modes are isolated by large gaps of low probability. Simulated Tempering (Marinari and Parisi, 1992; Geyer and Thompson, 1995) and Parallel Tempering (Swendsen and Wang, 1986; Geyer, 1991; Hukushima and Nemoto, 1996), are MCMC variants designed to ease the challenges of multi-modal distributions by incorporating an auxiliary temperature parameter to overcome the prohibitively low probability regions which otherwise trap samplers in local modes (Zhang and Ma, 2008). Sampling occurs at different temperatures, balancing short distance within-mode steps and longer distance between-mode steps.

Both Parallel Tempering (PT) and Simulated Tempering (ST) define a sequence of distributions for the vector of data $\boldsymbol{Y}\in\mathbb{R}^{N}$ and parameter $\boldsymbol{\theta}\in\mathbb{R}^{d}$, indexed by an inverse temperature $\tau\in[0,1]$, along a path between the prior and the target distribution usually defined by

[TABLE]

At the extremes of the distribution sequence, $P(\boldsymbol{\theta}\mid\tau=0,\boldsymbol{Y})=P(\boldsymbol{\theta})$ is the prior, and $P(\boldsymbol{\theta}\mid\tau=1,\boldsymbol{Y})\propto P(\boldsymbol{Y}\mid\boldsymbol{\theta})P(\boldsymbol{\theta})$ is the usual target distribution. Consequently, when $\tau=0$ it is easy for the sampler to explore the parameter space, but at $\tau=1$ the sampler explores the more challenging target distribution. Markov Chains in both PT and ST therefore avoid becoming trapped in a single mode. In PT, T independent chains are run at different temperatures where information about $\boldsymbol{\theta}$ is allowed to flow between them. ST samples the system at different temperatures using a single chain by augmenting the state space with the temperature parameter as a dynamic variable. In both standard ST and PT, samples of interest correspond to samples obtained at $\tau=1$ i.e., $P(\boldsymbol{\theta}\mid\boldsymbol{Y},\tau=1)$.

Marginal likelihood estimates, which are crucial to Baye-sian model comparison and selection, can be obtained using thermodynamic integration (TI) based on the samples obtained in all T chains of PT (Friel and Pettitt, 2008; Calderhead and Girolami, 2009). TI approximates the log marginal likelihood, $\log{P(\boldsymbol{Y})}$, by numerically solving the integral $\int\limits_{0}^{1}E_{\boldsymbol{\theta}\mid\tau,\boldsymbol{Y}}[\log(P(\boldsymbol{Y}\mid\boldsymbol{\theta}))]d\tau$. As with any numerical integration scheme, accuracy of the integral depends on the discretization, here the number and location of $\tau$ values used in PT. Approximating the thermodynamic integral using PT approximations produces biased marginal likelihood estima-tes (Calderhead and Girolami, 2009). In this paper, we propose a new ST algorithm variant that also solves the Thermodynamic Integral without needing the user to select and tune the number and location of discrete $\tau$ values.

In order to update the inverse temperature parameter, standard ST requires the user to select a prior for $\tau$. Choice of $P(\tau)$ typically requires estimates of the normalizing constant, $z(\boldsymbol{Y}\mid\tau)=\int\limits_{\Theta}P(\boldsymbol{Y}\mid\boldsymbol{\theta})^{\tau}P(\boldsymbol{\theta})d\boldsymbol{\theta}$, which reduces to a finite dimensional problem if ST is performed over a fixed discretized sequence of $\tau$ values rather than treating $\tau$ as a continuous variable. Consequently, in standard ST, the temperature schedule and the normalizing constants must be estimated through preliminary runs.

Several other methods have been proposed for obtaining normalizing constants (Geyer and Thompson, 1995): iterative adjustment, Metropolis-Coupled Markov Chain Monte Carlo (MCMCMC) (Geyer, 1991), and stochastic approximation (Wasan, 1969; Wang and Landau, 2001). Stochastic approximation could be embedded within standard ST to adaptively estimate the normalizing constants, but this still requires preliminary runs to learn the temperature schedule (Atchade and Liu, 2004). Extending this further, the Parallel Adaptive Wang-Landau (PAWL) can be embedded in ST (Bornn et al, 2013; Bornn, 2014) to automatically learn the temperature schedule using an adaptive binning strategy. This approach is well suited to parameter estimation but does not resolve the discretization concerns when applied to TI.

Instead of performing preliminary runs to obtain the most suitable discretized temperature schedule and normalizing constants thereof, we propose a new ST algorithm with a continuous inverse temperature variable defined on $[0,1]$. We remove the requirement for calculating normalizing constants through imposition of the constraint that the distribution of the $\tau$ when profiling over $\boldsymbol{\theta}$ is Uniform. This constraint defines a formula for $P(\tau)$. In the rest of the paper we refer to the new ST algorithm as ‘Simulated Tempering Without Normalizing Constants’ (STWNC). By sampling across a continuous temperature scale, samples from STWNC can be used for thermodynamic integral estimates of the marginal likelihood.

In PT and ST, samples from the target distribution are retained whenever $\tau=1$ which is possible (though typically inefficient in standard ST) because of the discretization of the domain of $\tau$. However, in STWNC the temperature parameter is a continuous variable and therefore $P(\tau=1)=0$. To sample from the target distribution while maintaining a continuously valued $\tau$, we embed STWNC within the PT framework. The resulting PT and STWNC hybrid algorithm named ‘Parallel Tempering - Simulated Tempering Without Normalizing Constants’ (PT-STWNC), runs PT with $\text{T}=2$ chains, each targeting a different goal. The first chain draws samples at continuous temperatures via STWNC. We refer to the first chain as a ‘tempered’ chain. The second chain draws samples from the target distribution where $\tau=1$. We refer to the second chain as a ‘target’ chain in the rest of our presentation.

The remainder of the paper is organized as follows. Section 2 gives an overview of tempering methods with emphasis on PT and standard ST as well as a review of modern variants of stochastic approximation. Section 3 proposes the new STWNC algorithm and provides detailed explanation of how it removes the need for temperature dependent normalizing constants. In Section 4 the hybrid PT-STWNC, which embeds STWNC within a Parallel Tempering algorithm is proposed. Section 5 presents an overview of marginal likelihood approximation via thermodynamic integration and describes marginal likelihood estimation via STWNC. Sections 6 - 8 illustrate the PT-STWNC algorithm using several examples: a mixture model of two Gaussian distributions, a mixture model applied to Galaxy velocity data, a Susceptible-Infected-Recovered (SIR) epidemiological ordinary differential Equations (ODE) model. Section 9 provides discussion.

2 Tempering Methods

This section outlines Parallel and Simulated Tempering methods.

2.1 Parallel tempering

Parallel Tempering (PT), also known as replica exchange or Population MCMC, is a sampling algorithm designed to improve the dynamic properties of the MCMC samplers especially when it comes to exploring a posterior with many isolated modes.

Given the likelihood $P(\mathbf{Y}\mid\boldsymbol{\theta})$, and prior distribution $P(\theta)$, the inverse temperature sequence $0\leq\tau_{1}<,...,<\tau_{\text{T}}=1$ defines the T approximations to the target posterior distribution:

[TABLE]

where

[TABLE]

is the temperature dependent normalizing constant of the $t^{th}$ approximation to the target posterior distribution. The parameter $\tau$ controls the contribution of the likelihood to the equation (1), thus enabling the sampler to move easily and explore the parameter space when $\tau$ is low-valued and to remain further within the basin of attraction of a local mode when $\tau$ is high-valued.

At each iteration PT performs one of two steps: a mutation step where, for example, Metropolis-Hastings (MH) is used to jitter $\boldsymbol{\theta}$ from each of the T chains independently, and an exchange step where, with some probability two chains $t$ and $l$ are randomly chosen to exchange their parameters $\boldsymbol{\theta}_{t}$ and $\boldsymbol{\theta}_{l}$. The exchange is accepted with probability:

[TABLE]

Note the cancellation of normalizing constants $z(\mathbf{Y}\mid\tau_{l})$ and $z(\mathbf{Y}\mid\tau_{t})$ when swapping parameters between two different temperature based chains. The cancellation of normalizing constants and the ease of movement between modes at low values of $\tau$ has made PT into a widely used algorithm. Samples from the $\text{T}^{th}$ chain correspond to samples from the target distribution.

2.2 Simulated tempering

Simulated tempering (ST), also known as serial tempering, is a single chain sampling method where the posterior parameter space is augmented by including the temperature, $\tau$, as a random variable. Consequently ST requires specification of the prior $P(\tau)$. As with PT, $\tau$ controls the influence of the likelihood on $P(\boldsymbol{\theta}\mid\mathbf{Y},\tau)$. However, in ST different temperatures are explored in a random walk through the joint distribution of $\boldsymbol{\theta}$ and $\tau$ :

[TABLE]

where the normalizing constant for the joint density is:

[TABLE]

Similar to PT, samples that arise at $\tau=1$ are samples from target distribution, however, when $\tau$ is continuous, $P(\tau=1)=0$.

2.2.1 Prior for $\tau$

In order for standard ST to mix well, a carefully chosen prior $P(\tau)$ needs to be defined. Following Geyer and Thompson (1995), the prior for the inverse temperature can be found by examining the marginal distribution of $\tau$:

[TABLE]

where $z(\boldsymbol{Y}\mid\tau)=\int\limits_{\boldsymbol{\Theta}}P(\boldsymbol{Y}\mid\boldsymbol{\theta})^{\tau}P(\boldsymbol{\theta})d\boldsymbol{\theta}$ is the $\tau$ dependent normalizing constant of the conditional posterior distribution $P(\boldsymbol{\theta}\mid\tau,\boldsymbol{Y})$. If the prior for $\tau$ is chosen to be approximately proportional to the inverse normalizing constant, that is, if

[TABLE]

then the marginal distribution of $\tau\mid\boldsymbol{Y}\sim\mbox{Uniform(0,1)}$. Following Geyer and Thompson (1995), the standard ST algorithm starts with $n$ tempered distributions at $\tau_{i}=\frac{i}{T}$ where $i=0,..,n\leq T$, and iterates between adjusting the prior of $\tau$ and adjusting the inverse temperature spacing until a desired rate for transitions between the interpolating distributions is met. Afterwards, new inverse temperatures are added, and a new iterative cycle of adjusting prior and inverse temperature spacing is started for the newly added temperature set.

2.3 Outline of the Standard Simulated Tempering

A Markov Chain using standard ST updates parameters through Gibbs steps, updating $\boldsymbol{\theta}\mid\tau,\boldsymbol{Y}$ and $\tau\mid\boldsymbol{\theta},\boldsymbol{Y}$ in turn. In the first kind of transition, a fixed temperature mutation step is typically a Metropolis Hastings step identical to what would be used to sample one of the fixed temperature chains $P_{t}(\boldsymbol{\theta}\mid\boldsymbol{Y})$ in PT.

Updating $\tau\mid\boldsymbol{\theta}^{(i+1)},\boldsymbol{Y}$ through Metropolis Hastings occurs by proposing a value $\tau^{*}$ sampled from a symmetric distribution (for expositional simplicity). Using the prior in (5), the value $\tau^{*}$ is accepted (setting $\tau^{(i+1)}=\tau^{*})$ with probability:

[TABLE]

Consequently, the acceptance probability relies on temperature dependent normalizing constant $z(\boldsymbol{Y}\mid\tau)$. Normalizing constants, also referred to as weights of the sequence of distributions (Neal, 1996), are generally unknown and finding suitable values requires pilot runs – defeating the purpose of using standard ST. Information from pilot runs also cannot be recycled when additional data becomes available because $P(\tau)$ depends on $\boldsymbol{Y}$. To make the normalizing constant problem simpler, $\tau$ is discretized instead of being continuously valued.

2.4 Related algorithms

The Wang-Landau (WL) algorithm (Wang and Landau, 2001), a stochastic approximation algorithm, could be used within standard ST to automatically obtain $P(\tau)$ (Geyer and Thompson, 1995; Atchade and Liu, 2004). Standard WL algorithm requires a predefined temperature schedule $0\leq\tau_{1}<,..,<\tau_{\text{T}}=1$ corresponding to partitioning the temperature state space $\mathbb{T}$ into T different regions i.e., bins $\text{E}_{1},..,\text{E}_{\text{T}}$. The goal is to construct a chain that could spend the same time in each $\text{E}_{t}$. The moves inside the $\text{E}_{t}$ are performed with a standard Metropolis-Hastings (MH) algorithm with target distribution $\pi$. The WL algorithm recursively re-weights $\pi$ in $\text{E}_{t}$ by a factor $\phi_{i}(t)$. Hence, given $\tau^{(i)}$ and some unnormalized weights $\phi_{i}$, the $\tau^{(i+1)}$ can be sampled using a MH algorithm with invariant density proportional to $\sum\limits_{t=1}^{\text{T}}\frac{\pi(\tau)}{\phi_{i}(t)}\mathbb{I}_{\text{E}_{t}}(\tau)$, where $\mathbb{I}$ is the indicator function. The normalizing constants at the $i$-th iteration $\phi_{i}(t)=z(\textbf{Y}\mid\tau_{t}^{(i)})$ are updated for each of the bins $t\in\{1,..,\text{T}\}$ until a predefined criterion is met. This criterion, also known as ’flat histogram’, ensures that proportions of visits of the chain in each of the bins $\text{E}_{1},..,\text{E}_{\text{T}}$ are approximately equal to $\text{T}^{-1}$. The updating rule for the normalizing constants yields $\phi_{i+1}(t)=\phi_{i}(t)(1+\gamma_{i}\mathbb{I}_{\text{E}_{t}}(\tau^{(i+1)}))$, where $\gamma_{i}$ is a learning rate which decreases stochastically until the criterion of ’flat histogram’ of visit frequency to the bins $\text{E}_{1},..,\text{E}_{\text{T}}$ is met.

The Parallel Adaptive WL (PAWL) algorithm (Bornn et al, 2013) removes the need for preliminary runs of WL to learn the optimal partitioning of the state space by exploiting an adaptive binning strategy. The adaptive binning strategy requires initial bins and bin range to be specified by the user. Since $\phi_{i}(t)$ represent normalizing constant for the $t$-th bin, the adaptive binning strategy maintains uniformity within a bin to allow within-bin movement. This is achieved by determining presence of heavy tails in the distribution of the samples within each bin. If the distribution is skewed towards the left side, then the sampler will have difficulty moving to the neighboring bin on the left. Hence, the binning strategy divides the bin into two chains by the middle point of the bin, and then it measures discrepancy between the chains using a ratio of the number of points in any of the two chains and the number of points within the bin. If this ratio is close to $50\%$ then the histogram of the within-bin distribution is close to uniform. Hence, if the ratio is below some threshold, for example $25\%$, two new bins are created using the middle point of the former bin, and otherwise, the bin remains unchanged. One can specify the threshold to be $50\%$, but then the number of newly created bins will be larger. PAWL checks if the bins have to be split until the ’flat histogram’ criterion is met, and afterwards the bin splitting stops since the sampler can move easily between the bins. Embedding PAWL within standard ST has been explored with a purpose to automate the two input requirements for ST: choice of temperature schedule and calculation of $P(\tau)$ (Bornn, 2014).

The Equi-Energy (EE) sampler (Kou et al, 2006), which utilizes temperature-energy duality, also allows wide moves by performing jumps between the states with similar energy levels. The EE, which is a powerful sampling and estimation methodology that addresses multi-modality in high-dimensional target distributions, provides estimates of expectations under any fixed temperature. However, the discrete nature of the temperature in EE sampler requires careful tuning of the temperature schedule in order for the EE samples to be applicable to thermodynamic integration.

3 Simulated Tempering Without Normalizing Constants

The standard ST is not widely used in practice because of the challenges that arise from finding suitable values of the unknown normalizing constants. The proposed Simulated Tempering Without Normalizing Constants (STWNC) algorithm removes the dependence on normalizing constants in the acceptance ratio in (6) by the way we define the prior for $\tau$. As with ST, STWNC is still moving through two kinds of transitions: updating $(\boldsymbol{\theta}\mid\tau,\boldsymbol{Y})$; and updating $(\tau\mid\boldsymbol{\theta},\boldsymbol{Y})$. The nuisance parameter $\tau$ is not of inferential interest, consequently its prior $P(\tau)$ can be selected for algorithmic convenience. Our algorithm therefore chooses $P(\tau$) by imposing a constraint that the profile posterior distribution of $\tau$ while maintaining $\boldsymbol{\theta}\mid\tau,\boldsymbol{Y}$ at its maximum value results in a uniform distribution. Using this constraint we derive a formula for the prior of $\tau$ which is computationally inexpensive and does not require preliminary runs of the algorithm or discretization of $\tau$. In the remainder of this section we describe the derivation of formula for the prior of $\tau$.

We impose the constraint that the posterior distribution of $\tau$ while profiling over $\boldsymbol{\theta}$ is uniform:

[TABLE]

[TABLE]

Equation (8) implies that the joint posterior distribution of $\tau$ and ${\theta}_{\max}(\tau)$, $P(\tau,\boldsymbol{\theta}_{\max}(\tau)\mid\boldsymbol{Y})$ has a ridge of constant maximum height from $\tau=0$ to $\tau=1$. In other words for any $\tau_{i},\tau_{j}\in[0,1]$,

[TABLE]

To derive the prior of $\tau$, we expand the profile posterior distribution given by the equation (8):

[TABLE]

where $P(\boldsymbol{Y}\mid\boldsymbol{\theta}_{\max}(\tau))^{\tau}$ is tempered profile likelihood, $P_{\text{prof}}(\boldsymbol{Y})=\int\limits_{0}^{1}P(\boldsymbol{Y}\mid\boldsymbol{\theta}_{\max}(\tau))^{\tau}P(\boldsymbol{\theta}_{\max}(\tau))P(\tau)d\tau$ is the normalizing constant of the profile posterior distribution of $\tau$ over $\boldsymbol{\theta}$, $P(\tau)$ is the prior of the inverse temperature and $P(\boldsymbol{\theta}_{\max}(\tau))$ is prior of $\boldsymbol{\theta}$ evaluated at $\boldsymbol{\theta}_{\max}(\tau)$. Expressing $P(\tau)$ by rearranging (LABEL:EQ:joint) gives

[TABLE]

As a direct consequence of the constraint that profile posterior distribution of $\tau$ while maintaining $\boldsymbol{\theta}\mid\tau,\boldsymbol{Y}$ at its maximum value is uniform, i.e., (8) and (9), the numerator in the formula for prior of $\tau$ in (11) is constant with respect to $\tau$.

Using the equation (11) as prior for $\tau$, the STWNC acceptance ratio for a proposed $\tau^{*}$ alters (6) into:

[TABLE]

The acceptance ratio in equation (LABEL:eq:Updtau) does not depend on the temperature dependent normalizing constant $z(\boldsymbol{Y}\mid\tau)$ thus eliminating the need for discrete temperatures, normalizing constant estimates, and tuning of bin widths.

The formula for the prior of $\tau$ imposes the property that the maximal contour of the joint posterior distribution is continuous for $\tau\in[0,1]$. However, calculation of $\boldsymbol{\theta}_{\max}(\tau^{*})$ requires optimization for each proposed value of $\tau^{*}$. While this initially may be prohibitive, the optimization cost decays with increasing iterations because of the improved starting points available for the optimizer. At iteration $i$, with proposed $\tau^{*}$ we initialize the optimizer at the point $\boldsymbol{\theta}_{\max}(\tau^{(k)})$ for the $k<i$ that minimizes: $||\tau^{(k)}-\tau^{*}||$. As $i\rightarrow\infty$ the optimization is initialized arbitrarily close to the optimum. The computational cost of evaluating the prior for $\tau$ is reduced, but some computational overhead remains by instead searching through previously obtained values for initialization. This approach was implemented in the examples in this paper and is incorporated into the Pseudo-code of STWNC is given in Algorithm 1.

For long running Markov Chains, further computational reduction may be possible within the MCMC iterations by approximating $P(\tau)$ with an iteratively updated Gaussian process approximation (Conrad et al, 2016), or by approximating the manifold of $\theta_{\max}(\tau)$ based on offline evaluations thereof. In this latter approach, before beginning MCMC iterations, optimization of $\theta(\tau)$ is performed over a grid of values of $\tau$. An interpolator is defined so as to replace the optimization of $\theta(\tau)$ within each MCMC iteration with evaluation of the interpolator instead. After this initial computational cost, the time per iteration is equivalent to that of parallel tempering with the same number of parallel chains. This approach is explored as proof of concept and not in a general sense because the precision of the approximation depends on the quality of the interpolator which in turn will be impacted by smoothness and continuity of the manifold, properties which are difficult to characterize or guarantee given the multi-modal nature of the likelihoods in this manuscript. This latter strategy is demonstrated as proof of concept in Section 8.2.

4 Parallel Tempering via Simulated Tempering without Normalizing Constants algorithm

In Parallel Tempering via Simulated Tempering without Normalizing Constants (PT-STWNC), the inverse temperature is a continuous parameter and $P(\boldsymbol{\theta}\mid\boldsymbol{Y})=\int\limits_{0}^{1}P(\boldsymbol{\theta},\tau\mid\boldsymbol{Y})d\tau$ is a marginal distribution, which does not coincide with the target distribution. Instead, the target distribution is $P(\boldsymbol{\theta}\mid\boldsymbol{Y},\tau=1)$. To obtain samples from the target distribution we run PT with $T=2$ defined as follows:

• •

The first PT chain, called the tempered chain, updates $\boldsymbol{\theta}$ and $\tau$ via STWNC using $P(\boldsymbol{\theta},\tau\mid\boldsymbol{Y})$

• •

The second PT chain, called the target chain, updates $\boldsymbol{\theta}$ via standard Metropolis-Hastings or Gibbs sampling from the target distribution $P(\boldsymbol{\theta}\mid\boldsymbol{Y},\tau=1)$.

As in standard PT, the PT-STWNC algorithm goes through two transitions: exchange and mutation. In the exchange step, the exchange between the two chains is proposed. If the exchange was accepted, then the two chains swap the parameter values between each other, and otherwise the two chains go through mutation steps. In the mutation step, the first chain updates the parameters of interest at different temperatures via STWNC; the second chain updates the parameters of interest at $\tau=1$ via Metropolis-Hastings. Pseudo-code of PT-STWNC is given in Algorithm 2.

5 Estimation of Marginal Likelihoods via thermodynamic integration

Posterior model probabilities $P(M\mid\boldsymbol{Y})$ provide an intuitive framework for evaluating model $M$ within a model class. Expanding the model class requires rescaling all posterior model probabilities. Consequently, comparing models $M_{1}$ and $M_{2}$ is typically performed through the posterior odds,

[TABLE]

The ratio of posterior and prior odds,

[TABLE]

is Bayes Factor of $\text{M}_{1}$ against $\text{M}_{2}$ (Kass and Raftery, 1995).

When there are no prior preferences for models, $\text{B}_{12}$ is equal to the posterior odds. The marginal likelihoods for models $\text{M}_{j},j=1,2$, in (15), are obtained by integrating over the parameter space,

[TABLE]

where $\boldsymbol{\theta}_{j}$ are the parameters corresponding to the $j$-th model. For expositional simplicity we will assume no prior preference for models throughout this paper. Computing meaningful Bayes Factors requires accurate estimates of the marginal likelihood in (16).

The Posterior Harmonic Mean estimator (PHM), uses importance sampling to integrate (16) (Newton and Raftery, 1994; Raftery et al, 2006) resulting in unbiased marginal likelihood estimates but potentially infinite variance. Steppingstone sampling (SS) (Xie et al, 2011), which uses importance sampling to estimate each ratio of normalizing constants of the sequence of interpolating distributions between prior and target distribution, provides reliable marginal likelihood estimates. Alternatively, thermodynamic integration (TI) (Friel and Pettitt, 2008) builds on ideas from path sampling (Gelman and Meng, 1998) to estimate the marginal likelihood via,

[TABLE]

where the expectation in the integrand is with respect to the tempered posterior distribution in (1) (Friel and Pettitt, 2008; Calderhead and Girolami, 2009). A numerical approximation to the thermodynamic integral in (17) is possible through discretization of $\tau$. In Friel and Pettitt (2008), samples from the discretized tempered posteriors were used from parallel chains in PT. At each discretized value of $\tau$,

$\operatorname{\mathbb{E}}_{\boldsymbol{\theta}\mid\boldsymbol{Y},\tau,\text{M}_{j}}[\log(P(\boldsymbol{Y}\mid\boldsymbol{\theta},\text{M}_{j}))]$ is evaluated and the marginal likelihood is approximated by applying a trapezoid rule to numerically integrate over $\tau$,

[TABLE]

where $\operatorname{\mathbb{E}}_{t,\text{M}_{j}}=\operatorname{\mathbb{E}}_{\boldsymbol{\theta}\mid\boldsymbol{Y},\tau_{t},{\text{M}_{j}}}[\log(P(\boldsymbol{Y}\mid\boldsymbol{\theta},\text{M}_{j}))]$ and $\Delta\tau_{t}=\tau_{t}-\tau_{t-1}$. The discretized trapezoidal rule in (18) was improved by Calderhead and Girolami (2009) by correcting for integration bias in terms of Kullberg-Leibler divergence, $\text{KL}(p_{t-1,{\text{M}_{j}}}\|p_{t,{\text{M}_{j}}})$, of $p_{t,{\text{M}_{j}}}$ from $p_{t-1,{\text{M}_{j}}}$ for model $\text{M}_{j}$,

[TABLE]

where $p_{t,{\text{M}_{j}}}$ is the tempered posterior distribution for the model $\text{M}_{j}$ given by the equation (1).

The thermodynamic integration via PT applies a numerical integration approximation to Monte Carlo approximations. The approximation error should decrease with number of chains and number of samples in each chain. As with any numerical integration, discretization over $\tau$ determines the accuracy of the result. Determining the optimal temperature schedule requires preliminary experimentation which contributes to unpopularity of the thermodynamic integration in practice. Based on the idea from path sampling, (Gelman and Meng, 1998), Calderhead and Girolami (2009) proposed that the temperature schedule could be chosen such that the Monte Carlo variance of the marginal likelihood estimates is minimized. The authors’ numerical simulations suggest that the optimal temperature schedule is $\tau_{i}=(\frac{i}{\text{T}})^{5}$ in their situations, which puts more emphasis on values closer to the prior.

5.1 Computing the marginal likelihood via STWNC

Following Friel and Pettitt (2008), in standard ST, samples $\{(\boldsymbol{\theta}^{1},\tau_{1}),..,(\boldsymbol{\theta}^{n},\tau_{n})\}$ drawn from $P(\boldsymbol{\theta},\tau\mid\boldsymbol{Y},\text{M}_{j})$ can be used to estimate the marginal likelihood by first obtaining Monte Carlo approximation $\operatorname{\mathbb{E}}_{\boldsymbol{\theta}\mid\boldsymbol{Y},\tau,\text{M}_{j}}[\log{P(\boldsymbol{Y}\mid\boldsymbol{\theta},\text{M}_{j})}]$, and then solving the thermodynamic integral in (18) via quadrature. This is based on the assumption that the prior of $\tau$ is proportional to the temperature-dependent normalizing constant, $P(\tau)\propto z(\boldsymbol{Y}\mid\tau,\text{M}_{j})$. According to Friel and Pettitt (2008), in single chain methods such as ST, the normalizing constant $z(\boldsymbol{Y}\mid\tau,\text{M}_{j})$ varies by orders of magnitude with $\tau$ which leads to poor estimation of the log untempered likelihood $\log{P(\boldsymbol{Y}\mid\boldsymbol{\theta},\text{M}_{j})}$. Thus in standard ST, small values of $\tau$ do not tend to be sampled with high frequencies. However, in STWNC, where $\tau$ is continuous, the marginal distribution of $\tau$ tends to spend a lot of time at near zero values The marginal distribution of $\tau$ in STWNC coincides with that of the recommended geometric temperature schedule for thermodynamic integration via PT suggested by Calderhead and Girolami (2009).

Samples $\{(\boldsymbol{\theta}^{1},\tau_{1}),..,(\boldsymbol{\theta}^{n},\tau_{n})\}$ from PT-STWNC can be used to solve the marginal likelihood integral in (17), by first ordering the samples with respect to $\tau$, and solving the integral numerically,

[TABLE]

where $\operatorname{\mathbb{E}}_{t,\text{M}_{j}}=\operatorname{\mathbb{E}}_{\boldsymbol{\theta}\mid\boldsymbol{Y},\tau_{t},{\text{M}_{j}}}[\log(P(\boldsymbol{Y}\mid\boldsymbol{\theta},\text{M}_{j}))]$ and $\Delta\tau_{t}=\tau_{t}-\tau_{t-1}$.

The variance of the marginal likelihood estimator in (20) asymptotically disappears with the number of samples, N. The variance estimator is:

[TABLE]

If all the samples are unique then $\text{T}=\text{N}$ and $\Delta\tau_{t}\approx\frac{1}{\text{N}}$ giving

[TABLE]

Furthermore, $\operatorname{\mathbb{E}}_{t,{\text{M}_{j}}}$ is a single sample from $P(\boldsymbol{\theta}\mid\boldsymbol{Y},\tau_{t},{\text{M}_{j}})$ giving the bounds:

[TABLE]

where upper and lower bounds go to zero as $\text{N}\rightarrow\infty$.

The TI via PT relies on two layers of approximation to produce marginal likelihood estimates. The first layer of approximation in TI via PT corresponds to the MCMC integration (i.e., obtaining PT samples), and the second layer occurs when $\tau$ is numerically integrated out from the thermodynamic integrals in (18) and (19). Similarly, in the TI via PT-STWNC, the first approximation layer corresponds to obtaining samples from PT-STWNC, and the second layer of approximation is a result of solving the marginal likelihood integral in (20) numerically. However, since the $\tau$ is continuous, the TI via PT-STWNC removes the need for optimal temperature discretization schedule. Evaluations of the log likelihood can be saved within Algorithm 1 prior to applying the temperature and consequently PT-STWNC provides marginal likelihood estimate with negligible additional computational cost.

6 Example: Bimodal model

The likelihood is bimodal with respect to $\mu$, but is unimodal with respect to $\boldsymbol{Y}$,

[TABLE]

and the posterior distribution of P($\mu\mid\boldsymbol{Y},\sigma^{2}$) is bimodal. The data were simulated from (24) with $n=25$, $\mu=1.5$ and $\sigma^{2}=1$.

The sampling distribution of the PT-STWNC is the tempered joint posterior distribution,

[TABLE]

Conjugate priors on $\mu$ and $\sigma^{2}$ were assigned, $P(\mu)\sim\mbox{N}(0,1)$; $P(\sigma^{2})\sim\mbox{InverseGamma}(1,1)$.

6.1 Results

The PT-STWNC algorithm was run for 50,000 iterations with the first 15,000 samples discarded as burn-in. Figure 1 shows the sampled joint posterior distribution of $\mu$ and $\tau$ obtained from the PT-STWNC ‘tempered’ chain. The perspective and contour plots in Figure 1 illustrate uniform profile as seen in the two ridges of the posterior surface having approximately constant maximum height along the $\tau$ axis. The last observation is a direct consequence of the constraint that the profile posterior distribution of $\tau$ while profiling over $\mu$ is uniform on the interval [0,1].

The wide contours of Figure 1B at low values of $\tau$ demonstrate that the ’tempered’ chain spends a lot of time sampling at low values of the inverse temperature, thus taking large steps to move between the two modes. Similarly, the marginal distribution of sampled $\tau$ demonstrates that low values of $\tau$ are sampled with higher frequencies (see Figure 2, gray color). The PT-STWNC ‘target’ chain updates the parameters of interest at $\tau=1$, which results in drawing samples from the target posterior distribution (Figure 2, diagonal, blue color).

PAWL within ST was also run on this toy example and the marginal distributions of $\mu$ and $\tau$ were compared to those from the PT-STWNC. The plot of the marginal distribution of $\mu$ is similar to the marginal distribution of $\mu$ obtained from the ’tempered’ PT-STWNC chain (Figure 2, green color). The marginal distribution of the discrete $\tau$ from the PAWL within ST has similar shape as that of the continuous $\tau$ from the PT-STWNC. Table 1 shows the posterior means of $\mu$ and $\sigma^{2}$ compared to their theoretical values. Full implementation details are given in the Appendix 10.2.

6.2 Marginal likelihood estimation

The PT-STWNC was used to estimate the marginal likelihood of the bimodal model introduced in the Section 6. The PT-STWNC marginal likelihood estimates were compared to the following three approaches: $i).$ analytical solution to the marginal likelihood integral in (16), $ii).$ the thermodynamic integration via PT with bias correction by Calderhead and Girolami (2009) (TI-PT-B) and $iii).$ the thermodynamic integration via PT without bias correction by Friel and Pettitt (2008) (TI-PT-NB).

This is a toy example and one of the very few cases where a closed form of the marginal likelihood integral exists (see Appendix 10.1). The PT-STWNC marginal likelihood estimates were obtained directly from the ’tempered’ and ’target’ chain using the equation (20). TI via PT estimates were obtained by running PT with T=30 chains. The inverse temperature schedule was chosen as the geometric schedule, $t_{i}=(\frac{i}{\text{T}})^{5}$.

The TI-PT-NB and the TI-PT-B estimates in Table 2, were obtained using samples from the same run of the PT. The only difference between the TI-PT-NB and the TI-PT-B is the bias term as per equation (19). The results in Table 2 show that the thermodynamic integral bias term has a small effect on the marginal likelihood estimate. Table 2 also found that increased number of PT chains (T=60 and T=100) did not yield substantially better TI-PT-NB and TI-PT-B estimates. This result complies with the finding by Ahlers and Engel (2008), who demonstrated that while thermodynamic integration via PT performs well in unimodal case, the method exhibits substantiative bias in a bimodal case with a tractable marginal likelihood. The results from PT and our PT-STWNC both exhibit bias, this is consistent with Ahlers and Engel (2008), who also found that the estimates do not improve with increased number of parallel chains. Those authors trace back this problem to the incomplete equilibrium between the two modes, which leads to failure to reproduce the exact mixture probabilities. In addition, the bias in PT-STWNC, TI-PT-NB, and TI-PT-B did not reduce when the integral in (20) was obtained by combining all the samples from the 20 replicate runs.

7 Example: Galaxy data

The Galaxy data comprises velocities of $82$ galaxies that diverge from our galaxy studied by Postman et al (1986); Carlin and Chib (1995); Neal (1999). The data are univariate identically and independently distributed samples from mixture of $K$ Gaussian components and denoted as $\boldsymbol{Y}=(y_{1},y_{2},..,y_{n})^{{}^{\prime}}$, with $n=82$. Parameters of the model are given as $\boldsymbol{\theta}=\left(\boldsymbol{\mu},\boldsymbol{\sigma^{2}},\boldsymbol{p}\right)^{{}^{\prime}}$ where $\boldsymbol{\mu}=\left(\mu_{1},..,\mu_{K}\right)^{{}^{\prime}}$ is a vector of mixture component means, $\boldsymbol{\sigma^{2}}=\left(\sigma^{2}_{1},..,\sigma^{2}_{\text{K}}\right)^{{}^{\prime}}$ is a vector of mixture component variances and $\boldsymbol{p}=\left(p_{1},..,p_{K}\right)^{{}^{\prime}}$ is a vector of mixture probabilities. The $k$-th mixture component has distribution $\mbox{N}(\textbf{Y}\mid\mu_{k},\sigma_{k}^{2})$. Then likelihood function is

[TABLE]

Conjugate priors were assigned, $P(\mu_{k})\sim\mbox{N}(20,100),$ $P(\sigma_{k}^{2})\sim\mbox{IGamma}(shape=3,scale=20)$ and $P(p_{1},..,p_{K})\sim\mbox{Dirichlet}(\alpha_{1}=1,..,\alpha_{K}=1)$.

A latent variable Z such that $P(Z_{ik}=1\mid p_{k})=p_{k}$ was introduced to help derive the necessary sampling distributions. Z is a $n\times K$ matrix of indicator variables in which $\text{Z}_{i,k}=1$ indicates the data point $i$ belongs to the mixture component $k$. Then likelihood of data point $\{y_{i}\}$ conditional on $\text{Z}_{ik}$ is,

[TABLE]

and the joint distribution of $\{y_{i}\}$ and $Z_{ik}$ is,

[TABLE]

Each row of Z is a multinomial variable with probabilities,

[TABLE]

The joint posterior distribution of the parameters of interest $\boldsymbol{\theta}=\left(\boldsymbol{\mu},\boldsymbol{\sigma^{2}},\boldsymbol{p}\right)^{{}^{\prime}}$, latent variable Z and the inverse temperature parameter $\tau$ is,

[TABLE]

7.1 Results

The PT-STWNC algorithm was run for 35,000 iterations with the first 1000 generated samples being removed as burn-in. Multi-modality in the Galaxy data with three components is illustrated by marginal distributions of $\mu_{1},\mu_{3},\sigma^{2}_{1}\mbox{ and }\tau$ and bivariate joint posterior distribution of $(\mu_{1},\mu_{3})$ (Figure 3). Plot of the marginal distribution of $\tau$ shows that $\tau$ spends a lot of time at values close to zero, which is crucial for the sampler to move easily between the isolated modes and produce accurate estimates (Figure  3).

7.2 Marginal likelihood estimation for the Galaxy data

The PT-STWNC was used to perform model selection on the Galaxy data set via the choice of $K$. The marginal likelihood was calculated for the Galaxy data to compare the following models: two components with equal variances, three components with unequal variances, three components with equal variances, four components with unequal variances and five components with unequal variances. A closed form expression of the marginal likelihood integral is not available.

The two thermodynamic integration approaches from the Section 5, were used to compare model selection abilities with PT-STWNC. Log marginal likelihood estimates and Bayes Factors demonstrate that all three estimation techniques agree that the model with five components is the best model (Tables 3 and 4). In addition, all the three estimation methods find support for the models with $3-5$ components, while the worst model is the model with $2$ components and equal variances. Findings from our model selection study comply with the results from previous studies. For instance, Steele and Raftery (2010) argued that the number of components in the Galaxy data is not known, and concluded that the models with $3$ and $6$ components are reasonable fit to the data. Gibbs sampling was used to show that the model with $3$ components is the best model Chib (1995). Evolutionary Monte Carlo was used in combination with bridge sampling to demonstrate that the model with 5 components is the best model and $2$ components with equal variances is the worst model Liang and Wong (2001). In addition, Liang and Wong (2001) found support for the models with $3-5$ components. Reversible Jump MCMC (RJMCMC) was used to conclude that the models with $5-7$ components are good fit to the data Richardson and Green (1997).

The TI-PT-NB and the TI-PT-B estimates in Table 3 were obtained using samples from the same run of the PT. Hence, the only difference between the TI-PT-NB and the TI-PT-B is in the bias term as per equation (19). The results in the Table 3 show that the estimates from the TI-PT-NB and the TI-PT-B are nearly the same, which suggests that the thermodynamic integral bias term has near zero effect in this example.

8 Example: Susceptible-Infected-Recovered (SIR) epidemiological model with real data

We illustrate the PT-STWNC on a Susceptible-Infected-Recovered (SIR) epidemiological model for number of daily deaths due to the black plague. The data were collected by the grave digger during the second black plague outbreak in the village of Eyam, UK, from June 19, 1666 to November 1, 1666 (Massad et al, 2004). The village had quarantined itself to avoid spreading the disease to the neighboring villages. Therefore, the population size is fixed to $\text{N}=261$, and the population is stratified into groups of susceptible $\text{S}(t)$, infected $\text{I}(t)$ and removed $\text{R}(t)$ individuals, $\text{N}=\text{S}(t)+\text{I}(t)+\text{R}(t)$. Since there is no recovery from the plague, the number of deaths correspond to the number of removed individuals up to time $t$, $\text{R}(t)$ (Campbell and Lele, 2014; Golchi and Campbell, 2016)

The disease spread dynamics can be described by the system of ordinary differential equations (ODE),

[TABLE]

where $\alpha$ describes the rate of death once the individual is infected and $\beta$ describes the plague transmission. Additional to the model parameters $\boldsymbol{\theta}=\left(\alpha,\beta\right)^{{}^{\prime}}$, the ordinary differential equations model requires estimates of the initial states $(\text{S}(0),\text{I}(0),\text{R}(0))^{{}^{\prime}}$. At the initial time the population consists of susceptible and infected individuals, and therefore $\text{R}(0)=0$ and $\text{S}(0)=\text{N}-\text{I}(0)$. Consequently, the only initial state parameter is $\text{I}(0)$ so that the unknown parameters of the model are $\boldsymbol{\theta}=\left(\alpha,\beta,\text{I}(0)\right)^{{}^{\prime}}$. The data denoted as $\boldsymbol{Y}=(y_{1},..,y_{n})^{{}^{\prime}}$ with $n=136$, represent cumulative number of deaths up to times $t_{1},..,t_{n}$. The data points $\boldsymbol{Y}$ were modeled by a Binomial distribution with expected value equal to the solution to the system (27), $R_{(\alpha,\beta,\text{I}(0))}(t)$, where $t\in\{t_{1},..,t_{n}\}$ .

The number of susceptible $\text{S}(t)$ and infected $\text{I}(t)$ are not observed. However, the number of infected at the end of the plague is 0, and the number of infected at time one before the end of the plague must therefore equal 1 (Campbell and Lele, 2014). Two additional data points on number of infected individuals $\boldsymbol{X}=(x_{n-1}=1,x_{n}=0)^{{}^{\prime}}$ at times $(t_{n-1},t_{n})^{{}^{\prime}}$ were modeled using Binomial distribution with expected value equal to $\text{I}_{(\alpha,\beta,\text{I}(0))}(t)$ for $t\in(t_{n-1},t_{n})^{{}^{\prime}}$,

[TABLE]

Prior distributions for $\boldsymbol{\theta}=\left(\alpha,\beta,\text{I}(0)\right)^{{}^{\prime}}$ were chosen to be: $\alpha,\beta\sim\mbox{Gamma}(1,1),$ $\text{I}(0)\sim\mbox{Binomial}(\text{N},\frac{5}{\text{N}}).$

Parameters $\alpha$ and $\beta$ are continuous and $\text{I}(0)$ is discrete. The discrete nature of the $\text{I}(0)$ induces multi-modality in the likelihood surface. This mixture of discrete and continuous parameters in the model imposes difficulties in sampling from the posterior distribution. Standard MCMC could get easily trapped in local modes of the posterior of the parameters of interest.

8.1 Results

The PT-STWNC was run on the SIR model for 35,000 iterations with 1000 burn-in samples. We applied a conjugate-gradient optimization routine for continuous parameters in parallel conditional on discrete parameters allowing a follow-up line search over discrete parameters. Multi-modality and topological challenges of the model are illustrated by the marginal distribution plots (Figure 4, diagonal) and by the bivariate joint posterior distributions plots (Figure 4, off-diagonal) of $\boldsymbol{\theta}=\left(\alpha,\beta,\text{I}(0)\right)^{{}^{\prime}}$. The marginal distributions of $\alpha$ and $\beta$ exhibit structures with three isolated modes. In Figure 4 (off-diagonal), clouds in the joint posterior distribution of $\alpha$ and $\beta$ represent the modes which correspond to the discrete samples of $\text{I}(0)=\{6,5,4,3\}$ from left to right.

Histograms of the marginal distributions of $\alpha,\beta,\text{I}(0)$ and $\tau$ obtained from the ‘tempered’ chain (Figure 5), illustrate the complexity and topological challenges of the model as well as the need for exploring the diffuse prior parameter space in order for PT-STWNC to draw samples from the target distribution. Figure 6 demonstrates that the prior parameter space (the grey contour lines) is much more diffuse than that of the joint posterior distribution of $\alpha$ and $\beta$ (the red dots, which when zoomed-in assumes the shape of the target posterior distribution). Consequently, the algorithm spends much time sampling at near-zero values of $\tau$ thus exploring the prior parameter space.

For implementation details we refer the reader to Appendix 10.3.

8.2 Computational Acceleration by Approximating $P(\tau)$

A promising direction for accelerating the implementation is through optimal manifold approximation. Prior to beginning MCMC iterations, optimization was performed over a grid of 301 values of $\tau$ and an interpolator is defined so as to replace the optimization within MCMC iterations with evaluation of the interpolator instead. In this example, $\theta_{\max}(\tau)$ is optimized over a $\log_{10}$-uniform grid of 301 values of $\tau\in[10^{-15},1]$. The resulting values are then interpolated using a $4^{th}$ order b-spline basis with 60 $\log_{10}$-uniformly spaced knots across $[10^{-15},1]$. The computational time for the initial 301 optimizations and spline interpolator set up was under one minute. After that initial computational cost, the time per iteration of this approximation to PT-STWNC is equivalent to that of parallel tempering also operating with the same number of chains. This is explored only as proof of concept because the degree of the approximation depends on the interpolator which in turn will be impacted by smoothness and continuity of the manifold which is not guaranteed in a general sense for multi-modal problems.

8.3 Comparison with PT and Complexity

Parallel tempering with 5 chains was run alongside PT-STWNC on the model from section 8 so as to compare the algorithms. The rule of thumb spacing for temperatures was used $\tau_{t}=(t/5)^{5}$ for $t\in\{1,\ldots,5\}$ following (Calderhead and Girolami, 2009), although using fewer temperatures than would be required for low bias thermodynamic integration. Both PT-STWNC and PT resulted in similar point and modal importance estimates with the latter values differing by at most 5% across the 4 modes. In both cases the lag 1 autocorrelation for the target chain was driven almost entirely by which mode it was sampling from and as a result in both methods that value was $\sim 0.97$ but differed slightly in the third decimal place.

When the optimum values are tracked, ST an additional diagnostic which is unavailable to PT in that posterior samples should be in the vicinity of the optimal value for all $\tau\in[0,1]$ if the chain is sampling correctly. Figure 7 shows the samples of $\beta$ and $\alpha$ plotted against $\tau$ for the STWNC chain. This figure includes the optimum for the profile at each value of $\tau$. The optimum line moves slowly with decreasing $\tau$ but curves dramatically near zero values as the priors (both exponential with mean 1) over-ride the likelihood.

9 Discussion

In this paper we presented a solution to the implementation challenge that has left Simulated Tempering unuseable as a Monte Carlo method. Our approach allows ST to be used with a continuous temperature and without pilot runs or normalizing constant approximations. STWNC is more computationally complex than PT but maintains the inherent capability of tempering algorithms to mix across challenging posterior regions. Furthermore STWNC provides easy access to Bayes Factors through thermodynamic integration and provides an additional heuristic diagnostic of assessing sampling in the neighbourhood of the optima along the full range of $\tau$.

Results from PT and PT-STWNC were comparable in all of our examples in terms of mixing of the target chain, estimation of points, and modal importance. However, STWNC is more computationally complex than PT because of the optimization stage. In the way the SIR example was set up in section 8, one would need to include 50 PT chains to produce approximately the same time per iteration as the 2 chain PT-STWNC. Results will depend on the complexity of the model in question and at the other extreme, the 2 chain PT and 2 chain PT-STWNC were approximately equivalent in speed in section 6 because the optima are analytically available. Optimization of the posterior distribution with respect to the parameters for a fixed $\tau$ may be initially computationally expensive in complex models. However, the optimization time decays with the number of MCMC iterations if previous optimal values of $\theta$ are used to initialize the optimizer in the next Markov Chain iteration. Further reduction in the cost as iterations progress could be obtained by replacing the optimizer with a function evaluation. Such a replacement would make the STWNC chain as fast as a PT chain but with much more capacity to overcome distant modes. The novelty of the STWNC is in the mathematical nature of the solution rather than its computational speed, but ongoing work shows considerable promise to accelerate the sampler. Additionally, having an implementable ST method allows additional work on exploiting the geometry of the problem, producing alternative pathways between the prior and posterior, and exploiting the now obtainable $\tau$ distribution as part of a modified Hamiltonian Monte Carlo method.

COMPUTER CODE

Examples in this paper are implemented in (R Core Team, 2017), using following packages: utils, graphics, parallel, deSolve (Soetaert et al, 2010), MASS (Venables and Ripley, 2002), gtools (Warnes et al, 2015), MCMCpack (Martin et al, 2011), mvtnorm (Genz and Bretz, 2009), truncnorm (Mersmann et al, 2018), optimx (Nash, 2014), coda (Plummer et al, 2006). Relevant R code and data is provided in the public github repository https://github.com/BiljanaJSJ/PT-STWNC.

10 Apendices

10.1 Analytical calculation of the marginal likelihood

In this section we provide the details on analytic calculation of the marginal likelihood in (16) for the bimodal model example given in the Section (6).

The posterior distribution of the unknown parameter $\mu$ is,

[TABLE]

where the variance was fixed at $\sigma^{2}=1$. The prior of $\mu$ was Gaussian: $P(\mu)\sim\text{N}(\lambda=0,\beta=1)$.

The Thermodynamic Integral is:

[TABLE]

where $a=\frac{n}{\sigma^{2}}+\frac{1}{\beta}$, $b=\frac{\sum_{i}y_{i}}{\sigma^{2}}+\frac{\lambda}{\beta}$, $b_{1}=\frac{-\sum_{i}y_{i}}{\sigma^{2}}+\frac{\lambda}{\beta}$ and $c=\frac{\sum_{i}y_{i}^{2}}{\sigma^{2}}+\frac{\lambda^{2}}{\beta}$.

Plug in the values of $\{\boldsymbol{Y},\lambda,\beta,\sigma^{2},n\}$ in the solution of the integral given by the equation (LABEL:eq:thInt) and take a log to obtain the analytical marginal likelihood reported in Table 2.

10.2 Implementation, bimodal model

The transition kernel of $\mu$ was updated using the optimal symmetric jumping kernel for Gaussian target distributions by (Gelman et al, 1996). The proposal distribution for of $\mu$ at the the $i$-th iteration is:

[TABLE]

where $\mathrm{Var}_{P(\mu\mid\sigma^{2},\textbf{Y},\tau)}(\mu)$ is the target variance of $\mu$ with respect to the target posterior distribution $P(\mu\mid\sigma^{2},\boldsymbol{Y},\tau)$, $[2.4/\sqrt{d}]^{2}$ is the optimal scale factor of the target variance found by (Gelman et al, 1996) with d being the dimension of the parameters updated in the MCMC step. The variance parameter $\sigma^{2}$ was sampled from a log normal proposal distribution. The transition step was tuned so that the acceptance rate is 44 $\%$. The inverse temperature parameter was updated by drawing independent samples from the standard uniform proposal distribution.

In order to evaluate the prior $P(\tau)$ in (11), $\mu$ and $\sigma^{2}$ were optimized using closed forms of the conditional posterior mean of $P(\mu\mid\boldsymbol{Y},\sigma^{2},\tau)$ and the conditional posterior mode of $P(\sigma^{2}\mid\boldsymbol{Y},\mu,\tau)$, respectively. In particular, $\mu$ and $\sigma^{2}$ were maximized in a conditional iterative manner by optimizing each of them conditional on the last optimized value of the other. Iterations were repeated until the optimized values of the both parameters stopped changing within a tolerance level of $10^{-3}$. Using explicit formulae of posterior means (or modes) avoids numerical issues that are usually associated with optimization routines.

Marginal likelihood estimation using thermodynamic integration via Parallel Tempering – bimodal model

Convergence of each of the PT chains used to obtain the TI-PT-NB and the TI-PT-B estimates was assessed using the Potential Scale Reduction Factor (PSRF) or $\hat{\text{R}}$ statistics by (Gelman and Rubin, 1992), which compares the in-chain and between-chain variances of the chains for each of the $20$ runs. The observed $\hat{R}<1.1$ in our runs indicates that the chains have converged.

10.3 Implementation PT-STWNC, SIR model

The proposal distribution of $\tau$ was a truncated standard normal at [0,1]. The proposal distributions of $\alpha$ and $\beta$ were log-normal. The proposal distribution of the $I(0)$ was Binomial, set up as follows: the proposed at the $i-th$ iteration is: $\text{I}(0)^{(i)}=\mbox{Binomial}(\text{N},\frac{\text{I}(0)^{(i-1)}}{\text{N}})$, where $\text{N}=261$ is the total population.

Optimization of $\boldsymbol{\theta}=\left(\alpha,\beta,\text{I}(0)\right)^{{}^{\prime}}$, which is needed to evaluate the prior of the inverse temperature $\tau$, was carried out by first optimizing the continuous $(\alpha,\beta)^{{}^{\prime}}$ conditional on fixed discrete values of $\text{I}(0)=\{1,..,8\}$ using the Nelder-Mead optimization routine. Then, out of the eight optimized values $\left(\alpha_{\max},\beta_{\max},I(0)\mid\text{I}(0)\ \in\{1,2,..,8\}\right)^{{}^{\prime}}$, the one that maximizes the posterior distribution $P(\alpha,\beta,\text{I}(0)\mid\boldsymbol{Y},\tau)$ was chosen as a maximum.