\subsection{Motorcycle Accident Data}
\label{sec:moto}

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

%\iffalse
The Motorcycle Accident Dataset \cite{silv:1985} is a classic
nonstationary data set used in recent literature
\cite{rasm:ghah:nips:2002} to demonstrate the success of nonstationary
models.  The data consists of measurements of the acceleration of the
head of a motorcycle rider as a function of time in the first moments
after an impact.  In addition to being nonstationary, the data has
input--dependent noise (heteroskedasticity) which makes it useful for
illustrating how the treed GP model handles this nuance.  There are at
least two---perhaps three---three regions where the response exhibits
different behavior both in terms of the correlation structure and
noise level.

The data is
%\else 
%In this section we return to the motivating Motorcycle Accident
%Dataset~\cite{silv:1985}, which is  
%\fi 
included as part of the {\tt MASS} library in
{\sf R}.
<<>>= 
library(MASS)
X <- data.frame(times=mcycle[,1])
Z <- data.frame(accel=mcycle[,2])
@ 
Figure \ref{f:moto:bgp} shows how a stationary GP is able to capture
the nonlinearity in the response but fails to capture the input
dependent noise and increased smoothness (perhaps linearity) in
parts of the input space.
<<echo=TRUE,results=hide>>=
moto.bgp <- bgp(X=X, Z=Z, verb=0)
@ 
Progress indicators are suppressed.
\begin{figure}[ht!]
\centering
<<label=moto-bgp,fig=TRUE,echo=TRUE,width=7,height=5,include=FALSE>>=
plot(moto.bgp, main='GP,', layout='surf')
@
<<echo=false,results=hide>>=
graphics.off()
@
\includegraphics[trim=0 25 0 25]{tgp-moto-bgp}
%\vspace{-0.5cm}
\caption{Posterior predictive distribution using {\tt
    bgp} on the motorcycle accident data: mean and 90\% credible
  interval}
\label{f:moto:bgp}
\end{figure}

A Bayesian Linear CART model is able to capture the input dependent
noise but fails to capture the waviness of the
``whiplash''---center--- segment of the response.
<<echo=TRUE,results=hide>>=
moto.btlm <- btlm(X=X, Z=Z, verb=0)
@ 
Figure \ref{f:moto:btlm} shows the resulting piecewise linear
predictive surface and MAP partition ($\hat{\mathcal{T}}$).
\begin{figure}[ht!]
\centering
<<label=moto-btlm,fig=TRUE,echo=TRUE,width=7,height=5,include=FALSE>>=
plot(moto.btlm, main='Bayesian CART,', layout='surf')
@
<<echo=false,results=hide>>=
graphics.off()
@
\includegraphics[trim=0 25 0 25]{tgp-moto-btlm}
%\vspace{-0.5cm}
\caption{Posterior predictive distribution using {\tt
    btlm} on the motorcycle accident data: mean and 90\% credible
  interval}
\label{f:moto:btlm}
\end{figure}

A treed GP model seems appropriate because it can model input
dependent smoothness {\em and} noise.  A treed GP LLM is probably
most appropriate since the left-hand part of the input space is likely
linear.  One might further hypothesize that the right--hand region is
also linear, perhaps with the same mean as the left--hand region, only
with higher noise.  The {\tt b*} functions can force
an i.i.d.~hierarchical linear model by setting \verb!bprior="b0"!.
<<echo=TRUE,results=hide>>=
moto.btgpllm <- btgpllm(X=X, Z=Z, bprior="b0", verb=0)
moto.btgpllm.p <- predict(moto.btgpllm) ## using MAP
@ 
The {\tt predict.tgp} function obtains posterior predictive estimates
from the MAP parameterization (a.k.a., {\em kriging}).
\begin{figure}[ht!]
\centering
<<label=moto-btgp,fig=TRUE,echo=TRUE,width=8,height=4,include=FALSE>>=
par(mfrow=c(1,2))
plot(moto.btgpllm, main='treed GP LLM,', layout='surf')
plot(moto.btgpllm.p, center='km', layout='surf')
@
<<echo=false,results=hide>>=
graphics.off()
@

\includegraphics[trim=50 25 50 20]{tgp-moto-btgp}

<<label=moto-btgpq,fig=TRUE,echo=TRUE,width=8,height=4,include=FALSE>>=
par(mfrow=c(1,2))
plot(moto.btgpllm, main='treed GP LLM,', layout='as')
plot(moto.btgpllm.p, as='ks2', layout='as')
@
<<echo=false,results=hide>>=
graphics.off()
@

\includegraphics[trim=50 25 50 20]{tgp-moto-btgpq}

%\vspace{-0.5cm}
\caption{{\em Top}: Posterior predictive distribution
  using treed GP LLM on the motorcycle accident data. The {\em
    left}--hand panes how mean and 90\% credible interval; {\em
    bottom}: Quantile-norm (90\%-5\%) showing input-dependent noise.
    The {\em right}--hand panes show similar {\em kriging} surfaces
  for the MAP parameterization.}
\label{f:moto:tgp}
\end{figure}
The resulting posterior predictive surface is shown in the {\em
  top--left} of Figure \ref{f:moto:tgp}.  The {\em bottom--left} of
the figure shows the norm (difference) in predictive quantiles,
clearly illustrating the treed GP's ability to capture input-specific
noise in the posterior predictive distribution.  The {\em right}--hand
side of the figure shows the MAP surfaces obtained from the output of
the {\tt predict.tgp} function.

The {\tt tgp}--default \verb!bprior="bflat"! implies an improper prior
on the regression coefficients $\bm{\beta}$.  It essentially forces
$\mb{W}=\mb{\infty}$, thus eliminating the need to specify priors on
$\bm{\beta}_0$ and $\mb{W}^{-1}$ in (\ref{eq:model}).  This was chosen
as the default because it works well in many examples, and leads to a
simpler overall model and a faster implementation.  However, the
Motorcycle data is an exception. Moreover, when the response data is
very noisy (i.e., low signal--to--noise ratio), {\tt tgp} can be
expected to partition heavily under the \verb!bprior="bflat"! prior.
In such cases, one of the other proper priors like the full
hierarchical \verb!bprior="b0"!  or \verb!bprior="bmzt"! might be preferred.

An anonymous reviewer pointed out a shortcoming of the treed GP model
on this data.  The sharp spike in predictive variance near the first
regime shift suggests that the symmetric Gaussian noise model may be
inappropriate.  A log Gaussian process might offer an improvement, at
least locally.  Running the treed GP MCMC for longer will eventually
result in the finding of a partition near time=17, just after the
first regime change.  The variance is still poorly modeled in this
region.  Since it is isolated by the tree it could potentially be fit
with a different noise model.
