+ +Recommended Reading: + +
First we study and compare the OLS and Ridge results. The next code compares all three methods. + +
+ + +
import os
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+
+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.
+
+X = np.array( [ [ 2, 0], [0, 1], [0,0]])
+y = np.array( [4, 2, 3])
+
+
+# matrix inversion to find beta
+OLSbeta = np.linalg.inv(X.T @ X) @ X.T @ y
+print(OLSbeta)
+# and then make the prediction
+ytildeOLS = X @ OLSbeta
+print("Training MSE for OLS")
+print(MSE(y,ytildeOLS))
+ypredictOLS = X @ OLSbeta
+
+# Repeat now for Ridge regression and various values of the regularization parameter
+I = np.eye(2,2)
+# Decide which values of lambda to use
+nlambdas = 100
+MSEPredict = np.zeros(nlambdas)
+lambdas = np.logspace(-4, 4, nlambdas)
+for i in range(nlambdas):
+ lmb = lambdas[i]
+ Ridgebeta = np.linalg.inv(X.T @ X+lmb*I) @ X.T @ y
+# print(Ridgebeta)
+ # and then make the prediction
+ ypredictRidge = X @ Ridgebeta
+ MSEPredict[i] = MSE(y,ypredictRidge)
+# print(MSEPredict[i])
+ # Now plot the results
+plt.figure()
+plt.plot(np.log10(lambdas), MSEPredict, 'r--', label = 'MSE Ridge Test')
+plt.xlabel('log10(lambda)')
+plt.ylabel('MSE')
+plt.legend()
+plt.show()
+
+We see here that we reach a plateau. What is actually happening?
+
+
+
+
+
+
+
+
+
+With Lasso Regression
+
+import os
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+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.
+
+X = np.array( [ [ 2, 0], [0, 1], [0,0]])
+y = np.array( [4, 2, 3])
+
+
+# matrix inversion to find beta
+OLSbeta = np.linalg.inv(X.T @ X) @ X.T @ y
+print(OLSbeta)
+# and then make the prediction
+ytildeOLS = X @ OLSbeta
+print("Training MSE for OLS")
+print(MSE(y,ytildeOLS))
+ypredictOLS = X @ OLSbeta
+
+# Repeat now for Ridge regression and various values of the regularization parameter
+I = np.eye(2,2)
+# Decide which values of lambda to use
+nlambdas = 100
+MSERidgePredict = np.zeros(nlambdas)
+MSELassoPredict = np.zeros(nlambdas)
+lambdas = np.logspace(-4, 4, nlambdas)
+for i in range(nlambdas):
+ lmb = lambdas[i]
+ Ridgebeta = np.linalg.inv(X.T @ X+lmb*I) @ X.T @ y
+# print(Ridgebeta)
+ # and then make the prediction
+ ypredictRidge = X @ Ridgebeta
+ MSERidgePredict[i] = MSE(y,ypredictRidge)
+# print(MSEPredict[i])
+ RegLasso = linear_model.Lasso(lmb)
+ RegLasso.fit(X,y)
+ ypredictLasso = RegLasso.predict(X)
+ MSELassoPredict[i] = MSE(y,ypredictLasso)
+# Now plot the results
+plt.figure()
+plt.plot(np.log10(lambdas), MSERidgePredict, 'r--', label = 'MSE Ridge Train')
+plt.plot(np.log10(lambdas), MSELassoPredict, 'r--', label = 'MSE Lasso Train')
+plt.xlabel('log10(lambda)')
+plt.ylabel('MSE')
+plt.legend()
+plt.show()
+
Another Example, now with a polynomial fit
+
+import os
+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)
+
+x = np.random.rand(100)
+y = 2.0+5*x*x+0.1*np.random.randn(100)
+
+# number of features p (here degree of polynomial
+p = 3
+# The design matrix now as function of a given polynomial
+X = np.zeros((len(x),p))
+X[:,0] = 1.0
+X[:,1] = x
+X[:,2] = x*x
+# 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.inv(X_train.T @ X_train) @ X_train.T @ y_train
+print(OLSbeta)
+# and then make the prediction
+ytildeOLS = X_train @ OLSbeta
+print("Training R2 for OLS")
+print(R2(y_train,ytildeOLS))
+print("Training MSE for OLS")
+print(MSE(y_train,ytildeOLS))
+ypredictOLS = X_test @ OLSbeta
+print("Test R2 for OLS")
+print(R2(y_test,ypredictOLS))
+print("Test MSE OLS")
+print(MSE(y_test,ypredictOLS))
+
+# Repeat now for Lasso and Ridge regression and various values of the regularization parameter
+I = np.eye(p,p)
+# Decide which values of lambda to use
+nlambdas = 100
+MSEPredict = np.zeros(nlambdas)
+MSETrain = np.zeros(nlambdas)
+MSELassoPredict = np.zeros(nlambdas)
+MSELassoTrain = np.zeros(nlambdas)
+lambdas = np.logspace(-4, 4, nlambdas)
+for i in range(nlambdas):
+ lmb = lambdas[i]
+ Ridgebeta = np.linalg.inv(X_train.T @ X_train+lmb*I) @ X_train.T @ y_train
+ # include lasso using Scikit-Learn
+ RegLasso = linear_model.Lasso(lmb)
+ RegLasso.fit(X_train,y_train)
+ # and then make the prediction
+ ytildeRidge = X_train @ Ridgebeta
+ ypredictRidge = X_test @ Ridgebeta
+ ytildeLasso = RegLasso.predict(X_train)
+ ypredictLasso = RegLasso.predict(X_test)
+ MSEPredict[i] = MSE(y_test,ypredictRidge)
+ MSETrain[i] = MSE(y_train,ytildeRidge)
+ MSELassoPredict[i] = MSE(y_test,ypredictLasso)
+ MSELassoTrain[i] = MSE(y_train,ytildeLasso)
+
+# Now plot the results
+plt.figure()
+plt.plot(np.log10(lambdas), MSETrain, label = 'MSE Ridge train')
+plt.plot(np.log10(lambdas), MSEPredict, 'r--', label = 'MSE Ridge Test')
+plt.plot(np.log10(lambdas), MSELassoTrain, label = 'MSE Lasso train')
+plt.plot(np.log10(lambdas), MSELassoPredict, 'r--', label = 'MSE Lasso Test')
+
+plt.xlabel('log10(lambda)')
+plt.ylabel('MSE')
+plt.legend()
+plt.show()
+
First we study and compare the OLS and Ridge results. The next code compares all three methods. +
+ + +
import os
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+
+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.
+
+X = np.array( [ [ 2, 0], [0, 1], [0,0]])
+y = np.array( [4, 2, 3])
+
+
+# matrix inversion to find beta
+OLSbeta = np.linalg.inv(X.T @ X) @ X.T @ y
+print(OLSbeta)
+# and then make the prediction
+ytildeOLS = X @ OLSbeta
+print("Training MSE for OLS")
+print(MSE(y,ytildeOLS))
+ypredictOLS = X @ OLSbeta
+
+# Repeat now for Ridge regression and various values of the regularization parameter
+I = np.eye(2,2)
+# Decide which values of lambda to use
+nlambdas = 100
+MSEPredict = np.zeros(nlambdas)
+lambdas = np.logspace(-4, 4, nlambdas)
+for i in range(nlambdas):
+ lmb = lambdas[i]
+ Ridgebeta = np.linalg.inv(X.T @ X+lmb*I) @ X.T @ y
+# print(Ridgebeta)
+ # and then make the prediction
+ ypredictRidge = X @ Ridgebeta
+ MSEPredict[i] = MSE(y,ypredictRidge)
+# print(MSEPredict[i])
+ # Now plot the results
+plt.figure()
+plt.plot(np.log10(lambdas), MSEPredict, 'r--', label = 'MSE Ridge Test')
+plt.xlabel('log10(lambda)')
+plt.ylabel('MSE')
+plt.legend()
+plt.show()
++We see here that we reach a plateau. What is actually happening? + +
+
+
+
+ + +
import os
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+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.
+
+X = np.array( [ [ 2, 0], [0, 1], [0,0]])
+y = np.array( [4, 2, 3])
+
+
+# matrix inversion to find beta
+OLSbeta = np.linalg.inv(X.T @ X) @ X.T @ y
+print(OLSbeta)
+# and then make the prediction
+ytildeOLS = X @ OLSbeta
+print("Training MSE for OLS")
+print(MSE(y,ytildeOLS))
+ypredictOLS = X @ OLSbeta
+
+# Repeat now for Ridge regression and various values of the regularization parameter
+I = np.eye(2,2)
+# Decide which values of lambda to use
+nlambdas = 100
+MSERidgePredict = np.zeros(nlambdas)
+MSELassoPredict = np.zeros(nlambdas)
+lambdas = np.logspace(-4, 4, nlambdas)
+for i in range(nlambdas):
+ lmb = lambdas[i]
+ Ridgebeta = np.linalg.inv(X.T @ X+lmb*I) @ X.T @ y
+# print(Ridgebeta)
+ # and then make the prediction
+ ypredictRidge = X @ Ridgebeta
+ MSERidgePredict[i] = MSE(y,ypredictRidge)
+# print(MSEPredict[i])
+ RegLasso = linear_model.Lasso(lmb)
+ RegLasso.fit(X,y)
+ ypredictLasso = RegLasso.predict(X)
+ MSELassoPredict[i] = MSE(y,ypredictLasso)
+# Now plot the results
+plt.figure()
+plt.plot(np.log10(lambdas), MSERidgePredict, 'r--', label = 'MSE Ridge Train')
+plt.plot(np.log10(lambdas), MSELassoPredict, 'r--', label = 'MSE Lasso Train')
+plt.xlabel('log10(lambda)')
+plt.ylabel('MSE')
+plt.legend()
+plt.show()
+
+
+
+
+ + +
import os
+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)
+
+x = np.random.rand(100)
+y = 2.0+5*x*x+0.1*np.random.randn(100)
+
+# number of features p (here degree of polynomial
+p = 3
+# The design matrix now as function of a given polynomial
+X = np.zeros((len(x),p))
+X[:,0] = 1.0
+X[:,1] = x
+X[:,2] = x*x
+# 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.inv(X_train.T @ X_train) @ X_train.T @ y_train
+print(OLSbeta)
+# and then make the prediction
+ytildeOLS = X_train @ OLSbeta
+print("Training R2 for OLS")
+print(R2(y_train,ytildeOLS))
+print("Training MSE for OLS")
+print(MSE(y_train,ytildeOLS))
+ypredictOLS = X_test @ OLSbeta
+print("Test R2 for OLS")
+print(R2(y_test,ypredictOLS))
+print("Test MSE OLS")
+print(MSE(y_test,ypredictOLS))
+
+# Repeat now for Lasso and Ridge regression and various values of the regularization parameter
+I = np.eye(p,p)
+# Decide which values of lambda to use
+nlambdas = 100
+MSEPredict = np.zeros(nlambdas)
+MSETrain = np.zeros(nlambdas)
+MSELassoPredict = np.zeros(nlambdas)
+MSELassoTrain = np.zeros(nlambdas)
+lambdas = np.logspace(-4, 4, nlambdas)
+for i in range(nlambdas):
+ lmb = lambdas[i]
+ Ridgebeta = np.linalg.inv(X_train.T @ X_train+lmb*I) @ X_train.T @ y_train
+ # include lasso using Scikit-Learn
+ RegLasso = linear_model.Lasso(lmb)
+ RegLasso.fit(X_train,y_train)
+ # and then make the prediction
+ ytildeRidge = X_train @ Ridgebeta
+ ypredictRidge = X_test @ Ridgebeta
+ ytildeLasso = RegLasso.predict(X_train)
+ ypredictLasso = RegLasso.predict(X_test)
+ MSEPredict[i] = MSE(y_test,ypredictRidge)
+ MSETrain[i] = MSE(y_train,ytildeRidge)
+ MSELassoPredict[i] = MSE(y_test,ypredictLasso)
+ MSELassoTrain[i] = MSE(y_train,ytildeLasso)
+
+# Now plot the results
+plt.figure()
+plt.plot(np.log10(lambdas), MSETrain, label = 'MSE Ridge train')
+plt.plot(np.log10(lambdas), MSEPredict, 'r--', label = 'MSE Ridge Test')
+plt.plot(np.log10(lambdas), MSELassoTrain, label = 'MSE Lasso train')
+plt.plot(np.log10(lambdas), MSELassoPredict, 'r--', label = 'MSE Lasso Test')
+
+plt.xlabel('log10(lambda)')
+plt.ylabel('MSE')
+plt.legend()
+plt.show()
+diff --git a/doc/pub/week36/html/week36.html b/doc/pub/week36/html/week36.html index 0b3d76a4b..ba6797338 100644 --- a/doc/pub/week36/html/week36.html +++ b/doc/pub/week36/html/week36.html @@ -123,6 +123,11 @@ div { text-align: justify; text-justify: inter-word; } 2, None, 'simple-code-for-solving-the-above-problem'), + ('With Lasso Regression', 2, None, 'with-lasso-regression'), + ('Another Example, now with a polynomial fit', + 2, + None, + 'another-example-now-with-a-polynomial-fit'), ('Linking the regression analysis with a statistical ' 'interpretation', 2, @@ -299,6 +304,14 @@ MathJax.Hub.Config({
First we study and compare the OLS and Ridge results. The next code compares all three methods. +
+ + +
import os
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+
+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.
+
+X = np.array( [ [ 2, 0], [0, 1], [0,0]])
+y = np.array( [4, 2, 3])
+
+
+# matrix inversion to find beta
+OLSbeta = np.linalg.inv(X.T @ X) @ X.T @ y
+print(OLSbeta)
+# and then make the prediction
+ytildeOLS = X @ OLSbeta
+print("Training MSE for OLS")
+print(MSE(y,ytildeOLS))
+ypredictOLS = X @ OLSbeta
+
+# Repeat now for Ridge regression and various values of the regularization parameter
+I = np.eye(2,2)
+# Decide which values of lambda to use
+nlambdas = 100
+MSEPredict = np.zeros(nlambdas)
+lambdas = np.logspace(-4, 4, nlambdas)
+for i in range(nlambdas):
+ lmb = lambdas[i]
+ Ridgebeta = np.linalg.inv(X.T @ X+lmb*I) @ X.T @ y
+# print(Ridgebeta)
+ # and then make the prediction
+ ypredictRidge = X @ Ridgebeta
+ MSEPredict[i] = MSE(y,ypredictRidge)
+# print(MSEPredict[i])
+ # Now plot the results
+plt.figure()
+plt.plot(np.log10(lambdas), MSEPredict, 'r--', label = 'MSE Ridge Test')
+plt.xlabel('log10(lambda)')
+plt.ylabel('MSE')
+plt.legend()
+plt.show()
++We see here that we reach a plateau. What is actually happening? + +
+
+
+
+ + +
import os
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+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.
+
+X = np.array( [ [ 2, 0], [0, 1], [0,0]])
+y = np.array( [4, 2, 3])
+
+
+# matrix inversion to find beta
+OLSbeta = np.linalg.inv(X.T @ X) @ X.T @ y
+print(OLSbeta)
+# and then make the prediction
+ytildeOLS = X @ OLSbeta
+print("Training MSE for OLS")
+print(MSE(y,ytildeOLS))
+ypredictOLS = X @ OLSbeta
+
+# Repeat now for Ridge regression and various values of the regularization parameter
+I = np.eye(2,2)
+# Decide which values of lambda to use
+nlambdas = 100
+MSERidgePredict = np.zeros(nlambdas)
+MSELassoPredict = np.zeros(nlambdas)
+lambdas = np.logspace(-4, 4, nlambdas)
+for i in range(nlambdas):
+ lmb = lambdas[i]
+ Ridgebeta = np.linalg.inv(X.T @ X+lmb*I) @ X.T @ y
+# print(Ridgebeta)
+ # and then make the prediction
+ ypredictRidge = X @ Ridgebeta
+ MSERidgePredict[i] = MSE(y,ypredictRidge)
+# print(MSEPredict[i])
+ RegLasso = linear_model.Lasso(lmb)
+ RegLasso.fit(X,y)
+ ypredictLasso = RegLasso.predict(X)
+ MSELassoPredict[i] = MSE(y,ypredictLasso)
+# Now plot the results
+plt.figure()
+plt.plot(np.log10(lambdas), MSERidgePredict, 'r--', label = 'MSE Ridge Train')
+plt.plot(np.log10(lambdas), MSELassoPredict, 'r--', label = 'MSE Lasso Train')
+plt.xlabel('log10(lambda)')
+plt.ylabel('MSE')
+plt.legend()
+plt.show()
+
+
+
+
+ + +
import os
+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)
+
+x = np.random.rand(100)
+y = 2.0+5*x*x+0.1*np.random.randn(100)
+
+# number of features p (here degree of polynomial
+p = 3
+# The design matrix now as function of a given polynomial
+X = np.zeros((len(x),p))
+X[:,0] = 1.0
+X[:,1] = x
+X[:,2] = x*x
+# 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.inv(X_train.T @ X_train) @ X_train.T @ y_train
+print(OLSbeta)
+# and then make the prediction
+ytildeOLS = X_train @ OLSbeta
+print("Training R2 for OLS")
+print(R2(y_train,ytildeOLS))
+print("Training MSE for OLS")
+print(MSE(y_train,ytildeOLS))
+ypredictOLS = X_test @ OLSbeta
+print("Test R2 for OLS")
+print(R2(y_test,ypredictOLS))
+print("Test MSE OLS")
+print(MSE(y_test,ypredictOLS))
+
+# Repeat now for Lasso and Ridge regression and various values of the regularization parameter
+I = np.eye(p,p)
+# Decide which values of lambda to use
+nlambdas = 100
+MSEPredict = np.zeros(nlambdas)
+MSETrain = np.zeros(nlambdas)
+MSELassoPredict = np.zeros(nlambdas)
+MSELassoTrain = np.zeros(nlambdas)
+lambdas = np.logspace(-4, 4, nlambdas)
+for i in range(nlambdas):
+ lmb = lambdas[i]
+ Ridgebeta = np.linalg.inv(X_train.T @ X_train+lmb*I) @ X_train.T @ y_train
+ # include lasso using Scikit-Learn
+ RegLasso = linear_model.Lasso(lmb)
+ RegLasso.fit(X_train,y_train)
+ # and then make the prediction
+ ytildeRidge = X_train @ Ridgebeta
+ ypredictRidge = X_test @ Ridgebeta
+ ytildeLasso = RegLasso.predict(X_train)
+ ypredictLasso = RegLasso.predict(X_test)
+ MSEPredict[i] = MSE(y_test,ypredictRidge)
+ MSETrain[i] = MSE(y_train,ytildeRidge)
+ MSELassoPredict[i] = MSE(y_test,ypredictLasso)
+ MSELassoTrain[i] = MSE(y_train,ytildeLasso)
+
+# Now plot the results
+plt.figure()
+plt.plot(np.log10(lambdas), MSETrain, label = 'MSE Ridge train')
+plt.plot(np.log10(lambdas), MSEPredict, 'r--', label = 'MSE Ridge Test')
+plt.plot(np.log10(lambdas), MSELassoTrain, label = 'MSE Lasso train')
+plt.plot(np.log10(lambdas), MSELassoPredict, 'r--', label = 'MSE Lasso Test')
+
+plt.xlabel('log10(lambda)')
+plt.ylabel('MSE')
+plt.legend()
+plt.show()
+diff --git a/doc/pub/week36/ipynb/ipynb-week36-src.tar.gz b/doc/pub/week36/ipynb/ipynb-week36-src.tar.gz index 264bd23e0..c07b5e271 100644 Binary files a/doc/pub/week36/ipynb/ipynb-week36-src.tar.gz and b/doc/pub/week36/ipynb/ipynb-week36-src.tar.gz differ diff --git a/doc/pub/week36/ipynb/week36.ipynb b/doc/pub/week36/ipynb/week36.ipynb index e5e167cdd..9395e12c4 100644 --- a/doc/pub/week36/ipynb/week36.ipynb +++ b/doc/pub/week36/ipynb/week36.ipynb @@ -23,6 +23,13 @@ "\n", "* Friday: Linear Regression and links with Statistics, Resampling methods and presentation of first project.\n", "\n", + "Recommended Reading:\n", + "1. Lectures on Regression\n", + "\n", + "2. Bishop 1.1, 1.2, 2.1, 2.2, 2.3 and 3.1\n", + "\n", + "3. Hastie et al chapter 3\n", + "\n", "## Thursday September 9\n", "\n", "\n", @@ -1092,9 +1099,252 @@ "\n", "Here we set up the OLS, Ridge and Lasso functionality in order to study the above example. Note that here we have opted for a set of values of $\\lambda$, meaning that we need to perform a search in order to find the optimal values.\n", "\n", - "First we study and compare the OLS and Ridge results. The next code compares all three methods.\n", + "First we study and compare the OLS and Ridge results. The next code compares all three methods." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false, + "editable": true + }, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "\n", + "import os\n", + "import numpy as np\n", + "import pandas as pd\n", + "import matplotlib.pyplot as plt\n", + "\n", + "def R2(y_data, y_model):\n", + " return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2)\n", + "def MSE(y_data,y_model):\n", + " n = np.size(y_model)\n", + " return np.sum((y_data-y_model)**2)/n\n", "\n", "\n", + "# A seed just to ensure that the random numbers are the same for every run.\n", + "# Useful for eventual debugging.\n", + "\n", + "X = np.array( [ [ 2, 0], [0, 1], [0,0]])\n", + "y = np.array( [4, 2, 3])\n", + "\n", + "\n", + "# matrix inversion to find beta\n", + "OLSbeta = np.linalg.inv(X.T @ X) @ X.T @ y\n", + "print(OLSbeta)\n", + "# and then make the prediction\n", + "ytildeOLS = X @ OLSbeta\n", + "print(\"Training MSE for OLS\")\n", + "print(MSE(y,ytildeOLS))\n", + "ypredictOLS = X @ OLSbeta\n", + "\n", + "# Repeat now for Ridge regression and various values of the regularization parameter\n", + "I = np.eye(2,2)\n", + "# Decide which values of lambda to use\n", + "nlambdas = 100\n", + "MSEPredict = np.zeros(nlambdas)\n", + "lambdas = np.logspace(-4, 4, nlambdas)\n", + "for i in range(nlambdas):\n", + " lmb = lambdas[i]\n", + " Ridgebeta = np.linalg.inv(X.T @ X+lmb*I) @ X.T @ y\n", + "# print(Ridgebeta)\n", + " # and then make the prediction\n", + " ypredictRidge = X @ Ridgebeta\n", + " MSEPredict[i] = MSE(y,ypredictRidge)\n", + "# print(MSEPredict[i])\n", + " # Now plot the results\n", + "plt.figure()\n", + "plt.plot(np.log10(lambdas), MSEPredict, 'r--', label = 'MSE Ridge Test')\n", + "plt.xlabel('log10(lambda)')\n", + "plt.ylabel('MSE')\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We see here that we reach a plateau. What is actually happening?\n", + "\n", + "\n", + "## With Lasso Regression" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false, + "editable": true + }, + "outputs": [], + "source": [ + "import os\n", + "import numpy as np\n", + "import pandas as pd\n", + "import matplotlib.pyplot as plt\n", + "from sklearn import linear_model\n", + "\n", + "def R2(y_data, y_model):\n", + " return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2)\n", + "def MSE(y_data,y_model):\n", + " n = np.size(y_model)\n", + " return np.sum((y_data-y_model)**2)/n\n", + "\n", + "\n", + "# A seed just to ensure that the random numbers are the same for every run.\n", + "# Useful for eventual debugging.\n", + "\n", + "X = np.array( [ [ 2, 0], [0, 1], [0,0]])\n", + "y = np.array( [4, 2, 3])\n", + "\n", + "\n", + "# matrix inversion to find beta\n", + "OLSbeta = np.linalg.inv(X.T @ X) @ X.T @ y\n", + "print(OLSbeta)\n", + "# and then make the prediction\n", + "ytildeOLS = X @ OLSbeta\n", + "print(\"Training MSE for OLS\")\n", + "print(MSE(y,ytildeOLS))\n", + "ypredictOLS = X @ OLSbeta\n", + "\n", + "# Repeat now for Ridge regression and various values of the regularization parameter\n", + "I = np.eye(2,2)\n", + "# Decide which values of lambda to use\n", + "nlambdas = 100\n", + "MSERidgePredict = np.zeros(nlambdas)\n", + "MSELassoPredict = np.zeros(nlambdas)\n", + "lambdas = np.logspace(-4, 4, nlambdas)\n", + "for i in range(nlambdas):\n", + " lmb = lambdas[i]\n", + " Ridgebeta = np.linalg.inv(X.T @ X+lmb*I) @ X.T @ y\n", + "# print(Ridgebeta)\n", + " # and then make the prediction\n", + " ypredictRidge = X @ Ridgebeta\n", + " MSERidgePredict[i] = MSE(y,ypredictRidge)\n", + "# print(MSEPredict[i])\n", + " RegLasso = linear_model.Lasso(lmb)\n", + " RegLasso.fit(X,y)\n", + " ypredictLasso = RegLasso.predict(X)\n", + " MSELassoPredict[i] = MSE(y,ypredictLasso)\n", + "# Now plot the results\n", + "plt.figure()\n", + "plt.plot(np.log10(lambdas), MSERidgePredict, 'r--', label = 'MSE Ridge Train')\n", + "plt.plot(np.log10(lambdas), MSELassoPredict, 'r--', label = 'MSE Lasso Train')\n", + "plt.xlabel('log10(lambda)')\n", + "plt.ylabel('MSE')\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Another Example, now with a polynomial fit" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false, + "editable": true + }, + "outputs": [], + "source": [ + "import os\n", + "import numpy as np\n", + "import pandas as pd\n", + "import matplotlib.pyplot as plt\n", + "from sklearn.model_selection import train_test_split\n", + "from sklearn import linear_model\n", + "\n", + "def R2(y_data, y_model):\n", + " return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2)\n", + "def MSE(y_data,y_model):\n", + " n = np.size(y_model)\n", + " return np.sum((y_data-y_model)**2)/n\n", + "\n", + "\n", + "# A seed just to ensure that the random numbers are the same for every run.\n", + "# Useful for eventual debugging.\n", + "np.random.seed(3155)\n", + "\n", + "x = np.random.rand(100)\n", + "y = 2.0+5*x*x+0.1*np.random.randn(100)\n", + "\n", + "# number of features p (here degree of polynomial\n", + "p = 3\n", + "# The design matrix now as function of a given polynomial\n", + "X = np.zeros((len(x),p))\n", + "X[:,0] = 1.0\n", + "X[:,1] = x\n", + "X[:,2] = x*x\n", + "# We split the data in test and training data\n", + "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)\n", + "\n", + "# matrix inversion to find beta\n", + "OLSbeta = np.linalg.inv(X_train.T @ X_train) @ X_train.T @ y_train\n", + "print(OLSbeta)\n", + "# and then make the prediction\n", + "ytildeOLS = X_train @ OLSbeta\n", + "print(\"Training R2 for OLS\")\n", + "print(R2(y_train,ytildeOLS))\n", + "print(\"Training MSE for OLS\")\n", + "print(MSE(y_train,ytildeOLS))\n", + "ypredictOLS = X_test @ OLSbeta\n", + "print(\"Test R2 for OLS\")\n", + "print(R2(y_test,ypredictOLS))\n", + "print(\"Test MSE OLS\")\n", + "print(MSE(y_test,ypredictOLS))\n", + "\n", + "# Repeat now for Lasso and Ridge regression and various values of the regularization parameter\n", + "I = np.eye(p,p)\n", + "# Decide which values of lambda to use\n", + "nlambdas = 100\n", + "MSEPredict = np.zeros(nlambdas)\n", + "MSETrain = np.zeros(nlambdas)\n", + "MSELassoPredict = np.zeros(nlambdas)\n", + "MSELassoTrain = np.zeros(nlambdas)\n", + "lambdas = np.logspace(-4, 4, nlambdas)\n", + "for i in range(nlambdas):\n", + " lmb = lambdas[i]\n", + " Ridgebeta = np.linalg.inv(X_train.T @ X_train+lmb*I) @ X_train.T @ y_train\n", + " # include lasso using Scikit-Learn\n", + " RegLasso = linear_model.Lasso(lmb)\n", + " RegLasso.fit(X_train,y_train)\n", + " # and then make the prediction\n", + " ytildeRidge = X_train @ Ridgebeta\n", + " ypredictRidge = X_test @ Ridgebeta\n", + " ytildeLasso = RegLasso.predict(X_train)\n", + " ypredictLasso = RegLasso.predict(X_test)\n", + " MSEPredict[i] = MSE(y_test,ypredictRidge)\n", + " MSETrain[i] = MSE(y_train,ytildeRidge)\n", + " MSELassoPredict[i] = MSE(y_test,ypredictLasso)\n", + " MSELassoTrain[i] = MSE(y_train,ytildeLasso)\n", + "\n", + "# Now plot the results\n", + "plt.figure()\n", + "plt.plot(np.log10(lambdas), MSETrain, label = 'MSE Ridge train')\n", + "plt.plot(np.log10(lambdas), MSEPredict, 'r--', label = 'MSE Ridge Test')\n", + "plt.plot(np.log10(lambdas), MSELassoTrain, label = 'MSE Lasso train')\n", + "plt.plot(np.log10(lambdas), MSELassoPredict, 'r--', label = 'MSE Lasso Test')\n", + "\n", + "plt.xlabel('log10(lambda)')\n", + "plt.ylabel('MSE')\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ "\n", "## Linking the regression analysis with a statistical interpretation\n", "\n", @@ -1961,8 +2211,6 @@ }, "outputs": [], "source": [ - "%matplotlib inline\n", - "\n", "from numpy import *\n", "from numpy.random import randint, randn\n", "from time import time\n", diff --git a/doc/src/week36/week36.do.txt b/doc/src/week36/week36.do.txt index 8316ee482..201bbac0f 100644 --- a/doc/src/week36/week36.do.txt +++ b/doc/src/week36/week36.do.txt @@ -9,6 +9,11 @@ DATE: today * Thursday: Summary from last week on SVD, Statistics, probability theory and linear regression * Friday: Linear Regression and links with Statistics, Resampling methods and presentation of first project. +Recommended Reading: +o Lectures on Regression +o Bishop 1.1, 1.2, 2.1, 2.2, 2.3 and 3.1 +o Hastie et al chapter 3 + !split ===== Thursday September 9 ===== @@ -594,6 +599,218 @@ Here we set up the OLS, Ridge and Lasso functionality in order to study the abov First we study and compare the OLS and Ridge results. The next code compares all three methods. +!bc pycod +import os +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt + +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. + +X = np.array( [ [ 2, 0], [0, 1], [0,0]]) +y = np.array( [4, 2, 3]) + + +# matrix inversion to find beta +OLSbeta = np.linalg.inv(X.T @ X) @ X.T @ y +print(OLSbeta) +# and then make the prediction +ytildeOLS = X @ OLSbeta +print("Training MSE for OLS") +print(MSE(y,ytildeOLS)) +ypredictOLS = X @ OLSbeta + +# Repeat now for Ridge regression and various values of the regularization parameter +I = np.eye(2,2) +# Decide which values of lambda to use +nlambdas = 100 +MSEPredict = np.zeros(nlambdas) +lambdas = np.logspace(-4, 4, nlambdas) +for i in range(nlambdas): + lmb = lambdas[i] + Ridgebeta = np.linalg.inv(X.T @ X+lmb*I) @ X.T @ y +# print(Ridgebeta) + # and then make the prediction + ypredictRidge = X @ Ridgebeta + MSEPredict[i] = MSE(y,ypredictRidge) +# print(MSEPredict[i]) + # Now plot the results +plt.figure() +plt.plot(np.log10(lambdas), MSEPredict, 'r--', label = 'MSE Ridge Test') +plt.xlabel('log10(lambda)') +plt.ylabel('MSE') +plt.legend() +plt.show() + +!ec + +We see here that we reach a plateau. What is actually happening? + + +!split +===== With Lasso Regression ===== + +!bc pycod +import os +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +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. + +X = np.array( [ [ 2, 0], [0, 1], [0,0]]) +y = np.array( [4, 2, 3]) + + +# matrix inversion to find beta +OLSbeta = np.linalg.inv(X.T @ X) @ X.T @ y +print(OLSbeta) +# and then make the prediction +ytildeOLS = X @ OLSbeta +print("Training MSE for OLS") +print(MSE(y,ytildeOLS)) +ypredictOLS = X @ OLSbeta + +# Repeat now for Ridge regression and various values of the regularization parameter +I = np.eye(2,2) +# Decide which values of lambda to use +nlambdas = 100 +MSERidgePredict = np.zeros(nlambdas) +MSELassoPredict = np.zeros(nlambdas) +lambdas = np.logspace(-4, 4, nlambdas) +for i in range(nlambdas): + lmb = lambdas[i] + Ridgebeta = np.linalg.inv(X.T @ X+lmb*I) @ X.T @ y +# print(Ridgebeta) + # and then make the prediction + ypredictRidge = X @ Ridgebeta + MSERidgePredict[i] = MSE(y,ypredictRidge) +# print(MSEPredict[i]) + RegLasso = linear_model.Lasso(lmb) + RegLasso.fit(X,y) + ypredictLasso = RegLasso.predict(X) + MSELassoPredict[i] = MSE(y,ypredictLasso) +# Now plot the results +plt.figure() +plt.plot(np.log10(lambdas), MSERidgePredict, 'r--', label = 'MSE Ridge Train') +plt.plot(np.log10(lambdas), MSELassoPredict, 'r--', label = 'MSE Lasso Train') +plt.xlabel('log10(lambda)') +plt.ylabel('MSE') +plt.legend() +plt.show() + +!ec + +!split +===== Another Example, now with a polynomial fit ===== + +!bc pycod +import os +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) + +x = np.random.rand(100) +y = 2.0+5*x*x+0.1*np.random.randn(100) + +# number of features p (here degree of polynomial +p = 3 +# The design matrix now as function of a given polynomial +X = np.zeros((len(x),p)) +X[:,0] = 1.0 +X[:,1] = x +X[:,2] = x*x +# 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.inv(X_train.T @ X_train) @ X_train.T @ y_train +print(OLSbeta) +# and then make the prediction +ytildeOLS = X_train @ OLSbeta +print("Training R2 for OLS") +print(R2(y_train,ytildeOLS)) +print("Training MSE for OLS") +print(MSE(y_train,ytildeOLS)) +ypredictOLS = X_test @ OLSbeta +print("Test R2 for OLS") +print(R2(y_test,ypredictOLS)) +print("Test MSE OLS") +print(MSE(y_test,ypredictOLS)) + +# Repeat now for Lasso and Ridge regression and various values of the regularization parameter +I = np.eye(p,p) +# Decide which values of lambda to use +nlambdas = 100 +MSEPredict = np.zeros(nlambdas) +MSETrain = np.zeros(nlambdas) +MSELassoPredict = np.zeros(nlambdas) +MSELassoTrain = np.zeros(nlambdas) +lambdas = np.logspace(-4, 4, nlambdas) +for i in range(nlambdas): + lmb = lambdas[i] + Ridgebeta = np.linalg.inv(X_train.T @ X_train+lmb*I) @ X_train.T @ y_train + # include lasso using Scikit-Learn + RegLasso = linear_model.Lasso(lmb) + RegLasso.fit(X_train,y_train) + # and then make the prediction + ytildeRidge = X_train @ Ridgebeta + ypredictRidge = X_test @ Ridgebeta + ytildeLasso = RegLasso.predict(X_train) + ypredictLasso = RegLasso.predict(X_test) + MSEPredict[i] = MSE(y_test,ypredictRidge) + MSETrain[i] = MSE(y_train,ytildeRidge) + MSELassoPredict[i] = MSE(y_test,ypredictLasso) + MSELassoTrain[i] = MSE(y_train,ytildeLasso) + +# Now plot the results +plt.figure() +plt.plot(np.log10(lambdas), MSETrain, label = 'MSE Ridge train') +plt.plot(np.log10(lambdas), MSEPredict, 'r--', label = 'MSE Ridge Test') +plt.plot(np.log10(lambdas), MSELassoTrain, label = 'MSE Lasso train') +plt.plot(np.log10(lambdas), MSELassoPredict, 'r--', label = 'MSE Lasso Test') + +plt.xlabel('log10(lambda)') +plt.ylabel('MSE') +plt.legend() +plt.show() + + + + +!ec + + !split ===== Linking the regression analysis with a statistical interpretation =====