\subsection{Parameter traces}
\label{sec:traces}

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

Traces of (almost) all parameters to the {\tt tgp} model can be
collected by supplying {\tt trace=TRUE} to the {\tt b*} functions.
In the current version, traces for the linear prior correlation matrix
($\mb{W}$) are not provided.  I shall illustrate the gathering and
analyzing of traces through example.  But first, a few notes and
cautions.

Models which involve treed partitioning may have more than one base
model (GP or LM).  The process governing a particular input
$\mb{x}$ depends on the coordinates of $\mb{x}$.  As such, {\tt tgp}
records region--specific traces of parameters to GP (and linear)
models at the locations enumerated in the {\tt XX} argument.  Even
traces of single--parameter Markov chains can require hefty amounts of
storage, so recording traces at each of the {\tt XX} locations can be
an enormous memory hog.  A related warning will be given if the
product of $|${\tt XX}$|$, \verb!(BTE[2]-BTE[1])/BTE[3]! and {\sf R}
is beyond a threshold.  The easiest way to keep the storage
requirements for traces down is to control the size of {\tt XX} and
the thinning level {\tt BTE[3]}.  Finally, traces for most of the
parameters are stored in output files.  The contents of the trace
files are read into {\sf R} and stored as {\tt data.frame} objects,
and the files are removed. The existence of partially written trace
files in the current working directory (CWD)---while the {\tt C} code
is executing---means that not more than one {\tt tgp} run (with
\verb!trace = TRUE!) should be active in the CWD at one time.

Consider again the exponential data.  For illustrative purposes I
chose {\tt XX} locations (where traces are gathered) to be (1) in the
interior of the interesting region, (2) on/near the plausible
intersection of partition boundaries, and (3) in the interior of the
flat region.  The hierarchical prior \verb!bprior = "b0"! is used to
leverage a (prior) belief the most of the input domain is
uninteresting.
<<>>=
exp2d.data <- exp2d.rand(n2=150, lh=0, dopt=10)
X <- exp2d.data$X
Z <- exp2d.data$Z
XX <- rbind(c(0,0),c(2,2),c(4,4))
@ 
We now fit a treed GP LLM and gather traces, and also gather
EI and ALC statistics for the purposes of illustration.  Prediction
at the input locations {\tt X} is turned off to be thrifty.
<<>>=
out <- btgpllm(X=X, Z=Z, XX=XX, corr="exp", bprior="b0", 
               pred.n=FALSE, Ds2x=TRUE, R=10, 
               trace=TRUE, verb=0)
@ 
\begin{figure}[hp]
\centering
<<>>=
out$trace
@
\caption{Listing the contents of {\tt "tgptraces"}--class objects.}
\label{f:tgptraces}
\end{figure} 
Figure \ref{f:tgptraces} shows a dump of \verb!out$trace! which is a
\verb!"tgptraces"!--class object.  It depicts the full set of parameter
traces broken down into the elements of a \verb!list!: \verb!$XX!
with GP/LLM parameter traces for each {\tt XX} location (the
parameters are listed); \verb!$hier! with traces for
(non--input--dependent) hierarchical parameters (listed);
\verb!$linarea!  recording proportions of the input space under the
LLM; \verb!$parts!  with the boundaries of all partitions visited;
\verb!$post! containing (log) posterior probabilities; \verb!preds!
containing traces of samples from the posterior predictive
distribution and adaptive sampling statistics.

\begin{figure}[ht!]
\centering
<<label=traces-XXd,fig=TRUE,echo=TRUE,include=FALSE,width=8,height=5>>=
trXX <- out$trace$XX; ltrXX <- length(trXX)
y <- trXX[[1]]$d
for(i in 2:ltrXX) y <- c(y, trXX[[i]]$d)
plot(log(trXX[[1]]$d), type="l", ylim=range(log(y)), ylab="log(d)",
     main="range (d) parameter traces")
names <- "XX[1,]"
for(i in 2:ltrXX) {
  lines(log(trXX[[i]]$d), col=i, lty=i)
  names <- c(names, paste("XX[", i, ",]", sep=""))
}
legend("bottomleft", names, col=1:ltrXX, lty=1:ltrXX)
@
<<echo=false,results=hide>>=
graphics.off()
@
\includegraphics[trim=55 25 65 20]{tgp-traces-XXd}
\caption{Traces of the (log of the) first range
  parameter for each of the three {\tt XX} locations}
\label{f:XXd}
\end{figure} 

Plots of traces are useful for assessing the mixing of the Markov
chain.  For example, Figure \ref{f:XXd} plots traces of the range
parameter ($d$) %in the first input dimension ($d_1$)
for each of the
\Sexpr{length(out$trace$XX)} predictive locations {\tt XX}.  It is
easy to see which of the locations is in the same partition with
others, and which have smaller range parameters than others.

The mean area under the LLM can be calculated as
<<>>=
linarea <- mean(out$trace$linarea$la)
linarea
@ 
\begin{figure}[ht!]
\centering
<<label=traces-la,fig=TRUE,echo=TRUE,include=FALSE,width=8,height=6>>=
hist(out$trace$linarea$la)
@
<<echo=false,results=hide>>=
graphics.off()
@
\includegraphics[trim=0 0 0 20]{tgp-traces-la}
\vspace{-0.5cm}
\caption{Histogram of proportions of the area of the input 
  domain under the LLM}
\label{f:la}
\end{figure} 
This means that the expected proportion of the input domain under the
full LLM is \Sexpr{signif(linarea[1], 3)}.  Figure \ref{f:la} shows a
histogram of areas under the LLM.  The clumps near 0, 0.25, 0.5, and
0.75 can be thought of as representing quadrants (none, one, two, and
tree) under the LLM.  Similarly, we can calculate the probability that
each of the {\tt XX} locations is governed by the
LLM. % (in total, and by dimension)
<<>>=
m <- matrix(0, nrow=length(trXX), ncol=3)#ncol=5)
for(i in 1:length(trXX))
  m[i,] <- as.double(c(out$XX[i,], mean(trXX[[i]]$b)))
m <- data.frame(cbind(m, 1-m[,3]))
names(m)=c("XX1","XX2","b","pllm")
m
@ 
The final column above represents the probability that the
corresponding {\tt XX} location is under the LLM (which is equal to
{\tt 1-b}).

\begin{figure}[ht!]
\centering
<<label=traces-alc,fig=TRUE,echo=TRUE,include=FALSE,width=8,height=5>>=
trALC <- out$trace$preds$Ds2x
y <- trALC[,1]
for(i in 2:ncol(trALC)) y <- c(y, trALC[,i])
plot(log(trALC[,1]), type="l", ylim=range(log(y)), ylab="Ds2x",
     main="ALC: samples from Ds2x")
names <- "XX[1,]"
for(i in 2:ncol(trALC)) {
  lines(log(trALC[,i]), col=i, lty=i)
  names <- c(names, paste("XX[", i, ",]", sep=""))
}
legend("bottomright", names, col=1:ltrXX, lty=1:ltrXX)
@
<<echo=false,results=hide>>=
graphics.off()
@
\includegraphics[trim=55 25 65 20]{tgp-traces-alc}
\caption{Traces of the (log of the) samples
         for the ALC statistic  $\Delta \sigma^2(\tilde{\mb{x}})$
         at for each of the three {\tt XX} locations}
\label{f:preds}
\end{figure} 
Traces of posterior predictive and adaptive sampling statistics are
contained in the \verb!$preds! field.  For example, Figure
\ref{f:preds} shows samples of the ALC statistic $\Delta
\sigma^2(\tilde{\mb{x}})$.  We can see from the trace that statistic
is generally lowest for {\tt XX[3,]} which is in the uninteresting
region, and that there is some competition between {\tt XX[2,]} which
lies on the boundary between the regions, and {\tt XX[1,]} which is in
the interior of the interesting region.  Similar plots can be made for
the other adaptive sampling statistics (i.e., ALM \& EI).
