-
@@ -438,7 +452,7 @@ MathJax.Hub.Config({
-Yes, it could be a bad idea to include the intercept column for the -exact reason you stated. If no transformation is applied to your data, -the intercept can be interpreted as the expected value of your target -variable when all your predictors are put to zero. Therefore, whenever -you cannot assume that the expected target variable is zero when all -your predictors are zero, it could be a bad idea to apply a model -which penalizes the intercept. Also, the analytical solution to the -ridge regression coefficients (when not shrinking $$\beta_0$$) is -derived under the assumption that both y and X are zero centered (mean -subtracted). What you are doing is correct, but you should also zero -center X (subtracting the mean of each column from the corresponding -column). +The intercept can be interpreted as the expected value of our +target/output variables when all other predictors are set to zero. +Thus, if we cannot assume that the expected outputs/targets are zero +when all predictors are zero (the columns in the design matrix), it +may be a bad idea to implement a model which penalizes the intercept. +Furthermore, in for example Ridge and Lasso regression, the solutions +(when not shrinking $$\beta_0$$) for the unknown parameters +\( \boldsymbol{\beta} \) are derived under the assumption that both \( \boldsymbol{y} \) and +\( \boldsymbol{X} \) are zero centered, that is we subtract the mean values. -
-If your predictors are of different scales, I would advice you to -standardize X by subtracting the mean of each column from the -corresponding column and dividing the column with its standard -deviation. If you dont do this, you will give an "unfair" penalization -of the parameters since their magnitude depends on the scale of their -corresponding predictor. Suppose that you have an input variable -"height". Human height might be measured in inches or meters or -kilometers. If measured in kilometers, a standard linear regression -model with this predictor would probably give a much bigger -coefficient term, than if measured in millimeters. You may see how -this could become a problem when considering the loss function for -ridge regression. - -
-Remember that when you do any transformation to your dataset before -training, the exact same transformation has to be applied to new data -before making a prediction. In your case, this means: - -
- - -
#Model training:
-y_train_mean = np.mean(y_train)
-X_train_mean = np.mean(X_train,axis=0)
-X_train = X_train - X_train_mean
-y_train = y_train - y_train_mean
-
-trained_model = some_model.fit(X_train,y_train)
-
-#Model prediction:
-X_test = X_test - X_train_mean #Use mean from training data
-y_pred = trained_model(X_test)
-y_pred = y_pred + y_train_mean
--Here is a mathematical explanation of the zero centering: - -
-The cost/loss function for Ridge regression is: - -$$ -C(\beta_0, \beta_1, ... , \beta_P) = \sum_{i=1}^{n} (y_i - \beta_0 - \sum_{p=1}^P X_{ip}\beta_p)^2 + \lambda \sum_{p=1}^P \beta_p^2. -$$ - -
-Notice that the intercept is left out of the \( L_2 \) regularization term. The design matrix -\( X \) does in this case not contain any intercept column. We want - -$$ -\frac{\partial L}{\partial \beta_j} = 0, -$$ - -
-for all \( j \), so lets start with \( \beta_0 \). This means that we have - -$$ -\frac{\partial L}{\partial \beta_0} = -2\sum_{i=1}^{n} (y_i - \beta_0 - \sum_{p=1}^P X_{ip} \beta_p). -$$ - -
-We want to solve -$$ --2\sum_{i=1}^{n} (y_i - \beta_0 - \sum_{p=1}^P X_{ip} \beta_p) = 0, -$$ - -
-which gives -$$ -\sum_{i=1}^{n} \beta_0 = \sum_{i=1}^{n}y_i - \sum_{i=1}^{n} \sum_{p=1}^P X_{ip} \beta_p, -$$ - -
-or -$ n\beta_0 = \sum_{i=1}^{n} y_i - \sum_{p=1}^P\beta_p \sum_{i=1}^{n} X_{ip}$. - -
-If we assume that every column of \( X \) is centered, whic we can do by subtracting the mean, -
- - -
X = X - np.mean(X,axis=0)
--the sum $ \sum_{i=1}^{n} X_{ip} $ - -
-can be rewritten as -$$ -\sum_{i=1}^{n} (X_{ip} - \frac{1}{n}\sum_{i=1}^{n} X_{ip}) = \sum_{i=1}^{n} X_{ip} - \sum_{i=1}^{n} \frac{1}{n} \sum_{i=1}^{n}X_{ip}, -$$ - -resulting in -$$ -\sum_{i=1}^{n} X_{ip} - n \frac{1}{n} \sum_{i=1}^{n}X_{ip} = 0. -$$ - -
-Finally we have -$$ -n\beta_0 = \sum_{i=1}^{n} y_i - \sum_{p=1}^P\beta_p \sum_{i=1}^{n} X_{ip}, -$$ - -or -$$ -\beta_0 = \frac{1}{n}\sum_{i=1}^{n} y_i = y_{average}. -$$ - -
-Replacing \( y_i \) with \( y_i - \beta_0 = y_i - y_{average} \) in the loss function will give us (in vector-matrix disguise) -$$ -C(\boldsymbol{\beta}) = (\boldsymbol{\tilde{y}} - \tilde{X}\boldsymbol{\beta})^T(\boldsymbol{\tilde{y}} - \tilde{X}\boldsymbol{\beta}) + \lambda \boldsymbol{\beta}^T\boldsymbol{\beta}, -$$ - -
-which has the solution - -
-\( \beta = (\tilde{X}^T\tilde{X} + \lambda I)^{-1}\tilde{X}^T\boldsymbol{\tilde{y}} \). -where \( \boldsymbol{\tilde{y}} = \boldsymbol{y} - y_{average} \) -and \( \tilde{X}_{ij} = X_{ij} - \frac{1}{n}\sum_{k=1}^{n-1}X_{kj} \). - -
- - -
import numpy as np
-import pandas as pd
-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
-
-
-# A seed just to ensure that the random numbers are the same for every run.
-# Useful for eventual debugging.
-np.random.seed(3155)
-
-n = 100
-x = np.random.rand(n)
-y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)
-
-Maxpolydegree = 20
-X = np.zeros((n,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)
-# and then make the prediction
-ytildeOLS = X_train @ OLSbeta
-print("Training MSE for OLS")
-print(MSE(y_train,ytildeOLS))
-ypredictOLS = X_test @ OLSbeta
-print("Test MSE OLS")
-print(MSE(y_test,ypredictOLS))
-
-p = len(OLSbeta)
-I = np.eye(p,p)
-# Decide which values of lambda to use
-nlambdas = 4
-MSEOwnRidgePredict = np.zeros(nlambdas)
-MSEOwnRidgeTrain = np.zeros(nlambdas)
-MSERidgePredict = np.zeros(nlambdas)
-MSERidgeTrain = np.zeros(nlambdas)
-
-lambdas = np.logspace(-4, 4, nlambdas)
-for i in range(nlambdas):
- lmb = lambdas[i]
- OwnRidgeBeta = np.linalg.pinv(X_train.T @ X_train+lmb*I) @ X_train.T @ y_train
- # include lasso using Scikit-Learn
- # Note: we include the intercept
- RegRidge = linear_model.Ridge(lmb,fit_intercept=False)
- RegRidge.fit(X_train,y_train)
- # and then make the prediction
- ytildeOwnRidge = X_train @ OwnRidgeBeta
- ypredictOwnRidge = X_test @ OwnRidgeBeta
- ytildeRidge = RegRidge.predict(X_train)
- ypredictRidge = RegRidge.predict(X_test)
- MSEOwnRidgePredict[i] = MSE(y_test,ypredictOwnRidge)
- MSEOwnRidgeTrain[i] = MSE(y_train,ytildeOwnRidge)
- MSERidgePredict[i] = MSE(y_test,ypredictRidge)
- MSERidgeTrain[i] = MSE(y_train,ytildeRidge)
- print("Beta values for own Ridge implementation")
- print(OwnRidgeBeta)
- print("Beta values for Scikit-Learn Ridge implementation")
- print(RegRidge.coef_)
-# Now plot the results
-plt.figure()
-plt.plot(np.log10(lambdas), MSEOwnRidgeTrain, 'b', label = 'MSE Ridge train')
-plt.plot(np.log10(lambdas), MSEOwnRidgePredict, 'r', label = 'MSE Ridge Test')
-plt.plot(np.log10(lambdas), MSERidgeTrain, 'y', label = 'MSE Ridge train')
-plt.plot(np.log10(lambdas), MSERidgePredict, 'g', label = 'MSE Ridge Test')
-
-plt.xlabel('log10(lambda)')
-plt.ylabel('MSE')
-plt.legend()
-plt.show()
-- - -
import numpy as np
-import pandas as pd
-import matplotlib.pyplot as plt
-from sklearn.model_selection import train_test_split
-from sklearn import linear_model
-from sklearn.preprocessing import StandardScaler
-
-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
-
-
-# A seed just to ensure that the random numbers are the same for every run.
-# Useful for eventual debugging.
-np.random.seed(315)
-
-n = 100
-x = np.random.rand(n)
-y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)
-
-Maxpolydegree = 5
-X = np.zeros((n,Maxpolydegree-1))
-
-for degree in range(1,Maxpolydegree): #No intercept column
- X[:,degree-1] = 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)
-
-
-
-
-
-#For our own implementation, we will need to deal with the intercept by centering the design matrix and the target variable
-X_train_mean = np.mean(X_train,axis=0)
-X_train_scaled = X_train - X_train_mean #Center by removing mean from each feature
-X_test_scaled = X_test - X_train_mean
-
-y_scaler = np.mean(y_train) #The model intercept (called y_scaler) is given by the mean of target variable (IF X is centered)
-y_train_scaled = y_train - y_scaler #Remove the intercept from the training data.
-
-
-p = Maxpolydegree-1
-I = np.eye(p,p)
-# Decide which values of lambda to use
-nlambdas = 4
-MSEOwnRidgePredict = np.zeros(nlambdas)
-MSERidgePredict = np.zeros(nlambdas)
-
-lambdas = np.logspace(-4, 1, nlambdas)
-for i in range(nlambdas):
- lmb = lambdas[i]
- OwnRidgeBeta = np.linalg.pinv(X_train_scaled.T @ X_train_scaled+lmb*I) @ X_train_scaled.T @ (y_train_scaled)
- intercept_ = y_scaler - X_train_mean@OwnRidgeBeta #The intercept can be shifted so the model can predict on uncentered data
-
- ypredictOwnRidge = X_test @ OwnRidgeBeta + intercept_ #Add intercept to prediction
- #EQUIVALENT PREDICTION:
- ypredictOwnRidge = X_test_scaled @ OwnRidgeBeta + y_scaler #Add intercept to prediction
- print("Values for own Ridge prediction")
- print(ypredictOwnRidge)
-
-
-
- RegRidge = linear_model.Ridge(lmb)
- RegRidge.fit(X_train,y_train)
- ypredictRidge = RegRidge.predict(X_test)
- print("Values for SL Ridge prediction")
- print(ypredictRidge)
-
-
- MSEOwnRidgePredict[i] = MSE(y_test,ypredictOwnRidge)
- MSERidgePredict[i] = MSE(y_test,ypredictRidge)
-
- print("Beta values for own Ridge implementation")
- print(OwnRidgeBeta) #Intercept is given by mean of target variable
- print("Beta values for Scikit-Learn Ridge implementation")
- print(RegRidge.coef_)
- print('Intercept from own implementation:')
- print(intercept_)
- print('Intercept from Scikit-Learn Ridge implementation')
- print(RegRidge.intercept_)
-
-
-
-# Now plot the results
-
-plt.figure()
-plt.plot(np.log10(lambdas), MSEOwnRidgePredict, 'b--', label = 'MSE own Ridge Test')
-plt.plot(np.log10(lambdas), MSERidgePredict, 'g--', label = 'MSE SL Ridge Test')
-
-plt.xlabel('log10(lambda)')
-plt.ylabel('MSE')
-plt.legend()
-plt.show()
-- - -
import numpy as np
-import matplotlib.pyplot as plt
-
-from sklearn.linear_model import LinearRegression
-
-
-np.random.seed(2021)
-
-
-def fit_beta(X, y):
- return np.linalg.pinv(X.T @ X) @ X.T @ y
-
-
-true_beta = [2, 0.5, 3.7]
-
-x = np.linspace(0, 1, 11)
-y = np.sum(
- np.asarray([x ** p * b for p, b in enumerate(true_beta)]), axis=0
-) + 0.1 * np.random.normal(size=len(x))
-
-degree = 3
-X = np.zeros((len(x), degree))
-
-# Include the intercept in the design matrix
-for p in range(degree):
- X[:, p] = x ** p
-
-beta = fit_beta(X, y)
-
-# Intercept is included in the design matrix
-clf = LinearRegression(fit_intercept=False).fit(X, y)
-
-print(f"True beta: {true_beta}")
-print(f"Fitted beta: {beta}")
-print(f"Sklearn fitted beta: {clf.coef_}")
-
-
-plt.figure()
-plt.scatter(x, y, label="Data")
-plt.plot(x, X @ beta, label="Fit")
-plt.plot(x, clf.predict(X), label="Sklearn (fit_intercept=False)")
-
-
-# Do not include the intercept in the design matrix
-X = np.zeros((len(x), degree - 1))
-
-for p in range(degree - 1):
- X[:, p] = x ** (p + 1)
-
-# Intercept is not included in the design matrix
-clf = LinearRegression(fit_intercept=True).fit(X, y)
-
-# Use centered values for X and y when computing coefficients
-y_offset = np.average(y, axis=0)
-X_offset = np.average(X, axis=0)
-
-beta = fit_beta(X - X_offset, y - y_offset)
-intercept = np.mean(y_offset - X_offset @ beta)
-
-print(f"Manual intercept: {intercept}")
-print(f"Fitted beta (sans intercept): {beta}")
-print(f"Sklearn intercept: {clf.intercept_}")
-print(f"Sklearn fitted beta (sans intercept): {clf.coef_}")
-
-plt.plot(x, X @ beta + intercept, "--", label="Fit (manual intercept)")
-plt.plot(x, clf.predict(X), "--", label="Sklearn (fit_intercept=True)")
-plt.grid()
-plt.legend()
-
-plt.show()
-
@@ -833,7 +456,7 @@ plt.show()
-The one-dimensional Ising model with nearest neighbor interaction, no -external field and a constant coupling constant \( J \) is given by - -$$ -\begin{align} - H = -J \sum_{k}^L s_k s_{k + 1}, -\tag{1} -\end{align} -$$ +If our predictors represent different scales, then it is important to +standardize the design matrix \( \boldsymbol{X} \) by subtracting the mean of each +column from the corresponding column and dividing the column with its +standard deviation.
-where \( s_i \in \{-1, 1\} \) and \( s_{N + 1} = s_1 \). The number of spins -in the system is determined by \( L \). For the one-dimensional system -there is no phase transition. +The +Standadscaler +function in Scikit-Learn does this for us. For the data sets we +have been studying in our various examples, the data are in many cases +already scaled and there is no need to scale them.
-We will look at a system of \( L = 40 \) spins with a coupling constant of -\( J = 1 \). To get enough training data we will generate 10000 states -with their respective energies. +If you need to scale the data, not doing so will give an unfair +penalization of the parameters since their magnitude depends on the +scale of their corresponding predictor.
- - -
import numpy as np
-import matplotlib.pyplot as plt
-from mpl_toolkits.axes_grid1 import make_axes_locatable
-import seaborn as sns
-import scipy.linalg as scl
-from sklearn.model_selection import train_test_split
-import tqdm
-sns.set(color_codes=True)
-cmap_args=dict(vmin=-1., vmax=1., cmap='seismic')
-
-L = 40
-n = int(1e4)
-
-spins = np.random.choice([-1, 1], size=(n, L))
-J = 1.0
-
-energies = np.zeros(n)
-
-for i in range(n):
- energies[i] = - J * np.dot(spins[i], np.roll(spins[i], 1))
--Here we use ordinary least squares -regression to predict the energy for the nearest neighbor -one-dimensional Ising model on a ring, i.e., the endpoints wrap -around. We will use linear regression to fit a value for -the coupling constant to achieve this. +Suppose as an example that you +you have an input variable given by the heights of different persons. +Human height might be measured in inches or meters or +kilometers. If measured in kilometers, a standard linear regression +model with this predictor would probably give a much bigger +coefficient term, than if measured in millimeters. +This can clearly lead to problems in evaluating the cost/loss functions.
@@ -474,7 +463,7 @@ the coupling constant to achieve this.
-A more general form for the one-dimensional Ising model is - -$$ -\begin{align} - H = - \sum_j^L \sum_k^L s_j s_k J_{jk}. -\tag{2} -\end{align} -$$ - -
-Here we allow for interactions beyond the nearest neighbors and a state dependent -coupling constant. This latter expression can be formulated as -a matrix-product -$$ -\begin{align} - \boldsymbol{H} = \boldsymbol{X} J, -\tag{3} -\end{align} -$$ - -
-where \( X_{jk} = s_j s_k \) and \( J \) is a matrix which consists of the -elements \( -J_{jk} \). This form of writing the energy fits perfectly -with the form utilized in linear regression, that is - -$$ -\begin{align} - \boldsymbol{y} = \boldsymbol{X}\boldsymbol{\beta} + \boldsymbol{\epsilon}, -\tag{4} -\end{align} -$$ - -
-We split the data in training and test data as discussed in the previous example +Keep in mind that when you transform your data set before training a model, the same transformation needs to be done +on your eventual new data set before making a prediction. If we translate this into a Python code, it would could be implemented as follows
-
X = np.zeros((n, L ** 2))
-for i in range(n):
- X[i] = np.outer(spins[i], spins[i]).ravel()
-y = energies
-X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)
+#Model training, we compute the mean value of y and X
+y_train_mean = np.mean(y_train)
+X_train_mean = np.mean(X_train,axis=0)
+X_train = X_train - X_train_mean
+y_train = y_train - y_train_mean
+
+# The we fit our model with the training data
+trained_model = some_model.fit(X_train,y_train)
+
+
+#Model prediction, here we need also to transform our data set used for the prediction.
+X_test = X_test - X_train_mean #Use mean from training data
+y_pred = trained_model(X_test)
+y_pred = y_pred + y_train_mean
@@ -468,7 +459,7 @@ X_train, X_test, y_train, y_test = train_tes
-In the ordinary least squares method we choose the cost function - -$$ -\begin{align} - C(\boldsymbol{X}, \boldsymbol{\beta})= \frac{1}{n}\left\{(\boldsymbol{X}\boldsymbol{\beta} - \boldsymbol{y})^T(\boldsymbol{X}\boldsymbol{\beta} - \boldsymbol{y})\right\}. -\tag{5} -\end{align} -$$ - -
-We then find the extremal point of \( C \) by taking the derivative with respect to \( \boldsymbol{\beta} \) as discussed above. -This yields the expression for \( \boldsymbol{\beta} \) to be - -$$ - \boldsymbol{\beta} = \frac{\boldsymbol{X}^T \boldsymbol{y}}{\boldsymbol{X}^T \boldsymbol{X}}, -$$ - -
-which immediately imposes some requirements on \( \boldsymbol{X} \) as there must exist -an inverse of \( \boldsymbol{X}^T \boldsymbol{X} \). If the expression we are modeling contains an -intercept, i.e., a constant term, we must make sure that the -first column of \( \boldsymbol{X} \) consists of \( 1 \). We do this here +This code shows a simple first-order fit to a data set using the above transformed data, where we consider the role of the intercept first, by either excluding it or including it (code example thanks to Øyvind Sigmundson Schøyen)
-
X_train_own = np.concatenate(
- (np.ones(len(X_train))[:, np.newaxis], X_train),
- axis=1
-)
-X_test_own = np.concatenate(
- (np.ones(len(X_test))[:, np.newaxis], X_test),
- axis=1
-)
-+
import numpy as np
+import matplotlib.pyplot as plt
-
-def ols_inv(x: np.ndarray, y: np.ndarray) -> np.ndarray:
- return scl.inv(x.T @ x) @ (x.T @ y)
-beta = ols_inv(X_train_own, y_train)
+from sklearn.linear_model import LinearRegression
+
+
+np.random.seed(2021)
+
+
+def fit_beta(X, y):
+ return np.linalg.pinv(X.T @ X) @ X.T @ y
+
+
+true_beta = [2, 0.5, 3.7]
+
+x = np.linspace(0, 1, 11)
+y = np.sum(
+ np.asarray([x ** p * b for p, b in enumerate(true_beta)]), axis=0
+) + 0.1 * np.random.normal(size=len(x))
+
+degree = 3
+X = np.zeros((len(x), degree))
+
+# Include the intercept in the design matrix
+for p in range(degree):
+ X[:, p] = x ** p
+
+beta = fit_beta(X, y)
+
+# Intercept is included in the design matrix
+clf = LinearRegression(fit_intercept=False).fit(X, y)
+
+print(f"True beta: {true_beta}")
+print(f"Fitted beta: {beta}")
+print(f"Sklearn fitted beta: {clf.coef_}")
+
+
+plt.figure()
+plt.scatter(x, y, label="Data")
+plt.plot(x, X @ beta, label="Fit")
+plt.plot(x, clf.predict(X), label="Sklearn (fit_intercept=False)")
+
+
+# Do not include the intercept in the design matrix
+X = np.zeros((len(x), degree - 1))
+
+for p in range(degree - 1):
+ X[:, p] = x ** (p + 1)
+
+# Intercept is not included in the design matrix
+clf = LinearRegression(fit_intercept=True).fit(X, y)
+
+# Use centered values for X and y when computing coefficients
+y_offset = np.average(y, axis=0)
+X_offset = np.average(X, axis=0)
+
+beta = fit_beta(X - X_offset, y - y_offset)
+intercept = np.mean(y_offset - X_offset @ beta)
+
+print(f"Manual intercept: {intercept}")
+print(f"Fitted beta (sans intercept): {beta}")
+print(f"Sklearn intercept: {clf.intercept_}")
+print(f"Sklearn fitted beta (sans intercept): {clf.coef_}")
+
+plt.plot(x, X @ beta + intercept, "--", label="Fit (manual intercept)")
+plt.plot(x, clf.predict(X), "--", label="Sklearn (fit_intercept=True)")
+plt.grid()
+plt.legend()
+
+plt.show()
@@ -466,7 +514,7 @@ beta = ols_inv(X_train_own, y_train)
-Doing the inversion directly turns out to be a bad idea since the matrix -\( \boldsymbol{X}^T\boldsymbol{X} \) is singular. An alternative approach is to use the singular -value decomposition. Using the definition of the Moore-Penrose -pseudoinverse we can write the equation for \( \boldsymbol{\beta} \) as +The cost/loss function for Ridge regression is: $$ - \boldsymbol{\beta} = \boldsymbol{X}^{+}\boldsymbol{y}, +C(\beta_0, \beta_1, ... , \beta_P) = \sum_{i=1}^{n} (y_i - \beta_0 - \sum_{p=1}^P X_{ip}\beta_p)^2 + \lambda \sum_{p=1}^P \beta_p^2. $$
-where the pseudoinverse of \( \boldsymbol{X} \) is given by +Notice that the intercept is left out of the \( L_2 \) regularization term. The design matrix +\( X \) does in this case not contain any intercept column. We want $$ - \boldsymbol{X}^{+} = \frac{\boldsymbol{X}^T}{\boldsymbol{X}^T\boldsymbol{X}}. +\frac{\partial L}{\partial \beta_j} = 0, $$
-Using singular value decomposition we can decompose the matrix \( \boldsymbol{X} = \boldsymbol{U}\boldsymbol{\Sigma} \boldsymbol{V}^T \), -where \( \boldsymbol{U} \) and \( \boldsymbol{V} \) are orthogonal(unitary) matrices and \( \boldsymbol{\Sigma} \) contains the singular values (more details below). -where \( X^{+} = V\Sigma^{+} U^T \). This reduces the equation for -\( \omega \) to +for all \( j \), so lets start with \( \beta_0 \). This means that we have + $$ -\begin{align} - \boldsymbol{\beta} = \boldsymbol{V}\boldsymbol{\Sigma}^{+} \boldsymbol{U}^T \boldsymbol{y}. -\tag{6} -\end{align} +\frac{\partial L}{\partial \beta_0} = -2\sum_{i=1}^{n} (y_i - \beta_0 - \sum_{p=1}^P X_{ip} \beta_p). $$
-Note that solving this equation by actually doing the pseudoinverse -(which is what we will do) is not a good idea as this operation scales -as \( \mathcal{O}(n^3) \), where \( n \) is the number of elements in a -general matrix. Instead, doing \( QR \)-factorization and solving the -linear system as an equation would reduce this down to -\( \mathcal{O}(n^2) \) operations. +We want to solve +$$ +-2\sum_{i=1}^{n} (y_i - \beta_0 - \sum_{p=1}^P X_{ip} \beta_p) = 0, +$$ +
+which gives +$$ +\sum_{i=1}^{n} \beta_0 = \sum_{i=1}^{n}y_i - \sum_{i=1}^{n} \sum_{p=1}^P X_{ip} \beta_p, +$$ + +
+or +$ n\beta_0 = \sum_{i=1}^{n} y_i - \sum_{p=1}^P\beta_p \sum_{i=1}^{n} X_{ip}$. + +
+If we assume that every column of \( X \) is centered, whic we can do by subtracting the mean,
-
def ols_svd(x: np.ndarray, y: np.ndarray) -> np.ndarray:
- u, s, v = scl.svd(x)
- return v.T @ scl.pinv(scl.diagsvd(s, u.shape[0], v.shape[0])) @ u.T @ y
+X = X - np.mean(X,axis=0)
+the sum $ \sum_{i=1}^{n} X_{ip} $
-
-
beta = ols_svd(X_train_own,y_train)
-
-When extracting the \( J \)-matrix we need to make sure that we remove the intercept, as is done here
+can be rewritten as
+$$
+\sum_{i=1}^{n} (X_{ip} - \frac{1}{n}\sum_{i=1}^{n} X_{ip}) = \sum_{i=1}^{n} X_{ip} - \sum_{i=1}^{n} \frac{1}{n} \sum_{i=1}^{n}X_{ip},
+$$
+
+resulting in
+$$
+\sum_{i=1}^{n} X_{ip} - n \frac{1}{n} \sum_{i=1}^{n}X_{ip} = 0.
+$$
+
+
+Finally we have
+$$
+n\beta_0 = \sum_{i=1}^{n} y_i - \sum_{p=1}^P\beta_p \sum_{i=1}^{n} X_{ip},
+$$
+
+or
+$$
+\beta_0 = \frac{1}{n}\sum_{i=1}^{n} y_i = y_{average}.
+$$
+
+
+Replacing \( y_i \) with \( y_i - \beta_0 = y_i - y_{average} \) in the loss function will give us (in vector-matrix disguise)
+$$
+C(\boldsymbol{\beta}) = (\boldsymbol{\tilde{y}} - \tilde{X}\boldsymbol{\beta})^T(\boldsymbol{\tilde{y}} - \tilde{X}\boldsymbol{\beta}) + \lambda \boldsymbol{\beta}^T\boldsymbol{\beta},
+$$
+
+
+which has the solution
+
+
+\( \beta = (\tilde{X}^T\tilde{X} + \lambda I)^{-1}\tilde{X}^T\boldsymbol{\tilde{y}} \).
+where \( \boldsymbol{\tilde{y}} = \boldsymbol{y} - y_{average} \)
+and \( \tilde{X}_{ij} = X_{ij} - \frac{1}{n}\sum_{k=1}^{n-1}X_{kj} \).
-
J = beta[1:].reshape(L, L)
-
-
-A way of looking at the coefficients in \( J \) is to plot the matrices as images.
+
import numpy as np
+import pandas as pd
+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
-
-
fig = plt.figure(figsize=(20, 14))
-im = plt.imshow(J, **cmap_args)
-plt.title("OLS", fontsize=18)
-plt.xticks(fontsize=18)
-plt.yticks(fontsize=18)
-cb = fig.colorbar(im)
-cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=18)
+
+# A seed just to ensure that the random numbers are the same for every run.
+# Useful for eventual debugging.
+np.random.seed(3155)
+
+n = 100
+x = np.random.rand(n)
+y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)
+
+Maxpolydegree = 20
+X = np.zeros((n,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)
+# and then make the prediction
+ytildeOLS = X_train @ OLSbeta
+print("Training MSE for OLS")
+print(MSE(y_train,ytildeOLS))
+ypredictOLS = X_test @ OLSbeta
+print("Test MSE OLS")
+print(MSE(y_test,ypredictOLS))
+
+p = len(OLSbeta)
+I = np.eye(p,p)
+# Decide which values of lambda to use
+nlambdas = 4
+MSEOwnRidgePredict = np.zeros(nlambdas)
+MSEOwnRidgeTrain = np.zeros(nlambdas)
+MSERidgePredict = np.zeros(nlambdas)
+MSERidgeTrain = np.zeros(nlambdas)
+
+lambdas = np.logspace(-4, 4, nlambdas)
+for i in range(nlambdas):
+ lmb = lambdas[i]
+ OwnRidgeBeta = np.linalg.pinv(X_train.T @ X_train+lmb*I) @ X_train.T @ y_train
+ # include lasso using Scikit-Learn
+ # Note: we include the intercept
+ RegRidge = linear_model.Ridge(lmb,fit_intercept=False)
+ RegRidge.fit(X_train,y_train)
+ # and then make the prediction
+ ytildeOwnRidge = X_train @ OwnRidgeBeta
+ ypredictOwnRidge = X_test @ OwnRidgeBeta
+ ytildeRidge = RegRidge.predict(X_train)
+ ypredictRidge = RegRidge.predict(X_test)
+ MSEOwnRidgePredict[i] = MSE(y_test,ypredictOwnRidge)
+ MSEOwnRidgeTrain[i] = MSE(y_train,ytildeOwnRidge)
+ MSERidgePredict[i] = MSE(y_test,ypredictRidge)
+ MSERidgeTrain[i] = MSE(y_train,ytildeRidge)
+ print("Beta values for own Ridge implementation")
+ print(OwnRidgeBeta)
+ print("Beta values for Scikit-Learn Ridge implementation")
+ print(RegRidge.coef_)
+# Now plot the results
+plt.figure()
+plt.plot(np.log10(lambdas), MSEOwnRidgeTrain, 'b', label = 'MSE Ridge train')
+plt.plot(np.log10(lambdas), MSEOwnRidgePredict, 'r', label = 'MSE Ridge Test')
+plt.plot(np.log10(lambdas), MSERidgeTrain, 'y', label = 'MSE Ridge train')
+plt.plot(np.log10(lambdas), MSERidgePredict, 'g', label = 'MSE Ridge Test')
+
+plt.xlabel('log10(lambda)')
+plt.ylabel('MSE')
+plt.legend()
plt.show()
-It is interesting to note that OLS
-considers both \( J_{j, j + 1} = -0.5 \) and \( J_{j, j - 1} = -0.5 \) as
-valid matrix elements for \( J \).
-In our discussion below on hyperparameters and Ridge and Lasso regression we will see that
-this problem can be removed, partly and only with Lasso regression.
-
-In this case our matrix inversion was actually possible. The obvious question now is what is the mathematics behind the SVD?
+
+
import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+from sklearn.model_selection import train_test_split
+from sklearn import linear_model
+from sklearn.preprocessing import StandardScaler
+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
+
+
+# A seed just to ensure that the random numbers are the same for every run.
+# Useful for eventual debugging.
+np.random.seed(315)
+
+n = 100
+x = np.random.rand(n)
+y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)
+
+Maxpolydegree = 5
+X = np.zeros((n,Maxpolydegree-1))
+
+for degree in range(1,Maxpolydegree): #No intercept column
+ X[:,degree-1] = 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)
+
+
+
+
+
+#For our own implementation, we will need to deal with the intercept by centering the design matrix and the target variable
+X_train_mean = np.mean(X_train,axis=0)
+X_train_scaled = X_train - X_train_mean #Center by removing mean from each feature
+X_test_scaled = X_test - X_train_mean
+
+y_scaler = np.mean(y_train) #The model intercept (called y_scaler) is given by the mean of target variable (IF X is centered)
+y_train_scaled = y_train - y_scaler #Remove the intercept from the training data.
+
+
+p = Maxpolydegree-1
+I = np.eye(p,p)
+# Decide which values of lambda to use
+nlambdas = 4
+MSEOwnRidgePredict = np.zeros(nlambdas)
+MSERidgePredict = np.zeros(nlambdas)
+
+lambdas = np.logspace(-4, 1, nlambdas)
+for i in range(nlambdas):
+ lmb = lambdas[i]
+ OwnRidgeBeta = np.linalg.pinv(X_train_scaled.T @ X_train_scaled+lmb*I) @ X_train_scaled.T @ (y_train_scaled)
+ intercept_ = y_scaler - X_train_mean@OwnRidgeBeta #The intercept can be shifted so the model can predict on uncentered data
+
+ ypredictOwnRidge = X_test @ OwnRidgeBeta + intercept_ #Add intercept to prediction
+ #EQUIVALENT PREDICTION:
+ ypredictOwnRidge = X_test_scaled @ OwnRidgeBeta + y_scaler #Add intercept to prediction
+ print("Values for own Ridge prediction")
+ print(ypredictOwnRidge)
+
+
+
+ RegRidge = linear_model.Ridge(lmb)
+ RegRidge.fit(X_train,y_train)
+ ypredictRidge = RegRidge.predict(X_test)
+ print("Values for SL Ridge prediction")
+ print(ypredictRidge)
+
+
+ MSEOwnRidgePredict[i] = MSE(y_test,ypredictOwnRidge)
+ MSERidgePredict[i] = MSE(y_test,ypredictRidge)
+
+ print("Beta values for own Ridge implementation")
+ print(OwnRidgeBeta) #Intercept is given by mean of target variable
+ print("Beta values for Scikit-Learn Ridge implementation")
+ print(RegRidge.coef_)
+ print('Intercept from own implementation:')
+ print(intercept_)
+ print('Intercept from Scikit-Learn Ridge implementation')
+ print(RegRidge.intercept_)
+
+
+
+# Now plot the results
+
+plt.figure()
+plt.plot(np.log10(lambdas), MSEOwnRidgePredict, 'b--', label = 'MSE own Ridge Test')
+plt.plot(np.log10(lambdas), MSERidgePredict, 'g--', label = 'MSE SL Ridge Test')
+
+plt.xlabel('log10(lambda)')
+plt.ylabel('MSE')
+plt.legend()
+plt.show()
+
@@ -505,7 +714,7 @@ In this case our matrix inversion was actually possible. The obvious question no
21
22
...
- 83
+ 87
»
diff --git a/doc/pub/week38/html/._week38-bs013.html b/doc/pub/week38/html/._week38-bs013.html
index c3f414954..e62d0af57 100644
--- a/doc/pub/week38/html/._week38-bs013.html
+++ b/doc/pub/week38/html/._week38-bs013.html
@@ -68,6 +68,16 @@ Automatically generated HTML file from DocOnce source
2,
None,
'to-think-about-first-part'),
+ ('More thinking', 2, None, 'more-thinking'),
+ ('Still thinking', 2, None, 'still-thinking'),
+ ('Linear Regression code, Intercept handling first',
+ 2,
+ None,
+ 'linear-regression-code-intercept-handling-first'),
+ ('What does centering mean mathematically?',
+ 2,
+ None,
+ 'what-does-centering-mean-mathematically'),
('More complicated Example: The Ising model',
2,
None,
@@ -306,80 +316,84 @@ MathJax.Hub.Config({
Cross-validation in brief
Code Example for Cross-validation and \( k \)-fold Cross-validation
To think about, first part
- More complicated Example: The Ising model
- Reformulating the problem to suit regression
- Linear regression
- Singular Value decomposition
- The one-dimensional Ising model
- Ridge regression
- LASSO regression
- Performance as function of the regularization parameter
- Finding the optimal value of \( \lambda \)
- Logistic Regression
- Classification problems
- Optimization and Deep learning
- Basics
- Linear classifier
- Some selected properties
- Simple example
- Plotting the mean value for each group
- The logistic function
- Examples of likelihood functions used in logistic regression and nueral networks
- Two parameters
- Maximum likelihood
- The cost function rewritten
- Minimizing the cross entropy
- A more compact expression
- Extending to more predictors
- Including more classes
- More classes
- Friday September 24
- Wisconsin Cancer Data
- Using the correlation matrix
- Discussing the correlation data
- Other measures in classification studies: Cancer Data again
- Optimization, the central part of any Machine Learning algortithm
- Revisiting our Logistic Regression case
- The equations to solve
- Solving using Newton-Raphson's method
- Brief reminder on Newton-Raphson's method
- The equations
- Simple geometric interpretation
- Extending to more than one variable
- Steepest descent
- More on Steepest descent
- The ideal
- The sensitiveness of the gradient descent
- Convex functions
- Convex function
- Conditions on convex functions
- More on convex functions
- Some simple problems
- Friday September 25
- Standard steepest descent
- Gradient method
- Steepest descent method
- Steepest descent method
- Final expressions
- Steepest descent example
- Conjugate gradient method
- Conjugate gradient method
- Conjugate gradient method
- Conjugate gradient method
- Conjugate gradient method and iterations
- Conjugate gradient method
- Conjugate gradient method
- Conjugate gradient method
- Revisiting some of our first Linear Regression Encounters
- Gradient descent example
- The derivative of the cost/loss function
- The Hessian matrix
- Simple program
- Gradient Descent Example
- And a corresponding example using scikit-learn
- Gradient descent and Ridge
- Program example for gradient descent with Ridge Regression
- Using gradient descent methods, limitations
+ More thinking
+ Still thinking
+ Linear Regression code, Intercept handling first
+ What does centering mean mathematically?
+ More complicated Example: The Ising model
+ Reformulating the problem to suit regression
+ Linear regression
+ Singular Value decomposition
+ The one-dimensional Ising model
+ Ridge regression
+ LASSO regression
+ Performance as function of the regularization parameter
+ Finding the optimal value of \( \lambda \)
+ Logistic Regression
+ Classification problems
+ Optimization and Deep learning
+ Basics
+ Linear classifier
+ Some selected properties
+ Simple example
+ Plotting the mean value for each group
+ The logistic function
+ Examples of likelihood functions used in logistic regression and nueral networks
+ Two parameters
+ Maximum likelihood
+ The cost function rewritten
+ Minimizing the cross entropy
+ A more compact expression
+ Extending to more predictors
+ Including more classes
+ More classes
+ Friday September 24
+ Wisconsin Cancer Data
+ Using the correlation matrix
+ Discussing the correlation data
+ Other measures in classification studies: Cancer Data again
+ Optimization, the central part of any Machine Learning algortithm
+ Revisiting our Logistic Regression case
+ The equations to solve
+ Solving using Newton-Raphson's method
+ Brief reminder on Newton-Raphson's method
+ The equations
+ Simple geometric interpretation
+ Extending to more than one variable
+ Steepest descent
+ More on Steepest descent
+ The ideal
+ The sensitiveness of the gradient descent
+ Convex functions
+ Convex function
+ Conditions on convex functions
+ More on convex functions
+ Some simple problems
+ Friday September 25
+ Standard steepest descent
+ Gradient method
+ Steepest descent method
+ Steepest descent method
+ Final expressions
+ Steepest descent example
+ Conjugate gradient method
+ Conjugate gradient method
+ Conjugate gradient method
+ Conjugate gradient method
+ Conjugate gradient method and iterations
+ Conjugate gradient method
+ Conjugate gradient method
+ Conjugate gradient method
+ Revisiting some of our first Linear Regression Encounters
+ Gradient descent example
+ The derivative of the cost/loss function
+ The Hessian matrix
+ Simple program
+ Gradient Descent Example
+ And a corresponding example using scikit-learn
+ Gradient descent and Ridge
+ Program example for gradient descent with Ridge Regression
+ Using gradient descent methods, limitations
@@ -395,27 +409,28 @@ MathJax.Hub.Config({
-The one-dimensional Ising model
+More complicated Example: The Ising model
-Let us bring back the Ising model again, but now with an additional
-focus on Ridge and Lasso regression as well. We repeat some of the
-basic parts of the Ising model and the setup of the training and test
-data. The one-dimensional Ising model with nearest neighbor
-interaction, no external field and a constant coupling constant \( J \) is
-given by
+The one-dimensional Ising model with nearest neighbor interaction, no
+external field and a constant coupling constant \( J \) is given by
$$
\begin{align}
H = -J \sum_{k}^L s_k s_{k + 1},
-\tag{7}
+\tag{1}
\end{align}
$$
-where \( s_i \in \{-1, 1\} \) and \( s_{N + 1} = s_1 \). The number of spins in the system is determined by \( L \). For the one-dimensional system there is no phase transition.
+
+where \( s_i \in \{-1, 1\} \) and \( s_{N + 1} = s_1 \). The number of spins
+in the system is determined by \( L \). For the one-dimensional system
+there is no phase transition.
-We will look at a system of \( L = 40 \) spins with a coupling constant of \( J = 1 \). To get enough training data we will generate 10000 states with their respective energies.
+We will look at a system of \( L = 40 \) spins with a coupling constant of
+\( J = 1 \). To get enough training data we will generate 10000 states
+with their respective energies.
@@ -426,7 +441,6 @@ We will look at a system of \( L = 40 \) spins with a coupling constant of \( J
import seaborn as sns
import scipy.linalg as scl
from sklearn.model_selection import train_test_split
-import sklearn.linear_model as skl
import tqdm
sns.set(color_codes=True)
cmap_args=dict(vmin=-1., vmax=1., cmap='seismic')
@@ -443,88 +457,11 @@ energies = np.<
energies[i] = - J * np.dot(spins[i], np.roll(spins[i], 1))
-A more general form for the one-dimensional Ising model is
-
-$$
-\begin{align}
- H = - \sum_j^L \sum_k^L s_j s_k J_{jk}.
-\tag{8}
-\end{align}
-$$
-
-
-Here we allow for interactions beyond the nearest neighbors and a more
-adaptive coupling matrix. This latter expression can be formulated as
-a matrix-product on the form
-$$
-\begin{align}
- H = X J,
-\tag{9}
-\end{align}
-$$
-
-
-where \( X_{jk} = s_j s_k \) and \( J \) is the matrix consisting of the
-elements \( -J_{jk} \). This form of writing the energy fits perfectly
-with the form utilized in linear regression, viz.
-$$
-\begin{align}
- \boldsymbol{y} = \boldsymbol{X}\boldsymbol{\beta} + \boldsymbol{\epsilon}.
-\tag{10}
-\end{align}
-$$
-
-We organize the data as we did above
-
-
-
-
X = np.zeros((n, L ** 2))
-for i in range(n):
- X[i] = np.outer(spins[i], spins[i]).ravel()
-y = energies
-X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.96)
-
-X_train_own = np.concatenate(
- (np.ones(len(X_train))[:, np.newaxis], X_train),
- axis=1
-)
-
-X_test_own = np.concatenate(
- (np.ones(len(X_test))[:, np.newaxis], X_test),
- axis=1
-)
-
-
-We will do all fitting with Scikit-Learn,
-
-
-
-
-
clf = skl.LinearRegression().fit(X_train, y_train)
-
-
-When extracting the \( J \)-matrix we make sure to remove the intercept
-
-
-
-
J_sk = clf.coef_.reshape(L, L)
-
-
-And then we plot the results
-
-
-
-
fig = plt.figure(figsize=(20, 14))
-im = plt.imshow(J_sk, **cmap_args)
-plt.title("LinearRegression from Scikit-learn", fontsize=18)
-plt.xticks(fontsize=18)
-plt.yticks(fontsize=18)
-cb = fig.colorbar(im)
-cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=18)
-plt.show()
-
-
-The results perfectly with our previous discussion where we used our own code.
+Here we use ordinary least squares
+regression to predict the energy for the nearest neighbor
+one-dimensional Ising model on a ring, i.e., the endpoints wrap
+around. We will use linear regression to fit a value for
+the coupling constant to achieve this.
@@ -552,7 +489,7 @@ The results perfectly with our previous discussion where we used our own code.
-Having explored the ordinary least squares we move on to ridge -regression. In ridge regression we include a regularizer. This -involves a new cost function which leads to a new estimate for the -weights \( \boldsymbol{\beta} \). This results in a penalized regression problem. The -cost function is given by +A more general form for the one-dimensional Ising model is $$ \begin{align} - C(\boldsymbol{X}, \boldsymbol{\beta}; \lambda) = (\boldsymbol{X}\boldsymbol{\beta} - \boldsymbol{y})^T(\boldsymbol{X}\boldsymbol{\beta} - \boldsymbol{y}) + \lambda \boldsymbol{\beta}^T\boldsymbol{\beta}. -\tag{11} + H = - \sum_j^L \sum_k^L s_j s_k J_{jk}. +\tag{2} \end{align} $$ +
+Here we allow for interactions beyond the nearest neighbors and a state dependent +coupling constant. This latter expression can be formulated as +a matrix-product +$$ +\begin{align} + \boldsymbol{H} = \boldsymbol{X} J, +\tag{3} +\end{align} +$$ + +
+where \( X_{jk} = s_j s_k \) and \( J \) is a matrix which consists of the +elements \( -J_{jk} \). This form of writing the energy fits perfectly +with the form utilized in linear regression, that is + +$$ +\begin{align} + \boldsymbol{y} = \boldsymbol{X}\boldsymbol{\beta} + \boldsymbol{\epsilon}, +\tag{4} +\end{align} +$$ + +
+We split the data in training and test data as discussed in the previous example +
-
_lambda = 0.1
-clf_ridge = skl.Ridge(alpha=_lambda).fit(X_train, y_train)
-J_ridge_sk = clf_ridge.coef_.reshape(L, L)
-fig = plt.figure(figsize=(20, 14))
-im = plt.imshow(J_ridge_sk, **cmap_args)
-plt.title("Ridge from Scikit-learn", fontsize=18)
-plt.xticks(fontsize=18)
-plt.yticks(fontsize=18)
-cb = fig.colorbar(im)
-cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=18)
-
-plt.show()
+X = np.zeros((n, L ** 2))
+for i in range(n):
+ X[i] = np.outer(spins[i], spins[i]).ravel()
+y = energies
+X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)
@@ -453,7 +482,7 @@ plt.show()
-In the Least Absolute Shrinkage and Selection Operator (LASSO)-method we get a third cost function. +In the ordinary least squares method we choose the cost function $$ \begin{align} - C(\boldsymbol{X}, \boldsymbol{\beta}; \lambda) = (\boldsymbol{X}\boldsymbol{\beta} - \boldsymbol{y})^T(\boldsymbol{X}\boldsymbol{\beta} - \boldsymbol{y}) + \lambda \sqrt{\boldsymbol{\beta}^T\boldsymbol{\beta}}. -\tag{12} + C(\boldsymbol{X}, \boldsymbol{\beta})= \frac{1}{n}\left\{(\boldsymbol{X}\boldsymbol{\beta} - \boldsymbol{y})^T(\boldsymbol{X}\boldsymbol{\beta} - \boldsymbol{y})\right\}. +\tag{5} \end{align} $$
-Finding the extremal point of this cost function is not so straight-forward as in least squares and ridge. We will therefore rely solely on the function ``Lasso`` from Scikit-Learn. +We then find the extremal point of \( C \) by taking the derivative with respect to \( \boldsymbol{\beta} \) as discussed above. +This yields the expression for \( \boldsymbol{\beta} \) to be + +$$ + \boldsymbol{\beta} = \frac{\boldsymbol{X}^T \boldsymbol{y}}{\boldsymbol{X}^T \boldsymbol{X}}, +$$ + +
+which immediately imposes some requirements on \( \boldsymbol{X} \) as there must exist +an inverse of \( \boldsymbol{X}^T \boldsymbol{X} \). If the expression we are modeling contains an +intercept, i.e., a constant term, we must make sure that the +first column of \( \boldsymbol{X} \) consists of \( 1 \). We do this here
-
clf_lasso = skl.Lasso(alpha=_lambda).fit(X_train, y_train)
-J_lasso_sk = clf_lasso.coef_.reshape(L, L)
-fig = plt.figure(figsize=(20, 14))
-im = plt.imshow(J_lasso_sk, **cmap_args)
-plt.title("Lasso from Scikit-learn", fontsize=18)
-plt.xticks(fontsize=18)
-plt.yticks(fontsize=18)
-cb = fig.colorbar(im)
-cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=18)
-
-plt.show()
+X_train_own = np.concatenate(
+ (np.ones(len(X_train))[:, np.newaxis], X_train),
+ axis=1
+)
+X_test_own = np.concatenate(
+ (np.ones(len(X_test))[:, np.newaxis], X_test),
+ axis=1
+)
-It is quite striking how LASSO breaks the symmetry of the coupling
-constant as opposed to ridge and OLS. We get a sparse solution with
-\( J_{j, j + 1} = -1 \).
+
+
def ols_inv(x: np.ndarray, y: np.ndarray) -> np.ndarray:
+ return scl.inv(x.T @ x) @ (x.T @ y)
+beta = ols_inv(X_train_own, y_train)
+
@@ -456,7 +480,7 @@ constant as opposed to ridge and OLS. We get a sparse solution with
-We see how the different models perform for a different set of values for \( \lambda \). +Doing the inversion directly turns out to be a bad idea since the matrix +\( \boldsymbol{X}^T\boldsymbol{X} \) is singular. An alternative approach is to use the singular +value decomposition. Using the definition of the Moore-Penrose +pseudoinverse we can write the equation for \( \boldsymbol{\beta} \) as + +$$ + \boldsymbol{\beta} = \boldsymbol{X}^{+}\boldsymbol{y}, +$$ + +
+where the pseudoinverse of \( \boldsymbol{X} \) is given by + +$$ + \boldsymbol{X}^{+} = \frac{\boldsymbol{X}^T}{\boldsymbol{X}^T\boldsymbol{X}}. +$$ + +
+Using singular value decomposition we can decompose the matrix \( \boldsymbol{X} = \boldsymbol{U}\boldsymbol{\Sigma} \boldsymbol{V}^T \), +where \( \boldsymbol{U} \) and \( \boldsymbol{V} \) are orthogonal(unitary) matrices and \( \boldsymbol{\Sigma} \) contains the singular values (more details below). +where \( X^{+} = V\Sigma^{+} U^T \). This reduces the equation for +\( \omega \) to +$$ +\begin{align} + \boldsymbol{\beta} = \boldsymbol{V}\boldsymbol{\Sigma}^{+} \boldsymbol{U}^T \boldsymbol{y}. +\tag{6} +\end{align} +$$ + +
+Note that solving this equation by actually doing the pseudoinverse +(which is what we will do) is not a good idea as this operation scales +as \( \mathcal{O}(n^3) \), where \( n \) is the number of elements in a +general matrix. Instead, doing \( QR \)-factorization and solving the +linear system as an equation would reduce this down to +\( \mathcal{O}(n^2) \) operations.
-
lambdas = np.logspace(-4, 5, 10)
+def ols_svd(x: np.ndarray, y: np.ndarray) -> np.ndarray:
+ u, s, v = scl.svd(x)
+ return v.T @ scl.pinv(scl.diagsvd(s, u.shape[0], v.shape[0])) @ u.T @ y
+
+
-train_errors = {
- "ols_sk": np.zeros(lambdas.size),
- "ridge_sk": np.zeros(lambdas.size),
- "lasso_sk": np.zeros(lambdas.size)
-}
+
+
beta = ols_svd(X_train_own,y_train)
+
+
+When extracting the \( J \)-matrix we need to make sure that we remove the intercept, as is done here
-test_errors = {
- "ols_sk": np.zeros(lambdas.size),
- "ridge_sk": np.zeros(lambdas.size),
- "lasso_sk": np.zeros(lambdas.size)
-}
+
-plot_counter = 1
+
+
J = beta[1:].reshape(L, L)
+
+
+A way of looking at the coefficients in \( J \) is to plot the matrices as images.
-fig = plt.figure(figsize=(32, 54))
-
-for i, _lambda in enumerate(tqdm.tqdm(lambdas)):
- for key, method in zip(
- ["ols_sk", "ridge_sk", "lasso_sk"],
- [skl.LinearRegression(), skl.Ridge(alpha=_lambda), skl.Lasso(alpha=_lambda)]
- ):
- method = method.fit(X_train, y_train)
-
- train_errors[key][i] = method.score(X_train, y_train)
- test_errors[key][i] = method.score(X_test, y_test)
-
- omega = method.coef_.reshape(L, L)
-
- plt.subplot(10, 5, plot_counter)
- plt.imshow(omega, **cmap_args)
- plt.title(r"%s, $\lambda = %.4f$" % (key, _lambda))
- plot_counter += 1
+
+
+
fig = plt.figure(figsize=(20, 14))
+im = plt.imshow(J, **cmap_args)
+plt.title("OLS", fontsize=18)
+plt.xticks(fontsize=18)
+plt.yticks(fontsize=18)
+cb = fig.colorbar(im)
+cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=18)
plt.show()
-We see that LASSO reaches a good solution for low
-values of \( \lambda \), but will "wither" when we increase \( \lambda \) too
-much. Ridge is more stable over a larger range of values for
-\( \lambda \), but eventually also fades away.
+It is interesting to note that OLS
+considers both \( J_{j, j + 1} = -0.5 \) and \( J_{j, j - 1} = -0.5 \) as
+valid matrix elements for \( J \).
+In our discussion below on hyperparameters and Ridge and Lasso regression we will see that
+this problem can be removed, partly and only with Lasso regression.
+
+
+In this case our matrix inversion was actually possible. The obvious question now is what is the mathematics behind the SVD?
@@ -472,7 +519,7 @@ much. Ridge is more stable over a larger range of values for
-To determine which value of \( \lambda \) is best we plot the accuracy of -the models when predicting the training and the testing set. We expect -the accuracy of the training set to be quite good, but if the accuracy -of the testing set is much lower this tells us that we might be -subject to an overfit model. The ideal scenario is an accuracy on the -testing set that is close to the accuracy of the training set. +Let us bring back the Ising model again, but now with an additional +focus on Ridge and Lasso regression as well. We repeat some of the +basic parts of the Ising model and the setup of the training and test +data. The one-dimensional Ising model with nearest neighbor +interaction, no external field and a constant coupling constant \( J \) is +given by + +$$ +\begin{align} + H = -J \sum_{k}^L s_k s_{k + 1}, +\tag{7} +\end{align} +$$ + +where \( s_i \in \{-1, 1\} \) and \( s_{N + 1} = s_1 \). The number of spins in the system is determined by \( L \). For the one-dimensional system there is no phase transition. + +
+We will look at a system of \( L = 40 \) spins with a coupling constant of \( J = 1 \). To get enough training data we will generate 10000 states with their respective energies.
+ +
import numpy as np
+import matplotlib.pyplot as plt
+from mpl_toolkits.axes_grid1 import make_axes_locatable
+import seaborn as sns
+import scipy.linalg as scl
+from sklearn.model_selection import train_test_split
+import sklearn.linear_model as skl
+import tqdm
+sns.set(color_codes=True)
+cmap_args=dict(vmin=-1., vmax=1., cmap='seismic')
+
+L = 40
+n = int(1e4)
+
+spins = np.random.choice([-1, 1], size=(n, L))
+J = 1.0
+
+energies = np.zeros(n)
+
+for i in range(n):
+ energies[i] = - J * np.dot(spins[i], np.roll(spins[i], 1))
++A more general form for the one-dimensional Ising model is + +$$ +\begin{align} + H = - \sum_j^L \sum_k^L s_j s_k J_{jk}. +\tag{8} +\end{align} +$$ + +
+Here we allow for interactions beyond the nearest neighbors and a more +adaptive coupling matrix. This latter expression can be formulated as +a matrix-product on the form +$$ +\begin{align} + H = X J, +\tag{9} +\end{align} +$$ + +
+where \( X_{jk} = s_j s_k \) and \( J \) is the matrix consisting of the +elements \( -J_{jk} \). This form of writing the energy fits perfectly +with the form utilized in linear regression, viz. +$$ +\begin{align} + \boldsymbol{y} = \boldsymbol{X}\boldsymbol{\beta} + \boldsymbol{\epsilon}. +\tag{10} +\end{align} +$$ + +We organize the data as we did above +
+ + +
X = np.zeros((n, L ** 2))
+for i in range(n):
+ X[i] = np.outer(spins[i], spins[i]).ravel()
+y = energies
+X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.96)
+
+X_train_own = np.concatenate(
+ (np.ones(len(X_train))[:, np.newaxis], X_train),
+ axis=1
+)
+
+X_test_own = np.concatenate(
+ (np.ones(len(X_test))[:, np.newaxis], X_test),
+ axis=1
+)
++We will do all fitting with Scikit-Learn, + +
+ + +
clf = skl.LinearRegression().fit(X_train, y_train)
++When extracting the \( J \)-matrix we make sure to remove the intercept +
+ + +
J_sk = clf.coef_.reshape(L, L)
++And then we plot the results +
+
fig = plt.figure(figsize=(20, 14))
-
-colors = {
- "ols_sk": "r",
- "ridge_sk": "y",
- "lasso_sk": "c"
-}
-
-for key in train_errors:
- plt.semilogx(
- lambdas,
- train_errors[key],
- colors[key],
- label="Train {0}".format(key),
- linewidth=4.0
- )
-
-for key in test_errors:
- plt.semilogx(
- lambdas,
- test_errors[key],
- colors[key] + "--",
- label="Test {0}".format(key),
- linewidth=4.0
- )
-plt.legend(loc="best", fontsize=18)
-plt.xlabel(r"$\lambda$", fontsize=18)
-plt.ylabel(r"$R^2$", fontsize=18)
-plt.tick_params(labelsize=18)
+im = plt.imshow(J_sk, **cmap_args)
+plt.title("LinearRegression from Scikit-learn", fontsize=18)
+plt.xticks(fontsize=18)
+plt.yticks(fontsize=18)
+cb = fig.colorbar(im)
+cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=18)
plt.show()
-From the above figure we can see that LASSO with \( \lambda = 10^{-2} \) -achieves a very good accuracy on the test set. This by far surpasses the -other models for all values of \( \lambda \). +The results perfectly with our previous discussion where we used our own code.
@@ -470,7 +566,7 @@ other models for all values of \( \lambda \).
- + -
-In linear regression our main interest was centered on learning the -coefficients of a functional fit (say a polynomial) in order to be -able to predict the response of a continuous variable on some unseen -data. The fit to the continuous variable \( y_i \) is based on some -independent variables \( \hat{x}_i \). Linear regression resulted in -analytical expressions for standard ordinary Least Squares or Ridge -regression (in terms of matrices to invert) for several quantities, -ranging from the variance and thereby the confidence intervals of the -parameters \( \hat{\beta} \) to the mean squared error. If we can invert -the product of the design matrices, linear regression gives then a -simple recipe for fitting our data. +Having explored the ordinary least squares we move on to ridge +regression. In ridge regression we include a regularizer. This +involves a new cost function which leads to a new estimate for the +weights \( \boldsymbol{\beta} \). This results in a penalized regression problem. The +cost function is given by +$$ +\begin{align} + C(\boldsymbol{X}, \boldsymbol{\beta}; \lambda) = (\boldsymbol{X}\boldsymbol{\beta} - \boldsymbol{y})^T(\boldsymbol{X}\boldsymbol{\beta} - \boldsymbol{y}) + \lambda \boldsymbol{\beta}^T\boldsymbol{\beta}. +\tag{11} +\end{align} +$$ + +
+ + +
_lambda = 0.1
+clf_ridge = skl.Ridge(alpha=_lambda).fit(X_train, y_train)
+J_ridge_sk = clf_ridge.coef_.reshape(L, L)
+fig = plt.figure(figsize=(20, 14))
+im = plt.imshow(J_ridge_sk, **cmap_args)
+plt.title("Ridge from Scikit-learn", fontsize=18)
+plt.xticks(fontsize=18)
+plt.yticks(fontsize=18)
+cb = fig.colorbar(im)
+cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=18)
+
+plt.show()
+
@@ -436,7 +467,7 @@ simple recipe for fitting our data.
- + -
-Classification problems, however, are concerned with outcomes taking -the form of discrete variables (i.e. categories). We may for example, -on the basis of DNA sequencing for a number of patients, like to find -out which mutations are important for a certain disease; or based on -scans of various patients' brains, figure out if there is a tumor or -not; or given a specific physical system, we'd like to identify its -state, say whether it is an ordered or disordered system (typical -situation in solid state physics); or classify the status of a -patient, whether she/he has a stroke or not and many other similar -situations. +In the Least Absolute Shrinkage and Selection Operator (LASSO)-method we get a third cost function. + +$$ +\begin{align} + C(\boldsymbol{X}, \boldsymbol{\beta}; \lambda) = (\boldsymbol{X}\boldsymbol{\beta} - \boldsymbol{y})^T(\boldsymbol{X}\boldsymbol{\beta} - \boldsymbol{y}) + \lambda \sqrt{\boldsymbol{\beta}^T\boldsymbol{\beta}}. +\tag{12} +\end{align} +$$
-The most common situation we encounter when we apply logistic -regression is that of two possible outcomes, normally denoted as a -binary outcome, true or false, positive or negative, success or -failure etc. +Finding the extremal point of this cost function is not so straight-forward as in least squares and ridge. We will therefore rely solely on the function ``Lasso`` from Scikit-Learn. + +
+ + +
clf_lasso = skl.Lasso(alpha=_lambda).fit(X_train, y_train)
+J_lasso_sk = clf_lasso.coef_.reshape(L, L)
+fig = plt.figure(figsize=(20, 14))
+im = plt.imshow(J_lasso_sk, **cmap_args)
+plt.title("Lasso from Scikit-learn", fontsize=18)
+plt.xticks(fontsize=18)
+plt.yticks(fontsize=18)
+cb = fig.colorbar(im)
+cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=18)
+
+plt.show()
++It is quite striking how LASSO breaks the symmetry of the coupling +constant as opposed to ridge and OLS. We get a sparse solution with +\( J_{j, j + 1} = -1 \).
@@ -441,7 +470,7 @@ failure etc.
-Logistic regression will also serve as our stepping stone towards -neural network algorithms and supervised deep learning. For logistic -learning, the minimization of the cost function leads to a non-linear -equation in the parameters \( \hat{\beta} \). The optimization of the -problem calls therefore for minimization algorithms. This forms the -bottle neck of all machine learning algorithms, namely how to find -reliable minima of a multi-variable function. This leads us to the -family of gradient descent methods. The latter are the working horses -of basically all modern machine learning algorithms. +We see how the different models perform for a different set of values for \( \lambda \).
-We note also that many of the topics discussed here on logistic -regression are also commonly used in modern supervised Deep Learning -models, as we will see later. + + +
lambdas = np.logspace(-4, 5, 10)
+
+train_errors = {
+ "ols_sk": np.zeros(lambdas.size),
+ "ridge_sk": np.zeros(lambdas.size),
+ "lasso_sk": np.zeros(lambdas.size)
+}
+
+test_errors = {
+ "ols_sk": np.zeros(lambdas.size),
+ "ridge_sk": np.zeros(lambdas.size),
+ "lasso_sk": np.zeros(lambdas.size)
+}
+
+plot_counter = 1
+
+fig = plt.figure(figsize=(32, 54))
+
+for i, _lambda in enumerate(tqdm.tqdm(lambdas)):
+ for key, method in zip(
+ ["ols_sk", "ridge_sk", "lasso_sk"],
+ [skl.LinearRegression(), skl.Ridge(alpha=_lambda), skl.Lasso(alpha=_lambda)]
+ ):
+ method = method.fit(X_train, y_train)
+
+ train_errors[key][i] = method.score(X_train, y_train)
+ test_errors[key][i] = method.score(X_test, y_test)
+
+ omega = method.coef_.reshape(L, L)
+
+ plt.subplot(10, 5, plot_counter)
+ plt.imshow(omega, **cmap_args)
+ plt.title(r"%s, $\lambda = %.4f$" % (key, _lambda))
+ plot_counter += 1
+
+plt.show()
++We see that LASSO reaches a good solution for low +values of \( \lambda \), but will "wither" when we increase \( \lambda \) too +much. Ridge is more stable over a larger range of values for +\( \lambda \), but eventually also fades away.
@@ -439,7 +486,7 @@ models, as we will see later.
- + -
-We consider the case where the dependent variables, also called the -responses or the outcomes, \( y_i \) are discrete and only take values -from \( k=0,\dots,K-1 \) (i.e. \( K \) classes). +To determine which value of \( \lambda \) is best we plot the accuracy of +the models when predicting the training and the testing set. We expect +the accuracy of the training set to be quite good, but if the accuracy +of the testing set is much lower this tells us that we might be +subject to an overfit model. The ideal scenario is an accuracy on the +testing set that is close to the accuracy of the training set.
-The goal is to predict the -output classes from the design matrix \( \hat{X}\in\mathbb{R}^{n\times p} \) -made of \( n \) samples, each of which carries \( p \) features or predictors. The -primary goal is to identify the classes to which new unseen samples -belong. + +
fig = plt.figure(figsize=(20, 14))
+
+colors = {
+ "ols_sk": "r",
+ "ridge_sk": "y",
+ "lasso_sk": "c"
+}
+
+for key in train_errors:
+ plt.semilogx(
+ lambdas,
+ train_errors[key],
+ colors[key],
+ label="Train {0}".format(key),
+ linewidth=4.0
+ )
+
+for key in test_errors:
+ plt.semilogx(
+ lambdas,
+ test_errors[key],
+ colors[key] + "--",
+ label="Test {0}".format(key),
+ linewidth=4.0
+ )
+plt.legend(loc="best", fontsize=18)
+plt.xlabel(r"$\lambda$", fontsize=18)
+plt.ylabel(r"$R^2$", fontsize=18)
+plt.tick_params(labelsize=18)
+plt.show()
+-Let us specialize to the case of two classes only, with outputs -\( y_i=0 \) and \( y_i=1 \). Our outcomes could represent the status of a -credit card user that could default or not on her/his credit card -debt. That is - -$$ -y_i = \begin{bmatrix} 0 & \mathrm{no}\\ 1 & \mathrm{yes} \end{bmatrix}. -$$ +From the above figure we can see that LASSO with \( \lambda = 10^{-2} \) +achieves a very good accuracy on the test set. This by far surpasses the +other models for all values of \( \lambda \).
@@ -445,7 +484,7 @@ $$
- + -
-Before moving to the logistic model, let us try to use our linear -regression model to classify these two outcomes. We could for example -fit a linear model to the default case if \( y_i > 0.5 \) and the no -default case \( y_i \leq 0.5 \). - -
-We would then have our -weighted linear combination, namely -$$ -\begin{equation} -\hat{y} = \hat{X}^T\hat{\beta} + \hat{\epsilon}, -\tag{13} -\end{equation} -$$ - -where \( \hat{y} \) is a vector representing the possible outcomes, \( \hat{X} \) is our -\( n\times p \) design matrix and \( \hat{\beta} \) represents our estimators/predictors. +In linear regression our main interest was centered on learning the +coefficients of a functional fit (say a polynomial) in order to be +able to predict the response of a continuous variable on some unseen +data. The fit to the continuous variable \( y_i \) is based on some +independent variables \( \hat{x}_i \). Linear regression resulted in +analytical expressions for standard ordinary Least Squares or Ridge +regression (in terms of matrices to invert) for several quantities, +ranging from the variance and thereby the confidence intervals of the +parameters \( \hat{\beta} \) to the mean squared error. If we can invert +the product of the design matrices, linear regression gives then a +simple recipe for fitting our data.
@@ -442,7 +450,7 @@ where \( \hat{y} \) is a vector representing the possible outcomes, \( \hat{X} \
- + -
-The main problem with our function is that it takes values on the -entire real axis. In the case of logistic regression, however, the -labels \( y_i \) are discrete variables. A typical example is the credit -card data discussed below here, where we can set the state of -defaulting the debt to \( y_i=1 \) and not to \( y_i=0 \) for one the persons -in the data set (see the full example below). +Classification problems, however, are concerned with outcomes taking +the form of discrete variables (i.e. categories). We may for example, +on the basis of DNA sequencing for a number of patients, like to find +out which mutations are important for a certain disease; or based on +scans of various patients' brains, figure out if there is a tumor or +not; or given a specific physical system, we'd like to identify its +state, say whether it is an ordered or disordered system (typical +situation in solid state physics); or classify the status of a +patient, whether she/he has a stroke or not and many other similar +situations.
-One simple way to get a discrete output is to have sign -functions that map the output of a linear regressor to values \( \{0,1\} \), -\( f(s_i)=sign(s_i)=1 \) if \( s_i\ge 0 \) and 0 if otherwise. -We will encounter this model in our first demonstration of neural networks. Historically it is called the ``perceptron" model in the machine learning -literature. This model is extremely simple. However, in many cases it is more -favorable to use a ``soft" classifier that outputs -the probability of a given category. This leads us to the logistic function. +The most common situation we encounter when we apply logistic +regression is that of two possible outcomes, normally denoted as a +binary outcome, true or false, positive or negative, success or +failure etc.
@@ -440,7 +455,7 @@ the probability of a given category. This leads us to the logistic function.
-The following example on data for coronary heart disease (CHD) as function of age may serve as an illustration. In the code here we read and plot whether a person has had CHD (output = 1) or not (output = 0). This ouput is plotted the person's against age. Clearly, the figure shows that attempting to make a standard linear regression fit may not be very meaningful. +Logistic regression will also serve as our stepping stone towards +neural network algorithms and supervised deep learning. For logistic +learning, the minimization of the cost function leads to a non-linear +equation in the parameters \( \hat{\beta} \). The optimization of the +problem calls therefore for minimization algorithms. This forms the +bottle neck of all machine learning algorithms, namely how to find +reliable minima of a multi-variable function. This leads us to the +family of gradient descent methods. The latter are the working horses +of basically all modern machine learning algorithms.
+We note also that many of the topics discussed here on logistic +regression are also commonly used in modern supervised Deep Learning +models, as we will see later. - -
# Common imports
-import os
-import numpy as np
-import pandas as pd
-import matplotlib.pyplot as plt
-from sklearn.linear_model import LinearRegression, Ridge, Lasso
-from sklearn.model_selection import train_test_split
-from sklearn.utils import resample
-from sklearn.metrics import mean_squared_error
-from IPython.display import display
-from pylab import plt, mpl
-plt.style.use('seaborn')
-mpl.rcParams['font.family'] = 'serif'
-
-# Where to save the figures and data files
-PROJECT_ROOT_DIR = "Results"
-FIGURE_ID = "Results/FigureFiles"
-DATA_ID = "DataFiles/"
-
-if not os.path.exists(PROJECT_ROOT_DIR):
- os.mkdir(PROJECT_ROOT_DIR)
-
-if not os.path.exists(FIGURE_ID):
- os.makedirs(FIGURE_ID)
-
-if not os.path.exists(DATA_ID):
- os.makedirs(DATA_ID)
-
-def image_path(fig_id):
- return os.path.join(FIGURE_ID, fig_id)
-
-def data_path(dat_id):
- return os.path.join(DATA_ID, dat_id)
-
-def save_fig(fig_id):
- plt.savefig(image_path(fig_id) + ".png", format='png')
-
-infile = open(data_path("chddata.csv"),'r')
-
-# Read the chd data as csv file and organize the data into arrays with age group, age, and chd
-chd = pd.read_csv(infile, names=('ID', 'Age', 'Agegroup', 'CHD'))
-chd.columns = ['ID', 'Age', 'Agegroup', 'CHD']
-output = chd['CHD']
-age = chd['Age']
-agegroup = chd['Agegroup']
-numberID = chd['ID']
-display(chd)
-
-plt.scatter(age, output, marker='o')
-plt.axis([18,70.0,-0.1, 1.2])
-plt.xlabel(r'Age')
-plt.ylabel(r'CHD')
-plt.title(r'Age distribution and Coronary heart disease')
-plt.show()
-
@@ -484,7 +453,7 @@ plt.show()
- + -
-What we could attempt however is to plot the mean value for each group. +We consider the case where the dependent variables, also called the +responses or the outcomes, \( y_i \) are discrete and only take values +from \( k=0,\dots,K-1 \) (i.e. \( K \) classes).
+The goal is to predict the +output classes from the design matrix \( \hat{X}\in\mathbb{R}^{n\times p} \) +made of \( n \) samples, each of which carries \( p \) features or predictors. The +primary goal is to identify the classes to which new unseen samples +belong. - -
agegroupmean = np.array([0.1, 0.133, 0.250, 0.333, 0.462, 0.625, 0.765, 0.800])
-group = np.array([1, 2, 3, 4, 5, 6, 7, 8])
-plt.plot(group, agegroupmean, "r-")
-plt.axis([0,9,0, 1.0])
-plt.xlabel(r'Age group')
-plt.ylabel(r'CHD mean values')
-plt.title(r'Mean values for each age group')
-plt.show()
--We are now trying to find a function \( f(y\vert x) \), that is a function which gives us an expected value for the output \( y \) with a given input \( x \). -In standard linear regression with a linear dependence on \( x \), we would write this in terms of our model +Let us specialize to the case of two classes only, with outputs +\( y_i=0 \) and \( y_i=1 \). Our outcomes could represent the status of a +credit card user that could default or not on her/his credit card +debt. That is + $$ -f(y_i\vert x_i)=\beta_0+\beta_1 x_i. +y_i = \begin{bmatrix} 0 & \mathrm{no}\\ 1 & \mathrm{yes} \end{bmatrix}. $$ -
-This expression implies however that \( f(y_i\vert x_i) \) could take any -value from minus infinity to plus infinity. If we however let -\( f(y\vert y) \) be represented by the mean value, the above example -shows us that we can constrain the function to take values between -zero and one, that is we have \( 0 \le f(y_i\vert x_i) \le 1 \). Looking -at our last curve we see also that it has an S-shaped form. This leads -us to a very popular model for the function \( f \), namely the so-called -Sigmoid function or logistic model. We will consider this function as -representing the probability for finding a value of \( y_i \) with a given -\( x_i \). -
@@ -457,7 +459,7 @@ representing the probability for finding a value of \( y_i \) with a given
-Another widely studied model, is the so-called -perceptron model, which is an example of a "hard classification" model. We -will encounter this model when we discuss neural networks as -well. Each datapoint is deterministically assigned to a category (i.e -\( y_i=0 \) or \( y_i=1 \)). In many cases, and the coronary heart disease data forms one of many such examples, it is favorable to have a "soft" -classifier that outputs the probability of a given category rather -than a single value. For example, given \( x_i \), the classifier -outputs the probability of being in a category \( k \). Logistic regression -is the most common example of a so-called soft classifier. In logistic -regression, the probability that a data point \( x_i \) -belongs to a category \( y_i=\{0,1\} \) is given by the so-called logit function (or Sigmoid) which is meant to represent the likelihood for a given event, +Before moving to the logistic model, let us try to use our linear +regression model to classify these two outcomes. We could for example +fit a linear model to the default case if \( y_i > 0.5 \) and the no +default case \( y_i \leq 0.5 \). + +
+We would then have our +weighted linear combination, namely $$ -p(t) = \frac{1}{1+\mathrm \exp{-t}}=\frac{\exp{t}}{1+\mathrm \exp{t}}. +\begin{equation} +\hat{y} = \hat{X}^T\hat{\beta} + \hat{\epsilon}, +\tag{13} +\end{equation} $$ -Note that \( 1-p(t)= p(-t) \). +where \( \hat{y} \) is a vector representing the possible outcomes, \( \hat{X} \) is our +\( n\times p \) design matrix and \( \hat{\beta} \) represents our estimators/predictors.
@@ -441,7 +456,7 @@ Note that \( 1-p(t)= p(-t) \).
-The following code plots the logistic function, the step function and other functions we will encounter from here and on. +The main problem with our function is that it takes values on the +entire real axis. In the case of logistic regression, however, the +labels \( y_i \) are discrete variables. A typical example is the credit +card data discussed below here, where we can set the state of +defaulting the debt to \( y_i=1 \) and not to \( y_i=0 \) for one the persons +in the data set (see the full example below).
+One simple way to get a discrete output is to have sign +functions that map the output of a linear regressor to values \( \{0,1\} \), +\( f(s_i)=sign(s_i)=1 \) if \( s_i\ge 0 \) and 0 if otherwise. +We will encounter this model in our first demonstration of neural networks. Historically it is called the ``perceptron" model in the machine learning +literature. This model is extremely simple. However, in many cases it is more +favorable to use a ``soft" classifier that outputs +the probability of a given category. This leads us to the logistic function. - -
"""The sigmoid function (or the logistic curve) is a
-function that takes any real number, z, and outputs a number (0,1).
-It is useful in neural networks for assigning weights on a relative scale.
-The value z is the weighted sum of parameters involved in the learning algorithm."""
-
-import numpy
-import matplotlib.pyplot as plt
-import math as mt
-
-z = numpy.arange(-5, 5, .1)
-sigma_fn = numpy.vectorize(lambda z: 1/(1+numpy.exp(-z)))
-sigma = sigma_fn(z)
-
-fig = plt.figure()
-ax = fig.add_subplot(111)
-ax.plot(z, sigma)
-ax.set_ylim([-0.1, 1.1])
-ax.set_xlim([-5,5])
-ax.grid(True)
-ax.set_xlabel('z')
-ax.set_title('sigmoid function')
-
-plt.show()
-
-"""Step Function"""
-z = numpy.arange(-5, 5, .02)
-step_fn = numpy.vectorize(lambda z: 1.0 if z >= 0.0 else 0.0)
-step = step_fn(z)
-
-fig = plt.figure()
-ax = fig.add_subplot(111)
-ax.plot(z, step)
-ax.set_ylim([-0.5, 1.5])
-ax.set_xlim([-5,5])
-ax.grid(True)
-ax.set_xlabel('z')
-ax.set_title('step function')
-
-plt.show()
-
-"""tanh Function"""
-z = numpy.arange(-2*mt.pi, 2*mt.pi, 0.1)
-t = numpy.tanh(z)
-
-fig = plt.figure()
-ax = fig.add_subplot(111)
-ax.plot(z, t)
-ax.set_ylim([-1.0, 1.0])
-ax.set_xlim([-2*mt.pi,2*mt.pi])
-ax.grid(True)
-ax.set_xlabel('z')
-ax.set_title('tanh function')
-
-plt.show()
-
@@ -484,7 +454,7 @@ plt.show()
-We assume now that we have two classes with \( y_i \) either \( 0 \) or \( 1 \). Furthermore we assume also that we have only two parameters \( \beta \) in our fitting of the Sigmoid function, that is we define probabilities -$$ -\begin{align*} -p(y_i=1|x_i,\hat{\beta}) &= \frac{\exp{(\beta_0+\beta_1x_i)}}{1+\exp{(\beta_0+\beta_1x_i)}},\nonumber\\ -p(y_i=0|x_i,\hat{\beta}) &= 1 - p(y_i=1|x_i,\hat{\beta}), -\end{align*} -$$ - -where \( \hat{\beta} \) are the weights we wish to extract from data, in our case \( \beta_0 \) and \( \beta_1 \). +The following example on data for coronary heart disease (CHD) as function of age may serve as an illustration. In the code here we read and plot whether a person has had CHD (output = 1) or not (output = 0). This ouput is plotted the person's against age. Clearly, the figure shows that attempting to make a standard linear regression fit may not be very meaningful.
-Note that we used -$$ -p(y_i=0\vert x_i, \hat{\beta}) = 1-p(y_i=1\vert x_i, \hat{\beta}). -$$ + +
# Common imports
+import os
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+from sklearn.linear_model import LinearRegression, Ridge, Lasso
+from sklearn.model_selection import train_test_split
+from sklearn.utils import resample
+from sklearn.metrics import mean_squared_error
+from IPython.display import display
+from pylab import plt, mpl
+plt.style.use('seaborn')
+mpl.rcParams['font.family'] = 'serif'
+
+# Where to save the figures and data files
+PROJECT_ROOT_DIR = "Results"
+FIGURE_ID = "Results/FigureFiles"
+DATA_ID = "DataFiles/"
+
+if not os.path.exists(PROJECT_ROOT_DIR):
+ os.mkdir(PROJECT_ROOT_DIR)
+
+if not os.path.exists(FIGURE_ID):
+ os.makedirs(FIGURE_ID)
+
+if not os.path.exists(DATA_ID):
+ os.makedirs(DATA_ID)
+
+def image_path(fig_id):
+ return os.path.join(FIGURE_ID, fig_id)
+
+def data_path(dat_id):
+ return os.path.join(DATA_ID, dat_id)
+
+def save_fig(fig_id):
+ plt.savefig(image_path(fig_id) + ".png", format='png')
+
+infile = open(data_path("chddata.csv"),'r')
+
+# Read the chd data as csv file and organize the data into arrays with age group, age, and chd
+chd = pd.read_csv(infile, names=('ID', 'Age', 'Agegroup', 'CHD'))
+chd.columns = ['ID', 'Age', 'Agegroup', 'CHD']
+output = chd['CHD']
+age = chd['Age']
+agegroup = chd['Agegroup']
+numberID = chd['ID']
+display(chd)
+
+plt.scatter(age, output, marker='o')
+plt.axis([18,70.0,-0.1, 1.2])
+plt.xlabel(r'Age')
+plt.ylabel(r'CHD')
+plt.title(r'Age distribution and Coronary heart disease')
+plt.show()
+
@@ -440,7 +498,7 @@ $$
- + -
-In order to define the total likelihood for all possible outcomes from a -dataset \( \mathcal{D}=\{(y_i,x_i)\} \), with the binary labels -\( y_i\in\{0,1\} \) and where the data points are drawn independently, we use the so-called Maximum Likelihood Estimation (MLE) principle. -We aim thus at maximizing -the probability of seeing the observed data. We can then approximate the -likelihood in terms of the product of the individual probabilities of a specific outcome \( y_i \), that is +What we could attempt however is to plot the mean value for each group. + +
+ + +
agegroupmean = np.array([0.1, 0.133, 0.250, 0.333, 0.462, 0.625, 0.765, 0.800])
+group = np.array([1, 2, 3, 4, 5, 6, 7, 8])
+plt.plot(group, agegroupmean, "r-")
+plt.axis([0,9,0, 1.0])
+plt.xlabel(r'Age group')
+plt.ylabel(r'CHD mean values')
+plt.title(r'Mean values for each age group')
+plt.show()
++We are now trying to find a function \( f(y\vert x) \), that is a function which gives us an expected value for the output \( y \) with a given input \( x \). +In standard linear regression with a linear dependence on \( x \), we would write this in terms of our model $$ -\begin{align*} -P(\mathcal{D}|\hat{\beta})& = \prod_{i=1}^n \left[p(y_i=1|x_i,\hat{\beta})\right]^{y_i}\left[1-p(y_i=1|x_i,\hat{\beta}))\right]^{1-y_i}\nonumber \\ -\end{align*} +f(y_i\vert x_i)=\beta_0+\beta_1 x_i. $$ -from which we obtain the log-likelihood and our cost/loss function -$$ -\mathcal{C}(\hat{\beta}) = \sum_{i=1}^n \left( y_i\log{p(y_i=1|x_i,\hat{\beta})} + (1-y_i)\log\left[1-p(y_i=1|x_i,\hat{\beta}))\right]\right). -$$ +
+This expression implies however that \( f(y_i\vert x_i) \) could take any +value from minus infinity to plus infinity. If we however let +\( f(y\vert y) \) be represented by the mean value, the above example +shows us that we can constrain the function to take values between +zero and one, that is we have \( 0 \le f(y_i\vert x_i) \le 1 \). Looking +at our last curve we see also that it has an S-shaped form. This leads +us to a very popular model for the function \( f \), namely the so-called +Sigmoid function or logistic model. We will consider this function as +representing the probability for finding a value of \( y_i \) with a given +\( x_i \).
@@ -441,7 +471,7 @@ $$
-Reordering the logarithms, we can rewrite the cost/loss function as +Another widely studied model, is the so-called +perceptron model, which is an example of a "hard classification" model. We +will encounter this model when we discuss neural networks as +well. Each datapoint is deterministically assigned to a category (i.e +\( y_i=0 \) or \( y_i=1 \)). In many cases, and the coronary heart disease data forms one of many such examples, it is favorable to have a "soft" +classifier that outputs the probability of a given category rather +than a single value. For example, given \( x_i \), the classifier +outputs the probability of being in a category \( k \). Logistic regression +is the most common example of a so-called soft classifier. In logistic +regression, the probability that a data point \( x_i \) +belongs to a category \( y_i=\{0,1\} \) is given by the so-called logit function (or Sigmoid) which is meant to represent the likelihood for a given event, $$ -\mathcal{C}(\hat{\beta}) = \sum_{i=1}^n \left(y_i(\beta_0+\beta_1x_i) -\log{(1+\exp{(\beta_0+\beta_1x_i)})}\right). +p(t) = \frac{1}{1+\mathrm \exp{-t}}=\frac{\exp{t}}{1+\mathrm \exp{t}}. $$ -
-The maximum likelihood estimator is defined as the set of parameters that maximize the log-likelihood where we maximize with respect to \( \beta \). -Since the cost (error) function is just the negative log-likelihood, for logistic regression we have that -$$ -\mathcal{C}(\hat{\beta})=-\sum_{i=1}^n \left(y_i(\beta_0+\beta_1x_i) -\log{(1+\exp{(\beta_0+\beta_1x_i)})}\right). -$$ - -This equation is known in statistics as the cross entropy. Finally, we note that just as in linear regression, -in practice we often supplement the cross-entropy with additional regularization terms, usually \( L_1 \) and \( L_2 \) regularization as we did for Ridge and Lasso regression. +Note that \( 1-p(t)= p(-t) \).
@@ -439,7 +455,7 @@ in practice we often supplement the cross-entropy with additional regularization
-The cross entropy is a convex function of the weights \( \hat{\beta} \) and, -therefore, any local minimizer is a global minimizer. +The following code plots the logistic function, the step function and other functions we will encounter from here and on.
-Minimizing this -cost function with respect to the two parameters \( \beta_0 \) and \( \beta_1 \) we obtain -$$ -\frac{\partial \mathcal{C}(\hat{\beta})}{\partial \beta_0} = -\sum_{i=1}^n \left(y_i -\frac{\exp{(\beta_0+\beta_1x_i)}}{1+\exp{(\beta_0+\beta_1x_i)}}\right), -$$ + +
"""The sigmoid function (or the logistic curve) is a
+function that takes any real number, z, and outputs a number (0,1).
+It is useful in neural networks for assigning weights on a relative scale.
+The value z is the weighted sum of parameters involved in the learning algorithm."""
-and
-$$
-\frac{\partial \mathcal{C}(\hat{\beta})}{\partial \beta_1} = -\sum_{i=1}^n \left(y_ix_i -x_i\frac{\exp{(\beta_0+\beta_1x_i)}}{1+\exp{(\beta_0+\beta_1x_i)}}\right).
-$$
+import numpy
+import matplotlib.pyplot as plt
+import math as mt
+z = numpy.arange(-5, 5, .1)
+sigma_fn = numpy.vectorize(lambda z: 1/(1+numpy.exp(-z)))
+sigma = sigma_fn(z)
+
+fig = plt.figure()
+ax = fig.add_subplot(111)
+ax.plot(z, sigma)
+ax.set_ylim([-0.1, 1.1])
+ax.set_xlim([-5,5])
+ax.grid(True)
+ax.set_xlabel('z')
+ax.set_title('sigmoid function')
+
+plt.show()
+
+"""Step Function"""
+z = numpy.arange(-5, 5, .02)
+step_fn = numpy.vectorize(lambda z: 1.0 if z >= 0.0 else 0.0)
+step = step_fn(z)
+
+fig = plt.figure()
+ax = fig.add_subplot(111)
+ax.plot(z, step)
+ax.set_ylim([-0.5, 1.5])
+ax.set_xlim([-5,5])
+ax.grid(True)
+ax.set_xlabel('z')
+ax.set_title('step function')
+
+plt.show()
+
+"""tanh Function"""
+z = numpy.arange(-2*mt.pi, 2*mt.pi, 0.1)
+t = numpy.tanh(z)
+
+fig = plt.figure()
+ax = fig.add_subplot(111)
+ax.plot(z, t)
+ax.set_ylim([-1.0, 1.0])
+ax.set_xlim([-2*mt.pi,2*mt.pi])
+ax.grid(True)
+ax.set_xlabel('z')
+ax.set_title('tanh function')
+
+plt.show()
+
@@ -440,7 +498,7 @@ $$
-Let us now define a vector \( \hat{y} \) with \( n \) elements \( y_i \), an -\( n\times p \) matrix \( \hat{X} \) which contains the \( x_i \) values and a -vector \( \hat{p} \) of fitted probabilities \( p(y_i\vert x_i,\hat{\beta}) \). We can rewrite in a more compact form the first -derivative of cost function as +We assume now that we have two classes with \( y_i \) either \( 0 \) or \( 1 \). Furthermore we assume also that we have only two parameters \( \beta \) in our fitting of the Sigmoid function, that is we define probabilities +$$ +\begin{align*} +p(y_i=1|x_i,\hat{\beta}) &= \frac{\exp{(\beta_0+\beta_1x_i)}}{1+\exp{(\beta_0+\beta_1x_i)}},\nonumber\\ +p(y_i=0|x_i,\hat{\beta}) &= 1 - p(y_i=1|x_i,\hat{\beta}), +\end{align*} +$$ -$$ -\frac{\partial \mathcal{C}(\hat{\beta})}{\partial \hat{\beta}} = -\hat{X}^T\left(\hat{y}-\hat{p}\right). -$$ +where \( \hat{\beta} \) are the weights we wish to extract from data, in our case \( \beta_0 \) and \( \beta_1 \).
-If we in addition define a diagonal matrix \( \hat{W} \) with elements -\( p(y_i\vert x_i,\hat{\beta})(1-p(y_i\vert x_i,\hat{\beta}) \), we can obtain a compact expression of the second derivative as - +Note that we used $$ -\frac{\partial^2 \mathcal{C}(\hat{\beta})}{\partial \hat{\beta}\partial \hat{\beta}^T} = \hat{X}^T\hat{W}\hat{X}. +p(y_i=0\vert x_i, \hat{\beta}) = 1-p(y_i=1\vert x_i, \hat{\beta}). $$
@@ -441,7 +454,7 @@ $$
- + -
-Within a binary classification problem, we can easily expand our model to include multiple predictors. Our ratio between likelihoods is then with \( p \) predictors +In order to define the total likelihood for all possible outcomes from a +dataset \( \mathcal{D}=\{(y_i,x_i)\} \), with the binary labels +\( y_i\in\{0,1\} \) and where the data points are drawn independently, we use the so-called Maximum Likelihood Estimation (MLE) principle. +We aim thus at maximizing +the probability of seeing the observed data. We can then approximate the +likelihood in terms of the product of the individual probabilities of a specific outcome \( y_i \), that is $$ -\log{ \frac{p(\hat{\beta}\hat{x})}{1-p(\hat{\beta}\hat{x})}} = \beta_0+\beta_1x_1+\beta_2x_2+\dots+\beta_px_p. +\begin{align*} +P(\mathcal{D}|\hat{\beta})& = \prod_{i=1}^n \left[p(y_i=1|x_i,\hat{\beta})\right]^{y_i}\left[1-p(y_i=1|x_i,\hat{\beta}))\right]^{1-y_i}\nonumber \\ +\end{align*} $$ -Here we defined \( \hat{x}=[1,x_1,x_2,\dots,x_p] \) and \( \hat{\beta}=[\beta_0, \beta_1, \dots, \beta_p] \) leading to +from which we obtain the log-likelihood and our cost/loss function $$ -p(\hat{\beta}\hat{x})=\frac{ \exp{(\beta_0+\beta_1x_1+\beta_2x_2+\dots+\beta_px_p)}}{1+\exp{(\beta_0+\beta_1x_1+\beta_2x_2+\dots+\beta_px_p)}}. +\mathcal{C}(\hat{\beta}) = \sum_{i=1}^n \left( y_i\log{p(y_i=1|x_i,\hat{\beta})} + (1-y_i)\log\left[1-p(y_i=1|x_i,\hat{\beta}))\right]\right). $$
@@ -434,7 +455,7 @@ $$
-Till now we have mainly focused on two classes, the so-called binary -system. Suppose we wish to extend to \( K \) classes. Let us for the sake -of simplicity assume we have only two predictors. We have then following model - +Reordering the logarithms, we can rewrite the cost/loss function as $$ -\log{\frac{p(C=1\vert x)}{p(K\vert x)}} = \beta_{10}+\beta_{11}x_1, -$$ - -and -$$ -\log{\frac{p(C=2\vert x)}{p(K\vert x)}} = \beta_{20}+\beta_{21}x_1, -$$ - -and so on till the class \( C=K-1 \) class -$$ -\log{\frac{p(C=K-1\vert x)}{p(K\vert x)}} = \beta_{(K-1)0}+\beta_{(K-1)1}x_1, +\mathcal{C}(\hat{\beta}) = \sum_{i=1}^n \left(y_i(\beta_0+\beta_1x_i) -\log{(1+\exp{(\beta_0+\beta_1x_i)})}\right). $$
-and the model is specified in term of \( K-1 \) so-called log-odds or -logit transformations. +The maximum likelihood estimator is defined as the set of parameters that maximize the log-likelihood where we maximize with respect to \( \beta \). +Since the cost (error) function is just the negative log-likelihood, for logistic regression we have that +$$ +\mathcal{C}(\hat{\beta})=-\sum_{i=1}^n \left(y_i(\beta_0+\beta_1x_i) -\log{(1+\exp{(\beta_0+\beta_1x_i)})}\right). +$$ + +This equation is known in statistics as the cross entropy. Finally, we note that just as in linear regression, +in practice we often supplement the cross-entropy with additional regularization terms, usually \( L_1 \) and \( L_2 \) regularization as we did for Ridge and Lasso regression.
@@ -446,7 +453,7 @@ and the model is specified in term of \( K-1 \) so-called log-odds or
-In our discussion of neural networks we will encounter the above again -in terms of a slightly modified function, the so-called Softmax function. +The cross entropy is a convex function of the weights \( \hat{\beta} \) and, +therefore, any local minimizer is a global minimizer.
-The softmax function is used in various multiclass classification -methods, such as multinomial logistic regression (also known as -softmax regression), multiclass linear discriminant analysis, naive -Bayes classifiers, and artificial neural networks. Specifically, in -multinomial logistic regression and linear discriminant analysis, the -input to the function is the result of \( K \) distinct linear functions, -and the predicted probability for the \( k \)-th class given a sample -vector \( \hat{x} \) and a weighting vector \( \hat{\beta} \) is (with two -predictors): +Minimizing this +cost function with respect to the two parameters \( \beta_0 \) and \( \beta_1 \) we obtain $$ -p(C=k\vert \mathbf {x} )=\frac{\exp{(\beta_{k0}+\beta_{k1}x_1)}}{1+\sum_{l=1}^{K-1}\exp{(\beta_{l0}+\beta_{l1}x_1)}}. +\frac{\partial \mathcal{C}(\hat{\beta})}{\partial \beta_0} = -\sum_{i=1}^n \left(y_i -\frac{\exp{(\beta_0+\beta_1x_i)}}{1+\exp{(\beta_0+\beta_1x_i)}}\right), $$ -It is easy to extend to more predictors. The final class is +and $$ -p(C=K\vert \mathbf {x} )=\frac{1}{1+\sum_{l=1}^{K-1}\exp{(\beta_{l0}+\beta_{l1}x_1)}}, +\frac{\partial \mathcal{C}(\hat{\beta})}{\partial \beta_1} = -\sum_{i=1}^n \left(y_ix_i -x_i\frac{\exp{(\beta_0+\beta_1x_i)}}{1+\exp{(\beta_0+\beta_1x_i)}}\right). $$ -
-and they sum to one. Our earlier discussions were all specialized to -the case with two classes only. It is easy to see from the above that -what we derived earlier is compatible with these equations. - -
-To find the optimal parameters we would typically use a gradient -descent method. Newton's method and gradient descent methods are -discussed in the material on optimization -methods. -
@@ -458,7 +454,7 @@ methods.
+Let us now define a vector \( \hat{y} \) with \( n \) elements \( y_i \), an +\( n\times p \) matrix \( \hat{X} \) which contains the \( x_i \) values and a +vector \( \hat{p} \) of fitted probabilities \( p(y_i\vert x_i,\hat{\beta}) \). We can rewrite in a more compact form the first +derivative of cost function as + +$$ +\frac{\partial \mathcal{C}(\hat{\beta})}{\partial \hat{\beta}} = -\hat{X}^T\left(\hat{y}-\hat{p}\right). +$$ + +
+If we in addition define a diagonal matrix \( \hat{W} \) with elements +\( p(y_i\vert x_i,\hat{\beta})(1-p(y_i\vert x_i,\hat{\beta}) \), we can obtain a compact expression of the second derivative as + +$$ +\frac{\partial^2 \mathcal{C}(\hat{\beta})}{\partial \hat{\beta}\partial \hat{\beta}^T} = \hat{X}^T\hat{W}\hat{X}. +$$
@@ -423,7 +455,7 @@ MathJax.Hub.Config({
-We show here how we can use a simple regression case on the breast -cancer data using Logistic regression as our algorithm for -classification. +Within a binary classification problem, we can easily expand our model to include multiple predictors. Our ratio between likelihoods is then with \( p \) predictors +$$ +\log{ \frac{p(\hat{\beta}\hat{x})}{1-p(\hat{\beta}\hat{x})}} = \beta_0+\beta_1x_1+\beta_2x_2+\dots+\beta_px_p. +$$ -
+Here we defined \( \hat{x}=[1,x_1,x_2,\dots,x_p] \) and \( \hat{\beta}=[\beta_0, \beta_1, \dots, \beta_p] \) leading to +$$ +p(\hat{\beta}\hat{x})=\frac{ \exp{(\beta_0+\beta_1x_1+\beta_2x_2+\dots+\beta_px_p)}}{1+\exp{(\beta_0+\beta_1x_1+\beta_2x_2+\dots+\beta_px_p)}}. +$$ - -
import matplotlib.pyplot as plt
-import numpy as np
-from sklearn.model_selection import train_test_split
-from sklearn.datasets import load_breast_cancer
-from sklearn.linear_model import LogisticRegression
-
-# Load the data
-cancer = load_breast_cancer()
-
-X_train, X_test, y_train, y_test = train_test_split(cancer.data,cancer.target,random_state=0)
-print(X_train.shape)
-print(X_test.shape)
-# Logistic Regression
-logreg = LogisticRegression(solver='lbfgs')
-logreg.fit(X_train, y_train)
-print("Test set accuracy with Logistic Regression: {:.2f}".format(logreg.score(X_test,y_test)))
-#now scale the data
-from sklearn.preprocessing import StandardScaler
-scaler = StandardScaler()
-scaler.fit(X_train)
-X_train_scaled = scaler.transform(X_train)
-X_test_scaled = scaler.transform(X_test)
-# Logistic Regression
-logreg.fit(X_train_scaled, y_train)
-print("Test set accuracy Logistic Regression with scaled data: {:.2f}".format(logreg.score(X_test_scaled,y_test)))
-
@@ -457,7 +448,7 @@ logreg.fit(X_train_scaled, y_train)
-In addition to the above scores, we could also study the covariance (and the correlation matrix). -We use Pandas to compute the correlation matrix. +Till now we have mainly focused on two classes, the so-called binary +system. Suppose we wish to extend to \( K \) classes. Let us for the sake +of simplicity assume we have only two predictors. We have then following model + +$$ +\log{\frac{p(C=1\vert x)}{p(K\vert x)}} = \beta_{10}+\beta_{11}x_1, +$$ + +and +$$ +\log{\frac{p(C=2\vert x)}{p(K\vert x)}} = \beta_{20}+\beta_{21}x_1, +$$ + +and so on till the class \( C=K-1 \) class +$$ +\log{\frac{p(C=K-1\vert x)}{p(K\vert x)}} = \beta_{(K-1)0}+\beta_{(K-1)1}x_1, +$$ +
+and the model is specified in term of \( K-1 \) so-called log-odds or +logit transformations. - -
import matplotlib.pyplot as plt
-import numpy as np
-from sklearn.model_selection import train_test_split
-from sklearn.datasets import load_breast_cancer
-from sklearn.linear_model import LogisticRegression
-cancer = load_breast_cancer()
-import pandas as pd
-# Making a data frame
-cancerpd = pd.DataFrame(cancer.data, columns=cancer.feature_names)
-
-fig, axes = plt.subplots(15,2,figsize=(10,20))
-malignant = cancer.data[cancer.target == 0]
-benign = cancer.data[cancer.target == 1]
-ax = axes.ravel()
-
-for i in range(30):
- _, bins = np.histogram(cancer.data[:,i], bins =50)
- ax[i].hist(malignant[:,i], bins = bins, alpha = 0.5)
- ax[i].hist(benign[:,i], bins = bins, alpha = 0.5)
- ax[i].set_title(cancer.feature_names[i])
- ax[i].set_yticks(())
-ax[0].set_xlabel("Feature magnitude")
-ax[0].set_ylabel("Frequency")
-ax[0].legend(["Malignant", "Benign"], loc ="best")
-fig.tight_layout()
-plt.show()
-
-import seaborn as sns
-correlation_matrix = cancerpd.corr().round(1)
-# use the heatmap function from seaborn to plot the correlation matrix
-# annot = True to print the values inside the square
-plt.figure(figsize=(15,8))
-sns.heatmap(data=correlation_matrix, annot=True)
-plt.show()
-
@@ -464,7 +460,7 @@ plt.show()
-In the above example we note two things. In the first plot we display -the overlap of benign and malignant tumors as functions of the various -features in the Wisconsing breast cancer data set. We see that for -some of the features we can distinguish clearly the benign and -malignant cases while for other features we cannot. This can point to -us which features may be of greater interest when we wish to classify -a benign or not benign tumour. +In our discussion of neural networks we will encounter the above again +in terms of a slightly modified function, the so-called Softmax function.
-In the second figure we have computed the so-called correlation -matrix, which in our case with thirty features becomes a \( 30\times 30 \) -matrix. +The softmax function is used in various multiclass classification +methods, such as multinomial logistic regression (also known as +softmax regression), multiclass linear discriminant analysis, naive +Bayes classifiers, and artificial neural networks. Specifically, in +multinomial logistic regression and linear discriminant analysis, the +input to the function is the result of \( K \) distinct linear functions, +and the predicted probability for the \( k \)-th class given a sample +vector \( \hat{x} \) and a weighting vector \( \hat{\beta} \) is (with two +predictors): + +$$ +p(C=k\vert \mathbf {x} )=\frac{\exp{(\beta_{k0}+\beta_{k1}x_1)}}{1+\sum_{l=1}^{K-1}\exp{(\beta_{l0}+\beta_{l1}x_1)}}. +$$ + +It is easy to extend to more predictors. The final class is +$$ +p(C=K\vert \mathbf {x} )=\frac{1}{1+\sum_{l=1}^{K-1}\exp{(\beta_{l0}+\beta_{l1}x_1)}}, +$$
-We constructed this matrix using pandas via the statements -
+and they sum to one. Our earlier discussions were all specialized to +the case with two classes only. It is easy to see from the above that +what we derived earlier is compatible with these equations. - -
cancerpd = pd.DataFrame(cancer.data, columns=cancer.feature_names)
--and then -
- - -
correlation_matrix = cancerpd.corr().round(1)
--Diagonalizing this matrix we can in turn say something about which -features are of relevance and which are not. This leads us to -the classical Principal Component Analysis (PCA) theorem with -applications. This will be discussed later this semester (week 43). +To find the optimal parameters we would typically use a gradient +descent method. Newton's method and gradient descent methods are +discussed in the material on optimization +methods.
@@ -457,7 +472,7 @@ applications. This will be discussed later this semester (48
+
import matplotlib.pyplot as plt
-import numpy as np
-from sklearn.model_selection import train_test_split
-from sklearn.datasets import load_breast_cancer
-from sklearn.linear_model import LogisticRegression
-
-# Load the data
-cancer = load_breast_cancer()
-
-X_train, X_test, y_train, y_test = train_test_split(cancer.data,cancer.target,random_state=0)
-print(X_train.shape)
-print(X_test.shape)
-# Logistic Regression
-logreg = LogisticRegression(solver='lbfgs')
-logreg.fit(X_train, y_train)
-print("Test set accuracy with Logistic Regression: {:.2f}".format(logreg.score(X_test,y_test)))
-#now scale the data
-from sklearn.preprocessing import StandardScaler
-scaler = StandardScaler()
-scaler.fit(X_train)
-X_train_scaled = scaler.transform(X_train)
-X_test_scaled = scaler.transform(X_test)
-# Logistic Regression
-logreg.fit(X_train_scaled, y_train)
-print("Test set accuracy Logistic Regression with scaled data: {:.2f}".format(logreg.score(X_test_scaled,y_test)))
-
-
-from sklearn.preprocessing import LabelEncoder
-from sklearn.model_selection import cross_validate
-#Cross validation
-accuracy = cross_validate(logreg,X_test_scaled,y_test,cv=10)['test_score']
-print(accuracy)
-print("Test set accuracy with Logistic Regression and scaled data: {:.2f}".format(logreg.score(X_test_scaled,y_test)))
-
-
-import scikitplot as skplt
-y_pred = logreg.predict(X_test_scaled)
-skplt.metrics.plot_confusion_matrix(y_test, y_pred, normalize=True)
-plt.show()
-y_probas = logreg.predict_proba(X_test_scaled)
-skplt.metrics.plot_roc(y_test, y_probas)
-plt.show()
-skplt.metrics.plot_cumulative_gain(y_test, y_probas)
-plt.show()
-
@@ -470,7 +437,7 @@ plt.show()
-Overview Video, why do we care about gradient methods? +We show here how we can use a simple regression case on the breast +cancer data using Logistic regression as our algorithm for +classification.
-Almost every problem in machine learning and data science starts with -a dataset \( X \), a model \( g(\beta) \), which is a function of the -parameters \( \beta \) and a cost function \( C(X, g(\beta)) \) that allows -us to judge how well the model \( g(\beta) \) explains the observations -\( X \). The model is fit by finding the values of \( \beta \) that minimize -the cost function. Ideally we would be able to solve for \( \beta \) -analytically, however this is not possible in general and we must use -some approximative/numerical method to compute the minimum. + +
import matplotlib.pyplot as plt
+import numpy as np
+from sklearn.model_selection import train_test_split
+from sklearn.datasets import load_breast_cancer
+from sklearn.linear_model import LogisticRegression
+
+# Load the data
+cancer = load_breast_cancer()
+
+X_train, X_test, y_train, y_test = train_test_split(cancer.data,cancer.target,random_state=0)
+print(X_train.shape)
+print(X_test.shape)
+# Logistic Regression
+logreg = LogisticRegression(solver='lbfgs')
+logreg.fit(X_train, y_train)
+print("Test set accuracy with Logistic Regression: {:.2f}".format(logreg.score(X_test,y_test)))
+#now scale the data
+from sklearn.preprocessing import StandardScaler
+scaler = StandardScaler()
+scaler.fit(X_train)
+X_train_scaled = scaler.transform(X_train)
+X_test_scaled = scaler.transform(X_test)
+# Logistic Regression
+logreg.fit(X_train_scaled, y_train)
+print("Test set accuracy Logistic Regression with scaled data: {:.2f}".format(logreg.score(X_test_scaled,y_test)))
+
@@ -436,7 +471,7 @@ some approximative/numerical method to compute the minimum.
-In our discussion on Logistic Regression we studied the -case of -two classes, with \( y_i \) either -\( 0 \) or \( 1 \). Furthermore we assumed also that we have only two -parameters \( \beta \) in our fitting, that is we -defined probabilities +In addition to the above scores, we could also study the covariance (and the correlation matrix). +We use Pandas to compute the correlation matrix. +
-$$ -\begin{align*} -p(y_i=1|x_i,\boldsymbol{\beta}) &= \frac{\exp{(\beta_0+\beta_1x_i)}}{1+\exp{(\beta_0+\beta_1x_i)}},\nonumber\\ -p(y_i=0|x_i,\boldsymbol{\beta}) &= 1 - p(y_i=1|x_i,\boldsymbol{\beta}), -\end{align*} -$$ + +
import matplotlib.pyplot as plt
+import numpy as np
+from sklearn.model_selection import train_test_split
+from sklearn.datasets import load_breast_cancer
+from sklearn.linear_model import LogisticRegression
+cancer = load_breast_cancer()
+import pandas as pd
+# Making a data frame
+cancerpd = pd.DataFrame(cancer.data, columns=cancer.feature_names)
-where \( \boldsymbol{\beta} \) are the weights we wish to extract from data, in our case \( \beta_0 \) and \( \beta_1 \).
+fig, axes = plt.subplots(15,2,figsize=(10,20))
+malignant = cancer.data[cancer.target == 0]
+benign = cancer.data[cancer.target == 1]
+ax = axes.ravel()
+for i in range(30):
+ _, bins = np.histogram(cancer.data[:,i], bins =50)
+ ax[i].hist(malignant[:,i], bins = bins, alpha = 0.5)
+ ax[i].hist(benign[:,i], bins = bins, alpha = 0.5)
+ ax[i].set_title(cancer.feature_names[i])
+ ax[i].set_yticks(())
+ax[0].set_xlabel("Feature magnitude")
+ax[0].set_ylabel("Frequency")
+ax[0].legend(["Malignant", "Benign"], loc ="best")
+fig.tight_layout()
+plt.show()
+
+import seaborn as sns
+correlation_matrix = cancerpd.corr().round(1)
+# use the heatmap function from seaborn to plot the correlation matrix
+# annot = True to print the values inside the square
+plt.figure(figsize=(15,8))
+sns.heatmap(data=correlation_matrix, annot=True)
+plt.show()
+
@@ -440,7 +478,7 @@ where \( \boldsymbol{\beta} \) are the weights we wish to extract from data, in
-
@@ -438,7 +452,7 @@ MathJax.Hub.Config({
-
@@ -437,57 +437,165 @@ set up the design matrix, what scaling we should use and other topics which may confuse us.
-Yes, it could be a bad idea to include the intercept column for the -exact reason you stated. If no transformation is applied to your data, -the intercept can be interpreted as the expected value of your target -variable when all your predictors are put to zero. Therefore, whenever -you cannot assume that the expected target variable is zero when all -your predictors are zero, it could be a bad idea to apply a model -which penalizes the intercept. Also, the analytical solution to the -ridge regression coefficients (when not shrinking
+The intercept can be interpreted as the expected value of our
+target/output variables when all other predictors are set to zero.
+Thus, if we cannot assume that the expected outputs/targets are zero
+when all predictors are zero (the columns in the design matrix), it
+may be a bad idea to implement a model which penalizes the intercept.
+Furthermore, in for example Ridge and Lasso regression, the solutions
+(when not shrinking
$$\beta_0$$
-
) is
-derived under the assumption that both y and X are zero centered (mean
-subtracted). What you are doing is correct, but you should also zero
-center X (subtracting the mean of each column from the corresponding
-column).
+
-If your predictors are of different scales, I would advice you to
-standardize X by subtracting the mean of each column from the
-corresponding column and dividing the column with its standard
-deviation. If you dont do this, you will give an "unfair" penalization
-of the parameters since their magnitude depends on the scale of their
-corresponding predictor. Suppose that you have an input variable
-"height". Human height might be measured in inches or meters or
+If our predictors represent different scales, then it is important to
+standardize the design matrix \( \boldsymbol{X} \) by subtracting the mean of each
+column from the corresponding column and dividing the column with its
+standard deviation.
+
+
+The
+Standadscaler
+function in Scikit-Learn does this for us. For the data sets we
+have been studying in our various examples, the data are in many cases
+already scaled and there is no need to scale them.
+
+
+If you need to scale the data, not doing so will give an unfair
+penalization of the parameters since their magnitude depends on the
+scale of their corresponding predictor.
+
+
+Suppose as an example that you
+you have an input variable given by the heights of different persons.
+Human height might be measured in inches or meters or
kilometers. If measured in kilometers, a standard linear regression
model with this predictor would probably give a much bigger
-coefficient term, than if measured in millimeters. You may see how
-this could become a problem when considering the loss function for
-ridge regression.
+coefficient term, than if measured in millimeters.
+This can clearly lead to problems in evaluating the cost/loss functions.
+
-Remember that when you do any transformation to your dataset before
-training, the exact same transformation has to be applied to new data
-before making a prediction. In your case, this means:
+Keep in mind that when you transform your data set before training a model, the same transformation needs to be done
+on your eventual new data set before making a prediction. If we translate this into a Python code, it would could be implemented as follows
-
+This code shows a simple first-order fit to a data set using the above transformed data, where we consider the role of the intercept first, by either excluding it or including it (code example thanks to Øyvind Sigmundson Schøyen)
+
+
+
+
+
@@ -784,80 +892,6 @@ plt.plot(np.log10(lambdas), MSERidgePredict, 'g
plt.xlabel('log10(lambda)')
plt.ylabel('MSE')
plt.legend()
-plt.show()
-
-
-
-
) for the unknown parameters
+\( \boldsymbol{\beta} \) are derived under the assumption that both \( \boldsymbol{y} \) and
+\( \boldsymbol{X} \) are zero centered, that is we subtract the mean values.
+
+
+
+More thinking
Still thinking
#Model training:
+
#Model training, we compute the mean value of y and X
y_train_mean = np.mean(y_train)
X_train_mean = np.mean(X_train,axis=0)
X_train = X_train - X_train_mean
y_train = y_train - y_train_mean
+# The we fit our model with the training data
trained_model = some_model.fit(X_train,y_train)
-#Model prediction:
+
+#Model prediction, here we need also to transform our data set used for the prediction.
X_test = X_test - X_train_mean #Use mean from training data
y_pred = trained_model(X_test)
y_pred = y_pred + y_train_mean
Linear Regression code, Intercept handling first
+
import numpy as np
+import matplotlib.pyplot as plt
+
+from sklearn.linear_model import LinearRegression
+
+
+np.random.seed(2021)
+
+
+def fit_beta(X, y):
+ return np.linalg.pinv(X.T @ X) @ X.T @ y
+
+
+true_beta = [2, 0.5, 3.7]
+
+x = np.linspace(0, 1, 11)
+y = np.sum(
+ np.asarray([x ** p * b for p, b in enumerate(true_beta)]), axis=0
+) + 0.1 * np.random.normal(size=len(x))
+
+degree = 3
+X = np.zeros((len(x), degree))
+
+# Include the intercept in the design matrix
+for p in range(degree):
+ X[:, p] = x ** p
+
+beta = fit_beta(X, y)
+
+# Intercept is included in the design matrix
+clf = LinearRegression(fit_intercept=False).fit(X, y)
+
+print(f"True beta: {true_beta}")
+print(f"Fitted beta: {beta}")
+print(f"Sklearn fitted beta: {clf.coef_}")
+
+
+plt.figure()
+plt.scatter(x, y, label="Data")
+plt.plot(x, X @ beta, label="Fit")
+plt.plot(x, clf.predict(X), label="Sklearn (fit_intercept=False)")
+
+
+# Do not include the intercept in the design matrix
+X = np.zeros((len(x), degree - 1))
+
+for p in range(degree - 1):
+ X[:, p] = x ** (p + 1)
+
+# Intercept is not included in the design matrix
+clf = LinearRegression(fit_intercept=True).fit(X, y)
+
+# Use centered values for X and y when computing coefficients
+y_offset = np.average(y, axis=0)
+X_offset = np.average(X, axis=0)
+
+beta = fit_beta(X - X_offset, y - y_offset)
+intercept = np.mean(y_offset - X_offset @ beta)
+
+print(f"Manual intercept: {intercept}")
+print(f"Fitted beta (sans intercept): {beta}")
+print(f"Sklearn intercept: {clf.intercept_}")
+print(f"Sklearn fitted beta (sans intercept): {clf.coef_}")
+
+plt.plot(x, X @ beta + intercept, "--", label="Fit (manual intercept)")
+plt.plot(x, clf.predict(X), "--", label="Sklearn (fit_intercept=True)")
+plt.grid()
+plt.legend()
+
+plt.show()
+
What does centering mean mathematically?
Here is a mathematical explanation of the zero centering:
import numpy as np
-import matplotlib.pyplot as plt
-
-from sklearn.linear_model import LinearRegression
-
-
-np.random.seed(2021)
-
-
-def fit_beta(X, y):
- return np.linalg.pinv(X.T @ X) @ X.T @ y
-
-
-true_beta = [2, 0.5, 3.7]
-
-x = np.linspace(0, 1, 11)
-y = np.sum(
- np.asarray([x ** p * b for p, b in enumerate(true_beta)]), axis=0
-) + 0.1 * np.random.normal(size=len(x))
-
-degree = 3
-X = np.zeros((len(x), degree))
-
-# Include the intercept in the design matrix
-for p in range(degree):
- X[:, p] = x ** p
-
-beta = fit_beta(X, y)
-
-# Intercept is included in the design matrix
-clf = LinearRegression(fit_intercept=False).fit(X, y)
-
-print(f"True beta: {true_beta}")
-print(f"Fitted beta: {beta}")
-print(f"Sklearn fitted beta: {clf.coef_}")
-
-
-plt.figure()
-plt.scatter(x, y, label="Data")
-plt.plot(x, X @ beta, label="Fit")
-plt.plot(x, clf.predict(X), label="Sklearn (fit_intercept=False)")
-
-
-# Do not include the intercept in the design matrix
-X = np.zeros((len(x), degree - 1))
-
-for p in range(degree - 1):
- X[:, p] = x ** (p + 1)
-
-# Intercept is not included in the design matrix
-clf = LinearRegression(fit_intercept=True).fit(X, y)
-
-# Use centered values for X and y when computing coefficients
-y_offset = np.average(y, axis=0)
-X_offset = np.average(X, axis=0)
-
-beta = fit_beta(X - X_offset, y - y_offset)
-intercept = np.mean(y_offset - X_offset @ beta)
-
-print(f"Manual intercept: {intercept}")
-print(f"Fitted beta (sans intercept): {beta}")
-print(f"Sklearn intercept: {clf.intercept_}")
-print(f"Sklearn fitted beta (sans intercept): {clf.coef_}")
-
-plt.plot(x, X @ beta + intercept, "--", label="Fit (manual intercept)")
-plt.plot(x, clf.predict(X), "--", label="Sklearn (fit_intercept=True)")
-plt.grid()
-plt.legend()
-
plt.show()
-
@@ -587,55 +597,161 @@ set up the design matrix, what scaling we should use and other topics
which may confuse us.
-Yes, it could be a bad idea to include the intercept column for the -exact reason you stated. If no transformation is applied to your data, -the intercept can be interpreted as the expected value of your target -variable when all your predictors are put to zero. Therefore, whenever -you cannot assume that the expected target variable is zero when all -your predictors are zero, it could be a bad idea to apply a model -which penalizes the intercept. Also, the analytical solution to the -ridge regression coefficients (when not shrinking $$\beta_0$$) is -derived under the assumption that both y and X are zero centered (mean -subtracted). What you are doing is correct, but you should also zero -center X (subtracting the mean of each column from the corresponding -column). +The intercept can be interpreted as the expected value of our +target/output variables when all other predictors are set to zero. +Thus, if we cannot assume that the expected outputs/targets are zero +when all predictors are zero (the columns in the design matrix), it +may be a bad idea to implement a model which penalizes the intercept. +Furthermore, in for example Ridge and Lasso regression, the solutions +(when not shrinking $$\beta_0$$) for the unknown parameters +\( \boldsymbol{\beta} \) are derived under the assumption that both \( \boldsymbol{y} \) and +\( \boldsymbol{X} \) are zero centered, that is we subtract the mean values.
-If your predictors are of different scales, I would advice you to
-standardize X by subtracting the mean of each column from the
-corresponding column and dividing the column with its standard
-deviation. If you dont do this, you will give an "unfair" penalization
-of the parameters since their magnitude depends on the scale of their
-corresponding predictor. Suppose that you have an input variable
-"height". Human height might be measured in inches or meters or
+
+
+
+If our predictors represent different scales, then it is important to +standardize the design matrix \( \boldsymbol{X} \) by subtracting the mean of each +column from the corresponding column and dividing the column with its +standard deviation. + +
+The +Standadscaler +function in Scikit-Learn does this for us. For the data sets we +have been studying in our various examples, the data are in many cases +already scaled and there is no need to scale them. + +
+If you need to scale the data, not doing so will give an unfair +penalization of the parameters since their magnitude depends on the +scale of their corresponding predictor. + +
+Suppose as an example that you +you have an input variable given by the heights of different persons. +Human height might be measured in inches or meters or kilometers. If measured in kilometers, a standard linear regression model with this predictor would probably give a much bigger -coefficient term, than if measured in millimeters. You may see how -this could become a problem when considering the loss function for -ridge regression. +coefficient term, than if measured in millimeters. +This can clearly lead to problems in evaluating the cost/loss functions.
-Remember that when you do any transformation to your dataset before
-training, the exact same transformation has to be applied to new data
-before making a prediction. In your case, this means:
+
+
+
+Keep in mind that when you transform your data set before training a model, the same transformation needs to be done +on your eventual new data set before making a prediction. If we translate this into a Python code, it would could be implemented as follows
-
#Model training:
+#Model training, we compute the mean value of y and X
y_train_mean = np.mean(y_train)
X_train_mean = np.mean(X_train,axis=0)
X_train = X_train - X_train_mean
y_train = y_train - y_train_mean
+# The we fit our model with the training data
trained_model = some_model.fit(X_train,y_train)
-#Model prediction:
+
+#Model prediction, here we need also to transform our data set used for the prediction.
X_test = X_test - X_train_mean #Use mean from training data
y_pred = trained_model(X_test)
y_pred = y_pred + y_train_mean
+
+
+
Linear Regression code, Intercept handling first
+
+
+This code shows a simple first-order fit to a data set using the above transformed data, where we consider the role of the intercept first, by either excluding it or including it (code example thanks to Øyvind Sigmundson Schøyen)
+
+
+
+
+
import numpy as np
+import matplotlib.pyplot as plt
+
+from sklearn.linear_model import LinearRegression
+
+
+np.random.seed(2021)
+
+
+def fit_beta(X, y):
+ return np.linalg.pinv(X.T @ X) @ X.T @ y
+
+
+true_beta = [2, 0.5, 3.7]
+
+x = np.linspace(0, 1, 11)
+y = np.sum(
+ np.asarray([x ** p * b for p, b in enumerate(true_beta)]), axis=0
+) + 0.1 * np.random.normal(size=len(x))
+
+degree = 3
+X = np.zeros((len(x), degree))
+
+# Include the intercept in the design matrix
+for p in range(degree):
+ X[:, p] = x ** p
+
+beta = fit_beta(X, y)
+
+# Intercept is included in the design matrix
+clf = LinearRegression(fit_intercept=False).fit(X, y)
+
+print(f"True beta: {true_beta}")
+print(f"Fitted beta: {beta}")
+print(f"Sklearn fitted beta: {clf.coef_}")
+
+
+plt.figure()
+plt.scatter(x, y, label="Data")
+plt.plot(x, X @ beta, label="Fit")
+plt.plot(x, clf.predict(X), label="Sklearn (fit_intercept=False)")
+
+
+# Do not include the intercept in the design matrix
+X = np.zeros((len(x), degree - 1))
+
+for p in range(degree - 1):
+ X[:, p] = x ** (p + 1)
+
+# Intercept is not included in the design matrix
+clf = LinearRegression(fit_intercept=True).fit(X, y)
+
+# Use centered values for X and y when computing coefficients
+y_offset = np.average(y, axis=0)
+X_offset = np.average(X, axis=0)
+
+beta = fit_beta(X - X_offset, y - y_offset)
+intercept = np.mean(y_offset - X_offset @ beta)
+
+print(f"Manual intercept: {intercept}")
+print(f"Fitted beta (sans intercept): {beta}")
+print(f"Sklearn intercept: {clf.intercept_}")
+print(f"Sklearn fitted beta (sans intercept): {clf.coef_}")
+
+plt.plot(x, X @ beta + intercept, "--", label="Fit (manual intercept)")
+plt.plot(x, clf.predict(X), "--", label="Sklearn (fit_intercept=True)")
+plt.grid()
+plt.legend()
+
+plt.show()
+
+
+
+
+
What does centering mean mathematically?
Here is a mathematical explanation of the zero centering:
@@ -912,80 +1028,6 @@ plt.plot(np.log10(lambdas), MSERidgePredict, 'g
plt.xlabel('log10(lambda)')
plt.ylabel('MSE')
plt.legend()
-plt.show()
-
- - -
import numpy as np
-import matplotlib.pyplot as plt
-
-from sklearn.linear_model import LinearRegression
-
-
-np.random.seed(2021)
-
-
-def fit_beta(X, y):
- return np.linalg.pinv(X.T @ X) @ X.T @ y
-
-
-true_beta = [2, 0.5, 3.7]
-
-x = np.linspace(0, 1, 11)
-y = np.sum(
- np.asarray([x ** p * b for p, b in enumerate(true_beta)]), axis=0
-) + 0.1 * np.random.normal(size=len(x))
-
-degree = 3
-X = np.zeros((len(x), degree))
-
-# Include the intercept in the design matrix
-for p in range(degree):
- X[:, p] = x ** p
-
-beta = fit_beta(X, y)
-
-# Intercept is included in the design matrix
-clf = LinearRegression(fit_intercept=False).fit(X, y)
-
-print(f"True beta: {true_beta}")
-print(f"Fitted beta: {beta}")
-print(f"Sklearn fitted beta: {clf.coef_}")
-
-
-plt.figure()
-plt.scatter(x, y, label="Data")
-plt.plot(x, X @ beta, label="Fit")
-plt.plot(x, clf.predict(X), label="Sklearn (fit_intercept=False)")
-
-
-# Do not include the intercept in the design matrix
-X = np.zeros((len(x), degree - 1))
-
-for p in range(degree - 1):
- X[:, p] = x ** (p + 1)
-
-# Intercept is not included in the design matrix
-clf = LinearRegression(fit_intercept=True).fit(X, y)
-
-# Use centered values for X and y when computing coefficients
-y_offset = np.average(y, axis=0)
-X_offset = np.average(X, axis=0)
-
-beta = fit_beta(X - X_offset, y - y_offset)
-intercept = np.mean(y_offset - X_offset @ beta)
-
-print(f"Manual intercept: {intercept}")
-print(f"Fitted beta (sans intercept): {beta}")
-print(f"Sklearn intercept: {clf.intercept_}")
-print(f"Sklearn fitted beta (sans intercept): {clf.coef_}")
-
-plt.plot(x, X @ beta + intercept, "--", label="Fit (manual intercept)")
-plt.plot(x, clf.predict(X), "--", label="Sklearn (fit_intercept=True)")
-plt.grid()
-plt.legend()
-
plt.show()
diff --git a/doc/pub/week38/html/week38.html b/doc/pub/week38/html/week38.html index 2b01b38a6..d98d373ee 100644 --- a/doc/pub/week38/html/week38.html +++ b/doc/pub/week38/html/week38.html @@ -93,6 +93,16 @@ div { text-align: justify; text-justify: inter-word; } 2, None, 'to-think-about-first-part'), + ('More thinking', 2, None, 'more-thinking'), + ('Still thinking', 2, None, 'still-thinking'), + ('Linear Regression code, Intercept handling first', + 2, + None, + 'linear-regression-code-intercept-handling-first'), + ('What does centering mean mathematically?', + 2, + None, + 'what-does-centering-mean-mathematically'), ('More complicated Example: The Ising model', 2, None, @@ -327,7 +337,7 @@ MathJax.Hub.Config({
-
@@ -592,55 +602,161 @@ set up the design matrix, what scaling we should use and other topics
which may confuse us.
-Yes, it could be a bad idea to include the intercept column for the -exact reason you stated. If no transformation is applied to your data, -the intercept can be interpreted as the expected value of your target -variable when all your predictors are put to zero. Therefore, whenever -you cannot assume that the expected target variable is zero when all -your predictors are zero, it could be a bad idea to apply a model -which penalizes the intercept. Also, the analytical solution to the -ridge regression coefficients (when not shrinking $$\beta_0$$) is -derived under the assumption that both y and X are zero centered (mean -subtracted). What you are doing is correct, but you should also zero -center X (subtracting the mean of each column from the corresponding -column). +The intercept can be interpreted as the expected value of our +target/output variables when all other predictors are set to zero. +Thus, if we cannot assume that the expected outputs/targets are zero +when all predictors are zero (the columns in the design matrix), it +may be a bad idea to implement a model which penalizes the intercept. +Furthermore, in for example Ridge and Lasso regression, the solutions +(when not shrinking $$\beta_0$$) for the unknown parameters +\( \boldsymbol{\beta} \) are derived under the assumption that both \( \boldsymbol{y} \) and +\( \boldsymbol{X} \) are zero centered, that is we subtract the mean values.
-If your predictors are of different scales, I would advice you to
-standardize X by subtracting the mean of each column from the
-corresponding column and dividing the column with its standard
-deviation. If you dont do this, you will give an "unfair" penalization
-of the parameters since their magnitude depends on the scale of their
-corresponding predictor. Suppose that you have an input variable
-"height". Human height might be measured in inches or meters or
+
+
+
+If our predictors represent different scales, then it is important to +standardize the design matrix \( \boldsymbol{X} \) by subtracting the mean of each +column from the corresponding column and dividing the column with its +standard deviation. + +
+The +Standadscaler +function in Scikit-Learn does this for us. For the data sets we +have been studying in our various examples, the data are in many cases +already scaled and there is no need to scale them. + +
+If you need to scale the data, not doing so will give an unfair +penalization of the parameters since their magnitude depends on the +scale of their corresponding predictor. + +
+Suppose as an example that you +you have an input variable given by the heights of different persons. +Human height might be measured in inches or meters or kilometers. If measured in kilometers, a standard linear regression model with this predictor would probably give a much bigger -coefficient term, than if measured in millimeters. You may see how -this could become a problem when considering the loss function for -ridge regression. +coefficient term, than if measured in millimeters. +This can clearly lead to problems in evaluating the cost/loss functions.
-Remember that when you do any transformation to your dataset before
-training, the exact same transformation has to be applied to new data
-before making a prediction. In your case, this means:
+
+
+
+Keep in mind that when you transform your data set before training a model, the same transformation needs to be done +on your eventual new data set before making a prediction. If we translate this into a Python code, it would could be implemented as follows
-
#Model training:
+#Model training, we compute the mean value of y and X
y_train_mean = np.mean(y_train)
X_train_mean = np.mean(X_train,axis=0)
X_train = X_train - X_train_mean
y_train = y_train - y_train_mean
+# The we fit our model with the training data
trained_model = some_model.fit(X_train,y_train)
-#Model prediction:
+
+#Model prediction, here we need also to transform our data set used for the prediction.
X_test = X_test - X_train_mean #Use mean from training data
y_pred = trained_model(X_test)
y_pred = y_pred + y_train_mean
+
+
+
Linear Regression code, Intercept handling first
+
+
+This code shows a simple first-order fit to a data set using the above transformed data, where we consider the role of the intercept first, by either excluding it or including it (code example thanks to Øyvind Sigmundson Schøyen)
+
+
+
+
+
import numpy as np
+import matplotlib.pyplot as plt
+
+from sklearn.linear_model import LinearRegression
+
+
+np.random.seed(2021)
+
+
+def fit_beta(X, y):
+ return np.linalg.pinv(X.T @ X) @ X.T @ y
+
+
+true_beta = [2, 0.5, 3.7]
+
+x = np.linspace(0, 1, 11)
+y = np.sum(
+ np.asarray([x ** p * b for p, b in enumerate(true_beta)]), axis=0
+) + 0.1 * np.random.normal(size=len(x))
+
+degree = 3
+X = np.zeros((len(x), degree))
+
+# Include the intercept in the design matrix
+for p in range(degree):
+ X[:, p] = x ** p
+
+beta = fit_beta(X, y)
+
+# Intercept is included in the design matrix
+clf = LinearRegression(fit_intercept=False).fit(X, y)
+
+print(f"True beta: {true_beta}")
+print(f"Fitted beta: {beta}")
+print(f"Sklearn fitted beta: {clf.coef_}")
+
+
+plt.figure()
+plt.scatter(x, y, label="Data")
+plt.plot(x, X @ beta, label="Fit")
+plt.plot(x, clf.predict(X), label="Sklearn (fit_intercept=False)")
+
+
+# Do not include the intercept in the design matrix
+X = np.zeros((len(x), degree - 1))
+
+for p in range(degree - 1):
+ X[:, p] = x ** (p + 1)
+
+# Intercept is not included in the design matrix
+clf = LinearRegression(fit_intercept=True).fit(X, y)
+
+# Use centered values for X and y when computing coefficients
+y_offset = np.average(y, axis=0)
+X_offset = np.average(X, axis=0)
+
+beta = fit_beta(X - X_offset, y - y_offset)
+intercept = np.mean(y_offset - X_offset @ beta)
+
+print(f"Manual intercept: {intercept}")
+print(f"Fitted beta (sans intercept): {beta}")
+print(f"Sklearn intercept: {clf.intercept_}")
+print(f"Sklearn fitted beta (sans intercept): {clf.coef_}")
+
+plt.plot(x, X @ beta + intercept, "--", label="Fit (manual intercept)")
+plt.plot(x, clf.predict(X), "--", label="Sklearn (fit_intercept=True)")
+plt.grid()
+plt.legend()
+
+plt.show()
+
+
+
+
+
What does centering mean mathematically?
Here is a mathematical explanation of the zero centering:
@@ -917,80 +1033,6 @@ plt.plot(np..xlabel('log10(lambda)')
plt.ylabel('MSE')
plt.legend()
-plt.show()
-
- - -
import numpy as np
-import matplotlib.pyplot as plt
-
-from sklearn.linear_model import LinearRegression
-
-
-np.random.seed(2021)
-
-
-def fit_beta(X, y):
- return np.linalg.pinv(X.T @ X) @ X.T @ y
-
-
-true_beta = [2, 0.5, 3.7]
-
-x = np.linspace(0, 1, 11)
-y = np.sum(
- np.asarray([x ** p * b for p, b in enumerate(true_beta)]), axis=0
-) + 0.1 * np.random.normal(size=len(x))
-
-degree = 3
-X = np.zeros((len(x), degree))
-
-# Include the intercept in the design matrix
-for p in range(degree):
- X[:, p] = x ** p
-
-beta = fit_beta(X, y)
-
-# Intercept is included in the design matrix
-clf = LinearRegression(fit_intercept=False).fit(X, y)
-
-print(f"True beta: {true_beta}")
-print(f"Fitted beta: {beta}")
-print(f"Sklearn fitted beta: {clf.coef_}")
-
-
-plt.figure()
-plt.scatter(x, y, label="Data")
-plt.plot(x, X @ beta, label="Fit")
-plt.plot(x, clf.predict(X), label="Sklearn (fit_intercept=False)")
-
-
-# Do not include the intercept in the design matrix
-X = np.zeros((len(x), degree - 1))
-
-for p in range(degree - 1):
- X[:, p] = x ** (p + 1)
-
-# Intercept is not included in the design matrix
-clf = LinearRegression(fit_intercept=True).fit(X, y)
-
-# Use centered values for X and y when computing coefficients
-y_offset = np.average(y, axis=0)
-X_offset = np.average(X, axis=0)
-
-beta = fit_beta(X - X_offset, y - y_offset)
-intercept = np.mean(y_offset - X_offset @ beta)
-
-print(f"Manual intercept: {intercept}")
-print(f"Fitted beta (sans intercept): {beta}")
-print(f"Sklearn intercept: {clf.intercept_}")
-print(f"Sklearn fitted beta (sans intercept): {clf.coef_}")
-
-plt.plot(x, X @ beta + intercept, "--", label="Fit (manual intercept)")
-plt.plot(x, clf.predict(X), "--", label="Sklearn (fit_intercept=True)")
-plt.grid()
-plt.legend()
-
plt.show()
diff --git a/doc/pub/week38/ipynb/ipynb-week38-src.tar.gz b/doc/pub/week38/ipynb/ipynb-week38-src.tar.gz index 19e55c41b..534605b92 100644 Binary files a/doc/pub/week38/ipynb/ipynb-week38-src.tar.gz and b/doc/pub/week38/ipynb/ipynb-week38-src.tar.gz differ diff --git a/doc/pub/week38/ipynb/week38.ipynb b/doc/pub/week38/ipynb/week38.ipynb index dc2d08352..29c201f3d 100644 --- a/doc/pub/week38/ipynb/week38.ipynb +++ b/doc/pub/week38/ipynb/week38.ipynb @@ -10,7 +10,7 @@ " \n", "**Morten Hjorth-Jensen**, Department of Physics, University of Oslo and Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University\n", "\n", - "Date: **Sep 21, 2021**\n", + "Date: **Sep 22, 2021**\n", "\n", "Copyright 1999-2021, Morten Hjorth-Jensen. Released under CC Attribution-NonCommercial 4.0 license\n", "\n", @@ -365,38 +365,48 @@ "\n", "\n", "\n", - "Yes, it could be a bad idea to include the intercept column for the\n", - "exact reason you stated. If no transformation is applied to your data,\n", - "the intercept can be interpreted as the expected value of your target\n", - "variable when all your predictors are put to zero. Therefore, whenever\n", - "you cannot assume that the expected target variable is zero when all\n", - "your predictors are zero, it could be a bad idea to apply a model\n", - "which penalizes the intercept. Also, the analytical solution to the\n", - "ridge regression coefficients (when not shrinking $$\\beta_0$$) is\n", - "derived under the assumption that both y and X are zero centered (mean\n", - "subtracted). What you are doing is correct, but you should also zero\n", - "center X (subtracting the mean of each column from the corresponding\n", - "column). \n", + "The intercept can be interpreted as the expected value of our\n", + "target/output variables when all other predictors are set to zero.\n", + "Thus, if we cannot assume that the expected outputs/targets are zero\n", + "when all predictors are zero (the columns in the design matrix), it\n", + "may be a bad idea to implement a model which penalizes the intercept.\n", + "Furthermore, in for example Ridge and Lasso regression, the solutions\n", + "(when not shrinking $$\\beta_0$$) for the unknown parameters\n", + "$\\boldsymbol{\\beta}$ are derived under the assumption that both $\\boldsymbol{y}$ and\n", + "$\\boldsymbol{X}$ are zero centered, that is we subtract the mean values.\n", "\n", - "If your predictors are of different scales, I would advice you to\n", - "standardize X by subtracting the mean of each column from the\n", - "corresponding column and dividing the column with its standard\n", - "deviation. If you dont do this, you will give an \"unfair\" penalization\n", - "of the parameters since their magnitude depends on the scale of their\n", - "corresponding predictor. Suppose that you have an input variable\n", - "\"height\". Human height might be measured in inches or meters or\n", + "\n", + "## More thinking\n", + "\n", + "\n", + "If our predictors represent different scales, then it is important to\n", + "standardize the design matrix $\\boldsymbol{X}$ by subtracting the mean of each\n", + "column from the corresponding column and dividing the column with its\n", + "standard deviation.\n", + "\n", + "The\n", + "[Standadscaler](https://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.StandardScaler.html)\n", + "function in **Scikit-Learn** does this for us. For the data sets we\n", + "have been studying in our various examples, the data are in many cases\n", + "already scaled and there is no need to scale them.\n", + "\n", + "If you need to scale the data, not doing so will give an *unfair*\n", + "penalization of the parameters since their magnitude depends on the\n", + "scale of their corresponding predictor.\n", + "\n", + "Suppose as an example that you \n", + "you have an input variable given by the heights of different persons.\n", + "Human height might be measured in inches or meters or\n", "kilometers. If measured in kilometers, a standard linear regression\n", "model with this predictor would probably give a much bigger\n", - "coefficient term, than if measured in millimeters. You may see how\n", - "this could become a problem when considering the loss function for\n", - "ridge regression.\n", + "coefficient term, than if measured in millimeters.\n", + "This can clearly lead to problems in evaluating the cost/loss functions.\n", "\n", "\n", + "## Still thinking\n", "\n", - "\n", - "Remember that when you do any transformation to your dataset before\n", - "training, the exact same transformation has to be applied to new data\n", - "before making a prediction. In your case, this means:" + "Keep in mind that when you transform your data set before training a model, the same transformation needs to be done\n", + "on your eventual new data set before making a prediction. If we translate this into a Python code, it would could be implemented as follows" ] }, { @@ -408,15 +418,17 @@ }, "outputs": [], "source": [ - "#Model training:\n", + "#Model training, we compute the mean value of y and X\n", "y_train_mean = np.mean(y_train)\n", "X_train_mean = np.mean(X_train,axis=0)\n", "X_train = X_train - X_train_mean\n", "y_train = y_train - y_train_mean\n", "\n", + "# The we fit our model with the training data\n", "trained_model = some_model.fit(X_train,y_train)\n", "\n", - "#Model prediction:\n", + "\n", + "#Model prediction, here we need also to transform our data set used for the prediction.\n", "X_test = X_test - X_train_mean #Use mean from training data\n", "y_pred = trained_model(X_test)\n", "y_pred = y_pred + y_train_mean" @@ -426,6 +438,97 @@ "cell_type": "markdown", "metadata": {}, "source": [ + "## Linear Regression code, Intercept handling first\n", + "\n", + "This code shows a simple first-order fit to a data set using the above transformed data, where we consider the role of the intercept first, by either excluding it or including it (*code example thanks to Øyvind Sigmundson Schøyen*)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false, + "editable": true + }, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "\n", + "from sklearn.linear_model import LinearRegression\n", + "\n", + "\n", + "np.random.seed(2021)\n", + "\n", + "\n", + "def fit_beta(X, y):\n", + " return np.linalg.pinv(X.T @ X) @ X.T @ y\n", + "\n", + "\n", + "true_beta = [2, 0.5, 3.7]\n", + "\n", + "x = np.linspace(0, 1, 11)\n", + "y = np.sum(\n", + " np.asarray([x ** p * b for p, b in enumerate(true_beta)]), axis=0\n", + ") + 0.1 * np.random.normal(size=len(x))\n", + "\n", + "degree = 3\n", + "X = np.zeros((len(x), degree))\n", + "\n", + "# Include the intercept in the design matrix\n", + "for p in range(degree):\n", + " X[:, p] = x ** p\n", + "\n", + "beta = fit_beta(X, y)\n", + "\n", + "# Intercept is included in the design matrix\n", + "clf = LinearRegression(fit_intercept=False).fit(X, y)\n", + "\n", + "print(f\"True beta: {true_beta}\")\n", + "print(f\"Fitted beta: {beta}\")\n", + "print(f\"Sklearn fitted beta: {clf.coef_}\")\n", + "\n", + "\n", + "plt.figure()\n", + "plt.scatter(x, y, label=\"Data\")\n", + "plt.plot(x, X @ beta, label=\"Fit\")\n", + "plt.plot(x, clf.predict(X), label=\"Sklearn (fit_intercept=False)\")\n", + "\n", + "\n", + "# Do not include the intercept in the design matrix\n", + "X = np.zeros((len(x), degree - 1))\n", + "\n", + "for p in range(degree - 1):\n", + " X[:, p] = x ** (p + 1)\n", + "\n", + "# Intercept is not included in the design matrix\n", + "clf = LinearRegression(fit_intercept=True).fit(X, y)\n", + "\n", + "# Use centered values for X and y when computing coefficients\n", + "y_offset = np.average(y, axis=0)\n", + "X_offset = np.average(X, axis=0)\n", + "\n", + "beta = fit_beta(X - X_offset, y - y_offset)\n", + "intercept = np.mean(y_offset - X_offset @ beta)\n", + "\n", + "print(f\"Manual intercept: {intercept}\")\n", + "print(f\"Fitted beta (sans intercept): {beta}\")\n", + "print(f\"Sklearn intercept: {clf.intercept_}\")\n", + "print(f\"Sklearn fitted beta (sans intercept): {clf.coef_}\")\n", + "\n", + "plt.plot(x, X @ beta + intercept, \"--\", label=\"Fit (manual intercept)\")\n", + "plt.plot(x, clf.predict(X), \"--\", label=\"Sklearn (fit_intercept=True)\")\n", + "plt.grid()\n", + "plt.legend()\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## What does centering mean mathematically?\n", "Here is a mathematical explanation of the zero centering:\n", "\n", "\n", @@ -832,87 +935,6 @@ "plt.show()" ] }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false, - "editable": true - }, - "outputs": [], - "source": [ - "import numpy as np\n", - "import matplotlib.pyplot as plt\n", - "\n", - "from sklearn.linear_model import LinearRegression\n", - "\n", - "\n", - "np.random.seed(2021)\n", - "\n", - "\n", - "def fit_beta(X, y):\n", - " return np.linalg.pinv(X.T @ X) @ X.T @ y\n", - "\n", - "\n", - "true_beta = [2, 0.5, 3.7]\n", - "\n", - "x = np.linspace(0, 1, 11)\n", - "y = np.sum(\n", - " np.asarray([x ** p * b for p, b in enumerate(true_beta)]), axis=0\n", - ") + 0.1 * np.random.normal(size=len(x))\n", - "\n", - "degree = 3\n", - "X = np.zeros((len(x), degree))\n", - "\n", - "# Include the intercept in the design matrix\n", - "for p in range(degree):\n", - " X[:, p] = x ** p\n", - "\n", - "beta = fit_beta(X, y)\n", - "\n", - "# Intercept is included in the design matrix\n", - "clf = LinearRegression(fit_intercept=False).fit(X, y)\n", - "\n", - "print(f\"True beta: {true_beta}\")\n", - "print(f\"Fitted beta: {beta}\")\n", - "print(f\"Sklearn fitted beta: {clf.coef_}\")\n", - "\n", - "\n", - "plt.figure()\n", - "plt.scatter(x, y, label=\"Data\")\n", - "plt.plot(x, X @ beta, label=\"Fit\")\n", - "plt.plot(x, clf.predict(X), label=\"Sklearn (fit_intercept=False)\")\n", - "\n", - "\n", - "# Do not include the intercept in the design matrix\n", - "X = np.zeros((len(x), degree - 1))\n", - "\n", - "for p in range(degree - 1):\n", - " X[:, p] = x ** (p + 1)\n", - "\n", - "# Intercept is not included in the design matrix\n", - "clf = LinearRegression(fit_intercept=True).fit(X, y)\n", - "\n", - "# Use centered values for X and y when computing coefficients\n", - "y_offset = np.average(y, axis=0)\n", - "X_offset = np.average(X, axis=0)\n", - "\n", - "beta = fit_beta(X - X_offset, y - y_offset)\n", - "intercept = np.mean(y_offset - X_offset @ beta)\n", - "\n", - "print(f\"Manual intercept: {intercept}\")\n", - "print(f\"Fitted beta (sans intercept): {beta}\")\n", - "print(f\"Sklearn intercept: {clf.intercept_}\")\n", - "print(f\"Sklearn fitted beta (sans intercept): {clf.coef_}\")\n", - "\n", - "plt.plot(x, X @ beta + intercept, \"--\", label=\"Fit (manual intercept)\")\n", - "plt.plot(x, clf.predict(X), \"--\", label=\"Sklearn (fit_intercept=True)\")\n", - "plt.grid()\n", - "plt.legend()\n", - "\n", - "plt.show()" - ] - }, { "cell_type": "markdown", "metadata": {}, diff --git a/doc/src/week38/week38.do.txt b/doc/src/week38/week38.do.txt index 965ada8e8..34f5c0c87 100644 --- a/doc/src/week38/week38.do.txt +++ b/doc/src/week38/week38.do.txt @@ -253,56 +253,149 @@ which may confuse us. -Yes, it could be a bad idea to include the intercept column for the -exact reason you stated. If no transformation is applied to your data, -the intercept can be interpreted as the expected value of your target -variable when all your predictors are put to zero. Therefore, whenever -you cannot assume that the expected target variable is zero when all -your predictors are zero, it could be a bad idea to apply a model -which penalizes the intercept. Also, the analytical solution to the -ridge regression coefficients (when not shrinking $$\beta_0$$) is -derived under the assumption that both y and X are zero centered (mean -subtracted). What you are doing is correct, but you should also zero -center X (subtracting the mean of each column from the corresponding -column). +The intercept can be interpreted as the expected value of our +target/output variables when all other predictors are set to zero. +Thus, if we cannot assume that the expected outputs/targets are zero +when all predictors are zero (the columns in the design matrix), it +may be a bad idea to implement a model which penalizes the intercept. +Furthermore, in for example Ridge and Lasso regression, the solutions +(when not shrinking $$\beta_0$$) for the unknown parameters +$\bm{\beta}$ are derived under the assumption that both $\bm{y}$ and +$\bm{X}$ are zero centered, that is we subtract the mean values. -If your predictors are of different scales, I would advice you to -standardize X by subtracting the mean of each column from the -corresponding column and dividing the column with its standard -deviation. If you dont do this, you will give an "unfair" penalization -of the parameters since their magnitude depends on the scale of their -corresponding predictor. Suppose that you have an input variable -"height". Human height might be measured in inches or meters or + +!split +===== More thinking ===== + + +If our predictors represent different scales, then it is important to +standardize the design matrix $\bm{X}$ by subtracting the mean of each +column from the corresponding column and dividing the column with its +standard deviation. + +The +"Standadscaler":"https://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.StandardScaler.html" +function in _Scikit-Learn_ does this for us. For the data sets we +have been studying in our various examples, the data are in many cases +already scaled and there is no need to scale them. + +If you need to scale the data, not doing so will give an *unfair* +penalization of the parameters since their magnitude depends on the +scale of their corresponding predictor. + +Suppose as an example that you +you have an input variable given by the heights of different persons. +Human height might be measured in inches or meters or kilometers. If measured in kilometers, a standard linear regression model with this predictor would probably give a much bigger -coefficient term, than if measured in millimeters. You may see how -this could become a problem when considering the loss function for -ridge regression. +coefficient term, than if measured in millimeters. +This can clearly lead to problems in evaluating the cost/loss functions. +!split +===== Still thinking ===== - -Remember that when you do any transformation to your dataset before -training, the exact same transformation has to be applied to new data -before making a prediction. In your case, this means: +Keep in mind that when you transform your data set before training a model, the same transformation needs to be done +on your eventual new data set before making a prediction. If we translate this into a Python code, it would could be implemented as follows !bc pycod -#Model training: +#Model training, we compute the mean value of y and X y_train_mean = np.mean(y_train) X_train_mean = np.mean(X_train,axis=0) X_train = X_train - X_train_mean y_train = y_train - y_train_mean +# The we fit our model with the training data trained_model = some_model.fit(X_train,y_train) -#Model prediction: + +#Model prediction, here we need also to transform our data set used for the prediction. X_test = X_test - X_train_mean #Use mean from training data y_pred = trained_model(X_test) y_pred = y_pred + y_train_mean !ec +!split +===== Linear Regression code, Intercept handling first ===== + +This code shows a simple first-order fit to a data set using the above transformed data, where we consider the role of the intercept first, by either excluding it or including it (*code example thanks to Øyvind Sigmundson Schøyen*) + +!bc pycod +import numpy as np +import matplotlib.pyplot as plt + +from sklearn.linear_model import LinearRegression +np.random.seed(2021) + + +def fit_beta(X, y): + return np.linalg.pinv(X.T @ X) @ X.T @ y + + +true_beta = [2, 0.5, 3.7] + +x = np.linspace(0, 1, 11) +y = np.sum( + np.asarray([x ** p * b for p, b in enumerate(true_beta)]), axis=0 +) + 0.1 * np.random.normal(size=len(x)) + +degree = 3 +X = np.zeros((len(x), degree)) + +# Include the intercept in the design matrix +for p in range(degree): + X[:, p] = x ** p + +beta = fit_beta(X, y) + +# Intercept is included in the design matrix +clf = LinearRegression(fit_intercept=False).fit(X, y) + +print(f"True beta: {true_beta}") +print(f"Fitted beta: {beta}") +print(f"Sklearn fitted beta: {clf.coef_}") + + +plt.figure() +plt.scatter(x, y, label="Data") +plt.plot(x, X @ beta, label="Fit") +plt.plot(x, clf.predict(X), label="Sklearn (fit_intercept=False)") + + +# Do not include the intercept in the design matrix +X = np.zeros((len(x), degree - 1)) + +for p in range(degree - 1): + X[:, p] = x ** (p + 1) + +# Intercept is not included in the design matrix +clf = LinearRegression(fit_intercept=True).fit(X, y) + +# Use centered values for X and y when computing coefficients +y_offset = np.average(y, axis=0) +X_offset = np.average(X, axis=0) + +beta = fit_beta(X - X_offset, y - y_offset) +intercept = np.mean(y_offset - X_offset @ beta) + +print(f"Manual intercept: {intercept}") +print(f"Fitted beta (sans intercept): {beta}") +print(f"Sklearn intercept: {clf.intercept_}") +print(f"Sklearn fitted beta (sans intercept): {clf.coef_}") + +plt.plot(x, X @ beta + intercept, "--", label="Fit (manual intercept)") +plt.plot(x, clf.predict(X), "--", label="Sklearn (fit_intercept=True)") +plt.grid() +plt.legend() + +plt.show() + +!ec + +!split +===== What does centering mean mathematically? ===== Here is a mathematical explanation of the zero centering: @@ -602,79 +695,6 @@ plt.show() !ec -!bc pycod -import numpy as np -import matplotlib.pyplot as plt - -from sklearn.linear_model import LinearRegression - - -np.random.seed(2021) - - -def fit_beta(X, y): - return np.linalg.pinv(X.T @ X) @ X.T @ y - - -true_beta = [2, 0.5, 3.7] - -x = np.linspace(0, 1, 11) -y = np.sum( - np.asarray([x ** p * b for p, b in enumerate(true_beta)]), axis=0 -) + 0.1 * np.random.normal(size=len(x)) - -degree = 3 -X = np.zeros((len(x), degree)) - -# Include the intercept in the design matrix -for p in range(degree): - X[:, p] = x ** p - -beta = fit_beta(X, y) - -# Intercept is included in the design matrix -clf = LinearRegression(fit_intercept=False).fit(X, y) - -print(f"True beta: {true_beta}") -print(f"Fitted beta: {beta}") -print(f"Sklearn fitted beta: {clf.coef_}") - - -plt.figure() -plt.scatter(x, y, label="Data") -plt.plot(x, X @ beta, label="Fit") -plt.plot(x, clf.predict(X), label="Sklearn (fit_intercept=False)") - - -# Do not include the intercept in the design matrix -X = np.zeros((len(x), degree - 1)) - -for p in range(degree - 1): - X[:, p] = x ** (p + 1) - -# Intercept is not included in the design matrix -clf = LinearRegression(fit_intercept=True).fit(X, y) - -# Use centered values for X and y when computing coefficients -y_offset = np.average(y, axis=0) -X_offset = np.average(X, axis=0) - -beta = fit_beta(X - X_offset, y - y_offset) -intercept = np.mean(y_offset - X_offset @ beta) - -print(f"Manual intercept: {intercept}") -print(f"Fitted beta (sans intercept): {beta}") -print(f"Sklearn intercept: {clf.intercept_}") -print(f"Sklearn fitted beta (sans intercept): {clf.coef_}") - -plt.plot(x, X @ beta + intercept, "--", label="Fit (manual intercept)") -plt.plot(x, clf.predict(X), "--", label="Sklearn (fit_intercept=True)") -plt.grid() -plt.legend() - -plt.show() - -!ec !split