\section{Importance tempering}
\label{sec:it}

<<echo=false,results=hide>>=
library(tgp)
library(maptree)
#options(width=65)
seed <- 0; set.seed(seed)
@ 

It is well--known that MCMC inference in Bayesian treed methods
suffers from poor mixing.  For example, Chipman et
al.~\cite{chip:geor:mccu:1998,chip:geor:mccu:2002} recommend
periodically restarting the MCMC to avoid chains becoming stuck in
local modes of the posterior distribution (particularly in tree
space).  The treed GP models are or no exception, although it is worth
remarking that using flexible GP models at the leaves of the tree
typically results in shallower trees, and thus less pathalogical
mixing in tree space.  Version 1.x provided some crude tools to help
mitigate the effects of poor mixing in tree space.  For example, the
{\tt R} argument to the {\tt b*} functions facilitates the restarts
suggested by Chipman et al.

A modern Monte Carlo technique for dealing with poor mixing in Markov
chain methods is to employ {\em tempering} to flatten the
peaks and raise the troughs in the posterior distribution so that
movements between modes is more fluid.  One such method, called {\em
  simulated tempering} (ST) \cite{geyer:1995}, is essentially the MCMC
analogue of the popular simulated annealing algorithm for
optimization.  The ST algorithm helps obtain samples from a multimodal
density $\pi(\theta)$ where standard methods, such as
Metropolis--Hastings (MH) \cite{met:1953,hast:1970} and Gibbs Sampling
(GS) \cite{geman:1984}, fail.

As will be shown in our examples, ST can guard against becoming stuck
in local modes of the {\tt tgp} posterior by encouraging better mixing
{\em between modes} via in increase in the acceptance rate of tree
modification proposals, particularly {\em prunes}.  However, as we
will see, ST suffers from inefficiency because it discards the lions
share of the samples it collects.  The discarded samples can be
recycled if they are given appropriate importance sampling (IS)
\cite{liu:2001} weights.  These weights, if combined carefully, can be
used to construct meta-estimators of expectations under the {\tt tgp}
posterior that have much lower variance compared to ST alone.  This
combined application of ST and IT is dubbed {\em importance tempering}
\cite{gra:samw:king:2009}.

\subsection{Simulated Tempering and related methods}
\label{sec:st}

ST is an application of the MH algorithm on the product space of
parameters and inverse temperatures $k\in [0,1]$.  That is, ST uses
MH to sample from the joint chain $\pi(\theta,k) \propto \pi(\theta)^k
p(k)$.  The inverse temperature is allowed to take on a discrete set
of values $k \in \{k_1,\dots,k_m: k_1 = 1, \; k_i > k_{i+1} \geq 0\}$,
called the {\em temperature ladder}.  Typically, ST calls for sampling
$(\theta,k)^{(t+1)}$ by first updating $\theta^{(t+1)}$ conditional on
$k^{(t)}$ and (possibly) on $\theta^{(t)}$, using MH or GS.  Then, for
a proposed $k' \sim q(k^{(t)} \rightarrow k')$, usually giving equal
probability to the nearest inverse temperatures greater and less than
$k^{(t)}$, an acceptance ratio is calculated:
\[
A^{(t+1)} = \frac{\pi(\theta^{(t+1)})^{k'} p(k') q(k' \rightarrow
  k^{(t)})}{\pi(\theta^{(t+1)})^{k^{(t)}} p(k^{(t)})
  q(k^{(t)}\rightarrow k')}.
\]
Finally, $k^{(t+1)}$ is determined according to the MH accept/reject
rule: set $k^{(t+1)} = k'$ with probability $\alpha^{(t+1)} =
\min\{1,A^{(t+1)}\}$, or $k^{(t+1)} = k^{(t)}$ otherwise.  Standard
theory for MH and GS gives that samples from the marginals
$\pi_{k_i}$ can be obtained by collecting samples
$\theta^{(t)}$ where $k^{(t)} = k_i$.  Samples from $\pi(\theta)$ are
obtained when $k^{(t)} = 1$.

The success of ST depends crucially on the ability of the Markov chain
frequently to: (a) visit high temperatures (low $k$) where the
probability of escaping local modes is increased; (b) visit $k=1$ to
obtain samples from $\pi$.  The algorithm can be tuned by:
(i.)~adjusting the number and location of the rungs of the temperature
ladder; or (ii.)~setting the pseudo-prior $p(k)$ for the inverse
temperature.

Geyer \& Thompson \cite{geyer:1995} give ways
of adjusting the spacing of the rungs of the ladder so that the ST
algorithm achieves between--temperature acceptance rates of 20--40\%.
More recently, authors have preferred to rely on defaults, e.g.,
\begin{equation}
 \;\;\;\;\;
k_i = \left\{ \begin{array}{cl}
(1+\Delta_k)^{1-i} & \mbox{geometric spacing}\\
\{1+\Delta_k (i-1)\}^{-1} & \mbox{harmonic spacing}
\end{array} \right. \;\;\;\;\ i=1,\dots,m.
\label{eq:ladder}
\end{equation}
Motivation for such default spacings is outlined by Liu
\cite{liu:2001}.  Geometric spacing, or uniform spacing of
$\log(k_i)$, is also advocated by Neal \cite{neal:1996,neal:2001} to
encourage the Markov chain to rapidly traverse the breadth of the
temperature ladder.  Harmonic spacing is more often used by a related
method called Metropolis coupled Markov chain Monte Carlo (MC$^3$)
\cite{geyer:1991}.  Both defaults are implemented in the {\tt tgp}
package, through the provided {\tt default.itemps} function.  A new
``sigmoidal'' option is also implemented, as discussed below.  The
rate parameter $\Delta_k>0$ can be problem specific.  Rather than work
with $\Delta_k$ the {\tt default.itemps} function allows the ladder to
be specified via $m$ and the hottest temperature $k_m$, thus fixing
$\Delta_k$ implicitly.  I.e., for the geometric ladder $\Delta_k =
(k_m)^{1/(1-m)}-1$, and for the harmonic ladder $\Delta_k =
\frac{(k_m)^{-1}-1}{m-1}$.

A sigmoidal ladder can provide a higher concentration of temperatures
near $k=1$ without sacrificing the other nice properties of the
geometric and harmonic ladders.  It is specified by first situating
$m$ indices $j_i\in \mathbb{R}$ so that $k_1 = k(j_1) = 1$ and $k_m =
k(j_m) = k_{\mbox{\tiny m}}$ under
\[
k(j_i) = 1.01 - \frac{1}{1+e^{j_i}}.
\]
The remaining $j_i, i=2,\dots,(m-1)$ are spaced evenly between $j_1$
and $j_m$ to fill out the ladder $k_i = k(j_i), i=1,\dots,(m-1)$.

By way of comparison, consider generating the three different types of
ladder with identical minimum inverse temperature $k_{\mbox{\tiny m}}
= 0.1$, the default setting in {\tt tgp}.
<<>>=
geo <- default.itemps(type="geometric")
har <- default.itemps(type="harmonic")
sig <- default.itemps(type="sigmoidal")
@ 
The plots in Figure \ref{f:itemps} show the resulting 
inverse temperature ladders, and their logarithms.
\begin{figure}[ht!]
<<label=it-itemps,fig=TRUE,echo=TRUE,width=6,height=9,include=FALSE>>=
par(mfrow=c(2,1))
all <- cbind(geo$k, har$k, sig$k)
matplot(all, pch=21:23,
        main="inv-temp ladders", xlab="indx", ylab="itemp")
legend("topright", pch=21:23, 
       c("geometric","harmonic","sigmoidal"), col=1:3)
matplot(log(all), pch=21:23,
        main="log(inv-temp) ladders", xlab="indx", ylab="itemp")
@
<<echo=false,results=hide>>=
graphics.off()
@
\centering
\includegraphics[height=5.9in,width=4.5in,trim=0 20 0 20]{tgp2-it-itemps}
\caption{Three different inverse temperature ladders, each with
$m=40$ temperatures starting at $k_1=1$ and ending at $k_m=0.1$}
\label{f:itemps}
\end{figure}
Observe how, relative to the geometric ladder, the harmonic ladder has a
higher concentration of inverse temperatures near zero, whereas the
sigmoidal ladder has a higher concentration near one.  

Once a suitable ladder has been chosen, the {\tt tgp} package
implementation of ST follows the suggestions of Geyer \& Thompson
\cite{geyer:1995} in setting the pseudo--prior, starting from a
uniform $p_0$.  First, $p_0$ is adjusted by {\em stochastic
  approximation}: add $c_0/[m(t+n_0)]$ to $\log p_0(k)$ for each $k_i
\ne k^{(t)}$ and subtract $c_0/(t+n_0)$ from $\log p_0(k^{(t)})$ over
$t=1,\dots,B$ {\em burn--in} MCMC rounds sampling from the joint
posterior of $(\theta, k)$.  Then, $p_0$ is normalized to obtain
$p_1$.  Before subsequent runs, specified via an {\tt R >= 2}
argument, {\em occupation numbers} $o(k_i) = \sum_{t=1}^B 1_{\{k^{(t)}
  = k_i\}}$, are used update $p(k_i) \propto p_1(k_i)/o(k_i)$.  Note
that, in this setting, the {\tt R} argument is used to update the
pseudo--prior only, not to restart the Markov chain.  

\subsection{Importance sampling from tempered distributions}
\label{sec:temp}

ST provides us with $\{(\theta^{(t)},k^{(t)}): t = 1,\ldots,T\}$,
where $\theta^{(t)}$ is an observation from $\pi_{k^{(t)}}$.  It is
convenient to write $\mathcal{T}_i = \{t: k^{(t)} = k_i\}$ for the
index set of observations at the $i^{\mbox{\tiny th}}$ temperature,
and let $T_i = |\mathcal{T}_i|$.  Let the vector of observations at
the $i^{\mbox{\tiny th}}$ temperature collect in $\bm{\theta}_i =
(\theta_{i1},\dots,\theta_{iT_i})$, so that
$\{\theta_{ij}\}_{j=1}^{T_i}\sim \pi_{k_i}$.  Each vector
$\bm{\theta}_i$ can be used to construct an IS estimator of
$E_{\pi}\{h(\theta)\}$ by setting
\[
\hat{h}_i = \frac{\sum_{j=1}^{T_i} w_i(\theta_{ij}) h(\theta_{ij})}
{\sum_{j=1}^{T_i} w_i(\theta_{ij})} 
\equiv \frac{\sum_{j=1}^{T_i} w_{ij}h(\theta_{ij})}{W_i},
\]
say.  That is, rather than obtain one estimator from ST (at the cold
temperature), we can obtain $m$ estimators (one at each temperature)
via IS. The efficiency of each estimator, $i=1,\dots,m$ can be measured
through its variance, but unfortunately this can be difficult to
calculate in general.  As a result, the notion of {\em effective
  sample size} \cite{liu:2001} (ESS) plays an important role in the study of
IS estimators.  Denote the vector of IS weights at the $i^{\mbox{\tiny
    th}}$ temperature as $\mathbf{w}_i = \mathbf{w}_i(\bm{\theta}_i) =
(w_i(\theta_{i1}),\ldots,w_i(\theta_{iT_i}))$, where $w_i(\theta) =
\pi(\theta)/\pi_{k_i}(\theta)$.  The ESS of $\hat{h}_i$ is defined by
\begin{equation}
\mathrm{ESS}(\mb{w}_i) = \frac{T}{1 +
  \mathrm{cv^2}(\mathbf{w}_i)}, \label{eq:essw}
\end{equation}
where $\mathrm{cv}(\mathbf{w}_i)$ is the \emph{coefficient of variation}
of the weights (in the $i^{\mbox{\tiny th}}$ temperature), given by
\begin{align*}
\mathrm{cv^2}(\mathbf{w}_i) &= \frac{\sum_{t=1}^T(w(\theta^{(t)}) -
  \bar{w})^2}{(T-1) \bar{w}^2}, &\mbox{where} &&
\bar{w} &= T^{-1} \sum_{t=1}^T w(\theta^{(t)}).  
\end{align*}
In {\sf R}:
<<>>=
ESS <- function(w)
{
  mw <- mean(w)
  cv2 <- sum((w-mw)^2)/((length(w)-1)*mw^2)
  ess <- length(w)/(1+cv2)
  return(ess)
}
@ 
This should not be confused with the concept of \emph{effective
  sample size due to autocorrelation} \cite{kass:1998} (due to
serially correlated samples coming from a Markov chain as in MCMC) as
implemented by the {\tt effectiveSize} function in the {\tt coda}
package \cite{coda:R} for {\sf R}.

Before attempting to combine $m$ IS estimators it is fruitful
backtrack briefly to obtain some perspective on the topic of applying
IS with a {\em single} tempered proposal distribution.  Jennison
\cite{jennison:1993} put this idea forward more than a decade ago,
although the question of how to choose the best temperature was neither
posed or resolved.  It is clear that larger $k$ leads to lower
variance estimators (and larger ESS), but at the expense of poorer
mixing in the Markov chain. It can be shown that the optimal inverse
temperature $k^*$ for IS, in the sense of constructing a minimum
variance estimator, may be significantly lower than one
\cite{gra:samw:king:2009}.  However, the variance of such an estimator
will indeed become unbounded as $k\rightarrow 0$, just as
ESS~$\rightarrow 0$.  Needless to say, the choice of how to best pick
the best temperatures (for ST or IS) is still an open problem.  But in
the context of the family of tempered distributions used by ST for
mixing considerations, this means that the discarded samples obtained
when $k^{(t)} < 1$ may actually lead to more efficient estimators than
the ones saved from the cold distribution.  So ST is wastefull indeed.

However, when combining IS estimators from the multiple temperatures
used in ST, the deleterious effect of the high variance ones obtained
at high temperature must be mitigated.  The possible strategies
involved in developing such a meta-estimator comprise the {\em
  importance tempering} (IT) family of methods.  The idea is that
small ESS will indicate high variance IS estimators which should be
relegated to having only a small influence on the overall estimator.

\subsection{An optimal way to combine IS estimators}
\label{sec:lambdas}

It is natural to consider an overall meta-estimator of
$E_{\pi}\{h(\theta)\}$ defined by a convex combination:
\begin{align}
\label{eq:hhatlambda}
\hat{h}_{\lambda} &= \sum_{i=1}^m \lambda_i \hat{h}_i,&
\mbox{where} && 0 \leq \lambda_i \leq \sum_{i=1}^m \lambda_i = 1.
\end{align}
Unfortunately, if $\lambda_1,\dots,\lambda_m$ are not chosen
carefully, $\mbox{Var}(\hat{h}_\lambda)$, can be nearly as large as
the largest $\mbox{Var}(\hat{h}_i)$ \cite{owen:2000}, due to the
considerations alluded to in Section \ref{sec:temp}.  Notice that ST
is recovered as a special case when $\lambda_1=1$ and
$\lambda_2,\dots,\lambda_m = 0$.  It may be tempting to choose
$\lambda_i = W_i/W$, where $W = \sum_{i=1}^m W_i$.  The resulting
estimator is equivalent to
\begin{align}
\label{Eq:hath}
\hat{h} &= W^{-1} \sum_{t=1}^T w(\theta^{(t)},k^{(t)})h(\theta^{(t)}), 
& \mbox{where} &&
W = \sum_{t=1}^T w(\theta^{(t)},k^{(t)}),
\end{align}
and $w(\theta,k) = \pi(\theta)/\pi(\theta)^k = \pi(\theta)^{1-k}$.  It
can lead to a very poor estimator, even compared to ST, as will be
demonstrated empirically in the examples to follow shortly.

Observe that we can equivalently write
\begin{align}
\hat{h}_{\lambda} &= \sum_{i=1}^m \sum_{j=1}^{T_i}
w_{ij}^{\lambda}h(\theta_{ij}), 
&& \mbox{where} & w_{ij}^{\lambda} &= \lambda_iw_{ij}/W_i.  
\label{eq:wlambda}
\end{align}
Let $\mathbf{w}^{\lambda} =
(w_{11}^\lambda,\ldots,w_{1T_1}^\lambda,w_{21}^\lambda,\ldots,w_{2T_2}^\lambda,
\ldots,w_{m1}^\lambda,\ldots,w_{mT_m}^\lambda)$.  Attempting to choose
$\lambda_1,\dots,\lambda_m$ to minimize $\mbox{Var}(\hat{h}_\lambda)$
directly can be difficult.  Moreover, for the applications that we
have in mind, it is important that our estimator can be constructed
without knowledge of the normalizing constants of
$\pi_{k_1},\ldots,\pi_{k_m}$, and without evaluating the MH transition
kernels $\mathcal{K}_{\pi_{k_i}}(\cdot,\cdot)$.  It is for this reason
that methods like the \emph{balance heuristic} \cite{veach:1995}, MCV
\cite{owen:2000}, or population Monte Carlo (PMC)
\cite{douc:etal:2007} cannot be applied.  Instead, we seek maximize
the effective sample size of $\hat{h}_\lambda$ in
(\ref{eq:hhatlambda}), and look for an $O(T)$ operation to determine
the optimal $\lambda^*$.

%\begin{thm}
%\label{thm:lambdastar}
Among estimators of the form~(\ref{eq:hhatlambda}), it can be shown
\cite{gra:samw:king:2009} that $\mathrm{ESS}(\mathbf{w}^{\lambda})$ is
maximized by $\lambda = \lambda^*$, where, for $i=1,\ldots,m$,
\begin{align*}
\lambda_i^* &= \frac{\ell_i}{\sum_{i=1}^m \ell_i}, & \mbox{and} && \ell_i &=
  \frac{W_i^2}{\sum_{j=1}^{T_i} w_{ij}^2}.
\end{align*}
The efficiency of each IS estimator $\hat{h}_i$ can be measured
through $\mathrm{ESS}(\mathbf{w}_i)$.  Intuitively, we hope that with
a good choice of $\lambda$, the ESS (\ref{eq:essw}) of
$\hat{h}_{\lambda}$, would be close to the sum over $i$ of the
effective sample sizes each of $\hat{h}_i$.  This is indeed the case
for $\hat{h}_{\lambda^*}$, because it can be shown
\cite{gra:samw:king:2009} that
\[
\mathrm{ESS}(\mathbf{w}^{\lambda^*}) \geq \sum_{i=1}^m \mathrm{ESS}(\mathbf{w}_i) 
- \frac{1}{4} - \frac{1}{T}.
\]
In practice we have found that this bound is conservative and that in
fact $\mathrm{ESS}(\mathbf{w}^{\lambda^*}) \geq \sum_{i=1}^m
\mathrm{ESS}(\mathbf{w}_i)$, as will be shown empirically in the
examples that follow.  Thus our optimally--combined IS estimator has a
highly desirable and intuitive property in terms of its effective
sample size: that the whole is greater than the sum of its parts.

$\mathrm{ESS}(\mathbf{w}^{\lambda^*})$ depends on
$\mathrm{ESS}(\mathbf{w}_i)$ which in turn depend on the $k_i$.
Smaller $k_i$ will lead to better mixing in the Markov chain, but
lower $\mathrm{ESS}(\mathbf{w}_i)$.  Therefore, we can expect that the
geometric and sigmoidal ladders will fare better than the harmonic
ones, so long as the desired improvements in mixing are achieved.  In
the examples to follow, we shall see that the sigmoidal ladder does
indeed leader to higher $\mathrm{ESS}(\mathbf{w}^{\lambda^*})$.


\subsection{Examples}
\label{sec:examples}

Here the IT method is shown in action for {\tt tgp} models.  IT is
controlled in {\tt b*} functions via the {\tt itemps} argument: a {\tt
  data.frame} coinciding with the output of the {\tt default.itemps}
function.  The {\tt lambda} argument to {\tt default.itemps} can be
used to base posterior predictive inference the other IT heuristics:
ST and the na\"ive approach (\ref{Eq:hath}).  Whenever the argument
{\tt m = 1} is used with {\tt k.min != 1} the resulting estimator is
constructed via tempered importance sampling at the single inverse
temperature {\tt k.min}, in the style of Jennison~\cite{jennison:1993}
as outlined in Section \ref{sec:temp}.  The parameters $c_0$ and $n_0$
for stochastic approximation of the pseudo--prior can be specified as
a 2--vector {\tt c0n0} argument to {\tt default.itemps}.  In the
examples which follow we simply use the default configuration of the
IT method, adjusting only the minimum inverse temperature via the {\tt
  k.min} argument.

Before delving into more involved examples, we illustrate the stages
involved in a small run of importance tempering (IT) on the
exponential data from Section 3.3 of \cite{gramacy:2007}.  The data
can be obtained as:
<<>>=
exp2d.data<-exp2d.rand() 
X<-exp2d.data$X 
Z<-exp2d.data$Z 
@ 
Now, consider applying IT to the Bayesian treed LM with a small
geometric ladder.  A warning will be given if the default setting of
\verb!bprior="bflat"! is used, as this (numerically) improper prior
can lead to improper posterior inference at high temperatures.
<<>>=
its <- default.itemps(m=10)
exp.btlm <- btlm(X=X,Z=Z, bprior="b0", R=2, itemps=its, pred.n=FALSE,
                 BTE=c(1000,3000,2)) 
@ 
Notice how the MCMC inference procedure starts with
$B+T=\Sexpr{exp.btlm$BTE[1] + exp.btlm$BTE[2]}$ rounds of stochastic
approximation (initial adjustment of the pseudo--prior) in place of
typical (default) the $B=\Sexpr{exp.btlm$BTE[1]}$ burn--in
rounds. Then, the first round of sampling from the posterior
commences, over $T=\Sexpr{exp.btlm$BTE[2]-exp.btlm$BTE[1]}$ 
rounds, during which the
observation counts in each temperature are tallied.  The progress
meter shows the current temperature the chain is in, say {\tt
  k=0.629961}, after each of 1000 sampling rounds.  The first repeat
starts with a pseudo--prior that has been adjusted by the observation
counts, which continue to be accumulated throughout the entire
procedure (i.e., they are never reset).  
Any subsequent repeats begin after a similar (re-)adjustment.

Before finishing, the routine summarizes the sample size and effective
sample sizes in each rung of the temperature ladder.  The number of
samples is given by {\tt len}, and the ESS by {\tt ess}.  These
quantities can also be recovered via {\tt traces}, as shown later. The
ESS of the optimal combined IT sample is the last quantity printed.
This, along with the ESS and total numbers of samples in each
temperature, can also be obtained via the {\tt tgp}-class output
object.
<<>>=
exp.btlm$ess
@ 


\subsubsection{Motorcycle accident data}
\label{sec:moto}

Recall the motorcycle accident data of Section 3.4 of the first {\tt
  tgp} vignette \cite{gramacy:2007}.  Consider using IT to sample from
the posterior distribution of the treed GP LLM model using the
geometric temperature ladder.
<<>>=
library(MASS)
moto.it <- btgpllm(X=mcycle[,1], Z=mcycle[,2], BTE=c(2000,52000,10),
        bprior="b0", R=3, itemps=geo, trace=TRUE, pred.n=FALSE, verb=0)
@ 
Out of a total of $\Sexpr{moto.it$R*moto.it$BTE[2]/moto.it$BTE[3]}$
samples from the joint chain, the resulting (optimally combined) ESS was:
<<>>=
moto.it$ess$combined
@ 
Alternatively, $\mb{w}^{\lambda^*}$ can be extracted from the
traces, and used to make the ESS calculation directly.
<<>>=
p <- moto.it$trace$post
ESS(p$wlambda)
@ 
The unadjusted weights $\mb{w}$ are also available from {\tt
  trace}.  We can see that the na\"{i}ve choice of $\lambda_i =
W_i/W$, leading to the estimator in (\ref{Eq:hath}), has a clearly
inferior effective sample size.
<<>>=
ESS(p$w)
@ 
To see the benefit of IT over ST we can simply count the number of
samples obtained when $k^{(t)} = 1$.  This can be accomplished in
several ways: either via the traces or through the output object.
<<>>=
as.numeric(c(sum(p$itemp == 1), moto.it$ess$each[1,2:3]))
@ 
That is, (optimal) IT gives effectively
$\Sexpr{signif(moto.it$ess$combined/sum(p$itemp==1), 3)}$ times more samples.
The na\"{i}ve combination, leading to the estimator in (\ref{Eq:hath}),
yields an estimator with an effective sample size that is
$\Sexpr{round(100*ESS(p$w)/sum(p$itemp==1))}$\% of the number of
samples obtained under ST.

Now, we should like to compare to the MCMC samples obtained
under the same model, without IT.
<<>>=
moto.reg <- btgpllm(X=mcycle[,1], Z=mcycle[,2], BTE=c(2000,52000,10),
        R=3, bprior="b0", trace=TRUE, pred.n=FALSE, verb=0)
@ 
The easiest comparison to make is to look at the heights explored
under the three chains: the regular one, the chain of heights visited
at all temperatures (combined), and those obtained after applying IT
via re-weighting under the optimal combination $\lambda^*$.
<<>>=
L <- length(p$height)
hw <- suppressWarnings(sample(p$height, L, prob=p$wlambda, replace=TRUE))
b <- hist2bar(cbind(moto.reg$trace$post$height, p$height, hw))
@
\begin{figure}[ht!]
<<label=it-moto-height,fig=TRUE,echo=TRUE,width=11,height=7,include=FALSE>>=
barplot(b, beside=TRUE, col=1:3, xlab="tree height", ylab="counts", 
         main="tree heights encountered")
legend("topright", c("reg MCMC", "All Temps", "IT"), fill=1:3)
@
<<echo=false,results=hide>>=
graphics.off()
@
\centering
\includegraphics[trim=0 25 0 25]{tgp2-it-moto-height}
\caption{Barplots indicating the counts of the number of times the 
  Markov chains (for regular MCMC, combining all temperatures in the
  inverse temperature ladder, and those re-weighted via IT) were in
  trees of various heights for the motorcycle data.}
\label{f:moto:it:heights}
\end{figure}
Figure \ref{f:moto:it:heights} shows barplots indicating the count of
the number of times the Markov chains were in trees of various heights
after burn--in.  Notice how the tempered chain (denoted ``All Temps''
in the figure) frequently visits trees of height one, whereas the
non--tempered chain (denoted ``reg MCMC'') never does.  The result is
that the non--tempered chain underestimates the probability of height
two trees and produces a corresponding overestimate of height four
trees---which are clearly not supported by the data---even visiting
trees of height five.  The IT estimator appropriately down--weights
height one trees and provides correspondingly more realistic estimates
of the probability of height two and four trees.

Whenever introducing another parameter into the model, like the
inverse temperature $k$, it is important to check that the marginal
posterior chain for that parameter is mixing well.  For ST it is
crucial that the chain makes rapid excursions between the cold
temperature, the hottest temperatures, and visits each temperature
roughly the same number of times.
\begin{figure}[ht!]
<<label=it-moto-ktrace,fig=TRUE,echo=TRUE,width=11,height=6,include=FALSE>>=
plot(log(moto.it$trace$post$itemp), type="l", ylab="log(k)", xlab="samples",
     main="trace of log(k)")
@
<<echo=false,results=hide>>=
graphics.off()
@
\centering
\includegraphics[trim=0 25 0 25]{tgp2-it-moto-ktrace}
\caption{A trace of the MCMC samples from the marginal posterior
  distribution of the inverse temperature parameter, $k$, in 
  the motorcycle experiment}
  \label{f:ktrace}
\end{figure}

Figure \ref{f:ktrace} shows a trace of the posterior samples for $k$
in the motorcycle experiment.  Arguably, the mixing in $k$--space
leaves something to be desired.  Since it can be very difficult to
tune the pseudo--prior and MH proposal mechanism to get good mixing in
$k$--space, it is fortunate that the IT methodology does not rely on
the same mixing properties as ST does.  Since samples can be obtained
from the posterior distribution of the parameters of interest by
re-weighting samples obtained when $k < 1$ it is only important that
the chain frequently visit low temperatures to obtain good sampling,
and high temperatures to obtain good mixing.  The actual time spent in
specific temperatures, i.e., $k=1$ is less important.
%%ylim <- c(0, 1.25*max(c(b[,1], moto.it$itemps$counts)))
%, ylim=ylim)
\begin{figure}[ht!]
<<label=it-moto-khist,fig=TRUE,echo=TRUE,width=10,height=6,include=FALSE>>=
b <- itemps.barplot(moto.it, plot.it=FALSE)
barplot(t(cbind(moto.it$itemps$counts, b)), col=1:2,
        beside=TRUE, ylab="counts", xlab="itemps", 
        main="inv-temp observation counts")
legend("topleft", c("observation counts", "posterior samples"), fill=1:2)
@
<<echo=false,results=hide>>=
graphics.off()
@
\centering
\includegraphics[trim=0 25 0 25]{tgp2-it-moto-khist}
  \caption{Comparing (thinned) samples from the posterior distribution for
  the inverse temperature parameter, $k$, (posterior samples), to
  the observation counts used to update the pseudo--prior, in the
  motorcycle experiment}
  \label{f:khist}
\end{figure}
Figure \ref{f:khist} shows the histogram of the inverse temperatures
visited in the Markov chain for the motorcycle experiment.  Also
plotted is a histogram of the {\em observation counts} in each
temperature.  The two histograms should have similar shape but
different totals.  Observation counts are tallied during every MCMC
sample after burn--in, whereas the posterior samples of $k$ are
thinned (at a rate specified in {\tt BTE[3]}).  When the default {\tt
  trace=FALSE} argument is used only the observation counts will be
available in the {\tt tgp}--class object, and these can be used as a
surrogate for a trace of $k$.

The compromise IT approach obtained using the sigmoidal ladder can
yield an increase in ESS.
<<>>=
moto.it.sig <- btgpllm(X=mcycle[,1], Z=mcycle[,2], BTE=c(2000,52000,10),
                      R=3, bprior="b0", krige=FALSE, itemps=sig, verb=0)
@ 
Compare the resulting ESS to the one given for the geometric
ladder above.
<<>>=
moto.it.sig$ess$combined
@ 
\begin{figure}[ht!]
<<label=it-moto-pred,fig=TRUE,echo=TRUE,width=10,height=5,include=FALSE>>=
plot(moto.it.sig)
@
<<echo=false,results=hide>>=
graphics.off()
@
\centering
\includegraphics[trim=0 25 0 25]{tgp2-it-moto-pred}
\caption{Posterior predictive surface for the motorcycle data, with
  90\% quantile errorbars, obtained under IT with the sigmoidal ladder.}
  \label{f:moto:pred}
\end{figure}
Plots of the resulting predictive surface is shown in Figure
\ref{f:moto:pred} for comparison with those in Section 1.1 of the
first {\tt tgp} vignette \cite{gramacy:2007}.  In particular, observe
that the transition from the middle region to the right one is much
less stark in this tempered version than than in the original---which
very likely spent a disproportionate amount of time stuck in a 
posterior mode with trees of depth three or greater.

\subsubsection{Synthetic 2--d Exponential Data}
\label{sec:exp}

Recall the synthetic 2--d exponential data of Section 3.4 of the tgp
vignette \cite{gramacy:2007}, where the true response is given by
\[
z(\mb{x}) = x_1 \exp(-x_1^2 - x_2^2).
\]
Here, we will take $\mb{x} \in [-6,6]\times [-6,6]$ with a
$D$--optimal design
<<>>=
Xcand <- lhs(10000, rbind(c(-6,6),c(-6,6)))
X <- dopt.gp(400, X=NULL, Xcand)$XX
Z <- exp2d.Z(X)$Z
@ 

Consider a treed GP LLM model fit to this data using the standard
MCMC.
<<>>=
exp.reg <- btgpllm(X=X, Z=Z, BTE=c(2000,52000,10), bprior="b0", 
                   trace=TRUE, krige=FALSE, R=10, verb=0)
@ 
\begin{figure}[ht!]
<<label=it-exp-pred,fig=TRUE,echo=TRUE,width=12.5,height=7,include=FALSE>>=
plot(exp.reg)
@
<<echo=false,results=hide>>=
graphics.off()
@
\centering
\includegraphics[trim=0 25 0 25]{tgp2-it-exp-pred}
\caption{Posterior predictive surface for the 2--d exponential data:
  mean surface {\em (left)} and 90\% quantile difference {\em (right)}}
  \label{f:exp:pred}
\end{figure}
Figure \ref{f:exp:pred} shows the resulting posterior predictive
surface. The maximum {\em a' posteriori} (MAP) tree is drawn over the
error surface in the {\em right--hand} plot.  The height of this tree
can be obtained from the {\tt tgp}-class object.
<<>>=
h <- exp.reg$post$height[which.max(exp.reg$posts$lpost)]
h
@ 
It is easy to see that many fewer partitions are actually necessary
to separate the interesting, central, region from the surrounding flat
region.
\begin{figure}[ht!]
<<label=it-exp-mapt,fig=TRUE,echo=TRUE,width=11,height=7,include=FALSE>>=
tgp.trees(exp.reg, "map")
@
<<echo=false,results=hide>>=
graphics.off()
@
\centering
\includegraphics[trim=0 100 0 25]{tgp2-it-exp-mapt}
\caption{Diagrammatic depiction of the maximum {\em a' posteriori} (MAP)
  tree for the 2--d exponential data under standard MCMC sampling }
  \label{f:exp:mapt}
\end{figure}
Figure \ref{f:exp:mapt} shows a diagrammatic representation of the MAP
tree.  Given the apparent over--partitioning in this height \Sexpr{h}
tree it would be surprising to find much posterior support for trees
of greater height.  One might indeed suspect that there are trees with
fewer partitions which would have higher posterior probability, and
thus guess that the Markov chain for the trees plotted in these
figures possibly became stuck in a local mode of tree space while on
an excursion into deeper trees.

Now consider using IT.  It will be important in this case to have a
$k_{\mbox{\tiny m}}$ small enough to ensure that the tree occasionally
prunes back to the root.  We shall therefore use a smaller
$k_{\mbox{\tiny m}}$.  % with an extra 10 rungs.
Generally speaking, some pilot tuning may be necessary to choose an
appropriate $k_{\mbox{\tiny m}}$ and number of rungs $m$, although the
defaults should give adequate performance in most cases.
<<>>=
its <- default.itemps(k.min=0.02)
exp.it <- btgpllm(X=X, Z=Z, BTE=c(2000,52000,10), bprior="b0", 
               trace=TRUE, krige=FALSE, itemps=its, R=10, verb=0)
@ 
As expected, the tempered chain moves more rapidly throughout tree
space by accepting more tree proposals.  The acceptance rates of tree
operations can be accessed from the {\tt tgp}--class object.
<<>>=
exp.it$gpcs
exp.reg$gpcs
@ 
The increased rate of {\em prune} operations explains how the
tempered distributions helped the chain escape the local modes of deep
trees.

We can quickly compare the effective sample sizes of the three possible
estimators: ST, na\"{i}ve IT, and optimal IT.
<<>>=
p <- exp.it$trace$post
data.frame(ST=sum(p$itemp == 1), nIT=ESS(p$w), oIT=exp.it$ess$combined)
@ 
Due to the thinning in the Markov chain ({\tt BTE[3] = 10}) and the
traversal between $m=10$ temperatures in the ladder, we can be
reasonably certain that the \Sexpr{round(exp.it$ess$combined)} samples obtained
via IT from the total of
\Sexpr{round(exp.it$R*(exp.it$BTE[2]-exp.it$BTE[1])/exp.it$BTE[3])} 
samples obtained from the posterior are far less correlated than the 
ones obtained via standard MCMC.

As with the motorcycle data, we can compare the tree heights visited
by the two chains.
<<>>=
L <- length(p$height)
hw <- suppressWarnings(sample(p$height, L, prob=p$wlambda, replace=TRUE))
b <- hist2bar(cbind(exp.reg$trace$post$height, p$height, hw))
@
\begin{figure}[ht!]
<<label=it-exp-height,fig=TRUE,echo=TRUE,width=11,height=7,include=FALSE>>=
barplot(b, beside=TRUE, col=1:3, xlab="tree height", ylab="counts", 
         main="tree heights encountered")
legend("topright", c("reg MCMC", "All Temps", "IT"), fill=1:3)
@
<<echo=false,results=hide>>=
graphics.off()
@
\centering
\includegraphics[trim=0 25 0 25]{tgp2-it-exp-height}
\caption{Barplots indicating the counts of the number of times the 
  Markov chains (for regular MCMC, combining all temperatures in the
  inverse temperature ladder, and those re-weighted via IT) were in
  trees of various heights for the 2--d exponential data.}
\label{f:exp:it:heights}
\end{figure}
Figure \ref{f:exp:it:heights} shows a barplot of {\tt b}, which
illustrates that the tempered chain frequently visited shallow trees.
IT with the optimal weights shows that the standard MCMC chain
missed many trees of height three and four with considerable posterior
support.

\begin{figure}[ht!]
<<label=it-exp-trace-height,fig=TRUE,echo=TRUE,width=11,height=7,include=FALSE>>=
ylim <- range(p$height, exp.reg$trace$post$height)
plot(p$height, type="l", main="trace of tree heights", 
     xlab="t", ylab="height", ylim=ylim)
lines(exp.reg$trace$post$height, col=2)
legend("topright", c("tempered", "reg MCMC"), lty=c(1,1), col=1:2)
@
<<echo=false,results=hide>>=
graphics.off()
@
\centering
\includegraphics[trim=0 25 0 25]{tgp2-it-exp-trace-height}
\caption{Traces of the tree heights obtained under the two Markov
  chains (for regular MCMC, combining all temperatures in the inverse
  temperature ladder) on the 2--d exponential data.}
\label{f:exp:trace:height}
\end{figure}
To more directly compare the mixing in tree space between the ST and
tempered chains, consider the trace plots of the heights of the trees
explored by the chains shown in Figure \ref{f:exp:trace:height}.
Despite being restarted \Sexpr{exp.reg$R} times, the regular MCMC
chain (almost) never visits trees of height less than five after burn--in and
instead makes rather lengthy excursions into deeper trees, exploring
a local mode in the posterior.  In contrast, the tempered chain
frequently prunes back to the tree root, and consequently discovers
posterior modes in tree heights three and four.

\begin{figure}[ht!]
<<label=it-expit-pred,fig=TRUE,echo=TRUE,width=14,height=8,include=FALSE>>=
plot(exp.it)
@
\vspace{-0.7cm}
<<label=it-expit-trees,fig=TRUE,echo=TRUE,width=12,height=8,include=FALSE>>=
tgp.trees(exp.it, "map")
@
<<echo=false,results=hide>>=
graphics.off()
@
\centering
\includegraphics[trim=0 15 0 0]{tgp2-it-expit-pred}
\includegraphics[trim=0 100 0 0]{tgp2-it-expit-trees}
\caption{2--d exponential data fit with IT. {\em Top:} Posterior
  predictive mean surface for the 2d--exponential, with the MAP tree
  overlayed.  {\em Bottom:} diagrammatic representation of the MAP tree. }
  \label{f:exp-it:pred}
\end{figure}
To conclude,  a plot of the posterior predictive surface is given in
Figure \ref{f:exp-it:pred}, where the MAP tree is shown both graphically and
diagrammatically.
