-
@@ -326,7 +321,7 @@ MathJax.Hub.Config({
+Note that we kept the intercept column in the fitting here. This means that we need to set the intercept in the call to the Scikit-Learn function as False. Alternatively, we could have set up the design matrix \( X \) without the first column of ones. +
@@ -384,7 +382,6 @@ plt.show()
+In this example we keep the intercept column again but add cross-validation in order to estimate the best possible value of the means squared error.
@@ -344,7 +341,7 @@ kfold = KFold(n_splits = polydegree for degree in range(polydegree): X[:,degree] = Density**(degree/3.0) - OLS = LinearRegression() + OLS = LinearRegression(fit_intercept=False) # loop over trials in order to estimate the expectation value of the MSE estimated_mse_folds = cross_val_score(OLS, X, Energies, scoring='neg_mean_squared_error', cv=kfold) #[:, np.newaxis] @@ -357,6 +354,7 @@ plt.legend() plt.show()
+
diff --git a/doc/pub/week37/html/week37-bs.html b/doc/pub/week37/html/week37-bs.html index f984145ed..2ae1b501d 100644 --- a/doc/pub/week37/html/week37-bs.html +++ b/doc/pub/week37/html/week37-bs.html @@ -177,11 +177,7 @@ Automatically generated HTML file from DocOnce source ('The same example but now with cross-validation', 2, None, - 'the-same-example-but-now-with-cross-validation'), - ('Cross-validation with Ridge', - 2, - None, - 'cross-validation-with-ridge')]} + 'the-same-example-but-now-with-cross-validation')]} end of tocinfo --> @@ -267,7 +263,6 @@ MathJax.Hub.Config({
-
@@ -326,7 +321,7 @@ MathJax.Hub.Config({
-
@@ -1845,7 +1845,7 @@ trials = 100 trainingerror[polydegree] = 0.0 for samples in range(trials): x_train, x_test, y_train, y_test = train_test_split(X, Energies, test_size=0.2) - model = LinearRegression(fit_intercept=True).fit(x_train, y_train) + model = LinearRegression(fit_intercept=False).fit(x_train, y_train) ypred = model.predict(x_train) ytilde = model.predict(x_test) testerror[polydegree] += mean_squared_error(y_test, ytilde) @@ -1864,12 +1864,16 @@ plt.ylabel('log10[MSE]') plt.legend() plt.show() +
+Note that we kept the intercept column in the fitting here. This means that we need to set the intercept in the call to the Scikit-Learn function as False. Alternatively, we could have set up the design matrix \( X \) without the first column of ones.
+In this example we keep the intercept column again but add cross-validation in order to estimate the best possible value of the means squared error.
@@ -1929,7 +1933,7 @@ kfold = KFold(n_splits = k)
polynomial[polydegree] = polydegree
for degree in range(polydegree):
X[:,degree] = Density**(degree/3.0)
- OLS = LinearRegression()
+ OLS = LinearRegression(fit_intercept=False)
# loop over trials in order to estimate the expectation value of the MSE
estimated_mse_folds = cross_val_score(OLS, X, Energies, scoring='neg_mean_squared_error', cv=kfold)
#[:, np.newaxis]
@@ -1944,50 +1948,6 @@ plt.show()
-
-
-The same example but now with cross-validation
+Cross-validation with Ridge
-import numpy as np
-import matplotlib.pyplot as plt
-from sklearn.model_selection import KFold
-from sklearn.linear_model import Ridge
-from sklearn.model_selection import cross_val_score
-from sklearn.preprocessing import PolynomialFeatures
-
-# A seed just to ensure that the random numbers are the same for every run.
-np.random.seed(3155)
-# Generate the data.
-n = 100
-x = np.linspace(-3, 3, n).reshape(-1, 1)
-y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)+ np.random.normal(0, 0.1, x.shape)
-# Decide degree on polynomial to fit
-poly = PolynomialFeatures(degree = 10)
-
-# Decide which values of lambda to use
-nlambdas = 500
-lambdas = np.logspace(-3, 5, nlambdas)
-# Initialize a KFold instance
-k = 5
-kfold = KFold(n_splits = k)
-estimated_mse_sklearn = np.zeros(nlambdas)
-i = 0
-for lmb in lambdas:
- ridge = Ridge(alpha = lmb)
- estimated_mse_folds = cross_val_score(ridge, x, y, scoring='neg_mean_squared_error', cv=kfold)
- estimated_mse_sklearn[i] = np.mean(-estimated_mse_folds)
- i += 1
-plt.figure()
-plt.plot(np.log10(lambdas), estimated_mse_sklearn, label = 'cross_val_score')
-plt.xlabel('log10(lambda)')
-plt.ylabel('MSE')
-plt.legend()
-plt.show()
-
-
@@ -1844,7 +1840,7 @@ trials = 100
trainingerror[polydegree] = 0.0
for samples in range(trials):
x_train, x_test, y_train, y_test = train_test_split(X, Energies, test_size=0.2)
- model = LinearRegression(fit_intercept=True).fit(x_train, y_train)
+ model = LinearRegression(fit_intercept=False).fit(x_train, y_train)
ypred = model.predict(x_train)
ytilde = model.predict(x_test)
testerror[polydegree] += mean_squared_error(y_test, ytilde)
@@ -1863,11 +1859,16 @@ plt.ylabel('log10[MSE]')
plt.legend()
plt.show()
+
+Note that we kept the intercept column in the fitting here. This means that we need to set the intercept in the call to the Scikit-Learn function as False. Alternatively, we could have set up the design matrix \( X \) without the first column of ones. +
+In this example we keep the intercept column again but add cross-validation in order to estimate the best possible value of the means squared error.
@@ -1927,7 +1928,7 @@ kfold = KFold(n_splits = k) polynomial[polydegree] = polydegree for degree in range(polydegree): X[:,degree] = Density**(degree/3.0) - OLS = LinearRegression() + OLS = LinearRegression(fit_intercept=False) # loop over trials in order to estimate the expectation value of the MSE estimated_mse_folds = cross_val_score(OLS, X, Energies, scoring='neg_mean_squared_error', cv=kfold) #[:, np.newaxis] @@ -1940,49 +1941,6 @@ plt.legend() plt.show()
-
-
-
- - -
import numpy as np
-import matplotlib.pyplot as plt
-from sklearn.model_selection import KFold
-from sklearn.linear_model import Ridge
-from sklearn.model_selection import cross_val_score
-from sklearn.preprocessing import PolynomialFeatures
-
-# A seed just to ensure that the random numbers are the same for every run.
-np.random.seed(3155)
-# Generate the data.
-n = 100
-x = np.linspace(-3, 3, n).reshape(-1, 1)
-y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)+ np.random.normal(0, 0.1, x.shape)
-# Decide degree on polynomial to fit
-poly = PolynomialFeatures(degree = 10)
-
-# Decide which values of lambda to use
-nlambdas = 500
-lambdas = np.logspace(-3, 5, nlambdas)
-# Initialize a KFold instance
-k = 5
-kfold = KFold(n_splits = k)
-estimated_mse_sklearn = np.zeros(nlambdas)
-i = 0
-for lmb in lambdas:
- ridge = Ridge(alpha = lmb)
- estimated_mse_folds = cross_val_score(ridge, x, y, scoring='neg_mean_squared_error', cv=kfold)
- estimated_mse_sklearn[i] = np.mean(-estimated_mse_folds)
- i += 1
-plt.figure()
-plt.plot(np.log10(lambdas), estimated_mse_sklearn, label = 'cross_val_score')
-plt.xlabel('log10(lambda)')
-plt.ylabel('MSE')
-plt.legend()
-plt.show()
-diff --git a/doc/pub/week37/html/week37.html b/doc/pub/week37/html/week37.html index cf8db4dde..148d17190 100644 --- a/doc/pub/week37/html/week37.html +++ b/doc/pub/week37/html/week37.html @@ -202,11 +202,7 @@ div { text-align: justify; text-justify: inter-word; } ('The same example but now with cross-validation', 2, None, - 'the-same-example-but-now-with-cross-validation'), - ('Cross-validation with Ridge', - 2, - None, - 'cross-validation-with-ridge')]} + 'the-same-example-but-now-with-cross-validation')]} end of tocinfo -->
@@ -248,7 +244,7 @@ MathJax.Hub.Config({-
+Note that we kept the intercept column in the fitting here. This means that we need to set the intercept in the call to the Scikit-Learn function as False. Alternatively, we could have set up the design matrix \( X \) without the first column of ones.
+
+In this example we keep the intercept column again but add cross-validation in order to estimate the best possible value of the means squared error.
@@ -1932,7 +1933,7 @@ kfold = KFold(n_splits = polydegree
for degree in range(polydegree):
X[:,degree] = Density**(degree/3.0)
- OLS = LinearRegression()
+ OLS = LinearRegression(fit_intercept=False)
# loop over trials in order to estimate the expectation value of the MSE
estimated_mse_folds = cross_val_score(OLS, X, Energies, scoring='neg_mean_squared_error', cv=kfold)
#[:, np.newaxis]
@@ -1945,49 +1946,6 @@ plt.legend()
plt.show()
-
-
-
-
diff --git a/doc/pub/week37/ipynb/ipynb-week37-src.tar.gz b/doc/pub/week37/ipynb/ipynb-week37-src.tar.gz
index ed785793e..d63efe22f 100644
Binary files a/doc/pub/week37/ipynb/ipynb-week37-src.tar.gz and b/doc/pub/week37/ipynb/ipynb-week37-src.tar.gz differ
diff --git a/doc/pub/week37/ipynb/week37.ipynb b/doc/pub/week37/ipynb/week37.ipynb
index ab0f83f40..15af1b5f5 100644
--- a/doc/pub/week37/ipynb/week37.ipynb
+++ b/doc/pub/week37/ipynb/week37.ipynb
@@ -10,7 +10,7 @@
" \n",
"**Morten Hjorth-Jensen**, Department of Physics, University of Oslo and Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University\n",
"\n",
- "Date: **Sep 17, 2021**\n",
+ "Date: **Sep 28, 2021**\n",
"\n",
"Copyright 1999-2021, Morten Hjorth-Jensen. Released under CC Attribution-NonCommercial 4.0 license\n",
"\n",
@@ -1984,7 +1984,7 @@
" trainingerror[polydegree] = 0.0\n",
" for samples in range(trials):\n",
" x_train, x_test, y_train, y_test = train_test_split(X, Energies, test_size=0.2)\n",
- " model = LinearRegression(fit_intercept=True).fit(x_train, y_train)\n",
+ " model = LinearRegression(fit_intercept=False).fit(x_train, y_train)\n",
" ypred = model.predict(x_train)\n",
" ytilde = model.predict(x_test)\n",
" testerror[polydegree] += mean_squared_error(y_test, ytilde)\n",
@@ -2008,8 +2008,12 @@
"cell_type": "markdown",
"metadata": {},
"source": [
+ "Note that we kept the intercept column in the fitting here. This means that we need to set the **intercept** in the call to the **Scikit-Learn** function as **False**. Alternatively, we could have set up the design matrix $X$ without the first column of ones.\n",
+ "\n",
"\n",
- "## The same example but now with cross-validation"
+ "## The same example but now with cross-validation\n",
+ "\n",
+ "In this example we keep the intercept column again but add cross-validation in order to estimate the best possible value of the means squared error."
]
},
{
@@ -2077,7 +2081,7 @@
" polynomial[polydegree] = polydegree\n",
" for degree in range(polydegree):\n",
" X[:,degree] = Density**(degree/3.0)\n",
- " OLS = LinearRegression()\n",
+ " OLS = LinearRegression(fit_intercept=False)\n",
"# loop over trials in order to estimate the expectation value of the MSE\n",
" estimated_mse_folds = cross_val_score(OLS, X, Energies, scoring='neg_mean_squared_error', cv=kfold)\n",
"#[:, np.newaxis]\n",
@@ -2089,59 +2093,6 @@
"plt.legend()\n",
"plt.show()"
]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Cross-validation with Ridge"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "collapsed": false,
- "editable": true
- },
- "outputs": [],
- "source": [
- "import numpy as np\n",
- "import matplotlib.pyplot as plt\n",
- "from sklearn.model_selection import KFold\n",
- "from sklearn.linear_model import Ridge\n",
- "from sklearn.model_selection import cross_val_score\n",
- "from sklearn.preprocessing import PolynomialFeatures\n",
- "\n",
- "# A seed just to ensure that the random numbers are the same for every run.\n",
- "np.random.seed(3155)\n",
- "# Generate the data.\n",
- "n = 100\n",
- "x = np.linspace(-3, 3, n).reshape(-1, 1)\n",
- "y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)+ np.random.normal(0, 0.1, x.shape)\n",
- "# Decide degree on polynomial to fit\n",
- "poly = PolynomialFeatures(degree = 10)\n",
- "\n",
- "# Decide which values of lambda to use\n",
- "nlambdas = 500\n",
- "lambdas = np.logspace(-3, 5, nlambdas)\n",
- "# Initialize a KFold instance\n",
- "k = 5\n",
- "kfold = KFold(n_splits = k)\n",
- "estimated_mse_sklearn = np.zeros(nlambdas)\n",
- "i = 0\n",
- "for lmb in lambdas:\n",
- " ridge = Ridge(alpha = lmb)\n",
- " estimated_mse_folds = cross_val_score(ridge, x, y, scoring='neg_mean_squared_error', cv=kfold)\n",
- " estimated_mse_sklearn[i] = np.mean(-estimated_mse_folds)\n",
- " i += 1\n",
- "plt.figure()\n",
- "plt.plot(np.log10(lambdas), estimated_mse_sklearn, label = 'cross_val_score')\n",
- "plt.xlabel('log10(lambda)')\n",
- "plt.ylabel('MSE')\n",
- "plt.legend()\n",
- "plt.show()"
- ]
}
],
"metadata": {},
diff --git a/doc/src/week37/week37.do.txt b/doc/src/week37/week37.do.txt
index 4cf333222..c16dd66b7 100644
--- a/doc/src/week37/week37.do.txt
+++ b/doc/src/week37/week37.do.txt
@@ -1485,7 +1485,7 @@ for polydegree in range(1, Maxpolydegree):
trainingerror[polydegree] = 0.0
for samples in range(trials):
x_train, x_test, y_train, y_test = train_test_split(X, Energies, test_size=0.2)
- model = LinearRegression(fit_intercept=True).fit(x_train, y_train)
+ model = LinearRegression(fit_intercept=False).fit(x_train, y_train)
ypred = model.predict(x_train)
ytilde = model.predict(x_test)
testerror[polydegree] += mean_squared_error(y_test, ytilde)
@@ -1506,10 +1506,12 @@ plt.show()
!ec
+Note that we kept the intercept column in the fitting here. This means that we need to set the _intercept_ in the call to the _Scikit-Learn_ function as _False_. Alternatively, we could have set up the design matrix $X$ without the first column of ones.
!split
===== The same example but now with cross-validation =====
+In this example we keep the intercept column again but add cross-validation in order to estimate the best possible value of the means squared error.
!bc pycod
# Common imports
import os
@@ -1567,7 +1569,7 @@ for polydegree in range(1, Maxpolydegree):
polynomial[polydegree] = polydegree
for degree in range(polydegree):
X[:,degree] = Density**(degree/3.0)
- OLS = LinearRegression()
+ OLS = LinearRegression(fit_intercept=False)
# loop over trials in order to estimate the expectation value of the MSE
estimated_mse_folds = cross_val_score(OLS, X, Energies, scoring='neg_mean_squared_error', cv=kfold)
#[:, np.newaxis]
@@ -1581,52 +1583,3 @@ plt.show()
!ec
-!split
-===== Cross-validation with Ridge =====
-!bc pycod
-import numpy as np
-import matplotlib.pyplot as plt
-from sklearn.model_selection import KFold
-from sklearn.linear_model import Ridge
-from sklearn.model_selection import cross_val_score
-from sklearn.preprocessing import PolynomialFeatures
-
-# A seed just to ensure that the random numbers are the same for every run.
-np.random.seed(3155)
-# Generate the data.
-n = 100
-x = np.linspace(-3, 3, n).reshape(-1, 1)
-y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)+ np.random.normal(0, 0.1, x.shape)
-# Decide degree on polynomial to fit
-poly = PolynomialFeatures(degree = 10)
-
-# Decide which values of lambda to use
-nlambdas = 500
-lambdas = np.logspace(-3, 5, nlambdas)
-# Initialize a KFold instance
-k = 5
-kfold = KFold(n_splits = k)
-estimated_mse_sklearn = np.zeros(nlambdas)
-i = 0
-for lmb in lambdas:
- ridge = Ridge(alpha = lmb)
- estimated_mse_folds = cross_val_score(ridge, x, y, scoring='neg_mean_squared_error', cv=kfold)
- estimated_mse_sklearn[i] = np.mean(-estimated_mse_folds)
- i += 1
-plt.figure()
-plt.plot(np.log10(lambdas), estimated_mse_sklearn, label = 'cross_val_score')
-plt.xlabel('log10(lambda)')
-plt.ylabel('MSE')
-plt.legend()
-plt.show()
-
-
-!ec
-
-
-
-
-
-
-
-
diff --git a/doc/src/week39/test1.py b/doc/src/week39/test1.py
new file mode 100644
index 000000000..1a9ead338
--- /dev/null
+++ b/doc/src/week39/test1.py
@@ -0,0 +1,16 @@
+import numpy as np
+from sklearn import datasets
+from sklearn.linear_model import Ridge
+from sklearn.model_selection import GridSearchCV
+# load the diabetes datasets
+dataset = datasets.load_diabetes()
+# prepare a range of alpha values to test
+alphas = np.array([1,0.1,0.01,0.001,0.0001,0])
+# create and fit a ridge regression model, testing each alpha
+model = Ridge()
+grid = GridSearchCV(estimator=model, param_grid=dict(alpha=alphas))
+grid.fit(dataset.data, dataset.target)
+print(grid)
+# summarize the results of the grid search
+print(grid.best_score_)
+print(grid.best_estimator_.alpha)
diff --git a/doc/src/week39/test2.py b/doc/src/week39/test2.py
new file mode 100644
index 000000000..8a58d4ce8
--- /dev/null
+++ b/doc/src/week39/test2.py
@@ -0,0 +1,17 @@
+import numpy as np
+from scipy.stats import uniform as sp_rand
+from sklearn import datasets
+from sklearn.linear_model import Ridge
+from sklearn.model_selection import RandomizedSearchCV
+# load the diabetes datasets
+dataset = datasets.load_diabetes()
+# prepare a uniform distribution to sample for the alpha parameter
+param_grid = {'alpha': sp_rand()}
+# create and fit a ridge regression model, testing random alpha values
+model = Ridge()
+rsearch = RandomizedSearchCV(estimator=model, param_distributions=param_grid, n_iter=100)
+rsearch.fit(dataset.data, dataset.target)
+print(rsearch)
+# summarize the results of the random parameter search
+print(rsearch.best_score_)
+print(rsearch.best_estimator_.alpha)
diff --git a/doc/src/week39/test3.py b/doc/src/week39/test3.py
new file mode 100644
index 000000000..dc4d67de8
--- /dev/null
+++ b/doc/src/week39/test3.py
@@ -0,0 +1,30 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from sklearn.model_selection import KFold
+from sklearn.linear_model import Ridge
+from sklearn.model_selection import cross_val_score
+from sklearn.preprocessing import PolynomialFeatures
+from sklearn.model_selection import GridSearchCV
+
+# A seed just to ensure that the random numbers are the same for every run.
+np.random.seed(3155)
+# Generate the data.
+n = 100
+x = np.linspace(-3, 3, n).reshape(-1, 1)
+y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)+ np.random.normal(0, 0.1, x.shape)
+# Decide degree on polynomial to fit
+poly = PolynomialFeatures(degree = 10)
+
+# Decide which values of lambda to use
+nlambdas = 10
+lambdas = np.logspace(-3, 3, nlambdas)
+
+# create and fit a ridge regression model, testing each alpha
+model = Ridge()
+grid = GridSearchCV(estimator=model, param_grid=dict(alpha=lambdas))
+grid.fit(x, y)
+print(grid)
+# summarize the results of the grid search
+print(grid.best_score_)
+print(grid.best_estimator_.alpha)
+
diff --git a/doc/src/week39/test4.py b/doc/src/week39/test4.py
new file mode 100644
index 000000000..82e34dc98
--- /dev/null
+++ b/doc/src/week39/test4.py
@@ -0,0 +1,35 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from sklearn.model_selection import KFold
+from sklearn.linear_model import Ridge
+from sklearn.model_selection import cross_val_score
+from sklearn.preprocessing import PolynomialFeatures
+
+# A seed just to ensure that the random numbers are the same for every run.
+np.random.seed(3155)
+# Generate the data.
+n = 100
+x = np.linspace(-3, 3, n).reshape(-1, 1)
+y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)+ np.random.normal(0, 0.1, x.shape)
+# Decide degree on polynomial to fit
+#poly = PolynomialFeatures(degree = 10)
+
+# Decide which values of lambda to use
+nlambdas = 500
+lambdas = np.logspace(-3, 5, nlambdas)
+# Initialize a KFold instance
+k = 5
+kfold = KFold(n_splits = k)
+estimated_mse_sklearn = np.zeros(nlambdas)
+i = 0
+for lmb in lambdas:
+ ridge = Ridge(alpha = lmb)
+ estimated_mse_folds = cross_val_score(ridge, x, y, scoring='neg_mean_squared_error', cv=kfold)
+ estimated_mse_sklearn[i] = np.mean(-estimated_mse_folds)
+ i += 1
+plt.figure()
+plt.plot(np.log10(lambdas), estimated_mse_sklearn, label = 'cross_val_score')
+plt.xlabel('log10(lambda)')
+plt.ylabel('MSE')
+plt.legend()
+plt.show()
diff --git a/doc/src/week39/test5.py b/doc/src/week39/test5.py
new file mode 100644
index 000000000..430ab655c
--- /dev/null
+++ b/doc/src/week39/test5.py
@@ -0,0 +1,61 @@
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+from sklearn.model_selection import train_test_split
+from sklearn.linear_model import Ridge
+from sklearn.preprocessing import StandardScaler
+from sklearn.model_selection import cross_val_score
+from sklearn.model_selection import GridSearchCV
+
+
+
+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)
+#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)
+#Remove the intercept from the training data.
+y_scaler = np.mean(y_train)
+y_train_scaled = y_train - y_scaler
+
+p = Maxpolydegree-1
+I = np.eye(p,p)
+# Decide which values of lambda to use
+nlambdas = 10
+MSEOwnRidgePredict = np.zeros(nlambdas)
+MSERidgePredict = np.zeros(nlambdas)
+
+lambdas = np.logspace(-4, 2, nlambdas)
+
+# create and fit a ridge regression model, testing each alpha
+model = Ridge()
+grid = GridSearchCV(estimator=model, param_grid=dict(alpha=lambdas))
+grid.fit(X_train_scaled, y_train_scaled)
+print(grid)
+ypredictRidge = grid.predict(X_test_scaled)
+# summarize the results of the grid search
+print(grid.best_score_)
+print(grid.best_estimator_.alpha)
+print(MSE(y_test,ypredictRidge))
diff --git a/doc/src/week39/test6.py b/doc/src/week39/test6.py
new file mode 100644
index 000000000..b5515f583
--- /dev/null
+++ b/doc/src/week39/test6.py
@@ -0,0 +1,72 @@
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+from sklearn.model_selection import train_test_split
+from sklearn.linear_model import Ridge
+from sklearn.preprocessing import StandardScaler
+from sklearn.model_selection import cross_val_score
+from sklearn.model_selection import GridSearchCV
+from scipy.stats import uniform as sp_rand
+from sklearn.model_selection import RandomizedSearchCV
+
+
+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)
+#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)
+#Remove the intercept from the training data.
+y_scaler = np.mean(y_train)
+y_train_scaled = y_train - y_scaler
+y_test_scaled = y_test - y_scaler
+
+p = Maxpolydegree-1
+I = np.eye(p,p)
+# Decide which values of lambda to use
+nlambdas = 10
+MSEOwnRidgePredict = np.zeros(nlambdas)
+MSERidgePredict = np.zeros(nlambdas)
+
+
+param_grid = {'alpha': sp_rand()}
+# create and fit a ridge regression model, testing each alpha
+model = Ridge()
+grid = RandomizedSearchCV(estimator=model, param_distributions=param_grid, n_iter=100)
+#GridSearchCV(estimator=model, param_grid=dict(alpha=lambdas))
+grid.fit(X_train_scaled, y_train_scaled)
+print(grid)
+ypredictRidge = grid.predict(X_test_scaled)
+# summarize the results of the grid search
+print(grid.best_score_)
+print(grid.best_estimator_.alpha)
+print(MSE(y_test_scaled,ypredictRidge))
+print(R2(y_test_scaled,ypredictRidge))
+
+
+
diff --git a/doc/src/week39/week39.do.txt b/doc/src/week39/week39.do.txt
index 8fc273586..25f3525dc 100644
--- a/doc/src/week39/week39.do.txt
+++ b/doc/src/week39/week39.do.txt
@@ -31,41 +31,6 @@ in evaluating the MSE as function of different $\lambda$ values.
Based on these calculations, one tries then to determine the value of the hyperparameter $\lambda$
which results in optimal scores (for example the smallest MSE or an $R2=1$).
!bc pycod
-import numpy as np
-import matplotlib.pyplot as plt
-from sklearn.model_selection import KFold
-from sklearn.linear_model import Ridge
-from sklearn.model_selection import cross_val_score
-from sklearn.preprocessing import PolynomialFeatures
-
-# A seed just to ensure that the random numbers are the same for every run.
-np.random.seed(3155)
-# Generate the data.
-n = 100
-x = np.linspace(-3, 3, n).reshape(-1, 1)
-y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)+ np.random.normal(0, 0.1, x.shape)
-# Decide degree on polynomial to fit
-poly = PolynomialFeatures(degree = 10)
-
-# Decide which values of lambda to use
-nlambdas = 500
-lambdas = np.logspace(-3, 5, nlambdas)
-# Initialize a KFold instance
-k = 5
-kfold = KFold(n_splits = k)
-estimated_mse_sklearn = np.zeros(nlambdas)
-i = 0
-for lmb in lambdas:
- ridge = Ridge(alpha = lmb)
- estimated_mse_folds = cross_val_score(ridge, x, y, scoring='neg_mean_squared_error', cv=kfold)
- estimated_mse_sklearn[i] = np.mean(-estimated_mse_folds)
- i += 1
-plt.figure()
-plt.plot(np.log10(lambdas), estimated_mse_sklearn, label = 'cross_val_score')
-plt.xlabel('log10(lambda)')
-plt.ylabel('MSE')
-plt.legend()
-plt.show()
!ec
We see from this plot that the optimal MSE occurs for a value of $\lambda\in [10,100]$.
@@ -80,41 +45,7 @@ included with the library _Scikit-Learn_, as demonstrated for the same
example here.
!bc pycod
-import numpy as np
-import matplotlib.pyplot as plt
-from sklearn.model_selection import KFold
-from sklearn.linear_model import Ridge
-from sklearn.model_selection import cross_val_score
-from sklearn.preprocessing import PolynomialFeatures
-# A seed just to ensure that the random numbers are the same for every run.
-np.random.seed(3155)
-# Generate the data.
-n = 100
-x = np.linspace(-3, 3, n).reshape(-1, 1)
-y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)+ np.random.normal(0, 0.1, x.shape)
-# Decide degree on polynomial to fit
-poly = PolynomialFeatures(degree = 10)
-
-# Decide which values of lambda to use
-nlambdas = 500
-lambdas = np.logspace(-3, 5, nlambdas)
-# Initialize a KFold instance
-k = 5
-kfold = KFold(n_splits = k)
-estimated_mse_sklearn = np.zeros(nlambdas)
-i = 0
-for lmb in lambdas:
- ridge = Ridge(alpha = lmb)
- estimated_mse_folds = cross_val_score(ridge, x, y, scoring='neg_mean_squared_error', cv=kfold)
- estimated_mse_sklearn[i] = np.mean(-estimated_mse_folds)
- i += 1
-plt.figure()
-plt.plot(np.log10(lambdas), estimated_mse_sklearn, label = 'cross_val_score')
-plt.xlabel('log10(lambda)')
-plt.ylabel('MSE')
-plt.legend()
-plt.show()
!ec
@@ -1849,7 +1845,7 @@ trials = 100= 0.0
for samples in range(trials):
x_train, x_test, y_train, y_test = train_test_split(X, Energies, test_size=0.2)
- model = LinearRegression(fit_intercept=True).fit(x_train, y_train)
+ model = LinearRegression(fit_intercept=False).fit(x_train, y_train)
ypred = model.predict(x_train)
ytilde = model.predict(x_test)
testerror[polydegree] += mean_squared_error(y_test, ytilde)
@@ -1868,11 +1864,16 @@ plt.ylabel('
plt.legend()
plt.show()
+The same example but now with cross-validation
+
-
-Cross-validation with Ridge
-import numpy as np
-import matplotlib.pyplot as plt
-from sklearn.model_selection import KFold
-from sklearn.linear_model import Ridge
-from sklearn.model_selection import cross_val_score
-from sklearn.preprocessing import PolynomialFeatures
-
-# A seed just to ensure that the random numbers are the same for every run.
-np.random.seed(3155)
-# Generate the data.
-n = 100
-x = np.linspace(-3, 3, n).reshape(-1, 1)
-y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)+ np.random.normal(0, 0.1, x.shape)
-# Decide degree on polynomial to fit
-poly = PolynomialFeatures(degree = 10)
-
-# Decide which values of lambda to use
-nlambdas = 500
-lambdas = np.logspace(-3, 5, nlambdas)
-# Initialize a KFold instance
-k = 5
-kfold = KFold(n_splits = k)
-estimated_mse_sklearn = np.zeros(nlambdas)
-i = 0
-for lmb in lambdas:
- ridge = Ridge(alpha = lmb)
- estimated_mse_folds = cross_val_score(ridge, x, y, scoring='neg_mean_squared_error', cv=kfold)
- estimated_mse_sklearn[i] = np.mean(-estimated_mse_folds)
- i += 1
-plt.figure()
-plt.plot(np.log10(lambdas), estimated_mse_sklearn, label = 'cross_val_score')
-plt.xlabel('log10(lambda)')
-plt.ylabel('MSE')
-plt.legend()
-plt.show()
-