|
|
|
@@ -11,12 +11,913 @@ DATE: today
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
!split
|
|
|
|
|
===== Thursday: =====
|
|
|
|
|
===== Thursday September 17 =====
|
|
|
|
|
|
|
|
|
|
Material will added during Wednesday 16
|
|
|
|
|
|
|
|
|
|
!split
|
|
|
|
|
===== Friday: Intro to Logistic Regression =====
|
|
|
|
|
===== Ridge and LASSO Regression, reminder =====
|
|
|
|
|
|
|
|
|
|
The expression for the standard Mean Squared Error (MSE) which we used to define our cost function and the equations for the ordinary least squares (OLS) method, that is
|
|
|
|
|
our optimization problem is
|
|
|
|
|
!bt
|
|
|
|
|
\[
|
|
|
|
|
{\displaystyle \min_{\bm{\beta}\in {\mathbb{R}}^{p}}}\frac{1}{n}\left\{\left(\bm{y}-\bm{X}\bm{\beta}\right)^T\left(\bm{y}-\bm{X}\bm{\beta}\right)\right\}.
|
|
|
|
|
\]
|
|
|
|
|
!et
|
|
|
|
|
or we can state it as
|
|
|
|
|
!bt
|
|
|
|
|
\[
|
|
|
|
|
{\displaystyle \min_{\bm{\beta}\in
|
|
|
|
|
{\mathbb{R}}^{p}}}\frac{1}{n}\sum_{i=0}^{n-1}\left(y_i-\tilde{y}_i\right)^2=\frac{1}{n}\vert\vert \bm{y}-\bm{X}\bm{\beta}\vert\vert_2^2,
|
|
|
|
|
\]
|
|
|
|
|
!et
|
|
|
|
|
where we have used the definition of a norm-2 vector, that is
|
|
|
|
|
!bt
|
|
|
|
|
\[
|
|
|
|
|
\vert\vert \bm{x}\vert\vert_2 = \sqrt{\sum_i x_i^2}.
|
|
|
|
|
\]
|
|
|
|
|
!et
|
|
|
|
|
|
|
|
|
|
By minimizing the above equation with respect to the parameters
|
|
|
|
|
$\bm{\beta}$ we could then obtain an analytical expression for the
|
|
|
|
|
parameters $\bm{\beta}$. We can add a regularization parameter $\lambda$ by
|
|
|
|
|
defining a new cost function to be optimized, that is
|
|
|
|
|
|
|
|
|
|
!bt
|
|
|
|
|
\[
|
|
|
|
|
{\displaystyle \min_{\bm{\beta}\in
|
|
|
|
|
{\mathbb{R}}^{p}}}\frac{1}{n}\vert\vert \bm{y}-\bm{X}\bm{\beta}\vert\vert_2^2+\lambda\vert\vert \bm{\beta}\vert\vert_2^2
|
|
|
|
|
\]
|
|
|
|
|
!et
|
|
|
|
|
|
|
|
|
|
which leads to the Ridge regression minimization problem where we
|
|
|
|
|
require that $\vert\vert \bm{\beta}\vert\vert_2^2\le t$, where $t$ is
|
|
|
|
|
a finite number larger than zero. By defining
|
|
|
|
|
|
|
|
|
|
!bt
|
|
|
|
|
\[
|
|
|
|
|
C(\bm{X},\bm{\beta})=\frac{1}{n}\vert\vert \bm{y}-\bm{X}\bm{\beta}\vert\vert_2^2+\lambda\vert\vert \bm{\beta}\vert\vert_1,
|
|
|
|
|
\]
|
|
|
|
|
!et
|
|
|
|
|
|
|
|
|
|
we have a new optimization equation
|
|
|
|
|
!bt
|
|
|
|
|
\[
|
|
|
|
|
{\displaystyle \min_{\bm{\beta}\in
|
|
|
|
|
{\mathbb{R}}^{p}}}\frac{1}{n}\vert\vert \bm{y}-\bm{X}\bm{\beta}\vert\vert_2^2+\lambda\vert\vert \bm{\beta}\vert\vert_1
|
|
|
|
|
\]
|
|
|
|
|
!et
|
|
|
|
|
which leads to Lasso regression. Lasso stands for least absolute shrinkage and selection operator.
|
|
|
|
|
|
|
|
|
|
Here we have defined the norm-1 as
|
|
|
|
|
!bt
|
|
|
|
|
\[
|
|
|
|
|
\vert\vert \bm{x}\vert\vert_1 = \sum_i \vert x_i\vert.
|
|
|
|
|
\]
|
|
|
|
|
!et
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
!split
|
|
|
|
|
===== Various steps in cross-validation =====
|
|
|
|
|
|
|
|
|
|
When the repetitive splitting of the data set is done randomly,
|
|
|
|
|
samples may accidently end up in a fast majority of the splits in
|
|
|
|
|
either training or test set. Such samples may have an unbalanced
|
|
|
|
|
influence on either model building or prediction evaluation. To avoid
|
|
|
|
|
this $k$-fold cross-validation structures the data splitting. The
|
|
|
|
|
samples are divided into $k$ more or less equally sized exhaustive and
|
|
|
|
|
mutually exclusive subsets. In turn (at each split) one of these
|
|
|
|
|
subsets plays the role of the test set while the union of the
|
|
|
|
|
remaining subsets constitutes the training set. Such a splitting
|
|
|
|
|
warrants a balanced representation of each sample in both training and
|
|
|
|
|
test set over the splits. Still the division into the $k$ subsets
|
|
|
|
|
involves a degree of randomness. This may be fully excluded when
|
|
|
|
|
choosing $k=n$. This particular case is referred to as leave-one-out
|
|
|
|
|
cross-validation (LOOCV).
|
|
|
|
|
|
|
|
|
|
!split
|
|
|
|
|
===== How to set up the cross-validation for Ridge and/or Lasso =====
|
|
|
|
|
|
|
|
|
|
* Define a range of interest for the penalty parameter.
|
|
|
|
|
|
|
|
|
|
* Divide the data set into training and test set comprising samples $\{1, \ldots, n\} \setminus i$ and $\{ i \}$, respectively.
|
|
|
|
|
|
|
|
|
|
* Fit the linear regression model by means of ridge estimation for each $\lambda$ in the grid using the training set, and the corresponding estimate of the error variance $\bm{\sigma}_{-i}^2(\lambda)$, as
|
|
|
|
|
!bt
|
|
|
|
|
\begin{align*}
|
|
|
|
|
\bm{\beta}_{-i}(\lambda) & = ( \bm{X}_{-i, \ast}^{T}
|
|
|
|
|
\bm{X}_{-i, \ast} + \lambda \bm{I}_{pp})^{-1}
|
|
|
|
|
\bm{X}_{-i, \ast}^{T} \bm{y}_{-i}
|
|
|
|
|
\end{align*}
|
|
|
|
|
!et
|
|
|
|
|
|
|
|
|
|
* Evaluate the prediction performance of these models on the test set by $\log\{L[y_i, \bm{X}_{i, \ast}; \bm{\beta}_{-i}(\lambda), \bm{\sigma}_{-i}^2(\lambda)]\}$. Or, by the prediction error $|y_i - \bm{X}_{i, \ast} \bm{\beta}_{-i}(\lambda)|$, the relative error, the error squared or the R2 score function.
|
|
|
|
|
|
|
|
|
|
* Repeat the first three steps such that each sample plays the role of the test set once.
|
|
|
|
|
|
|
|
|
|
* Average the prediction performances of the test sets at each grid point of the penalty bias/parameter. It is an estimate of the prediction performance of the model corresponding to this value of the penalty parameter on novel data. It is defined as
|
|
|
|
|
!bt
|
|
|
|
|
\begin{align*}
|
|
|
|
|
\frac{1}{n} \sum_{i = 1}^n \log\{L[y_i, \mathbf{X}_{i, \ast}; \bm{\beta}_{-i}(\lambda), \bm{\sigma}_{-i}^2(\lambda)]\}.
|
|
|
|
|
\end{align*}
|
|
|
|
|
!et
|
|
|
|
|
|
|
|
|
|
!split
|
|
|
|
|
===== Cross-validation in brief =====
|
|
|
|
|
|
|
|
|
|
For the various values of $k$
|
|
|
|
|
|
|
|
|
|
o shuffle the dataset randomly.
|
|
|
|
|
o Split the dataset into $k$ groups.
|
|
|
|
|
o For each unique group:
|
|
|
|
|
o Decide which group to use as set for test data
|
|
|
|
|
o Take the remaining groups as a training data set
|
|
|
|
|
o Fit a model on the training set and evaluate it on the test set
|
|
|
|
|
o Retain the evaluation score and discard the model
|
|
|
|
|
o Summarize the model using the sample of model evaluation scores
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
!split
|
|
|
|
|
===== Code Example for Cross-validation and $k$-fold Cross-validation =====
|
|
|
|
|
|
|
|
|
|
The code here uses Ridge regression with cross-validation (CV) resampling and $k$-fold CV in order to fit a specific polynomial.
|
|
|
|
|
!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.
|
|
|
|
|
# Useful for eventual debugging.
|
|
|
|
|
np.random.seed(3155)
|
|
|
|
|
|
|
|
|
|
# Generate the data.
|
|
|
|
|
nsamples = 100
|
|
|
|
|
x = np.random.randn(nsamples)
|
|
|
|
|
y = 3*x**2 + np.random.randn(nsamples)
|
|
|
|
|
|
|
|
|
|
## Cross-validation on Ridge regression using KFold only
|
|
|
|
|
|
|
|
|
|
# Decide degree on polynomial to fit
|
|
|
|
|
poly = PolynomialFeatures(degree = 6)
|
|
|
|
|
|
|
|
|
|
# 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)
|
|
|
|
|
|
|
|
|
|
# Perform the cross-validation to estimate MSE
|
|
|
|
|
scores_KFold = np.zeros((nlambdas, k))
|
|
|
|
|
|
|
|
|
|
i = 0
|
|
|
|
|
for lmb in lambdas:
|
|
|
|
|
ridge = Ridge(alpha = lmb)
|
|
|
|
|
j = 0
|
|
|
|
|
for train_inds, test_inds in kfold.split(x):
|
|
|
|
|
xtrain = x[train_inds]
|
|
|
|
|
ytrain = y[train_inds]
|
|
|
|
|
|
|
|
|
|
xtest = x[test_inds]
|
|
|
|
|
ytest = y[test_inds]
|
|
|
|
|
|
|
|
|
|
Xtrain = poly.fit_transform(xtrain[:, np.newaxis])
|
|
|
|
|
ridge.fit(Xtrain, ytrain[:, np.newaxis])
|
|
|
|
|
|
|
|
|
|
Xtest = poly.fit_transform(xtest[:, np.newaxis])
|
|
|
|
|
ypred = ridge.predict(Xtest)
|
|
|
|
|
|
|
|
|
|
scores_KFold[i,j] = np.sum((ypred - ytest[:, np.newaxis])**2)/np.size(ypred)
|
|
|
|
|
|
|
|
|
|
j += 1
|
|
|
|
|
i += 1
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
estimated_mse_KFold = np.mean(scores_KFold, axis = 1)
|
|
|
|
|
|
|
|
|
|
## Cross-validation using cross_val_score from sklearn along with KFold
|
|
|
|
|
|
|
|
|
|
# kfold is an instance initialized above as:
|
|
|
|
|
# kfold = KFold(n_splits = k)
|
|
|
|
|
|
|
|
|
|
estimated_mse_sklearn = np.zeros(nlambdas)
|
|
|
|
|
i = 0
|
|
|
|
|
for lmb in lambdas:
|
|
|
|
|
ridge = Ridge(alpha = lmb)
|
|
|
|
|
|
|
|
|
|
X = poly.fit_transform(x[:, np.newaxis])
|
|
|
|
|
estimated_mse_folds = cross_val_score(ridge, X, y[:, np.newaxis], scoring='neg_mean_squared_error', cv=kfold)
|
|
|
|
|
|
|
|
|
|
# cross_val_score return an array containing the estimated negative mse for every fold.
|
|
|
|
|
# we have to the the mean of every array in order to get an estimate of the mse of the model
|
|
|
|
|
estimated_mse_sklearn[i] = np.mean(-estimated_mse_folds)
|
|
|
|
|
|
|
|
|
|
i += 1
|
|
|
|
|
|
|
|
|
|
## Plot and compare the slightly different ways to perform cross-validation
|
|
|
|
|
|
|
|
|
|
plt.figure()
|
|
|
|
|
|
|
|
|
|
plt.plot(np.log10(lambdas), estimated_mse_sklearn, label = 'cross_val_score')
|
|
|
|
|
plt.plot(np.log10(lambdas), estimated_mse_KFold, 'r--', label = 'KFold')
|
|
|
|
|
|
|
|
|
|
plt.xlabel('log10(lambda)')
|
|
|
|
|
plt.ylabel('mse')
|
|
|
|
|
|
|
|
|
|
plt.legend()
|
|
|
|
|
|
|
|
|
|
plt.show()
|
|
|
|
|
|
|
|
|
|
!ec
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
!split
|
|
|
|
|
===== Bias-Variance tradeoff with Bootstrap =====
|
|
|
|
|
!bc pycod
|
|
|
|
|
import matplotlib.pyplot as plt
|
|
|
|
|
import numpy as np
|
|
|
|
|
from sklearn.linear_model import LinearRegression, Ridge, Lasso
|
|
|
|
|
from sklearn.preprocessing import PolynomialFeatures
|
|
|
|
|
from sklearn.model_selection import train_test_split
|
|
|
|
|
from sklearn.pipeline import make_pipeline
|
|
|
|
|
from sklearn.utils import resample
|
|
|
|
|
|
|
|
|
|
np.random.seed(2018)
|
|
|
|
|
|
|
|
|
|
n = 40
|
|
|
|
|
n_boostraps = 100
|
|
|
|
|
maxdegree = 14
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
# Make data set.
|
|
|
|
|
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)
|
|
|
|
|
error = np.zeros(maxdegree)
|
|
|
|
|
bias = np.zeros(maxdegree)
|
|
|
|
|
variance = np.zeros(maxdegree)
|
|
|
|
|
polydegree = np.zeros(maxdegree)
|
|
|
|
|
x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.2)
|
|
|
|
|
|
|
|
|
|
for degree in range(maxdegree):
|
|
|
|
|
model = make_pipeline(PolynomialFeatures(degree=degree), LinearRegression(fit_intercept=False))
|
|
|
|
|
y_pred = np.empty((y_test.shape[0], n_boostraps))
|
|
|
|
|
for i in range(n_boostraps):
|
|
|
|
|
x_, y_ = resample(x_train, y_train)
|
|
|
|
|
y_pred[:, i] = model.fit(x_, y_).predict(x_test).ravel()
|
|
|
|
|
|
|
|
|
|
polydegree[degree] = degree
|
|
|
|
|
error[degree] = np.mean( np.mean((y_test - y_pred)**2, axis=1, keepdims=True) )
|
|
|
|
|
bias[degree] = np.mean( (y_test - np.mean(y_pred, axis=1, keepdims=True))**2 )
|
|
|
|
|
variance[degree] = np.mean( np.var(y_pred, axis=1, keepdims=True) )
|
|
|
|
|
print('Polynomial degree:', degree)
|
|
|
|
|
print('Error:', error[degree])
|
|
|
|
|
print('Bias^2:', bias[degree])
|
|
|
|
|
print('Var:', variance[degree])
|
|
|
|
|
print('{} >= {} + {} = {}'.format(error[degree], bias[degree], variance[degree], bias[degree]+variance[degree]))
|
|
|
|
|
|
|
|
|
|
plt.plot(polydegree, error, label='Error')
|
|
|
|
|
plt.plot(polydegree, bias, label='bias')
|
|
|
|
|
plt.plot(polydegree, variance, label='Variance')
|
|
|
|
|
plt.legend()
|
|
|
|
|
plt.show()
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
!ec
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
!split
|
|
|
|
|
===== Another Example from Scikit-Learn's Repository =====
|
|
|
|
|
!bc pycod
|
|
|
|
|
"""
|
|
|
|
|
============================
|
|
|
|
|
Underfitting vs. Overfitting
|
|
|
|
|
============================
|
|
|
|
|
|
|
|
|
|
This example demonstrates the problems of underfitting and overfitting and
|
|
|
|
|
how we can use linear regression with polynomial features to approximate
|
|
|
|
|
nonlinear functions. The plot shows the function that we want to approximate,
|
|
|
|
|
which is a part of the cosine function. In addition, the samples from the
|
|
|
|
|
real function and the approximations of different models are displayed. The
|
|
|
|
|
models have polynomial features of different degrees. We can see that a
|
|
|
|
|
linear function (polynomial with degree 1) is not sufficient to fit the
|
|
|
|
|
training samples. This is called **underfitting**. A polynomial of degree 4
|
|
|
|
|
approximates the true function almost perfectly. However, for higher degrees
|
|
|
|
|
the model will **overfit** the training data, i.e. it learns the noise of the
|
|
|
|
|
training data.
|
|
|
|
|
We evaluate quantitatively **overfitting** / **underfitting** by using
|
|
|
|
|
cross-validation. We calculate the mean squared error (MSE) on the validation
|
|
|
|
|
set, the higher, the less likely the model generalizes correctly from the
|
|
|
|
|
training data.
|
|
|
|
|
"""
|
|
|
|
|
|
|
|
|
|
print(__doc__)
|
|
|
|
|
|
|
|
|
|
import numpy as np
|
|
|
|
|
import matplotlib.pyplot as plt
|
|
|
|
|
from sklearn.pipeline import Pipeline
|
|
|
|
|
from sklearn.preprocessing import PolynomialFeatures
|
|
|
|
|
from sklearn.linear_model import LinearRegression
|
|
|
|
|
from sklearn.model_selection import cross_val_score
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
def true_fun(X):
|
|
|
|
|
return np.cos(1.5 * np.pi * X)
|
|
|
|
|
|
|
|
|
|
np.random.seed(0)
|
|
|
|
|
|
|
|
|
|
n_samples = 30
|
|
|
|
|
degrees = [1, 4, 15]
|
|
|
|
|
|
|
|
|
|
X = np.sort(np.random.rand(n_samples))
|
|
|
|
|
y = true_fun(X) + np.random.randn(n_samples) * 0.1
|
|
|
|
|
|
|
|
|
|
plt.figure(figsize=(14, 5))
|
|
|
|
|
for i in range(len(degrees)):
|
|
|
|
|
ax = plt.subplot(1, len(degrees), i + 1)
|
|
|
|
|
plt.setp(ax, xticks=(), yticks=())
|
|
|
|
|
|
|
|
|
|
polynomial_features = PolynomialFeatures(degree=degrees[i],
|
|
|
|
|
include_bias=False)
|
|
|
|
|
linear_regression = LinearRegression()
|
|
|
|
|
pipeline = Pipeline([("polynomial_features", polynomial_features),
|
|
|
|
|
("linear_regression", linear_regression)])
|
|
|
|
|
pipeline.fit(X[:, np.newaxis], y)
|
|
|
|
|
|
|
|
|
|
# Evaluate the models using crossvalidation
|
|
|
|
|
scores = cross_val_score(pipeline, X[:, np.newaxis], y,
|
|
|
|
|
scoring="neg_mean_squared_error", cv=10)
|
|
|
|
|
|
|
|
|
|
X_test = np.linspace(0, 1, 100)
|
|
|
|
|
plt.plot(X_test, pipeline.predict(X_test[:, np.newaxis]), label="Model")
|
|
|
|
|
plt.plot(X_test, true_fun(X_test), label="True function")
|
|
|
|
|
plt.scatter(X, y, edgecolor='b', s=20, label="Samples")
|
|
|
|
|
plt.xlabel("x")
|
|
|
|
|
plt.ylabel("y")
|
|
|
|
|
plt.xlim((0, 1))
|
|
|
|
|
plt.ylim((-2, 2))
|
|
|
|
|
plt.legend(loc="best")
|
|
|
|
|
plt.title("Degree {}\nMSE = {:.2e}(+/- {:.2e})".format(
|
|
|
|
|
degrees[i], -scores.mean(), scores.std()))
|
|
|
|
|
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
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
!split
|
|
|
|
|
===== The Ising model =====
|
|
|
|
|
|
|
|
|
|
The one-dimensional Ising model with nearest neighbor interaction, no
|
|
|
|
|
external field and a constant coupling constant $J$ is given by
|
|
|
|
|
|
|
|
|
|
!bt
|
|
|
|
|
\begin{align}
|
|
|
|
|
H = -J \sum_{k}^L s_k s_{k + 1},
|
|
|
|
|
\end{align}
|
|
|
|
|
!et
|
|
|
|
|
|
|
|
|
|
where $s_i \in \{-1, 1\}$ and $s_{N + 1} = s_1$. The number of spins
|
|
|
|
|
in the system is determined by $L$. For the one-dimensional system
|
|
|
|
|
there is no phase transition.
|
|
|
|
|
|
|
|
|
|
We will look at a system of $L = 40$ spins with a coupling constant of
|
|
|
|
|
$J = 1$. To get enough training data we will generate 10000 states
|
|
|
|
|
with their respective energies.
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
!bc pycod
|
|
|
|
|
import numpy as np
|
|
|
|
|
import matplotlib.pyplot as plt
|
|
|
|
|
from mpl_toolkits.axes_grid1 import make_axes_locatable
|
|
|
|
|
import seaborn as sns
|
|
|
|
|
import scipy.linalg as scl
|
|
|
|
|
from sklearn.model_selection import train_test_split
|
|
|
|
|
import tqdm
|
|
|
|
|
sns.set(color_codes=True)
|
|
|
|
|
cmap_args=dict(vmin=-1., vmax=1., cmap='seismic')
|
|
|
|
|
|
|
|
|
|
L = 40
|
|
|
|
|
n = int(1e4)
|
|
|
|
|
|
|
|
|
|
spins = np.random.choice([-1, 1], size=(n, L))
|
|
|
|
|
J = 1.0
|
|
|
|
|
|
|
|
|
|
energies = np.zeros(n)
|
|
|
|
|
|
|
|
|
|
for i in range(n):
|
|
|
|
|
energies[i] = - J * np.dot(spins[i], np.roll(spins[i], 1))
|
|
|
|
|
!ec
|
|
|
|
|
|
|
|
|
|
Here we use ordinary least squares
|
|
|
|
|
regression to predict the energy for the nearest neighbor
|
|
|
|
|
one-dimensional Ising model on a ring, i.e., the endpoints wrap
|
|
|
|
|
around. We will use linear regression to fit a value for
|
|
|
|
|
the coupling constant to achieve this.
|
|
|
|
|
|
|
|
|
|
!split
|
|
|
|
|
===== Reformulating the problem to suit regression =====
|
|
|
|
|
|
|
|
|
|
A more general form for the one-dimensional Ising model is
|
|
|
|
|
|
|
|
|
|
!bt
|
|
|
|
|
\begin{align}
|
|
|
|
|
H = - \sum_j^L \sum_k^L s_j s_k J_{jk}.
|
|
|
|
|
\end{align}
|
|
|
|
|
!et
|
|
|
|
|
|
|
|
|
|
Here we allow for interactions beyond the nearest neighbors and a state dependent
|
|
|
|
|
coupling constant. This latter expression can be formulated as
|
|
|
|
|
a matrix-product
|
|
|
|
|
!bt
|
|
|
|
|
\begin{align}
|
|
|
|
|
\bm{H} = \bm{X} J,
|
|
|
|
|
\end{align}
|
|
|
|
|
!et
|
|
|
|
|
|
|
|
|
|
where $X_{jk} = s_j s_k$ and $J$ is a matrix which consists of the
|
|
|
|
|
elements $-J_{jk}$. This form of writing the energy fits perfectly
|
|
|
|
|
with the form utilized in linear regression, that is
|
|
|
|
|
|
|
|
|
|
!bt
|
|
|
|
|
\begin{align}
|
|
|
|
|
\bm{y} = \bm{X}\bm{\beta} + \bm{\epsilon},
|
|
|
|
|
\end{align}
|
|
|
|
|
!et
|
|
|
|
|
|
|
|
|
|
We split the data in training and test data as discussed in the previous example
|
|
|
|
|
|
|
|
|
|
!bc pycod
|
|
|
|
|
X = np.zeros((n, L ** 2))
|
|
|
|
|
for i in range(n):
|
|
|
|
|
X[i] = np.outer(spins[i], spins[i]).ravel()
|
|
|
|
|
y = energies
|
|
|
|
|
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)
|
|
|
|
|
!ec
|
|
|
|
|
|
|
|
|
|
!split
|
|
|
|
|
===== Linear regression =====
|
|
|
|
|
|
|
|
|
|
In the ordinary least squares method we choose the cost function
|
|
|
|
|
|
|
|
|
|
!bt
|
|
|
|
|
\begin{align}
|
|
|
|
|
C(\bm{X}, \bm{\beta})= \frac{1}{n}\left\{(\bm{X}\bm{\beta} - \bm{y})^T(\bm{X}\bm{\beta} - \bm{y})\right\}.
|
|
|
|
|
\end{align}
|
|
|
|
|
!et
|
|
|
|
|
|
|
|
|
|
We then find the extremal point of $C$ by taking the derivative with respect to $\bm{\beta}$ as discussed above.
|
|
|
|
|
This yields the expression for $\bm{\beta}$ to be
|
|
|
|
|
|
|
|
|
|
!bt
|
|
|
|
|
\[
|
|
|
|
|
\bm{\beta} = \frac{\bm{X}^T \bm{y}}{\bm{X}^T \bm{X}},
|
|
|
|
|
\]
|
|
|
|
|
!et
|
|
|
|
|
|
|
|
|
|
which immediately imposes some requirements on $\bm{X}$ as there must exist
|
|
|
|
|
an inverse of $\bm{X}^T \bm{X}$. If the expression we are modeling contains an
|
|
|
|
|
intercept, i.e., a constant term, we must make sure that the
|
|
|
|
|
first column of $\bm{X}$ consists of $1$. We do this here
|
|
|
|
|
|
|
|
|
|
!bc pycod
|
|
|
|
|
X_train_own = np.concatenate(
|
|
|
|
|
(np.ones(len(X_train))[:, np.newaxis], X_train),
|
|
|
|
|
axis=1
|
|
|
|
|
)
|
|
|
|
|
X_test_own = np.concatenate(
|
|
|
|
|
(np.ones(len(X_test))[:, np.newaxis], X_test),
|
|
|
|
|
axis=1
|
|
|
|
|
)
|
|
|
|
|
!ec
|
|
|
|
|
|
|
|
|
|
!bc pycod
|
|
|
|
|
def ols_inv(x: np.ndarray, y: np.ndarray) -> np.ndarray:
|
|
|
|
|
return scl.inv(x.T @ x) @ (x.T @ y)
|
|
|
|
|
beta = ols_inv(X_train_own, y_train)
|
|
|
|
|
!ec
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
!split
|
|
|
|
|
===== Singular Value decomposition =====
|
|
|
|
|
|
|
|
|
|
Doing the inversion directly turns out to be a bad idea since the matrix
|
|
|
|
|
$\bm{X}^T\bm{X}$ is singular. An alternative approach is to use the _singular
|
|
|
|
|
value decomposition_. Using the definition of the Moore-Penrose
|
|
|
|
|
pseudoinverse we can write the equation for $\bm{\beta}$ as
|
|
|
|
|
|
|
|
|
|
!bt
|
|
|
|
|
\[
|
|
|
|
|
\bm{\beta} = \bm{X}^{+}\bm{y},
|
|
|
|
|
\]
|
|
|
|
|
!et
|
|
|
|
|
|
|
|
|
|
where the pseudoinverse of $\bm{X}$ is given by
|
|
|
|
|
|
|
|
|
|
!bt
|
|
|
|
|
\[
|
|
|
|
|
\bm{X}^{+} = \frac{\bm{X}^T}{\bm{X}^T\bm{X}}.
|
|
|
|
|
\]
|
|
|
|
|
!et
|
|
|
|
|
|
|
|
|
|
Using singular value decomposition we can decompose the matrix $\bm{X} = \bm{U}\bm{\Sigma} \bm{V}^T$,
|
|
|
|
|
where $\bm{U}$ and $\bm{V}$ are orthogonal(unitary) matrices and $\bm{\Sigma}$ contains the singular values (more details below).
|
|
|
|
|
where $X^{+} = V\Sigma^{+} U^T$. This reduces the equation for
|
|
|
|
|
$\omega$ to
|
|
|
|
|
!bt
|
|
|
|
|
\begin{align}
|
|
|
|
|
\bm{\beta} = \bm{V}\bm{\Sigma}^{+} \bm{U}^T \bm{y}.
|
|
|
|
|
\end{align}
|
|
|
|
|
!et
|
|
|
|
|
|
|
|
|
|
Note that solving this equation by actually doing the pseudoinverse
|
|
|
|
|
(which is what we will do) is not a good idea as this operation scales
|
|
|
|
|
as $\mathcal{O}(n^3)$, where $n$ is the number of elements in a
|
|
|
|
|
general matrix. Instead, doing $QR$-factorization and solving the
|
|
|
|
|
linear system as an equation would reduce this down to
|
|
|
|
|
$\mathcal{O}(n^2)$ operations.
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
!bc pycod
|
|
|
|
|
def ols_svd(x: np.ndarray, y: np.ndarray) -> np.ndarray:
|
|
|
|
|
u, s, v = scl.svd(x)
|
|
|
|
|
return v.T @ scl.pinv(scl.diagsvd(s, u.shape[0], v.shape[0])) @ u.T @ y
|
|
|
|
|
!ec
|
|
|
|
|
|
|
|
|
|
!bc pycod
|
|
|
|
|
beta = ols_svd(X_train_own,y_train)
|
|
|
|
|
!ec
|
|
|
|
|
|
|
|
|
|
When extracting the $J$-matrix we need to make sure that we remove the intercept, as is done here
|
|
|
|
|
|
|
|
|
|
!bc pycod
|
|
|
|
|
J = beta[1:].reshape(L, L)
|
|
|
|
|
!ec
|
|
|
|
|
|
|
|
|
|
A way of looking at the coefficients in $J$ is to plot the matrices as images.
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
!bc pycod
|
|
|
|
|
fig = plt.figure(figsize=(20, 14))
|
|
|
|
|
im = plt.imshow(J, **cmap_args)
|
|
|
|
|
plt.title("OLS", fontsize=18)
|
|
|
|
|
plt.xticks(fontsize=18)
|
|
|
|
|
plt.yticks(fontsize=18)
|
|
|
|
|
cb = fig.colorbar(im)
|
|
|
|
|
cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=18)
|
|
|
|
|
plt.show()
|
|
|
|
|
!ec
|
|
|
|
|
It is interesting to note that OLS
|
|
|
|
|
considers both $J_{j, j + 1} = -0.5$ and $J_{j, j - 1} = -0.5$ as
|
|
|
|
|
valid matrix elements for $J$.
|
|
|
|
|
In our discussion below on hyperparameters and Ridge and Lasso regression we will see that
|
|
|
|
|
this problem can be removed, partly and only with Lasso regression.
|
|
|
|
|
|
|
|
|
|
In this case our matrix inversion was actually possible. The obvious question now is what is the mathematics behind the SVD?
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
!split
|
|
|
|
|
===== The one-dimensional Ising model =====
|
|
|
|
|
|
|
|
|
|
Let us bring back the Ising model again, but now with an additional
|
|
|
|
|
focus on Ridge and Lasso regression as well. We repeat some of the
|
|
|
|
|
basic parts of the Ising model and the setup of the training and test
|
|
|
|
|
data. The one-dimensional Ising model with nearest neighbor
|
|
|
|
|
interaction, no external field and a constant coupling constant $J$ is
|
|
|
|
|
given by
|
|
|
|
|
|
|
|
|
|
!bt
|
|
|
|
|
\begin{align}
|
|
|
|
|
H = -J \sum_{k}^L s_k s_{k + 1},
|
|
|
|
|
\end{align}
|
|
|
|
|
!et
|
|
|
|
|
where $s_i \in \{-1, 1\}$ and $s_{N + 1} = s_1$. The number of spins in the system is determined by $L$. For the one-dimensional system there is no phase transition.
|
|
|
|
|
|
|
|
|
|
We will look at a system of $L = 40$ spins with a coupling constant of $J = 1$. To get enough training data we will generate 10000 states with their respective energies.
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
!bc pycod
|
|
|
|
|
import numpy as np
|
|
|
|
|
import matplotlib.pyplot as plt
|
|
|
|
|
from mpl_toolkits.axes_grid1 import make_axes_locatable
|
|
|
|
|
import seaborn as sns
|
|
|
|
|
import scipy.linalg as scl
|
|
|
|
|
from sklearn.model_selection import train_test_split
|
|
|
|
|
import sklearn.linear_model as skl
|
|
|
|
|
import tqdm
|
|
|
|
|
sns.set(color_codes=True)
|
|
|
|
|
cmap_args=dict(vmin=-1., vmax=1., cmap='seismic')
|
|
|
|
|
|
|
|
|
|
L = 40
|
|
|
|
|
n = int(1e4)
|
|
|
|
|
|
|
|
|
|
spins = np.random.choice([-1, 1], size=(n, L))
|
|
|
|
|
J = 1.0
|
|
|
|
|
|
|
|
|
|
energies = np.zeros(n)
|
|
|
|
|
|
|
|
|
|
for i in range(n):
|
|
|
|
|
energies[i] = - J * np.dot(spins[i], np.roll(spins[i], 1))
|
|
|
|
|
!ec
|
|
|
|
|
|
|
|
|
|
A more general form for the one-dimensional Ising model is
|
|
|
|
|
|
|
|
|
|
!bt
|
|
|
|
|
\begin{align}
|
|
|
|
|
H = - \sum_j^L \sum_k^L s_j s_k J_{jk}.
|
|
|
|
|
\end{align}
|
|
|
|
|
!et
|
|
|
|
|
|
|
|
|
|
Here we allow for interactions beyond the nearest neighbors and a more
|
|
|
|
|
adaptive coupling matrix. This latter expression can be formulated as
|
|
|
|
|
a matrix-product on the form
|
|
|
|
|
!bt
|
|
|
|
|
\begin{align}
|
|
|
|
|
H = X J,
|
|
|
|
|
\end{align}
|
|
|
|
|
!et
|
|
|
|
|
|
|
|
|
|
where $X_{jk} = s_j s_k$ and $J$ is the matrix consisting of the
|
|
|
|
|
elements $-J_{jk}$. This form of writing the energy fits perfectly
|
|
|
|
|
with the form utilized in linear regression, viz.
|
|
|
|
|
!bt
|
|
|
|
|
\begin{align}
|
|
|
|
|
\bm{y} = \bm{X}\bm{\beta} + \bm{\epsilon}.
|
|
|
|
|
\end{align}
|
|
|
|
|
!et
|
|
|
|
|
We organize the data as we did above
|
|
|
|
|
!bc pycod
|
|
|
|
|
X = np.zeros((n, L ** 2))
|
|
|
|
|
for i in range(n):
|
|
|
|
|
X[i] = np.outer(spins[i], spins[i]).ravel()
|
|
|
|
|
y = energies
|
|
|
|
|
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.96)
|
|
|
|
|
|
|
|
|
|
X_train_own = np.concatenate(
|
|
|
|
|
(np.ones(len(X_train))[:, np.newaxis], X_train),
|
|
|
|
|
axis=1
|
|
|
|
|
)
|
|
|
|
|
|
|
|
|
|
X_test_own = np.concatenate(
|
|
|
|
|
(np.ones(len(X_test))[:, np.newaxis], X_test),
|
|
|
|
|
axis=1
|
|
|
|
|
)
|
|
|
|
|
!ec
|
|
|
|
|
|
|
|
|
|
We will do all fitting with _Scikit-Learn_,
|
|
|
|
|
|
|
|
|
|
!bc pycod
|
|
|
|
|
clf = skl.LinearRegression().fit(X_train, y_train)
|
|
|
|
|
!ec
|
|
|
|
|
When extracting the $J$-matrix we make sure to remove the intercept
|
|
|
|
|
!bc pycod
|
|
|
|
|
J_sk = clf.coef_.reshape(L, L)
|
|
|
|
|
!ec
|
|
|
|
|
And then we plot the results
|
|
|
|
|
!bc pycod
|
|
|
|
|
fig = plt.figure(figsize=(20, 14))
|
|
|
|
|
im = plt.imshow(J_sk, **cmap_args)
|
|
|
|
|
plt.title("LinearRegression from Scikit-learn", fontsize=18)
|
|
|
|
|
plt.xticks(fontsize=18)
|
|
|
|
|
plt.yticks(fontsize=18)
|
|
|
|
|
cb = fig.colorbar(im)
|
|
|
|
|
cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=18)
|
|
|
|
|
plt.show()
|
|
|
|
|
!ec
|
|
|
|
|
The results perfectly with our previous discussion where we used our own code.
|
|
|
|
|
|
|
|
|
|
!split
|
|
|
|
|
===== Ridge regression =====
|
|
|
|
|
|
|
|
|
|
Having explored the ordinary least squares we move on to ridge
|
|
|
|
|
regression. In ridge regression we include a _regularizer_. This
|
|
|
|
|
involves a new cost function which leads to a new estimate for the
|
|
|
|
|
weights $\bm{\beta}$. This results in a penalized regression problem. The
|
|
|
|
|
cost function is given by
|
|
|
|
|
|
|
|
|
|
!bt
|
|
|
|
|
\begin{align}
|
|
|
|
|
C(\bm{X}, \bm{\beta}; \lambda) = (\bm{X}\bm{\beta} - \bm{y})^T(\bm{X}\bm{\beta} - \bm{y}) + \lambda \bm{\beta}^T\bm{\beta}.
|
|
|
|
|
\end{align}
|
|
|
|
|
!et
|
|
|
|
|
!bc pycod
|
|
|
|
|
_lambda = 0.1
|
|
|
|
|
clf_ridge = skl.Ridge(alpha=_lambda).fit(X_train, y_train)
|
|
|
|
|
J_ridge_sk = clf_ridge.coef_.reshape(L, L)
|
|
|
|
|
fig = plt.figure(figsize=(20, 14))
|
|
|
|
|
im = plt.imshow(J_ridge_sk, **cmap_args)
|
|
|
|
|
plt.title("Ridge from Scikit-learn", fontsize=18)
|
|
|
|
|
plt.xticks(fontsize=18)
|
|
|
|
|
plt.yticks(fontsize=18)
|
|
|
|
|
cb = fig.colorbar(im)
|
|
|
|
|
cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=18)
|
|
|
|
|
|
|
|
|
|
plt.show()
|
|
|
|
|
!ec
|
|
|
|
|
|
|
|
|
|
!split
|
|
|
|
|
===== LASSO regression =====
|
|
|
|
|
|
|
|
|
|
In the _Least Absolute Shrinkage and Selection Operator_ (LASSO)-method we get a third cost function.
|
|
|
|
|
|
|
|
|
|
!bt
|
|
|
|
|
\begin{align}
|
|
|
|
|
C(\bm{X}, \bm{\beta}; \lambda) = (\bm{X}\bm{\beta} - \bm{y})^T(\bm{X}\bm{\beta} - \bm{y}) + \lambda \sqrt{\bm{\beta}^T\bm{\beta}}.
|
|
|
|
|
\end{align}
|
|
|
|
|
!et
|
|
|
|
|
|
|
|
|
|
Finding the extremal point of this cost function is not so straight-forward as in least squares and ridge. We will therefore rely solely on the function ``Lasso`` from _Scikit-Learn_.
|
|
|
|
|
|
|
|
|
|
!bc pycod
|
|
|
|
|
clf_lasso = skl.Lasso(alpha=_lambda).fit(X_train, y_train)
|
|
|
|
|
J_lasso_sk = clf_lasso.coef_.reshape(L, L)
|
|
|
|
|
fig = plt.figure(figsize=(20, 14))
|
|
|
|
|
im = plt.imshow(J_lasso_sk, **cmap_args)
|
|
|
|
|
plt.title("Lasso from Scikit-learn", fontsize=18)
|
|
|
|
|
plt.xticks(fontsize=18)
|
|
|
|
|
plt.yticks(fontsize=18)
|
|
|
|
|
cb = fig.colorbar(im)
|
|
|
|
|
cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=18)
|
|
|
|
|
|
|
|
|
|
plt.show()
|
|
|
|
|
!ec
|
|
|
|
|
|
|
|
|
|
It is quite striking how LASSO breaks the symmetry of the coupling
|
|
|
|
|
constant as opposed to ridge and OLS. We get a sparse solution with
|
|
|
|
|
$J_{j, j + 1} = -1$.
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
!split
|
|
|
|
|
===== Performance as function of the regularization parameter =====
|
|
|
|
|
|
|
|
|
|
We see how the different models perform for a different set of values for $\lambda$.
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
!bc pycod
|
|
|
|
|
lambdas = np.logspace(-4, 5, 10)
|
|
|
|
|
|
|
|
|
|
train_errors = {
|
|
|
|
|
"ols_sk": np.zeros(lambdas.size),
|
|
|
|
|
"ridge_sk": np.zeros(lambdas.size),
|
|
|
|
|
"lasso_sk": np.zeros(lambdas.size)
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
test_errors = {
|
|
|
|
|
"ols_sk": np.zeros(lambdas.size),
|
|
|
|
|
"ridge_sk": np.zeros(lambdas.size),
|
|
|
|
|
"lasso_sk": np.zeros(lambdas.size)
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
plot_counter = 1
|
|
|
|
|
|
|
|
|
|
fig = plt.figure(figsize=(32, 54))
|
|
|
|
|
|
|
|
|
|
for i, _lambda in enumerate(tqdm.tqdm(lambdas)):
|
|
|
|
|
for key, method in zip(
|
|
|
|
|
["ols_sk", "ridge_sk", "lasso_sk"],
|
|
|
|
|
[skl.LinearRegression(), skl.Ridge(alpha=_lambda), skl.Lasso(alpha=_lambda)]
|
|
|
|
|
):
|
|
|
|
|
method = method.fit(X_train, y_train)
|
|
|
|
|
|
|
|
|
|
train_errors[key][i] = method.score(X_train, y_train)
|
|
|
|
|
test_errors[key][i] = method.score(X_test, y_test)
|
|
|
|
|
|
|
|
|
|
omega = method.coef_.reshape(L, L)
|
|
|
|
|
|
|
|
|
|
plt.subplot(10, 5, plot_counter)
|
|
|
|
|
plt.imshow(omega, **cmap_args)
|
|
|
|
|
plt.title(r"%s, $\lambda = %.4f$" % (key, _lambda))
|
|
|
|
|
plot_counter += 1
|
|
|
|
|
|
|
|
|
|
plt.show()
|
|
|
|
|
!ec
|
|
|
|
|
|
|
|
|
|
We see that LASSO reaches a good solution for low
|
|
|
|
|
values of $\lambda$, but will "wither" when we increase $\lambda$ too
|
|
|
|
|
much. Ridge is more stable over a larger range of values for
|
|
|
|
|
$\lambda$, but eventually also fades away.
|
|
|
|
|
|
|
|
|
|
!split
|
|
|
|
|
===== Finding the optimal value of $\lambda$ =====
|
|
|
|
|
|
|
|
|
|
To determine which value of $\lambda$ is best we plot the accuracy of
|
|
|
|
|
the models when predicting the training and the testing set. We expect
|
|
|
|
|
the accuracy of the training set to be quite good, but if the accuracy
|
|
|
|
|
of the testing set is much lower this tells us that we might be
|
|
|
|
|
subject to an overfit model. The ideal scenario is an accuracy on the
|
|
|
|
|
testing set that is close to the accuracy of the training set.
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
!bc pycod
|
|
|
|
|
fig = plt.figure(figsize=(20, 14))
|
|
|
|
|
|
|
|
|
|
colors = {
|
|
|
|
|
"ols_sk": "r",
|
|
|
|
|
"ridge_sk": "y",
|
|
|
|
|
"lasso_sk": "c"
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
for key in train_errors:
|
|
|
|
|
plt.semilogx(
|
|
|
|
|
lambdas,
|
|
|
|
|
train_errors[key],
|
|
|
|
|
colors[key],
|
|
|
|
|
label="Train {0}".format(key),
|
|
|
|
|
linewidth=4.0
|
|
|
|
|
)
|
|
|
|
|
|
|
|
|
|
for key in test_errors:
|
|
|
|
|
plt.semilogx(
|
|
|
|
|
lambdas,
|
|
|
|
|
test_errors[key],
|
|
|
|
|
colors[key] + "--",
|
|
|
|
|
label="Test {0}".format(key),
|
|
|
|
|
linewidth=4.0
|
|
|
|
|
)
|
|
|
|
|
plt.legend(loc="best", fontsize=18)
|
|
|
|
|
plt.xlabel(r"$\lambda$", fontsize=18)
|
|
|
|
|
plt.ylabel(r"$R^2$", fontsize=18)
|
|
|
|
|
plt.tick_params(labelsize=18)
|
|
|
|
|
plt.show()
|
|
|
|
|
!ec
|
|
|
|
|
|
|
|
|
|
From the above figure we can see that LASSO with $\lambda = 10^{-2}$
|
|
|
|
|
achieves a very good accuracy on the test set. This by far surpasses the
|
|
|
|
|
other models for all values of $\lambda$.
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
!split
|
|
|
|
|
===== Friday September 18: Intro to Logistic Regression =====
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
!split
|
|
|
|
@@ -550,3 +1451,4 @@ plt.show()
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Morten Hjorth-Jensen