update
This commit is contained in:
Morten Hjorth-Jensen
File diff suppressed because one or more lines are too long
+22
-135
@@ -7,9 +7,22 @@ DATE: today
|
||||
|
||||
===== This note contains code examples with a simple scaling =====
|
||||
|
||||
The programs here use both ordinrary least squares and Ridge regression with one value only for
|
||||
the hyperparameter $\lambda$. The first example has no scaling and includes the intercept as well and we are trying to fit a second-order
|
||||
polynomial.
|
||||
The programs here use both ordinrary least squares (OLS) and Ridge
|
||||
regression with one value only for the hyperparameter $\lambda$. The
|
||||
first example has no scaling and includes the intercept as well and we
|
||||
are trying to fit a second-order polynomial. The second code takes out
|
||||
the intercept and subtracts the mean values of each column of the
|
||||
design matrix and the mean value of the outputs.
|
||||
|
||||
The third and final code uses _Scikit-Learn_ as library in order to
|
||||
calculate the optimal parameters for OLS and Ridge regression. Note
|
||||
that it is highly recommended to not include the intercept in Ridge
|
||||
and Lasso regression, in order to avoid penalizing the optimization by
|
||||
the intercept. The second and third codes do thus not include the
|
||||
intercept. In the second code we do the scaling ourselves while the
|
||||
last code uses the standard scaler option included in _Scikit-Learn_, known as centering (where
|
||||
we subtract the mean values).
|
||||
|
||||
|
||||
!bc pycod
|
||||
import matplotlib.pyplot as plt
|
||||
@@ -34,14 +47,14 @@ def Ridge_fit_beta(X, y,L,d):
|
||||
np.random.seed(2018)
|
||||
n = 100
|
||||
d = 3
|
||||
# hyperparameter lambda
|
||||
Lambda = 0.01
|
||||
true_beta = [2, 0.5, 3.7]
|
||||
|
||||
# Make data set.
|
||||
# Make data set, simple second-order polynomial
|
||||
x = np.linspace(-3, 3, n)
|
||||
y = 2 + 0.5*x + 3.7*x**2
|
||||
y = 2.0 + 0.5*x + 5.0*(x**2)+ np.random.randn(n)
|
||||
|
||||
#Design matrix X includes the intercept and scaling is made
|
||||
# The design matrix X includes the intercept and no scaling is made
|
||||
X = np.zeros((len(x), d))
|
||||
for p in range(d):
|
||||
X[:, p] = x ** (p)
|
||||
@@ -74,17 +87,17 @@ plt.plot(x, X @ beta_Ridge, label="Ridge_Fit")
|
||||
plt.grid()
|
||||
plt.legend()
|
||||
plt.show()
|
||||
|
||||
!ec
|
||||
|
||||
|
||||
In this example we do not include the intercept and we scale the data by subtracting the mean values. This follows the discussion in the "lecture material":"https://compphysics.github.io/MachineLearning/doc/LectureNotes/_build/html/chapter3.html#more-on-rescaling-data".
|
||||
see also the weekly slides "for week 36":"https://compphysics.github.io/MachineLearning/doc/pub/week36/html/._week36-bs029.html".
|
||||
It is recommended whrn we use Ridge and Lasso regression to not include the intercept in the optimization process.
|
||||
|
||||
Before we discuss the code, we repeat some of the basic math from the slides of week 36.
|
||||
|
||||
Let us try to understand what this may imply mathematically when we
|
||||
subtract the mean values, also known as *zero centering*. For
|
||||
subtract the mean values, also known as *zero centering* or simply *centering*. For
|
||||
simplicity, we will focus on ordinary regression, as done in the above example.
|
||||
|
||||
The cost/loss function for regression is
|
||||
@@ -208,137 +221,14 @@ Now we try to implement this.
|
||||
|
||||
!bc pycod
|
||||
|
||||
np.random.seed(2018)
|
||||
n = 100
|
||||
# we do not include the intercept
|
||||
d = 2
|
||||
Lambda = 0.01
|
||||
true_beta = [2, 0.5, 3.7]
|
||||
|
||||
# Make data set.
|
||||
x = np.linspace(-3, 3, n)
|
||||
y = 2 + 0.5*x + 3.7*x**2
|
||||
|
||||
#Design matrix X does not include the intercept.
|
||||
X = np.zeros((len(x), d))
|
||||
for p in range(d):
|
||||
X[:, p] = x ** (p+1)
|
||||
|
||||
|
||||
#Split data in train and test
|
||||
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)
|
||||
|
||||
# Scale data by subtracting mean value,own implementation
|
||||
#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)
|
||||
#Center by removing mean from each feature
|
||||
X_train_scaled = X_train - X_train_mean
|
||||
X_test_scaled = X_test - X_train_mean
|
||||
#The model intercept (called y_scaler) is given by the mean of the target variable (IF X is centered, note)
|
||||
y_scaler = np.mean(y_train)
|
||||
y_train_scaled = y_train - y_scaler
|
||||
|
||||
|
||||
#Calculate beta
|
||||
beta_OLS = OLS_fit_beta(X_train_scaled, y_train_scaled)
|
||||
beta_Ridge = Ridge_fit_beta(X_train_scaled, y_train_scaled,Lambda,d)
|
||||
print(beta_OLS)
|
||||
print(beta_Ridge)
|
||||
# calculate intercepts and print them
|
||||
interceptOLS = y_scaler - X_train_mean @ beta_OLS
|
||||
interceptRidge = y_scaler - X_train_mean @ beta_Ridge
|
||||
print(interceptOLS)
|
||||
print(interceptRidge)
|
||||
|
||||
#predict value with intercept
|
||||
ytilde_test_OLS = X_test_scaled @ beta_OLS+y_scaler
|
||||
ytilde_test_Ridge = X_test_scaled @ beta_Ridge+y_scaler
|
||||
|
||||
|
||||
#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, X @ beta_OLS+interceptOLS,'*', label="OLS_Fit")
|
||||
plt.plot(x, X @ beta_Ridge+interceptRidge, label="Ridge_Fit")
|
||||
plt.grid()
|
||||
plt.legend()
|
||||
plt.show()
|
||||
|
||||
!ec
|
||||
|
||||
We see that we get the same values for the parameters! As it should be. The MSE may however change (not the case here).
|
||||
|
||||
Finally, instead of using our own function we repeat the same example
|
||||
using the _standardscaler_ functionality of the library
|
||||
_Scikit-Learn_. Here we limit ourselves to Ridge regression only.
|
||||
|
||||
!bc pycod
|
||||
|
||||
from sklearn import linear_model
|
||||
np.random.seed(2018)
|
||||
n = 100
|
||||
d = 2
|
||||
Lambda = 0.01
|
||||
true_beta = [2, 0.5, 3.7]
|
||||
|
||||
# Make data set.
|
||||
x = np.linspace(-3, 3, n)
|
||||
y = (2 + 0.5*x + 3.7*x**2)
|
||||
|
||||
#Design matrix X does not include the intercept.
|
||||
X = np.zeros((n, d))
|
||||
for p in range(d):
|
||||
X[:, p] = x ** (p+1)
|
||||
|
||||
#Split data in train and test
|
||||
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)
|
||||
|
||||
|
||||
# Scale data by subtracting mean value using scikit-learn
|
||||
from sklearn.preprocessing import StandardScaler
|
||||
scaler = StandardScaler()
|
||||
#scaler.fit(X_train)
|
||||
#scaler.fit(y_train)
|
||||
#X_train_scaled = scaler.transform(X_train)
|
||||
#X_test_scaled = scaler.transform(X_test)
|
||||
#y_train_scaled = scaler.transform(y_train)
|
||||
|
||||
#Calculate beta
|
||||
OLS = LinearRegression()
|
||||
OLS.fit(X_train,y_train)
|
||||
ypredictOLS = OLS.predict(X_test)
|
||||
RegRidge = linear_model.Ridge(Lambda)
|
||||
RegRidge.fit(X_train,y_train)
|
||||
ypredictRidge = RegRidge.predict(X_test)
|
||||
print(OLS.coef_)
|
||||
print(RegRidge.coef_)
|
||||
print(OLS.intercept_)
|
||||
interceptRidge = RegRidge.intercept_
|
||||
print(RegRidge.intercept_)
|
||||
#predict value without intercept
|
||||
ytilde_test_Ridge = X_test @ RegRidge.coef_+ RegRidge.intercept_
|
||||
ytilde_test_OLS = X_test @ OLS.coef_+ OLS.intercept_
|
||||
|
||||
#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, X @ RegRidge.coef_ + RegRidge.intercept_ , label="Ridge_Fit")
|
||||
plt.grid()
|
||||
plt.legend()
|
||||
plt.show()
|
||||
!ec
|
||||
|
||||
|
||||
@@ -346,6 +236,3 @@ plt.show()
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
Reference in New Issue
Block a user