\documentclass[]{article}
\usepackage{lmodern}
\usepackage{amssymb,amsmath}
\usepackage{ifxetex,ifluatex}
\usepackage{fixltx2e} % provides \textsubscript
\ifnum 0\ifxetex 1\fi\ifluatex 1\fi=0 % if pdftex
  \usepackage[T1]{fontenc}
  \usepackage[utf8]{inputenc}
\else % if luatex or xelatex
  \ifxetex
    \usepackage{mathspec}
  \else
    \usepackage{fontspec}
  \fi
  \defaultfontfeatures{Ligatures=TeX,Scale=MatchLowercase}
\fi
% use upquote if available, for straight quotes in verbatim environments
\IfFileExists{upquote.sty}{\usepackage{upquote}}{}
% use microtype if available
\IfFileExists{microtype.sty}{%
\usepackage{microtype}
\UseMicrotypeSet[protrusion]{basicmath} % disable protrusion for tt fonts
}{}
\usepackage[margin=1in]{geometry}
\usepackage{hyperref}
\hypersetup{unicode=true,
            pdftitle={Solutions to Recommended Exercises},
            pdfauthor={Thea Roksvåg, Julia Debik and Mette Langaas, Department of Mathematical Sciences, NTNU},
            pdfborder={0 0 0},
            breaklinks=true}
\urlstyle{same}  % don't use monospace font for urls
\usepackage{color}
\usepackage{fancyvrb}
\newcommand{\VerbBar}{|}
\newcommand{\VERB}{\Verb[commandchars=\\\{\}]}
\DefineVerbatimEnvironment{Highlighting}{Verbatim}{commandchars=\\\{\}}
% Add ',fontsize=\small' for more characters per line
\usepackage{framed}
\definecolor{shadecolor}{RGB}{248,248,248}
\newenvironment{Shaded}{\begin{snugshade}}{\end{snugshade}}
\newcommand{\KeywordTok}[1]{\textcolor[rgb]{0.13,0.29,0.53}{\textbf{#1}}}
\newcommand{\DataTypeTok}[1]{\textcolor[rgb]{0.13,0.29,0.53}{#1}}
\newcommand{\DecValTok}[1]{\textcolor[rgb]{0.00,0.00,0.81}{#1}}
\newcommand{\BaseNTok}[1]{\textcolor[rgb]{0.00,0.00,0.81}{#1}}
\newcommand{\FloatTok}[1]{\textcolor[rgb]{0.00,0.00,0.81}{#1}}
\newcommand{\ConstantTok}[1]{\textcolor[rgb]{0.00,0.00,0.00}{#1}}
\newcommand{\CharTok}[1]{\textcolor[rgb]{0.31,0.60,0.02}{#1}}
\newcommand{\SpecialCharTok}[1]{\textcolor[rgb]{0.00,0.00,0.00}{#1}}
\newcommand{\StringTok}[1]{\textcolor[rgb]{0.31,0.60,0.02}{#1}}
\newcommand{\VerbatimStringTok}[1]{\textcolor[rgb]{0.31,0.60,0.02}{#1}}
\newcommand{\SpecialStringTok}[1]{\textcolor[rgb]{0.31,0.60,0.02}{#1}}
\newcommand{\ImportTok}[1]{#1}
\newcommand{\CommentTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textit{#1}}}
\newcommand{\DocumentationTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textbf{\textit{#1}}}}
\newcommand{\AnnotationTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textbf{\textit{#1}}}}
\newcommand{\CommentVarTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textbf{\textit{#1}}}}
\newcommand{\OtherTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{#1}}
\newcommand{\FunctionTok}[1]{\textcolor[rgb]{0.00,0.00,0.00}{#1}}
\newcommand{\VariableTok}[1]{\textcolor[rgb]{0.00,0.00,0.00}{#1}}
\newcommand{\ControlFlowTok}[1]{\textcolor[rgb]{0.13,0.29,0.53}{\textbf{#1}}}
\newcommand{\OperatorTok}[1]{\textcolor[rgb]{0.81,0.36,0.00}{\textbf{#1}}}
\newcommand{\BuiltInTok}[1]{#1}
\newcommand{\ExtensionTok}[1]{#1}
\newcommand{\PreprocessorTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textit{#1}}}
\newcommand{\AttributeTok}[1]{\textcolor[rgb]{0.77,0.63,0.00}{#1}}
\newcommand{\RegionMarkerTok}[1]{#1}
\newcommand{\InformationTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textbf{\textit{#1}}}}
\newcommand{\WarningTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textbf{\textit{#1}}}}
\newcommand{\AlertTok}[1]{\textcolor[rgb]{0.94,0.16,0.16}{#1}}
\newcommand{\ErrorTok}[1]{\textcolor[rgb]{0.64,0.00,0.00}{\textbf{#1}}}
\newcommand{\NormalTok}[1]{#1}
\usepackage{graphicx,grffile}
\makeatletter
\def\maxwidth{\ifdim\Gin@nat@width>\linewidth\linewidth\else\Gin@nat@width\fi}
\def\maxheight{\ifdim\Gin@nat@height>\textheight\textheight\else\Gin@nat@height\fi}
\makeatother
% Scale images if necessary, so that they will not overflow the page
% margins by default, and it is still possible to overwrite the defaults
% using explicit options in \includegraphics[width, height, ...]{}
\setkeys{Gin}{width=\maxwidth,height=\maxheight,keepaspectratio}
\IfFileExists{parskip.sty}{%
\usepackage{parskip}
}{% else
\setlength{\parindent}{0pt}
\setlength{\parskip}{6pt plus 2pt minus 1pt}
}
\setlength{\emergencystretch}{3em}  % prevent overfull lines
\providecommand{\tightlist}{%
  \setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}
\setcounter{secnumdepth}{0}
% Redefines (sub)paragraphs to behave more like sections
\ifx\paragraph\undefined\else
\let\oldparagraph\paragraph
\renewcommand{\paragraph}[1]{\oldparagraph{#1}\mbox{}}
\fi
\ifx\subparagraph\undefined\else
\let\oldsubparagraph\subparagraph
\renewcommand{\subparagraph}[1]{\oldsubparagraph{#1}\mbox{}}
\fi

%%% Use protect on footnotes to avoid problems with footnotes in titles
\let\rmarkdownfootnote\footnote%
\def\footnote{\protect\rmarkdownfootnote}

%%% Change title format to be more compact
\usepackage{titling}

% Create subtitle command for use in maketitle
\newcommand{\subtitle}[1]{
  \posttitle{
    \begin{center}\large#1\end{center}
    }
}

\setlength{\droptitle}{-2em}

  \title{Solutions to Recommended Exercises}
    \pretitle{\vspace{\droptitle}\centering\huge}
  \posttitle{\par}
  \subtitle{TMA4268 Statistical Learning V2019. Module 2: STATISTICAL LEARNING}
  \author{Thea Roksvåg, Julia Debik and Mette Langaas, Department of Mathematical
Sciences, NTNU}
    \preauthor{\centering\large\emph}
  \postauthor{\par}
      \predate{\centering\large\emph}
  \postdate{\par}
    \date{14.01.2019}


\begin{document}
\maketitle

{
\setcounter{tocdepth}{2}
\tableofcontents
}
Last changes: (14.01.2019: first version)

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\section{Problem 1}\label{problem-1}

\begin{enumerate}
\def\labelenumi{\alph{enumi}.}
\tightlist
\item
  Observation 1: \(\sqrt{(3-1)^2+(3-2)^2}= \sqrt{5}\)\\
  Observation 2: \(\sqrt{(2-1)^2+(0-2)^2}= \sqrt{5}\)\\
  Observation 3: \(\sqrt{(1-1)^2+(1-2)^2}= 1\)\\
  Observation 4: \(\sqrt{(0-1)^2+(1-2)^2}= \sqrt{2}\)\\
  Observation 5: \(\sqrt{(-1-1)^2+(0-2)^2}= \sqrt{8}\)\\
  Observation 6: \(\sqrt{(2-1)^2+(1-2)^2}= \sqrt{2}\)\\
  Observation 7: \(\sqrt{(1-1)^2+(0-2)^2} = 2\)
\item
  The closest point is observation 3. This observation belongs to class
  A hence the predicted class membership for our test observation is A.
\item
  The four closest points are observations 3, 7, 4 and 6 with
  corresponding classes \{A, B, B, B\}. The predited class for our test
  observation is thus B.
\item
  If the Bayes decision boundary is highly non-linear, we would choose a
  \(K\) that is not to big, as the decision boundary becomes
  approximately linear when \(K\) tends to infinity
\item
  We start by installing and loading the \texttt{ggplot2} packages.
\end{enumerate}

\begin{Shaded}
\begin{Highlighting}[]
\CommentTok{#install.packages(ggplot2)}
\KeywordTok{library}\NormalTok{(ggplot2)}
\end{Highlighting}
\end{Shaded}

We make a data frame with our observations

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{knnframe =}\StringTok{ }\KeywordTok{data.frame}\NormalTok{(}\DataTypeTok{x1 =} \KeywordTok{c}\NormalTok{(}\DecValTok{3}\NormalTok{, }\DecValTok{2}\NormalTok{, }\DecValTok{1}\NormalTok{, }\DecValTok{0}\NormalTok{, }\OperatorTok{-}\DecValTok{1}\NormalTok{, }\DecValTok{2}\NormalTok{, }\DecValTok{1}\NormalTok{), }\DataTypeTok{x2 =} \KeywordTok{c}\NormalTok{(}\DecValTok{3}\NormalTok{, }\DecValTok{0}\NormalTok{, }\DecValTok{1}\NormalTok{, }\DecValTok{1}\NormalTok{, }\DecValTok{0}\NormalTok{, }\DecValTok{1}\NormalTok{, }\DecValTok{0}\NormalTok{),  }\DataTypeTok{y=}\KeywordTok{as.factor}\NormalTok{(}\KeywordTok{c}\NormalTok{(}\StringTok{"A"}\NormalTok{, }\StringTok{"A"}\NormalTok{, }\StringTok{"A"}\NormalTok{, }\StringTok{"B"}\NormalTok{, }\StringTok{"B"}\NormalTok{, }\StringTok{"B"}\NormalTok{, }\StringTok{"B"}\NormalTok{)))}
\end{Highlighting}
\end{Shaded}

We plot the observations unsing the \texttt{ggplot} function. We set
\texttt{color=y} to obtain a colored response.

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{ggplot}\NormalTok{(knnframe, }\KeywordTok{aes}\NormalTok{(}\DataTypeTok{x=}\NormalTok{x1, }\DataTypeTok{y=}\NormalTok{x2, }\DataTypeTok{color=}\NormalTok{y))}\OperatorTok{+}\KeywordTok{geom_point}\NormalTok{()}\OperatorTok{+}\KeywordTok{theme_minimal}\NormalTok{()}
\end{Highlighting}
\end{Shaded}

\includegraphics{2RecEx-sol_files/figure-latex/unnamed-chunk-3-1.pdf}

\begin{enumerate}
\def\labelenumi{\alph{enumi}.}
\setcounter{enumi}{5}
\item
\end{enumerate}

\begin{Shaded}
\begin{Highlighting}[]
\CommentTok{#install.packages(class)}
\KeywordTok{library}\NormalTok{(class)}
\KeywordTok{knn}\NormalTok{(}\DataTypeTok{train =}\NormalTok{ knnframe[,}\DecValTok{1}\OperatorTok{:}\DecValTok{2}\NormalTok{], }\DataTypeTok{cl =}\NormalTok{ knnframe[,}\DecValTok{3}\NormalTok{], }\DataTypeTok{test =} \KeywordTok{c}\NormalTok{(}\DecValTok{1}\NormalTok{, }\DecValTok{2}\NormalTok{), }\DataTypeTok{k=}\DecValTok{1}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## [1] A
## attr(,"nn.index")
##      [,1]
## [1,]    3
## attr(,"nn.dist")
##      [,1]
## [1,]    1
## Levels: A
\end{verbatim}

\begin{enumerate}
\def\labelenumi{\alph{enumi}.}
\setcounter{enumi}{6}
\item
\end{enumerate}

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{knn}\NormalTok{(}\DataTypeTok{train =}\NormalTok{ knnframe[,}\DecValTok{1}\OperatorTok{:}\DecValTok{2}\NormalTok{], }\DataTypeTok{cl =}\NormalTok{ knnframe[,}\DecValTok{3}\NormalTok{], }\DataTypeTok{test =} \KeywordTok{c}\NormalTok{(}\DecValTok{1}\NormalTok{, }\DecValTok{2}\NormalTok{), }\DataTypeTok{k=}\DecValTok{4}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## [1] B
## attr(,"nn.index")
##      [,1] [,2] [,3] [,4]
## [1,]    3    6    4    7
## attr(,"nn.dist")
##      [,1]     [,2]     [,3] [,4]
## [1,]    1 1.414214 1.414214    2
## Levels: B
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{knn}\NormalTok{(}\DataTypeTok{train =}\NormalTok{ knnframe[,}\DecValTok{1}\OperatorTok{:}\DecValTok{2}\NormalTok{], }\DataTypeTok{cl =}\NormalTok{ knnframe[,}\DecValTok{3}\NormalTok{], }\DataTypeTok{test =} \KeywordTok{c}\NormalTok{(}\DecValTok{1}\NormalTok{, }\DecValTok{2}\NormalTok{), }\DataTypeTok{k=}\DecValTok{7}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## [1] B
## attr(,"nn.index")
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,]    3    6    4    7    1    2    5
## attr(,"nn.dist")
##      [,1]     [,2]     [,3] [,4]     [,5]     [,6]     [,7]
## [1,]    1 1.414214 1.414214    2 2.236068 2.236068 2.828427
## Levels: B
\end{verbatim}

\section{Problem 2}\label{problem-2}

We consider a regression problem, where the true underlying curve is
\(f(x)=-x+x^2+x^3\) and we are considering \(x \in [-3,3]\).

This non-linear curve is only observed with added noise (either a random
phenomenon, or unobservable variables influence the observations), that
is, we observe \(y=f(x)+\varepsilon\). In our example the error is
sampled from \(\varepsilon\sim N(0,2^2)\).

\begin{figure}
\centering
\includegraphics{Prob1f1.png}
\caption{Figure 1}
\end{figure}

In real life we are presented with a data set of pairs \((x_i,y_i)\),
\(i=1,\ldots,n\), and asked to provide a prediction at a value \(x\). We
will use the method of K nearest neighbour regression to do this here.

We have a training set of \(n=61\) observations \((x_i,y_i)\),
\(i=1,\ldots,n\). The KNN regression method provides a prediction at a
value \(x\) by finding the closes \(K\) points and calculating the
average of the observed \(y\) values at these points.

In addition we have a test set of \(n=61\) observations (at the same
grid points as for the training set), but now with new observed values
\(y\).

We have considered \(K=1,\ldots,25\) in the KNN method. Our experiment
has been repeated \(M=1000\) times (that is, \(M\) versions of training
and test set).

\subsection{a) Training and test MSE}\label{a-training-and-test-mse}

\begin{figure}
\centering
\includegraphics{Prob1f2.png}
\caption{Figure 2}
\end{figure}

In the Figure 2 (above) you see the result of applying the KNN method
with \(K=1,2,10,25\) to our training data, repeated for \(M\) different
training sets (blue lines). The black lines show the true underlying
curve.

\begin{itemize}
\tightlist
\item
  Comment briefly on what you see.
\item
  Does a high or low value of \(K\) give the most flexible fit?
\end{itemize}

\subsubsection{Answers:}\label{answers}

{Here \(K=1\) gives predicted values that on average (over the training
sets) are not far from the true curve, but the large variability.
Increasing the number of neighbours to \(K=2\) gives less variance in
the preditions, and still on average are not far from the true curve.
Increasing the number of neighbours further to \(K=10\), the variance is
even more reduced, but the valued on the edge of our training set (at
=-3 and x=3) are not well estimated, and \(K=25\) gives little
flexibility to the curve and doesn't fit well, especially on the edges
of the data set. Hence, a low value of \(K\) gives the most flexible
(complex) fit.}

In Figure 3 (below) you see mean-squared errors (mean of squared
differences between observed and fitted values) for the training set and
for the test set (right panel for one training and one test set, and
left panel for \(M\)).

\begin{itemize}
\tightlist
\item
  Comment on what you see.
\item
  What do you think is the ``best'' choice for K?
\end{itemize}

\subsubsection{Answers:}\label{answers-1}

{ We see that larger values of \(K\) give larger variations in the MSE
for both the test and the train set than for the smaller choices of
\(K\). For the training set, the MSE increases linearly with \(K\),
i.e.~worse fit, but for the test set, there seems to be a small decrease
at \(K=5\) before it increases again. Based on this, we would choose
\(K=5\) as the best choice for \(K\). }

Remark: in real life we do not know the true curve, and need to use the
test data to decide on model flexibility (choosing \(K\)).

\begin{figure}
\centering
\includegraphics{Prob1f3.png}
\caption{Figure 3}
\end{figure}

\subsection{b) Bias-variance trade-off}\label{b-bias-variance-trade-off}

Now we leave the real world situation, and assume we know the truth
(this is to focus on bias-variance trade-off). You will not observe
these curves in real life - but the understanding of the bias-variance
trade-off is a core skill in this course!

In the Figure 4 (below) you see a plot of estimated squared bias,
estimated variance, true irreducible error and the sum of these
(labelled total) and averaged over all values of \(x\)

The the squared bias and the variance is calculated based on the
predicted values and the ``true'' values (without the added noise) at
each \(x\).

\begin{itemize}
\tightlist
\item
  Explain how that is done. Hint: this is what the \(M\) repeated
  training data sets are used for.
\item
  Focus on Figure 4. As the flexibility of the model increases (\(K\)
  decreases), what happens with

  \begin{itemize}
  \tightlist
  \item
    the squared bias,\\
  \item
    the variance, and\\
  \item
    the irreducible error?
  \end{itemize}
\item
  What would you recommend is the optimal value of \(K\)? Is this in
  agreement with what you found in a)?
\end{itemize}

\subsubsection{Answers:}\label{answers-2}

\begin{itemize}
\item
  {How this is done? The bias and variance of \(\hat{f(x_0)}\) - when
  the true value is \(f(x_0)\) is defined by \[
  \text{Bias}(\hat f(x_0)) =\text{E}[\hat f(x_0) - f(x_0)]
  \] \[
  \text{Var}(\hat f(x_0)) = \text{E}[\hat f(x_0)^2] - \text{E}[\hat f(x_0)]^2
  \] Using the M resamples of the test data, we estimate the mean of
  \(\hat f(x_0)\) by the average over the \(M\) predicted values at
  \(x_0\), \(\text{Ave}(\hat{f}(x_0))\), and the variance by the
  emipircal variance
  \(\frac{1}{M-1} \sum_{i=1}^M (\hat{f}(x_0)-\text{Ave}(\hat{f}(x_0)))^2\).}
  \[
  \text{E}\left[\left(Y-\hat{f}\left(x_0\right)\right)^2\right]=\underbrace{\text{E}(\epsilon)}_{\text{Irredusible error}}+\underbrace{\text{Var}\left(\hat{f}\left(x_0\right)\right)}_{\text{Variance}}+\underbrace{\left[\text{E}\left(\hat{f}\left(x_0\right)\right)-f(x_0))\right]^2}_{\text{Squared bias}}
  \]
\item
  { From Figure 4, we se that as K increases the squared bias increases
  rapidly, the variance is slightly reduced, but the irreducible error
  stays constant. Looking at the total MSE (squared bias plus variance)
  we see that the optimal K-value lies somewhere between 3 and 5.}
\item
  { Extra: We observe that the ``optimal'' \(K\) varies for different
  \(x_0\)s. For \(x_0=-2\) the best \(K\) is between 5 and 10. For
  \(x_0=0\) most values of \(K\) is ok, and for \(x_0=1\) the best is
  around \(K=10\). For \(x_0=2.5\) \(K\) needs to be below 15. The
  reason for showing this is to see some of the \(x_0\) values in our
  interval \([-3,3]\) that is given equal weight in deciding the optimal
  \(K\).}
\end{itemize}

Extra: We have chosen to also plot curves at four values of \(x\) -
Figure 5 (below). Based on these four curves, that would you recommend
is the optimal value of \(K\)? Is this in agreement with what you found
previously (averaged over \(x\))?

\begin{figure}
\centering
\includegraphics{Prob1f4.png}
\caption{Figure 4}
\end{figure}

\begin{figure}
\centering
\includegraphics{Prob1f5.png}
\caption{Figure 5}
\end{figure}

\section{Problem 3}\label{problem-3}

\subsection{a)}\label{a}

We sample from the model \(y=x^2+\epsilon\) where
\(\epsilon \sim \mathcal{N}(0,2^2)\) and
\(x\in \{-2,-1.9,-1.8,...,3.8,3.9,4\}\). This means that
\(y \sim \mathcal{N}(x^2,2^2)\). A total of 100 samples from this model
are generated for each of the 61 \(x\)'s.

See comments in code for further explanations.

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{library}\NormalTok{(ggplot2)}
\KeywordTok{library}\NormalTok{(ggpubr)}
\KeywordTok{set.seed}\NormalTok{(}\DecValTok{2}\NormalTok{) }\CommentTok{# to reproduce}

\NormalTok{M=}\DecValTok{100} \CommentTok{# repeated samplings, x fixed }
\NormalTok{nord=}\DecValTok{20} \CommentTok{# order of polynoms}


\NormalTok{x =}\StringTok{ }\KeywordTok{seq}\NormalTok{(}\OperatorTok{-}\DecValTok{2}\NormalTok{, }\DecValTok{4}\NormalTok{, }\FloatTok{0.1}\NormalTok{) }\CommentTok{#We make a sequence of 61 points, x. These are the points for which we evaluate the function f(x).}
\NormalTok{truefunc=}\ControlFlowTok{function}\NormalTok{(x) }\KeywordTok{return}\NormalTok{(x}\OperatorTok{^}\DecValTok{2}\NormalTok{) }\CommentTok{#The true f(x)=x^2. }
\NormalTok{true_y =}\StringTok{ }\KeywordTok{truefunc}\NormalTok{(x) }\CommentTok{#We find f(x) for each element in vector x.}

\NormalTok{error =}\StringTok{ }\KeywordTok{matrix}\NormalTok{(}\KeywordTok{rnorm}\NormalTok{(}\KeywordTok{length}\NormalTok{(x)}\OperatorTok{*}\NormalTok{M, }\DataTypeTok{mean=}\DecValTok{0}\NormalTok{, }\DataTypeTok{sd=}\DecValTok{2}\NormalTok{),}\DataTypeTok{nrow=}\NormalTok{M,}\DataTypeTok{byrow=}\OtherTok{TRUE}\NormalTok{) }\CommentTok{#Noise (epsilon) is sampled from a normal distribution and stored in this matrix. Each column corresponds to one value of x.}
\NormalTok{ymat =}\StringTok{ }\KeywordTok{matrix}\NormalTok{(}\KeywordTok{rep}\NormalTok{(true_y,M),}\DataTypeTok{byrow=}\NormalTok{T,}\DataTypeTok{nrow=}\NormalTok{M) }\OperatorTok{+}\StringTok{ }\NormalTok{error }\CommentTok{#The 100 samples or the observations are stored in this matrix.}

\NormalTok{predarray=}\KeywordTok{array}\NormalTok{(}\OtherTok{NA}\NormalTok{,}\DataTypeTok{dim=}\KeywordTok{c}\NormalTok{(M,}\KeywordTok{length}\NormalTok{(x),nord))}
\ControlFlowTok{for}\NormalTok{ (i }\ControlFlowTok{in} \DecValTok{1}\OperatorTok{:}\NormalTok{M)}
\NormalTok{\{}
  \ControlFlowTok{for}\NormalTok{ (j }\ControlFlowTok{in} \DecValTok{1}\OperatorTok{:}\NormalTok{nord)}
\NormalTok{  \{}
\NormalTok{    predarray[i,,j]=}\KeywordTok{predict}\NormalTok{(}\KeywordTok{lm}\NormalTok{(ymat[i,]}\OperatorTok{~}\KeywordTok{poly}\NormalTok{(x, j,}\DataTypeTok{raw=}\OtherTok{TRUE}\NormalTok{)))}
    \CommentTok{#Based on the response y_i and the x_i's, we fit a polynomial model of degre 1,...,20. This means that we assume y_i~Normal(x_i^j,0). }
\NormalTok{  \}}
\NormalTok{\}}
\CommentTok{# M matrices of size length(x) times nord}
\CommentTok{# first, only look at variablity in the M fits and plot M curves where we had 1.}

\CommentTok{# for plotting need to stack the matrices underneath eachother and make new variable "rep"}
\NormalTok{stackmat=}\OtherTok{NULL}
\ControlFlowTok{for}\NormalTok{ (i }\ControlFlowTok{in} \DecValTok{1}\OperatorTok{:}\NormalTok{M) stackmat=}\KeywordTok{rbind}\NormalTok{(stackmat,}\KeywordTok{cbind}\NormalTok{(x,}\KeywordTok{rep}\NormalTok{(i,}\KeywordTok{length}\NormalTok{(x)),predarray[i,,]))}
\CommentTok{#dim(stackmat)}
\KeywordTok{colnames}\NormalTok{(stackmat)=}\KeywordTok{c}\NormalTok{(}\StringTok{"x"}\NormalTok{,}\StringTok{"rep"}\NormalTok{,}\KeywordTok{paste}\NormalTok{(}\StringTok{"poly"}\NormalTok{,}\DecValTok{1}\OperatorTok{:}\DecValTok{20}\NormalTok{,}\DataTypeTok{sep=}\StringTok{""}\NormalTok{))}
\NormalTok{sdf=}\KeywordTok{as.data.frame}\NormalTok{(stackmat) }\CommentTok{#NB have poly1-20 now - but first only use 1,2,20}
\CommentTok{# to add true curve using stat_function - easiest solution}
\NormalTok{true_x=x}
\NormalTok{yrange=}\KeywordTok{range}\NormalTok{(}\KeywordTok{apply}\NormalTok{(sdf,}\DecValTok{2}\NormalTok{,range)[,}\DecValTok{3}\OperatorTok{:}\DecValTok{22}\NormalTok{])}
\NormalTok{p1=}\KeywordTok{ggplot}\NormalTok{(}\DataTypeTok{data=}\NormalTok{sdf,}\KeywordTok{aes}\NormalTok{(}\DataTypeTok{x=}\NormalTok{x,}\DataTypeTok{y=}\NormalTok{poly1,}\DataTypeTok{group=}\NormalTok{rep,}\DataTypeTok{colour=}\NormalTok{rep))}\OperatorTok{+}\KeywordTok{scale_y_continuous}\NormalTok{(}\DataTypeTok{limits=}\NormalTok{yrange)}\OperatorTok{+}\KeywordTok{geom_line}\NormalTok{()}
\NormalTok{p1=p1}\OperatorTok{+}\KeywordTok{stat_function}\NormalTok{(}\DataTypeTok{fun=}\NormalTok{truefunc,}\DataTypeTok{lwd=}\FloatTok{1.3}\NormalTok{,}\DataTypeTok{colour=}\StringTok{"black"}\NormalTok{)}\OperatorTok{+}\KeywordTok{ggtitle}\NormalTok{(}\StringTok{"poly1"}\NormalTok{)}\OperatorTok{+}\KeywordTok{theme_minimal}\NormalTok{()}
\NormalTok{p2=}\KeywordTok{ggplot}\NormalTok{(}\DataTypeTok{data=}\NormalTok{sdf,}\KeywordTok{aes}\NormalTok{(}\DataTypeTok{x=}\NormalTok{x,}\DataTypeTok{y=}\NormalTok{poly2,}\DataTypeTok{group=}\NormalTok{rep,}\DataTypeTok{colour=}\NormalTok{rep))}\OperatorTok{+}\KeywordTok{scale_y_continuous}\NormalTok{(}\DataTypeTok{limits=}\NormalTok{yrange)}\OperatorTok{+}\KeywordTok{geom_line}\NormalTok{()}
\NormalTok{p2=p2}\OperatorTok{+}\KeywordTok{stat_function}\NormalTok{(}\DataTypeTok{fun=}\NormalTok{truefunc,}\DataTypeTok{lwd=}\FloatTok{1.3}\NormalTok{,}\DataTypeTok{colour=}\StringTok{"black"}\NormalTok{)}\OperatorTok{+}\KeywordTok{ggtitle}\NormalTok{(}\StringTok{"poly2"}\NormalTok{)}\OperatorTok{+}\KeywordTok{theme_minimal}\NormalTok{()}
\NormalTok{p10=}\KeywordTok{ggplot}\NormalTok{(}\DataTypeTok{data=}\NormalTok{sdf,}\KeywordTok{aes}\NormalTok{(}\DataTypeTok{x=}\NormalTok{x,}\DataTypeTok{y=}\NormalTok{poly10,}\DataTypeTok{group=}\NormalTok{rep,}\DataTypeTok{colour=}\NormalTok{rep))}\OperatorTok{+}\KeywordTok{scale_y_continuous}\NormalTok{(}\DataTypeTok{limits=}\NormalTok{yrange)}\OperatorTok{+}\KeywordTok{geom_line}\NormalTok{()}
\NormalTok{p10=p10}\OperatorTok{+}\KeywordTok{stat_function}\NormalTok{(}\DataTypeTok{fun=}\NormalTok{truefunc,}\DataTypeTok{lwd=}\FloatTok{1.3}\NormalTok{,}\DataTypeTok{colour=}\StringTok{"black"}\NormalTok{)}\OperatorTok{+}\KeywordTok{ggtitle}\NormalTok{(}\StringTok{"poly10"}\NormalTok{)}\OperatorTok{+}\KeywordTok{theme_minimal}\NormalTok{()}
\NormalTok{p20=}\KeywordTok{ggplot}\NormalTok{(}\DataTypeTok{data=}\NormalTok{sdf,}\KeywordTok{aes}\NormalTok{(}\DataTypeTok{x=}\NormalTok{x,}\DataTypeTok{y=}\NormalTok{poly20,}\DataTypeTok{group=}\NormalTok{rep,}\DataTypeTok{colour=}\NormalTok{rep))}\OperatorTok{+}\KeywordTok{scale_y_continuous}\NormalTok{(}\DataTypeTok{limits=}\NormalTok{yrange)}\OperatorTok{+}\KeywordTok{geom_line}\NormalTok{()}
\NormalTok{p20=p20}\OperatorTok{+}\KeywordTok{stat_function}\NormalTok{(}\DataTypeTok{fun=}\NormalTok{truefunc,}\DataTypeTok{lwd=}\FloatTok{1.3}\NormalTok{,}\DataTypeTok{colour=}\StringTok{"black"}\NormalTok{)}\OperatorTok{+}\KeywordTok{ggtitle}\NormalTok{(}\StringTok{"poly20"}\NormalTok{)}\OperatorTok{+}\KeywordTok{theme_minimal}\NormalTok{()}
\KeywordTok{ggarrange}\NormalTok{(p1,p2,p10,p20)}
\end{Highlighting}
\end{Shaded}

\includegraphics{2RecEx-sol_files/figure-latex/unnamed-chunk-7-1.pdf}

The upper left plot shows 100 predictions when we assume that \(y\) is a
linear function of \(x\), the upper right plot hows 100 predictions when
we assume that \(y\) is function of poynomials up to \(x^2\), the lower
left plot shows 100 predictions when we assume \(y\) is a function of
polynomials up to \(x^{10}\) and the lower right plot shows 100
predictions when assuming \(y\) is a function of polynomials up to
\(x^{20}\).

\subsection{b)}\label{b}

Run the code attached and consider the following plots:

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{set.seed}\NormalTok{(}\DecValTok{2}\NormalTok{) }\CommentTok{# to reproduce}

\NormalTok{M=}\DecValTok{100} \CommentTok{# repeated samplings,x fixed but new errors}
\NormalTok{nord=}\DecValTok{20}
\NormalTok{x =}\StringTok{ }\KeywordTok{seq}\NormalTok{(}\OperatorTok{-}\DecValTok{2}\NormalTok{, }\DecValTok{4}\NormalTok{, }\FloatTok{0.1}\NormalTok{)}
\NormalTok{truefunc=}\ControlFlowTok{function}\NormalTok{(x) }\KeywordTok{return}\NormalTok{(x}\OperatorTok{^}\DecValTok{2}\NormalTok{)}
\NormalTok{true_y =}\StringTok{ }\KeywordTok{truefunc}\NormalTok{(x)}

\NormalTok{error =}\StringTok{ }\KeywordTok{matrix}\NormalTok{(}\KeywordTok{rnorm}\NormalTok{(}\KeywordTok{length}\NormalTok{(x)}\OperatorTok{*}\NormalTok{M, }\DataTypeTok{mean=}\DecValTok{0}\NormalTok{, }\DataTypeTok{sd=}\DecValTok{2}\NormalTok{),}\DataTypeTok{nrow=}\NormalTok{M,}\DataTypeTok{byrow=}\OtherTok{TRUE}\NormalTok{)}
\NormalTok{testerror =}\StringTok{ }\KeywordTok{matrix}\NormalTok{(}\KeywordTok{rnorm}\NormalTok{(}\KeywordTok{length}\NormalTok{(x)}\OperatorTok{*}\NormalTok{M, }\DataTypeTok{mean=}\DecValTok{0}\NormalTok{, }\DataTypeTok{sd=}\DecValTok{2}\NormalTok{),}\DataTypeTok{nrow=}\NormalTok{M,}\DataTypeTok{byrow=}\OtherTok{TRUE}\NormalTok{)}
\NormalTok{ymat =}\StringTok{ }\KeywordTok{matrix}\NormalTok{(}\KeywordTok{rep}\NormalTok{(true_y,M),}\DataTypeTok{byrow=}\NormalTok{T,}\DataTypeTok{nrow=}\NormalTok{M) }\OperatorTok{+}\StringTok{ }\NormalTok{error}
\NormalTok{testymat =}\StringTok{ }\KeywordTok{matrix}\NormalTok{(}\KeywordTok{rep}\NormalTok{(true_y,M),}\DataTypeTok{byrow=}\NormalTok{T,}\DataTypeTok{nrow=}\NormalTok{M) }\OperatorTok{+}\StringTok{ }\NormalTok{testerror}

\NormalTok{predarray=}\KeywordTok{array}\NormalTok{(}\OtherTok{NA}\NormalTok{,}\DataTypeTok{dim=}\KeywordTok{c}\NormalTok{(M,}\KeywordTok{length}\NormalTok{(x),nord))}
\ControlFlowTok{for}\NormalTok{ (i }\ControlFlowTok{in} \DecValTok{1}\OperatorTok{:}\NormalTok{M)}
\NormalTok{\{}
  \ControlFlowTok{for}\NormalTok{ (j }\ControlFlowTok{in} \DecValTok{1}\OperatorTok{:}\NormalTok{nord)}
\NormalTok{  \{}
\NormalTok{    predarray[i,,j]=}\KeywordTok{predict}\NormalTok{(}\KeywordTok{lm}\NormalTok{(ymat[i,]}\OperatorTok{~}\KeywordTok{poly}\NormalTok{(x, j,}\DataTypeTok{raw=}\OtherTok{TRUE}\NormalTok{)))}
\NormalTok{  \}}
\NormalTok{\}  }
\NormalTok{trainMSE=}\KeywordTok{matrix}\NormalTok{(}\DataTypeTok{ncol=}\NormalTok{nord,}\DataTypeTok{nrow=}\NormalTok{M)}
\ControlFlowTok{for}\NormalTok{ (i }\ControlFlowTok{in} \DecValTok{1}\OperatorTok{:}\NormalTok{M) trainMSE[i,]=}\KeywordTok{apply}\NormalTok{((predarray[i,,]}\OperatorTok{-}\NormalTok{ymat[i,])}\OperatorTok{^}\DecValTok{2}\NormalTok{,}\DecValTok{2}\NormalTok{,mean)}
\NormalTok{testMSE=}\KeywordTok{matrix}\NormalTok{(}\DataTypeTok{ncol=}\NormalTok{nord,}\DataTypeTok{nrow=}\NormalTok{M)}
\ControlFlowTok{for}\NormalTok{ (i }\ControlFlowTok{in} \DecValTok{1}\OperatorTok{:}\NormalTok{M) testMSE[i,]=}\KeywordTok{apply}\NormalTok{((predarray[i,,]}\OperatorTok{-}\NormalTok{testymat[i,])}\OperatorTok{^}\DecValTok{2}\NormalTok{,}\DecValTok{2}\NormalTok{,mean)}

\KeywordTok{library}\NormalTok{(ggplot2)}
\KeywordTok{library}\NormalTok{(ggpubr)}

\CommentTok{# format suitable for plotting }
\NormalTok{stackmat=}\OtherTok{NULL}
\ControlFlowTok{for}\NormalTok{ (i }\ControlFlowTok{in} \DecValTok{1}\OperatorTok{:}\NormalTok{M) stackmat=}\KeywordTok{rbind}\NormalTok{(stackmat,}\KeywordTok{cbind}\NormalTok{(}\KeywordTok{rep}\NormalTok{(i,nord),}\DecValTok{1}\OperatorTok{:}\NormalTok{nord,trainMSE[i,],testMSE[i,]))}
\KeywordTok{colnames}\NormalTok{(stackmat)=}\KeywordTok{c}\NormalTok{(}\StringTok{"rep"}\NormalTok{,}\StringTok{"poly"}\NormalTok{,}\StringTok{"trainMSE"}\NormalTok{,}\StringTok{"testMSE"}\NormalTok{)}
\NormalTok{sdf=}\KeywordTok{as.data.frame}\NormalTok{(stackmat) }
\NormalTok{yrange=}\KeywordTok{range}\NormalTok{(sdf[,}\DecValTok{3}\OperatorTok{:}\DecValTok{4}\NormalTok{])}
\NormalTok{p1=}\KeywordTok{ggplot}\NormalTok{(}\DataTypeTok{data=}\NormalTok{sdf[}\DecValTok{1}\OperatorTok{:}\NormalTok{nord,],}\KeywordTok{aes}\NormalTok{(}\DataTypeTok{x=}\NormalTok{poly,}\DataTypeTok{y=}\NormalTok{trainMSE))}\OperatorTok{+}\KeywordTok{scale_y_continuous}\NormalTok{(}\DataTypeTok{limits=}\NormalTok{yrange)}\OperatorTok{+}\KeywordTok{geom_line}\NormalTok{()}\OperatorTok{+}\KeywordTok{theme_minimal}\NormalTok{()}
\NormalTok{pall=}\StringTok{ }\KeywordTok{ggplot}\NormalTok{(}\DataTypeTok{data=}\NormalTok{sdf,}\KeywordTok{aes}\NormalTok{(}\DataTypeTok{x=}\NormalTok{poly,}\DataTypeTok{group=}\NormalTok{rep,}\DataTypeTok{y=}\NormalTok{trainMSE,}\DataTypeTok{colour=}\NormalTok{rep))}\OperatorTok{+}\KeywordTok{scale_y_continuous}\NormalTok{(}\DataTypeTok{limits=}\NormalTok{yrange)}\OperatorTok{+}\KeywordTok{geom_line}\NormalTok{()}\OperatorTok{+}\KeywordTok{theme_minimal}\NormalTok{()}
\NormalTok{testp1=}\KeywordTok{ggplot}\NormalTok{(}\DataTypeTok{data=}\NormalTok{sdf[}\DecValTok{1}\OperatorTok{:}\NormalTok{nord,],}\KeywordTok{aes}\NormalTok{(}\DataTypeTok{x=}\NormalTok{poly,}\DataTypeTok{y=}\NormalTok{testMSE))}\OperatorTok{+}\KeywordTok{scale_y_continuous}\NormalTok{(}\DataTypeTok{limits=}\NormalTok{yrange)}\OperatorTok{+}\KeywordTok{geom_line}\NormalTok{()}\OperatorTok{+}\KeywordTok{theme_minimal}\NormalTok{()}
\NormalTok{testpall=}\StringTok{ }\KeywordTok{ggplot}\NormalTok{(}\DataTypeTok{data=}\NormalTok{sdf,}\KeywordTok{aes}\NormalTok{(}\DataTypeTok{x=}\NormalTok{poly,}\DataTypeTok{group=}\NormalTok{rep,}\DataTypeTok{y=}\NormalTok{testMSE,}\DataTypeTok{colour=}\NormalTok{rep))}\OperatorTok{+}\KeywordTok{scale_y_continuous}\NormalTok{(}\DataTypeTok{limits=}\NormalTok{yrange)}\OperatorTok{+}\KeywordTok{geom_line}\NormalTok{()}\OperatorTok{+}\KeywordTok{theme_minimal}\NormalTok{()}
\KeywordTok{ggarrange}\NormalTok{(p1,pall,testp1,testpall)}
\end{Highlighting}
\end{Shaded}

\includegraphics{2RecEx-sol_files/figure-latex/unnamed-chunk-8-1.pdf}

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{library}\NormalTok{(reshape2)}
\NormalTok{df=}\KeywordTok{melt}\NormalTok{(sdf,}\DataTypeTok{id=}\KeywordTok{c}\NormalTok{(}\StringTok{"poly"}\NormalTok{,}\StringTok{"rep"}\NormalTok{))[,}\OperatorTok{-}\DecValTok{2}\NormalTok{]}
\KeywordTok{colnames}\NormalTok{(df)[}\DecValTok{2}\NormalTok{]=}\StringTok{"MSEtype"}
\KeywordTok{ggplot}\NormalTok{(}\DataTypeTok{data=}\NormalTok{df,}\KeywordTok{aes}\NormalTok{(}\DataTypeTok{x=}\KeywordTok{as.factor}\NormalTok{(poly),}\DataTypeTok{y=}\NormalTok{value))}\OperatorTok{+}\KeywordTok{geom_boxplot}\NormalTok{(}\KeywordTok{aes}\NormalTok{(}\DataTypeTok{fill=}\NormalTok{MSEtype))}\OperatorTok{+}\KeywordTok{theme_minimal}\NormalTok{()}
\end{Highlighting}
\end{Shaded}

\includegraphics{2RecEx-sol_files/figure-latex/unnamed-chunk-8-2.pdf}

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{trainMSEmean=}\KeywordTok{apply}\NormalTok{(trainMSE,}\DecValTok{2}\NormalTok{,mean)}
\NormalTok{testMSEmean=}\KeywordTok{apply}\NormalTok{(testMSE,}\DecValTok{2}\NormalTok{,mean)}
\NormalTok{meandf=}\KeywordTok{melt}\NormalTok{(}\KeywordTok{data.frame}\NormalTok{(}\KeywordTok{cbind}\NormalTok{(}\StringTok{"poly"}\NormalTok{=}\DecValTok{1}\OperatorTok{:}\NormalTok{nord,trainMSEmean,testMSEmean)),}\DataTypeTok{id=}\StringTok{"poly"}\NormalTok{)}
\KeywordTok{ggplot}\NormalTok{(}\DataTypeTok{data=}\NormalTok{meandf,}\KeywordTok{aes}\NormalTok{(}\DataTypeTok{x=}\NormalTok{poly,}\DataTypeTok{y=}\NormalTok{value,}\DataTypeTok{colour=}\NormalTok{variable))}\OperatorTok{+}\KeywordTok{geom_line}\NormalTok{()}\OperatorTok{+}\KeywordTok{theme_minimal}\NormalTok{()}
\end{Highlighting}
\end{Shaded}

\includegraphics{2RecEx-sol_files/figure-latex/unnamed-chunk-8-3.pdf}

The plots show that the test MSE in general is larger than the train
MSE. This is reasonable. The fitted model is fitted based on the
training set. Thus, the error will be smaller for the train data than
for the test data. Furthermore, the plots show that the difference
between the MSE for the test set and the training set increases when the
degree of the polynomial increases. When the degree of the polynomial
increases, we get a more flexible model. The fitted curve will try to
pass through the training data if possible, which typically leads to an
overfitted model that performs bad for test data.

\subsection{c)}\label{c}

Run the code and consider the following plots:

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{meanmat=}\KeywordTok{matrix}\NormalTok{(}\DataTypeTok{ncol=}\KeywordTok{length}\NormalTok{(x),}\DataTypeTok{nrow=}\NormalTok{nord)}
\NormalTok{varmat=}\KeywordTok{matrix}\NormalTok{(}\DataTypeTok{ncol=}\KeywordTok{length}\NormalTok{(x),}\DataTypeTok{nrow=}\NormalTok{nord)}
\ControlFlowTok{for}\NormalTok{ (j }\ControlFlowTok{in} \DecValTok{1}\OperatorTok{:}\NormalTok{nord)}
\NormalTok{\{}
\NormalTok{  meanmat[j,]=}\KeywordTok{apply}\NormalTok{(predarray[,,j],}\DecValTok{2}\NormalTok{,mean) }\CommentTok{# we now take the mean over the M simulations - to mimic E and Var at each x value and each poly model}
\NormalTok{  varmat[j,]=}\KeywordTok{apply}\NormalTok{(predarray[,,j],}\DecValTok{2}\NormalTok{,var)}
\NormalTok{\}}
\CommentTok{# nord times length(x)}
\NormalTok{bias2mat=(meanmat}\OperatorTok{-}\KeywordTok{matrix}\NormalTok{(}\KeywordTok{rep}\NormalTok{(true_y,nord),}\DataTypeTok{byrow=}\OtherTok{TRUE}\NormalTok{,}\DataTypeTok{nrow=}\NormalTok{nord))}\OperatorTok{^}\DecValTok{2} \CommentTok{#here the truth is}

\NormalTok{df=}\KeywordTok{data.frame}\NormalTok{(}\KeywordTok{rep}\NormalTok{(x,}\DataTypeTok{each=}\NormalTok{nord),}\KeywordTok{rep}\NormalTok{(}\DecValTok{1}\OperatorTok{:}\NormalTok{nord,}\KeywordTok{length}\NormalTok{(x)),}\KeywordTok{c}\NormalTok{(bias2mat),}\KeywordTok{c}\NormalTok{(varmat),}\KeywordTok{rep}\NormalTok{(}\DecValTok{4}\NormalTok{,}\KeywordTok{prod}\NormalTok{(}\KeywordTok{dim}\NormalTok{(varmat)))) }\CommentTok{#irr is 2^2.}
\KeywordTok{colnames}\NormalTok{(df)=}\KeywordTok{c}\NormalTok{(}\StringTok{"x"}\NormalTok{,}\StringTok{"poly"}\NormalTok{,}\StringTok{"bias2"}\NormalTok{,}\StringTok{"variance"}\NormalTok{,}\StringTok{"irreducible error"}\NormalTok{) }\CommentTok{#suitable for plotting}
\NormalTok{df}\OperatorTok{$}\NormalTok{total=df}\OperatorTok{$}\NormalTok{bias2}\OperatorTok{+}\NormalTok{df}\OperatorTok{$}\NormalTok{variance}\OperatorTok{+}\NormalTok{df}\OperatorTok{$}\StringTok{`}\DataTypeTok{irreducible error}\StringTok{`}
\NormalTok{hdf=}\KeywordTok{melt}\NormalTok{(df,}\DataTypeTok{id=}\KeywordTok{c}\NormalTok{(}\StringTok{"x"}\NormalTok{,}\StringTok{"poly"}\NormalTok{))}
\NormalTok{hdf1=hdf[hdf}\OperatorTok{$}\NormalTok{poly}\OperatorTok{==}\DecValTok{1}\NormalTok{,]}
\NormalTok{hdf2=hdf[hdf}\OperatorTok{$}\NormalTok{poly}\OperatorTok{==}\DecValTok{2}\NormalTok{,]}
\NormalTok{hdf10=hdf[hdf}\OperatorTok{$}\NormalTok{poly}\OperatorTok{==}\DecValTok{10}\NormalTok{,]}
\NormalTok{hdf20=hdf[hdf}\OperatorTok{$}\NormalTok{poly}\OperatorTok{==}\DecValTok{20}\NormalTok{,]}

\NormalTok{p1=}\KeywordTok{ggplot}\NormalTok{(}\DataTypeTok{data=}\NormalTok{hdf1,}\KeywordTok{aes}\NormalTok{(}\DataTypeTok{x=}\NormalTok{x,}\DataTypeTok{y=}\NormalTok{value,}\DataTypeTok{colour=}\NormalTok{variable))}\OperatorTok{+}\KeywordTok{geom_line}\NormalTok{()}\OperatorTok{+}\KeywordTok{ggtitle}\NormalTok{(}\StringTok{"poly1"}\NormalTok{)}\OperatorTok{+}\KeywordTok{theme_minimal}\NormalTok{()}
\NormalTok{p2=}\KeywordTok{ggplot}\NormalTok{(}\DataTypeTok{data=}\NormalTok{hdf2,}\KeywordTok{aes}\NormalTok{(}\DataTypeTok{x=}\NormalTok{x,}\DataTypeTok{y=}\NormalTok{value,}\DataTypeTok{colour=}\NormalTok{variable))}\OperatorTok{+}\KeywordTok{geom_line}\NormalTok{()}\OperatorTok{+}\KeywordTok{ggtitle}\NormalTok{(}\StringTok{"poly2"}\NormalTok{)}\OperatorTok{+}\KeywordTok{theme_minimal}\NormalTok{()}
\NormalTok{p10=}\KeywordTok{ggplot}\NormalTok{(}\DataTypeTok{data=}\NormalTok{hdf10,}\KeywordTok{aes}\NormalTok{(}\DataTypeTok{x=}\NormalTok{x,}\DataTypeTok{y=}\NormalTok{value,}\DataTypeTok{colour=}\NormalTok{variable))}\OperatorTok{+}\KeywordTok{geom_line}\NormalTok{()}\OperatorTok{+}\KeywordTok{ggtitle}\NormalTok{(}\StringTok{"poly10"}\NormalTok{)}\OperatorTok{+}\KeywordTok{theme_minimal}\NormalTok{()}
\NormalTok{p20=}\KeywordTok{ggplot}\NormalTok{(}\DataTypeTok{data=}\NormalTok{hdf20,}\KeywordTok{aes}\NormalTok{(}\DataTypeTok{x=}\NormalTok{x,}\DataTypeTok{y=}\NormalTok{value,}\DataTypeTok{colour=}\NormalTok{variable))}\OperatorTok{+}\KeywordTok{geom_line}\NormalTok{()}\OperatorTok{+}\KeywordTok{ggtitle}\NormalTok{(}\StringTok{"poly20"}\NormalTok{)}\OperatorTok{+}\KeywordTok{theme_minimal}\NormalTok{()}
\KeywordTok{ggarrange}\NormalTok{(p1,p2,p10,p20)}
\end{Highlighting}
\end{Shaded}

\includegraphics{2RecEx-sol_files/figure-latex/unnamed-chunk-9-1.pdf}

We see that the variance (green) increases with the complexity of the
model. A linear model gives variance close to zero, while a polynomial
of degree 20 gives variance close to 1 (larger at the borders). A more
complex model is more flexible as it can turn up and down and change
direction fast. This leads to larger variance. (Look at the plot in 2a,
there is a larger variety of curves you can make when the degree is 20
compared to if the degree is 1.)

Further, we see that the bias is large for poly1, the linear model. The
linear model is quite rigid, so if the true underlying model is
non-linear, we typically get large deviations between the fitted line
and the training data. If we study the first plot, it seems like the
fitted line goes through the training data in \(x=-1\) and \(x=3\) as
the bias is close to zero here (this is confirmed by looking at the
upper left plot in 2a).

The polynomial models with degree larger than one lead to lower bias.
Recall that this is the training bias: The polynomial models will try to
pass through the training points if possible, and when the degree of the
polynomial is large they are able to do so because they have large
flexibility compared to a linear model.

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{hdfatxa=hdf[hdf}\OperatorTok{$}\NormalTok{x}\OperatorTok{==-}\DecValTok{1}\NormalTok{,]}
\NormalTok{hdfatxb=hdf[hdf}\OperatorTok{$}\NormalTok{x}\OperatorTok{==}\FloatTok{0.5}\NormalTok{,]}
\NormalTok{hdfatxc=hdf[hdf}\OperatorTok{$}\NormalTok{x}\OperatorTok{==}\DecValTok{2}\NormalTok{,]}
\NormalTok{hdfatxd=hdf[hdf}\OperatorTok{$}\NormalTok{x}\OperatorTok{==}\FloatTok{3.5}\NormalTok{,]}
\NormalTok{pa=}\KeywordTok{ggplot}\NormalTok{(}\DataTypeTok{data=}\NormalTok{hdfatxa,}\KeywordTok{aes}\NormalTok{(}\DataTypeTok{x=}\NormalTok{poly,}\DataTypeTok{y=}\NormalTok{value,}\DataTypeTok{colour=}\NormalTok{variable))}\OperatorTok{+}\KeywordTok{geom_line}\NormalTok{()}\OperatorTok{+}\KeywordTok{ggtitle}\NormalTok{(}\StringTok{"x0=-1"}\NormalTok{)}\OperatorTok{+}\KeywordTok{theme_minimal}\NormalTok{()}
\NormalTok{pb=}\KeywordTok{ggplot}\NormalTok{(}\DataTypeTok{data=}\NormalTok{hdfatxb,}\KeywordTok{aes}\NormalTok{(}\DataTypeTok{x=}\NormalTok{poly,}\DataTypeTok{y=}\NormalTok{value,}\DataTypeTok{colour=}\NormalTok{variable))}\OperatorTok{+}\KeywordTok{geom_line}\NormalTok{()}\OperatorTok{+}\KeywordTok{ggtitle}\NormalTok{(}\StringTok{"x0=0.5"}\NormalTok{)}\OperatorTok{+}\KeywordTok{theme_minimal}\NormalTok{()}
\NormalTok{pc=}\KeywordTok{ggplot}\NormalTok{(}\DataTypeTok{data=}\NormalTok{hdfatxc,}\KeywordTok{aes}\NormalTok{(}\DataTypeTok{x=}\NormalTok{poly,}\DataTypeTok{y=}\NormalTok{value,}\DataTypeTok{colour=}\NormalTok{variable))}\OperatorTok{+}\KeywordTok{geom_line}\NormalTok{()}\OperatorTok{+}\KeywordTok{ggtitle}\NormalTok{(}\StringTok{"x0=2"}\NormalTok{)}\OperatorTok{+}\KeywordTok{theme_minimal}\NormalTok{()}
\NormalTok{pd=}\KeywordTok{ggplot}\NormalTok{(}\DataTypeTok{data=}\NormalTok{hdfatxd,}\KeywordTok{aes}\NormalTok{(}\DataTypeTok{x=}\NormalTok{poly,}\DataTypeTok{y=}\NormalTok{value,}\DataTypeTok{colour=}\NormalTok{variable))}\OperatorTok{+}\KeywordTok{geom_line}\NormalTok{()}\OperatorTok{+}\KeywordTok{ggtitle}\NormalTok{(}\StringTok{"x0=3.5"}\NormalTok{)}\OperatorTok{+}\KeywordTok{theme_minimal}\NormalTok{()}
\KeywordTok{ggarrange}\NormalTok{(pa,pb,pc,pd)}
\end{Highlighting}
\end{Shaded}

\includegraphics{2RecEx-sol_files/figure-latex/2bbiasvariance2-1.pdf}

Compare to Figures in 2.12 on page 36 in ISL (our textbook).

\subsection{d)}\label{d}

To change \(f(x)\), replace

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{truefunc=}\ControlFlowTok{function}\NormalTok{(x) }\KeywordTok{return}\NormalTok{(x}\OperatorTok{^}\DecValTok{2}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

by for example

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{truefunc=}\ControlFlowTok{function}\NormalTok{(x) }\KeywordTok{return}\NormalTok{(x}\OperatorTok{^}\DecValTok{4}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

or

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{truefunc=}\ControlFlowTok{function}\NormalTok{(x) }\KeywordTok{return}\NormalTok{(}\KeywordTok{exp}\NormalTok{(}\DecValTok{2}\OperatorTok{*}\NormalTok{x))}
\end{Highlighting}
\end{Shaded}

and rerun the code. Study the results.

If you want to set the variance to 1 for example, set \(sd=1\) in these
parts of the code in 2a and 2b:

\begin{verbatim}
rnorm(length(x)*M, mean=0, sd=1)
\end{verbatim}

Also change the following part in 2c:

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{df=}\KeywordTok{data.frame}\NormalTok{(}\KeywordTok{rep}\NormalTok{(x,}\DataTypeTok{each=}\NormalTok{nord),}\KeywordTok{rep}\NormalTok{(}\DecValTok{1}\OperatorTok{:}\NormalTok{nord,}\KeywordTok{length}\NormalTok{(x)),}\KeywordTok{c}\NormalTok{(bias2mat),}\KeywordTok{c}\NormalTok{(varmat),}\KeywordTok{rep}\NormalTok{(}\DecValTok{1}\NormalTok{,}\KeywordTok{prod}\NormalTok{(}\KeywordTok{dim}\NormalTok{(varmat)))) }\CommentTok{#irr is 1^2.}
\end{Highlighting}
\end{Shaded}

to get correct plots of the irreducible error. Here,
\(rep(4,prod(dim(varmat)))\) is replaced by
\(rep(1,prod(dim(varmat)))\).

\section{ R packages}\label{r-packages}

If you want to look at the .Rmd file and \texttt{knit} it, you need to
first install the following packages (only once).

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{install.packages}\NormalTok{(}\StringTok{"ggplot2"}\NormalTok{)}
\KeywordTok{install.packages}\NormalTok{(}\StringTok{"gamlss.data"}\NormalTok{)}
\KeywordTok{install.packages}\NormalTok{(}\StringTok{"tidyverse"}\NormalTok{)}
\KeywordTok{install.packages}\NormalTok{(}\StringTok{"GGally"}\NormalTok{)}
\KeywordTok{install.packages}\NormalTok{(}\StringTok{"Matrix"}\NormalTok{)}
\KeywordTok{install.packages}\NormalTok{(}\StringTok{"ggpubr"}\NormalTok{)}
\end{Highlighting}
\end{Shaded}


\end{document}
