1. Introduction

In this practical, you will implement linear regression and get to see it work on real datasets. To get started with the exercise, you will need to download the starter code and unzip its contents to the directory where you wish to complete the exercise. You will be required to modify several parts of the code, by following the instructions in this assignment. Files needed for this practical include:


2. Linear regression with one variable

In this part of this exercise, you will implement linear regression with one variable to predict profits for a new laboratory. Suppose you work for a pharma company and are considering different cities for opening a new outstation. The company has already laboratories in various cities and you have data for profits and populations from the cities. You would like to use this data to select which city to expand to next. The file data1.txt contains the dataset for this linear regression problem. The first column is the population of a city and the second column is the profit of your laboratory in that city. A negative value for profit indicates a loss.


2.1. Plot the data

Before starting on any task, it is often useful to understand the data by visualizing it. For this dataset, you can use a scatterplot to visualize the data, since it has only two variables to plot (profit and population). Note that most problems that you will encounter in real life are multi-dimensional and can’t be plotted on a 2D plot.

data1 <- read.table("data1.txt", sep = ",")  # read comma-separated data
X <- data1[, 1]   
y <- data1[, 2] 
m <- length(y)  # Number of training examples

plot(X, y, xlab = "City population in 10,000s", ylab = "Profit in $10,000s",
     col = "blue", pch = 19)


2.2. Gradient descent

In this part, you will fit the linear regression parameters to our dataset using gradient descent.


2.2.1. Update equations

The objective of linear regression is to minimize the cost function: \[ J(\theta) = \frac{1}{2m}\sum\limits_{i=1}^m(h_{\theta}(x^{(i)}) - y^{(i)})^2 \] where the hypothesis \(h_\theta(x)\) is given by the linear model: \[ h_\theta(x) = \theta^Tx = \theta_0 + \theta_1x_1 \]

Recall that the parameters of your model are the \(\theta_j\) values. These are the values you will adjust to minimize cost \(J(\theta)\). One way to do this is to use the batch gradient descent algorithm. In batch gradient descent, each iteration performs the update (simultaneously for all \(j\)): \[ \theta_j := \theta_j - \alpha\frac{1}{m}\sum\limits_{i=1}^m(h_{\theta}(x^{(i)}) - y^{(i)})x_j^{(i)} \] With each step of gradient descent, your parameters \(\theta_j\) come closer to the optimal values that will achieve the lowest cost \(J(\theta)\).


2.2.2. Implementation

We store each example as a row in the the \(X\) matrix. To take into account the intercept term (\(\theta_0\)), we add an additional first column to \(X\) and set it to all ones. This allows us to treat \(\theta_0\) as simply another feature. We also initialize the initial parameters to 0 and the learning rate alpha to 0.01.

X <- cbind(rep(1, m), X)    # Add a column of ones to X
colnames(X) <- NULL
theta <- c(0, 0)            # Initialize fitting parameters
iterations <- 1500
alpha <- 0.01


2.2.3. Computing the cost \(J(\theta)\)

As you perform gradient descent to learn minimize the cost function \(J(\theta)\), it is helpful to monitor the convergence by computing the cost. In this section, you will implement a function to calculate \(J(\theta)\) so you can check the convergence of your gradient descent implementation. Your next task is to complete the code in the file computeCost.R, which is a function that computes \(J(\theta)\). As you are doing this, remember that the variables \(X\) and \(y\) are not scalar values, but matrices whose rows represent the examples from the training set. Once you have completed the function, the next code will run computeCost once using \(\theta\) initialized to zeros, and you will see the cost printed to the screen. You should expect to see a cost of 32.07.

source("computeCost.R")         # Load and run the computeCost function
J <- computeCost(X, y, theta)


2.2.4. Gradient descent

Next, you will implement gradient descent in the file gradientDescent.R The loop structure has been written for you, and you only need to supply the updates to \(\theta\) within each iteration. As you program, make sure you understand what you are trying to optimize and what is being updated. Keep in mind that the cost \(J(\theta)\) is parameterized by the vector \(\theta\), not \(X\) and \(y\). That is, we minimize the value of \(J(\theta)\) by changing the values of the vector \(\theta\), not by changing \(X\) or \(y\).

A good way to verify that gradient descent is working correctly is to look at the value of \(J(\theta)\) and check that it is decreasing with each step. The starter code for gradientDescent calls computeCost on every iteration and prints the cost. Assuming you have implemented both functions correctly, your value of \(J(\theta)\) should never increase, and should converge to a steady value by the end of the algorithm.

source("gradientDescent.R")         # Load and run the gradientDescent function
gd <- gradientDescent(X, y, theta, alpha, iterations)   
theta <- gd$theta
J_history <- gd$J_history
# Plot the convergence graph
plot(J_history, ylab = "Cost", xlab = "No. of iterations", pch = 19)

After you are finished, we will use your final parameters to plot the linear fit. The result should look something like this (see below). Your final values for \(\theta\) will also be used to make predictions on profits in areas of 70,000 and 120,000 people. Note the way that the following code lines use matrix multiplication, rather than explicit summation or looping, to calculate the predictions. This is an example of code vectorization.

# Plot linear fit
plot(X[, 2], y, xlab = "City population in 10,000s", ylab = "Profit in $10,000s", 
     col = "blue", pch = 19)
lines(X[, 2], X%*%theta, col = "red")

# Predict values for population sizes of 70,000 and 120,000
predict1 <- t(theta) %*% c(1, 7) * 10000
predict2 <- t(theta) %*% c(1, 12) * 10000


2.2.5. Visualizing \(J(\theta)\)

To understand the cost function \(J(\theta)\) better, you will now plot the cost over a 2-dimensional grid of \(\theta_0\) and \(\theta_1\) values. You will not need to code anything new for this part, but you should understand how the code you have written already is creating these images.

# Grid over which we will calculate 
theta0_vals <- seq(-10, 10, 0.2)
theta1_vals <- seq(-1, 4, 0.05)

# Initialize J_vals to a matrix of 0's
J_vals <- matrix(0, nrow = length(theta0_vals), 
                    ncol = length(theta1_vals))

# Fill out J_vals
for (i in 1:length(theta0_vals)){
  for (j in 1:length(theta1_vals)){
     t <- c(theta0_vals[i], theta1_vals[j])
     J_vals[i, j] <- computeCost(X, y, t)
  }
}

# Plot
plot(x = theta[1, ], y = theta[2, ], pch = 4, col = "red",  
     xlim = c(-10,10), ylim = c(-1,4), xlab = "theta0", ylab = "theta1")
contour(theta0_vals, theta1_vals, J_vals, col = "blue", 
        levels = seq(0, 500, 20), drawlabels = F, add = TRUE)

library(plotly)
plot_ly(z = J_vals, x = theta0_vals, y = theta1_vals) %>% 
  add_surface() %>% 
  layout(scene = list(xaxis = list(title = "theta0"), 
                      yaxis = list(title = "theta1"),
                      zaxis = list(title = "J(theta)")))

The purpose of these graphs is to show you that how \(J(\theta)\) varies with changes in \(\theta_0\) and \(\theta_1\). The cost function \(J(\theta)\) is bowl-shaped and has a global mininum (this is easier to see in the contour plot than in the 3D surface plot). This minimum is the optimal point for \(\theta_0\) and \(\theta_1\), and each step of gradient descent moves closer to this point.


3. Linear regression with multiple variables

In this part, you will implement linear regression with multiple variables. The dataset for this example, available in the file data2.txt comes from a study by Stamey et al. that examined the correlation between the level of prostate specific antigen (PSA) and a number of clinical measures, in 97 men who were about to receive a radical prostatectomy. The goal is to predict the log of PSA (lpsa) from a number of measurements including log cancer volume (lcavol), log prostate weight (lweight), age, log of benign prostatic hyperplasia amount (bph), seminal vesicle invasion (svi), log capsular penetration (lcp), Gleason score (gleason), and percent of Gleason scores 4 or 5 (pgg45). The next figure is a scatterplot matrix of the variables. Some correlations with lpsa are evident.

data2 <- read.table("data2.txt", sep = "\t", header = T)  # read tab-separated data
X <- as.matrix(data2[, -9])   
y <- data2[, 9] 
m <- length(y)  # Number of training examples

pairs(data2, pch = 20, col = "blue")


3.1. Feature scaling

By looking at the values, note that some variables are several times larger than others. When features differ by orders of magnitude, performing feature scaling can make gradient descent converge much more quickly. Your task here is to complete the code in featureScaling.R to i) subtract the mean value of each feature from the dataset, and ii) scale (divide) the feature values by their respective standard deviations. You will do this for all the features and your code should work with datasets of all sizes (any number of features / examples). Note that each column of the matrix \(X\) corresponds to one feature. After learning the parameters from the model, we often want to predict new responses that we have not seen before. Given a new \(X\) value, note that we must first normalize using the mean and standard deviation that we had previously computed from the training set. Thus, it is important to store those values.

# Scale and center features 
source("featureScaling.R")
fs <- featureScaling(X)

Xnorm <- fs$Xnorm
mu <- fs$mu
sigma <- fs$sigma

# Add a column of ones to X
Xnorm <- cbind(rep(1, m), Xnorm)


3.2. Cost function and gradient descent

Previously, you implemented gradient descent on a simple regression problem. The only difference now is that there are more features in the matrix \(X\). The hypothesis function and the batch gradient descent update rule remain unchanged. You should update the code in computeCost.R and gradientDescent.R to implement the cost function and gradient descent for linear regression with multiple variables. If your code in the previous part (single variable) already supports multiple variables, you can use it here directly. Make sure your code is well-vectorized.

The following starter code runs gradient descent with a particular learning rate (alpha). After making sure that your functions computeCost and gradientDescent already work with this starter code and support multiple variables, try running gradient descent with different values of alpha and see which one gives you the best result. Notice the changes in the convergence curves as the learning rate changes. With a small learning rate, you should find that gradient descent takes a very long time to converge to the optimal value. Conversely, with a large learning rate, gradient descent might not converge or might even diverge. Using the best learning rate that you found, run gradient descent until convergence to find the final values of \(\theta\). Next,use this value of \(\theta\) to predict the PSA levels for a new patient.

# Initialize parameters
alpha <- 0.01
iterations <- 1500
theta <- rep(0, ncol(Xnorm))

# Run gradient descent (try different values of alpha)
source("computeCost.R")
source("gradientDescent.R")
gd <- gradientDescent(Xnorm, y, theta, alpha, iterations);
theta <- gd$theta
J_history <- gd$J_history

# Plot the convergence graph
plot(J_history, ylab = "Cost", xlab = "No. of iterations", pch = 19)

# Display gradient descent's result 
theta

# Now predict the level of prostate specific antigen (PSA) for a patient with: 
#
#      lcavol   lweight   age  lbph  svi    lcp gleason pgg45
#        1.63      3.54    78  1.95    1  -1.49       7     0
#
# Hint: At prediction, make sure you do the same feature normalization


3.3. Normal equations

The closed-form solution to linear regression is \[ \theta = (X^TX)^{-1}X^Ty \]

Using this formula does not require any feature scaling, and you will get an exact solution in one calculation: there is no “loop until convergence” like in gradient descent. Complete the code in normalEqn.R to use the formula above to calculate \(\theta\). Remember that while you don’t need to scale your features, you still need to add a column of 1’s to the \(X\) matrix to have an intercept term (\(\theta_0\)). Once you have found \(\theta\), use it to predict the log level of prostate specific antigen (lpsa) for the patient above and compare the results with those obtained using gradient descent.

# Add a column of ones to X
X <- cbind(rep(1, m), X)   

source("normalEqn.R")
theta <- normalEqn(X, y)