corrected beta--> theta
This commit is contained in:
Morten Hjorth-Jensen
@@ -348,8 +348,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 <b>Scikit-Learn</b> (when not shrinking \( \beta_0 \)) for the unknown parameters
|
||||
\( \boldsymbol{\beta} \), are derived under the assumption that both \( \boldsymbol{y} \) and
|
||||
from the library <b>Scikit-Learn</b> (when not shrinking \( \theta_0 \)) for the unknown parameters
|
||||
\( \boldsymbol{\theta} \), are derived under the assumption that both \( \boldsymbol{y} \) and
|
||||
\( \boldsymbol{X} \) are zero centered, that is we subtract the mean values.
|
||||
</p>
|
||||
|
||||
@@ -433,7 +433,7 @@ simplicity, we will focus on ordinary regression, as done in the above example.
|
||||
|
||||
<p>The cost/loss function for regression is</p>
|
||||
$$
|
||||
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,.
|
||||
$$
|
||||
|
||||
<p>Recall also that we use the squared value. This expression can lead to an
|
||||
@@ -441,37 +441,37 @@ increased penalty for higher differences between predicted and
|
||||
output/target values.
|
||||
</p>
|
||||
|
||||
<p>What we have done is to single out the \( \beta_0 \) term in the
|
||||
<p>What we have done is to single out the \( \theta_0 \) term in the
|
||||
definition of the mean squared error (MSE). The design matrix \( X \)
|
||||
does in this case not contain any intercept column. When we take the
|
||||
derivative with respect to \( \beta_0 \), we want the derivative to obey
|
||||
derivative with respect to \( \theta_0 \), we want the derivative to obey
|
||||
</p>
|
||||
|
||||
$$
|
||||
\frac{\partial C}{\partial \beta_j} = 0,
|
||||
\frac{\partial C}{\partial \theta_j} = 0,
|
||||
$$
|
||||
|
||||
<p>for all \( j \). For \( \beta_0 \) we have</p>
|
||||
<p>for all \( j \). For \( \theta_0 \) we have</p>
|
||||
|
||||
$$
|
||||
\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).
|
||||
$$
|
||||
|
||||
<p>Multiplying away the constant \( 2/n \), we obtain</p>
|
||||
$$
|
||||
\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.
|
||||
$$
|
||||
|
||||
<p>Let us specialize first to the case where we have only two parameters \( \beta_0 \) and \( \beta_1 \).
|
||||
Our result for \( \beta_0 \) simplifies then to
|
||||
<p>Let us specialize first to the case where we have only two parameters \( \theta_0 \) and \( \theta_1 \).
|
||||
Our result for \( \theta_0 \) simplifies then to
|
||||
</p>
|
||||
$$
|
||||
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.
|
||||
$$
|
||||
|
||||
<p>We obtain then</p>
|
||||
$$
|
||||
\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}.
|
||||
$$
|
||||
|
||||
<p>If we define</p>
|
||||
@@ -486,17 +486,17 @@ $$
|
||||
|
||||
<p>we have</p>
|
||||
$$
|
||||
\beta_0 = \mu_y - \beta_1\mu_{\boldsymbol{x}_1}.
|
||||
\theta_0 = \mu_y - \theta_1\mu_{\boldsymbol{x}_1}.
|
||||
$$
|
||||
|
||||
<p>In the general case with more parameters than \( \beta_0 \) and \( \beta_1 \), we have</p>
|
||||
<p>In the general case with more parameters than \( \theta_0 \) and \( \theta_1 \), we have</p>
|
||||
$$
|
||||
\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.
|
||||
$$
|
||||
|
||||
<p>We can rewrite the latter equation as</p>
|
||||
$$
|
||||
\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,
|
||||
$$
|
||||
|
||||
<p>where we have defined</p>
|
||||
@@ -508,22 +508,22 @@ $$
|
||||
|
||||
<p>Replacing \( y_i \) with \( y_i - y_i - \overline{\boldsymbol{y}} \) and centering also our design matrix results in a cost function (in vector-matrix disguise)</p>
|
||||
$$
|
||||
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}).
|
||||
$$
|
||||
|
||||
<p>If we minimize with respect to \( \boldsymbol{\beta} \) we have then</p>
|
||||
<p>If we minimize with respect to \( \boldsymbol{\theta} \) we have then</p>
|
||||
|
||||
$$
|
||||
\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}},
|
||||
$$
|
||||
|
||||
<p>where \( \boldsymbol{\tilde{y}} = \boldsymbol{y} - \overline{\boldsymbol{y}} \)
|
||||
and \( \tilde{X}_{ij} = X_{ij} - \frac{1}{n}\sum_{k=0}^{n-1}X_{kj} \).
|
||||
</p>
|
||||
|
||||
<p>For Ridge regression we need to add \( \lambda \boldsymbol{\beta}^T\boldsymbol{\beta} \) to the cost function and get then</p>
|
||||
<p>For Ridge regression we need to add \( \lambda \boldsymbol{\theta}^T\boldsymbol{\theta} \) to the cost function and get then</p>
|
||||
$$
|
||||
\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}}.
|
||||
$$
|
||||
|
||||
<p>What does this mean? And why do we insist on all this? Let us look at some examples.</p>
|
||||
@@ -552,15 +552,15 @@ np<span style="color: #666666">.</span>random<span style="color: #666666">.</spa
|
||||
<span style="color: #008000; font-weight: bold">return</span> np<span style="color: #666666">.</span>sum((y_data<span style="color: #666666">-</span>y_model)<span style="color: #666666">**2</span>)<span style="color: #666666">/</span>n
|
||||
|
||||
|
||||
<span style="color: #008000; font-weight: bold">def</span> <span style="color: #0000FF">fit_beta</span>(X, y):
|
||||
<span style="color: #008000; font-weight: bold">def</span> <span style="color: #0000FF">fit_theta</span>(X, y):
|
||||
<span style="color: #008000; font-weight: bold">return</span> np<span style="color: #666666">.</span>linalg<span style="color: #666666">.</span>pinv(X<span style="color: #666666">.</span>T <span style="color: #666666">@</span> X) <span style="color: #666666">@</span> X<span style="color: #666666">.</span>T <span style="color: #666666">@</span> y
|
||||
|
||||
|
||||
true_beta <span style="color: #666666">=</span> [<span style="color: #666666">2</span>, <span style="color: #666666">0.5</span>, <span style="color: #666666">3.7</span>]
|
||||
true_theta <span style="color: #666666">=</span> [<span style="color: #666666">2</span>, <span style="color: #666666">0.5</span>, <span style="color: #666666">3.7</span>]
|
||||
|
||||
x <span style="color: #666666">=</span> np<span style="color: #666666">.</span>linspace(<span style="color: #666666">0</span>, <span style="color: #666666">1</span>, <span style="color: #666666">11</span>)
|
||||
y <span style="color: #666666">=</span> np<span style="color: #666666">.</span>sum(
|
||||
np<span style="color: #666666">.</span>asarray([x <span style="color: #666666">**</span> p <span style="color: #666666">*</span> b <span style="color: #008000; font-weight: bold">for</span> p, b <span style="color: #AA22FF; font-weight: bold">in</span> <span style="color: #008000">enumerate</span>(true_beta)]), axis<span style="color: #666666">=0</span>
|
||||
np<span style="color: #666666">.</span>asarray([x <span style="color: #666666">**</span> p <span style="color: #666666">*</span> b <span style="color: #008000; font-weight: bold">for</span> p, b <span style="color: #AA22FF; font-weight: bold">in</span> <span style="color: #008000">enumerate</span>(true_theta)]), axis<span style="color: #666666">=0</span>
|
||||
) <span style="color: #666666">+</span> <span style="color: #666666">0.1</span> <span style="color: #666666">*</span> np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>normal(size<span style="color: #666666">=</span><span style="color: #008000">len</span>(x))
|
||||
|
||||
degree <span style="color: #666666">=</span> <span style="color: #666666">3</span>
|
||||
@@ -570,15 +570,15 @@ X <span style="color: #666666">=</span> np<span style="color: #666666">.</span>z
|
||||
<span style="color: #008000; font-weight: bold">for</span> p <span style="color: #AA22FF; font-weight: bold">in</span> <span style="color: #008000">range</span>(degree):
|
||||
X[:, p] <span style="color: #666666">=</span> x <span style="color: #666666">**</span> p
|
||||
|
||||
beta <span style="color: #666666">=</span> fit_beta(X, y)
|
||||
theta <span style="color: #666666">=</span> fit_theta(X, y)
|
||||
|
||||
<span style="color: #408080; font-style: italic"># Intercept is included in the design matrix</span>
|
||||
skl <span style="color: #666666">=</span> LinearRegression(fit_intercept<span style="color: #666666">=</span><span style="color: #008000; font-weight: bold">False</span>)<span style="color: #666666">.</span>fit(X, y)
|
||||
|
||||
<span style="color: #008000">print</span>(<span style="color: #BA2121">f"True beta: </span><span style="color: #BB6688; font-weight: bold">{</span>true_beta<span style="color: #BB6688; font-weight: bold">}</span><span style="color: #BA2121">"</span>)
|
||||
<span style="color: #008000">print</span>(<span style="color: #BA2121">f"Fitted beta: </span><span style="color: #BB6688; font-weight: bold">{</span>beta<span style="color: #BB6688; font-weight: bold">}</span><span style="color: #BA2121">"</span>)
|
||||
<span style="color: #008000">print</span>(<span style="color: #BA2121">f"Sklearn fitted beta: </span><span style="color: #BB6688; font-weight: bold">{</span>skl<span style="color: #666666">.</span>coef_<span style="color: #BB6688; font-weight: bold">}</span><span style="color: #BA2121">"</span>)
|
||||
ypredictOwn <span style="color: #666666">=</span> X <span style="color: #666666">@</span> beta
|
||||
<span style="color: #008000">print</span>(<span style="color: #BA2121">f"True theta: </span><span style="color: #BB6688; font-weight: bold">{</span>true_theta<span style="color: #BB6688; font-weight: bold">}</span><span style="color: #BA2121">"</span>)
|
||||
<span style="color: #008000">print</span>(<span style="color: #BA2121">f"Fitted theta: </span><span style="color: #BB6688; font-weight: bold">{</span>theta<span style="color: #BB6688; font-weight: bold">}</span><span style="color: #BA2121">"</span>)
|
||||
<span style="color: #008000">print</span>(<span style="color: #BA2121">f"Sklearn fitted theta: </span><span style="color: #BB6688; font-weight: bold">{</span>skl<span style="color: #666666">.</span>coef_<span style="color: #BB6688; font-weight: bold">}</span><span style="color: #BA2121">"</span>)
|
||||
ypredictOwn <span style="color: #666666">=</span> X <span style="color: #666666">@</span> theta
|
||||
ypredictSKL <span style="color: #666666">=</span> skl<span style="color: #666666">.</span>predict(X)
|
||||
<span style="color: #008000">print</span>(<span style="color: #BA2121">f"MSE with intercept column"</span>)
|
||||
<span style="color: #008000">print</span>(MSE(y,ypredictOwn))
|
||||
@@ -588,7 +588,7 @@ ypredictSKL <span style="color: #666666">=</span> skl<span style="color: #666666
|
||||
|
||||
plt<span style="color: #666666">.</span>figure()
|
||||
plt<span style="color: #666666">.</span>scatter(x, y, label<span style="color: #666666">=</span><span style="color: #BA2121">"Data"</span>)
|
||||
plt<span style="color: #666666">.</span>plot(x, X <span style="color: #666666">@</span> beta, label<span style="color: #666666">=</span><span style="color: #BA2121">"Fit"</span>)
|
||||
plt<span style="color: #666666">.</span>plot(x, X <span style="color: #666666">@</span> theta, label<span style="color: #666666">=</span><span style="color: #BA2121">"Fit"</span>)
|
||||
plt<span style="color: #666666">.</span>plot(x, skl<span style="color: #666666">.</span>predict(X), label<span style="color: #666666">=</span><span style="color: #BA2121">"Sklearn (fit_intercept=False)"</span>)
|
||||
|
||||
|
||||
@@ -605,21 +605,21 @@ skl <span style="color: #666666">=</span> LinearRegression(fit_intercept<span st
|
||||
y_offset <span style="color: #666666">=</span> np<span style="color: #666666">.</span>average(y, axis<span style="color: #666666">=0</span>)
|
||||
X_offset <span style="color: #666666">=</span> np<span style="color: #666666">.</span>average(X, axis<span style="color: #666666">=0</span>)
|
||||
|
||||
beta <span style="color: #666666">=</span> fit_beta(X <span style="color: #666666">-</span> X_offset, y <span style="color: #666666">-</span> y_offset)
|
||||
intercept <span style="color: #666666">=</span> np<span style="color: #666666">.</span>mean(y_offset <span style="color: #666666">-</span> X_offset <span style="color: #666666">@</span> beta)
|
||||
theta <span style="color: #666666">=</span> fit_theta(X <span style="color: #666666">-</span> X_offset, y <span style="color: #666666">-</span> y_offset)
|
||||
intercept <span style="color: #666666">=</span> np<span style="color: #666666">.</span>mean(y_offset <span style="color: #666666">-</span> X_offset <span style="color: #666666">@</span> theta)
|
||||
|
||||
<span style="color: #008000">print</span>(<span style="color: #BA2121">f"Manual intercept: </span><span style="color: #BB6688; font-weight: bold">{</span>intercept<span style="color: #BB6688; font-weight: bold">}</span><span style="color: #BA2121">"</span>)
|
||||
<span style="color: #008000">print</span>(<span style="color: #BA2121">f"Fitted beta (without intercept): </span><span style="color: #BB6688; font-weight: bold">{</span>beta<span style="color: #BB6688; font-weight: bold">}</span><span style="color: #BA2121">"</span>)
|
||||
<span style="color: #008000">print</span>(<span style="color: #BA2121">f"Fitted theta (without intercept): </span><span style="color: #BB6688; font-weight: bold">{</span>theta<span style="color: #BB6688; font-weight: bold">}</span><span style="color: #BA2121">"</span>)
|
||||
<span style="color: #008000">print</span>(<span style="color: #BA2121">f"Sklearn intercept: </span><span style="color: #BB6688; font-weight: bold">{</span>skl<span style="color: #666666">.</span>intercept_<span style="color: #BB6688; font-weight: bold">}</span><span style="color: #BA2121">"</span>)
|
||||
<span style="color: #008000">print</span>(<span style="color: #BA2121">f"Sklearn fitted beta (without intercept): </span><span style="color: #BB6688; font-weight: bold">{</span>skl<span style="color: #666666">.</span>coef_<span style="color: #BB6688; font-weight: bold">}</span><span style="color: #BA2121">"</span>)
|
||||
ypredictOwn <span style="color: #666666">=</span> X <span style="color: #666666">@</span> beta
|
||||
<span style="color: #008000">print</span>(<span style="color: #BA2121">f"Sklearn fitted theta (without intercept): </span><span style="color: #BB6688; font-weight: bold">{</span>skl<span style="color: #666666">.</span>coef_<span style="color: #BB6688; font-weight: bold">}</span><span style="color: #BA2121">"</span>)
|
||||
ypredictOwn <span style="color: #666666">=</span> X <span style="color: #666666">@</span> theta
|
||||
ypredictSKL <span style="color: #666666">=</span> skl<span style="color: #666666">.</span>predict(X)
|
||||
<span style="color: #008000">print</span>(<span style="color: #BA2121">f"MSE with Manual intercept"</span>)
|
||||
<span style="color: #008000">print</span>(MSE(y,ypredictOwn<span style="color: #666666">+</span>intercept))
|
||||
<span style="color: #008000">print</span>(<span style="color: #BA2121">f"MSE with Sklearn intercept"</span>)
|
||||
<span style="color: #008000">print</span>(MSE(y,ypredictSKL))
|
||||
|
||||
plt<span style="color: #666666">.</span>plot(x, X <span style="color: #666666">@</span> beta <span style="color: #666666">+</span> intercept, <span style="color: #BA2121">"--"</span>, label<span style="color: #666666">=</span><span style="color: #BA2121">"Fit (manual intercept)"</span>)
|
||||
plt<span style="color: #666666">.</span>plot(x, X <span style="color: #666666">@</span> theta <span style="color: #666666">+</span> intercept, <span style="color: #BA2121">"--"</span>, label<span style="color: #666666">=</span><span style="color: #BA2121">"Fit (manual intercept)"</span>)
|
||||
plt<span style="color: #666666">.</span>plot(x, skl<span style="color: #666666">.</span>predict(X), <span style="color: #BA2121">"--"</span>, label<span style="color: #666666">=</span><span style="color: #BA2121">"Sklearn (fit_intercept=True)"</span>)
|
||||
plt<span style="color: #666666">.</span>grid()
|
||||
plt<span style="color: #666666">.</span>legend()
|
||||
@@ -649,23 +649,23 @@ 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 \( \beta_0 \). For different values
|
||||
regularization term does not include \( \theta_0 \). For different values
|
||||
of \( \lambda \), 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>
|
||||
$$
|
||||
\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,
|
||||
$$
|
||||
|
||||
<p>but when we take out the intercept, this equation becomes</p>
|
||||
$$
|
||||
\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.
|
||||
$$
|
||||
|
||||
<p>For Lasso regression we have</p>
|
||||
$$
|
||||
\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.
|
||||
$$
|
||||
|
||||
<p>It means that, when scaling the design matrix and the outputs/targets,
|
||||
@@ -721,20 +721,20 @@ MSERidgePredict <span style="color: #666666">=</span> np<span style="color: #666
|
||||
lambdas <span style="color: #666666">=</span> np<span style="color: #666666">.</span>logspace(<span style="color: #666666">-4</span>, <span style="color: #666666">2</span>, nlambdas)
|
||||
<span style="color: #008000; font-weight: bold">for</span> i <span style="color: #AA22FF; font-weight: bold">in</span> <span style="color: #008000">range</span>(nlambdas):
|
||||
lmb <span style="color: #666666">=</span> lambdas[i]
|
||||
OwnRidgeBeta <span style="color: #666666">=</span> np<span style="color: #666666">.</span>linalg<span style="color: #666666">.</span>pinv(X_train<span style="color: #666666">.</span>T <span style="color: #666666">@</span> X_train<span style="color: #666666">+</span>lmb<span style="color: #666666">*</span>I) <span style="color: #666666">@</span> X_train<span style="color: #666666">.</span>T <span style="color: #666666">@</span> y_train
|
||||
OwnRidgeTheta <span style="color: #666666">=</span> np<span style="color: #666666">.</span>linalg<span style="color: #666666">.</span>pinv(X_train<span style="color: #666666">.</span>T <span style="color: #666666">@</span> X_train<span style="color: #666666">+</span>lmb<span style="color: #666666">*</span>I) <span style="color: #666666">@</span> X_train<span style="color: #666666">.</span>T <span style="color: #666666">@</span> y_train
|
||||
<span style="color: #408080; font-style: italic"># Note: we include the intercept column and no scaling</span>
|
||||
RegRidge <span style="color: #666666">=</span> linear_model<span style="color: #666666">.</span>Ridge(lmb,fit_intercept<span style="color: #666666">=</span><span style="color: #008000; font-weight: bold">False</span>)
|
||||
RegRidge<span style="color: #666666">.</span>fit(X_train,y_train)
|
||||
<span style="color: #408080; font-style: italic"># and then make the prediction</span>
|
||||
ytildeOwnRidge <span style="color: #666666">=</span> X_train <span style="color: #666666">@</span> OwnRidgeBeta
|
||||
ypredictOwnRidge <span style="color: #666666">=</span> X_test <span style="color: #666666">@</span> OwnRidgeBeta
|
||||
ytildeOwnRidge <span style="color: #666666">=</span> X_train <span style="color: #666666">@</span> OwnRidgeTheta
|
||||
ypredictOwnRidge <span style="color: #666666">=</span> X_test <span style="color: #666666">@</span> OwnRidgeTheta
|
||||
ytildeRidge <span style="color: #666666">=</span> RegRidge<span style="color: #666666">.</span>predict(X_train)
|
||||
ypredictRidge <span style="color: #666666">=</span> RegRidge<span style="color: #666666">.</span>predict(X_test)
|
||||
MSEOwnRidgePredict[i] <span style="color: #666666">=</span> MSE(y_test,ypredictOwnRidge)
|
||||
MSERidgePredict[i] <span style="color: #666666">=</span> MSE(y_test,ypredictRidge)
|
||||
<span style="color: #008000">print</span>(<span style="color: #BA2121">"Beta values for own Ridge implementation"</span>)
|
||||
<span style="color: #008000">print</span>(OwnRidgeBeta)
|
||||
<span style="color: #008000">print</span>(<span style="color: #BA2121">"Beta values for Scikit-Learn Ridge implementation"</span>)
|
||||
<span style="color: #008000">print</span>(<span style="color: #BA2121">"Theta values for own Ridge implementation"</span>)
|
||||
<span style="color: #008000">print</span>(OwnRidgeTheta)
|
||||
<span style="color: #008000">print</span>(<span style="color: #BA2121">"Theta values for Scikit-Learn Ridge implementation"</span>)
|
||||
<span style="color: #008000">print</span>(RegRidge<span style="color: #666666">.</span>coef_)
|
||||
<span style="color: #008000">print</span>(<span style="color: #BA2121">"MSE values for own Ridge implementation"</span>)
|
||||
<span style="color: #008000">print</span>(MSEOwnRidgePredict[i])
|
||||
@@ -825,18 +825,18 @@ MSERidgePredict <span style="color: #666666">=</span> np<span style="color: #666
|
||||
lambdas <span style="color: #666666">=</span> np<span style="color: #666666">.</span>logspace(<span style="color: #666666">-4</span>, <span style="color: #666666">2</span>, nlambdas)
|
||||
<span style="color: #008000; font-weight: bold">for</span> i <span style="color: #AA22FF; font-weight: bold">in</span> <span style="color: #008000">range</span>(nlambdas):
|
||||
lmb <span style="color: #666666">=</span> lambdas[i]
|
||||
OwnRidgeBeta <span style="color: #666666">=</span> np<span style="color: #666666">.</span>linalg<span style="color: #666666">.</span>pinv(X_train_scaled<span style="color: #666666">.</span>T <span style="color: #666666">@</span> X_train_scaled<span style="color: #666666">+</span>lmb<span style="color: #666666">*</span>I) <span style="color: #666666">@</span> X_train_scaled<span style="color: #666666">.</span>T <span style="color: #666666">@</span> (y_train_scaled)
|
||||
intercept_ <span style="color: #666666">=</span> y_scaler <span style="color: #666666">-</span> X_train_mean<span style="color: #AA22FF">@OwnRidgeBeta</span> <span style="color: #408080; font-style: italic">#The intercept can be shifted so the model can predict on uncentered data</span>
|
||||
OwnRidgeTheta <span style="color: #666666">=</span> np<span style="color: #666666">.</span>linalg<span style="color: #666666">.</span>pinv(X_train_scaled<span style="color: #666666">.</span>T <span style="color: #666666">@</span> X_train_scaled<span style="color: #666666">+</span>lmb<span style="color: #666666">*</span>I) <span style="color: #666666">@</span> X_train_scaled<span style="color: #666666">.</span>T <span style="color: #666666">@</span> (y_train_scaled)
|
||||
intercept_ <span style="color: #666666">=</span> y_scaler <span style="color: #666666">-</span> X_train_mean<span style="color: #AA22FF">@OwnRidgeTheta</span> <span style="color: #408080; font-style: italic">#The intercept can be shifted so the model can predict on uncentered data</span>
|
||||
<span style="color: #408080; font-style: italic">#Add intercept to prediction</span>
|
||||
ypredictOwnRidge <span style="color: #666666">=</span> X_test_scaled <span style="color: #666666">@</span> OwnRidgeBeta <span style="color: #666666">+</span> y_scaler
|
||||
ypredictOwnRidge <span style="color: #666666">=</span> X_test_scaled <span style="color: #666666">@</span> OwnRidgeTheta <span style="color: #666666">+</span> y_scaler
|
||||
RegRidge <span style="color: #666666">=</span> linear_model<span style="color: #666666">.</span>Ridge(lmb)
|
||||
RegRidge<span style="color: #666666">.</span>fit(X_train,y_train)
|
||||
ypredictRidge <span style="color: #666666">=</span> RegRidge<span style="color: #666666">.</span>predict(X_test)
|
||||
MSEOwnRidgePredict[i] <span style="color: #666666">=</span> MSE(y_test,ypredictOwnRidge)
|
||||
MSERidgePredict[i] <span style="color: #666666">=</span> MSE(y_test,ypredictRidge)
|
||||
<span style="color: #008000">print</span>(<span style="color: #BA2121">"Beta values for own Ridge implementation"</span>)
|
||||
<span style="color: #008000">print</span>(OwnRidgeBeta) <span style="color: #408080; font-style: italic">#Intercept is given by mean of target variable</span>
|
||||
<span style="color: #008000">print</span>(<span style="color: #BA2121">"Beta values for Scikit-Learn Ridge implementation"</span>)
|
||||
<span style="color: #008000">print</span>(<span style="color: #BA2121">"Theta values for own Ridge implementation"</span>)
|
||||
<span style="color: #008000">print</span>(OwnRidgeTheta) <span style="color: #408080; font-style: italic">#Intercept is given by mean of target variable</span>
|
||||
<span style="color: #008000">print</span>(<span style="color: #BA2121">"Theta values for Scikit-Learn Ridge implementation"</span>)
|
||||
<span style="color: #008000">print</span>(RegRidge<span style="color: #666666">.</span>coef_)
|
||||
<span style="color: #008000">print</span>(<span style="color: #BA2121">'Intercept from own implementation:'</span>)
|
||||
<span style="color: #008000">print</span>(intercept_)
|
||||
@@ -874,7 +874,7 @@ plt<span style="color: #666666">.</span>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
|
||||
\( \beta_0 \) in the fitting. This applies to Lasso regularization as
|
||||
\( \theta_0 \) 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.
|
||||
</p>
|
||||
|
||||
@@ -3016,8 +3016,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 <b>Scikit-Learn</b> (when not shrinking \( \beta_0 \)) for the unknown parameters
|
||||
\( \boldsymbol{\beta} \), are derived under the assumption that both \( \boldsymbol{y} \) and
|
||||
from the library <b>Scikit-Learn</b> (when not shrinking \( \theta_0 \)) for the unknown parameters
|
||||
\( \boldsymbol{\theta} \), are derived under the assumption that both \( \boldsymbol{y} \) and
|
||||
\( \boldsymbol{X} \) are zero centered, that is we subtract the mean values.
|
||||
</p>
|
||||
|
||||
@@ -3102,7 +3102,7 @@ simplicity, we will focus on ordinary regression, as done in the above example.
|
||||
<p>The cost/loss function for regression is</p>
|
||||
<p> <br>
|
||||
$$
|
||||
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,.
|
||||
$$
|
||||
<p> <br>
|
||||
|
||||
@@ -3111,46 +3111,46 @@ increased penalty for higher differences between predicted and
|
||||
output/target values.
|
||||
</p>
|
||||
|
||||
<p>What we have done is to single out the \( \beta_0 \) term in the
|
||||
<p>What we have done is to single out the \( \theta_0 \) term in the
|
||||
definition of the mean squared error (MSE). The design matrix \( X \)
|
||||
does in this case not contain any intercept column. When we take the
|
||||
derivative with respect to \( \beta_0 \), we want the derivative to obey
|
||||
derivative with respect to \( \theta_0 \), we want the derivative to obey
|
||||
</p>
|
||||
|
||||
<p> <br>
|
||||
$$
|
||||
\frac{\partial C}{\partial \beta_j} = 0,
|
||||
\frac{\partial C}{\partial \theta_j} = 0,
|
||||
$$
|
||||
<p> <br>
|
||||
|
||||
<p>for all \( j \). For \( \beta_0 \) we have</p>
|
||||
<p>for all \( j \). For \( \theta_0 \) we have</p>
|
||||
|
||||
<p> <br>
|
||||
$$
|
||||
\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).
|
||||
$$
|
||||
<p> <br>
|
||||
|
||||
<p>Multiplying away the constant \( 2/n \), we obtain</p>
|
||||
<p> <br>
|
||||
$$
|
||||
\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.
|
||||
$$
|
||||
<p> <br>
|
||||
|
||||
<p>Let us specialize first to the case where we have only two parameters \( \beta_0 \) and \( \beta_1 \).
|
||||
Our result for \( \beta_0 \) simplifies then to
|
||||
<p>Let us specialize first to the case where we have only two parameters \( \theta_0 \) and \( \theta_1 \).
|
||||
Our result for \( \theta_0 \) simplifies then to
|
||||
</p>
|
||||
<p> <br>
|
||||
$$
|
||||
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.
|
||||
$$
|
||||
<p> <br>
|
||||
|
||||
<p>We obtain then</p>
|
||||
<p> <br>
|
||||
$$
|
||||
\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}.
|
||||
$$
|
||||
<p> <br>
|
||||
|
||||
@@ -3171,21 +3171,21 @@ $$
|
||||
<p>we have</p>
|
||||
<p> <br>
|
||||
$$
|
||||
\beta_0 = \mu_y - \beta_1\mu_{\boldsymbol{x}_1}.
|
||||
\theta_0 = \mu_y - \theta_1\mu_{\boldsymbol{x}_1}.
|
||||
$$
|
||||
<p> <br>
|
||||
|
||||
<p>In the general case with more parameters than \( \beta_0 \) and \( \beta_1 \), we have</p>
|
||||
<p>In the general case with more parameters than \( \theta_0 \) and \( \theta_1 \), we have</p>
|
||||
<p> <br>
|
||||
$$
|
||||
\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.
|
||||
$$
|
||||
<p> <br>
|
||||
|
||||
<p>We can rewrite the latter equation as</p>
|
||||
<p> <br>
|
||||
$$
|
||||
\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,
|
||||
$$
|
||||
<p> <br>
|
||||
|
||||
@@ -3201,15 +3201,15 @@ $$
|
||||
<p>Replacing \( y_i \) with \( y_i - y_i - \overline{\boldsymbol{y}} \) and centering also our design matrix results in a cost function (in vector-matrix disguise)</p>
|
||||
<p> <br>
|
||||
$$
|
||||
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}).
|
||||
$$
|
||||
<p> <br>
|
||||
|
||||
<p>If we minimize with respect to \( \boldsymbol{\beta} \) we have then</p>
|
||||
<p>If we minimize with respect to \( \boldsymbol{\theta} \) we have then</p>
|
||||
|
||||
<p> <br>
|
||||
$$
|
||||
\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}},
|
||||
$$
|
||||
<p> <br>
|
||||
|
||||
@@ -3217,10 +3217,10 @@ $$
|
||||
and \( \tilde{X}_{ij} = X_{ij} - \frac{1}{n}\sum_{k=0}^{n-1}X_{kj} \).
|
||||
</p>
|
||||
|
||||
<p>For Ridge regression we need to add \( \lambda \boldsymbol{\beta}^T\boldsymbol{\beta} \) to the cost function and get then</p>
|
||||
<p>For Ridge regression we need to add \( \lambda \boldsymbol{\theta}^T\boldsymbol{\theta} \) to the cost function and get then</p>
|
||||
<p> <br>
|
||||
$$
|
||||
\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}}.
|
||||
$$
|
||||
<p> <br>
|
||||
|
||||
@@ -3250,15 +3250,15 @@ np.random.seed(<span style="color: #B452CD">2021</span>)
|
||||
<span style="color: #8B008B; font-weight: bold">return</span> np.sum((y_data-y_model)**<span style="color: #B452CD">2</span>)/n
|
||||
|
||||
|
||||
<span style="color: #8B008B; font-weight: bold">def</span> <span style="color: #008b45">fit_beta</span>(X, y):
|
||||
<span style="color: #8B008B; font-weight: bold">def</span> <span style="color: #008b45">fit_theta</span>(X, y):
|
||||
<span style="color: #8B008B; font-weight: bold">return</span> np.linalg.pinv(X.T @ X) @ X.T @ y
|
||||
|
||||
|
||||
true_beta = [<span style="color: #B452CD">2</span>, <span style="color: #B452CD">0.5</span>, <span style="color: #B452CD">3.7</span>]
|
||||
true_theta = [<span style="color: #B452CD">2</span>, <span style="color: #B452CD">0.5</span>, <span style="color: #B452CD">3.7</span>]
|
||||
|
||||
x = np.linspace(<span style="color: #B452CD">0</span>, <span style="color: #B452CD">1</span>, <span style="color: #B452CD">11</span>)
|
||||
y = np.sum(
|
||||
np.asarray([x ** p * b <span style="color: #8B008B; font-weight: bold">for</span> p, b <span style="color: #8B008B">in</span> <span style="color: #658b00">enumerate</span>(true_beta)]), axis=<span style="color: #B452CD">0</span>
|
||||
np.asarray([x ** p * b <span style="color: #8B008B; font-weight: bold">for</span> p, b <span style="color: #8B008B">in</span> <span style="color: #658b00">enumerate</span>(true_theta)]), axis=<span style="color: #B452CD">0</span>
|
||||
) + <span style="color: #B452CD">0.1</span> * np.random.normal(size=<span style="color: #658b00">len</span>(x))
|
||||
|
||||
degree = <span style="color: #B452CD">3</span>
|
||||
@@ -3268,15 +3268,15 @@ X = np.zeros((<span style="color: #658b00">len</span>(x), degree))
|
||||
<span style="color: #8B008B; font-weight: bold">for</span> p <span style="color: #8B008B">in</span> <span style="color: #658b00">range</span>(degree):
|
||||
X[:, p] = x ** p
|
||||
|
||||
beta = fit_beta(X, y)
|
||||
theta = fit_theta(X, y)
|
||||
|
||||
<span style="color: #228B22"># Intercept is included in the design matrix</span>
|
||||
skl = LinearRegression(fit_intercept=<span style="color: #8B008B; font-weight: bold">False</span>).fit(X, y)
|
||||
|
||||
<span style="color: #658b00">print</span>(<span style="color: #CD5555">f"True beta: {</span>true_beta<span style="color: #CD5555">}"</span>)
|
||||
<span style="color: #658b00">print</span>(<span style="color: #CD5555">f"Fitted beta: {</span>beta<span style="color: #CD5555">}"</span>)
|
||||
<span style="color: #658b00">print</span>(<span style="color: #CD5555">f"Sklearn fitted beta: {</span>skl.coef_<span style="color: #CD5555">}"</span>)
|
||||
ypredictOwn = X @ beta
|
||||
<span style="color: #658b00">print</span>(<span style="color: #CD5555">f"True theta: {</span>true_theta<span style="color: #CD5555">}"</span>)
|
||||
<span style="color: #658b00">print</span>(<span style="color: #CD5555">f"Fitted theta: {</span>theta<span style="color: #CD5555">}"</span>)
|
||||
<span style="color: #658b00">print</span>(<span style="color: #CD5555">f"Sklearn fitted theta: {</span>skl.coef_<span style="color: #CD5555">}"</span>)
|
||||
ypredictOwn = X @ theta
|
||||
ypredictSKL = skl.predict(X)
|
||||
<span style="color: #658b00">print</span>(<span style="color: #CD5555">f"MSE with intercept column"</span>)
|
||||
<span style="color: #658b00">print</span>(MSE(y,ypredictOwn))
|
||||
@@ -3286,7 +3286,7 @@ ypredictSKL = skl.predict(X)
|
||||
|
||||
plt.figure()
|
||||
plt.scatter(x, y, label=<span style="color: #CD5555">"Data"</span>)
|
||||
plt.plot(x, X @ beta, label=<span style="color: #CD5555">"Fit"</span>)
|
||||
plt.plot(x, X @ theta, label=<span style="color: #CD5555">"Fit"</span>)
|
||||
plt.plot(x, skl.predict(X), label=<span style="color: #CD5555">"Sklearn (fit_intercept=False)"</span>)
|
||||
|
||||
|
||||
@@ -3303,21 +3303,21 @@ skl = LinearRegression(fit_intercept=<span style="color: #8B008B; font-weight: b
|
||||
y_offset = np.average(y, axis=<span style="color: #B452CD">0</span>)
|
||||
X_offset = np.average(X, axis=<span style="color: #B452CD">0</span>)
|
||||
|
||||
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)
|
||||
|
||||
<span style="color: #658b00">print</span>(<span style="color: #CD5555">f"Manual intercept: {</span>intercept<span style="color: #CD5555">}"</span>)
|
||||
<span style="color: #658b00">print</span>(<span style="color: #CD5555">f"Fitted beta (without intercept): {</span>beta<span style="color: #CD5555">}"</span>)
|
||||
<span style="color: #658b00">print</span>(<span style="color: #CD5555">f"Fitted theta (without intercept): {</span>theta<span style="color: #CD5555">}"</span>)
|
||||
<span style="color: #658b00">print</span>(<span style="color: #CD5555">f"Sklearn intercept: {</span>skl.intercept_<span style="color: #CD5555">}"</span>)
|
||||
<span style="color: #658b00">print</span>(<span style="color: #CD5555">f"Sklearn fitted beta (without intercept): {</span>skl.coef_<span style="color: #CD5555">}"</span>)
|
||||
ypredictOwn = X @ beta
|
||||
<span style="color: #658b00">print</span>(<span style="color: #CD5555">f"Sklearn fitted theta (without intercept): {</span>skl.coef_<span style="color: #CD5555">}"</span>)
|
||||
ypredictOwn = X @ theta
|
||||
ypredictSKL = skl.predict(X)
|
||||
<span style="color: #658b00">print</span>(<span style="color: #CD5555">f"MSE with Manual intercept"</span>)
|
||||
<span style="color: #658b00">print</span>(MSE(y,ypredictOwn+intercept))
|
||||
<span style="color: #658b00">print</span>(<span style="color: #CD5555">f"MSE with Sklearn intercept"</span>)
|
||||
<span style="color: #658b00">print</span>(MSE(y,ypredictSKL))
|
||||
|
||||
plt.plot(x, X @ beta + intercept, <span style="color: #CD5555">"--"</span>, label=<span style="color: #CD5555">"Fit (manual intercept)"</span>)
|
||||
plt.plot(x, X @ theta + intercept, <span style="color: #CD5555">"--"</span>, label=<span style="color: #CD5555">"Fit (manual intercept)"</span>)
|
||||
plt.plot(x, skl.predict(X), <span style="color: #CD5555">"--"</span>, label=<span style="color: #CD5555">"Sklearn (fit_intercept=True)"</span>)
|
||||
plt.grid()
|
||||
plt.legend()
|
||||
@@ -3347,28 +3347,28 @@ 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 \( \beta_0 \). For different values
|
||||
regularization term does not include \( \theta_0 \). For different values
|
||||
of \( \lambda \), 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>
|
||||
<p> <br>
|
||||
$$
|
||||
\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,
|
||||
$$
|
||||
<p> <br>
|
||||
|
||||
<p>but when we take out the intercept, this equation becomes</p>
|
||||
<p> <br>
|
||||
$$
|
||||
\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.
|
||||
$$
|
||||
<p> <br>
|
||||
|
||||
<p>For Lasso regression we have</p>
|
||||
<p> <br>
|
||||
$$
|
||||
\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.
|
||||
$$
|
||||
<p> <br>
|
||||
|
||||
@@ -3425,20 +3425,20 @@ MSERidgePredict = np.zeros(nlambdas)
|
||||
lambdas = np.logspace(-<span style="color: #B452CD">4</span>, <span style="color: #B452CD">2</span>, nlambdas)
|
||||
<span style="color: #8B008B; font-weight: bold">for</span> i <span style="color: #8B008B">in</span> <span style="color: #658b00">range</span>(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
|
||||
<span style="color: #228B22"># Note: we include the intercept column and no scaling</span>
|
||||
RegRidge = linear_model.Ridge(lmb,fit_intercept=<span style="color: #8B008B; font-weight: bold">False</span>)
|
||||
RegRidge.fit(X_train,y_train)
|
||||
<span style="color: #228B22"># and then make the prediction</span>
|
||||
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)
|
||||
<span style="color: #658b00">print</span>(<span style="color: #CD5555">"Beta values for own Ridge implementation"</span>)
|
||||
<span style="color: #658b00">print</span>(OwnRidgeBeta)
|
||||
<span style="color: #658b00">print</span>(<span style="color: #CD5555">"Beta values for Scikit-Learn Ridge implementation"</span>)
|
||||
<span style="color: #658b00">print</span>(<span style="color: #CD5555">"Theta values for own Ridge implementation"</span>)
|
||||
<span style="color: #658b00">print</span>(OwnRidgeTheta)
|
||||
<span style="color: #658b00">print</span>(<span style="color: #CD5555">"Theta values for Scikit-Learn Ridge implementation"</span>)
|
||||
<span style="color: #658b00">print</span>(RegRidge.coef_)
|
||||
<span style="color: #658b00">print</span>(<span style="color: #CD5555">"MSE values for own Ridge implementation"</span>)
|
||||
<span style="color: #658b00">print</span>(MSEOwnRidgePredict[i])
|
||||
@@ -3529,18 +3529,18 @@ MSERidgePredict = np.zeros(nlambdas)
|
||||
lambdas = np.logspace(-<span style="color: #B452CD">4</span>, <span style="color: #B452CD">2</span>, nlambdas)
|
||||
<span style="color: #8B008B; font-weight: bold">for</span> i <span style="color: #8B008B">in</span> <span style="color: #658b00">range</span>(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<span style="color: #707a7c">@OwnRidgeBeta</span> <span style="color: #228B22">#The intercept can be shifted so the model can predict on uncentered data</span>
|
||||
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<span style="color: #707a7c">@OwnRidgeTheta</span> <span style="color: #228B22">#The intercept can be shifted so the model can predict on uncentered data</span>
|
||||
<span style="color: #228B22">#Add intercept to prediction</span>
|
||||
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)
|
||||
<span style="color: #658b00">print</span>(<span style="color: #CD5555">"Beta values for own Ridge implementation"</span>)
|
||||
<span style="color: #658b00">print</span>(OwnRidgeBeta) <span style="color: #228B22">#Intercept is given by mean of target variable</span>
|
||||
<span style="color: #658b00">print</span>(<span style="color: #CD5555">"Beta values for Scikit-Learn Ridge implementation"</span>)
|
||||
<span style="color: #658b00">print</span>(<span style="color: #CD5555">"Theta values for own Ridge implementation"</span>)
|
||||
<span style="color: #658b00">print</span>(OwnRidgeTheta) <span style="color: #228B22">#Intercept is given by mean of target variable</span>
|
||||
<span style="color: #658b00">print</span>(<span style="color: #CD5555">"Theta values for Scikit-Learn Ridge implementation"</span>)
|
||||
<span style="color: #658b00">print</span>(RegRidge.coef_)
|
||||
<span style="color: #658b00">print</span>(<span style="color: #CD5555">'Intercept from own implementation:'</span>)
|
||||
<span style="color: #658b00">print</span>(intercept_)
|
||||
@@ -3578,7 +3578,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
|
||||
\( \beta_0 \) in the fitting. This applies to Lasso regularization as
|
||||
\( \theta_0 \) 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.
|
||||
</p>
|
||||
|
||||
@@ -2775,8 +2775,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 <b>Scikit-Learn</b> (when not shrinking \( \beta_0 \)) for the unknown parameters
|
||||
\( \boldsymbol{\beta} \), are derived under the assumption that both \( \boldsymbol{y} \) and
|
||||
from the library <b>Scikit-Learn</b> (when not shrinking \( \theta_0 \)) for the unknown parameters
|
||||
\( \boldsymbol{\theta} \), are derived under the assumption that both \( \boldsymbol{y} \) and
|
||||
\( \boldsymbol{X} \) are zero centered, that is we subtract the mean values.
|
||||
</p>
|
||||
|
||||
@@ -2860,7 +2860,7 @@ simplicity, we will focus on ordinary regression, as done in the above example.
|
||||
|
||||
<p>The cost/loss function for regression is</p>
|
||||
$$
|
||||
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,.
|
||||
$$
|
||||
|
||||
<p>Recall also that we use the squared value. This expression can lead to an
|
||||
@@ -2868,37 +2868,37 @@ increased penalty for higher differences between predicted and
|
||||
output/target values.
|
||||
</p>
|
||||
|
||||
<p>What we have done is to single out the \( \beta_0 \) term in the
|
||||
<p>What we have done is to single out the \( \theta_0 \) term in the
|
||||
definition of the mean squared error (MSE). The design matrix \( X \)
|
||||
does in this case not contain any intercept column. When we take the
|
||||
derivative with respect to \( \beta_0 \), we want the derivative to obey
|
||||
derivative with respect to \( \theta_0 \), we want the derivative to obey
|
||||
</p>
|
||||
|
||||
$$
|
||||
\frac{\partial C}{\partial \beta_j} = 0,
|
||||
\frac{\partial C}{\partial \theta_j} = 0,
|
||||
$$
|
||||
|
||||
<p>for all \( j \). For \( \beta_0 \) we have</p>
|
||||
<p>for all \( j \). For \( \theta_0 \) we have</p>
|
||||
|
||||
$$
|
||||
\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).
|
||||
$$
|
||||
|
||||
<p>Multiplying away the constant \( 2/n \), we obtain</p>
|
||||
$$
|
||||
\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.
|
||||
$$
|
||||
|
||||
<p>Let us specialize first to the case where we have only two parameters \( \beta_0 \) and \( \beta_1 \).
|
||||
Our result for \( \beta_0 \) simplifies then to
|
||||
<p>Let us specialize first to the case where we have only two parameters \( \theta_0 \) and \( \theta_1 \).
|
||||
Our result for \( \theta_0 \) simplifies then to
|
||||
</p>
|
||||
$$
|
||||
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.
|
||||
$$
|
||||
|
||||
<p>We obtain then</p>
|
||||
$$
|
||||
\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}.
|
||||
$$
|
||||
|
||||
<p>If we define</p>
|
||||
@@ -2913,17 +2913,17 @@ $$
|
||||
|
||||
<p>we have</p>
|
||||
$$
|
||||
\beta_0 = \mu_y - \beta_1\mu_{\boldsymbol{x}_1}.
|
||||
\theta_0 = \mu_y - \theta_1\mu_{\boldsymbol{x}_1}.
|
||||
$$
|
||||
|
||||
<p>In the general case with more parameters than \( \beta_0 \) and \( \beta_1 \), we have</p>
|
||||
<p>In the general case with more parameters than \( \theta_0 \) and \( \theta_1 \), we have</p>
|
||||
$$
|
||||
\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.
|
||||
$$
|
||||
|
||||
<p>We can rewrite the latter equation as</p>
|
||||
$$
|
||||
\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,
|
||||
$$
|
||||
|
||||
<p>where we have defined</p>
|
||||
@@ -2935,22 +2935,22 @@ $$
|
||||
|
||||
<p>Replacing \( y_i \) with \( y_i - y_i - \overline{\boldsymbol{y}} \) and centering also our design matrix results in a cost function (in vector-matrix disguise)</p>
|
||||
$$
|
||||
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}).
|
||||
$$
|
||||
|
||||
<p>If we minimize with respect to \( \boldsymbol{\beta} \) we have then</p>
|
||||
<p>If we minimize with respect to \( \boldsymbol{\theta} \) we have then</p>
|
||||
|
||||
$$
|
||||
\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}},
|
||||
$$
|
||||
|
||||
<p>where \( \boldsymbol{\tilde{y}} = \boldsymbol{y} - \overline{\boldsymbol{y}} \)
|
||||
and \( \tilde{X}_{ij} = X_{ij} - \frac{1}{n}\sum_{k=0}^{n-1}X_{kj} \).
|
||||
</p>
|
||||
|
||||
<p>For Ridge regression we need to add \( \lambda \boldsymbol{\beta}^T\boldsymbol{\beta} \) to the cost function and get then</p>
|
||||
<p>For Ridge regression we need to add \( \lambda \boldsymbol{\theta}^T\boldsymbol{\theta} \) to the cost function and get then</p>
|
||||
$$
|
||||
\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}}.
|
||||
$$
|
||||
|
||||
<p>What does this mean? And why do we insist on all this? Let us look at some examples.</p>
|
||||
@@ -2979,15 +2979,15 @@ np.random.seed(<span style="color: #B452CD">2021</span>)
|
||||
<span style="color: #8B008B; font-weight: bold">return</span> np.sum((y_data-y_model)**<span style="color: #B452CD">2</span>)/n
|
||||
|
||||
|
||||
<span style="color: #8B008B; font-weight: bold">def</span> <span style="color: #008b45">fit_beta</span>(X, y):
|
||||
<span style="color: #8B008B; font-weight: bold">def</span> <span style="color: #008b45">fit_theta</span>(X, y):
|
||||
<span style="color: #8B008B; font-weight: bold">return</span> np.linalg.pinv(X.T @ X) @ X.T @ y
|
||||
|
||||
|
||||
true_beta = [<span style="color: #B452CD">2</span>, <span style="color: #B452CD">0.5</span>, <span style="color: #B452CD">3.7</span>]
|
||||
true_theta = [<span style="color: #B452CD">2</span>, <span style="color: #B452CD">0.5</span>, <span style="color: #B452CD">3.7</span>]
|
||||
|
||||
x = np.linspace(<span style="color: #B452CD">0</span>, <span style="color: #B452CD">1</span>, <span style="color: #B452CD">11</span>)
|
||||
y = np.sum(
|
||||
np.asarray([x ** p * b <span style="color: #8B008B; font-weight: bold">for</span> p, b <span style="color: #8B008B">in</span> <span style="color: #658b00">enumerate</span>(true_beta)]), axis=<span style="color: #B452CD">0</span>
|
||||
np.asarray([x ** p * b <span style="color: #8B008B; font-weight: bold">for</span> p, b <span style="color: #8B008B">in</span> <span style="color: #658b00">enumerate</span>(true_theta)]), axis=<span style="color: #B452CD">0</span>
|
||||
) + <span style="color: #B452CD">0.1</span> * np.random.normal(size=<span style="color: #658b00">len</span>(x))
|
||||
|
||||
degree = <span style="color: #B452CD">3</span>
|
||||
@@ -2997,15 +2997,15 @@ X = np.zeros((<span style="color: #658b00">len</span>(x), degree))
|
||||
<span style="color: #8B008B; font-weight: bold">for</span> p <span style="color: #8B008B">in</span> <span style="color: #658b00">range</span>(degree):
|
||||
X[:, p] = x ** p
|
||||
|
||||
beta = fit_beta(X, y)
|
||||
theta = fit_theta(X, y)
|
||||
|
||||
<span style="color: #228B22"># Intercept is included in the design matrix</span>
|
||||
skl = LinearRegression(fit_intercept=<span style="color: #8B008B; font-weight: bold">False</span>).fit(X, y)
|
||||
|
||||
<span style="color: #658b00">print</span>(<span style="color: #CD5555">f"True beta: {</span>true_beta<span style="color: #CD5555">}"</span>)
|
||||
<span style="color: #658b00">print</span>(<span style="color: #CD5555">f"Fitted beta: {</span>beta<span style="color: #CD5555">}"</span>)
|
||||
<span style="color: #658b00">print</span>(<span style="color: #CD5555">f"Sklearn fitted beta: {</span>skl.coef_<span style="color: #CD5555">}"</span>)
|
||||
ypredictOwn = X @ beta
|
||||
<span style="color: #658b00">print</span>(<span style="color: #CD5555">f"True theta: {</span>true_theta<span style="color: #CD5555">}"</span>)
|
||||
<span style="color: #658b00">print</span>(<span style="color: #CD5555">f"Fitted theta: {</span>theta<span style="color: #CD5555">}"</span>)
|
||||
<span style="color: #658b00">print</span>(<span style="color: #CD5555">f"Sklearn fitted theta: {</span>skl.coef_<span style="color: #CD5555">}"</span>)
|
||||
ypredictOwn = X @ theta
|
||||
ypredictSKL = skl.predict(X)
|
||||
<span style="color: #658b00">print</span>(<span style="color: #CD5555">f"MSE with intercept column"</span>)
|
||||
<span style="color: #658b00">print</span>(MSE(y,ypredictOwn))
|
||||
@@ -3015,7 +3015,7 @@ ypredictSKL = skl.predict(X)
|
||||
|
||||
plt.figure()
|
||||
plt.scatter(x, y, label=<span style="color: #CD5555">"Data"</span>)
|
||||
plt.plot(x, X @ beta, label=<span style="color: #CD5555">"Fit"</span>)
|
||||
plt.plot(x, X @ theta, label=<span style="color: #CD5555">"Fit"</span>)
|
||||
plt.plot(x, skl.predict(X), label=<span style="color: #CD5555">"Sklearn (fit_intercept=False)"</span>)
|
||||
|
||||
|
||||
@@ -3032,21 +3032,21 @@ skl = LinearRegression(fit_intercept=<span style="color: #8B008B; font-weight: b
|
||||
y_offset = np.average(y, axis=<span style="color: #B452CD">0</span>)
|
||||
X_offset = np.average(X, axis=<span style="color: #B452CD">0</span>)
|
||||
|
||||
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)
|
||||
|
||||
<span style="color: #658b00">print</span>(<span style="color: #CD5555">f"Manual intercept: {</span>intercept<span style="color: #CD5555">}"</span>)
|
||||
<span style="color: #658b00">print</span>(<span style="color: #CD5555">f"Fitted beta (without intercept): {</span>beta<span style="color: #CD5555">}"</span>)
|
||||
<span style="color: #658b00">print</span>(<span style="color: #CD5555">f"Fitted theta (without intercept): {</span>theta<span style="color: #CD5555">}"</span>)
|
||||
<span style="color: #658b00">print</span>(<span style="color: #CD5555">f"Sklearn intercept: {</span>skl.intercept_<span style="color: #CD5555">}"</span>)
|
||||
<span style="color: #658b00">print</span>(<span style="color: #CD5555">f"Sklearn fitted beta (without intercept): {</span>skl.coef_<span style="color: #CD5555">}"</span>)
|
||||
ypredictOwn = X @ beta
|
||||
<span style="color: #658b00">print</span>(<span style="color: #CD5555">f"Sklearn fitted theta (without intercept): {</span>skl.coef_<span style="color: #CD5555">}"</span>)
|
||||
ypredictOwn = X @ theta
|
||||
ypredictSKL = skl.predict(X)
|
||||
<span style="color: #658b00">print</span>(<span style="color: #CD5555">f"MSE with Manual intercept"</span>)
|
||||
<span style="color: #658b00">print</span>(MSE(y,ypredictOwn+intercept))
|
||||
<span style="color: #658b00">print</span>(<span style="color: #CD5555">f"MSE with Sklearn intercept"</span>)
|
||||
<span style="color: #658b00">print</span>(MSE(y,ypredictSKL))
|
||||
|
||||
plt.plot(x, X @ beta + intercept, <span style="color: #CD5555">"--"</span>, label=<span style="color: #CD5555">"Fit (manual intercept)"</span>)
|
||||
plt.plot(x, X @ theta + intercept, <span style="color: #CD5555">"--"</span>, label=<span style="color: #CD5555">"Fit (manual intercept)"</span>)
|
||||
plt.plot(x, skl.predict(X), <span style="color: #CD5555">"--"</span>, label=<span style="color: #CD5555">"Sklearn (fit_intercept=True)"</span>)
|
||||
plt.grid()
|
||||
plt.legend()
|
||||
@@ -3076,23 +3076,23 @@ 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 \( \beta_0 \). For different values
|
||||
regularization term does not include \( \theta_0 \). For different values
|
||||
of \( \lambda \), 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>
|
||||
$$
|
||||
\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,
|
||||
$$
|
||||
|
||||
<p>but when we take out the intercept, this equation becomes</p>
|
||||
$$
|
||||
\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.
|
||||
$$
|
||||
|
||||
<p>For Lasso regression we have</p>
|
||||
$$
|
||||
\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.
|
||||
$$
|
||||
|
||||
<p>It means that, when scaling the design matrix and the outputs/targets,
|
||||
@@ -3148,20 +3148,20 @@ MSERidgePredict = np.zeros(nlambdas)
|
||||
lambdas = np.logspace(-<span style="color: #B452CD">4</span>, <span style="color: #B452CD">2</span>, nlambdas)
|
||||
<span style="color: #8B008B; font-weight: bold">for</span> i <span style="color: #8B008B">in</span> <span style="color: #658b00">range</span>(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
|
||||
<span style="color: #228B22"># Note: we include the intercept column and no scaling</span>
|
||||
RegRidge = linear_model.Ridge(lmb,fit_intercept=<span style="color: #8B008B; font-weight: bold">False</span>)
|
||||
RegRidge.fit(X_train,y_train)
|
||||
<span style="color: #228B22"># and then make the prediction</span>
|
||||
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)
|
||||
<span style="color: #658b00">print</span>(<span style="color: #CD5555">"Beta values for own Ridge implementation"</span>)
|
||||
<span style="color: #658b00">print</span>(OwnRidgeBeta)
|
||||
<span style="color: #658b00">print</span>(<span style="color: #CD5555">"Beta values for Scikit-Learn Ridge implementation"</span>)
|
||||
<span style="color: #658b00">print</span>(<span style="color: #CD5555">"Theta values for own Ridge implementation"</span>)
|
||||
<span style="color: #658b00">print</span>(OwnRidgeTheta)
|
||||
<span style="color: #658b00">print</span>(<span style="color: #CD5555">"Theta values for Scikit-Learn Ridge implementation"</span>)
|
||||
<span style="color: #658b00">print</span>(RegRidge.coef_)
|
||||
<span style="color: #658b00">print</span>(<span style="color: #CD5555">"MSE values for own Ridge implementation"</span>)
|
||||
<span style="color: #658b00">print</span>(MSEOwnRidgePredict[i])
|
||||
@@ -3252,18 +3252,18 @@ MSERidgePredict = np.zeros(nlambdas)
|
||||
lambdas = np.logspace(-<span style="color: #B452CD">4</span>, <span style="color: #B452CD">2</span>, nlambdas)
|
||||
<span style="color: #8B008B; font-weight: bold">for</span> i <span style="color: #8B008B">in</span> <span style="color: #658b00">range</span>(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<span style="color: #707a7c">@OwnRidgeBeta</span> <span style="color: #228B22">#The intercept can be shifted so the model can predict on uncentered data</span>
|
||||
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<span style="color: #707a7c">@OwnRidgeTheta</span> <span style="color: #228B22">#The intercept can be shifted so the model can predict on uncentered data</span>
|
||||
<span style="color: #228B22">#Add intercept to prediction</span>
|
||||
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)
|
||||
<span style="color: #658b00">print</span>(<span style="color: #CD5555">"Beta values for own Ridge implementation"</span>)
|
||||
<span style="color: #658b00">print</span>(OwnRidgeBeta) <span style="color: #228B22">#Intercept is given by mean of target variable</span>
|
||||
<span style="color: #658b00">print</span>(<span style="color: #CD5555">"Beta values for Scikit-Learn Ridge implementation"</span>)
|
||||
<span style="color: #658b00">print</span>(<span style="color: #CD5555">"Theta values for own Ridge implementation"</span>)
|
||||
<span style="color: #658b00">print</span>(OwnRidgeTheta) <span style="color: #228B22">#Intercept is given by mean of target variable</span>
|
||||
<span style="color: #658b00">print</span>(<span style="color: #CD5555">"Theta values for Scikit-Learn Ridge implementation"</span>)
|
||||
<span style="color: #658b00">print</span>(RegRidge.coef_)
|
||||
<span style="color: #658b00">print</span>(<span style="color: #CD5555">'Intercept from own implementation:'</span>)
|
||||
<span style="color: #658b00">print</span>(intercept_)
|
||||
@@ -3301,7 +3301,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
|
||||
\( \beta_0 \) in the fitting. This applies to Lasso regularization as
|
||||
\( \theta_0 \) 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.
|
||||
</p>
|
||||
|
||||
@@ -2852,8 +2852,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 <b>Scikit-Learn</b> (when not shrinking \( \beta_0 \)) for the unknown parameters
|
||||
\( \boldsymbol{\beta} \), are derived under the assumption that both \( \boldsymbol{y} \) and
|
||||
from the library <b>Scikit-Learn</b> (when not shrinking \( \theta_0 \)) for the unknown parameters
|
||||
\( \boldsymbol{\theta} \), are derived under the assumption that both \( \boldsymbol{y} \) and
|
||||
\( \boldsymbol{X} \) are zero centered, that is we subtract the mean values.
|
||||
</p>
|
||||
|
||||
@@ -2937,7 +2937,7 @@ simplicity, we will focus on ordinary regression, as done in the above example.
|
||||
|
||||
<p>The cost/loss function for regression is</p>
|
||||
$$
|
||||
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,.
|
||||
$$
|
||||
|
||||
<p>Recall also that we use the squared value. This expression can lead to an
|
||||
@@ -2945,37 +2945,37 @@ increased penalty for higher differences between predicted and
|
||||
output/target values.
|
||||
</p>
|
||||
|
||||
<p>What we have done is to single out the \( \beta_0 \) term in the
|
||||
<p>What we have done is to single out the \( \theta_0 \) term in the
|
||||
definition of the mean squared error (MSE). The design matrix \( X \)
|
||||
does in this case not contain any intercept column. When we take the
|
||||
derivative with respect to \( \beta_0 \), we want the derivative to obey
|
||||
derivative with respect to \( \theta_0 \), we want the derivative to obey
|
||||
</p>
|
||||
|
||||
$$
|
||||
\frac{\partial C}{\partial \beta_j} = 0,
|
||||
\frac{\partial C}{\partial \theta_j} = 0,
|
||||
$$
|
||||
|
||||
<p>for all \( j \). For \( \beta_0 \) we have</p>
|
||||
<p>for all \( j \). For \( \theta_0 \) we have</p>
|
||||
|
||||
$$
|
||||
\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).
|
||||
$$
|
||||
|
||||
<p>Multiplying away the constant \( 2/n \), we obtain</p>
|
||||
$$
|
||||
\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.
|
||||
$$
|
||||
|
||||
<p>Let us specialize first to the case where we have only two parameters \( \beta_0 \) and \( \beta_1 \).
|
||||
Our result for \( \beta_0 \) simplifies then to
|
||||
<p>Let us specialize first to the case where we have only two parameters \( \theta_0 \) and \( \theta_1 \).
|
||||
Our result for \( \theta_0 \) simplifies then to
|
||||
</p>
|
||||
$$
|
||||
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.
|
||||
$$
|
||||
|
||||
<p>We obtain then</p>
|
||||
$$
|
||||
\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}.
|
||||
$$
|
||||
|
||||
<p>If we define</p>
|
||||
@@ -2990,17 +2990,17 @@ $$
|
||||
|
||||
<p>we have</p>
|
||||
$$
|
||||
\beta_0 = \mu_y - \beta_1\mu_{\boldsymbol{x}_1}.
|
||||
\theta_0 = \mu_y - \theta_1\mu_{\boldsymbol{x}_1}.
|
||||
$$
|
||||
|
||||
<p>In the general case with more parameters than \( \beta_0 \) and \( \beta_1 \), we have</p>
|
||||
<p>In the general case with more parameters than \( \theta_0 \) and \( \theta_1 \), we have</p>
|
||||
$$
|
||||
\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.
|
||||
$$
|
||||
|
||||
<p>We can rewrite the latter equation as</p>
|
||||
$$
|
||||
\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,
|
||||
$$
|
||||
|
||||
<p>where we have defined</p>
|
||||
@@ -3012,22 +3012,22 @@ $$
|
||||
|
||||
<p>Replacing \( y_i \) with \( y_i - y_i - \overline{\boldsymbol{y}} \) and centering also our design matrix results in a cost function (in vector-matrix disguise)</p>
|
||||
$$
|
||||
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}).
|
||||
$$
|
||||
|
||||
<p>If we minimize with respect to \( \boldsymbol{\beta} \) we have then</p>
|
||||
<p>If we minimize with respect to \( \boldsymbol{\theta} \) we have then</p>
|
||||
|
||||
$$
|
||||
\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}},
|
||||
$$
|
||||
|
||||
<p>where \( \boldsymbol{\tilde{y}} = \boldsymbol{y} - \overline{\boldsymbol{y}} \)
|
||||
and \( \tilde{X}_{ij} = X_{ij} - \frac{1}{n}\sum_{k=0}^{n-1}X_{kj} \).
|
||||
</p>
|
||||
|
||||
<p>For Ridge regression we need to add \( \lambda \boldsymbol{\beta}^T\boldsymbol{\beta} \) to the cost function and get then</p>
|
||||
<p>For Ridge regression we need to add \( \lambda \boldsymbol{\theta}^T\boldsymbol{\theta} \) to the cost function and get then</p>
|
||||
$$
|
||||
\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}}.
|
||||
$$
|
||||
|
||||
<p>What does this mean? And why do we insist on all this? Let us look at some examples.</p>
|
||||
@@ -3056,15 +3056,15 @@ np<span style="color: #666666">.</span>random<span style="color: #666666">.</spa
|
||||
<span style="color: #008000; font-weight: bold">return</span> np<span style="color: #666666">.</span>sum((y_data<span style="color: #666666">-</span>y_model)<span style="color: #666666">**2</span>)<span style="color: #666666">/</span>n
|
||||
|
||||
|
||||
<span style="color: #008000; font-weight: bold">def</span> <span style="color: #0000FF">fit_beta</span>(X, y):
|
||||
<span style="color: #008000; font-weight: bold">def</span> <span style="color: #0000FF">fit_theta</span>(X, y):
|
||||
<span style="color: #008000; font-weight: bold">return</span> np<span style="color: #666666">.</span>linalg<span style="color: #666666">.</span>pinv(X<span style="color: #666666">.</span>T <span style="color: #666666">@</span> X) <span style="color: #666666">@</span> X<span style="color: #666666">.</span>T <span style="color: #666666">@</span> y
|
||||
|
||||
|
||||
true_beta <span style="color: #666666">=</span> [<span style="color: #666666">2</span>, <span style="color: #666666">0.5</span>, <span style="color: #666666">3.7</span>]
|
||||
true_theta <span style="color: #666666">=</span> [<span style="color: #666666">2</span>, <span style="color: #666666">0.5</span>, <span style="color: #666666">3.7</span>]
|
||||
|
||||
x <span style="color: #666666">=</span> np<span style="color: #666666">.</span>linspace(<span style="color: #666666">0</span>, <span style="color: #666666">1</span>, <span style="color: #666666">11</span>)
|
||||
y <span style="color: #666666">=</span> np<span style="color: #666666">.</span>sum(
|
||||
np<span style="color: #666666">.</span>asarray([x <span style="color: #666666">**</span> p <span style="color: #666666">*</span> b <span style="color: #008000; font-weight: bold">for</span> p, b <span style="color: #AA22FF; font-weight: bold">in</span> <span style="color: #008000">enumerate</span>(true_beta)]), axis<span style="color: #666666">=0</span>
|
||||
np<span style="color: #666666">.</span>asarray([x <span style="color: #666666">**</span> p <span style="color: #666666">*</span> b <span style="color: #008000; font-weight: bold">for</span> p, b <span style="color: #AA22FF; font-weight: bold">in</span> <span style="color: #008000">enumerate</span>(true_theta)]), axis<span style="color: #666666">=0</span>
|
||||
) <span style="color: #666666">+</span> <span style="color: #666666">0.1</span> <span style="color: #666666">*</span> np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>normal(size<span style="color: #666666">=</span><span style="color: #008000">len</span>(x))
|
||||
|
||||
degree <span style="color: #666666">=</span> <span style="color: #666666">3</span>
|
||||
@@ -3074,15 +3074,15 @@ X <span style="color: #666666">=</span> np<span style="color: #666666">.</span>z
|
||||
<span style="color: #008000; font-weight: bold">for</span> p <span style="color: #AA22FF; font-weight: bold">in</span> <span style="color: #008000">range</span>(degree):
|
||||
X[:, p] <span style="color: #666666">=</span> x <span style="color: #666666">**</span> p
|
||||
|
||||
beta <span style="color: #666666">=</span> fit_beta(X, y)
|
||||
theta <span style="color: #666666">=</span> fit_theta(X, y)
|
||||
|
||||
<span style="color: #408080; font-style: italic"># Intercept is included in the design matrix</span>
|
||||
skl <span style="color: #666666">=</span> LinearRegression(fit_intercept<span style="color: #666666">=</span><span style="color: #008000; font-weight: bold">False</span>)<span style="color: #666666">.</span>fit(X, y)
|
||||
|
||||
<span style="color: #008000">print</span>(<span style="color: #BA2121">f"True beta: </span><span style="color: #BB6688; font-weight: bold">{</span>true_beta<span style="color: #BB6688; font-weight: bold">}</span><span style="color: #BA2121">"</span>)
|
||||
<span style="color: #008000">print</span>(<span style="color: #BA2121">f"Fitted beta: </span><span style="color: #BB6688; font-weight: bold">{</span>beta<span style="color: #BB6688; font-weight: bold">}</span><span style="color: #BA2121">"</span>)
|
||||
<span style="color: #008000">print</span>(<span style="color: #BA2121">f"Sklearn fitted beta: </span><span style="color: #BB6688; font-weight: bold">{</span>skl<span style="color: #666666">.</span>coef_<span style="color: #BB6688; font-weight: bold">}</span><span style="color: #BA2121">"</span>)
|
||||
ypredictOwn <span style="color: #666666">=</span> X <span style="color: #666666">@</span> beta
|
||||
<span style="color: #008000">print</span>(<span style="color: #BA2121">f"True theta: </span><span style="color: #BB6688; font-weight: bold">{</span>true_theta<span style="color: #BB6688; font-weight: bold">}</span><span style="color: #BA2121">"</span>)
|
||||
<span style="color: #008000">print</span>(<span style="color: #BA2121">f"Fitted theta: </span><span style="color: #BB6688; font-weight: bold">{</span>theta<span style="color: #BB6688; font-weight: bold">}</span><span style="color: #BA2121">"</span>)
|
||||
<span style="color: #008000">print</span>(<span style="color: #BA2121">f"Sklearn fitted theta: </span><span style="color: #BB6688; font-weight: bold">{</span>skl<span style="color: #666666">.</span>coef_<span style="color: #BB6688; font-weight: bold">}</span><span style="color: #BA2121">"</span>)
|
||||
ypredictOwn <span style="color: #666666">=</span> X <span style="color: #666666">@</span> theta
|
||||
ypredictSKL <span style="color: #666666">=</span> skl<span style="color: #666666">.</span>predict(X)
|
||||
<span style="color: #008000">print</span>(<span style="color: #BA2121">f"MSE with intercept column"</span>)
|
||||
<span style="color: #008000">print</span>(MSE(y,ypredictOwn))
|
||||
@@ -3092,7 +3092,7 @@ ypredictSKL <span style="color: #666666">=</span> skl<span style="color: #666666
|
||||
|
||||
plt<span style="color: #666666">.</span>figure()
|
||||
plt<span style="color: #666666">.</span>scatter(x, y, label<span style="color: #666666">=</span><span style="color: #BA2121">"Data"</span>)
|
||||
plt<span style="color: #666666">.</span>plot(x, X <span style="color: #666666">@</span> beta, label<span style="color: #666666">=</span><span style="color: #BA2121">"Fit"</span>)
|
||||
plt<span style="color: #666666">.</span>plot(x, X <span style="color: #666666">@</span> theta, label<span style="color: #666666">=</span><span style="color: #BA2121">"Fit"</span>)
|
||||
plt<span style="color: #666666">.</span>plot(x, skl<span style="color: #666666">.</span>predict(X), label<span style="color: #666666">=</span><span style="color: #BA2121">"Sklearn (fit_intercept=False)"</span>)
|
||||
|
||||
|
||||
@@ -3109,21 +3109,21 @@ skl <span style="color: #666666">=</span> LinearRegression(fit_intercept<span st
|
||||
y_offset <span style="color: #666666">=</span> np<span style="color: #666666">.</span>average(y, axis<span style="color: #666666">=0</span>)
|
||||
X_offset <span style="color: #666666">=</span> np<span style="color: #666666">.</span>average(X, axis<span style="color: #666666">=0</span>)
|
||||
|
||||
beta <span style="color: #666666">=</span> fit_beta(X <span style="color: #666666">-</span> X_offset, y <span style="color: #666666">-</span> y_offset)
|
||||
intercept <span style="color: #666666">=</span> np<span style="color: #666666">.</span>mean(y_offset <span style="color: #666666">-</span> X_offset <span style="color: #666666">@</span> beta)
|
||||
theta <span style="color: #666666">=</span> fit_theta(X <span style="color: #666666">-</span> X_offset, y <span style="color: #666666">-</span> y_offset)
|
||||
intercept <span style="color: #666666">=</span> np<span style="color: #666666">.</span>mean(y_offset <span style="color: #666666">-</span> X_offset <span style="color: #666666">@</span> theta)
|
||||
|
||||
<span style="color: #008000">print</span>(<span style="color: #BA2121">f"Manual intercept: </span><span style="color: #BB6688; font-weight: bold">{</span>intercept<span style="color: #BB6688; font-weight: bold">}</span><span style="color: #BA2121">"</span>)
|
||||
<span style="color: #008000">print</span>(<span style="color: #BA2121">f"Fitted beta (without intercept): </span><span style="color: #BB6688; font-weight: bold">{</span>beta<span style="color: #BB6688; font-weight: bold">}</span><span style="color: #BA2121">"</span>)
|
||||
<span style="color: #008000">print</span>(<span style="color: #BA2121">f"Fitted theta (without intercept): </span><span style="color: #BB6688; font-weight: bold">{</span>theta<span style="color: #BB6688; font-weight: bold">}</span><span style="color: #BA2121">"</span>)
|
||||
<span style="color: #008000">print</span>(<span style="color: #BA2121">f"Sklearn intercept: </span><span style="color: #BB6688; font-weight: bold">{</span>skl<span style="color: #666666">.</span>intercept_<span style="color: #BB6688; font-weight: bold">}</span><span style="color: #BA2121">"</span>)
|
||||
<span style="color: #008000">print</span>(<span style="color: #BA2121">f"Sklearn fitted beta (without intercept): </span><span style="color: #BB6688; font-weight: bold">{</span>skl<span style="color: #666666">.</span>coef_<span style="color: #BB6688; font-weight: bold">}</span><span style="color: #BA2121">"</span>)
|
||||
ypredictOwn <span style="color: #666666">=</span> X <span style="color: #666666">@</span> beta
|
||||
<span style="color: #008000">print</span>(<span style="color: #BA2121">f"Sklearn fitted theta (without intercept): </span><span style="color: #BB6688; font-weight: bold">{</span>skl<span style="color: #666666">.</span>coef_<span style="color: #BB6688; font-weight: bold">}</span><span style="color: #BA2121">"</span>)
|
||||
ypredictOwn <span style="color: #666666">=</span> X <span style="color: #666666">@</span> theta
|
||||
ypredictSKL <span style="color: #666666">=</span> skl<span style="color: #666666">.</span>predict(X)
|
||||
<span style="color: #008000">print</span>(<span style="color: #BA2121">f"MSE with Manual intercept"</span>)
|
||||
<span style="color: #008000">print</span>(MSE(y,ypredictOwn<span style="color: #666666">+</span>intercept))
|
||||
<span style="color: #008000">print</span>(<span style="color: #BA2121">f"MSE with Sklearn intercept"</span>)
|
||||
<span style="color: #008000">print</span>(MSE(y,ypredictSKL))
|
||||
|
||||
plt<span style="color: #666666">.</span>plot(x, X <span style="color: #666666">@</span> beta <span style="color: #666666">+</span> intercept, <span style="color: #BA2121">"--"</span>, label<span style="color: #666666">=</span><span style="color: #BA2121">"Fit (manual intercept)"</span>)
|
||||
plt<span style="color: #666666">.</span>plot(x, X <span style="color: #666666">@</span> theta <span style="color: #666666">+</span> intercept, <span style="color: #BA2121">"--"</span>, label<span style="color: #666666">=</span><span style="color: #BA2121">"Fit (manual intercept)"</span>)
|
||||
plt<span style="color: #666666">.</span>plot(x, skl<span style="color: #666666">.</span>predict(X), <span style="color: #BA2121">"--"</span>, label<span style="color: #666666">=</span><span style="color: #BA2121">"Sklearn (fit_intercept=True)"</span>)
|
||||
plt<span style="color: #666666">.</span>grid()
|
||||
plt<span style="color: #666666">.</span>legend()
|
||||
@@ -3153,23 +3153,23 @@ 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 \( \beta_0 \). For different values
|
||||
regularization term does not include \( \theta_0 \). For different values
|
||||
of \( \lambda \), 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>
|
||||
$$
|
||||
\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,
|
||||
$$
|
||||
|
||||
<p>but when we take out the intercept, this equation becomes</p>
|
||||
$$
|
||||
\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.
|
||||
$$
|
||||
|
||||
<p>For Lasso regression we have</p>
|
||||
$$
|
||||
\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.
|
||||
$$
|
||||
|
||||
<p>It means that, when scaling the design matrix and the outputs/targets,
|
||||
@@ -3225,20 +3225,20 @@ MSERidgePredict <span style="color: #666666">=</span> np<span style="color: #666
|
||||
lambdas <span style="color: #666666">=</span> np<span style="color: #666666">.</span>logspace(<span style="color: #666666">-4</span>, <span style="color: #666666">2</span>, nlambdas)
|
||||
<span style="color: #008000; font-weight: bold">for</span> i <span style="color: #AA22FF; font-weight: bold">in</span> <span style="color: #008000">range</span>(nlambdas):
|
||||
lmb <span style="color: #666666">=</span> lambdas[i]
|
||||
OwnRidgeBeta <span style="color: #666666">=</span> np<span style="color: #666666">.</span>linalg<span style="color: #666666">.</span>pinv(X_train<span style="color: #666666">.</span>T <span style="color: #666666">@</span> X_train<span style="color: #666666">+</span>lmb<span style="color: #666666">*</span>I) <span style="color: #666666">@</span> X_train<span style="color: #666666">.</span>T <span style="color: #666666">@</span> y_train
|
||||
OwnRidgeTheta <span style="color: #666666">=</span> np<span style="color: #666666">.</span>linalg<span style="color: #666666">.</span>pinv(X_train<span style="color: #666666">.</span>T <span style="color: #666666">@</span> X_train<span style="color: #666666">+</span>lmb<span style="color: #666666">*</span>I) <span style="color: #666666">@</span> X_train<span style="color: #666666">.</span>T <span style="color: #666666">@</span> y_train
|
||||
<span style="color: #408080; font-style: italic"># Note: we include the intercept column and no scaling</span>
|
||||
RegRidge <span style="color: #666666">=</span> linear_model<span style="color: #666666">.</span>Ridge(lmb,fit_intercept<span style="color: #666666">=</span><span style="color: #008000; font-weight: bold">False</span>)
|
||||
RegRidge<span style="color: #666666">.</span>fit(X_train,y_train)
|
||||
<span style="color: #408080; font-style: italic"># and then make the prediction</span>
|
||||
ytildeOwnRidge <span style="color: #666666">=</span> X_train <span style="color: #666666">@</span> OwnRidgeBeta
|
||||
ypredictOwnRidge <span style="color: #666666">=</span> X_test <span style="color: #666666">@</span> OwnRidgeBeta
|
||||
ytildeOwnRidge <span style="color: #666666">=</span> X_train <span style="color: #666666">@</span> OwnRidgeTheta
|
||||
ypredictOwnRidge <span style="color: #666666">=</span> X_test <span style="color: #666666">@</span> OwnRidgeTheta
|
||||
ytildeRidge <span style="color: #666666">=</span> RegRidge<span style="color: #666666">.</span>predict(X_train)
|
||||
ypredictRidge <span style="color: #666666">=</span> RegRidge<span style="color: #666666">.</span>predict(X_test)
|
||||
MSEOwnRidgePredict[i] <span style="color: #666666">=</span> MSE(y_test,ypredictOwnRidge)
|
||||
MSERidgePredict[i] <span style="color: #666666">=</span> MSE(y_test,ypredictRidge)
|
||||
<span style="color: #008000">print</span>(<span style="color: #BA2121">"Beta values for own Ridge implementation"</span>)
|
||||
<span style="color: #008000">print</span>(OwnRidgeBeta)
|
||||
<span style="color: #008000">print</span>(<span style="color: #BA2121">"Beta values for Scikit-Learn Ridge implementation"</span>)
|
||||
<span style="color: #008000">print</span>(<span style="color: #BA2121">"Theta values for own Ridge implementation"</span>)
|
||||
<span style="color: #008000">print</span>(OwnRidgeTheta)
|
||||
<span style="color: #008000">print</span>(<span style="color: #BA2121">"Theta values for Scikit-Learn Ridge implementation"</span>)
|
||||
<span style="color: #008000">print</span>(RegRidge<span style="color: #666666">.</span>coef_)
|
||||
<span style="color: #008000">print</span>(<span style="color: #BA2121">"MSE values for own Ridge implementation"</span>)
|
||||
<span style="color: #008000">print</span>(MSEOwnRidgePredict[i])
|
||||
@@ -3329,18 +3329,18 @@ MSERidgePredict <span style="color: #666666">=</span> np<span style="color: #666
|
||||
lambdas <span style="color: #666666">=</span> np<span style="color: #666666">.</span>logspace(<span style="color: #666666">-4</span>, <span style="color: #666666">2</span>, nlambdas)
|
||||
<span style="color: #008000; font-weight: bold">for</span> i <span style="color: #AA22FF; font-weight: bold">in</span> <span style="color: #008000">range</span>(nlambdas):
|
||||
lmb <span style="color: #666666">=</span> lambdas[i]
|
||||
OwnRidgeBeta <span style="color: #666666">=</span> np<span style="color: #666666">.</span>linalg<span style="color: #666666">.</span>pinv(X_train_scaled<span style="color: #666666">.</span>T <span style="color: #666666">@</span> X_train_scaled<span style="color: #666666">+</span>lmb<span style="color: #666666">*</span>I) <span style="color: #666666">@</span> X_train_scaled<span style="color: #666666">.</span>T <span style="color: #666666">@</span> (y_train_scaled)
|
||||
intercept_ <span style="color: #666666">=</span> y_scaler <span style="color: #666666">-</span> X_train_mean<span style="color: #AA22FF">@OwnRidgeBeta</span> <span style="color: #408080; font-style: italic">#The intercept can be shifted so the model can predict on uncentered data</span>
|
||||
OwnRidgeTheta <span style="color: #666666">=</span> np<span style="color: #666666">.</span>linalg<span style="color: #666666">.</span>pinv(X_train_scaled<span style="color: #666666">.</span>T <span style="color: #666666">@</span> X_train_scaled<span style="color: #666666">+</span>lmb<span style="color: #666666">*</span>I) <span style="color: #666666">@</span> X_train_scaled<span style="color: #666666">.</span>T <span style="color: #666666">@</span> (y_train_scaled)
|
||||
intercept_ <span style="color: #666666">=</span> y_scaler <span style="color: #666666">-</span> X_train_mean<span style="color: #AA22FF">@OwnRidgeTheta</span> <span style="color: #408080; font-style: italic">#The intercept can be shifted so the model can predict on uncentered data</span>
|
||||
<span style="color: #408080; font-style: italic">#Add intercept to prediction</span>
|
||||
ypredictOwnRidge <span style="color: #666666">=</span> X_test_scaled <span style="color: #666666">@</span> OwnRidgeBeta <span style="color: #666666">+</span> y_scaler
|
||||
ypredictOwnRidge <span style="color: #666666">=</span> X_test_scaled <span style="color: #666666">@</span> OwnRidgeTheta <span style="color: #666666">+</span> y_scaler
|
||||
RegRidge <span style="color: #666666">=</span> linear_model<span style="color: #666666">.</span>Ridge(lmb)
|
||||
RegRidge<span style="color: #666666">.</span>fit(X_train,y_train)
|
||||
ypredictRidge <span style="color: #666666">=</span> RegRidge<span style="color: #666666">.</span>predict(X_test)
|
||||
MSEOwnRidgePredict[i] <span style="color: #666666">=</span> MSE(y_test,ypredictOwnRidge)
|
||||
MSERidgePredict[i] <span style="color: #666666">=</span> MSE(y_test,ypredictRidge)
|
||||
<span style="color: #008000">print</span>(<span style="color: #BA2121">"Beta values for own Ridge implementation"</span>)
|
||||
<span style="color: #008000">print</span>(OwnRidgeBeta) <span style="color: #408080; font-style: italic">#Intercept is given by mean of target variable</span>
|
||||
<span style="color: #008000">print</span>(<span style="color: #BA2121">"Beta values for Scikit-Learn Ridge implementation"</span>)
|
||||
<span style="color: #008000">print</span>(<span style="color: #BA2121">"Theta values for own Ridge implementation"</span>)
|
||||
<span style="color: #008000">print</span>(OwnRidgeTheta) <span style="color: #408080; font-style: italic">#Intercept is given by mean of target variable</span>
|
||||
<span style="color: #008000">print</span>(<span style="color: #BA2121">"Theta values for Scikit-Learn Ridge implementation"</span>)
|
||||
<span style="color: #008000">print</span>(RegRidge<span style="color: #666666">.</span>coef_)
|
||||
<span style="color: #008000">print</span>(<span style="color: #BA2121">'Intercept from own implementation:'</span>)
|
||||
<span style="color: #008000">print</span>(intercept_)
|
||||
@@ -3378,7 +3378,7 @@ plt<span style="color: #666666">.</span>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
|
||||
\( \beta_0 \) in the fitting. This applies to Lasso regularization as
|
||||
\( \theta_0 \) 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.
|
||||
</p>
|
||||
|
||||
Binary file not shown.
+502
-624
File diff suppressed because one or more lines are too long
@@ -23,10 +23,6 @@ o Goodfellow, Bengio and Courville, Deep Learning, chapter 2 on linear algebra
|
||||
o Raschka et al on preprocessing of data, relevant for exercise 3 this week, see chapter 4.
|
||||
o For exercise 1 of week 35, the book by A. Aldo Faisal, Cheng Soon Ong, and Marc Peter Deisenroth on the Mathematics of Machine Learning, may be very relevant. In particular chapter 5 at URL"https://mml-book.github.io/" (section 5.5 on derivatives) is very useful for exercise 1 this coming week.
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
!split
|
||||
===== Reminder from last week =====
|
||||
|
||||
@@ -2342,8 +2338,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 _Scikit-Learn_ (when not shrinking $\beta_0$) for the unknown parameters
|
||||
$\bm{\beta}$, are derived under the assumption that both $\bm{y}$ and
|
||||
from the library _Scikit-Learn_ (when not shrinking $\theta_0$) for the unknown parameters
|
||||
$\bm{\theta}$, are derived under the assumption that both $\bm{y}$ and
|
||||
$\bm{X}$ are zero centered, that is we subtract the mean values.
|
||||
|
||||
|
||||
@@ -2407,7 +2403,7 @@ simplicity, we will focus on ordinary regression, as done in the above example.
|
||||
The cost/loss function for regression is
|
||||
!bt
|
||||
\[
|
||||
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,.
|
||||
\]
|
||||
!et
|
||||
|
||||
@@ -2415,42 +2411,42 @@ 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.
|
||||
|
||||
What we have done is to single out the $\beta_0$ term in the
|
||||
What we have done is to single out the $\theta_0$ term in the
|
||||
definition of the mean squared error (MSE). The design matrix $X$
|
||||
does in this case not contain any intercept column. When we take the
|
||||
derivative with respect to $\beta_0$, we want the derivative to obey
|
||||
derivative with respect to $\theta_0$, we want the derivative to obey
|
||||
|
||||
!bt
|
||||
\[
|
||||
\frac{\partial C}{\partial \beta_j} = 0,
|
||||
\frac{\partial C}{\partial \theta_j} = 0,
|
||||
\]
|
||||
!et
|
||||
|
||||
for all $j$. For $\beta_0$ we have
|
||||
for all $j$. For $\theta_0$ we have
|
||||
|
||||
!bt
|
||||
\[
|
||||
\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).
|
||||
\]
|
||||
!et
|
||||
Multiplying away the constant $2/n$, we obtain
|
||||
!bt
|
||||
\[
|
||||
\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.
|
||||
\]
|
||||
!et
|
||||
|
||||
Let us specialize first to the case where we have only two parameters $\beta_0$ and $\beta_1$.
|
||||
Our result for $\beta_0$ simplifies then to
|
||||
Let us specialize first to the case where we have only two parameters $\theta_0$ and $\theta_1$.
|
||||
Our result for $\theta_0$ simplifies then to
|
||||
!bt
|
||||
\[
|
||||
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.
|
||||
\]
|
||||
!et
|
||||
We obtain then
|
||||
!bt
|
||||
\[
|
||||
\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}.
|
||||
\]
|
||||
!et
|
||||
If we define
|
||||
@@ -2468,20 +2464,20 @@ and the mean value of the outputs as
|
||||
we have
|
||||
!bt
|
||||
\[
|
||||
\beta_0 = \mu_y - \beta_1\mu_{\bm{x}_1}.
|
||||
\theta_0 = \mu_y - \theta_1\mu_{\bm{x}_1}.
|
||||
\]
|
||||
!et
|
||||
In the general case with more parameters than $\beta_0$ and $\beta_1$, we have
|
||||
In the general case with more parameters than $\theta_0$ and $\theta_1$, we have
|
||||
!bt
|
||||
\[
|
||||
\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.
|
||||
\]
|
||||
!et
|
||||
|
||||
We can rewrite the latter equation as
|
||||
!bt
|
||||
\[
|
||||
\beta_0 = \frac{1}{n}\sum_{i=0}^{n-1}y_i - \sum_{j=1}^{p-1} \mu_{\bm{x}_j}\beta_j,
|
||||
\theta_0 = \frac{1}{n}\sum_{i=0}^{n-1}y_i - \sum_{j=1}^{p-1} \mu_{\bm{x}_j}\theta_j,
|
||||
\]
|
||||
!et
|
||||
where we have defined
|
||||
@@ -2497,27 +2493,27 @@ the mean value for all elements of the column vector $\bm{x}_j$.
|
||||
Replacing $y_i$ with $y_i - y_i - \overline{\bm{y}}$ and centering also our design matrix results in a cost function (in vector-matrix disguise)
|
||||
!bt
|
||||
\[
|
||||
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}).
|
||||
\]
|
||||
!et
|
||||
|
||||
|
||||
|
||||
If we minimize with respect to $\bm{\beta}$ we have then
|
||||
If we minimize with respect to $\bm{\theta}$ we have then
|
||||
|
||||
!bt
|
||||
\[
|
||||
\hat{\bm{\beta}} = (\tilde{X}^T\tilde{X})^{-1}\tilde{X}^T\boldsymbol{\tilde{y}},
|
||||
\hat{\bm{\theta}} = (\tilde{X}^T\tilde{X})^{-1}\tilde{X}^T\boldsymbol{\tilde{y}},
|
||||
\]
|
||||
!et
|
||||
|
||||
where $\boldsymbol{\tilde{y}} = \boldsymbol{y} - \overline{\bm{y}}$
|
||||
and $\tilde{X}_{ij} = X_{ij} - \frac{1}{n}\sum_{k=0}^{n-1}X_{kj}$.
|
||||
|
||||
For Ridge regression we need to add $\lambda \boldsymbol{\beta}^T\boldsymbol{\beta}$ to the cost function and get then
|
||||
For Ridge regression we need to add $\lambda \boldsymbol{\theta}^T\boldsymbol{\theta}$ to the cost function and get then
|
||||
!bt
|
||||
\[
|
||||
\hat{\bm{\beta}} = (\tilde{X}^T\tilde{X} + \lambda I)^{-1}\tilde{X}^T\boldsymbol{\tilde{y}}.
|
||||
\hat{\bm{\theta}} = (\tilde{X}^T\tilde{X} + \lambda I)^{-1}\tilde{X}^T\boldsymbol{\tilde{y}}.
|
||||
\]
|
||||
!et
|
||||
|
||||
@@ -2541,15 +2537,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
|
||||
@@ -2559,15 +2555,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"True beta: {true_beta}")
|
||||
print(f"Fitted beta: {beta}")
|
||||
print(f"Sklearn fitted beta: {skl.coef_}")
|
||||
ypredictOwn = X @ beta
|
||||
print(f"True theta: {true_theta}")
|
||||
print(f"Fitted theta: {theta}")
|
||||
print(f"Sklearn fitted theta: {skl.coef_}")
|
||||
ypredictOwn = X @ theta
|
||||
ypredictSKL = skl.predict(X)
|
||||
print(f"MSE with intercept column")
|
||||
print(MSE(y,ypredictOwn))
|
||||
@@ -2577,7 +2573,7 @@ print(MSE(y,ypredictSKL))
|
||||
|
||||
plt.figure()
|
||||
plt.scatter(x, y, label="Data")
|
||||
plt.plot(x, X @ beta, label="Fit")
|
||||
plt.plot(x, X @ theta, label="Fit")
|
||||
plt.plot(x, skl.predict(X), label="Sklearn (fit_intercept=False)")
|
||||
|
||||
|
||||
@@ -2594,21 +2590,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"Manual intercept: {intercept}")
|
||||
print(f"Fitted beta (without intercept): {beta}")
|
||||
print(f"Fitted theta (without intercept): {theta}")
|
||||
print(f"Sklearn intercept: {skl.intercept_}")
|
||||
print(f"Sklearn fitted beta (without intercept): {skl.coef_}")
|
||||
ypredictOwn = X @ beta
|
||||
print(f"Sklearn fitted theta (without intercept): {skl.coef_}")
|
||||
ypredictOwn = X @ theta
|
||||
ypredictSKL = skl.predict(X)
|
||||
print(f"MSE with Manual intercept")
|
||||
print(MSE(y,ypredictOwn+intercept))
|
||||
print(f"MSE with Sklearn intercept")
|
||||
print(MSE(y,ypredictSKL))
|
||||
|
||||
plt.plot(x, X @ beta + intercept, "--", label="Fit (manual intercept)")
|
||||
plt.plot(x, X @ theta + intercept, "--", label="Fit (manual intercept)")
|
||||
plt.plot(x, skl.predict(X), "--", label="Sklearn (fit_intercept=True)")
|
||||
plt.grid()
|
||||
plt.legend()
|
||||
@@ -2625,26 +2621,26 @@ 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 $\beta_0$. For different values
|
||||
regularization term does not include $\theta_0$. For different values
|
||||
of $\lambda$, this may lead to different MSE values.
|
||||
|
||||
To remind the reader, the regularization term, with the intercept in Ridge regression, is given by
|
||||
!bt
|
||||
\[
|
||||
\lambda \vert\vert \bm{\beta} \vert\vert_2^2 = \lambda \sum_{j=0}^{p-1}\beta_j^2,
|
||||
\lambda \vert\vert \bm{\theta} \vert\vert_2^2 = \lambda \sum_{j=0}^{p-1}\theta_j^2,
|
||||
\]
|
||||
!et
|
||||
but when we take out the intercept, this equation becomes
|
||||
!bt
|
||||
\[
|
||||
\lambda \vert\vert \bm{\beta} \vert\vert_2^2 = \lambda \sum_{j=1}^{p-1}\beta_j^2.
|
||||
\lambda \vert\vert \bm{\theta} \vert\vert_2^2 = \lambda \sum_{j=1}^{p-1}\theta_j^2.
|
||||
\]
|
||||
!et
|
||||
|
||||
For Lasso regression we have
|
||||
!bt
|
||||
\[
|
||||
\lambda \vert\vert \bm{\beta} \vert\vert_1 = \lambda \sum_{j=1}^{p-1}\vert\beta_j\vert.
|
||||
\lambda \vert\vert \bm{\theta} \vert\vert_1 = \lambda \sum_{j=1}^{p-1}\vert\theta_j\vert.
|
||||
\]
|
||||
!et
|
||||
|
||||
@@ -2695,20 +2691,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("Beta values for own Ridge implementation")
|
||||
print(OwnRidgeBeta)
|
||||
print("Beta values for Scikit-Learn Ridge implementation")
|
||||
print("Theta values for own Ridge implementation")
|
||||
print(OwnRidgeTheta)
|
||||
print("Theta values for Scikit-Learn Ridge implementation")
|
||||
print(RegRidge.coef_)
|
||||
print("MSE values for own Ridge implementation")
|
||||
print(MSEOwnRidgePredict[i])
|
||||
@@ -2780,18 +2776,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("Beta values for own Ridge implementation")
|
||||
print(OwnRidgeBeta) #Intercept is given by mean of target variable
|
||||
print("Beta values for Scikit-Learn Ridge implementation")
|
||||
print("Theta values for own Ridge implementation")
|
||||
print(OwnRidgeTheta) #Intercept is given by mean of target variable
|
||||
print("Theta values for Scikit-Learn Ridge implementation")
|
||||
print(RegRidge.coef_)
|
||||
print('Intercept from own implementation:')
|
||||
print(intercept_)
|
||||
@@ -2815,7 +2811,7 @@ plt.show()
|
||||
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
|
||||
$\beta_0$ in the fitting. This applies to Lasso regularization as
|
||||
$\theta_0$ 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.
|
||||
|
||||
|
||||
Reference in New Issue
Block a user