From f4d8b8ceca945576ea359178c4075048d08b3c48 Mon Sep 17 00:00:00 2001 From: Morten Hjorth-Jensen Date: Thu, 9 Sep 2021 09:24:06 +0200 Subject: [PATCH] updating lectures --- doc/pub/week36/html/week36-bs.html | 89 ++++--- doc/pub/week36/html/week36-reveal.html | 223 ++++++++++++++++ doc/pub/week36/html/week36-solarized.html | 225 ++++++++++++++++ doc/pub/week36/html/week36.html | 225 ++++++++++++++++ doc/pub/week36/ipynb/ipynb-week36-src.tar.gz | Bin 191 -> 192 bytes doc/pub/week36/ipynb/week36.ipynb | 254 ++++++++++++++++++- doc/src/week36/week36.do.txt | 217 ++++++++++++++++ 7 files changed, 1189 insertions(+), 44 deletions(-) diff --git a/doc/pub/week36/html/week36-bs.html b/doc/pub/week36/html/week36-bs.html index a1cd812e4..9def9c7a3 100644 --- a/doc/pub/week36/html/week36-bs.html +++ b/doc/pub/week36/html/week36-bs.html @@ -98,6 +98,11 @@ Automatically generated HTML file from DocOnce source 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, @@ -281,46 +286,48 @@ MathJax.Hub.Config({
  • Lasso case
  • The first Case
  • Simple code for solving the above problem
    • -
  • Linking the regression analysis with a statistical interpretation
    • -
  • Assumptions made
    • -
  • Expectation value and variance
    • -
  • Expectation value and variance for \( \boldsymbol{\beta} \)
    • -
  • Deriving OLS from a probability distribution
    • -
  • Independent and Identically Distrubuted (iid)
    • -
  • Maximum Likelihood Estimation (MLE)
    • -
  • A new Cost Function
    • -
  • More basic Statistics and Bayes' theorem
    • -
  • Marginal Probability
    • -
  • Conditional Probability
    • -
  • Bayes' Theorem
    • -
  • Frliday September 10
    • -
  • Why resampling methods
    • -
  • Resampling methods
    • -
  • Resampling approaches can be computationally expensive
    • -
  • Why resampling methods ?
    • -
  • Statistical analysis
    • -
  • Resampling methods
    • -
  • Resampling methods: Jackknife and Bootstrap
    • -
  • Resampling methods: Jackknife
    • -
  • Jackknife code example
    • -
  • Resampling methods: Bootstrap
    • -
  • Resampling methods: Bootstrap background
    • -
  • Resampling methods: More Bootstrap background
    • -
  • Resampling methods: Bootstrap approach
    • -
  • Resampling methods: Bootstrap steps
    • -
  • Code example for the Bootstrap method
    • -
  • Various steps in cross-validation
    • -
  • How to set up the cross-validation for Ridge and/or Lasso
    • -
  • Cross-validation in brief
    • -
  • Code Example for Cross-validation and \( k \)-fold Cross-validation
    • -
  • The bias-variance tradeoff
    • -
  • Example code for Bias-Variance tradeoff
    • -
  • Understanding what happens
    • -
  • Summing up
    • -
  • Another Example from Scikit-Learn's Repository
    • -
  • More examples on bootstrap and cross-validation and errors
    • -
  • The same example but now with cross-validation
    • -
  • Cross-validation with Ridge
    • +
  • With Lasso Regression
    • +
  • Another Example, now with a polynomial fit
    • +
  • Linking the regression analysis with a statistical interpretation
    • +
  • Assumptions made
    • +
  • Expectation value and variance
    • +
  • Expectation value and variance for \( \boldsymbol{\beta} \)
    • +
  • Deriving OLS from a probability distribution
    • +
  • Independent and Identically Distrubuted (iid)
    • +
  • Maximum Likelihood Estimation (MLE)
    • +
  • A new Cost Function
    • +
  • More basic Statistics and Bayes' theorem
    • +
  • Marginal Probability
    • +
  • Conditional Probability
    • +
  • Bayes' Theorem
    • +
  • Frliday September 10
    • +
  • Why resampling methods
    • +
  • Resampling methods
    • +
  • Resampling approaches can be computationally expensive
    • +
  • Why resampling methods ?
    • +
  • Statistical analysis
    • +
  • Resampling methods
    • +
  • Resampling methods: Jackknife and Bootstrap
    • +
  • Resampling methods: Jackknife
    • +
  • Jackknife code example
    • +
  • Resampling methods: Bootstrap
    • +
  • Resampling methods: Bootstrap background
    • +
  • Resampling methods: More Bootstrap background
    • +
  • Resampling methods: Bootstrap approach
    • +
  • Resampling methods: Bootstrap steps
    • +
  • Code example for the Bootstrap method
    • +
  • Various steps in cross-validation
    • +
  • How to set up the cross-validation for Ridge and/or Lasso
    • +
  • Cross-validation in brief
    • +
  • Code Example for Cross-validation and \( k \)-fold Cross-validation
    • +
  • The bias-variance tradeoff
    • +
  • Example code for Bias-Variance tradeoff
    • +
  • Understanding what happens
    • +
  • Summing up
    • +
  • Another Example from Scikit-Learn's Repository
    • +
  • More examples on bootstrap and cross-validation and errors
    • +
  • The same example but now with cross-validation
    • +
  • Cross-validation with Ridge
    • @@ -379,7 +386,7 @@ MathJax.Hub.Config({
  • 9
  • 10
  • ...
    • -
  • 64
    • +
  • 66
  • »
    • diff --git a/doc/pub/week36/html/week36-reveal.html b/doc/pub/week36/html/week36-reveal.html index 343fc90ab..2eff9beef 100644 --- a/doc/pub/week36/html/week36-reveal.html +++ b/doc/pub/week36/html/week36-reveal.html @@ -165,6 +165,15 @@ MathJax.Hub.Config({

  • 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: + +

      +

    1. Lectures on Regression
      2. +

    2. Bishop 1.1, 1.2, 2.1, 2.2, 2.3 and 3.1
      3. +

    3. Hastie et al chapter 3
      4. +
    @@ -864,6 +873,220 @@ 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. + +

    + + +

    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()
+
    diff --git a/doc/pub/week36/html/week36-solarized.html b/doc/pub/week36/html/week36-solarized.html index f7eca64eb..fc95d417d 100644 --- a/doc/pub/week36/html/week36-solarized.html +++ b/doc/pub/week36/html/week36-solarized.html @@ -118,6 +118,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, @@ -294,6 +299,14 @@ MathJax.Hub.Config({
  • Friday: Linear Regression and links with Statistics, Resampling methods and presentation of first project.
    • +Recommended Reading: + +
      +
    1. Lectures on Regression
      2. +
    2. Bishop 1.1, 1.2, 2.1, 2.2, 2.3 and 3.1
      3. +
    3. Hastie et al chapter 3
      4. +
    +









    Thursday September 9

    @@ -879,6 +892,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. +

    + + +

    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()
+

    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({

  • Friday: Linear Regression and links with Statistics, Resampling methods and presentation of first project.
    • +Recommended Reading: + +
      +
    1. Lectures on Regression
      2. +
    2. Bishop 1.1, 1.2, 2.1, 2.2, 2.3 and 3.1
      3. +
    3. Hastie et al chapter 3
      4. +
    +









    Thursday September 9

    @@ -884,6 +897,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. +

    + + +

    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()
+

    diff --git a/doc/pub/week36/ipynb/ipynb-week36-src.tar.gz b/doc/pub/week36/ipynb/ipynb-week36-src.tar.gz index 264bd23e0fb4a71a090aa16a65fddd923556d1d7..c07b5e271a7db591519421c96b786d43f0383c4e 100644 GIT binary patch literal 192 zcmV;x06+g9iwFRVwmD${1MSbv3c@f92k@Qu6nTQtu3J3|?%+WX@dY}TxjNU*wnO*! z?gR9sco`z}cli?%LUPE~n_U*Uy9*XW2;oTrW3G5QC5h*HLTL&#p_DVmV#_Wl|}V#ZWt@e4}0cU;F*8oSV;@peeWu*Kxv1$)-~J^YsMql_9}-$ uqa9md@Y+d(AXE>cD5R5GiA&fTeKMl4QTXd;JkRsIuRQ<-`4^@D2mkuHia9|^nulw_E?fvAULd8ZjkQTlQna_X z573q3rihSl^E1pa%p9`yW|t-Y?xV#ZggB)z=8Cf^Q99QXj2U2xF()adi~+)g5flMh z?xdH_d0~fBn(8dIlltAwd?K@bGN*B;VgYXtxZ006s5SIYnZ 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 =====