Finished exercise2.

ex2/ex2.py

@@ -5,7 +5,6 @@ import scipy.optimize
 import math
 
 
-
 # PLOTDATA Plots the data points X and y into a new figure
 # PLOTDATA(x,y) plots the data points with + for the positive examples
 # and o for the negative examples. X is assumed to be a Mx2 matrix.
@@ -16,6 +15,7 @@ def plotdata(X, y):
  plt.plot(X[:, :1][pos], X[:, 1:2][pos], 'k+', linewidth=2, markersize=7)
  plt.plot(X[:, :1][neg], X[:, 1:2][neg], 'ko', markerfacecolor='y', markersize=7)
 
+
 # SIGMOID Compute sigmoid functoon
 # J = SIGMOID(z) computes the sigmoid of z.
 def sigmoid(z):
@@ -36,7 +36,8 @@ def cost_function(theta, X, y):
  h = sigmoid(X.dot(theta))
  J = np.sum(-y * np.log(h) - (1 - y) * np.log(1 - h)) / m
  grad = (h - y).T.dot(X) / m
- return J.flatten()[0], grad.flatten()
+ return J, grad.flatten()
+
 
 # MAPFEATURE Feature mapping function to polynomial features
 #
@@ -49,12 +50,11 @@ def cost_function(theta, X, y):
 # Inputs X1, X2 must be the same size
 def map_feature(X1, X2):
  degree = 6
- out = np.ones(X1[:,0].size)
- for i in range(1, degree+1):
- for j in range(0, i):
- pass
- # TODO out[:, end+1] = (X1 ** (i-j)).*(X2.^j)
- pass
+ out = np.ones(X1[:, :1].shape)
+ for i in range(1, degree + 1):
+ for j in range(0, i + 1):
+ out = np.hstack((out, (X1 ** (i - j)) * (X2 ** j)))
+ return out
 
 
 # PLOTDECISIONBOUNDARY Plots the data points X and y into a new figure with
@@ -67,13 +67,13 @@ def map_feature(X1, X2):
 # 2) MxN, N>3 matrix, where the first column is all-ones
 def plot_decision_boundary(theta, X, y):
  # Plot Data
- plotdata(X[:,2:4], y)
+ plotdata(X[:, 1:3], y)
  if theta.ndim == 2 and (theta.shape[0] == 1 or theta.shape[1] == 1):
  theta = theta.flatten()
 
  if X.shape[1] <= 3:
  # Only need 2 points to define a line, so choose two endpoints
- plot_x = np.array([np.min(X[:,2])[0]-2,  np.max(X[:,2])[0]+2])
+ plot_x = np.array([np.min(X[:, 1]) - 2, np.max(X[:, 1]) + 2])
 
  # Calculate the decision boundary line
  plot_y = (-1 / theta[2]) * (theta[1] * plot_x + theta[0])
@@ -86,88 +86,118 @@ def plot_decision_boundary(theta, X, y):
  plt.axis([30, 100, 30, 100])
  else:
  # Here is the grid range
- u = np.linspace(-1, 1.5, 50)
- v = np.linspace(-1, 1.5, 50)
+ u = np.linspace(-1, 1.5, 50).reshape(-1, 1)
+ v = np.linspace(-1, 1.5, 50).reshape(-1, 1)
 
- z = np.zeros(u.size, v.size)
+ z = np.zeros((u.size, v.size))
  # Evaluate z = theta*x over the grid
  for i in range(0, u.size):
  for j in range(0, v.size):
- z[i,j] = map_feature(u[i], v[j]) * theta
+ z[i, j] = map_feature(u[i:i + 1, :1], v[j:j + 1, :1]).dot(theta)
+
+ z = z.T # important to transpose z before calling contour
+
+ # Plot z = 0
+ # Notice you need to specify the range [0, 0]
+ plt.contour(u.flatten(), v.flatten(), z, [0], linewidth=2)
+
+
+# PREDICT Predict whether the label is 0 or 1 using learned logistic
+# regression parameters theta
+# p = PREDICT(theta, X) computes the predictions for X using a
+# threshold at 0.5 (i.e., if sigmoid(theta'*x) >= 0.5, predict 1)
+def predict(theta, X):
+ m = X.shape[0]
+ p = sigmoid(X.dot(theta)) >= 0.5
+ return p
+
+
+## ===================== Start Main ======================
+if __name__ == "__main__":
+ ## Load Data
+ # The first two columns contains the exam scores and the third column
+ # contains the label.
+ data = np.loadtxt('textdata/ex2data1.txt', delimiter=',')
+ X = data[:, :2]
+ y = data[:, 2:3]
+
+ ## ==================== Part 1: Plotting ====================
+ print('Plotting data with + indicating (y = 1) examples and o indicating (y = 0) examples.')
 
- z = z.T # important to transpose z before calling contour
+ plotdata(X, y)
 
- # Plot z = 0
- # Notice you need to specify the range [0, 0]
- plt.contour(u, v, z, [0, 0], 'LineWidth', 2)
- pass
+ # Put some labels
+ # Labels and Legend
+ plt.xlabel('Exam 1 score')
+ plt.ylabel('Exam 2 score')
 
+ # Specified in plot order
+ plt.legend(['Admitted', 'Not admitted'])
+ plt.show()
 
-## Load Data
-# The first two columns contains the exam scores and the third column
-# contains the label.
-data = np.loadtxt('textdata/ex2data1.txt', delimiter=',')
-X = data[:, :2]
-y = data[:, 2:3]
+ print('\nProgram paused. Press enter to continue.\n')
+ input()
 
-## ==================== Part 1: Plotting ====================
-print('Plotting data with + indicating (y = 1) examples and o indicating (y = 0) examples.')
+ ## ============ Part 2: Compute Cost and Gradient ============
+ # Setup the data matrix appropriately, and add ones for the intercept term
+ m, n = X.shape
 
-plotdata(X, y)
+ # Add intercept term to x and X_test
+ X = np.hstack((np.ones((m, 1)), X))
 
-# Put some labels
-# Labels and Legend
-plt.xlabel('Exam 1 score')
-plt.ylabel('Exam 2 score')
+ # Initialize fitting parameters
+ initial_theta = np.zeros((n + 1, 1))
 
-# Specified in plot order
-plt.legend(['Admitted', 'Not admitted'])
-# plt.show()
+ # Compute and display initial cost and gradient
+ cost, grad = cost_function(initial_theta, X, y)
 
-print('\nProgram paused. Press enter to continue.\n')
-input()
+ print('Cost at initial theta (zeros): {0:f}'.format(cost))
+ print('Gradient at initial theta (zeros): ')
+ print(np.array_str(grad.reshape(-1, 1)).replace('[', ' ').replace(']', ' '))
 
-## ============ Part 2: Compute Cost and Gradient ============
-# Setup the data matrix appropriately, and add ones for the intercept term
-m, n = X.shape
+ print('\nProgram paused. Press enter to continue.\n')
+ input()
 
-# Add intercept term to x and X_test
-X = np.concatenate((np.ones((m, 1)), X), axis=1)
+ ## ============= Part 3: Optimizing using fmin =============
 
-# Initialize fitting parameters
-initial_theta = np.zeros((n + 1, 1))
+ # Run fmin to obtain the optimal theta
+ # This function will return theta and the cost
+ initial_theta += 0.01
+ res = scipy.optimize.minimize(cost_function, initial_theta, args=(X, y), method='BFGS', jac=True,
+ options={'maxiter': 400})
+ theta = res['x']
 
-# Compute and display initial cost and gradient
-cost, grad = cost_function(initial_theta, X, y)
+ # Print theta to screen
+ print('Cost at theta found by fminunc: %f\n', res['fun'])
+ print('theta: \n')
+ print(np.array_str(theta).replace('[', ' ').replace(']', ' '))
 
-print('Cost at initial theta (zeros): {0:f}'.format(cost))
-print('Gradient at initial theta (zeros): ')
-print(np.array_str(grad.reshape(-1,1)).replace('[', ' ').replace(']', ' '))
+ # Plot Boundary
+ plot_decision_boundary(theta, X, y)
 
-print('\nProgram paused. Press enter to continue.\n')
-input()
+ # Put some labels
+ # Labels and Legend
+ plt.xlabel('Exam 1 score')
+ plt.ylabel('Exam 2 score')
 
-## ============= Part 3: Optimizing using fmin =============
+ # Specified in plot order
+ plt.legend(['Admitted', 'Not admitted'])
+ plt.show()
 
-# Run fmin to obtain the optimal theta
-# This function will return theta and the cost
-res = scipy.optimize.minimize(cost_function, initial_theta, args=(X, y), method='L-BFGS-B', jac=True, options={'maxiter': 400})
+ print('\nProgram paused. Press enter to continue.')
+ input()
 
-# Print theta to screen
-print('Cost at theta found by fminunc: %f\n', res['fun'])
-print('theta: \n')
-print(np.array_str(res['x']).replace('[', ' ').replace(']', ' '))
+ ## ============== Part 4: Predict and Accuracies ==============
 
-# Plot Boundary
-plot_decision_boundary(theta, X, y)
+ # Predict probability for a student with score 45 on exam 1
+ # and score 85 on exam 2
+ prob = sigmoid(np.array([1, 45, 85]).dot(theta))
+ print('For a student with scores 45 and 85, we predict an admission probability of {0:f}\n'.format(prob))
 
-# Put some labels
-# Labels and Legend
-plt.xlabel('Exam 1 score')
-plt.ylabel('Exam 2 score')
+ # Compute accuracy on our training set
+ p = predict(theta, X)
 
-# Specified in plot order
-plt.legend(['Admitted', 'Not admitted'])
+ print('Train Accuracy: {0:f}\n'.format(np.mean((p == y.flatten())) * 100))
 
-print('\nProgram paused. Press enter to continue.')
-input()
+ print('\nProgram paused. Press enter to continue.\n')
+ input()

ex2/ex2_reg.py

@@ -1,2 +1,92 @@
 import numpy as np
-import matplotlib.pyplot as plt
+import matplotlib.pyplot as plt
+import scipy
+
+from ex2 import plotdata, map_feature, sigmoid, plot_decision_boundary, predict
+
+
+#COSTFUNCTIONREG Compute cost and gradient for logistic regression with regularization
+# J = COSTFUNCTIONREG(theta, X, y, lambda) computes the cost of using
+# theta as the parameter for regularized logistic regression and the
+# gradient of the cost w.r.t. to the parameters.
+# TODO Also suffers from unstablity
+def cost_function_reg(theta, X, y, ld):
+ m = y.shape[0]
+ d1 = X.shape[1]
+ d2 = y.shape[1]
+ theta = theta.reshape((d1, d2))
+ h = sigmoid(X.dot(theta))
+ if theta.size > 0:
+ theta2 = theta
+ theta2[0, :] = 0
+ J = np.sum(-y * np.log(h) - (1 - y) * np.log(1 - h)) / m + ld/m/2 * np.sum(theta2 ** 2)
+ grad = (h - y).T.dot(X).T / m + ld / m * theta2
+ return J, grad.flatten()
+
+
+## Load Data
+# The first two columns contains the X values and the third column
+# contains the label (y).
+data = np.loadtxt('textdata/ex2data2.txt', delimiter=',')
+X = data[:, :2]
+y = data[:, 2:3]
+
+plotdata(X, y)
+
+# Put some labels
+# Labels and Legend
+plt.xlabel('Microchip Test 1')
+plt.ylabel('Microchip Test 2')
+
+# Specified in plot order
+plt.legend(['y = 1', 'y = 0'])
+plt.show()
+
+## =========== Part 1: Regularized Logistic Regression ============
+# Add Polynomial Features
+
+# Note that mapFeature also adds a column of ones for us, so the intercept
+# term is handled
+X = map_feature(X[:,:1], X[:,1:2])
+
+# Initialize fitting parameters
+initial_theta = np.zeros((X.shape[1], 1))
+
+# Set regularization parameter lambda to 1
+ld = 1
+
+# Compute and display initial cost and gradient for regularized logistic
+# regression
+[cost, grad] = cost_function_reg(initial_theta, X, y, ld)
+
+print('Cost at initial theta (zeros): {0:f}\n'.format(cost))
+
+print('\nProgram paused. Press enter to continue.\n')
+input()
+
+## ============= Part 2: Regularization and Accuracies =============
+# Initialize fitting parameters
+initial_theta = np.zeros((X.shape[1], 1)) + 0.01
+
+# Set regularization parameter lambda to 1 (you should vary this)
+ld = 1
+
+# Optimize
+res = scipy.optimize.minimize(cost_function_reg, initial_theta, args=(X, y, ld), method='BFGS', jac=True, options={'maxiter': 400})
+theta = res['x']
+
+# Plot Boundary
+plot_decision_boundary(theta, X, y)
+plt.title(print('lambda = {0:g}'.format(ld)))
+
+# Labels and Legend
+plt.xlabel('Microchip Test 1')
+plt.ylabel('Microchip Test 2')
+
+plt.legend(['y = 1', 'y = 0', 'Decision boundary'])
+plt.show()
+
+# Compute accuracy on our training set
+p = predict(theta, X)
+
+print('Train Accuracy: {0:f}\n'.format(np.mean(p == y.flatten()) * 100))