Morten Hjorth-Jensen
107 lines
3.0 KiB
Python
107 lines
3.0 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
|
|
|
|
|
|
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 = 5
|
|
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 = 5
|
|
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()
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|