{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "import matplotlib as mpl\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import numpy.linalg as la" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Dimension Reduction" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "np.random.seed(123)\n", "np.set_printoptions(3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### PCA from scratch\n", "\n", "Principal Components Analysis (PCA) basically means to find and rank all the eigenvalues and eigenvectors of a covariance matrix. This is useful because high-dimensional data (with $p$ features) may have nearly all their variation in a small number of dimensions $k$, i.e. in the subspace spanned by the eigenvectors of the covariance matrix that have the $k$ largest eigenvalues. If we project the original data into this subspace, we can have a dimension reduction (from $p$ to $k$) with hopefully little loss of information.\n", "\n", "For zero-centered vectors,\n", "\n", "\\begin{align}\n", "\\text{Cov}(X, Y) &= \\frac{\\sum_{i=1}^n(X_i - \\bar{X})(Y_i - \\bar{Y})}{n-1} \\\\\n", " &= \\frac{\\sum_{i=1}^nX_iY_i}{n-1} \\\\\n", " &= \\frac{XY^T}{n-1}\n", "\\end{align}\n", "\n", "and so the covariance matrix for a data set X that has zero mean in each feature vector is just $XX^T/(n-1)$. \n", "\n", "In other words, we can also get the eigendecomposition of the covariance matrix from the positive semi-definite matrix $XX^T$.\n", "\n", "We will take advantage of this when we cover the SVD later in the course.\n", "\n", "Note: Here we are using a matrix of **row** vectors" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "n = 100\n", "x1, x2, x3 = np.random.normal(0, 10, (3, n))\n", "x4 = x1 + np.random.normal(size=x1.shape)\n", "x5 = (x1 + x2)/2 + np.random.normal(size=x1.shape)\n", "x6 = (x1 + x2 + x3)/3 + np.random.normal(size=x1.shape)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### For PCA calculations, each column is an observation" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "xs = np.c_[x1, x2, x3, x4, x5, x6].T" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[-10.856, 9.973, 2.83 , -15.063, -5.786, 16.514, -24.267,\n", " -4.289, 12.659, -8.667],\n", " [ 6.421, -19.779, 7.123, 25.983, -0.246, 0.341, 1.795,\n", " -18.62 , 4.261, -16.054],\n", " [ 7.033, -5.981, 22.007, 6.883, -0.063, -2.067, -0.865,\n", " -9.153, -0.952, 2.787],\n", " [-10.091, 9.144, 2.171, -14.452, -5.93 , 17.831, -24.971,\n", " -3.539, 13.002, -8.794],\n", " [ -0.684, -5.433, 4.485, 4.151, -3.025, 9.405, -12.987,\n", " -12.12 , 8.496, -11.511],\n", " [ 1.618, -5.193, 10.388, 6.864, -0.771, 6.267, -8.769,\n", " -11.222, 3.621, -7.55 ]])" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "xs[:, :10]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Center each observation" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "xc = xs - np.mean(xs, 1)[:, np.newaxis]" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[-1.113e+01, 9.702e+00, 2.559e+00, -1.533e+01, -6.057e+00,\n", " 1.624e+01, -2.454e+01, -4.560e+00, 1.239e+01, -8.938e+00],\n", " [ 6.616e+00, -1.958e+01, 7.318e+00, 2.618e+01, -5.090e-02,\n", " 5.368e-01, 1.991e+00, -1.842e+01, 4.457e+00, -1.586e+01],\n", " [ 7.984e+00, -5.030e+00, 2.296e+01, 7.834e+00, 8.882e-01,\n", " -1.115e+00, 8.609e-02, -8.202e+00, -7.123e-04, 3.738e+00],\n", " [-1.023e+01, 9.007e+00, 2.033e+00, -1.459e+01, -6.068e+00,\n", " 1.769e+01, -2.511e+01, -3.676e+00, 1.286e+01, -8.931e+00],\n", " [-7.496e-01, -5.498e+00, 4.419e+00, 4.085e+00, -3.091e+00,\n", " 9.339e+00, -1.305e+01, -1.219e+01, 8.431e+00, -1.158e+01],\n", " [ 1.801e+00, -5.009e+00, 1.057e+01, 7.048e+00, -5.873e-01,\n", " 6.451e+00, -8.585e+00, -1.104e+01, 3.805e+00, -7.366e+00]])" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "xc[:, :10]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Covariance\n", "\n", "Remember the formula for covariance\n", "\n", "$$\n", "\\text{Cov}(X, Y) = \\frac{\\sum_{i=1}^n(X_i - \\bar{X})(Y_i - \\bar{Y})}{n-1}\n", "$$\n", "\n", "where $\\text{Cov}(X, X)$ is the sample variance of $X$." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "cov = (xc @ xc.T)/(n-1)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[128.578, -2.133, -6.661, 127.395, 63.27 , 42.369],\n", " [ -2.133, 95.051, 2.269, -1.349, 44.885, 31.637],\n", " [ -6.661, 2.269, 94.934, -5.048, -1.141, 30.275],\n", " [127.395, -1.349, -5.048, 126.993, 63.196, 42.723],\n", " [ 63.27 , 44.885, -1.141, 63.196, 54.408, 36.838],\n", " [ 42.369, 31.637, 30.275, 42.723, 36.838, 36.563]])" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "cov" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Check" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[128.578, -2.133, -6.661, 127.395, 63.27 , 42.369],\n", " [ -2.133, 95.051, 2.269, -1.349, 44.885, 31.637],\n", " [ -6.661, 2.269, 94.934, -5.048, -1.141, 30.275],\n", " [127.395, -1.349, -5.048, 126.993, 63.196, 42.723],\n", " [ 63.27 , 44.885, -1.141, 63.196, 54.408, 36.838],\n", " [ 42.369, 31.637, 30.275, 42.723, 36.838, 36.563]])" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.cov(xs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Eigendecomposition" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "e, v = la.eigh(cov)" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "idx = np.argsort(e)[::-1]" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "e = e[idx]\n", "v = v[:, idx]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Explain the magnitude of the eigenvalues\n", "\n", "Note that $x4, x5, x6$ are linear combinations of $x1, x2, x3$ with some added noise, and hence the last 3 eigenvalues are small. " ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAERpJREFUeJzt3X+s3XV9x/Hna6XgFX9Uxh0rt9US\nx7rgmK25QYxmcTotMrNWshlcpmwhqckw08zUUf8ZJiO41F8z2cyqEHFTkEgF4oyVIcawKHhLkfLD\nzg5Beqn0+qP+yDqF+t4f91u8uLb33B+np/dzn4/k5H7P5/v5nu/7E8LrfPs5n3O+qSokSe36tUEX\nIEnqL4Nekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1LiTBl0AwOmnn16rVq0adBmS\ntKDs2LHje1U1PF2/EyLoV61axdjY2KDLkKQFJckjvfRz6kaSGmfQS1LjDHpJapxBL0mNM+glqXEn\nxKqb2bhp5zhbtu/msQMHOXPZEJvWrWbD2pFBlyVJJ5wFGfQ37Rxn87ZdHHziEADjBw6yedsuAMNe\nkn7Fgpy62bJ991Mhf9jBJw6xZfvuAVUkSSeuBRn0jx04OKN2SVrMFmTQn7lsaEbtkrSYLcig37Ru\nNUNLlzytbWjpEjatWz2giiTpxDVt0Cd5RpK7knwjyf1J3tO1n5XkziR7knw6ycld+ynd8z3d/lXz\nXfSGtSNcddG5nLxksvyRZUNcddG5fhArSUfQyxX9z4BXVdWLgTXABUnOB/4B+GBV/RbwQ+DSrv+l\nwA+79g92/ebdhrUjrH3+Ml561mn85+WvMuQl6SimDfqa9NPu6dLuUcCrgM907dcCG7rt9d1zuv2v\nTpJ5q1iSNCM9zdEnWZLkHmA/cCvw38CBqnqy67IXOHxJPQI8CtDt/xHw60d4zY1JxpKMTUxMzG0U\nkqSj6inoq+pQVa0BVgDnAb8z1xNX1daqGq2q0eHhaX83X5I0SzNadVNVB4DbgZcBy5Ic/mbtCmC8\n2x4HVgJ0+58LfH9eqpUkzVgvq26GkyzrtoeA1wAPMhn4f9J1uwS4udu+pXtOt/9LVVXzWbQkqXe9\n/NbNcuDaJEuYfGO4oao+l+QB4Pokfw/sBK7u+l8N/GuSPcAPgIv7ULckqUfTBn1V3QusPUL7Q0zO\n1/9q+/8Cfzov1UmS5mxBfjNWktQ7g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLU\nOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z\n6CWpcdMGfZKVSW5P8kCS+5O8vWu/Isl4knu6x4VTjtmcZE+S3UnW9XMAkqRjO6mHPk8C76yqu5M8\nG9iR5NZu3wer6n1TOyc5B7gYeBFwJvAfSX67qg7NZ+GSpN5Me0VfVfuq6u5u+yfAg8DIMQ5ZD1xf\nVT+rqm8De4Dz5qNYSdLMzWiOPskqYC1wZ9f0tiT3JrkmyfO6thHg0SmH7eXYbwySpD7qOeiTPAu4\nEXhHVf0Y+AjwQmANsA94/0xOnGRjkrEkYxMTEzM5VJI0Az0FfZKlTIb8J6tqG0BVPV5Vh6rqF8BH\n+eX0zDiwcsrhK7q2p6mqrVU1WlWjw8PDcxmDJOkYell1E+Bq4MGq+sCU9uVTur0BuK/bvgW4OMkp\nSc4Czgbumr+SJUkz0cuqm5cDbwZ2Jbmna3s38KYka4ACHgbeClBV9ye5AXiAyRU7l7niRpIGZ9qg\nr6o7gBxh1+ePccyVwJVzqEuSNE/8ZqwkNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWp\ncQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn\n0EtS4wx6SWqcQS9JjZs26JOsTHJ7kgeS3J/k7V37aUluTfKt7u/zuvYk+XCSPUnuTfKSfg9CknR0\nvVzRPwm8s6rOAc4HLktyDnA5cFtVnQ3c1j0HeB1wdvfYCHxk3quWJPVs2qCvqn1VdXe3/RPgQWAE\nWA9c23W7FtjQba8HPlGTvgYsS7J83iuXJPVkRnP0SVYBa4E7gTOqal+367vAGd32CPDolMP2dm2S\npAHoOeiTPAu4EXhHVf146r6qKqBmcuIkG5OMJRmbmJiYyaGSpBnoKeiTLGUy5D9ZVdu65scPT8l0\nf/d37ePAyimHr+janqaqtlbVaFWNDg8Pz7Z+SdI0ell1E+Bq4MGq+sCUXbcAl3TblwA3T2l/S7f6\n5nzgR1OmeCRJx9lJPfR5OfBmYFeSe7q2dwPvBW5IcinwCPDGbt/ngQuBPcD/AH85rxVLkmZk2qCv\nqjuAHGX3q4/Qv4DL5liXJGme+M1YSWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMM\neklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCX\npMYZ9JLUOINekho3bdAnuSbJ/iT3TWm7Isl4knu6x4VT9m1OsifJ7iTr+lW4JKk3vVzRfxy44Ajt\nH6yqNd3j8wBJzgEuBl7UHfPPSZbMV7GSpJmbNuir6ivAD3p8vfXA9VX1s6r6NrAHOG8O9UmS5mgu\nc/RvS3JvN7XzvK5tBHh0Sp+9XZskaUBmG/QfAV4IrAH2Ae+f6Qsk2ZhkLMnYxMTELMuQJE1nVkFf\nVY9X1aGq+gXwUX45PTMOrJzSdUXXdqTX2FpVo1U1Ojw8PJsyJEk9mFXQJ1k+5ekbgMMrcm4BLk5y\nSpKzgLOBu+ZWoiRpLk6arkOS64BXAqcn2Qv8HfDKJGuAAh4G3gpQVfcnuQF4AHgSuKyqDvWndElS\nL6YN+qp60xGarz5G/yuBK+dSlCRp/vjNWElqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGjft\nOnqdOG7aOc6W7bt57MBBzlw2xKZ1q9mw1t+Mk3RsBv0CcdPOcTZv28XBJya/aDx+4CCbt+0CMOwl\nHZNTNwvElu27nwr5ww4+cYgt23cPqCJJC4VBv0A8duDgjNol6TCDfoE4c9nQjNol6TCDfoHYtG41\nQ0uffvvdoaVL2LRu9YAqkrRQ+GHsAnH4A9d3feZefn7oF4y46kZSjwz6BWTD2hGuu+s7AHz6rS8b\ncDWSFgqnbiSpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNmzbok1yTZH+S+6a0nZbk\n1iTf6v4+r2tPkg8n2ZPk3iQv6WfxkqTp9XJF/3Hggl9puxy4rarOBm7rngO8Dji7e2wEPjI/ZUqS\nZmvan0Coqq8kWfUrzeuBV3bb1wJfBv62a/9EVRXwtSTLkiyvqn3zVbAWF++qJc3dbH/r5owp4f1d\n4IxuewR4dEq/vV2bQa8Z865a0vyY84ex3dV7zfS4JBuTjCUZm5iYmGsZapB31ZLmx2yD/vEkywG6\nv/u79nFg5ZR+K7q2/6eqtlbVaFWNDg8Pz7IMtcy7aknzY7ZBfwtwSbd9CXDzlPa3dKtvzgd+5Py8\nZsu7aknzo5flldcBXwVWJ9mb5FLgvcBrknwL+MPuOcDngYeAPcBHgb/qS9VaFLyrljQ/ell186aj\n7Hr1EfoWcNlci5LAu2pJ88U7TOmE5l21pLnzJxAkqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6\nSWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJek\nxhn0ktQ4g16SGmfQS1LjTprLwUkeBn4CHAKerKrRJKcBnwZWAQ8Db6yqH86tTEnSbM3HFf0fVNWa\nqhrtnl8O3FZVZwO3dc8lSQPSj6mb9cC13fa1wIY+nEOS1KO5Bn0BX0yyI8nGru2MqtrXbX8XOGOO\n55AkzcGc5uiBV1TVeJLfAG5N8s2pO6uqktSRDuzeGDYCPP/5z59jGZKko5nTFX1VjXd/9wOfBc4D\nHk+yHKD7u/8ox26tqtGqGh0eHp5LGZKkY5h10Cc5NcmzD28DrwXuA24BLum6XQLcPNciJUmzN5ep\nmzOAzyY5/DqfqqovJPk6cEOSS4FHgDfOvUxJ0mzNOuir6iHgxUdo/z7w6rkUJUmaP34zVpIaZ9BL\nUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1\nzqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGndSvF05yAfCPwBLgY1X1\n3n6dS2rJTTvH2bJ9N48dOMiZy4bYtG41G9aODLqsvlqMYz6e+hL0SZYA/wS8BtgLfD3JLVX1QD/O\nJ7Xipp3jbN62i4NPHAJg/MBBNm/bBdBs8C3GMcPxfXPr19TNecCeqnqoqn4OXA+s79O5pGZs2b77\nqcA77OATh9iyffeAKuq/xTjmw29u4wcOUvzyze2mneN9OV+/pm5GgEenPN8LvHS+T3LBlz/Fb048\nyiN3PGe+X/qE9Rf7fgzgmBv11w99/6j7Hrn/48evkONoMY6Z7xzgiicn39weeu4I//J76596c+vH\nVf3APoxNsjHJWJKxiYmJWb3GaaeewjNPXjLPlZ3YnnnyEsfcsFNOOvI4j9begsU45p89eeiI7Y8d\nONiX8/Xrin4cWDnl+Yqu7SlVtRXYCjA6OlqzOcn6j71vtvUtWC8YdAEDsJjGvHPnOFdMma8GGFq6\nhKsuOpcXNDpfvRjH/Gfv/RLjRwj1M5cN9eV8/bqi/zpwdpKzkpwMXAzc0qdzSc3YsHaEqy46l5Fl\nQwQYWTbEVRed2/SHkotxzJvWrWZo6dP/xTK0dAmb1q3uy/lSNauL6elfOLkQ+BCTyyuvqaorj9Z3\ndHS0xsbG+lKHJJ2I5mPVTZIdVTU6bb9+Bf1MGPSSNHO9Br3fjJWkxhn0ktQ4g16SGmfQS1LjDHpJ\natwJseomyQTwyCwPPx343jyWsxA45sXBMS8OcxnzC6pqeLpOJ0TQz0WSsV6WF7XEMS8OjnlxOB5j\ndupGkhpn0EtS41oI+q2DLmAAHPPi4JgXh76PecHP0UuSjq2FK3pJ0jEs6KBPckGS3Un2JLl80PX0\nW5JrkuxPct+gazlekqxMcnuSB5Lcn+Ttg66p35I8I8ldSb7Rjfk9g67peEiyJMnOJJ8bdC3HQ5KH\nk+xKck+Svv6q44KduuluQP5fTLkBOfCmlm9AnuT3gZ8Cn6iq3x10PcdDkuXA8qq6O8mzgR3Ahsb/\nOwc4tap+mmQpcAfw9qr62oBL66skfwOMAs+pqtcPup5+S/IwMFpVff/ewEK+ol90NyCvqq8APxh0\nHcdTVe2rqru77Z8ADzJ5T+Jm1aSfdk+Xdo+FeUXWoyQrgD8CPjboWlq0kIP+SDcgbzoAFrskq4C1\nwJ2DraT/ummMe4D9wK1V1fqYPwS8C/jFoAs5jgr4YpIdSTb280QLOei1iCR5FnAj8I6q+vGg6+m3\nqjpUVWuYvN/yeUmanapL8npgf1XtGHQtx9krquolwOuAy7qp2b5YyEE/7Q3I1YZunvpG4JNVtW3Q\n9RxPVXUAuB24YNC19NHLgT/u5qyvB16V5N8GW1L/VdV493c/8Fkmp6P7YiEHvTcgXwS6DyavBh6s\nqg8Mup7jIclwkmXd9hCTCw6+Odiq+qeqNlfViqpaxeT/x1+qqj8fcFl9leTUbnEBSU4FXgv0bTXd\ngg36qnoSeBuwnckP6G6oqvsHW1V/JbkO+CqwOsneJJcOuqbj4OXAm5m8yrune1w46KL6bDlwe5J7\nmbygubWqFsWSw0XkDOCOJN8A7gL+vaq+0K+TLdjllZKk3izYK3pJUm8MeklqnEEvSY0z6CWpcQa9\nJDXOoJekxhn0ktQ4g16SGvd/iCoTQGnhs44AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.stem(e)\n", "pass" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### The eigenvalues and eigenvectors give a factorization of the covariance matrix" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[128.578, -2.133, -6.661, 127.395, 63.27 , 42.369],\n", " [ -2.133, 95.051, 2.269, -1.349, 44.885, 31.637],\n", " [ -6.661, 2.269, 94.934, -5.048, -1.141, 30.275],\n", " [127.395, -1.349, -5.048, 126.993, 63.196, 42.723],\n", " [ 63.27 , 44.885, -1.141, 63.196, 54.408, 36.838],\n", " [ 42.369, 31.637, 30.275, 42.723, 36.838, 36.563]])" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "v @ np.diag(e) @ v.T" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Geometry of PCA\n", "\n", "![Geometry off PCA](https://i.stack.imgur.com/AaF1w.jpg)\n", "\n", "### Algebra of PCA\n", "\n", "![Commuative diagram](figs/spectral.png)\n", "\n", "Note that $Q^T X$ results in a new data set that is uncorrelated." ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "m = np.zeros(2)\n", "s = np.array([[1, 0.8], [0.8, 1]])\n", "x = np.random.multivariate_normal(m, s, n).T" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(2, 100)" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x.shape" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Calculate covariance matrix from centered observations" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [], "source": [ "xc = (x - x.mean(1)[:, np.newaxis])" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [], "source": [ "cov = (xc @ xc.T)/(n-1)" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[1.014, 0.881],\n", " [0.881, 1.163]])" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "cov" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Find eigendecoposition" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [], "source": [ "e, v = la.eigh(cov)\n", "idx = np.argsort(e)[::-1]\n", "e = e[idx]\n", "v = v[:, idx]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### In original coordinates" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAARMAAAEOCAYAAACq4fP9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xt83HWZ6PHPM5fMTDLJpLm0TS/p\nxdKSUmEpkcK64o1KV1SOF1zc3SPetrpHWQQ8XgDX3dUFXNdFZF2le/SsIojrC0FWkduKW3Ep0Pa0\n2DYQe6dtmiZNM8kkM5O5fM8fkwmTNJeZyS/5/WbyvF+vvkiamckTYJ58L8/3+YoxBqWUmi6X3QEo\npcqDJhOllCU0mSilLKHJRCllCU0mSilLaDJRSlnCMclERPwi8ryI7BaRvSLyt3bHpJTKnzilzkRE\nBKgyxkRExAs8A1xvjNlmc2hKqTx47A4gy2SyWmT4U+/wH2dkOqXUlByTTABExA3sAFYB3zLGPDfO\nYzYDmwGqqqouOvfcc2c3SKXmkB07dnQbYxrzeaxjpjm5RKQWeAi4zhizZ6LHtba2mu3bt89eYErN\nMSKywxjTms9jHbMAm8sY0ws8DWyyOxalVH4ck0xEpHF4RIKIBICNwEv2RqWUypeT1kyagO8Pr5u4\ngH83xvzc5piUUnlyTDIxxrwIXGh3HEqp4jhmmqOUKm2aTJRSltBkopSyhCYTpZQlHLMAq9Rc1NYR\n5rE9nRzvjbK4NsCmdQtoaQrZHVZRdGSilE3aOsJs2XqIcDRBU8hPOJpgy9ZDtHWE7Q6tKJpMlLLJ\nY3s6CQW8hAJeXCIjHz+2p9Pu0IqiyUQpmxzvjVLtH73SUO33cLw3alNE06PJRCmbLK4N0B9Ljvq7\n/liSxbUBmyKaHk0mStlk07oFhKMJwtEEaWNGPt60boHdoRVFk4lSNmlpCrH5shWEAl46wjFCAS+b\nL1tRsrs5ujWslI1amkIlmzzG0pGJUsoSmkyUUpbQZKKUsoQmE6WUJXQBVil1luyZIU/twuX5PkdH\nJkqpUXLPDJlUcijf52kyUUqNkntmqBCaTJRSo4x3ZigfmkyUUqOMd2YoH5pMlFKj5J4ZKoQmE6XU\nKLlnhsTtqcj3ebo1rNQMK8XWjNkzQzf2njyc73M0mShVgEITQ3abNRTwjmrNONHp4FJMPFk6zVEq\nT8X0bC2kNWOp94TVkYmak4oZAYytv8j+87E9nRM+93hvlKaQf9TfTdSaMfv6Q8kUzx3qIxJL4nUL\n9z57hNvec34xP+ascszIRESWisjTIrJPRPaKyPV2x6TKU7EjgGJ6thbSmvF4b5RYIsnOo73EEymC\nPjcYwzP7T5fE6MQxyQRIAjcZY9YClwCfFJG1NsekylCxXeGL6dlaSGvGxbUB2jr68Xlc+L1uRARE\nmFdZGh3rHZNMjDEdxpidwx/3A23AYnujUuWo2K7wxfRsLaQ146Z1CzgzmMAYgzGGWCJFPJlm7aLq\nkuhY78g1ExFZDlwIPGdvJKocLa4NEI4mRp09yacrfDYx5K61/Mnrlky51pJva8aWphBvWFXPnhN9\n9MeT1Pi9rFtcg9ftZn51Yedk7OC4ZCIiQeBB4NPGmL5xvr4Z2AzQ3Nw8y9GpcrBp3QK2bD0EZEYk\n/bEk4WiCP3ndkimfO9M9W//80mUjW8mFxmY3McbYHcMIEfECPwceN8b801SPb21tNdu3b5/5wFTZ\ncXI9h5NiE5EdxpjWfB7rmJGJiAjwXaAtn0Si1HQ4uSu8k2ObjGMWYIHXA/8TeIuI7Br+83a7g1JK\n5ccxIxNjzDOA2B2HmpucNLUoVY5JJkrZpdDzM05ld0J00jRHKVsUW8TmJE4416PJRM15xRaxOYkT\nEqImEzXnFVMm7zROSIi6ZqLmpNz1BZ9bOBGOsay+quQKxbKKreq1ko5M1Jwzdn3B63HjEmEomZry\n/IxTFXNuyGo6MlFzzrh9SeoqCQW83LBxtc3RFafYc0NW0mSi5pxCGhZNxu6t2LHsrpzVaY6ac3IX\nXLv6Yzx78DS//N1JjvYM5r2V6oStWKfRZKLmnOz6wqGuCDuP9NIXTeB2QVONL++E4IStWKfRZKLm\nnOz6wsn+OIl0mlDAS+vyeSxvCOadEJywFes0umai5qSWphDNdZVsWFGHS149EpZvQsjdiu2OxNh/\naoDuSJz6oI+2jnBJ7QRZRUcmas5p6whz55Pt7D0RZmt7F92R2MjX8q3NyE6VDndH2H74DOFoAq/L\nxcLq/KdK5UaTiSpZ2aTwmZ/s5s4n2/N6A+cunF6wJEQkluTZAz2c6o9OXZtx113wla+AMSNTpY6+\nOKk01AS8rF9Wy4rG/KdK5UaTiSpJxe6m5C6cLqgJsGFlHdV+D7tf6Zu8WO2HP4RPfxq++EV4/nng\n1anSH792IZeurKexOrPdPFfXTnTNZA5wWj2EFYq5EAvOrjFprPZz2WofHeHYxAVrv/gFfPjDmY//\n8R9hw4aRLzmhjN0pdGRS5sq1HqLY3ZSCD/X99rdw9dWQTMLnPgc33TTqy04oY3cKTSZlrlzrIYo9\n6VvQm/93v4N3vAOiUfjoR+H22896SCH34pQ7neaUOatKx52m2Osq8j7DcvAgXHEF9PbCu98N3/kO\nyPhdRe0uY3cKTSZlrlzn9NM52Dblm//kSXjb26CjA978Zrj/fvDoW2Uq+m+ozE3nwimnm5ERQTgM\nmzbBgQOwfj08/DD4/VM/bwrluAg+lq6ZlDmd0xcgGoV3vQt274ZzzoFf/hJqaqb9suW6CD6Wjkzm\nACfM6R3/mzmZhGuuga1bYdEiePJJmD/fkpcudhu71OjIRM04J/9mbusIc+cTL/PCxvfCI4+Qqq2F\nJ56AZcss+x5z5VCgjkzUjCvkN/NsjmDaOsJs+a+DXP3ju3jdrx9hyBfgX278BlfULaHFwu9Trovg\nY+nIRM24fH8zz/YI5rE9nbz9l/fyhw9/n5Tbw8//+pv0XtBqeQ3OXCls02SiZly+BWazXWC3+Kf3\ns/EHd2JEePyzX+XI6y6bkenHXFkE12mOmnH5bk/PaoHdQw/xvi1fBuDpT36Rl9/8DmDmph9OWASf\naY4amYjI90TklIjssTsWZZ18fzPP2mVYTz8N11yDK53msfd+nK0b31/W04/Z4rSRyb8B/wz8wOY4\nlMXy+c08KwV2O3fCVVfB0BB88pMsu/krtO09Zdv1EOXEUcnEGLNVRJbbHYcqnBW7MDN+90t7e6a6\ntb8/U1PyzW/S4nLRsqjWmtef4xyVTPIhIpuBzQDNzc02R6Pg1V2YUMA7ahemmEXGGVtbOH4cNm6E\nrq7MAb7vfx9cjprll7ySSybGmC3AFoDW1lZjcziKEqjw7OnJHNw7ehQuuQQefBAqKiZ9ip0Vu46v\nFp6ApmY1bY6u8BwYgCuvhH37YO3aTNe0qqpJn2Jnxa6Tq4WnoslETdus7cIUamgI3vte2LYtUx7/\nxBNQVzfl0+xsKFXKzawclUxE5EfAs8AaETkmIh+1OyY1NTsrPCfsUJ9Ow7XXwuOPQ2NjJpEsXpzX\na9o50nL0KG8KjlozMcZ8wO4Y1Gj5zN8L2YXJvt6+jjDhaJIav4fzFoWKWheYcOH3Dctpue1WeOAB\nqK7OtBJYPUGz6HHYeZamlM/xOCqZKGcpZJcmn12Y7Oul02mOnh5EROgbTFDpdbNl62DBuz8TLfz2\nfv6v4d5vZRZZf/YzuOiign5uOxtKlXIzK0dNc5SzWD1/z77eyb44fq+bUMCLz+viZH+8qNcdb0rw\nhid+zKX33p3Z9n3ggUzbxQLZeZamlM/x6MhETWiiszJ7T2TWKQrdusy+Xl8sQbUv87+ez+MiEksW\ntS4wdkqw5umf85Z/+Urmi/fck2kEXSQ7z9KU6jkeHZmoCY23S3Oke4BjZ6JFbV1mX6/G7yWeTAMQ\nT6YJDg/nC10XyF34XfrCVq74h88hxnDq5r+Bj32soNdS06fJRE1ovF2a9s4IaxYEi5r6ZF9vYY2P\nWCJFOJognkizsNpX1O5Pdkqw5vBe3vm31+FOJXl8059x3xuvKYm6jHKjyURNaLz5+9L6AM31o4u+\n8p2iZF9veUOQ5vpKagJeltQFWNEYLHpdoKXnGO//0sfxDcXY9Zar2PPpWwnHkiVT6FVOdM1ETWrs\n/P3OJ9tHrVN09cfYe6KPRMpw55PtU66fWLoecOQIXHEF7t5eXrr4zfz6f9+Gy+UiFMj8jpyonL9U\ny9WdTkcmqiC5U5/OvijPHeyhP5bkgqU1s1v6fepU5rzN8eMcaLmIx2+9E+P20B2Jse3gabYd7ObJ\nfZ1nxVLK5epOp8lEFSR36rP7WJig38Olr6ljfnVgVkq/2zrCfOvhnRy75E3Q3k54zVr+bvPtPNJ+\nhl+1dfLb/aeJJVL43C68bhmVKNo6wnzpkX28eKyXto4+egbiJVWu7nQ6zVEFy05Vslu9rpw7eGey\n9LutI8z3nnqJv7jjUyw51Max+kX8+ZVfpN7rx2sMHeEYyVSaCrcLl0tY3xyiwuMeSRRbth6iJzJE\nXaWXWCLFjiO9XLSslroqX0mUqzudJhNVtHxLv61ao3h893H+/O4vsHrfdnpq6vj0h++gN1hHf0+M\ni5fPo3sgjriEyFCSN61upLHaT9oYjvdGRwrm6oIVxBMp/F43APtPDdDS5C6JcnWn02mOKlo+B/ws\nW6Mwhku/disXvPA0/YEg13/oDn7nqyeRSpFIpTk9mKCpxo/XJQzEkxzoGqA7EuNI9wBHewZ5eNdx\n9p0IU1+ZqXGJJVJUuIXuSFz7vlpERyaqaGMP+FW4hUqvi+8+c3hkBGJZ46Sbb2bDrx4i5vVx/Z9+\nmfb5K/AOJUmk0qRNmuM9A3g9bmKJNIEKF9GhJL9+qYuUMWxYUQfG0BdL0hdLsqK+ktODCXoiQ9QH\nfSVTru50mkzUtGTXT3IPBdYFPSMjkP5Ygpam0Zd/F7yu8vWvwx13YNxuPvf+W9m5dC2Vrkwp/lAy\njdslhGNJlsyroC5YQZXPQzINybQhFPCyvCFI0O9hx5FeAE4PDLF2UYhwNKGJxEKaTJQlJhqBHO+N\n0h9LFn2k/sRd32bRZz6T+R433c7BJZfg7o0STaQJeN0EQ5lS/GgiRY3fw/rmWhqrM+eJnth7EkOm\ns2dD0M9Fy2r5fWeEV85EQYRQwDOyOKsJZfp0zURZYqKmPqGAp+jGSa98/wEW3HgdAL/6xM089/q3\nk0gaXrs4xOoF1YQqvSwMBXjDOQ28dkktaxeFRhIJQIXHhc/jHvncGAhHEwzGk8QTKeqrKrTOxEI6\nMlGWmGhnZ21TaGTtJLub87rltTy2p3PU2spZI4OtW2na/CHc6RTb/vQv2f2eawkBqxcEaT8V4ZKV\n9aP6fVx7aTNPtXUBr/YBqauqwCVCOJoglkjywqEznB4YYn61D4Bdr4RZ31xLKODl3meP0Fjt16rY\nadBkoooydrt39YKqs97M2aY+uSX0eTVc2rUL3vlOPENxdr/9T3j22utHvu+yhioGEylCAe9ZXd1W\nNgZHxfTZTWuAzBTsuUOnCfo9pDHMq6pAhmtj9ncNsLKhkv8+0MNbzp0/7as65jJNJnPUdGo/xksI\nT7V1cXlLI+2dA5O2bpxobSU7MujatYfPfuUvqO3r49n1b+aBD9zE8pyiuP5YkvMWhbhh49ltGCc6\n95NbYPf8oR5iw3Um2V4qbR39zKt08FUdJUKTyRw03UuzJkoI7Z0D477Jc43XcCmWSPLfB3p4Y3CI\n6+74FLV9PexcfRHf/Yu/4eXj/RiXm2UNVdNqYZidhq2aXzWyq4MxeN3CmcEEr181umt9qTRxdhJd\ngJ2DptuOcaLF1mwHtrM6xecYr+FSW0c/S4jxsds+SVNPB/uXtfDVzbfRJ17WN9dysj8+7RaG2QI7\nr9vNhc2Z55+JJjlvUQ1vWFWPzzP65ymVJs5OoiOTOWiidoz5/iaucAtb27tIpAxBv4f6Si8Hugbo\n7I/R1R+npal6wtHOeA2TB3v7ueu+m1l2/ADHFi7jjuu+jqkKEoklaV5Rh9fj5h+vvmBaP/PoArsk\nb1ozf2Rqlx2p5cZUKk2cnUSTyRw0nesU2jrCdPbF6Y8lCfrcnInE2HcijMBIgsrdJRm77jC2anZp\n0MM/P3Q7y1/aRfe8BXzpf32dSLCWeCJVdDvHiUy2pjKjF6bPEZpM5qDpXKfw2J5OltZVsjDkY/+p\nAU6EY/g9bpLpNLWVmV0SdypJT/shzll/7rijnZE3dfairJ2/IVJdyz03f4tj7gYkmgADy+oqZ22E\nUKpNnJ1Ek8kcNJ3fxK+2HfDSEPTTH0tSVeHi0OlB4sk0Abdw/Y++ymvbnufHX/5XFresG/+FjIEb\nb4Qf/hCqquh64EEqvUtozrmca0VjUOs9Sogmkzmq2N/EY6dIQb+HvmiChTV+4sk0FUNDzOvppK6v\nhw994Vo67/sJMM4Oz223wV13gdcLDz/Misvfwg3T/JmUvXQ3RxVkbNuBhdU+BuJJXtNYxYXNIYYq\n/Nx47e3s2fBWAoMRln/g3fCLX4x+kXvugVtvBRG4/364/HJ7fhhlKUclExHZJCIvi8h+Efm83fHM\nFRNe/j2OsR3rVzQG+fwfr2F5Q5BECt60Zj53f+QS1v328czdNbEYXHUV3HsvAMf+9Qek//IvAbj7\nfTdyi2u1nospE2KMye+BIg8D/wd41BiTtjwQETfQDmwEjgEvAB8wxuyb6Dmtra1m+/btVocyp+QW\nsOUuxlpSSm4M3HIL3H47AD0f/Ag19/0ATyrJD6/8GP9+xQeJxFOsaKjiM1es1rURBxKRHcaY1nwe\nW8iayQDwYyAsIv8G/F9jzO+LiG8iFwP7jTEHAUTkAeAqYMJkMhdZfU2DZc2LxiOSWRtpaICbbqLu\nB98D4JE3vo//eMeHCYggkul2pqXrpS/vaY4x5s+AJuDLwOXAyyKyVUQ+KCJWFAIsBl7J+fzY8N+N\nIiKbRWS7iGzv6uqy4NuWjpm4pmGialZLS8lvvBE+8QkM0FHXxP3vuy6TaMg0OIonU1q6XgYK2s0x\nxvQB3wa+LSLnAR8D7gG+KSI/Br5hjGmzPsxRMWwBtkBmmjOT38tpZmIUMZ0CtoJ8+9s8VdHElmV/\nyBAusvW38WQan6f4hs56oZZzFLUAKyKLyExB3gEkgQeBpcCLIvKZImM5PvwaWUuG/04Nm4lRxNjd\nmUNdEbYdOM2+4UXZ3FFPIQu141ny+esJzqshEksSG0oSHUrSH0vSEPQV1dBZL9RylryTiYh4ReR9\nIvIocAT4H8A/AE3GmI8aY94OvBe4tchYXgDOEZEVIlIBXAM8UuRrlaXxDslNdxSRuzvT1tFH+6kI\naxYGOXfh6Bv6rHjjtjSF+OymNWxYWcdQ2pBMw6Ur64pefJ3ugUVlrUKmOR2AAPcDnzfGvDjOY7YC\nZ4oJxBiTFJFPAY8DbuB7xpi9xbxWuZpOGfxksgVsdz7ZzpJ5leNOo7KfT3eK1dIU4rb3nD+teLOm\ne2BRWauQZHID8BNjTGyiBxhjeoEVxQZjjHkUeLTY55e7mT6QNtWbs9g3br7rGoWuf8zaeo/KS951\nJk6kdSbWuvPJ9rPenLmf537cHYmx53gfQ6k0b1u7cNIEkU8dS/ZxqVSak30xegYSeFzCdW99DVee\nf9am3rivfaR7gPbOCEvrAyO9Z3UxdnoKqTNxVAWsstdkN/Tlfu1Uf5RnD/QQiSW5YElo0vWTfNc1\nHtvTSSqVpv1UhHgyTV2VFwTu/tWBCddl8l3vUbNDk4kaMbZUPrezWe7Xdr/SR7Xfw4aVdSyoCUy6\n8JnvDtTx3ign+2L4PC78XjciQo3fQyKVnnRBtaUp0w/2vEUhLllZz/KGoC7G2kRPDatRJjtNnP3a\nq20IXm30PF6CaOsIc7RnkP939AwNQR+r5leNtC0Yu66xuDbAcwdPk0ynSaQMFR4XwQoP9VUVea3L\n6GKs/TSZqIJNtfDZ1hHm3meP8Mz+0/g8QiqdmTJtP3yGcxdW43K5ztqBWr2ginA0gUuEygoX8USa\nSCzGolp/XguquhhrP53mqLxli9b2ngiz7eBpDndHzlpbyS6K7j3Rx7yAh8oKz0gX+FQaOvri4x4i\nbO8cYP3SEG6XEE2k8XmE+mAFJ/vieRW0Tbbeo2aHjkxUXnJ3Tlqaaqj0utl9rJffHQ/j87q5cGkt\n8OqCayJlCPoyax9UVuD3unnTmjo6wrFxp1HHe6O8duk8Fs2rZH/XAJHhHrM1AW9eOzLax9V+mkxU\nXsaeCwr6PbjFRdDv4bLVjfTHkmzZeoj+WIKWphqCfg/xnMuu+mKJSacd2WlKY7V/5L7gsdOWqWgf\nV3tpMikjM3nobewC5/6uAYI+N0Op9MjuSfZx/bEkqxqr2Hk0c9mVMYYKt2ukWne8OGequlcPAs4e\nXTMpEzN96G3suaBILIkBavyvjhyq/R5CAQ/haIIKj5s/WJp50/ZGk6xbVMPmyzLF0ePFCUy4LV0s\nPQg4u3RkUiZmtMkRZ58L8rqF/liS1y559bX7Y8mRytPH9nQSiY++7AoyVbYTxXnDRmu7rc30vxM1\nmiaTMtDWEeaJfSfBQHXAy6rGKhqr/ZbWWYxd4DxvUQ2dfXG8bjdpY0ZNSyZbuyi2HqSY6YrWnswu\nTSYlLjuUr3C7MMYQT6TYebSX9c21VEyj6dB4xiaJsW/wfHZPiqkHKfaida09mV2aTEpcdii/bnEN\nO4704vMIPrew90QfKxuDM3obXjG7J8UstBY7XZmpRV01Pk0mJS73hr2LltWy/9QA4egQLnFZ0mF+\n7Ohj9YIq2jsHit4dKaYepNjpitaezC5NJiUudyjfEPTTEPRzuDtCR1+c7z5zeFrboWOnF4e7I/x0\n5zEuXFrLsoaqvKcbYxU6opnOdEVrT2aPbg2XuLFl5Ie7I+w82svCah8eF/z65VN8/N6d3PLTFwve\nEh3bPuBkX5wqn4eT/fFZPZmrpfKlQZNJiRvbNqCjL86FS2sJ+j3seiVMLJEiNpTk0T0nueHHu/nF\ni/n36B7bPqAvlqDa5yaSU28yG7sjk7VGUM6h05wykDuU/8xPdtMU8vPcoR6MMXRHhnC7wD3cLuDu\n/zzAysZgXm/EsdOLGr+XcDRBjQ27IzpdcT4dmZSZbKVqJJakP5ZpfSgIFV4XNX4PybTJe1py1iXl\nNZlLyhdW+3S6oc6iI5Myk90O9bqFwaEUPo+LZBoaq30j7RDHa2I0XkHY2N2Q5Q1B3nbeglG7Obo7\norI0mZSZbAK499kjHD49SCoNC2t8uF1CPJlmWV3lqGnJVAVh400vrpztH0qVBJ3mlKHs3TS3vfs8\nGoI+4kmDz+Ni9fwgbrdr1LREL7JSVtGRSRm78vzFrGwMTnqmRc+vKKtoMilzU+2CWHl+RXuHzG06\nzZnjrCoI094hSpPJHGdVQZiuvShHTHNE5Grgb4AW4GJjjN75OYusKAjTtRfliGQC7AHeA9xjdyD5\ncMLagBNiyKW9Q5QjpjnGmDZjzMt2x5EPJ6wNOCGGsfQwnnJEMimEiGwWke0isr2rq2vWv78T1gac\nEMNYehhPzdo0R0SeAhaO86VbjDE/y/d1jDFbgC0Ara2txqLw8uaEtQE7Y5hseqWH8ea2WUsmxpjL\nZ+t7zSQnrA3YFUOxvVjV3FBy0xy7OWFtwK4YnDi9Us7hiGQiIu8WkWPApcAvRORxu2OaiBPWBuyK\nYWyzJNDtX/UqR2wNG2MeAh6yO458OWFtYDZjyK6T7D0R5ved/axbXENDMLNmo9u/KssRyURZy8qO\n8rnrJBcsCfHCoTM8e6CHDSvn4fN49OoINUKTSZnJt6P85S2NeSWY3HWSUMDLhpWZO3l2v9LHxrUL\ntDmSGqHJpMz88NkjHOyKMJRKU+P3MhhPjnSUX9EYJBTwcmYgzt3/eYBLXlM/5a5Mdhu6qz/G/q4B\nIrEkQZ+bxmovN2xcbdNPqZzIEQuwyhptHWF+s/80xhiqfR5iiRRHegZxiyESS9IdibHt4Gme2d9N\nZ3+MRCo15a7M4toAR7oH2Hm0l3giRdDnpi+W5NiZqJ4IVqNoMikjj+3pZF6lFxFBRPB73QS8bjr7\nMh3qdxzpJZZIkU6D1y3sONJLdyQGTLwrs2ndAto7IwD4PC7iyTQAaxYEdUtYjaLJpIwc743S0lRN\nPJkmlkhhjKEm4GEwkSKWSFHhzlx34XJlRiM+j4v9pwaAiXdlWppCLK0PUOP30B9P4ve6uWhZLc31\nVbolrEbRNRObzMSp32xl7Prm2pH1jcoKD288p4HDPYMkUobqgIeLl8/j0OlBMIZwdGik6G2iXZm1\nTaGzKm7D0YRuCatRdGRig5k69ZutjK3wuNmwoo6LV9SxsjHIX11+Dm9bu5ANK+u5dGU9qxfWsL65\nFkRwiWvKojcnVP0q59NkYoOZKkufrDJ2bEKo8LhZ2Rjk6+8/nxs2rp50VOSEql/lfDrNscFMnvqd\nqDJ27IVahV6g5YSqX+VsmkxsYNepX00IaibpNMcGugahypEmExvoGoQqRzrNsYlOOVS50WQyhzmt\nw70qbZpMSpAVSWC6LRg1EamxdM2kxFhV8DadWhcnXrWh7KcjE4cbOwLo7o+NvPGBkX8+tqezoJHB\ndGpdchPRdGJQ5UWTiYONNxX5zf7TvH5VHZB5A3dHYvy+M0Jnfxwg7+nGdGpdnHDdh3IeneY42HhT\nkXmVXvad6AcyiWTHkV76YkkWVPsKmm5Mp9ZlcW2A/lhy1N9pL1ilycTBxusG39JUzZnBzBv/98N9\nRgBWzQ8WtO4xnVoXLbpT49FpjoONNxXxez380ap6QgEvnf1xFlT7WDU/SGN1ZtpRyHSj2FqX6Z7z\nUeVJk4mDbVq3gC1bDwGZJNEfSxKOJkaNIOy6XVCL7tRYOs1xsKmmIjrdUE6iIxOHm2wEoNMN5SSa\nTEqcTjeUU+g0RyllCUckExH5moi8JCIvishDIlJrd0xKqcI4IpkATwLrjDHnA+3AF2yORylVIEck\nE2PME8aYbEnlNkBvwlaqxDgimYzxEeCXdgehlCrMrO3miMhTwMJxvnSLMeZnw4+5BUgC903yOpuB\nzQDNzc0zEGmG9utQqjBijLE7+FLwAAAFlklEQVQ7BgBE5EPAx4G3GmMG83lOa2ur2b59u+Wx5J7W\nnajyVKm5QER2GGNa83msI6Y5IrIJ+CzwrnwTyUyaqUuylCpnjkgmwD8D1cCTIrJLRL5jZzDjndbV\nfh1KTc4RFbDGmFV2x5DLrkuydJ1GlTKnjEwcxY4DdNpXVZU6TSbjsOOSLF2nUaXOEdMcJ5rtA3Ta\nV1WVOh2ZOIT2VVWlTpOJQ2ijI1XqNJk4hF5mrkqdrpk4iDY6UqVMRyZKKUtoMlFKWUKTiVLKEppM\nlFKW0GSilLKEJhOllCU0mSilLKHJRCllCU0mSilLaDJRSllCk4lSyhKaTJRSltBkopSyhCYTpZQl\nNJkopSyhyUQpZQlNJkopS2gyUUpZQpOJUsoSmkyUUpbQZKKUsoQjkomIfFlEXhSRXSLyhIgssjsm\npVRhHJFMgK8ZY843xvwB8HPgr+0OSClVGEckE2NMX86nVYCxKxalVHEccwmXiPw98EEgDLx5ksdt\nBjYPfxoXkT2zEF6+GoBuu4PI4bR4wHkxaTyTW5PvA8WY2RkEiMhTwMJxvnSLMeZnOY/7AuA3xnwp\nj9fcboxptTDMadF4pua0mDSeyRUSz6yNTIwxl+f50PuAR4Epk4lSyjkcsWYiIufkfHoV8JJdsSil\niuOUNZM7RGQNkAaOAJ/I83lbZi6komg8U3NaTBrP5PKOZ9bWTJRS5c0R0xylVOnTZKKUskTJJxOn\nleKLyNdE5KXhmB4SkVqb47laRPaKSFpEbNtyFJFNIvKyiOwXkc/bFUdOPN8TkVNOqVMSkaUi8rSI\n7Bv+73W9zfH4ReR5Edk9HM/fTvkkY0xJ/wFqcj7+K+A7NsfzNsAz/PFXga/aHE8LmcKjXwOtNsXg\nBg4AK4EKYDew1uZ/L5cB64E9dsaRE08TsH7442qg3c5/R4AAweGPvcBzwCWTPafkRybGYaX4xpgn\njDHJ4U+3AUtsjqfNGPOynTEAFwP7jTEHjTFDwANkSgBsY4zZCvTYGUMuY0yHMWbn8Mf9QBuw2MZ4\njDEmMvypd/jPpO+tkk8mkCnFF5FXgD/DWYcEPwL80u4gHGAx8ErO58ew8Y3idCKyHLiQzGjAzjjc\nIrILOAU8aYyZNJ6SSCYi8pSI7Bnnz1UAxphbjDFLyVTPfsrueIYfcwuQHI7J9nhUaRCRIPAg8Okx\no+5ZZ4xJmcxJ/iXAxSKybrLHO6VobVLGYaX4U8UjIh8C3gG81QxPOu2MxwGOA0tzPl8y/Hcqh4h4\nySSS+4wxP7U7nixjTK+IPA1sAiZcsC6JkclknFaKLyKbgM8C7zLGDNoZi4O8AJwjIitEpAK4BnjE\n5pgcRUQE+C7QZoz5JwfE05jdiRSRALCRKd5bJV8BKyIPktmtGCnFN8bY9ltPRPYDPuD08F9tM8bk\nezxgJuJ5N3A30Aj0AruMMVfYEMfbgW+Q2dn5njHm72c7hjHx/Ah4E5kj/53Al4wx37Uxnj8CfgP8\njsz/ywA3G2MetSme84Hvk/nv5QL+3Rjzd5M+p9STiVLKGUp+mqOUcgZNJkopS2gyUUpZQpOJUsoS\nmkyUUpbQZKKUsoQmE6WUJTSZKKUsoclEzajhsuwOEflSzt+dLyIxEbnaztiUtbQCVs04EbkC+A/g\njcAuYDvwvDHmw7YGpiylyUTNChH5BvAu4L+ANwB/kNN8R5UBTSZqVoiIj0y7xnOAP5yq0Y4qPbpm\nombLcjI9TQyZXrCqzOjIRM244aY/28g0SX6OTPOqC4wxR20NTFlKk4macSJyB/CnwPlAmExfXD/w\nFmNMerLnqtKh0xw1o0TkjcBNwAeNMb3DbSw/BKwFPmdnbMpaOjJRSllCRyZKKUtoMlFKWUKTiVLK\nEppMlFKW0GSilLKEJhOllCU0mSilLKHJRCllif8PdnXBp/oNSoEAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.scatter(x[0], x[1], alpha=0.5)\n", "for e_, v_ in zip(e, v.T):\n", " plt.plot([0, e_*v_[0]], [0, e_*v_[1]], 'r-', lw=2)\n", "plt.xlabel('x', fontsize=14)\n", "plt.ylabel('y', fontsize=14)\n", "plt.axis('square')\n", "plt.axis([-3,3,-3,3])\n", "pass" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### After change of basis" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [], "source": [ "yc = v.T @ xc" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAARMAAAEOCAYAAACq4fP9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3X10XHd95/H3dx40I2v0EMuyLT/G\njhPHxk1JEJAcKLCBFLdQKOzh0JZDy3bPettuw8OWw1KyuxzawznQcJp2U/Zw3IVCF2g2LKRwCASc\nhK5JSQJOSIJtESPZiYkty7JsjTSyZjQP3/1jZpTRZCTNw2/m3pG/r3N07BmNRt+5uvdzf/d3f/d3\nRVUxxphGBbwuwBizOliYGGOcsDAxxjhhYWKMccLCxBjjhIWJMcYJ34SJiERF5Mci8rSIHBORT3hd\nkzGmeuKXcSYiIkCXqiZEJAw8AnxAVR/zuDRjTBVCXhdQpPlUSxQehgtf/kg6Y8yKfBMmACISBJ4A\ndgGfVdXHK7zmAHAAoKur6xXXX399a4s05gryxBNPXFDVgWpe65vDnFIi0gfcB9yuqkeXet3Q0JAe\nOXKkdYUZc4URkSdUdaia1/qmA7aUqk4BPwD2e12LMaY6vgkTERkotEgQkU7gNuDn3lZljKmWn/pM\nBoEvFfpNAsC9qvptj2syxlTJN2Giqs8AN3pdhzGmPr45zDHGtDcLE2OMExYmxhgnLEyMMU5YmBhj\nnLAwMcY4YWFijHHCwsQY44SFiTHGCQsTY4wTFibGGCcsTIwxTliYGGOcsDAxxjhhYWKMccLCxBjj\nhIWJMcYJCxNjjBMWJsYYJyxMjDFOWJgYY5ywMDHGOGFhYoxxwsLEGOOEhYkxxgkLE2OMExYmxhgn\nfBMmIrJVRH4gIsdF5JiIfMDrmowx1fPNjcuBDPBnqvqkiHQDT4jIIVU97nVhxpiV+aZloqpjqvpk\n4f8zwDCw2duqjDHV8k2YlBKRq4Ebgce9rcQYUy3fhYmIxICvAx9U1ekK3z8gIkdE5MjExETrCzTG\nVOSrMBGRMPkg+YqqfqPSa1T1oKoOqerQwMBAaws0xizJN2EiIgJ8HhhW1b/2uh5jTG18EybAa4D3\nAreKyFOFr9/0uihjTHV8c2pYVR8BxOs6jDH18VPLxBjTxixMjDFOWJgYY5ywMDHGOGFhYoxxwsLE\nGOOEhYkxxgnfjDMxV67hsTgPHB3nzNQcm/s62b9vA3sGe70uy9TIWibGU8NjcQ4ePkV8Ls1gb5T4\nXJqDh08xPBb3ujRTIwsT46kHjo7T2xmmtzNMQGTh/w8cHfe6NFMjCxPjqTNTc3RHFx9td0dDnJma\n86giUy8LE+OpzX2dzCQzi56bSWbY3NfpUUWmXhYmxlP7920gPpcmPpcmp7rw//37NnhdmqmRhYnx\n1J7BXg68bge9nWHG4kl6O8MceN0OO5vThuzUsPHcnsFeC49VwFomxhgnLEyMMU5YmBhjnLAwMcY4\nYWFijHHCwsQY44SFiTHGCQsTY4wTFibGGCcsTIwxTliYGGOcsDAxxjhhYWKMccJXYSIiXxCR8yJy\n1OtajDG18VWYAF8E9ntdhDGmdr4KE1U9DFz0ug5jTO18FSbVEJEDInJERI5MTEx4XY4xpqDtZlpT\n1YPAQYChoSH1uBxfsZtZGS+1XZi0QjtulMWbWfV2hhfdzMpv86m247I11Wm7w5xma9c7zLXDzaza\nddma6viqZSIi/wS8AVgnIi8AH1fVz7eyhtKNElj494Gj477eg56ZmmOwN7roOb/dzKrWZWutmPbi\nq5aJqv6uqg6qalhVt7Q6SKB97zDXDjezqmXZWium/fiqZeIHm/s6ic+lF/aa4L+NspL9+zZw8PAp\nIL+BziQzxOfSvPuVWzyu7EW1LNt6WojWkvGWr1omftCud5hrxs2shsfi3HXoBB/+2tPcdehEw62C\n4rI9NZHg0dEL3P/MGI+NTnLdhq6XvLbWFqK1ZLxnLZMyxY2ydA/37lduaYs9nMubWTXj7NCewV7e\ntGeAux8eJZ3N0d/VwWBvlAeHJ9g5EFv0vrW2ENu1r2s1sTCpwO4w17yN88T4LDfv7F8UEvG59Eve\nt9bDtnbogF7tLExMRc3aOKt931pbiO3a11Wq3ft8LExMRc3aOGt531paiO3QAb2cdhl0uBzrgDUV\nNasjulnv24wO6FZqh0GHK7GWiamoWR3RxU7YLz16mvHpJBt6ovzBLducbPTt3Ne1Gvp8LEzMkpqx\ncQ6PxXlweIK9gz28esdaZpKZimdz/KJV/Riroc/HwsQsaMWG0+pTuI18plb2Y7R7nw9Yn0lbcj2Y\nrPierRj01crLFRr9TK3sx2j3Ph+wlklbKN27RoLC2XiS7f1dTveWrWox1Nucr6eF0ehnanU/Rjv3\n+YCFie+VN7UPn5ggkcww2BslIGFnG72LDWe5Db74veNjcX45Ocd1G2JsX9dVVXO+uAxyuRxj8SQ/\nPX2J7x07x+23XsNbbtjctM+0GvoxWsnCxOfK967z2RyxSJCRiVkGuvMbiou9ZaMbznL9C8DC967f\n2MOacJBnxxNcTmd52abeFc8SPXB0nFwux7PjCSKhAP1dHUwnM9z90OiyHbeNfiYX/RjVBGy9/Tl+\nG+BmfSY+V97H0BMNo0CiZLoBF3vLesZ/lPbdfPxbx8lmcxX7F8r7Hq5eF+Pmnf28bFMvH7rtuhU3\ngjNTc4zFk0RCAaLhICJCTzREJqfL9l80Oqal0X6M5fpsGunP8etFjdYy8bnyveuu9V08OnqR7miI\nnKqTXv/iXi6RSnNmao6eaGjFFkN5S+Snpy8xfTlNLBqq2GJa6nCjmj3s5r5Ofnr6Ev1dHQvPpTI5\n1naFl22RuRgr00g/xnJ9NsXH9fTn+PWiRgsTnytvaoeDQXas62JDT4SxeLLhwWSloXD9xp6FcFqp\n2Vy+Qq+LRYjPpRcdfpW2mCodbkSCUtWp1/37NvC9Y+eYTmboiYZIZXKkMjm2r12zYovMy07Nlfps\n6u3P8esANwsTn1tqxOhyHY+1qHcvV75C93eFGTk/w5lLc6DKxp4owWBgocVUqe+hMxyo6nfvGezl\n9luv4e6HRrk4m2ZtV5jta9cQDAZ8Pc/MSn029fbn+LVj2MKkAj91bjV7xGi9e7nSFXpiJsnJC5fp\n6QwzM5fmF+cTPDue4DU71wJLH258/pHn6I9VN+bkLTdsZudAzDd/l3Ll68x1G7q4MJPkhyOTXLUm\nzJ7BbqLh0KJD0no7d/06wM3CpIzfrt584Og42WyO42PTJJIZYtEQG7sjzo6P693Lla7QI+cTAKhC\ndzRMT2cYVDkTTy5adpX6Qmr53X4dh1G+zpyaSPCNJ1/gpm19vGbXWo6fneFHoxd57a7+RetRvf05\nfp3Ay8KkTDM7t+pp8Rw7G2fkfIK5+QyZnDKZEM5PJ7mczjZUS9Fye7nl6i1docdnUmzojhAOCIGA\nEA0HUVUSqezCGZ1Kn9Ove9hala8z52ZSdEVCnJtOcfPOftbvfjE0yw/fGpm1zuvwKFf1qWEReY2I\n3CkiHxORrWXfu0pEHnZfXus1a7h3vafzxuNJJhPzgBAJBQFhMjHPeDzZUD1FS53+BFasd89g/tTu\nb798M3s39ZJRJRLKr1KpTI5YNLTsslsNQ8jhpetMIpmhOxJkOpleeM4PHaTNVlXLRER+C7gPeALo\nBv6LiPyuqn6n8JIO4PXNKbG1mtW5VW+LZ3Y+Q1AWPxeU/POuVNrL3XXoRNX1FlsYHcEAyXQWESGV\nyfGyTT2Ll50UPoi+eFdXP+5ha1W+zsSiIaZ90kHayv6/alsmdwB/oaqvVtW9wMeAe0XkHU2pykPN\nmryn3hZPJBxkoCdCMCjMZ3MEg8JAT4RIONhQPSuppd5iC2Pfph6m5vIh9/KtvXSEgksuu2ZcrOiV\n8nVmY3eE2VSGjT0RT+9w0OrBbdX2mewFfq/4QFU/KyLngC+LyO8DjzSjOC80q3Or3hbPjVv7ePzk\nRQZiESKhAKlMjkQyw41X9zVUz3KGx+KcvniZn56+xLpYhF3ru1gXi67YOfrJd96waE/Y2xmuuOxa\n2cndij1z+TqzYyDGm/dt4MT4rKcdpK0e3FZtmCSBtcDJ4hOq+nUREeAfgY86r8xDzWh619vZ+N5b\ntjMWT3Jxdp6ZZIaOUIBt/Wt47y3bndZXVNzQN3ZHmL6c36Meee4S12/sJhAIrFhvNcuuVSt5K0Or\n0ud+i9PfULtWD26rNkx+CtwKHCl9UlX/r4gEgC+7Lmy1qbfFs2ewl4/s392UvWulvXbphh6LhhiZ\nmOViYp6x6RSfeNteT8e21Mqvw85bpdWD26oNk8+xRAerqt5bCJT/2GgxIrIf+FsgCPwvVf1Uo+/p\nJ/W2eJo1fWKlvXYileb6jT0ADHRHGeiOklNlLJ50VkNxJZ/PZBmZmCWRzBAOCi/b1OPk/Yv8Ouy8\nVVp96r2qMFHV+8ifzVnq+/cA9zRSiIgEgc8CtwEvAD8RkW+p6vFG3nc1ctEPsNRe+8zUHDPJTFP3\nZvv3beCvHniW05OXiUWChAL53zE+nWJ4LO48tFx+Fj+Njl5Jqwe3VXtqeAD4I+BvVXW67Hu9wPuB\n/6mqkw3U8ipgRFVPFt73HuDtgIVJCVf9AKV77YmZJCMTs8zMpZlLZ+iOhGDtmqbM4QGwZ1Mf/1D9\nR67bh5rwnnsKX83y4XufchpSrTz1Xu1hzgeAq8uDBEBV4yJyLfBB4L81UMtm4Jclj18AXl3+IhE5\nABwA2LZtWwO/rj3V0g+w3EZdeqjx5OkpIqEA4aDQEeogp0o6k2Usnql5b1YMu2w2x7npJE+dnuJ7\nR89x+xuv8bxDsh344RKOelUbJr9FPlCW8gXgf9BYmFRFVQ8CBwGGhoZ0hZf7gsumcbX9ACu1YIrH\n0ycnEkQKo+Lms8ortvcSDgbp7Qzzoduuq7m+4rVEJ87nZ0Vb2xXOz4r28Cg7z04tfO67Dp0gPpdm\neGyaZDpLNBwkmc4SCQfZO9jD8bFp9g72vOSexL2dYX4xPs2PRi8iAqGAkM0p85kcV6/r4o/fsJMH\nhyfo7Qwvalm52DA//LWnC9NlvjiKsNif9Jl3/Wrd71tcFr2dYQK0pqO4GYdr1Q5auwYYXeb7J4Ed\nDVUCZ4DSYfpbCs+1lOvBVPc/c4YP/Z+n+fYzZzk9OctzFxINDRza3NfJTHLx6NdK/QArzaxePJ5O\nZ5VUNkc0HOQV2/tYF4tW1Um51HI6MzXHuemXzoqWzuYWzYpWHOh1IZGiIygk01lSmRy7BrrojoYY\nn05WHDR37Gycfx2dJFAIkulkhtn5LOFg/nff/dAouVzlGd8aVe2yr1UrZ+yHwjp5z9Pc/8wYpydn\nOTXR2DpZVG2YpFm8oZfbAjQ6vvsnwLUiskNEOoDfAb7V4HvWxPWIweGxOHc/lM/g/q4OUpn8PKa5\nXK7ulbvaEbqlK+iFRJLHTk7y2MkLHDo+vvB59gz2ctveDdy8cx037+xnXezFSY0iQVkyVJdbTpv7\nOrk4m164Rgfy1+n0d3Us2jiKYdYfi3DpcoZIOMhN2/oY6M4PjtvQE6244U4nMwhCqBBAoYAQCkAq\nqwQDQiaXbymUcrVhNmN0dHGA4Hd+NsZjJye5kMjX3qxTuMNjce5+eBQE1naFSWXyrchstv51sqja\nMHkSWG7o/L8lPxalbqqaAf4U+B4wDNyrqscaec9aub5PygNHx8nklJ5oCJH81bS5nPLj5y7yz0+d\nqavlU+3FccW96IVEkieen8ofQgTz/SKlAVlpA/nlxcucjScXwuLURIIP3fM0/+Eff8Jdh07w5Uef\nX3I57d+3YaHFoKoLLY7B3uhLNo49g7184m17+ZUtvewd7KG/MFtbfC7NH9yyreKG29sZYnNflPms\nMp/NAYoCmayyoSdKNCz84nyC7x8/t7BxutowXV+YWDpAMBwILAwQfO5ComnD7x84Ok46m1u0TkZC\nAc5NJxsO3Gr7TD5L/lqcF4C/U9UsgIiEyAfA+8m3JBpSuHDwOyu+sEmW6o84djbfpK/1+PLM1NxC\n+kfDQWZTGSZn50lnc2y9ak3dHW3V9NCX9ol0FPpEUlnlpm35a2aKx+OVTh/O90ToCAUXJj46cT4B\nAvHCiNgfjkzyml1rgRf7M4p7/z2Dvdz+xmu4++FRJmfn6e/q4Or+NQQClWdFW+70ZflkSO9+5RYe\nODrOmsKyTGeyZHPFOFEuzaaYnssQEOgICHPzGR4dvciOdV28+8219//Uu+yr1YoBguXOTM0ttJKj\nheu7IqEAF2fT3HJNY4Fb7TiTb4jIp4G7gL8UkWL/yTVAF3Cnqn69oUp8oNK4hOcvzPLCpTm2XLWm\n5p72zX2dzKez+Y0RmJxNkc0p4WCAXetjTe1oK26k//nep0GhuzPEyzb1LAxCKz/kKP39H/7a0wsz\noI1MzBIJBYiEAsyk8uNPrloT5vjZGdbvfnHlK9371zor2lIb6FLPHzx8mV/d2sfo+QQnJ2dJZ2Bj\nT4RUVgkEhK6OADkgm8uH3IaeiC/PipTuvJo1QLDc5r5O0pksz47n18lIKMB0MkMoIA23hKqeHElV\n7xCRfwbeA1wLCPD/gK+q6o8bqsInKo0YPDGeYPfGFzf8WgIg/36XuW59jHPTSRLJLB0h4dU7rnJ6\nz5ul7Bns5df3bqx54FZpqCaSGWKRIKlMjp5ouPC+3fxo9CLxufSSY1GaNb6htCUTDgVBhFQmSzgY\nZG4mydarOgkF852/N+/sX9g4/ciLuVyL6+TuDTHG4kkmZ+cJBwPc/sZrGv57VTtobQ1wJ/Db5Nu2\nDwG3q+qFhn67z1Rqcm/t72Rbf9ei11UbAKXv1xHOr/hrwgEmZ9OcOn6OnmiYjT0Rrl4Xa9ZHqmtI\ndenPxCJBpgsdofs254e7R8MhXrurn97OsCdXxZYGVenp2kdPTpJKZwt72/zERH6YaHkpXsw0Vx7G\nt1yzztkAOVFdeaiGiNwJ/AnwFWCO/HQE/6Kq72q4ggYMDQ3pkSNHVn5hA+46dILnLiQ4N51iOple\nFAC1jsO4/5kzfOq7z9IVCdEdCXIhMc9kYp69m7p59U53f9Ry9YwpKP7MsbNxXrg0x+4NMbb1d604\nbqPVw81Lx2hMzCT50cgFZlJZRPJ7/rVdHXxk/272DPbWXVszP5Pfh+eLyBOqOlTVa6sMk1HgjsI1\nOIjIq4B/BaLFzlgvtCJMygNgJpVlNpXho7+xu+bbTdx16ASnJhKcm0lxYSZJfC4/vd/6nih7N/U6\nG1zlWrUrfOlAOdcDxparrfg7k+kMj/xikkQqw1VrQvR0dhAJ5funzs+kagpFLz+Tn9QSJtX2mWwF\nflh8oKo/FpEMsInFQ+A95zrpT4zPcuPWPs7NpEgkM/R0hrlufYwT47M1Dw8/MzXH9nVd7BiI8ejJ\nSbqj2UUdm9Del8f/70ef5+REgnRWiUVD7BroWnZCaRdKm+2Pn5qkP9bB63evY10syoVEkkdHL5LK\nTC+Me3l2PEEsGloYU7NSbVf6NAa1qDZMgsB82XOZGn6+JZoxGU5pABSVnw2p1kodm368PL7aZTo8\nFueRkUmu6gzlP1c6f83Py7f2kki5m6+2kmIfSvHsSHG4+8j5WWKRIOmsks5mF+4GOHJ+tuqRvlf6\nNAa1qDYMhPwUjamS56LA34vI5eITqvo2l8XVqhl7EZc97it1bPqxs7CaZTo8Fufj3zrOTDJNKp1l\noDtCVyRU+N4Mb9i9viW1lv+tppNpOgJCrDASuJ7OWb/ePc+Pqh0B+yXgLDBZ8vVl8oc4pc95qhnX\nOLgcQl06grKnsHLu3hBjbVekpvdt5WTMKy3TYsvlYmKeTb0RUpkcZy7NkUilUVUuXW7dRMrXbeji\nsdFJ7n9mjEdHL5DJ5kiksuwa6GLXQBepTI7pZH6KhWqXd7MmGF+Nqh209u+aXYgLzdiLuJ5gpvS0\nZnn/zlKTL5ffdrJ4VexKhx0u+o5WWqbFlsvaWAepdJatazs5P51ifDrF5r5Ofm1Xf0v6Foq3Ud29\nMT9+Yiw+x0wyS0CUY2en2TPYXejrStC7JrzkZNflmj3BkN/P5tSiqrM5flV+Nme19bxX+jyPjU6y\ne2Ns0diU4sZePFXtcjms9F7FcR6TidTCvCgdQeHi5TQ3bOlr2bIvPUVcvB4JIBSAYCDApctpXrur\nn/fest0360I7rK+1nM2p+o5+7cD1hVheq3ThYTVXxbq8YHGlZVq8oHCgO8pN2/qIhINcupyhPxap\natm7OmQrPRwbOZ+/BKAnGiKnwht2r+fW69cz0B311brg+sJSr/nqbIwLrZymbiWNNGGHx+J8//i5\nwnU1YXYNdDHQHWVtV5jJ2cUn1soP5VyfgVhumZZ2KvfHIgs33qo2SFydfSs9HJtOpumOhBZuUQr+\nPAOz2s4UraqWiZ80MjdK8Wc7ClMGFE+zTswk2dgTJRwMLNsh2KxJfCpppDXocs9c2lHaHQkxXThk\nmE1lOHR8nMMnJhaunvaLVv6dWmHVtUz8opHT1MWf3be5hyeenyISEiJB4djZaXYOxLj91muWvVtc\nq6/5qLc16HLPXNpR2rsmzPmZJJmMEoqCSnNmv2+UF9fmlHLd+Wth0iSNbCgvDr4K84rtfYycnyU+\nN09AAgt7/eVG3zb7DIQrrs++lYbaHd94hqNnp5nP5gcF/sqW/Ny2fhq56uXfqRkDPC1MmqSRDaX0\nZ9fFoqyLRRceV/uH9lPfUSXDY3EuzCT54cgkV60Js2ewm2g45GzPnMoqr7tu4CWTP/utP8Krv1Mz\nBnhan0mTNDLYabUPlCruFcOhYGHGNvjR6EXmM1lnZ99WW3+Ea80Y4GktkyZppAnbLocpUN9x9+K9\nYpj1uztrbnmtxOv+CL9rxgBPC5MmaqQJ6/fDFKj/uLsVp0TbKZC90IywtTAxdav3uLtVF8+1QyB7\npRlha2Fi6lZvC8MOQfzBddhaB6ypW72dnKvtsgeTZy0TD7X7FaONtDDsEGT1sZaJR1zfitQL1sIw\npaxl4pHVMrdorS2Mdm+NmaVZmHhktV0xWo1mDOFuV6sxVH1xmCMi7xKRYyKSE5GqJmJpd1fiCM3V\nNn9HvVbDIW4lvggT4CjwTuCw14W0ymofMl9JM4Zwt6PVGqq+CBNVHVbVZ72uo5WuxM7LK7E1Vslq\nDdW26zMRkQPAAYBt27Z5XE1jrrTTozZYLW+13j6jZS0TEXlQRI5W+Hp7Le+jqgdVdUhVhwYGBppV\nrmmCK7E1VslqPcRtWctEVd/Uqt9l/OtKa41VslovQmy7wxxjVoPVGKq+6IAVkXeIyAvALcD9IvI9\nr2syxtTGFy0TVb0PuM/rOowx9fNFy8QY0/4sTIwxTliYGGOcsDAxxjhhYWKMccLCxBjjhIWJMcYJ\nCxNjjBMWJsYYJyxMjDFOWJgYY5ywMDHGOGFhYoxxwsLEGOOEhYkxxgkLE2OMExYmxhgnLEyMMU5Y\nmBhjnLAwMcY4YWFijHHCwsQY44SFiTHGCQsTY4wTFibGGCcsTIwxTliYGGOc8EWYiMidIvJzEXlG\nRO4TkT6vazLG1MYXYQIcAvap6g3ACeDPPa7HGFMjX4SJqn5fVTOFh48BW7ysxxhTO1+ESZk/BL7r\ndRHGmNqEWvWLRORBYGOFb92hqt8svOYOIAN8ZZn3OQAcANi2bVsTKjXG1KNlYaKqb1ru+yLyPuCt\nwBtVVZd5n4PAQYChoaElX2eMaa2WhclyRGQ/8BHg9ap62et6jDG180ufyd8B3cAhEXlKRD7ndUHG\nmNr4omWiqru8rsEY0xi/tEyMMW3OwsQY44SFiTHGCQsTY4wTFibGGCcsTIwxTliYGGOcsDAxxjhh\nYWKMccLCxBjjhIWJMcYJCxNjjBMWJsYYJyxMjDFOWJgYY5ywMDHGOGFhYoxxwsLEGOOEhYkxxgkL\nE2OMExYmxhgnLEyMMU5YmBhjnLAwMcY4YWFijHHCwsQY44SFiTHGCQsTY4wTvggTEflLEXlGRJ4S\nke+LyCavazLG1MYXYQLcqao3qOrLgW8D/93rgowxtfFFmKjqdMnDLkC9qsUYU5+Q1wUUicgngd8H\n4sC/WeZ1B4ADhYcpETnagvKqtQ644HURJfxWD/ivJqtneburfaGotqYRICIPAhsrfOsOVf1myev+\nHIiq6sereM8jqjrksMyGWD0r81tNVs/yaqmnZS0TVX1TlS/9CvAdYMUwMcb4hy/6TETk2pKHbwd+\n7lUtxpj6+KXP5FMishvIAc8Df1Tlzx1sXkl1sXpW5rearJ7lVV1Py/pMjDGrmy8Oc4wx7c/CxBjj\nRNuHid+G4ovInSLy80JN94lIn8f1vEtEjolITkQ8O+UoIvtF5FkRGRGRj3pVR0k9XxCR834ZpyQi\nW0XkByJyvPD3+oDH9URF5Mci8nShnk+s+EOq2tZfQE/J/98PfM7jen4dCBX+/2ng0x7Xs4f8wKN/\nAYY8qiEIjAI7gQ7gaWCvx8vldcBNwFEv6yipZxC4qfD/buCEl8sIECBW+H8YeBy4ebmfafuWifps\nKL6qfl9VM4WHjwFbPK5nWFWf9bIG4FXAiKqeVNV54B7yQwA8o6qHgYte1lBKVcdU9cnC/2eAYWCz\nh/WoqiYKD8OFr2W3rbYPE8gPxReRXwLvwV8XCf4h8F2vi/CBzcAvSx6/gIcbit+JyNXAjeRbA17W\nERSRp4DzwCFVXbaetggTEXlQRI5W+Ho7gKreoapbyY+e/VOv6ym85g4gU6jJ83pMexCRGPB14INl\nre6WU9Ws5q/k3wK8SkT2Lfd6vwxaW5b6bCj+SvWIyPuAtwJv1MJBp5f1+MAZYGvJ4y2F50wJEQmT\nD5KvqOo3vK6nSFWnROQHwH5gyQ7rtmiZLMdvQ/FFZD/wEeBtqnrZy1p85CfAtSKyQ0Q6gN8BvuVx\nTb4iIgJ8HhhW1b/2QT0DxTORItIJ3MYK21bbj4AVka+TP1uxMBRfVT3b64nICBABJgtPPaaq1V4e\n0Ix63gHcDQwAU8BTqvpmD+r4TeBvyJ/Z+YKqfrLVNZTV80/AG8hf8j8OfFxVP+9hPa8Ffgj8jPy6\nDPAxVf2OR/XcAHyJ/N8rANywAV4CAAAByElEQVSrqn+x7M+0e5gYY/yh7Q9zjDH+YGFijHHCwsQY\n44SFiTHGCQsTY4wTFibGGCcsTEzdROSLIqKFr7SInBSRz4hIV8lr3ikiD4vIlIjMisjPCtdSrS98\nf1BEvlqYtiErIl/07AOZhliYmEY9SP7y+Z3AfwX+BPgMLNwL6WvAU+QvL9gLfADYAfxx4ecj5O8T\n8yk8vrDNNMYGrZm6FVoR61T1rSXP/T354Hg7+XD4s0rDw0WkT1Wnyp77NnBBVd/XzLpNc1jLxLg2\nR37ui/cAs+SH8r9EeZCY9mdhYpwRkVcBvwc8BFwLjKpq2tuqTKtYmJhG7ReRhIgkgUeBw8Dt5Kf9\nM1eQtpjPxPjaYfI3kk8DZ4stERE5AfyaiHQUpmo0q5y1TEyjLqvqiKo+X3ZI81Xyc/JWnPnO61n7\njXvWMjFNoaqPi8hfAXeKyBbyM4i9QP608L8HRoBPAIjIyws/1gPkCo/nVfV46ys39bJTw6ZulU4N\nV3jNu4D/RH6C5BBwCvgm8DeqOlF4TaWV8HlVvdp1zaZ5LEyMMU5Yn4kxxgkLE2OMExYmxhgnLEyM\nMU5YmBhjnLAwMcY4YWFijHHCwsQY48T/B8mijQl9oGL/AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.scatter(yc[0,:], yc[1,:], alpha=0.5)\n", "for e_, v_ in zip(e, np.eye(2)):\n", " plt.plot([0, e_*v_[0]], [0, e_*v_[1]], 'r-', lw=2)\n", "plt.xlabel('PC1', fontsize=14)\n", "plt.ylabel('PC2', fontsize=14)\n", "plt.axis('square')\n", "plt.axis([-3,3,-3,3])\n", "pass" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Eigenvectors from PCA" ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[ 0.677, 0.736],\n", " [-0.736, 0.677]])" ] }, "execution_count": 35, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pca.components_" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Eigenvalues from PCA" ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([1.404, 0.452])" ] }, "execution_count": 36, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pca.explained_variance_" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Explained variance\n", "\n", "This is just a consequence of the invariance of the trace under change of basis. Since the original diagnonal entries in the covariance matrix are the variances of the featrues, the sum of the eigenvalues must also be the sum of the orignal varainces. In other words, the cumulateive proportion of the top $n$ eigenvaluee is the \"explained variance\" of the first $n$ principal components. " ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([0.906, 0.094])" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "e/e.sum()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Check" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [], "source": [ "from sklearn.decomposition import PCA" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [], "source": [ "pca = PCA()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Note that the PCA from scikit-learn works with feature vectors, not observation vectors" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [], "source": [ "z = pca.fit_transform(x.T)" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([0.906, 0.094])" ] }, "execution_count": 30, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pca.explained_variance_ratio_" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### The principal components are identical to our home-brew version, up to a flip in direction of eigenvectors" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAARMAAAEOCAYAAACq4fP9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3X10XHd95/H3dx40I2v0EMuyLT/GjhPHxk1JEJAcKLCBFLdQKOzh0JZDy3bPettuw8OWw1KyuxzawznQcJp2U/Zw3IVCF2g2LKRwCASchK5JSQJOSIJtESPZiYkty7JsjTSyZjQP3/1jZpTRZCTNw2/m3pG/r3N07BmNRt+5uvdzf/d3f/d3RVUxxphGBbwuwBizOliYGGOcsDAxxjhhYWKMccLCxBjjhIWJMcYJ34SJiERF5Mci8rSIHBORT3hdkzGmeuKXcSYiIkCXqiZEJAw8AnxAVR/zuDRjTBVCXhdQpPlUSxQehgtf/kg6Y8yKfBMmACISBJ4AdgGfVdXHK7zmAHAAoKur6xXXX399a4s05gryxBNPXFDVgWpe65vDnFIi0gfcB9yuqkeXet3Q0JAeOXKkdYUZc4URkSdUdaia1/qmA7aUqk4BPwD2e12LMaY6vgkTERkotEgQkU7gNuDn3lZljKmWn/pMBoEvFfpNAsC9qvptj2syxlTJN2Giqs8AN3pdhzGmPr45zDHGtDcLE2OMExYmxhgnLEyMMU5YmBhjnLAwMcY4YWFijHHCwsQY44SFiTHGCQsTY4wTFibGGCcsTIwxTliYGGOcsDAxxjhhYWKMccLCxBjjhIWJMcYJCxNjjBMWJsYYJyxMjDFOWJgYY5ywMDHGOGFhYoxxwsLEGOOEhYkxxgkLE2OMExYmxhgnfBMmIrJVRH4gIsdF5JiIfMDrmowx1fPNjcuBDPBnqvqkiHQDT4jIIVU97nVhxpiV+aZloqpjqvpk4f8zwDCw2duqjDHV8k2YlBKRq4Ebgce9rcQYUy3fhYmIxICvAx9U1ekK3z8gIkdE5MjExETrCzTGVOSrMBGRMPkg+YqqfqPSa1T1oKoOqerQwMBAaws0xizJN2EiIgJ8HhhW1b/2uh5jTG18EybAa4D3AreKyFOFr9/0uihjTHV8c2pYVR8BxOs6jDH18VPLxBjTxixMjDFOWJgYY5ywMDHGOGFhYoxxwsLEGOOEhYkxxgnfjDMxV67hsTgPHB3nzNQcm/s62b9vA3sGe70uy9TIWibGU8NjcQ4ePkV8Ls1gb5T4XJqDh08xPBb3ujRTIwsT46kHjo7T2xmmtzNMQGTh/w8cHfe6NFMjCxPjqTNTc3RHFx9td0dDnJma86giUy8LE+OpzX2dzCQzi56bSWbY3NfpUUWmXhYmxlP7920gPpcmPpcmp7rw//37NnhdmqmRhYnx1J7BXg68bge9nWHG4kl6O8MceN0OO5vThuzUsPHcnsFeC49VwFomxhgnLEyMMU5YmBhjnLAwMcY4YWFijHHCwsQY44SFiTHGCQsTY4wTFibGGCcsTIwxTliYGGOcsDAxxjhhYWKMccJXYSIiXxCR8yJy1OtajDG18VWYAF8E9ntdhDGmdr4KE1U9DFz0ug5jTO18FSbVEJEDInJERI5MTEx4XY4xpqDtZlpT1YPAQYChoSH1uBxfsZtZGS+1XZi0QjtulMWbWfV2hhfdzMpv86m247I11Wm7w5xma9c7zLXDzazaddma6viqZSIi/wS8AVgnIi8AH1fVz7eyhtKNElj494Gj477eg56ZmmOwN7roOb/dzKrWZWutmPbiq5aJqv6uqg6qalhVt7Q6SKB97zDXDjezqmXZWium/fiqZeIHm/s6ic+lF/aa4L+NspL9+zZw8PApIL+BziQzxOfSvPuVWzyu7EW1LNt6WojWkvGWr1omftCud5hrxs2shsfi3HXoBB/+2tPcdehEw62C4rI9NZHg0dEL3P/MGI+NTnLdhq6XvLbWFqK1ZLxnLZMyxY2ydA/37lduaYs9nMubWTXj7NCewV7etGeAux8eJZ3N0d/VwWBvlAeHJ9g5EFv0vrW2ENu1r2s1sTCpwO4w17yN88T4LDfv7F8UEvG59Evet9bDtnbogF7tLExMRc3aOKt931pbiO3a11Wq3ft8LExMRc3aOGt531paiO3QAb2cdhl0uBzrgDUVNasjulnv24wO6FZqh0GHK7GWiamoWR3RxU7YLz16mvHpJBt6ovzBLducbPTt3Ne1Gvp8LEzMkpqxcQ6PxXlweIK9gz28esdaZpKZimdz/KJV/Riroc/HwsQsaMWG0+pTuI18plb2Y7R7nw9Yn0lbcj2YrPierRj01crLFRr9TK3sx2j3Ph+wlklbKN27RoLC2XiS7f1dTveWrWox1Nucr6eF0ehnanU/Rjv3+YCFie+VN7UPn5ggkcww2BslIGFnG72LDWe5Db74veNjcX45Ocd1G2JsX9dVVXO+uAxyuRxj8SQ/PX2J7x07x+23XsNbbtjctM+0GvoxWsnCxOfK967z2RyxSJCRiVkGuvMbiou9ZaMbznL9C8DC967f2MOacJBnxxNcTmd52abeFc8SPXB0nFwux7PjCSKhAP1dHUwnM9z90OiyHbeNfiYX/RjVBGy9/Tl+G+BmfSY+V97H0BMNo0CiZLoBF3vLesZ/lPbdfPxbx8lmcxX7F8r7Hq5eF+Pmnf28bFMvH7rtuhU3gjNTc4zFk0RCAaLhICJCTzREJqfL9l80Oqal0X6M5fpsGunP8etFjdYy8bnyveuu9V08OnqR7miInKqTXv/iXi6RSnNmao6eaGjFFkN5S+Snpy8xfTlNLBqq2GJa6nCjmj3s5r5Ofnr6Ev1dHQvPpTI51naFl22RuRgr00g/xnJ9NsXH9fTn+PWiRgsTnytvaoeDQXas62JDT4SxeLLhwWSloXD9xp6FcFqp2Vy+Qq+LRYjPpRcdfpW2mCodbkSCUtWp1/37NvC9Y+eYTmboiYZIZXKkMjm2r12zYovMy07Nlfps6u3P8esANwsTn1tqxOhyHY+1qHcvV75C93eFGTk/w5lLc6DKxp4owWBgocVUqe+hMxyo6nfvGezl9luv4e6HRrk4m2ZtV5jta9cQDAZ8Pc/MSn029fbn+LVj2MKkAj91bjV7xGi9e7nSFXpiJsnJC5fp6QwzM5fmF+cTPDue4DU71wJLH258/pHn6I9VN+bkLTdsZudAzDd/l3Ll68x1G7q4MJPkhyOTXLUmzJ7BbqLh0KJD0no7d/06wM3CpIzfrt584Og42WyO42PTJJIZYtEQG7sjzo6P693Lla7QI+cTAKhCdzRMT2cYVDkTTy5adpX6Qmr53X4dh1G+zpyaSPCNJ1/gpm19vGbXWo6fneFHoxd57a7+RetRvf05fp3Ay8KkTDM7t+pp8Rw7G2fkfIK5+QyZnDKZEM5PJ7mczjZUS9Fye7nl6i1docdnUmzojhAOCIGAEA0HUVUSqezCGZ1Kn9Ove9hala8z52ZSdEVCnJtOcfPOftbvfjE0yw/fGpm1zuvwKFf1qWEReY2I3CkiHxORrWXfu0pEHnZfXus1a7h3vafzxuNJJhPzgBAJBQFhMjHPeDzZUD1FS53+BFasd89g/tTub798M3s39ZJRJRLKr1KpTI5YNLTsslsNQ8jhpetMIpmhOxJkOpleeM4PHaTNVlXLRER+C7gPeALoBv6LiPyuqn6n8JIO4PXNKbG1mtW5VW+LZ3Y+Q1AWPxeU/POuVNrL3XXoRNX1FlsYHcEAyXQWESGVyfGyTT2Ll50UPoi+eFdXP+5ha1W+zsSiIaZ90kHayv6/alsmdwB/oaqvVtW9wMeAe0XkHU2pykPNmryn3hZPJBxkoCdCMCjMZ3MEg8JAT4RIONhQPSuppd5iC2Pfph6m5vIh9/KtvXSEgksuu2ZcrOiV8nVmY3eE2VSGjT0RT+9w0OrBbdX2mewFfq/4QFU/KyLngC+LyO8DjzSjOC80q3Or3hbPjVv7ePzkRQZiESKhAKlMjkQyw41X9zVUz3KGx+KcvniZn56+xLpYhF3ru1gXi67YOfrJd96waE/Y2xmuuOxa2cndij1z+TqzYyDGm/dt4MT4rKcdpK0e3FZtmCSBtcDJ4hOq+nUREeAfgY86r8xDzWh619vZ+N5btjMWT3Jxdp6ZZIaOUIBt/Wt47y3bndZXVNzQN3ZHmL6c36Meee4S12/sJhAIrFhvNcuuVSt5K0Or0ud+i9PfULtWD26rNkx+CtwKHCl9UlX/r4gEgC+7Lmy1qbfFs2ewl4/s392UvWulvXbphh6LhhiZmOViYp6x6RSfeNteT8e21Mqvw85bpdWD26oNk8+xRAerqt5bCJT/2GgxIrIf+FsgCPwvVf1Uo+/pJ/W2eJo1fWKlvXYileb6jT0ADHRHGeiOklNlLJ50VkNxJZ/PZBmZmCWRzBAOCi/b1OPk/Yv8Ouy8VVp96r2qMFHV+8ifzVnq+/cA9zRSiIgEgc8CtwEvAD8RkW+p6vFG3nc1ctEPsNRe+8zUHDPJTFP3Zvv3beCvHniW05OXiUWChAL53zE+nWJ4LO48tFx+Fj+Njl5Jqwe3VXtqeAD4I+BvVXW67Hu9wPuB/6mqkw3U8ipgRFVPFt73HuDtgIVJCVf9AKV77YmZJCMTs8zMpZlLZ+iOhGDtmqbM4QGwZ1Mf/1D9R67bh5rwnnsKX83y4XufchpSrTz1Xu1hzgeAq8uDBEBV4yJyLfBB4L81UMtm4Jclj18AXl3+IhE5ABwA2LZtWwO/rj3V0g+w3EZdeqjx5OkpIqEA4aDQEeogp0o6k2Usnql5b1YMu2w2x7npJE+dnuJ7R89x+xuv8bxDsh344RKOelUbJr9FPlCW8gXgf9BYmFRFVQ8CBwGGhoZ0hZf7gsumcbX9ACu1YIrH0ycnEkQKo+Lms8ortvcSDgbp7Qzzoduuq7m+4rVEJ87nZ0Vb2xXOz4r28Cg7z04tfO67Dp0gPpdmeGyaZDpLNBwkmc4SCQfZO9jD8bFp9g72vOSexL2dYX4xPs2PRi8iAqGAkM0p85kcV6/r4o/fsJMHhyfo7Qwvalm52DA//LWnC9NlvjiKsNif9Jl3/Wrd71tcFr2dYQK0pqO4GYdr1Q5auwYYXeb7J4EdDVUCZ4DSYfpbCs+1lOvBVPc/c4YP/Z+n+fYzZzk9OctzFxINDRza3NfJTHLx6NdK/QArzaxePJ5OZ5VUNkc0HOQV2/tYF4tW1Um51HI6MzXHuemXzoqWzuYWzYpWHOh1IZGiIygk01lSmRy7BrrojoYYn05WHDR37Gycfx2dJFAIkulkhtn5LOFg/nff/dAouVzlGd8aVe2yr1UrZ+yHwjp5z9Pc/8wYpydnOTXR2DpZVG2YpFm8oZfbAjQ6vvsnwLUiskNEOoDfAb7V4HvWxPWIweGxOHc/lM/g/q4OUpn8PKa5XK7ulbvaEbqlK+iFRJLHTk7y2MkLHDo+vvB59gz2ctveDdy8cx037+xnXezFSY0iQVkyVJdbTpv7Ork4m164Rgfy1+n0d3Us2jiKYdYfi3DpcoZIOMhN2/oY6M4PjtvQE6244U4nMwhCqBBAoYAQCkAqqwQDQiaXbymUcrVhNmN0dHGA4Hd+NsZjJye5kMjX3qxTuMNjce5+eBQE1naFSWXyrchstv51sqjaMHkSWG7o/L8lPxalbqqaAf4U+B4wDNyrqscaec9aub5PygNHx8nklJ5oCJH81bS5nPLj5y7yz0+dqavlU+3FccW96IVEkieen8ofQgTz/SKlAVlpA/nlxcucjScXwuLURIIP3fM0/+Eff8Jdh07w5UefX3I57d+3YaHFoKoLLY7B3uhLNo49g7184m17+ZUtvewd7KG/MFtbfC7NH9yyreKG29sZYnNflPmsMp/NAYoCmayyoSdKNCz84nyC7x8/t7BxutowXV+YWDpAMBwILAwQfO5ComnD7x84Ok46m1u0TkZCAc5NJxsO3Gr7TD5L/lqcF4C/U9UsgIiEyAfA+8m3JBpSuHDwOyu+sEmW6o84djbfpK/1+PLM1NxC+kfDQWZTGSZn50lnc2y9ak3dHW3V9NCX9ol0FPpEUlnlpm35a2aKx+OVTh/O90ToCAUXJj46cT4BAvHCiNgfjkzyml1rgRf7M4p7/z2Dvdz+xmu4++FRJmfn6e/q4Or+NQQClWdFW+70ZflkSO9+5RYeODrOmsKyTGeyZHPFOFEuzaaYnssQEOgICHPzGR4dvciOdV28+8219//Uu+yr1YoBguXOTM0ttJKjheu7IqEAF2fT3HJNY4Fb7TiTb4jIp4G7gL8UkWL/yTVAF3Cnqn69oUp8oNK4hOcvzPLCpTm2XLWm5p72zX2dzKez+Y0RmJxNkc0p4WCAXetjTe1oK26k//nep0GhuzPEyzb1LAxCKz/kKP39H/7a0wszoI1MzBIJBYiEAsyk8uNPrloT5vjZGdbvfnHlK9371zor2lIb6FLPHzx8mV/d2sfo+QQnJ2dJZ2BjT4RUVgkEhK6OADkgm8uH3IaeiC/PipTuvJo1QLDc5r5O0pksz47n18lIKMB0MkMoIA23hKqeHElV7xCRfwbeA1wLCPD/gK+q6o8bqsInKo0YPDGeYPfGFzf8WgIg/36XuW59jHPTSRLJLB0h4dU7rnJ6z5ul7Bns5df3bqx54FZpqCaSGWKRIKlMjp5ouPC+3fxo9CLxufSSY1GaNb6htCUTDgVBhFQmSzgYZG4mydarOgkF852/N+/sX9g4/ciLuVyL6+TuDTHG4kkmZ+cJBwPc/sZrGv57VTtobQ1wJ/Db5Nu2DwG3q+qFhn67z1Rqcm/t72Rbf9ei11UbAKXv1xHOr/hrwgEmZ9OcOn6OnmiYjT0Rrl4Xa9ZHqmtIdenPxCJBpgsdofs254e7R8MhXrurn97OsCdXxZYGVenp2kdPTpJKZwt72/zERH6YaHkpXsw0Vx7Gt1yzztkAOVFdeaiGiNwJ/AnwFWCO/HQE/6Kq72q4ggYMDQ3pkSNHVn5hA+46dILnLiQ4N51iOpleFAC1jsO4/5kzfOq7z9IVCdEdCXIhMc9kYp69m7p59U53f9Ry9YwpKP7MsbNxXrg0x+4NMbb1d604bqPVw81Lx2hMzCT50cgFZlJZRPJ7/rVdHXxk/272DPbWXVszP5Pfh+eLyBOqOlTVa6sMk1HgjsI1OIjIq4B/BaLFzlgvtCJMygNgJpVlNpXho7+xu+bbTdx16ASnJhKcm0lxYSZJfC4/vd/6nih7N/U6G1zlWrUrfOlAOdcDxparrfg7k+kMj/xikkQqw1VrQvR0dhAJ5funzs+kagpFLz+Tn9QSJtX2mWwFflh8oKo/FpEMsInFQ+A95zrpT4zPcuPWPs7NpEgkM/R0hrlufYwT47M1Dw8/MzXH9nVd7BiI8ejJSbqj2UUdm9Del8f/70ef5+REgnRWiUVD7BroWnZCaRdKm+2Pn5qkP9bB63evY10syoVEkkdHL5LKTC+Me3l2PEEsGloYU7NSbVf6NAa1qDZMgsB82XOZGn6+JZoxGU5pABSVnw2p1kodm368PL7aZTo8FueRkUmu6gzlP1c6f83Py7f2kki5m6+2kmIfSvHsSHG4+8j5WWKRIOmsks5mF+4GOHJ+tuqRvlf6NAa1qDYMhPwUjamS56LA34vI5eITqvo2l8XVqhl7EZc97it1bPqxs7CaZTo8Fufj3zrOTDJNKp1loDtCVyRU+N4Mb9i9viW1lv+tppNpOgJCrDASuJ7OWb/ePc+Pqh0B+yXgLDBZ8vVl8oc4pc95qhnXOLgcQl06grKnsHLu3hBjbVekpvdt5WTMKy3TYsvlYmKeTb0RUpkcZy7NkUilUVUuXW7dRMrXbejisdFJ7n9mjEdHL5DJ5kiksuwa6GLXQBepTI7pZH6KhWqXd7MmGF+Nqh209u+aXYgLzdiLuJ5gpvS0Znn/zlKTL5ffdrJ4VexKhx0u+o5WWqbFlsvaWAepdJatazs5P51ifDrF5r5Ofm1Xf0v6Foq3Ud29MT9+Yiw+x0wyS0CUY2en2TPYXejrStC7JrzkZNflmj3BkN/P5tSiqrM5flV+Nme19bxX+jyPjU6ye2Ns0diU4sZePFXtcjms9F7FcR6TidTCvCgdQeHi5TQ3bOlr2bIvPUVcvB4JIBSAYCDApctpXrurn/fest0360I7rK+1nM2p+o5+7cD1hVheq3ThYTVXxbq8YHGlZVq8oHCgO8pN2/qIhINcupyhPxapatm7OmQrPRwbOZ+/BKAnGiKnwht2r+fW69cz0B311brg+sJSr/nqbIwLrZymbiWNNGGHx+J8//i5wnU1YXYNdDHQHWVtV5jJ2cUn1soP5VyfgVhumZZ2KvfHIgs33qo2SFydfSs9HJtOpumOhBZuUQr+PAOz2s4UraqWiZ80MjdK8Wc7ClMGFE+zTswk2dgTJRwMLNsh2KxJfCpppDXocs9c2lHaHQkxXThkmE1lOHR8nMMnJhaunvaLVv6dWmHVtUz8opHT1MWf3be5hyeenyISEiJB4djZaXYOxLj91muWvVtcq6/5qLc16HLPXNpR2rsmzPmZJJmMEoqCSnNmv2+UF9fmlHLd+Wth0iSNbCgvDr4K84rtfYycnyU+N09AAgt7/eVG3zb7DIQrrs++lYbaHd94hqNnp5nP5gcF/sqW/Ny2fhq56uXfqRkDPC1MmqSRDaX0Z9fFoqyLRRceV/uH9lPfUSXDY3EuzCT54cgkV60Js2ewm2g45GzPnMoqr7tu4CWTP/utP8Krv1MzBnhan0mTNDLYabUPlCruFcOhYGHGNvjR6EXmM1lnZ99WW3+Ea80Y4GktkyZppAnbLocpUN9x9+K9Ypj1uztrbnmtxOv+CL9rxgBPC5MmaqQJ6/fDFKj/uLsVp0TbKZC90IywtTAxdav3uLtVF8+1QyB7pRlha2Fi6lZvC8MOQfzBddhaB6ypW72dnKvtsgeTZy0TD7X7FaONtDDsEGT1sZaJR1zfitQL1sIwpaxl4pHVMrdorS2Mdm+NmaVZmHhktV0xWo1mDOFuV6sxVH1xmCMi7xKRYyKSE5GqJmJpd1fiCM3VNn9HvVbDIW4lvggT4CjwTuCw14W0ymofMl9JM4Zwt6PVGqq+CBNVHVbVZ72uo5WuxM7LK7E1VslqDdW26zMRkQPAAYBt27Z5XE1jrrTTozZYLW+13j6jZS0TEXlQRI5W+Hp7Le+jqgdVdUhVhwYGBppVrmmCK7E1VslqPcRtWctEVd/Uqt9l/OtKa41VslovQmy7wxxjVoPVGKq+6IAVkXeIyAvALcD9IvI9r2syxtTGFy0TVb0PuM/rOowx9fNFy8QY0/4sTIwxTliYGGOcsDAxxjhhYWKMccLCxBjjhIWJMcYJCxNjjBMWJsYYJyxMjDFOWJgYY5ywMDHGOGFhYoxxwsLEGOOEhYkxxgkLE2OMExYmxhgnLEyMMU5YmBhjnLAwMcY4YWFijHHCwsQY44SFiTHGCQsTY4wTFibGGCcsTIwxTliYGGOc8EWYiMidIvJzEXlGRO4TkT6vazLG1MYXYQIcAvap6g3ACeDPPa7HGFMjX4SJqn5fVTOFh48BW7ysxxhTO1+ESZk/BL7rdRHGmNqEWvWLRORBYGOFb92hqt8svOYOIAN8ZZn3OQAcANi2bVsTKjXG1KNlYaKqb1ru+yLyPuCtwBtVVZd5n4PAQYChoaElX2eMaa2WhclyRGQ/8BHg9ap62et6jDG180ufyd8B3cAhEXlKRD7ndUHGmNr4omWiqru8rsEY0xi/tEyMMW3OwsQY44SFiTHGCQsTY4wTFibGGCcsTIwxTliYGGOcsDAxxjhhYWKMccLCxBjjhIWJMcYJCxNjjBMWJsYYJyxMjDFOWJgYY5ywMDHGOGFhYoxxwsLEGOOEhYkxxgkLE2OMExYmxhgnLEyMMU5YmBhjnLAwMcY4YWFijHHCwsQY44SFiTHGCQsTY4wTvggTEflLEXlGRJ4Ske+LyCavazLG1MYXYQLcqao3qOrLgW8D/93rgowxtfFFmKjqdMnDLkC9qsUYU5+Q1wUUicgngd8H4sC/WeZ1B4ADhYcpETnagvKqtQ644HURJfxWD/ivJqtneburfaGotqYRICIPAhsrfOsOVf1myev+HIiq6sereM8jqjrksMyGWD0r81tNVs/yaqmnZS0TVX1TlS/9CvAdYMUwMcb4hy/6TETk2pKHbwd+7lUtxpj6+KXP5FMishvIAc8Df1Tlzx1sXkl1sXpW5rearJ7lVV1Py/pMjDGrmy8Oc4wx7c/CxBjjRNuHid+G4ovInSLy80JN94lIn8f1vEtEjolITkQ8O+UoIvtF5FkRGRGRj3pVR0k9XxCR834ZpyQiW0XkByJyvPD3+oDH9URF5Mci8nShnk+s+EOq2tZfQE/J/98PfM7jen4dCBX+/2ng0x7Xs4f8wKN/AYY8qiEIjAI7gQ7gaWCvx8vldcBNwFEv6yipZxC4qfD/buCEl8sIECBW+H8YeBy4ebmfafuWifpsKL6qfl9VM4WHjwFbPK5nWFWf9bIG4FXAiKqeVNV54B7yQwA8o6qHgYte1lBKVcdU9cnC/2eAYWCzh/WoqiYKD8OFr2W3rbYPE8gPxReRXwLvwV8XCf4h8F2vi/CBzcAvSx6/gIcbit+JyNXAjeRbA17WERSRp4DzwCFVXbaetggTEXlQRI5W+Ho7gKreoapbyY+e/VOv6ym85g4gU6jJ83pMexCRGPB14INlre6WU9Ws5q/k3wK8SkT2Lfd6vwxaW5b6bCj+SvWIyPuAtwJv1MJBp5f1+MAZYGvJ4y2F50wJEQmTD5KvqOo3vK6nSFWnROQHwH5gyQ7rtmiZLMdvQ/FFZD/wEeBtqnrZy1p85CfAtSKyQ0Q6gN8BvuVxTb4iIgJ8HhhW1b/2QT0DxTORItIJ3MYK21bbj4AVka+TP1uxMBRfVT3b64nICBABJgtPPaaq1V4e0Ix63gHcDQwAU8BTqvpmD+r4TeBvyJ/Z+YKqfrLVNZTV80/AG8hf8j8OfFxVP+9hPa8Ffgj8jPy6DPAxVf2OR/XcAHyJ/N8rANywAV4CAAAByElEQVSrqn+x7M+0e5gYY/yh7Q9zjDH+YGFijHHCwsQY44SFiTHGCQsTY4wTFibGGCcsTEzdROSLIqKFr7SInBSRz4hIV8lr3ikiD4vIlIjMisjPCtdSrS98f1BEvlqYtiErIl/07AOZhliYmEY9SP7y+Z3AfwX+BPgMLNwL6WvAU+QvL9gLfADYAfxx4ecj5O8T8yk8vrDNNMYGrZm6FVoR61T1rSXP/T354Hg7+XD4s0rDw0WkT1Wnyp77NnBBVd/XzLpNc1jLxLg2R37ui/cAs+SH8r9EeZCY9mdhYpwRkVcBvwc8BFwLjKpq2tuqTKtYmJhG7ReRhIgkgUeBw8Dt5Kf9M1eQtpjPxPjaYfI3kk8DZ4stERE5AfyaiHQUpmo0q5y1TEyjLqvqiKo+X3ZI81Xyc/JWnPnO61n7jXvWMjFNoaqPi8hfAXeKyBbyM4i9QP608L8HRoBPAIjIyws/1gPkCo/nVfV46ys39bJTw6ZulU4NV3jNu4D/RH6C5BBwCvgm8DeqOlF4TaWV8HlVvdp1zaZ5LEyMMU5Yn4kxxgkLE2OMExYmxhgnLEyMMU5YmBhjnLAwMcY4YWFijHHCwsQY48T/B8mijQl9oGL/AAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.scatter(z[:, 0], z[:, 1], alpha=0.5)\n", "for e_, v_ in zip(e, np.eye(2)):\n", " plt.plot([0, e_*v_[0]], [0, e_*v_[1]], 'r-', lw=2)\n", "plt.xlabel('PC1', fontsize=14)\n", "plt.ylabel('PC2', fontsize=14)\n", "plt.axis('square')\n", "plt.axis([-3,3,-3,3])\n", "pass" ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAADdCAYAAADuKuYGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJztvXl8XGd56P99ZtGMpNFiWbIsO94Sx44TkxAwpClcSFnThkJpC3RjaWnzo/eWAhdKF+7vtrSXXsreQre09EILlLIE6GUxJNCQBkI2Z8GxY8exE2+yLFmWRiPNPs/945wzHskj6YxG0pyZPN/Pxx9r5pw585wz7/M+7/O8z/u8oqoYhmEYRtAINVoAwzAMw6iGGSjDMAwjkJiBMgzDMAKJGSjDMAwjkJiBMgzDMAKJGSjDMAwjkJiBAkTkURG5wf37T0TkM/Oc90ci8o81XltFZPsyiNl0iMjLReSrFa+fJyKPi0hKRH5ORO4Qkd90j/2qiHynATJudX+jiPv6yyLy06sthxEcFuoDjNWlpQyUiDxfRH4oIpMiMi4iPxCR5yz2OVW9SlXv8HHen6uq16HO6tiCjojcICIl1zhMicghEfn1iuNtrmI+LiLTIvKkiPyTiGydc51PiUhBRIZ8fO37gPdXvP5T4BOqmlDVr1aeqKqfVdWXLf0Ol42/AP5Xo4V4uuG2t5fMee9NInJXo2RqNO4ALuPq7JiI3FqpdyLyXBH5pohMuP3dvZU67Z6zzdX7v139O6ifljFQItINfB34ONAHbATeC2QbKddKIA5L+e1Oq2oC6AZ+H/gHEbnSPfYl4JXArwA9wDXAA8CLK763E/gFYBL4tUVkfA7Qo6o/qnh7C/DoEuReNVT1XqBbRPY0WhYj+KzCAPV3XJ3dAfQCH3W/93rge8D3ge3AWuC3gbne/xuA88DrRCS2wrIuOy1joHB+QFT1X1W1qKppVf2Oqj7inSAivyUiB10P4oCIPMt9/6LRm/t+VET+1Q37eB6G5/rf6f4/4Y5wrl9MQBGJiciHROS4iIyIyN+JSLt7bI2IfF1ERkXkvPv3JRWfvUNE3iciPwBmgEvd9/7M9RSnROQ7ItK/mBzq8FWchnule+8vBV6lqvepakFVJ1X1r1X1kxUf/QVgAscTeuMiX/PTOMrjyf8EcCnwf93nNUtZ5o6WXe/0d0XkqDt6/KBnlN1zfyAin3C95cdEpNKQ9ojIJ0VkWEROicj/EpGweyzs/gZjInIUuKmK7HfM877RQERkl9vmJ8QJy7+y4tinRORvRORbbvv6gYisF5GPufr0mIhcW3H+BlevR0XkmIj87pyvi4vIv7l6tU9Erqn47JMi8vsi8ggwLSKR+a4nInERSXt6KSLvEScC0e2+/jMR+dhi966q48CXgd3uWx8EPq2qf6GqY65OP6Cqr62QU3AM1P8A8sDP1vK8g0ArGajDQFFEPi0iPy0iayoPishrgD/B+cG6cbyFc/NdzDUcX8XxwF6rqrk5p7zA/b/XDVnd7UPG9+MY0mfijHo2Av/TPRYC/g+Ol7EZSAOfmPP51wM3A13AU+57vwL8OrAOaAPetZgQIhISkVfjjMh+DLwEuFdVTyzy0TcC/wp8HrhCRJ69wLnPAA55L1T1MuA48LPu8/Lj2b4a2AM8C3gV8BsVx64DngD6gT8GbhWRPvfYp4ACzjO+FngZ8Jvusd8CXuG+vwf4xSrfexDHgzQCgohEgf8LfAenrb8V+KyI7Kw47bU4nXE/jt7eDexzX38J+Ih7rZB7rYdxdPDFwNtF5OUV13oV8EWcaMzngK+6Mnj8Ms4gphcozXc9Vc0A9wEvdD/3QhzdfV7F6++zCK6B+wXgQRHpAK5372khng9cgqOvX2DxQWXgaBkDpapJnB9EgX8ARkXk30Vk0D3lN4EPuB6CquoRVX1qnst1A3txOsBfV9VivfK5o5mbgXeo6riqTgF/DvySK/85Vf2yqs64x97HhUbt8SlVfdT1cPLue/9HVQ+rahqnET5zATE2iMgEMIbTqb9eVQ/hhAeGF5F/M/BTwOdUdQT4Lo6xn49eYGqha/rgL9xndRz4GE6n4HEW+Jiq5lX133CM4U3u7/0zwNtVdVpVz+KERX7J/dxr3c+dcEel/7vK90658hury1dd72jCbad/U3HsJ4AE8H5Vzanq93BC+pVt4iuuF5EBvgJkVPWfXf39N5xBCcBzgAFV/VP3Wkdx+oxfqrjWA6r6JVfPPgLEXRk8/sptQ2kf1/s+8EJxwoFXA3/lvo67n72T+fkr91k8jKOj/x1Yg9N3L6izOAbpW6p6HsfI3igi6xb5TKBoigl+v6jqQeBNACJyBfAZLnRsm3AMjh9+AogCv6w1VNMVkUdxPCCAn1bV/6w4PAB0AA84tsr5COCFnjpwOtIbcRogQJeIhCsMZDUP50zF3zM4Sjwfp1X1kirvn8MNkS7A64GDqvqQ+/qzwIdF5F0VxrKS8zieXj1U3u9TwIaK16fm/Dbe8S04v91wxXMOVVxrQ5XrzqULJ5RprC4/p6q3ey9E5E1c8Hw3ACdUtVRx/lM4HovHSMXf6SqvPd3YwoXBmkcYqNTXchtR1ZKInGR2+6tsQ4td7/s4Ru5ZOBGL24BP4vQzR1R13kgO8LuqOitz2O0rSsAQ8Fi1D7kRoNfgPj9VvVtEjuNEXBYNKQaFlvGg5qKqj+GEeryY7QngMp8f/w7OyPq7FR7YRV9R5TuvcsNXiTnGCRyvJQ1cpaq97r8edwIU4J3ATuA6Ve3mQghRKq6xUqXnbweeWznnVYU34Mx7nRGRMzgK14/jrVTjERY3eouxqeLvzcDpitcbpcICVRw/gRPe6a94zt2qepV73nCV685lF86I1QgOp4FNMjs5aDNwagnXOgEcq2gfvarapaqVbbncRtzvvITZ7a9SFxe73g9xdPvVwPdV9YAr+8/gI7w3F1WdwQlf/sICp70aJxL0NxU6u5EmC/O1jIESkStE5J1eJysim3A8Jy+L7B+Bd4nIs8Vhu4hsme96qvoBHLf4u1I98WAUZxRzqR/53JHfPwAf9dxsEdlYEffuwjFgE+5cyh/7ue5y4I5abwO+4j6fiIh0ichbROQ3xEkAuQx4Lk4I8Zk4hv9zzB/m+yYXhyhr5ffESR7ZBLwNJ0zjsQ74XXESWV6DY1S+qarDOAOMD4tItzvfdpmIeLJ8wf3cJe485R9U+d4XAt+qU3ZjebkHJ0Lwbvc3vwFn0v/zS7jWvcCUm+jQLk7izG6ZvSTl2SLy825Y7u04g54fVb3aItdzDcoDwH/jgkH6IfAWlmCgXN4NvElEfk9E1gKIyDUi4j2PNwL/hDMX7Ons84BrROQZS/zOVadlDBTOvMF1wD0iMo3TmPbjeCao6hdx5nU+5577VZwJ0HlR1T9zz7u9YgLeOzbjXu8Hbsz8J6pdYw6/DxwBfiQiSRzPxZvk/RjQjuNp/QhnDmw1+UUco/JvOGnk+3GSCG7HaexfU9Ufq+oZ7x/wl8Ar5j4bAFXdB0yKyHV1yPQ1HMV+CPgGTljE4x7gcpzn9T7gFytCJW/ASRg5gBNq/BJOOAScQcK3cTykfcCtlV/odiopddLNjYDgJin9LE526BjO/NQb3EhJrdcq4iTKPBM45l7vH3GWV3h8DXgdTvt5PfDz84Sy/V7v+zih53srXnex8PzTQvfwQ+BF7r+jIjIO3AJ8U0S8RI2PVeqrqj6A0680jRclNUyxGEZNiMjLgP+qqj+3hM8qcLmqHqly7E3Ab6rq8+uX8qJrfxn4pKp+c7mvbRhGbbRUkoQRLFT1OzjhtqZBVReK6xuGsYoEJsQnzoK2e0XkYXEW4b230TIZRjNiumS0CoEJ8bkZWZ2qmnIXxN0FvE1nl8oxDGMRTJeMViEwIT53TUvKfRl1/wXDehpGE2G6ZLQKgTFQ4NRJw8na2g78tareU+Wcm3EqMtDZ2fnsK664YnWFNJ7WPPDAA2OqOtBoORbDdMkIMn71KDAhvkpEpBenVMlbVXX/fOft2bNH77///tUTzHjaIyIPqGrTVDo3XTKCiF89CkySRCWqOgH8B07ZH8MwlojpktHMBMZAiciAO9rz6ki9lHnqTBmGMT+mS0arEKQ5qCHg027sPAR8QVW/3mCZDKMZMV0yWoLAGCh1Nha8dtETDcNYENMlo1UITIjPMAzDMCoxA2UYhmEEEjNQhmEYRiAxA2UYhmEEEjNQhmEYRiAxA2UYhmEEEjNQhmEYRiAxA2UYhmEEEjNQhmEYRiAxA2UYhmEEEjNQhmEYRiAxA2UYhmEEEjNQhmEYRiAxA2UYhmEEEjNQhmEYRiAxA2UYhmEEEjNQhmEYRiAJjIESkU0i8h8ickBEHhWRtzVaJsNoRkyXjFYhMFu+AwXgnaq6T0S6gAdE5DZVPdBowQyjyTBdMlqCwHhQqjqsqvvcv6eAg8DGxkplGM2H6ZLRKgTGQFUiIluBa4F7qhy7WUTuF5H7R0dHV1s0w2gqTJeMZiZwBkpEEsCXgberanLucVW9RVX3qOqegYGB1RfQMJoE0yWj2QmUgRKRKI5CfVZVb220PIbRrJguGa1AYAyUiAjwSeCgqn6k0fIYRrNiumS0CoExUMDzgNcDLxKRh9x/P9NooQyjCTFdMlqCwKSZq+pdgDRaDsNodkyXjFYhSB6UYRiGYZQxA2UYhmEEEjNQhmEYRiAJzByUYRhGM3NweJK9+0c4NZFmY287N+4eZNdQT6PFamrMgzIMw6iTg8OT3HLnMSbTeYZ64kym89xy5zEODk82WrSmxgyUYRhGnezdP0JPe5Se9ighkfLfe/ePNFq0psYMlGEYRp2cmkjTFZ89Y9IVj3BqIt0giVoDM1CGYRh1srG3nalMYdZ7U5kCG3vbGyRRa2AGyjAMo05u3D3IZDrPZDpPSbX89427BxstWlNjBsowDKNOdg31cPMLttHTHmV4MkNPe5SbX7DNsvjqxNLMDcMwloFdQz1mkJYZ86AMwzCMQGIGyjAMwwgkZqAMwzCMQGIGyjAMwwgkZqAMwzCMQGIGyjAMwwgkgTJQIvJPInJWRPY3WhbDaFZMj4xWIWjroD4FfAL45wbLcRFWSt9oIj5FQPXIMGohUAZKVe8Uka2NlmMuXin9nvborFL6zbRS3Azs04eg6pFh1EqgDJQfRORm4GaAzZs3r8p3VpbSB8r/790/0hSdfLMZWDOmq0MjdGmlsDbTmgRqDsoPqnqLqu5R1T0DAwOr8p3NXkq/mfaqsY3fVo/l0qWDw5N89LbDvOuLD/PR2w6v+m/VyDbT6HtvdZrOQDWCZi+l30wGtpmMqRGMAUWj2kwQ7r3VMQPlg2Yvpd9MBraZjKkRjAFFo9qMn3s3D6s+AmWgRORfgbuBnSJyUkTe3GiZoPlL6fsxsEFRpGYypkFlNfUoCAOKRrWZxe7dj4cVFL0LKoFKklDVX260DPPRzKX0PQNbOYn8uudcUr6fICVR3Lh7kFvuPAY4yj6VKTCZzvO651yyqnI0M6upRxt725lM58uJQ7D6A4p62kw9yRUbe9s5NprizFSWVKZAIh5hfVeMbQMJYPHkqiDpXVAJlAdlLD/eCO2Tdz0JwJufv5V3vHTHLAUIQpjGo9m91acbQQh/L7XN1DuHtGOwkwdPTJBM5+lsC5FM53nwxAQ7BjuBxT2sIOldUAmUB1UPrZZmuhz343eEdmoizVBPfNZnGznv08ze6tONxbzz1ZSj1u+sd/nI4ZFpnrW5lzPJLMmM40XuHExweGSam1jcuwyS3gW1/2wJA9VqrvJy3Y9fBQxCmMZoXpp1QFGvgTg1kWbz2k629ifK75VUy59fLPQYFL0Lcv/ZEgaq2RfSzmW57ufA8CSTM3mmsgW641G2r+ukrzN2kQI2ct6nlpFbUEd5RnNSr4FYbA5qMe8yKPOtQe4/W8JArbSrvNod43Lcz8HhSU6cS4NAdzxCJl/kgacm2DmYmDXig8aFaWoZuQV5lGc0J/UaiB2Dndy67ySdsQhdsTDJdJ7hiTQvr5h/W8i7rKZ3z9nay979I3zyridXbcAWpFDjXFrCQK2kq9yIjtG7n3yxyJGz04wk06RzJTpjYT5622FfDXHv/hF2DCY4fDZFtlAiFgmRLZQ4NJLiLTdcdtH5yxmmmU9x5r7/+EiSUxMZcsVS2cPzJonnyhLkUZ7RXFS2w45oiFyhyPBkoeaB2WJzUH6o1LtGDdiCEmqsRksYqJV0lRvRMd64e5APffswx8amCYsyMVOgVFLi0RBPjqW45c6ZRRviqYk0W/o7ScQjHBmdJpUp0B2P0N0erdrYl8tDnE9xXrJrgNsPjpbff3IsxR2Hx9jYE2dNZ1vZw7t2cw+nJgoXXTfIozxjdahsp7GwoECuqDW12bnt0+srltKxLzYHVfmdn7n7KR48MYmiXLupl9dfv6WuQdhy9ktBCTVWoyUM1EqGqBrRMe4a6mGwO8ZYKsupiTRtEWFdVzvhkHAmmWXXUPeiDdEbFQ10xRnocuSfO0oC+MYjp/j4d5+gUFL6OqPk8kVfBnA+5lOcT999nCuHusuvzySztEfDTGby9CVixKNhAA6cnuKGnevmvZ/VGOXZXFfwqDQs0TDcfXQcAZ6zbU1N3sNSOvZq7QHg+PgMDx2foC/RxvaBTga64he1yYPDk+XBZiIWRhDuOTrO8GSGd9+4c8nZtMvZLwUlE7MaLWGgYOUyiRrh/h4cnuShE5OUtATAuq4YnbEIqkoyk/fVEP2Mig4OT/Lx7z0BAn2dUbKFEofPptixLlFei1FrRz2f4owkM1y3ra/8XjKTZ7CrjVOTWTL5IrFICFQ5ny5UXUOzWqM8m+sKJpWG5UdHk+X1RUfHZrj+0rXlc5baPqt5PXv3j3DAncvdMZhgS38nk+k8H9h7iJAIQ90xJmdyJNN59j01wY7BBOFwaFab3Lt/hLFUlq54pDwIQ4Tx6Vxd2bRL6ZfmGtodg50cHpkOdIJSyxiolWI5Oka/P64XCvjPI+fI5ot0x8OUSsrjZ1O0R8PEImEGumK+DKSfUdFn7n6KkWSGkMBUJkxfZxuxSIgzyQzpQpHj4zM1d9TzKc5gtzO69N7vjkeZTOe5ZE07sWiYVKZANCw8f/vaqvNVN+4e9D3Kq0eZbK4rmFQalmQmT1fM6bpSbokjv95DW1i48/Ao+aKSiEfYPtBJWyR8kdfjDVImZ/IgcPhsikQ8wkBXnPHpHAC7N65zQuhnpxlLZTk2Ns0zLumZleBwaiJNtlCkO35BH2KREFOZQl3ZtLX2S3MHXk+Opbh130mu3dRbNryrmaAk0ZivEb4ZqEWo1/31++N65x0dTdHbHiETCXFqYoZSSVGFTL6IKkzM5DgxPsPrXr7Dl+wLjYj+88g5oiFBBApF5fREhqGeGOPTRRBhY29HzR31fIrzxus3c/vB0fL767tjnJ5IlxXEO+/112+5KJxzx6GzfOXBUzx/+9qqsfulPO/5sLmuYLKxt50nx1KcSWYZncoynsrR2xGlp6MN8BfVODg8yUgyy1SmQCIWJpsrcM/RcTav7eDdN+4sn1c5SDk7lSFXUNL5InccHuWGHQPkCiUUBaA/Eac/EefsVJofHBmnLRJmbSJSbnft0RCxSJhsoVT2oLKFEm2R0EXy1tLX1NovzR14nUlm6YxFODOVZdtAYtXmuzxC8a5eP+eZgfJBPeFDvz+ud16uWKIrFqG9LcJYKkumUKQ9FCJXLLE20UZXPMpgd6zu0fze/SOs6YiSjoQ4N50jEoJwCM5MZulPxOiOR5ZUBHQhxbl0IFF+f2t/gpddNTgrxOCd99HbDtPTHiVfLPLg8UlikRBr2iM8ejq5qLHZu3+EYrHEgeHkrLUpfpUpyBlNQWU1wj+VKd3ruto4OZ5hOldgW39HubzSYlGNvftH2NTXwfqeGEfOTpPM5EnEI2zoiVedCxqdyjCVKSJAWOBcKss3fnyGsEB/Ijbr2gdOT7Gm42I9zxeK9CdiHBubRlURIJUtsnltR9VQdi19zWIDUO83aQsLdx89RzwSpqs9yvaBTtcLDZc9UFid+S4PCUdii59Vg4ESkecBPwecB/5FVU9UHFsDfFlVX1SroK2O3x/XO687HiWTLxKPhhGBjrYIG3vbiUfDbF/XyeMjKX7wxDnf6eYLybVrqIuHTkyytrONqWyeTE4porz1xZdxeGR6yR31fIpT7f1q6bjes7j3WJJiqcRYKk+24MzHbenrWNDYPHp6kpPjaWLRkDNKzhc5PJJiJl9cVG5Y+bmuVtOj1ZqzOzwyzbpElMNnp8kUSkRCQk97hOPjGXau7/EV1fDaVUii9CccnRxJptl3YoJ3ffHhsnH1BilHRqfpT7RxZjJDOl8kEhYiIcjmSyQzeZ4cS7F5reP9n5/J85OX9c36vq54hOHJAu96+Y5ZWXzXXdq3aCSgHip/k0gI7jk6zvnpPOu6hGy+yL7jE0RCwlS2SPcKzXcthhYLWT/n+TJQIvKzwFeAB4Au4PdF5JdV9ZvuKW3AC5ciaDNQb8VjPz+ud972dZ088NQEAOGQkCuUOHk+TTwa4sjZFN3xCINdsbo7Au/7nrW5lyOj00TCIaJdwlUburnp6o1c6mYfPZjKkkznmM4VCYnwoivWcXB4ctnnfqrJNpJMM5UpEgm5VY2FRY1NMlMAoRxOiUed8Eoyc3HqejVWMqOpFfWo3vCP3zZzz9ExHj87TSwSJhGLkC8qM9kiW9eGecdLFw93w8W6ODqV4b5j50nEI1WXRIyncqzpiDCWEqQgtIVDxCIhetqj7N7Qw3AyS9Sdv3r+9rW0RcKzvs/T811DPbzv56/2JeNyUPmbHBhOkohHCIeEc9M5LmmLEAsLmUKJmVyRHesSlFSXdb7LD6XM1ISf8/xWM38P8Keqep2qXgn8EfAFEXn1UgVcKZZ7f5V6Kx77rfbsnRcNh7l2s6OgxRLkCiW62yOUSkqxpJybzrG2s23plY9FQKT8fW2RMNdt6+O52/q4dCDB87av5aO3HeYjtx3mibNTTMzkODedp1iCeCREKlOoev/LubuoJ1s6VwIUEIrqJFaMpbI8dGJi3t+2p93JdnTm7Jz/07kCY1MZ321i11AP73jpDj70mmsuqvxeJ02jR36pZz+oWtrM8fNpQqEQsWiIUEjc/0McP+8/zDRXFx89nUSB3Ru7Z1UTPzwyzc0v2EZfoo3xmTwiwuXrOtm5vpuBrjj9XXG29Hfyhbf8JB967TN5x0t38PrrtzS8qrtH5W+SyhSIRUL0dkTpbncyCbPFEuFQiD/46Z1sG0gsWgHeG7TNVy1+KX2u5rO+fji/Ib4rgV8pX1z1r0XkDPAZEXkDcJfP66wo33jkFB//3hPkiyXWdraRL9S3pgfqHyH6HZHPPq/ADTvXMTqVIZUpcGYqy6EzSaIhoaTCD4+e4+FTztxMPBpesqfSEQ1x77Hx8uLBysW0kzN54m0RRqeyDHbH6OuMkckXOTOV5coq67DKc2iFIvccS5az8v7l7qf48xpHj96zePjEeUanCmhIKeGsOxFgbaJtXg/yyqEeOqLh8ur+sAjFktLb0RaEtPGm0KNaqCf8U4tuqULITeYJh5zfNCTO+36Zq4v5onLdpWvK4T64YFx3DfXw3ldeWU5c8gY72UKJqzZ0z9og0fMCU9k8pybSdMcjXLXBX9hxJaj8TRLxCFk34rCuK85PXLq2fOymqzcuqeJFJSsd4vVroDJAH3DUe0NVvywiAvwz8Ad1S1InB4cn+fh3nTU9azvbymV9dg4mlpxtcnB4ku8cOANKeXJxoCte8wSh34nPuef95qfvI5nOk8oWaQuHSOeLiCpTuVJZQbpiET6w9xC/+OyNsxIOFlvjcMudxyiVSiRiYc5N59h3fILz0zkGe9rpaY+6BWad9UvT2SJ9nU56bCpTIJMvcM+xc7OufWoiTSQED51wDKc3/3PXkXPzhgQXexY/e81GHj4+zo9PT5HO5omEhLZIiGxByRWKVcsiOeGIGXYNddMVj3Dn4VFikTBXbbgwSoaGpY0HXo9qpZ7wTy2T75v6OhhNZsgWS+TcLLiOtigD3fGLzl2ISh17z62PsP90koeKk+VSW9HwhZRzz6B5Sz/WdER55qYe2iJhJtP58jW9DvqK9d3l+1/q9jjLESKv/E0u7e/gvmPnUeDKDV2+E0r8stLLMvwaqAeBFwH3V76pql8SkRDwmbolAUTkRuAvgTDwj6r6fr+f3bt/pFwNQUTKcxDDkxmic2LDfvBGBm3hEKpanlx81ubei9ZNLHSNpTa4g8OTnHTDF46hUNK5IkV3xChASSGdL3DgdJIPjqa44YpBhnriHBt11jhc1t/BTL7EQ8cn+Pb+M7z1xZeVR0ylkmPAY5EQazvbSGYK/OCJc/zMM9bT0x4tJ2u0Rx3DCE56rAhV4/bt0RCPnk6WvTpHSGFNxwUjUuvzuHH3IN/ef4aBrhigREJCoeQMQI6MTnPdtr6LOrO5o+RcscRztq0pV9OAhqaNr4oeQX26VAv1zNnV4n298frNvP9bh+iOR+mKhZnKFpnOFnjj9Zt9yzo3s+3I2ZST6RkLk84VuPuJcbb1d85awuHNH1V+tqc9OquDX44Oejk9kcrfJJUtcN2lfQiQLSrruqLL6tmt9LIMvwbq75hn8lZVv+Aq1/9XjyAiEgb+GngpcBK4T0T+XVUP+Pn8qYl0uRqC10HG3BTq6y/rr1keb2Swe2M3Dzw1QSwixMLCo6eTXDqQWHQEUm+D27t/hJ2DCQ6NOMVeBXFCG0UnVTUcEmIhoaQwlc2TK4QvrHGYyhIW+PHpKS5Z005fZ5RkpsDHv/dE2UANT2ZmGZPueISRkHBweIrB7vZyskbMTTKYTOdBnVpjCmzsjXPvsXGSmTxt4RAbe+Ocn8nT684BZQslsoWSW1svvaTnsWuoh01r25mcyTOWEkAY6onR0eakx87XmVWOkj962+FZo11oaNr4iusR1K9LtbLUZRi1eF83Xb0RcEpmjSQzdMci7N7QzX8cGuPwyPSig5257e/Ow6NMZQpcsT7B2HSeVKZAVzwy7xKOhe6xKx5hdCpTrnmZiIVnZcf5YSWjsvX+AAAcv0lEQVSSTfwmj9TDSi/L8GWgVPUrONlH8x3/PPD5OmV5LnBEVY8CiMjngVcBvpRqY287uXyRw2dTgGOckpkC0XBoSROVlSmpz97Sy5Gz00ymc4QktGI1vyrx9nLKFYokM3kKpRKCE4ePR4RwKERJcdNei3S0XfgpU5kCuWKJkuosA3TOXQEPlJMtPLKFUtnITKbz9HXGygby8nUxprJFZnIFxqfzDCSiHDidpKejja6Ys5XH/tNT7N7QxamJTHn/qd0bu4mGw6zrii75eVw51MNkOs/lgwkeeGqCcMhJlY2GxVeoIkiFMFdJj6BOXVotavW+brp6IzddvXGWsemKR3wNdubOkT55boZYWDg+nubFu5z+oaTK8GSm6uerev/usafGpjl8NkWppExl85w8XyQaDvH333+cmZxeFDH4xiOnyoZ2sDvOG6/fXJcnshze11KjPSutX37TzAeAtwB/qarJOcd6gN8F/kZVz9Uhy0bgRMXrk8B1fj/szT3sWJfgTDLD+LQzZ/HWF19WV6pzT3u0vFrce73SFQkq93LqT8QcDyYcoljKQqFEUUHdChPxaIh8SOjruGBsEvEIpyZm6GwLu0Ylx0yuSDx6IWnzu+/6qRqeRuN4h98T55ktX8m08VpZJT2COnUp6Cw22KnW2c6dI+1sC5Nzy3mNpTL0Jy4u9OrhGYBiscSZZKYcMt/rHv+rX312TfLfxMXr//wmK1RjF/DROj7vXWMX8OFvP1aTgZtPv8CJXtQ7n+Y3xPc2YOtcpQJQ1UkRuRx4O/D/1yxBjYjIzcDNAJs3X4g/Vz6otmiY6y+7UHV4KQ+q3pFBvdlNc/dyikfDdMWd8Nlk2lnrEwk5Ya/L1yUY6HIMaJdbOeGxYedapyfSgBASIb6EubhWoJ5KIMtMYPQI5telWqhnnnWpWbcLDf7m8ybmzpH2dbZxciJNWzjE4yMpouFwWb/n3tPYVIZi0SmkHIuEyiHzVmQpiURz9Ws5d0gQ9ZGnKSIPA29T1TvmOX4D8FequuTVaCJyPfAnqvpy9/UfAqjq/57vM3v27NH7779/vsMXhQJq3ful3iSHpX73u774MEM9cc6lshfFtd/5sh38y91P8eCJCQTh2k09/Nr1W4DZlcc72oRb7nySfNHZ6LArFiUUEnasS7BtIME7XrrD9/158oREAPjuwRHOJjNki8rO9V3lgps97dF54971PI8gbX8hIg+o6p4lfnbF9ci9zrLrUjXq/U3f8fmHyzs+e3OW3o7PC82fePOKlYO/ytfe32OpTLmQa0c0zGgqy7quWHnx9uRMjkhYGE3l6euM0tHm7Iw7kS6wczBRrhLxvcfO0h0PEw5dmLNVddYkvuLqDQB8/ZHTrO1sQ1wdOTwyRXs0RG9HjJde6QyU735ijH3HJ+htj5QXvitKJlcCgaHuOCrOesOueJTdG7p9LfBd6HksNg91cHiSd37hEUpaoqe9je3rOulPxMvhzg+95ppFv3/u9d7xbw8Ds3/Xyn4H/OuRXw/qMuCJBY4fBbb5vNZ83AdcLiLbgFPAL1GxZmQpLMcapqV2hMuR3VRtL6ddQz3zriuae+37npwop6l7lZvXJmLlMKPf+5tbAXpzXzszuSJD8QjXbevz5V3OfR7TmTwnz6d5y2f2lePw3kR4JS22/cVq6BGsgC5Vox79mpt1Wywp56dz3HXEqRix0CBkoejGJ+96kqGeOGOpjJvcFKKvI8r5mQLt0TCZfIl8SemOR9m4oZvHzqQY7IoRcmV4/GyK7vYoh0ac6uX9iThrOqIcH5/h8nUXNibMFhyv79REmjc/fyvffvQMSXdT0GyhREicyhOJikXM49N54lFhOudkxeaLJQolZ41fCGdJ+oAb0r9yQxfZor9FXkuN9ni6FQ0LaKi8aeizt/TOSrevhb37R8oecWU29ZlkhrZo7REcvwYqD2xidly7kkuAunxeVS2IyO8A38ZJjf0nVX20nmuuRmXqxXb5XEomzXJNPF61oafqyMpreH48k2oVoA9NZelPtLF9nbMK3a/x9QziNx45xfu/dcgp/JloI5nO8/5vHQK4yEitVimdVWLF9QhWRpeqUU2/qq2Rq/a8K7NuiyVntB4JQbGkRMOy4CBkocFfuY7e2elyOC+TL9KXaGN9V4zDZ1Nct21teY2cANFIiEjI6UxHkhl3q5soR85O05+Is2uoi2Nj07MMULZQYuvajnIpo7e+6DI+/t0nGJ92vLFnbOjiibEZ1nfFyqWEIiFh52AXD52cREtOUZeSt2xEYSZXIB7tAObfuLMaSx0Me7p11YZu9h13MpXbwsL+U/4ylatxaiJdXodamU09Pp3n+stqN3h+DdQ+4NXAD+c5/gs4azzqwq1J9s1FT/TJSqdALtcun3NZron9hQzd3IKS3pYW/2X7Wn6topDlfBWgL1+XWHJ9sU/ffZzOWIRISDg1mSFXKIHC395x9CID1ejspmVmVfQIll+XquG3tl21512ZdTs6lSGdK1AoOcsnNvbGqy7CrmQ+799r82OpLH0d0VnVH9YmYszknQXelWvkHjnpzE0BtEfDpN0NNJMZZ3lCPBrheZet5di5mXL269a1HYRCFzKEb7p646xq/Rt723nNczbNWiz/1hdfxu0HR3l8JMV0rojnH8XDoAj5YglVXXDjzvlYSrSnMlPZq8k5lc6DaF01PvOFIodGZmdTR0KypGxqvwbqr3Fqhp0EPqGqRQARiQC/g5N99Es1f/sKs9IpkMu1y2c16p3Y9zyHqUzeXVwY4cqhC+VXvC0tcoViObOptz3C/jlbWlSrAL1QOq4fRpIZOttCDE9mneoQ4RCFYpFjY9MXVZ1YrVI6q0RT6tF8zNWvR08nyRaKtJfC3H5whO54lPXd1bc68bJu13e1cXTU6cwEiIaEu46cY/OadobW1D6Q9AZ3f/zvBxhP5ehLtHHVhu5yEtFVG3rKUQ1v7sYrBxSPhumMXSgu3O2msE+m87zimiF+cOQcD56YcELmscisgZz33XPvc2523qUDCR46cR5J5YhHwxRKRQpFZ8F9rqAcHZumr6ONq4a62Lt/ZNbmh/OVGlqOQtbedEItmcrV8H7XnYNOdOXcdI5oOLTkbGpfxWJV9VbgL3CyGc+LyIMi8iAwDnwY+LCqfrnmb19hvMa6nEUOK6ksypjM5IlFQuVyQNC4igWVRTh3DXVz5VA3iVh0VuP1ZD8yeiEUUlJnzcYjJyf4438/wMHhSTb2ts+qOwZw/Nw0x8dnlvzcBrvjnE3miISESFjcUIeQiEUuKn7rt9huNeopZLoSNKsezcdc/ZrK5AmHnAXl3vq4x85McaBK+/A+O1NQOtsiJGIR4m1h2tsixKMhziSznDiXXlKx4V1DTh29Z1zSw5VDjudUrd14bWt9l7N1+5GzKUaSWRKxMJlcgR53fyevRmVbJMxLdg3y3G19zORLS35mr7xmIy/cMcCW/k7i0QjpfJG2cJj+RIyueIST52e498nz3HHoLNEw8xbRXa1C1rXe380v2MbWfifJ5BVXb+Cjr7um6vyyH3zvB6Wq7xGRrwK/ClyOM+D5PvA5Vb13Sd9eJxkfe/zM54nUGv6ZO1LZMdjJ8fEZHjx+nv5EjLBIec8ib2K0URUL5vMc/uXupxjoinNqIs3x8Rly+SJjUxky+SIzuSLpvFObr68jyngqN2vrAXA69+Pnptl3fIJrN/UuOWz2xus383tf+jHt0RAhCZEvKrliiZ/cunjpopUqpbNaBFGPKql1RF6pX6/9+7tJpvMXbXUyma4+rbZrqIfNfR1c2t/BbQfOEhUhHIJCCXKlEjvqqKPpp9145/zL3U+RKZQIh5yK+d3tbfQnYvx3t5L9H936CEdHU7O2iV8sBLkQnpdx5VA3qNIWdsJgoZCQTBdoi4QIu67Dg8cnefaW3qrft1qFrGtlOZd1+F2o2wF8EGejtSjwXeCtqjq2LFIskbklbGqhlh93rjGrrHWXDIWYTOfJ5ArkS0o8EmbX0PIXZazGweHJqinn801e//CJcV50xTqGeuLkC0XuOTbulOOPOjv2qiqFopLMFOhLtM3aesBrxMPJLNdu6mXbQGLR5zYfN129kX/90XEePDlBMlMgHglxzSU9rO/pmGVMPFajlM5qEFQ98qh3zq47HiE5kyfjzuFk3bnF7vj83Yw3iPCSD3LFEmERNq/pYEt/Z13erp92s2uoh4GuODfuHroomWjv/hGOjqb4+o+HCSPE20IUiiX2pfM8c1MPqezS8lkqDcPIVJbBrhg/uX0tR85Ok8kXOTWRpli6UAXmyNlpnlul7uRyJIEFaI1gVfx6UO8F3gR8FkjjpKz+LfCalRHLH3mfaZjVqOXHnWvMzkxl6YxFSBeUZ21xJhcLRaUrHuIZl/SQK2q5oORK/fgHhyf5wN5DHD83QyIWRnGSNM4ks6zvjjGVKcxSuIPDs7ek3tqf4ODwFLlCiULJyZ5yOhJhLJXj2VvWzNp6wLuP3/rn+ziTzPD4aKpcAbqvM+Y7aWHv/hEePT3JqckMvR1tTjoqMDyZJRqZ5t037ly2ZxSkKhIugdQjj3pH5FdtcLc6mco6a/fiEbb0dZQHM9XwBhHd7VGnGogbibhmU88sb7daBKNatf6lzMnM1xc8enqSb+8/Q1soRCgEpdKFEmEHh/1n2VXD06nRqQyPnk5y9xPnGJ7M0BYOkS+VyovqvWSNhTY5DVKEYLnxa6B+HnizWysMEfks8AMRCXsTvY0gGpYlf7aWH3duA05lCnTFwiQzF9YqeYkD861RqkVx/Jy7d/8I49M5EvFIeaQlIoylnP2bPO/S8xzOz+R53vbZW1KHQs580LWbe7nj8CjpnFO9vC0i5QnTyudRWYKpOx4pr5vwFlcuROXoPJnOE3NlLikU1QmdbOiJL7vxCNgIMZB65FHviLwydFXpsS40p1EZZrvL3dLi2s09F1V2qBbBeNbmXjav7bxoJ9xaPcD5+oJkpkChpKzvibkJPU7ty4mZHJlCpO7NCL0lHOdSzjIOVSWdLxIWyBaKnJ/OEY+GaAuHqnr+yxEhCNgyjIvwu6PuJuA/vRdurLwAbFgJofxSLRzkl1omCOcmCiTiEaayzloJj4VGLrVMZvo999REmpxbBsnDCasUyRV11uR1rlCkJx7h3mPn+dHRc4ylMu75YdoiIQa64tywY4DB7jjd7VHWucZp7vPwSjAB5RJMAIdGUosqa+XoPJUt0h2P0N0epTMW4WVXrucFOwZ8L0xsYgKpRx7VEmJqGZEvlpRUSWWC0t79I7z++i38/eufxQ0715EvMuuzlW0nJFKOYJxJZmfthPvpu4/POs/7e7Fdp+frC3raI/R1RgmHQmzojRMJC4WSki8p/2X72ln3tZSEK28JR29HG/FomHg0jAi0tzmDtcl0nol0gd0buqs+x1qedzWWcxfslcKvBxUGcnPeK9Tw+RUhvoSVyR61hH/mjlTWd8UYnkizczBRXoS30MilltCJ33M39rbz+MjUrAVxjtEIlxcPeiGPW+48xrb+Tg6PpJhM57n/yfNcsb6L/kSsrJBrEzF2rEtweCRVzl6a+zxOTaTZ0t9Jws3+S7mpuN0+0lIrR+deWm/lWpNWC03MQyD1yGM5RuR+PNaF5rqqLWx/9PTkrIooo8kM/Ym2ctvx5B1JZrhu2+wogR8PcL6+YO/+EY6Npso1+Db2tpdr8HnlxRa7n4WehacTRVW2rO1ARJjO5jmTzBIOh1ibiPDh117tO0mlVgK4DOMi/CqG4GxLna14Lw78g4jMeG+o6iuXU7iVZrEft9L9bY+GyBeKDE8W2DaQ4OW7B2fFwBea26gldOL33Bt3D/LIyQmOn5sBd4+mVLbItv7Oi7werxF6hmU8lWM4meW9r7yyfM6piTTbBhL89k/Nv15hoRJMi1EZRtk+0Mm+4xNkC6VZa00albywigRaj1Zrzq7WBKXKjTuz+SJJNxw22HNhQDOVKTDYHb9o7tXvwGe+vsDPDglL7egvJIg4C4qLJeXslNM0vF2gV9JQrHSlneUIH/o1UJ+u8t6y7f4ZROaOiqoVwvRbIr+W+S6/5+4a6uHdN+6clcV3/aV9Fy0erGyEc+fLvPP8Npp6RtiVn/XjrbUogdej1ZizqzVBqXLjzlgkVN7b7JpNvbMiGG+8fvOsJRH1Zm1WGuzKHRLmPp+5Ht7cmpfz4enE+u4YD52Y4Px0nlBIGEhEmcoUGElmL1q4vpysZJLFclVx8bth4a8vWdImZTnd31o69lrO3TU0f+FYj4297Tw5luJMMksyky+v7l8sqaGS+TzJWkbYc0fni3lrrcjTUY+qUWuC0ua1bljZLbU10BVjsCfO1v7ERZ7e3HJD9Q58/ERZ5np4+45PlKt3L3ZtTyd+fGqS9rYw7W0h+hJxtq/rJBoOr2i4bSWXYSxX/xmI2HcQWU73t5bQSbVzn7O111fZk2rsGOzk1n0n6Yw5WwlMpvOcnkjzsqv8ZSD58SRreQ5PJ4NkVKeWjtEzZt6mobDwVhLL2cb8ZtPO9fCyhRKHR1L89k9dtuh3ePJeKCl2ITPZq+yyUqxkSHe5+k8zUPOw3O5vLYpTea4fV3khRTo8Ms21m3rLa1O626NueG3aV4iyGSZSjeaingSl1Vps7TdEVc3D645H6OmorZ5do9Y0rdSgcbnuxwzUPASlCoGf7a0XUiQv864y3FDLyGw1tiwxnn747Rgbtdi6lmza+Ty8WghKf7NcLNf9mIGah6BUIVjMQCymSPXOQTX7avWgL0Q0FqcRoeFasmmXoyMOSn+zXCzX/ZiBWoAgzJksZiAWU6R656DqUcBGG4cA7gdlLEKj24xHLdm0y2VYgtDfLCfLcT9+K0kYDWKxiheLrf735qC626NM50p0t0e5dlMvh0emfX3/UlerB2GV+twKBH4rCxiNIQhtxsNPpRmvesQn73oSgDc/fyvvcCugG8uDeVABZ7ER2mIeTr1zUJ4MtSpdEJIrbP6suQhCm/FYTO9awTsPire6EIEwUCLyGuBPgF3Ac1X1/sZKFCwWMhCLKVKj5pCCYByaff5sKTSzLgWhzVSykN4FyZguhWYxsIEwUMB+nErPf99oQZqRhRSpUdlBQTAOrZYZ5ZOm1aUgtBm/BM2Y1kqzGNhAzEGp6kFVPdRoOVqRpc4h1ctKbCddK42690bSzLoUhDbjl3orvzeaUxNpuuZsJBlEAxsUD8o3InIzcDPA5s2bGyxNc9CI7KCgpM22WmbUchI0XQpKm/FDs3vnzeKtrpqBEpHbgfVVDr1HVb/m9zqqegtwC8CePXtafgOhZsaMw8rQyrrULG2mmYxpNZrFwK6agVLVl6zWdxlGK2O6FAyaxZhWo1kMbNOF+AzDMIz6aQYDG4gkCRF5tYicBK4HviEi3260TIbRjJguGa1EIDwoVf0K8JVGy2EYzY7pktFKBMKDMgzDMIy5mIEyDMMwAokZKMMwDCOQmIEyDMMwAokZKMMwDCOQmIEyDMMwAkkg0syN5acZ9noxDMNYCPOgWpAg7UxqGIaxVMxAtSC21blhGK2AGagWpFn2ejEMw1gIM1AtSLNvpmYYhgFmoFqSZtqZ1DAMYz7MQLUgT8etzg3DaD0szbxFaYa9XgzDMBbCPCjDMAwjkJiBMgzDMAKJGSjDMAwjkATCQInIB0XkMRF5RES+IiK9jZbJMJoR0yWjlQiEgQJuA3ar6tXAYeAPGyyPYTQrpktGyxAIA6Wq31FVb2Xpj4BLGimPYTQrpktGKxEIAzWH3wC+Nd9BEblZRO4XkftHR0dXUSzDaDpMl4ymZtXWQYnI7cD6Kofeo6pfc895D1AAPjvfdVT1FuAWgD179ugKiGoYgcZ0yXi6sGoGSlVfstBxEXkT8ArgxapqymIY82C6ZDxdCEQlCRG5EXg38EJVnWm0PIbRrJguGa1EUOagPgF0AbeJyEMi8neNFsgwmhTTJaNlCIQHparbGy2DYbQCpktGKxEUD8owDMMwZmEGyjAMwwgkZqAMwzCMQGIGyjAMwwgkZqAMwzCMQGIGyjAMwwgkZqAMwzCMQGIGyjAMwwgkZqAMwzCMQGIGyjAMwwgkZqAMwzCMQGIGyjAMwwgkZqAMwzCMQGIGyjAMwwgkZqAMwzCMQGIGyjAMwwgkZqAMwzCMQBIIAyUifyYij7hbVH9HRDY0WibDaEZMl4xWIhAGCvigql6tqs8Evg78z0YLZBhNiumS0TIEwkCparLiZSegjZLFMJoZ0yWjlYg0WgAPEXkf8AZgEvipBc67GbjZfZkSkUOrIF4/MLYK37Na2P0snS2r9D1LxnRpVbH7WRq+9EhUV2eAJSK3A+urHHqPqn6t4rw/BOKq+serIpgPROR+Vd3TaDmWC7uf5sZ0KTjY/awsq+ZBqepLfJ76WeCbQGCUyjCChOmS8XQhEHNQInJ5xctXAY81ShbDaGZMl4xWIihzUO8XkZ1ACXgKeEuD5ZnLLY0WYJmx+2ldTJdWF7ufFWTV5qAMwzAMoxYCEeIzDMMwjLmYgTIMwzACiRkon4jIB0XkMbeMzFdEpLfRMi0FEblRRA6JyBER+YNGy1MPIrJJRP5DRA6IyKMi8rZGy2QsjOlRMAmqLtkclE9E5GXA91S1ICJ/AaCqv99gsWpCRMLAYeClwEngPuCXVfVAQwVbIiIyBAyp6j4R6QIeAH6uWe/n6YDpUTAJqi6ZB+UTVf2Oqhbclz8CLmmkPEvkucARVT2qqjng8zipyE2Jqg6r6j737yngILCxsVIZC2F6FEyCqktmoJbGbwDfarQQS2AjcKLi9UkC0AiXAxHZClwL3NNYSYwaMD0KIEHSpaCsgwoEfkrIiMh7gALOKn0jAIhIAvgy8PY5xVKNBmB61LwETZfMQFWwWAkZEXkT8Argxdqck3engE0Vry9x32taRCSKo1CfVdVbGy2PYXrUrARRlyxJwiciciPwEeCFqjraaHmWgohEcCZ3X4yjUPcBv6KqjzZUsCUiIgJ8GhhX1bc3Wh5jcUyPgklQdckMlE9E5AgQA865b/1IVYNWRmZRRORngI8BYeCfVPV9DRZpyYjI84H/BH6MU9oH4I9U9ZuNk8pYCNOjYBJUXTIDZRiGYQQSy+IzDMMwAokZKMMwDCOQmIEyDMMwAokZKMMwDCOQmIEyDMMwAokZKMMwDCOQmIFqAUTkUyKi7r+8iBwVkQ+JSGfFOT8vIt8TkQkRmRaRH4vI+0RknXt8SEQ+526FUBSRTzXshgyjQZguBQszUK3D7cAQcCnwP4D/CnwIQETeB3wReAinxMyVwNuAbcBvu5+PAWPA+wlAkUjDaCCmSwHBFuq2AO4IrV9VX1Hx3j/gKNCrcJTknar6kSqf7VXViTnvfR0YU9U3raTchhE0TJeChXlQrUsaiAK/CkwDH6920lyFMgzjIkyXGoQZqBZERJ4L/ArwXeBy4AlVzTdWKsNoPkyXGosZqNbhRhFJiUgGuBu4E3grII0VyzCaDtOlgGD7QbUOdwI3A3ngtDfKE5HDwH8RkTZ3e2rDMBbGdCkgmAfVOsyo6hFVfWpOCOJzQCfwO9U+JCK9qyKdYTQPpksBwTyoFkdV7xGRDwAfFJFLcHbMPImTFvtm4AjwXgAReab7sW6g5L7OqeqB1ZfcMIKF6dLqY2nmLUC11Ngq57wG+G/AtTgDk2PA14CPeTubiki1xvCUqm5dbpkNI4iYLgULM1CGYRhGILE5KMMwDCOQmIEyDMMwAokZKMMwDCOQmIEyDMMwAokZKMMwDCOQmIEyDMMwAokZKMMwDCOQmIEyDMMwAsn/A7mKEF9ZkNvtAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.subplot(121)\n", "plt.scatter(-z[:, 0], -z[:, 1], alpha=0.5)\n", "for e_, v_ in zip(e, np.eye(2)):\n", " plt.plot([0, e_*v_[0]], [0, e_*v_[1]], 'r-', lw=2)\n", "plt.xlabel('PC1', fontsize=14)\n", "plt.ylabel('PC2', fontsize=14)\n", "plt.axis('square')\n", "plt.axis([-3,3,-3,3])\n", "plt.title('Scikit-learn PCA (flipped)')\n", "\n", "plt.subplot(122)\n", "plt.scatter(yc[0,:], yc[1,:], alpha=0.5)\n", "for e_, v_ in zip(e, np.eye(2)):\n", " plt.plot([0, e_*v_[0]], [0, e_*v_[1]], 'r-', lw=2)\n", "plt.xlabel('PC1', fontsize=14)\n", "plt.ylabel('PC2', fontsize=14)\n", "plt.axis('square')\n", "plt.axis([-3,3,-3,3])\n", "plt.title('Homebrew PCA')\n", "plt.tight_layout()\n", "pass" ] } ], "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.8.5" } }, "nbformat": 4, "nbformat_minor": 2 }