2
$\begingroup$

I am trying to reproduce this figure from Bishop:

enter image description here

Residual vs. Alpha (lambda in figure)

The code is pasted below:

import numpy as np
import matplotlib.pyplot as plt
from sklearn import linear_model

from sklearn.linear_model import Ridge
from sklearn.preprocessing import PolynomialFeatures
from sklearn.pipeline import make_pipeline
from sklearn.metrics import mean_squared_error

x_train = np.linspace(0, 1, 10)
noise_10 = np.random.normal(0, 0.3, 10)
y_train = np.sin(2*np.pi*x_train) + noise_10
x_train = x_train[:, np.newaxis]

x_test = np.linspace(0, 1, 100)
noise_100 = np.random.normal(0, 0.3, 100)
y_test = np.sin(2*np.pi*x_test) + noise_100
x_test = x_test[:, np.newaxis]


n_alphas = 200
alphas = np.logspace(-40, -18, n_alphas)

errors = []
for a in alphas:
    ridge = make_pipeline(PolynomialFeatures(degree = 9), 
              Ridge(alpha=a))
    ridge.fit(x_train, y_train)
    mse = mean_squared_error(y_test, ridge.predict(x_test))
    errors.append(np.sqrt(mse))    

print(errors)

However, the errors array has the same value for all values of alpha. It's taking the first value of alpha = np.exp(-40) and all the other values seem to be the same for all future iterations of the for loop. How can I correct this error?

$\endgroup$
1
$\begingroup$

I still haven't figured out what the previously posted code is doing wrong. However, manually populating the alphas array gives me results that are close to the original figure.

import numpy as np
import matplotlib.pyplot as plt 
from sklearn import linear_model
from sklearn.linear_model import Ridge
from sklearn.preprocessing import PolynomialFeatures 
from sklearn.pipeline import make_pipeline
from sklearn.metrics import mean_squared_error
np.random.seed(12344)
x_1000 = np.linspace(0, 1, 1000)
x_train = np.linspace(0, 1, 10)
noise_10 = np.random.normal(0, 0.3, 10)
y_train = np.sin(2*np.pi*x_train) + noise_10
x_train = x_train[:, np.newaxis]
x_test = np.linspace(0, 1, 100)
noise_100 = np.random.normal(0, 0.3, 100)
y_test = np.sin(2*np.pi*x_test) + noise_100
x_test = x_test[:, np.newaxis]
low_alpha = make_pipeline(PolynomialFeatures(degree = 9), Ridge(alpha=np.exp(-18)))
low_alpha.fit(x_test, y_test)
plt.figure(1)
plt.plot(x_train, low_alpha.predict(x_train), label = 'alpha = ln(-18)')
plt.plot(x_1000, np.sin(2*np.pi*x_1000))
plt.legend()
plt.show
high_alpha = make_pipeline(PolynomialFeatures(degree = 9),
              Ridge(alpha=np.exp(0)))
high_alpha.fit(x_test, y_test)
plt.figure(2)
plt.plot(x_train, high_alpha.predict(x_train), label = 'alpha = 1')
plt.plot(x_1000, np.sin(2*np.pi*x_1000))
plt.legend()
plt.show
alphas = np.array([np.exp(-30), np.exp(-29), np.exp(-28), np.exp(-27), np.exp(-26), np.exp(-25), np.exp(-24), np.exp(-23), np.exp(-22), np.exp(-21), np.exp(-20), np.exp(-19), np.exp(-18), np.exp(-17), np.exp(-16), np.exp(-15), np.exp(-14), np.exp(-13), np.exp(-12), np.exp(-11), np.exp(-10), np.exp(-9), np.exp(-8), np.exp(-7), np.exp(-6), np.exp(-5), np.exp(-4), np.exp(-3), np.exp(-2), np.exp(-1), np.exp(-1), np.exp(0), np.exp(1)])

test_errors = []
train_errors = []

for a in np.nditer(alphas):
    ridge = make_pipeline(PolynomialFeatures(degree = 9),
              Ridge(alpha=a))
    ridge.fit(x_train, y_train)
    mse_train = mean_squared_error(y_train, ridge.predict(x_train))
    mse_test = mean_squared_error(y_test, ridge.predict(x_test))
    train_errors.append(np.sqrt(mse_train))
    test_errors.append(np.sqrt(mse_test))
plt.figure(3)
plt.plot(alphas, test_errors, 'g^', label = 'Test Error')
plt.plot(alphas, train_errors, 'bs', label = 'Train Error')
plt.xscale('log')
plt.xlabel('Regression coefficient Lambda')
plt.ylabel('Residuals')
plt.legend()
plt.show()

The output (only for the third figure) is shown below: Reproduction of Figure from Bishop

| improve this answer | |
$\endgroup$
0
$\begingroup$

I have taken a look on your code. You obtain same errors results for each alpha value because your regularization strength is too small. Replacing :

alphas = np.logspace(-40, -18, n_alphas)

with :

alphas = np.logspace(-40, -1, n_alphas)

will yields different errors values for alpha values large enough. Are you sure about figure alpha values? Do you have a link to this hands-on?

Also, I would like to highlight the fact that you have to standardize your features before using regularization. Reason is, by creating polynomial features, polynomial features will have different magnitudes. Therefore when fitting model, coefficients to be estimated won't have the same magnitudes neither and so regularization will highly penalize coefficients with large values. Standardization / Normalization is a strong prerequisite to regularization.

| improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.