update week 38

2021-09-20 18:03:32 +02:00
parent 1751d37f48
commit 03ff1cd657
47 changed files with 2908 additions and 2202 deletions
@@ -0,0 +1,86 @@
+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, label = 'MSE Ridge train')
+plt.plot(np.log10(lambdas), MSEOwnRidgePredict, 'r--', label = 'MSE Ridge Test')
+plt.plot(np.log10(lambdas), MSERidgeTrain, label = 'MSE Ridge train')
+plt.plot(np.log10(lambdas), MSERidgePredict, 'r--', label = 'MSE Ridge Test')
+
+plt.xlabel('log10(lambda)')
+plt.ylabel('MSE')
+plt.legend()
+plt.show()
+
@@ -233,6 +233,100 @@ plt.show()
 !ec


+!split
+===== To think about =====
+
+When you are comparing your own code with for example _Scikit-Learn_'s library, there are some minor things to keep in mind.
+The example here shows how one can leave out or keep the intercept.
+
+!bc pycod
+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, label = 'MSE Ridge train')
+plt.plot(np.log10(lambdas), MSEOwnRidgePredict, 'r--', label = 'MSE Ridge Test')
+plt.plot(np.log10(lambdas), MSERidgeTrain, label = 'MSE Ridge train')
+plt.plot(np.log10(lambdas), MSERidgePredict, 'r--', label = 'MSE Ridge Test')
+
+plt.xlabel('log10(lambda)')
+plt.ylabel('MSE')
+plt.legend()
+plt.show()
+
+!ec

 !split
 =====  More complicated Example: The Ising model =====