\subsection{Friedman data}
\label{sec:fried}

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

This Friedman data set is the first one of a suite that was used to
illustrate MARS (Multivariate Adaptive Regression Splines)
\cite{freid:1991}.  There are 10 covariates in the data ($\mb{x} =
\{x_1,x_2,\dots,x_{10}\}$).  The function that describes the
responses ($Z$), observed with standard Normal noise, has mean
\begin{equation}
E(Z|\mb{x}) = \mu = 10 \sin(\pi x_1 x_2) + 20(x_3 - 0.5)^2 + 10x_4 + 5 x_5,
\label{eq:f1}
\end{equation}
but depends only on $\{x_1,\dots,x_5\}$, thus combining nonlinear, linear,
and irrelevant effects.  Comparisons are made on this data to results
provided for several other models in recent literature.  Chipman et
al.~\cite{chip:geor:mccu:2002} used this data to compare
their treed LM algorithm to four other methods of varying
parameterization: linear regression, greedy tree, MARS, and neural
networks.  The statistic they use for comparison is root mean-square
error (RMSE)
\begin{align*}
\mbox{MSE} &= \textstyle \sum_{i=1}^n (\mu_i - \hat{z}_i)^2/n 
& \mbox{RMSE} &= \sqrt{\mbox{MSE}}
\end{align*}
where $\hat{z}_i$ is the model--predicted response for input
$\mb{x}_i$.  The $\mb{x}$'s are randomly distributed on the unit
interval.

Input data, responses, and predictive locations of size $N=200$ and
$N'=1000$, respectively, can be obtained by a function included in the
{\tt tgp} package.
<<>>=
f <- friedman.1.data(200)
ff <- friedman.1.data(1000)
X <- f[,1:10]; Z <- f$Y
XX <- ff[,1:10]
@ 
This example compares Bayesian treed LMs with Bayesian GP LLM (not
treed), following the RMSE experiments of Chipman et al.  It
helps to scale the responses so that they have a mean of zero and a
range of one.  First, fit the Bayesian treed LM, and obtain
the RMSE.
<<>>=
fr.btlm <- btlm(X=X, Z=Z, XX=XX, tree=c(0.95,2), pred.n=FALSE, verb=0)
fr.btlm.mse <- sqrt(mean((fr.btlm$ZZ.mean - ff$Ytrue)^2))
fr.btlm.mse
@ 
Next, fit the GP LLM, and obtain its RMSE.
<<>>=
fr.bgpllm <- bgpllm(X=X, Z=Z, XX=XX, pred.n=FALSE, verb=0)
fr.bgpllm.mse <- sqrt(mean((fr.bgpllm$ZZ.mean - ff$Ytrue)^2))
fr.bgpllm.mse
@ 
So, the GP LLM is \Sexpr{signif(fr.btlm.mse/fr.bgpllm.mse,4)} times
better than Bayesian treed LM on this data, in terms of RMSE (in
terms of MSE the GP LLM is
\Sexpr{signif(sqrt(fr.btlm.mse)/sqrt(fr.bgpllm.mse),4)} times better).

Parameter traces need to be gathered in order to judge the ability
of the GP LLM model to identify linear and irrelevant effects.
<<>>=
XX1 <- matrix(rep(0,10), nrow=1)
fr.bgpllm.tr <- bgpllm(X=X, Z=Z, XX=XX1, pred.n=FALSE, trace=TRUE, 
                       m0r1=FALSE, verb=0)
@ 
Here, \verb!m0r1 = FALSE! has been specified instead so that the
$\bm{\beta}$ estimates provided below will be on the original 
scale.\footnote{The default setting of {\tt m0r1 = TRUE} causes the 
{\tt Z}--values to undergo pre-processing so that they have a mean of 
zero and a range of one.  The default prior specification has been
tuned so as to work well this range.} 
A summary of the parameter traces show that the
Markov chain had the following (average) configuration for the booleans.
<<>>=
trace <- fr.bgpllm.tr$trace$XX[[1]]
apply(trace[,27:36], 2, mean)
@ 
Therefore the GP LLM model correctly identified that only the first three
input variables interact only linearly with
the response.  This agrees with dimension--wise estimate of the total
area of the input domain under the LLM (out of a total of 10 input variables).
<<>>=
mean(fr.bgpllm.tr$trace$linarea$ba)
@ 
A similar summary of the parameter traces for $\bm{\beta}$ shows
that the GP LLM correctly identified the linear regression
coefficients associated with the fourth and fifth input covariates
(from (\ref{eq:f1}))
<<>>=
summary(trace[,9:10])
@ 
and that the rest are much closer to zero.
<<>>=
apply(trace[,11:15], 2, mean)
@ 
