- -
We will now couple the discussions of ordinary least squares, Ridge -and Lasso regression with a statistical interpretation, that is we -move from a linear algebra analysis to a statistical analysis. In -particular, we will focus on what the regularization terms can result -in. We will amongst other things show that the regularization -parameter can reduce considerably the variance of the parameters -\( \beta \). +
Hitherto we have discussed Ridge and Lasso regression in terms of a +linear analysis. This may to many of you feel rather technical and +perhaps not that intuitive. The question is whether we can develop a +more intuitive way of understanding what Ridge and Lasso express.
-The -advantage of doing linear regression is that we actually end up with -analytical expressions for several statistical quantities. -Standard least squares and Ridge regression allow us to -derive quantities like the variance and other expectation values in a -rather straightforward way. +
Before we proceed let us perform a Ridge, Lasso and OLS analysis of a polynomial fit.
+ +We will play around with a study of the values for the optimal +parameters \( \boldsymbol{\beta} \) using OLS, Ridge and Lasso regression. For +OLS, you will notice as function of the noise and polynomial degree, +that the parameters \( \beta \) will fluctuate from order to order in the +polynomial fit and that for larger and larger polynomial degrees of freedom, the parameters will tend to increase in value for OLS.
-It is assumed that \( \varepsilon_i -\sim \mathcal{N}(0, \sigma^2) \) and the \( \varepsilon_{i} \) are -independent, i.e.: -
-$$ -\begin{align*} -\mbox{Cov}(\varepsilon_{i_1}, -\varepsilon_{i_2}) & = \left\{ \begin{array}{lcc} \sigma^2 & \mbox{if} -& i_1 = i_2, \\ 0 & \mbox{if} & i_1 \not= i_2. \end{array} \right. -\end{align*} -$$ +For Ridge and Lasso regression, the higher order parameters will typically be reduced, providing thereby less fluctuations from one order to another one.
-The randomness of \( \varepsilon_i \) implies that -\( \mathbf{y}_i \) is also a random variable. In particular, -\( \mathbf{y}_i \) is normally distributed, because \( \varepsilon_i \sim -\mathcal{N}(0, \sigma^2) \) and \( \mathbf{X}_{i,\ast} \, \boldsymbol{\beta} \) is a -non-random scalar. To specify the parameters of the distribution of -\( \mathbf{y}_i \) we need to calculate its first two moments. -
-Recall that \( \boldsymbol{X} \) is a matrix of dimensionality \( n\times p \). The -notation above \( \mathbf{X}_{i,\ast} \) means that we are looking at the -row number \( i \) and perform a sum over all values \( p \). -
+ +import numpy as np
+import matplotlib.pyplot as plt
+from sklearn.model_selection import train_test_split
+from sklearn import linear_model
+
+def R2(y_data, y_model):
+ return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2)
+def MSE(y_data,y_model):
+ n = np.size(y_model)
+ return np.sum((y_data-y_model)**2)/n
+
+# Make data set.
+n = 10000
+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):
+ X[:,polydegree] = x**(polydegree)
+
+# 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)
+
+# matrix inversion to find beta
+OLSbeta = np.linalg.pinv(X_train.T @ X_train) @ X_train.T @ y_train
+print(OLSbeta)
+ypredictOLS = X_test @ OLSbeta
+print("Test MSE OLS")
+print(MSE(y_test,ypredictOLS))
+# Repeat now for Lasso and Ridge regression and various values of the regularization parameter using Scikit-Learn
+# Decide which values of lambda to use
+nlambdas = 4
+MSERidgePredict = np.zeros(nlambdas)
+MSELassoPredict = np.zeros(nlambdas)
+lambdas = np.logspace(-3, 1, nlambdas)
+for i in range(nlambdas):
+ lmb = lambdas[i]
+ # Make the fit using Ridge and Lasso
+ RegRidge = linear_model.Ridge(lmb,fit_intercept=False)
+ RegRidge.fit(X_train,y_train)
+ RegLasso = linear_model.Lasso(lmb,fit_intercept=False)
+ RegLasso.fit(X_train,y_train)
+ # and then make the prediction
+ ypredictRidge = RegRidge.predict(X_test)
+ ypredictLasso = RegLasso.predict(X_test)
+ # Compute the MSE and print it
+ MSERidgePredict[i] = MSE(y_test,ypredictRidge)
+ MSELassoPredict[i] = MSE(y_test,ypredictLasso)
+ print(lmb,RegRidge.coef_)
+ print(lmb,RegLasso.coef_)
+# Now plot the results
+plt.figure()
+plt.plot(np.log10(lambdas), MSERidgePredict, 'b', label = 'MSE Ridge Test')
+plt.plot(np.log10(lambdas), MSELassoPredict, 'r', label = 'MSE Lasso Test')
+plt.xlabel('log10(lambda)')
+plt.ylabel('MSE')
+plt.legend()
+plt.show()
+
+How can we understand this?
@@ -334,7 +393,7 @@ row number \( i \) and perform a sum over all values \( p \).
- -
The assumption we have made here can be summarized as (and this is going to be useful when we discuss the bias-variance trade off) -that there exists a function \( f(\boldsymbol{x}) \) and a normal distributed error \( \boldsymbol{\varepsilon}\sim \mathcal{N}(0, \sigma^2) \) -which describe our data +
We will now couple the discussions of ordinary least squares, Ridge +and Lasso regression with a statistical interpretation, that is we +move from a linear algebra analysis to a statistical analysis. In +particular, we will focus on what the regularization terms can result +in. We will amongst other things show that the regularization +parameter can reduce considerably the variance of the parameters +\( \beta \). +
+ +The +advantage of doing linear regression is that we actually end up with +analytical expressions for several statistical quantities. +Standard least squares and Ridge regression allow us to +derive quantities like the variance and other expectation values in a +rather straightforward way. +
+ +It is assumed that \( \varepsilon_i +\sim \mathcal{N}(0, \sigma^2) \) and the \( \varepsilon_{i} \) are +independent, i.e.:
$$ -\boldsymbol{y} = f(\boldsymbol{x})+\boldsymbol{\varepsilon} +\begin{align*} +\mbox{Cov}(\varepsilon_{i_1}, +\varepsilon_{i_2}) & = \left\{ \begin{array}{lcc} \sigma^2 & \mbox{if} +& i_1 = i_2, \\ 0 & \mbox{if} & i_1 \not= i_2. \end{array} \right. +\end{align*} $$ -We approximate this function with our model from the solution of the linear regression equations, that is our -function \( f \) is approximated by \( \boldsymbol{\tilde{y}} \) where we want to minimize \( (\boldsymbol{y}-\boldsymbol{\tilde{y}})^2 \), our MSE, with +
The randomness of \( \varepsilon_i \) implies that +\( \mathbf{y}_i \) is also a random variable. In particular, +\( \mathbf{y}_i \) is normally distributed, because \( \varepsilon_i \sim +\mathcal{N}(0, \sigma^2) \) and \( \mathbf{X}_{i,\ast} \, \boldsymbol{\beta} \) is a +non-random scalar. To specify the parameters of the distribution of +\( \mathbf{y}_i \) we need to calculate its first two moments.
-$$ -\boldsymbol{\tilde{y}} = \boldsymbol{X}\boldsymbol{\beta}. -$$ +Recall that \( \boldsymbol{X} \) is a matrix of dimensionality \( n\times p \). The +notation above \( \mathbf{X}_{i,\ast} \) means that we are looking at the +row number \( i \) and perform a sum over all values \( p \). +
@@ -308,7 +332,7 @@ $$
-
We can calculate the expectation value of \( \boldsymbol{y} \) for a given element \( i \)
-$$ -\begin{align*} -\mathbb{E}(y_i) & = -\mathbb{E}(\mathbf{X}_{i, \ast} \, \boldsymbol{\beta}) + \mathbb{E}(\varepsilon_i) -\, \, \, = \, \, \, \mathbf{X}_{i, \ast} \, \beta, -\end{align*} -$$ - -while -its variance is +
The assumption we have made here can be summarized as (and this is going to be useful when we discuss the bias-variance trade off) +that there exists a function \( f(\boldsymbol{x}) \) and a normal distributed error \( \boldsymbol{\varepsilon}\sim \mathcal{N}(0, \sigma^2) \) +which describe our data
$$ -\begin{align*} \mbox{Var}(y_i) & = \mathbb{E} \{ [y_i -- \mathbb{E}(y_i)]^2 \} \, \, \, = \, \, \, \mathbb{E} ( y_i^2 ) - -[\mathbb{E}(y_i)]^2 \\ & = \mathbb{E} [ ( \mathbf{X}_{i, \ast} \, -\beta + \varepsilon_i )^2] - ( \mathbf{X}_{i, \ast} \, \boldsymbol{\beta})^2 \\ & -= \mathbb{E} [ ( \mathbf{X}_{i, \ast} \, \boldsymbol{\beta})^2 + 2 \varepsilon_i -\mathbf{X}_{i, \ast} \, \boldsymbol{\beta} + \varepsilon_i^2 ] - ( \mathbf{X}_{i, -\ast} \, \beta)^2 \\ & = ( \mathbf{X}_{i, \ast} \, \boldsymbol{\beta})^2 + 2 -\mathbb{E}(\varepsilon_i) \mathbf{X}_{i, \ast} \, \boldsymbol{\beta} + -\mathbb{E}(\varepsilon_i^2 ) - ( \mathbf{X}_{i, \ast} \, \boldsymbol{\beta})^2 -\\ & = \mathbb{E}(\varepsilon_i^2 ) \, \, \, = \, \, \, -\mbox{Var}(\varepsilon_i) \, \, \, = \, \, \, \sigma^2. -\end{align*} +\boldsymbol{y} = f(\boldsymbol{x})+\boldsymbol{\varepsilon} $$ -Hence, \( y_i \sim \mathcal{N}( \mathbf{X}_{i, \ast} \, \boldsymbol{\beta}, \sigma^2) \), that is \( \boldsymbol{y} \) follows a normal distribution with -mean value \( \boldsymbol{X}\boldsymbol{\beta} \) and variance \( \sigma^2 \) (not be confused with the singular values of the SVD). +
We approximate this function with our model from the solution of the linear regression equations, that is our +function \( f \) is approximated by \( \boldsymbol{\tilde{y}} \) where we want to minimize \( (\boldsymbol{y}-\boldsymbol{\tilde{y}})^2 \), our MSE, with
+$$ +\boldsymbol{\tilde{y}} = \boldsymbol{X}\boldsymbol{\beta}. +$$ +@@ -323,7 +306,7 @@ mean value \( \boldsymbol{X}\boldsymbol{\beta} \) and variance \( \sigma^2 \) (n
-
With the OLS expressions for the optimal parameters \( \boldsymbol{\hat{\beta}} \) we can evaluate the expectation value
+We can calculate the expectation value of \( \boldsymbol{y} \) for a given element \( i \)
$$ -\mathbb{E}(\boldsymbol{\hat{\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}. +\begin{align*} +\mathbb{E}(y_i) & = +\mathbb{E}(\mathbf{X}_{i, \ast} \, \boldsymbol{\beta}) + \mathbb{E}(\varepsilon_i) +\, \, \, = \, \, \, \mathbf{X}_{i, \ast} \, \beta, +\end{align*} $$ -This means that the estimator of the regression parameters is unbiased.
- -We can also calculate the variance
- -The variance of the optimal value \( \boldsymbol{\hat{\beta}} \) is
-$$ -\begin{eqnarray*} -\mbox{Var}(\boldsymbol{\hat{\beta}}) & = & \mathbb{E} \{ [\boldsymbol{\beta} - \mathbb{E}(\boldsymbol{\beta})] [\boldsymbol{\beta} - \mathbb{E}(\boldsymbol{\beta})]^{T} \} -\\ -& = & \mathbb{E} \{ [(\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \mathbf{Y} - \boldsymbol{\beta}] \, [(\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \mathbf{Y} - \boldsymbol{\beta}]^{T} \} -\\ -% & = & \mathbb{E} \{ [(\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \mathbf{Y}] \, [(\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \mathbf{Y}]^{T} \} - \boldsymbol{\beta} \, \boldsymbol{\beta}^{T} -% \\ -% & = & \mathbb{E} \{ (\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \mathbf{Y} \, \mathbf{Y}^{T} \, \mathbf{X} \, (\mathbf{X}^{T} \mathbf{X})^{-1} \} - \boldsymbol{\beta} \, \boldsymbol{\beta}^{T} -% \\ -& = & (\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \, \mathbb{E} \{ \mathbf{Y} \, \mathbf{Y}^{T} \} \, \mathbf{X} \, (\mathbf{X}^{T} \mathbf{X})^{-1} - \boldsymbol{\beta} \, \boldsymbol{\beta}^{T} -\\ -& = & (\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \, \{ \mathbf{X} \, \boldsymbol{\beta} \, \boldsymbol{\beta}^{T} \, \mathbf{X}^{T} + \sigma^2 \} \, \mathbf{X} \, (\mathbf{X}^{T} \mathbf{X})^{-1} - \boldsymbol{\beta} \, \boldsymbol{\beta}^{T} -% \\ -% & = & (\mathbf{X}^T \mathbf{X})^{-1} \, \mathbf{X}^T \, \mathbf{X} \, \boldsymbol{\beta} \, \boldsymbol{\beta}^T \, \mathbf{X}^T \, \mathbf{X} \, (\mathbf{X}^T % \mathbf{X})^{-1} -% \\ -% & & + \, \, \sigma^2 \, (\mathbf{X}^T \mathbf{X})^{-1} \, \mathbf{X}^T \, \mathbf{X} \, (\mathbf{X}^T \mathbf{X})^{-1} - \boldsymbol{\beta} \boldsymbol{\beta}^T -\\ -& = & \boldsymbol{\beta} \, \boldsymbol{\beta}^{T} + \sigma^2 \, (\mathbf{X}^{T} \mathbf{X})^{-1} - \boldsymbol{\beta} \, \boldsymbol{\beta}^{T} -\, \, \, = \, \, \, \sigma^2 \, (\mathbf{X}^{T} \mathbf{X})^{-1}, -\end{eqnarray*} -$$ - -where we have used that \( \mathbb{E} (\mathbf{Y} \mathbf{Y}^{T}) = -\mathbf{X} \, \boldsymbol{\beta} \, \boldsymbol{\beta}^{T} \, \mathbf{X}^{T} + -\sigma^2 \, \mathbf{I}_{nn} \). From \( \mbox{Var}(\boldsymbol{\beta}) = \sigma^2 -\, (\mathbf{X}^{T} \mathbf{X})^{-1} \), one obtains an estimate of the -variance of the estimate of the \( j \)-th regression coefficient: -\( \boldsymbol{\sigma}^2 (\boldsymbol{\beta}_j ) = \boldsymbol{\sigma}^2 [(\mathbf{X}^{T} \mathbf{X})^{-1}]_{jj} \). This may be used to -construct a confidence interval for the estimates. +
while +its variance is
- -In a similar way, we can obtain analytical expressions for say the -expectation values of the parameters \( \boldsymbol{\beta} \) and their variance -when we employ Ridge regression, allowing us again to define a confidence interval. -
- -It is rather straightforward to show that
$$ -\mathbb{E} \big[ \boldsymbol{\beta}^{\mathrm{Ridge}} \big]=(\mathbf{X}^{T} \mathbf{X} + \lambda \mathbf{I}_{pp})^{-1} (\mathbf{X}^{\top} \mathbf{X})\boldsymbol{\beta}^{\mathrm{OLS}}. +\begin{align*} \mbox{Var}(y_i) & = \mathbb{E} \{ [y_i +- \mathbb{E}(y_i)]^2 \} \, \, \, = \, \, \, \mathbb{E} ( y_i^2 ) - +[\mathbb{E}(y_i)]^2 \\ & = \mathbb{E} [ ( \mathbf{X}_{i, \ast} \, +\beta + \varepsilon_i )^2] - ( \mathbf{X}_{i, \ast} \, \boldsymbol{\beta})^2 \\ & += \mathbb{E} [ ( \mathbf{X}_{i, \ast} \, \boldsymbol{\beta})^2 + 2 \varepsilon_i +\mathbf{X}_{i, \ast} \, \boldsymbol{\beta} + \varepsilon_i^2 ] - ( \mathbf{X}_{i, +\ast} \, \beta)^2 \\ & = ( \mathbf{X}_{i, \ast} \, \boldsymbol{\beta})^2 + 2 +\mathbb{E}(\varepsilon_i) \mathbf{X}_{i, \ast} \, \boldsymbol{\beta} + +\mathbb{E}(\varepsilon_i^2 ) - ( \mathbf{X}_{i, \ast} \, \boldsymbol{\beta})^2 +\\ & = \mathbb{E}(\varepsilon_i^2 ) \, \, \, = \, \, \, +\mbox{Var}(\varepsilon_i) \, \, \, = \, \, \, \sigma^2. +\end{align*} $$ -We see clearly that -\( \mathbb{E} \big[ \boldsymbol{\beta}^{\mathrm{Ridge}} \big] \not= \boldsymbol{\beta}^{\mathrm{OLS}} \) for any \( \lambda > 0 \). We say then that the ridge estimator is biased. -
- -We can also compute the variance as
- -$$ -\mbox{Var}[\boldsymbol{\beta}^{\mathrm{Ridge}}]=\sigma^2[ \mathbf{X}^{T} \mathbf{X} + \lambda \mathbf{I} ]^{-1} \mathbf{X}^{T} \mathbf{X} \{ [ \mathbf{X}^{\top} \mathbf{X} + \lambda \mathbf{I} ]^{-1}\}^{T}, -$$ - -and it is easy to see that if the parameter \( \lambda \) goes to infinity then the variance of Ridge parameters \( \boldsymbol{\beta} \) goes to zero.
- -With this, we can compute the difference
- -$$ -\mbox{Var}[\boldsymbol{\beta}^{\mathrm{OLS}}]-\mbox{Var}(\boldsymbol{\beta}^{\mathrm{Ridge}})=\sigma^2 [ \mathbf{X}^{T} \mathbf{X} + \lambda \mathbf{I} ]^{-1}[ 2\lambda\mathbf{I} + \lambda^2 (\mathbf{X}^{T} \mathbf{X})^{-1} ] \{ [ \mathbf{X}^{T} \mathbf{X} + \lambda \mathbf{I} ]^{-1}\}^{T}. -$$ - -The difference is non-negative definite since each component of the -matrix product is non-negative definite. -This means the variance we obtain with the standard OLS will always for \( \lambda > 0 \) be larger than the variance of \( \boldsymbol{\beta} \) obtained with the Ridge estimator. This has interesting consequences when we discuss the so-called bias-variance trade-off below. +
Hence, \( y_i \sim \mathcal{N}( \mathbf{X}_{i, \ast} \, \boldsymbol{\beta}, \sigma^2) \), that is \( \boldsymbol{y} \) follows a normal distribution with +mean value \( \boldsymbol{X}\boldsymbol{\beta} \) and variance \( \sigma^2 \) (not be confused with the singular values of the SVD).
@@ -367,7 +321,7 @@ This means the variance we obtain with the standard OLS will always for \( \lamb
-
Our basic assumption when we derived the OLS equations was to assume -that our output is determined by a given continuous function -\( f(\boldsymbol{x}) \) and a random noise \( \boldsymbol{\epsilon} \) given by the normal -distribution with zero mean value and an undetermined variance -\( \sigma^2 \). -
- -We found above that the outputs \( \boldsymbol{y} \) have a mean value given by -\( \boldsymbol{X}\hat{\boldsymbol{\beta}} \) and variance \( \sigma^2 \). Since the entries to -the design matrix are not stochastic variables, we can assume that the -probability distribution of our targets is also a normal distribution -but now with mean value \( \boldsymbol{X}\hat{\boldsymbol{\beta}} \). This means that a -single output \( y_i \) is given by the Gaussian distribution -
+With the OLS expressions for the optimal parameters \( \boldsymbol{\hat{\beta}} \) we can evaluate the expectation value
$$ -y_i\sim \mathcal{N}(\boldsymbol{X}_{i,*}\boldsymbol{\beta}, \sigma^2)=\frac{1}{\sqrt{2\pi\sigma^2}}\exp{\left[-\frac{(y_i-\boldsymbol{X}_{i,*}\boldsymbol{\beta})^2}{2\sigma^2}\right]}. +\mathbb{E}(\boldsymbol{\hat{\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}. $$ +This means that the estimator of the regression parameters is unbiased.
+ +We can also calculate the variance
+ +The variance of the optimal value \( \boldsymbol{\hat{\beta}} \) is
+$$ +\begin{eqnarray*} +\mbox{Var}(\boldsymbol{\hat{\beta}}) & = & \mathbb{E} \{ [\boldsymbol{\beta} - \mathbb{E}(\boldsymbol{\beta})] [\boldsymbol{\beta} - \mathbb{E}(\boldsymbol{\beta})]^{T} \} +\\ +& = & \mathbb{E} \{ [(\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \mathbf{Y} - \boldsymbol{\beta}] \, [(\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \mathbf{Y} - \boldsymbol{\beta}]^{T} \} +\\ +% & = & \mathbb{E} \{ [(\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \mathbf{Y}] \, [(\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \mathbf{Y}]^{T} \} - \boldsymbol{\beta} \, \boldsymbol{\beta}^{T} +% \\ +% & = & \mathbb{E} \{ (\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \mathbf{Y} \, \mathbf{Y}^{T} \, \mathbf{X} \, (\mathbf{X}^{T} \mathbf{X})^{-1} \} - \boldsymbol{\beta} \, \boldsymbol{\beta}^{T} +% \\ +& = & (\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \, \mathbb{E} \{ \mathbf{Y} \, \mathbf{Y}^{T} \} \, \mathbf{X} \, (\mathbf{X}^{T} \mathbf{X})^{-1} - \boldsymbol{\beta} \, \boldsymbol{\beta}^{T} +\\ +& = & (\mathbf{X}^{T} \mathbf{X})^{-1} \, \mathbf{X}^{T} \, \{ \mathbf{X} \, \boldsymbol{\beta} \, \boldsymbol{\beta}^{T} \, \mathbf{X}^{T} + \sigma^2 \} \, \mathbf{X} \, (\mathbf{X}^{T} \mathbf{X})^{-1} - \boldsymbol{\beta} \, \boldsymbol{\beta}^{T} +% \\ +% & = & (\mathbf{X}^T \mathbf{X})^{-1} \, \mathbf{X}^T \, \mathbf{X} \, \boldsymbol{\beta} \, \boldsymbol{\beta}^T \, \mathbf{X}^T \, \mathbf{X} \, (\mathbf{X}^T % \mathbf{X})^{-1} +% \\ +% & & + \, \, \sigma^2 \, (\mathbf{X}^T \mathbf{X})^{-1} \, \mathbf{X}^T \, \mathbf{X} \, (\mathbf{X}^T \mathbf{X})^{-1} - \boldsymbol{\beta} \boldsymbol{\beta}^T +\\ +& = & \boldsymbol{\beta} \, \boldsymbol{\beta}^{T} + \sigma^2 \, (\mathbf{X}^{T} \mathbf{X})^{-1} - \boldsymbol{\beta} \, \boldsymbol{\beta}^{T} +\, \, \, = \, \, \, \sigma^2 \, (\mathbf{X}^{T} \mathbf{X})^{-1}, +\end{eqnarray*} +$$ + +where we have used that \( \mathbb{E} (\mathbf{Y} \mathbf{Y}^{T}) = +\mathbf{X} \, \boldsymbol{\beta} \, \boldsymbol{\beta}^{T} \, \mathbf{X}^{T} + +\sigma^2 \, \mathbf{I}_{nn} \). From \( \mbox{Var}(\boldsymbol{\beta}) = \sigma^2 +\, (\mathbf{X}^{T} \mathbf{X})^{-1} \), one obtains an estimate of the +variance of the estimate of the \( j \)-th regression coefficient: +\( \boldsymbol{\sigma}^2 (\boldsymbol{\beta}_j ) = \boldsymbol{\sigma}^2 [(\mathbf{X}^{T} \mathbf{X})^{-1}]_{jj} \). This may be used to +construct a confidence interval for the estimates. +
+ +In a similar way, we can obtain analytical expressions for say the +expectation values of the parameters \( \boldsymbol{\beta} \) and their variance +when we employ Ridge regression, allowing us again to define a confidence interval. +
+ +It is rather straightforward to show that
+$$ +\mathbb{E} \big[ \boldsymbol{\beta}^{\mathrm{Ridge}} \big]=(\mathbf{X}^{T} \mathbf{X} + \lambda \mathbf{I}_{pp})^{-1} (\mathbf{X}^{\top} \mathbf{X})\boldsymbol{\beta}^{\mathrm{OLS}}. +$$ + +We see clearly that +\( \mathbb{E} \big[ \boldsymbol{\beta}^{\mathrm{Ridge}} \big] \not= \boldsymbol{\beta}^{\mathrm{OLS}} \) for any \( \lambda > 0 \). We say then that the ridge estimator is biased. +
+ +We can also compute the variance as
+ +$$ +\mbox{Var}[\boldsymbol{\beta}^{\mathrm{Ridge}}]=\sigma^2[ \mathbf{X}^{T} \mathbf{X} + \lambda \mathbf{I} ]^{-1} \mathbf{X}^{T} \mathbf{X} \{ [ \mathbf{X}^{\top} \mathbf{X} + \lambda \mathbf{I} ]^{-1}\}^{T}, +$$ + +and it is easy to see that if the parameter \( \lambda \) goes to infinity then the variance of Ridge parameters \( \boldsymbol{\beta} \) goes to zero.
+ +With this, we can compute the difference
+ +$$ +\mbox{Var}[\boldsymbol{\beta}^{\mathrm{OLS}}]-\mbox{Var}(\boldsymbol{\beta}^{\mathrm{Ridge}})=\sigma^2 [ \mathbf{X}^{T} \mathbf{X} + \lambda \mathbf{I} ]^{-1}[ 2\lambda\mathbf{I} + \lambda^2 (\mathbf{X}^{T} \mathbf{X})^{-1} ] \{ [ \mathbf{X}^{T} \mathbf{X} + \lambda \mathbf{I} ]^{-1}\}^{T}. +$$ + +The difference is non-negative definite since each component of the +matrix product is non-negative definite. +This means the variance we obtain with the standard OLS will always for \( \lambda > 0 \) be larger than the variance of \( \boldsymbol{\beta} \) obtained with the Ridge estimator. This has interesting consequences when we discuss the so-called bias-variance trade-off below. +
@@ -312,7 +365,7 @@ $$
-
We assume now that the various \( y_i \) values are stochastically distributed according to the above Gaussian distribution. -We define this distribution as +
Our basic assumption when we derived the OLS equations was to assume +that our output is determined by a given continuous function +\( f(\boldsymbol{x}) \) and a random noise \( \boldsymbol{\epsilon} \) given by the normal +distribution with zero mean value and an undetermined variance +\( \sigma^2 \).
-$$ -p(y_i, \boldsymbol{X}\vert\boldsymbol{\beta})=\frac{1}{\sqrt{2\pi\sigma^2}}\exp{\left[-\frac{(y_i-\boldsymbol{X}_{i,*}\boldsymbol{\beta})^2}{2\sigma^2}\right]}, -$$ -which reads as finding the likelihood of an event \( y_i \) with the input variables \( \boldsymbol{X} \) given the parameters (to be determined) \( \boldsymbol{\beta} \).
- -Since these events are assumed to be independent and identicall distributed we can build the probability distribution function (PDF) for all possible event \( \boldsymbol{y} \) as the product of the single events, that is we have
- -$$ -p(\boldsymbol{y},\boldsymbol{X}\vert\boldsymbol{\beta})=\prod_{i=0}^{n-1}\frac{1}{\sqrt{2\pi\sigma^2}}\exp{\left[-\frac{(y_i-\boldsymbol{X}_{i,*}\boldsymbol{\beta})^2}{2\sigma^2}\right]}=\prod_{i=0}^{n-1}p(y_i,\boldsymbol{X}\vert\boldsymbol{\beta}). -$$ - -We will write this in a more compact form reserving \( \boldsymbol{D} \) for the domain of events, including the ouputs (targets) and the inputs. That is -in case we have a simple one-dimensional input and output case +
We found above that the outputs \( \boldsymbol{y} \) have a mean value given by +\( \boldsymbol{X}\hat{\boldsymbol{\beta}} \) and variance \( \sigma^2 \). Since the entries to +the design matrix are not stochastic variables, we can assume that the +probability distribution of our targets is also a normal distribution +but now with mean value \( \boldsymbol{X}\hat{\boldsymbol{\beta}} \). This means that a +single output \( y_i \) is given by the Gaussian distribution
+ $$ -\boldsymbol{D}=[(x_0,y_0), (x_1,y_1),\dots, (x_{n-1},y_{n-1})]. +y_i\sim \mathcal{N}(\boldsymbol{X}_{i,*}\boldsymbol{\beta}, \sigma^2)=\frac{1}{\sqrt{2\pi\sigma^2}}\exp{\left[-\frac{(y_i-\boldsymbol{X}_{i,*}\boldsymbol{\beta})^2}{2\sigma^2}\right]}. $$ -In the more general case the various inputs should be replaced by the possible features represented by the input data set \( \boldsymbol{X} \). -We can now rewrite the above probability as -
-$$ -p(\boldsymbol{D}\vert\boldsymbol{\beta})=\prod_{i=0}^{n-1}\frac{1}{\sqrt{2\pi\sigma^2}}\exp{\left[-\frac{(y_i-\boldsymbol{X}_{i,*}\boldsymbol{\beta})^2}{2\sigma^2}\right]}. -$$ - -It is a conditional probability (see below) and reads as the likelihood of a domain of events \( \boldsymbol{D} \) given a set of parameters \( \boldsymbol{\beta} \).
@@ -323,7 +310,7 @@ $$
-
In statistics, maximum likelihood estimation (MLE) is a method of -estimating the parameters of an assumed probability distribution, -given some observed data. This is achieved by maximizing a likelihood -function so that, under the assumed statistical model, the observed -data is the most probable. +
We assume now that the various \( y_i \) values are stochastically distributed according to the above Gaussian distribution. +We define this distribution as
+$$ +p(y_i, \boldsymbol{X}\vert\boldsymbol{\beta})=\frac{1}{\sqrt{2\pi\sigma^2}}\exp{\left[-\frac{(y_i-\boldsymbol{X}_{i,*}\boldsymbol{\beta})^2}{2\sigma^2}\right]}, +$$ -We will assume here that our events are given by the above Gaussian -distribution and we will determine the optimal parameters \( \beta \) by -maximizing the above PDF. However, computing the derivatives of a -product function is cumbersome and can easily lead to overflow and/or -underflowproblems, with potentials for loss of numerical precision. -
+which reads as finding the likelihood of an event \( y_i \) with the input variables \( \boldsymbol{X} \) given the parameters (to be determined) \( \boldsymbol{\beta} \).
-In practice, it is more convenient to maximize the logarithm of the -PDF because it is a monotonically increasing function of the argument. -Alternatively, and this will be our option, we will minimize the -negative of the logarithm since this is a monotonically decreasing -function. -
+Since these events are assumed to be independent and identicall distributed we can build the probability distribution function (PDF) for all possible event \( \boldsymbol{y} \) as the product of the single events, that is we have
-Note also that maximization/minimization of the logarithm of the PDF -is equivalent to the maximization/minimization of the function itself. +$$ +p(\boldsymbol{y},\boldsymbol{X}\vert\boldsymbol{\beta})=\prod_{i=0}^{n-1}\frac{1}{\sqrt{2\pi\sigma^2}}\exp{\left[-\frac{(y_i-\boldsymbol{X}_{i,*}\boldsymbol{\beta})^2}{2\sigma^2}\right]}=\prod_{i=0}^{n-1}p(y_i,\boldsymbol{X}\vert\boldsymbol{\beta}). +$$ + +
We will write this in a more compact form reserving \( \boldsymbol{D} \) for the domain of events, including the ouputs (targets) and the inputs. That is +in case we have a simple one-dimensional input and output case
+$$ +\boldsymbol{D}=[(x_0,y_0), (x_1,y_1),\dots, (x_{n-1},y_{n-1})]. +$$ + +In the more general case the various inputs should be replaced by the possible features represented by the input data set \( \boldsymbol{X} \). +We can now rewrite the above probability as +
+$$ +p(\boldsymbol{D}\vert\boldsymbol{\beta})=\prod_{i=0}^{n-1}\frac{1}{\sqrt{2\pi\sigma^2}}\exp{\left[-\frac{(y_i-\boldsymbol{X}_{i,*}\boldsymbol{\beta})^2}{2\sigma^2}\right]}. +$$ + +It is a conditional probability (see below) and reads as the likelihood of a domain of events \( \boldsymbol{D} \) given a set of parameters \( \boldsymbol{\beta} \).
@@ -317,7 +321,7 @@ is equivalent to the maximization/minimization of the function itself.
-
We could now define a new cost function to minimize, namely the negative logarithm of the above PDF
+In statistics, maximum likelihood estimation (MLE) is a method of +estimating the parameters of an assumed probability distribution, +given some observed data. This is achieved by maximizing a likelihood +function so that, under the assumed statistical model, the observed +data is the most probable. +
-$$ -C(\boldsymbol{\beta}=-\log{\prod_{i=0}^{n-1}p(y_i,\boldsymbol{X}\vert\boldsymbol{\beta})}=-\sum_{i=0}^{n-1}\log{p(y_i,\boldsymbol{X}\vert\boldsymbol{\beta})}, -$$ +We will assume here that our events are given by the above Gaussian +distribution and we will determine the optimal parameters \( \beta \) by +maximizing the above PDF. However, computing the derivatives of a +product function is cumbersome and can easily lead to overflow and/or +underflowproblems, with potentials for loss of numerical precision. +
-which becomes
-$$ -C(\boldsymbol{\beta}=\frac{n}{2}\log{2\pi\sigma^2}+\frac{\vert\vert (\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta})\vert\vert_2^2}{2\sigma^2}. -$$ +In practice, it is more convenient to maximize the logarithm of the +PDF because it is a monotonically increasing function of the argument. +Alternatively, and this will be our option, we will minimize the +negative of the logarithm since this is a monotonically decreasing +function. +
-Taking the derivative of the new cost function with respect to the parameters \( \beta \) we recognize our familiar OLS equation, namely
- -$$ -\boldsymbol{X}^T\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right) =0, -$$ - -which leads to the well-known OLS equation for the optimal paramters \( \beta \)
-$$ -\hat{\boldsymbol{\beta}}^{\mathrm{OLS}}=\left(\boldsymbol{X}^T\boldsymbol{X}\right)^{-1}\boldsymbol{X}^T\boldsymbol{y}! -$$ - -Before we make a similar analysis for Ridge and Lasso regression, we need a short reminder on statistics.
+Note also that maximization/minimization of the logarithm of the PDF +is equivalent to the maximization/minimization of the function itself. +
@@ -316,7 +315,7 @@ $$
-
A central theorem in statistics is Bayes' theorem. This theorem plays a similar role as the good old Pythagoras' theorem in geometry. -Bayes' theorem is extremely simple to derive. But to do so we need some basic axioms from statistics. -
+We could now define a new cost function to minimize, namely the negative logarithm of the above PDF
-Assume we have two domains of events \( X=[x_0,x_1,\dots,x_{n-1}] \) and \( Y=[y_0,y_1,\dots,y_{n-1}] \).
- -We define also the likelihood for \( X \) and \( Y \) as \( p(X) \) and \( p(Y) \) respectively. -The likelihood of a specific event \( x_i \) (or \( y_i \)) is then written as \( p(X=x_i) \) or just \( p(x_i)=p_i \). -
- -where we read \( p(X\vert Y) \) as the likelihood of obtaining \( X \) given \( Y \).
-which becomes
+$$ +C(\boldsymbol{\beta}=\frac{n}{2}\log{2\pi\sigma^2}+\frac{\vert\vert (\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta})\vert\vert_2^2}{2\sigma^2}. +$$ +Taking the derivative of the new cost function with respect to the parameters \( \beta \) we recognize our familiar OLS equation, namely
-If we have independent events then \( p(X,Y)=p(X)p(Y) \).
+$$ +\boldsymbol{X}^T\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right) =0, +$$ + +which leads to the well-known OLS equation for the optimal paramters \( \beta \)
+$$ +\hat{\boldsymbol{\beta}}^{\mathrm{OLS}}=\left(\boldsymbol{X}^T\boldsymbol{X}\right)^{-1}\boldsymbol{X}^T\boldsymbol{y}! +$$ + +Before we make a similar analysis for Ridge and Lasso regression, we need a short reminder on statistics.
@@ -326,7 +314,7 @@ $$
-
A central theorem in statistics is Bayes' theorem. This theorem plays a similar role as the good old Pythagoras' theorem in geometry. +Bayes' theorem is extremely simple to derive. But to do so we need some basic axioms from statistics. +
+ +Assume we have two domains of events \( X=[x_0,x_1,\dots,x_{n-1}] \) and \( Y=[y_0,y_1,\dots,y_{n-1}] \).
+ +We define also the likelihood for \( X \) and \( Y \) as \( p(X) \) and \( p(Y) \) respectively. +The likelihood of a specific event \( x_i \) (or \( y_i \)) is then written as \( p(X=x_i) \) or just \( p(x_i)=p_i \). +
-The marginal probability is defined in terms of only one of the set of variables \( X,Y \). For a discrete probability we have
where we read \( p(X\vert Y) \) as the likelihood of obtaining \( X \) given \( Y \).
+If we have independent events then \( p(X,Y)=p(X)p(Y) \).
+diff --git a/doc/pub/week36/html/._week36-bs040.html b/doc/pub/week36/html/._week36-bs040.html index e01284206..2e3a0fd4b 100644 --- a/doc/pub/week36/html/._week36-bs040.html +++ b/doc/pub/week36/html/._week36-bs040.html @@ -116,6 +116,10 @@ doconce format html week36.do.txt --html_style=bootstrap --pygments_html_style=d 2, None, 'important-technicalities-more-on-rescaling-data'), + ('Test Function for what happens with OLS, Ridge and Lasso', + 2, + None, + 'test-function-for-what-happens-with-ols-ridge-and-lasso'), ('Linking the regression analysis with a statistical ' 'interpretation', 2, @@ -163,11 +167,6 @@ doconce format html week36.do.txt --html_style=bootstrap --pygments_html_style=d 2, None, 'bayes-theorem-and-ridge-and-lasso-regression'), - ('Test Function for what happens with OLS, Ridge and Lasso', - 2, - None, - 'test-function-for-what-happens-with-ols-ridge-and-lasso'), - ("Invoking Bayes' theorem", 2, None, 'invoking-bayes-theorem'), ('Ridge and Bayes', 2, None, 'ridge-and-bayes'), ('Lasso and Bayes', 2, None, 'lasso-and-bayes')]} end of tocinfo --> @@ -234,26 +233,25 @@ MathJax.Hub.Config({
-
The conditional probability, if \( p(Y) > 0 \), is
+The marginal probability is defined in terms of only one of the set of variables \( X,Y \). For a discrete probability we have
-
If we combine the conditional probability with the marginal probability and the standard product rule, we have
+The conditional probability, if \( p(Y) > 0 \), is
+which we can rewrite as
- -$$ -p(X\vert Y)= \frac{p(X,Y)}{\sum_{i=0}^{n-1}p(Y\vert X=x_i)p(x_i)}=\frac{p(Y\vert X)p(X)}{\sum_{i=0}^{n-1}p(Y\vert X=x_i)p(x_i)}, -$$ - -which is Bayes' theorem. It allows us to evaluate the uncertainty in in \( X \) after we have observed \( Y \). We can easily interchange \( X \) with \( Y \).
@@ -302,7 +298,6 @@ $$
-
The quantity \( p(Y\vert X) \) on the right-hand side of the theorem is -evaluated for the observed data \( Y \) and can be viewed as a function of -the parameter space represented by \( X \). This function is not -necesseraly normalized and is normally called the likelihood function. -
+If we combine the conditional probability with the marginal probability and the standard product rule, we have
+$$ +p(X\vert Y)= \frac{p(X,Y)}{p(Y)}, +$$ -The function \( p(X) \) on the right hand side is called the prior while the function on the left hand side is the called the posterior probability. The denominator on the right hand side serves as a normalization factor for the posterior distribution.
+which we can rewrite as
-Let us try to illustrate Bayes' theorem through an example.
+$$ +p(X\vert Y)= \frac{p(X,Y)}{\sum_{i=0}^{n-1}p(Y\vert X=x_i)p(x_i)}=\frac{p(Y\vert X)p(X)}{\sum_{i=0}^{n-1}p(Y\vert X=x_i)p(x_i)}, +$$ + +which is Bayes' theorem. It allows us to evaluate the uncertainty in in \( X \) after we have observed \( Y \). We can easily interchange \( X \) with \( Y \).
@@ -298,7 +299,6 @@ necesseraly normalized and is normally called the likelihood function.
-
Let us suppose that you are undergoing a series of mammography scans in -order to rule out possible breast cancer cases. We define the -sensitivity for a positive event by the variable \( X \). It takes binary -values with \( X=1 \) representing a positive event and \( X=0 \) being a -negative event. We reserve \( Y \) as a classification parameter for -either a negative or a positive breast cancer confirmation. (Short note on wordings: positive here means having breast cancer, although none of us would consider this being a positive thing). +
The quantity \( p(Y\vert X) \) on the right-hand side of the theorem is +evaluated for the observed data \( Y \) and can be viewed as a function of +the parameter space represented by \( X \). This function is not +necesseraly normalized and is normally called the likelihood function.
-We let \( Y=1 \) represent the the case of having breast cancer and \( Y=0 \) as not.
+The function \( p(X) \) on the right hand side is called the prior while the function on the left hand side is the called the posterior probability. The denominator on the right hand side serves as a normalization factor for the posterior distribution.
-Let us assume that if you have breast cancer, the test will be positive with a probability of \( 0.8 \), that is we have
- -$$ -p(X=1\vert Y=1) =0.8. -$$ - -This obviously sounds scary since many would conclude that if the test is positive, there is a likelihood of \( 80\% \) for having cancer. -It is however not correct, as the following Bayesian analysis shows. -
+Let us try to illustrate Bayes' theorem through an example.
@@ -307,7 +295,6 @@ It is however not correct, as the following Bayesian analysis shows.
-
If we look at various national surveys on breast cancer, the general likelihood of developing breast cancer is a very small number. -Let us assume that the prior probability in the population as a whole is +
Let us suppose that you are undergoing a series of mammography scans in +order to rule out possible breast cancer cases. We define the +sensitivity for a positive event by the variable \( X \). It takes binary +values with \( X=1 \) representing a positive event and \( X=0 \) being a +negative event. We reserve \( Y \) as a classification parameter for +either a negative or a positive breast cancer confirmation. (Short note on wordings: positive here means having breast cancer, although none of us would consider this being a positive thing).
-$$ -p(Y=1) =0.004. -$$ +We let \( Y=1 \) represent the the case of having breast cancer and \( Y=0 \) as not.
-We need also to account for the fact that the test may produce a false positive result (false alarm). Let us here assume that we have
-$$ -p(X=1\vert Y=0) =0.1. -$$ - -Using Bayes' theorem we can then find the posterior probability that the person has breast cancer in case of a positive test, that is we can compute
+Let us assume that if you have breast cancer, the test will be positive with a probability of \( 0.8 \), that is we have
$$ -p(Y=1\vert X=1)=\frac{p(X=1\vert Y=1)p(Y=1)}{p(X=1\vert Y=1)p(Y=1)+p(X=1\vert Y=0)p(Y=0)}=\frac{0.8\times 0.004}{0.8\times 0.004+0.1\times 0.996}=0.031. +p(X=1\vert Y=1) =0.8. $$ -That is, in case of a positive test, there is only a \( 3\% \) chance of having breast cancer!
+This obviously sounds scary since many would conclude that if the test is positive, there is a likelihood of \( 80\% \) for having cancer. +It is however not correct, as the following Bayesian analysis shows. +
@@ -307,7 +304,6 @@ $$
-
Hitherto we have discussed Ridge and Lasso regression in terms of a -linear analysis. This may to many of you feel rather technical and -perhaps not that intuitive. The question is whether we can develop a -more intuitive way of understanding what Ridge and Lasso express. +
If we look at various national surveys on breast cancer, the general likelihood of developing breast cancer is a very small number. +Let us assume that the prior probability in the population as a whole is
-Before we proceed let us perform a Ridge, Lasso and OLS analysis of a polynomial fit.
+$$ +p(Y=1) =0.004. +$$ + +We need also to account for the fact that the test may produce a false positive result (false alarm). Let us here assume that we have
+$$ +p(X=1\vert Y=0) =0.1. +$$ + +Using Bayes' theorem we can then find the posterior probability that the person has breast cancer in case of a positive test, that is we can compute
+ +$$ +p(Y=1\vert X=1)=\frac{p(X=1\vert Y=1)p(Y=1)}{p(X=1\vert Y=1)p(Y=1)+p(X=1\vert Y=0)p(Y=0)}=\frac{0.8\times 0.004}{0.8\times 0.004+0.1\times 0.996}=0.031. +$$ + +That is, in case of a positive test, there is only a \( 3\% \) chance of having breast cancer!
@@ -293,7 +304,6 @@ more intuitive way of understanding what Ridge and Lasso express.
-
We will play around with a study of the values for the optimal -parameters \( \boldsymbol{\beta} \) using OLS, Ridge and Lasso regression. For -OLS, you will notice as function of the noise and polynomial degree, -that the parameters \( \beta \) will fluctuate from order to order in the -polynomial fit and that for larger and larger polynomial degrees of freedom, the parameters will tend to increase in value for OLS. -
+Using Bayes' theorem we can gain a better intuition about Ridge and Lasso regression.
-For Ridge and Lasso regression, the higher order parameters will typically be reduced, providing thereby less fluctuations from one order to another one.
+For ordinary least squares we postulated that the maximum likelihood for the doamin of events \( \boldsymbol{D} \) (one-dimensional case)
+$$ +\boldsymbol{D}=[(x_0,y_0), (x_1,y_1),\dots, (x_{n-1},y_{n-1})], +$$ +is given by
+$$ +p(\boldsymbol{D}\vert\boldsymbol{\beta})=\prod_{i=0}^{n-1}\frac{1}{\sqrt{2\pi\sigma^2}}\exp{\left[-\frac{(y_i-\boldsymbol{X}_{i,*}\boldsymbol{\beta})^2}{2\sigma^2}\right]}. +$$ - -import numpy as np
-import matplotlib.pyplot as plt
-from sklearn.model_selection import train_test_split
-from sklearn import linear_model
+In Bayes' theorem this function plays the role of the so-called likelihood. We could now ask the question what is the posterior probability of a parameter set \( \boldsymbol{\beta} \) given a domain of events \( \boldsymbol{D} \)? That is, how can we define the posterior probability
-def R2(y_data, y_model):
- return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2)
-def MSE(y_data,y_model):
- n = np.size(y_model)
- return np.sum((y_data-y_model)**2)/n
+$$
+p(\boldsymbol{\beta}\vert\boldsymbol{D}).
+$$
-# Make data set.
-n = 10000
-x = np.random.rand(n)
-y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)+ np.random.randn(n)
+Bayes' theorem comes to our rescue here since (omitting the normalization constant)
+$$
+p(\boldsymbol{\beta}\vert\boldsymbol{D})\propto p(\boldsymbol{D}\vert\boldsymbol{\beta})p(\boldsymbol{\beta}).
+$$
-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)
-
-# matrix inversion to find beta
-OLSbeta = np.linalg.pinv(X_train.T @ X_train) @ X_train.T @ y_train
-print(OLSbeta)
-ypredictOLS = X_test @ OLSbeta
-print("Test MSE OLS")
-print(MSE(y_test,ypredictOLS))
-# Repeat now for Lasso and Ridge regression and various values of the regularization parameter using Scikit-Learn
-# Decide which values of lambda to use
-nlambdas = 4
-MSERidgePredict = np.zeros(nlambdas)
-MSELassoPredict = np.zeros(nlambdas)
-lambdas = np.logspace(-3, 1, nlambdas)
-for i in range(nlambdas):
- lmb = lambdas[i]
- # Make the fit using Ridge and Lasso
- RegRidge = linear_model.Ridge(lmb,fit_intercept=False)
- RegRidge.fit(X_train,y_train)
- RegLasso = linear_model.Lasso(lmb,fit_intercept=False)
- RegLasso.fit(X_train,y_train)
- # and then make the prediction
- ypredictRidge = RegRidge.predict(X_test)
- ypredictLasso = RegLasso.predict(X_test)
- # Compute the MSE and print it
- MSERidgePredict[i] = MSE(y_test,ypredictRidge)
- MSELassoPredict[i] = MSE(y_test,ypredictLasso)
- print(lmb,RegRidge.coef_)
- print(lmb,RegLasso.coef_)
-# Now plot the results
-plt.figure()
-plt.plot(np.log10(lambdas), MSERidgePredict, 'b', label = 'MSE Ridge Test')
-plt.plot(np.log10(lambdas), MSELassoPredict, 'r', label = 'MSE Lasso Test')
-plt.xlabel('log10(lambda)')
-plt.ylabel('MSE')
-plt.legend()
-plt.show()
-
-How can we understand this?
+We have a model for \( p(\boldsymbol{D}\vert\boldsymbol{\beta}) \) but need one for the prior \( p(\boldsymbol{\beta} \)!
@@ -380,7 +307,6 @@ plt.show()
-
Using Bayes' theorem we can gain a better intuition about Ridge and Lasso regression.
+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. +
-For ordinary least squares we postulated that the maximum likelihood for the doamin of events \( \boldsymbol{D} \) (one-dimensional case)
-$$ -\boldsymbol{D}=[(x_0,y_0), (x_1,y_1),\dots, (x_{n-1},y_{n-1})], -$$ - -is given by
-$$ -p(\boldsymbol{D}\vert\boldsymbol{\beta})=\prod_{i=0}^{n-1}\frac{1}{\sqrt{2\pi\sigma^2}}\exp{\left[-\frac{(y_i-\boldsymbol{X}_{i,*}\boldsymbol{\beta})^2}{2\sigma^2}\right]}. -$$ - -In Bayes' theorem this function plays the role of the so-called likelihood. We could now ask the question what is the posterior probability of a parameter set \( \boldsymbol{\beta} \) given a domain of events \( \boldsymbol{D} \)? That is, how can we define the posterior probability
+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}\vert\boldsymbol{D}). +p(\boldsymbol{\beta})=\prod_{j=0}^{p-1}\exp{\left(-\frac{\beta_j^2}{2\tau^2}\right)}. $$ -Bayes' theorem comes to our rescue here since (omitting the normalization constant)
+Our posterior probability becomes then (omitting the normalization factor which is just a constant)
$$ -p(\boldsymbol{\beta}\vert\boldsymbol{D})\propto p(\boldsymbol{D}\vert\boldsymbol{\beta})p(\boldsymbol{\beta}). +p(\boldsymbol{\beta\vert\boldsymbol{D})}=\prod_{i=0}^{n-1}\frac{1}{\sqrt{2\pi\sigma^2}}\exp{\left[-\frac{(y_i-\boldsymbol{X}_{i,*}\boldsymbol{\beta})^2}{2\sigma^2}\right]}\prod_{j=0}^{p-1}\exp{\left(-\frac{\beta_j^2}{2\tau^2}\right)}. $$ -We have a model for \( p(\boldsymbol{D}\vert\boldsymbol{\beta}) \) but need one for the prior \( p(\boldsymbol{\beta} \)!
+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 +
+ +$$ +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, +$$ + +and replacing \( 1/2\tau^2 \) with \( \lambda \) we have
+ +$$ +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_2^2, +$$ + +which is our Ridge cost function! Nice, isn't it?
@@ -308,7 +315,6 @@ $$
-
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
+To derive the Lasso cost function, we simply replace the Gaussian prior with an exponential distribution (Laplace in this case) with zero mean value, that is
$$ -p(\boldsymbol{\beta})=\prod_{j=0}^{p-1}\exp{\left(-\frac{\beta_j^2}{2\tau^2}\right)}. +p(\boldsymbol{\beta})=\prod_{j=0}^{p-1}\exp{\left(-\frac{\vert\beta_j\vert}{\tau}\right)}. $$Our posterior probability becomes then (omitting the normalization factor which is just a constant)
$$ -p(\boldsymbol{\beta\vert\boldsymbol{D})}=\prod_{i=0}^{n-1}\frac{1}{\sqrt{2\pi\sigma^2}}\exp{\left[-\frac{(y_i-\boldsymbol{X}_{i,*}\boldsymbol{\beta})^2}{2\sigma^2}\right]}\prod_{j=0}^{p-1}\exp{\left(-\frac{\beta_j^2}{2\tau^2}\right)}. +p(\boldsymbol{\beta}\vert\boldsymbol{D})=\prod_{i=0}^{n-1}\frac{1}{\sqrt{2\pi\sigma^2}}\exp{\left[-\frac{(y_i-\boldsymbol{X}_{i,*}\boldsymbol{\beta})^2}{2\sigma^2}\right]}\prod_{j=0}^{p-1}\exp{\left(-\frac{\vert\beta_j\vert}{\tau}\right)}. $$ -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 +
Taking the negative +logarithm of the posterior probability and leaving out the constants 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, +C(\boldsymbol{\beta}=\frac{\vert\vert (\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta})\vert\vert_2^2}{2\sigma^2}+\frac{1}{\tau}\vert\vert\boldsymbol{\beta}\vert\vert_1, $$ -and replacing \( 1/2\tau^2 \) with \( \lambda \) we have
+and replacing \( 1/\tau \) with \( \lambda \) we have
$$ -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_2^2, +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 Ridge cost function! Nice, isn't it?
+which is our Lasso cost function!
@@ -316,8 +308,6 @@ $$
Hitherto we have discussed Ridge and Lasso regression in terms of a +linear analysis. This may to many of you feel rather technical and +perhaps not that intuitive. The question is whether we can develop a +more intuitive way of understanding what Ridge and Lasso express. +
+ +Before we proceed let us perform a Ridge, Lasso and OLS analysis of a polynomial fit.
+ +We will play around with a study of the values for the optimal +parameters \( \boldsymbol{\beta} \) using OLS, Ridge and Lasso regression. For +OLS, you will notice as function of the noise and polynomial degree, +that the parameters \( \beta \) will fluctuate from order to order in the +polynomial fit and that for larger and larger polynomial degrees of freedom, the parameters will tend to increase in value for OLS. +
+ +For Ridge and Lasso regression, the higher order parameters will typically be reduced, providing thereby less fluctuations from one order to another one.
+ + + +import numpy as np
+import matplotlib.pyplot as plt
+from sklearn.model_selection import train_test_split
+from sklearn import linear_model
+
+def R2(y_data, y_model):
+ return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2)
+def MSE(y_data,y_model):
+ n = np.size(y_model)
+ return np.sum((y_data-y_model)**2)/n
+
+# Make data set.
+n = 10000
+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):
+ X[:,polydegree] = x**(polydegree)
+
+# 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)
+
+# matrix inversion to find beta
+OLSbeta = np.linalg.pinv(X_train.T @ X_train) @ X_train.T @ y_train
+print(OLSbeta)
+ypredictOLS = X_test @ OLSbeta
+print("Test MSE OLS")
+print(MSE(y_test,ypredictOLS))
+# Repeat now for Lasso and Ridge regression and various values of the regularization parameter using Scikit-Learn
+# Decide which values of lambda to use
+nlambdas = 4
+MSERidgePredict = np.zeros(nlambdas)
+MSELassoPredict = np.zeros(nlambdas)
+lambdas = np.logspace(-3, 1, nlambdas)
+for i in range(nlambdas):
+ lmb = lambdas[i]
+ # Make the fit using Ridge and Lasso
+ RegRidge = linear_model.Ridge(lmb,fit_intercept=False)
+ RegRidge.fit(X_train,y_train)
+ RegLasso = linear_model.Lasso(lmb,fit_intercept=False)
+ RegLasso.fit(X_train,y_train)
+ # and then make the prediction
+ ypredictRidge = RegRidge.predict(X_test)
+ ypredictLasso = RegLasso.predict(X_test)
+ # Compute the MSE and print it
+ MSERidgePredict[i] = MSE(y_test,ypredictRidge)
+ MSELassoPredict[i] = MSE(y_test,ypredictLasso)
+ print(lmb,RegRidge.coef_)
+ print(lmb,RegLasso.coef_)
+# Now plot the results
+plt.figure()
+plt.plot(np.log10(lambdas), MSERidgePredict, 'b', label = 'MSE Ridge Test')
+plt.plot(np.log10(lambdas), MSELassoPredict, 'r', label = 'MSE Lasso Test')
+plt.xlabel('log10(lambda)')
+plt.ylabel('MSE')
+plt.legend()
+plt.show()
+
+How can we understand this?
+Hitherto we have discussed Ridge and Lasso regression in terms of a -linear analysis. This may to many of you feel rather technical and -perhaps not that intuitive. The question is whether we can develop a -more intuitive way of understanding what Ridge and Lasso express. -
- -Before we proceed let us perform a Ridge, Lasso and OLS analysis of a polynomial fit.
-We will play around with a study of the values for the optimal -parameters \( \boldsymbol{\beta} \) using OLS, Ridge and Lasso regression. For -OLS, you will notice as function of the noise and polynomial degree, -that the parameters \( \beta \) will fluctuate from order to order in the -polynomial fit and that for larger and larger polynomial degrees of freedom, the parameters will tend to increase in value for OLS. -
- -For Ridge and Lasso regression, the higher order parameters will typically be reduced, providing thereby less fluctuations from one order to another one.
- - - -import numpy as np
-import matplotlib.pyplot as plt
-from sklearn.model_selection import train_test_split
-from sklearn import linear_model
-
-def R2(y_data, y_model):
- return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2)
-def MSE(y_data,y_model):
- n = np.size(y_model)
- return np.sum((y_data-y_model)**2)/n
-
-# Make data set.
-n = 10000
-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)
-
-# matrix inversion to find beta
-OLSbeta = np.linalg.pinv(X_train.T @ X_train) @ X_train.T @ y_train
-print(OLSbeta)
-ypredictOLS = X_test @ OLSbeta
-print("Test MSE OLS")
-print(MSE(y_test,ypredictOLS))
-# Repeat now for Lasso and Ridge regression and various values of the regularization parameter using Scikit-Learn
-# Decide which values of lambda to use
-nlambdas = 4
-MSERidgePredict = np.zeros(nlambdas)
-MSELassoPredict = np.zeros(nlambdas)
-lambdas = np.logspace(-3, 1, nlambdas)
-for i in range(nlambdas):
- lmb = lambdas[i]
- # Make the fit using Ridge and Lasso
- RegRidge = linear_model.Ridge(lmb,fit_intercept=False)
- RegRidge.fit(X_train,y_train)
- RegLasso = linear_model.Lasso(lmb,fit_intercept=False)
- RegLasso.fit(X_train,y_train)
- # and then make the prediction
- ypredictRidge = RegRidge.predict(X_test)
- ypredictLasso = RegLasso.predict(X_test)
- # Compute the MSE and print it
- MSERidgePredict[i] = MSE(y_test,ypredictRidge)
- MSELassoPredict[i] = MSE(y_test,ypredictLasso)
- print(lmb,RegRidge.coef_)
- print(lmb,RegLasso.coef_)
-# Now plot the results
-plt.figure()
-plt.plot(np.log10(lambdas), MSERidgePredict, 'b', label = 'MSE Ridge Test')
-plt.plot(np.log10(lambdas), MSELassoPredict, 'r', label = 'MSE Lasso Test')
-plt.xlabel('log10(lambda)')
-plt.ylabel('MSE')
-plt.legend()
-plt.show()
-
-How can we understand this?
-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)
diff --git a/doc/pub/week36/html/week36-solarized.html b/doc/pub/week36/html/week36-solarized.html index dbd9a99a5..24579e5da 100644 --- a/doc/pub/week36/html/week36-solarized.html +++ b/doc/pub/week36/html/week36-solarized.html @@ -143,6 +143,10 @@ div.toc p,a { 2, None, 'important-technicalities-more-on-rescaling-data'), + ('Test Function for what happens with OLS, Ridge and Lasso', + 2, + None, + 'test-function-for-what-happens-with-ols-ridge-and-lasso'), ('Linking the regression analysis with a statistical ' 'interpretation', 2, @@ -190,11 +194,6 @@ div.toc p,a { 2, None, 'bayes-theorem-and-ridge-and-lasso-regression'), - ('Test Function for what happens with OLS, Ridge and Lasso', - 2, - None, - 'test-function-for-what-happens-with-ols-ridge-and-lasso'), - ("Invoking Bayes' theorem", 2, None, 'invoking-bayes-theorem'), ('Ridge and Bayes', 2, None, 'ridge-and-bayes'), ('Lasso and Bayes', 2, None, 'lasso-and-bayes')]} end of tocinfo --> @@ -252,7 +251,7 @@ MathJax.Hub.Config({Hitherto we have discussed Ridge and Lasso regression in terms of a +linear analysis. This may to many of you feel rather technical and +perhaps not that intuitive. The question is whether we can develop a +more intuitive way of understanding what Ridge and Lasso express. +
+ +Before we proceed let us perform a Ridge, Lasso and OLS analysis of a polynomial fit.
+ +We will play around with a study of the values for the optimal +parameters \( \boldsymbol{\beta} \) using OLS, Ridge and Lasso regression. For +OLS, you will notice as function of the noise and polynomial degree, +that the parameters \( \beta \) will fluctuate from order to order in the +polynomial fit and that for larger and larger polynomial degrees of freedom, the parameters will tend to increase in value for OLS. +
+ +For Ridge and Lasso regression, the higher order parameters will typically be reduced, providing thereby less fluctuations from one order to another one.
+ + + +import numpy as np
+import matplotlib.pyplot as plt
+from sklearn.model_selection import train_test_split
+from sklearn import linear_model
+
+def R2(y_data, y_model):
+ return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2)
+def MSE(y_data,y_model):
+ n = np.size(y_model)
+ return np.sum((y_data-y_model)**2)/n
+
+# Make data set.
+n = 10000
+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):
+ X[:,polydegree] = x**(polydegree)
+
+# 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)
+
+# matrix inversion to find beta
+OLSbeta = np.linalg.pinv(X_train.T @ X_train) @ X_train.T @ y_train
+print(OLSbeta)
+ypredictOLS = X_test @ OLSbeta
+print("Test MSE OLS")
+print(MSE(y_test,ypredictOLS))
+# Repeat now for Lasso and Ridge regression and various values of the regularization parameter using Scikit-Learn
+# Decide which values of lambda to use
+nlambdas = 4
+MSERidgePredict = np.zeros(nlambdas)
+MSELassoPredict = np.zeros(nlambdas)
+lambdas = np.logspace(-3, 1, nlambdas)
+for i in range(nlambdas):
+ lmb = lambdas[i]
+ # Make the fit using Ridge and Lasso
+ RegRidge = linear_model.Ridge(lmb,fit_intercept=False)
+ RegRidge.fit(X_train,y_train)
+ RegLasso = linear_model.Lasso(lmb,fit_intercept=False)
+ RegLasso.fit(X_train,y_train)
+ # and then make the prediction
+ ypredictRidge = RegRidge.predict(X_test)
+ ypredictLasso = RegLasso.predict(X_test)
+ # Compute the MSE and print it
+ MSERidgePredict[i] = MSE(y_test,ypredictRidge)
+ MSELassoPredict[i] = MSE(y_test,ypredictLasso)
+ print(lmb,RegRidge.coef_)
+ print(lmb,RegLasso.coef_)
+# Now plot the results
+plt.figure()
+plt.plot(np.log10(lambdas), MSERidgePredict, 'b', label = 'MSE Ridge Test')
+plt.plot(np.log10(lambdas), MSELassoPredict, 'r', label = 'MSE Lasso Test')
+plt.xlabel('log10(lambda)')
+plt.ylabel('MSE')
+plt.legend()
+plt.show()
+
+How can we understand this?
+Hitherto we have discussed Ridge and Lasso regression in terms of a -linear analysis. This may to many of you feel rather technical and -perhaps not that intuitive. The question is whether we can develop a -more intuitive way of understanding what Ridge and Lasso express. -
- -Before we proceed let us perform a Ridge, Lasso and OLS analysis of a polynomial fit.
- -We will play around with a study of the values for the optimal -parameters \( \boldsymbol{\beta} \) using OLS, Ridge and Lasso regression. For -OLS, you will notice as function of the noise and polynomial degree, -that the parameters \( \beta \) will fluctuate from order to order in the -polynomial fit and that for larger and larger polynomial degrees of freedom, the parameters will tend to increase in value for OLS. -
- -For Ridge and Lasso regression, the higher order parameters will typically be reduced, providing thereby less fluctuations from one order to another one.
- - - -import numpy as np
-import matplotlib.pyplot as plt
-from sklearn.model_selection import train_test_split
-from sklearn import linear_model
-
-def R2(y_data, y_model):
- return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2)
-def MSE(y_data,y_model):
- n = np.size(y_model)
- return np.sum((y_data-y_model)**2)/n
-
-# Make data set.
-n = 10000
-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)
-
-# matrix inversion to find beta
-OLSbeta = np.linalg.pinv(X_train.T @ X_train) @ X_train.T @ y_train
-print(OLSbeta)
-ypredictOLS = X_test @ OLSbeta
-print("Test MSE OLS")
-print(MSE(y_test,ypredictOLS))
-# Repeat now for Lasso and Ridge regression and various values of the regularization parameter using Scikit-Learn
-# Decide which values of lambda to use
-nlambdas = 4
-MSERidgePredict = np.zeros(nlambdas)
-MSELassoPredict = np.zeros(nlambdas)
-lambdas = np.logspace(-3, 1, nlambdas)
-for i in range(nlambdas):
- lmb = lambdas[i]
- # Make the fit using Ridge and Lasso
- RegRidge = linear_model.Ridge(lmb,fit_intercept=False)
- RegRidge.fit(X_train,y_train)
- RegLasso = linear_model.Lasso(lmb,fit_intercept=False)
- RegLasso.fit(X_train,y_train)
- # and then make the prediction
- ypredictRidge = RegRidge.predict(X_test)
- ypredictLasso = RegLasso.predict(X_test)
- # Compute the MSE and print it
- MSERidgePredict[i] = MSE(y_test,ypredictRidge)
- MSELassoPredict[i] = MSE(y_test,ypredictLasso)
- print(lmb,RegRidge.coef_)
- print(lmb,RegLasso.coef_)
-# Now plot the results
-plt.figure()
-plt.plot(np.log10(lambdas), MSERidgePredict, 'b', label = 'MSE Ridge Test')
-plt.plot(np.log10(lambdas), MSELassoPredict, 'r', label = 'MSE Lasso Test')
-plt.xlabel('log10(lambda)')
-plt.ylabel('MSE')
-plt.legend()
-plt.show()
-
-How can we understand this?
- -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)
diff --git a/doc/pub/week36/html/week36.html b/doc/pub/week36/html/week36.html index 79997f0cf..626ad64eb 100644 --- a/doc/pub/week36/html/week36.html +++ b/doc/pub/week36/html/week36.html @@ -220,6 +220,10 @@ div.toc p,a { 2, None, 'important-technicalities-more-on-rescaling-data'), + ('Test Function for what happens with OLS, Ridge and Lasso', + 2, + None, + 'test-function-for-what-happens-with-ols-ridge-and-lasso'), ('Linking the regression analysis with a statistical ' 'interpretation', 2, @@ -267,11 +271,6 @@ div.toc p,a { 2, None, 'bayes-theorem-and-ridge-and-lasso-regression'), - ('Test Function for what happens with OLS, Ridge and Lasso', - 2, - None, - 'test-function-for-what-happens-with-ols-ridge-and-lasso'), - ("Invoking Bayes' theorem", 2, None, 'invoking-bayes-theorem'), ('Ridge and Bayes', 2, None, 'ridge-and-bayes'), ('Lasso and Bayes', 2, None, 'lasso-and-bayes')]} end of tocinfo --> @@ -329,7 +328,7 @@ MathJax.Hub.Config({Hitherto we have discussed Ridge and Lasso regression in terms of a +linear analysis. This may to many of you feel rather technical and +perhaps not that intuitive. The question is whether we can develop a +more intuitive way of understanding what Ridge and Lasso express. +
+ +Before we proceed let us perform a Ridge, Lasso and OLS analysis of a polynomial fit.
+ +We will play around with a study of the values for the optimal +parameters \( \boldsymbol{\beta} \) using OLS, Ridge and Lasso regression. For +OLS, you will notice as function of the noise and polynomial degree, +that the parameters \( \beta \) will fluctuate from order to order in the +polynomial fit and that for larger and larger polynomial degrees of freedom, the parameters will tend to increase in value for OLS. +
+ +For Ridge and Lasso regression, the higher order parameters will typically be reduced, providing thereby less fluctuations from one order to another one.
+ + + +import numpy as np
+import matplotlib.pyplot as plt
+from sklearn.model_selection import train_test_split
+from sklearn import linear_model
+
+def R2(y_data, y_model):
+ return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2)
+def MSE(y_data,y_model):
+ n = np.size(y_model)
+ return np.sum((y_data-y_model)**2)/n
+
+# Make data set.
+n = 10000
+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):
+ X[:,polydegree] = x**(polydegree)
+
+# 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)
+
+# matrix inversion to find beta
+OLSbeta = np.linalg.pinv(X_train.T @ X_train) @ X_train.T @ y_train
+print(OLSbeta)
+ypredictOLS = X_test @ OLSbeta
+print("Test MSE OLS")
+print(MSE(y_test,ypredictOLS))
+# Repeat now for Lasso and Ridge regression and various values of the regularization parameter using Scikit-Learn
+# Decide which values of lambda to use
+nlambdas = 4
+MSERidgePredict = np.zeros(nlambdas)
+MSELassoPredict = np.zeros(nlambdas)
+lambdas = np.logspace(-3, 1, nlambdas)
+for i in range(nlambdas):
+ lmb = lambdas[i]
+ # Make the fit using Ridge and Lasso
+ RegRidge = linear_model.Ridge(lmb,fit_intercept=False)
+ RegRidge.fit(X_train,y_train)
+ RegLasso = linear_model.Lasso(lmb,fit_intercept=False)
+ RegLasso.fit(X_train,y_train)
+ # and then make the prediction
+ ypredictRidge = RegRidge.predict(X_test)
+ ypredictLasso = RegLasso.predict(X_test)
+ # Compute the MSE and print it
+ MSERidgePredict[i] = MSE(y_test,ypredictRidge)
+ MSELassoPredict[i] = MSE(y_test,ypredictLasso)
+ print(lmb,RegRidge.coef_)
+ print(lmb,RegLasso.coef_)
+# Now plot the results
+plt.figure()
+plt.plot(np.log10(lambdas), MSERidgePredict, 'b', label = 'MSE Ridge Test')
+plt.plot(np.log10(lambdas), MSELassoPredict, 'r', label = 'MSE Lasso Test')
+plt.xlabel('log10(lambda)')
+plt.ylabel('MSE')
+plt.legend()
+plt.show()
+
+How can we understand this?
+Hitherto we have discussed Ridge and Lasso regression in terms of a -linear analysis. This may to many of you feel rather technical and -perhaps not that intuitive. The question is whether we can develop a -more intuitive way of understanding what Ridge and Lasso express. -
- -Before we proceed let us perform a Ridge, Lasso and OLS analysis of a polynomial fit.
- -We will play around with a study of the values for the optimal -parameters \( \boldsymbol{\beta} \) using OLS, Ridge and Lasso regression. For -OLS, you will notice as function of the noise and polynomial degree, -that the parameters \( \beta \) will fluctuate from order to order in the -polynomial fit and that for larger and larger polynomial degrees of freedom, the parameters will tend to increase in value for OLS. -
- -For Ridge and Lasso regression, the higher order parameters will typically be reduced, providing thereby less fluctuations from one order to another one.
- - - -import numpy as np
-import matplotlib.pyplot as plt
-from sklearn.model_selection import train_test_split
-from sklearn import linear_model
-
-def R2(y_data, y_model):
- return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2)
-def MSE(y_data,y_model):
- n = np.size(y_model)
- return np.sum((y_data-y_model)**2)/n
-
-# Make data set.
-n = 10000
-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)
-
-# matrix inversion to find beta
-OLSbeta = np.linalg.pinv(X_train.T @ X_train) @ X_train.T @ y_train
-print(OLSbeta)
-ypredictOLS = X_test @ OLSbeta
-print("Test MSE OLS")
-print(MSE(y_test,ypredictOLS))
-# Repeat now for Lasso and Ridge regression and various values of the regularization parameter using Scikit-Learn
-# Decide which values of lambda to use
-nlambdas = 4
-MSERidgePredict = np.zeros(nlambdas)
-MSELassoPredict = np.zeros(nlambdas)
-lambdas = np.logspace(-3, 1, nlambdas)
-for i in range(nlambdas):
- lmb = lambdas[i]
- # Make the fit using Ridge and Lasso
- RegRidge = linear_model.Ridge(lmb,fit_intercept=False)
- RegRidge.fit(X_train,y_train)
- RegLasso = linear_model.Lasso(lmb,fit_intercept=False)
- RegLasso.fit(X_train,y_train)
- # and then make the prediction
- ypredictRidge = RegRidge.predict(X_test)
- ypredictLasso = RegLasso.predict(X_test)
- # Compute the MSE and print it
- MSERidgePredict[i] = MSE(y_test,ypredictRidge)
- MSELassoPredict[i] = MSE(y_test,ypredictLasso)
- print(lmb,RegRidge.coef_)
- print(lmb,RegLasso.coef_)
-# Now plot the results
-plt.figure()
-plt.plot(np.log10(lambdas), MSERidgePredict, 'b', label = 'MSE Ridge Test')
-plt.plot(np.log10(lambdas), MSELassoPredict, 'r', label = 'MSE Lasso Test')
-plt.xlabel('log10(lambda)')
-plt.ylabel('MSE')
-plt.legend()
-plt.show()
-
-How can we understand this?
- -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)
diff --git a/doc/pub/week36/ipynb/ipynb-week36-src.tar.gz b/doc/pub/week36/ipynb/ipynb-week36-src.tar.gz index db754fea3..69318c705 100644 Binary files a/doc/pub/week36/ipynb/ipynb-week36-src.tar.gz and b/doc/pub/week36/ipynb/ipynb-week36-src.tar.gz differ diff --git a/doc/pub/week36/ipynb/week36.ipynb b/doc/pub/week36/ipynb/week36.ipynb index efb2e069f..d17233c27 100644 --- a/doc/pub/week36/ipynb/week36.ipynb +++ b/doc/pub/week36/ipynb/week36.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "markdown", - "id": "e06b0e94", + "id": "ff34ab4f", "metadata": { "editable": true }, @@ -14,7 +14,7 @@ }, { "cell_type": "markdown", - "id": "145c9d1c", + "id": "a88eeeaf", "metadata": { "editable": true }, @@ -27,7 +27,7 @@ }, { "cell_type": "markdown", - "id": "035eed6a", + "id": "49b54282", "metadata": { "editable": true }, @@ -48,14 +48,14 @@ "\n", " * Linear Regression and links with Statistics\n", "\n", - " * Recommended Reading: Goodfellow et al chapter 3 on probability theory, see URL:\"\"\n", + " * [Recommended Reading: Goodfellow et al chapter 3 on probability theory](https://www.deeplearningbook.org/)\n", "\n", " * See also Murphy, sections 2.4 (Gaussian distributions) and 3.2 (Bayesian Statistics, basis)" ] }, { "cell_type": "markdown", - "id": "1f3dd303", + "id": "cb61e293", "metadata": { "editable": true }, @@ -67,7 +67,7 @@ }, { "cell_type": "markdown", - "id": "c81143e1", + "id": "0b12f4e7", "metadata": { "editable": true }, @@ -79,7 +79,7 @@ }, { "cell_type": "markdown", - "id": "622d1ab3", + "id": "b77ff400", "metadata": { "editable": true }, @@ -91,7 +91,7 @@ }, { "cell_type": "markdown", - "id": "7746de71", + "id": "94203cb0", "metadata": { "editable": true }, @@ -101,7 +101,7 @@ }, { "cell_type": "markdown", - "id": "d2cb17f0", + "id": "9adad1a7", "metadata": { "editable": true }, @@ -113,7 +113,7 @@ }, { "cell_type": "markdown", - "id": "6d3a1838", + "id": "8d44ef1e", "metadata": { "editable": true }, @@ -123,7 +123,7 @@ }, { "cell_type": "markdown", - "id": "c20c462c", + "id": "2b65e6c1", "metadata": { "editable": true }, @@ -140,7 +140,7 @@ }, { "cell_type": "markdown", - "id": "58a5ee08", + "id": "4bec5357", "metadata": { "editable": true }, @@ -150,7 +150,7 @@ }, { "cell_type": "markdown", - "id": "f9ca7d34", + "id": "4240b8b5", "metadata": { "editable": true }, @@ -162,7 +162,7 @@ }, { "cell_type": "markdown", - "id": "02ef6562", + "id": "bf4a51f1", "metadata": { "editable": true }, @@ -172,7 +172,7 @@ }, { "cell_type": "markdown", - "id": "d401a468", + "id": "fbec17ea", "metadata": { "editable": true }, @@ -184,7 +184,7 @@ }, { "cell_type": "markdown", - "id": "d6dc5c03", + "id": "722612e1", "metadata": { "editable": true }, @@ -198,7 +198,7 @@ }, { "cell_type": "markdown", - "id": "7e399fe1", + "id": "916f14e8", "metadata": { "editable": true }, @@ -210,7 +210,7 @@ }, { "cell_type": "markdown", - "id": "7157ab81", + "id": "34de5c7a", "metadata": { "editable": true }, @@ -232,7 +232,7 @@ }, { "cell_type": "markdown", - "id": "16086935", + "id": "4db4d16d", "metadata": { "editable": true }, @@ -244,7 +244,7 @@ }, { "cell_type": "markdown", - "id": "fb93fadb", + "id": "4fcefcd0", "metadata": { "editable": true }, @@ -259,7 +259,7 @@ }, { "cell_type": "markdown", - "id": "c587b168", + "id": "3f3197db", "metadata": { "editable": true }, @@ -271,7 +271,7 @@ }, { "cell_type": "markdown", - "id": "ef8c785a", + "id": "0a6eaa09", "metadata": { "editable": true }, @@ -283,7 +283,7 @@ }, { "cell_type": "markdown", - "id": "4f2998ee", + "id": "c7d43a29", "metadata": { "editable": true }, @@ -293,7 +293,7 @@ }, { "cell_type": "markdown", - "id": "451ce7b4", + "id": "33e14cc4", "metadata": { "editable": true }, @@ -305,7 +305,7 @@ }, { "cell_type": "markdown", - "id": "172c78b9", + "id": "8ad076ff", "metadata": { "editable": true }, @@ -315,7 +315,7 @@ }, { "cell_type": "markdown", - "id": "38fef505", + "id": "b2b1aeac", "metadata": { "editable": true }, @@ -327,7 +327,7 @@ }, { "cell_type": "markdown", - "id": "9a010e65", + "id": "8deb30d4", "metadata": { "editable": true }, @@ -337,7 +337,7 @@ }, { "cell_type": "markdown", - "id": "97bd739b", + "id": "f1c7bce6", "metadata": { "editable": true }, @@ -349,7 +349,7 @@ }, { "cell_type": "markdown", - "id": "55a422e9", + "id": "d90876f4", "metadata": { "editable": true }, @@ -366,7 +366,7 @@ }, { "cell_type": "markdown", - "id": "5c781ed1", + "id": "85930202", "metadata": { "editable": true }, @@ -380,7 +380,7 @@ { "cell_type": "code", "execution_count": 1, - "id": "4efabee8", + "id": "540328fb", "metadata": { "collapsed": false, "editable": true @@ -394,7 +394,7 @@ }, { "cell_type": "markdown", - "id": "4116bb0c", + "id": "58c2ee4a", "metadata": { "editable": true }, @@ -405,7 +405,7 @@ { "cell_type": "code", "execution_count": 2, - "id": "f2b03369", + "id": "ec5c2fe8", "metadata": { "collapsed": false, "editable": true @@ -445,7 +445,7 @@ }, { "cell_type": "markdown", - "id": "7d91bf0b", + "id": "15898f2d", "metadata": { "editable": true }, @@ -465,7 +465,7 @@ }, { "cell_type": "markdown", - "id": "2d3dab62", + "id": "82a70c91", "metadata": { "editable": true }, @@ -477,7 +477,7 @@ }, { "cell_type": "markdown", - "id": "35031c2d", + "id": "c60476a5", "metadata": { "editable": true }, @@ -488,7 +488,7 @@ { "cell_type": "code", "execution_count": 3, - "id": "4b48f5b7", + "id": "114372da", "metadata": { "collapsed": false, "editable": true @@ -522,7 +522,7 @@ }, { "cell_type": "markdown", - "id": "e5671fb7", + "id": "f08b1538", "metadata": { "editable": true }, @@ -532,7 +532,7 @@ }, { "cell_type": "markdown", - "id": "5d94606d", + "id": "515bf446", "metadata": { "editable": true }, @@ -545,7 +545,7 @@ }, { "cell_type": "markdown", - "id": "b7b42f2a", + "id": "8e63f91c", "metadata": { "editable": true }, @@ -557,7 +557,7 @@ }, { "cell_type": "markdown", - "id": "5cf829d1", + "id": "37cc5c42", "metadata": { "editable": true }, @@ -567,7 +567,7 @@ }, { "cell_type": "markdown", - "id": "cf7e3359", + "id": "3cb8a97a", "metadata": { "editable": true }, @@ -580,7 +580,7 @@ }, { "cell_type": "markdown", - "id": "6ca412cc", + "id": "d00c65cc", "metadata": { "editable": true }, @@ -590,7 +590,7 @@ }, { "cell_type": "markdown", - "id": "a497f398", + "id": "1c4684ff", "metadata": { "editable": true }, @@ -602,7 +602,7 @@ }, { "cell_type": "markdown", - "id": "95a2876e", + "id": "899a6755", "metadata": { "editable": true }, @@ -617,7 +617,7 @@ }, { "cell_type": "markdown", - "id": "dca0e875", + "id": "35858a3a", "metadata": { "editable": true }, @@ -630,7 +630,7 @@ }, { "cell_type": "markdown", - "id": "c09ec13f", + "id": "fcbba7d8", "metadata": { "editable": true }, @@ -644,7 +644,7 @@ }, { "cell_type": "markdown", - "id": "c8043e52", + "id": "26cb747b", "metadata": { "editable": true }, @@ -656,7 +656,7 @@ }, { "cell_type": "markdown", - "id": "4b3098ad", + "id": "ae7c22be", "metadata": { "editable": true }, @@ -666,7 +666,7 @@ }, { "cell_type": "markdown", - "id": "2a4fd104", + "id": "0ef7ea7a", "metadata": { "editable": true }, @@ -679,7 +679,7 @@ }, { "cell_type": "markdown", - "id": "38f204d5", + "id": "285f0e32", "metadata": { "editable": true }, @@ -691,7 +691,7 @@ }, { "cell_type": "markdown", - "id": "09b2359a", + "id": "48dc1686", "metadata": { "editable": true }, @@ -703,7 +703,7 @@ }, { "cell_type": "markdown", - "id": "029e4332", + "id": "94bea820", "metadata": { "editable": true }, @@ -715,7 +715,7 @@ }, { "cell_type": "markdown", - "id": "b95acbbd", + "id": "79dbe6d8", "metadata": { "editable": true }, @@ -727,7 +727,7 @@ }, { "cell_type": "markdown", - "id": "59312592", + "id": "299d4456", "metadata": { "editable": true }, @@ -741,7 +741,7 @@ }, { "cell_type": "markdown", - "id": "5b886126", + "id": "28f1259b", "metadata": { "editable": true }, @@ -753,7 +753,7 @@ }, { "cell_type": "markdown", - "id": "2e89e4e5", + "id": "ef07fce1", "metadata": { "editable": true }, @@ -763,7 +763,7 @@ }, { "cell_type": "markdown", - "id": "02b4539c", + "id": "567e07ca", "metadata": { "editable": true }, @@ -775,7 +775,7 @@ }, { "cell_type": "markdown", - "id": "c1090a64", + "id": "ab200db0", "metadata": { "editable": true }, @@ -785,7 +785,7 @@ }, { "cell_type": "markdown", - "id": "db24e71d", + "id": "e2d8ca37", "metadata": { "editable": true }, @@ -797,7 +797,7 @@ }, { "cell_type": "markdown", - "id": "15f9d2c6", + "id": "604b3117", "metadata": { "editable": true }, @@ -809,7 +809,7 @@ }, { "cell_type": "markdown", - "id": "530ace6f", + "id": "da04983b", "metadata": { "editable": true }, @@ -819,7 +819,7 @@ }, { "cell_type": "markdown", - "id": "25ef4bcc", + "id": "2bebf877", "metadata": { "editable": true }, @@ -830,7 +830,7 @@ }, { "cell_type": "markdown", - "id": "3068fc1c", + "id": "6152d137", "metadata": { "editable": true }, @@ -842,7 +842,7 @@ }, { "cell_type": "markdown", - "id": "35b818a4", + "id": "eccd4760", "metadata": { "editable": true }, @@ -858,7 +858,7 @@ }, { "cell_type": "markdown", - "id": "42a61f8c", + "id": "73aef7e2", "metadata": { "editable": true }, @@ -871,7 +871,7 @@ }, { "cell_type": "markdown", - "id": "bf7c24e1", + "id": "37b86951", "metadata": { "editable": true }, @@ -883,7 +883,7 @@ }, { "cell_type": "markdown", - "id": "488834ec", + "id": "9c638256", "metadata": { "editable": true }, @@ -893,7 +893,7 @@ }, { "cell_type": "markdown", - "id": "b094f27c", + "id": "a67ccade", "metadata": { "editable": true }, @@ -905,7 +905,7 @@ }, { "cell_type": "markdown", - "id": "9a2bf498", + "id": "d60a0acb", "metadata": { "editable": true }, @@ -915,7 +915,7 @@ }, { "cell_type": "markdown", - "id": "69e176f5", + "id": "7d16dcf3", "metadata": { "editable": true }, @@ -927,7 +927,7 @@ }, { "cell_type": "markdown", - "id": "d550d346", + "id": "3b162883", "metadata": { "editable": true }, @@ -939,7 +939,7 @@ }, { "cell_type": "markdown", - "id": "b4f33917", + "id": "fa81d3cf", "metadata": { "editable": true }, @@ -955,7 +955,7 @@ }, { "cell_type": "markdown", - "id": "3301a74c", + "id": "c35848c1", "metadata": { "editable": true }, @@ -967,7 +967,7 @@ }, { "cell_type": "markdown", - "id": "7b9baeb9", + "id": "f75f3dea", "metadata": { "editable": true }, @@ -979,7 +979,7 @@ }, { "cell_type": "markdown", - "id": "a5b1552b", + "id": "2b407fea", "metadata": { "editable": true }, @@ -989,7 +989,7 @@ }, { "cell_type": "markdown", - "id": "4b1ff318", + "id": "5fad7e1b", "metadata": { "editable": true }, @@ -1001,7 +1001,7 @@ }, { "cell_type": "markdown", - "id": "a3fc6e1a", + "id": "de011fe9", "metadata": { "editable": true }, @@ -1011,7 +1011,7 @@ }, { "cell_type": "markdown", - "id": "128ed808", + "id": "ad61d5a1", "metadata": { "editable": true }, @@ -1023,7 +1023,7 @@ }, { "cell_type": "markdown", - "id": "0c72c2fe", + "id": "eff9b775", "metadata": { "editable": true }, @@ -1040,7 +1040,7 @@ }, { "cell_type": "markdown", - "id": "a23c0f35", + "id": "6c4f7a24", "metadata": { "editable": true }, @@ -1052,7 +1052,7 @@ }, { "cell_type": "markdown", - "id": "b19d4f06", + "id": "f4deb050", "metadata": { "editable": true }, @@ -1064,7 +1064,7 @@ }, { "cell_type": "markdown", - "id": "12422bbc", + "id": "9d6cfa3e", "metadata": { "editable": true }, @@ -1074,7 +1074,7 @@ }, { "cell_type": "markdown", - "id": "5fd9439c", + "id": "c2303440", "metadata": { "editable": true }, @@ -1086,7 +1086,7 @@ }, { "cell_type": "markdown", - "id": "a9930104", + "id": "7bf0b1fc", "metadata": { "editable": true }, @@ -1096,7 +1096,7 @@ }, { "cell_type": "markdown", - "id": "0d210eff", + "id": "25fe5980", "metadata": { "editable": true }, @@ -1108,7 +1108,7 @@ }, { "cell_type": "markdown", - "id": "4aefc549", + "id": "38a46583", "metadata": { "editable": true }, @@ -1118,7 +1118,7 @@ }, { "cell_type": "markdown", - "id": "dfed14eb", + "id": "45666dad", "metadata": { "editable": true }, @@ -1130,7 +1130,7 @@ }, { "cell_type": "markdown", - "id": "41b7e154", + "id": "d72cd292", "metadata": { "editable": true }, @@ -1140,7 +1140,7 @@ }, { "cell_type": "markdown", - "id": "df64b9cc", + "id": "1cdd45ed", "metadata": { "editable": true }, @@ -1155,7 +1155,7 @@ }, { "cell_type": "markdown", - "id": "54d3cef5", + "id": "bf327b4f", "metadata": { "editable": true }, @@ -1167,7 +1167,7 @@ }, { "cell_type": "markdown", - "id": "6c297d55", + "id": "235e1172", "metadata": { "editable": true }, @@ -1177,7 +1177,7 @@ }, { "cell_type": "markdown", - "id": "2d0ee0d8", + "id": "c0df9f8a", "metadata": { "editable": true }, @@ -1189,7 +1189,7 @@ }, { "cell_type": "markdown", - "id": "81d831f1", + "id": "88642f9a", "metadata": { "editable": true }, @@ -1201,7 +1201,7 @@ }, { "cell_type": "markdown", - "id": "663ad00e", + "id": "c48f87e2", "metadata": { "editable": true }, @@ -1213,7 +1213,7 @@ }, { "cell_type": "markdown", - "id": "748acfc4", + "id": "ae22178d", "metadata": { "editable": true }, @@ -1223,7 +1223,7 @@ }, { "cell_type": "markdown", - "id": "6c7eca03", + "id": "3c8c1324", "metadata": { "editable": true }, @@ -1235,7 +1235,7 @@ }, { "cell_type": "markdown", - "id": "ad13171b", + "id": "f9f58a15", "metadata": { "editable": true }, @@ -1247,7 +1247,7 @@ }, { "cell_type": "markdown", - "id": "538976df", + "id": "865aa159", "metadata": { "editable": true }, @@ -1259,7 +1259,7 @@ }, { "cell_type": "markdown", - "id": "dc880394", + "id": "9e6565ec", "metadata": { "editable": true }, @@ -1269,7 +1269,7 @@ }, { "cell_type": "markdown", - "id": "322aa407", + "id": "8b2950b4", "metadata": { "editable": true }, @@ -1281,7 +1281,7 @@ }, { "cell_type": "markdown", - "id": "96430ed6", + "id": "23a7c1a6", "metadata": { "editable": true }, @@ -1291,7 +1291,7 @@ }, { "cell_type": "markdown", - "id": "f3d0d853", + "id": "4a733817", "metadata": { "editable": true }, @@ -1305,7 +1305,7 @@ }, { "cell_type": "markdown", - "id": "827a8646", + "id": "dc517d14", "metadata": { "editable": true }, @@ -1315,7 +1315,7 @@ }, { "cell_type": "markdown", - "id": "f33aa235", + "id": "ed79ab9a", "metadata": { "editable": true }, @@ -1327,7 +1327,7 @@ }, { "cell_type": "markdown", - "id": "a2ba7f8e", + "id": "116b5106", "metadata": { "editable": true }, @@ -1339,7 +1339,7 @@ }, { "cell_type": "markdown", - "id": "dfcc2679", + "id": "aaa68ab6", "metadata": { "editable": true }, @@ -1349,7 +1349,7 @@ }, { "cell_type": "markdown", - "id": "c373a2ec", + "id": "17fa6c43", "metadata": { "editable": true }, @@ -1361,7 +1361,7 @@ }, { "cell_type": "markdown", - "id": "c2ee229a", + "id": "01264686", "metadata": { "editable": true }, @@ -1371,7 +1371,7 @@ }, { "cell_type": "markdown", - "id": "6ad191aa", + "id": "1b53fa9e", "metadata": { "editable": true }, @@ -1383,7 +1383,7 @@ }, { "cell_type": "markdown", - "id": "461735f1", + "id": "0a7cdc65", "metadata": { "editable": true }, @@ -1395,7 +1395,7 @@ }, { "cell_type": "markdown", - "id": "071db765", + "id": "6fac325b", "metadata": { "editable": true }, @@ -1405,7 +1405,7 @@ }, { "cell_type": "markdown", - "id": "002ab2d0", + "id": "2befb748", "metadata": { "editable": true }, @@ -1417,7 +1417,7 @@ }, { "cell_type": "markdown", - "id": "b54920c3", + "id": "c67fea47", "metadata": { "editable": true }, @@ -1427,7 +1427,7 @@ }, { "cell_type": "markdown", - "id": "5bb8e57d", + "id": "65e47e22", "metadata": { "editable": true }, @@ -1439,7 +1439,7 @@ }, { "cell_type": "markdown", - "id": "f7d0b833", + "id": "221a5457", "metadata": { "editable": true }, @@ -1451,7 +1451,7 @@ }, { "cell_type": "markdown", - "id": "d9a21c66", + "id": "5517307a", "metadata": { "editable": true }, @@ -1461,7 +1461,7 @@ }, { "cell_type": "markdown", - "id": "d35bc697", + "id": "b20940ce", "metadata": { "editable": true }, @@ -1473,7 +1473,7 @@ }, { "cell_type": "markdown", - "id": "1a5578d7", + "id": "4110bef8", "metadata": { "editable": true }, @@ -1486,7 +1486,7 @@ }, { "cell_type": "markdown", - "id": "a2cf0dc8", + "id": "0dad3074", "metadata": { "editable": true }, @@ -1498,7 +1498,7 @@ }, { "cell_type": "markdown", - "id": "6740771b", + "id": "350a0faa", "metadata": { "editable": true }, @@ -1510,7 +1510,7 @@ }, { "cell_type": "markdown", - "id": "9147be96", + "id": "24cd86b7", "metadata": { "editable": true }, @@ -1522,7 +1522,7 @@ }, { "cell_type": "markdown", - "id": "e32c6fc0", + "id": "17c643f0", "metadata": { "editable": true }, @@ -1532,7 +1532,7 @@ }, { "cell_type": "markdown", - "id": "71b8da53", + "id": "4cd2f888", "metadata": { "editable": true }, @@ -1544,7 +1544,7 @@ }, { "cell_type": "markdown", - "id": "ed90c422", + "id": "13ce9a3c", "metadata": { "editable": true }, @@ -1554,7 +1554,7 @@ }, { "cell_type": "markdown", - "id": "f803b922", + "id": "9c52a280", "metadata": { "editable": true }, @@ -1566,7 +1566,7 @@ }, { "cell_type": "markdown", - "id": "7ed0142e", + "id": "4aa82b09", "metadata": { "editable": true }, @@ -1576,7 +1576,7 @@ }, { "cell_type": "markdown", - "id": "38c06ddd", + "id": "7369b16e", "metadata": { "editable": true }, @@ -1588,7 +1588,7 @@ }, { "cell_type": "markdown", - "id": "c6a1c64a", + "id": "1fc52972", "metadata": { "editable": true }, @@ -1598,7 +1598,7 @@ }, { "cell_type": "markdown", - "id": "427819e1", + "id": "40baac3a", "metadata": { "editable": true }, @@ -1611,7 +1611,7 @@ }, { "cell_type": "markdown", - "id": "8c59625a", + "id": "c038e0d3", "metadata": { "editable": true }, @@ -1623,7 +1623,7 @@ }, { "cell_type": "markdown", - "id": "571cb692", + "id": "4ce575fa", "metadata": { "editable": true }, @@ -1635,7 +1635,7 @@ }, { "cell_type": "markdown", - "id": "986c666e", + "id": "1e5a661f", "metadata": { "editable": true }, @@ -1645,7 +1645,7 @@ }, { "cell_type": "markdown", - "id": "21174b3b", + "id": "0a170322", "metadata": { "editable": true }, @@ -1657,7 +1657,7 @@ }, { "cell_type": "markdown", - "id": "2db9fd19", + "id": "2443ec27", "metadata": { "editable": true }, @@ -1674,7 +1674,7 @@ }, { "cell_type": "markdown", - "id": "fc01dbbd", + "id": "64723ae2", "metadata": { "editable": true }, @@ -1686,7 +1686,7 @@ }, { "cell_type": "markdown", - "id": "a3645253", + "id": "d5723b13", "metadata": { "editable": true }, @@ -1698,7 +1698,7 @@ }, { "cell_type": "markdown", - "id": "0cfa4a22", + "id": "b64de210", "metadata": { "editable": true }, @@ -1708,7 +1708,7 @@ }, { "cell_type": "markdown", - "id": "562f4c4d", + "id": "207377d9", "metadata": { "editable": true }, @@ -1720,7 +1720,7 @@ }, { "cell_type": "markdown", - "id": "338354d4", + "id": "38658c17", "metadata": { "editable": true }, @@ -1730,7 +1730,7 @@ }, { "cell_type": "markdown", - "id": "8a32486c", + "id": "038520f6", "metadata": { "editable": true }, @@ -1742,7 +1742,7 @@ }, { "cell_type": "markdown", - "id": "f76d2153", + "id": "342f1a76", "metadata": { "editable": true }, @@ -1752,7 +1752,7 @@ }, { "cell_type": "markdown", - "id": "77d83d13", + "id": "6addac68", "metadata": { "editable": true }, @@ -1764,7 +1764,7 @@ }, { "cell_type": "markdown", - "id": "f8d54178", + "id": "78cc3bce", "metadata": { "editable": true }, @@ -1774,7 +1774,7 @@ }, { "cell_type": "markdown", - "id": "af88a9fc", + "id": "420ca9dd", "metadata": { "editable": true }, @@ -1789,7 +1789,7 @@ { "cell_type": "code", "execution_count": 4, - "id": "e039011c", + "id": "585576f1", "metadata": { "collapsed": false, "editable": true @@ -1851,7 +1851,7 @@ }, { "cell_type": "markdown", - "id": "621b2394", + "id": "9d887ba8", "metadata": { "editable": true }, @@ -1861,7 +1861,7 @@ }, { "cell_type": "markdown", - "id": "6ca32c73", + "id": "8c6667cd", "metadata": { "editable": true }, @@ -1872,7 +1872,7 @@ { "cell_type": "code", "execution_count": 5, - "id": "9e7fca6e", + "id": "469fcfa1", "metadata": { "collapsed": false, "editable": true @@ -1939,7 +1939,7 @@ }, { "cell_type": "markdown", - "id": "52dca90c", + "id": "62e61bf5", "metadata": { "editable": true }, @@ -1950,7 +1950,7 @@ { "cell_type": "code", "execution_count": 6, - "id": "386780bb", + "id": "b0302344", "metadata": { "collapsed": false, "editable": true @@ -2039,7 +2039,7 @@ }, { "cell_type": "markdown", - "id": "dd7db85c", + "id": "6e699b2f", "metadata": { "editable": true }, @@ -2049,7 +2049,7 @@ }, { "cell_type": "markdown", - "id": "e8a44276", + "id": "d39bb025", "metadata": { "editable": true }, @@ -2106,7 +2106,7 @@ { "cell_type": "code", "execution_count": 7, - "id": "ced43b74", + "id": "6f772122", "metadata": { "collapsed": false, "editable": true @@ -2133,7 +2133,7 @@ }, { "cell_type": "markdown", - "id": "78931069", + "id": "38ff7f71", "metadata": { "editable": true }, @@ -2147,7 +2147,7 @@ }, { "cell_type": "markdown", - "id": "87469805", + "id": "7dba297c", "metadata": { "editable": true }, @@ -2159,7 +2159,7 @@ }, { "cell_type": "markdown", - "id": "d2f82e8f", + "id": "772f2dc0", "metadata": { "editable": true }, @@ -2176,7 +2176,7 @@ }, { "cell_type": "markdown", - "id": "17c43864", + "id": "1908cc29", "metadata": { "editable": true }, @@ -2188,7 +2188,7 @@ }, { "cell_type": "markdown", - "id": "4d119d27", + "id": "fb31721c", "metadata": { "editable": true }, @@ -2198,7 +2198,7 @@ }, { "cell_type": "markdown", - "id": "0775ec0e", + "id": "01d332e2", "metadata": { "editable": true }, @@ -2210,7 +2210,7 @@ }, { "cell_type": "markdown", - "id": "2d4ee560", + "id": "fd2cc608", "metadata": { "editable": true }, @@ -2220,7 +2220,7 @@ }, { "cell_type": "markdown", - "id": "d458b594", + "id": "abe97943", "metadata": { "editable": true }, @@ -2232,7 +2232,7 @@ }, { "cell_type": "markdown", - "id": "cf214030", + "id": "4b0c9c12", "metadata": { "editable": true }, @@ -2243,7 +2243,7 @@ }, { "cell_type": "markdown", - "id": "bdeb1a3a", + "id": "addd3d97", "metadata": { "editable": true }, @@ -2255,7 +2255,7 @@ }, { "cell_type": "markdown", - "id": "db41b9ba", + "id": "a96063de", "metadata": { "editable": true }, @@ -2265,7 +2265,7 @@ }, { "cell_type": "markdown", - "id": "d7c3855d", + "id": "61d186e9", "metadata": { "editable": true }, @@ -2277,7 +2277,7 @@ }, { "cell_type": "markdown", - "id": "b3f77327", + "id": "9b976247", "metadata": { "editable": true }, @@ -2287,7 +2287,7 @@ }, { "cell_type": "markdown", - "id": "ded3e3c5", + "id": "c935d0d7", "metadata": { "editable": true }, @@ -2299,7 +2299,7 @@ }, { "cell_type": "markdown", - "id": "2f96f46b", + "id": "9f4e09e3", "metadata": { "editable": true }, @@ -2309,7 +2309,7 @@ }, { "cell_type": "markdown", - "id": "67e05813", + "id": "9438c6ec", "metadata": { "editable": true }, @@ -2321,7 +2321,7 @@ }, { "cell_type": "markdown", - "id": "9efbc0ca", + "id": "f773450a", "metadata": { "editable": true }, @@ -2331,7 +2331,7 @@ }, { "cell_type": "markdown", - "id": "6d85e760", + "id": "1b530b85", "metadata": { "editable": true }, @@ -2343,7 +2343,7 @@ }, { "cell_type": "markdown", - "id": "f087fc8d", + "id": "4fa771d8", "metadata": { "editable": true }, @@ -2353,7 +2353,7 @@ }, { "cell_type": "markdown", - "id": "13efce6b", + "id": "09116fa4", "metadata": { "editable": true }, @@ -2365,7 +2365,7 @@ }, { "cell_type": "markdown", - "id": "e9b11189", + "id": "5011c21f", "metadata": { "editable": true }, @@ -2375,7 +2375,7 @@ }, { "cell_type": "markdown", - "id": "3c8c288e", + "id": "b1d51013", "metadata": { "editable": true }, @@ -2387,7 +2387,7 @@ }, { "cell_type": "markdown", - "id": "50a413a7", + "id": "f12877c0", "metadata": { "editable": true }, @@ -2397,7 +2397,7 @@ }, { "cell_type": "markdown", - "id": "9f98e4d5", + "id": "0b27f009", "metadata": { "editable": true }, @@ -2409,7 +2409,7 @@ }, { "cell_type": "markdown", - "id": "44803cbe", + "id": "dc8e6f9a", "metadata": { "editable": true }, @@ -2421,7 +2421,7 @@ }, { "cell_type": "markdown", - "id": "9d71963d", + "id": "c76d277c", "metadata": { "editable": true }, @@ -2433,7 +2433,7 @@ }, { "cell_type": "markdown", - "id": "ff6865ff", + "id": "3986fa93", "metadata": { "editable": true }, @@ -2443,7 +2443,7 @@ }, { "cell_type": "markdown", - "id": "03e7606a", + "id": "07298cf7", "metadata": { "editable": true }, @@ -2455,7 +2455,7 @@ }, { "cell_type": "markdown", - "id": "65288609", + "id": "84f96f80", "metadata": { "editable": true }, @@ -2468,7 +2468,7 @@ }, { "cell_type": "markdown", - "id": "8df4447a", + "id": "54ba61b3", "metadata": { "editable": true }, @@ -2480,7 +2480,7 @@ }, { "cell_type": "markdown", - "id": "696034eb", + "id": "e23b4f5e", "metadata": { "editable": true }, @@ -2494,7 +2494,7 @@ { "cell_type": "code", "execution_count": 8, - "id": "cc398bad", + "id": "3a6e4b8d", "metadata": { "collapsed": false, "editable": true @@ -2591,7 +2591,7 @@ }, { "cell_type": "markdown", - "id": "1122d69c", + "id": "d938305a", "metadata": { "editable": true }, @@ -2612,7 +2612,7 @@ }, { "cell_type": "markdown", - "id": "354c1f91", + "id": "9797b358", "metadata": { "editable": true }, @@ -2624,7 +2624,7 @@ }, { "cell_type": "markdown", - "id": "dfa71b80", + "id": "4962007b", "metadata": { "editable": true }, @@ -2634,7 +2634,7 @@ }, { "cell_type": "markdown", - "id": "f9a2bc3c", + "id": "84738217", "metadata": { "editable": true }, @@ -2646,7 +2646,7 @@ }, { "cell_type": "markdown", - "id": "186d5010", + "id": "340cdb6e", "metadata": { "editable": true }, @@ -2656,7 +2656,7 @@ }, { "cell_type": "markdown", - "id": "591e3c2f", + "id": "ea0d39ea", "metadata": { "editable": true }, @@ -2668,7 +2668,7 @@ }, { "cell_type": "markdown", - "id": "5f86c0da", + "id": "e18a34e9", "metadata": { "editable": true }, @@ -2686,7 +2686,7 @@ { "cell_type": "code", "execution_count": 9, - "id": "7ee195c4", + "id": "c7c1b283", "metadata": { "collapsed": false, "editable": true @@ -2762,7 +2762,7 @@ }, { "cell_type": "markdown", - "id": "21e60356", + "id": "26406e58", "metadata": { "editable": true }, @@ -2776,7 +2776,7 @@ { "cell_type": "code", "execution_count": 10, - "id": "34cc4689", + "id": "ce41d0b7", "metadata": { "collapsed": false, "editable": true @@ -2865,7 +2865,7 @@ }, { "cell_type": "markdown", - "id": "39f93f17", + "id": "3a37bf34", "metadata": { "editable": true }, @@ -2880,962 +2880,19 @@ }, { "cell_type": "markdown", - "id": "6b489532", + "id": "f33e21aa", "metadata": { "editable": true }, "source": [ - "## Linking the regression analysis with a statistical interpretation\n", - "\n", - "We will now couple the discussions of ordinary least squares, Ridge\n", - "and Lasso regression with a statistical interpretation, that is we\n", - "move from a linear algebra analysis to a statistical analysis. In\n", - "particular, we will focus on what the regularization terms can result\n", - "in. We will amongst other things show that the regularization\n", - "parameter can reduce considerably the variance of the parameters\n", - "$\\beta$.\n", - "\n", - "The\n", - "advantage of doing linear regression is that we actually end up with\n", - "analytical expressions for several statistical quantities. \n", - "Standard least squares and Ridge regression allow us to\n", - "derive quantities like the variance and other expectation values in a\n", - "rather straightforward way.\n", - "\n", - "It is assumed that $\\varepsilon_i\n", - "\\sim \\mathcal{N}(0, \\sigma^2)$ and the $\\varepsilon_{i}$ are\n", - "independent, i.e.:" - ] - }, - { - "cell_type": "markdown", - "id": "8dc63362", - "metadata": { - "editable": true - }, - "source": [ - "$$\n", - "\\begin{align*} \n", - "\\mbox{Cov}(\\varepsilon_{i_1},\n", - "\\varepsilon_{i_2}) & = \\left\\{ \\begin{array}{lcc} \\sigma^2 & \\mbox{if}\n", - "& i_1 = i_2, \\\\ 0 & \\mbox{if} & i_1 \\not= i_2. \\end{array} \\right.\n", - "\\end{align*}\n", - "$$" - ] - }, - { - "cell_type": "markdown", - "id": "8f7d73e3", - "metadata": { - "editable": true - }, - "source": [ - "The randomness of $\\varepsilon_i$ implies that\n", - "$\\mathbf{y}_i$ is also a random variable. In particular,\n", - "$\\mathbf{y}_i$ is normally distributed, because $\\varepsilon_i \\sim\n", - "\\mathcal{N}(0, \\sigma^2)$ and $\\mathbf{X}_{i,\\ast} \\, \\boldsymbol{\\beta}$ is a\n", - "non-random scalar. To specify the parameters of the distribution of\n", - "$\\mathbf{y}_i$ we need to calculate its first two moments. \n", - "\n", - "Recall that $\\boldsymbol{X}$ is a matrix of dimensionality $n\\times p$. The\n", - "notation above $\\mathbf{X}_{i,\\ast}$ means that we are looking at the\n", - "row number $i$ and perform a sum over all values $p$." - ] - }, - { - "cell_type": "markdown", - "id": "987be257", - "metadata": { - "editable": true - }, - "source": [ - "## Assumptions made\n", - "\n", - "The assumption we have made here can be summarized as (and this is going to be useful when we discuss the bias-variance trade off)\n", - "that there exists a function $f(\\boldsymbol{x})$ and a normal distributed error $\\boldsymbol{\\varepsilon}\\sim \\mathcal{N}(0, \\sigma^2)$\n", - "which describe our data" - ] - }, - { - "cell_type": "markdown", - "id": "3c6d8a07", - "metadata": { - "editable": true - }, - "source": [ - "$$\n", - "\\boldsymbol{y} = f(\\boldsymbol{x})+\\boldsymbol{\\varepsilon}\n", - "$$" - ] - }, - { - "cell_type": "markdown", - "id": "344d5c27", - "metadata": { - "editable": true - }, - "source": [ - "We approximate this function with our model from the solution of the linear regression equations, that is our\n", - "function $f$ is approximated by $\\boldsymbol{\\tilde{y}}$ where we want to minimize $(\\boldsymbol{y}-\\boldsymbol{\\tilde{y}})^2$, our MSE, with" - ] - }, - { - "cell_type": "markdown", - "id": "14d6b4ec", - "metadata": { - "editable": true - }, - "source": [ - "$$\n", - "\\boldsymbol{\\tilde{y}} = \\boldsymbol{X}\\boldsymbol{\\beta}.\n", - "$$" - ] - }, - { - "cell_type": "markdown", - "id": "f1c0ce92", - "metadata": { - "editable": true - }, - "source": [ - "## Expectation value and variance\n", - "\n", - "We can calculate the expectation value of $\\boldsymbol{y}$ for a given element $i$" - ] - }, - { - "cell_type": "markdown", - "id": "2528cc25", - "metadata": { - "editable": true - }, - "source": [ - "$$\n", - "\\begin{align*} \n", - "\\mathbb{E}(y_i) & =\n", - "\\mathbb{E}(\\mathbf{X}_{i, \\ast} \\, \\boldsymbol{\\beta}) + \\mathbb{E}(\\varepsilon_i)\n", - "\\, \\, \\, = \\, \\, \\, \\mathbf{X}_{i, \\ast} \\, \\beta, \n", - "\\end{align*}\n", - "$$" - ] - }, - { - "cell_type": "markdown", - "id": "6cd0ff77", - "metadata": { - "editable": true - }, - "source": [ - "while\n", - "its variance is" - ] - }, - { - "cell_type": "markdown", - "id": "740af370", - "metadata": { - "editable": true - }, - "source": [ - "$$\n", - "\\begin{align*} \\mbox{Var}(y_i) & = \\mathbb{E} \\{ [y_i\n", - "- \\mathbb{E}(y_i)]^2 \\} \\, \\, \\, = \\, \\, \\, \\mathbb{E} ( y_i^2 ) -\n", - "[\\mathbb{E}(y_i)]^2 \\\\ & = \\mathbb{E} [ ( \\mathbf{X}_{i, \\ast} \\,\n", - "\\beta + \\varepsilon_i )^2] - ( \\mathbf{X}_{i, \\ast} \\, \\boldsymbol{\\beta})^2 \\\\ &\n", - "= \\mathbb{E} [ ( \\mathbf{X}_{i, \\ast} \\, \\boldsymbol{\\beta})^2 + 2 \\varepsilon_i\n", - "\\mathbf{X}_{i, \\ast} \\, \\boldsymbol{\\beta} + \\varepsilon_i^2 ] - ( \\mathbf{X}_{i,\n", - "\\ast} \\, \\beta)^2 \\\\ & = ( \\mathbf{X}_{i, \\ast} \\, \\boldsymbol{\\beta})^2 + 2\n", - "\\mathbb{E}(\\varepsilon_i) \\mathbf{X}_{i, \\ast} \\, \\boldsymbol{\\beta} +\n", - "\\mathbb{E}(\\varepsilon_i^2 ) - ( \\mathbf{X}_{i, \\ast} \\, \\boldsymbol{\\beta})^2 \n", - "\\\\ & = \\mathbb{E}(\\varepsilon_i^2 ) \\, \\, \\, = \\, \\, \\,\n", - "\\mbox{Var}(\\varepsilon_i) \\, \\, \\, = \\, \\, \\, \\sigma^2. \n", - "\\end{align*}\n", - "$$" - ] - }, - { - "cell_type": "markdown", - "id": "c806e1e9", - "metadata": { - "editable": true - }, - "source": [ - "Hence, $y_i \\sim \\mathcal{N}( \\mathbf{X}_{i, \\ast} \\, \\boldsymbol{\\beta}, \\sigma^2)$, that is $\\boldsymbol{y}$ follows a normal distribution with \n", - "mean value $\\boldsymbol{X}\\boldsymbol{\\beta}$ and variance $\\sigma^2$ (not be confused with the singular values of the SVD)." - ] - }, - { - "cell_type": "markdown", - "id": "72059494", - "metadata": { - "editable": true - }, - "source": [ - "## Expectation value and variance for $\\boldsymbol{\\beta}$\n", - "\n", - "With the OLS expressions for the optimal parameters $\\boldsymbol{\\hat{\\beta}}$ we can evaluate the expectation value" - ] - }, - { - "cell_type": "markdown", - "id": "11212568", - "metadata": { - "editable": true - }, - "source": [ - "$$\n", - "\\mathbb{E}(\\boldsymbol{\\hat{\\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}.\n", - "$$" - ] - }, - { - "cell_type": "markdown", - "id": "a8949807", - "metadata": { - "editable": true - }, - "source": [ - "This means that the estimator of the regression parameters is unbiased.\n", - "\n", - "We can also calculate the variance\n", - "\n", - "The variance of the optimal value $\\boldsymbol{\\hat{\\beta}}$ is" - ] - }, - { - "cell_type": "markdown", - "id": "bf8deaae", - "metadata": { - "editable": true - }, - "source": [ - "$$\n", - "\\begin{eqnarray*}\n", - "\\mbox{Var}(\\boldsymbol{\\hat{\\beta}}) & = & \\mathbb{E} \\{ [\\boldsymbol{\\beta} - \\mathbb{E}(\\boldsymbol{\\beta})] [\\boldsymbol{\\beta} - \\mathbb{E}(\\boldsymbol{\\beta})]^{T} \\}\n", - "\\\\\n", - "& = & \\mathbb{E} \\{ [(\\mathbf{X}^{T} \\mathbf{X})^{-1} \\, \\mathbf{X}^{T} \\mathbf{Y} - \\boldsymbol{\\beta}] \\, [(\\mathbf{X}^{T} \\mathbf{X})^{-1} \\, \\mathbf{X}^{T} \\mathbf{Y} - \\boldsymbol{\\beta}]^{T} \\}\n", - "\\\\\n", - "% & = & \\mathbb{E} \\{ [(\\mathbf{X}^{T} \\mathbf{X})^{-1} \\, \\mathbf{X}^{T} \\mathbf{Y}] \\, [(\\mathbf{X}^{T} \\mathbf{X})^{-1} \\, \\mathbf{X}^{T} \\mathbf{Y}]^{T} \\} - \\boldsymbol{\\beta} \\, \\boldsymbol{\\beta}^{T}\n", - "% \\\\\n", - "% & = & \\mathbb{E} \\{ (\\mathbf{X}^{T} \\mathbf{X})^{-1} \\, \\mathbf{X}^{T} \\mathbf{Y} \\, \\mathbf{Y}^{T} \\, \\mathbf{X} \\, (\\mathbf{X}^{T} \\mathbf{X})^{-1} \\} - \\boldsymbol{\\beta} \\, \\boldsymbol{\\beta}^{T}\n", - "% \\\\\n", - "& = & (\\mathbf{X}^{T} \\mathbf{X})^{-1} \\, \\mathbf{X}^{T} \\, \\mathbb{E} \\{ \\mathbf{Y} \\, \\mathbf{Y}^{T} \\} \\, \\mathbf{X} \\, (\\mathbf{X}^{T} \\mathbf{X})^{-1} - \\boldsymbol{\\beta} \\, \\boldsymbol{\\beta}^{T}\n", - "\\\\\n", - "& = & (\\mathbf{X}^{T} \\mathbf{X})^{-1} \\, \\mathbf{X}^{T} \\, \\{ \\mathbf{X} \\, \\boldsymbol{\\beta} \\, \\boldsymbol{\\beta}^{T} \\, \\mathbf{X}^{T} + \\sigma^2 \\} \\, \\mathbf{X} \\, (\\mathbf{X}^{T} \\mathbf{X})^{-1} - \\boldsymbol{\\beta} \\, \\boldsymbol{\\beta}^{T}\n", - "% \\\\\n", - "% & = & (\\mathbf{X}^T \\mathbf{X})^{-1} \\, \\mathbf{X}^T \\, \\mathbf{X} \\, \\boldsymbol{\\beta} \\, \\boldsymbol{\\beta}^T \\, \\mathbf{X}^T \\, \\mathbf{X} \\, (\\mathbf{X}^T % \\mathbf{X})^{-1}\n", - "% \\\\\n", - "% & & + \\, \\, \\sigma^2 \\, (\\mathbf{X}^T \\mathbf{X})^{-1} \\, \\mathbf{X}^T \\, \\mathbf{X} \\, (\\mathbf{X}^T \\mathbf{X})^{-1} - \\boldsymbol{\\beta} \\boldsymbol{\\beta}^T\n", - "\\\\\n", - "& = & \\boldsymbol{\\beta} \\, \\boldsymbol{\\beta}^{T} + \\sigma^2 \\, (\\mathbf{X}^{T} \\mathbf{X})^{-1} - \\boldsymbol{\\beta} \\, \\boldsymbol{\\beta}^{T}\n", - "\\, \\, \\, = \\, \\, \\, \\sigma^2 \\, (\\mathbf{X}^{T} \\mathbf{X})^{-1},\n", - "\\end{eqnarray*}\n", - "$$" - ] - }, - { - "cell_type": "markdown", - "id": "9b017145", - "metadata": { - "editable": true - }, - "source": [ - "where we have used that $\\mathbb{E} (\\mathbf{Y} \\mathbf{Y}^{T}) =\n", - "\\mathbf{X} \\, \\boldsymbol{\\beta} \\, \\boldsymbol{\\beta}^{T} \\, \\mathbf{X}^{T} +\n", - "\\sigma^2 \\, \\mathbf{I}_{nn}$. From $\\mbox{Var}(\\boldsymbol{\\beta}) = \\sigma^2\n", - "\\, (\\mathbf{X}^{T} \\mathbf{X})^{-1}$, one obtains an estimate of the\n", - "variance of the estimate of the $j$-th regression coefficient:\n", - "$\\boldsymbol{\\sigma}^2 (\\boldsymbol{\\beta}_j ) = \\boldsymbol{\\sigma}^2 [(\\mathbf{X}^{T} \\mathbf{X})^{-1}]_{jj} $. This may be used to\n", - "construct a confidence interval for the estimates.\n", - "\n", - "In a similar way, we can obtain analytical expressions for say the\n", - "expectation values of the parameters $\\boldsymbol{\\beta}$ and their variance\n", - "when we employ Ridge regression, allowing us again to define a confidence interval. \n", - "\n", - "It is rather straightforward to show that" - ] - }, - { - "cell_type": "markdown", - "id": "9a3b13f5", - "metadata": { - "editable": true - }, - "source": [ - "$$\n", - "\\mathbb{E} \\big[ \\boldsymbol{\\beta}^{\\mathrm{Ridge}} \\big]=(\\mathbf{X}^{T} \\mathbf{X} + \\lambda \\mathbf{I}_{pp})^{-1} (\\mathbf{X}^{\\top} \\mathbf{X})\\boldsymbol{\\beta}^{\\mathrm{OLS}}.\n", - "$$" - ] - }, - { - "cell_type": "markdown", - "id": "9451c68a", - "metadata": { - "editable": true - }, - "source": [ - "We see clearly that \n", - "$\\mathbb{E} \\big[ \\boldsymbol{\\beta}^{\\mathrm{Ridge}} \\big] \\not= \\boldsymbol{\\beta}^{\\mathrm{OLS}}$ for any $\\lambda > 0$. We say then that the ridge estimator is biased.\n", - "\n", - "We can also compute the variance as" - ] - }, - { - "cell_type": "markdown", - "id": "52b40833", - "metadata": { - "editable": true - }, - "source": [ - "$$\n", - "\\mbox{Var}[\\boldsymbol{\\beta}^{\\mathrm{Ridge}}]=\\sigma^2[ \\mathbf{X}^{T} \\mathbf{X} + \\lambda \\mathbf{I} ]^{-1} \\mathbf{X}^{T} \\mathbf{X} \\{ [ \\mathbf{X}^{\\top} \\mathbf{X} + \\lambda \\mathbf{I} ]^{-1}\\}^{T},\n", - "$$" - ] - }, - { - "cell_type": "markdown", - "id": "a15fb2e5", - "metadata": { - "editable": true - }, - "source": [ - "and it is easy to see that if the parameter $\\lambda$ goes to infinity then the variance of Ridge parameters $\\boldsymbol{\\beta}$ goes to zero. \n", - "\n", - "With this, we can compute the difference" - ] - }, - { - "cell_type": "markdown", - "id": "9fc98887", - "metadata": { - "editable": true - }, - "source": [ - "$$\n", - "\\mbox{Var}[\\boldsymbol{\\beta}^{\\mathrm{OLS}}]-\\mbox{Var}(\\boldsymbol{\\beta}^{\\mathrm{Ridge}})=\\sigma^2 [ \\mathbf{X}^{T} \\mathbf{X} + \\lambda \\mathbf{I} ]^{-1}[ 2\\lambda\\mathbf{I} + \\lambda^2 (\\mathbf{X}^{T} \\mathbf{X})^{-1} ] \\{ [ \\mathbf{X}^{T} \\mathbf{X} + \\lambda \\mathbf{I} ]^{-1}\\}^{T}.\n", - "$$" - ] - }, - { - "cell_type": "markdown", - "id": "fb9ca21d", - "metadata": { - "editable": true - }, - "source": [ - "The difference is non-negative definite since each component of the\n", - "matrix product is non-negative definite. \n", - "This means the variance we obtain with the standard OLS will always for $\\lambda > 0$ be larger than the variance of $\\boldsymbol{\\beta}$ obtained with the Ridge estimator. This has interesting consequences when we discuss the so-called bias-variance trade-off below." - ] - }, - { - "cell_type": "markdown", - "id": "aa65553b", - "metadata": { - "editable": true - }, - "source": [ - "## Deriving OLS from a probability distribution\n", - "\n", - "Our basic assumption when we derived the OLS equations was to assume\n", - "that our output is determined by a given continuous function\n", - "$f(\\boldsymbol{x})$ and a random noise $\\boldsymbol{\\epsilon}$ given by the normal\n", - "distribution with zero mean value and an undetermined variance\n", - "$\\sigma^2$.\n", - "\n", - "We found above that the outputs $\\boldsymbol{y}$ have a mean value given by\n", - "$\\boldsymbol{X}\\hat{\\boldsymbol{\\beta}}$ and variance $\\sigma^2$. Since the entries to\n", - "the design matrix are not stochastic variables, we can assume that the\n", - "probability distribution of our targets is also a normal distribution\n", - "but now with mean value $\\boldsymbol{X}\\hat{\\boldsymbol{\\beta}}$. This means that a\n", - "single output $y_i$ is given by the Gaussian distribution" - ] - }, - { - "cell_type": "markdown", - "id": "50f53b28", - "metadata": { - "editable": true - }, - "source": [ - "$$\n", - "y_i\\sim \\mathcal{N}(\\boldsymbol{X}_{i,*}\\boldsymbol{\\beta}, \\sigma^2)=\\frac{1}{\\sqrt{2\\pi\\sigma^2}}\\exp{\\left[-\\frac{(y_i-\\boldsymbol{X}_{i,*}\\boldsymbol{\\beta})^2}{2\\sigma^2}\\right]}.\n", - "$$" - ] - }, - { - "cell_type": "markdown", - "id": "256983f1", - "metadata": { - "editable": true - }, - "source": [ - "## Independent and Identically Distrubuted (iid)\n", - "\n", - "We assume now that the various $y_i$ values are stochastically distributed according to the above Gaussian distribution. \n", - "We define this distribution as" - ] - }, - { - "cell_type": "markdown", - "id": "2ad84e01", - "metadata": { - "editable": true - }, - "source": [ - "$$\n", - "p(y_i, \\boldsymbol{X}\\vert\\boldsymbol{\\beta})=\\frac{1}{\\sqrt{2\\pi\\sigma^2}}\\exp{\\left[-\\frac{(y_i-\\boldsymbol{X}_{i,*}\\boldsymbol{\\beta})^2}{2\\sigma^2}\\right]},\n", - "$$" - ] - }, - { - "cell_type": "markdown", - "id": "b66da8f8", - "metadata": { - "editable": true - }, - "source": [ - "which reads as finding the likelihood of an event $y_i$ with the input variables $\\boldsymbol{X}$ given the parameters (to be determined) $\\boldsymbol{\\beta}$.\n", - "\n", - "Since these events are assumed to be independent and identicall distributed we can build the probability distribution function (PDF) for all possible event $\\boldsymbol{y}$ as the product of the single events, that is we have" - ] - }, - { - "cell_type": "markdown", - "id": "1488805d", - "metadata": { - "editable": true - }, - "source": [ - "$$\n", - "p(\\boldsymbol{y},\\boldsymbol{X}\\vert\\boldsymbol{\\beta})=\\prod_{i=0}^{n-1}\\frac{1}{\\sqrt{2\\pi\\sigma^2}}\\exp{\\left[-\\frac{(y_i-\\boldsymbol{X}_{i,*}\\boldsymbol{\\beta})^2}{2\\sigma^2}\\right]}=\\prod_{i=0}^{n-1}p(y_i,\\boldsymbol{X}\\vert\\boldsymbol{\\beta}).\n", - "$$" - ] - }, - { - "cell_type": "markdown", - "id": "1b8081db", - "metadata": { - "editable": true - }, - "source": [ - "We will write this in a more compact form reserving $\\boldsymbol{D}$ for the domain of events, including the ouputs (targets) and the inputs. That is\n", - "in case we have a simple one-dimensional input and output case" - ] - }, - { - "cell_type": "markdown", - "id": "236b5d6e", - "metadata": { - "editable": true - }, - "source": [ - "$$\n", - "\\boldsymbol{D}=[(x_0,y_0), (x_1,y_1),\\dots, (x_{n-1},y_{n-1})].\n", - "$$" - ] - }, - { - "cell_type": "markdown", - "id": "a3a7c862", - "metadata": { - "editable": true - }, - "source": [ - "In the more general case the various inputs should be replaced by the possible features represented by the input data set $\\boldsymbol{X}$. \n", - "We can now rewrite the above probability as" - ] - }, - { - "cell_type": "markdown", - "id": "917b0b1d", - "metadata": { - "editable": true - }, - "source": [ - "$$\n", - "p(\\boldsymbol{D}\\vert\\boldsymbol{\\beta})=\\prod_{i=0}^{n-1}\\frac{1}{\\sqrt{2\\pi\\sigma^2}}\\exp{\\left[-\\frac{(y_i-\\boldsymbol{X}_{i,*}\\boldsymbol{\\beta})^2}{2\\sigma^2}\\right]}.\n", - "$$" - ] - }, - { - "cell_type": "markdown", - "id": "b99286e3", - "metadata": { - "editable": true - }, - "source": [ - "It is a conditional probability (see below) and reads as the likelihood of a domain of events $\\boldsymbol{D}$ given a set of parameters $\\boldsymbol{\\beta}$." - ] - }, - { - "cell_type": "markdown", - "id": "80550e52", - "metadata": { - "editable": true - }, - "source": [ - "## Maximum Likelihood Estimation (MLE)\n", - "\n", - "In statistics, maximum likelihood estimation (MLE) is a method of\n", - "estimating the parameters of an assumed probability distribution,\n", - "given some observed data. This is achieved by maximizing a likelihood\n", - "function so that, under the assumed statistical model, the observed\n", - "data is the most probable. \n", - "\n", - "We will assume here that our events are given by the above Gaussian\n", - "distribution and we will determine the optimal parameters $\\beta$ by\n", - "maximizing the above PDF. However, computing the derivatives of a\n", - "product function is cumbersome and can easily lead to overflow and/or\n", - "underflowproblems, with potentials for loss of numerical precision.\n", - "\n", - "In practice, it is more convenient to maximize the logarithm of the\n", - "PDF because it is a monotonically increasing function of the argument.\n", - "Alternatively, and this will be our option, we will minimize the\n", - "negative of the logarithm since this is a monotonically decreasing\n", - "function.\n", - "\n", - "Note also that maximization/minimization of the logarithm of the PDF\n", - "is equivalent to the maximization/minimization of the function itself." - ] - }, - { - "cell_type": "markdown", - "id": "8897e74b", - "metadata": { - "editable": true - }, - "source": [ - "## A new Cost Function\n", - "\n", - "We could now define a new cost function to minimize, namely the negative logarithm of the above PDF" - ] - }, - { - "cell_type": "markdown", - "id": "8c3051a7", - "metadata": { - "editable": true - }, - "source": [ - "$$\n", - "C(\\boldsymbol{\\beta}=-\\log{\\prod_{i=0}^{n-1}p(y_i,\\boldsymbol{X}\\vert\\boldsymbol{\\beta})}=-\\sum_{i=0}^{n-1}\\log{p(y_i,\\boldsymbol{X}\\vert\\boldsymbol{\\beta})},\n", - "$$" - ] - }, - { - "cell_type": "markdown", - "id": "6894e38a", - "metadata": { - "editable": true - }, - "source": [ - "which becomes" - ] - }, - { - "cell_type": "markdown", - "id": "a726009e", - "metadata": { - "editable": true - }, - "source": [ - "$$\n", - "C(\\boldsymbol{\\beta}=\\frac{n}{2}\\log{2\\pi\\sigma^2}+\\frac{\\vert\\vert (\\boldsymbol{y}-\\boldsymbol{X}\\boldsymbol{\\beta})\\vert\\vert_2^2}{2\\sigma^2}.\n", - "$$" - ] - }, - { - "cell_type": "markdown", - "id": "3128e8f2", - "metadata": { - "editable": true - }, - "source": [ - "Taking the derivative of the *new* cost function with respect to the parameters $\\beta$ we recognize our familiar OLS equation, namely" - ] - }, - { - "cell_type": "markdown", - "id": "d0517e13", - "metadata": { - "editable": true - }, - "source": [ - "$$\n", - "\\boldsymbol{X}^T\\left(\\boldsymbol{y}-\\boldsymbol{X}\\boldsymbol{\\beta}\\right) =0,\n", - "$$" - ] - }, - { - "cell_type": "markdown", - "id": "6c551c61", - "metadata": { - "editable": true - }, - "source": [ - "which leads to the well-known OLS equation for the optimal paramters $\\beta$" - ] - }, - { - "cell_type": "markdown", - "id": "d402c6a3", - "metadata": { - "editable": true - }, - "source": [ - "$$\n", - "\\hat{\\boldsymbol{\\beta}}^{\\mathrm{OLS}}=\\left(\\boldsymbol{X}^T\\boldsymbol{X}\\right)^{-1}\\boldsymbol{X}^T\\boldsymbol{y}!\n", - "$$" - ] - }, - { - "cell_type": "markdown", - "id": "62495095", - "metadata": { - "editable": true - }, - "source": [ - "Before we make a similar analysis for Ridge and Lasso regression, we need a short reminder on statistics." - ] - }, - { - "cell_type": "markdown", - "id": "8fb09b93", - "metadata": { - "editable": true - }, - "source": [ - "## More basic Statistics and Bayes' theorem\n", - "\n", - "A central theorem in statistics is Bayes' theorem. This theorem plays a similar role as the good old Pythagoras' theorem in geometry.\n", - "Bayes' theorem is extremely simple to derive. But to do so we need some basic axioms from statistics.\n", - "\n", - "Assume we have two domains of events $X=[x_0,x_1,\\dots,x_{n-1}]$ and $Y=[y_0,y_1,\\dots,y_{n-1}]$.\n", - "\n", - "We define also the likelihood for $X$ and $Y$ as $p(X)$ and $p(Y)$ respectively.\n", - "The likelihood of a specific event $x_i$ (or $y_i$) is then written as $p(X=x_i)$ or just $p(x_i)=p_i$. \n", - "\n", - "**Union of events is given by.**" - ] - }, - { - "cell_type": "markdown", - "id": "745faaaa", - "metadata": { - "editable": true - }, - "source": [ - "$$\n", - "p(X \\cup Y)= p(X)+p(Y)-p(X \\cap Y).\n", - "$$" - ] - }, - { - "cell_type": "markdown", - "id": "8d1b2f69", - "metadata": { - "editable": true - }, - "source": [ - "**The product rule (aka joint probability) is given by.**" - ] - }, - { - "cell_type": "markdown", - "id": "af234eb7", - "metadata": { - "editable": true - }, - "source": [ - "$$\n", - "p(X \\cup Y)= p(X,Y)= p(X\\vert Y)p(Y)=p(Y\\vert X)p(X),\n", - "$$" - ] - }, - { - "cell_type": "markdown", - "id": "8ead38c0", - "metadata": { - "editable": true - }, - "source": [ - "where we read $p(X\\vert Y)$ as the likelihood of obtaining $X$ given $Y$.\n", - "\n", - "If we have independent events then $p(X,Y)=p(X)p(Y)$." - ] - }, - { - "cell_type": "markdown", - "id": "35aa26f0", - "metadata": { - "editable": true - }, - "source": [ - "## Marginal Probability\n", - "\n", - "The marginal probability is defined in terms of only one of the set of variables $X,Y$. For a discrete probability we have" - ] - }, - { - "cell_type": "markdown", - "id": "262cf2ba", - "metadata": { - "editable": true - }, - "source": [ - "$$\n", - "p(X)=\\sum_{i=0}^{n-1}p(X,Y=y_i)=\\sum_{i=0}^{n-1}p(X\\vert Y=y_i)p(Y=y_i)=\\sum_{i=0}^{n-1}p(X\\vert y_i)p(y_i).\n", - "$$" - ] - }, - { - "cell_type": "markdown", - "id": "4c322386", - "metadata": { - "editable": true - }, - "source": [ - "## Conditional Probability\n", - "\n", - "The conditional probability, if $p(Y) > 0$, is" - ] - }, - { - "cell_type": "markdown", - "id": "67c4e730", - "metadata": { - "editable": true - }, - "source": [ - "$$\n", - "p(X\\vert Y)= \\frac{p(X,Y)}{p(Y)}=\\frac{p(X,Y)}{\\sum_{i=0}^{n-1}p(Y\\vert X=x_i)p(x_i)}.\n", - "$$" - ] - }, - { - "cell_type": "markdown", - "id": "06a9c929", - "metadata": { - "editable": true - }, - "source": [ - "## Bayes' Theorem\n", - "\n", - "If we combine the conditional probability with the marginal probability and the standard product rule, we have" - ] - }, - { - "cell_type": "markdown", - "id": "6ae0b543", - "metadata": { - "editable": true - }, - "source": [ - "$$\n", - "p(X\\vert Y)= \\frac{p(X,Y)}{p(Y)},\n", - "$$" - ] - }, - { - "cell_type": "markdown", - "id": "ea757688", - "metadata": { - "editable": true - }, - "source": [ - "which we can rewrite as" - ] - }, - { - "cell_type": "markdown", - "id": "cacc60d2", - "metadata": { - "editable": true - }, - "source": [ - "$$\n", - "p(X\\vert Y)= \\frac{p(X,Y)}{\\sum_{i=0}^{n-1}p(Y\\vert X=x_i)p(x_i)}=\\frac{p(Y\\vert X)p(X)}{\\sum_{i=0}^{n-1}p(Y\\vert X=x_i)p(x_i)},\n", - "$$" - ] - }, - { - "cell_type": "markdown", - "id": "c8f425c7", - "metadata": { - "editable": true - }, - "source": [ - "which is Bayes' theorem. It allows us to evaluate the uncertainty in in $X$ after we have observed $Y$. We can easily interchange $X$ with $Y$." - ] - }, - { - "cell_type": "markdown", - "id": "0e2f3195", - "metadata": { - "editable": true - }, - "source": [ - "## Interpretations of Bayes' Theorem\n", - "\n", - "The quantity $p(Y\\vert X)$ on the right-hand side of the theorem is\n", - "evaluated for the observed data $Y$ and can be viewed as a function of\n", - "the parameter space represented by $X$. This function is not\n", - "necesseraly normalized and is normally called the likelihood function.\n", - "\n", - "The function $p(X)$ on the right hand side is called the prior while the function on the left hand side is the called the posterior probability. The denominator on the right hand side serves as a normalization factor for the posterior distribution.\n", - "\n", - "Let us try to illustrate Bayes' theorem through an example." - ] - }, - { - "cell_type": "markdown", - "id": "f0501701", - "metadata": { - "editable": true - }, - "source": [ - "## Example of Usage of Bayes' theorem\n", - "\n", - "Let us suppose that you are undergoing a series of mammography scans in\n", - "order to rule out possible breast cancer cases. We define the\n", - "sensitivity for a positive event by the variable $X$. It takes binary\n", - "values with $X=1$ representing a positive event and $X=0$ being a\n", - "negative event. We reserve $Y$ as a classification parameter for\n", - "either a negative or a positive breast cancer confirmation. (Short note on wordings: positive here means having breast cancer, although none of us would consider this being a positive thing).\n", - "\n", - "We let $Y=1$ represent the the case of having breast cancer and $Y=0$ as not.\n", - "\n", - "Let us assume that if you have breast cancer, the test will be positive with a probability of $0.8$, that is we have" - ] - }, - { - "cell_type": "markdown", - "id": "c4114c79", - "metadata": { - "editable": true - }, - "source": [ - "$$\n", - "p(X=1\\vert Y=1) =0.8.\n", - "$$" - ] - }, - { - "cell_type": "markdown", - "id": "8dcf2619", - "metadata": { - "editable": true - }, - "source": [ - "This obviously sounds scary since many would conclude that if the test is positive, there is a likelihood of $80\\%$ for having cancer.\n", - "It is however not correct, as the following Bayesian analysis shows." - ] - }, - { - "cell_type": "markdown", - "id": "c0ebfeb6", - "metadata": { - "editable": true - }, - "source": [ - "## Doing it correctly\n", - "\n", - "If we look at various national surveys on breast cancer, the general likelihood of developing breast cancer is a very small number.\n", - "Let us assume that the prior probability in the population as a whole is" - ] - }, - { - "cell_type": "markdown", - "id": "2a174072", - "metadata": { - "editable": true - }, - "source": [ - "$$\n", - "p(Y=1) =0.004.\n", - "$$" - ] - }, - { - "cell_type": "markdown", - "id": "5084d033", - "metadata": { - "editable": true - }, - "source": [ - "We need also to account for the fact that the test may produce a false positive result (false alarm). Let us here assume that we have" - ] - }, - { - "cell_type": "markdown", - "id": "c011b70a", - "metadata": { - "editable": true - }, - "source": [ - "$$\n", - "p(X=1\\vert Y=0) =0.1.\n", - "$$" - ] - }, - { - "cell_type": "markdown", - "id": "fb7da2bd", - "metadata": { - "editable": true - }, - "source": [ - "Using Bayes' theorem we can then find the posterior probability that the person has breast cancer in case of a positive test, that is we can compute" - ] - }, - { - "cell_type": "markdown", - "id": "492beb10", - "metadata": { - "editable": true - }, - "source": [ - "$$\n", - "p(Y=1\\vert X=1)=\\frac{p(X=1\\vert Y=1)p(Y=1)}{p(X=1\\vert Y=1)p(Y=1)+p(X=1\\vert Y=0)p(Y=0)}=\\frac{0.8\\times 0.004}{0.8\\times 0.004+0.1\\times 0.996}=0.031.\n", - "$$" - ] - }, - { - "cell_type": "markdown", - "id": "753962c7", - "metadata": { - "editable": true - }, - "source": [ - "That is, in case of a positive test, there is only a $3\\%$ chance of having breast cancer!" - ] - }, - { - "cell_type": "markdown", - "id": "fa58071d", - "metadata": { - "editable": true - }, - "source": [ - "## Bayes' Theorem and Ridge and Lasso Regression\n", + "## Test Function for what happens with OLS, Ridge and Lasso\n", "\n", "Hitherto we have discussed Ridge and Lasso regression in terms of a\n", "linear analysis. This may to many of you feel rather technical and\n", "perhaps not that intuitive. The question is whether we can develop a\n", "more intuitive way of understanding what Ridge and Lasso express.\n", "\n", - "Before we proceed let us perform a Ridge, Lasso and OLS analysis of a polynomial fit." - ] - }, - { - "cell_type": "markdown", - "id": "6e12d6bc", - "metadata": { - "editable": true - }, - "source": [ - "## Test Function for what happens with OLS, Ridge and Lasso\n", + "Before we proceed let us perform a Ridge, Lasso and OLS analysis of a polynomial fit. \n", "\n", "We will play around with a study of the values for the optimal\n", "parameters $\\boldsymbol{\\beta}$ using OLS, Ridge and Lasso regression. For\n", @@ -3849,7 +2906,7 @@ { "cell_type": "code", "execution_count": 11, - "id": "e58e929b", + "id": "94659a7f", "metadata": { "collapsed": false, "editable": true @@ -3876,10 +2933,9 @@ "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", + "for polydegree in range(1,Maxpolydegree):\n", + " X[:,polydegree] = x**(polydegree)\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", @@ -3923,7 +2979,7 @@ }, { "cell_type": "markdown", - "id": "dfcc4646", + "id": "2056e242", "metadata": { "editable": true }, @@ -3933,43 +2989,475 @@ }, { "cell_type": "markdown", - "id": "18a51f22", + "id": "3a93dae3", "metadata": { "editable": true }, "source": [ - "## Invoking Bayes' theorem\n", + "## Linking the regression analysis with a statistical interpretation\n", "\n", - "Using Bayes' theorem we can gain a better intuition about Ridge and Lasso regression. \n", + "We will now couple the discussions of ordinary least squares, Ridge\n", + "and Lasso regression with a statistical interpretation, that is we\n", + "move from a linear algebra analysis to a statistical analysis. In\n", + "particular, we will focus on what the regularization terms can result\n", + "in. We will amongst other things show that the regularization\n", + "parameter can reduce considerably the variance of the parameters\n", + "$\\beta$.\n", "\n", - "For ordinary least squares we postulated that the maximum likelihood for the doamin of events $\\boldsymbol{D}$ (one-dimensional case)" + "The\n", + "advantage of doing linear regression is that we actually end up with\n", + "analytical expressions for several statistical quantities. \n", + "Standard least squares and Ridge regression allow us to\n", + "derive quantities like the variance and other expectation values in a\n", + "rather straightforward way.\n", + "\n", + "It is assumed that $\\varepsilon_i\n", + "\\sim \\mathcal{N}(0, \\sigma^2)$ and the $\\varepsilon_{i}$ are\n", + "independent, i.e.:" ] }, { "cell_type": "markdown", - "id": "eac138af", + "id": "265c2329", "metadata": { "editable": true }, "source": [ "$$\n", - "\\boldsymbol{D}=[(x_0,y_0), (x_1,y_1),\\dots, (x_{n-1},y_{n-1})],\n", + "\\begin{align*} \n", + "\\mbox{Cov}(\\varepsilon_{i_1},\n", + "\\varepsilon_{i_2}) & = \\left\\{ \\begin{array}{lcc} \\sigma^2 & \\mbox{if}\n", + "& i_1 = i_2, \\\\ 0 & \\mbox{if} & i_1 \\not= i_2. \\end{array} \\right.\n", + "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", - "id": "79162e39", + "id": "1eb62f34", "metadata": { "editable": true }, "source": [ - "is given by" + "The randomness of $\\varepsilon_i$ implies that\n", + "$\\mathbf{y}_i$ is also a random variable. In particular,\n", + "$\\mathbf{y}_i$ is normally distributed, because $\\varepsilon_i \\sim\n", + "\\mathcal{N}(0, \\sigma^2)$ and $\\mathbf{X}_{i,\\ast} \\, \\boldsymbol{\\beta}$ is a\n", + "non-random scalar. To specify the parameters of the distribution of\n", + "$\\mathbf{y}_i$ we need to calculate its first two moments. \n", + "\n", + "Recall that $\\boldsymbol{X}$ is a matrix of dimensionality $n\\times p$. The\n", + "notation above $\\mathbf{X}_{i,\\ast}$ means that we are looking at the\n", + "row number $i$ and perform a sum over all values $p$." ] }, { "cell_type": "markdown", - "id": "a386a495", + "id": "93523fed", + "metadata": { + "editable": true + }, + "source": [ + "## Assumptions made\n", + "\n", + "The assumption we have made here can be summarized as (and this is going to be useful when we discuss the bias-variance trade off)\n", + "that there exists a function $f(\\boldsymbol{x})$ and a normal distributed error $\\boldsymbol{\\varepsilon}\\sim \\mathcal{N}(0, \\sigma^2)$\n", + "which describe our data" + ] + }, + { + "cell_type": "markdown", + "id": "5b2fe793", + "metadata": { + "editable": true + }, + "source": [ + "$$\n", + "\\boldsymbol{y} = f(\\boldsymbol{x})+\\boldsymbol{\\varepsilon}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "645cfe7a", + "metadata": { + "editable": true + }, + "source": [ + "We approximate this function with our model from the solution of the linear regression equations, that is our\n", + "function $f$ is approximated by $\\boldsymbol{\\tilde{y}}$ where we want to minimize $(\\boldsymbol{y}-\\boldsymbol{\\tilde{y}})^2$, our MSE, with" + ] + }, + { + "cell_type": "markdown", + "id": "acd315af", + "metadata": { + "editable": true + }, + "source": [ + "$$\n", + "\\boldsymbol{\\tilde{y}} = \\boldsymbol{X}\\boldsymbol{\\beta}.\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "3e5bcdc9", + "metadata": { + "editable": true + }, + "source": [ + "## Expectation value and variance\n", + "\n", + "We can calculate the expectation value of $\\boldsymbol{y}$ for a given element $i$" + ] + }, + { + "cell_type": "markdown", + "id": "445b0b0d", + "metadata": { + "editable": true + }, + "source": [ + "$$\n", + "\\begin{align*} \n", + "\\mathbb{E}(y_i) & =\n", + "\\mathbb{E}(\\mathbf{X}_{i, \\ast} \\, \\boldsymbol{\\beta}) + \\mathbb{E}(\\varepsilon_i)\n", + "\\, \\, \\, = \\, \\, \\, \\mathbf{X}_{i, \\ast} \\, \\beta, \n", + "\\end{align*}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "9ba5bb36", + "metadata": { + "editable": true + }, + "source": [ + "while\n", + "its variance is" + ] + }, + { + "cell_type": "markdown", + "id": "cd53bf61", + "metadata": { + "editable": true + }, + "source": [ + "$$\n", + "\\begin{align*} \\mbox{Var}(y_i) & = \\mathbb{E} \\{ [y_i\n", + "- \\mathbb{E}(y_i)]^2 \\} \\, \\, \\, = \\, \\, \\, \\mathbb{E} ( y_i^2 ) -\n", + "[\\mathbb{E}(y_i)]^2 \\\\ & = \\mathbb{E} [ ( \\mathbf{X}_{i, \\ast} \\,\n", + "\\beta + \\varepsilon_i )^2] - ( \\mathbf{X}_{i, \\ast} \\, \\boldsymbol{\\beta})^2 \\\\ &\n", + "= \\mathbb{E} [ ( \\mathbf{X}_{i, \\ast} \\, \\boldsymbol{\\beta})^2 + 2 \\varepsilon_i\n", + "\\mathbf{X}_{i, \\ast} \\, \\boldsymbol{\\beta} + \\varepsilon_i^2 ] - ( \\mathbf{X}_{i,\n", + "\\ast} \\, \\beta)^2 \\\\ & = ( \\mathbf{X}_{i, \\ast} \\, \\boldsymbol{\\beta})^2 + 2\n", + "\\mathbb{E}(\\varepsilon_i) \\mathbf{X}_{i, \\ast} \\, \\boldsymbol{\\beta} +\n", + "\\mathbb{E}(\\varepsilon_i^2 ) - ( \\mathbf{X}_{i, \\ast} \\, \\boldsymbol{\\beta})^2 \n", + "\\\\ & = \\mathbb{E}(\\varepsilon_i^2 ) \\, \\, \\, = \\, \\, \\,\n", + "\\mbox{Var}(\\varepsilon_i) \\, \\, \\, = \\, \\, \\, \\sigma^2. \n", + "\\end{align*}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "3d7101bc", + "metadata": { + "editable": true + }, + "source": [ + "Hence, $y_i \\sim \\mathcal{N}( \\mathbf{X}_{i, \\ast} \\, \\boldsymbol{\\beta}, \\sigma^2)$, that is $\\boldsymbol{y}$ follows a normal distribution with \n", + "mean value $\\boldsymbol{X}\\boldsymbol{\\beta}$ and variance $\\sigma^2$ (not be confused with the singular values of the SVD)." + ] + }, + { + "cell_type": "markdown", + "id": "3670204e", + "metadata": { + "editable": true + }, + "source": [ + "## Expectation value and variance for $\\boldsymbol{\\beta}$\n", + "\n", + "With the OLS expressions for the optimal parameters $\\boldsymbol{\\hat{\\beta}}$ we can evaluate the expectation value" + ] + }, + { + "cell_type": "markdown", + "id": "09f3c5ea", + "metadata": { + "editable": true + }, + "source": [ + "$$\n", + "\\mathbb{E}(\\boldsymbol{\\hat{\\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}.\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "814246ef", + "metadata": { + "editable": true + }, + "source": [ + "This means that the estimator of the regression parameters is unbiased.\n", + "\n", + "We can also calculate the variance\n", + "\n", + "The variance of the optimal value $\\boldsymbol{\\hat{\\beta}}$ is" + ] + }, + { + "cell_type": "markdown", + "id": "e831728d", + "metadata": { + "editable": true + }, + "source": [ + "$$\n", + "\\begin{eqnarray*}\n", + "\\mbox{Var}(\\boldsymbol{\\hat{\\beta}}) & = & \\mathbb{E} \\{ [\\boldsymbol{\\beta} - \\mathbb{E}(\\boldsymbol{\\beta})] [\\boldsymbol{\\beta} - \\mathbb{E}(\\boldsymbol{\\beta})]^{T} \\}\n", + "\\\\\n", + "& = & \\mathbb{E} \\{ [(\\mathbf{X}^{T} \\mathbf{X})^{-1} \\, \\mathbf{X}^{T} \\mathbf{Y} - \\boldsymbol{\\beta}] \\, [(\\mathbf{X}^{T} \\mathbf{X})^{-1} \\, \\mathbf{X}^{T} \\mathbf{Y} - \\boldsymbol{\\beta}]^{T} \\}\n", + "\\\\\n", + "% & = & \\mathbb{E} \\{ [(\\mathbf{X}^{T} \\mathbf{X})^{-1} \\, \\mathbf{X}^{T} \\mathbf{Y}] \\, [(\\mathbf{X}^{T} \\mathbf{X})^{-1} \\, \\mathbf{X}^{T} \\mathbf{Y}]^{T} \\} - \\boldsymbol{\\beta} \\, \\boldsymbol{\\beta}^{T}\n", + "% \\\\\n", + "% & = & \\mathbb{E} \\{ (\\mathbf{X}^{T} \\mathbf{X})^{-1} \\, \\mathbf{X}^{T} \\mathbf{Y} \\, \\mathbf{Y}^{T} \\, \\mathbf{X} \\, (\\mathbf{X}^{T} \\mathbf{X})^{-1} \\} - \\boldsymbol{\\beta} \\, \\boldsymbol{\\beta}^{T}\n", + "% \\\\\n", + "& = & (\\mathbf{X}^{T} \\mathbf{X})^{-1} \\, \\mathbf{X}^{T} \\, \\mathbb{E} \\{ \\mathbf{Y} \\, \\mathbf{Y}^{T} \\} \\, \\mathbf{X} \\, (\\mathbf{X}^{T} \\mathbf{X})^{-1} - \\boldsymbol{\\beta} \\, \\boldsymbol{\\beta}^{T}\n", + "\\\\\n", + "& = & (\\mathbf{X}^{T} \\mathbf{X})^{-1} \\, \\mathbf{X}^{T} \\, \\{ \\mathbf{X} \\, \\boldsymbol{\\beta} \\, \\boldsymbol{\\beta}^{T} \\, \\mathbf{X}^{T} + \\sigma^2 \\} \\, \\mathbf{X} \\, (\\mathbf{X}^{T} \\mathbf{X})^{-1} - \\boldsymbol{\\beta} \\, \\boldsymbol{\\beta}^{T}\n", + "% \\\\\n", + "% & = & (\\mathbf{X}^T \\mathbf{X})^{-1} \\, \\mathbf{X}^T \\, \\mathbf{X} \\, \\boldsymbol{\\beta} \\, \\boldsymbol{\\beta}^T \\, \\mathbf{X}^T \\, \\mathbf{X} \\, (\\mathbf{X}^T % \\mathbf{X})^{-1}\n", + "% \\\\\n", + "% & & + \\, \\, \\sigma^2 \\, (\\mathbf{X}^T \\mathbf{X})^{-1} \\, \\mathbf{X}^T \\, \\mathbf{X} \\, (\\mathbf{X}^T \\mathbf{X})^{-1} - \\boldsymbol{\\beta} \\boldsymbol{\\beta}^T\n", + "\\\\\n", + "& = & \\boldsymbol{\\beta} \\, \\boldsymbol{\\beta}^{T} + \\sigma^2 \\, (\\mathbf{X}^{T} \\mathbf{X})^{-1} - \\boldsymbol{\\beta} \\, \\boldsymbol{\\beta}^{T}\n", + "\\, \\, \\, = \\, \\, \\, \\sigma^2 \\, (\\mathbf{X}^{T} \\mathbf{X})^{-1},\n", + "\\end{eqnarray*}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "522aaab4", + "metadata": { + "editable": true + }, + "source": [ + "where we have used that $\\mathbb{E} (\\mathbf{Y} \\mathbf{Y}^{T}) =\n", + "\\mathbf{X} \\, \\boldsymbol{\\beta} \\, \\boldsymbol{\\beta}^{T} \\, \\mathbf{X}^{T} +\n", + "\\sigma^2 \\, \\mathbf{I}_{nn}$. From $\\mbox{Var}(\\boldsymbol{\\beta}) = \\sigma^2\n", + "\\, (\\mathbf{X}^{T} \\mathbf{X})^{-1}$, one obtains an estimate of the\n", + "variance of the estimate of the $j$-th regression coefficient:\n", + "$\\boldsymbol{\\sigma}^2 (\\boldsymbol{\\beta}_j ) = \\boldsymbol{\\sigma}^2 [(\\mathbf{X}^{T} \\mathbf{X})^{-1}]_{jj} $. This may be used to\n", + "construct a confidence interval for the estimates.\n", + "\n", + "In a similar way, we can obtain analytical expressions for say the\n", + "expectation values of the parameters $\\boldsymbol{\\beta}$ and their variance\n", + "when we employ Ridge regression, allowing us again to define a confidence interval. \n", + "\n", + "It is rather straightforward to show that" + ] + }, + { + "cell_type": "markdown", + "id": "0dfec07a", + "metadata": { + "editable": true + }, + "source": [ + "$$\n", + "\\mathbb{E} \\big[ \\boldsymbol{\\beta}^{\\mathrm{Ridge}} \\big]=(\\mathbf{X}^{T} \\mathbf{X} + \\lambda \\mathbf{I}_{pp})^{-1} (\\mathbf{X}^{\\top} \\mathbf{X})\\boldsymbol{\\beta}^{\\mathrm{OLS}}.\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "c0b31a8b", + "metadata": { + "editable": true + }, + "source": [ + "We see clearly that \n", + "$\\mathbb{E} \\big[ \\boldsymbol{\\beta}^{\\mathrm{Ridge}} \\big] \\not= \\boldsymbol{\\beta}^{\\mathrm{OLS}}$ for any $\\lambda > 0$. We say then that the ridge estimator is biased.\n", + "\n", + "We can also compute the variance as" + ] + }, + { + "cell_type": "markdown", + "id": "287d622e", + "metadata": { + "editable": true + }, + "source": [ + "$$\n", + "\\mbox{Var}[\\boldsymbol{\\beta}^{\\mathrm{Ridge}}]=\\sigma^2[ \\mathbf{X}^{T} \\mathbf{X} + \\lambda \\mathbf{I} ]^{-1} \\mathbf{X}^{T} \\mathbf{X} \\{ [ \\mathbf{X}^{\\top} \\mathbf{X} + \\lambda \\mathbf{I} ]^{-1}\\}^{T},\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "0afd8826", + "metadata": { + "editable": true + }, + "source": [ + "and it is easy to see that if the parameter $\\lambda$ goes to infinity then the variance of Ridge parameters $\\boldsymbol{\\beta}$ goes to zero. \n", + "\n", + "With this, we can compute the difference" + ] + }, + { + "cell_type": "markdown", + "id": "bb005c68", + "metadata": { + "editable": true + }, + "source": [ + "$$\n", + "\\mbox{Var}[\\boldsymbol{\\beta}^{\\mathrm{OLS}}]-\\mbox{Var}(\\boldsymbol{\\beta}^{\\mathrm{Ridge}})=\\sigma^2 [ \\mathbf{X}^{T} \\mathbf{X} + \\lambda \\mathbf{I} ]^{-1}[ 2\\lambda\\mathbf{I} + \\lambda^2 (\\mathbf{X}^{T} \\mathbf{X})^{-1} ] \\{ [ \\mathbf{X}^{T} \\mathbf{X} + \\lambda \\mathbf{I} ]^{-1}\\}^{T}.\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "2d9cb9a1", + "metadata": { + "editable": true + }, + "source": [ + "The difference is non-negative definite since each component of the\n", + "matrix product is non-negative definite. \n", + "This means the variance we obtain with the standard OLS will always for $\\lambda > 0$ be larger than the variance of $\\boldsymbol{\\beta}$ obtained with the Ridge estimator. This has interesting consequences when we discuss the so-called bias-variance trade-off below." + ] + }, + { + "cell_type": "markdown", + "id": "4760d2c7", + "metadata": { + "editable": true + }, + "source": [ + "## Deriving OLS from a probability distribution\n", + "\n", + "Our basic assumption when we derived the OLS equations was to assume\n", + "that our output is determined by a given continuous function\n", + "$f(\\boldsymbol{x})$ and a random noise $\\boldsymbol{\\epsilon}$ given by the normal\n", + "distribution with zero mean value and an undetermined variance\n", + "$\\sigma^2$.\n", + "\n", + "We found above that the outputs $\\boldsymbol{y}$ have a mean value given by\n", + "$\\boldsymbol{X}\\hat{\\boldsymbol{\\beta}}$ and variance $\\sigma^2$. Since the entries to\n", + "the design matrix are not stochastic variables, we can assume that the\n", + "probability distribution of our targets is also a normal distribution\n", + "but now with mean value $\\boldsymbol{X}\\hat{\\boldsymbol{\\beta}}$. This means that a\n", + "single output $y_i$ is given by the Gaussian distribution" + ] + }, + { + "cell_type": "markdown", + "id": "5af957ae", + "metadata": { + "editable": true + }, + "source": [ + "$$\n", + "y_i\\sim \\mathcal{N}(\\boldsymbol{X}_{i,*}\\boldsymbol{\\beta}, \\sigma^2)=\\frac{1}{\\sqrt{2\\pi\\sigma^2}}\\exp{\\left[-\\frac{(y_i-\\boldsymbol{X}_{i,*}\\boldsymbol{\\beta})^2}{2\\sigma^2}\\right]}.\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "fd6a18ae", + "metadata": { + "editable": true + }, + "source": [ + "## Independent and Identically Distrubuted (iid)\n", + "\n", + "We assume now that the various $y_i$ values are stochastically distributed according to the above Gaussian distribution. \n", + "We define this distribution as" + ] + }, + { + "cell_type": "markdown", + "id": "401d6e0a", + "metadata": { + "editable": true + }, + "source": [ + "$$\n", + "p(y_i, \\boldsymbol{X}\\vert\\boldsymbol{\\beta})=\\frac{1}{\\sqrt{2\\pi\\sigma^2}}\\exp{\\left[-\\frac{(y_i-\\boldsymbol{X}_{i,*}\\boldsymbol{\\beta})^2}{2\\sigma^2}\\right]},\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "48e126ea", + "metadata": { + "editable": true + }, + "source": [ + "which reads as finding the likelihood of an event $y_i$ with the input variables $\\boldsymbol{X}$ given the parameters (to be determined) $\\boldsymbol{\\beta}$.\n", + "\n", + "Since these events are assumed to be independent and identicall distributed we can build the probability distribution function (PDF) for all possible event $\\boldsymbol{y}$ as the product of the single events, that is we have" + ] + }, + { + "cell_type": "markdown", + "id": "34eda977", + "metadata": { + "editable": true + }, + "source": [ + "$$\n", + "p(\\boldsymbol{y},\\boldsymbol{X}\\vert\\boldsymbol{\\beta})=\\prod_{i=0}^{n-1}\\frac{1}{\\sqrt{2\\pi\\sigma^2}}\\exp{\\left[-\\frac{(y_i-\\boldsymbol{X}_{i,*}\\boldsymbol{\\beta})^2}{2\\sigma^2}\\right]}=\\prod_{i=0}^{n-1}p(y_i,\\boldsymbol{X}\\vert\\boldsymbol{\\beta}).\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "90714c71", + "metadata": { + "editable": true + }, + "source": [ + "We will write this in a more compact form reserving $\\boldsymbol{D}$ for the domain of events, including the ouputs (targets) and the inputs. That is\n", + "in case we have a simple one-dimensional input and output case" + ] + }, + { + "cell_type": "markdown", + "id": "74731195", + "metadata": { + "editable": true + }, + "source": [ + "$$\n", + "\\boldsymbol{D}=[(x_0,y_0), (x_1,y_1),\\dots, (x_{n-1},y_{n-1})].\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "1bc063bc", + "metadata": { + "editable": true + }, + "source": [ + "In the more general case the various inputs should be replaced by the possible features represented by the input data set $\\boldsymbol{X}$. \n", + "We can now rewrite the above probability as" + ] + }, + { + "cell_type": "markdown", + "id": "a9de140e", "metadata": { "editable": true }, @@ -3981,7 +3469,508 @@ }, { "cell_type": "markdown", - "id": "8684f5af", + "id": "b8c38a5a", + "metadata": { + "editable": true + }, + "source": [ + "It is a conditional probability (see below) and reads as the likelihood of a domain of events $\\boldsymbol{D}$ given a set of parameters $\\boldsymbol{\\beta}$." + ] + }, + { + "cell_type": "markdown", + "id": "2b9c970d", + "metadata": { + "editable": true + }, + "source": [ + "## Maximum Likelihood Estimation (MLE)\n", + "\n", + "In statistics, maximum likelihood estimation (MLE) is a method of\n", + "estimating the parameters of an assumed probability distribution,\n", + "given some observed data. This is achieved by maximizing a likelihood\n", + "function so that, under the assumed statistical model, the observed\n", + "data is the most probable. \n", + "\n", + "We will assume here that our events are given by the above Gaussian\n", + "distribution and we will determine the optimal parameters $\\beta$ by\n", + "maximizing the above PDF. However, computing the derivatives of a\n", + "product function is cumbersome and can easily lead to overflow and/or\n", + "underflowproblems, with potentials for loss of numerical precision.\n", + "\n", + "In practice, it is more convenient to maximize the logarithm of the\n", + "PDF because it is a monotonically increasing function of the argument.\n", + "Alternatively, and this will be our option, we will minimize the\n", + "negative of the logarithm since this is a monotonically decreasing\n", + "function.\n", + "\n", + "Note also that maximization/minimization of the logarithm of the PDF\n", + "is equivalent to the maximization/minimization of the function itself." + ] + }, + { + "cell_type": "markdown", + "id": "9f18acb3", + "metadata": { + "editable": true + }, + "source": [ + "## A new Cost Function\n", + "\n", + "We could now define a new cost function to minimize, namely the negative logarithm of the above PDF" + ] + }, + { + "cell_type": "markdown", + "id": "db679b3a", + "metadata": { + "editable": true + }, + "source": [ + "$$\n", + "C(\\boldsymbol{\\beta}=-\\log{\\prod_{i=0}^{n-1}p(y_i,\\boldsymbol{X}\\vert\\boldsymbol{\\beta})}=-\\sum_{i=0}^{n-1}\\log{p(y_i,\\boldsymbol{X}\\vert\\boldsymbol{\\beta})},\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "6b743abf", + "metadata": { + "editable": true + }, + "source": [ + "which becomes" + ] + }, + { + "cell_type": "markdown", + "id": "43048bc0", + "metadata": { + "editable": true + }, + "source": [ + "$$\n", + "C(\\boldsymbol{\\beta}=\\frac{n}{2}\\log{2\\pi\\sigma^2}+\\frac{\\vert\\vert (\\boldsymbol{y}-\\boldsymbol{X}\\boldsymbol{\\beta})\\vert\\vert_2^2}{2\\sigma^2}.\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "53939c7a", + "metadata": { + "editable": true + }, + "source": [ + "Taking the derivative of the *new* cost function with respect to the parameters $\\beta$ we recognize our familiar OLS equation, namely" + ] + }, + { + "cell_type": "markdown", + "id": "22af3bc5", + "metadata": { + "editable": true + }, + "source": [ + "$$\n", + "\\boldsymbol{X}^T\\left(\\boldsymbol{y}-\\boldsymbol{X}\\boldsymbol{\\beta}\\right) =0,\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "89e9f36e", + "metadata": { + "editable": true + }, + "source": [ + "which leads to the well-known OLS equation for the optimal paramters $\\beta$" + ] + }, + { + "cell_type": "markdown", + "id": "10eef97b", + "metadata": { + "editable": true + }, + "source": [ + "$$\n", + "\\hat{\\boldsymbol{\\beta}}^{\\mathrm{OLS}}=\\left(\\boldsymbol{X}^T\\boldsymbol{X}\\right)^{-1}\\boldsymbol{X}^T\\boldsymbol{y}!\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "6f491288", + "metadata": { + "editable": true + }, + "source": [ + "Before we make a similar analysis for Ridge and Lasso regression, we need a short reminder on statistics." + ] + }, + { + "cell_type": "markdown", + "id": "aa38d555", + "metadata": { + "editable": true + }, + "source": [ + "## More basic Statistics and Bayes' theorem\n", + "\n", + "A central theorem in statistics is Bayes' theorem. This theorem plays a similar role as the good old Pythagoras' theorem in geometry.\n", + "Bayes' theorem is extremely simple to derive. But to do so we need some basic axioms from statistics.\n", + "\n", + "Assume we have two domains of events $X=[x_0,x_1,\\dots,x_{n-1}]$ and $Y=[y_0,y_1,\\dots,y_{n-1}]$.\n", + "\n", + "We define also the likelihood for $X$ and $Y$ as $p(X)$ and $p(Y)$ respectively.\n", + "The likelihood of a specific event $x_i$ (or $y_i$) is then written as $p(X=x_i)$ or just $p(x_i)=p_i$. \n", + "\n", + "**Union of events is given by.**" + ] + }, + { + "cell_type": "markdown", + "id": "e9b2c3a4", + "metadata": { + "editable": true + }, + "source": [ + "$$\n", + "p(X \\cup Y)= p(X)+p(Y)-p(X \\cap Y).\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "68cb51e6", + "metadata": { + "editable": true + }, + "source": [ + "**The product rule (aka joint probability) is given by.**" + ] + }, + { + "cell_type": "markdown", + "id": "e0325fa4", + "metadata": { + "editable": true + }, + "source": [ + "$$\n", + "p(X \\cup Y)= p(X,Y)= p(X\\vert Y)p(Y)=p(Y\\vert X)p(X),\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "9e0da087", + "metadata": { + "editable": true + }, + "source": [ + "where we read $p(X\\vert Y)$ as the likelihood of obtaining $X$ given $Y$.\n", + "\n", + "If we have independent events then $p(X,Y)=p(X)p(Y)$." + ] + }, + { + "cell_type": "markdown", + "id": "6af44809", + "metadata": { + "editable": true + }, + "source": [ + "## Marginal Probability\n", + "\n", + "The marginal probability is defined in terms of only one of the set of variables $X,Y$. For a discrete probability we have" + ] + }, + { + "cell_type": "markdown", + "id": "71162259", + "metadata": { + "editable": true + }, + "source": [ + "$$\n", + "p(X)=\\sum_{i=0}^{n-1}p(X,Y=y_i)=\\sum_{i=0}^{n-1}p(X\\vert Y=y_i)p(Y=y_i)=\\sum_{i=0}^{n-1}p(X\\vert y_i)p(y_i).\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "622c0483", + "metadata": { + "editable": true + }, + "source": [ + "## Conditional Probability\n", + "\n", + "The conditional probability, if $p(Y) > 0$, is" + ] + }, + { + "cell_type": "markdown", + "id": "157eb8eb", + "metadata": { + "editable": true + }, + "source": [ + "$$\n", + "p(X\\vert Y)= \\frac{p(X,Y)}{p(Y)}=\\frac{p(X,Y)}{\\sum_{i=0}^{n-1}p(Y\\vert X=x_i)p(x_i)}.\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "a73c711d", + "metadata": { + "editable": true + }, + "source": [ + "## Bayes' Theorem\n", + "\n", + "If we combine the conditional probability with the marginal probability and the standard product rule, we have" + ] + }, + { + "cell_type": "markdown", + "id": "2d83fce3", + "metadata": { + "editable": true + }, + "source": [ + "$$\n", + "p(X\\vert Y)= \\frac{p(X,Y)}{p(Y)},\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "fa09b6e5", + "metadata": { + "editable": true + }, + "source": [ + "which we can rewrite as" + ] + }, + { + "cell_type": "markdown", + "id": "c21e7b34", + "metadata": { + "editable": true + }, + "source": [ + "$$\n", + "p(X\\vert Y)= \\frac{p(X,Y)}{\\sum_{i=0}^{n-1}p(Y\\vert X=x_i)p(x_i)}=\\frac{p(Y\\vert X)p(X)}{\\sum_{i=0}^{n-1}p(Y\\vert X=x_i)p(x_i)},\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "b0f6f2ec", + "metadata": { + "editable": true + }, + "source": [ + "which is Bayes' theorem. It allows us to evaluate the uncertainty in in $X$ after we have observed $Y$. We can easily interchange $X$ with $Y$." + ] + }, + { + "cell_type": "markdown", + "id": "fdb50d53", + "metadata": { + "editable": true + }, + "source": [ + "## Interpretations of Bayes' Theorem\n", + "\n", + "The quantity $p(Y\\vert X)$ on the right-hand side of the theorem is\n", + "evaluated for the observed data $Y$ and can be viewed as a function of\n", + "the parameter space represented by $X$. This function is not\n", + "necesseraly normalized and is normally called the likelihood function.\n", + "\n", + "The function $p(X)$ on the right hand side is called the prior while the function on the left hand side is the called the posterior probability. The denominator on the right hand side serves as a normalization factor for the posterior distribution.\n", + "\n", + "Let us try to illustrate Bayes' theorem through an example." + ] + }, + { + "cell_type": "markdown", + "id": "80ea6e5e", + "metadata": { + "editable": true + }, + "source": [ + "## Example of Usage of Bayes' theorem\n", + "\n", + "Let us suppose that you are undergoing a series of mammography scans in\n", + "order to rule out possible breast cancer cases. We define the\n", + "sensitivity for a positive event by the variable $X$. It takes binary\n", + "values with $X=1$ representing a positive event and $X=0$ being a\n", + "negative event. We reserve $Y$ as a classification parameter for\n", + "either a negative or a positive breast cancer confirmation. (Short note on wordings: positive here means having breast cancer, although none of us would consider this being a positive thing).\n", + "\n", + "We let $Y=1$ represent the the case of having breast cancer and $Y=0$ as not.\n", + "\n", + "Let us assume that if you have breast cancer, the test will be positive with a probability of $0.8$, that is we have" + ] + }, + { + "cell_type": "markdown", + "id": "bc467488", + "metadata": { + "editable": true + }, + "source": [ + "$$\n", + "p(X=1\\vert Y=1) =0.8.\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "e2f5ad4f", + "metadata": { + "editable": true + }, + "source": [ + "This obviously sounds scary since many would conclude that if the test is positive, there is a likelihood of $80\\%$ for having cancer.\n", + "It is however not correct, as the following Bayesian analysis shows." + ] + }, + { + "cell_type": "markdown", + "id": "488981c2", + "metadata": { + "editable": true + }, + "source": [ + "## Doing it correctly\n", + "\n", + "If we look at various national surveys on breast cancer, the general likelihood of developing breast cancer is a very small number.\n", + "Let us assume that the prior probability in the population as a whole is" + ] + }, + { + "cell_type": "markdown", + "id": "82e6140d", + "metadata": { + "editable": true + }, + "source": [ + "$$\n", + "p(Y=1) =0.004.\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "a4453800", + "metadata": { + "editable": true + }, + "source": [ + "We need also to account for the fact that the test may produce a false positive result (false alarm). Let us here assume that we have" + ] + }, + { + "cell_type": "markdown", + "id": "346f9f17", + "metadata": { + "editable": true + }, + "source": [ + "$$\n", + "p(X=1\\vert Y=0) =0.1.\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "a16b35ad", + "metadata": { + "editable": true + }, + "source": [ + "Using Bayes' theorem we can then find the posterior probability that the person has breast cancer in case of a positive test, that is we can compute" + ] + }, + { + "cell_type": "markdown", + "id": "4cd7eeb4", + "metadata": { + "editable": true + }, + "source": [ + "$$\n", + "p(Y=1\\vert X=1)=\\frac{p(X=1\\vert Y=1)p(Y=1)}{p(X=1\\vert Y=1)p(Y=1)+p(X=1\\vert Y=0)p(Y=0)}=\\frac{0.8\\times 0.004}{0.8\\times 0.004+0.1\\times 0.996}=0.031.\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "0bad6e84", + "metadata": { + "editable": true + }, + "source": [ + "That is, in case of a positive test, there is only a $3\\%$ chance of having breast cancer!" + ] + }, + { + "cell_type": "markdown", + "id": "de649672", + "metadata": { + "editable": true + }, + "source": [ + "## Bayes' Theorem and 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)" + ] + }, + { + "cell_type": "markdown", + "id": "247ef7ea", + "metadata": { + "editable": true + }, + "source": [ + "$$\n", + "\\boldsymbol{D}=[(x_0,y_0), (x_1,y_1),\\dots, (x_{n-1},y_{n-1})],\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "e668ae58", + "metadata": { + "editable": true + }, + "source": [ + "is given by" + ] + }, + { + "cell_type": "markdown", + "id": "503efa2d", + "metadata": { + "editable": true + }, + "source": [ + "$$\n", + "p(\\boldsymbol{D}\\vert\\boldsymbol{\\beta})=\\prod_{i=0}^{n-1}\\frac{1}{\\sqrt{2\\pi\\sigma^2}}\\exp{\\left[-\\frac{(y_i-\\boldsymbol{X}_{i,*}\\boldsymbol{\\beta})^2}{2\\sigma^2}\\right]}.\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "b2b875a0", "metadata": { "editable": true }, @@ -3991,7 +3980,7 @@ }, { "cell_type": "markdown", - "id": "32c55602", + "id": "55a767c5", "metadata": { "editable": true }, @@ -4003,7 +3992,7 @@ }, { "cell_type": "markdown", - "id": "9af380d7", + "id": "14bbfa0a", "metadata": { "editable": true }, @@ -4013,7 +4002,7 @@ }, { "cell_type": "markdown", - "id": "455577eb", + "id": "240ed0c3", "metadata": { "editable": true }, @@ -4025,7 +4014,7 @@ }, { "cell_type": "markdown", - "id": "bd3923f7", + "id": "f458005f", "metadata": { "editable": true }, @@ -4035,7 +4024,7 @@ }, { "cell_type": "markdown", - "id": "8e361210", + "id": "b8b753eb", "metadata": { "editable": true }, @@ -4051,7 +4040,7 @@ }, { "cell_type": "markdown", - "id": "19e3c297", + "id": "5e4534e0", "metadata": { "editable": true }, @@ -4063,7 +4052,7 @@ }, { "cell_type": "markdown", - "id": "3665066b", + "id": "c601b314", "metadata": { "editable": true }, @@ -4073,7 +4062,7 @@ }, { "cell_type": "markdown", - "id": "469f484f", + "id": "fe20511b", "metadata": { "editable": true }, @@ -4085,7 +4074,7 @@ }, { "cell_type": "markdown", - "id": "08f13d2f", + "id": "254593d7", "metadata": { "editable": true }, @@ -4098,7 +4087,7 @@ }, { "cell_type": "markdown", - "id": "2203fa0e", + "id": "a62f7507", "metadata": { "editable": true }, @@ -4110,7 +4099,7 @@ }, { "cell_type": "markdown", - "id": "01d952e5", + "id": "7c1a1880", "metadata": { "editable": true }, @@ -4120,7 +4109,7 @@ }, { "cell_type": "markdown", - "id": "46cbe057", + "id": "5719b77c", "metadata": { "editable": true }, @@ -4132,7 +4121,7 @@ }, { "cell_type": "markdown", - "id": "761243e1", + "id": "976134e0", "metadata": { "editable": true }, @@ -4142,7 +4131,7 @@ }, { "cell_type": "markdown", - "id": "02efcbd8", + "id": "ddc710a7", "metadata": { "editable": true }, @@ -4154,7 +4143,7 @@ }, { "cell_type": "markdown", - "id": "05414e8f", + "id": "9eba4148", "metadata": { "editable": true }, @@ -4166,7 +4155,7 @@ }, { "cell_type": "markdown", - "id": "fc664941", + "id": "1e028fb0", "metadata": { "editable": true }, @@ -4176,7 +4165,7 @@ }, { "cell_type": "markdown", - "id": "027a1818", + "id": "4ed079fb", "metadata": { "editable": true }, @@ -4188,7 +4177,7 @@ }, { "cell_type": "markdown", - "id": "7161d112", + "id": "879dc4c0", "metadata": { "editable": true }, @@ -4200,7 +4189,7 @@ }, { "cell_type": "markdown", - "id": "2f1fb6e3", + "id": "570d880d", "metadata": { "editable": true }, @@ -4212,7 +4201,7 @@ }, { "cell_type": "markdown", - "id": "54457246", + "id": "9b9e59f3", "metadata": { "editable": true }, @@ -4222,7 +4211,7 @@ }, { "cell_type": "markdown", - "id": "a3def508", + "id": "9bed769e", "metadata": { "editable": true }, @@ -4234,7 +4223,7 @@ }, { "cell_type": "markdown", - "id": "5769759a", + "id": "0d2b3bbb", "metadata": { "editable": true }, diff --git a/doc/src/week36/week36.do.txt b/doc/src/week36/week36.do.txt index b865938b3..d630ca7c8 100644 --- a/doc/src/week36/week36.do.txt +++ b/doc/src/week36/week36.do.txt @@ -1417,6 +1417,97 @@ well. It means that our optimization is now done only with the centered matrix and/or vector that enter the fitting procedure. + + +!split +===== Test Function for what happens with OLS, Ridge and Lasso ===== + +Hitherto we have discussed Ridge and Lasso regression in terms of a +linear analysis. This may to many of you feel rather technical and +perhaps not that intuitive. The question is whether we can develop a +more intuitive way of understanding what Ridge and Lasso express. + +Before we proceed let us perform a Ridge, Lasso and OLS analysis of a polynomial fit. + +We will play around with a study of the values for the optimal +parameters $\bm{\beta}$ using OLS, Ridge and Lasso regression. For +OLS, you will notice as function of the noise and polynomial degree, +that the parameters $\beta$ will fluctuate from order to order in the +polynomial fit and that for larger and larger polynomial degrees of freedom, the parameters will tend to increase in value for OLS. + +For Ridge and Lasso regression, the higher order parameters will typically be reduced, providing thereby less fluctuations from one order to another one. + +!bc pycod +import numpy as np +import matplotlib.pyplot as plt +from sklearn.model_selection import train_test_split +from sklearn import linear_model + +def R2(y_data, y_model): + return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2) +def MSE(y_data,y_model): + n = np.size(y_model) + return np.sum((y_data-y_model)**2)/n + +# Make data set. +n = 10000 +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): + X[:,polydegree] = x**(polydegree) + +# 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) + +# matrix inversion to find beta +OLSbeta = np.linalg.pinv(X_train.T @ X_train) @ X_train.T @ y_train +print(OLSbeta) +ypredictOLS = X_test @ OLSbeta +print("Test MSE OLS") +print(MSE(y_test,ypredictOLS)) +# Repeat now for Lasso and Ridge regression and various values of the regularization parameter using Scikit-Learn +# Decide which values of lambda to use +nlambdas = 4 +MSERidgePredict = np.zeros(nlambdas) +MSELassoPredict = np.zeros(nlambdas) +lambdas = np.logspace(-3, 1, nlambdas) +for i in range(nlambdas): + lmb = lambdas[i] + # Make the fit using Ridge and Lasso + RegRidge = linear_model.Ridge(lmb,fit_intercept=False) + RegRidge.fit(X_train,y_train) + RegLasso = linear_model.Lasso(lmb,fit_intercept=False) + RegLasso.fit(X_train,y_train) + # and then make the prediction + ypredictRidge = RegRidge.predict(X_test) + ypredictLasso = RegLasso.predict(X_test) + # Compute the MSE and print it + MSERidgePredict[i] = MSE(y_test,ypredictRidge) + MSELassoPredict[i] = MSE(y_test,ypredictLasso) + print(lmb,RegRidge.coef_) + print(lmb,RegLasso.coef_) +# Now plot the results +plt.figure() +plt.plot(np.log10(lambdas), MSERidgePredict, 'b', label = 'MSE Ridge Test') +plt.plot(np.log10(lambdas), MSELassoPredict, 'r', label = 'MSE Lasso Test') +plt.xlabel('log10(lambda)') +plt.ylabel('MSE') +plt.legend() +plt.show() + + +!ec + +How can we understand this? + + + !split ===== Linking the regression analysis with a statistical interpretation ===== @@ -1853,98 +1944,6 @@ That is, in case of a positive test, there is only a $3\%$ chance of having brea !split ===== Bayes' Theorem and Ridge and Lasso Regression ===== -Hitherto we have discussed Ridge and Lasso regression in terms of a -linear analysis. This may to many of you feel rather technical and -perhaps not that intuitive. The question is whether we can develop a -more intuitive way of understanding what Ridge and Lasso express. - -Before we proceed let us perform a Ridge, Lasso and OLS analysis of a polynomial fit. - -!split -===== Test Function for what happens with OLS, Ridge and Lasso ===== - -We will play around with a study of the values for the optimal -parameters $\bm{\beta}$ using OLS, Ridge and Lasso regression. For -OLS, you will notice as function of the noise and polynomial degree, -that the parameters $\beta$ will fluctuate from order to order in the -polynomial fit and that for larger and larger polynomial degrees of freedom, the parameters will tend to increase in value for OLS. - -For Ridge and Lasso regression, the higher order parameters will typically be reduced, providing thereby less fluctuations from one order to another one. - -!bc pycod -import numpy as np -import matplotlib.pyplot as plt -from sklearn.model_selection import train_test_split -from sklearn import linear_model - -def R2(y_data, y_model): - return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2) -def MSE(y_data,y_model): - n = np.size(y_model) - return np.sum((y_data-y_model)**2)/n - -# Make data set. -n = 10000 -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) - -# matrix inversion to find beta -OLSbeta = np.linalg.pinv(X_train.T @ X_train) @ X_train.T @ y_train -print(OLSbeta) -ypredictOLS = X_test @ OLSbeta -print("Test MSE OLS") -print(MSE(y_test,ypredictOLS)) -# Repeat now for Lasso and Ridge regression and various values of the regularization parameter using Scikit-Learn -# Decide which values of lambda to use -nlambdas = 4 -MSERidgePredict = np.zeros(nlambdas) -MSELassoPredict = np.zeros(nlambdas) -lambdas = np.logspace(-3, 1, nlambdas) -for i in range(nlambdas): - lmb = lambdas[i] - # Make the fit using Ridge and Lasso - RegRidge = linear_model.Ridge(lmb,fit_intercept=False) - RegRidge.fit(X_train,y_train) - RegLasso = linear_model.Lasso(lmb,fit_intercept=False) - RegLasso.fit(X_train,y_train) - # and then make the prediction - ypredictRidge = RegRidge.predict(X_test) - ypredictLasso = RegLasso.predict(X_test) - # Compute the MSE and print it - MSERidgePredict[i] = MSE(y_test,ypredictRidge) - MSELassoPredict[i] = MSE(y_test,ypredictLasso) - print(lmb,RegRidge.coef_) - print(lmb,RegLasso.coef_) -# Now plot the results -plt.figure() -plt.plot(np.log10(lambdas), MSERidgePredict, 'b', label = 'MSE Ridge Test') -plt.plot(np.log10(lambdas), MSELassoPredict, 'r', label = 'MSE Lasso Test') -plt.xlabel('log10(lambda)') -plt.ylabel('MSE') -plt.legend() -plt.show() - - -!ec - -How can we understand this? - - -!split -===== Invoking Bayes' theorem ===== - 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)