<p><div><div><div><div data-width="100" data-show-faces="false" data-send="true" data-layout="button_count" data-height="21" data-href="https://www.r-bloggers.com/analytical-and-numerical-solutions-to-linear-regression-problems/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://datascienceplus.com/analytical-and-numerical-solutions-to-linear-regression-problems/"> DataScience+</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>This exercise focuses on linear regression with both analytical (normal equation) and numerical (gradient descent) methods. We will start with linear regression with one variable. From this part of the exercise, we will create plots that help to visualize how gradient descent gets the coefficient of the predictor and the intercept. In the second part, we will implement linear regression with multiple variables.</p>
<h2>Linear regression with one variable</h2>
<p>In this part of this exercise, we will implement linear regression with one variable to predict profits for a food truck. Suppose you are the CEO of a restaurant franchise and are considering different cities for opening a new outlet. The chain already has trucks in various cities and you have data for profits and populations from cities. We would like to use this data to help us select which city to expand to next.<br>
The file <a target="_blank" rel="nofollow" href="http://datascience-enthusiast.com/data/ex1data1.txt">ex1data1.txt</a> contains the dataset for our linear regression problem. The first column is the population of a city and the second column is the profit of a food truck in that city. A negative value for profit indicates a loss. The first column refers to the population size in 10,000s and the second column refers to the profit in $10,000s.</p>
<pre>library(ggplot2)
library(data.table)
library(magrittr)
library(caret)
library(fields)
library(plot3D)</pre>
<p><strong>Load the data and display first 6 observations</strong></p>
<pre>ex1data1 &lt;- fread(&quot;ex1data1.txt&quot;,col.names=c(&quot;population&quot;,&quot;profit&quot;))    
head(ex1data1)
<em> population profit
1: 6.1101 17.5920
2: 5.5277 9.1302
3: 8.5186 13.6620
4: 7.0032 11.8540
5: 5.8598 6.8233
6: 8.3829 11.8860</em></pre>
<h3>Plotting the Data</h3>
<p>Before starting on any task, it is often useful to understand the data by visualizing it. For this dataset, we can use a scatter plot to visualize the data, since it has only two properties to plot (profit and population). Many other problems that we encounter in real life are multi-dimensional and can’t be plotted on a 2-d plot.</p>
<pre>ex1data1%>%ggplot(aes(x=population, y=profit))+
geom_point(color=&quot;blue&quot;,size=4,alpha=0.5)+
ylab('Profit in $10,000s')+          
xlab('Population of City in 10,000s')+ggtitle ('Figure 1: Scatter plot of training data')+
theme(plot.title = element_text(size = 16,colour=&quot;red&quot;))</pre>
<p>Here it is the plot:<br>
<a title="" data-rel="lightbox-0" target="_blank" rel="nofollow" href="https://i1.wp.com/datascienceplus.com/wp-content/uploads/2017/02/fig2.png?ssl=1"><img sizes="(max-width: 490px) 100vw, 490px" alt="" class="wp-image-5555" src="https://i2.wp.com/datascienceplus.com/wp-content/uploads/2017/02/fig2-490x306.png?resize=456%2C285&ssl=1" data-recalc-dims="1"></a></p>
<h3>Gradient Descent</h3>
<p>In this part, we will fit linear regression parameters θ to our dataset using gradient descent.<br>
To implement gradient descent, we need to calculate the cost, which is given by:<br>
<a title="" data-rel="lightbox-1" target="_blank" rel="nofollow" href="https://i1.wp.com/datascienceplus.com/wp-content/uploads/2017/02/fig3-1.png?ssl=1"><img sizes="(max-width: 490px) 100vw, 490px" alt="" class="wp-image-5558" src="https://i2.wp.com/datascienceplus.com/wp-content/uploads/2017/02/fig3-1-490x361.png?resize=456%2C336&ssl=1" data-recalc-dims="1"></a></p>
<h3>Computing the cost J(θ)</h3>
<p>As we perform gradient descent to learn minimize the cost function J(θ), it is helpful to monitor the convergence by computing the cost. In this section, we will implement a function to calculate J(θ) so we can check the convergence of our gradient descent implementation.<br>
But remember, before we go any further we need to add the θ intercept term.</p>
<pre>X=cbind(1,ex1data1$population)
y=ex1data1$profit
head(X)
<em> [,1] [,2]
[1,] 1 6.1101
[2,] 1 5.5277
[3,] 1 8.5186
[4,] 1 7.0032
[5,] 1 5.8598
[6,] 1 8.3829</em></pre>
<p>The function below calcuates cost based on the equation given above.</p>
<pre>computeCost=function(X,y,theta){
    z=((X%*%theta)-y)^2
    return(sum(z)/(2*nrow(X)))
}</pre>
<p>Now, we can calculate the initial cost by initilizating the initial parameters to 0.</p>
<pre>theta=matrix(rep(0,ncol(X)))
round(computeCost(X,y,theta),2)
<em>32.07</em></pre>
<h3>Gradient descent</h3>
<p>Next, we will implement gradient descent by calling the <em>computeCost</em> function above. A good way to verify that gradient descent is working correctly is to look at the value of J(θ) and check that it is decreasing with each step. Our value of J(θ) should never increase, and should converge to a steady value by the end of the algorithm.<br>
The gradient descent function below returns the cost in every iteration and the optimal parameters for the number of iterations and learning rate alpha specified.</p>
<pre>gradientDescent=function(X, y, theta, alpha, iters){
gd=list()
cost=rep(0,iters)
for(k in 1:iters){
  z=rep(0,ncol(X))
  for(i in 1:ncol(X)){
      for(j in 1:nrow(X)){
          z[i]=z[i]+(((X[j,]%*%theta)-y[j])*X[j,i])
          }
      }
theta= theta-((alpha/nrow(X))*z)
cost[k]=computeCost(X,y,theta)
    }
    gd$theta= theta
    gd$cost=cost
    gd 
}</pre>
<p>Now, let’s use the <em>gradientDescent</em> function to find the parameters and we have to make sure that our cost never increases. Let’s use 1500 iterations and a learning rate alpha of 0.04 for now. Later, we will see the effect of these values in our application.</p>
<pre>iterations = 1500
alpha = 0.01
theta= matrix(rep(0, ncol(X)))
gradientDescent_results=gradientDescent(X,y,theta,alpha,iterations)
theta=gradientDescent_results$theta
theta
<em> [,1]
[1,] -3.630291
[2,] 1.166362</em></pre>
<h3>Ploting the cost function as a function of the number of iterations</h3>
<p>We expect that the cost should be decreasing or at least should not be at all increasing if our implementaion is correct. Let’s plot the cost fucntion as a function of number of iterations and make sure our implemenation makes sense.</p>
<pre>data.frame(Cost=gradientDescent_results$cost,Iterations=1:iterations)%>%
ggplot(aes(x=Iterations,y=Cost))+geom_line(color=&quot;blue&quot;)+
ggtitle(&quot;Cost as a function of number of iteration&quot;)+
theme(plot.title = element_text(size = 16,colour=&quot;red&quot;))</pre>
<p>Here it is the plot:<br>
<a title="" data-rel="lightbox-2" target="_blank" rel="nofollow" href="https://i0.wp.com/datascienceplus.com/wp-content/uploads/2017/02/fig6.png?ssl=1"><img sizes="(max-width: 490px) 100vw, 490px" alt="" class="wp-image-5564" src="https://i2.wp.com/datascienceplus.com/wp-content/uploads/2017/02/fig6-490x306.png?resize=456%2C285&ssl=1" data-recalc-dims="1"></a><br>
As the plot above shows, the cost decrases with number of iterations and gets almost close to convergence with 1500 iterations.</p>
<h3>Plot the linear fit</h3>
<p>Now, since we know the parameters (slope and intercept), we can plot the linear fit on the scatter plot.</p>
<pre>ex1data1%>%ggplot(aes(x=population, y=profit))+
geom_point(color=&quot;blue&quot;,size=3,alpha=0.5)+
ylab('Profit in $10,000s')+          
xlab('Population of City in 10,000s')+
ggtitle ('Figure 1: Scatter plot of training data') +
geom_abline(intercept = theta[1], slope = theta[2],col=&quot;red&quot;,show.legend=TRUE)+
theme(plot.title = element_text(size = 16,colour=&quot;red&quot;))+
annotate(&quot;text&quot;, x = 12, y = 20, label = paste0(&quot;Profit = &quot;,round(theta[1],4),&quot;+&quot;,round(theta[2],4),&quot; * Population&quot;))</pre>
<p>Here it is the plot:<br>
<a title="" data-rel="lightbox-3" target="_blank" rel="nofollow" href="https://i2.wp.com/datascienceplus.com/wp-content/uploads/2017/02/fig7.png?ssl=1"><img sizes="(max-width: 490px) 100vw, 490px" alt="" class="wp-image-5565" src="https://i0.wp.com/datascienceplus.com/wp-content/uploads/2017/02/fig7-490x306.png?resize=456%2C285&ssl=1" data-recalc-dims="1"></a></p>
<h3>Visualizing J(θ)</h3>
<p>To understand the cost function J(θ) better, let’s now plot the cost over a 2-dimensional grid of θ0 and θ1 values. The global minimum is the optimal point for θ0 and θ1, and each step of gradient descent moves closer to this point.</p>
<pre>
Intercept=seq(from=-10,to=10,length=100)
Slope=seq(from=-1,to=4,length=100)
# initialize cost values to a matrix of 0's
Cost = matrix(0,length(Intercept), length(Slope));
for(i in 1:length(Intercept)){
    for(j in 1:length(Slope)){
        t = c(Intercept[i],Slope[j])
        Cost[i,j]= computeCost(X, y, t)
    }
}
persp3D(Intercept,Slope,Cost,theta=-45, phi=25, expand=0.75,lighting = TRUE,
        ticktype=&quot;detailed&quot;, xlab=&quot;Intercept&quot;, ylab=&quot;Slope&quot;,
        zlab=&quot;&quot;,axes=TRUE, main=&quot;Surface&quot;)
image.plot(Intercept,Slope,Cost, main=&quot;Contour, showing minimum&quot;)
contour(Intercept,Slope,Cost, add = TRUE,n=30,labels='')
points(theta[1],theta[2],col='red',pch=4,lwd=6)</pre>
<p>Here it is the plot:<br>
<a title="" data-rel="lightbox-4" target="_blank" rel="nofollow" href="https://i0.wp.com/datascienceplus.com/wp-content/uploads/2017/02/fig8.png?ssl=1"><img sizes="(max-width: 490px) 100vw, 490px" alt="" class="wp-image-5566" src="https://i2.wp.com/datascienceplus.com/wp-content/uploads/2017/02/fig8-490x368.png?resize=456%2C342&ssl=1" data-recalc-dims="1"></a><br>
<a title="" data-rel="lightbox-5" target="_blank" rel="nofollow" href="https://i0.wp.com/datascienceplus.com/wp-content/uploads/2017/02/fig9.png?ssl=1"><img sizes="(max-width: 490px) 100vw, 490px" alt="" class="wp-image-5567" src="https://i2.wp.com/datascienceplus.com/wp-content/uploads/2017/02/fig9-490x368.png?resize=456%2C342&ssl=1" data-recalc-dims="1"></a></p>
<h3>Normal Equation</h3>
<p>Since linear regression has closed-form solution, we can solve it analytically and it is called normal equation. It is given by the formula below. we do not need to iterate or choose learning curve. However, we need to calculate inverse of a matrix , which make it slow if the number of records is very large. Gradient descent is applicable to other machine learning techniques as well. Further, gradient descent method is more appropriate even to linear regression when the number of observations is very large.<br>
θ=(X<sup>T</sup>X)<sup>−1</sup>X<sup>T</sup>y</p>
<pre>theta2=solve((t(X)%*%X))%*%t(X)%*%y
theta2
<em> [,1]
[1,] -3.895781
[2,] 1.193034</em></pre>
<p>There is very small difference between the parameters we got from normal equation and using gradient descent. Let’s increase the number of iteration and see if they get closer to each other. I increased the number of iterations from 1500 to 15000.</p>
<pre>iterations = 15000
alpha = 0.01
theta= matrix(rep(0, ncol(X)))
gradientDescent_results=gradientDescent(X,y,theta,alpha,iterations)
theta=gradientDescent_results$theta
theta
<em> [,1]
[1,] -3.895781
[2,] 1.193034</em></pre>
<p>As you can see, now the results from normal equation and gradient descent are the same.</p>
<h3>Using caret package</h3>
<p>By the way, we can use packages develoed by experts in the field and perform our machine learning tasks. There are many machine learning packages in R for differnt types of machine learning tasks. To verify that we get the same results, let’s use the caret package, which is among the most commonly used machine learning packages in R.</p>
<pre>my_lm &lt;- train(profit~population, data=ex1data1,method = &quot;lm&quot;)
my_lm$finalModel$coefficients
<em>(Intercept) -3.89578087831187
population 1.1930336441896</em></pre>
<h2>Linear regression with multiple variables</h2>
<p>In this part of the exercise, we will implement linear regression with multiple variables to predict the prices of houses. Suppose you are selling your house and you want to know what a good market price would be. One way to do this is to first collect information on recent houses sold and make a model of housing prices. The file <a target="_blank" rel="nofollow" href="http://datascience-enthusiast.com/data/ex1data2.txt">ex1data2.txt</a> contains a training set of housing prices in Port- land, Oregon. The first column is the size of the house (in square feet), the second column is the number of bedrooms, and the third column is the price of the house.</p>
<pre>ex1data2 &lt;- fread(&quot;ex1data2.txt&quot;,col.names=c(&quot;size&quot;,&quot;bedrooms&quot;,&quot;price&quot;))    
head(ex1data2)
<em> size bedrooms price
1: 2104 3 399900
2: 1600 3 329900
3: 2400 3 369000
4: 1416 2 232000
5: 3000 4 539900
6: 1985 4 299900</em></pre>
<h3>Feature Normalization</h3>
<p>House sizes are about 1000 times the number of bedrooms. When features differ by orders of magnitude, first performing feature scaling can make gradient descent converge much more quickly.</p>
<pre>ex1data2=as.data.frame(ex1data2)
for(i in 1:(ncol(ex1data2)-1)){
    ex1data2[,i]=(ex1data2[,i]-mean(ex1data2[,i]))/sd(ex1data2[,i])
}</pre>
<h3>Gradient Descent</h3>
<p>Previously, we implemented gradient descent on a univariate regression problem. The only difference now is that there is one more feature in the matrix X. The hypothesis function and the batch gradient descent update rule remain unchanged.</p>
<pre>X=cbind(1,ex1data2$size,ex1data2$bedrooms)
y=ex1data2$price
head(X)

iterations = 6000
alpha = 0.01
theta= matrix(rep(0, ncol(X)))
gradientDescent_results=gradientDescent(X,y,theta,alpha,iterations)
theta=gradientDescent_results$theta
theta
<em> [,1] [,2] [,3]
[1,] 1 0.13000987 -0.2236752
[2,] 1 -0.50418984 -0.2236752
[3,] 1 0.50247636 -0.2236752
[4,] 1 -0.73572306 -1.5377669
[5,] 1 1.25747602 1.0904165
[6,] 1 -0.01973173 1.0904165
[,1]
[1,] 340412.660
[2,] 110631.050
[3,] -6649.474</em></pre>
<h3>Normal Equation</h3>
<pre>theta2=solve((t(X)%*%X))%*%t(X)%*%y
theta2
<em> [,1]
[1,] 340412.660
[2,] 110631.050
[3,] -6649.474</em></pre>
<h3>Using caret package</h3>
<pre>ex1data2 &lt;- fread(&quot;ex1data2.txt&quot;,col.names=c(&quot;size&quot;,&quot;bedrooms&quot;,&quot;price&quot;))    
my_lm &lt;- train(price~size+bedrooms, data=ex1data2,method = &quot;lm&quot;, 
               preProcess = c(&quot;center&quot;,&quot;scale&quot;))
my_lm$finalModel$coefficients
<em>(Intercept) 340412.659574468
size 110631.050278846
bedrooms -6649.4742708198</em></pre>
<h2>Summary</h2>
<p>In this post, we saw how to implement numerical and analytical solutions to linear regression problems using R. We also used <code>caret</code> -the famous R machine learning package- to verify our results. The data sets are from the Coursera machine learning course offered by Andrew Ng. The course is offered with Matlab/Octave. I am doing the exercises in that course with R. You can get the code from <a target="_blank" rel="nofollow" href="https://github.com/fissehab/Andrew_Machine_Learning_R">this Github repository</a>.</p>
<ul></ul>
<p><strong>
<p>Related Post</p>
<p>
<ol>
<li><a title="Permanent Link toHow to create a loop to run multiple regression models" target="_blank" rel="nofollow" href="https://datascienceplus.com/how-to-create-a-loop-to-run-multiple-regression-models/">How to create a loop to run multiple regression models</a></li>
<li><a title="Permanent Link toRegression model with auto correlated errors – Part 3, some astrology" target="_blank" rel="nofollow" href="https://datascienceplus.com/regression-model-with-auto-correlated-errors-part-3-some-astrology/">Regression model with auto correlated errors – Part 3, some astrology</a></li>
<li><a title="Permanent Link toRegression model with auto correlated errors – Part 2, the models" target="_blank" rel="nofollow" href="https://datascienceplus.com/regression-model-with-auto-correlated-errors-part-2-the-models/">Regression model with auto correlated errors – Part 2, the models</a></li>
<li><a title="Permanent Link toRegression model with auto correlated errors – Part 1, the data" target="_blank" rel="nofollow" href="https://datascienceplus.com/fitting-a-simple-regression-model-with-auto-correlated-errors-an-astrological-example-part-1-the-data/">Regression model with auto correlated errors – Part 1, the data</a></li>
<li><a title="Permanent Link toLinear Regression from Scratch in R" target="_blank" rel="nofollow" href="https://datascienceplus.com/linear-regression-from-scratch-in-r/">Linear Regression from Scratch in R</a></li>
</ol></strong>
<div>
<a target="_blank" rel="nofollow" href="http://feeds.feedburner.com/~ff/datascienceplus?a=EHFUveK8Uis:yxaHcb_9-IU:D7DqB2pKExk"><img src="http://feeds.feedburner.com/~ff/datascienceplus?i=EHFUveK8Uis:yxaHcb_9-IU:D7DqB2pKExk" border="0"></a> <a target="_blank" rel="nofollow" href="http://feeds.feedburner.com/~ff/datascienceplus?a=EHFUveK8Uis:yxaHcb_9-IU:qj6IDK7rITs"><img src="http://feeds.feedburner.com/~ff/datascienceplus?d=qj6IDK7rITs" border="0"></a> <a target="_blank" rel="nofollow" href="http://feeds.feedburner.com/~ff/datascienceplus?a=EHFUveK8Uis:yxaHcb_9-IU:yIl2AUoC8zA"><img src="http://feeds.feedburner.com/~ff/datascienceplus?d=yIl2AUoC8zA" border="0"></a>
</div>
<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="https://datascienceplus.com/analytical-and-numerical-solutions-to-linear-regression-problems/"> DataScience+</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a title="The R Project for Statistical Computing" rel="nofollow" href="https://www.r-project.org/">R</a> news and <a title="R tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="https://www.r-bloggers.com/search/data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data">Big Data, <a title="R jobs" rel="nofollow" href="https://www.r-users.com/">R jobs</a>, visualization (<a title="ggplot and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="https://www.r-bloggers.com/search/map">maps</a>, <a title="Animation in R" rel="nofollow" href="https://www.r-bloggers.com/search/animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="https://www.r-bloggers.com/search/sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="https://www.r-bloggers.com/search/SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="https://www.r-bloggers.com/search/git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="https://www.r-bloggers.com/search/hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="https://www.r-bloggers.com/search/PCA">PCA</a>, <a title="Time series" rel="nofollow" href="https://www.r-bloggers.com/search/time+series">time series</a>, <a title="finance trading" rel="nofollow" href="https://www.r-bloggers.com/search/trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/b0gbephKqRk"></a></div>