PCA from scratch.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Principal Component Analysis from scratch"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import packages\n",
    "\n",
    "import pandas as pd \n",
    "import numpy as np\n",
    "import cvxpy as cp "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read in file with data\n",
    "\n",
    "cars = pd.read_excel(r'cars.xlsx', header = 3)\n",
    "cars = cars.iloc[:,1:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(32, 11)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "carsdata = cars.iloc[:,1:].values\n",
    "carsdata.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "X (rounded to two digits)\n",
      "[[ 0.15 -0.11 -0.58 -0.54  0.58 -0.62 -0.79 -0.88  1.21  0.43  0.75]\n",
      " [ 0.15 -0.11 -0.58 -0.54  0.58 -0.36 -0.47 -0.88  1.21  0.43  0.75]\n",
      " [ 0.46 -1.24 -1.01 -0.8   0.48 -0.93  0.43  1.13  1.21  0.43 -1.14]\n",
      " [ 0.22 -0.11  0.22 -0.54 -0.98 -0.    0.9   1.13 -0.83 -0.95 -1.14]\n",
      " [-0.23  1.03  1.06  0.42 -0.85  0.23 -0.47 -0.88 -0.83 -0.95 -0.51]\n",
      " [-0.34 -0.11 -0.05 -0.62 -1.59  0.25  1.35  1.13 -0.83 -0.95 -1.14]\n",
      " [-0.98  1.03  1.06  1.46 -0.73  0.37 -1.14 -0.88 -0.83 -0.95  0.75]\n",
      " [ 0.73 -1.24 -0.69 -1.25  0.18 -0.03  1.22  1.13 -0.83  0.43 -0.51]\n",
      " [ 0.46 -1.24 -0.74 -0.77  0.61 -0.07  2.87  1.13 -0.83  0.43 -0.51]\n",
      " [-0.15 -0.11 -0.52 -0.35  0.61  0.23  0.26  1.13 -0.83  0.43  0.75]\n",
      " [-0.39 -0.11 -0.52 -0.35  0.61  0.23  0.6   1.13 -0.83  0.43  0.75]\n",
      " [-0.62  1.03  0.37  0.49 -1.    0.89 -0.26 -0.88 -0.83 -0.95  0.12]\n",
      " [-0.47  1.03  0.37  0.49 -1.    0.53 -0.14 -0.88 -0.83 -0.95  0.12]\n",
      " [-0.82  1.03  0.37  0.49 -1.    0.58  0.09 -0.88 -0.83 -0.95  0.12]\n",
      " [-1.63  1.03  1.98  0.86 -1.27  2.11  0.07 -0.88 -0.83 -0.95  0.75]\n",
      " [-1.63  1.03  1.88  1.01 -1.13  2.29 -0.02 -0.88 -0.83 -0.95  0.75]\n",
      " [-0.91  1.03  1.72  1.23 -0.7   2.21 -0.24 -0.88 -0.83 -0.95  0.75]\n",
      " [ 2.08 -1.24 -1.25 -1.2   0.92 -1.06  0.92  1.13  1.21  0.43 -1.14]\n",
      " [ 1.74 -1.24 -1.27 -1.4   2.53 -1.66  0.38  1.13  1.21  0.43 -0.51]\n",
      " [ 2.33 -1.24 -1.31 -1.21  1.18 -1.44  1.17  1.13  1.21  0.43 -1.14]\n",
      " [ 0.24 -1.24 -0.91 -0.74  0.2  -0.78  1.23  1.13 -0.83 -0.95 -1.14]\n",
      " [-0.77  1.03  0.72  0.05 -1.59  0.31 -0.56 -0.88 -0.83 -0.95 -0.51]\n",
      " [-0.82  1.03  0.6   0.05 -0.85  0.23 -0.31 -0.88 -0.83 -0.95 -0.51]\n",
      " [-1.14  1.03  0.98  1.46  0.25  0.65 -1.39 -0.88 -0.83 -0.95  0.75]\n",
      " [-0.15  1.03  1.39  0.42 -0.98  0.65 -0.45 -0.88 -0.83 -0.95 -0.51]\n",
      " [ 1.22 -1.24 -1.24 -1.2   0.92 -1.33  0.6   1.13  1.21  0.43 -1.14]\n",
      " [ 1.   -1.24 -0.91 -0.83  1.58 -1.12 -0.65 -0.88  1.21  1.81 -0.51]\n",
      " [ 1.74 -1.24 -1.11 -0.5   0.33 -1.77 -0.54  1.13  1.21  1.81 -0.51]\n",
      " [-0.72  1.03  0.99  1.74  1.18 -0.05 -1.9  -0.88  1.21  1.81  0.75]\n",
      " [-0.07 -0.11 -0.7   0.42  0.04 -0.46 -1.34 -0.88  1.21  1.81  2.01]\n",
      " [-0.86  1.03  0.58  2.79 -0.11  0.37 -1.85 -0.88  1.21  1.81  3.26]\n",
      " [ 0.22 -1.24 -0.9  -0.56  0.98 -0.45  0.43  1.13  1.21  0.43 -0.51]]\n",
      "Shat (rounded to two digits)\n",
      "[[ 1.03 -0.88 -0.87 -0.8   0.7  -0.9   0.43  0.69  0.62  0.5  -0.57]\n",
      " [-0.88  1.03  0.93  0.86 -0.72  0.81 -0.61 -0.84 -0.54 -0.51  0.54]\n",
      " [-0.87  0.93  1.03  0.82 -0.73  0.92 -0.45 -0.73 -0.61 -0.57  0.41]\n",
      " [-0.8   0.86  0.82  1.03 -0.46  0.68 -0.73 -0.75 -0.25 -0.13  0.77]\n",
      " [ 0.7  -0.72 -0.73 -0.46  1.03 -0.74  0.09  0.45  0.74  0.72 -0.09]\n",
      " [-0.9   0.81  0.92  0.68 -0.74  1.03 -0.18 -0.57 -0.71 -0.6   0.44]\n",
      " [ 0.43 -0.61 -0.45 -0.73  0.09 -0.18  1.03  0.77 -0.24 -0.22 -0.68]\n",
      " [ 0.69 -0.84 -0.73 -0.75  0.45 -0.57  0.77  1.03  0.17  0.21 -0.59]\n",
      " [ 0.62 -0.54 -0.61 -0.25  0.74 -0.71 -0.24  0.17  1.03  0.82  0.06]\n",
      " [ 0.5  -0.51 -0.57 -0.13  0.72 -0.6  -0.22  0.21  0.82  1.03  0.28]\n",
      " [-0.57  0.54  0.41  0.77 -0.09  0.44 -0.68 -0.59  0.06  0.28  1.03]]\n",
      "Eigenvectors (rounded to two digits):\n",
      "[[ 0.36 -0.02  0.23  0.12  0.14  0.24  0.75 -0.37 -0.02  0.1   0.11]\n",
      " [-0.37 -0.04  0.18  0.14 -0.85  0.05  0.23 -0.06 -0.    0.06 -0.17]\n",
      " [-0.37  0.05  0.06 -0.66  0.05  0.2  -0.   -0.21  0.26  0.39  0.34]\n",
      " [-0.33 -0.25 -0.14  0.26  0.25 -0.58  0.22  0.   -0.07  0.54 -0.07]\n",
      " [ 0.29 -0.27 -0.16  0.04 -0.1  -0.05 -0.03 -0.02  0.85  0.08 -0.24]\n",
      " [-0.35  0.14 -0.34  0.57  0.09  0.36  0.01  0.02  0.25 -0.08  0.46]\n",
      " [ 0.2   0.46 -0.4  -0.18 -0.27 -0.53  0.23 -0.05  0.07 -0.16  0.33]\n",
      " [ 0.31  0.23 -0.43 -0.01 -0.16  0.36 -0.03  0.27 -0.21  0.6  -0.19]\n",
      " [ 0.23 -0.43  0.21 -0.03 -0.18 -0.05  0.06  0.59 -0.03  0.09  0.57]\n",
      " [ 0.21 -0.46 -0.29  0.05 -0.21 -0.   -0.34 -0.61 -0.26  0.05  0.24]\n",
      " [-0.21 -0.41 -0.53 -0.32  0.07  0.17  0.4   0.17 -0.13 -0.36 -0.18]]\n",
      "Eigenvalues (sorted descending and rounded to two digits):\n",
      "[6.82, 2.74, 0.65, 0.28, 0.23, 0.22, 0.14, 0.13, 0.08, 0.05, 0.02]\n",
      "Z matrix (rounded to two digits)\n",
      "[[ 0.66 -1.74  0.6  -0.18 -0.07  0.15  0.01  0.43  0.12 -0.96  0.02]\n",
      " [ 0.63 -1.55  0.38 -0.09 -0.13  0.08  0.09  0.42  0.2  -1.03  0.25]\n",
      " [ 2.78  0.15  0.24  0.1   0.05 -0.13 -0.59  0.62 -0.25  0.41  0.35]\n",
      " [ 0.31  2.36  0.14 -0.15 -0.06  0.22  0.05  0.04 -0.51  0.56 -0.02]\n",
      " [-1.97  0.75  1.13 -0.15 -0.06 -0.02  0.16 -0.39  0.08  0.21 -0.15]\n",
      " [ 0.06  2.79 -0.16 -0.02 -0.2  -0.03 -0.26  0.3  -0.99  0.22  0.25]\n",
      " [-3.   -0.33  0.36 -0.17  0.37 -0.18  0.17  0.14 -0.05  0.35 -0.72]\n",
      " [ 2.06  1.47 -0.94  0.02  0.3   0.38 -0.17 -0.65 -0.14 -0.32  0.  ]\n",
      " [ 2.29  1.98 -1.8  -0.16 -0.12 -0.88  0.11 -0.63  0.29 -0.34  0.34]\n",
      " [ 0.53  0.16 -1.49  0.13 -0.22  0.55 -0.1  -0.16  0.07 -0.07 -0.83]\n",
      " [ 0.51  0.32 -1.68  0.04 -0.35  0.31 -0.2  -0.09  0.1  -0.15 -0.74]\n",
      " [-2.25  0.68  0.38  0.4  -0.07 -0.06  0.2   0.02 -0.13 -0.38 -0.13]\n",
      " [-2.05  0.68  0.48  0.2  -0.11 -0.21  0.33 -0.05 -0.21 -0.36 -0.24]\n",
      " [-2.15  0.8   0.3   0.14 -0.22 -0.39  0.12  0.07 -0.18 -0.44 -0.18]\n",
      " [-3.9   0.83 -0.65 -0.27  0.06  0.2  -0.14  0.17  0.3  -0.05  0.9 ]\n",
      " [-3.95  0.73 -0.72 -0.04  0.12  0.2  -0.13  0.2   0.41  0.    0.88]\n",
      " [-3.59  0.42 -0.55  0.23  0.28  0.28  0.4  -0.04  0.68  0.21  0.66]\n",
      " [ 3.86  0.3   0.42  0.21 -0.03  0.11  0.64  0.04  0.06  0.22  0.48]\n",
      " [ 4.25 -0.69  0.21 -0.25 -0.16  0.24  0.4   0.26  1.19  0.1  -0.53]\n",
      " [ 4.23  0.28  0.47  0.03 -0.13 -0.11  0.87 -0.07  0.19  0.23  0.32]\n",
      " [ 1.9   2.12 -0.16 -0.06  0.52 -0.45 -0.21  0.28  0.05  0.04 -0.74]\n",
      " [-2.18  1.01  1.17 -0.05 -0.14  0.11 -0.32 -0.1  -0.59 -0.23 -0.11]\n",
      " [-1.86  0.91  0.96 -0.05 -0.3  -0.12 -0.33 -0.08  0.01 -0.26 -0.29]\n",
      " [-2.89 -0.68  0.16  0.1   0.33 -0.05 -0.04  0.21  0.83  0.4  -0.96]\n",
      " [-2.25  0.87  1.04 -0.12  0.02  0.21  0.24 -0.48  0.15  0.3   0.2 ]\n",
      " [ 3.57  0.12  0.45  0.01 -0.09 -0.02 -0.09  0.36 -0.01  0.21  0.15]\n",
      " [ 2.65 -2.05  0.83  0.11  0.36 -0.13 -0.65 -0.83  0.58 -0.61  0.34]\n",
      " [ 3.39 -1.38  0.45 -0.03  0.24  0.3  -0.   -0.52 -1.17  0.71 -0.02]\n",
      " [-1.37 -3.5   0.14  0.04 -0.51 -0.24 -0.62 -0.42  0.6   1.12  0.18]\n",
      " [ 0.   -3.22 -0.4  -0.04  0.15  0.1  -0.01 -0.03 -0.95 -0.86  0.01]\n",
      " [-2.67 -4.38 -1.35 -0.06  0.05 -0.35  0.57  0.19 -0.89  0.46  0.02]\n",
      " [ 2.42 -0.23 -0.41  0.15  0.12 -0.05 -0.48  0.79  0.23  0.33  0.33]]\n"
     ]
    }
   ],
   "source": [
    "# PCA implementation\n",
    "\n",
    "# 1. Shifting\n",
    "\n",
    "# calculate the mean of each column\n",
    "M = np.mean(carsdata.T, axis=1)\n",
    "# center columns by subtracting column means\n",
    "C = carsdata - M\n",
    "\n",
    "# 2. Scaling + 3. Denote X\n",
    "\n",
    "# calculate standard deviation\n",
    "SD = np.std(carsdata.T, axis = 1)\n",
    "# divide columns by standard deviation\n",
    "X = C/SD\n",
    "print(\"X (rounded to two digits)\")\n",
    "print(np.round(X, 2))\n",
    "\n",
    "# 4. Covariance matrix\n",
    "\n",
    "# calculate covariance matrix of centered matrix\n",
    "S = np.cov(X.T)\n",
    "print(\"Shat (rounded to two digits)\")\n",
    "print(np.round(S, 2))\n",
    "\n",
    "# 5. Eigendecomposition\n",
    "\n",
    "# eigendecomposition of covariance matrix\n",
    "values, vectors = np.linalg.eig(S)\n",
    "print(\"Eigenvectors (rounded to two digits):\")\n",
    "print(np.round(vectors, 2))\n",
    "print(\"Eigenvalues (sorted descending and rounded to two digits):\")\n",
    "print(sorted(np.round(values, 2), reverse = True))\n",
    "\n",
    "# 6./7. Define Q and Z\n",
    "# project data\n",
    "Z = vectors.T.dot(X.T)\n",
    "print(\"Z matrix (rounded to two digits)\")\n",
    "print(np.round(Z.T, 2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Percentages of variation captured by each PC in ascending order of PC index (rounded to two digits)\n",
      "[60.08, 24.1, 5.7, 2.45, 2.03, 1.92, 1.23, 1.12, 0.7, 0.47, 0.2]\n"
     ]
    }
   ],
   "source": [
    "percentages = sorted(np.round(100 * (values/np.sum(values)),2), reverse = True)\n",
    "print(\"Percentages of variation captured by each PC in ascending order of PC index (rounded to two digits)\")\n",
    "print(percentages)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cummulative percentages of variation captured by each PC in ascending order of PC index (rounded to two digits)\n",
      "[ 60.08  84.18  89.88  92.33  94.36  96.28  97.51  98.63  99.33  99.8\n",
      " 100.  ]\n"
     ]
    }
   ],
   "source": [
    "cumm_percentages = np.empty(len(percentages))\n",
    "for i in range(len(percentages)):\n",
    "    cumm_percentages[i] = sum(percentages[:i+1])\n",
    "print(\"Cummulative percentages of variation captured by each PC in ascending order of PC index (rounded to two digits)\")\n",
    "print(cumm_percentages)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Standard deviation of each PC in ascending order of PC index (rounded to two digits)\n",
      "[2.57, 1.63, 0.79, 0.52, 0.47, 0.46, 0.37, 0.35, 0.28, 0.23, 0.15]\n"
     ]
    }
   ],
   "source": [
    "SDs = sorted(np.round(np.std(Z.T, axis = 0),2), reverse = True)\n",
    "print(\"Standard deviation of each PC in ascending order of PC index (rounded to two digits)\")\n",
    "print(SDs)"
   ]
  }
 ],
 "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
}