added file

doc/src/week36/franke.py

@@ -0,0 +1,106 @@
+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()
+
+
+
+
+
+
+

doc/src/week36/test.py

@@ -74,10 +74,10 @@ for i in range(nlambdas):
     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.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')