updated txt

This commit is contained in:
Morten Hjorth-Jensen
2021-09-17 07:29:28 +02:00
parent 2eba427d80
commit 1800607d04
32 changed files with 2426 additions and 267 deletions
+91 -6
View File
@@ -2571,16 +2571,90 @@ plt.show()
!ec
How can we understand this? _More text to be added_.
How can we understand this?
Let us write out the values of the coefficients $\beta_i$ as functions
of the polynomial degree and noise. We will focus only on the Ridge
results and some few selected values of the hyperparameter $\lambda$.
If we don't include any noise and run this code for different values
of the polynomial degree, we notice that the results for $\beta_i$ do
not show great changes from one order to the next. This is an
indication that for higher polynomial orders, our parameters become
less important.
If we however add noise, what happens is that the polynomial fit is
trying to adjust the fit to traverse in the best possible way all data
points. This can lead to large fluctuations in the parameters
$\beta_i$ as functions of polynomial order. It will also be reflected
in a larger value of the variance of each parameter $\beta_i$. What
Ridge regression (and Lasso as well) are doing then is to try to
quench the fluctuations in the parameters of $\beta_i$ which have a
large variance (normally for higher orders in the polynomial).
!bc pycod
import numpy as np
import pandas as pd
from IPython.display import display
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn import linear_model
# Make data set.
n = 1000
x = np.random.rand(n)
y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)+ np.random.randn(n)
Maxpolydegree = 5
X = np.zeros((len(x),Maxpolydegree))
X[:,0] = 1.0
for polydegree in range(1, Maxpolydegree):
for degree in range(polydegree):
X[:,degree] = x**(degree)
# We split the data in test and training data
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)
# Decide which values of lambda to use
nlambdas = 5
lambdas = np.logspace(-3, 2, nlambdas)
for i in range(nlambdas):
lmb = lambdas[i]
# Make the fit using Ridge only
RegRidge = linear_model.Ridge(lmb,fit_intercept=False)
RegRidge.fit(X_train,y_train)
# and then make the prediction
ypredictRidge = RegRidge.predict(X_test)
Coeffs = np.array(RegRidge.coef_)
BetaValues = pd.DataFrame(Coeffs)
BetaValues.columns = ['beta']
display(BetaValues)
!ec
As an exercise, repeat these calculations with ordinary least squares
only with and without noise. Calculate thereafter the variance of the
parameters $\beta_j$ as function of polynomial order and of the added
noise. Here we recommend to use $\sigma^2=1$ as variance for the
added noise (which follows a normal distribution with mean value zero).
Comment your results. If you have a large noise term, do the parameters $\beta_j$ vary more as function
model complexity? And what about their variance?
===== Linking Bayes' Theorem with Ridge and Lasso Regression =====
We have seen that Ridge regression suppresses those features which
have a small singular value. This corresponds to a feature which exhibits
a large variance in the parameters $\beta_j$.
Our analysis hitherto has been based on linear algebra. To add to our intuition, we will use
Bayes' theorem in order to deepen our understanding of Ridge and Lasso regression.
Using Bayes' theorem we can gain a better intuition about Ridge and Lasso regression.
For ordinary least squares we postulated that the maximum likelihood for the doamin of events $\bm{D}$ (one-dimensional case)
For ordinary least squares we postulated that the maximum likelihood for the domain of events $\bm{D}$ (one-dimensional case)
!bt
\[
\bm{D}=[(x_0,y_0), (x_1,y_1),\dots, (x_{n-1},y_{n-1})],
@@ -2616,7 +2690,10 @@ With the posterior probability defined by a likelihood which we have
already modeled and an unknown prior, we are now ready to make
additional models for the prior.
We can, based on our discussions of the variance of $\bm{\beta}$ and the mean value, assume that the prior for the values $\bm{\beta}$ is given by a Gaussian with mean value zero and variance $\tau^2$, that is
We can, based on our discussions of the variance of $\bm{\beta}$ and
the mean value, assume that the prior for the values $\bm{\beta}$ is
given by a Gaussian with mean value zero and variance $\tau^2$, that
is
!bt
\[
@@ -2635,7 +2712,7 @@ p(\bm{\beta\vert\bm{D})}=\prod_{i=0}^{n-1}\frac{1}{\sqrt{2\pi\sigma^2}}\exp{\lef
We can now optimize this quantity with respect to $\bm{\beta}$. As we
did for OLS, this is most conveniently done by taking the negative
logarithm of the posterior probability. Doing so and leaving out the
constants terms that do not depend on $\beta$, we have
terms that do not depend on $\beta$, we have
!bt
@@ -2689,4 +2766,12 @@ C(\bm{\beta}=\frac{\vert\vert (\bm{y}-\bm{X}\bm{\beta})\vert\vert_2^2}{2\sigma^2
which is our Lasso cost function!
Plotting these prior functions shows us that we can use the parameter
$\lambda$ to shrink or increase the role of a given parameter
$\beta_j$. The variance for the Laplace distribution is
$2\tau^2=1/\lambda$ while for the Gaussian distribution it is
$\sigma^2=1/(2\lambda)$. Thus, increasing the variance means
decreasing $\lambda$ and shrinking the variance means increasing
$\lamdbda$. When we increase $\lambda$, this corresponds to shrinking the role of less important features (small singular values).
+144 -2
View File
@@ -134,7 +134,7 @@ With the OLS expressions for the parameters $\bm{\beta}$ we can evaluate the exp
\]
!et
This means that the estimator of the regression parameters is unbiased.
v
We can also calculate the variance
The variance of $\bm{\beta}$ is
@@ -363,7 +363,148 @@ histogram of the relative frequency of $\widehat{\beta}^*$. Instead
you use the estimators corresponding to the statistic of interest. For
example, if you are interested in estimating the variance of $\widehat
\beta$, apply the etsimator $\widehat \sigma^2$ to the values
$\widehat \beta ^*$.
$\widehat \beta^*$.
Before we proceed however, we need to remind ourselves about a central
theorem in statistics, namely the so-called _central limit theorem_.
This theorem plays a central role in understanding why the Bootstrap
(and other resampling methods) work so well on independent and
identically distributed variables.
Suppose we have a PDF $p(x)$ from which we generate a series $N$
of averages $\langle x_i \rangle$. Each mean value $\langle x_i \rangle$
is viewed as the average of a specific measurement, e.g., throwing
dice 100 times and then taking the average value, or producing a certain
amount of random numbers.
For notational ease, we set $\langle x_i \rangle=x_i$ in the discussion
which follows.
If we compute the mean $z$ of $m$ such mean values $x_i$
!bt
\[
z=\frac{x_1+x_2+\dots+x_m}{m},
\]
!et
the question we pose is which is the PDF of the new variable $z$.
The probability of obtaining an average value $z$ is the product of the
probabilities of obtaining arbitrary individual mean values $x_i$,
but with the constraint that the average is $z$. We can express this through
the following expression
!bt
\[
\tilde{p}(z)=\int dx_1p(x_1)\int dx_2p(x_2)\dots\int dx_mp(x_m)
\delta(z-\frac{x_1+x_2+\dots+x_m}{m}),
\]
!et
where the $\delta$-function enbodies the constraint that the mean is $z$.
All measurements that lead to each individual $x_i$ are expected to
be independent, which in turn means that we can express $\tilde{p}$ as the
product of individual $p(x_i)$. The independence assumption is important in the derivation of the central limit theorem.
If we use the integral expression for the $\delta$-function
!bt
\[
\delta(z-\frac{x_1+x_2+\dots+x_m}{m})=\frac{1}{2\pi}\int_{-\infty}^{\infty}
dq\exp{\left(iq(z-\frac{x_1+x_2+\dots+x_m}{m})\right)},
\]
!et
and inserting $e^{i\mu q-i\mu q}$ where $\mu$ is the mean value
we arrive at
!bt
\[
\tilde{p}(z)=\frac{1}{2\pi}\int_{-\infty}^{\infty}
dq\exp{\left(iq(z-\mu)\right)}\left[\int_{-\infty}^{\infty}
dxp(x)\exp{\left(iq(\mu-x)/m\right)}\right]^m,
\]
!et
with the integral over $x$ resulting in
!bt
\[
\int_{-\infty}^{\infty}dxp(x)\exp{\left(iq(\mu-x)/m\right)}=
\int_{-\infty}^{\infty}dxp(x)
\left[1+\frac{iq(\mu-x)}{m}-\frac{q^2(\mu-x)^2}{2m^2}+\dots\right].
\]
!et
The second term on the rhs disappears since this is just the mean and
employing the definition of $\sigma^2$ we have
!bt
\[
\int_{-\infty}^{\infty}dxp(x)e^{\left(iq(\mu-x)/m\right)}=
1-\frac{q^2\sigma^2}{2m^2}+\dots,
\]
!et
resulting in
!bt
\[
\left[\int_{-\infty}^{\infty}dxp(x)\exp{\left(iq(\mu-x)/m\right)}\right]^m\approx
\left[1-\frac{q^2\sigma^2}{2m^2}+\dots \right]^m,
\]
!et
and in the limit $m\rightarrow \infty$ we obtain
!bt
\[
\tilde{p}(z)=\frac{1}{\sqrt{2\pi}(\sigma/\sqrt{m})}
\exp{\left(-\frac{(z-\mu)^2}{2(\sigma/\sqrt{m})^2}\right)},
\]
!et
which is the normal distribution with variance
$\sigma^2_m=\sigma^2/m$, where $\sigma$ is the variance of the PDF $p(x)$
and $\mu$ is also the mean of the PDF $p(x)$.
Thus, the central limit theorem states that the PDF $\tilde{p}(z)$ of
the average of $m$ random values corresponding to a PDF $p(x)$
is a normal distribution whose mean is the
mean value of the PDF $p(x)$ and whose variance is the variance
of the PDF $p(x)$ divided by $m$, the number of values used to compute $z$.
The central limit theorem leads to the well-known expression for the
standard deviation, given by
!bt
\[
\sigma_m=
\frac{\sigma}{\sqrt{m}}.
\]
!et
The latter is true only if the average value is known exactly. This is obtained in the limit
$m\rightarrow \infty$ only. Because the mean and the variance are measured quantities we obtain
the familiar expression in statistics
!bt
\[
\sigma_m\approx
\frac{\sigma}{\sqrt{m-1}}.
\]
!et
In many cases however the above estimate for the standard deviation,
in particular if correlations are strong, may be too simplistic. Keep
in mind that we have assumed that the variables $x$ are independent
and identically distributed. This is obviously not always the
case. For example, the random numbers (or better pseudorandom numbers)
we generate in various calculations do always exhibit some
correlations.
The theorem is satisfied by a large class of PDFs. Note however that for a
finite $m$, it is not always possible to find a closed form /analytic expression for
$\tilde{p}(x)$.
@@ -1337,3 +1478,4 @@ At the end, you should present a critical evaluation of your results
and discuss the applicability of these regression methods to the type
of data presented here (either the terrain data we propose or other data sets).
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.

After

Width:  |  Height:  |  Size: 23 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 16 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 37 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 16 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 14 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 12 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 10 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 32 KiB

@@ -4079,16 +4079,101 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"How can we understand this? **More text to be added**.\n",
"How can we understand this?\n",
"\n",
"Let us write out the values of the coefficients $\\beta_i$ as functions\n",
"of the polynomial degree and noise. We will focus only on the Ridge\n",
"results and some few selected values of the hyperparameter $\\lambda$.\n",
"\n",
"If we don't include any noise and run this code for different values\n",
"of the polynomial degree, we notice that the results for $\\beta_i$ do\n",
"not show great changes from one order to the next. This is an\n",
"indication that for higher polynomial orders, our parameters become\n",
"less important.\n",
"\n",
"If we however add noise, what happens is that the polynomial fit is\n",
"trying to adjust the fit to traverse in the best possible way all data\n",
"points. This can lead to large fluctuations in the parameters\n",
"$\\beta_i$ as functions of polynomial order. It will also be reflected\n",
"in a larger value of the variance of each parameter $\\beta_i$. What\n",
"Ridge regression (and Lasso as well) are doing then is to try to\n",
"quench the fluctuations in the parameters of $\\beta_i$ which have a\n",
"large variance (normally for higher orders in the polynomial)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false,
"editable": true
},
"outputs": [],
"source": [
"import numpy as np\n",
"import pandas as pd\n",
"from IPython.display import display\n",
"import matplotlib.pyplot as plt\n",
"from sklearn.model_selection import train_test_split\n",
"from sklearn import linear_model\n",
"\n",
"# Make data set.\n",
"n = 1000\n",
"x = np.random.rand(n)\n",
"y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)+ np.random.randn(n)\n",
"\n",
"Maxpolydegree = 5\n",
"X = np.zeros((len(x),Maxpolydegree))\n",
"X[:,0] = 1.0\n",
"\n",
"for polydegree in range(1, Maxpolydegree):\n",
" for degree in range(polydegree):\n",
" X[:,degree] = x**(degree)\n",
"\n",
"\n",
"# We split the data in test and training data\n",
"X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)\n",
"\n",
"# Decide which values of lambda to use\n",
"nlambdas = 5\n",
"lambdas = np.logspace(-3, 2, nlambdas)\n",
"for i in range(nlambdas):\n",
" lmb = lambdas[i]\n",
" # Make the fit using Ridge only\n",
" RegRidge = linear_model.Ridge(lmb,fit_intercept=False)\n",
" RegRidge.fit(X_train,y_train)\n",
" # and then make the prediction\n",
" ypredictRidge = RegRidge.predict(X_test)\n",
" Coeffs = np.array(RegRidge.coef_)\n",
" BetaValues = pd.DataFrame(Coeffs)\n",
" BetaValues.columns = ['beta']\n",
" display(BetaValues)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As an exercise, repeat these calculations with ordinary least squares\n",
"only with and without noise. Calculate thereafter the variance of the\n",
"parameters $\\beta_j$ as function of polynomial order and of the added\n",
"noise. Here we recommend to use $\\sigma^2=1$ as variance for the\n",
"added noise (which follows a normal distribution with mean value zero).\n",
"Comment your results. If you have a large noise term, do the parameters $\\beta_j$ vary more as function\n",
"model complexity? And what about their variance? \n",
"\n",
"\n",
"\n",
"\n",
"## Linking Bayes' Theorem with Ridge and Lasso Regression\n",
"\n",
"We have seen that Ridge regression suppresses those features which\n",
"have a small singular value. This corresponds to a feature which exhibits\n",
"a large variance in the parameters $\\beta_j$.\n",
"Our analysis hitherto has been based on linear algebra. To add to our intuition, we will use\n",
"Bayes' theorem in order to deepen our understanding of Ridge and Lasso regression. \n",
"\n",
"Using Bayes' theorem we can gain a better intuition about Ridge and Lasso regression. \n",
"\n",
"For ordinary least squares we postulated that the maximum likelihood for the doamin of events $\\boldsymbol{D}$ (one-dimensional case)"
"For ordinary least squares we postulated that the maximum likelihood for the domain of events $\\boldsymbol{D}$ (one-dimensional case)"
]
},
{
@@ -4160,7 +4245,10 @@
"already modeled and an unknown prior, we are now ready to make\n",
"additional models for the prior.\n",
"\n",
"We can, based on our discussions of the variance of $\\boldsymbol{\\beta}$ and the mean value, assume that the prior for the values $\\boldsymbol{\\beta}$ is given by a Gaussian with mean value zero and variance $\\tau^2$, that is"
"We can, based on our discussions of the variance of $\\boldsymbol{\\beta}$ and\n",
"the mean value, assume that the prior for the values $\\boldsymbol{\\beta}$ is\n",
"given by a Gaussian with mean value zero and variance $\\tau^2$, that\n",
"is"
]
},
{
@@ -4195,7 +4283,7 @@
"We can now optimize this quantity with respect to $\\boldsymbol{\\beta}$. As we\n",
"did for OLS, this is most conveniently done by taking the negative\n",
"logarithm of the posterior probability. Doing so and leaving out the\n",
"constants terms that do not depend on $\\beta$, we have"
"terms that do not depend on $\\beta$, we have"
]
},
{
@@ -4296,7 +4384,16 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"which is our Lasso cost function!"
"which is our Lasso cost function! \n",
"\n",
"\n",
"Plotting these prior functions shows us that we can use the parameter\n",
"$\\lambda$ to shrink or increase the role of a given parameter\n",
"$\\beta_j$. The variance for the Laplace distribution is\n",
"$2\\tau^2=1/\\lambda$ while for the Gaussian distribution it is\n",
"$\\sigma^2=1/(2\\lambda)$. Thus, increasing the variance means\n",
"decreasing $\\lambda$ and shrinking the variance means increasing\n",
"$\\lamdbda$. When we increase $\\lambda$, this corresponds to shrinking the role of less important features (small singular values)."
]
}
],
@@ -211,7 +211,7 @@
"metadata": {},
"source": [
"This means that the estimator of the regression parameters is unbiased.\n",
"\n",
"v\n",
"We can also calculate the variance\n",
"\n",
"The variance of $\\boldsymbol{\\beta}$ is"
@@ -503,7 +503,239 @@
"you use the estimators corresponding to the statistic of interest. For\n",
"example, if you are interested in estimating the variance of $\\widehat\n",
"\\beta$, apply the etsimator $\\widehat \\sigma^2$ to the values\n",
"$\\widehat \\beta ^*$.\n",
"$\\widehat \\beta^*$.\n",
"\n",
"Before we proceed however, we need to remind ourselves about a central\n",
"theorem in statistics, namely the so-called **central limit theorem**.\n",
"This theorem plays a central role in understanding why the Bootstrap\n",
"(and other resampling methods) work so well on independent and\n",
"identically distributed variables.\n",
"\n",
"\n",
"Suppose we have a PDF $p(x)$ from which we generate a series $N$\n",
"of averages $\\langle x_i \\rangle$. Each mean value $\\langle x_i \\rangle$\n",
"is viewed as the average of a specific measurement, e.g., throwing \n",
"dice 100 times and then taking the average value, or producing a certain\n",
"amount of random numbers. \n",
"For notational ease, we set $\\langle x_i \\rangle=x_i$ in the discussion\n",
"which follows. \n",
"\n",
"If we compute the mean $z$ of $m$ such mean values $x_i$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"z=\\frac{x_1+x_2+\\dots+x_m}{m},\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"the question we pose is which is the PDF of the new variable $z$.\n",
"\n",
"\n",
"The probability of obtaining an average value $z$ is the product of the \n",
"probabilities of obtaining arbitrary individual mean values $x_i$,\n",
"but with the constraint that the average is $z$. We can express this through\n",
"the following expression"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\tilde{p}(z)=\\int dx_1p(x_1)\\int dx_2p(x_2)\\dots\\int dx_mp(x_m)\n",
" \\delta(z-\\frac{x_1+x_2+\\dots+x_m}{m}),\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"where the $\\delta$-function enbodies the constraint that the mean is $z$.\n",
"All measurements that lead to each individual $x_i$ are expected to\n",
"be independent, which in turn means that we can express $\\tilde{p}$ as the \n",
"product of individual $p(x_i)$. The independence assumption is important in the derivation of the central limit theorem.\n",
"\n",
"\n",
"\n",
"If we use the integral expression for the $\\delta$-function"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\delta(z-\\frac{x_1+x_2+\\dots+x_m}{m})=\\frac{1}{2\\pi}\\int_{-\\infty}^{\\infty}\n",
" dq\\exp{\\left(iq(z-\\frac{x_1+x_2+\\dots+x_m}{m})\\right)},\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"and inserting $e^{i\\mu q-i\\mu q}$ where $\\mu$ is the mean value\n",
"we arrive at"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\tilde{p}(z)=\\frac{1}{2\\pi}\\int_{-\\infty}^{\\infty}\n",
" dq\\exp{\\left(iq(z-\\mu)\\right)}\\left[\\int_{-\\infty}^{\\infty}\n",
" dxp(x)\\exp{\\left(iq(\\mu-x)/m\\right)}\\right]^m,\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"with the integral over $x$ resulting in"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\int_{-\\infty}^{\\infty}dxp(x)\\exp{\\left(iq(\\mu-x)/m\\right)}=\n",
" \\int_{-\\infty}^{\\infty}dxp(x)\n",
" \\left[1+\\frac{iq(\\mu-x)}{m}-\\frac{q^2(\\mu-x)^2}{2m^2}+\\dots\\right].\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The second term on the rhs disappears since this is just the mean and \n",
"employing the definition of $\\sigma^2$ we have"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\int_{-\\infty}^{\\infty}dxp(x)e^{\\left(iq(\\mu-x)/m\\right)}=\n",
" 1-\\frac{q^2\\sigma^2}{2m^2}+\\dots,\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"resulting in"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\left[\\int_{-\\infty}^{\\infty}dxp(x)\\exp{\\left(iq(\\mu-x)/m\\right)}\\right]^m\\approx\n",
" \\left[1-\\frac{q^2\\sigma^2}{2m^2}+\\dots \\right]^m,\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"and in the limit $m\\rightarrow \\infty$ we obtain"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\tilde{p}(z)=\\frac{1}{\\sqrt{2\\pi}(\\sigma/\\sqrt{m})}\n",
" \\exp{\\left(-\\frac{(z-\\mu)^2}{2(\\sigma/\\sqrt{m})^2}\\right)},\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"which is the normal distribution with variance\n",
"$\\sigma^2_m=\\sigma^2/m$, where $\\sigma$ is the variance of the PDF $p(x)$\n",
"and $\\mu$ is also the mean of the PDF $p(x)$. \n",
"\n",
"\n",
"Thus, the central limit theorem states that the PDF $\\tilde{p}(z)$ of\n",
"the average of $m$ random values corresponding to a PDF $p(x)$ \n",
"is a normal distribution whose mean is the \n",
"mean value of the PDF $p(x)$ and whose variance is the variance\n",
"of the PDF $p(x)$ divided by $m$, the number of values used to compute $z$.\n",
"\n",
"The central limit theorem leads to the well-known expression for the\n",
"standard deviation, given by"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\sigma_m=\n",
"\\frac{\\sigma}{\\sqrt{m}}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The latter is true only if the average value is known exactly. This is obtained in the limit\n",
"$m\\rightarrow \\infty$ only. Because the mean and the variance are measured quantities we obtain \n",
"the familiar expression in statistics"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\sigma_m\\approx \n",
"\\frac{\\sigma}{\\sqrt{m-1}}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In many cases however the above estimate for the standard deviation,\n",
"in particular if correlations are strong, may be too simplistic. Keep\n",
"in mind that we have assumed that the variables $x$ are independent\n",
"and identically distributed. This is obviously not always the\n",
"case. For example, the random numbers (or better pseudorandom numbers)\n",
"we generate in various calculations do always exhibit some\n",
"correlations.\n",
"\n",
"\n",
"\n",
"The theorem is satisfied by a large class of PDFs. Note however that for a\n",
"finite $m$, it is not always possible to find a closed form /analytic expression for\n",
"$\\tilde{p}(x)$.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
+371 -63
View File
@@ -157,6 +157,11 @@
6. Logistic Regression
</a>
</li>
<li class="toctree-l1">
<a class="reference internal" href="chapteroptimization.html">
7. Optimization, the central part of any Machine Learning algortithm
</a>
</li>
<li class="toctree-l1">
<a class="reference internal" href="chapter5.html">
8. Support Vector Machines, overarching aims
@@ -1101,10 +1106,10 @@ covariance matrix through the <strong>np.linalg.eig()</strong> function.</p>
</div>
</div>
<div class="cell_output docutils container">
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>-0.1831277634928002
3.3041320306136366
[[1.03929932 3.07670437]
[3.07670437 9.7701384 ]]
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>-0.055529303095955385
3.9154537458386387
[[1.0003451 2.94327895]
[2.94327895 9.80021057]]
</pre></div>
</div>
</div>
@@ -1141,10 +1146,10 @@ a more brute force way. Here we scale the mean values for each column of the des
</div>
</div>
<div class="cell_output docutils container">
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>0.09456011349477329
1.6600855222624895
[[1. 0.70183798]
[0.70183798 1. ]]
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>0.095786152583691
1.5864689451163851
[[1. 0.6713619]
[0.6713619 1. ]]
</pre></div>
</div>
</div>
@@ -1174,30 +1179,30 @@ this matrix we easily see that it is a positive definite matrix.</p>
</div>
</div>
<div class="cell_output docutils container">
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>[[ 0.54316925 0.40968888]
[-1.8393131 -5.8788431 ]
[ 0.24613822 0.9957744 ]
[ 1.38892672 3.93679587]
[ 0.06540809 -0.84008474]
[ 0.49841285 2.31081134]
[-0.49323032 -1.25263916]
[ 0.08380813 -0.7925146 ]
[ 0.64078247 2.92814088]
[-1.13410232 -1.81712976]]
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>[[-0.10979444 1.25912505]
[ 2.17273709 7.48017006]
[ 0.31520842 0.66152576]
[-0.07076926 -0.44655382]
[ 0.24348403 0.38561052]
[-0.57984245 -0.86011441]
[ 0.42375621 -1.21244261]
[-0.2407922 0.56831157]
[-2.23691315 -6.6936767 ]
[ 0.08292574 -1.14195542]]
0 1
0 0.543169 0.409689
1 -1.839313 -5.878843
2 0.246138 0.995774
3 1.388927 3.936796
4 0.065408 -0.840085
5 0.498413 2.310811
6 -0.493230 -1.252639
7 0.083808 -0.792515
8 0.640782 2.928141
9 -1.134102 -1.817130
0 -0.109794 1.259125
1 2.172737 7.480170
2 0.315208 0.661526
3 -0.070769 -0.446554
4 0.243484 0.385611
5 -0.579842 -0.860114
6 0.423756 -1.212443
7 -0.240792 0.568312
8 -2.236913 -6.693677
9 0.082926 -1.141955
0 1
0 1.000000 0.942726
1 0.942726 1.000000
0 1.000000 0.931066
1 0.931066 1.000000
</pre></div>
</div>
</div>
@@ -1254,37 +1259,37 @@ this matrix we easily see that it is a positive definite matrix.</p>
<div class="cell_output docutils container">
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span> 0 1 2 3 4 5 6 7 \
0 0.0 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
1 0.0 0.080633 0.085811 0.085748 0.084826 0.083221 0.080341 0.078342
2 0.0 0.085811 0.092659 0.093248 0.092932 0.091659 0.088414 0.086561
3 0.0 0.085748 0.093248 0.096639 0.096802 0.095941 0.093928 0.092288
4 0.0 0.084826 0.092932 0.096802 0.097346 0.096763 0.094774 0.093328
5 0.0 0.083221 0.091659 0.095941 0.096763 0.096406 0.094513 0.093236
6 0.0 0.080341 0.088414 0.093928 0.094774 0.094513 0.093573 0.092375
7 0.0 0.078342 0.086561 0.092288 0.093328 0.093236 0.092375 0.091318
8 0.0 0.076213 0.084461 0.090355 0.091535 0.091578 0.090817 0.089881
9 0.0 0.074037 0.082239 0.088258 0.089540 0.089693 0.089039 0.088209
10 0.0 0.073546 0.081474 0.088137 0.089326 0.089437 0.089352 0.088484
11 0.0 0.071453 0.079337 0.086052 0.087334 0.087544 0.087526 0.086757
12 0.0 0.069371 0.077163 0.083905 0.085253 0.085544 0.085598 0.084916
13 0.0 0.067326 0.074996 0.081744 0.083138 0.083495 0.083622 0.083015
14 0.0 0.065334 0.072862 0.079597 0.081024 0.081435 0.081632 0.081092
1 0.0 0.090107 0.080265 0.095180 0.087148 0.079801 0.089928 0.082666
2 0.0 0.080265 0.073052 0.087032 0.080565 0.074557 0.083797 0.077645
3 0.0 0.095180 0.087032 0.106946 0.099079 0.091756 0.104887 0.097150
4 0.0 0.087148 0.080565 0.099079 0.092378 0.086076 0.098028 0.091253
5 0.0 0.079801 0.074557 0.091756 0.086076 0.080678 0.091542 0.085627
6 0.0 0.089928 0.083797 0.104887 0.098028 0.091542 0.105387 0.098187
7 0.0 0.082666 0.077645 0.097150 0.091253 0.085627 0.098187 0.091855
8 0.0 0.076121 0.072046 0.090099 0.085041 0.080173 0.091566 0.086006
9 0.0 0.070212 0.066949 0.083672 0.079351 0.075151 0.085484 0.080612
10 0.0 0.083337 0.078717 0.099609 0.093725 0.088085 0.101781 0.095275
11 0.0 0.076774 0.072995 0.092276 0.087205 0.082303 0.094708 0.088983
12 0.0 0.070866 0.067811 0.085627 0.081270 0.077021 0.088260 0.083228
13 0.0 0.065545 0.063114 0.079598 0.075869 0.072197 0.082382 0.077968
14 0.0 0.060748 0.058856 0.074130 0.070955 0.067793 0.077026 0.073162
8 9 10 11 12 13 14
0 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
1 0.076213 0.074037 0.073546 0.071453 0.069371 0.067326 0.065334
2 0.084461 0.082239 0.081474 0.079337 0.077163 0.074996 0.072862
3 0.090355 0.088258 0.088137 0.086052 0.083905 0.081744 0.079597
4 0.091535 0.089540 0.089326 0.087334 0.085253 0.083138 0.081024
5 0.091578 0.089693 0.089437 0.087544 0.085544 0.083495 0.081435
6 0.090817 0.089039 0.089352 0.087526 0.085598 0.083622 0.081632
7 0.089881 0.088209 0.088484 0.086757 0.084916 0.083015 0.081092
8 0.088556 0.086985 0.087242 0.085610 0.083855 0.082032 0.080182
9 0.086985 0.085510 0.085763 0.084221 0.082549 0.080806 0.079029
10 0.087242 0.085763 0.086430 0.084864 0.083177 0.081425 0.079643
11 0.085610 0.084221 0.084864 0.083384 0.081778 0.080101 0.078390
12 0.083855 0.082549 0.083177 0.081778 0.080249 0.078646 0.077005
13 0.082032 0.080806 0.081425 0.080101 0.078646 0.077115 0.075542
14 0.080182 0.079029 0.079643 0.078390 0.077005 0.075542 0.074036
1 0.076121 0.070212 0.083337 0.076774 0.070866 0.065545 0.060748
2 0.072046 0.066949 0.078717 0.072995 0.067811 0.063114 0.058856
3 0.090099 0.083672 0.099609 0.092276 0.085627 0.079598 0.074130
4 0.085041 0.079351 0.093725 0.087205 0.081270 0.075869 0.070955
5 0.080173 0.075151 0.088085 0.082303 0.077021 0.072197 0.067793
6 0.091566 0.085484 0.101781 0.094708 0.088260 0.082382 0.077026
7 0.086006 0.080612 0.095275 0.088983 0.083228 0.077968 0.073162
8 0.080846 0.076068 0.089250 0.083659 0.078529 0.073827 0.069519
9 0.076068 0.071841 0.083682 0.078720 0.074154 0.069956 0.066098
10 0.089250 0.083682 0.099504 0.092936 0.086920 0.081414 0.076377
11 0.083659 0.078720 0.092936 0.087096 0.081732 0.076810 0.072296
12 0.078529 0.074154 0.086920 0.081732 0.076954 0.072557 0.068515
13 0.073827 0.069956 0.081414 0.076810 0.072557 0.068633 0.065016
14 0.069519 0.066098 0.076377 0.072296 0.068515 0.065016 0.061783
</pre></div>
</div>
</div>
@@ -2864,12 +2869,305 @@ Test MSE OLS
<img alt="_images/chapter2_316_1.png" src="_images/chapter2_316_1.png" />
</div>
</div>
<p>How can we understand this? <strong>More text to be added</strong>.</p>
<p>How can we understand this?</p>
<p>Let us write out the values of the coefficients <span class="math notranslate nohighlight">\(\beta_i\)</span> as functions
of the polynomial degree and noise. We will focus only on the Ridge
results and some few selected values of the hyperparameter <span class="math notranslate nohighlight">\(\lambda\)</span>.</p>
<p>If we dont include any noise and run this code for different values
of the polynomial degree, we notice that the results for <span class="math notranslate nohighlight">\(\beta_i\)</span> do
not show great changes from one order to the next. This is an
indication that for higher polynomial orders, our parameters become
less important.</p>
<p>If we however add noise, what happens is that the polynomial fit is
trying to adjust the fit to traverse in the best possible way all data
points. This can lead to large fluctuations in the parameters
<span class="math notranslate nohighlight">\(\beta_i\)</span> as functions of polynomial order. It will also be reflected
in a larger value of the variance of each parameter <span class="math notranslate nohighlight">\(\beta_i\)</span>. What
Ridge regression (and Lasso as well) are doing then is to try to
quench the fluctuations in the parameters of <span class="math notranslate nohighlight">\(\beta_i\)</span> which have a
large variance (normally for higher orders in the polynomial).</p>
<div class="cell docutils container">
<div class="cell_input docutils container">
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
<span class="kn">import</span> <span class="nn">pandas</span> <span class="k">as</span> <span class="nn">pd</span>
<span class="kn">from</span> <span class="nn">IPython.display</span> <span class="kn">import</span> <span class="n">display</span>
<span class="kn">import</span> <span class="nn">matplotlib.pyplot</span> <span class="k">as</span> <span class="nn">plt</span>
<span class="kn">from</span> <span class="nn">sklearn.model_selection</span> <span class="kn">import</span> <span class="n">train_test_split</span>
<span class="kn">from</span> <span class="nn">sklearn</span> <span class="kn">import</span> <span class="n">linear_model</span>
<span class="c1"># Make data set.</span>
<span class="n">n</span> <span class="o">=</span> <span class="mi">1000</span>
<span class="n">x</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">random</span><span class="o">.</span><span class="n">rand</span><span class="p">(</span><span class="n">n</span><span class="p">)</span>
<span class="n">y</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="o">-</span><span class="n">x</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span> <span class="o">+</span> <span class="mf">1.5</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="o">-</span><span class="p">(</span><span class="n">x</span><span class="o">-</span><span class="mi">2</span><span class="p">)</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span><span class="o">+</span> <span class="n">np</span><span class="o">.</span><span class="n">random</span><span class="o">.</span><span class="n">randn</span><span class="p">(</span><span class="n">n</span><span class="p">)</span>
<span class="n">Maxpolydegree</span> <span class="o">=</span> <span class="mi">5</span>
<span class="n">X</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">((</span><span class="nb">len</span><span class="p">(</span><span class="n">x</span><span class="p">),</span><span class="n">Maxpolydegree</span><span class="p">))</span>
<span class="n">X</span><span class="p">[:,</span><span class="mi">0</span><span class="p">]</span> <span class="o">=</span> <span class="mf">1.0</span>
<span class="k">for</span> <span class="n">polydegree</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="mi">1</span><span class="p">,</span> <span class="n">Maxpolydegree</span><span class="p">):</span>
<span class="k">for</span> <span class="n">degree</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">polydegree</span><span class="p">):</span>
<span class="n">X</span><span class="p">[:,</span><span class="n">degree</span><span class="p">]</span> <span class="o">=</span> <span class="n">x</span><span class="o">**</span><span class="p">(</span><span class="n">degree</span><span class="p">)</span>
<span class="c1"># We split the data in test and training data</span>
<span class="n">X_train</span><span class="p">,</span> <span class="n">X_test</span><span class="p">,</span> <span class="n">y_train</span><span class="p">,</span> <span class="n">y_test</span> <span class="o">=</span> <span class="n">train_test_split</span><span class="p">(</span><span class="n">X</span><span class="p">,</span> <span class="n">y</span><span class="p">,</span> <span class="n">test_size</span><span class="o">=</span><span class="mf">0.2</span><span class="p">)</span>
<span class="c1"># Decide which values of lambda to use</span>
<span class="n">nlambdas</span> <span class="o">=</span> <span class="mi">5</span>
<span class="n">lambdas</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">logspace</span><span class="p">(</span><span class="o">-</span><span class="mi">3</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="n">nlambdas</span><span class="p">)</span>
<span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">nlambdas</span><span class="p">):</span>
<span class="n">lmb</span> <span class="o">=</span> <span class="n">lambdas</span><span class="p">[</span><span class="n">i</span><span class="p">]</span>
<span class="c1"># Make the fit using Ridge only</span>
<span class="n">RegRidge</span> <span class="o">=</span> <span class="n">linear_model</span><span class="o">.</span><span class="n">Ridge</span><span class="p">(</span><span class="n">lmb</span><span class="p">,</span><span class="n">fit_intercept</span><span class="o">=</span><span class="kc">False</span><span class="p">)</span>
<span class="n">RegRidge</span><span class="o">.</span><span class="n">fit</span><span class="p">(</span><span class="n">X_train</span><span class="p">,</span><span class="n">y_train</span><span class="p">)</span>
<span class="c1"># and then make the prediction</span>
<span class="n">ypredictRidge</span> <span class="o">=</span> <span class="n">RegRidge</span><span class="o">.</span><span class="n">predict</span><span class="p">(</span><span class="n">X_test</span><span class="p">)</span>
<span class="n">Coeffs</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">(</span><span class="n">RegRidge</span><span class="o">.</span><span class="n">coef_</span><span class="p">)</span>
<span class="n">BetaValues</span> <span class="o">=</span> <span class="n">pd</span><span class="o">.</span><span class="n">DataFrame</span><span class="p">(</span><span class="n">Coeffs</span><span class="p">)</span>
<span class="n">BetaValues</span><span class="o">.</span><span class="n">columns</span> <span class="o">=</span> <span class="p">[</span><span class="s1">&#39;beta&#39;</span><span class="p">]</span>
<span class="n">display</span><span class="p">(</span><span class="n">BetaValues</span><span class="p">)</span>
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<div class="output text_html"><div>
<style scoped>
.dataframe tbody tr th:only-of-type {
vertical-align: middle;
}
.dataframe tbody tr th {
vertical-align: top;
}
.dataframe thead th {
text-align: right;
}
</style>
<table border="1" class="dataframe">
<thead>
<tr style="text-align: right;">
<th></th>
<th>beta</th>
</tr>
</thead>
<tbody>
<tr>
<th>0</th>
<td>0.986699</td>
</tr>
<tr>
<th>1</th>
<td>-0.606760</td>
</tr>
<tr>
<th>2</th>
<td>1.280573</td>
</tr>
<tr>
<th>3</th>
<td>-0.850164</td>
</tr>
<tr>
<th>4</th>
<td>0.000000</td>
</tr>
</tbody>
</table>
</div></div><div class="output text_html"><div>
<style scoped>
.dataframe tbody tr th:only-of-type {
vertical-align: middle;
}
.dataframe tbody tr th {
vertical-align: top;
}
.dataframe thead th {
text-align: right;
}
</style>
<table border="1" class="dataframe">
<thead>
<tr style="text-align: right;">
<th></th>
<th>beta</th>
</tr>
</thead>
<tbody>
<tr>
<th>0</th>
<td>0.978553</td>
</tr>
<tr>
<th>1</th>
<td>-0.511888</td>
</tr>
<tr>
<th>2</th>
<td>1.051418</td>
</tr>
<tr>
<th>3</th>
<td>-0.701370</td>
</tr>
<tr>
<th>4</th>
<td>0.000000</td>
</tr>
</tbody>
</table>
</div></div><div class="output text_html"><div>
<style scoped>
.dataframe tbody tr th:only-of-type {
vertical-align: middle;
}
.dataframe tbody tr th {
vertical-align: top;
}
.dataframe thead th {
text-align: right;
}
</style>
<table border="1" class="dataframe">
<thead>
<tr style="text-align: right;">
<th></th>
<th>beta</th>
</tr>
</thead>
<tbody>
<tr>
<th>0</th>
<td>0.946957</td>
</tr>
<tr>
<th>1</th>
<td>-0.162246</td>
</tr>
<tr>
<th>2</th>
<td>0.221921</td>
</tr>
<tr>
<th>3</th>
<td>-0.167787</td>
</tr>
<tr>
<th>4</th>
<td>0.000000</td>
</tr>
</tbody>
</table>
</div></div><div class="output text_html"><div>
<style scoped>
.dataframe tbody tr th:only-of-type {
vertical-align: middle;
}
.dataframe tbody tr th {
vertical-align: top;
}
.dataframe thead th {
text-align: right;
}
</style>
<table border="1" class="dataframe">
<thead>
<tr style="text-align: right;">
<th></th>
<th>beta</th>
</tr>
</thead>
<tbody>
<tr>
<th>0</th>
<td>0.906747</td>
</tr>
<tr>
<th>1</th>
<td>0.017665</td>
</tr>
<tr>
<th>2</th>
<td>-0.029483</td>
</tr>
<tr>
<th>3</th>
<td>-0.053849</td>
</tr>
<tr>
<th>4</th>
<td>0.000000</td>
</tr>
</tbody>
</table>
</div></div><div class="output text_html"><div>
<style scoped>
.dataframe tbody tr th:only-of-type {
vertical-align: middle;
}
.dataframe tbody tr th {
vertical-align: top;
}
.dataframe thead th {
text-align: right;
}
</style>
<table border="1" class="dataframe">
<thead>
<tr style="text-align: right;">
<th></th>
<th>beta</th>
</tr>
</thead>
<tbody>
<tr>
<th>0</th>
<td>0.718165</td>
</tr>
<tr>
<th>1</th>
<td>0.156956</td>
</tr>
<tr>
<th>2</th>
<td>0.040102</td>
</tr>
<tr>
<th>3</th>
<td>-0.001880</td>
</tr>
<tr>
<th>4</th>
<td>0.000000</td>
</tr>
</tbody>
</table>
</div></div></div>
</div>
<p>As an exercise, repeat these calculations with ordinary least squares
only with and without noise. Calculate thereafter the variance of the
parameters <span class="math notranslate nohighlight">\(\beta_j\)</span> as function of polynomial order and of the added
noise. Here we recommend to use <span class="math notranslate nohighlight">\(\sigma^2=1\)</span> as variance for the
added noise (which follows a normal distribution with mean value zero).
Comment your results. If you have a large noise term, do the parameters <span class="math notranslate nohighlight">\(\beta_j\)</span> vary more as function
model complexity? And what about their variance?</p>
</div>
<div class="section" id="linking-bayes-theorem-with-ridge-and-lasso-regression">
<h2><span class="section-number">4.14. </span>Linking Bayes Theorem with Ridge and Lasso Regression<a class="headerlink" href="#linking-bayes-theorem-with-ridge-and-lasso-regression" title="Permalink to this headline"></a></h2>
<p>Using Bayes theorem we can gain a better intuition about Ridge and Lasso regression.</p>
<p>For ordinary least squares we postulated that the maximum likelihood for the doamin of events <span class="math notranslate nohighlight">\(\boldsymbol{D}\)</span> (one-dimensional case)</p>
<p>We have seen that Ridge regression suppresses those features which
have a small singular value. This corresponds to a feature which exhibits
a large variance in the parameters <span class="math notranslate nohighlight">\(\beta_j\)</span>.
Our analysis hitherto has been based on linear algebra. To add to our intuition, we will use
Bayes theorem in order to deepen our understanding of Ridge and Lasso regression.</p>
<p>For ordinary least squares we postulated that the maximum likelihood for the domain of events <span class="math notranslate nohighlight">\(\boldsymbol{D}\)</span> (one-dimensional case)</p>
<div class="math notranslate nohighlight">
\[
\boldsymbol{D}=[(x_0,y_0), (x_1,y_1),\dots, (x_{n-1},y_{n-1})],
@@ -2893,7 +3191,10 @@ p(\boldsymbol{\beta}\vert\boldsymbol{D})\propto p(\boldsymbol{D}\vert\boldsymbol
<p>With the posterior probability defined by a likelihood which we have
already modeled and an unknown prior, we are now ready to make
additional models for the prior.</p>
<p>We can, based on our discussions of the variance of <span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span> and the mean value, assume that the prior for the values <span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span> is given by a Gaussian with mean value zero and variance <span class="math notranslate nohighlight">\(\tau^2\)</span>, that is</p>
<p>We can, based on our discussions of the variance of <span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span> and
the mean value, assume that the prior for the values <span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span> is
given by a Gaussian with mean value zero and variance <span class="math notranslate nohighlight">\(\tau^2\)</span>, that
is</p>
<div class="math notranslate nohighlight">
\[
p(\boldsymbol{\beta})=\prod_{j=0}^{p-1}\exp{\left(-\frac{\beta_j^2}{2\tau^2}\right)}.
@@ -2906,7 +3207,7 @@ p(\boldsymbol{\beta\vert\boldsymbol{D})}=\prod_{i=0}^{n-1}\frac{1}{\sqrt{2\pi\si
<p>We can now optimize this quantity with respect to <span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span>. As we
did for OLS, this is most conveniently done by taking the negative
logarithm of the posterior probability. Doing so and leaving out the
constants terms that do not depend on <span class="math notranslate nohighlight">\(\beta\)</span>, we have</p>
terms that do not depend on <span class="math notranslate nohighlight">\(\beta\)</span>, we have</p>
<div class="math notranslate nohighlight">
\[
C(\boldsymbol{\beta})=\frac{\vert\vert (\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta})\vert\vert_2^2}{2\sigma^2}+\frac{1}{2\tau^2}\vert\vert\boldsymbol{\beta}\vert\vert_2^2,
@@ -2940,6 +3241,13 @@ C(\boldsymbol{\beta}=\frac{\vert\vert (\boldsymbol{y}-\boldsymbol{X}\boldsymbol{
C(\boldsymbol{\beta}=\frac{\vert\vert (\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta})\vert\vert_2^2}{2\sigma^2}+\lambda\vert\vert\boldsymbol{\beta}\vert\vert_1,
\]</div>
<p>which is our Lasso cost function!</p>
<p>Plotting these prior functions shows us that we can use the parameter
<span class="math notranslate nohighlight">\(\lambda\)</span> to shrink or increase the role of a given parameter
<span class="math notranslate nohighlight">\(\beta_j\)</span>. The variance for the Laplace distribution is
<span class="math notranslate nohighlight">\(2\tau^2=1/\lambda\)</span> while for the Gaussian distribution it is
<span class="math notranslate nohighlight">\(\sigma^2=1/(2\lambda)\)</span>. Thus, increasing the variance means
decreasing <span class="math notranslate nohighlight">\(\lambda\)</span> and shrinking the variance means increasing
<span class="math notranslate nohighlight">\(\lamdbda\)</span>. When we increase <span class="math notranslate nohighlight">\(\lambda\)</span>, this corresponds to shrinking the role of less important features (small singular values).</p>
</div>
</div>
+157 -50
View File
@@ -495,8 +495,9 @@ mean value <span class="math notranslate nohighlight">\(\boldsymbol{X}\boldsymbo
\[
\mathbb{E}(\boldsymbol{\beta}) = \mathbb{E}[ (\mathbf{X}^{\top} \mathbf{X})^{-1}\mathbf{X}^{T} \mathbf{Y}]=(\mathbf{X}^{T} \mathbf{X})^{-1}\mathbf{X}^{T} \mathbb{E}[ \mathbf{Y}]=(\mathbf{X}^{T} \mathbf{X})^{-1} \mathbf{X}^{T}\mathbf{X}\boldsymbol{\beta}=\boldsymbol{\beta}.
\]</div>
<p>This means that the estimator of the regression parameters is unbiased.</p>
<p>We can also calculate the variance</p>
<p>This means that the estimator of the regression parameters is unbiased.
v
We can also calculate the variance</p>
<p>The variance of <span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span> is</p>
<div class="math notranslate nohighlight">
\[\begin{split}
@@ -629,10 +630,10 @@ number <span class="math notranslate nohighlight">\(i\)</span> is left out. Usin
</div>
</div>
<div class="cell_output docutils container">
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>Runtime: 0.137546 sec
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>Runtime: 0.141725 sec
Jackknife Statistics :
original bias std. error
99.9031 99.8931 0.149233
100.148 100.138 0.1492
</pre></div>
</div>
</div>
@@ -698,7 +699,111 @@ histogram of the relative frequency of <span class="math notranslate nohighlight
you use the estimators corresponding to the statistic of interest. For
example, if you are interested in estimating the variance of <span class="math notranslate nohighlight">\(\widehat
\beta\)</span>, apply the etsimator <span class="math notranslate nohighlight">\(\widehat \sigma^2\)</span> to the values
<span class="math notranslate nohighlight">\(\widehat \beta ^*\)</span>.</p>
<span class="math notranslate nohighlight">\(\widehat \beta^*\)</span>.</p>
<p>Before we proceed however, we need to remind ourselves about a central
theorem in statistics, namely the so-called <strong>central limit theorem</strong>.
This theorem plays a central role in understanding why the Bootstrap
(and other resampling methods) work so well on independent and
identically distributed variables.</p>
<p>Suppose we have a PDF <span class="math notranslate nohighlight">\(p(x)\)</span> from which we generate a series <span class="math notranslate nohighlight">\(N\)</span>
of averages <span class="math notranslate nohighlight">\(\langle x_i \rangle\)</span>. Each mean value <span class="math notranslate nohighlight">\(\langle x_i \rangle\)</span>
is viewed as the average of a specific measurement, e.g., throwing
dice 100 times and then taking the average value, or producing a certain
amount of random numbers.
For notational ease, we set <span class="math notranslate nohighlight">\(\langle x_i \rangle=x_i\)</span> in the discussion
which follows.</p>
<p>If we compute the mean <span class="math notranslate nohighlight">\(z\)</span> of <span class="math notranslate nohighlight">\(m\)</span> such mean values <span class="math notranslate nohighlight">\(x_i\)</span></p>
<div class="math notranslate nohighlight">
\[
z=\frac{x_1+x_2+\dots+x_m}{m},
\]</div>
<p>the question we pose is which is the PDF of the new variable <span class="math notranslate nohighlight">\(z\)</span>.</p>
<p>The probability of obtaining an average value <span class="math notranslate nohighlight">\(z\)</span> is the product of the
probabilities of obtaining arbitrary individual mean values <span class="math notranslate nohighlight">\(x_i\)</span>,
but with the constraint that the average is <span class="math notranslate nohighlight">\(z\)</span>. We can express this through
the following expression</p>
<div class="math notranslate nohighlight">
\[
\tilde{p}(z)=\int dx_1p(x_1)\int dx_2p(x_2)\dots\int dx_mp(x_m)
\delta(z-\frac{x_1+x_2+\dots+x_m}{m}),
\]</div>
<p>where the <span class="math notranslate nohighlight">\(\delta\)</span>-function enbodies the constraint that the mean is <span class="math notranslate nohighlight">\(z\)</span>.
All measurements that lead to each individual <span class="math notranslate nohighlight">\(x_i\)</span> are expected to
be independent, which in turn means that we can express <span class="math notranslate nohighlight">\(\tilde{p}\)</span> as the
product of individual <span class="math notranslate nohighlight">\(p(x_i)\)</span>. The independence assumption is important in the derivation of the central limit theorem.</p>
<p>If we use the integral expression for the <span class="math notranslate nohighlight">\(\delta\)</span>-function</p>
<div class="math notranslate nohighlight">
\[
\delta(z-\frac{x_1+x_2+\dots+x_m}{m})=\frac{1}{2\pi}\int_{-\infty}^{\infty}
dq\exp{\left(iq(z-\frac{x_1+x_2+\dots+x_m}{m})\right)},
\]</div>
<p>and inserting <span class="math notranslate nohighlight">\(e^{i\mu q-i\mu q}\)</span> where <span class="math notranslate nohighlight">\(\mu\)</span> is the mean value
we arrive at</p>
<div class="math notranslate nohighlight">
\[
\tilde{p}(z)=\frac{1}{2\pi}\int_{-\infty}^{\infty}
dq\exp{\left(iq(z-\mu)\right)}\left[\int_{-\infty}^{\infty}
dxp(x)\exp{\left(iq(\mu-x)/m\right)}\right]^m,
\]</div>
<p>with the integral over <span class="math notranslate nohighlight">\(x\)</span> resulting in</p>
<div class="math notranslate nohighlight">
\[
\int_{-\infty}^{\infty}dxp(x)\exp{\left(iq(\mu-x)/m\right)}=
\int_{-\infty}^{\infty}dxp(x)
\left[1+\frac{iq(\mu-x)}{m}-\frac{q^2(\mu-x)^2}{2m^2}+\dots\right].
\]</div>
<p>The second term on the rhs disappears since this is just the mean and
employing the definition of <span class="math notranslate nohighlight">\(\sigma^2\)</span> we have</p>
<div class="math notranslate nohighlight">
\[
\int_{-\infty}^{\infty}dxp(x)e^{\left(iq(\mu-x)/m\right)}=
1-\frac{q^2\sigma^2}{2m^2}+\dots,
\]</div>
<p>resulting in</p>
<div class="math notranslate nohighlight">
\[
\left[\int_{-\infty}^{\infty}dxp(x)\exp{\left(iq(\mu-x)/m\right)}\right]^m\approx
\left[1-\frac{q^2\sigma^2}{2m^2}+\dots \right]^m,
\]</div>
<p>and in the limit <span class="math notranslate nohighlight">\(m\rightarrow \infty\)</span> we obtain</p>
<div class="math notranslate nohighlight">
\[
\tilde{p}(z)=\frac{1}{\sqrt{2\pi}(\sigma/\sqrt{m})}
\exp{\left(-\frac{(z-\mu)^2}{2(\sigma/\sqrt{m})^2}\right)},
\]</div>
<p>which is the normal distribution with variance
<span class="math notranslate nohighlight">\(\sigma^2_m=\sigma^2/m\)</span>, where <span class="math notranslate nohighlight">\(\sigma\)</span> is the variance of the PDF <span class="math notranslate nohighlight">\(p(x)\)</span>
and <span class="math notranslate nohighlight">\(\mu\)</span> is also the mean of the PDF <span class="math notranslate nohighlight">\(p(x)\)</span>.</p>
<p>Thus, the central limit theorem states that the PDF <span class="math notranslate nohighlight">\(\tilde{p}(z)\)</span> of
the average of <span class="math notranslate nohighlight">\(m\)</span> random values corresponding to a PDF <span class="math notranslate nohighlight">\(p(x)\)</span>
is a normal distribution whose mean is the
mean value of the PDF <span class="math notranslate nohighlight">\(p(x)\)</span> and whose variance is the variance
of the PDF <span class="math notranslate nohighlight">\(p(x)\)</span> divided by <span class="math notranslate nohighlight">\(m\)</span>, the number of values used to compute <span class="math notranslate nohighlight">\(z\)</span>.</p>
<p>The central limit theorem leads to the well-known expression for the
standard deviation, given by</p>
<div class="math notranslate nohighlight">
\[
\sigma_m=
\frac{\sigma}{\sqrt{m}}.
\]</div>
<p>The latter is true only if the average value is known exactly. This is obtained in the limit
<span class="math notranslate nohighlight">\(m\rightarrow \infty\)</span> only. Because the mean and the variance are measured quantities we obtain
the familiar expression in statistics</p>
<div class="math notranslate nohighlight">
\[
\sigma_m\approx
\frac{\sigma}{\sqrt{m-1}}.
\]</div>
<p>In many cases however the above estimate for the standard deviation,
in particular if correlations are strong, may be too simplistic. Keep
in mind that we have assumed that the variables <span class="math notranslate nohighlight">\(x\)</span> are independent
and identically distributed. This is obviously not always the
case. For example, the random numbers (or better pseudorandom numbers)
we generate in various calculations do always exhibit some
correlations.</p>
<p>The theorem is satisfied by a large class of PDFs. Note however that for a
finite <span class="math notranslate nohighlight">\(m\)</span>, it is not always possible to find a closed form /analytic expression for
<span class="math notranslate nohighlight">\(\tilde{p}(x)\)</span>.</p>
</div>
</div>
<div class="section" id="the-bias-variance-tradeoff">
@@ -818,7 +923,7 @@ Var: 0.0010479246398391328
0.01312157412031145 &gt;= 0.012073649480472317 + 0.0010479246398391328 = 0.01312157412031145
</pre></div>
</div>
<img alt="_images/chapter3_37_1.png" src="_images/chapter3_37_1.png" />
<img alt="_images/chapter3_57_1.png" src="_images/chapter3_57_1.png" />
</div>
</div>
<div class="cell docutils container">
@@ -910,10 +1015,9 @@ Error: 0.03781367141738898
Bias^2: 0.03365768507152761
Var: 0.004155986345861379
0.03781367141738898 &gt;= 0.03365768507152761 + 0.004155986345861379 = 0.03781367141738899
Polynomial degree:
</pre></div>
</div>
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span> 7
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>Polynomial degree: 7
Error: 0.027609773491022498
Bias^2: 0.02299949826036597
Var: 0.004610275230656537
@@ -923,15 +1027,14 @@ Error: 0.017355848195591973
Bias^2: 0.010331721306655588
Var: 0.007024126888936384
0.017355848195591973 &gt;= 0.010331721306655588 + 0.007024126888936384 = 0.017355848195591973
Polynomial degree: 9
</pre></div>
</div>
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>Polynomial degree: 9
Error: 0.026605727637189085
Bias^2: 0.010018312644140933
Var: 0.016587414993048166
0.026605727637189085 &gt;= 0.010018312644140933 + 0.016587414993048166 = 0.0266057276371891
Polynomial degree:
</pre></div>
</div>
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span> 10
Polynomial degree: 10
Error: 0.021592704588043153
Bias^2: 0.010516485576652981
Var: 0.011076219011390184
@@ -946,14 +1049,16 @@ Error: 0.1154777721897675
Bias^2: 0.01628578269590588
Var: 0.09919198949386163
0.1154777721897675 &gt;= 0.01628578269590588 + 0.09919198949386163 = 0.11547777218976751
Polynomial degree: 13
</pre></div>
</div>
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>Polynomial degree: 13
Error: 0.22842468702166951
Bias^2: 0.01975416527163567
Var: 0.20867052175003387
0.22842468702166951 &gt;= 0.01975416527163567 + 0.20867052175003387 = 0.22842468702166954
</pre></div>
</div>
<img alt="_images/chapter3_38_4.png" src="_images/chapter3_38_4.png" />
<img alt="_images/chapter3_58_5.png" src="_images/chapter3_58_5.png" />
</div>
</div>
<p>The bias-variance tradeoff summarizes the fundamental tension in
@@ -1078,7 +1183,7 @@ set, the higher, the less likely the model generalizes correctly from the
training data.
</pre></div>
</div>
<img alt="_images/chapter3_40_1.png" src="_images/chapter3_40_1.png" />
<img alt="_images/chapter3_60_1.png" src="_images/chapter3_60_1.png" />
</div>
</div>
<div class="cell docutils container">
@@ -1185,84 +1290,86 @@ Mean squared error on test data: 5.98822371
Degree of polynomial: 6
Mean squared error on training data: 3.66204648
Mean squared error on test data: 8.14812206
Degree of polynomial: 7
Mean squared error on training data: 0.47075725
Mean squared error on test data: 2.00607783
</pre></div>
</div>
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>Degree of polynomial: 8
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>Degree of polynomial: 7
Mean squared error on training data: 0.47075725
Mean squared error on test data: 2.00607783
Degree of polynomial: 8
Mean squared error on training data: 0.04912436
Mean squared error on test data: 0.21596432
Degree of polynomial: 9
Mean squared error on training data: 0.02522069
Mean squared error on test data: 0.08576932
Degree of polynomial: 10
Mean squared error on training data: 0.02511518
Mean squared error on test data: 1.20015436
</pre></div>
</div>
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>Degree of polynomial: 11
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>Degree of polynomial: 10
Mean squared error on training data: 0.02511518
Mean squared error on test data: 1.20015436
Degree of polynomial: 11
Mean squared error on training data: 0.01640891
Mean squared error on test data: 1.35533773
Degree of polynomial: 12
Mean squared error on training data: 0.00813803
Mean squared error on test data: 0.17446471
Degree of polynomial: 13
</pre></div>
</div>
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>Degree of polynomial: 13
Mean squared error on training data: 0.00759119
Mean squared error on test data: 1.08131003
Degree of polynomial: 14
Mean squared error on training data: 0.00472199
Mean squared error on test data: 0.81333793
</pre></div>
</div>
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>Degree of polynomial: 15
Degree of polynomial: 15
Mean squared error on training data: 0.00410478
Mean squared error on test data: 92.09145189
Degree of polynomial: 16
</pre></div>
</div>
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>Degree of polynomial: 16
Mean squared error on training data: 0.00315593
Mean squared error on test data: 234.39716546
Degree of polynomial: 17
Mean squared error on training data: 0.00242998
Mean squared error on test data: 1271.05295709
</pre></div>
</div>
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>Degree of polynomial: 18
Degree of polynomial: 18
Mean squared error on training data: 0.00228740
Mean squared error on test data: 108.42208194
Degree of polynomial: 19
</pre></div>
</div>
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>Degree of polynomial: 19
Mean squared error on training data: 0.00156372
Mean squared error on test data: 1388.41078073
Degree of polynomial: 20
Mean squared error on training data: 0.00137982
Mean squared error on test data: 1761.43341615
</pre></div>
</div>
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>Degree of polynomial: 21
Degree of polynomial: 21
Mean squared error on training data: 0.00118170
Mean squared error on test data: 15061.31603087
Degree of polynomial: 22
</pre></div>
</div>
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>Degree of polynomial: 22
Mean squared error on training data: 0.00092354
Mean squared error on test data: 890.63488525
Degree of polynomial: 23
Mean squared error on training data: 0.00085887
Mean squared error on test data: 5483.16796929
</pre></div>
</div>
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>Degree of polynomial: 24
Degree of polynomial: 24
Mean squared error on training data: 0.00084589
Mean squared error on test data: 1695.57143061
Degree of polynomial: 25
</pre></div>
</div>
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>Degree of polynomial: 25
Mean squared error on training data: 0.00078806
Mean squared error on test data: 131343.30655001
Degree of polynomial: 26
Mean squared error on training data: 0.00076916
Mean squared error on test data: 17709.14370264
</pre></div>
</div>
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>Degree of polynomial: 27
Degree of polynomial: 27
Mean squared error on training data: 0.00068970
Mean squared error on test data: 2975.38903780
Degree of polynomial: 28
</pre></div>
</div>
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>Degree of polynomial: 28
Mean squared error on training data: 0.00062588
Mean squared error on test data: 3848.64522721
Degree of polynomial: 29
@@ -1276,7 +1383,7 @@ Mean squared error on test data: 2988.64001211
plt.plot(polynomial, np.log10(testerror), label=&#39;Test Error&#39;)
</pre></div>
</div>
<img alt="_images/chapter3_41_10.png" src="_images/chapter3_41_10.png" />
<img alt="_images/chapter3_61_11.png" src="_images/chapter3_61_11.png" />
</div>
</div>
</div>
@@ -1429,7 +1536,7 @@ cross-validation (LOOCV).</p>
</div>
</div>
<div class="cell_output docutils container">
<img alt="_images/chapter3_47_0.png" src="_images/chapter3_47_0.png" />
<img alt="_images/chapter3_67_0.png" src="_images/chapter3_67_0.png" />
</div>
</div>
<p>More examples of the application of cross-validation follow here.</p>
@@ -1510,7 +1617,7 @@ cross-validation (LOOCV).</p>
plt.plot(polynomial, np.log10(estimated_mse_sklearn), label=&#39;Test Error&#39;)
</pre></div>
</div>
<img alt="_images/chapter3_49_1.png" src="_images/chapter3_49_1.png" />
<img alt="_images/chapter3_69_1.png" src="_images/chapter3_69_1.png" />
</div>
</div>
<div class="cell docutils container">
@@ -1554,7 +1661,7 @@ cross-validation (LOOCV).</p>
</div>
</div>
<div class="cell_output docutils container">
<img alt="_images/chapter3_50_0.png" src="_images/chapter3_50_0.png" />
<img alt="_images/chapter3_70_0.png" src="_images/chapter3_70_0.png" />
</div>
</div>
</div>
@@ -1643,7 +1750,7 @@ which polynomial fits the data best.</p>
</div>
</div>
<div class="cell_output docutils container">
<img alt="_images/chapter3_54_0.png" src="_images/chapter3_54_0.png" />
<img alt="_images/chapter3_74_0.png" src="_images/chapter3_74_0.png" />
</div>
</div>
<div class="section" id="exercise-ordinary-least-square-ols-on-the-franke-function">
File diff suppressed because one or more lines are too long
@@ -1470,10 +1470,10 @@
"name": "stdout",
"output_type": "stream",
"text": [
"-0.1831277634928002\n",
"3.3041320306136366\n",
"[[1.03929932 3.07670437]\n",
" [3.07670437 9.7701384 ]]\n"
"-0.055529303095955385\n",
"3.9154537458386387\n",
"[[1.0003451 2.94327895]\n",
" [2.94327895 9.80021057]]\n"
]
}
],
@@ -1513,10 +1513,10 @@
"name": "stdout",
"output_type": "stream",
"text": [
"0.09456011349477329\n",
"1.6600855222624895\n",
"[[1. 0.70183798]\n",
" [0.70183798 1. ]]\n"
"0.095786152583691\n",
"1.5864689451163851\n",
"[[1. 0.6713619]\n",
" [0.6713619 1. ]]\n"
]
}
],
@@ -1571,30 +1571,30 @@
"name": "stdout",
"output_type": "stream",
"text": [
"[[ 0.54316925 0.40968888]\n",
" [-1.8393131 -5.8788431 ]\n",
" [ 0.24613822 0.9957744 ]\n",
" [ 1.38892672 3.93679587]\n",
" [ 0.06540809 -0.84008474]\n",
" [ 0.49841285 2.31081134]\n",
" [-0.49323032 -1.25263916]\n",
" [ 0.08380813 -0.7925146 ]\n",
" [ 0.64078247 2.92814088]\n",
" [-1.13410232 -1.81712976]]\n",
"[[-0.10979444 1.25912505]\n",
" [ 2.17273709 7.48017006]\n",
" [ 0.31520842 0.66152576]\n",
" [-0.07076926 -0.44655382]\n",
" [ 0.24348403 0.38561052]\n",
" [-0.57984245 -0.86011441]\n",
" [ 0.42375621 -1.21244261]\n",
" [-0.2407922 0.56831157]\n",
" [-2.23691315 -6.6936767 ]\n",
" [ 0.08292574 -1.14195542]]\n",
" 0 1\n",
"0 0.543169 0.409689\n",
"1 -1.839313 -5.878843\n",
"2 0.246138 0.995774\n",
"3 1.388927 3.936796\n",
"4 0.065408 -0.840085\n",
"5 0.498413 2.310811\n",
"6 -0.493230 -1.252639\n",
"7 0.083808 -0.792515\n",
"8 0.640782 2.928141\n",
"9 -1.134102 -1.817130\n",
"0 -0.109794 1.259125\n",
"1 2.172737 7.480170\n",
"2 0.315208 0.661526\n",
"3 -0.070769 -0.446554\n",
"4 0.243484 0.385611\n",
"5 -0.579842 -0.860114\n",
"6 0.423756 -1.212443\n",
"7 -0.240792 0.568312\n",
"8 -2.236913 -6.693677\n",
"9 0.082926 -1.141955\n",
" 0 1\n",
"0 1.000000 0.942726\n",
"1 0.942726 1.000000\n"
"0 1.000000 0.931066\n",
"1 0.931066 1.000000\n"
]
}
],
@@ -1636,37 +1636,37 @@
"text": [
" 0 1 2 3 4 5 6 7 \\\n",
"0 0.0 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 \n",
"1 0.0 0.080633 0.085811 0.085748 0.084826 0.083221 0.080341 0.078342 \n",
"2 0.0 0.085811 0.092659 0.093248 0.092932 0.091659 0.088414 0.086561 \n",
"3 0.0 0.085748 0.093248 0.096639 0.096802 0.095941 0.093928 0.092288 \n",
"4 0.0 0.084826 0.092932 0.096802 0.097346 0.096763 0.094774 0.093328 \n",
"5 0.0 0.083221 0.091659 0.095941 0.096763 0.096406 0.094513 0.093236 \n",
"6 0.0 0.080341 0.088414 0.093928 0.094774 0.094513 0.093573 0.092375 \n",
"7 0.0 0.078342 0.086561 0.092288 0.093328 0.093236 0.092375 0.091318 \n",
"8 0.0 0.076213 0.084461 0.090355 0.091535 0.091578 0.090817 0.089881 \n",
"9 0.0 0.074037 0.082239 0.088258 0.089540 0.089693 0.089039 0.088209 \n",
"10 0.0 0.073546 0.081474 0.088137 0.089326 0.089437 0.089352 0.088484 \n",
"11 0.0 0.071453 0.079337 0.086052 0.087334 0.087544 0.087526 0.086757 \n",
"12 0.0 0.069371 0.077163 0.083905 0.085253 0.085544 0.085598 0.084916 \n",
"13 0.0 0.067326 0.074996 0.081744 0.083138 0.083495 0.083622 0.083015 \n",
"14 0.0 0.065334 0.072862 0.079597 0.081024 0.081435 0.081632 0.081092 \n",
"1 0.0 0.090107 0.080265 0.095180 0.087148 0.079801 0.089928 0.082666 \n",
"2 0.0 0.080265 0.073052 0.087032 0.080565 0.074557 0.083797 0.077645 \n",
"3 0.0 0.095180 0.087032 0.106946 0.099079 0.091756 0.104887 0.097150 \n",
"4 0.0 0.087148 0.080565 0.099079 0.092378 0.086076 0.098028 0.091253 \n",
"5 0.0 0.079801 0.074557 0.091756 0.086076 0.080678 0.091542 0.085627 \n",
"6 0.0 0.089928 0.083797 0.104887 0.098028 0.091542 0.105387 0.098187 \n",
"7 0.0 0.082666 0.077645 0.097150 0.091253 0.085627 0.098187 0.091855 \n",
"8 0.0 0.076121 0.072046 0.090099 0.085041 0.080173 0.091566 0.086006 \n",
"9 0.0 0.070212 0.066949 0.083672 0.079351 0.075151 0.085484 0.080612 \n",
"10 0.0 0.083337 0.078717 0.099609 0.093725 0.088085 0.101781 0.095275 \n",
"11 0.0 0.076774 0.072995 0.092276 0.087205 0.082303 0.094708 0.088983 \n",
"12 0.0 0.070866 0.067811 0.085627 0.081270 0.077021 0.088260 0.083228 \n",
"13 0.0 0.065545 0.063114 0.079598 0.075869 0.072197 0.082382 0.077968 \n",
"14 0.0 0.060748 0.058856 0.074130 0.070955 0.067793 0.077026 0.073162 \n",
"\n",
" 8 9 10 11 12 13 14 \n",
"0 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 \n",
"1 0.076213 0.074037 0.073546 0.071453 0.069371 0.067326 0.065334 \n",
"2 0.084461 0.082239 0.081474 0.079337 0.077163 0.074996 0.072862 \n",
"3 0.090355 0.088258 0.088137 0.086052 0.083905 0.081744 0.079597 \n",
"4 0.091535 0.089540 0.089326 0.087334 0.085253 0.083138 0.081024 \n",
"5 0.091578 0.089693 0.089437 0.087544 0.085544 0.083495 0.081435 \n",
"6 0.090817 0.089039 0.089352 0.087526 0.085598 0.083622 0.081632 \n",
"7 0.089881 0.088209 0.088484 0.086757 0.084916 0.083015 0.081092 \n",
"8 0.088556 0.086985 0.087242 0.085610 0.083855 0.082032 0.080182 \n",
"9 0.086985 0.085510 0.085763 0.084221 0.082549 0.080806 0.079029 \n",
"10 0.087242 0.085763 0.086430 0.084864 0.083177 0.081425 0.079643 \n",
"11 0.085610 0.084221 0.084864 0.083384 0.081778 0.080101 0.078390 \n",
"12 0.083855 0.082549 0.083177 0.081778 0.080249 0.078646 0.077005 \n",
"13 0.082032 0.080806 0.081425 0.080101 0.078646 0.077115 0.075542 \n",
"14 0.080182 0.079029 0.079643 0.078390 0.077005 0.075542 0.074036 \n"
"1 0.076121 0.070212 0.083337 0.076774 0.070866 0.065545 0.060748 \n",
"2 0.072046 0.066949 0.078717 0.072995 0.067811 0.063114 0.058856 \n",
"3 0.090099 0.083672 0.099609 0.092276 0.085627 0.079598 0.074130 \n",
"4 0.085041 0.079351 0.093725 0.087205 0.081270 0.075869 0.070955 \n",
"5 0.080173 0.075151 0.088085 0.082303 0.077021 0.072197 0.067793 \n",
"6 0.091566 0.085484 0.101781 0.094708 0.088260 0.082382 0.077026 \n",
"7 0.086006 0.080612 0.095275 0.088983 0.083228 0.077968 0.073162 \n",
"8 0.080846 0.076068 0.089250 0.083659 0.078529 0.073827 0.069519 \n",
"9 0.076068 0.071841 0.083682 0.078720 0.074154 0.069956 0.066098 \n",
"10 0.089250 0.083682 0.099504 0.092936 0.086920 0.081414 0.076377 \n",
"11 0.083659 0.078720 0.092936 0.087096 0.081732 0.076810 0.072296 \n",
"12 0.078529 0.074154 0.086920 0.081732 0.076954 0.072557 0.068515 \n",
"13 0.073827 0.069956 0.081414 0.076810 0.072557 0.068633 0.065016 \n",
"14 0.069519 0.066098 0.076377 0.072296 0.068515 0.065016 0.061783 \n"
]
}
],
@@ -4546,16 +4546,407 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"How can we understand this? **More text to be added**.\n",
"How can we understand this?\n",
"\n",
"Let us write out the values of the coefficients $\\beta_i$ as functions\n",
"of the polynomial degree and noise. We will focus only on the Ridge\n",
"results and some few selected values of the hyperparameter $\\lambda$.\n",
"\n",
"If we don't include any noise and run this code for different values\n",
"of the polynomial degree, we notice that the results for $\\beta_i$ do\n",
"not show great changes from one order to the next. This is an\n",
"indication that for higher polynomial orders, our parameters become\n",
"less important.\n",
"\n",
"If we however add noise, what happens is that the polynomial fit is\n",
"trying to adjust the fit to traverse in the best possible way all data\n",
"points. This can lead to large fluctuations in the parameters\n",
"$\\beta_i$ as functions of polynomial order. It will also be reflected\n",
"in a larger value of the variance of each parameter $\\beta_i$. What\n",
"Ridge regression (and Lasso as well) are doing then is to try to\n",
"quench the fluctuations in the parameters of $\\beta_i$ which have a\n",
"large variance (normally for higher orders in the polynomial)."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"collapsed": false,
"editable": true
},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>beta</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>0.986699</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>-0.606760</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>1.280573</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>-0.850164</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>0.000000</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" beta\n",
"0 0.986699\n",
"1 -0.606760\n",
"2 1.280573\n",
"3 -0.850164\n",
"4 0.000000"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>beta</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>0.978553</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>-0.511888</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>1.051418</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>-0.701370</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>0.000000</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" beta\n",
"0 0.978553\n",
"1 -0.511888\n",
"2 1.051418\n",
"3 -0.701370\n",
"4 0.000000"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>beta</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>0.946957</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>-0.162246</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>0.221921</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>-0.167787</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>0.000000</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" beta\n",
"0 0.946957\n",
"1 -0.162246\n",
"2 0.221921\n",
"3 -0.167787\n",
"4 0.000000"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>beta</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>0.906747</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>0.017665</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>-0.029483</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>-0.053849</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>0.000000</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" beta\n",
"0 0.906747\n",
"1 0.017665\n",
"2 -0.029483\n",
"3 -0.053849\n",
"4 0.000000"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>beta</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>0.718165</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>0.156956</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>0.040102</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>-0.001880</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>0.000000</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" beta\n",
"0 0.718165\n",
"1 0.156956\n",
"2 0.040102\n",
"3 -0.001880\n",
"4 0.000000"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"import numpy as np\n",
"import pandas as pd\n",
"from IPython.display import display\n",
"import matplotlib.pyplot as plt\n",
"from sklearn.model_selection import train_test_split\n",
"from sklearn import linear_model\n",
"\n",
"# Make data set.\n",
"n = 1000\n",
"x = np.random.rand(n)\n",
"y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)+ np.random.randn(n)\n",
"\n",
"Maxpolydegree = 5\n",
"X = np.zeros((len(x),Maxpolydegree))\n",
"X[:,0] = 1.0\n",
"\n",
"for polydegree in range(1, Maxpolydegree):\n",
" for degree in range(polydegree):\n",
" X[:,degree] = x**(degree)\n",
"\n",
"\n",
"# We split the data in test and training data\n",
"X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)\n",
"\n",
"# Decide which values of lambda to use\n",
"nlambdas = 5\n",
"lambdas = np.logspace(-3, 2, nlambdas)\n",
"for i in range(nlambdas):\n",
" lmb = lambdas[i]\n",
" # Make the fit using Ridge only\n",
" RegRidge = linear_model.Ridge(lmb,fit_intercept=False)\n",
" RegRidge.fit(X_train,y_train)\n",
" # and then make the prediction\n",
" ypredictRidge = RegRidge.predict(X_test)\n",
" Coeffs = np.array(RegRidge.coef_)\n",
" BetaValues = pd.DataFrame(Coeffs)\n",
" BetaValues.columns = ['beta']\n",
" display(BetaValues)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As an exercise, repeat these calculations with ordinary least squares\n",
"only with and without noise. Calculate thereafter the variance of the\n",
"parameters $\\beta_j$ as function of polynomial order and of the added\n",
"noise. Here we recommend to use $\\sigma^2=1$ as variance for the\n",
"added noise (which follows a normal distribution with mean value zero).\n",
"Comment your results. If you have a large noise term, do the parameters $\\beta_j$ vary more as function\n",
"model complexity? And what about their variance? \n",
"\n",
"\n",
"\n",
"\n",
"## Linking Bayes' Theorem with Ridge and Lasso Regression\n",
"\n",
"We have seen that Ridge regression suppresses those features which\n",
"have a small singular value. This corresponds to a feature which exhibits\n",
"a large variance in the parameters $\\beta_j$.\n",
"Our analysis hitherto has been based on linear algebra. To add to our intuition, we will use\n",
"Bayes' theorem in order to deepen our understanding of Ridge and Lasso regression. \n",
"\n",
"Using Bayes' theorem we can gain a better intuition about Ridge and Lasso regression. \n",
"\n",
"For ordinary least squares we postulated that the maximum likelihood for the doamin of events $\\boldsymbol{D}$ (one-dimensional case)"
"For ordinary least squares we postulated that the maximum likelihood for the domain of events $\\boldsymbol{D}$ (one-dimensional case)"
]
},
{
@@ -4627,7 +5018,10 @@
"already modeled and an unknown prior, we are now ready to make\n",
"additional models for the prior.\n",
"\n",
"We can, based on our discussions of the variance of $\\boldsymbol{\\beta}$ and the mean value, assume that the prior for the values $\\boldsymbol{\\beta}$ is given by a Gaussian with mean value zero and variance $\\tau^2$, that is"
"We can, based on our discussions of the variance of $\\boldsymbol{\\beta}$ and\n",
"the mean value, assume that the prior for the values $\\boldsymbol{\\beta}$ is\n",
"given by a Gaussian with mean value zero and variance $\\tau^2$, that\n",
"is"
]
},
{
@@ -4662,7 +5056,7 @@
"We can now optimize this quantity with respect to $\\boldsymbol{\\beta}$. As we\n",
"did for OLS, this is most conveniently done by taking the negative\n",
"logarithm of the posterior probability. Doing so and leaving out the\n",
"constants terms that do not depend on $\\beta$, we have"
"terms that do not depend on $\\beta$, we have"
]
},
{
@@ -4763,7 +5157,16 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"which is our Lasso cost function!"
"which is our Lasso cost function! \n",
"\n",
"\n",
"Plotting these prior functions shows us that we can use the parameter\n",
"$\\lambda$ to shrink or increase the role of a given parameter\n",
"$\\beta_j$. The variance for the Laplace distribution is\n",
"$2\\tau^2=1/\\lambda$ while for the Gaussian distribution it is\n",
"$\\sigma^2=1/(2\\lambda)$. Thus, increasing the variance means\n",
"decreasing $\\lambda$ and shrinking the variance means increasing\n",
"$\\lamdbda$. When we increase $\\lambda$, this corresponds to shrinking the role of less important features (small singular values)."
]
}
],
@@ -2418,16 +2418,86 @@ plt.ylabel('MSE')
plt.legend()
plt.show()
How can we understand this? **More text to be added**.
How can we understand this?
Let us write out the values of the coefficients $\beta_i$ as functions
of the polynomial degree and noise. We will focus only on the Ridge
results and some few selected values of the hyperparameter $\lambda$.
If we don't include any noise and run this code for different values
of the polynomial degree, we notice that the results for $\beta_i$ do
not show great changes from one order to the next. This is an
indication that for higher polynomial orders, our parameters become
less important.
If we however add noise, what happens is that the polynomial fit is
trying to adjust the fit to traverse in the best possible way all data
points. This can lead to large fluctuations in the parameters
$\beta_i$ as functions of polynomial order. It will also be reflected
in a larger value of the variance of each parameter $\beta_i$. What
Ridge regression (and Lasso as well) are doing then is to try to
quench the fluctuations in the parameters of $\beta_i$ which have a
large variance (normally for higher orders in the polynomial).
import numpy as np
import pandas as pd
from IPython.display import display
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn import linear_model
# Make data set.
n = 1000
x = np.random.rand(n)
y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)+ np.random.randn(n)
Maxpolydegree = 5
X = np.zeros((len(x),Maxpolydegree))
X[:,0] = 1.0
for polydegree in range(1, Maxpolydegree):
for degree in range(polydegree):
X[:,degree] = x**(degree)
# We split the data in test and training data
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)
# Decide which values of lambda to use
nlambdas = 5
lambdas = np.logspace(-3, 2, nlambdas)
for i in range(nlambdas):
lmb = lambdas[i]
# Make the fit using Ridge only
RegRidge = linear_model.Ridge(lmb,fit_intercept=False)
RegRidge.fit(X_train,y_train)
# and then make the prediction
ypredictRidge = RegRidge.predict(X_test)
Coeffs = np.array(RegRidge.coef_)
BetaValues = pd.DataFrame(Coeffs)
BetaValues.columns = ['beta']
display(BetaValues)
As an exercise, repeat these calculations with ordinary least squares
only with and without noise. Calculate thereafter the variance of the
parameters $\beta_j$ as function of polynomial order and of the added
noise. Here we recommend to use $\sigma^2=1$ as variance for the
added noise (which follows a normal distribution with mean value zero).
Comment your results. If you have a large noise term, do the parameters $\beta_j$ vary more as function
model complexity? And what about their variance?
## Linking Bayes' Theorem with Ridge and Lasso Regression
We have seen that Ridge regression suppresses those features which
have a small singular value. This corresponds to a feature which exhibits
a large variance in the parameters $\beta_j$.
Our analysis hitherto has been based on linear algebra. To add to our intuition, we will use
Bayes' theorem in order to deepen our understanding of Ridge and Lasso regression.
Using Bayes' theorem we can gain a better intuition about Ridge and Lasso regression.
For ordinary least squares we postulated that the maximum likelihood for the doamin of events $\boldsymbol{D}$ (one-dimensional case)
For ordinary least squares we postulated that the maximum likelihood for the domain of events $\boldsymbol{D}$ (one-dimensional case)
$$
\boldsymbol{D}=[(x_0,y_0), (x_1,y_1),\dots, (x_{n-1},y_{n-1})],
@@ -2459,7 +2529,10 @@ With the posterior probability defined by a likelihood which we have
already modeled and an unknown prior, we are now ready to make
additional models for the prior.
We can, based on our discussions of the variance of $\boldsymbol{\beta}$ and the mean value, assume that the prior for the values $\boldsymbol{\beta}$ is given by a Gaussian with mean value zero and variance $\tau^2$, that is
We can, based on our discussions of the variance of $\boldsymbol{\beta}$ and
the mean value, assume that the prior for the values $\boldsymbol{\beta}$ is
given by a Gaussian with mean value zero and variance $\tau^2$, that
is
$$
p(\boldsymbol{\beta})=\prod_{j=0}^{p-1}\exp{\left(-\frac{\beta_j^2}{2\tau^2}\right)}.
@@ -2474,7 +2547,7 @@ $$
We can now optimize this quantity with respect to $\boldsymbol{\beta}$. As we
did for OLS, this is most conveniently done by taking the negative
logarithm of the posterior probability. Doing so and leaving out the
constants terms that do not depend on $\beta$, we have
terms that do not depend on $\beta$, we have
$$
C(\boldsymbol{\beta})=\frac{\vert\vert (\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta})\vert\vert_2^2}{2\sigma^2}+\frac{1}{2\tau^2}\vert\vert\boldsymbol{\beta}\vert\vert_2^2,
@@ -2515,4 +2588,13 @@ $$
C(\boldsymbol{\beta}=\frac{\vert\vert (\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta})\vert\vert_2^2}{2\sigma^2}+\lambda\vert\vert\boldsymbol{\beta}\vert\vert_1,
$$
which is our Lasso cost function!
which is our Lasso cost function!
Plotting these prior functions shows us that we can use the parameter
$\lambda$ to shrink or increase the role of a given parameter
$\beta_j$. The variance for the Laplace distribution is
$2\tau^2=1/\lambda$ while for the Gaussian distribution it is
$\sigma^2=1/(2\lambda)$. Thus, increasing the variance means
decreasing $\lambda$ and shrinking the variance means increasing
$\lamdbda$. When we increase $\lambda$, this corresponds to shrinking the role of less important features (small singular values).
@@ -211,7 +211,7 @@
"metadata": {},
"source": [
"This means that the estimator of the regression parameters is unbiased.\n",
"\n",
"v\n",
"We can also calculate the variance\n",
"\n",
"The variance of $\\boldsymbol{\\beta}$ is"
@@ -397,10 +397,10 @@
"name": "stdout",
"output_type": "stream",
"text": [
"Runtime: 0.137546 sec\n",
"Runtime: 0.141725 sec\n",
"Jackknife Statistics :\n",
"original bias std. error\n",
" 99.9031 99.8931 0.149233\n"
" 100.148 100.138 0.1492\n"
]
}
],
@@ -514,7 +514,239 @@
"you use the estimators corresponding to the statistic of interest. For\n",
"example, if you are interested in estimating the variance of $\\widehat\n",
"\\beta$, apply the etsimator $\\widehat \\sigma^2$ to the values\n",
"$\\widehat \\beta ^*$.\n",
"$\\widehat \\beta^*$.\n",
"\n",
"Before we proceed however, we need to remind ourselves about a central\n",
"theorem in statistics, namely the so-called **central limit theorem**.\n",
"This theorem plays a central role in understanding why the Bootstrap\n",
"(and other resampling methods) work so well on independent and\n",
"identically distributed variables.\n",
"\n",
"\n",
"Suppose we have a PDF $p(x)$ from which we generate a series $N$\n",
"of averages $\\langle x_i \\rangle$. Each mean value $\\langle x_i \\rangle$\n",
"is viewed as the average of a specific measurement, e.g., throwing \n",
"dice 100 times and then taking the average value, or producing a certain\n",
"amount of random numbers. \n",
"For notational ease, we set $\\langle x_i \\rangle=x_i$ in the discussion\n",
"which follows. \n",
"\n",
"If we compute the mean $z$ of $m$ such mean values $x_i$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"z=\\frac{x_1+x_2+\\dots+x_m}{m},\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"the question we pose is which is the PDF of the new variable $z$.\n",
"\n",
"\n",
"The probability of obtaining an average value $z$ is the product of the \n",
"probabilities of obtaining arbitrary individual mean values $x_i$,\n",
"but with the constraint that the average is $z$. We can express this through\n",
"the following expression"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\tilde{p}(z)=\\int dx_1p(x_1)\\int dx_2p(x_2)\\dots\\int dx_mp(x_m)\n",
" \\delta(z-\\frac{x_1+x_2+\\dots+x_m}{m}),\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"where the $\\delta$-function enbodies the constraint that the mean is $z$.\n",
"All measurements that lead to each individual $x_i$ are expected to\n",
"be independent, which in turn means that we can express $\\tilde{p}$ as the \n",
"product of individual $p(x_i)$. The independence assumption is important in the derivation of the central limit theorem.\n",
"\n",
"\n",
"\n",
"If we use the integral expression for the $\\delta$-function"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\delta(z-\\frac{x_1+x_2+\\dots+x_m}{m})=\\frac{1}{2\\pi}\\int_{-\\infty}^{\\infty}\n",
" dq\\exp{\\left(iq(z-\\frac{x_1+x_2+\\dots+x_m}{m})\\right)},\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"and inserting $e^{i\\mu q-i\\mu q}$ where $\\mu$ is the mean value\n",
"we arrive at"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\tilde{p}(z)=\\frac{1}{2\\pi}\\int_{-\\infty}^{\\infty}\n",
" dq\\exp{\\left(iq(z-\\mu)\\right)}\\left[\\int_{-\\infty}^{\\infty}\n",
" dxp(x)\\exp{\\left(iq(\\mu-x)/m\\right)}\\right]^m,\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"with the integral over $x$ resulting in"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\int_{-\\infty}^{\\infty}dxp(x)\\exp{\\left(iq(\\mu-x)/m\\right)}=\n",
" \\int_{-\\infty}^{\\infty}dxp(x)\n",
" \\left[1+\\frac{iq(\\mu-x)}{m}-\\frac{q^2(\\mu-x)^2}{2m^2}+\\dots\\right].\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The second term on the rhs disappears since this is just the mean and \n",
"employing the definition of $\\sigma^2$ we have"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\int_{-\\infty}^{\\infty}dxp(x)e^{\\left(iq(\\mu-x)/m\\right)}=\n",
" 1-\\frac{q^2\\sigma^2}{2m^2}+\\dots,\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"resulting in"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\left[\\int_{-\\infty}^{\\infty}dxp(x)\\exp{\\left(iq(\\mu-x)/m\\right)}\\right]^m\\approx\n",
" \\left[1-\\frac{q^2\\sigma^2}{2m^2}+\\dots \\right]^m,\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"and in the limit $m\\rightarrow \\infty$ we obtain"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\tilde{p}(z)=\\frac{1}{\\sqrt{2\\pi}(\\sigma/\\sqrt{m})}\n",
" \\exp{\\left(-\\frac{(z-\\mu)^2}{2(\\sigma/\\sqrt{m})^2}\\right)},\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"which is the normal distribution with variance\n",
"$\\sigma^2_m=\\sigma^2/m$, where $\\sigma$ is the variance of the PDF $p(x)$\n",
"and $\\mu$ is also the mean of the PDF $p(x)$. \n",
"\n",
"\n",
"Thus, the central limit theorem states that the PDF $\\tilde{p}(z)$ of\n",
"the average of $m$ random values corresponding to a PDF $p(x)$ \n",
"is a normal distribution whose mean is the \n",
"mean value of the PDF $p(x)$ and whose variance is the variance\n",
"of the PDF $p(x)$ divided by $m$, the number of values used to compute $z$.\n",
"\n",
"The central limit theorem leads to the well-known expression for the\n",
"standard deviation, given by"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\sigma_m=\n",
"\\frac{\\sigma}{\\sqrt{m}}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The latter is true only if the average value is known exactly. This is obtained in the limit\n",
"$m\\rightarrow \\infty$ only. Because the mean and the variance are measured quantities we obtain \n",
"the familiar expression in statistics"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\sigma_m\\approx \n",
"\\frac{\\sigma}{\\sqrt{m-1}}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In many cases however the above estimate for the standard deviation,\n",
"in particular if correlations are strong, may be too simplistic. Keep\n",
"in mind that we have assumed that the variables $x$ are independent\n",
"and identically distributed. This is obviously not always the\n",
"case. For example, the random numbers (or better pseudorandom numbers)\n",
"we generate in various calculations do always exhibit some\n",
"correlations.\n",
"\n",
"\n",
"\n",
"The theorem is satisfied by a large class of PDFs. Note however that for a\n",
"finite $m$, it is not always possible to find a closed form /analytic expression for\n",
"$\\tilde{p}(x)$.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
@@ -668,7 +900,7 @@
},
"metadata": {
"filenames": {
"image/png": "/Users/mhjensen/Teaching/MachineLearning/doc/LectureNotes/_build/jupyter_execute/chapter3_37_1.png"
"image/png": "/Users/mhjensen/Teaching/MachineLearning/doc/LectureNotes/_build/jupyter_execute/chapter3_57_1.png"
},
"needs_background": "light"
},
@@ -786,15 +1018,14 @@
"Error: 0.03781367141738898\n",
"Bias^2: 0.03365768507152761\n",
"Var: 0.004155986345861379\n",
"0.03781367141738898 >= 0.03365768507152761 + 0.004155986345861379 = 0.03781367141738899\n",
"Polynomial degree:"
"0.03781367141738898 >= 0.03365768507152761 + 0.004155986345861379 = 0.03781367141738899\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" 7\n",
"Polynomial degree: 7\n",
"Error: 0.027609773491022498\n",
"Bias^2: 0.02299949826036597\n",
"Var: 0.004610275230656537\n",
@@ -803,20 +1034,19 @@
"Error: 0.017355848195591973\n",
"Bias^2: 0.010331721306655588\n",
"Var: 0.007024126888936384\n",
"0.017355848195591973 >= 0.010331721306655588 + 0.007024126888936384 = 0.017355848195591973\n",
"Polynomial degree: 9\n",
"Error: 0.026605727637189085\n",
"Bias^2: 0.010018312644140933\n",
"Var: 0.016587414993048166\n",
"0.026605727637189085 >= 0.010018312644140933 + 0.016587414993048166 = 0.0266057276371891\n",
"Polynomial degree:"
"0.017355848195591973 >= 0.010331721306655588 + 0.007024126888936384 = 0.017355848195591973\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" 10\n",
"Polynomial degree: 9\n",
"Error: 0.026605727637189085\n",
"Bias^2: 0.010018312644140933\n",
"Var: 0.016587414993048166\n",
"0.026605727637189085 >= 0.010018312644140933 + 0.016587414993048166 = 0.0266057276371891\n",
"Polynomial degree: 10\n",
"Error: 0.021592704588043153\n",
"Bias^2: 0.010516485576652981\n",
"Var: 0.011076219011390184\n",
@@ -830,7 +1060,13 @@
"Error: 0.1154777721897675\n",
"Bias^2: 0.01628578269590588\n",
"Var: 0.09919198949386163\n",
"0.1154777721897675 >= 0.01628578269590588 + 0.09919198949386163 = 0.11547777218976751\n",
"0.1154777721897675 >= 0.01628578269590588 + 0.09919198949386163 = 0.11547777218976751\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Polynomial degree: 13\n",
"Error: 0.22842468702166951\n",
"Bias^2: 0.01975416527163567\n",
@@ -847,7 +1083,7 @@
},
"metadata": {
"filenames": {
"image/png": "/Users/mhjensen/Teaching/MachineLearning/doc/LectureNotes/_build/jupyter_execute/chapter3_38_4.png"
"image/png": "/Users/mhjensen/Teaching/MachineLearning/doc/LectureNotes/_build/jupyter_execute/chapter3_58_5.png"
},
"needs_background": "light"
},
@@ -984,7 +1220,7 @@
},
"metadata": {
"filenames": {
"image/png": "/Users/mhjensen/Teaching/MachineLearning/doc/LectureNotes/_build/jupyter_execute/chapter3_40_1.png"
"image/png": "/Users/mhjensen/Teaching/MachineLearning/doc/LectureNotes/_build/jupyter_execute/chapter3_60_1.png"
}
},
"output_type": "display_data"
@@ -1099,112 +1335,118 @@
"Mean squared error on test data: 5.98822371\n",
"Degree of polynomial: 6\n",
"Mean squared error on training data: 3.66204648\n",
"Mean squared error on test data: 8.14812206\n",
"Degree of polynomial: 7\n",
"Mean squared error on training data: 0.47075725\n",
"Mean squared error on test data: 2.00607783\n"
"Mean squared error on test data: 8.14812206\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Degree of polynomial: 7\n",
"Mean squared error on training data: 0.47075725\n",
"Mean squared error on test data: 2.00607783\n",
"Degree of polynomial: 8\n",
"Mean squared error on training data: 0.04912436\n",
"Mean squared error on test data: 0.21596432\n",
"Degree of polynomial: 9\n",
"Mean squared error on training data: 0.02522069\n",
"Mean squared error on test data: 0.08576932\n",
"Degree of polynomial: 10\n",
"Mean squared error on training data: 0.02511518\n",
"Mean squared error on test data: 1.20015436\n"
"Mean squared error on test data: 0.08576932\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Degree of polynomial: 10\n",
"Mean squared error on training data: 0.02511518\n",
"Mean squared error on test data: 1.20015436\n",
"Degree of polynomial: 11\n",
"Mean squared error on training data: 0.01640891\n",
"Mean squared error on test data: 1.35533773\n",
"Degree of polynomial: 12\n",
"Mean squared error on training data: 0.00813803\n",
"Mean squared error on test data: 0.17446471\n",
"Mean squared error on test data: 0.17446471\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Degree of polynomial: 13\n",
"Mean squared error on training data: 0.00759119\n",
"Mean squared error on test data: 1.08131003\n",
"Degree of polynomial: 14\n",
"Mean squared error on training data: 0.00472199\n",
"Mean squared error on test data: 0.81333793\n"
"Mean squared error on test data: 0.81333793\n",
"Degree of polynomial: 15\n",
"Mean squared error on training data: 0.00410478\n",
"Mean squared error on test data: 92.09145189\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Degree of polynomial: 15\n",
"Mean squared error on training data: 0.00410478\n",
"Mean squared error on test data: 92.09145189\n",
"Degree of polynomial: 16\n",
"Mean squared error on training data: 0.00315593\n",
"Mean squared error on test data: 234.39716546\n",
"Degree of polynomial: 17\n",
"Mean squared error on training data: 0.00242998\n",
"Mean squared error on test data: 1271.05295709\n"
"Mean squared error on test data: 1271.05295709\n",
"Degree of polynomial: 18\n",
"Mean squared error on training data: 0.00228740\n",
"Mean squared error on test data: 108.42208194\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Degree of polynomial: 18\n",
"Mean squared error on training data: 0.00228740\n",
"Mean squared error on test data: 108.42208194\n",
"Degree of polynomial: 19\n",
"Mean squared error on training data: 0.00156372\n",
"Mean squared error on test data: 1388.41078073\n",
"Degree of polynomial: 20\n",
"Mean squared error on training data: 0.00137982\n",
"Mean squared error on test data: 1761.43341615\n"
"Mean squared error on test data: 1761.43341615\n",
"Degree of polynomial: 21\n",
"Mean squared error on training data: 0.00118170\n",
"Mean squared error on test data: 15061.31603087\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Degree of polynomial: 21\n",
"Mean squared error on training data: 0.00118170\n",
"Mean squared error on test data: 15061.31603087\n",
"Degree of polynomial: 22\n",
"Mean squared error on training data: 0.00092354\n",
"Mean squared error on test data: 890.63488525\n",
"Degree of polynomial: 23\n",
"Mean squared error on training data: 0.00085887\n",
"Mean squared error on test data: 5483.16796929\n"
"Mean squared error on test data: 5483.16796929\n",
"Degree of polynomial: 24\n",
"Mean squared error on training data: 0.00084589\n",
"Mean squared error on test data: 1695.57143061\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Degree of polynomial: 24\n",
"Mean squared error on training data: 0.00084589\n",
"Mean squared error on test data: 1695.57143061\n",
"Degree of polynomial: 25\n",
"Mean squared error on training data: 0.00078806\n",
"Mean squared error on test data: 131343.30655001\n",
"Degree of polynomial: 26\n",
"Mean squared error on training data: 0.00076916\n",
"Mean squared error on test data: 17709.14370264\n"
"Mean squared error on test data: 17709.14370264\n",
"Degree of polynomial: 27\n",
"Mean squared error on training data: 0.00068970\n",
"Mean squared error on test data: 2975.38903780\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Degree of polynomial: 27\n",
"Mean squared error on training data: 0.00068970\n",
"Mean squared error on test data: 2975.38903780\n",
"Degree of polynomial: 28\n",
"Mean squared error on training data: 0.00062588\n",
"Mean squared error on test data: 3848.64522721\n",
@@ -1232,7 +1474,7 @@
},
"metadata": {
"filenames": {
"image/png": "/Users/mhjensen/Teaching/MachineLearning/doc/LectureNotes/_build/jupyter_execute/chapter3_41_10.png"
"image/png": "/Users/mhjensen/Teaching/MachineLearning/doc/LectureNotes/_build/jupyter_execute/chapter3_61_11.png"
},
"needs_background": "light"
},
@@ -1427,7 +1669,7 @@
},
"metadata": {
"filenames": {
"image/png": "/Users/mhjensen/Teaching/MachineLearning/doc/LectureNotes/_build/jupyter_execute/chapter3_47_0.png"
"image/png": "/Users/mhjensen/Teaching/MachineLearning/doc/LectureNotes/_build/jupyter_execute/chapter3_67_0.png"
},
"needs_background": "light"
},
@@ -1558,7 +1800,7 @@
},
"metadata": {
"filenames": {
"image/png": "/Users/mhjensen/Teaching/MachineLearning/doc/LectureNotes/_build/jupyter_execute/chapter3_49_1.png"
"image/png": "/Users/mhjensen/Teaching/MachineLearning/doc/LectureNotes/_build/jupyter_execute/chapter3_69_1.png"
},
"needs_background": "light"
},
@@ -1652,7 +1894,7 @@
},
"metadata": {
"filenames": {
"image/png": "/Users/mhjensen/Teaching/MachineLearning/doc/LectureNotes/_build/jupyter_execute/chapter3_50_0.png"
"image/png": "/Users/mhjensen/Teaching/MachineLearning/doc/LectureNotes/_build/jupyter_execute/chapter3_70_0.png"
},
"needs_background": "light"
},
@@ -1780,7 +2022,7 @@
},
"metadata": {
"filenames": {
"image/png": "/Users/mhjensen/Teaching/MachineLearning/doc/LectureNotes/_build/jupyter_execute/chapter3_54_0.png"
"image/png": "/Users/mhjensen/Teaching/MachineLearning/doc/LectureNotes/_build/jupyter_execute/chapter3_74_0.png"
},
"needs_background": "light"
},
@@ -145,7 +145,7 @@ $$
$$
This means that the estimator of the regression parameters is unbiased.
v
We can also calculate the variance
The variance of $\boldsymbol{\beta}$ is
@@ -372,7 +372,139 @@ histogram of the relative frequency of $\widehat{\beta}^*$. Instead
you use the estimators corresponding to the statistic of interest. For
example, if you are interested in estimating the variance of $\widehat
\beta$, apply the etsimator $\widehat \sigma^2$ to the values
$\widehat \beta ^*$.
$\widehat \beta^*$.
Before we proceed however, we need to remind ourselves about a central
theorem in statistics, namely the so-called **central limit theorem**.
This theorem plays a central role in understanding why the Bootstrap
(and other resampling methods) work so well on independent and
identically distributed variables.
Suppose we have a PDF $p(x)$ from which we generate a series $N$
of averages $\langle x_i \rangle$. Each mean value $\langle x_i \rangle$
is viewed as the average of a specific measurement, e.g., throwing
dice 100 times and then taking the average value, or producing a certain
amount of random numbers.
For notational ease, we set $\langle x_i \rangle=x_i$ in the discussion
which follows.
If we compute the mean $z$ of $m$ such mean values $x_i$
$$
z=\frac{x_1+x_2+\dots+x_m}{m},
$$
the question we pose is which is the PDF of the new variable $z$.
The probability of obtaining an average value $z$ is the product of the
probabilities of obtaining arbitrary individual mean values $x_i$,
but with the constraint that the average is $z$. We can express this through
the following expression
$$
\tilde{p}(z)=\int dx_1p(x_1)\int dx_2p(x_2)\dots\int dx_mp(x_m)
\delta(z-\frac{x_1+x_2+\dots+x_m}{m}),
$$
where the $\delta$-function enbodies the constraint that the mean is $z$.
All measurements that lead to each individual $x_i$ are expected to
be independent, which in turn means that we can express $\tilde{p}$ as the
product of individual $p(x_i)$. The independence assumption is important in the derivation of the central limit theorem.
If we use the integral expression for the $\delta$-function
$$
\delta(z-\frac{x_1+x_2+\dots+x_m}{m})=\frac{1}{2\pi}\int_{-\infty}^{\infty}
dq\exp{\left(iq(z-\frac{x_1+x_2+\dots+x_m}{m})\right)},
$$
and inserting $e^{i\mu q-i\mu q}$ where $\mu$ is the mean value
we arrive at
$$
\tilde{p}(z)=\frac{1}{2\pi}\int_{-\infty}^{\infty}
dq\exp{\left(iq(z-\mu)\right)}\left[\int_{-\infty}^{\infty}
dxp(x)\exp{\left(iq(\mu-x)/m\right)}\right]^m,
$$
with the integral over $x$ resulting in
$$
\int_{-\infty}^{\infty}dxp(x)\exp{\left(iq(\mu-x)/m\right)}=
\int_{-\infty}^{\infty}dxp(x)
\left[1+\frac{iq(\mu-x)}{m}-\frac{q^2(\mu-x)^2}{2m^2}+\dots\right].
$$
The second term on the rhs disappears since this is just the mean and
employing the definition of $\sigma^2$ we have
$$
\int_{-\infty}^{\infty}dxp(x)e^{\left(iq(\mu-x)/m\right)}=
1-\frac{q^2\sigma^2}{2m^2}+\dots,
$$
resulting in
$$
\left[\int_{-\infty}^{\infty}dxp(x)\exp{\left(iq(\mu-x)/m\right)}\right]^m\approx
\left[1-\frac{q^2\sigma^2}{2m^2}+\dots \right]^m,
$$
and in the limit $m\rightarrow \infty$ we obtain
$$
\tilde{p}(z)=\frac{1}{\sqrt{2\pi}(\sigma/\sqrt{m})}
\exp{\left(-\frac{(z-\mu)^2}{2(\sigma/\sqrt{m})^2}\right)},
$$
which is the normal distribution with variance
$\sigma^2_m=\sigma^2/m$, where $\sigma$ is the variance of the PDF $p(x)$
and $\mu$ is also the mean of the PDF $p(x)$.
Thus, the central limit theorem states that the PDF $\tilde{p}(z)$ of
the average of $m$ random values corresponding to a PDF $p(x)$
is a normal distribution whose mean is the
mean value of the PDF $p(x)$ and whose variance is the variance
of the PDF $p(x)$ divided by $m$, the number of values used to compute $z$.
The central limit theorem leads to the well-known expression for the
standard deviation, given by
$$
\sigma_m=
\frac{\sigma}{\sqrt{m}}.
$$
The latter is true only if the average value is known exactly. This is obtained in the limit
$m\rightarrow \infty$ only. Because the mean and the variance are measured quantities we obtain
the familiar expression in statistics
$$
\sigma_m\approx
\frac{\sigma}{\sqrt{m-1}}.
$$
In many cases however the above estimate for the standard deviation,
in particular if correlations are strong, may be too simplistic. Keep
in mind that we have assumed that the variables $x$ are independent
and identically distributed. This is obviously not always the
case. For example, the random numbers (or better pseudorandom numbers)
we generate in various calculations do always exhibit some
correlations.
The theorem is satisfied by a large class of PDFs. Note however that for a
finite $m$, it is not always possible to find a closed form /analytic expression for
$\tilde{p}(x)$.
Binary file not shown.

After

Width:  |  Height:  |  Size: 23 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 16 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 37 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 16 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 14 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 12 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 10 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 32 KiB

+104 -7
View File
@@ -4079,16 +4079,101 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"How can we understand this? **More text to be added**.\n",
"How can we understand this?\n",
"\n",
"Let us write out the values of the coefficients $\\beta_i$ as functions\n",
"of the polynomial degree and noise. We will focus only on the Ridge\n",
"results and some few selected values of the hyperparameter $\\lambda$.\n",
"\n",
"If we don't include any noise and run this code for different values\n",
"of the polynomial degree, we notice that the results for $\\beta_i$ do\n",
"not show great changes from one order to the next. This is an\n",
"indication that for higher polynomial orders, our parameters become\n",
"less important.\n",
"\n",
"If we however add noise, what happens is that the polynomial fit is\n",
"trying to adjust the fit to traverse in the best possible way all data\n",
"points. This can lead to large fluctuations in the parameters\n",
"$\\beta_i$ as functions of polynomial order. It will also be reflected\n",
"in a larger value of the variance of each parameter $\\beta_i$. What\n",
"Ridge regression (and Lasso as well) are doing then is to try to\n",
"quench the fluctuations in the parameters of $\\beta_i$ which have a\n",
"large variance (normally for higher orders in the polynomial)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false,
"editable": true
},
"outputs": [],
"source": [
"import numpy as np\n",
"import pandas as pd\n",
"from IPython.display import display\n",
"import matplotlib.pyplot as plt\n",
"from sklearn.model_selection import train_test_split\n",
"from sklearn import linear_model\n",
"\n",
"# Make data set.\n",
"n = 1000\n",
"x = np.random.rand(n)\n",
"y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)+ np.random.randn(n)\n",
"\n",
"Maxpolydegree = 5\n",
"X = np.zeros((len(x),Maxpolydegree))\n",
"X[:,0] = 1.0\n",
"\n",
"for polydegree in range(1, Maxpolydegree):\n",
" for degree in range(polydegree):\n",
" X[:,degree] = x**(degree)\n",
"\n",
"\n",
"# We split the data in test and training data\n",
"X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)\n",
"\n",
"# Decide which values of lambda to use\n",
"nlambdas = 5\n",
"lambdas = np.logspace(-3, 2, nlambdas)\n",
"for i in range(nlambdas):\n",
" lmb = lambdas[i]\n",
" # Make the fit using Ridge only\n",
" RegRidge = linear_model.Ridge(lmb,fit_intercept=False)\n",
" RegRidge.fit(X_train,y_train)\n",
" # and then make the prediction\n",
" ypredictRidge = RegRidge.predict(X_test)\n",
" Coeffs = np.array(RegRidge.coef_)\n",
" BetaValues = pd.DataFrame(Coeffs)\n",
" BetaValues.columns = ['beta']\n",
" display(BetaValues)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As an exercise, repeat these calculations with ordinary least squares\n",
"only with and without noise. Calculate thereafter the variance of the\n",
"parameters $\\beta_j$ as function of polynomial order and of the added\n",
"noise. Here we recommend to use $\\sigma^2=1$ as variance for the\n",
"added noise (which follows a normal distribution with mean value zero).\n",
"Comment your results. If you have a large noise term, do the parameters $\\beta_j$ vary more as function\n",
"model complexity? And what about their variance? \n",
"\n",
"\n",
"\n",
"\n",
"## Linking Bayes' Theorem with Ridge and Lasso Regression\n",
"\n",
"We have seen that Ridge regression suppresses those features which\n",
"have a small singular value. This corresponds to a feature which exhibits\n",
"a large variance in the parameters $\\beta_j$.\n",
"Our analysis hitherto has been based on linear algebra. To add to our intuition, we will use\n",
"Bayes' theorem in order to deepen our understanding of Ridge and Lasso regression. \n",
"\n",
"Using Bayes' theorem we can gain a better intuition about Ridge and Lasso regression. \n",
"\n",
"For ordinary least squares we postulated that the maximum likelihood for the doamin of events $\\boldsymbol{D}$ (one-dimensional case)"
"For ordinary least squares we postulated that the maximum likelihood for the domain of events $\\boldsymbol{D}$ (one-dimensional case)"
]
},
{
@@ -4160,7 +4245,10 @@
"already modeled and an unknown prior, we are now ready to make\n",
"additional models for the prior.\n",
"\n",
"We can, based on our discussions of the variance of $\\boldsymbol{\\beta}$ and the mean value, assume that the prior for the values $\\boldsymbol{\\beta}$ is given by a Gaussian with mean value zero and variance $\\tau^2$, that is"
"We can, based on our discussions of the variance of $\\boldsymbol{\\beta}$ and\n",
"the mean value, assume that the prior for the values $\\boldsymbol{\\beta}$ is\n",
"given by a Gaussian with mean value zero and variance $\\tau^2$, that\n",
"is"
]
},
{
@@ -4195,7 +4283,7 @@
"We can now optimize this quantity with respect to $\\boldsymbol{\\beta}$. As we\n",
"did for OLS, this is most conveniently done by taking the negative\n",
"logarithm of the posterior probability. Doing so and leaving out the\n",
"constants terms that do not depend on $\\beta$, we have"
"terms that do not depend on $\\beta$, we have"
]
},
{
@@ -4296,7 +4384,16 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"which is our Lasso cost function!"
"which is our Lasso cost function! \n",
"\n",
"\n",
"Plotting these prior functions shows us that we can use the parameter\n",
"$\\lambda$ to shrink or increase the role of a given parameter\n",
"$\\beta_j$. The variance for the Laplace distribution is\n",
"$2\\tau^2=1/\\lambda$ while for the Gaussian distribution it is\n",
"$\\sigma^2=1/(2\\lambda)$. Thus, increasing the variance means\n",
"decreasing $\\lambda$ and shrinking the variance means increasing\n",
"$\\lamdbda$. When we increase $\\lambda$, this corresponds to shrinking the role of less important features (small singular values)."
]
}
],
+234 -2
View File
@@ -211,7 +211,7 @@
"metadata": {},
"source": [
"This means that the estimator of the regression parameters is unbiased.\n",
"\n",
"v\n",
"We can also calculate the variance\n",
"\n",
"The variance of $\\boldsymbol{\\beta}$ is"
@@ -503,7 +503,239 @@
"you use the estimators corresponding to the statistic of interest. For\n",
"example, if you are interested in estimating the variance of $\\widehat\n",
"\\beta$, apply the etsimator $\\widehat \\sigma^2$ to the values\n",
"$\\widehat \\beta ^*$.\n",
"$\\widehat \\beta^*$.\n",
"\n",
"Before we proceed however, we need to remind ourselves about a central\n",
"theorem in statistics, namely the so-called **central limit theorem**.\n",
"This theorem plays a central role in understanding why the Bootstrap\n",
"(and other resampling methods) work so well on independent and\n",
"identically distributed variables.\n",
"\n",
"\n",
"Suppose we have a PDF $p(x)$ from which we generate a series $N$\n",
"of averages $\\langle x_i \\rangle$. Each mean value $\\langle x_i \\rangle$\n",
"is viewed as the average of a specific measurement, e.g., throwing \n",
"dice 100 times and then taking the average value, or producing a certain\n",
"amount of random numbers. \n",
"For notational ease, we set $\\langle x_i \\rangle=x_i$ in the discussion\n",
"which follows. \n",
"\n",
"If we compute the mean $z$ of $m$ such mean values $x_i$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"z=\\frac{x_1+x_2+\\dots+x_m}{m},\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"the question we pose is which is the PDF of the new variable $z$.\n",
"\n",
"\n",
"The probability of obtaining an average value $z$ is the product of the \n",
"probabilities of obtaining arbitrary individual mean values $x_i$,\n",
"but with the constraint that the average is $z$. We can express this through\n",
"the following expression"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\tilde{p}(z)=\\int dx_1p(x_1)\\int dx_2p(x_2)\\dots\\int dx_mp(x_m)\n",
" \\delta(z-\\frac{x_1+x_2+\\dots+x_m}{m}),\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"where the $\\delta$-function enbodies the constraint that the mean is $z$.\n",
"All measurements that lead to each individual $x_i$ are expected to\n",
"be independent, which in turn means that we can express $\\tilde{p}$ as the \n",
"product of individual $p(x_i)$. The independence assumption is important in the derivation of the central limit theorem.\n",
"\n",
"\n",
"\n",
"If we use the integral expression for the $\\delta$-function"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\delta(z-\\frac{x_1+x_2+\\dots+x_m}{m})=\\frac{1}{2\\pi}\\int_{-\\infty}^{\\infty}\n",
" dq\\exp{\\left(iq(z-\\frac{x_1+x_2+\\dots+x_m}{m})\\right)},\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"and inserting $e^{i\\mu q-i\\mu q}$ where $\\mu$ is the mean value\n",
"we arrive at"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\tilde{p}(z)=\\frac{1}{2\\pi}\\int_{-\\infty}^{\\infty}\n",
" dq\\exp{\\left(iq(z-\\mu)\\right)}\\left[\\int_{-\\infty}^{\\infty}\n",
" dxp(x)\\exp{\\left(iq(\\mu-x)/m\\right)}\\right]^m,\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"with the integral over $x$ resulting in"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\int_{-\\infty}^{\\infty}dxp(x)\\exp{\\left(iq(\\mu-x)/m\\right)}=\n",
" \\int_{-\\infty}^{\\infty}dxp(x)\n",
" \\left[1+\\frac{iq(\\mu-x)}{m}-\\frac{q^2(\\mu-x)^2}{2m^2}+\\dots\\right].\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The second term on the rhs disappears since this is just the mean and \n",
"employing the definition of $\\sigma^2$ we have"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\int_{-\\infty}^{\\infty}dxp(x)e^{\\left(iq(\\mu-x)/m\\right)}=\n",
" 1-\\frac{q^2\\sigma^2}{2m^2}+\\dots,\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"resulting in"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\left[\\int_{-\\infty}^{\\infty}dxp(x)\\exp{\\left(iq(\\mu-x)/m\\right)}\\right]^m\\approx\n",
" \\left[1-\\frac{q^2\\sigma^2}{2m^2}+\\dots \\right]^m,\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"and in the limit $m\\rightarrow \\infty$ we obtain"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\tilde{p}(z)=\\frac{1}{\\sqrt{2\\pi}(\\sigma/\\sqrt{m})}\n",
" \\exp{\\left(-\\frac{(z-\\mu)^2}{2(\\sigma/\\sqrt{m})^2}\\right)},\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"which is the normal distribution with variance\n",
"$\\sigma^2_m=\\sigma^2/m$, where $\\sigma$ is the variance of the PDF $p(x)$\n",
"and $\\mu$ is also the mean of the PDF $p(x)$. \n",
"\n",
"\n",
"Thus, the central limit theorem states that the PDF $\\tilde{p}(z)$ of\n",
"the average of $m$ random values corresponding to a PDF $p(x)$ \n",
"is a normal distribution whose mean is the \n",
"mean value of the PDF $p(x)$ and whose variance is the variance\n",
"of the PDF $p(x)$ divided by $m$, the number of values used to compute $z$.\n",
"\n",
"The central limit theorem leads to the well-known expression for the\n",
"standard deviation, given by"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\sigma_m=\n",
"\\frac{\\sigma}{\\sqrt{m}}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The latter is true only if the average value is known exactly. This is obtained in the limit\n",
"$m\\rightarrow \\infty$ only. Because the mean and the variance are measured quantities we obtain \n",
"the familiar expression in statistics"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\sigma_m\\approx \n",
"\\frac{\\sigma}{\\sqrt{m-1}}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In many cases however the above estimate for the standard deviation,\n",
"in particular if correlations are strong, may be too simplistic. Keep\n",
"in mind that we have assumed that the variables $x$ are independent\n",
"and identically distributed. This is obviously not always the\n",
"case. For example, the random numbers (or better pseudorandom numbers)\n",
"we generate in various calculations do always exhibit some\n",
"correlations.\n",
"\n",
"\n",
"\n",
"The theorem is satisfied by a large class of PDFs. Note however that for a\n",
"finite $m$, it is not always possible to find a closed form /analytic expression for\n",
"$\\tilde{p}(x)$.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",