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

\subsection{1-d Linear data}
\label{sec:ex:1dlinear}

Consider data sampled from a linear model.
\begin{equation} 
z_i = 1 + 2x_i + \epsilon_, \;\;\;\;\; \mbox{where} \;\;\;
\epsilon_i \stackrel{\mbox{\tiny iid}}{\sim} N(0,0.25^2) 
\label{eq:linear:sim}
\end{equation} 

The following {\sf R} code takes a sample $\{\mb{X}, \mb{Z}\}$ of size
$N=50$ from (\ref{eq:linear:sim}).  It also chooses $N'=99$ evenly spaced
predictive locations $\tilde{\mb{X}} = \mbox{\tt XX}$.
<<>>=
# 1-d linear data input and predictive data
X <- seq(0,1,length=50)  # inputs
XX <- seq(0,1,length=99) # predictive locations
Z <- 1 + 2*X + rnorm(length(X),sd=0.25) # responses
@ 

Using {\tt tgp} on this data with a Bayesian
hierarchical linear model goes as follows:
<<>>=
lin.blm <- blm(X=X, XX=XX, Z=Z)
@ 
\begin{figure}[ht!]
\centering
<<label=linear-blm,fig=TRUE,echo=TRUE,width=7,height=4.5,include=FALSE>>=
plot(lin.blm, main='Linear Model,', layout='surf')
abline(1,2,lty=3,col='blue')
@
<<echo=false,results=hide>>=
graphics.off()
@
\includegraphics[trim=0 25 0 25]{tgp-linear-blm}
%\vspace{-0.5cm}
\caption{Posterior predictive distribution using {\tt blm} on
  synthetic linear data: mean and 90\% credible interval.  The actual
  generating lines are shown as blue-dotted.}
\label{f:lin:blm}
\end{figure}

MCMC progress indicators are echoed every 1,000 rounds.  The linear
model is indicated by {\tt d=[0]}.  For {\tt btlm} the MCMC progress
indicators are boring, but we will see more interesting ones later.
In terminal versions, e.g. {\tt Unix}, the progress indicators can
give a sense of when the code will finish.  GUI versions of {\tt
  R}---{\tt Windows} or {\tt MacOS X}---can buffer {\tt stdout},
rendering this feature essentially useless as a real--time indicator of
progress.  Progress indicators can be turned off by providing the
argument {\tt verb=0}.  Further explanation on the verbosity of screen
output and interpretations is provided in Appendix \ref{sec:progress}.

The generic {\tt plot} method can be used to visualize the fitted
posterior predictive surface (with option {\tt layout = 'surf'}) 
in terms of means and credible intervals.
Figure \ref{f:lin:blm} shows how to do it, and what you get.  
The default option {\tt layout = 'both'} shows both a predictive
surface and error (or uncertainty) plot, side by side.
The error plot can be obtained alone via {\tt layout = 'as'}.
Examples of these layouts appear later.

If, say, you were unsure about the dubious ``linearness'' of this
data, you might try a GP LLM (using {\tt bgpllm}) and let a more
flexible model speak as to the linearity of the process.
<<>>=
lin.gpllm <- bgpllm(X=X, XX=XX, Z=Z)
@
\begin{figure}[ht!]
\centering
<<label=linear-gplm,fig=TRUE,echo=TRUE,width=7,height=4.5,include=FALSE>>=
plot(lin.gpllm, main='GP LLM,', layout='surf')
abline(1,2,lty=4,col='blue')
@
<<echo=false,results=hide>>=
graphics.off()
@
\includegraphics[trim=0 25 0 25]{tgp-linear-gplm}
%\vspace{-0.5cm}
\caption{Posterior predictive distribution using {\tt
    bgpllm} on synthetic linear data: mean and 90\% credible interval.
  The actual generating lines are shown as blue-dotted.}
\label{f:lin:gpllm}
\end{figure}
Whenever the progress indicators show {\tt d=[0]} the process is
under the LLM in that round, and the GP otherwise.  A plot of the
resulting surface is shown in Figure \ref{f:lin:gpllm} for comparison.
Since the data is linear, the resulting predictive surfaces should
look strikingly similar to one another.  On occasion, the GP LLM may
find some ``bendyness'' in the surface.  This happens rarely with samples
as large as $N=50$, but is quite a bit more common for $N<20$.

To see the proportion of time the Markov chain spent in the LLM
requires the gathering of traces (Appendix \ref{sec:traces}).  For
example
<<>>=
lin.gpllm.tr <- bgpllm(X=X, XX=0.5, Z=Z, pred.n=FALSE, trace=TRUE,
                       verb=0)
mla <- mean(lin.gpllm.tr$trace$linarea$la)
mla
@ 
shows that the average area under the LLM is \Sexpr{signif(mla,3)}.
Progress indicators are suppressed with \verb!verb=0!.  Alternatively,
the probability that input location {\tt xx} =
\Sexpr{lin.gpllm.tr$XX[1,]} is under the LLM is given by
<<>>=
1-mean(lin.gpllm.tr$trace$XX[[1]]$b1)
@ 
This is the same value as the area under the LLM since the process
is stationary (i.e., there is no treed partitioning).
