<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-2-5-pca/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://enhancedatascience.com/2018/05/15/create-machine-learning-library-r-pca/"> Enhance Data Science</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>This is this second post of the “Create your Machine Learning library from scratch with R !” series. Today, we will see how you can implement <strong>Principal components analysis (PCA)</strong> using only the linear algebra available in R. Previously, we managed to implement <a rel="nofollow" href="http://enhancedatascience.com/2018/01/30/your-own-machine-learning-library-from-scratch-with-r/" target="_blank">linear regression and logistic regression from scratch</a> and next time we will deal with K nearest neighbors (KNN).</p>
<h2> Principal components analysis </h2>
<p>The PCA is a <a rel="nofollow" href="http://enhancedatascience.com/2017/07/31/machine-learning-explained-dimensionality-reduction/" target="_blank">dimensionality reduction</a> method which seeks the vectors which explains most of the variance in the dataset. From a mathematical standpoint, the PCA is just a coordinates change to represent the points in a more appropriate basis. Picking few of these coordinates is enough to explain an important part of the variance in the dataset.</p>
<h2> The mathematics of PCA </h2>
<p>Let <img src="https://i2.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-3981a72eefc3ea28c988345ea7f9de44_l3.png?resize=60%2C12" alt="mathbf{x}_1, ... mathbf{x}_n" width="60" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="12"> be the observations of our datasets, the points are in <img src="https://i2.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-f2b0b0e5a2f7c38423baad3d27eb9c50_l3.png?resize=20%2C12" alt="mathbb{R}^p" width="20" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="12">. We assume that they are centered and of unit variance. We denote <img src="https://i1.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-b6c5ac5177c5d386dccf31b7ba54e1f8_l3.png?resize=132%2C19" alt="mathbf{X}=(mathbf{x}_1, ..., mathbf{x}_n)^T" width="132" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="19"> the matrix of observations.<br>
Then, <img src="https://i1.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-cfefd5f68b490a0a3e436e457ad82793_l3.png?resize=41%2C16" alt="mathbf{X}^T mathbf{X}" width="41" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="16"> can be diagonalized and has real and positive eigenvalues (it is a symmetric positive definite matrix).<br>
We denote <img src="https://i0.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-5eb9e7bef30b00f2943efdb55bebde64_l3.png?resize=97%2C18" alt="lambda_1 &gt; ... &gt; lambda_p" width="97" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="18"> its ordered eigenvalues and <img src="https://i0.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-bf7f69de868eff3efbb0c5387d90ef07_l3.png?resize=65%2C14" alt="u_1, ... , u_p" width="65" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="14"> the associated eigenvectors. It <a rel="nofollow" href="https://en.wikipedia.org/wiki/Principal_component_analysis" target="_blank">can be shown</a> that <img src="https://i0.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-99033ed455edbd42603cb26f41646d7c_l3.png?resize=67%2C33" alt="frac{sum_{1leq i leq k} lambda_i}{sum_{1leq i leq p} lambda_i}" width="67" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="33"> is the cumulative variance explained by <img src="https://i0.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-7af5b053c4524cde415505619e3ea70b_l3.png?resize=65%2C12" alt="u_1, ..., u_k" width="65" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="12">.<br>
It can also <a rel="nofollow" href="https://en.wikipedia.org/wiki/Principal_component_analysis" target="_blank">be shown</a> that <img src="https://i0.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-7af5b053c4524cde415505619e3ea70b_l3.png?resize=65%2C12" alt="u_1, ..., u_k" width="65" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="12"> is the orthonormal basis of size <img src="https://i1.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-d42bc2203d6f76ad01b27ac9acc0bee1_l3.png?resize=9%2C13" alt="k" width="9" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="13"> which explains the most variances. </p>
<p><strong>This is exactly what we wanted ! We have a smaller basis which explains as much variance as possible !</strong></p>
<h3> PCA in R </h3>
<p>The implementation in R has three-steps:</p>
<ol>
<li>We center the data and divide them by their deviations. Our data now comply with PCA hypothesis.</li>
<li>We diagonalise <img src="https://i1.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-cfefd5f68b490a0a3e436e457ad82793_l3.png?resize=41%2C16" alt="mathbf{X}^T mathbf{X}" width="41" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="16"> and store the eigenvectors and eigenvalues</li>
<li>The cumulative variance is computed and the required numbers of eigenvectors <img src="https://i1.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-d42bc2203d6f76ad01b27ac9acc0bee1_l3.png?resize=9%2C13" alt="k" width="9" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="13"> to reach the variance threshold is stored. We only keep the first <img src="https://i1.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-d42bc2203d6f76ad01b27ac9acc0bee1_l3.png?resize=9%2C13" alt="k" width="9" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="13"> eigenvectors</li>
</ol>
<pre>
###PCA
my_pca&lt;-function(x,variance_explained=0.9,center=T,scale=T)
{
  my_pca=list()
  ##Compute the mean of each variable
  if (center)
  {
    my_pca[['center']]=colMeans(x)
  }
  ## Otherwise, we set the mean to 0 
  else
    my_pca[['center']]=rep(0,dim(x)[2])
  ####Compute the standard dev of each variable
  if (scale)
  {
    my_pca[['std']]=apply(x,2,sd)
  }
  ## Otherwise, we set the sd to 0 
  else
    my_pca[['std']]=rep(1,dim(x)[2])

  ##Normalization
  ##Centering
  x_std=sweep(x,2,my_pca[['center']])
  ##Standardization
  x_std=x_std%*%diag(1/my_pca[['std']])
  
  ##Cov matrix
  eigen_cov=eigen(crossprod(x_std,x_std))
  ##Computing the cumulative variance
  my_pca[['cumulative_variance']] =cumsum(eigen_cov[['values']])
  ##Number of required components
  my_pca[['n_components']] =sum((my_pca[['cumulative_variance']]/sum(eigen_cov[['values']]))&lt;variance_explained)+1
  ##Selection of the principal components
  my_pca[['transform']] =eigen_cov[['vectors']][,1:my_pca[['n_components']]]
  attr(my_pca, &quot;class&quot;) &lt;- &quot;my_pca&quot;
  return(my_pca)
}
</pre>
<p>Now that we have the transformation matrix, we can perform the projection on the new basis.</p>
<pre>
predict.my_pca&lt;-function(pca,x,..)
{
  ##Centering
  x_std=sweep(x,2,pca[['center']])
  ##Standardization
  x_std=x_std%*%diag(1/pca[['std']])
  return(x_std%*%pca[['transform']])
}
</pre>
<p>The function applies <a rel="nofollow" href="https://en.wikipedia.org/wiki/Change_of_basis" target="_blank">the change of basis</a> formula and a projection on the <img src="https://i1.wp.com/enhancedatascience.com/wp-content/ql-cache/quicklatex.com-d42bc2203d6f76ad01b27ac9acc0bee1_l3.png?resize=9%2C13" alt="k" width="9" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="13"> principals components.</p>
<h2> Plot the PCA projection</h2>
<p>Using the predict function, we can now plot the projection of the observations on the two main components. As in the part 1, we used the Iris dataset.</p>
<pre>
library(ggplot2)
pca1=my_pca(as.matrix(iris[,1:4]),1,scale=TRUE,center = TRUE)
projected=predict(pca1,as.matrix(iris[,1:4]))
ggplot()+geom_point(aes(x=projected[,1],y=projected[,2],color=iris[,5]))
</pre>
<figure class="wp-caption"><img data-orig-file="https://i0.wp.com/enhancedatascience.com/wp-content/uploads/2018/05/PCA_iris.png?fit=833%2C560" data-orig-size="833,560" data-medium-file="https://i0.wp.com/enhancedatascience.com/wp-content/uploads/2018/05/PCA_iris.png?fit=300%2C202" data-attachment-id="2930" data-image-meta="{" data-image-description="" src="https://i0.wp.com/enhancedatascience.com/wp-content/uploads/2018/05/PCA_iris.png?resize=640%2C430" data-permalink="http://enhancedatascience.com/2018/05/15/create-machine-learning-library-r-pca/pca_iris/" data-large-file="https://i0.wp.com/enhancedatascience.com/wp-content/uploads/2018/05/PCA_iris.png?fit=640%2C430" alt="PCA iris" data-image-title="PCA_iris" sizes="(max-width: 640px) 100vw, 640px" width="450" data-comments-opened="1" class="wp-image-2930" data-recalc-dims="1"><figcaption class="wp-caption-text">Projection of the Iris dataset on the two mains PCA</figcaption></figure>
<h2> Comparison with the FactoMineR implementation</h2>
<p>We can now compare our implementation with the standard <a rel="nofollow" href="http://factominer.free.fr/" target="_blank">FactoMineR</a> implementation of Principal Component Analysis.</p>
<pre>
library(FactoMineR)
pca_stats= PCA(as.matrix(iris[,1:4]))
projected_stats=predict(pca_stats,as.matrix(iris[,1:4]))$coord[,1:2]
ggplot(data=iris)+geom_point(aes(x=projected_stats[,1],y=-projected_stats[,2],color=Species))+xlab('PC1')+ylab('PC2')+ggtitle('Iris dataset projected on the two mains PC (FactomineR)')
</pre>
<p>When running this, you should get a plot very similar to the previous one. This ensures the sanity of our implementation.</p>
<figure class="wp-caption"><img data-orig-file="https://i2.wp.com/enhancedatascience.com/wp-content/uploads/2018/05/FactominerPCA.png?fit=833%2C560" data-orig-size="833,560" data-medium-file="https://i2.wp.com/enhancedatascience.com/wp-content/uploads/2018/05/FactominerPCA.png?fit=300%2C202" data-attachment-id="2993" data-image-meta="{" data-image-description="" src="https://i2.wp.com/enhancedatascience.com/wp-content/uploads/2018/05/FactominerPCA.png?resize=640%2C430" data-permalink="http://enhancedatascience.com/2018/05/15/create-machine-learning-library-r-pca/factominerpca/" data-large-file="https://i2.wp.com/enhancedatascience.com/wp-content/uploads/2018/05/FactominerPCA.png?fit=640%2C430" alt="" data-image-title="FactominerPCA" sizes="(max-width: 640px) 100vw, 640px" width="450" data-comments-opened="1" class="wp-image-2993" data-recalc-dims="1"><figcaption class="wp-caption-text">Projection of the Iris dataset using the FactoMineR implementation</figcaption></figure>

<p>Thanks for reading ! To find more posts on Machine Learning, Python and R, you can follow us on <a rel="nofollow" href="https://www.facebook.com/EnhanceDataScience/" target="_blank">Facebook</a> or <a rel="nofollow" href="https://twitter.com/AntGuilllot" target="_blank">Twitter</a>.<br>
. </p>
<p>The post <a rel="nofollow" href="http://enhancedatascience.com/2018/05/15/create-machine-learning-library-r-pca/" target="_blank">Create your Machine Learning library from scratch with R ! (2/5) – PCA</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/05/15/create-machine-learning-library-r-pca/"> 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 alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/1LlMoeEPHNA" height="1"></a></div>