diff --git a/doc/src/week36/scale2.py b/doc/src/week36/scale2.py new file mode 100644 index 000000000..8f1004bfd --- /dev/null +++ b/doc/src/week36/scale2.py @@ -0,0 +1,94 @@ +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +from sklearn.model_selection import train_test_split +from sklearn import linear_model +from sklearn.preprocessing import StandardScaler + +def R2(y_data, y_model): + return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2) +def MSE(y_data,y_model): + n = np.size(y_model) + return np.sum((y_data-y_model)**2)/n + + +# A seed just to ensure that the random numbers are the same for every run. +# Useful for eventual debugging. +np.random.seed(315) + +n = 100 +x = np.random.rand(n) +y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2) + +Maxpolydegree = 5 +X = np.zeros((n,Maxpolydegree-1)) + +for degree in range(1,Maxpolydegree): #No intercept column + X[:,degree-1] = x**(degree) + + + + +# We split the data in test and training data +X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2) + + + + + +#For our own implementation, we will need to deal with the intercept by centering the design matrix and the target variable +X_train_mean = np.mean(X_train,axis=0) +X_train_scaled = X_train - X_train_mean #Center by removing mean from each feature +X_test_scaled = X_test - X_train_mean + +y_scaler = np.mean(y_train) #The model intercept (called y_scaler) is given by the mean of target variable (IF X is centered) +y_train_scaled = y_train - y_scaler #Remove the intercept from the training data. + + +p = Maxpolydegree-1 +I = np.eye(p,p) +# Decide which values of lambda to use +nlambdas = 1 +MSEOwnRidgePredict = np.zeros(nlambdas) +MSERidgePredict = np.zeros(nlambdas) + +lambdas = np.logspace(-4, 1, nlambdas) +for i in range(nlambdas): + lmb = lambdas[i] + OwnRidgeBeta = np.linalg.pinv(X_train_scaled.T @ X_train_scaled+lmb*I) @ X_train_scaled.T @ (y_train_scaled) + ypredictOwnRidge = X_test_scaled @ OwnRidgeBeta + y_scaler #Add intercept (y_scaler) to prediction + print("Values for own Ridge prediction") + print(ypredictOwnRidge) + + + + RegRidge = linear_model.Ridge(lmb) + RegRidge.fit(X_train,y_train) + ypredictRidge = RegRidge.predict(X_test) + print("Values for SL Ridge prediction") + print(ypredictRidge) + + + MSEOwnRidgePredict[i] = MSE(y_test,ypredictOwnRidge) + MSERidgePredict[i] = MSE(y_test,ypredictRidge) + + print("Beta values for own Ridge implementation") + print(OwnRidgeBeta) #Intercept is given by mean of target variable + print("Beta values for Scikit-Learn Ridge implementation") + print(RegRidge.coef_) + print('Intercept from own implementation:') + print(y_scaler) + print('Intercept from Scikit-Learn Ridge implementation') + print(RegRidge.intercept_) + +# Now plot the results + +plt.figure() +plt.plot(np.log10(lambdas), MSEOwnRidgePredict, 'b--', label = 'MSE own Ridge Test') +plt.plot(np.log10(lambdas), MSERidgePredict, 'g--', label = 'MSE SL Ridge Test') + +plt.xlabel('log10(lambda)') +plt.ylabel('MSE') +plt.legend() +plt.show() + diff --git a/doc/src/week36/test2.py b/doc/src/week36/test2.py new file mode 100644 index 000000000..b101077b5 --- /dev/null +++ b/doc/src/week36/test2.py @@ -0,0 +1,70 @@ +import numpy as np +import matplotlib.pyplot as plt + +from sklearn.linear_model import LinearRegression + + +np.random.seed(2021) + + +def fit_beta(X, y): + return np.linalg.pinv(X.T @ X) @ X.T @ y + + +true_beta = [2, 0.5, 3.7] + +x = np.linspace(0, 1, 11) +y = np.sum( + np.asarray([x ** p * b for p, b in enumerate(true_beta)]), axis=0 +) + 0.1 * np.random.normal(size=len(x)) + +degree = 3 +X = np.zeros((len(x), degree)) + +# Include the intercept in the design matrix +for p in range(degree): + X[:, p] = x ** p + +beta = fit_beta(X, y) + +# Intercept is included in the design matrix +clf = LinearRegression(fit_intercept=False).fit(X, y) + +print(f"True beta: {true_beta}") +print(f"Fitted beta: {beta}") +print(f"Sklearn fitted beta: {clf.coef_}") + + +plt.figure() +plt.scatter(x, y, label="Data") +plt.plot(x, X @ beta, label="Fit") +plt.plot(x, clf.predict(X), label="Sklearn (fit_intercept=False)") + + +# Do not include the intercept in the design matrix +X = np.zeros((len(x), degree - 1)) + +for p in range(degree - 1): + X[:, p] = x ** (p + 1) + +# Intercept is not included in the design matrix +clf = LinearRegression(fit_intercept=True).fit(X, y) + +# Use centered values for X and y when computing coefficients +y_offset = np.average(y, axis=0) +X_offset = np.average(X, axis=0) + +beta = fit_beta(X - X_offset, y - y_offset) +intercept = np.mean(y_offset - X_offset @ beta) + +print(f"Manual intercept: {intercept}") +print(f"Fitted beta (sans intercept): {beta}") +print(f"Sklearn intercept: {clf.intercept_}") +print(f"Sklearn fitted beta (sans intercept): {clf.coef_}") + +plt.plot(x, X @ beta + intercept, "--", label="Fit (manual intercept)") +plt.plot(x, clf.predict(X), "--", label="Sklearn (fit_intercept=True)") +plt.grid() +plt.legend() + +plt.show()