<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/create-your-machine-learning-library-from-scratch-with-r-1-3/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://enhancedatascience.com/2018/01/30/your-own-machine-learning-library-from-scratch-with-r/"> Enhance Data Science</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>When dealing with Machine Learning problems in R, most of the time you rely on already existing libraries. This fastens the analysis process, but do you really understand what is behind the algorithms? Could you implement a logistic regression from scratch with R?<br>
The goal of this post is to create our own basic machine learning library from scratch with R. We will only use the linear algebra tools available in R. There will be three posts:</p>
<ol>
<li>Linear and logistic regression (this one)</li>
<li>PCA and k-nearest neighbors classifiers and regressors</li>
<li>Tree-based methods and SVM</li>
</ol>
<h2>Linear Regression (Least-Square)</h2>
<p>The goal of liner regression is to estimate a continuous variable <img src="https://i0.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-82606c3098bb09002088b0f6f9ffbb2a_l3.png?resize=14%2C12" alt="Y" width="14" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="12"> given a matrix of observations <img src="https://i2.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-d4ee28752517d6062a3ca0314890342d_l3.png?resize=16%2C12" alt="X" width="16" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="12">. Before dealing with the code, we need to derive the solution of the linear regression.</p>
<h3>Solution derivation of linear regression</h3>
<p>Given a matrix of observations <img src="https://i2.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-d4ee28752517d6062a3ca0314890342d_l3.png?resize=16%2C12" alt="X" width="16" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="12"> and the target <img src="https://i0.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-82606c3098bb09002088b0f6f9ffbb2a_l3.png?resize=14%2C12" alt="Y" width="14" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="12">. The goal of the linear regression is to minimize the <img src="https://i1.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-b050896aa188721e07a72c4ee2e6e85a_l3.png?resize=20%2C12" alt="L2" width="20" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="12"> norm between <img src="https://i0.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-82606c3098bb09002088b0f6f9ffbb2a_l3.png?resize=14%2C12" alt="Y" width="14" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="12"> and a linear estimate of <img src="https://i0.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-82606c3098bb09002088b0f6f9ffbb2a_l3.png?resize=14%2C12" alt="Y" width="14" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="12"> : <img src="https://i1.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-9d1e5a4bc3f1c144a658d071491dcc08_l3.png?resize=73%2C17" alt="hat{Y}=WX" width="73" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="17">. Hence, linear regression can be rewritten as an optimization problem: <img src="https://i1.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-ffec1d8252a45658e70bb3924e63e451_l3.png?resize=178%2C19" alt="arg ; min_W ||Y-WX||^2" width="178" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="19">. <a rel="nofollow" href="https://en.wikipedia.org/wiki/Linear_least_squares_(mathematics)" target="_blank">A closed-form solution can easily be derived</a> and the optimal <img src="https://i1.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-4caed22919a1780df1b6310b338b904e_l3.png?resize=19%2C12" alt="W" width="19" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="12"> is <img src="https://i1.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-c46d5449aa19c859253857625077e301_l3.png?resize=159%2C21" alt="hat{W}=(X^T X)^{-1}X^T Y" width="159" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="21"></p>
<h3> Linear regression in R </h3>
<p>Using the closed-form solution, we can easily code the linear regression. Our linear model object will have three methods, an init method where the model is fitted, a predict method to work with new data and a plot method to visualize the residuals’ distribution.</p>
<pre>
###Linear model
fit_lm&lt;-function(x,y,intercept=TRUE,lambda=0)
{
  ##Conversion to matrix if required
  if (!is.matrix(x))
  {
    x=as.matrix(x)
  }
  if (!is.matrix(y))
  {
    y=as.matrix(y)
  }
  #Add the intercept coefficient
  if (intercept)
  {
    x=cbind(x,1)
  }
  my_lm=list(intercept=intercept)
  ##Compute coefficients estimates
  my_lm[['coeffs']]=solve(t(x) %*% x) %*% t(x) %*% y
  ##Compute estimates for the train dataset
  my_lm[['preds']]=x %*% my_lm[['coeffs']]
  my_lm[['residuals']]=my_lm[['preds']]-y
  my_lm[['mse']]=mean(my_lm[['residuals']]^2)
  attr(my_lm, &quot;class&quot;) &lt;- &quot;my_lm&quot;
  return(my_lm)
}
</pre>
<p>The fit function is simple, the first few lines transform the data to matrices and add an intercept if required. Then, the ‘my_lm’ object is created and the coefficients are computed. The <em>solve()</em> function is used to invert the matrix and <em>%*%</em> denotes matrix multiplication. A the end, the residuals and the estimates are computed and the class of the object is set as ‘my_lm’.<br>
Now let’s implement the predict and plot methods for the my_lm class:</p>
<pre>

predict.my_lm&lt;-function(my_lm,x,..)
{
  if (!is.matrix(x))
  {
    x=as.matrix(x)
  }
  if (my_lm[[&quot;intercept&quot;]])
  {
    x=cbind(x,1)
  }
  x%*%my_lm[[&quot;coeffs&quot;]]
}

plot.my_lm&lt;-function(my_lm,bins=30,..)
{
  library(ggplot2)
  qplot(my_lm[[&quot;residuals&quot;]], geom=&quot;histogram&quot;,bins=bins) + xlab('Residuals values') + ggtitle('Residual distribution')
  
}

</pre>
<p>You can test the code on some preinstalled R dataset such as the <em>car</em> one. The code will give you the same coefficients estimates as the <em>lm</em> function. For instance, on the <em>car</em> dataset:</p>
<pre>
my_lm1=fit_lm(cars[,1],cars[,2])
vanilla_lm=lm(dist~speed,cars)
print(vanilla_lm[['coefficients']])
print(my_lm1[['coeffs']])
</pre>
<h2>Logistic regression</h2>
<p>Previously, we worked on regression and the estimation of a continuous variable. Now, with logistic regression, we try to estimate a binary outcome (for instance, ill vs healthy, pass vs failed, …). Again, let’s deal with the maths first:</p>
<h3>The mathematics of logistic regression</h3>
<p>The goal is to estimate a binary outcome <img src="https://i0.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-82606c3098bb09002088b0f6f9ffbb2a_l3.png?resize=14%2C12" alt="Y" width="14" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="12"> given the observations <img src="https://i2.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-d4ee28752517d6062a3ca0314890342d_l3.png?resize=16%2C12" alt="X" width="16" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="12">. We assume that <img src="https://i0.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-e22e63fb3474ccf828f0065bcd6c31fe_l3.png?resize=35%2C18" alt="Y|X" width="35" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="18"> follows a Bernoulli distribution of parameter <img src="https://i2.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-7d0e09763bb4ab368f952f4ad094f3aa_l3.png?resize=197%2C27" alt="theta(X)=frac{e^{WX}}{1+e^{WX}}=sigma(WX)" width="197" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="27">. <img src="https://i2.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-1c9cc40f96a1492e298e7da85a2c1692_l3.png?resize=11%2C8" alt="sigma" width="11" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="8"> is called the <a rel="nofollow" href="https://en.wikipedia.org/wiki/Sigmoid_function" target="_blank">sigmoid function</a>.<br>
Hence, we have <img src="https://i0.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-4014eb657f984013749cd41c17bb7cea_l3.png?resize=231%2C19" alt="P(y|x)=theta(X)^y(1-theta(X))^{1-y}" width="231" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="19">.<br>
We want to maximize the log-likelihood of the observed sample(over <img src="https://i2.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-356a08e839ab6974a16448e16e56745d_l3.png?resize=9%2C12" alt="theta" width="9" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="12"> and hence over <img src="https://i1.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-4caed22919a1780df1b6310b338b904e_l3.png?resize=19%2C12" alt="W" width="19" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="12">):</p>
<p><span>   </span><span>   </span><img src="https://i1.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-76b5ed9b95691d7d97eb1be5ef315a33_l3.png?resize=361%2C40" width="361" alt="[l(theta)=sum_{x,y} ; y.log(theta(X)) + (1-y).log(1-theta(X)) ]" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="40"></p>
<p>This maximization will be done using <a rel="nofollow" href="https://en.wikipedia.org/wiki/Newton%27s_method_in_optimization" target="_blank">Newton’s Method</a>. Newton’s method is a variant of gradient descent in which we try to find the optimal curvature of the function to increase the speed of convergence. If you are not familiar with the Newton method, you can just see it as a variant of batch <a rel="nofollow" href="https://en.wikipedia.org/wiki/Gradient_descent" target="_blank">gradient descent.</a>. The weights updates has the following form:</p>
<p><span>   </span><span>   </span><img src="https://i1.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-e9aff433234353f5b6fa753162c7e339_l3.png?resize=261%2C21" width="261" alt="[ W_{t+1}=W_t+(nabla^2 l(W_t))^{-1} nabla l(W_t) ]" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="21"></p>
<p>with: The Hessian </p>
<p><span>   </span><span>   </span><img src="https://i2.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-167f0d0fbfb3b9064f228c660e208538_l3.png?resize=402%2C21" width="402" alt="[ nabla^2 l(W_t))=-X^T Diag(sigma( W^T X)_i(1-sigma( W^T X))_i) X]" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="21"></p>
<p> and the gradient </p>
<p><span>   </span><span>   </span><img src="https://i2.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-a6e1559ce670e8eddd57c3f010335849_l3.png?resize=227%2C21" width="227" alt="[nabla l(W_t)=X^T.(y-sigma( W^T X))]" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="21"></p>
<p>The algorithm in R will update the weights using this update until the termination criterion is met. Here, the termination criterion is met when the mean square error is below the user-defined tolerance.</p>
<h3>Logistic regression in R</h3>
<pre>
###Sigmoid function
sigmoid=function(x) {1/(1+exp(-x))}

###Fit logistic regression
fit_logit=function(x,y,intercept=T,tol=10e-5,max_it=100)
{
  ##Type conversion
  if (!is.matrix(x))
  {
    x=as.matrix(x)
  }
  if (!is.matrix(y))
  {
    y=as.matrix(y)
  }
   ##Add intercept is required
  if (intercept)
  {
    x=cbind(x,1)
  }
  ##Algorithm initialization
  iterations=0
  converged=F
  ##Weights are initialized to 1 
  coeffs=matrix(1,dim(x)[2])
  
  ##Updates the weight until the max number of iter
  ##Or the termination criterion is met
  while (iterations&lt;max_it&amp; !converged)
  {
    iterations=iterations+1
    nu&lt;-sigmoid(x %*% coeffs)
    old_pred=sigmoid(x %*% coeffs)
    nu_diag=diag(nu[,1])
    ##Weights update
    coeffs=coeffs + solve(t(x) %*% nu_diag %*% x)%*% t(x) %*% (y-nu)
    ##compute mse to check termination
    mse=mean((y-sigmoid(x%*%coeffs))^2)
    ##Stop computation if tolerance is reached
    if (mse&lt;tol)
    {
      converged=T
    }
      
  }
  ##Creates the logit objects 
  my_logit=list(intercept=intercept)
  my_logit[['coeffs']]=coeffs
  my_logit[['preds']]=sigmoid(x%*%coeffs)
  my_logit[['residuals']]=my_logit[['preds']]-y
  my_logit[['mse']]=mean(my_logit[['residuals']]^2)
  my_logit[['iteration']]=iterations
  attr(my_logit, &quot;class&quot;) &lt;- &quot;my_logit&quot;
  return(my_logit)

}

##Predict the outcome on new data
predict.my_logit&lt;-function(my_logit,x,probs=T,..)
{
  if (!is.matrix(x))
  {
    x=as.matrix(x)
  }
  if (my_logit[['intercept']])
  {
    x=cbind(x,1)
  }
  if (probs)
  {
    sigmoid(x %*% my_logit[['coeffs']])
  }
  else
  {
    sigmoid(x %*% my_logit[['coeffs']])>0.5
  }
}
</pre>
<p>The code is split into two part:</p>
<ul>
<li>
The fit part, where the logistic model is fitted using Newton’s method.<br>
This part has three main components. First, the data is put in the proper matrix format and, if required, the intercept is added. Then the algorithm parameters are updated (initially all the weights are set to one).<br>
Finally, the algorithm updates the weights until the MSE goes below the tolerance. The most important line is probably the weights update one where the update formula is used to update the weights of the model.
</li>
<li>
The predict method, where the outcome is estimated using the weights computed previously.
</li>
</ul>
<p>We can now use our logistic regression to predict the class of a flower from the iris dataset:</p>
<pre>
fit_logit(iris[,1:4],iris[,5]=='setosa')
</pre>
<p>As expected, the algorithm can predict efficiently if a flower is a setosa or not.</p>
<p>If you like this post, follow us to learn how to create your Machine Learning library from scratch with R!</p>
<p>The post <a rel="nofollow" href="http://enhancedatascience.com/2018/01/30/your-own-machine-learning-library-from-scratch-with-r/" target="_blank">Create your Machine Learning library from scratch with R ! (1/3)</a> appeared first on <a rel="nofollow" href="http://enhancedatascience.com/" target="_blank">Enhance Data Science</a>.</p>
	<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="http://enhancedatascience.com/2018/01/30/your-own-machine-learning-library-from-scratch-with-r/"> Enhance Data Science</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 rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img width="1" alt="" src="http://feeds.feedburner.com/~r/RBloggers/~4/XNz4-ERsIak" height="1"></a></div>