\documentclass[11pt]{article}

\evensidemargin 0cm
\oddsidemargin 0cm

\topmargin -.94cm
\textwidth 16.0cm
\textheight 24.6cm

\usepackage[dvips]{epsfig}
%\usepackage{nips}
\usepackage{wrapfig}
\usepackage{subfigure}
\usepackage{graphics}
\usepackage{url}
\usepackage{natbib}
\usepackage{algorithm}
\usepackage{algorithmic}
\usepackage{rotate}
\usepackage{latexsym}

\usepackage{bbold}
\newcommand{\real}{\mathbb R}
\newcommand{\iden}{\mathbb 1}

\input xy
\xyoption{all}
\xyoption{dvips}
\xyoption{color}
%\input{../xyps-col.tex}

\newcommand{\scalefac}{1}
\newcommand{\comment}[1]{}

%\newcommand{\proof}[1]{\vspace*{5mm} \noindent {{\small \bf Proof.} \small #1}%
%\vspace*{3mm}}
%\newcommand{\example}[2]{\vspace*{3mm} \noindent {\bf Example #1.} #2%
%\vspace*{1mm}}

\newcommand{\onder}[2]{#1_{\mbox{\rm \scriptsize #2}}}
\newcommand{\boven}[2]{#1^{\mbox{\rm \scriptsize #2}}}
\newcommand{\Xmap}{\boven{X}{map}}
\newcommand{\mycite}[1]{(\cite{#1})}
\newcommand{\myciteyear}[1]{(\citeyear{#1})}
\newcommand{\qold}{\boven{q}{old}}
\newcommand{\Pold}{\boven{P}{old}}
\newcommand{\qnew}{\boven{q}{new}}
\newcommand{\Pnew}{\boven{P}{new}}
\newcommand{\mnew}{\boven{\mu}{new}}
\newcommand{\bgamma}{\bar{\gamma}}
\newcommand{\blambda}{\bar{\lambda}}
\newcommand{\hlambda}{\hat{\lambda}}
\newcommand{\hbeta}{\hat{\beta}}
\newcommand{\hy}{\hat{y}}
\newcommand{\he}{\hat{e}}
\newcommand{\hsigma}{\hat{\sigma}}
\newcommand{\halpha}{\hat{\alpha}}
\newcommand{\blnew}{\boven{\bar{\lambda}}{new}}
\newcommand{\tlambda}{\tilde{\lambda}}
\newcommand{\bmu}{\bar{\mu}}
\newcommand{\tmu}{\tilde{\mu}}
\newcommand{\tm}{\tilde{m}}
\newcommand{\ttheta}{\tilde{\theta}}
\newcommand{\hm}{\hat{m}}
\newcommand{\htheta}{\hat{\theta}}
\newcommand{\KL}{\mbox{KL}}
\newcommand{\ev}[1]{\left<#1\right>}
\newcommand{\collapse}{\mathop{\rm Collapse}\limits}
\newcommand{\argmax}{\mathop{\rm argmax}\limits}
\newcommand{\argmin}{\mathop{\rm argmin}\limits}
\newcommand{\cmin}{\mathop{\rm min'}\limits}
\newcommand{\cmax}{\mathop{\rm max'}\limits}
\newcommand{\ninp}{\onder{n}{inp}}
\newcommand{\Tr}{\mbox{Tr~}}
\newcommand{\vc}[1]{{\bf #1}}
\newcommand{\vb}[1]{\mbox{\boldmath $#1$}}
\newcommand{\vs}[1]{\mbox{\boldmath \scriptsize $#1$}}
\newcommand{\sss}[1]{\underline{#1}}
\newcommand{\stheta}{\sss{\theta}}
\newcommand{\vx}{\vc{x}}
\newcommand{\sx}{\sss{x}}
\newcommand{\se}{\sss{e}}
\newcommand{\sy}{\sss{y}}
\newcommand{\sz}{\sss{z}}
\newcommand{\spi}{\sss{\pi}}
\newcommand{\svx}{\sss{\vx}}
\newcommand{\vX}{\vc{X}}
\newcommand{\sX}{\sss{X}}
\newcommand{\bX}{\bar{X}}
\newcommand{\by}{\bar{y}}
\newcommand{\bx}{\bar{x}}
\newcommand{\vy}{\vc{y}}
\newcommand{\vt}{\vc{t}}
\newcommand{\vY}{\vc{Y}}
\newcommand{\vz}{\vc{z}}
\newcommand{\vf}{\vc{f}}
\newcommand{\vg}{\vc{g}}
\newcommand{\vm}{\vc{m}}
\newcommand{\vepsilon}{\vc{\epsilon}}
\newcommand{\vmmin}{\vm^{-}}
\newcommand{\vmnul}{\vm^{0}}
\newcommand{\vmplus}{\vm^{+}}
\newcommand{\vM}{\vc{M}}
\newcommand{\hvm}{\hat{\vc{m}}}
\newcommand{\vtheta}{\mbox{\boldmath $\theta$}}
\newcommand{\vT}{\mbox{\boldmath $\Theta$}}
\newcommand{\vxi}{\mbox{\boldmath $\xi$}}
\newcommand{\eye}{{\mathbb 1}}
\newcommand{\zero}{{\mathbb 0}}
\newcommand{\tP}{\tilde{P}}
\newcommand{\hP}{\hat{P}}
\newcommand{\hp}{p}
%\newcommand{\hp}{\hat{p}}
\newcommand{\pold}{\boven{\hp}{old}}
\newcommand{\hV}{\hat{V}}
\newcommand{\ee}{\mbox{e}}
\newcommand{\evnul}[1]{\ev{#1}_{q_t}}
\newcommand{\evbac}[1]{\ev{#1}_{\hp_t}}
\newcommand{\evfor}[1]{\ev{#1}_{\hp_{t+1}}}
\newcommand{\valpha}{\vb{\alpha}}
\newcommand{\vmu}{\vb{\mu}}
\newcommand{\hmu}{\hat{\mu}}
\newcommand{\vbeta}{\vb{\beta}}
\newcommand{\hvbeta}{\hat{\vbeta}}
\newcommand{\hvtheta}{\hat{\vtheta}}
\newcommand{\vgamma}{\vb{\gamma}}
\newcommand{\vGamma}{\vb{\Gamma}}
\newcommand{\vdelta}{\vb{\delta}}
\newcommand{\vDelta}{\vb{\Delta}}
\newcommand{\vsalpha}{\vs{\alpha}}
\newcommand{\vsbeta}{\vs{\beta}}
\newcommand{\vsgamma}{\vs{\gamma}}
\newcommand{\vsdelta}{\vs{\delta}}
\newcommand{\talpha}{\tilde{\valpha}}
\newcommand{\tbeta}{\tilde{\vbeta}}
\newcommand{\tgamma}{\tilde{\vgamma}}
\newcommand{\tdelta}{\tilde{\vdelta}}
\newcommand{\dalpha}{\partial \alpha}
\newcommand{\dbeta}{\partial \beta}
\newcommand{\dgamma}{\partial \gamma}
\newcommand{\ddelta}{\partial \delta}

\newcommand{\evv}[1]{\ev{\ev{#1}}}
\newcommand{\evvnul}[1]{\evv{#1}_{q_t}}
\newcommand{\evvbac}[1]{\evv{#1}_{\hp_t}}
\newcommand{\evvfor}[1]{\evv{#1}_{\hp_{t+1}}}

\newcommand{\gnew}{\boven{\gamma}{new}}
\newcommand{\lnew}{\boven{\lambda}{new}}
\newcommand{\gold}{\boven{\gamma}{old}}
\newcommand{\lold}{\boven{\lambda}{old}}
\newcommand{\dnew}{\boven{\vdelta}{new}}
\newcommand{\dold}{\boven{\vdelta}{old}}
\newcommand{\anew}{\boven{\valpha}{new}}
\newcommand{\aold}{\boven{\valpha}{old}}
\newcommand{\bnew}{\boven{\vbeta}{new}}
\newcommand{\bold}{\boven{\vbeta}{old}}

\newcommand{\Var}{\mbox{Var}}
\newcommand{\RSS}{\mbox{RSS}}
\newcommand{\lik}{\mbox{lik}}
\newcommand{\TrainError}{\mbox{TE}}
\newcommand{\PredError}{\mbox{PE}}
\newcommand{\ValError}{\mbox{VE}}

\newenvironment{proof}{\begin{quotation} {\it Proof.\ }}
			{$\Box$ \end{quotation}}
\newtheorem{proposition}{Proposition}


\begin{document}

\section*{Exercise 5: Linear Least Squares}

As part of a nuclear safeguards program, the contents of a tank are
routinely measured. The determination of volume is made indirectly by
measuring the difference in pressure at the top and at the bottom of
the tank. The tank is cylindrical in shape, but its internal geometry
is complicated by various pipes and agitator paddles. Without these
complications, pressure and volume should have a linear relationship:
\[
p = \rho g h \: ,
\]
with $p$ the pressure, $\rho$ the liquid density, $g$ the gravitational
constant, and $h$ the height. So if the tank is purely cylindrical, height
(and thus pressure) is proportional to volume.

To calibrate pressure with respect to volume, known quantities ($x$)
of liquid are placed in the tank and pressure readings ($y$) are taken.
The data are shown in Table~\ref{tab1} and can be downloaded from
{\verb#http://www.mbfys.kun.nl/~tom/files/nuclear.dat#}.

\begin{table}[h]
\begin{center}
\begin{tabular}{cr}
\hline \\[-2mm]
Volume & \multicolumn{1}{c}{Pressure} \\ 
(kiloliters) & \multicolumn{1}{c}{(pascal)} \\ \\[-2mm] \hline \\[-1mm]
0.18941 &   215.25 \\
0.18949 &   218.28 \\
0.37880 &   632.62 \\
0.37884 &   627.14 \\
0.37884 &   628.71 \\
0.56840 &  1034.05 \\
0.56843 &  1033.34 \\
0.75755 &  1469.11 \\
0.75767 &  1474.11 \\
0.75769 &  1475.24 \\
0.94740 &  1924.52 \\
0.94746 &  1921.85 \\
1.13642 &  2372.04 \\
1.13646 &  2374.43 \\
1.13665 &  2377.76 \\
1.32640 &  2819.47 \\
1.32640 &  2815.70 \\
1.51523 &  3263.50 \\
1.51528 &  3261.94 \\
1.51561 &  3268.54 \\
1.70525 &  3711.64
\end{tabular}
\caption{Volume versus pressure.}
\label{tab1}
\end{center}
\end{table}

\begin{description}
\item[a.] Plot pressure versus volume. Does the relationship appear linear?
\end{description}

\begin{description}
\item[b.] Calculate the linear regression of pressure on volume, and plot the
residuals versus volume. What does the residual plot show?
\end{description}

\begin{description}
\item[c.] Try fitting pressure as a quadratic function of volume. What do
you think of the fit?
\end{description}

\begin{description}
\item[d.] Give an (unbiased) estimate of the noise variance and use this
to compute the standard error on the quadratic coefficient. Based on
this, would you expect the quadratic model to be ``better'' than the
linear model?
\end{description}

\begin{description}
\item[e.] Use leave-one-out cross-validation to estimate the prediction
error of the linear and quadratic model. Compare them with theoretical
estimates of the prediction error that follow from the assumptions of
the standard statistical model and discuss the difference (if any).
\end{description}

\begin{description}
\item[f.] Try fitting higher order polynomials. What calibration model
would you propose based on your results?
\end{description}

\footnotetext{Exercise and data adapted from Rice, Mathematical Statistics and
Data Analysis, chapter 14, problem 35.}

\newpage

\section*{Hints for Exercise~5.}

\begin{description}
\item[a.] Store the data locally and read it into {\sc Matlab} with
{\tt load nuclear.dat}. The inputs $x$ are in the first column of {\tt
nuclear}, the outputs $y$ in the second column.
\end{description}

\begin{description}
\item[b.] You should be able to program least squares estimation
yourself from the description in chapter~5 of the syllabus. In fact,
{\sc Matlab}'s standard (left) matrix division routines do most of the job
for you. See e.g.\ {\tt help mldivide}. Alternatively, you can build your
solution around {\sc Matlab}'s {\tt solvelin} and {\tt simulin}. These
functions are, however, obsolete. To avoid annoying warning messages add
{\tt nntwarn off} before calling them. With {\tt type solvelin} you can
see how it works. {\tt solvelin} and {\tt simulin} assume that the inputs are
of size $\onder{n}{input} \times n$ and the outputs of size $\onder{n}{output}
\times n$ (with $n$ the number of samples and $\onder{n}{output} = 1$ in
this exercise). For this you may have to ``transpose'' your variables, e.g.,
with ${\tt x = x';}$.
\end{description}

\begin{description}
\item[c.] Note that you can still use linear regression! The quadratic model
can be cast in the form
\[
\hat{y}_i = \beta_0 + \beta_1 x_{1i} + \beta_2 x_{2i} \mbox{~~with~~}
x_{1i} = x_i \mbox{~~and~~} x_{2i} = x_i^2 \: .
\]
That is, redefine the matrix of inputs that you use, for example, as input
to {\tt solvelin} and apply exactly the same routines. In other words, the
``linear'' in linear regression stands for linear in the parameters of the
model $\vbeta$, not necessarily linear in the input $x$.
\end{description}

\begin{description}
\item[d.] See section~5.3 in the syllabus.
\end{description}

\begin{description}
\item[e.] See section~5.5 in the syllabus for a description of leave-one-out
cross-validation as well as for theoretical estimates of the prediction error.
Implement the cross-validation procedure as generally as you can.
\end{description}

\begin{description}
\item[f.] Obviously, we can also incorporate terms $x^3$, $x^4$, and so on,
in our regression function and still apply linear least squares. Build a
``forward selection procedure'' that adds terms until the (estimated)
prediction error no longer decreases.
\end{description}

\end{document}