C1_W3_lecture_nb_03_pca.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Another explanation about PCA\n",
    "\n",
    "\n",
    "<img src = './images/pca.jpeg' width=\"width\" height=\"height\"/>\n",
    "<sub>photo credit: Raunak Joshi</sub>\n",
    "\n",
    "\n",
    "In this lab, we are going to view another explanation about Principal Component Analysis(PCA). PCA is a statistical technique invented in 1901 by Karl Pearson that uses orthogonal transformations to map a set of variables into a set of linearly uncorrelated variables called Principal Components. \n",
    "\n",
    "PCA is based on the Singular Value Decomposition (SVD) of the Covariance Matrix of the original dataset. The Eigenvectors of such decomposition are used as a rotation matrix.  The Eigenvectors are arranged in the rotation matrix in decreasing order according to its explained variance. This last term is related to the EigenValues of the SVD.\n",
    "\n",
    "PCA is a potent technique with applications ranging from simple space transformation, dimensionality reduction, and mixture separation from spectral information.\n",
    "\n",
    "Follow this lab to view another explanation for PCA. In this case, we are going to use the concept of rotation matrices applied to correlated random data, just as illustrated in the next picture.\n",
    "\n",
    "<img src=./images/GaussianScatterPCA.svg>\n",
    "\n",
    "Source: https://en.wikipedia.org/wiki/Principal_component_analysis\n",
    "\n",
    "As usual, we must import the libraries that will use in this lab."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np                         # Linear algebra library\n",
    "import matplotlib.pyplot as plt            # library for visualization\n",
    "from sklearn.decomposition import PCA      # PCA library\n",
    "import pandas as pd                        # Data frame library\n",
    "import math                                # Library for math functions\n",
    "import random                              # Library for pseudo random numbers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To start, let us consider a pair of random variables x, y. Consider the base case when y = n * x. The x and y variables will be perfectly correlated to each other since y is just a scaling of x."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAXEklEQVR4nO3df2xd5X3H8c+n13KipUmLR2IYkDoNriK2RjAsEhR1ZRuRUjZBpG4FRrQgVUUaqjTGKs3IUaP+iEpbNaVSK23QTqOlK7RodaPFTWtos0oVYRg5I6Jrm/BjEI+StN36U4Em/e6Pe8OMe22fx/fce8+55/2SrJx7z8M9X4zz8cNznvM8jggBAHrfa7pdAACgMwh8AKgIAh8AKoLAB4CKIPABoCL6ul3AfM4999wYGhrqdhkAUCqPP/74DyNidbNzhQ38oaEhTU1NdbsMACgV2/813zmGdACgIgh8AKgIAh8AKoLAB4CKIPABoCIKO0sHAJZq1/gRfeHR53UmQjVbN266SB/c/uZul9V1BD6AnnLTPY/o20/9+JXXZyJ036HnJKnyoc+QDoCeMTfsZzsb+lVGDx9AT9i4+4B++tKZbpdRaAQ+gNLbMDahU2fYzGkxDOkAKLWNuw8Q9hnRwwdQWik9+y3rB9pcTfHRwwdQSilhb0mff9eV7S2oBAh8AKWzac9k5rDvs/TMnX/S5orKgSEdAKUxPj2jv/3iYWUdst+yfoCe/SwEPoBS2DV+JGku/eDKfsJ+DoZ0ABTe+PRMctg/Ora1jRWVEz18AIU2Pj2j2x44nLn9Xddfqu2XXdC+gkqMwAdQWAstlTAXvfrFEfgACmnTnkm9+LOXM7UdXrNCk7df1d6CegBj+AAKJyXsJRH2GdHDB1AoF9+xX6cTVkp4ljn2mdHDB1AYKWG/vGbCPhGBD6AQhkazh/2qZTV9d8817S2oB+US+La32f6e7WO2Rxdo93bbYXskj+sC6A1Do/szt92yfkBPvG9bG6vpXS0Hvu2apE9JepukSyTdaPuSJu1WSvprSY+2ek0AvWF8eiYp7IfXrODp2Rbk0cO/QtKxiHg6Il6WdL+k65q0+4CkD0s6lcM1AZRc6gNVgyv7mY3TojwC/wJJz896fbzx3its/76kiyJiwV/ltm+xPWV76uTJkzmUBqCIdo0fSQr74TUreKgqB22/aWv7NZL2SvrbxdpGxN0RMRIRI6tXr253aQC6IHURtC3rB+jZ5ySPefgzki6a9frCxntnrZT0e5IO2pak8yTts31tREzlcH0AJZGyVILEujh5yyPwH5M0bHud6kF/g6S/OHsyIn4i6dyzr20flPQewh6olq17D+roiV9kbs8c+/y1PKQTEaclvVvS1yT9p6QvRsSTtt9v+9pWPx9A+e0aP5I57Hmgqn1yWVohIiYkTcx5773ztL0qj2sCKIeNuw/opy+dydSWFS/bi7V0ALRN6hx7bs62F0srAGiLdQlh32dWvOwEAh9A7jaMTSjrgpeWdOxDjNl3AoEPIFdDo/t16ky2uN+yfkDPcIO2YxjDB5CL1Aeqdmxeqw9uf3MbK8JcBD6AlqU+UDW4sp+w7wKGdAC0ZHx6JinsWRene+jhA2hJyiJoLJXQXQQ+gCUZn57R3ySEPU/Pdh+BDyBZyro4y2tmO8KCIPABJNkwNpE07ZIdqoqDwAeQWcpSCauW1Qj7gmGWDoBMUsK+z2Kj8QIi8AEsKmVdnMGV/SyVUFAEPoAFXXzH/szr4jDHvtgYwwcwr3WjaWHPipfFRuAD+A1LWSqBsC8+Ah/Aq6TsUCWxS1WZMIYP4BU33fNIUtjv2LyWsC8RevgAJKUtgtZnNi0pIwIfQNIwDjtUlReBD1RcygNVrItTbozhAxVG2FcLgQ9U1MV3ZA97S4R9DyDwgQraNX5EpzM+UTW4sp+NxnsEY/hAxaQ8PctG472FwAcqJGXMfnjNCsK+xzCkA1TEpj2TmduyLk5voocPVEDKPHuGcXoXgQ/0sPHpGd2WcaNxgr73EfhAj0pZ8dISYV8BuYzh295m+3u2j9kebXL+dtvfsf2E7YdtvyGP6wJobuPuA0nLGzPtshpaDnzbNUmfkvQ2SZdIutH2JXOaTUsaiYiNkh6U9JFWrwuguZTx+sGV/XqWsK+MPHr4V0g6FhFPR8TLku6XdN3sBhHxzYj4ZePlIUkX5nBdAHNsGJvIHParltVY2rhi8gj8CyQ9P+v18cZ783mnpK82O2H7FttTtqdOnjyZQ2lAdawb3a9TZ7I9UtVn6Yn3bWtzRSiajs7Dt71D0oikjzY7HxF3R8RIRIysXr26k6UBpbZx94HMT8+yln115RH4M5IumvX6wsZ7r2L7akljkq6NiJdyuC4ApY/ZE/bVlUfgPyZp2PY62/2SbpC0b3YD25dJ+gfVw/5EDtcEoPpSCYzZI6uW5+FHxGnb75b0NUk1Sf8YEU/afr+kqYjYp/oQzmslfcm2JD0XEde2em2gqnaNH9F9h57L3N5izB45PXgVEROSJua8995Zx1fncR0A6WHPujg4iydtgRIZn55JCvu7rr9U2y9baNIcqoTAB0oiZakEqb42DmGP2Qh8oARSZuJI9OzRHIEPFNyGsYnMD1RJYqkEzIsNUIAC27RnMnPYb1k/QNhjQfTwgYJKGcbZsn5An3/XlW2uCGVH4AMFlLLROGGPrAh8oGBSNhon7JGCMXygQFLCfsfmtYQ9khD4QEGsS+zZsyUhUhH4QAGkjNkPruynZ48lIfCBLksJ++E1K1jxEkvGTVugS5ayVALDOGgFgQ90wda9B3X0xC8yt2epBOSBwAc6bHx6JnPYW9IzPD2LnBD4QAdt2jOpF3/2cqa2hD3yRuADHZIyx56wRzswSwfogJSwX7WsRtijLQh8oM02jE0s3qhhy/oB9p5F2zCkA7RR6hx7HqhCOxH4QJukDOMsr5mNxtF2DOkAORufnkkK++E1K/TdPde0sSKgjh4+kKNd40d036HnMrdneWN0EoEP5ISlElB0BD6Qg5SlEphjj24h8IEWXXzHfp3OOhVHhD26h8AHWpAy7VKSniXs0UXM0gGWKCXsl9dM2KPrCHxgCVJ79ky7RBEQ+ECiDWMT9OxRSgQ+kGDd6H6dOpMt7ndsXkvPHoWSy01b29skfUJSTdKnI+LOOeeXSfqspMsl/UjS9RHxbB7X/g33Xis9829t+WhUW4T09LKMjS35sKTD7asHFTDyTulP9+b2cS338G3XJH1K0tskXSLpRtuXzGn2Tkn/ExEXS/q4pA+3et2mCHu0kZ3w1e1i0RumPiP96+25fVweQzpXSDoWEU9HxMuS7pd03Zw210m6t3H8oKQ/tp3/3wnCHkCvefyfcvuoPAL/AknPz3p9vPFe0zYRcVrSTyT99twPsn2L7SnbUydPnsyhNAAouTiT20cV6qZtRNwdESMRMbJ69epulwMA3edabh+VR+DPSLpo1usLG+81bWO7T9LrVL95m691b839IwGgqy6/ObePyiPwH5M0bHud7X5JN0jaN6fNPkk7G8d/JukbEZHy3Eo2O/cR+mhJqD4bJ+sX0FY5z9JpeVpmRJy2/W5JX1N9WuY/RsSTtt8vaSoi9kn6jKTP2T4m6ceq/1Joj51zf9cA2WzcfUA/fSnbeGmfpWMf4oEqlEsu8/AjYkLSxJz33jvr+JSkP8/jWkA7pKx4SdijrAp10xbohqHR7GG/Zf0AYY/SIvBRaSl7z+7YvJbtCFFqBD4qa13iRuNsR4iyI/BRSTfd80jmFS+H16zQ5O1XtbMcoCPY8QqVMj49o9seOJy5PRuNo5cQ+KiMTXsm9eLPXs7cnnXs0WsIfFRCyhx7Sbrr+kvbVwzQJYzho+elhv2W9QPaftnc9f+A8iPw0dM27ZlMCnumXqKXMaSDnpUyZr+8ZrYjRM8j8NGTUh6okkTYoxIY0kHPSQl7i9k4qA4CHz1lw9jE4o0aBlf26xnCHhVC4KNnbNozqVNnsj0/u2pZTY+ObW1zRUCxMIaP0ts1fkT3HXouc/vBlf2EPSqJwEep3XTPI/r2Uz/O3J6lElBlBD5Ka9f4kaSw5+Ysqo7ARyml9uxZKgHgpi1KKKVnb9XDnqUSAHr4KJmtew/q6IlfZGrLzVng1Qh8lMaGsYnM0y6X10zYA3MwpINSSAl7i6USgGbo4aPwUpZKYBE0YH708FFoqRuNE/bA/Ojho7C27j2YeaNxZuIAiyPwUUgpu1SxQxWQDYGPQhmfntFtDxzO3H7Vsho7VAEZEfgojNSnZ5lnD6Qh8FEIKQ9USSyCBiwFgY+u2zV+JHPYD69Zocnbr2pvQUCPamlapu0B25O2jzb+PKdJm0ttP2L7SdtP2L6+lWuit2zdezDzWvarltUIe6AFrc7DH5X0cEQMS3q48XquX0r6y4j4XUnbJN1l+/UtXhc94OI79icN4zzxvm1trAbofa0G/nWS7m0c3ytp+9wGEfH9iDjaOP5vSSckrW7xuii5i+/Yr9MZJ9mz0TiQj1YDfzAiXmgc/0DS4EKNbV8hqV/SU/Ocv8X2lO2pkydPtlgaiiol7CWx0TiQk0Vv2tp+SNJ5TU6NzX4REWF73r/Gts+X9DlJOyPi183aRMTdku6WpJGRkYRIQFls3H0gc9j3WTr2IcIeyMuigR8RV893zvaLts+PiBcagX5innarJO2XNBYRh5ZcLUpt057JzE/PMu0SyF+rQzr7JO1sHO+U9JW5DWz3S/qypM9GxIMtXg8lND49o6HR/XrxZy9nan/X9ZcS9kAbtBr4d0raavuopKsbr2V7xPanG23eIekPJN1s+3Dj69IWr4uS2DV+JGmphB2b17IuDtAmjijmUPnIyEhMTU11uwy0IHWphC3rB1gXB2iR7ccjYqTZOdbDR1ukhv2OzWsJe6DNWFoBuds1fiRz2FvSx1nLHugIAh+5SunZr1pW4+lZoIMIfORm057JzDNx+sxSCUCnEfjIRcoOVRIPVAHdwE1btGzr3oOZw355zayLA3QJPXy0ZMPYhE6dyTa1t8/Sd/dc0+aKAMyHwMeSDY3uz9zWYhgH6DaGdJDs7FIJWW1ZP8CKl0AB0MNHkvHpmaSlEu5ijj1QGPTwkVlq2G9ZP0DYAwVCDx+Z7Bo/knnvWam+2ThLJQDFQuBjUVv3Hkzae5a17IFiIvCxoJSnZy22IwSKjMDHvFLm2C+vmTn2QMER+GgqZaNxHqgCyoFZOvgNQ6PZw37VshoPVAElQeDjVVIeqGJ5Y6BcCHy8YuPuA5nbLq+ZsAdKhsCHpLTljRmzB8qJwEdS2A+vWcGYPVBSzNKpsNSlEljHHig3Ar+iUh6okuqLoAEoN4Z0KmgpYc8iaED50cOvmK17D2YOe6ZdAr2FwK+QlJ49YQ/0HgK/ItaN7lfGh2eZYw/0KAK/AlKenpWYYw/0Km7a9riUp2clpl4CvYwefg9LWd5YIuyBXkfg96iUMfs+i6dngQpoaUjH9oDtSdtHG3+es0DbVbaP2/5kK9fEwsanZzSUEPZb1g8Q9kBFtNrDH5X0cETcaXu08frv5mn7AUnfavF6WEDqRuMM4QDV0upN2+sk3ds4vlfS9maNbF8uaVDS11u8HuYxPj2TOexXLasR9kAFtRr4gxHxQuP4B6qH+qvYfo2kj0l6z2IfZvsW21O2p06ePNliadWRuggac+yBalp0SMf2Q5LOa3JqbPaLiAjbzYaOb5U0ERHHbS94rYi4W9LdkjQyMpJ9ekmF3XTPI/r2Uz/O3J5F0IDqWjTwI+Lq+c7ZftH2+RHxgu3zJZ1o0uxKSW+xfauk10rqt/3ziBhdctWQVB+zTw17FkEDqqvVm7b7JO2UdGfjz6/MbRARN509tn2zpBHCvnVb9x7U0RO/yNR2eM0KTd5+VXsLAlB4rY7h3ylpq+2jkq5uvJbtEdufbrU4NHfxHfszh/2OzWsJewCSJEcUc6h8ZGQkpqamul1G4aQ8Pbtl/YA+/64r21wRgCKx/XhEjDQ7x1o6JZIS9oMr+wl7AK9C4JfExt0HMof9qmU1PTq2tc0VASgbAr8ENu2Z1E9fOpOpLRuXAJgPgV9w7FIFIC8EfoGl7D+7Zf0AYQ9gQSyPXEApc+wlZuMAyIbAL5iUIRyJsAeQHYFfIBt3H8h8c3ZwZT8zcQAkIfALImWjccIewFJw07YAUsJ+eM0Kwh7AkhD4XbYuIeyX18y6OACWjMDvog1jE5n3npWk7+65pm21AOh9BH6X7Bo/knmpBIn9ZwG0jsDvki88+nymdn0m7AHkg8DvkjMZlqUeXNmvYx8i7AHkg2mZXVKzFwx9tiMEkDd6+F1y46aLmr6/vGY9e+efEPYAckfgd8kHt79ZOzavVc2WVO/x79i8lpk4ANqGLQ4BoIewxSEAgMAHgKog8AGgIgh8AKgIAh8AKqKws3Rsn5T0X92uo4lzJf2w20VkRK3tQa3tU6Z6i1rrGyJidbMThQ38orI9Nd+Up6Kh1vag1vYpU71lqvUshnQAoCIIfACoCAI/3d3dLiABtbYHtbZPmeotU62SGMMHgMqghw8AFUHgA0BFEPiLsD1ge9L20caf58zTbq3tr9v+T9vfsT3U4VIz19pou8r2cduf7GSNs66/aK22L7X9iO0nbT9h+/oO17jN9vdsH7M92uT8MtsPNM4/2o3/5rNqWazW2xs/l0/Yftj2G7pRZ6OWBWud1e7ttsN216Y+ZqnV9jsa39snbf9zp2tMEhF8LfAl6SOSRhvHo5I+PE+7g5K2No5fK+m3ilpr4/wnJP2zpE8W9fsq6U2ShhvHvyPpBUmv71B9NUlPSXqjpH5J/yHpkjltbpX0943jGyQ90KXvZZZa//Dsz6SkvypyrY12KyV9S9IhSSNFrVXSsKRpSec0Xq/pRq1Zv+jhL+46Sfc2ju+VtH1uA9uXSOqLiElJioifR8QvO1bh/1u0VkmyfbmkQUlf70xZTS1aa0R8PyKONo7/W9IJSU2fIGyDKyQdi4inI+JlSferXvNss/8dHpT0x3ZjR5vOWrTWiPjmrJ/JQ5Iu7HCNZ2X5vkrSByR9WNKpThY3R5Za3yXpUxHxP5IUESc6XGMSAn9xgxHxQuP4B6oH5VxvkvS/tv/F9rTtj9quda7EVyxaq+3XSPqYpPd0srAmsnxfX2H7CtV7WU+1u7CGCyQ9P+v18cZ7TdtExGlJP5H02x2pbp46GprVOts7JX21rRXNb9Fabf++pIsiYn8nC2siy/f1TZLeZPvbtg/Z3tax6paATcwl2X5I0nlNTo3NfhERYbvZPNY+SW+RdJmk5yQ9IOlmSZ/Jt9Jcar1V0kREHG93ZzSHWs9+zvmSPidpZ0T8Ot8qq8X2Dkkjkt7a7VqaaXRI9qr+96cM+lQf1rlK9f9r+pbtN0fE/3azqPkQ+JIi4ur5ztl+0fb5EfFCI3ia/S/bcUmHI+Lpxj8zLmmz2hD4OdR6paS32L5V9XsN/bZ/HhHz3jzrYq2yvUrSfkljEXEo7xoXMCNp9k7zFzbea9bmuO0+Sa+T9KPOlNe0jrOa1SrbV6v+y/atEfFSh2qba7FaV0r6PUkHGx2S8yTts31tRHR6z9Ms39fjkh6NiF9Jesb291X/BfBYZ0pMw5DO4vZJ2tk43inpK03aPCbp9bbPji//kaTvdKC2uRatNSJuioi1ETGk+rDOZ9sR9hksWqvtfklfVr3GBztYm1T/bzpse12jjhtUr3m22f8OfybpG9G4c9dhi9Zq+zJJ/yDp2i6PMy9Ya0T8JCLOjYihxs/oIdVr7sYG11l+BsZV793L9rmqD/E83cEa03T7rnHRv1Qfk31Y0lFJD0kaaLw/IunTs9ptlfSEpCOS/klSf1FrndX+ZnVvls6itUraIelXkg7P+rq0gzVeI+n7qt83GGu8937VA0iSlkv6kqRjkv5d0hu7+HO6WK0PSXpx1vdxX1FrndP2oLo0Syfj99WqD0F9p/F3/4Zu1Zrli6UVAKAiGNIBgIog8AGgIgh8AKgIAh8AKoLAB4CKIPABoCIIfACoiP8Df2Yz0l1HF9MAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "np.random.seed(1)\n",
    "n = 1  # The amount of the correlation\n",
    "x = np.random.uniform(1,2,1000) # Generate 1000 samples from a uniform random variable\n",
    "y = x.copy() * n # Make y = n * x\n",
    "\n",
    "# PCA works better if the data is centered\n",
    "x = x - np.mean(x) # Center x. Remove its mean\n",
    "y = y - np.mean(y) # Center y. Remove its mean\n",
    "\n",
    "data = pd.DataFrame({'x': x, 'y': y}) # Create a data frame with x and y\n",
    "plt.scatter(data.x, data.y) # Plot the original correlated data in blue\n",
    "\n",
    "pca = PCA(n_components=2) # Instantiate a PCA. Choose to get 2 output variables\n",
    "\n",
    "# Create the transformation model for this data. Internally, it gets the rotation \n",
    "# matrix and the explained variance\n",
    "pcaTr = pca.fit(data)\n",
    "\n",
    "rotatedData = pcaTr.transform(data) # Transform the data base on the rotation matrix of pcaTr\n",
    "# # Create a data frame with the new variables. We call these new variables PC1 and PC2\n",
    "dataPCA = pd.DataFrame(data = rotatedData, columns = ['PC1', 'PC2']) \n",
    "\n",
    "# Plot the transformed data in orange\n",
    "plt.scatter(dataPCA.PC1, dataPCA.PC2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, what is the direction in which the variables point?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Understanding the transformation model pcaTr\n",
    "\n",
    "As mentioned before, a PCA model is composed of a rotation matrix and its corresponding explained variance. In the next module, we will explain the details of the rotation matrices. \n",
    "\n",
    "* `pcaTr.components_` has the rotation matrix \n",
    "* `pcaTr.explained_variance_` has the explained variance of each principal component"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Eigenvectors or principal component: First row must be in the direction of [1, n]\n",
      "[[-0.70710678 -0.70710678]\n",
      " [-0.70710678  0.70710678]]\n",
      "\n",
      "Eigenvalues or explained variance\n",
      "[1.66465230e-01 4.66587934e-33]\n"
     ]
    }
   ],
   "source": [
    "print('Eigenvectors or principal component: First row must be in the direction of [1, n]')\n",
    "print(pcaTr.components_)\n",
    "\n",
    "print()\n",
    "print('Eigenvalues or explained variance')\n",
    "print(pcaTr.explained_variance_)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$cos(45^o) = 0.7071$\n",
    "\n",
    "The rotation matrix is equal to:\n",
    "\n",
    "$$R = \\begin{bmatrix} cos(45^o) & sin(45^o) \\\\ -sin(45^o) & cos(45^o) \\end{bmatrix}$$ \n",
    "\n",
    "And $45^o$ is the same angle that form the variables y = 1 * x.\n",
    "\n",
    "Then, PCA has identified the angle in which point the original variables.\n",
    "\n",
    "And the explained Variance is around [0.166 0]. Remember that the Variance of a uniform random variable x ~ U(1, 2), as our x and y, is equal to:\n",
    "\n",
    "$$Var(x) = \\frac {(2 - 1)^2}{12} = 0.083333$$\n",
    "    \n",
    "Then the explained variance given by the PCA can be interpret as\n",
    "\n",
    "$$[Var(x) + Var(y),  \\ 0] = [0.0833 + 0.0833, \\ 0] = [0.166, \\ 0]$$\n",
    "\n",
    "Which means that all the explained variance of our new system is explained by our first principal component. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Correlated Normal Random Variables.\n",
    "\n",
    "Now, we will use a controlled dataset composed of 2 random variables with different variances and with a specific Covariance among them. The only way I know to get such a dataset is, first, create two independent Normal random variables with the desired variances and then combine them using a rotation matrix. In this way, the new resulting variables will be a linear combination of the original random variables and thus be dependent and correlated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "angle:  45.0\n",
      "rotationMatrix\n",
      "[[ 0.70710678  0.70710678]\n",
      " [-0.70710678  0.70710678]]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAqSUlEQVR4nO2df4wc533en+/uzZF7tK091odWWomiqgQkRNMko4tNg/9UqisqkUlfKNuKIQdIU1Qo0AAhITCgYsMkCyEiQNhCAQdohToICrMOJZG+UKUDWoYIGJFL2cfc0dRZZGrZ1o+VUTMmT7J5S3Hv7u0fe7M3Ozvv/HxnZmf3+QAEdHu7M+/u6p75zvN+f4hSCoQQQopLKe8FEEIISQaFnBBCCg6FnBBCCg6FnBBCCg6FnBBCCs5QHif98Ic/rNavX5/HqQkhpLCcP3/+n5VSY+7HcxHy9evXY2pqKo9TE0JIYRGRN7wep7VCCCEFh0JOCCEFh0JOCCEFh0JOCCEFh0JOCCEFJ5esFUKIGSan6zh65jLemWvgtmoF+3duwMS2Wt7LIhlDISekoExO1/HEyYtoNBcBAPW5Bp44eREAKOYDBoWckIJy9MzltojbNJqLOHrmciGFnHcX8aGQE1JQ3plrRHq8lxmEu4s0L1Tc7CSkoNxWrUR6vJfxu7swzeR0HTuOvIS7DpzGjiMvYXK6bvwcXud84uRF1OcaUFi5UJk6N4WckIKyf+cGVKxyx2MVq4z9OzfktKL4ZHV3kbag6kj7QkUhJ6SgTGyr4ak9m1GrViAAatUKntqzuZBWRFZ3F1lG/k7SvlDRIyekwExsqxVSuN3s37mhwyMH0rm7yGtf4bZqBXWPc5i6UDEiJ4TkTlZ3F3ntK6RtgzEiJ4T0BFncXWQV+bux31daWSsUckLIwJC2oAadO63zUMgJIT1JWnnX/bKv4IRCTgjpOQahQMgk3OwkhPQceaUJFhVG5ISQyKTdFyXLAqF+6O9CISekgOQpQFnYHmnnXQP9Zd/QWiEkQ0z0+cirzNwmC9sji/YD/WTfMCInJCNMRYA6ATp0ajaTSNKk7aG7s8giTbCfukdSyAnJCFP9w3VCM9doYnK6nrqYm7I9gi5saacJZmHfZAWtFUIywlQE6Cc0WdgCpmyPvK2NfuoeSSEnJCNM9fnwE5osbAF3X5RqxcJqq4R9x2ci+f55Wxv91D2S1gohGWGqz8fEthoOvzCLa/PNrt9lZQvYtkcS378XrI1+qfJkRE5IRpiMAA/u2tQTtkASe6SfrI28YUROSIaYigDzbP7kJIk9EvY99EvRTppQyAkpKL1gCyS1R4LeQz8V7aRJYmtFRO4QkbMi8mMRmRWRPzOxMEJI75PUHgkqkMo7s6UomIjIFwA8rpT6RxH5IIDzIvKiUurHBo5NCOlhklg8YaLtvDNbikJiIVdK/QLAL5b/+9ci8hqAGgAKOSEDQFyLJ0yBVC9kthQBo1krIrIewDYAr3j87jERmRKRqStXrpg8LSGkgISJtpnZEg5jQi4iHwBwAsBepdR77t8rpZ5RSo0rpcbHxsZMnZYQUlDCFEj1U9FOmhjJWhERCy0RP6aUOmnimISQ/iZsgVQvZOf0OomFXEQEwNcBvKaU+mryJRFSPJjrHJ1eyYXvB0xE5DsA/BGAiyIys/zYXyilvm3g2IT0PP2U65z1BYnRthlMZK38AwAxsBZCComp9rR504sXJN7phIOVnYQkpF9ynbO4IEUR5l68sPQqFHJCElLkXGensCrNc0xdkKIKc7/c6WQBux8SkpBeyHWOMwt0crqO/c9faM/+1HFbtWJk1mjUcvt+udPJAkbkhCQk7+yLuBbE4Rdm0Vz0k/DWBem+jWNGLI6owlzkO52soZATYoA8sy/iWhBegylsBGhfkExZHFGF2dQgjkGAQk7IMkXNkEjDgvjZkYfa/73v+IyR40cV5rzvdIoEhZwQFDtDIq4FUa1YmGt0R+XVimXk+G6ChFl3Ie31z78XoJATguwzJExG/3EtiEO7N2H/cxfQXFrxya2S4NDuTUaO74VOmIt8Ie0FKOSEINsMCdOiFdeCCPu6LCwOphomg0JOCLLNkEhDtOJaEH4Rslu4Xz5wf6y1hYGphsmgkBOCdDMk3KLodcEA0hWtXq+oZKphMlgQRAjS63tti6JddFOfa2gbE6UlWl5reOLkRW1RTx5zMnuhqKrIMCInZJk0MiS8RFGhlaftLMVJU7SiWjl52BxMNUwGhZyQFNGJn0Ir6o8jWlEzXqIK88hwGddvLnY9Xh2xPJ5tDqYaxodCTkiK6LzfWrUSa/Mwjn8dxX+enK57ijgAKP9q/sQUtSCrF6BHTkiKmPZ+4/jXUdbgd5x3PYqHTBHVxyedMCInJEWier9BUWkc/zrKGvyOk2YGCfPIk0EhJyRl3N6v3RLWLapetsn+5y7g8AuzmJtv4rZqBdURy7PZVZDIhvWf/dIj7QjefbG5b+MYzl66ksgSYR55MijkhDhI26f187i9otLmkmoLd32uAasksMrS0X7WZL77/M0Fz9+VS62kyS9NXsSxc2+2M27qcw1849yb7efFzTlnHnkyRKW9g+HB+Pi4mpqayvy8hPjhFlmgJZIm8sltdhx5Sbv56Telx4kzdXHNcBlWuYR3G81YFx77wmXnt/udf3TEwtx8M9Qao27mZvHZ9wMicl4pNe5+nBE5IcvE9WmjRPF+FoKfreHEKaStDJOV6H7f8RnsPT6DWghRd4tnkED79S93E9USYR55MijkhCwTx6eNmg7oZyF4tQmIitPyCLI4vC5cpohjiTCPPD5MPyR9SZwZkzrx8RMlXRR/+IVZz+d7pQJaJcH8zQXsOz6D1VYJ1YoFQasvuFXWFfQHE5SWGCVqrljlrj7lfs9Nu7TexAzRfoJCTvqOuDnJcXK+dWJ4bb7peT53T5dqxQKk9Xy1/Lr3F5bw6PZ1WLNqKHCmZhAm0gntvjOHdm/q+nwEwI671xrvUeMHc8674WYn6Tv8NhSDNuCiZq3ozgW0RHrNqiHfY+leH7TxGBa/9+y1wejEKguOfmZLV+qk1+eTZVVmku+36HCzkwwMSXKSo/q0+3duwF7NTMu5RrM9Sk3nWfv1YklK0N2Ec4PRSxibi6pro9fr88m67S1zzruhkJO+I+2cZHf0WbFKaDSXAl/XaC7i8WcvtH8+euayEcG2GR2xMDLsfwfgxhbmuw6c9lxLGHHMuiqTOefdUMhJ35H2kAh39GmVBVZJOmZf6lhUCnuPz6BcEixqnq+zVcLYLXEtDZ04lkRw14HTscr604qQ0/x+iwo3O0nfYWpIhFdmhGf15aLCB1YPtc8XBp2I16oVPLp9nedxwuR57zs+gy9NXgy5ihW8NnqB1oVHt6Fofz66daUVIac1BKTIMCInxAOd76vbGLw238S1+SbKIliMmUAgQHuzzln2HgUF4Ni5NzF+59pIwjaxrYapN652lN+7cdolQRulaUfIzDnvhBE56TtMpKcdOjXr6fsGEVfEgZYI25F/TRPNliU45ldYaUcbJd/67KUrgVG/bZf4FRMxQs4eCjnpO5LOnJycrrezTbLG7ni4/l94C/n2fz3qaYG4eWeuEfmCFsbTtu0SXcqlfVdBEc8WCjkpNF4Rp06Q6sviFoTJIcMC4Asaz1tHc0nh+z+96vm7n/+qgaf2bA6MzBWAvcdnIl3Qgjxt2y6ZnK5nPkCa+GNEyEXkr0XklyLyqonjERIGXcTpN1syzGagyWyL26oVPDmxGY9uXxfpdTqHxo6Ev/K5LaEicy9070+34QkAFauE1VYJ+47P4PFnL2izagY5cyRPTG12/g2ArwH4n4aOR0ggOgtl1VAJFavs6eE6NwPtY7gHJJjK7XZGsMdeibd56cUTJy/iqT2b8dSeze1iniibrLqo2bZDDp2a7bKWGs2ldq687jxq+RicvZk9RiJypdT3AHjfCxKSErrI8t1GE0/t2ax9nUJLrNzR/DfOvRmqjWxY7DU8cfJi5MHFI1ZJGx03mos4dGq2LZa1agVf+dyWUMcNU+25ZlW8+K5WrbAPSk5k5pGLyGMiMiUiU1euXMnqtKSP8etWOLGtps38AFrl82m1cAVaojaxrRa7Vexf7vmo78VortHsEssgyiJtj/xLkxe12SxxrCX7ApF0o5nEIzMhV0o9o5QaV0qNj42NZXVa0scEdSvcv3NDpE1GUzjXkMRvjyJ+fheLElr+tW2JOO8+7AvB/ucutMU86oalM92QfVDygVkrpLAEVfhNbKt5VklWrDJGfTZEk1CSlqjuPT6DrYe/47vx6se+4zPGbJ4lBFeFNpcUDp1q9VD32/R0Y3cctD/zWzQ9y92Ps5+4WVjZSQpNUIXfkxObMX7n2q7NNwChpvG4Bx1bZQEUtH1VnA/PNZqxIyUTG65lEXyoMhR6RJu9wek1dk13UXFH2rqsSOfjWXdLHASMCLmIfBPAvwHwYRF5G8BBpdTXTRybEC+iZEb4ib2uBa3Nx9aP4vs/vdrerLRKgj333u5byu4kuCdieiwqFWnOphP3Z6brAe62YeY053M+nnW3xEHAiJArpT5v4jiE6HAK9y0VC9dvLrQj5bgR3dQbwYlWL7/e+Zz55hL+1yvhRDxv7KKhsGmJbrvJ/Zm77068MmDCtJilj24eeuSk53GntM01ml0j0KJmRkxO13EsZmOqEN1qOxhOMHczCYtKhRZxqyw4uGtT+2evzxyqJfZ+HQfDjMuLMxuV+EOPnPQ8YVP47BJ8v6jcjjJN5osHkXTuZlwErX7iXmIeNIbOs13vksLI8BCmv/yA9pxe/rr72Ownbh4KOel5otxyOy2Wyel6R5XicFlwM2NRrVas3BpwKXjbKlZJcGj3Jt8LXtjNTS+CNqDDiD2JBoWc9Dx+WRNunBbL/ucudGSXZC3iAPCpLbfG7i3upGKVsdoqxd687CDA6bGbYnl9WqbsD/YTNws9ctLzRMlrBlpR49Ezl0ONXkuT0RELZy8lr2K2/eiDuzbFbpTlxB6qrEM3S9SvKRbzwvOFQk4yI+4fu134U9UUm7ipjliRMyBq1QpKBvckK1YZB3dtSpyJ4e7vvWrIzJ+s37p0v1NoFSq5v7vJ6Tr2P3+hs1L0+QsU8wyhkJNMCGqmFCTyUZo5KRXPAjAVwDszOpJaEXZFpP35mfLb/dbl9zuv7+7wC7NdG7rNRYXDL8waWSsJhh45yQS/IhD3rEg7ovuLkz/C/HLr1Cibhu82mji0e1OXR+6HySyWX77XOtbkdB3zNxcSHeu9G03t0Ock2O11vTYcvbJK3DgLeHS+vRE/n4SCQk4ywW9qj1eVZHNRdUR5USJRu/vh1BtXjWw0RqW51KoYdRfQxGFJdW/aJsUu/Akqk7dFXndmFvD0DhRykhrOiE+Xzyyin4YTB+cwhxPn8/VoTeWPxxVxqySAoKsa8+CuTYFl8s6skqDyfN3dUtg9DZIceuQkFkGetnsDTDtVJoZGOasLd9y9tl2K3tqsVO1xZGn2G+91yiJ45GN3YM3wSqw2OmLFajcbVK15aPem1kXDgZ2rTrKBQk4iE2YKjNcGmCmUAp5+ZCvmby7g5devti8SS6o1kszvwjEoLCqFE+frHZHyjeZKC68oZfJh2gUf/eyWjt8f/ewW5olniKgc/ocfHx9XU1NTmZ+XmEF3q233pgaA9QdOJzqHriDFRjeTk7TQzfC0S/Prc42uz7hilT37p5DeQUTOK6XG3Y/TIyeRMd29bnTEwsjwUFe/cL8Ws3FE3Cq1xC3nOqHU8bvIzTWa7ShdYeWCWWOZfKGhtUIiE+a2PMpG1z23frDrsYltNSNTfMqOiQbNpf4XcdsH95tX6sQWcWfRkQ2rNYsDhZxEJmjzK+of/Pdfv+rptx/ctQnlhOWW5QH7P/y9RitvPUpbA687qTD7IKR3GLD/zYkXUSMvv80vXQWiX3TtDpKdja+WEobQeTTKShs/gV5Uqp0P7v6OdN+B1x2WX3oi6T3okQ84cecnenWvm5yu4/FnL3huso0MD2FkeCh0BaW9jv6T4WSURfDUns0dVa9ubMF12yXu7xrQ9wHnFJ9iwYh8wNFFXo8/eyGSN2qLhC7t7525huftvs44KYswK8VFxSrjK59rpfU1FvyngXoJblAaoRNO8SkWjMgHHF2EZQty2Ag9qBeIXTZvP9fOULlv4xhOnK93RYkUce9sHvszDMoa1glu2D7gnOJTLCjkA06YoQ1hJpwH3XLP31xoj2FzH2f8zrU4/MKso8lSfoZKbfnikkePFjfX5psYGR7C049sbX9m9tQjPwStC/COIy+FTin0aqD11J7NnOJTEFgQNOB4+aZeCICfHXlI+3tdkZCbasXqGjMWdg1ZsOPutfj+61d7ypu3c71HRyy8O9+Ev6ni/Vq/PHGdd87ioN5DVxBEj3zAcfumzrxrJ0HeaNh0t7lGsyuNzXSL1iS83GMiDqzcn1wLEHGvb87ZGliXPnj4hVlmqBQcWiukw+6IktngPgYQ3PoUaInE3uMz2Ht8BqMjFvtWGyLoAuRlkU1O17WfPzNUigOFnHSQZMK587Vh0wwp4ubQ9Vdx4hZnv6ibGSrFgUJOuog74TyJ1x3UJIv4Y5UFj/zuHV0ZQG7c4uwXdTNDpTjQIyfGSOJ12xtyxB+rLKhYnX+2oyMWjn5mC56c2NzRZ8XtmXtZZLqou1qxuNFZIBiREwDe6WdO3zyM1ZLEU+U0GT0lWRkoHWRzufc7gr43Xb44h0IUC6YfDjh2XrK7N4qdfgZAa5dUKxZEgLn5Jm6rVnD9/QXtbE07P/v4D97qGl1WEqBcSj7fssj4tdi1SpLqoIawF2qSP7r0Qwr5ABPkadu36FEmzJcEHWLkzEeenK67Cn+AEauEGwtLfd9e1g87xxsA9h2f8dwrcA7tIIMLB0uQLoI87ThWyZJqebZ2lG5Hd7qLhq7x06AwOmJ1CPQ+zTANpgISPyjkA0yQOIgrug7LyPAQpr/8QMdjvVT000vMudIvdS0TmApI/DAi5CLyIID/CqAM4H8opY6YOC5Jl6A+K3HtjnfmGl2+axR7ZpBwC3ScZlX0uEni9EMRKQP4KwC/B+AeAJ8XkXuSHpeYQzc4IsoUmSgotCwC53SZIqNrWxCFL2xf5ztVycZumeAcArFqSP9nykk+BDATkX8MwE+UUj8FABH5WwCfBvBjA8cmCQkaHDH1xtXInf7CVBD2097lkoGEgCcnNmP8zrWhI+cbjr0Duz8NgK7UQq+LZFC3Skbw/YcJIa8BeMvx89sAPu5+kog8BuAxAFi3bp2B05IwBI3sOnE+fOTmzJyYnK77TrnvJ0ohLlx+2NF12IpZv+/Mb+PYiW7/I+5EKNLbZFbZqZR6Rik1rpQaHxsby+q0A4/fyC7dBqRIq4Kw4zGs9Lj+0uTF9h//IBAk4u7Pys21+WakKfRBY9bCbBxXRyxPO42zOPsTE0JeB3CH4+fblx8jPYDfyC5t1ooCHvndOzq8YWc71GPn3mQGioPmosLoiAVBK/q2q1Sd8h7Fu9Z9ZyUR3HXgdOCeg1UW/ObGgqdvzlmc/YkJIf8hgN8WkbtEZBjAHwI4ZeC4xABeG5r2Jpu2z8aIhRPn69pItJ/8b1Mo1Rq8Mf3lBzBz8AHUqpWuzyls5HvfRu871kWlAj/7WrWCNcNDXdWz9rk5i7M/SSzkSqkFAH8K4AyA1wA8q5Tyn0VFMsNv4K5O5JUCI+6IuFsTBEW+ukwiADh76UqsNYyOWNi/cwPe1bRJ0A3A5izO4mMkj1wp9W0A3zZxrEHHdEaB3/F0vcd11YUkPH6FPUEbjnFtjmvzreyWWyqWZ88b3QBsZq0UH/Za6SFMz070Op5VEnxg9VBXCb2ToPmb7B3ezeiI1VHN6vfZ64Zp2FlBYeef6tJAR0cs3GgucQZnH8KZnQXAdEaB1/GaSwrX5pvtTbD9z1/A1sPf6bjF9ysUKotQxD146KO3dvzstrSqFQsQ/4lIdiQeplBLoM+mmZtvau000p8wIu8h7jpw2lMkgybYRz2eH872tXbBCSPwYPwi3snpOh5/9kJgGmNZBEtK4bbllr9nL13x/PwFwKPb17V/74adEvsXRuQFwGRGweR0HXEqy73uACjiwdgDpd0bl7bFEqagyM5Kqc81cOJ8687o50cewtOPbO2Irp9+ZCuenNjMjUvShhF5D2HKI08yO9N5XmauxMP5nYXxu3Ved5jImuX2gwX7kfcAQX90pjIKkraMLYtQxH0I6jXjLKcPykCpajJMgHBFOnEHZZP+gkKeEWF7XJj4w0xSpcdIPJgPrh7C9ZsLvqPp7O8gqIXvXKOp3YNgkQ4JCz3yjMiyx0VcAbCzG5wtVL2eMxzQW6TfmWs0AQXfz8n+Dvbv3BDYi0Uh3MR7QnRQyDNCFyXbjai8KvziEqfPuNOP/c2NBc/nWCXBfRvHcHOAhyTbNJcU3mt4f06Aq8w+xMelAKYLktjQWskI3S223VUQMNdS1H5tmJQ3oDP6O3rmclefDpvmksI3X3nL83eDiN9n+81X3mr3H9d9nk6YMkiSwIg8I7yiZC9vNI7d4tW3Y2JbDV/53JauW3absohn9Bfkryfpy11U4kwIWlQKe5enJAVBG4UkhRF5RnhlpOj+yKNsVoaZAHTs3JsdFwx30c++4zM4euYy7ts4FjhEIcx0oH5jUSnjm8DO4h+mDJKkMI88R3Q5xlFus8McwyvtEQD2P3ch1G2/TcUq40Ory/h/v74Z+jX9QG35M7M/w6QTg9j3hMSFeeQ9SJyJ6W78NlG3/Zfv4OCuTR0pjX6zHv0QALePrsb//eX1SK8rOl7fx4cqQ/jNjYVIF0GbGiNwkgIU8hwxUQCka1kKtBo07X/+QvtcX5q82GWzhEUBhRbxsP1idty9Fj//VaPr7sV5wb0234RVFlQrFt5tNENH6GURbmiSVKCQ50zYAiCdPfLeDX03PaA1huzQqVkcOjWrFfxB4OlHtrY/v1sqlmdBz4671+LYf/xE12t3HHmpu4vkosKaVUOYOfhA6JYIg7a3QLKDQl4AvDY09z9/AVBAmLv7QRZwoGVnOC+YO4685PmZ/PxX3XbT5HQ9cFPafWcl4v291FipSVKCQl4APPuK92BRTi9mtHh53GEHENsXUB0KrYuCbYc59yGS7n0QEgXmkReAokw4//zH79DmredBtWJ5ZoeEbRccpvmYc0K9jd+cVELSgBF5AQhqvKSjJMDqoRLmm0sprKqbuBupphG0PHGncDr3GKojFqySdGSdRInc3Ti7HbrPxTxxkgWMyAuAV1WoVZbAL29JAe9naMH0gogD3iL+xMmLqM81oNDKOlnCSqOqsggevrd70zlK8zFb9N3n8orYCTENhbzHsaO7RnOxXSpeq1Zw9DNbcItP9z2bxRi5zkVGpGVtONsWPP7shS6LZHFJtS88i0rhxPl6l9jq2ip4YYt+ll0uCbGhtdJDuG/J79s4hhPn621hsDcS52+2uu7N+QzyHVQe/fi6rs3GMBuwbnvEZrVVah+nWrHwqS23dnwnQKctE3YjlRCTMCLvEbxuyY+de9Nzs+3afBNPnLyIkeForWr7nS9sX4cnJzbHnpDkFFv7+3BOvX9/YQnjd6713cg0OXeVkLAwIu8RvMTHL47kFJ9ObBEH4ke/TrH1s0hePnC/dvPSRNsFQqLCiLxHyPrWu2KVsePutZmeM01sEQf00a/dund0OWvFiVts41okTD0kecCIvEfwGzyRxnbl76y7xbOSsYi4KyZ1UbFTUINSBHXfRxiLhAORSdZQyDNGJyD3bRzz7Bv+O+tuwcuvXw19/LDVlVGO2csI0GVbhGlGFiS2tEhIkaCQZ4huCMTUG1dx4ny9Q8QFwMP31nD20pXQx69YZTx8by2XwhyBneGRTfGRjYL3WLykUbGJzpSEZAWFPEN0G2jffOWtrihaATh76Upo77ws0rYOvnHuTVNLDo3dXTBOBaoTeyDG5HQd+5+/ENhTJs1GVLRISFGgkGeITpR1VkjQSDgbqyxYMzzUHtm2ZriM6zezzWrZd3wm8V2A07pwR8TVEatrmAOtDkJaUMgzRCfKOl+7JIL6XMNzw3PNcBnzNxfbAme3ZU0aEcclrojXlj+TskhHBaQdDev6pdDqIGQFCnkCogqLbgPt4XtrXdWCwEqk7iWSS2rFzriWsMKzFqEpV9VnIlFUvrB9HcbvXOs7PNoJrQ5CvEmURy4inxWRWRFZEpGugaD9TJzmSLoc4ycnNnc8XgrRC9aOXpPmn9szJMPyqS23JvalRfyrMNmbhJBoJI3IXwWwB8B/N7CWQuEnQH5Roy6qtB+fnK5j7/GZUGsI66HrqFhl3LdxzHd4gpvTP/oF3k9YVbp6qIzxO1vFSOxNQkhyEkXkSqnXlFIDGTqlJUBRItHVVqndQCsqa4ZbBTJnL12JVO5/bb6ZuL+5M+JmbxJCkpNZib6IPCYiUyIydeVK+NzoXiUtAYpyIWg0l2L743byR16bo/b79GoVy2wUQqIRKOQi8l0RedXj36ejnEgp9YxSalwpNT42NhZ/xT1CWgKURiRq9zF30mgu4vALs8bPFZaSCO46cBpHz1zGw/fWUK2s9FZfbbEFECFRCPTIlVKfzGIhRSNs5V+czJawHnkYKlZZa50ERfOjI1bijBgddkZOfa6B4z94q2Nig92mF/Cu2iSEdML0w5DoBNlPaHQl+UB2AmWX+Ue1UATBQm+KpscUozAbx4SQFknTD/9ARN4G8AkAp0XkjJll9RZx5zDGSa3z+12cCfVnL13R2kBOO8NNlAIfL+vGjZ1uGQVmrhASjqRZK99SSt2ulFqllPqXSqmdphbWS8TNdY6T2eL3u1t8hNfveLr89UO7N3UJfFQEwOtP/T5+fuQh3+fdVq3gneUKzrAwc4WQcNBaCUHcVMM4Pa39+pLHqai0z6WzgabeuOrZtCvq8QF9hahgJTsm7HmYuUJIeJgesIxz6vqOIy912CZxUg0np+ueOd5BArV/5wZY5e6oNY7MBp1rcrqOE+frsUXcfXzd1PmoR+dUHUKiwYgcwZuSUYcMuI9nU61YOLR7U2Dl56FTs4n7mYxYJfxlgBjGGVJcW7ZIbnOU9u848lK7Q6FTtkesUuTiIbuNLSEkPBRyBJfbRx0yoBPINauGQkWZ7xpoSuUe8OCVdRM1k8Utsu4LljvLxU/ER0cs3GgucQIPIQagkCOcBx6l817S8v0k/VNsFIC9y/3J79s41tFd0b7jKMlKhWcQVkm6RDZORA+0BPvgrk3tY8RpS8uWtoSsQCFHskG7aRzvvo1jxqb81OcanqPfggTYWQyks4Tipgc6/e844tsL+fmE9BIUcoQftBs2CvQ6nlUWXH9/AXcdOB0YQermdMbZOESM14T1qePcOdSqlcRiG7fzJCH9ysBnrdji3GgutnOcvbImohQFufO2R0csQLXSB8MUFOkiXdMDlasVy7dfjF8mD+CdpeKHKQ+crW8J6WSghdwpzkArx9kWG3dkF7UoaGJbDS8fuB8/O/IQRoaHusrQ/V6bRiGMO6GxYpVxaPcmz0Ihuy960IVrYlsND99b8xyEUatW8IXt6zyPnRS2viWkk4GyVtzWyPX3F0LfoieJAqO+1suaCUKWfZdbKhau31zomD5vj5M7e+mKpy3kJa5h7IvJ6TqO//Ctrg3T0vJ7SMvmiJoOSki/MzBC7rVBpsNLYJNsYEZ9rS2AUbogigK++sjWdjQdpiujnf/t9ZwwF5+jZy53XDBslgAcOjWbmpBHTQclpN8ZGCGPkirnJbBJosA4r41aGOQUTxNdGcNcfPzuRkwNaNbBQcyErDAwQh52I0wnsGGiQDsSri83h1pUqj3c+OF7a+2eJmURPHxvsNhejzjGLax4hrFNwlx8TOS7E0KSMzBCrhOd0RELI8NDoW7R/aJAd5TrHJyw/7kLgKw8tqgUTpyvY/zOtb7VoV62hQnCFkDZ69B9Nn5DMEZHondqJITEY2CEXBdhHtzl3/skLH7WTZzBCXFS6cKKZ1jPPsi+mNhWw9QbV7uKl6yytCs3CSHpMzBCnvYGWRzh9XtNVNsiSDydG6DVEQtWSTouMHGzPp6c2IzxO9dy45GQHBEVs4VpEsbHx9XU1FTm502THUdeStyEysnkdB37js+EKgIqCfCh1RbebTS13r1Xpema4SHtawghvYeInFdKjbsfH5iIPG38cr+tkgCCrtzuoKyVqTeuevZJcbOkVjY6vTJQvGyf5qLCmlVDmDn4QJi3RwjpYQpT2RlULp43zrJ8AB3l/kc/uwVHP7MlcpXjkxOb8fQjWzteZ1dLAvoZnu6qUZa0E9LfFCIiL0q3uzCbg26Cind0xwyycpwibbq7IyGktyhERB53+HGvE6URl5sgP94p0l7NrVjSTkj/UAgh71drQHeBevzZC4EWkt80erdIu7sxciYmIf1FIayVIlsDftaJ7kLkLCbyspAmp/0HJnuJNEvaCelfChGRF9UaCLJOwlyI3BbS5HS9VSmqwWtwg8mN4l7fdCZkECmEkBfVGgjy9sMOZnBG7odOzXpWigL6qUZxfXg3Jo9FCDFHIawVoDjWgNNK0ZkftjC7q01Ly4223Dgjd7/GWF4XN5Nj0ThijZDepDBCXgS8Kii9cAqz8wLl9fooFpKXmJrcKO7XTWdCik4hrJWiEKbnuZ8wh7GQdI2xdI+bHIvGEWuE9CYUcoP4RaZhvX3nrM+XD9zf9dyDuzbBKnemHvo1zDK5UVzUTWdC+h1aKwbRpUn6NceKStQujia7PnLEGiG9CbsfGkTncRchw8ZJmJmfhJDsYffDDOiHiLUofW0IISskishF5CiAXQBuAngdwL9XSs0Fva4XIvK0o848oloT59Q14zJpDxFC4qGLyJNudr4I4CNKqY8C+CcATyQ8XiakXdiSR+GMqXMyxZCQ4pFIyJVS31FK2aPezwG4PfmS0iftbop5dGs0dc5bKtHSGAkh+WMy/fBPAPy9weOlRtpRZx5RrYlzTk7Xcf3mQtfjVkmYYkhIDxMo5CLyXRF51ePfpx3P+SKABQDHfI7zmIhMicjUlStXzKw+JmkXtuRROGPinEfPXO4YR2fzgdVD3OgkpIcJFHKl1CeVUh/x+Pd3ACAifwzgUwAeVT47p0qpZ5RS40qp8bGxMWNvIA5pF7bkUThj4py66H1uXt/fhRCSP4msFRF5EMCfA9itlJo3s6T0sUvhnWXtq4bifxTu1q4AMu/WaKJDJEvwCSkmSfPIvwZgFYAXpTWx5pxS6j8lXlUAplL7bjSX2v8912jGypfW5V0/tWdz5ul6STtE7t+5IVHTLkJIPiQScqXUb5laSFhMFayYasnaT61d+6GgiZBBpHCVnaaE01RmSb/lXRel7zshZIXCdT80JZym/GD6yoSQvCmckJsSTlOZJWztSgjJm8IJuSnhNDUHtKjzRAkh/UMh29iyzSohZBDpqza23JAjhJAVCmetEEII6YRCTgghBYdCTgghBYdCTgghBYdCTgghBSeX9EMRuQLgjcxPHI4PA/jnvBeRMv3+Hvv9/QF8j/1C1Pd4p1Kqqw94LkLey4jIlFeeZj/R7++x398fwPfYL5h6j7RWCCGk4FDICSGk4FDIu3km7wVkQL+/x35/fwDfY79g5D3SIyeEkILDiJwQQgoOhZwQQgoOhdwDETkqIpdE5Eci8i0Rqea9JpOIyGdFZFZElkSkr9K7RORBEbksIj8RkQN5r8c0IvLXIvJLEXk177WkhYjcISJnReTHy/+f/lneazKJiKwWkR+IyIXl93c46TEp5N68COAjSqmPAvgnAE/kvB7TvApgD4Dv5b0Qk4hIGcBfAfg9APcA+LyI3JPvqozzNwAezHsRKbMA4HGl1D0AtgP4z332Pb4P4H6l1BYAWwE8KCLbkxyQQu6BUuo7SqmF5R/PAbg9z/WYRin1mlLqct7rSIGPAfiJUuqnSqmbAP4WwKdzXpNRlFLfA3A173WkiVLqF0qpf1z+718DeA1A3wwgUC1+s/yjtfwvUdYJhTyYPwHw93kvgoSiBuAtx89vo48EYBARkfUAtgF4JeelGEVEyiIyA+CXAF5USiV6f4WcEGQCEfkugH/l8asvKqX+bvk5X0TrNu9YlmszQZj3R0gvIyIfAHACwF6l1Ht5r8ckSqlFAFuX99++JSIfUUrF3vcYWCFXSn3S7/ci8scAPgXg36oCJtsHvb8+pQ7gDsfPty8/RgqGiFhoifgxpdTJvNeTFkqpORE5i9a+R2whp7XigYg8CODPAexWSs3nvR4Smh8C+G0RuUtEhgH8IYBTOa+JREREBMDXAbymlPpq3usxjYiM2ZlwIlIB8O8AXEpyTAq5N18D8EEAL4rIjIj8t7wXZBIR+QMReRvAJwCcFpEzea/JBMsb1H8K4AxaG2TPKqVm812VWUTkmwD+D4ANIvK2iPyHvNeUAjsA/BGA+5f//mZE5PfzXpRBbgVwVkR+hFbw8aJS6n8nOSBL9AkhpOAwIieEkIJDISeEkIJDISeEkIJDISeEkIJDISeEkIJDISeEkIJDISeEkILz/wHCcpAV0njHngAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.lines as mlines\n",
    "import matplotlib.transforms as mtransforms\n",
    "\n",
    "np.random.seed(100)\n",
    "\n",
    "std1 = 1     # The desired standard deviation of our first random variable\n",
    "std2 = 0.333 # The desired standard deviation of our second random variable\n",
    "\n",
    "x = np.random.normal(0, std1, 1000) # Get 1000 samples from x ~ N(0, std1)\n",
    "y = np.random.normal(0, std2, 1000)  # Get 1000 samples from y ~ N(0, std2)\n",
    "#y = y + np.random.normal(0,1,1000)*noiseLevel * np.sin(0.78)\n",
    "\n",
    "# PCA works better if the data is centered\n",
    "x = x - np.mean(x) # Center x \n",
    "y = y - np.mean(y) # Center y\n",
    "\n",
    "#Define a pair of dependent variables with a desired amount of covariance\n",
    "n = 1 # Magnitude of covariance. \n",
    "angle = np.arctan(1 / n) # Convert the covariance to and angle\n",
    "print('angle: ',  angle * 180 / math.pi)\n",
    "\n",
    "# Create a rotation matrix using the given angle\n",
    "rotationMatrix = np.array([[np.cos(angle), np.sin(angle)],\n",
    "                 [-np.sin(angle), np.cos(angle)]])\n",
    "\n",
    "\n",
    "print('rotationMatrix')\n",
    "print(rotationMatrix)\n",
    "\n",
    "xy = np.concatenate(([x] , [y]), axis=0).T # Create a matrix with columns x and y\n",
    "\n",
    "# Transform the data using the rotation matrix. It correlates the two variables\n",
    "data = np.dot(xy, rotationMatrix) # Return a nD array\n",
    "\n",
    "# Print the rotated data\n",
    "plt.scatter(data[:,0], data[:,1])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us print the original and the resulting transformed system using the result of the PCA in the same plot alongside with the 2 Principal Component vectors in red and blue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Eigenvectors or principal component: First row must be in the direction of [1, n]\n",
      "[[ 0.70827652  0.7059351 ]\n",
      " [-0.7059351   0.70827652]]\n",
      "\n",
      "Eigenvalues or explained variance\n",
      "[1.09488457 0.11243237]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABA0UlEQVR4nO29eZgU5bn3/3m6p2emB5SRJVEGcEkMeV0QIjHmlRiVuEXBESNEXOJ2PDma1y0B0UQFkwg6OW5v4i8xJm+OW4AojrgF93jwxESQRVGJcYkymASQwcD0zPR0P78/aqqnl6rqqurqbfr+XBcXTHd11dND1bfuuleltUYQBEGoXkLlXoAgCIJQGCLkgiAIVY4IuSAIQpUjQi4IglDliJALgiBUOXXlOOjIkSP1PvvsU45DC4IgVC2rV6/eqrUelf16WYR8n332YdWqVeU4tCAIQtWilPqb1eviWhEEQahyRMgFQRCqHBFyQRCEKkeEXBAEocoRIRcEQahyypK1IghCJu1rOmhbsZHNnTFGN0eZc/x4Wie1lHtZQpUgQi4IZaZ9TQdXL3uNWDwBQEdnjKuXvQYgYi64QoRcEMpM24qNKRE3icUTtK3YWHYhlyeF6kCEXBDKzObOmKfXS0U1PSnU+g1Hgp2CUGZGN0c9vV4qnJ4U/NK+poMjFj3HvvMe54hFz9G+pqPQZaZuOB2dMTQDN5wg9l0tiJALQpmZc/x4opFwxmvRSJg5x48v04oMgn5SKJbgFuOGU22IkAtCmWmd1MLCGQfT0hxFAS3NURbOOLjsroGgnxSKJbiV6poqJeIjF4QKoHVSS9mFO5s5x4/P8JFDYU8KxRLc0c1ROiz2UW7XVCkRi1wQBEuCflIoViygUl1TpUQsckEQbAnySSFoC9/EXF8tZ62IkAuCUBKKKbiV6JoqJSLkgiC4Iohc7VoX3GIhQi4IQl6qqTioFpFgpyAIeZFc7cpGLHJBGMQEVbpejNTBWi+rDxIRckEoI8UUsyDdIUHnaourJljEtSIIAeG1j0ixe4QE6Q4JOldbXDXBIkIuCAHgR5TtxGz+8g2BrKkQd0j2TQkItDhIyuqDRVwrghAAfnqK24lWZyxO+5qOgl0Mft0hdm6PhTMO5qV5xxS0pkLXJlgjFrkgBIAfC9NJtIJwMfh1h5TC7SFl9cEiQi4IAeCnj4iTaAXhYsjuldIcjdAYCXHFkrWOPvxSuD0qteNjtSKuFUEIAD99RFontbDg0Q1s74rnvBeUi8GspPSSJVIqt4dUeQaHWOSCEAB+Lczrpx1YEheDF3eJJ7fHwoWwYEGgaxW8Ixa5IASEHwuzVJ37vLhLnNaUynvf3sW1q5Zy/nP3wuzZkExCSOzCciFCLghlphQuBq/uEqs1pdwzvX3Me+H/cf6fl/HQIccRvmIhrSLiZaXg375SaqxS6nml1BtKqQ1KqcuCWJggCMHhNUvEqripbcVGunvjXP/sXXz7z8u4d9LX+d7x36Htmb+W4isIDgRhkfcB39Vav6qU2g1YrZR6Wmv9RgD7FgQhALy4cOwCo929cX684mfMXreCX36xlR8ffQEoJUU8FUDBQq61/gj4qP/f/1JKvQm0ACLkglBBuHXhWAVGe3t6ueWJ2zl1w3P83y/P4j+/chYoBUgRTyUQqI9cKbUPMAn4k8V7FwEXAYwbNy7IwwqCECDZFnZdoo/bHvtPTn7rv7n9q+dw6+EzU+9JEU9lEFiEQik1FHgIuFxr/Un2+1rru7TWk7XWk0eNGhXUYQVBCJh0C7u+L86djyzi5Lf+m//79X9n71tvlCKeCiQQi1wpFcEQ8fu11suC2KcgCOXBLG5KdnXxi4dv5Kj3VvPDEy7m4B9dI0U8FUrBQq6UUsCvgDe11rcUviRBqExqZRBC66QWwrEu9jpnFl94by2LZlzJwT+4clB+18FCEBb5EcDZwGtKqbX9r12jtX4igH0LQkVQyYMQAr/BfPIJ0+aeB++tg3vvYd5ZZwW3WKEoBJG1shJQAaxFECoWP21qS0HgN5jt2+GEE+DVV2HxYjj99Jzj1cJTSbUh5ViC4IJKHYQQZMvZJ557jb8ceBg9q1/lqjOup/2z/zvj/WJPNBL8I0IuCC7w06a2mJiVl1Zl9+D9BvPk02vY/5vTGLflA/5txrUsGT0pR6RlPFvlIkIuCC4oxSAEtzM/29d0MOfBdbYiDjAsGnE/P3TTJg6cPZ3Rnf/gvG/M58X9DgVyRbpSn0oEaZolCK4odpdCL77uBY9uIJ7QtvuKhBS7evvojMXz7ov334djjmGPT7ZxzswbWD3mgIy300VaxrNVLiLkguCSYuZQewmmWg2iMGlpjtLV25ezjeW+3n4bpk6Ff/2LKy78Cat32ztnf+ki7Wd4hlAaxLUi1DRu3RnFJii3xUvzjqHTRugz9vXGG/DVr0IsBs8/z8kXnpLXdSTj2SoXsciFmqWScsO9uC2ao5GU2yT7dVf7Wr8evvY1CIfhhRfgwANp7d/GcZiEpBxWLCLkQs1SjNxwv6LnxW0xf/qBzPndOuLJAT95JKSYP/3A/PtavRqOOw6amuDZZ+Fzn0tt4zhMogJudoI94loRapagszAKybP24rZondRC2+mHZGzbdvohqW1t99X9ARxzDOy+O7z4YoaI2yEph9WBWORCzRJ0FkahFr6XYGr2tqavP/1J4KV5xwx84A9/gJNOgtGjDUt87FhXx5GUw+pALHKhZgkqNzzo4hy7/dsFZPM+CTz1FJx4IowbZwi6SxGHyiuEEqwRIRdqliCyMNJF1I5CRM+Nu8bR/fHYYzBtGuy/vxHY3GsvT8cvRSGUUDjiWhFqmkJzw61ENJ1CRc+Nu8bO4j/4T8/AtW0wcSKsWAHDh3s+frELoYRgECEXhAJwcpu0uBC9fFkubnzUzU2RnAKg6W+8wC2P3cK2CV9gxDPPwLBhbr9SDjJMovIRIReEArALmLY0RzODjRa4Se3LF5BtX9PBzu6+jPdOX/80Nz15B38adxDXnjqfZwoQcckhrw7ERy4IBVCID9lNal++/bet2JiRT37Wmidoe/J2Vu4zkfO+cT3vFBBnlba11YNY5IJQAPl8yE4WrRu3Sb79p297wSvtXPvc3Tz92cO45JSr6a2L0FJAoLVSh2kIuSit7buoFYvJkyfrVatWlfy4glBMskX76M+P4qHVHRliGAkphjbW0dkVJ6QUCYvrz41bxsRMe7z4j0uZ++I9PD7+CC6f9j3i4QiRkGLWYWN5/q0tbO6MMSwaQSno7Iq7cpPsO+9xrNRBAe8tOsnV+oRgUUqt1lpPzn5dLHKhpgnKB2zl777/5Q9yhDCe1KnApJWIe81yOXr8SD71k4Vc+j+LefiAo/jeSVeQCIVT+7/v5Q9S26b3Z3FTai9ta6sH8ZELNUuQPmArN4SfZ92GOneXZPuaDiYtWMGYRQu49H8Ws3jCcXw3TcQBknkWkK/UXnLIqwexyIWaxYsP2G+aoFc6Y3GuWLKWy5estU1fbF/TwTUPrWPO73/Oeasf5Z5JJ3H9sf+OVt7tMqd1Sw559SBCLtQsbvuIFJImqPBumZvb27k/fvLkm/zgsTuYvW4Fd33xVG48+nxQyuNRSK3bCckhrw7EtSIMStwMjHDbR8TOcl/w6IbUz1ZuiEhI0VQ/8Fo04v1yy3F/9PVx5QM3MnvdCu748qyCRDwoN0mlDOeoZUTIhUGHW9+3Wx+wneW+vSue2md235bmaAQU7OoduAF0x5O+vk/q+PE4zJ7NjA3P0/aVs7nlyLNtRXyPJiNrJZsh9eFAp/tIrnllIEIuDDrc9tB22zTLyf3w3aXrUpYoGKPW3lt0EkMa6nIGJPtN9B3dHIWeHjj9dPjd73j9imv59VdnW24bCStumzWRNdcdl9Oz/KzDx9HcVO9zFdZIv/LKQHzkwqDDSw9tNz7gOceP5/Ilay3fM1MIs/3ZQQU/o5EwV311HLS2wu9/Dz/9KQddcgkL13Tw3aXrclIY4wmdCtamf7diTfqRfuWVgVjkwqAjqB7apu/3iiVrceOFjsUTXLNsPUcses639Q0QViplRd984meY/oN/N7oX3n03XHIJYIhv0qaYz0pEi2U5S7/yykCEXBh0BJH/nO37dSvMXfGkY29yN+werePWWRN56eLJTJt7ntFH/J574IILMrazE0sNOUHHYlnOkmteGYhrRRh0+Ml/zs4T7+rts+wzHrYpqw+S7V1xfnjvSo58ehHD33odliyBb3wjZzurIcsm6a4TwLYdQKGWs+SaVwYi5MKgw2vZvZX/2A4/Iu41l3x41w7uWXItQ7d9CA8vMyb8WNA6qYVVf/vYshUAGK6T+cs30NOXDKQdgB2B5ZqvXwrP3gA7NsGwMTD1Opgws/D91gAi5MKgwk9QL9+Un0IJKUVjJJSRimjHqJ3buX/x9xm34+9cOOMHzBjzBdqyhiqnf4/n39rieJNI76+STlgp/+mHxRDc9Uvh0Ush3n8T3fGh8TOImLtAfOTCoMJPUK9Qn3Y+Elq7EvE9P9nK4t/Oo+WTf3Lu6fN5cb9D8+Zo+/VxJ7X2L+KPXmoILXpAcNcv9bWOFM/eMCDiJvGY8bqQFxFyoepJryy0E+WOzphtkUrYZ2VkoTRFQqljj9nxD5Y+cBWf2vkx58y8gZfHTQDIe1Ny8nFHI2GG1Ict32tuivhbdLEEd8cmb68LGQQi5EqpXyul/qmUej2I/QmCW7KzS5y4YslaftD+Ws7rhQYv/ZTeg5HhcsaXxvLZHZtZcv88du/ZxZnf/DGrxxzg+Ll0K9wqa8QkFk/QZfMk4PsrF0twh43x9rqQQVA+8t8APwXuCWh/guAKL/5tDdz38gc8vv6j1HCFoz8/qqBMlEhI+S69B1j7+/9h2dJrSSR6OeOMG3nzU/vl/UzKCl+/lNYXbuCU8CY61Ahu7pvJ8uSUjG3tvlVnLM4RDr53W4aN6XerWLxeCFOvy/SRA0SixutCXgKxyLXWLwIfB7EvQfCCHx/x9q54yud838sfFGaRK/+l9//rn+/ym3uvoru3j5lnLOSTzx2Y9zOpTJM0X7VCMya0lUWRu5keWun4+emhlaysv5R3G2azpOvfmBZa6a0/ytTrDIFNJwjBnTATpt0Bw8YCyvh72h0S6HRJYKPelFL7AI9prQ+yef8i4CKAcePGHfq3v/0tkOMKtY056iybUuR7F3KMgz96m3uXXktXpJHZ3/wx7w93F3hsioSorwvzWOI/GBPamvP+puRIpvTeYfnZ6aGVLIrcTZPqTb3WpeuZF7+Q5ckp7kfMSZpg2bAb9VayYKfW+i6t9WSt9eRRo0aV6rBCtbF+Kdx6EMxvNv7Okw1hV1l4xpfGuiqr90s0EvYt4l/oeJP7F3+fTxqGMHP2ItciDoZfvTMWZ7TKFXGA0Wqb7Wfn1i3NEHGAJtXL3Drjd+z66WbCTLjidZjfafwtIl52JI+8VqgGK8pHLnG+ykK7YplC8Zt3/qUPXuPXDy7gH0OHc+Y3f8xHu/szajbrkYyxEPPNeoTtZ0aHrEXeFH/T9x7UHFOhdEj6YS3gNvfXozUcOH5S29YvpfWF43mpewbvffoqXvr61pTo/Kj1YG6dNTGjR3gkXJ5UQ4C579/Lb393DU3NvTR+q44vNvtvWHVz30y6dGZL2hgNtPXZ35w3J61FfrMegcJ4ugm0v3i5z6caIqj0w98CfwTGK6U2KaUuyPcZoYS4EchiFXp4wTa17UPrdaxfSt8j/ydjzX2P/J+MbVsntaR6hK+9/jjavnGI6+UEKflXv/Mb/uPBJYRGKNS3mhg9rNNVcNKO5ckpzItfyKbkSJLaCA5GZ/yU229cSItNbvnd9WfliH+XrufmvplojN9VYF0S1y+F9oszz6f2i0XMi0RQWStnaK330lpHtNZjtNa/CmK/QkC4yf2thMo6pxS2tJuKWQC06cGrqUt0Z2xWl+im60n7DIpVf3OfXBWUS+aEjS9x0bIHUZ8OwbeGwBDjskv3T/theXIKU3rv4NDwUrjiddoTR6SCv9k3oWgkzMSTLuLmyMUp8d+UHJkR6IQAuyQ+eRUks9oDJOPG60LgiI+8FnCT+1sJlXVWucQm/TeV9sQRqV4qoxusA36Nsb9bvt6+poP7X/4gyBUDRjbI3LqljFZb2axHZuRzT3/jBW557BZoCcOZTdCYKbFOwUk3RMKK66cdmNNjRjPQrKslw899Eccu+zKxtEKh9OZZdkOkPXdJjNncMO1eFwpChLwWmHqd8VibbiGFIpm5v8Uq9PCCGdBc9m/W7+/YlPHobxfw25WsZ7cFw0EnQIXh0HN5JTmeL66+mXcacsW2ELJT+sYoI5/70MRfmLb+RfZ4dCs9+zSyY1YzIxq7cj5vF5x0ujmYpAv0EYuey3GJmCKenlKYLzhs1RpX+otXPiLktUJ2P5Hsnyulsm7CzP7sGuubyuZ/DKzv5r6ZOXnRca0YqnoG/CI6gX7lV0wiRJ0yKjBNsSVOwWJul9J39ponCD3RDZ8J0zgrQqiulx4dpkENCKTpn87G7uZgrjcSUrSdfkhGJoldjxmv4+0C6y8eHW5tfUeHe9uP4AoR8lrg2RsgkSk2JHqN100r2Py70BRFpzTH1HsfGpayThgVfNnHcbipjH5i4NF/eXIKxOm3XLexWY9gL7XN8p5VR2YZvemfXt5bmJBb5nP/sYfQUz3wuTo4PQp1inr6+FgPZUuyMbVWu6cCp3zv5b1TcqKw7Ws6bHue+xkcEUh/8RNvgkcuyTzvwvXG60LgiJDXAm793xNmFpZbnpUHntjxAWEzDxwyxVn3W6ZWueION5U5iQ5WPnwnl7OY0Wor2/XQfuE2ZCzkIUSZzz/txr2R49757x54rgcOqIMZUUhLd2xmF1/ovQsFnHn4OJbb+OvzFfukD1gGw3q2+tZmSmE6JcsRD8owEFwhQl4LlMr/nZb58lvi3EEvT8RhDzPzxSqI2f/635ddw5cfGJImLtY3ldbwS3w9/EvqdQ8AI9TO1Htj1FZPXf0UmvcaZtNhIdJW7o3bIndyO3dmbJ9y79ADL/TAi73ogyOo1kYIZZrO6f7wyXsP5z4bIXdT7JPuMrHLKNGQShtsndRC+5oO5jy4jnjC+CV1dMaY8+C61PuBU6hhILhGhLwWKJX/O83CHwK8SoJj2MVTOz5gVJ6s7E/prRkFKCY51uMLN6RE3AqljBatblqMm9tY+cyt3BuhtO1vqr+bIdTx2+4p0Kv58R9+zm5/7KV7UhORk8KEs0S8V9el/OGNkVDGd8zGyvef7U9Pd5nYZZpA5u9zwaMbUiJuEk9oFjy6QSo3qxwR8lqgGI+5Vr7wNMt/OhEeQXEqXRzFLp6liT0dyhaSKN5tmJ1yYVyzTBOLJ1Mug47OGJcvWcv0hg+zDV1L3Iq5SbbP3M69YRKllyvDi1msD+MLT29gt1d38F9fOIlpJ73K8NCunO17CaduErE8bW+tfP/ZTwy7evrYd97jqVa8D63uyMlaSXcN/fORUXyl93QIkesu6io8e0coL4F1P/TC5MmT9apVq0p+XCEgsnuigGHhHzIb1j2Q8fpz9DGNLsYQ4lmaGGMh5tmi26vr+F78IstA4Mr6Sy27/qXTpwcyVLygNezb8wDTQyu5JfLzvPtIJmHJI4dzxvqn+MVhM1h41Hm813im5Q3E3HcxiEbCnHZoC8+/tSVlmVt1OuzVdWh0TubMvPiF3HHjwqKsTQiWsnc/FKqcVN+MYUaet1UV6NtPpfWUNjiGOp6iiY9IciS7eD8re8TKcq5XfcyPWM8oublvJj3aeiIOGMIUwruIm9wT+TGLInfnvxEkNbH2JGesf4o7vjyLhUed5+0RoJ/0/uAr6y91LNm32zYWT/D8W1t4ad4xqQpNK9dQverLEHEwnkTm1UvZfLUz+IW80hv3VPr6IKsPiwPm+1e8TnqO3BHU8QxD2I7mSHbxNgNiYqd9ezAQxEwXsBsjv6aerMIXDUlt9OK+pu9CNuuRlvtMahxvAkrBV0IbcgQwh4RGPxRjyGs7afvK2dxy5NmpL/KxHmp9bFSGSE8PreTVhou4PXInY0JbCSkch0NMD63kpvq7bbc1A55mW998rqF09qKw6lKh/AxuH7mPtqi+j+PH/1yq9RW6Tqs+LHb0r78ruidNsY9SLx9GmOcZwrF08VW6eIYmDsBeVIGUSKW7CIbSnbOdUrBdD00NVEhqaIv8Isf67KOOxYmjmBpaS4va6seAhj4ND8ZQG/v44dEX8KvDTs14e0HfOZbHDiudCqhmf6d0mlQv19fdw1wy/djzI/cQxT633Ax+mkHLfz4yij3Z4uorKZmLWfUMbou80EZQbqxlN10D7fbjdX1O68m31nzrdPp8PkvcYv03x2fldNo7RId5gSY0cBRdrMO+p7dS8JPIXVxfd09+CxnDgjeFf3lyCrvILYSpV31MDa21rKZ0RVzDki7Y2Mf2E0fkiLh57Dnxf6dP515aTaqXWyI/Z37E+TsNVzszLO+fRO7KeEJJZ7TallNC3zqphT1n3Jg7ki1cb7RmSKca5mJWw1NrmRncwc75zVjXuyljuokTVgG9cD3UD4XY9gGL1racfKzhYrALDE67A5Zd5H59TvsB+/dMi/vWg+zXaZeeaH5+wR6gvfmd9+1+gPl1v+bs8DMZWSZaw9sqwVS62IVmhR7CF5W9Ze4l+yR9zNm7DbMts1uSWrFZj7ANmGoNMepzhbZXw2+74P0EyWmNqEnGTWo7Q5kfPycnMGt3fK/fKR+bGcmfT3nROn3Q6gkMqqtIx+m8r+R1Fwm7YOfgFnIn8bridX+fTScSdXA59Iux0xrA/fr87EeFDQG2KwjKhynydk2sbFHMr7uMC3vvsxXMd3WSo1WMLSh+rxs40kbMvYheUiv267kfsM9uSWhFCG27T60NP3ez2kVY9V8b3Roe6IJNCXRrFDUhYvmZx5KHMzW0ltFqK4rgxNoOrWFB5HLm/2ABYFRtrn38Li7svY/RoW10R/ek6cQbghG8ck2YKuQaHoTUZtZKIRO/3bRvjccMsbTC9Ds6lcd7WZ/T0AU7kdYJUm4UP2MSdmzy2Y9cMzeyxLEEfj8V4r91lL3QfF3t4gX6bLfNdtHYsUvXpwKjo9VWkhY2SljZizgY4jsitJMEYSMwGtNw7y7oSJD8Rq6Ip3/mnPAzKZdIPhEPwn5KotjeZTw5tK8xWhfMjfcHT9E0xT7KGbThi3IOHSlGe+VB6KoZ3EI+YWZaOpwxRcX1I5nbAJBOOIux3X6GjfG2voIDUmaHag9E9/B9wTTF/k53056O24xTIVYSZRwhTqSLpyzE/GM9NGMSzsd6qK0IDlE9LIoMZHaElH/BrFd9dO2qp/e/etD/SLL19E+j/leuiKfjxQJ32jZ7zVpbf4+w0iyq/xWsX0rbio1czuIcl1Bdohse/nZholXOoSNO148fKmESVhEY3EIO/id+W1nLVpjiayfG+axut+tzux5HtP0ThB1+L5joHkb/kTzsRYgXaOJzhJhGF48x0DNda3gseXja1pou3Wi7LwU5QubbvbEzSfM9W4hs60XPakJ9zuPvrQCsujfafY8oPfDsDUYbA7uUw/QnMz+iVc6hI4U8VVtRCZOwisDgF3K/ZFvL0eH2EX8nMXZjdbt51MvYTwFo+0yRHGLbLS8kK4sx42cwelGn9aN2Mow/RYjnGcIEQpxKjIf6xVwpOCP8Am2RX2RkcRSdT5Lwmy5UZxJ1ZhOh/esYEdpZsM+7aOGoHZsY3Ry1zZ/PwI9oBW0Ve6GQp2orKmESVhEY3MHOoMkX8PETEPITlXcTiA2KYWNh/+OMqs0dm9iUHMGzyYn9Qb1tdDKEeuIM6be+CxG7HWi+Thd/IsE9RJmNvSsjOwiqNezUDewWyv8U4EhnEv5rl+EbP7MJxlZBqcWwsbQftYJdD1/GGeppV71omPFL4/xyc86WI3OkWMHVKg+e1mbWSinxe7L7ObEeuxJW/ZrgxgO7p5PduK73bJYnp1j28yiUnWim0cUfSHA3jZyPdaDTSshNfN9MtiXgni4j1fDsITDavTvFPH6xM1VyjgssUyfwYs9+LKq/O6NoyDEqYtMbx/ac9WPEQDCFck7r8kqVpzOKkBcbN4JsdbJ7ySWH4om4CrnOFTezQZL4a06Vjy40p9LFUyrBz3QjF9uIuVccUxm39It4Eji7Cfb05hNP6pz2447rAG+i7yTKvbqOf+lGRoSsiobsZgcxMKUpG6/WqdU5Ga7v752QNifWrWAW22ouVyplANgJeRU8N1YJ+XxvduX40T2sZxta+R/XLy2eJe6h4McUrEKaUznRhOIRmpipY1yiuunRmitoKHi/tsL59wTc22VEjM5tglHeA5teDHG3Aq618T+9WY9kL7VtIK89i3rVx3Cbyk/Hc8UuXuLFX2x3TmaPFoQB/3w+0Sy2H3sQDryQYGdQ5AsI2UXLwX1U/tkbKIc7pRw0oniQKKfpOq5UPdzoIgPGF5sThk+8DksRt0v9yyboQGiXruey+MVcHr8YMEbYBf7wnK8Gwg1ez0k3YlyMlMNBljeejQh5UORLk7I7gWPb3UXl1y91EeDsVxOvKYYVSj2KxUQ5U0f4vurhOrrRQd7IPuyDe3ZBo4Jzh8CI3N+bU+pfkJjHMIK2jcR0PbdH7kx1R8y3ju0MdV04lWKfKe6MCMc+PB6tZDdi7DXlMF8PokGYN56NuFaCIt8UHqe5mXaPeulT5109vGsjTRKs3TVVSB2K/6KRRg0/VL3ENNxMA8pPpWo67/cZZfe7heBbTbB7Zdg0SsEQ3c1QD8vp0yHmx8/h0NBfOCf8jPsbz/sr4dBzUxlJtkFMpw6dTu0fQpFcH7mb/G8vE63yrc8pb3wQuVeqR8irIUDh5Hubeh08ckmm7zBc72xlZETXXVqiVgJePwR6c8ePVQthFHfRSIOGn6heurXmdhoJ+RXzd/pgcRfsETICm7tVhoibeH0CiGkjTfP08IvePqsTsOpXxs1/xl3W524+IbRquIaCyecb/1z9G+M4Kmxkybi9Zt36sfOtb5DmjWdTHUJe6r7dfsl3s8lXSZOOlx7g+ahiETcJofgpjUS14j9VL90afk4jYa9i/pc4LI3BqBCc1QRDKkvE/TBU9XBVJHcikGtiH9tfT/mE0M56BmOfZkBVJ4xUx3GHB3vN5luf05PwIKI6zuJqKKtdvxTaL870xbVfnNl7PP0xE4yfrb6DK3947aFQtNHAD3Q9d6s459FNnxef+RtxWBKDT4fgnCGDQsQBUPmHRefF7npyE3i0qmwu1TWbb31Bl/hXKNVxJlfD49GTV1kL9ZNXGf92+x3MG4JgiULxQxr5oW7gXhVnNjHibsT8tTg8GIOWsCHi0RJX7hQRha/elrlYnaP5hNAu0FiqroVuehkFWeJfoVSHa6XQx6NS+Nftgouxj42Tzm2+uNUNQcjhBzQQ1fA91UOvjrGEKA12cramF5Z3wz5hOKMJ6ssn4kEOlQgcq+vJKfDo5PIM2qVhd6xpdxh/nK7vIPLGKzxGVx1CbjfBxs3jUSX413d8ODBmyy6Kb54ogyTbpBR8lwYateI7qptW3cUymohmi/krvfBEN3wmDLOaIOJdRYMSX63hnsTXODn0MsNDO4OxooPC6XpKF0LzPLUbNmK6Twq5Zq1wctV46Wrqh0rQkDwE4lpRSp2glNqolPqrUmpeEPvMoJDHo1L56sy0PzsSvdCwm/V3cDulXsjhEuq5WzeyggQn0cXOdDfLH3sMER9fB98cEHGvhTVxQoEV41zfdz7dqrF4Ih4dnr+OYPIFRtMsr9eT2/N0x6bB1bWwCmJ0BVvkSqkw8DPgWGAT8IpSarnW+o1C952B38ejUp0AJ96Um16YTWw7XPVe7utBZqjUIBdQTwPwLbo5gS4ep4lhL/bC8z1wQB3MiEJ4QDo99TjRECEZmDvk3YbZxXetOLUqViE4+Rbj33bXk50bwe15arpPgiyFd3LVFNvtUQUxuiAs8sOAv2qt39Va9wKLgVMC2G8wlKqX8oSZcMrPnPuF2x1zEFjiZei9lsFZ1LOEKH8iwbGf7GT7H3tgQgROyxRxrwRZ2an6pxYVVcdj253PwXw9dawqIZddBPOHuTtP/bpP8pXR2wU19z+u+JWb5ezH7pIghLwFSP8f3tT/WgZKqYuUUquUUqu2bNkSwGFtyD4h9j+udOlHZhrWjF+6P+b6pRT50i4JlRDA+wYRln3SwLohml9MD0Fro/uWhIOF6B7O53Y+t4ul1e3yLq3C/twnbsro7Vw1bz9VfLdHFaQwlizYqbW+C7gLjDa2RTmIVVBi3QNGRZlTGbKb/XrpteymxDij/F4Iimm7N/DK1hAHfT5cGXeXUtOzAx693P79Q891/rxfd0EhPb3dltFbuWqWXWS9zyDdHl5aBpSJIIS8A0h/lhvT/1rpsTsh3n7Kfx9jq5tD+8WGSJj+cKsotpN/0Kq5fRBEouJrByaMdB6SPKhJJiDpUMn79lPG+Wd3btqlyToxbGxhwlaID7pUlZsV3vo2CNfKK8D+Sql9lVL1wDeB5QHs1zvFCEpY3RyS8dygppfHuSevCkZwVSj3UbPQmZ7C4MbJh7x+KfT8y/2+VJhA3IKF+KCrwO1RCgq2yLXWfUqp7wArgDDwa631hoJX5odi3J293AR2bHI3EiuoXHGdtH7SKIa1L5SFohQQ2XX/s2oj4YSZHVNoXrWbnHO76ypot0eFF/7YMbhGvRVjHp+XQcfR4dAXcz5+oIOT+6/w7BNu/VJ4+NvOaWhCxdOjwyxOHM3U0Fpa1NaABd1ilOD8ZhwDm3aj4dLfP/Xn/q41JwEt1ZzNKpjnWTszO4O+o1r954YimT5yMP7D66I2ZfhpswbzXSzZuJ2laQ7TNYO6AQ5gKNdg4VpGa7gsfjHLk1MAeC/o/HOr+ZduZmXmM0SKIXzFnuFZ6uMUQO3M7Aw6KOHUptNykLIF6SeHUyN+K9zeaOMx1/M8Haer2yAiXnpui9zJXL2UP4UnBzvhz9ZtYXFeZm+bz9Xod2iDkwFWqoKcKij8sWPwCXkQWJ1UVndkKx+jpUirgUwBy0b8Dpj+fVfi7/Jq96jkIuL+KMS/rfoLh8aorYzRK4JbVHaGiWMGlcodBuHGEPEqfH6nEAWdmVLFvcuro41tMbCrJCtkxt/U67BWSD2Q0WIWNriZqxmKGPu0iswXQDmEudyVn+VAqaC+d4C/vOwGU45l99pw1aXj5lz0Knz5epmUKjOl2Mcp4hDo2hRyJ7EupEHOhJnYXnTpVsqEmfkLM6LDofXOAVdRelVbdLgh8lVErVr1Ffe9H7sy8+d81nP6++nXR8oQyfqCfoTPzRSiUvQUL+ZxijwEujZdK05iXaifbNjY/I9n65fCq/fkbuMU9c/2/fdfVLo/sFlpelFNVHSP8KBZ/RujaZYpyvmsffO8zXZ/6ERugN1vcoEbl0apCnKKdZwiD4GuTYvcSawLbZDjZqLKw9+2ztfVCfdFRRNm0n7UCg5ILKYjOdLdZwRLakbEwTjHXLdNVsY2tx5kXcSWXjWdPubNCT9TfgYDRQ6k1qaQO4l1oSeVlRtEhYxG/POHGX875eI6/cdmXQRrH7+LWDzBzX0z6dL17tZXALXo5x50qHAev7hK+7v/P3zHh/ZFbF6EyM69AIN/HFuROyjWpmvFqZIsiEox8/Fs/dL8PcqzseuvDDmR/bn6Tj4O9Rq5xnGYW7eU0Wobm/UI3tWf5ivhDb5dLlbuBjN4V1MW7GDj0HP701RtCEcg2eeudgG8CVE5p/yUm6AnJmVRm0KeT6yD8pM9e4M3EQ/XD/RXzk7FqsttiNWkeplbt5TlvVPSXjWsqAcTX+U9vRdnhp8jrJIoFYZ9psCmP2fuJ1xvuUYnse7TIcIkodi9tQWfKOOcHr4fvL/SeAJUYUPET76l36dt41bxcr56FaIqztMumCJ3UBx8lZ2VhJcqzuhwY8qQx9a2Sa24PP4fLIrcTZMauAi7dD3z4heyPDmFluYoL807xngj29rv3eWp94vWsCByOXMjS2js+gjIbPktFrsb0twWQeOmCtHPkyIY52j9EP9CVAWVk5VO7VR2VhJuiieyS5rtqkNt6G7ak2tiv6OJzIsy3Vrf3BnLrN5TYUwh0bGPPVnVcRXhqr7/j2hfT8ocT2rjnx16JC1qq6f1u6GULQKSuhSzKLR9dlMG2YKvYPL5Dha1sh9ekm0J1g/11rwtEjUMjUIsyCK7F2qZ2gx2loqp1xmuCzusgjp2PsfocMsgbNOJN7An1uI5Wm0D4FtD/5yZpZDWtc6LXZgEduoGovRkvB5ShohP6b2DDl14Bo3WkNCGqG5KjuSy+MVcFr+YHu2iiKoAkhpWJg8koVVOYFfrgIO9XR9DKN+Q5PPJdF5pY1CK1dQrU+SzhdYuwOhFxP1O/smmVPngNYgIeTEx53hGhw+8Fh1ujIKbv8MQ+mdvcJeKdeJN9heBjfhv1iOIRsLMjSyxzVIIkStQdoLVp+toZqfle+ZN4+a+mfRq9w96pkCaf0zh/kzPA+zX8wBTeu8AjEBuPYmiZc6YIj459DZhpVPWv9awLTmUy+IX296kElqxLTnU29riu4whEHaocH9VZdZOzZS/Q2YPFOWosCHi5lBlGMhwWvZv1gFGN5XFYJx7fjsaWmGOQ3Sbrii4Qnzk5cKpZSZ4C4pY7KtL13Nz5GImnnQRrY8ciJNPNtuvHdeKEBBWuZ+x84FvSo5Mie700Epuj9xp6woxXTHbMcRvD7WTzXokzyYnMjW0ltFqK5v1SG7uM75ztv8/SLSGDkbyMzWbS5IPMCaU+3RjfrfpoZW2sQiA+ZF7GK6sb3T+cPClhyKZtQjpLjq3E6iyJ0qZBT4bHh6w2M3YjQhuRVA7bWyrhaADP9l+0P2PG6i4UyHPvcmdgpbZ76UHVk1W1l9qKYp9OsSV8W8DuQJttd+YrmdEqHBxtOsTltSK5a0baJ3UYhucTmrFoeGldMeTHJv4Q0aaZ9FuNuakJy+dMs1zx03Pe7N5VrbB8MHLuV00K6wndy0jwc5KI+hUrPSUyezBEj4GTDgFFpUyBDmETolZuoiD4WJxyqRZWX9pjvBlH7NJ9RLFWhy1NhoTJFHUKZc5zxZ8ooYaIg50RfekKfZRzjab9Qg6u+LcOmsibSvqmdKZ+V2tvkuK6HDo3ektQ8RslgbeOmWa506+cyi9ZiK77YNVK+QAS8mF4iA+8mLgpstZMSu9Hr3cUryTBBe0C6HZr+d+pvTekSPiAMuTU5gXv5BNyZEktWJTcmSG1T66wOyWDj2S/Xru58r4t11VtSZtvnOi/432NR1ct+u0nH116Xpu7ptJqP8u89K8Y7ht1kSikQEfs/13UXDVe0acRLm81CJDcpulucU8d5zOIacAo1PvlVrI9a5ixCIPmny9lU2KmYoVt56irjTs2/OA4fYoUEg36xF5t1menJJVrJT++ZGu1rCLRkI6mWPZmy4NgJiuT1nuPURoIJ5j3dulFO4R2gm3HsT0HR9ysjIKnYynjSRJQkQx0jjpg6uXGWLcOqmFlg8fY+yrbXxKbyGpjO1zMAXV/H93M36vriHzPJkw011tQfq5s/9xsOpXudtMviAzIJqNk1gH1ZO7SmdiVjpikQeN2za45UjF6hczq94syf50vz6d/5TIFlI/WK3B6kkhQpzfJY7MsOx/lziSuXVLebdhNrdF7mREaKcxiEFBvYWIg30FqgLY8SEhoE4lUar/77Sfx4S2sihyN8cm/kDbio28svwXHPTqtezJFkIK6qxEPPumPGGmkf2RnsFkRezj3Cc4q0ymcH3/vizOnewe4iZ2r5vYirVNfrpXitzKtZaRYGfQ2FZzWgy7LRYL9rDslaFRHJBYTCyeYHpoZU7QbnlyimVmRo8Os4sozeyy9Yn7IXsNzepfDFU9OdtlZ8QEFVhM4s2SMdfxUsOl1oVPKmz83t1Ymk4BSauAtxdL1ukcnHGXtyHHZn66kyXvFqnsLBgJdpaKShgXdeh5lo/W7+49i4YPQsTiCZYnp/BYfEqO79iqAVdQwp1Ntuvl3YbZltuZOepgrKtQEdc6rQrVQxXnaLWN1vBLjLYpwEIn3d+sp15n5HhbYeXi8NL/x+4cjDT1Vw6ndTVMd/tZ9QMxs5/mNxfuCqnlXitFRlwrQVMJvZVPvoV39v4mfYTQGvoI8cfhp3Lyu6fSGRvIPbYLAC5PTmFK7x2OwcxisNmm4CbdH+8mSGr3vdIZrbaS8Hj6dzKEG+t+aZ/R4+VmPWGmvZvFzX6cAupW52Ao0h87sclISV+XWbAz9TqjkjQoV0iRW7nWMiLkQRO079vHnL/2NR2c/O6pfLb7PvbteYDPdt/H7M2nE4t7T0MsJVZ+815dR5Ru3m2Yzcr6S9muhzruo0vXc2/ia46ZOUoZwc86lbTdLvv1GA0Mqa+zfRroQxkNyOY3w037Gn/y/Z+deFP+ISR+5spanYMNu9n/Quws4kLGHlpRCUbOIEVcK8XAy2Owk+/TbQZMFm0rNuaIdjXMhMh263QyhCHEUgVBY9RWenUdca2IZFWdam1Uis6Pn8Py5BSmhta6yoqxsq5zi6EU0clnOfbxDqMHqiHT+5g4/Z/ZtTYF40Zgtx83Y8Oyz8H5zbZrt7WIi1HrAJK1UgREyMtJPqH2Oedvc6fLApIKJN1vvrL+UoZnVXXWqz4SOld9lYKuZCOrdz8WOmM8m5zIOeoZXx0Tcz/TP03eoZul42Gc/s+sinLsioAKmStru3aHjJRixHtKNXuzxhDXSjnJ9+jq0yIa3ZzdGa86sfOHh2yeL0aHtjHn+PFEI2GmhtYG2/Z2x6bCXABurVjHMWz4nytr5daw65jo9BlxhVQkYpGXmnRXSr4qOpcWUfuaDtpWbGRzZ4zRzVGO/vwolvz5Q+JpUb8QgHIXCKwU7IqGEoQsc7e7o3umyu1HP7It53032PYjHzbGELwnr/LWAjb9827IJ/imO8JrMZkft4a4QqoGEfJS4rYrnXnRu7hg29d0cPWy11I+8Y7OGEv+/GGOzCXBt6N8IN97oCthKTJZ7Pq1/C5xJKeHX8zJdW+iB+Y30zpsDET3sBFc1f/e9pzmYl3RPXmy+xBOTD6XGdRM/52feBO0X5zZeTAfXqxYp2Ekhc6V9ePWEFdIVSBCXkryPTZD5kXv4oK1CmzGAzS7swtwxiijypE4RRdzM/g5r34pe+rMnPbVyc+lgqLb9RB2C3UPCPeOD43KR6dWrxY0AaeBcwDa/NvRMjdNem0UCR0y270YWt28IbedrAiskIZUdpYSxxmeytej677zHi9qRopdO9r0astiEo2EWTjjYMC4aXV0xnK6dNtWWhY6YzIfbmeymn2+zbbC+dYi/UgEG6SysxKw9Xn7L1Ee3Rylo4AslYa6ED199m1g7QKO6dWWxSKsFLF4grYVG5lz/PjUAOnsmMDobpu1xLYb3QeLhZuZrGBY1+ntYfOlkYq1LXhEslZKiY8sgPY1HRyx6Dn2nfc4Ryx6jvY1HRnvm1kaGbsMKSJhdykbvQ4iDu6qLYPmiM8MJxoJk+h/WuzojHH1stdS3711UgsvzTuG9xadxEvzjkGVq2LQMhPEjjwVlYJQAAUJuVLqdKXUBqVUUimVY+4LWXis+jQDmR2dMTS5ggaGqC2ccTAtzVEU0NIcpe30Q2j7xiG0uEhDzOcYsKq2DKL7oR1nHT6O97fFcvz+pmVuSbnS5Kz+P/N1N0xHeowIAVGoa+V1YAbwiwDWUht4eGy2CmSagmam2YEh5uk/p7++z7zHC1puKZtoAUzeezj3v/yB5Xu2hU7lTJNzVdBjM3tTeowIAVGQkGut3wRQgVZeCCZ2wuWlcjOsVMpF4RenARFBM3/5Blu/v2OhU6X4le06CK57oDhDRASBEgY7lVIXARcBjBs3rlSHrWp8CVoWhYp4qemMxZk//cCM3HgwslfmHD++jCvzgNVNZdzhkokiFI28Qq6UegbY0+Kt72utH3F7IK31XcBdYKQful5hDTPn+PEFC1pzNJLRurYauGLJWqKREEoZDazCSnHaodbuo6qhUp4YhEFJXiHXWn+tFAsRcjGFKz3Vbs7x42md1JKTgme+nk1vX2W3rrVCA13xgWyahNY8tLqDyXsPr24xF4QiIXnkFYaVQJv50+nbZJflX73sNYAMoWtf05EhiH4YUh9mV2/5bwZWQV5BEAwKTT88VSm1Cfgy8LhSakUwy6pN3KQbgnM2S/Z2hRIJh3Ly1MtFNbfnFYRiUpCQa60f1lqP0Vo3aK0/rbU+PqiF1SJuBdptNksQwrcjFk/lqYOnEZcZBJHXNFja8wpC0IhrpYJwK9Bus1kKLd8395Gep96+poMFj25ge5e3AKq5Nr/rqaqsFUEoMVKiXyasSu/tLM7s14/+/KgcC9dK6OYcP951qb4VduLZ7dHvbu7Hqp2AWxbOOFj844Jggwh5GbDzhR/9+VE5Qpctpu1rOnhodUdGnaACy/S81kktDKn399DVFAlZiqeV+8cOs2XAaYe20LZiI5cvWUu3jwHQLf1PBYIgWCOulTJg5wt//q0tLJxxsGNaod1g5eff2mJ5rB0+c8hj/VZ3dhaNW9dIWCneWfj1nAwbrwUE4lIRhPyIkJcBJ1+4Xd8UN5+1wq+fXAOXL1mb0SXEy37O+NJYwJsFDxAJK4bU17EjFnfMjzdxm08vCIMZEfIyUEjpvdfPHv35Udxn04TKDV4t6LBSnPGlsfyo1RgG4TVzpu0bh7gWYrf59IIw2BEfeRmwCvqluxCcepBbfVZhiJhVv3I7l0vQtDRHeX/RSbyz8OspEQdvKYNefeFu0zUFYbAjQl5iTFdALJ4g3N81sqU5mgos5isKSu8/DuS4PrILiIpRRGOVMXP050dZ3nzcZqr48YUH0R1SEAYDIuQlJF2kweghYgpYel+VfFamOSGnpTma4/rI3jboIppoJMyZh4/LGGRx2qEtPLS6w/Lm0zqphdMObSGUpf5NkRB7NEVS+/CTXug2XVMQBjviIy8y6cG4kEVv8OweIl6sTDfbzjl+PFcsWeva190cjbCrt494YuATptXfYhNMPGLRc443nyWvfEgyawHd8SQ3zphQkC87iO6QgjAYEIu8iGS7Sex6g6cLrxcr0822rZNaXIt4CJg//cDUmDjTWr511kTeX3QSc44fT9uKjTnuE6cbStuKjRk3BZMkxhCJQrAacyeFQ0ItIhZ5EXGbepcuvF6sTLfbtrhMQTTFde31x+WIoVOGiFMmjZO/Oog+6fnSNQWhFhCLvIi4CbplC6+Tldm+poOJC55in3mPs8+8x1nw6Aa+MG5YKmhqN4Dh6M+Pcr1mO3F18t07ZeGIv1oQio9Y5EXEzlINK0VSa9sCFisrs31NB3N+t454mrN5e1ecl975OPWz3QCGIFIQ8xUxgfUADDAKi6zYoylS8LoEQRAhLyp2rg8/fty2FRszRNwOqwEMXtLx7MQ1XyGSnYujdVILq/72cU5RUiSsuH7aga7XJQiCPeJaKSJBBuO8iLFV21s3WImrWZzU0Rlz1XHRih+1HsxtsyZm/B68VHAKguCMWORFJqhgnJeeKdnC7TYFMZ7QqZTB9OKk9IZX+VIR7ZCgpCAUD7HIq4Q5x48nkl1VY4GVldw6qYUzDx/nakpPejGPXafFluYoL807RoRZECoEscg9Uq5ue+Yx5i/fkMos2aMpwkkT9uL5t7bkXc+PWg9m8t7DaVuxkY7OGGGL4iQT088uJfCCUB0obXMxF5PJkyfrVatWlfy4hZLtagD/wcti4vZmY/q+rVDYu3NMi1wQhNKilFqttZ6c/bq4VjxQDd328jXdSsfJ527eAPJNLBIEofyIkHugGlwNdjeby5eszWlzaxYSWWFa8VICLwiVj/jIPVDIQIigsXOfON1U0svqwb73Cwz45CXbRBAqHxFyD1RKtz0/fU9MYvEECx7dwM7uPtttzF7nhQR2ZQSbIJQOca14oFJcDV77nmSzvStuWyVq3pi8+NqzKeSzgiB4Ryxyj5TT1WBauXYWd3bfEz9Dl80bk1OP8Xzf3+lGI1a5IASPWORVQvZ0ISvS+568NO8Ybps10TLrxAk/Ay7cblNJQWFBGEyIkFcJ+Xqb21V0WrmC7Bpjpb9eyBg1GcEmCKVFXCtVgpM169T3xM4VNOfBdRmTe7IbZhUS2K2UoLAg1Aoi5FVCkFWW+fqHu92mkP0LghAcUqJfJVRyewBJNRSE0mBXoi8WeZVQqVauU057udcmCLVCQRa5UqoNmAb0Au8A52mtO/N9rlYt8qAs12JZwH72a9d4SxprCULwFKtp1tPAQVrrCcBfgKsL3N+gJagimWIV2/jdr6QaCkL5KUjItdZPaa3NWu+XgTGFL2lwElTnxGJ1YPS732FR+xmfgiCUhiDzyM8Hngxwf4OKoCzXYlnAfvbbvqaDXb25PVsiISWphoJQQvIKuVLqGaXU6xZ/Tknb5vtAH3C/w34uUkqtUkqt2rJlSzCrryKCKpIpVrGNn/22rdiYkYtuMrSxTgKdglBC8gq51vprWuuDLP48AqCUOhc4GThTO0ROtdZ3aa0na60njxo1KrAvUC0ENaShWMMe/OzXzlrv7IoXtBZBELxRUPqhUuoEYC7wVa11VzBLGpyYFuqCRzewvV/oGurcebays0lOO7TF1ZxOP+vzkrVSSf3ZBaGWKTSP/KdAA/C0MqbNvKy1/nbBq6oC/KYAdseTqX93xuJ5c66t8rQfWt1RlEIgr50dpRRfECqDgoRca/3ZoBZSTfgtgvHT3rWSW8JWapGSINQaUtnpA7/i6iczpNLztGUUnCCUH2lj6wO/4uonM0RawgqCkA8Rch/4FVc/mSHFylIRBGHwIELuA7/i6mfmZ6XMCRUEoXKRNrY+kdatgiCUGmljGzAS5BMEoVIQ14ogCEKVI0IuCIJQ5YiQC4IgVDki5IIgCFWOCLkgCEKVU5b0Q6XUFuBvRdr9SGBrkfZdDKppvbLW4iBrLQ6Dca17a61z+oCXRciLiVJqlVWeZaVSTeuVtRYHWWtxqKW1imtFEAShyhEhFwRBqHIGo5DfVe4FeKSa1itrLQ6y1uJQM2sddD5yQRCEWmMwWuSCIAg1hQi5IAhClTMohVwp9UOl1Hql1Fql1FNKqdHlXpMdSqk2pdRb/et9WCnVXO412aGUOl0ptUEplVRKVWRal1LqBKXURqXUX5VS88q9HieUUr9WSv1TKfV6udfihFJqrFLqeaXUG/3//5eVe012KKUalVJ/Vkqt61/rgnKvKR9KqbBSao1S6jG/+xiUQg60aa0naK0nAo8B15V5PU48DRyktZ4A/AW4uszrceJ1YAbwYrkXYoVSKgz8DDgROAA4Qyl1QHlX5chvgBPKvQgX9AHf1VofABwOXFLBv9ce4Bit9SHAROAEpdTh5V1SXi4D3ixkB4NSyLXWn6T9OASo2Iiu1voprXVf/48vA2PKuR4ntNZvaq03lnsdDhwG/FVr/a7WuhdYDJxS5jXZorV+Efi43OvIh9b6I631q/3//heG6FRkM35tsLP/x0j/n4q9/pVSY4CTgLsL2c+gFHIApdSPlVIfAmdS2RZ5OucDT5Z7EVVMC/Bh2s+bqFDBqVaUUvsAk4A/lXkptvS7KtYC/wSe1lpX7FqB24C5QLKQnVStkCulnlFKvW7x5xQArfX3tdZjgfuB71TyWvu3+T7GI+z95Vupu7UKtYlSaijwEHB51lNvRaG1TvS7VccAhymlDirzkixRSp0M/FNrvbrQfVXtqDet9ddcbno/8ARwfRGX40i+tSqlzgVOBqbqMif2e/i9ViIdwNi0n8f0vyYUiFIqgiHi92utl5V7PW7QWncqpZ7HiENUYkD5CGC6UurrQCOwu1LqPq31WV53VLUWuRNKqf3TfjwFeKtca8mHUuoEjEer6VrrrnKvp8p5BdhfKbWvUqoe+CawvMxrqnqUUgr4FfCm1vqWcq/HCaXUKDPzSykVBY6lQq9/rfXVWusxWut9MM7V5/yIOAxSIQcW9bsD1gPHYUSFK5WfArsBT/enS/683AuyQyl1qlJqE/Bl4HGl1Ipyrymd/qDxd4AVGAG5pVrrDeVdlT1Kqd8CfwTGK6U2KaUuKPeabDgCOBs4pv8cXdtvRVYiewHP91/7r2D4yH2n9VULUqIvCIJQ5QxWi1wQBKFmECEXBEGockTIBUEQqhwRckEQhCpHhFwQBKHKESEXBEGockTIBUEQqpz/H5m54zKhK2n6AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(data[:,0], data[:,1]) # Print the original data in blue\n",
    "\n",
    "# Apply PCA. In theory, the Eigenvector matrix must be the \n",
    "# inverse of the original rotationMatrix. \n",
    "pca = PCA(n_components=2)  # Instantiate a PCA. Choose to get 2 output variables\n",
    "\n",
    "# Create the transformation model for this data. Internally it gets the rotation \n",
    "# matrix and the explained variance\n",
    "pcaTr = pca.fit(data)\n",
    "\n",
    "# Create an array with the transformed data\n",
    "dataPCA = pcaTr.transform(data)\n",
    "\n",
    "print('Eigenvectors or principal component: First row must be in the direction of [1, n]')\n",
    "print(pcaTr.components_)\n",
    "\n",
    "print()\n",
    "print('Eigenvalues or explained variance')\n",
    "print(pcaTr.explained_variance_)\n",
    "\n",
    "# Print the rotated data\n",
    "plt.scatter(dataPCA[:,0], dataPCA[:,1])\n",
    "\n",
    "# Plot the first component axe. Use the explained variance to scale the vector\n",
    "plt.plot([0, rotationMatrix[0][0] * std1 * 3], [0, rotationMatrix[0][1] * std1 * 3], 'k-', color='red')\n",
    "# Plot the second component axe. Use the explained variance to scale the vector\n",
    "plt.plot([0, rotationMatrix[1][0] * std2 * 3], [0, rotationMatrix[1][1] * std2 * 3], 'k-', color='green')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The explanation of this chart is as follows:\n",
    "* The rotation matrix used to create our correlated variables took the original uncorrelated variables `x` and `y` and transformed them into the blue points.\n",
    "* The PCA transformation finds out the rotation matrix used to create our correlated variables (blue points). Using the PCA model to transform our data, puts back the variables as our original uncorrelated variables.\n",
    "* The explained Variance of the PCA is \n",
    "\n",
    "$$[1.0094, 0.1125] $$\n",
    "\n",
    "which is approximately\n",
    "\n",
    "$$[1, 0.333 * 0.333] = [std1^2, std2^2],$$\n",
    "\n",
    "the parameters of our original random variables x and y\n",
    "\n",
    "You can use the previous code to try with other standard deviations and correlations and convince yourself of this fact.   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## PCA as a strategy for dimensionality reduction\n",
    "\n",
    "The principal components contained in the rotation matrix, are decreasingly sorted depending on its explained Variance. It usually means that the first components retain most of the power of the data to explain the patterns that **generalize** the data. Nevertheless, for some applications, we are interested in the patterns that explain much less Variance, for example, in novelty detection. \n",
    "\n",
    "In the next figure, we can see the original data and its corresponding projection using dimenson axes as principal components. In other words, data comprised of a single variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAtHklEQVR4nO2dfZAc5Z3fv7+ZbaFZgTVS2AQYQJK5K1SWhRBsbHzKVQXsA84gkIVt7NhX5dxVqFRdKoiidBFnYkkuciilOnOknKqEil1O6ogtXuS1iOwS2FKVzxCwV+wKkJFsg3kbQ9BZWrC1A5rdffLHbM/29PTTr09PT/d8P1XYbM9Md88M8+1ff5/fiyilQAghJL+Usj4BQgghyaCQE0JIzqGQE0JIzqGQE0JIzqGQE0JIzhnK4qDnnnuuWrlyZRaHJoSQ3HL48OF/VEqNuLdnIuQrV67E+Ph4FocmhJDcIiKvem2ntUIIITmHQk4IITmHQk4IITmHQk4IITmHQk4IITknk6wVQtJg/8v7cf+z9+Ot02/hvCXn4fYrbscNH7wh69NKlbGJOnYfOI7fTDVwQbWCrdddik3ra1mfFukxFHJSCPa/vB87ntqB92bfAwC8efpN7HhqBwAUVszHJuq4a+/zaDRnAQD1qQbu2vs8AFDMBwxaK6QQ3P/s/W0Rt3lv9j3c/+z9GZ1R+uw+cLwt4jaN5ix2Hzie0RklY2yijg27DmLVtv3YsOsgxibqWZ9SbmBETgrBm6ffjLS9CPxmqhFpez8zCHcXadpgjMhJIRBIpO1F4IJqJdL2fqaXdxdZRP72hao+1YDCwoXK1LEp5KQQKHhPutJtLwJbr7sUFavcsa1ilbH1ukszOqP49OruIm1B1ZH2hYpCTkhO2bS+hns3r0WtWoEAqFUruHfz2lxaEb26u8hqXSHtCxU9ckJyzKb1tVwKt5ut113a4ZED6dxdZLWucEG1grrHMUxdqBiRk0KwdNHSSNtJf9Gru4us1hXStsEYkZNCcNdH78LdP7kbM2qmvW1IhnDXR+/K8KxIFHpxd9GryN+N/b7SylqhkJNCYBf9DFplJ4lG2oIadOy0jkMhJ4T0JWnlXRdlXcEJhZwUgkEs0S8yg1AgZBIudpJCMIgl+kWmaO0H0oYROSkEb51+K9J2koy0uy72skCoCN0jKeSkEFSGKpiemfbcXkSyFKBe2B5p510DxbJvaK2QQuAl4n7bs8JEn4+sysxtemF79KL9QJHsG0bkhPQIUxGgToB27Dvak0jSpO2hu7PoRZpgkbpHUsgJ6RF+EWAUgdIJzVSjibGJeupibsr2CLqwpZ0m2Av7plfQWiGkR5iKAP2Ephe2gCnbI2tro0jdIynkhPQIU30+/ISmF7aAuy9KtWJhsVXCHXsmI/n+WVsbReoeSWuFFAKBePYe76fBEqb6fGxaX8POx47i1HSz67Fe2QK27ZHE9+8Ha6MoVZ6MyEkhWFxeHGl7FpiMALdvXNMXtkASe6RI1kbWMCInhcBd1Rm0PStMRYBZNn9yksQeCfseilK0kyYUclIIPrDoA3jnzDue24tKP9gCSe2RoPdQpKKdNElsrYjIRSJySER+LiJHReR2EydGSBRENMOXNduJGZLaI0EFUllntuQFExH5DIA7lVLPisg5AA6LyBNKqZ8b2DchoZh6fyrSdmKGJBZPmGg768yWvJBYyJVSbwJ4c/7ffyciLwKoAaCQk55RkhLm1JzndpIucS2eMAVS/ZDZkgeM/lcuIisBrAfwjMdjt4nIuIiMnzhxwuRhCfEUcb/tJHvCRNvMbAmHMSEXkbMBPApgi1LqXffjSqkHlFKjSqnRkZERU4clBABw/pLzI20n2ROmQKpIRTtpYiRrRUQstET8QaXUXhP7JCQKt19xe8eEIKCVQ377FVx771fCFkj1Q3ZOv5NYyKWVFvANAC8qpb6W/JQIic4NH7wBE29P4OFfPIw5NYeSlHDzH9zcszFvzHWOTr/kwhcBExH5BgB/BuB5EZmc3/bXSqnvG9g3IaHY//J+7P3l3rYnPqfmsPeXe7H+n65PXcyLlOvc6wsSo20ziFLd/SnSZnR0VI2Pj/f8uKS4/PF3/tgz1bB6VhX/8Ll/SPXYG3Yd9MysqFUreHLbNake2yTuCxLQsjqy9KR5p9OJiBxWSo26t7OykxSCLPPIi5LrbKpfuh9RhLlIdzppQyEnJCF5znV2Cqvu3tzUBSmqMPfiwlIUWC1BSEL6Idc5zizQsYk6tj5ypD37U8cF1YqRWaNRy+2LcqfTCxiRE5KQrLMv4loQOx87iuas/xpZxSrj6tUjRiyOqMKc5zudXkMhJ8QAWWZfxLUgvAZT2AjQviCZsjiiCrOpQRyDAIWckHnymiGRhgXx610LKZt37Jk0sv+owpz1nU6eoJATgnxnSMS1IKoVC1ON7qi8WrGM7N9NkDDrLqT9/vn3AxRyUgjKKGMWs57bw9DrDAmT0X9cC2LHTWuw9eEjaM4t+ORWSbDjpjVG9u+FTpjzfCHtByjkpBB4ibjfdje9zJAwLVpxLYiwr+uFxcFUw2RQyAlBbzMk0hCtuBaEX4TsFu40q1SZapgMCjkhSDdDwi2KXhcMIF3R6veKSqYaJoMFQYQgvb7XtijaRTf1qQZ0U0TTEi2vc7hr7/Paop4s5mT2Q1FVnmFETsg8aWRIeImiQitP21mKk6ZoRbVysrA5mGqYDAo5ISmiEz+FVtQfR7SiZrxEFebhRWWcPtO9SFwdtjyebQ6mGsaHQk5Iiui837gtbuP411H857GJuqeIA0DaHa/zWpDVD9AjJyRFTHu/cfzrKOfgt593PIqHTBHVxyedMCInJEWier9BUWkc/zrKOfjtJ80MEuaRJ4NCTkjKuL1fuyWsW1S9bJOtDx/BzseOYmq6iQuqFVSHLc9mV0EiG9Z/9kuPtCN498Xm6tUjOHTsRCJLhHnkyaCQE+IgbZ/Wz+P2ikqbc6ot3PWpBqySwCpLR/tZk/nu02dmPB8rl1pJk3ePPY8Hn36tnXFTn2rg759+rf28uDnnzCNPBoWckHl6UQjjZyGEiT6bc6ojD33JojKscgl37JnE7gPHI1947AuXnd+uW8+cnVPtO4OgNc84lghb1iaDQk7IPHF92ihRvJ+F4GdrOHEKaSvDZOHCc8eeSWzZM4laiLsJ94UrSKD9+pe7iWqJMI88GRRyQuaJ49NGjeL9LASvqDQqTssj6G7C68JlijiWCPPI48P0Q1JolELoGZM68fETJV0Uv/Oxo57P90oFtEqC6TMzuGPPJBZbJVQrFgStvuBWWVfQH0xQWmKUqLlilbv6lPs9N21LxMQM0SJBISeFJ2xOcpycb50Ynppueh7P3dOlWrEAaT1fzb/u/Zk5fOGqi7HkrKHAmZpBmEgntPvO7LhpTdfnIwA2XLLceI8aP5hz3g2tFTIQhPG64/i0fr72jn1HtRNv7H1u2HWwa0pPoznbkRmSBD+xDrJyrLJg96fXdb1/r/fkXCew7wLSEnPmnHdDIScDQxgrIapPu/W6S7FFM9NyqtFsi7TOs/brxZKUoLsJ54XL62LUnFVd4uj1+fS67S1zzruhkJOBwVROsjtLpWKV0GjOBb6u0ZzFnQ8daf+9+8BxI4Jts2zYwvCiochTgjatr2HVtv2e5xJGHHsdITPnvBsKORkITBbNuKNPqyywStIx+1LHrFLYsmcS5ZJgVvN8XT63X563TdyUPZ04lkSwatv+WGX9aUXIzDnvhoudpPDEXYDzyozwrL6cVTh78VB7wS8MOhGvVSv4wlUXe+4nTJ73HXsmcffY8yHPYgGvhV6gdeHRLSjan4/uvNKKkNMaApJnGJET4oHO99UtDJ6abuLUdBNlEczG7PcqQLu1rbPsPQoKwINPv4bRFcsjCdum9TWMv3rSd5HVaZe4Px83aUfIzDnvhBE5KTxx0tN27Dvq6fsGEVfEgZYI25F/TRPNliU45ldYaEcbJd/60LETgVG/bZf4FRMxQu49FHIyEESZOTk2Ue9KCewVdsfDlf/EW8iv+uAyTwvEzW+mGpHzrcN42rZdoku5tO8qKOK9hUJOck3bpw0RCNfnxS0Ik0OGBcAXNZ63juacwlMvn/R87JXfNnDv5rWBkbkCsGXPZKQhFEGetm2XjE3Uez5AmvhjRMhF5Jsi8raIvGBif4SEwRlxhiXMYqDJbIsLqhXcs2ktvnDVxZFep7sw2e/1bz+7LlRk7oXu/ekWPAGgYpWw2Gp1WbzzoSParJpBzhzJElOLnd8C8HUA/8vQ/ggJJE7TJ+dioL0P94AEU7ndzgj2wWfiLV56cdfe53Hv5rW4d/PadjFPlEVWXdRs2yE79h31qDada+fK646j5vfB2Zu9x4iQK6V+LCIrTeyLkLDEjZwVWmL1/sxcR1ZK3EwRHfduXgugJbxR10CHrRIUxPNC1WjOYse+o1hyVqv4x25Zq6swdRKm2nP3geOx1ghq1UrPqzxJi56lH4rIbQBuA4CLL452m0mIF2H7d3uR9mJmrVrBpvU1bNh1MFar2L/ZfBkARC7/96Ms0vbIx189qR3PFucCaV8g2AclG3q22KmUekApNaqUGh0ZGenVYUmB8fN0s8QZ9Sbx26MsuvpdLEpo+de2JWLffTizWbY+fKS9EBx1wdKZbsg+KNnArBWSW5wVfjpE0JVhUbHKWDYcrrd2VErSEtUteyZx+c7HUY15nDv2TMa+23Azh+Cq0Oacwo59rR7qUS6QtWqlI91wqaZnuXs7+4mbhUJOcs2m9TU8ue0a+GXj3Xfr5V3l3Ns3dvfW9sI92MHuq6LDWXk/1WjinQjj0ZyYWHAti0S6YNlWjVcJvA53pK37Hpzb2U/cPEY8chH5NoB/CeBcEXkDwHal1DdM7JsQL9yZEThf/1y/cu6gBcKPrFyGp14+2V6stEqCzVdeGLpfeHBPxPSYVSrSnE0n7s9sw66DoToOTmmO59xOH908prJWPm9iP4TocAr30oqF02dm2tNz6lMNnH2ePhrUMf6qd9GNkydf6nzOdHMO//sZM0Mf0sYuGgqbluiO3t2fuVWWjolFXhkwYVrM0kc3D60V0ve4b8WnGs3EI9DGJup4MGa6YYhutR0sSjB3MwmzSoUWcass2L5xTftvr88cqiX2fh0Hw4zLizMblfjD7oek7zE57d2OMk0tJIYh6UUnLoJWP3EvMa9WrHYeulfRjme73jmF4UVDmPjKtdpjhhmXx37i5qGQk74n0S23AlZu2w+gFRmf6bGoVitWZg24FLxtFask2HHTGl8/WnehMzEuL85sVOIPhZz0PUkKf5wy1msRB4Ab151vpGK0YpWx2CrFXrzsIMDpsZtieX1apuwP9hM3Cz1y0vf0a+FPEMuGLRw6diLxfqKmTAZhD1XWoZsl6tcUi3nh2UIhJz0j7o/dzmuuaopNTFCrVuCTHh6ZilXG9o1rEmdiuPt7nzVk5ifrd166xxRahUru725soo6tjxzprBR95AjFvIdQyElPCCoCCRL5TetrWHJWuk5g1GwUHc6MjqRWhF0RaX9+pvx2v/Pye8zru9v52NGuBd3mrMLOx44aOVcSDD1y0hP8ikDcsyLtiO6v9z6H6fnWqWkvGprMYnn73da+xibqmD4zk2hf777X1A59ToLdXtdrwdErq8SNs4BH59sb8fNJKCjkpCfobtfrUw3PKsnmrOqI8rLK/IhDc65VMeouoInDnAK2PnwETVO3C1go/AlqN2uLvO7ILODpHyjkJDWcEZ8un1lEPw0n75jKH48r4lap1THMXY25feOawDJ5Z1ZJUHm+7m4pzTUN0gk9chKLIE/bvQCmnSoTQ6Oc1YUbLlkearL8oFEWwa0fuQhLFi3EasuGrVjtZoOqNXfctKarkZidq056A4WcRCZM9zqvBTBTKNXqaDh9ZgZPvnQydBn6IDGrFB49XO+IlN9rLrTwilIm79UN0Vmev2l9Dbs/s67j8d2fWcc88R4iKoMfwejoqBofH+/5cYkZdLfadm9qYKGaMi66ghSbilXusAbOXr3Ns2mWUsDvj+1KdC55RDfD0y7Nr081uj7jilX27J9C+gcROayUGnVvp0dOImO6e92yYQvDizr7fgD+LWbjZHBYpZa4GVw37EvcFzknzhFxCgsXzBrL5HMNrRUSmTC35VEWuj50/jld2zatrxmf4tOcK76I2z643zAIJ7aIO4uObFitmR8o5CQyQYtfUX/wT7100tNv375xDcomyy0HgHcbrbz1KG0NvO6kOMUnX1DISeTIy2/xS1eB6Bddu4NkOw0OAOaKHkLHwE+gZ5Vq54O7vyPdd+B1h+WXnkj6D3rkA44tvH6FIV54da8bm6jjzoeOeC6yDS8awvCiodAVlPZ5UMY7KYvg3s1rO6pe3diC67ZL3N81oO8Dzik++YIR+YCji7zufOhIJG/UFgldKuBvphqet/s646QsYrQkvQhUrDL+9rOttL7GjP80UC/BDUojdMIpPvmCEfmAo4uwbEEOG6EH9QK5oFrxHChw9eoRPHq43hUlUsS9s3nszzAoa1gnuGH7gHOKT76gkA84YYY2hJlwHnTLPX1mBmMTdU8hGV2xHDsfO+pospSdoVKbv7iYGAaRlFPTTQwvGsJ9t17e/szGJurYsc+/q6CgdQHesOtg6JRCrwZa925eyyk+OYEFQQOOl2/qhQD49a4btI/rioTcVCtW15ixsOfgh6mCoA2XLMdTL53sK2/ezvVeNmzhnekm/E0V79f65YnrvHMWB/UfuoIgeuQDjts31fUtCfJGw6a7TTWaXWlsplu0JuHJPhNxYOH+5FSAiHt9c87WwLr0wZ2PHWWGSs6htUI67I4omQ3ufQDBrU+Blkhs2TOJLXsmsWzYYt9qQwRdgLwssrGJuvbzZ4ZKfqCQkw6STDh3vjZsmiFF3By6/ipO3OLsF3UzQyU/UMhJF3EnnCfxuoOaZBF/rLLg1n9+UVcGkBu3OPtF3cxQyQ/0yIkxknjd9oIc8ccqCypW58922bCF3Z9eh3s2re3os+L2zL0sMl3UXa1YXOjMEYzICQDv9DOnbx7GakniqXKajJ7S/BSlMDaXe70j6HvT5YtzKES+YPrhgGPnJbt7o9jpZwC0dkm1YkEEmJpu4oJqBaffn9HO1rTzs/f89PWu0WUlAcqlZPMt896P3K/FrlWSVAc1hL1Qk+xhP3LShZ+n7Uw/8+ttbWMvbpYEHWLkzEcem6hj/3NvdixwDlslvDczl9o0oTxg53gDwB17JrsHUc+pwIKsJMRdEyH9A4V8gAnytONYJXOq5dnaUbod3ekuGrrGT4PCsmGrPVUJaAm5F0wFJH5QyAeYIHEQV3QdluFFQ5j4yrUd2/qp6KefmHKlX+paJjAVkPhhxCMXkesB3A+gDOB/KKV8TclYHvlzDwE/+irwzhtAZVlrW+MUsPRC4ONfAS77bPR9zLwPNE873kgJUHPA0ou69/ncQ8AP/gPQONn6u7IcWPMp4JePA++8DkgZULML+whDZTlw3lrglZ+0XttjVPt/UsD2q53zxFLkslUXQWeSP/fr19M9eBKkM7tE+52IvlOkcv6LtJ/eH0gZuPJLwI1fa/3d/g06fjM6zl0NnH574TcHtH4zf/qfW/9u/5ajaEDHOSTQkjg4jxvzWDqPPLGQi0gZwC8A/AmANwD8DMDnlVI/170mspA/9xDw2L8HmpoI0qoAG/+L/4cStA+/fT73EPC9vwRmz4Q/Z9JT1q7UC/nzr/SxkA8Ko38BXHxVtN+gjlIZQAmYc9zNhNEAwIyWxMHruDGOlWavlY8A+JVS6mWl1BkA3wFws4H9LvCjr/p/+c1G6zlJ9uG3zx99lSJOSBIOfyv6b1DH3GyniAPhNAAwoyVx8DquwWOZEPIaAGfI88b8tg5E5DYRGReR8RMnTkQ7wjtvJH9OmH3oXhPntYSQBdRs+r8jEzoR9jlR0e3T0LF6VtmplHpAKTWqlBodGRmJ9uKlFyZ/Tph96F4T57WEkAWknP7vyIROhH1OVHT7NHQsE0JeB3CR4+8L57eZ4+NfaflJOqxK6zlJ9uG3z49/BSgvCv9aQkgnV34p+m9QR6kMlFyVwGE0ADCjJXHwOq7BY5lIP/wZgD8UkVVoCfjnAPwrA/tdwF4MSLLS7LWPsFkr9v/nMGvFLwsCfo8VqD6nJ8XLBrJPAl8b8jwWdhLtsVQyXdxZK0C2WSsmtCQO7uMaPpap9MNPAvg7tNIPv6mU+k9+z2eJfm/RlWDreo+fNVTSltr3K/1Qov+KY4LSqm37tVr66103+JbFh5225GbZsIXtG9d4Vofax77v1ss5DSjHpFqir5T6PoDvm9jXoGO674Xf/nS9x3XVhSQ8foU97guoe8B13CrOU9Ot6UtLK5bnhVg3AJu9VfIPKzv7iKAfuIn9bX34CHY+drSjhN5ZIg4ED4YomPNihGXDnZ6tV1dBqySYPjODLR4XSuf0njADsQHvQRKN5iwWWyVUrLJ2yhN7qxQP9iPvI7zK2JPMTvTaX3NO4dR0Ewrzwv7IEVy+83Gs2rYfG3YdxNhE3Xf+ZlmEIu7BDZed3/G3exZqtWIB4j8RyY7Ew8w/FUA7DWhqutlx7Fq1Quuk4DAi7yN0t9Rxb7XDvK45q9q34fYdwL2bWwMK7MjcGYEHjRIbVB49XMfoiuUdYmlHvmMTddz50JHAz64kglXb9uOCagW3XFnDoWMnuj5/oCXiX7jq4vbjbmwLhcI9ODAi7yN0jZHiNEwam6h7Lv4F4XUHQOkOxh4obd/V2Nj2VpgL4KxS7TulRw+37oxe2XUD7rv18o7o+r5bL8c9m9Z6Ru5hBmWT4sGIvI/QTWuJ+sO0xSNO50JgITJnt8LouNc1wnR91Hndtmeui665cElsKOQ9JCgjxdQPM2nL2LIIRdyHoGn1ThEOsreqmgwTIJw1RguFABTynhE2I8XEDzPJEAJ3tgPp5pzFQzh9ZsZ3qpH9HQRloEw1mtosIPYgJ2GhR94jTGek+BFXAOzsBncqnfs5i8p90+06E6YaTUB1pxw6sb+DrdddCivg83K0EW9Dr5tEgULeI3RRcn2qgQ27Dnak/yUlTPqam1q10s4n//17M57PsUqCq1eP4MwAz9e0ac4pvNvw/pwA4OrVjsZwIT4uBTBdkMSG1kqP0N1iCxYGFyctALKxXxsm5Q3ojP52HzjeNeXepjmn8O1nOKTBxu+z/fYzr2N0xXLfz9OJ80JKSFQYkfcIryjZyxuNY7eMTdS7ovpN62v428+u0zZCKot4Rn9B/vog5pGXY+RxziqFLXsmQ1Vo0kYhSWFE3iO8MlJ0P/Ioi5VBi6jjr57Eg0+/1nHBsJsk2edzx55J7D5wHFevHkEpICMjKGOjiMwqZXwRuCyCOaWYMkiMQCHvIe6MFF2XuyiLlX6LqJvW13DPprXtW3xnSiMAbH34SPu2vz7VwN8//ZrvsSpWGR9YXMb/+91gjb2rzX9m9mcYdLELgt0GiWko5BliogDIbxF1/Vcfx/aNazouIHYue9Q2qQLgwmWL8cu3Twc+t0h4fR8fqAzh9+/NhPK+3dQYgZMUoJBniIkCIF3LUqDVoGnrI0fax7p77PkumyUsCsi1iIft2LjhkuV45beNrrsX5wX31HQTVllQrVh4p9EMHaGXRbigSVKBQp4xYQuAvKpCAeDd9/wHQDRnFXbsO4od+47mbliESe679fL257e0YnkW9Gy4ZDke/Dcf63rthl0Hu7tIziosOWsIk9uv9RzQ4cWgrS2Q3kEhzwGefcUfOQIohOqnMsgCDrTsDOcFc8Oug56fySu/7babxibqgYvS7jsrEe/vpcZKTZISFPIc4NlXvA+Lcvo1o8XtcYdtF2xfQHUotC4Kth3mXIcw0fyMkLAwjzwHJOmd0ks+/9GL0hngmwCR7uKqsO2CwzQfs9M9nRW57qESrNQkaUMhzwFxe6eUBBi2evcVx11I7QXOoqnpMzOwSp2XHK+IOewF1F3EZXruKiFBUMhzgFdVqFWWwC9vTgHv99CC6WcRv2vv86hPNaDQyjqZw0KjqrIIbrmye9E5ygXUFn33sbwidkJMQyHvc+zortGcbZeK16oV7P70Oiz16b5nMxt3ukRBUKrVc8ZtkczOqY7xdY8erneJra6tghe26PeyyyUhNlzs7CPct+RXrx7Bo4frbWGwFxKnz7S67k35DPIlC4RZgHVWwzpZbJXan3+1YuHGded3fCdApy1jeu4qIWFgRN4neN2SP/j0a56Lbaemm7hr7/MYXhStVS3xxym29vfhnHr//swcRlcs913INDl3lZCwMCLvE7xuyf3iSE7x6STOoGk3TrH1s0ie3HaNdvHS1NxVQqLAiLxP6PWtd8UqY8Mly3t6zKyxW/cuG7YCs1biWiRMPSRZwIi8T/AbPJHGcuUVFy/1rGTMI7VqBe8GPMfdcTAoRVD3fYSxSDgQmfQaCnmP0QnI1atHPPuGX3HxUjz50snQ+w9bXRlln/2MoGVn/Mfn9M/x6jgYJLa0SEieoLXSQ3Q5xnePPY9HD9c7RFwA3HJlLVLUXLHKmVVXCoBKD4uPbBT8x+KJwNfT1kGLhOQJRuQ9RLeA9u1nXu+KohWAQ8dOhPbOyyJtoQkaEJEGdnfBqH3O3dizK8cm6tj6yJHAnjJpNqKiRULyAoW8h+hEWWeFBI2Es7HKgiWLhtoj25YsKuP0md5mtdyxZzKxl++0LtwdBavDVtcwB1odhLSgkPcQnSjrfO2SCOpTDc8FzyWLypg+M9sWOLsta9KIOC5xRbw2/5mURToqIO1o2BkRs4cJId5QyBMQVVh0C2i3XFnrqhYEFiJ1L5GcUwt2xqmEFZ61EFG/TdVnIlFUvnjVxRhdsdx3eLQTWh2EeJNodUpEPiMiR0VkTkRGTZ1UHojTHEm3gHbPprUd20shVivt6DVp/rmd0RGWG9edn9iXFmmJ+D2b1rI3CSEGSBqRvwBgM4D/buBcckXQ9HoduqjS3j42UceWPZOhziGsh66jYpVx9eoR3+EJbvY/9ybeT1hVuniojNEVrWIk9iYhJDmJInKl1ItKqYEMndISoCiR6GKr1G6gFZUli1oFMoeOnYhU7n9quonp5lysY9o4I272JiEkOT1L/BWR20RkXETGT5w40avDpkZaAhTlQtBozsX2x+3kj6wWR+336dUqltkohEQjUMhF5Ici8oLHPzdHOZBS6gGl1KhSanRkZCT+GfcJaQlQGpFo2aOjVKM5i52PHTV+rLCURLBq237sPnAct1xZQ7Wy0Ft9cQaFRYTkmUCPXCn1iV6cSN5w5znrslbiZLaE9cjDULHKWuskKJpfNmwlzojRYWfk1Kca2PPT1zsmNthtegH/qk1CSAumH4ZEJ8h+QuOepu6XWpcWt1xZw6FjJyJbKIJgoTdF02OKUZiFY0JIi6Tph58SkTcAfAzAfhE5YOa0+ou4cxjjpNb5PRanh8qhYye0NpDTznATpcDHy7pxY6dbRoGZK4SEI2nWyneVUhcqpc5SSv0zpdR1pk6sn4ib6xwns8XvsaU+wuu3P13++o6b1nQJfFQEwEv3fhKv7LrB93kXVCv4zXwFZ1iYuUJIOGithCBuqmGcntZ+fcnjVFTax9LZQOOvnvRs2hV1/4C+QlSwkB0T9jjMXCEkPEwPmGdsoo4Nuw5i1bb92LDrYIdtEifVcGyi7pnjHSRQW6+7FFa5O2qNI7NBxxqbqOPRw/XYIu7ev27qfNS9s2UsIdFgRI7gRcmoQwbc+7OpVizsuGlNYOXnjn1HE/czGbZK+JsAMfSyjIKozVskFzhK+zfsOtjuUOiU7WGrFLl4yG5jSwgJD4UcweX2YVMN/fYHAEvOGgoVZb5joClVwyWgXlk3UTNZ3CLrvmC5s1z8RHzZsIX3mnOcwEOIASjkCOeBR+m8l7R8P0n/FBsFYMt8f/KrV490dFe07zhKslDhGYRVki6RjRPRAy3B3r5xTXsfcdrSui9MOD/yaRBSGCjkSDZoN439Xb16xNiUn/pUo2sWKIBAAXYWA+ksobjpgU7/O44P7mWFnX1eq6siIYMIhRzhB+2GrdL02p9VFpx+fwartu0PjD4PHfPuRRNn4RAxXhPWp45z51CrVhIvYsa9EyCkqAx81ootzo3mbDvH2StrIkpRkDtve9mwBahW+mCYgiJdpJt0lJqbasXy7Rfjl8kDeGep+GHKA2ehECGdDLSQO8UZaOU422LjjhqjFgVtWl/Dk9uuwa933YDhRUNdZeh+r02jEMbtOlSsMnbctEY7KT7MhWvT+hpuubLmOQijVq3gi1ddnMoUehYKEdLJQFkrbmvk9PszoYdDJFnAjPpaL2smCJn3XZZWLJw+M9Mxfd4eJ3fo2AlPW8hLXMMMzhibqGPPz17vWjAtzb+HtPLA43w+hBSZgRFyrwUyHV4Cm2QBM+prbQGM0gVRFPC1Wy9vR9NhujLa+d9ezwlz8dl94HjHBcNmDsCOfUdTE3KvdNB3udBJBpiBEfIoC2ReAhu1KCjpa6MWBjnF00RXxjAXH7+7EVMDmnW43+Pa/5nq4QjpawZGyMMukOkENkxRkB0J1+ebQ80q1R5ufMuVtXZPk7IIbrkyWGxPRxzjFlY8w9gmYS4+JvLdCSHJGRgh14nOsmELw4uGQhWl+EW67ijXOThh68NHAFnYNqsUHj1cx+iK5b7VoV62hQnCFkDZ56H7bPyGYCwbjt6pkRASj4ERcl2EuX2jf++TsPhZN3EGJ8RJsQsrnmE9+yCLZtP6GsZfPdlVvGSVpV25SQhJn4ER8qj9UqISR3j9XhPVtggST+cCaHXYglWSjgtM3BzvezatxeiK5al9roSQYAZGyIFo/VKiEscv9st42Xrdpbhjz2SoIqCSAEsWDeGO+d4qXt69u7mVVRZUKxbeaTQTi2+anyshJJiBEvI08ctttkoCCLpyu4OyVsZfPenZJ8XNnFpY6PTKQPGyfZqzCkvOGsLk9mvDvD1CSB+Tm8rOoHLxrHGW5QPoKPff/Zl12P3pdZGrHO/ZtBb33Xp5x+vsaklAP8PTXTWatBsjIaS/yUVE3g/T6MMQZnHQTVDxjm6fG3YdDF3UZLq7IyGkv8hFRB53+HG/E6URl5sgP94p0l7NrTjEgZDikAshL6o1oLtA3fnQkUALyW8avVuk3d0YOROTkGKRC2slz9aAn3WiuxA5i4m8LKSxCf+ByV4izcwSQopLLiLyvFoDQdZJmAuR20Iam6i3KkU1eA1uMLlQ3O+LzoQMIrkQ8rxaA0HeftjBDM7Ifce+o56VooB+qlFcH96NyX0RQsyRC2sFyI814LRSdOaHLczuatPSfKMtN87I3a8xltfFLUyDrLCY3BchxBy5EfI84E6T1OEUZucFyuv1USwkLzE1uVBc1EVnQvJOLqyVvBCm57mfMIexkHSNsXTbdT58nIVik/sihJiDQm4Qv8g0rLfvnPX55LZrup67feMaWOXO1EO/hlkmF4rzuuhMSNGhtWIQXZpkrVrBk9uuMXKMqF0cTXZ9TLuDZBLOX3I+3jz9pud2QoqOKJ985LQYHR1V4+PjPT9u2ug87jxk2DgJM/Oz39j/8n7c/ZO7MaMWpioNyRDu+Rf34IYP3pDhmRFiDhE5rJQadW9nRG6Qfo5Yw5KXvjZeiAicqULiU/1KSJFIFJGLyG4AGwGcAfASgH+tlJoKel0/RORpR51ZRLUmjqlrxmXSHkqDax+5VmutPP7pxzM4I0LMo4vIky52PgHgw0qpywD8AsBdCffXE9IubMmicMbUMfOaYugl4n7bCSkSiYRcKfW4Um1T8mkAFyY/pfRJu5tiFt0aTR1zaSVaGiMhJHtMph/+OYAfGNxfaqQddWYR1Zo45thEHafPzHRtt0rCFENC+phAIReRH4rICx7/3Ox4zpcBzAB40Gc/t4nIuIiMnzhxwszZxyTtwpYsCmdMHHP3geMd4+hszl481PcLnYQMMoFCrpT6hFLqwx7/fA8ARORLAG4E8AXls3KqlHpAKTWqlBodGRkx9gbikHZhSxaFMyaOqYvep6b1/V0IIdmTyFoRkesB/BWAm5RS02ZOKX3sUnhnWftZQ/E/CndrVwA979ZookNknkvwdYU/LAgig0DSPPKvAzgLwBPzObtPK6X+beKzCsBUat97zbn2v081mrHypXV51/duXtvzdL2kHSK3XndpoqZdWbLinBWeGSorzlmRwdkQ0lsSCblS6g9MnUhYTBWsmGrJWqTWrnkuaHrmrWcibSekSOSustOUcJrKLMlr3rWOvPR9d6M03d912wkpErnrfmhKOE35wXn2lQkhxSB3Qm5KOE1llrC1a39QKXt//7rthBSJ3Am5KeE0NQc0r/NEi8b2P9qOkus/5xJK2P5H2zM6I0J6Ry7b2OaxzSpJn/0v78f9z96Pt06/hfOWnIfbr7idLWxJodA1zcqlkBNCyCCSVvdDQgghGUMhJ4SQnEMhJ4SQnEMhJ4SQnEMhJ4SQnJNJ1oqInADwas8PHI5zAfxj1ieRMkV/j0V/fwDfY1GI+h5XKKW6+oBnIuT9jIiMe6X3FImiv8eivz+A77EomHqPtFYIISTnUMgJISTnUMi7eSDrE+gBRX+PRX9/AN9jUTDyHumRE0JIzmFETgghOYdCTgghOYdC7oGI7BaRYyLynIh8V0SqWZ+TSUTkMyJyVETmRKRQ6V0icr2IHBeRX4nItqzPxzQi8k0ReVtEXsj6XNJCRC4SkUMi8vP5/05vz/qcTCIii0XkpyJyZP797Uy6Twq5N08A+LBS6jIAvwBwV8bnY5oXAGwG8OOsT8QkIlIG8F8B/CmADwH4vIh8KNuzMs63AFyf9UmkzAyAO5VSHwJwFYC/LNj3+D6Aa5RS6wBcDuB6EbkqyQ4p5B4opR5XSs3M//k0gAuzPB/TKKVeVEodz/o8UuAjAH6llHpZKXUGwHcA3JzxORlFKfVjACezPo80UUq9qZR6dv7ffwfgRQCFmRyjWvx+/k9r/p9EWScU8mD+HMAPsj4JEooagNcdf7+BAgnAICIiKwGsB/BMxqdiFBEpi8gkgLcBPKGUSvT+hoycVQ4RkR8COM/joS8rpb43/5wvo3Wb92Avz80EYd4fIf2MiJwN4FEAW5RS72Z9PiZRSs0CuHx+/e27IvJhpVTsdY+BFXKl1Cf8HheRLwG4EcDHVQ6T7YPeX0GpA7jI8feF89tIzhARCy0Rf1AptTfr80kLpdSUiBxCa90jtpDTWvFARK4H8FcAblJKTWd9PiQ0PwPwhyKySkQWAfgcgH0ZnxOJiIgIgG8AeFEp9bWsz8c0IjJiZ8KJSAXAnwA4lmSfFHJvvg7gHABPiMikiPy3rE/IJCLyKRF5A8DHAOwXkQNZn5MJ5heo/x2AA2gtkD2klDqa7VmZRUS+DeD/ArhURN4Qkb/I+pxSYAOAPwNwzfzvb1JEPpn1SRnkfACHROQ5tIKPJ5RS/yfJDlmiTwghOYcROSGE5BwKOSGE5BwKOSGE5BwKOSGE5BwKOSGE5BwKOSGE5BwKOSGE5Jz/D7yOktHi1OKtAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "nPoints = len(data)\n",
    "\n",
    "# Plot the original data in blue\n",
    "plt.scatter(data[:,0], data[:,1])\n",
    "\n",
    "#Plot the projection along the first component in orange\n",
    "plt.scatter(data[:,0], np.zeros(nPoints))\n",
    "\n",
    "#Plot the projection along the second component in green\n",
    "plt.scatter(np.zeros(nPoints), data[:,1])\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## PCA as a strategy to plot complex data\n",
    "\n",
    "The next chart shows a sample diagram displaying a dataset of pictures of cats and dogs. Raw pictures are composed of hundreds or even thousands of features. However, PCA allows us to reduce that many features to only two. In that reduced space of uncorrelated variables, we can easily separate cats and dogs. \n",
    "\n",
    "<img src = './images/catdog.png'>\n",
    "\n",
    "You will learn how to generate a chart like this with word vectors in this week's programming assignment.\n",
    "\n",
    "### Optional Reading:\n",
    "\n",
    "To get more clarity on the Rotation Matrix, you can read the <i>optional reading item</i> which follows this notebook, <a href=\"https://www.coursera.org/learn/classification-vector-spaces-in-nlp/supplement/fwEUM/the-rotation-matrix-optional-reading\">The Rotation Matrix</a>."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}