Part 1 to 7 of exercise5. Function minimization is not fully working.

ex5/ex5.py

@@ -0,0 +1,255 @@
+## Machine Learning Online Class Exercise 5 | Regularized Linear Regression and Bias-Variance
+import numpy as np
+import matplotlib.pyplot as plt
+import scipy, scipy.io, scipy.optimize
+
+
+# LINEARREGCOSTFUNCTION Compute cost and gradient for regularized linear
+# regression with multiple variables
+# [J, grad] = LINEARREGCOSTFUNCTION(X, y, theta, lambda) computes the
+# cost of using theta as the parameter for linear regression to fit the
+# data points in X and y. Returns the cost in J and the gradient in grad
+def linear_reg_cost_function(X, y, theta, ld):
+ # Initialize some useful values
+ m = y.shape[0] # number of training examples
+ d1 = X.shape[1]
+ d2 = y.shape[1]
+
+ theta = theta.reshape((d1, d2))
+
+ J = 0.5 / m * (np.sum((X.dot(theta) - y) ** 2) + ld * np.sum(theta[1:, :] ** 2))
+ grad = 1 / m * X.T.dot((X.dot(theta) - y)) + ld / m * np.vstack((np.zeros((1, theta.shape[1])), theta[1:, :]))
+ return J, grad.flatten()
+
+
+# TRAINLINEARREG Trains linear regression given a dataset (X, y) and a
+# regularization parameter lambda
+# [theta] = TRAINLINEARREG (X, y, lambda) trains linear regression using
+# the dataset (X, y) and regularization parameter lambda. Returns the
+# trained parameters theta.
+def train_linear_reg(X, y, ld):
+ # Initialize Theta
+ initial_theta = np.zeros((X.shape[1], 1))
+
+ # Create "short hand" for the cost function to be minimized
+ cost_function = lambda t: linear_reg_cost_function(X, y, t, ld)
+
+ # Now, costFunction is a function that takes in only one argument
+ options = {'maxiter': 2000, 'norm': 3, 'disp': True}
+
+ # Minimize using minimize
+ # Because the implementation of Scipy, the result looks somehow different
+ theta = scipy.optimize.minimize(cost_function, initial_theta, options=options, jac=True, method='BFGS')['x']
+ return theta
+
+
+# LEARNINGCURVE Generates the train and cross validation set errors needed
+# to plot a learning curve
+# [error_train, error_val] = ...
+# LEARNINGCURVE(X, y, Xval, yval, lambda) returns the train and
+# cross validation set errors for a learning curve. In particular,
+# it returns two vectors of the same length - error_train and
+# error_val. Then, error_train(i) contains the training error for
+# i examples (and similarly for error_val(i)).
+#
+# In this function, you will compute the train and test errors for
+# dataset sizes from 1 up to m. In practice, when working with larger
+# datasets, you might want to do this in larger intervals.
+def learning_curve(X, y, Xval, yval, ld):
+ # Number of training examples
+ m = X.shape[0]
+
+ # You need to return these values correctly
+ error_train = np.zeros((m, 1))
+ error_val = np.zeros((m, 1))
+
+ for i in range(1, m + 1):
+ print(i)
+ theta = train_linear_reg(X[:i, :], y[:i, :], ld).reshape(-1, 1)
+ error_train[i - 1, 0] = 0.5 / i * np.sum((X[:i, :].dot(theta) - y[:i, :]) ** 2)
+ error_val[i - 1, 0] = 0.5 / yval.shape[0] * np.sum((Xval.dot(theta) - yval) ** 2)
+ return error_train, error_val
+
+
+# POLYFEATURES Maps X (1D vector) into the p-th power
+# [X_poly] = POLYFEATURES(X, p) takes a data matrix X (size m x 1) and
+# maps each example into its polynomial features where
+# X_poly(i, :) = [X(i) X(i).^2 X(i).^3 ... X(i).^p];
+def poly_features(X, p):
+ X_poly = X ** np.array(range(1, p + 1))
+ return X_poly
+
+
+# FEATURENORMALIZE Normalizes the features in X
+# FEATURENORMALIZE(X) returns a normalized version of X where
+# the mean value of each feature is 0 and the standard deviation
+# is 1. This is often a good preprocessing step to do when
+# working with learning algorithms.
+def feature_normalize(X):
+ mu = np.mean(X, axis=0)
+ X_norm = X - mu
+
+ sigma = np.std(X_norm, axis=0, ddof=1)
+ X_norm = X_norm / sigma
+ return X_norm, mu, sigma
+
+
+# PLOTFIT Plots a learned polynomial regression fit over an existing figure.
+# Also works with linear regression.
+# PLOTFIT(min_x, max_x, mu, sigma, theta, p) plots the learned polynomial
+# fit with power p and feature normalization (mu, sigma).
+def plotfit(min_x, max_x, mu, sigma, theta, p):
+ # We plot a range slightly bigger than the min and max values to get
+ # an idea of how the fit will vary outside the range of the data points
+ nnum = np.floor((max_x - min_x + 40) / 0.05)
+ x = np.linspace(min_x - 15, nnum * 0.05 + min_x - 15, nnum).reshape(-1, 1)
+
+ # Map the X values
+ X_poly = poly_features(x, p)
+ X_poly = X_poly - mu
+ X_poly = X_poly / sigma
+
+ # Add ones
+ X_poly = np.hstack((np.ones((x.shape[0], 1)), X_poly))
+
+ # Plot
+ plt.plot(x, X_poly.dot(theta), '--', linewidth=2)
+
+
+np.set_printoptions(formatter={'float': '{: 0.5f}'.format}, edgeitems=50, linewidth=150)
+## =========== Part 1: Loading and Visualizing Data =============
+# Load Training Data
+print('Loading and Visualizing Data ...\n')
+
+# Load from ex5data1:
+# You will have X, y, Xval, yval, Xtest, ytest in your environment
+data = scipy.io.loadmat('mat/ex5data1.mat', matlab_compatible=True)
+X = data['X']
+y = data['y']
+Xval = data['Xval']
+yval = data['yval']
+Xtest = data['Xtest']
+ytest = data['ytest']
+
+# m = Number of examples
+m = X.shape[0]
+
+# Plot training data
+plt.plot(X, y, 'rx', markersize=10, linewidth=1.5)
+plt.xlabel('Change in water level (x)')
+plt.ylabel('Water flowing out of the dam (y)')
+plt.show()
+
+print('Program paused. Press enter to continue.\n')
+input()
+
+## =========== Part 2: Regularized Linear Regression Cost =============
+theta = np.array([[1], [1]])
+J = linear_reg_cost_function(np.hstack((np.ones((m, 1)), X)), y, theta, 1)[0]
+
+print('Cost at theta = [1 ; 1]: {0:f} \n(this value should be about 303.993192)\n'.format(J))
+
+print('Program paused. Press enter to continue.\n')
+
+## =========== Part 3: Regularized Linear Regression Gradient =============
+J, grad = linear_reg_cost_function(np.hstack((np.ones((m, 1)), X)), y, theta, 1)
+
+print('Gradient at theta = [1 ; 1]: [{0:f}; {1:f}] \n(this value should be about [-15.303016; 598.250744])\n'.format(
+ grad[0], grad[1]))
+
+print('Program paused. Press enter to continue.\n')
+
+## =========== Part 4: Train Linear Regression =============
+# Train linear regression with ld = 0
+ld = 0
+
+theta = train_linear_reg(np.hstack((np.ones((m, 1)), X)), y, ld)
+
+# Plot fit over the data
+plt.plot(X, y, 'rx', markersize=10, linewidth=1.5)
+plt.xlabel('Change in water level (x)')
+plt.ylabel('Water flowing out of the dam (y)')
+plt.plot(X, np.hstack((np.ones((m, 1)), X)).dot(theta), '--', linewidth=2)
+plt.show()
+
+print('Program paused. Press enter to continue.\n')
+input()
+
+## =========== Part 5: Learning Curve for Linear Regression =============
+ld = 0
+
+error_train, error_val = learning_curve(np.hstack((np.ones((m, 1)), X)), y,
+ np.hstack((np.ones((Xval.shape[0], 1)), Xval)), yval, ld)
+
+plt.plot(np.array(range(1, m + 1)), error_train, np.array(range(1, m + 1)), error_val)
+plt.title('Learning curve for linear regression')
+plt.legend(['Train', 'Cross Validation'])
+plt.xlabel('Number of training examples')
+plt.ylabel('Error')
+plt.axis([0, 13, 0, 150])
+plt.show()
+
+print('# Training Examples\tTrain Error\tCross Validation Error\n')
+for i in range(0, m):
+ print(' \t{0:d}\t\t{1:f}\t{2:f}'.format(i, error_train[i, 0], error_val[i, 0]))
+
+print('Program paused. Press enter to continue.\n')
+input()
+
+## =========== Part 6: Feature Mapping for Polynomial Regression =============
+p = 8
+
+# Map X onto Polynomial Features and Normalize
+X_poly = poly_features(X, p)
+X_poly2 = np.copy(X_poly)
+X_poly, mu, sigma = feature_normalize(X_poly) # Normalize
+X_poly = np.hstack((np.ones((m, 1)), X_poly)) # Add Ones
+
+# Map X_poly_test and normalize (using mu and sigma)
+X_poly_test = poly_features(Xtest, p)
+X_poly_test = X_poly_test - mu
+X_poly_test = X_poly_test / sigma
+X_poly_test = np.hstack((np.ones((X_poly_test.shape[0], 1)), X_poly_test)) # Add Ones
+
+# Map X_poly_val and normalize (using mu and sigma)
+X_poly_val = poly_features(Xval, p)
+X_poly_val2 = np.copy(X_poly_val)
+X_poly_val = X_poly_val - mu
+X_poly_val = X_poly_val / sigma
+X_poly_val = np.hstack((np.ones((X_poly_val.shape[0], 1)), X_poly_val)) # Add Ones
+
+print('Normalized Training Example 1:\n')
+print(X_poly[0, :])
+
+print('\nProgram paused. Press enter to continue.\n')
+input()
+
+## =========== Part 7: Learning Curve for Polynomial Regression =============
+ld = 0
+theta = train_linear_reg(X_poly, y, ld)
+
+# Plot training data and fit
+plt.figure()
+plt.plot(X, y, 'rx', markersize=10, linewidth=1.5)
+plotfit(np.min(X), np.max(X), mu, sigma, theta, p)
+plt.xlabel('Change in water level (x)')
+plt.ylabel('Water flowing out of the dam (y)')
+plt.title('Polynomial Regression Fit (lambda = {0:f})'.format(ld))
+
+plt.figure()
+error_train, error_val = learning_curve(X_poly, y, X_poly_val, yval, ld)
+plt.plot(np.array(range(1, m + 1)), error_train, np.array(range(1, m + 1)), error_val)
+
+plt.title('Polynomial Regression Learning Curve (lambda = {0:f})'.format(ld))
+plt.xlabel('Number of training examples')
+plt.ylabel('Error')
+plt.axis([0, 13, 0, 100])
+plt.legend(['Train', 'Cross Validation'])
+plt.show()
+
+print('Polynomial Regression (lambda = {0:f})\n\n'.format(ld))
+print('# Training Examples\tTrain Error\tCross Validation Error\n')
+for i in range(0, m):
+ print(' \t{0:d}\t\t{1:f}\t{2:f}\n'.format(i + 1, error_train[i, 0], error_val[i, 0]))
+
+print('Program paused. Press enter to continue.\n')