diff --git a/doc/src/week37/programs/franke.py b/doc/src/week37/programs/franke.py new file mode 100644 index 000000000..fc2924700 --- /dev/null +++ b/doc/src/week37/programs/franke.py @@ -0,0 +1,128 @@ +from mpl_toolkits.mplot3d import Axes3D +from matplotlib import cm +from matplotlib.ticker import LinearLocator, FormatStrFormatter +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 MSE(y_data,y_model): + n = np.size(y_model) + return np.sum((y_data-y_model)**2)/n + +def FrankeFunction(x,y): + term1 = 0.75*np.exp(-(0.25*(9*x-2)**2) - 0.25*((9*y-2)**2)) + term2 = 0.75*np.exp(-((9*x+1)**2)/49.0 - 0.1*(9*y+1)) + term3 = 0.5*np.exp(-(9*x-7)**2/4.0 - 0.25*((9*y-3)**2)) + term4 = -0.2*np.exp(-(9*x-4)**2 - (9*y-7)**2) + return term1 + term2 + term3 + term4 + + +def create_X(x, y, n ): + if len(x.shape) > 1: + x = np.ravel(x) + y = np.ravel(y) + + N = len(x) + l = int((n+1)*(n+2)/2) # Number of elements in beta + X = np.ones((N,l)) + + for i in range(1,n+1): + q = int((i)*(i+1)/2) + for k in range(i+1): + X[:,q+k] = (x**(i-k))*(y**k) + + return X + + +# Making meshgrid of datapoints and compute Franke's function +n = 10 +N = 1000 +x = np.sort(np.random.uniform(0, 1, N)) +y = np.sort(np.random.uniform(0, 1, N)) +z = FrankeFunction(x, y) +X = create_X(x, y, n=n) +# We split the data in test and training data +X_train, X_test, y_train, y_test = train_test_split(X, z, 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 = 100 +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, 'r', label = 'MSE own Ridge train') +plt.plot(np.log10(lambdas), MSEOwnRidgePredict, 'b--', label = 'MSE own Ridge Test') +plt.plot(np.log10(lambdas), MSERidgeTrain, 'y', label = 'MSE SL Ridge train') +plt.plot(np.log10(lambdas), MSERidgePredict, 'g--', label = 'MSE SL Ridge Test') + +plt.xlabel('log10(lambda)') +plt.ylabel('MSE') +plt.legend() +plt.show() + +# Now plot the function +fig = plt.figure() +ax=fig.add_subplot(projection='3d') + +# Make data. +x, y = np.meshgrid(x,y) +z = FrankeFunction(x, y) +znew = X @ OLSbeta +# Plot the surface. +surf = ax.plot_surface(x, y, z, cmap=cm.coolwarm, + linewidth=0, antialiased=False) + + +# Customize the z axis. +ax.set_zlim(-0.10, 1.40) +ax.zaxis.set_major_locator(LinearLocator(10)) +ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f')) + +# Add a color bar which maps values to colors. +fig.colorbar(surf, shrink=0.5, aspect=5) + +plt.show() + + + + diff --git a/doc/src/week37/programs/noscaling.py b/doc/src/week37/programs/noscaling.py deleted file mode 100644 index 7eb9769e9..000000000 --- a/doc/src/week37/programs/noscaling.py +++ /dev/null @@ -1,72 +0,0 @@ -import matplotlib.pyplot as plt -import numpy as np -from sklearn.linear_model import LinearRegression -from sklearn.preprocessing import PolynomialFeatures -from sklearn.model_selection import train_test_split -from sklearn.preprocessing import StandardScaler - -def MSE(y_data,y_model): - n = np.size(y_model) - return np.sum((y_data-y_model)**2)/n - -def OLS_fit_beta(X, y): - return np.linalg.pinv(X.T @ X) @ X.T @ y - -def Ridge_fit_beta(X, y,L,d): - I = np.eye(d,d) - return np.linalg.pinv(X.T @ X + L*I) @ X.T @ y - - -np.random.seed(2018) -n = 100 -d = 3 -L = 0.001 -true_beta = [2, 0.5, 3.7] - -# Make data set. -x = np.linspace(-3, 3, n) -y_real = 2 + 0.5*x + 3.7*x**2 - -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)) - - -#Design matrix X including the intercept -X = np.zeros((len(x), d)) -for p in range(d): # (d-1) - X[:, p] = x ** (p) # (p+1 if not intercept included) - - -#Split datamatrix -X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2) - - -#Calculate beta, own code -beta_OLS = OLS_fit_beta(X_train, y_train) -beta_Ridge = Ridge_fit_beta(X_train, y_train,L,d) -print(beta_OLS) -print(beta_Ridge) - -#predict value -ytilde_test_OLS = X_test @ beta_OLS -ytilde_test_Ridge = X_test @ beta_Ridge - - -#Calculate MSE - -print(" ") -print("test MSE of OLS:") -print(MSE(y_test,ytilde_test_OLS)) -print(" ") -print("test MSE of Ridge") -print(MSE(y_test,ytilde_test_Ridge)) - - -plt.scatter(x,y,label='Data') -#plt.plot(x,y_real,label='no noise') -plt.plot(x, X @ beta_OLS,'*', label="OLS_Fit") -plt.plot(x, X @ beta_Ridge, label="Ridge_Fit") -plt.grid() -plt.legend() -plt.show() diff --git a/doc/src/week37/programs/plotfranke.py b/doc/src/week37/programs/plotfranke.py new file mode 100644 index 000000000..644d68208 --- /dev/null +++ b/doc/src/week37/programs/plotfranke.py @@ -0,0 +1,39 @@ +from mpl_toolkits.mplot3d import Axes3D +import matplotlib.pyplot as plt +from matplotlib import cm +from matplotlib.ticker import LinearLocator, FormatStrFormatter +import numpy as np +from random import random, seed + +fig = plt.figure() +ax=fig.add_subplot(projection='3d') + +# Make data. +x = np.arange(0, 1, 0.05) +y = np.arange(0, 1, 0.05) +x, y = np.meshgrid(x,y) + + +def FrankeFunction(x,y): + term1 = 0.75*np.exp(-(0.25*(9*x-2)**2) - 0.25*((9*y-2)**2)) + term2 = 0.75*np.exp(-((9*x+1)**2)/49.0 - 0.1*(9*y+1)) + term3 = 0.5*np.exp(-(9*x-7)**2/4.0 - 0.25*((9*y-3)**2)) + term4 = -0.2*np.exp(-(9*x-4)**2 - (9*y-7)**2) + return term1 + term2 + term3 + term4 + + +z = FrankeFunction(x, y) + +# Plot the surface. +surf = ax.plot_surface(x, y, z, cmap=cm.coolwarm, + linewidth=0, antialiased=False) + +# Customize the z axis. +ax.set_zlim(-0.10, 1.40) +ax.zaxis.set_major_locator(LinearLocator(10)) +ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f')) + +# Add a color bar which maps values to colors. +fig.colorbar(surf, shrink=0.5, aspect=5) + +plt.show() diff --git a/doc/src/week37/programs/sklearnscaling.py b/doc/src/week37/programs/sklearnscaling.py deleted file mode 100644 index 44c51e37b..000000000 --- a/doc/src/week37/programs/sklearnscaling.py +++ /dev/null @@ -1,88 +0,0 @@ -import matplotlib.pyplot as plt -import numpy as np -from sklearn.linear_model import LinearRegression -from sklearn.preprocessing import PolynomialFeatures -from sklearn.model_selection import train_test_split -from sklearn.preprocessing import StandardScaler - -def MSE(y_data,y_model): - n = np.size(y_model) - return np.sum((y_data-y_model)**2)/n - -def OLS_fit_beta(X, y): - return np.linalg.pinv(X.T @ X) @ X.T @ y - -def Ridge_fit_beta(X, y,L,d): - I = np.eye(d,d) - return np.linalg.pinv(X.T @ X + L*I) @ X.T @ y - - -np.random.seed(2018) -n = 1000 -d = 3 -L = 0.001 -true_beta = [2, 0.5, 3.7] - -# Make data set. -x = np.linspace(-3, 3, n) -y_real = 2 + 0.5*x + 3.7*x**2 - -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)) - - -#Design matrix X does include the intercept -X = np.zeros((len(x), d)) -for p in range(d): # (d-1) - X[:, p] = x ** (p+1) # (p+1 if not intercept included) - - -#Split datamatrix -X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2) - -scaler = StandardScaler() -yscaler = StandardScaler() -scaler.fit(X_train) -yscaler.fit(y_train) -X_train_scaled = scaler.transform(X_train) -X_test_scaled = scaler.transform(X_test) -y_train_scaled = yscaler.transform(y_train) -y_test_scaled = yscaler.transform(y_test) - - - -#Calculate beta -beta_OLS = OLS_fit_beta(X_train_scaled, y_train_scaled) -beta_Ridge = Ridge_fit_beta(X_train_scaled, y_train_scaled,L,d) -print(beta_OLS) -print(beta_Ridge) - -""" -interceptOLS = y_scaler - X_train_mean @ beta_OLS -interceptRidge = y_scaler - X_train_mean @ beta_Ridge -print(interceptOLS) -print(interceptRidge) -""" -#predict value -ytilde_test_OLS = X_test_scaled @ beta_OLS -ytilde_test_Ridge = X_test_scaled @ beta_Ridge - - -#Calculate MSE - -print(" ") -print("test MSE of OLS:") -print(MSE(y_test,ytilde_test_OLS)) -print(" ") -print("test MSE of Ridge") -print(MSE(y_test,ytilde_test_Ridge)) - - -plt.scatter(x,y,label='Data') -#plt.plot(x,y_real,label='no noise') -plt.plot(x, X @ beta_OLS,'*', label="OLS_Fit") -plt.plot(x, X @ beta_Ridge, label="Ridge_Fit") -plt.grid() -plt.legend() -plt.show()