Morten Hjorth-Jensen
101 lines
3.1 KiB
Python
101 lines
3.1 KiB
Python
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 = 4
|
|
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)
|
|
intercept_ = y_scaler - X_train_mean@OwnRidgeBeta #The intercept can be shifted so the model can predict on uncentered data
|
|
|
|
ypredictOwnRidge = X_test @ OwnRidgeBeta + intercept_ #Add intercept to prediction
|
|
#EQUIVALENT PREDICTION:
|
|
ypredictOwnRidge = X_test_scaled @ OwnRidgeBeta + y_scaler #Add intercept 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(intercept_)
|
|
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()
|
|
|