<h1>Table of Contents<span class="tocSkip"></span></h1>
<div class="toc"><ul class="toc-item"><li><span><a href="#Regression" data-toc-modified-id="Regression-1"><span class="toc-item-num">1&nbsp;&nbsp;</span>Regression</a></span><ul class="toc-item"><li><ul class="toc-item"><li><span><a href="#Bibliography:-Chapter-8-Ivezic+" data-toc-modified-id="Bibliography:-Chapter-8-Ivezic+-1.0.1"><span class="toc-item-num">1.0.1&nbsp;&nbsp;</span>Bibliography: Chapter 8 Ivezic+</a></span></li></ul></li></ul></li><li><span><a href="#Regression:" data-toc-modified-id="Regression:-2"><span class="toc-item-num">2&nbsp;&nbsp;</span>Regression:</a></span></li><li><span><a href="#Metrics:" data-toc-modified-id="Metrics:-3"><span class="toc-item-num">3&nbsp;&nbsp;</span>Metrics:</a></span><ul class="toc-item"><li><span><a href="#Mean-squared-error-(MSE)" data-toc-modified-id="Mean-squared-error-(MSE)-3.1"><span class="toc-item-num">3.1&nbsp;&nbsp;</span>Mean squared error (MSE)</a></span></li><li><span><a href="#$MSE-\equiv-\frac{1}{n}-\sum_i-(y_i---\hat{y}_i)^2$" data-toc-modified-id="$MSE-\equiv-\frac{1}{n}-\sum_i-(y_i---\hat{y}_i)^2$-3.2"><span class="toc-item-num">3.2&nbsp;&nbsp;</span>$MSE \equiv \frac{1}{n} \sum_i (y_i - \hat{y}_i)^2$</a></span></li><li><span><a href="#Coefficient-of-determination-(R$^2$)" data-toc-modified-id="Coefficient-of-determination-(R$^2$)-3.3"><span class="toc-item-num">3.3&nbsp;&nbsp;</span>Coefficient of determination (R$^2$)</a></span></li><li><span><a href="#$R^2-\equiv-1----\frac{SS_{res}}{SS_{tot}}$," data-toc-modified-id="$R^2-\equiv-1----\frac{SS_{res}}{SS_{tot}}$,-3.4"><span class="toc-item-num">3.4&nbsp;&nbsp;</span>$R^2 \equiv 1  - \frac{SS_{res}}{SS_{tot}}$,</a></span><ul class="toc-item"><li><span><a href="#$-SS_{res}-=-\sum_i-(y_i---\hat{y}_i)^2$" data-toc-modified-id="$-SS_{res}-=-\sum_i-(y_i---\hat{y}_i)^2$-3.4.1"><span class="toc-item-num">3.4.1&nbsp;&nbsp;</span>$ SS_{res} = \sum_i (y_i - \hat{y}_i)^2$</a></span></li><li><span><a href="#$-SS_{tot}-=-\sum_i-(y_i---\bar{y})^2$" data-toc-modified-id="$-SS_{tot}-=-\sum_i-(y_i---\bar{y})^2$-3.4.2"><span class="toc-item-num">3.4.2&nbsp;&nbsp;</span>$ SS_{tot} = \sum_i (y_i - \bar{y})^2$</a></span></li></ul></li><li><span><a href="#Variance-and-bias" data-toc-modified-id="Variance-and-bias-3.5"><span class="toc-item-num">3.5&nbsp;&nbsp;</span>Variance and bias</a></span></li><li><span><a href="#$Var_{\hat{\theta}}-=-E_{\hat{\theta}}\biggl[-\bigl(\hat{\theta}---E_{\hat{\theta}}[\hat{\theta}]\bigr)^2\biggr]$" data-toc-modified-id="$Var_{\hat{\theta}}-=-E_{\hat{\theta}}\biggl[-\bigl(\hat{\theta}---E_{\hat{\theta}}[\hat{\theta}]\bigr)^2\biggr]$-3.6"><span class="toc-item-num">3.6&nbsp;&nbsp;</span>$Var_{\hat{\theta}} = E_{\hat{\theta}}\biggl[ \bigl(\hat{\theta} - E_{\hat{\theta}}[\hat{\theta}]\bigr)^2\biggr]$</a></span></li><li><span><a href="#$Bias_{\hat{\theta}}-=-E_{\hat{\theta}}---\theta$" data-toc-modified-id="$Bias_{\hat{\theta}}-=-E_{\hat{\theta}}---\theta$-3.7"><span class="toc-item-num">3.7&nbsp;&nbsp;</span>$Bias_{\hat{\theta}} = E_{\hat{\theta}} - \theta$</a></span></li><li><span><a href="#Relation-between-MSE,-variance-and-bias" data-toc-modified-id="Relation-between-MSE,-variance-and-bias-3.8"><span class="toc-item-num">3.8&nbsp;&nbsp;</span>Relation between MSE, variance and bias</a></span></li></ul></li><li><span><a href="#Linear-regression" data-toc-modified-id="Linear-regression-4"><span class="toc-item-num">4&nbsp;&nbsp;</span>Linear regression</a></span><ul class="toc-item"><li><span><a href="#Least-squares-regression-for-linear-regression" data-toc-modified-id="Least-squares-regression-for-linear-regression-4.1"><span class="toc-item-num">4.1&nbsp;&nbsp;</span>Least squares regression for linear regression</a></span></li><li><span><a href="#1.-Ordinary-least-squares" data-toc-modified-id="1.-Ordinary-least-squares-4.2"><span class="toc-item-num">4.2&nbsp;&nbsp;</span>1. Ordinary least squares</a></span><ul class="toc-item"><li><span><a href="#$\min-RSS-=-\min-\sum\limits_{i=1}^{n}-(y_i---\theta_0---\theta_1-x_i)^2$" data-toc-modified-id="$\min-RSS-=-\min-\sum\limits_{i=1}^{n}-(y_i---\theta_0---\theta_1-x_i)^2$-4.2.1"><span class="toc-item-num">4.2.1&nbsp;&nbsp;</span>$\min RSS = \min \sum\limits_{i=1}^{n} (y_i - \theta_0 - \theta_1 x_i)^2$</a></span></li></ul></li><li><span><a href="#$\hat{\theta}_{1,-OLS}-=-\frac{\sum\limits_{i=1}^n-(x_i---\bar{x})-(y_i---\bar{y})}{\sum\limits_{i=1}^n-(x_i---\bar{x})^2}-=-\frac{S_{xy}}{S_{xx}}$" data-toc-modified-id="$\hat{\theta}_{1,-OLS}-=-\frac{\sum\limits_{i=1}^n-(x_i---\bar{x})-(y_i---\bar{y})}{\sum\limits_{i=1}^n-(x_i---\bar{x})^2}-=-\frac{S_{xy}}{S_{xx}}$-4.3"><span class="toc-item-num">4.3&nbsp;&nbsp;</span>$\hat{\theta}_{1, OLS} = \frac{\sum\limits_{i=1}^n (x_i - \bar{x}) (y_i - \bar{y})}{\sum\limits_{i=1}^n (x_i - \bar{x})^2} = \frac{S_{xy}}{S_{xx}}$</a></span></li><li><span><a href="#$\hat{\theta}_{0,-OLS}-=-\bar{y}---\hat{\theta}_{1,OLS}-~\bar{x}$" data-toc-modified-id="$\hat{\theta}_{0,-OLS}-=-\bar{y}---\hat{\theta}_{1,OLS}-~\bar{x}$-4.4"><span class="toc-item-num">4.4&nbsp;&nbsp;</span>$\hat{\theta}_{0, OLS} = \bar{y} - \hat{\theta}_{1,OLS} ~\bar{x}$</a></span></li><li><span><a href="#$\hat{\sigma}^2-=-S^2-=-\frac{RSS}{n---2}$" data-toc-modified-id="$\hat{\sigma}^2-=-S^2-=-\frac{RSS}{n---2}$-4.5"><span class="toc-item-num">4.5&nbsp;&nbsp;</span>$\hat{\sigma}^2 = S^2 = \frac{RSS}{n - 2}$</a></span></li><li><span><a href="#$\hat{\theta}_{1,-OLS}-\sim-N\biggl(\theta_1,-\frac{\sigma^2}{S_{XX}}\biggr)$" data-toc-modified-id="$\hat{\theta}_{1,-OLS}-\sim-N\biggl(\theta_1,-\frac{\sigma^2}{S_{XX}}\biggr)$-4.6"><span class="toc-item-num">4.6&nbsp;&nbsp;</span>$\hat{\theta}_{1, OLS} \sim N\biggl(\theta_1, \frac{\sigma^2}{S_{XX}}\biggr)$</a></span></li><li><span><a href="#$\hat{\theta}_{0,-OLS}-\sim-N\biggl(\theta_0,-\sigma^2-\biggl(\frac{1}{n}-+-\frac{\bar{x}^2}{S_{XX}}\biggr)\biggr)$" data-toc-modified-id="$\hat{\theta}_{0,-OLS}-\sim-N\biggl(\theta_0,-\sigma^2-\biggl(\frac{1}{n}-+-\frac{\bar{x}^2}{S_{XX}}\biggr)\biggr)$-4.7"><span class="toc-item-num">4.7&nbsp;&nbsp;</span>$\hat{\theta}_{0, OLS} \sim N\biggl(\theta_0, \sigma^2 \biggl(\frac{1}{n} + \frac{\bar{x}^2}{S_{XX}}\biggr)\biggr)$</a></span></li><li><span><a href="#$\hat{\sigma}^2-\sim-\biggl(\frac{\sigma^2}{n---2}-\biggr)-\chi_{n-2}^{2}$" data-toc-modified-id="$\hat{\sigma}^2-\sim-\biggl(\frac{\sigma^2}{n---2}-\biggr)-\chi_{n-2}^{2}$-4.8"><span class="toc-item-num">4.8&nbsp;&nbsp;</span>$\hat{\sigma}^2 \sim \biggl(\frac{\sigma^2}{n - 2} \biggr) \chi_{n-2}^{2}$</a></span></li><li><span><a href="#$-T-=-\frac{\hat{\theta_1}---\theta_1}{SE}$" data-toc-modified-id="$-T-=-\frac{\hat{\theta_1}---\theta_1}{SE}$-4.9"><span class="toc-item-num">4.9&nbsp;&nbsp;</span>$ T = \frac{\hat{\theta_1} - \theta_1}{SE}$</a></span></li><li><span><a href="#$y(x)-=-\hat{\theta}_0-+-\hat{\theta}_1-x-\pm-t_{\alpha/2,-n---2}-~-S-~-\sqrt{1-+-\frac{1}{n}-\frac{(x---\bar{x})^2}{S_{xx}}}$" data-toc-modified-id="$y(x)-=-\hat{\theta}_0-+-\hat{\theta}_1-x-\pm-t_{\alpha/2,-n---2}-~-S-~-\sqrt{1-+-\frac{1}{n}-\frac{(x---\bar{x})^2}{S_{xx}}}$-4.10"><span class="toc-item-num">4.10&nbsp;&nbsp;</span>$y(x) = \hat{\theta}_0 + \hat{\theta}_1 x \pm t_{\alpha/2, n - 2} ~ S ~ \sqrt{1 + \frac{1}{n} \frac{(x - \bar{x})^2}{S_{xx}}}$</a></span></li><li><span><a href="#2.-Maximum-likelihood-estimator-(MLE)" data-toc-modified-id="2.-Maximum-likelihood-estimator-(MLE)-4.11"><span class="toc-item-num">4.11&nbsp;&nbsp;</span>2. Maximum likelihood estimator (MLE)</a></span></li><li><span><a href="#$p(\lbrace-y_i-\rbrace-|-\lbrace-x_i-\rbrace,-\theta)-=-\prod\limits_{i-=-1}^{N}-\frac{1}{\sqrt{2-\pi}-\sigma_i}-\exp-\biggl(-\frac{--(y_i---(\theta_0-+-\theta_1-x_i))^2}{2-\sigma_i^2}-\biggr)$" data-toc-modified-id="$p(\lbrace-y_i-\rbrace-|-\lbrace-x_i-\rbrace,-\theta)-=-\prod\limits_{i-=-1}^{N}-\frac{1}{\sqrt{2-\pi}-\sigma_i}-\exp-\biggl(-\frac{--(y_i---(\theta_0-+-\theta_1-x_i))^2}{2-\sigma_i^2}-\biggr)$-4.12"><span class="toc-item-num">4.12&nbsp;&nbsp;</span>$p(\lbrace y_i \rbrace | \lbrace x_i \rbrace, \theta) = \prod\limits_{i = 1}^{N} \frac{1}{\sqrt{2 \pi} \sigma_i} \exp \biggl( \frac{- (y_i - (\theta_0 + \theta_1 x_i))^2}{2 \sigma_i^2} \biggr)$</a></span></li><li><span><a href="#$\log(L)-=-\log-p(\lbrace-y_i-\rbrace-|-\lbrace-x_i-\rbrace,-\theta)-=-\sum\limits_{i-=-1}^{N}--\bigl(-\frac{--(y_i---(\theta_0-+-\theta_1-x_i))^2}{2-\sigma_i^2}-\bigr)$" data-toc-modified-id="$\log(L)-=-\log-p(\lbrace-y_i-\rbrace-|-\lbrace-x_i-\rbrace,-\theta)-=-\sum\limits_{i-=-1}^{N}--\bigl(-\frac{--(y_i---(\theta_0-+-\theta_1-x_i))^2}{2-\sigma_i^2}-\bigr)$-4.13"><span class="toc-item-num">4.13&nbsp;&nbsp;</span>$\log(L) = \log p(\lbrace y_i \rbrace | \lbrace x_i \rbrace, \theta) = \sum\limits_{i = 1}^{N}  \bigl( \frac{- (y_i - (\theta_0 + \theta_1 x_i))^2}{2 \sigma_i^2} \bigr)$</a></span></li><li><span><a href="#$P(\theta-|-\lbrace-y_i-\rbrace,-\lbrace-x_i-\rbrace)-=-P(\lbrace-y_i-\rbrace-|-\lbrace-x_i-\rbrace,-\theta)-~-\frac{p(\theta)}{p(y_i)}$" data-toc-modified-id="$P(\theta-|-\lbrace-y_i-\rbrace,-\lbrace-x_i-\rbrace)-=-P(\lbrace-y_i-\rbrace-|-\lbrace-x_i-\rbrace,-\theta)-~-\frac{p(\theta)}{p(y_i)}$-4.14"><span class="toc-item-num">4.14&nbsp;&nbsp;</span>$P(\theta | \lbrace y_i \rbrace, \lbrace x_i \rbrace) = P(\lbrace y_i \rbrace | \lbrace x_i \rbrace, \theta) ~ \frac{p(\theta)}{p(y_i)}$</a></span></li><li><span><a href="#Online-estimation" data-toc-modified-id="Online-estimation-4.15"><span class="toc-item-num">4.15&nbsp;&nbsp;</span>Online estimation</a></span><ul class="toc-item"><li><span><a href="#Diabetes-dataset-example" data-toc-modified-id="Diabetes-dataset-example-4.15.1"><span class="toc-item-num">4.15.1&nbsp;&nbsp;</span>Diabetes dataset example</a></span></li></ul></li><li><span><a href="#3.-Least-squares-in-matrix-form" data-toc-modified-id="3.-Least-squares-in-matrix-form-4.16"><span class="toc-item-num">4.16&nbsp;&nbsp;</span>3. Least squares in matrix form</a></span></li><li><span><a href="#$y-=-{\bf-M}-~-\theta$" data-toc-modified-id="$y-=-{\bf-M}-~-\theta$-4.17"><span class="toc-item-num">4.17&nbsp;&nbsp;</span>$y = {\bf M} ~ \theta$</a></span></li><li><span><a href="#$y-=-\begin{bmatrix}-y_0-\\-y_1-\\-.-\\-y_{N-1}-\end{bmatrix}$" data-toc-modified-id="$y-=-\begin{bmatrix}-y_0-\\-y_1-\\-.-\\-y_{N-1}-\end{bmatrix}$-4.18"><span class="toc-item-num">4.18&nbsp;&nbsp;</span>$y = \begin{bmatrix} y_0 \\ y_1 \\ . \\ y_{N-1} \end{bmatrix}$</a></span></li><li><span><a href="#$\theta-=-\begin{bmatrix}-\theta_0-\\-\theta_1-\end{bmatrix}$" data-toc-modified-id="$\theta-=-\begin{bmatrix}-\theta_0-\\-\theta_1-\end{bmatrix}$-4.19"><span class="toc-item-num">4.19&nbsp;&nbsp;</span>$\theta = \begin{bmatrix} \theta_0 \\ \theta_1 \end{bmatrix}$</a></span></li><li><span><a href="#${\bf-M}-=-\begin{bmatrix}-1-&amp;-x_0-\\-1-&amp;-x_1-\\-...-&amp;-...-\\-1-&amp;-x_{N-1}-\\-\end{bmatrix}$" data-toc-modified-id="${\bf-M}-=-\begin{bmatrix}-1-&amp;-x_0-\\-1-&amp;-x_1-\\-...-&amp;-...-\\-1-&amp;-x_{N-1}-\\-\end{bmatrix}$-4.20"><span class="toc-item-num">4.20&nbsp;&nbsp;</span>${\bf M} = \begin{bmatrix} 1 &amp; x_0 \\ 1 &amp; x_1 \\ ... &amp; ... \\ 1 &amp; x_{N-1} \\ \end{bmatrix}$</a></span></li><li><span><a href="#$-C-=-\begin{bmatrix}-\sigma_0^2-&amp;-0-&amp;-.-&amp;-0-\\-0-&amp;-\sigma_1^2-&amp;-.-&amp;-0-\\-.-&amp;-.-&amp;-.-&amp;-.-\\-0-&amp;-0-&amp;-.-&amp;-\sigma_{N---1}^2\\-\end{bmatrix}$" data-toc-modified-id="$-C-=-\begin{bmatrix}-\sigma_0^2-&amp;-0-&amp;-.-&amp;-0-\\-0-&amp;-\sigma_1^2-&amp;-.-&amp;-0-\\-.-&amp;-.-&amp;-.-&amp;-.-\\-0-&amp;-0-&amp;-.-&amp;-\sigma_{N---1}^2\\-\end{bmatrix}$-4.21"><span class="toc-item-num">4.21&nbsp;&nbsp;</span>$ C = \begin{bmatrix} \sigma_0^2 &amp; 0 &amp; . &amp; 0 \\ 0 &amp; \sigma_1^2 &amp; . &amp; 0 \\ . &amp; . &amp; . &amp; . \\ 0 &amp; 0 &amp; . &amp; \sigma_{N - 1}^2\\ \end{bmatrix}$</a></span></li><li><span><a href="#${\bf-\hat{\theta}}-=-({\bf-M}^T-C^{-1}-{\bf-M})^{-1}-({\bf-M}^T-C^{-1}-{\bf-y})$" data-toc-modified-id="${\bf-\hat{\theta}}-=-({\bf-M}^T-C^{-1}-{\bf-M})^{-1}-({\bf-M}^T-C^{-1}-{\bf-y})$-4.22"><span class="toc-item-num">4.22&nbsp;&nbsp;</span>${\bf \hat{\theta}} = ({\bf M}^T C^{-1} {\bf M})^{-1} ({\bf M}^T C^{-1} {\bf y})$</a></span></li><li><span><a href="#$({\bf-y}---{\bf-M}-~-\theta)^{-1}-C^{-1}-({\bf-y}----{\bf-M}-~-\theta)$" data-toc-modified-id="$({\bf-y}---{\bf-M}-~-\theta)^{-1}-C^{-1}-({\bf-y}----{\bf-M}-~-\theta)$-4.23"><span class="toc-item-num">4.23&nbsp;&nbsp;</span>$({\bf y} - {\bf M} ~ \theta)^{-1} C^{-1} ({\bf y} -  {\bf M} ~ \theta)$</a></span></li><li><span><a href="#$\Sigma_\theta-=-\begin{bmatrix}-\sigma_{\theta_0}^2-&amp;-\dots-&amp;-\sigma_{\theta_0-\theta_i}-&amp;-\dots-\\-\vdots-&amp;-\ddots-&amp;-\vdots-&amp;-\vdots-\\-\sigma_{\theta_0-\theta_i}-&amp;-\dots-&amp;-\sigma_{\theta_i}^2-\\--\vdots-&amp;-\vdots-&amp;-\vdots-&amp;-\ddots-\\-\end{bmatrix}-=-[{\bf-M}^T-C^{-1}-{\bf-M}]^{-1}$" data-toc-modified-id="$\Sigma_\theta-=-\begin{bmatrix}-\sigma_{\theta_0}^2-&amp;-\dots-&amp;-\sigma_{\theta_0-\theta_i}-&amp;-\dots-\\-\vdots-&amp;-\ddots-&amp;-\vdots-&amp;-\vdots-\\-\sigma_{\theta_0-\theta_i}-&amp;-\dots-&amp;-\sigma_{\theta_i}^2-\\--\vdots-&amp;-\vdots-&amp;-\vdots-&amp;-\ddots-\\-\end{bmatrix}-=-[{\bf-M}^T-C^{-1}-{\bf-M}]^{-1}$-4.24"><span class="toc-item-num">4.24&nbsp;&nbsp;</span>$\Sigma_\theta = \begin{bmatrix} \sigma_{\theta_0}^2 &amp; \dots &amp; \sigma_{\theta_0 \theta_i} &amp; \dots \\ \vdots &amp; \ddots &amp; \vdots &amp; \vdots \\ \sigma_{\theta_0 \theta_i} &amp; \dots &amp; \sigma_{\theta_i}^2 \\  \vdots &amp; \vdots &amp; \vdots &amp; \ddots \\ \end{bmatrix} = [{\bf M}^T C^{-1} {\bf M}]^{-1}$</a></span></li><li><span><a href="#4.-Polynomial-regression" data-toc-modified-id="4.-Polynomial-regression-4.25"><span class="toc-item-num">4.25&nbsp;&nbsp;</span>4. Polynomial regression</a></span></li><li><span><a href="#$y_i-=-\theta_0-+-\theta_1-x_i-+-\theta_2-x_i^2-+-\theta_3-x_i^3$" data-toc-modified-id="$y_i-=-\theta_0-+-\theta_1-x_i-+-\theta_2-x_i^2-+-\theta_3-x_i^3$-4.26"><span class="toc-item-num">4.26&nbsp;&nbsp;</span>$y_i = \theta_0 + \theta_1 x_i + \theta_2 x_i^2 + \theta_3 x_i^3$</a></span></li><li><span><a href="#$-M-=-\begin{bmatrix}
1-&amp;-x_0-&amp;-x_0^2-&amp;-x_0^3-\\
1-&amp;-x_1-&amp;-x_1^2-&amp;-x_1^3-\\
.-&amp;-.-&amp;-.-&amp;-.-\\
1-&amp;-x_N-&amp;-x_N^2-&amp;-x_N^3-\\
\end{bmatrix}$" data-toc-modified-id="$-M-=-\begin{bmatrix}
1-&amp;-x_0-&amp;-x_0^2-&amp;-x_0^3-\\
1-&amp;-x_1-&amp;-x_1^2-&amp;-x_1^3-\\
.-&amp;-.-&amp;-.-&amp;-.-\\
1-&amp;-x_N-&amp;-x_N^2-&amp;-x_N^3-\\
\end{bmatrix}$-4.27"><span class="toc-item-num">4.27&nbsp;&nbsp;</span>$ M = \begin{bmatrix}
1 &amp; x_0 &amp; x_0^2 &amp; x_0^3 \\
1 &amp; x_1 &amp; x_1^2 &amp; x_1^3 \\
. &amp; . &amp; . &amp; . \\
1 &amp; x_N &amp; x_N^2 &amp; x_N^3 \\
\end{bmatrix}$</a></span></li></ul></li><li><span><a href="#Regularization" data-toc-modified-id="Regularization-5"><span class="toc-item-num">5&nbsp;&nbsp;</span>Regularization</a></span><ul class="toc-item"><li><span><a href="#1.-Ridge-regression" data-toc-modified-id="1.-Ridge-regression-5.1"><span class="toc-item-num">5.1&nbsp;&nbsp;</span>1. Ridge regression</a></span></li><li><span><a href="#$({\bf-y}---M~-\theta)^{-1}-C^{-1}({\bf-y}---M~-\theta)-+-\lambda~-|\theta^T-\theta|$" data-toc-modified-id="$({\bf-y}---M~-\theta)^{-1}-C^{-1}({\bf-y}---M~-\theta)-+-\lambda~-|\theta^T-\theta|$-5.2"><span class="toc-item-num">5.2&nbsp;&nbsp;</span>$({\bf y} - M~ \theta)^{-1} C^{-1}({\bf y} - M~ \theta) + \lambda~ |\theta^T \theta|$</a></span></li><li><span><a href="#${\bf-\hat{\theta}}-=-({\bf-M}^T-C^{-1}-{\bf-M}-+-\lambda-I)^{-1}-({\bf-M}^T-C^{-1}-{\bf-y})$" data-toc-modified-id="${\bf-\hat{\theta}}-=-({\bf-M}^T-C^{-1}-{\bf-M}-+-\lambda-I)^{-1}-({\bf-M}^T-C^{-1}-{\bf-y})$-5.3"><span class="toc-item-num">5.3&nbsp;&nbsp;</span>${\bf \hat{\theta}} = ({\bf M}^T C^{-1} {\bf M} + \lambda I)^{-1} ({\bf M}^T C^{-1} {\bf y})$</a></span><ul class="toc-item"><li><span><a href="#Bayes-equivalent" data-toc-modified-id="Bayes-equivalent-5.3.1"><span class="toc-item-num">5.3.1&nbsp;&nbsp;</span>Bayes equivalent</a></span></li></ul></li><li><span><a href="#$-p(\theta)-\propto-\exp-\biggl(-\frac{--\lambda-(\theta^T-\theta)}{2}-\biggr)$" data-toc-modified-id="$-p(\theta)-\propto-\exp-\biggl(-\frac{--\lambda-(\theta^T-\theta)}{2}-\biggr)$-5.4"><span class="toc-item-num">5.4&nbsp;&nbsp;</span>$ p(\theta) \propto \exp \biggl( \frac{- \lambda (\theta^T \theta)}{2} \biggr)$</a></span></li><li><span><a href="#$|{\bf-\theta}|^2-<-s$" data-toc-modified-id="$|{\bf-\theta}|^2-<-s$-5.5"><span class="toc-item-num">5.5&nbsp;&nbsp;</span>$|{\bf \theta}|^2 &lt; s$</a></span></li><li><span><a href="#2.-LASSO-regression" data-toc-modified-id="2.-LASSO-regression-5.6"><span class="toc-item-num">5.6&nbsp;&nbsp;</span>2. LASSO regression</a></span></li><li><span><a href="#$({\bf-y}---\theta-{\bf-X})^{-1}({\bf-y}---\theta-{\bf-X})-+-\lambda-|\theta|$" data-toc-modified-id="$({\bf-y}---\theta-{\bf-X})^{-1}({\bf-y}---\theta-{\bf-X})-+-\lambda-|\theta|$-5.7"><span class="toc-item-num">5.7&nbsp;&nbsp;</span>$({\bf y} - \theta {\bf X})^{-1}({\bf y} - \theta {\bf X}) + \lambda |\theta|$</a></span></li><li><span><a href="#$|{\bf-\theta}|-<-s$" data-toc-modified-id="$|{\bf-\theta}|-<-s$-5.8"><span class="toc-item-num">5.8&nbsp;&nbsp;</span>$|{\bf \theta}| &lt; s$</a></span><ul class="toc-item"><li><span><a href="#Understanding-Ridge-and-LASSO-regression" data-toc-modified-id="Understanding-Ridge-and-LASSO-regression-5.8.1"><span class="toc-item-num">5.8.1&nbsp;&nbsp;</span>Understanding Ridge and LASSO regression</a></span></li></ul></li></ul></li><li><span><a href="#Comparison-between-ordinary,-Ridge-and-LASSO-regression" data-toc-modified-id="Comparison-between-ordinary,-Ridge-and-LASSO-regression-6"><span class="toc-item-num">6&nbsp;&nbsp;</span>Comparison between ordinary, Ridge and LASSO regression</a></span><ul class="toc-item"><li><span><a href="#Choosing-the-regularization-parameter-$\lambda$" data-toc-modified-id="Choosing-the-regularization-parameter-$\lambda$-6.1"><span class="toc-item-num">6.1&nbsp;&nbsp;</span>Choosing the regularization parameter $\lambda$</a></span></li><li><span><a href="#$-Err(\lambda)-=-k^{1}-\sum_k-N_k^{-1}-\sum\limits_{i}^{N_k}-\frac{[y_i---f(x_i-|-\theta)]^2}{\sigma_i^2}$" data-toc-modified-id="$-Err(\lambda)-=-k^{1}-\sum_k-N_k^{-1}-\sum\limits_{i}^{N_k}-\frac{[y_i---f(x_i-|-\theta)]^2}{\sigma_i^2}$-6.2"><span class="toc-item-num">6.2&nbsp;&nbsp;</span>$ Err(\lambda) = k^{1} \sum_k N_k^{-1} \sum\limits_{i}^{N_k} \frac{[y_i - f(x_i | \theta)]^2}{\sigma_i^2}$</a></span></li></ul></li><li><span><a href="#Outliers" data-toc-modified-id="Outliers-7"><span class="toc-item-num">7&nbsp;&nbsp;</span>Outliers</a></span><ul class="toc-item"><li><span><a href="#Huber-loss" data-toc-modified-id="Huber-loss-7.1"><span class="toc-item-num">7.1&nbsp;&nbsp;</span>Huber loss</a></span><ul class="toc-item"><li><span><a href="#$-\sum\limits_{i=1}^N-e(y_i|y)$" data-toc-modified-id="$-\sum\limits_{i=1}^N-e(y_i|y)$-7.1.1"><span class="toc-item-num">7.1.1&nbsp;&nbsp;</span>$ \sum\limits_{i=1}^N e(y_i|y)$</a></span></li><li><span><a href="#$\begin{equation}
--\phi(t)=\left\{
--\begin{array}{@{}ll@{}}
----\frac{1}{2}-t^2,-&amp;-\text{if}\-|t|\le-c-\\
----c|t|---\frac{1}{2}-c^2,-&amp;-\text{otherwise}
--\end{array}\right.
\end{equation}$" data-toc-modified-id="$\begin{equation}
--\phi(t)=\left\{
--\begin{array}{@{}ll@{}}
----\frac{1}{2}-t^2,-&amp;-\text{if}\-|t|\le-c-\\
----c|t|---\frac{1}{2}-c^2,-&amp;-\text{otherwise}
--\end{array}\right.
\end{equation}$-7.1.2"><span class="toc-item-num">7.1.2&nbsp;&nbsp;</span>$\begin{equation}
  \phi(t)=\left\{
  \begin{array}{@{}ll@{}}
    \frac{1}{2} t^2, &amp; \text{if}\ |t|\le c \\
    c|t| - \frac{1}{2} c^2, &amp; \text{otherwise}
  \end{array}\right.
\end{equation}$</a></span></li></ul></li><li><span><a href="#Bayesian-Outlier-Methods" data-toc-modified-id="Bayesian-Outlier-Methods-7.2"><span class="toc-item-num">7.2&nbsp;&nbsp;</span>Bayesian Outlier Methods</a></span></li></ul></li><li><span><a href="#Cross-Validation-Testing" data-toc-modified-id="Cross-Validation-Testing-8"><span class="toc-item-num">8&nbsp;&nbsp;</span>Cross-Validation Testing</a></span></li></ul></div>


```python
%config InlineBackend.figure_format = 'retina'
%matplotlib inline

import numpy as np
import matplotlib.pyplot as plt

plt.style.use(['ggplot', 'assets/class.mplstyle'])

red = '#E24A33'
blue = '#348ABD'
purple = '#988ED5'
gray = '#777777'
yellow = '#FBC15E'
green = '#8EBA42'
pink = '#FFB5B8'
```

# Regression

### Bibliography: Chapter 8 Ivezic+

# Regression:

- Finding the relation between a dependent variable $y$ and a set of independent variables $\bf x$. It describes the expectation value of $y$ given $\bf x$, $E[y|{\bf x}]$.

- The dependent variable $y$ is also called **response variable**.

- The problem is to find the expectation of $y$ given $\bf x$:

$E[y|{\bf x}] = f({\bf x}, \theta) + \epsilon$

where $\epsilon$ is a random error.

Note that we estimate $E[y | {\bf x}]$, the conditional expectation of $y$ given $\bf x$, rather than $p(y, {\bf x})$, the joint distribution of $y$ and $\bf x$, which is a much harder problem to solve. Regression is a simpler problem: instead of determining the multidimensional pdf, we would like to infer the expectation of $y$ given $\bf x$.

**The goal of regression usually focuses on estimation of the parameters $\theta$ of the function and their confidence intervals.**

# Metrics:

## Mean squared error (MSE)

The mean squared error (MSE) or mean squared deviation (MSD) of an estimator measures the average of the squares of the errors:

It is defined by (estimation of $y$):

## $MSE \equiv \frac{1}{n} \sum_i (y_i - \hat{y}_i)^2$

## Coefficient of determination (R$^2$)

It is the proportion of the variance in the dependent variable that is predictable from the independent variable(s).

It is defined as:

## $R^2 \equiv 1  - \frac{SS_{res}}{SS_{tot}}$,

where

### $ SS_{res} = \sum_i (y_i - \hat{y}_i)^2$ 

and

### $ SS_{tot} = \sum_i (y_i - \bar{y})^2$

## Variance and bias

The **variance** of an estimator of the parameter $\theta$ will be:

## $Var_{\hat{\theta}} = E_{\hat{\theta}}\biggl[ \bigl(\hat{\theta} - E_{\hat{\theta}}[\hat{\theta}]\bigr)^2\biggr]$

The **bias** of an estimator will be: 

## $Bias_{\hat{\theta}} = E_{\hat{\theta}} - \theta$



## Relation between MSE, variance and bias

Note that when estimating a parameter the MSE equals the variance plus the bias squared



# Linear regression

The one dimensional linear regression model is defined by

$f(x, \theta) = \theta_0 + \theta_1 x$


## Least squares regression for linear regression

## 1. Ordinary least squares

Minimizes the residual sum of squares (RSS) between the observed response values $y_i$ and the model predictions

### $\min RSS = \min \sum\limits_{i=1}^{n} (y_i - \theta_0 - \theta_1 x_i)^2$

Differentiating w.r.t. to $\theta_0$ and $\theta_1$ and setting these derivatives to zero gives two equations which can be solved to give:



## $\hat{\theta}_{1, OLS} = \frac{\sum\limits_{i=1}^n (x_i - \bar{x}) (y_i - \bar{y})}{\sum\limits_{i=1}^n (x_i - \bar{x})^2} = \frac{S_{xy}}{S_{xx}}$

## $\hat{\theta}_{0, OLS} = \bar{y} - \hat{\theta}_{1,OLS} ~\bar{x}$

OLS = ordinary least squares

These estimators are also the maximum likelihood estimator (MLE) when the errors are normally distributed.

The OLS estimator of the error variance $\sigma^2$ is:

## $\hat{\sigma}^2 = S^2 = \frac{RSS}{n - 2}$

If the errors are i.i.d. random variables, the slope and intercept are unbiased and asymptotically normally distributed:

## $\hat{\theta}_{1, OLS} \sim N\biggl(\theta_1, \frac{\sigma^2}{S_{XX}}\biggr)$

## $\hat{\theta}_{0, OLS} \sim N\biggl(\theta_0, \sigma^2 \biggl(\frac{1}{n} + \frac{\bar{x}^2}{S_{XX}}\biggr)\biggr)$

If the errors are also normally distributed with zero mean and known variance $\sigma^2$, then the variance is $\chi^2$ distributed:


## $\hat{\sigma}^2 \sim \biggl(\frac{\sigma^2}{n - 2} \biggr) \chi_{n-2}^{2}$

However, $\sigma^2$ is rarely known, and only the sample variance $S^2$ is known, changing the statistical properties of the estimators.

Defining:

## $ T = \frac{\hat{\theta_1} - \theta_1}{SE}$

where $SE = \frac{S}{\sqrt{S_{XX}}}$ is also known as the **standard error**,

one can show that $T$ is t (t-Student) distributed with $n-2$ degrees of freedom.

Then, the $(1-\alpha)$ confidence interval for $\theta_1$ becomes $(\hat{\theta}_1 - t_{\alpha/2, n-2} ~ SE, \hat{\theta}_1 + t_{\alpha/2, n-2} ~ SE)$ and the confidence interval for $y$ at a given $x$ becomes:

## $y(x) = \hat{\theta}_0 + \hat{\theta}_1 x \pm t_{\alpha/2, n - 2} ~ S ~ \sqrt{1 + \frac{1}{n} \frac{(x - \bar{x})^2}{S_{xx}}}$

## 2. Maximum likelihood estimator (MLE)

Let us consider the probability of observing $\lbrace y_i \rbrace$ given the parameters $\theta$ and the set of independent variable $\lbrace x_i \rbrace$.

Here we assume that the observations are Gaussian with standard deviation $\sigma_i$:

## $p(\lbrace y_i \rbrace | \lbrace x_i \rbrace, \theta) = \prod\limits_{i = 1}^{N} \frac{1}{\sqrt{2 \pi} \sigma_i} \exp \biggl( \frac{- (y_i - (\theta_0 + \theta_1 x_i))^2}{2 \sigma_i^2} \biggr)$


Thus, we can maximize the log likelihood minus a constant to find the **maximum likelihood estimator (MLE)**

## $\log(L) = \log p(\lbrace y_i \rbrace | \lbrace x_i \rbrace, \theta) = \sum\limits_{i = 1}^{N}  \bigl( \frac{- (y_i - (\theta_0 + \theta_1 x_i))^2}{2 \sigma_i^2} \bigr)$

which is **equivalent to minimizing the RSS**

Using Bayes theorem $P(A|B) = P(B|A)~\dfrac{P(A)}{P(B)}$ we can also see that 

## $P(\theta | \lbrace y_i \rbrace, \lbrace x_i \rbrace) = P(\lbrace y_i \rbrace | \lbrace x_i \rbrace, \theta) ~ \frac{p(\theta)}{p(y_i)}$

Therefore, this is also the **maximum a posteriori (MAP)** estimator assuming a flat prior for $\theta$.

## Online estimation

This figure illustrates how the uncertainties in the parameters shrink as the number of observations grow.



### Diabetes dataset example
(https://scikit-learn.org/stable/auto_examples/linear_model/plot_ols.html#sphx-glr-auto-examples-linear-model-plot-ols-py)

We want to try and figure out how bmi affects how sick the person is after one year.


```python
import matplotlib.pyplot as plt
import numpy as np
from sklearn import datasets, linear_model
from sklearn.metrics import mean_squared_error, r2_score

# Load the diabetes dataset
diabetes = datasets.load_diabetes()
print(diabetes.DESCR)
```

    .. _diabetes_dataset:
    
    Diabetes dataset
    ----------------
    
    Ten baseline variables, age, sex, body mass index, average blood
    pressure, and six blood serum measurements were obtained for each of n =
    442 diabetes patients, as well as the response of interest, a
    quantitative measure of disease progression one year after baseline.
    
    **Data Set Characteristics:**
    
      :Number of Instances: 442
    
      :Number of Attributes: First 10 columns are numeric predictive values
    
      :Target: Column 11 is a quantitative measure of disease progression one year after baseline
    
      :Attribute Information:
          - age     age in years
          - sex
          - bmi     body mass index
          - bp      average blood pressure
          - s1      tc, total serum cholesterol
          - s2      ldl, low-density lipoproteins
          - s3      hdl, high-density lipoproteins
          - s4      tch, total cholesterol / HDL
          - s5      ltg, possibly log of serum triglycerides level
          - s6      glu, blood sugar level
    
    Note: Each of these 10 feature variables have been mean centered and scaled by the standard deviation times `n_samples` (i.e. the sum of squares of each column totals 1).
    
    Source URL:
    https://www4.stat.ncsu.edu/~boos/var.select/diabetes.html
    
    For more information see:
    Bradley Efron, Trevor Hastie, Iain Johnstone and Robert Tibshirani (2004) "Least Angle Regression," Annals of Statistics (with discussion), 407-499.
    (https://web.stanford.edu/~hastie/Papers/LARS/LeastAngle_2002.pdf)



```python
# Use only one feature (body mass index)
diabetes_X = diabetes.data[:, np.newaxis, 2]  # axis 2 = bmi

# Split the data into training/testing sets
diabetes_X_train = diabetes_X[:-20]  # last 20 points
diabetes_X_test = diabetes_X[-20:]  # we put them into the test set

# Split the targets into training/testing sets accordingly
diabetes_y_train = diabetes.target[:-20]
diabetes_y_test = diabetes.target[-20:]

# Create linear regression object
regr = linear_model.LinearRegression()

# Train the model using the training sets
regr.fit(diabetes_X_train, diabetes_y_train)

# Make predictions using the testing set
diabetes_y_pred = regr.predict(diabetes_X_test)

# The coefficients
print('Coefficients: \n', regr.coef_)
# The mean squared error
print("Mean squared error (MSE): %.2f"
      % mean_squared_error(diabetes_y_test, diabetes_y_pred))
# Explained variance score: 1 is perfect prediction
print('Variance score (coefficient of determination, R2): %.2f' % r2_score(diabetes_y_test, diabetes_y_pred))
```

    Coefficients: 
     [938.23786125]
    Mean squared error (MSE): 2548.07
    Variance score (coefficient of determination, R2): 0.47



```python
fig, ax = plt.subplots(figsize = (11, 7))

ax.scatter(diabetes_X_test, diabetes_y_test,  c='k')
ax.scatter(diabetes_X_train, diabetes_y_train, c='gray', alpha=0.3)
ax.plot(diabetes_X_test, diabetes_y_pred, color=blue, linewidth=3)
ax.set_xlabel("BMI (scaled)")
ax.set_ylabel("Diabetes disease progression")
```

## 3. Least squares in matrix form

The previous equations can be derived in matrix form

In general one can define regression in terms of a design matrix $M$:

## $y = {\bf M} ~ \theta$

where $y$ is the $N$–dimensional vector of values $y_i$.

## $y = \begin{bmatrix} y_0 \\ y_1 \\ . \\ y_{N-1} \end{bmatrix}$

For straight--line regression, $\theta$ is a two--dimensional vector:

## $\theta = \begin{bmatrix} \theta_0 \\ \theta_1 \end{bmatrix}$

and $M$ is a $2 \times N$ matrix:

## ${\bf M} = \begin{bmatrix} 1 & x_0 \\ 1 & x_1 \\ ... & ... \\ 1 & x_{N-1} \\ \end{bmatrix}$

For the case of heteroscedastic (different) uncertainties, we can define the covariance matrix $C$:

## $ C = \begin{bmatrix} \sigma_0^2 & 0 & . & 0 \\ 0 & \sigma_1^2 & . & 0 \\ . & . & . & . \\ 0 & 0 & . & \sigma_{N - 1}^2\\ \end{bmatrix}$

The MLE solution for this regression is:

## ${\bf \hat{\theta}} = ({\bf M}^T C^{-1} {\bf M})^{-1} ({\bf M}^T C^{-1} {\bf y})$

which is the also the solution that minimizes the equation

## $({\bf y} - {\bf M} ~ \theta)^{-1} C^{-1} ({\bf y} -  {\bf M} ~ \theta)$

The uncertainties on the regression coefficients $\bf \theta$ can be expressed as the following matrix (general case):

## $\Sigma_\theta = \begin{bmatrix} \sigma_{\theta_0}^2 & \dots & \sigma_{\theta_0 \theta_i} & \dots \\ \vdots & \ddots & \vdots & \vdots \\ \sigma_{\theta_0 \theta_i} & \dots & \sigma_{\theta_i}^2 \\  \vdots & \vdots & \vdots & \ddots \\ \end{bmatrix} = [{\bf M}^T C^{-1} {\bf M}]^{-1}$

Values outside of the diagonal appear when they are not fully independent.

In the general case, $C$ can be non--diagonal too, which indicates correlated noise.

We can extend it to...

## 4. Polynomial regression

Let us consider the case where we perform a polynomial representation of the data, e.g.:

## $y_i = \theta_0 + \theta_1 x_i + \theta_2 x_i^2 + \theta_3 x_i^3$

In this case the design matrix $M$ becomes:

## $ M = \begin{bmatrix}
1 & x_0 & x_0^2 & x_0^3 \\
1 & x_1 & x_1^2 & x_1^3 \\
. & . & . & . \\
1 & x_N & x_N^2 & x_N^3 \\
\end{bmatrix}$

And the previous formulae can be used to find the least squares solution since the coefficients are linear

# Regularization

Sometimes we may want to trade an increase in bias for a reduction in variance (e.g. if the matrix $M$ is ill-conditioned the variance could be very large)

## 1. Ridge regression

In this case we include the following regularization term in the least squares equation:

## $({\bf y} - M~ \theta)^{-1} C^{-1}({\bf y} - M~ \theta) + \lambda~ |\theta^T \theta|$

$\lambda$ is called the **regularization or smoothing parameter** and **$|\theta^T \theta|$** is called the **penalty function**

It can be shown that the new solution will be:

## ${\bf \hat{\theta}} = ({\bf M}^T C^{-1} {\bf M} + \lambda I)^{-1} ({\bf M}^T C^{-1} {\bf y})$

It is worth noting that even if $M^T C^{-1} M$ were singular, $M^T C^{-1} M + \lambda I$ may not be singular

=> $\lambda$ wants $\theta$ to have elements with less variance

"Continuous"

### Bayes equivalent

In the Bayesian approach, it can be shown that the MAP solution is the same as the ridge solution if we assume a Gaussian prior for the parameters $\theta$:


## $ p(\theta) \propto \exp \biggl( \frac{- \lambda (\theta^T \theta)}{2} \biggr)$

Ridge regression effectively constrains the norm of the parameters $\theta$:

## $|{\bf \theta}|^2 < s$

## 2. LASSO regression

Another common regularization technique

In this case we include a different regularization term in the least squares equation:

## $({\bf y} - \theta {\bf X})^{-1}({\bf y} - \theta {\bf X}) + \lambda |\theta|$

Now we use the $\ell_1$ norm instead of the $\ell_2$ norm.

In this case there is no closed form, but it can be shown that the MAP solution is the same if an exponential prior for the parameters $\theta$ is assumed.

LASSO regression effectively constrains the $\ell_1$ norm of the parameters $\theta$:

## $|{\bf \theta}| < s$

=> $\lambda$ wants $\theta$ to have less elements

"Discrete"

### Understanding Ridge and LASSO regression



# Comparison between ordinary, Ridge and LASSO regression


```python
from sklearn import linear_model

X_train = np.c_[.5, 1].T
y_train = [.5, 1]
X_test = np.c_[0, 2].T

np.random.seed(0)

classifiers = dict(ols=linear_model.LinearRegression(),
                   ridge=linear_model.Ridge(alpha=.05),
                    lasso=linear_model.Lasso(alpha=.02))

for name, clf in classifiers.items():
    fig, ax = plt.subplots(figsize=(4, 3))

    for _ in range(6):
        this_X = .1 * np.random.normal(size=(2, 1)) + X_train
        clf.fit(this_X, y_train)

        ax.plot(X_test, clf.predict(X_test), color='gray')
        ax.scatter(this_X, y_train, s=3, c='gray', marker='o', zorder=10)

    clf.fit(X_train, y_train)
    ax.plot(X_test, clf.predict(X_test), linewidth=2, color=blue)
    ax.scatter(X_train, y_train, s=30, c=red, marker='+', zorder=10)

    ax.set_title(name)
    ax.set_xlim(0, 2)
    ax.set_ylim((0, 1.6))
    ax.set_xlabel('X')
    ax.set_ylabel('y')

    fig.tight_layout()

plt.show()
```

## Choosing the regularization parameter $\lambda$

One way to choose $\lambda$ is to do k-fold cross-validation on $\lambda$, defining the error:

## $ Err(\lambda) = k^{1} \sum_k N_k^{-1} \sum\limits_{i}^{N_k} \frac{[y_i - f(x_i | \theta)]^2}{\sigma_i^2}$

and then choose the $\lambda$ with the smallest error.

# Outliers
## Huber loss

The Huber estimator minimizes

### $ \sum\limits_{i=1}^N e(y_i|y)$

where $e(y_i|y)$ is modeled as:

### $\begin{equation}
  \phi(t)=\left\{
  \begin{array}{@{}ll@{}}
    \frac{1}{2} t^2, & \text{if}\ |t|\le c \\
    c|t| - \frac{1}{2} c^2, & \text{otherwise}
  \end{array}\right.
\end{equation}$

With $t = y_i - \hat{y}_i$

c is a new parameter

If it's close to the line, it is treated as a Gaussian. If it is far, as an exponencial.





## Bayesian Outlier Methods

Alternative to the Huber loss, we can modify the likelihood function by assuming a mixture model, where some points belong to an outlier distribution, which is modeled as a Gaussian with mean $\mu_b$ and variance $V_b$, and others belong to the main population, which follow the linear model.

1. Assume that the probability of belonging to the outlier population is $p_b$



2. Assume that each point can be tagged as either outlier ($g_i = 0$) or not ($g_i = 1$)






```python

```


```python
from sklearn.datasets import make_regression
from sklearn.linear_model import HuberRegressor, Ridge

# Generate toy data.
rng = np.random.RandomState(0)
X, y = make_regression(n_samples=20, n_features=1, random_state=0, noise=4.0,
                       bias=100.0)

# Add four strong outliers to the dataset.
X_outliers = rng.normal(0, 0.5, size=(4, 1))
y_outliers = rng.normal(0, 2.0, size=4)
X_outliers[:2, :] += X.max() + X.mean() / 4.
X_outliers[2:, :] += X.min() - X.mean() / 4.
y_outliers[:2] += y.min() - y.mean() / 4.
y_outliers[2:] += y.max() + y.mean() / 4.
X = np.vstack((X, X_outliers))
y = np.concatenate((y, y_outliers))
plt.plot(X, y, 'b.')

# Fit the huber regressor over a series of epsilon values.
colors = ['r-', 'b-', 'y-', 'm-']

x = np.linspace(X.min(), X.max(), 7)
epsilon_values = [1.35, 1.5, 1.75, 1.9]
for k, epsilon in enumerate(epsilon_values):
    huber = HuberRegressor(fit_intercept=True, alpha=0.0, max_iter=100,
                           epsilon=epsilon)
    huber.fit(X, y)
    coef_ = huber.coef_ * x + huber.intercept_
    plt.plot(x, coef_, colors[k], label="huber loss, %s" % epsilon)

# Fit a ridge regressor to compare it to huber regressor.
ridge = Ridge(fit_intercept=True, alpha=0.0, random_state=0, normalize=True)
ridge.fit(X, y)
coef_ridge = ridge.coef_
coef_ = ridge.coef_ * x + ridge.intercept_
plt.plot(x, coef_, 'g-', label="ridge regression")

plt.title("Comparison of HuberRegressor vs Ridge")
plt.xlabel("X")
plt.ylabel("y")
plt.legend(loc=0)
plt.show()
```

# Cross-Validation Testing

Analogous to the classification problem, one can also do cross-validation during regression to determine the best hyperparameters. For example, choosing the regularization parameter for ridge or lasso regression, or choosing the polynomial degree of a polynomial fit.

=> **how to know when to stop increasing the degree of the polynomial?**

We want to get close to the data points, but not EXACTLY. We want the underlying distribution, not get this specific set.



We want to be somewhere in the good range that makes it better for both sets.




```python

```

ea2fed3bb9d8b3a046845d17113d099b967570d0
ipynb
Jupyter Notebook
05-Linear Regression.ipynb
NicoGalvarino/AS450-astroinformatica
2e40b6809287c231e367ec26be49cc043dfb8dcd
MIT
05-Linear Regression.ipynb
NicoGalvarino/AS450-astroinformatica
2e40b6809287c231e367ec26be49cc043dfb8dcd
MIT
05-Linear Regression.ipynb
NicoGalvarino/AS450-astroinformatica
2e40b6809287c231e367ec26be49cc043dfb8dcd
MIT
Qwen/Qwen-72B
1. YES
2. YES