lars/FYS-STK4155
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

updates

Browse Source
This commit is contained in:
Morten Hjorth-Jensen
2025-08-18 10:47:43 +02:00
parent d4db3104e6
commit 35c45696d6
37 changed files with 5815 additions and 5344 deletions
Download Patch File Download Diff File Expand all files Collapse all files
doc/LectureNotes/_build/html/chapter3.html
+129 -129
View File
@@ -479,7 +479,7 @@ independent, i.e.:</p>
<p>The randomness of <span class="math notranslate nohighlight">\(\varepsilon_i\)</span> implies that
<span class="math notranslate nohighlight">\(\mathbf{y}_i\)</span> is also a random variable. In particular,
<span class="math notranslate nohighlight">\(\mathbf{y}_i\)</span> is normally distributed, because <span class="math notranslate nohighlight">\(\varepsilon_i \sim
\mathcal{N}(0, \sigma^2)\)</span> and <span class="math notranslate nohighlight">\(\mathbf{X}_{i,\ast} \, \boldsymbol{\beta}\)</span> is a
\mathcal{N}(0, \sigma^2)\)</span> and <span class="math notranslate nohighlight">\(\mathbf{X}_{i,\ast} \, \boldsymbol{\theta}\)</span> is a
non-random scalar. To specify the parameters of the distribution of
<span class="math notranslate nohighlight">\(\mathbf{y}_i\)</span> we need to calculate its first two moments.</p>
<p>Recall that <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span> is a matrix of dimensionality <span class="math notranslate nohighlight">\(n\times p\)</span>. The
@@ -496,15 +496,15 @@ which describe our data</p>
function <span class="math notranslate nohighlight">\(f\)</span> is approximated by <span class="math notranslate nohighlight">\(\boldsymbol{\tilde{y}}\)</span> where we want to minimize <span class="math notranslate nohighlight">\((\boldsymbol{y}-\boldsymbol{\tilde{y}})^2\)</span>, our MSE, with</p>
<div class="math notranslate nohighlight">
\[
\boldsymbol{\tilde{y}} = \boldsymbol{X}\boldsymbol{\beta}.
\boldsymbol{\tilde{y}} = \boldsymbol{X}\boldsymbol{\theta}.
\]</div>
<p>We can calculate the expectation value of <span class="math notranslate nohighlight">\(\boldsymbol{y}\)</span> for a given element <span class="math notranslate nohighlight">\(i\)</span></p>
<div class="math notranslate nohighlight">
\[
\begin{align*}
\mathbb{E}(y_i) &amp; =
\mathbb{E}(\mathbf{X}_{i, \ast} \, \boldsymbol{\beta}) + \mathbb{E}(\varepsilon_i)
\, \, \, = \, \, \, \mathbf{X}_{i, \ast} \, \beta,
\mathbb{E}(\mathbf{X}_{i, \ast} \, \boldsymbol{\theta}) + \mathbb{E}(\varepsilon_i)
\, \, \, = \, \, \, \mathbf{X}_{i, \ast} \, \theta,
\end{align*}
\]</div>
<p>while
@@ -514,81 +514,81 @@ its variance is</p>
\begin{align*} \mbox{Var}(y_i) &amp; = \mathbb{E} \{ [y_i
- \mathbb{E}(y_i)]^2 \} \, \, \, = \, \, \, \mathbb{E} ( y_i^2 ) -
[\mathbb{E}(y_i)]^2 \\ &amp; = \mathbb{E} [ ( \mathbf{X}_{i, \ast} \,
\beta + \varepsilon_i )^2] - ( \mathbf{X}_{i, \ast} \, \boldsymbol{\beta})^2 \\ &amp;
= \mathbb{E} [ ( \mathbf{X}_{i, \ast} \, \boldsymbol{\beta})^2 + 2 \varepsilon_i
\mathbf{X}_{i, \ast} \, \boldsymbol{\beta} + \varepsilon_i^2 ] - ( \mathbf{X}_{i,
\ast} \, \beta)^2 \\ &amp; = ( \mathbf{X}_{i, \ast} \, \boldsymbol{\beta})^2 + 2
\mathbb{E}(\varepsilon_i) \mathbf{X}_{i, \ast} \, \boldsymbol{\beta} +
\mathbb{E}(\varepsilon_i^2 ) - ( \mathbf{X}_{i, \ast} \, \boldsymbol{\beta})^2
\theta + \varepsilon_i )^2] - ( \mathbf{X}_{i, \ast} \, \boldsymbol{\theta})^2 \\ &amp;
= \mathbb{E} [ ( \mathbf{X}_{i, \ast} \, \boldsymbol{\theta})^2 + 2 \varepsilon_i
\mathbf{X}_{i, \ast} \, \boldsymbol{\theta} + \varepsilon_i^2 ] - ( \mathbf{X}_{i,
\ast} \, \theta)^2 \\ &amp; = ( \mathbf{X}_{i, \ast} \, \boldsymbol{\theta})^2 + 2
\mathbb{E}(\varepsilon_i) \mathbf{X}_{i, \ast} \, \boldsymbol{\theta} +
\mathbb{E}(\varepsilon_i^2 ) - ( \mathbf{X}_{i, \ast} \, \boldsymbol{\theta})^2
\\ &amp; = \mathbb{E}(\varepsilon_i^2 ) \, \, \, = \, \, \,
\mbox{Var}(\varepsilon_i) \, \, \, = \, \, \, \sigma^2.
\end{align*}
\end{split}\]</div>
<p>Hence, <span class="math notranslate nohighlight">\(y_i \sim \mathcal{N}( \mathbf{X}_{i, \ast} \, \boldsymbol{\beta}, \sigma^2)\)</span>, that is <span class="math notranslate nohighlight">\(\boldsymbol{y}\)</span> follows a normal distribution with
mean value <span class="math notranslate nohighlight">\(\boldsymbol{X}\boldsymbol{\beta}\)</span> and variance <span class="math notranslate nohighlight">\(\sigma^2\)</span> (not be confused with the singular values of the SVD).</p>
<p>With the OLS expressions for the parameters <span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span> we can evaluate the expectation value</p>
<p>Hence, <span class="math notranslate nohighlight">\(y_i \sim \mathcal{N}( \mathbf{X}_{i, \ast} \, \boldsymbol{\theta}, \sigma^2)\)</span>, that is <span class="math notranslate nohighlight">\(\boldsymbol{y}\)</span> follows a normal distribution with
mean value <span class="math notranslate nohighlight">\(\boldsymbol{X}\boldsymbol{\theta}\)</span> and variance <span class="math notranslate nohighlight">\(\sigma^2\)</span> (not be confused with the singular values of the SVD).</p>
<p>With the OLS expressions for the parameters <span class="math notranslate nohighlight">\(\boldsymbol{\theta}\)</span> we can evaluate the expectation value</p>
<div class="math notranslate nohighlight">
\[
\mathbb{E}(\boldsymbol{\beta}) = \mathbb{E}[ (\mathbf{X}^{\top} \mathbf{X})^{-1}\mathbf{X}^{T} \mathbf{Y}]=(\mathbf{X}^{T} \mathbf{X})^{-1}\mathbf{X}^{T} \mathbb{E}[ \mathbf{Y}]=(\mathbf{X}^{T} \mathbf{X})^{-1} \mathbf{X}^{T}\mathbf{X}\boldsymbol{\beta}=\boldsymbol{\beta}.
\mathbb{E}(\boldsymbol{\theta}) = \mathbb{E}[ (\mathbf{X}^{\top} \mathbf{X})^{-1}\mathbf{X}^{T} \mathbf{Y}]=(\mathbf{X}^{T} \mathbf{X})^{-1}\mathbf{X}^{T} \mathbb{E}[ \mathbf{Y}]=(\mathbf{X}^{T} \mathbf{X})^{-1} \mathbf{X}^{T}\mathbf{X}\boldsymbol{\theta}=\boldsymbol{\theta}.
\]</div>
<p>This means that the estimator of the regression parameters is unbiased.</p>
<p>We can also calculate the variance</p>
<p>The variance of <span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span> is</p>
<p>The variance of <span class="math notranslate nohighlight">\(\boldsymbol{\theta}\)</span> is</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\begin{eqnarray*}
\mbox{Var}(\boldsymbol{\beta}) &amp; = &amp; \mathbb{E} \{ [\boldsymbol{\beta} - \mathbb{E}(\boldsymbol{\beta})] [\boldsymbol{\beta} - \mathbb{E}(\boldsymbol{\beta})]^{T} \}
\mbox{Var}(\boldsymbol{\theta}) &amp; = &amp; \mathbb{E} \{ [\boldsymbol{\theta} - \mathbb{E}(\boldsymbol{\theta})] [\boldsymbol{\theta} - \mathbb{E}(\boldsymbol{\theta})]^{T} \}
\\
&amp; = &amp; \mathbb{E} \{ [(\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \mathbf{Y} - \boldsymbol{\beta}] \, [(\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \mathbf{Y} - \boldsymbol{\beta}]^{T} \}
&amp; = &amp; \mathbb{E} \{ [(\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \mathbf{Y} - \boldsymbol{\theta}] \, [(\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \mathbf{Y} - \boldsymbol{\theta}]^{T} \}
\\
% &amp; = &amp; \mathbb{E} \{ [(\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \mathbf{Y}] \, [(\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \mathbf{Y}]^{T} \} - \boldsymbol{\beta} \, \boldsymbol{\beta}^{T}
% &amp; = &amp; \mathbb{E} \{ [(\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \mathbf{Y}] \, [(\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \mathbf{Y}]^{T} \} - \boldsymbol{\theta} \, \boldsymbol{\theta}^{T}
% \\
% &amp; = &amp; \mathbb{E} \{ (\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \mathbf{Y} \, \mathbf{Y}^{T} \, \mathbf{X} \, (\mathbf{X}^{T} \mathbf{X})^{-1} \} - \boldsymbol{\beta} \, \boldsymbol{\beta}^{T}
% &amp; = &amp; \mathbb{E} \{ (\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \mathbf{Y} \, \mathbf{Y}^{T} \, \mathbf{X} \, (\mathbf{X}^{T} \mathbf{X})^{-1} \} - \boldsymbol{\theta} \, \boldsymbol{\theta}^{T}
% \\
&amp; = &amp; (\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \, \mathbb{E} \{ \mathbf{Y} \, \mathbf{Y}^{T} \} \, \mathbf{X} \, (\mathbf{X}^{T} \mathbf{X})^{-1} - \boldsymbol{\beta} \, \boldsymbol{\beta}^{T}
&amp; = &amp; (\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \, \mathbb{E} \{ \mathbf{Y} \, \mathbf{Y}^{T} \} \, \mathbf{X} \, (\mathbf{X}^{T} \mathbf{X})^{-1} - \boldsymbol{\theta} \, \boldsymbol{\theta}^{T}
\\
&amp; = &amp; (\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \, \{ \mathbf{X} \, \boldsymbol{\beta} \, \boldsymbol{\beta}^{T} \, \mathbf{X}^{T} + \sigma^2 \} \, \mathbf{X} \, (\mathbf{X}^{T} \mathbf{X})^{-1} - \boldsymbol{\beta} \, \boldsymbol{\beta}^{T}
&amp; = &amp; (\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \, \{ \mathbf{X} \, \boldsymbol{\theta} \, \boldsymbol{\theta}^{T} \, \mathbf{X}^{T} + \sigma^2 \} \, \mathbf{X} \, (\mathbf{X}^{T} \mathbf{X})^{-1} - \boldsymbol{\theta} \, \boldsymbol{\theta}^{T}
% \\
% &amp; = &amp; (\mathbf{X}^T \mathbf{X})^{-1} \, \mathbf{X}^T \, \mathbf{X} \, \boldsymbol{\beta} \, \boldsymbol{\beta}^T \, \mathbf{X}^T \, \mathbf{X} \, (\mathbf{X}^T % \mathbf{X})^{-1}
% &amp; = &amp; (\mathbf{X}^T \mathbf{X})^{-1} \, \mathbf{X}^T \, \mathbf{X} \, \boldsymbol{\theta} \, \boldsymbol{\theta}^T \, \mathbf{X}^T \, \mathbf{X} \, (\mathbf{X}^T % \mathbf{X})^{-1}
% \\
% &amp; &amp; + \, \, \sigma^2 \, (\mathbf{X}^T \mathbf{X})^{-1} \, \mathbf{X}^T \, \mathbf{X} \, (\mathbf{X}^T \mathbf{X})^{-1} - \boldsymbol{\beta} \boldsymbol{\beta}^T
% &amp; &amp; + \, \, \sigma^2 \, (\mathbf{X}^T \mathbf{X})^{-1} \, \mathbf{X}^T \, \mathbf{X} \, (\mathbf{X}^T \mathbf{X})^{-1} - \boldsymbol{\theta} \boldsymbol{\theta}^T
\\
&amp; = &amp; \boldsymbol{\beta} \, \boldsymbol{\beta}^{T} + \sigma^2 \, (\mathbf{X}^{T} \mathbf{X})^{-1} - \boldsymbol{\beta} \, \boldsymbol{\beta}^{T}
&amp; = &amp; \boldsymbol{\theta} \, \boldsymbol{\theta}^{T} + \sigma^2 \, (\mathbf{X}^{T} \mathbf{X})^{-1} - \boldsymbol{\theta} \, \boldsymbol{\theta}^{T}
\, \, \, = \, \, \, \sigma^2 \, (\mathbf{X}^{T} \mathbf{X})^{-1},
\end{eqnarray*}
\end{split}\]</div>
<p>where we have used that <span class="math notranslate nohighlight">\(\mathbb{E} (\mathbf{Y} \mathbf{Y}^{T}) =
\mathbf{X} \, \boldsymbol{\beta} \, \boldsymbol{\beta}^{T} \, \mathbf{X}^{T} +
\sigma^2 \, \mathbf{I}_{nn}\)</span>. From <span class="math notranslate nohighlight">\(\mbox{Var}(\boldsymbol{\beta}) = \sigma^2
\mathbf{X} \, \boldsymbol{\theta} \, \boldsymbol{\theta}^{T} \, \mathbf{X}^{T} +
\sigma^2 \, \mathbf{I}_{nn}\)</span>. From <span class="math notranslate nohighlight">\(\mbox{Var}(\boldsymbol{\theta}) = \sigma^2
\, (\mathbf{X}^{T} \mathbf{X})^{-1}\)</span>, one obtains an estimate of the
variance of the estimate of the <span class="math notranslate nohighlight">\(j\)</span>-th regression coefficient:
<span class="math notranslate nohighlight">\(\boldsymbol{\sigma}^2 (\boldsymbol{\beta}_j ) = \boldsymbol{\sigma}^2 \sqrt{
<span class="math notranslate nohighlight">\(\boldsymbol{\sigma}^2 (\boldsymbol{\theta}_j ) = \boldsymbol{\sigma}^2 \sqrt{
[(\mathbf{X}^{T} \mathbf{X})^{-1}]_{jj} }\)</span>. This may be used to
construct a confidence interval for the estimates.</p>
<p>In a similar way, we can obtain analytical expressions for say the
expectation values of the parameters <span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span> and their variance
expectation values of the parameters <span class="math notranslate nohighlight">\(\boldsymbol{\theta}\)</span> and their variance
when we employ Ridge regression, allowing us again to define a confidence interval.</p>
<p>It is rather straightforward to show that</p>
<div class="math notranslate nohighlight">
\[
\mathbb{E} \big[ \boldsymbol{\beta}^{\mathrm{Ridge}} \big]=(\mathbf{X}^{T} \mathbf{X} + \lambda \mathbf{I}_{pp})^{-1} (\mathbf{X}^{\top} \mathbf{X})\boldsymbol{\beta}^{\mathrm{OLS}}.
\mathbb{E} \big[ \boldsymbol{\theta}^{\mathrm{Ridge}} \big]=(\mathbf{X}^{T} \mathbf{X} + \lambda \mathbf{I}_{pp})^{-1} (\mathbf{X}^{\top} \mathbf{X})\boldsymbol{\theta}^{\mathrm{OLS}}.
\]</div>
<p>We see clearly that
<span class="math notranslate nohighlight">\(\mathbb{E} \big[ \boldsymbol{\beta}^{\mathrm{Ridge}} \big] \not= \boldsymbol{\beta}^{\mathrm{OLS}}\)</span> for any <span class="math notranslate nohighlight">\(\lambda &gt; 0\)</span>. We say then that the ridge estimator is biased.</p>
<span class="math notranslate nohighlight">\(\mathbb{E} \big[ \boldsymbol{\theta}^{\mathrm{Ridge}} \big] \not= \boldsymbol{\theta}^{\mathrm{OLS}}\)</span> for any <span class="math notranslate nohighlight">\(\lambda &gt; 0\)</span>. We say then that the ridge estimator is biased.</p>
<p>We can also compute the variance as</p>
<div class="math notranslate nohighlight">
\[
\mbox{Var}[\boldsymbol{\beta}^{\mathrm{Ridge}}]=\sigma^2[ \mathbf{X}^{T} \mathbf{X} + \lambda \mathbf{I} ]^{-1} \mathbf{X}^{T} \mathbf{X} \{ [ \mathbf{X}^{\top} \mathbf{X} + \lambda \mathbf{I} ]^{-1}\}^{T},
\mbox{Var}[\boldsymbol{\theta}^{\mathrm{Ridge}}]=\sigma^2[ \mathbf{X}^{T} \mathbf{X} + \lambda \mathbf{I} ]^{-1} \mathbf{X}^{T} \mathbf{X} \{ [ \mathbf{X}^{\top} \mathbf{X} + \lambda \mathbf{I} ]^{-1}\}^{T},
\]</div>
<p>and it is easy to see that if the parameter <span class="math notranslate nohighlight">\(\lambda\)</span> goes to infinity then the variance of Ridge parameters <span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span> goes to zero.</p>
<p>and it is easy to see that if the parameter <span class="math notranslate nohighlight">\(\lambda\)</span> goes to infinity then the variance of Ridge parameters <span class="math notranslate nohighlight">\(\boldsymbol{\theta}\)</span> goes to zero.</p>
<p>With this, we can compute the difference</p>
<div class="math notranslate nohighlight">
\[
\mbox{Var}[\boldsymbol{\beta}^{\mathrm{OLS}}]-\mbox{Var}(\boldsymbol{\beta}^{\mathrm{Ridge}})=\sigma^2 [ \mathbf{X}^{T} \mathbf{X} + \lambda \mathbf{I} ]^{-1}[ 2\lambda\mathbf{I} + \lambda^2 (\mathbf{X}^{T} \mathbf{X})^{-1} ] \{ [ \mathbf{X}^{T} \mathbf{X} + \lambda \mathbf{I} ]^{-1}\}^{T}.
\mbox{Var}[\boldsymbol{\theta}^{\mathrm{OLS}}]-\mbox{Var}(\boldsymbol{\theta}^{\mathrm{Ridge}})=\sigma^2 [ \mathbf{X}^{T} \mathbf{X} + \lambda \mathbf{I} ]^{-1}[ 2\lambda\mathbf{I} + \lambda^2 (\mathbf{X}^{T} \mathbf{X})^{-1} ] \{ [ \mathbf{X}^{T} \mathbf{X} + \lambda \mathbf{I} ]^{-1}\}^{T}.
\]</div>
<p>The difference is non-negative definite since each component of the
matrix product is non-negative definite.
This means the variance we obtain with the standard OLS will always for <span class="math notranslate nohighlight">\(\lambda &gt; 0\)</span> be larger than the variance of <span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span> obtained with the Ridge estimator. This has interesting consequences when we discuss the so-called bias-variance trade-off below.</p>
This means the variance we obtain with the standard OLS will always for <span class="math notranslate nohighlight">\(\lambda &gt; 0\)</span> be larger than the variance of <span class="math notranslate nohighlight">\(\boldsymbol{\theta}\)</span> obtained with the Ridge estimator. This has interesting consequences when we discuss the so-called bias-variance trade-off below.</p>
</section>
<section id="id1">
<h2><span class="section-number">5.3. </span>Resampling methods<a class="headerlink" href="#id1" title="Link to this heading">#</a></h2>
@@ -620,7 +620,7 @@ satisfied, the methods will fail. Yet, it should be said that if the data are
independent, identically distributed, and we only want to estimate the
variance of <span class="math notranslate nohighlight">\(\overline{X}\)</span> (which often is the case), then there is no
need for bootstrapping.</p>
<p>The Jackknife works by making many replicas of the estimator <span class="math notranslate nohighlight">\(\widehat{\beta}\)</span>.
<p>The Jackknife works by making many replicas of the estimator <span class="math notranslate nohighlight">\(\widehat{\theta}\)</span>.
The jackknife is a resampling method where we systematically leave out one observation from the vector of observed values <span class="math notranslate nohighlight">\(\boldsymbol{x} = (x_1,x_2,\cdots,X_n)\)</span>.
Let <span class="math notranslate nohighlight">\(\boldsymbol{x}_i\)</span> denote the vector</p>
<div class="math notranslate nohighlight">
@@ -629,8 +629,8 @@ Let <span class="math notranslate nohighlight">\(\boldsymbol{x}_i\)</span> denot
\]</div>
<p>which equals the vector <span class="math notranslate nohighlight">\(\boldsymbol{x}\)</span> with the exception that observation
number <span class="math notranslate nohighlight">\(i\)</span> is left out. Using this notation, define
<span class="math notranslate nohighlight">\(\widehat{\beta}_i\)</span> to be the estimator
<span class="math notranslate nohighlight">\(\widehat{\beta}\)</span> computed using <span class="math notranslate nohighlight">\(\vec{X}_i\)</span>.</p>
<span class="math notranslate nohighlight">\(\widehat{\theta}_i\)</span> to be the estimator
<span class="math notranslate nohighlight">\(\widehat{\theta}\)</span> computed using <span class="math notranslate nohighlight">\(\vec{X}_i\)</span>.</p>
<div class="cell docutils container">
<div class="cell_input docutils container">
<div class="highlight-none notranslate"><div class="highlight"><pre><span></span>from numpy import *
@@ -677,27 +677,27 @@ advantages:</p>
<li><p>It is possible to apply the bootstrap to statistics with sampling distributions that are difficult to derive, even asymptotically.</p></li>
<li><p>It is relatively simple to apply the bootstrap to complex data-collection plans (such as stratified and clustered samples).</p></li>
</ol>
<p>Since <span class="math notranslate nohighlight">\(\widehat{\beta} = \widehat{\beta}(\boldsymbol{X})\)</span> is a function of random variables,
<span class="math notranslate nohighlight">\(\widehat{\beta}\)</span> itself must be a random variable. Thus it has
<p>Since <span class="math notranslate nohighlight">\(\widehat{\theta} = \widehat{\theta}(\boldsymbol{X})\)</span> is a function of random variables,
<span class="math notranslate nohighlight">\(\widehat{\theta}\)</span> itself must be a random variable. Thus it has
a pdf, call this function <span class="math notranslate nohighlight">\(p(\boldsymbol{t})\)</span>. The aim of the bootstrap is to
estimate <span class="math notranslate nohighlight">\(p(\boldsymbol{t})\)</span> by the relative frequency of
<span class="math notranslate nohighlight">\(\widehat{\beta}\)</span>. You can think of this as using a histogram
<span class="math notranslate nohighlight">\(\widehat{\theta}\)</span>. You can think of this as using a histogram
in the place of <span class="math notranslate nohighlight">\(p(\boldsymbol{t})\)</span>. If the relative frequency closely
resembles <span class="math notranslate nohighlight">\(p(\vec{t})\)</span>, then using numerics, it is straight forward to
estimate all the interesting parameters of <span class="math notranslate nohighlight">\(p(\boldsymbol{t})\)</span> using point
estimators.</p>
<p>In the case that <span class="math notranslate nohighlight">\(\widehat{\beta}\)</span> has
<p>In the case that <span class="math notranslate nohighlight">\(\widehat{\theta}\)</span> has
more than one component, and the components are independent, we use the
same estimator on each component separately. If the probability
density function of <span class="math notranslate nohighlight">\(X_i\)</span>, <span class="math notranslate nohighlight">\(p(x)\)</span>, had been known, then it would have
been straight forward to do this by:</p>
<ol class="arabic simple">
<li><p>Drawing lots of numbers from <span class="math notranslate nohighlight">\(p(x)\)</span>, suppose we call one such set of numbers <span class="math notranslate nohighlight">\((X_1^*, X_2^*, \cdots, X_n^*)\)</span>.</p></li>
<li><p>Then using these numbers, we could compute a replica of <span class="math notranslate nohighlight">\(\widehat{\beta}\)</span> called <span class="math notranslate nohighlight">\(\widehat{\beta}^*\)</span>.</p></li>
<li><p>Then using these numbers, we could compute a replica of <span class="math notranslate nohighlight">\(\widehat{\theta}\)</span> called <span class="math notranslate nohighlight">\(\widehat{\theta}^*\)</span>.</p></li>
</ol>
<p>By repeated use of (1) and (2), many
estimates of <span class="math notranslate nohighlight">\(\widehat{\beta}\)</span> could have been obtained. The
idea is to use the relative frequency of <span class="math notranslate nohighlight">\(\widehat{\beta}^*\)</span>
estimates of <span class="math notranslate nohighlight">\(\widehat{\theta}\)</span> could have been obtained. The
idea is to use the relative frequency of <span class="math notranslate nohighlight">\(\widehat{\theta}^*\)</span>
(think of a histogram) as an estimate of <span class="math notranslate nohighlight">\(p(\boldsymbol{t})\)</span>.</p>
<p>But
unless there is enough information available about the process that
@@ -715,18 +715,18 @@ frequency of the observation <span class="math notranslate nohighlight">\(X_i\)<
<ol class="arabic simple">
<li><p>Draw with replacement <span class="math notranslate nohighlight">\(n\)</span> numbers for the observed variables <span class="math notranslate nohighlight">\(\boldsymbol{x} = (x_1,x_2,\cdots,x_n)\)</span>.</p></li>
<li><p>Define a vector <span class="math notranslate nohighlight">\(\boldsymbol{x}^*\)</span> containing the values which were drawn from <span class="math notranslate nohighlight">\(\boldsymbol{x}\)</span>.</p></li>
<li><p>Using the vector <span class="math notranslate nohighlight">\(\boldsymbol{x}^*\)</span> compute <span class="math notranslate nohighlight">\(\widehat{\beta}^*\)</span> by evaluating <span class="math notranslate nohighlight">\(\widehat \beta\)</span> under the observations <span class="math notranslate nohighlight">\(\boldsymbol{x}^*\)</span>.</p></li>
<li><p>Using the vector <span class="math notranslate nohighlight">\(\boldsymbol{x}^*\)</span> compute <span class="math notranslate nohighlight">\(\widehat{\theta}^*\)</span> by evaluating <span class="math notranslate nohighlight">\(\widehat \theta\)</span> under the observations <span class="math notranslate nohighlight">\(\boldsymbol{x}^*\)</span>.</p></li>
<li><p>Repeat this process <span class="math notranslate nohighlight">\(k\)</span> times.</p></li>
</ol>
<p>When you are done, you can draw a histogram of the relative frequency
of <span class="math notranslate nohighlight">\(\widehat \beta^*\)</span>. This is your estimate of the probability
of <span class="math notranslate nohighlight">\(\widehat \theta^*\)</span>. This is your estimate of the probability
distribution <span class="math notranslate nohighlight">\(p(t)\)</span>. Using this probability distribution you can
estimate any statistics thereof. In principle you never draw the
histogram of the relative frequency of <span class="math notranslate nohighlight">\(\widehat{\beta}^*\)</span>. Instead
histogram of the relative frequency of <span class="math notranslate nohighlight">\(\widehat{\theta}^*\)</span>. Instead
you use the estimators corresponding to the statistic of interest. For
example, if you are interested in estimating the variance of <span class="math notranslate nohighlight">\(\widehat
\beta\)</span>, apply the estimator <span class="math notranslate nohighlight">\(\widehat \sigma^2\)</span> to the values
<span class="math notranslate nohighlight">\(\widehat \beta^*\)</span>.</p>
\theta\)</span>, apply the estimator <span class="math notranslate nohighlight">\(\widehat \sigma^2\)</span> to the values
<span class="math notranslate nohighlight">\(\widehat \theta^*\)</span>.</p>
<p>Before we proceed however, we need to remind ourselves about a central
theorem in statistics, namely the so-called <strong>central limit theorem</strong>.
This theorem plays a central role in understanding why the Bootstrap
@@ -911,12 +911,12 @@ tasks. Consider a dataset <span class="math notranslate nohighlight">\(\mathcal{
<p>where <span class="math notranslate nohighlight">\(\epsilon\)</span> is normally distributed with mean zero and standard deviation <span class="math notranslate nohighlight">\(\sigma^2\)</span>.</p>
<p>In our derivation of the ordinary least squares method we defined then
an approximation to the function <span class="math notranslate nohighlight">\(f\)</span> in terms of the parameters
<span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span> and the design matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span> which embody our model,
that is <span class="math notranslate nohighlight">\(\boldsymbol{\tilde{y}}=\boldsymbol{X}\boldsymbol{\beta}\)</span>.</p>
<p>Thereafter we found the parameters <span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span> by optimizing the means squared error via the so-called cost function</p>
<span class="math notranslate nohighlight">\(\boldsymbol{\theta}\)</span> and the design matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span> which embody our model,
that is <span class="math notranslate nohighlight">\(\boldsymbol{\tilde{y}}=\boldsymbol{X}\boldsymbol{\theta}\)</span>.</p>
<p>Thereafter we found the parameters <span class="math notranslate nohighlight">\(\boldsymbol{\theta}\)</span> by optimizing the means squared error via the so-called cost function</p>
<div class="math notranslate nohighlight">
\[
C(\boldsymbol{X},\boldsymbol{\beta}) =\frac{1}{n}\sum_{i=0}^{n-1}(y_i-\tilde{y}_i)^2=\mathbb{E}\left[(\boldsymbol{y}-\boldsymbol{\tilde{y}})^2\right].
C(\boldsymbol{X},\boldsymbol{\theta}) =\frac{1}{n}\sum_{i=0}^{n-1}(y_i-\tilde{y}_i)^2=\mathbb{E}\left[(\boldsymbol{y}-\boldsymbol{\tilde{y}})^2\right].
\]</div>
<p>We can rewrite this as</p>
<div class="math notranslate nohighlight">
@@ -1267,20 +1267,20 @@ cross-validation (LOOCV).</p>
<div class="math notranslate nohighlight">
\[
\begin{align*}
\boldsymbol{\beta}_{-i}(\lambda) &amp; = ( \boldsymbol{X}_{-i, \ast}^{T}
\boldsymbol{\theta}_{-i}(\lambda) &amp; = ( \boldsymbol{X}_{-i, \ast}^{T}
\boldsymbol{X}_{-i, \ast} + \lambda \boldsymbol{I}_{pp})^{-1}
\boldsymbol{X}_{-i, \ast}^{T} \boldsymbol{y}_{-i}
\end{align*}
\]</div>
<ul class="simple">
<li><p>Evaluate the prediction performance of these models on the test set by <span class="math notranslate nohighlight">\(\log\{L[y_i, \boldsymbol{X}_{i, \ast}; \boldsymbol{\beta}_{-i}(\lambda), \boldsymbol{\sigma}_{-i}^2(\lambda)]\}\)</span>. Or, by the prediction error <span class="math notranslate nohighlight">\(|y_i - \boldsymbol{X}_{i, \ast} \boldsymbol{\beta}_{-i}(\lambda)|\)</span>, the relative error, the error squared or the R2 score function.</p></li>
<li><p>Evaluate the prediction performance of these models on the test set by <span class="math notranslate nohighlight">\(\log\{L[y_i, \boldsymbol{X}_{i, \ast}; \boldsymbol{\theta}_{-i}(\lambda), \boldsymbol{\sigma}_{-i}^2(\lambda)]\}\)</span>. Or, by the prediction error <span class="math notranslate nohighlight">\(|y_i - \boldsymbol{X}_{i, \ast} \boldsymbol{\theta}_{-i}(\lambda)|\)</span>, the relative error, the error squared or the R2 score function.</p></li>
<li><p>Repeat the first three steps such that each sample plays the role of the test set once.</p></li>
<li><p>Average the prediction performances of the test sets at each grid point of the penalty bias/parameter. It is an estimate of the prediction performance of the model corresponding to this value of the penalty parameter on novel data. It is defined as</p></li>
</ul>
<div class="math notranslate nohighlight">
\[
\begin{align*}
\frac{1}{n} \sum_{i = 1}^n \log\{L[y_i, \mathbf{X}_{i, \ast}; \boldsymbol{\beta}_{-i}(\lambda), \boldsymbol{\sigma}_{-i}^2(\lambda)]\}.
\frac{1}{n} \sum_{i = 1}^n \log\{L[y_i, \mathbf{X}_{i, \ast}; \boldsymbol{\theta}_{-i}(\lambda), \boldsymbol{\sigma}_{-i}^2(\lambda)]\}.
\end{align*}
\]</div>
<p>For the various values of <span class="math notranslate nohighlight">\(k\)</span></p>
@@ -1483,8 +1483,8 @@ Thus, if we cannot assume that the expected outputs/targets are zero
when all predictors are zero (the columns in the design matrix), it
may be a bad idea to implement a model which penalizes the intercept.
Furthermore, in for example Ridge and Lasso regression, the default solutions
from the library <strong>Scikit-Learn</strong> (when not shrinking <span class="math notranslate nohighlight">\(\beta_0\)</span>) for the unknown parameters
<span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span>, are derived under the assumption that both <span class="math notranslate nohighlight">\(\boldsymbol{y}\)</span> and
from the library <strong>Scikit-Learn</strong> (when not shrinking <span class="math notranslate nohighlight">\(\theta_0\)</span>) for the unknown parameters
<span class="math notranslate nohighlight">\(\boldsymbol{\theta}\)</span>, are derived under the assumption that both <span class="math notranslate nohighlight">\(\boldsymbol{y}\)</span> and
<span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span> are zero centered, that is we subtract the mean values.</p>
<p>If our predictors represent different scales, then it is important to
standardize the design matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span> by subtracting the mean of each
@@ -1537,39 +1537,39 @@ simplicity, we will focus on ordinary regression, as done in the above example.
<p>The cost/loss function for regression is</p>
<div class="math notranslate nohighlight">
\[
C(\beta_0, \beta_1, ... , \beta_{p-1}) = \frac{1}{n}\sum_{i=0}^{n} \left(y_i - \beta_0 - \sum_{j=1}^{p-1} X_{ij}\beta_j\right)^2,.
C(\theta_0, \theta_1, ... , \theta_{p-1}) = \frac{1}{n}\sum_{i=0}^{n} \left(y_i - \theta_0 - \sum_{j=1}^{p-1} X_{ij}\theta_j\right)^2,.
\]</div>
<p>Recall also that we use the squared value. This expression can lead to an
increased penalty for higher differences between predicted and
output/target values.</p>
<p>What we have done is to single out the <span class="math notranslate nohighlight">\(\beta_0\)</span> term in the
<p>What we have done is to single out the <span class="math notranslate nohighlight">\(\theta_0\)</span> term in the
definition of the mean squared error (MSE). The design matrix <span class="math notranslate nohighlight">\(X\)</span>
does in this case not contain any intercept column. When we take the
derivative with respect to <span class="math notranslate nohighlight">\(\beta_0\)</span>, we want the derivative to obey</p>
derivative with respect to <span class="math notranslate nohighlight">\(\theta_0\)</span>, we want the derivative to obey</p>
<div class="math notranslate nohighlight">
\[
\frac{\partial C}{\partial \beta_j} = 0,
\frac{\partial C}{\partial \theta_j} = 0,
\]</div>
<p>for all <span class="math notranslate nohighlight">\(j\)</span>. For <span class="math notranslate nohighlight">\(\beta_0\)</span> we have</p>
<p>for all <span class="math notranslate nohighlight">\(j\)</span>. For <span class="math notranslate nohighlight">\(\theta_0\)</span> we have</p>
<div class="math notranslate nohighlight">
\[
\frac{\partial C}{\partial \beta_0} = -\frac{2}{n}\sum_{i=0}^{n-1} \left(y_i - \beta_0 - \sum_{j=1}^{p-1} X_{ij} \beta_j\right).
\frac{\partial C}{\partial \theta_0} = -\frac{2}{n}\sum_{i=0}^{n-1} \left(y_i - \theta_0 - \sum_{j=1}^{p-1} X_{ij} \theta_j\right).
\]</div>
<p>Multiplying away the constant <span class="math notranslate nohighlight">\(2/n\)</span>, we obtain</p>
<div class="math notranslate nohighlight">
\[
\sum_{i=0}^{n-1} \beta_0 = \sum_{i=0}^{n-1}y_i - \sum_{i=0}^{n-1} \sum_{j=1}^{p-1} X_{ij} \beta_j.
\sum_{i=0}^{n-1} \theta_0 = \sum_{i=0}^{n-1}y_i - \sum_{i=0}^{n-1} \sum_{j=1}^{p-1} X_{ij} \theta_j.
\]</div>
<p>Let us specialize first to the case where we have only two parameters <span class="math notranslate nohighlight">\(\beta_0\)</span> and <span class="math notranslate nohighlight">\(\beta_1\)</span>.
Our result for <span class="math notranslate nohighlight">\(\beta_0\)</span> simplifies then to</p>
<p>Let us specialize first to the case where we have only two parameters <span class="math notranslate nohighlight">\(\theta_0\)</span> and <span class="math notranslate nohighlight">\(\theta_1\)</span>.
Our result for <span class="math notranslate nohighlight">\(\theta_0\)</span> simplifies then to</p>
<div class="math notranslate nohighlight">
\[
n\beta_0 = \sum_{i=0}^{n-1}y_i - \sum_{i=0}^{n-1} X_{i1} \beta_1.
n\theta_0 = \sum_{i=0}^{n-1}y_i - \sum_{i=0}^{n-1} X_{i1} \theta_1.
\]</div>
<p>We obtain then</p>
<div class="math notranslate nohighlight">
\[
\beta_0 = \frac{1}{n}\sum_{i=0}^{n-1}y_i - \beta_1\frac{1}{n}\sum_{i=0}^{n-1} X_{i1}.
\theta_0 = \frac{1}{n}\sum_{i=0}^{n-1}y_i - \theta_1\frac{1}{n}\sum_{i=0}^{n-1} X_{i1}.
\]</div>
<p>If we define</p>
<div class="math notranslate nohighlight">
@@ -1584,17 +1584,17 @@ n\beta_0 = \sum_{i=0}^{n-1}y_i - \sum_{i=0}^{n-1} X_{i1} \beta_1.
<p>we have</p>
<div class="math notranslate nohighlight">
\[
\beta_0 = \mu_y - \beta_1\mu_{\boldsymbol{x}_1}.
\theta_0 = \mu_y - \theta_1\mu_{\boldsymbol{x}_1}.
\]</div>
<p>In the general case with more parameters than <span class="math notranslate nohighlight">\(\beta_0\)</span> and <span class="math notranslate nohighlight">\(\beta_1\)</span>, we have</p>
<p>In the general case with more parameters than <span class="math notranslate nohighlight">\(\theta_0\)</span> and <span class="math notranslate nohighlight">\(\theta_1\)</span>, we have</p>
<div class="math notranslate nohighlight">
\[
\beta_0 = \frac{1}{n}\sum_{i=0}^{n-1}y_i - \frac{1}{n}\sum_{i=0}^{n-1}\sum_{j=1}^{p-1} X_{ij}\beta_j.
\theta_0 = \frac{1}{n}\sum_{i=0}^{n-1}y_i - \frac{1}{n}\sum_{i=0}^{n-1}\sum_{j=1}^{p-1} X_{ij}\theta_j.
\]</div>
<p>We can rewrite the latter equation as</p>
<div class="math notranslate nohighlight">
\[
\beta_0 = \frac{1}{n}\sum_{i=0}^{n-1}y_i - \sum_{j=1}^{p-1} \mu_{\boldsymbol{x}_j}\beta_j,
\theta_0 = \frac{1}{n}\sum_{i=0}^{n-1}y_i - \sum_{j=1}^{p-1} \mu_{\boldsymbol{x}_j}\theta_j,
\]</div>
<p>where we have defined</p>
<div class="math notranslate nohighlight">
@@ -1605,19 +1605,19 @@ n\beta_0 = \sum_{i=0}^{n-1}y_i - \sum_{i=0}^{n-1} X_{i1} \beta_1.
<p>Replacing <span class="math notranslate nohighlight">\(y_i\)</span> with <span class="math notranslate nohighlight">\(y_i - y_i - \overline{\boldsymbol{y}}\)</span> and centering also our design matrix results in a cost function (in vector-matrix disguise)</p>
<div class="math notranslate nohighlight">
\[
C(\boldsymbol{\beta}) = (\boldsymbol{\tilde{y}} - \tilde{X}\boldsymbol{\beta})^T(\boldsymbol{\tilde{y}} - \tilde{X}\boldsymbol{\beta}).
C(\boldsymbol{\theta}) = (\boldsymbol{\tilde{y}} - \tilde{X}\boldsymbol{\theta})^T(\boldsymbol{\tilde{y}} - \tilde{X}\boldsymbol{\theta}).
\]</div>
<p>If we minimize with respect to <span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span> we have then</p>
<p>If we minimize with respect to <span class="math notranslate nohighlight">\(\boldsymbol{\theta}\)</span> we have then</p>
<div class="math notranslate nohighlight">
\[
\hat{\boldsymbol{\beta}} = (\tilde{X}^T\tilde{X})^{-1}\tilde{X}^T\boldsymbol{\tilde{y}},
\hat{\boldsymbol{\theta}} = (\tilde{X}^T\tilde{X})^{-1}\tilde{X}^T\boldsymbol{\tilde{y}},
\]</div>
<p>where <span class="math notranslate nohighlight">\(\boldsymbol{\tilde{y}} = \boldsymbol{y} - \overline{\boldsymbol{y}}\)</span>
and <span class="math notranslate nohighlight">\(\tilde{X}_{ij} = X_{ij} - \frac{1}{n}\sum_{k=0}^{n-1}X_{kj}\)</span>.</p>
<p>For Ridge regression we need to add <span class="math notranslate nohighlight">\(\lambda \boldsymbol{\beta}^T\boldsymbol{\beta}\)</span> to the cost function and get then</p>
<p>For Ridge regression we need to add <span class="math notranslate nohighlight">\(\lambda \boldsymbol{\theta}^T\boldsymbol{\theta}\)</span> to the cost function and get then</p>
<div class="math notranslate nohighlight">
\[
\hat{\boldsymbol{\beta}} = (\tilde{X}^T\tilde{X} + \lambda I)^{-1}\tilde{X}^T\boldsymbol{\tilde{y}}.
\hat{\boldsymbol{\theta}} = (\tilde{X}^T\tilde{X} + \lambda I)^{-1}\tilde{X}^T\boldsymbol{\tilde{y}}.
\]</div>
<p>What does this mean? And why do we insist on all this? Let us look at some examples.</p>
<p>This code shows a simple first-order fit to a data set using the above transformed data, where we consider the role of the intercept first, by either excluding it or including it (<em>code example thanks to Øyvind Sigmundson Schøyen</em>). Here our scaling of the data is done by subtracting the mean values only.
@@ -1637,15 +1637,15 @@ def MSE(y_data,y_model):
return np.sum((y_data-y_model)**2)/n
def fit_beta(X, y):
def fit_theta(X, y):
return np.linalg.pinv(X.T @ X) @ X.T @ y
true_beta = [2, 0.5, 3.7]
true_theta = [2, 0.5, 3.7]
x = np.linspace(0, 1, 11)
y = np.sum(
np.asarray([x ** p * b for p, b in enumerate(true_beta)]), axis=0
np.asarray([x ** p * b for p, b in enumerate(true_theta)]), axis=0
) + 0.1 * np.random.normal(size=len(x))
degree = 3
@@ -1655,15 +1655,15 @@ X = np.zeros((len(x), degree))
for p in range(degree):
X[:, p] = x ** p
beta = fit_beta(X, y)
theta = fit_theta(X, y)
# Intercept is included in the design matrix
skl = LinearRegression(fit_intercept=False).fit(X, y)
print(f&quot;True beta: {true_beta}&quot;)
print(f&quot;Fitted beta: {beta}&quot;)
print(f&quot;Sklearn fitted beta: {skl.coef_}&quot;)
ypredictOwn = X @ beta
print(f&quot;True theta: {true_theta}&quot;)
print(f&quot;Fitted theta: {theta}&quot;)
print(f&quot;Sklearn fitted theta: {skl.coef_}&quot;)
ypredictOwn = X @ theta
ypredictSKL = skl.predict(X)
print(f&quot;MSE with intercept column&quot;)
print(MSE(y,ypredictOwn))
@@ -1673,7 +1673,7 @@ print(MSE(y,ypredictSKL))
plt.figure()
plt.scatter(x, y, label=&quot;Data&quot;)
plt.plot(x, X @ beta, label=&quot;Fit&quot;)
plt.plot(x, X @ theta, label=&quot;Fit&quot;)
plt.plot(x, skl.predict(X), label=&quot;Sklearn (fit_intercept=False)&quot;)
@@ -1690,21 +1690,21 @@ skl = LinearRegression(fit_intercept=True).fit(X, y)
y_offset = np.average(y, axis=0)
X_offset = np.average(X, axis=0)
beta = fit_beta(X - X_offset, y - y_offset)
intercept = np.mean(y_offset - X_offset @ beta)
theta = fit_theta(X - X_offset, y - y_offset)
intercept = np.mean(y_offset - X_offset @ theta)
print(f&quot;Manual intercept: {intercept}&quot;)
print(f&quot;Fitted beta (wiothout intercept): {beta}&quot;)
print(f&quot;Fitted theta (wiothout intercept): {theta}&quot;)
print(f&quot;Sklearn intercept: {skl.intercept_}&quot;)
print(f&quot;Sklearn fitted beta (without intercept): {skl.coef_}&quot;)
ypredictOwn = X @ beta
print(f&quot;Sklearn fitted theta (without intercept): {skl.coef_}&quot;)
ypredictOwn = X @ theta
ypredictSKL = skl.predict(X)
print(f&quot;MSE with Manual intercept&quot;)
print(MSE(y,ypredictOwn+intercept))
print(f&quot;MSE with Sklearn intercept&quot;)
print(MSE(y,ypredictSKL))
plt.plot(x, X @ beta + intercept, &quot;--&quot;, label=&quot;Fit (manual intercept)&quot;)
plt.plot(x, X @ theta + intercept, &quot;--&quot;, label=&quot;Fit (manual intercept)&quot;)
plt.plot(x, skl.predict(X), &quot;--&quot;, label=&quot;Sklearn (fit_intercept=True)&quot;)
plt.grid()
plt.legend()
@@ -1721,22 +1721,22 @@ they should. However, when we move to for example Ridge regression,
the way we treat the intercept may give a larger or smaller MSE,
meaning that the MSE can be penalized by the value of the
intercept. Not including the intercept in the fit, means that the
regularization term does not include <span class="math notranslate nohighlight">\(\beta_0\)</span>. For different values
regularization term does not include <span class="math notranslate nohighlight">\(\theta_0\)</span>. For different values
of <span class="math notranslate nohighlight">\(\lambda\)</span>, this may lead to different MSE values.</p>
<p>To remind the reader, the regularization term, with the intercept in Ridge regression, is given by</p>
<div class="math notranslate nohighlight">
\[
\lambda \vert\vert \boldsymbol{\beta} \vert\vert_2^2 = \lambda \sum_{j=0}^{p-1}\beta_j^2,
\lambda \vert\vert \boldsymbol{\theta} \vert\vert_2^2 = \lambda \sum_{j=0}^{p-1}\theta_j^2,
\]</div>
<p>but when we take out the intercept, this equation becomes</p>
<div class="math notranslate nohighlight">
\[
\lambda \vert\vert \boldsymbol{\beta} \vert\vert_2^2 = \lambda \sum_{j=1}^{p-1}\beta_j^2.
\lambda \vert\vert \boldsymbol{\theta} \vert\vert_2^2 = \lambda \sum_{j=1}^{p-1}\theta_j^2.
\]</div>
<p>For Lasso regression we have</p>
<div class="math notranslate nohighlight">
\[
\lambda \vert\vert \boldsymbol{\beta} \vert\vert_1 = \lambda \sum_{j=1}^{p-1}\vert\beta_j\vert.
\lambda \vert\vert \boldsymbol{\theta} \vert\vert_1 = \lambda \sum_{j=1}^{p-1}\vert\theta_j\vert.
\]</div>
<p>It means that, when scaling the design matrix and the outputs/targets,
by subtracting the mean values, we have an optimization problem which
@@ -1783,20 +1783,20 @@ MSERidgePredict = np.zeros(nlambdas)
lambdas = np.logspace(-4, 2, nlambdas)
for i in range(nlambdas):
lmb = lambdas[i]
OwnRidgeBeta = np.linalg.pinv(X_train.T @ X_train+lmb*I) @ X_train.T @ y_train
OwnRidgeTheta = np.linalg.pinv(X_train.T @ X_train+lmb*I) @ X_train.T @ y_train
# Note: we include the intercept column and no scaling
RegRidge = linear_model.Ridge(lmb,fit_intercept=False)
RegRidge.fit(X_train,y_train)
# and then make the prediction
ytildeOwnRidge = X_train @ OwnRidgeBeta
ypredictOwnRidge = X_test @ OwnRidgeBeta
ytildeOwnRidge = X_train @ OwnRidgeTheta
ypredictOwnRidge = X_test @ OwnRidgeTheta
ytildeRidge = RegRidge.predict(X_train)
ypredictRidge = RegRidge.predict(X_test)
MSEOwnRidgePredict[i] = MSE(y_test,ypredictOwnRidge)
MSERidgePredict[i] = MSE(y_test,ypredictRidge)
print(&quot;Beta values for own Ridge implementation&quot;)
print(OwnRidgeBeta)
print(&quot;Beta values for Scikit-Learn Ridge implementation&quot;)
print(&quot;Theta values for own Ridge implementation&quot;)
print(OwnRidgeTheta)
print(&quot;Theta values for Scikit-Learn Ridge implementation&quot;)
print(RegRidge.coef_)
print(&quot;MSE values for own Ridge implementation&quot;)
print(MSEOwnRidgePredict[i])
@@ -1869,18 +1869,18 @@ MSERidgePredict = np.zeros(nlambdas)
lambdas = np.logspace(-4, 2, nlambdas)
for i in range(nlambdas):
lmb = lambdas[i]
OwnRidgeBeta = np.linalg.pinv(X_train_scaled.T @ X_train_scaled+lmb*I) @ X_train_scaled.T @ (y_train_scaled)
intercept_ = y_scaler - X_train_mean@OwnRidgeBeta #The intercept can be shifted so the model can predict on uncentered data
OwnRidgeTheta = np.linalg.pinv(X_train_scaled.T @ X_train_scaled+lmb*I) @ X_train_scaled.T @ (y_train_scaled)
intercept_ = y_scaler - X_train_mean@OwnRidgeTheta #The intercept can be shifted so the model can predict on uncentered data
#Add intercept to prediction
ypredictOwnRidge = X_test_scaled @ OwnRidgeBeta + y_scaler
ypredictOwnRidge = X_test_scaled @ OwnRidgeTheta + y_scaler
RegRidge = linear_model.Ridge(lmb)
RegRidge.fit(X_train,y_train)
ypredictRidge = RegRidge.predict(X_test)
MSEOwnRidgePredict[i] = MSE(y_test,ypredictOwnRidge)
MSERidgePredict[i] = MSE(y_test,ypredictRidge)
print(&quot;Beta values for own Ridge implementation&quot;)
print(OwnRidgeBeta) #Intercept is given by mean of target variable
print(&quot;Beta values for Scikit-Learn Ridge implementation&quot;)
print(&quot;Theta values for own Ridge implementation&quot;)
print(OwnRidgeTheta) #Intercept is given by mean of target variable
print(&quot;Theta values for Scikit-Learn Ridge implementation&quot;)
print(RegRidge.coef_)
print(&#39;Intercept from own implementation:&#39;)
print(intercept_)
@@ -1907,7 +1907,7 @@ plt.show()
<p>We see here, when compared to the code which includes explicitely the
intercept column, that our MSE value is actually smaller. This is
because the regularization term does not include the intercept value
<span class="math notranslate nohighlight">\(\beta_0\)</span> in the fitting. This applies to Lasso regularization as
<span class="math notranslate nohighlight">\(\theta_0\)</span> in the fitting. This applies to Lasso regularization as
well. It means that our optimization is now done only with the
centered matrix and/or vector that enter the fitting procedure. Note
also that the problem with the intercept occurs mainly in these type
@@ -1994,7 +1994,7 @@ with the form utilized in linear regression, that is</p>
<div class="math notranslate nohighlight">
\[
\begin{equation}
\boldsymbol{y} = \boldsymbol{X}\boldsymbol{\beta} + \boldsymbol{\epsilon},
\boldsymbol{y} = \boldsymbol{X}\boldsymbol{\theta} + \boldsymbol{\epsilon},
\label{_auto4} \tag{4}
\end{equation}
\]</div>
@@ -2016,15 +2016,15 @@ X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)
<div class="math notranslate nohighlight">
\[
\begin{equation}
C(\boldsymbol{X}, \boldsymbol{\beta})= \frac{1}{n}\left\{(\boldsymbol{X}\boldsymbol{\beta} - \boldsymbol{y})^T(\boldsymbol{X}\boldsymbol{\beta} - \boldsymbol{y})\right\}.
C(\boldsymbol{X}, \boldsymbol{\theta})= \frac{1}{n}\left\{(\boldsymbol{X}\boldsymbol{\theta} - \boldsymbol{y})^T(\boldsymbol{X}\boldsymbol{\theta} - \boldsymbol{y})\right\}.
\label{_auto5} \tag{5}
\end{equation}
\]</div>
<p>We then find the extremal point of <span class="math notranslate nohighlight">\(C\)</span> by taking the derivative with respect to <span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span> as discussed above.
This yields the expression for <span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span> to be</p>
<p>We then find the extremal point of <span class="math notranslate nohighlight">\(C\)</span> by taking the derivative with respect to <span class="math notranslate nohighlight">\(\boldsymbol{\theta}\)</span> as discussed above.
This yields the expression for <span class="math notranslate nohighlight">\(\boldsymbol{\theta}\)</span> to be</p>
<div class="math notranslate nohighlight">
\[
\boldsymbol{\beta} = \frac{\boldsymbol{X}^T \boldsymbol{y}}{\boldsymbol{X}^T \boldsymbol{X}},
\boldsymbol{\theta} = \frac{\boldsymbol{X}^T \boldsymbol{y}}{\boldsymbol{X}^T \boldsymbol{X}},
\]</div>
<p>which immediately imposes some requirements on <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span> as there must exist
an inverse of <span class="math notranslate nohighlight">\(\boldsymbol{X}^T \boldsymbol{X}\)</span>. If the expression we are modeling contains an
@@ -2047,10 +2047,10 @@ X_test_own = np.concatenate(
<p>Doing the inversion directly turns out to be a bad idea since the matrix
<span class="math notranslate nohighlight">\(\boldsymbol{X}^T\boldsymbol{X}\)</span> is singular. An alternative approach is to use the <strong>singular
value decomposition</strong>. Using the definition of the Moore-Penrose
pseudoinverse we can write the equation for <span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span> as</p>
pseudoinverse we can write the equation for <span class="math notranslate nohighlight">\(\boldsymbol{\theta}\)</span> as</p>
<div class="math notranslate nohighlight">
\[
\boldsymbol{\beta} = \boldsymbol{X}^{+}\boldsymbol{y},
\boldsymbol{\theta} = \boldsymbol{X}^{+}\boldsymbol{y},
\]</div>
<p>where the pseudoinverse of <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span> is given by</p>
<div class="math notranslate nohighlight">
@@ -2066,7 +2066,7 @@ where <span class="math notranslate nohighlight">\(X^{+} = V\Sigma^{+} U^T\)</sp
<div class="math notranslate nohighlight">
\[
\begin{equation}
\boldsymbol{\beta} = \boldsymbol{V}\boldsymbol{\Sigma}^{+} \boldsymbol{U}^T \boldsymbol{y}.
\boldsymbol{\theta} = \boldsymbol{V}\boldsymbol{\Sigma}^{+} \boldsymbol{U}^T \boldsymbol{y}.
\label{_auto6} \tag{6}
\end{equation}
\]</div>
@@ -2087,7 +2087,7 @@ linear system as an equation would reduce this down to
</div>
<div class="cell docutils container">
<div class="cell_input docutils container">
<div class="highlight-none notranslate"><div class="highlight"><pre><span></span>beta = ols_svd(X_train_own,y_train)
<div class="highlight-none notranslate"><div class="highlight"><pre><span></span>theta = ols_svd(X_train_own,y_train)
</pre></div>
</div>
</div>
@@ -2095,7 +2095,7 @@ linear system as an equation would reduce this down to
<p>When extracting the <span class="math notranslate nohighlight">\(J\)</span>-matrix we need to make sure that we remove the intercept, as is done here</p>
<div class="cell docutils container">
<div class="cell_input docutils container">
<div class="highlight-none notranslate"><div class="highlight"><pre><span></span>J = beta[1:].reshape(L, L)
<div class="highlight-none notranslate"><div class="highlight"><pre><span></span>J = theta[1:].reshape(L, L)
</pre></div>
</div>
</div>
@@ -2195,7 +2195,7 @@ with the form utilized in linear regression, viz.</p>
<div class="math notranslate nohighlight">
\[
\begin{equation}
\boldsymbol{y} = \boldsymbol{X}\boldsymbol{\beta} + \boldsymbol{\epsilon}.
\boldsymbol{y} = \boldsymbol{X}\boldsymbol{\theta} + \boldsymbol{\epsilon}.
\label{_auto10} \tag{10}
\end{equation}
\]</div>
@@ -2256,14 +2256,14 @@ plt.show()
<p>Having explored the ordinary least squares we move on to ridge
regression. In ridge regression we include a <strong>regularizer</strong>. This
involves a new cost function which leads to a new estimate for the
weights <span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span>. This results in a penalized regression problem. The
weights <span class="math notranslate nohighlight">\(\boldsymbol{\theta}\)</span>. This results in a penalized regression problem. The
cost function is given by</p>
<!-- Equation labels as ordinary links -->
<div id="_auto11"></div>
<div class="math notranslate nohighlight">
\[
\begin{equation}
C(\boldsymbol{X}, \boldsymbol{\beta}; \lambda) = (\boldsymbol{X}\boldsymbol{\beta} - \boldsymbol{y})^T(\boldsymbol{X}\boldsymbol{\beta} - \boldsymbol{y}) + \lambda \boldsymbol{\beta}^T\boldsymbol{\beta}.
C(\boldsymbol{X}, \boldsymbol{\theta}; \lambda) = (\boldsymbol{X}\boldsymbol{\theta} - \boldsymbol{y})^T(\boldsymbol{X}\boldsymbol{\theta} - \boldsymbol{y}) + \lambda \boldsymbol{\theta}^T\boldsymbol{\theta}.
\label{_auto11} \tag{11}
\end{equation}
\]</div>
@@ -2291,7 +2291,7 @@ plt.show()
<div class="math notranslate nohighlight">
\[
\begin{equation}
C(\boldsymbol{X}, \boldsymbol{\beta}; \lambda) = (\boldsymbol{X}\boldsymbol{\beta} - \boldsymbol{y})^T(\boldsymbol{X}\boldsymbol{\beta} - \boldsymbol{y}) + \lambda \sqrt{\boldsymbol{\beta}^T\boldsymbol{\beta}}.
C(\boldsymbol{X}, \boldsymbol{\theta}; \lambda) = (\boldsymbol{X}\boldsymbol{\theta} - \boldsymbol{y})^T(\boldsymbol{X}\boldsymbol{\theta} - \boldsymbol{y}) + \lambda \sqrt{\boldsymbol{\theta}^T\boldsymbol{\theta}}.
\label{_auto12} \tag{12}
\end{equation}
\]</div>
@@ -2505,7 +2505,7 @@ distribution <span class="math notranslate nohighlight">\(N(0,1)\)</span>.</p>
value decomposition from e.g., <strong>numpy</strong> ) or use your code from
homeworks 1 and 2 and perform a standard least square regression
analysis using polynomials in <span class="math notranslate nohighlight">\(x\)</span> and <span class="math notranslate nohighlight">\(y\)</span> up to fifth order. Find the
<a class="reference external" href="https://en.wikipedia.org/wiki/Confidence_interval">confidence intervals</a> of the parameters (estimators) <span class="math notranslate nohighlight">\(\beta\)</span> by computing their
<a class="reference external" href="https://en.wikipedia.org/wiki/Confidence_interval">confidence intervals</a> of the parameters (estimators) <span class="math notranslate nohighlight">\(\theta\)</span> by computing their
variances, evaluate the Mean Squared error (MSE)</p>
<div class="math notranslate nohighlight">
\[
@@ -2565,13 +2565,13 @@ dataset <span class="math notranslate nohighlight">\(\mathcal{L}\)</span> consis
deviation <span class="math notranslate nohighlight">\(\sigma^2\)</span>.</p>
<p>In our derivation of the ordinary least squares method we defined then
an approximation to the function <span class="math notranslate nohighlight">\(f\)</span> in terms of the parameters
<span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span> and the design matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span> which embody our model,
that is <span class="math notranslate nohighlight">\(\boldsymbol{\tilde{y}}=\boldsymbol{X}\boldsymbol{\beta}\)</span>.</p>
<p>The parameters <span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span> are in turn found by optimizing the means
<span class="math notranslate nohighlight">\(\boldsymbol{\theta}\)</span> and the design matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span> which embody our model,
that is <span class="math notranslate nohighlight">\(\boldsymbol{\tilde{y}}=\boldsymbol{X}\boldsymbol{\theta}\)</span>.</p>
<p>The parameters <span class="math notranslate nohighlight">\(\boldsymbol{\theta}\)</span> are in turn found by optimizing the means
squared error via the so-called cost function</p>
<div class="math notranslate nohighlight">
\[
C(\boldsymbol{X},\boldsymbol{\beta}) =\frac{1}{n}\sum_{i=0}^{n-1}(y_i-\tilde{y}_i)^2=\mathbb{E}\left[(\boldsymbol{y}-\boldsymbol{\tilde{y}})^2\right].
C(\boldsymbol{X},\boldsymbol{\theta}) =\frac{1}{n}\sum_{i=0}^{n-1}(y_i-\tilde{y}_i)^2=\mathbb{E}\left[(\boldsymbol{y}-\boldsymbol{\tilde{y}})^2\right].
\]</div>
<p>Here the expected value <span class="math notranslate nohighlight">\(\mathbb{E}\)</span> is the sample value.</p>
<p>Show that you can rewrite this as</p>