3 - PCA/2-MeanCovariance.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Mean/Covariance of a data set and effect of a linear transformation\n",
    "\n",
    "We are going to investigate how the mean and (co)variance of a dataset changes\n",
    "when we apply affine transformation to the dataset."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Learning objectives\n",
    "1. Get Farmiliar with basic programming using Python and Numpy/Scipy.\n",
    "2. Learn to appreciate implementing\n",
    "   functions to compute statistics of dataset in vectorized way.\n",
    "3. Understand the effects of affine transformations on a dataset.\n",
    "4. Understand the importance of testing in programming for machine learning."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, let's import the packages that we will use for the week"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PACKAGE: DO NOT EDIT THIS CELL\n",
    "import numpy as np\n",
    "import matplotlib\n",
    "matplotlib.use('Agg')\n",
    "import matplotlib.pyplot as plt\n",
    "matplotlib.style.use('fivethirtyeight')\n",
    "from sklearn.datasets import fetch_lfw_people, fetch_olivetti_faces\n",
    "import time\n",
    "import timeit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from ipywidgets import interact"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we are going to retrieve Olivetti faces dataset.\n",
    "\n",
    "When working with some datasets, before digging into further analysis, it is almost always\n",
    "useful to do a few things to understand your dataset. First of all, answer the following\n",
    "set of questions:\n",
    "\n",
    "1. What is the size of your dataset?\n",
    "2. What is the dimensionality of your data?\n",
    "\n",
    "The dataset we have are usually stored as 2D matrices, then it would be really important\n",
    "to know which dimension represents the dimension of the dataset, and which represents\n",
    "the data points in the dataset. \n",
    "\n",
    "__When you implement the functions for your assignment, make sure you read\n",
    "the docstring for what each dimension of your inputs represents the data points, and which \n",
    "represents the dimensions of the dataset!__. For this assignment, our data is organized as\n",
    "__(D,N)__, where D is the dimensionality of the samples and N is the number of samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Shape of the faces dataset: (4096, 400)\n",
      "400 data points\n"
     ]
    }
   ],
   "source": [
    "image_shape = (64, 64)\n",
    "# Load faces data\n",
    "dataset = fetch_olivetti_faces('./')\n",
    "faces = dataset.data.T\n",
    "\n",
    "print('Shape of the faces dataset: {}'.format(faces.shape))\n",
    "print('{} data points'.format(faces.shape[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When your dataset are images, it's a really good idea to see what they look like.\n",
    "\n",
    "One very\n",
    "convenient tool in Jupyter is the `interact` widget, which we use to visualize the images (faces). For more information on how to use interact, have a look at the documentation [here](http://ipywidgets.readthedocs.io/en/stable/examples/Using%20Interact.html).\n",
    "\n",
    "We have created two function which help you visuzlie the faces dataset. You do not need to modify them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def show_face(face):\n",
    "    plt.figure()\n",
    "    plt.imshow(face.reshape((64, 64)), cmap='gray')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6f140be6814c45568cf96c4910303c4d",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=0, description='n', max=399), Output()), _dom_classes=('widget-interact'…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "@interact(n=(0, faces.shape[1]-1))\n",
    "def display_faces(n=0):\n",
    "    plt.figure()\n",
    "    plt.imshow(faces[:,n].reshape((64, 64)), cmap='gray')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Mean and Covariance of a Dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this week, you will need to implement functions in the cell below which compute the mean and covariance of a dataset.\n",
    "\n",
    "You will implement both mean and covariance in two different ways. First, we will implement them using Python's for loops to iterate over the entire dataset. Later, you will learn to take advantage of Numpy and use its library routines. In the end, we will compare the speed differences between the different approaches."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# GRADED FUNCTION: DO NOT EDIT THIS LINE\n",
    "def mean_naive(X):\n",
    "    \"Compute the mean for a dataset X nby iterating over the data points\"\n",
    "    # X is of size (D,N) where D is the dimensionality and N the number of data points\n",
    "    D, N = X.shape\n",
    "    mean = np.zeros((D,1))\n",
    "    ### Edit the code; iterate over the dataset and compute the mean vector.\n",
    "    for d in range(D):\n",
    "        for n in range(N):\n",
    "            mean[d] = mean[d] + X[d,n]/N\n",
    "\n",
    "    return mean\n",
    "\n",
    "def cov_naive(X):\n",
    "    \"\"\"Compute the covariance for a dataset of size (D,N)\n",
    "    where D is the dimension and N is the number of data points\"\"\"\n",
    "    D, N = X.shape\n",
    "    ### Edit the code below to compute the covariance matrix by iterating over the dataset.\n",
    "    covariance = np.zeros((D, D))\n",
    "    ### Update covariance\n",
    "    X = (X - mean_naive(X)).T\n",
    "    covariance = np.dot(X.T,X)/N\n",
    "\n",
    "    return covariance\n",
    "\n",
    "\n",
    "def mean(X):\n",
    "    \"Compute the mean for a dataset of size (D,N) where D is the dimension and N is the number of data points\"\n",
    "    # given a dataset of size (D, N), the mean should be an array of size (D,1)\n",
    "    # you can use np.mean, but pay close attention to the shape of the mean vector you are returning.\n",
    "    D, N = X.shape\n",
    "    ### Edit the code to compute a (D,1) array `mean` for the mean of dataset.\n",
    "    mean = np.zeros((D,1))\n",
    "    ### Update mean here\n",
    "    mean = np.mean(X, axis=1, keepdims=True)\n",
    "\n",
    "    return mean\n",
    "\n",
    "def cov(X):\n",
    "    \"Compute the covariance for a dataset\"\n",
    "    # X is of size (D,N)\n",
    "    # It is possible to vectorize our code for computing the covariance with matrix multiplications,\n",
    "    # i.e., we do not need to explicitly\n",
    "    # iterate over the entire dataset as looping in Python tends to be slow\n",
    "    # We challenge you to give a vectorized implementation without using np.cov, but if you choose to use np.cov,\n",
    "    # be sure to pass in bias=True.\n",
    "    D, N = X.shape\n",
    "    ### Edit the code to compute the covariance matrix\n",
    "    covariance_matrix = np.zeros((D, D))\n",
    "    ### Update covariance_matrix here\n",
    "    covariance_matrix = np.cov(X, bias=True)\n",
    "\n",
    "    return covariance_matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's see whether our implementations are consistent"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "X:\n",
      " [[0 1 2]\n",
      " [3 4 5]]\n",
      "Expected mean:\n",
      " [[ 1.]\n",
      " [ 4.]]\n",
      "Expected covariance:\n",
      " [[ 0.66666667  0.66666667]\n",
      " [ 0.66666667  0.66666667]]\n"
     ]
    }
   ],
   "source": [
    "# Let's first test the functions on some hand-crafted dataset.\n",
    "\n",
    "X_test = np.arange(6).reshape(2,3)\n",
    "expected_test_mean = np.array([1., 4.]).reshape(-1, 1)\n",
    "expected_test_cov = np.array([[2/3., 2/3.], [2/3.,2/3.]])\n",
    "print('X:\\n', X_test)\n",
    "print('Expected mean:\\n', expected_test_mean)\n",
    "print('Expected covariance:\\n', expected_test_cov)\n",
    "\n",
    "np.testing.assert_almost_equal(mean(X_test), expected_test_mean)\n",
    "np.testing.assert_almost_equal(mean_naive(X_test), expected_test_mean)\n",
    "\n",
    "np.testing.assert_almost_equal(cov(X_test), expected_test_cov)\n",
    "np.testing.assert_almost_equal(cov_naive(X_test), expected_test_cov)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now test that both implementation should give identical results running on the faces dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.testing.assert_almost_equal(mean(faces), mean_naive(faces), decimal=6)\n",
    "np.testing.assert_almost_equal(cov(faces), cov_naive(faces))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the `mean` function implemented, let's take a look at the _mean_ face of our dataset!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPcAAAD3CAYAAADBjMJTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztnX+MHdWV57/H/cPdbbfbbkMTB1vYg5txEKxJhLL8SEYzMBllsyNAERmB0MasUDZKZhGrYbU4s9FKI+1qGa0yySgZZXYECY7EQFiGLBaazMRiyGwQCwQSjE0gNAskGBvb4O62u9v9k7t/dNXLqeM65533+vVrKM5Hslz1btWt+6rq9vuee849l1JKCIKgeqxa6QYEQbA8ROcOgooSnTsIKkp07iCoKNG5g6CiROcOgoqypM5NRJ8mol8S0StEtLtVjQqCYOlQs35uIuoA8DKATwE4BOCnAG5MKf0CAMbHx8OBHgRtZGBggPj+Un65Pw7glZTSqymlWQD3A7h2KY0LgqB1LKVznwvgDbZ/KPssCIL3AJ1LOJdKPiuV4iMjI0u4TBAEGsPDw2rZUjr3IQBb2P5mAIe1g++6667C/tlnn126DQCDg4Ol2+vXry8ct3bt2tp2b29voay7u7u23dXVBQA4fPgwtmzZUjiuo6OjdBsAVq1aVVrGP5f7RMW/eXw/396/fz927tx5xrEerDGSd999t6HzDh48iIsuusjdDu9x/FryuryN+fZLL72EHTt2lJbJ7XplCwsLpdvz8/OF4/h+ftypU6fQ39+Pubk59TxeNjMzU9uenp4uHDc1NVXbPnnyZKFsdHS0tn38+PHa9okTJwAAN998M+65557avjyHH/ujH/0IGkuR5T8FMExE24ioG8ANAPYuob4gCFpI07/cKaV5Ivr3AP4RQAeA76SUXmhZy4IgWBJLkeVIKf09gL/3HJtL49qFOztLt+WxvMySw94yeZwlt7Wypcry/Jxm5LAly3k75HHaeR0dHWo7GjEbtDZK2azdD2kSWdfm9Vv3m2M9M95GImrJe6WZdEDxnebmo+wTvB80cn9qbah7RBAE70uicwdBRVmSLG8EKSv4vpTlmsy1pIhXnlk0Mxpctr/U+r3S3mqXt02W3JPXsuSwdj3vd/F+r3pl2j3gI+dAUYrLc7Qyue99npa5Z8lyy2QMWR4EH2CicwdBRYnOHQQVpW02t7QZuH1hRYY1awdrUUzS9tLOke1oxt1lnTc7O+seW5AuNK3MwuvuamRcg6PZulZ0GY/+mp+fX1J0Wdk+3+aRZbIO3qaZmRm1jfWureF1k1k2txyXkvul13W1LgiC9x3RuYOgorRNllsRapZsseDyT8onrUwG+Hslr9dU8MryiYmJpqLoLJei1Q6tbGFhQTUJvG4mQDd9LNksnwuXztq23Lcmjnjle15HV1cXJicnC2Xye2rXtt4/y5WnvWOrVq0q7IcsD4KgRnTuIKgoKybLregbbbRVSisukbwjl+Pj42pZM1Ft9erQZNfx48dVuS0lV9nc9LJ9S9pz+LXm5ubU723Jcq/ctiQ1nw89Ojpa2J+dnS09TpZZkWdl0rvsu+QMDQ1hdHTUnBev1WmZKVKyazQSoaaZY5z45Q6CihKdOwgqSnTuIKgobbO5pR1p2QyabSdtL8tlUWYP9ff34+2331aPk2g2lRV1JSmzqwcHB3Ho0KGCvcy3V69eXaijp6endFsea03u1+z7mZkZ1ea23ExWxJdlL3NXZJ5nrKenB0ePHsXp06dLj5PuS77fyGwvTtlYyNDQEI4ePeqeGafVV+/a/N7xMjlmslR3cfxyB0FFic4dBBVlxWS5le+LSy0u8SyXgmcCQX9/P44ePVo4TptAIMusqCt+nuWyyO/B4OAg3njjjYKktqT3mjVratt9fX2FMp7SmZ9nRQTyNk1NTakSz5o0wZ8LoMtoLrUBYHJysnBtANi6dSveeuutQhk/T9bhdYVpedKAonuR348TJ064k0hwvFGVgG7eSClvRaiFKywIPsBE5w6CihKdOwgqyoola+BIe4jbVNK241juBs11derUqcJx1uwjmUxAO85yhZXl5t65cyeOHDlSsLm57Sztam5z822guKQSL5N2O7fBub05PT2t2pXye3K3lnRxcXuZ32NpL09MTJyxvXXrVhw5cqSwBI/mFgPsd4JjJQTh7efHTUxMmLPwON4klFboqzaOMTs7WzjOOxuQE7/cQVBR6nZuIvoOER0jooPss0Ei2kdEI9n/G5a3mUEQNIpHlt8D4FsAvsc+2w3g0ZTSnUS0O9u/w6rEiuCRLhdtRo1XBsnr8fOkpLNWdNRmOsk6LDcZP49L41OnTqnuIysiS15bc9dJKchlOr8fMzMzqltSynLLxcWluLYNFOU7l+hjY2OqLLfut0SbXcdNEUA3pebn580EHp46APs+au4vGdnn/Z4adX+5U0r/B8AJ8fG1APZk23sAXFf3SkEQtJVmbe5zUkpHACD7f6h1TQqCoBWQR+oS0VYAj6SULsr2x1JK61n5aEqpYHePj4/XKh4ZGWlVe4MgYAwPD9e2BwYGCrZvs66wo0S0KaV0hIg2AThW74R9+/YV9rnrR4bWaXZkI/nNy2zuyy+/HE888UThOO/MMt4mKxzSY199/vOfx/e+973C9+b2OL83QNHF1d/fXyhbt25dbXtwcLC2zV1kQNG9ltuf09PT6OnpKbSDf2dp6/Lvze1jYNFmLts+efJk4Thug+d1XH/99XjwwQcL7ilu38uxEMsm1paGljZ32XHXXHMN9u7day4vzWnm3QGK964sv37+fljvXJ5R6Lvf/a7avmZl+V4Au7LtXQAebrKeIAiWCY8r7D4A/xfAbxPRISK6BcCdAD5FRCMAPpXtB0HwHqKuLE8p3agUXd3IhayZX5b7yOuW8Ep2GbmltUnuc/lkufVkHTJHOK9PSvgc+V0seaklRpBRbt4kBtY9tpb44e3gEpK7vgDdrTc7O+u+31YEGb8/1iw5LbFkT09PoU7rWWtLVtUr42jvupVcAoi85UHwgSY6dxBUlBWbOGLJLm8dVp5uLerKG8QP6JJXjhTzEWBrkoOcwMK/N5eTlnvSSgbB65Aj7nySiiUZrWvzY+X31EbS5b2yZDm/x1oUF2DLcn5/+LUs2cwl+vz8vJl8Q5PY3uNkm7m81vKpefZL21T3iCAI3pdE5w6CihKdOwgqyorZ3N4lZvm2dxYOULSfuS0n85bz86R9ODo6Wtvma4zJ47Q1rgA9Yur48eOF47i9LKPQBgYGoMHvK7ezG5lJZblgOJbLT1sDzEqkKBM8cLvdWqZXzmrjaK5TaZvze7Vhw28ip998881CRKDMIc+vXZaIo147yo7NsaIxmyF+uYOgokTnDoKK0jZZ3oh7xxuVZiVQ0PJ4vfXWW4Xj+HlchgPAK6+8UnqeFYUm+fCHP1zb/tCHPlTbTimp8lJ+Z2tSiZbnzTJhpAtRi+bzLkkk0Sb+yDKZk42bO1qSC1n26quvFsr4RBX+XaR85+0/77zzAABf+MIX8Mwzz9T2gTNz1vG2cDNLPhe+73Xh8rp7enoKJky4woIgqBGdOwgqSttkuZSufHRYlmkRQtYIrZw3/M4775SWvfnmm4XjpOzi8BxfXLJb+cnkxA4+3/rcc88tbPN2cdNBTgrQlh0qOzbHmhCSnzMzM4POzk5VYsuIPWsVUU3OW8sayVVOeRv5feSmjfwuhw4dKpQdO/ab1AL8/ZARdfz+cG/E2NgYNm/eXNuXy0/xZ8/nzEuPBr8HchKPdr95mzo7O1UTwKqDE7/cQVBRonMHQUWJzh0EFWXFbG5reVLNFSaTG1jL+2pLEsmII24HSzuGR7NxV4pMQMDdHhs3biyUXXDBBbXt888/v7a9ffv2gj3HI+Ckm4Pb7evXry+U8e/Dt63ZUjJqjtt2VrQaHxewEkpwG9ObZEBG5fGxBW4DA8VxEjmGwm1rPu4ix0J4Hfy5nH/++di2bZtaP7fp+ZiMfK/4+JB8ntoMPU53d3fh3skxAzmWUUb8cgdBRYnOHQQVpW2yXMoxK7qHSygu/2SUEUemfuWyiMva7du3F47jslxKN24G8MkFUpZzicTTCwMoRDtxudrX11eo03KncSkuZTm/tibRZZ1yW1tqSEpGXqeVKrgsjXK9sqGhoYJrSXOLAcV7deGFF6pt5JNzpEnH7z2X5RdeeCHOOecc9dpcznNzT7q7eDukhOZlWj671atXF9596QL1rHQav9xBUFGicwdBRYnOHQQVpW02t2V3cJsYKNrg3B0g7QzLxpQuqZytW7cW9rkNJdu4Y8eO2vZZZ51V25ZuCW4fynZwO1KG1fJr8+8sxyC4m0guE6RdW9ah2fTd3d2Ffd5GGe5r5U/n9Vuz07hdyevYuHFjwW61bEpeNjRUXIOSt4PPDJSuUn4tPk6yefNm873iz4K3Q7q7rIQP/D3T7nd3d7cZxivvfxnxyx0EFcWznNAWInqMiF4koheI6Lbs80Ei2kdEI9n/G+rVFQRB+/DI8nkAt6eUfkZE/QCeJaJ9AG4G8GhK6U4i2g1gN4A7tEqkjOCyUboR+L4mYYCi1OKuDaDo+uDHnX322YXjrIgsPtOHSyspGb15wOXMKc0kkPeDS00pt+VSODlythu//zJCzZt3zFpVVZPlUg5rubn7+voK19PyvQN29Bc38aQJw9FMjDVr1hTKpAtK3n+tHVa0oJaLjpt73d3dhXsn+09LZHlK6UhK6WfZ9ikALwI4F8C1APZkh+0BcF3dqwVB0DYasrmJaCuAjwJ4CsA5KaUjwOIfAABD+plBELQbsnKbFQ4kWgvgnwH8t5TSQ0Q0llJaz8pHU0o1bTw+Pl6reGRkpIVNDoIgZ3h4uLY9MDBQsJVcrjAi6gLwdwDuTSk9lH18lIg2pZSOENEmAMf0GoADBw4U9rk9K+1lbi9ym8e7vCxQnju7p6en4B4BijahDFHk+9wGlMdpS88CRds0t+UuuOACvPzyy2pIqHSdWPYbr58fJ+1NbsfnZa+99hq2bdtWsPf5H3s5tsAzx/AZV0AxUw2fLWWNQeT3bXh4GCMjI6qtLu+319XGkWMEZVlOdu7cif379xdsbivnOD/OmiUnr60licxt7osvvhgHDhwo2OBjY2OFOvL3+I471GEu12g5AbgbwIsppb9gRXsB7Mq2dwF4uF5dQRC0D88v95UA/g2AA0T0XPbZnwK4E8ADRHQLgF8D+NzyNDEIgmao27lTSo8D0BJzX+29kJV32ZI0mnSVdVgzmHJZNzMzc0Y0nLX0DZeD/NrW0jzye2oJATdt2qQuNWR9Tyk7eVu4m0bWoeUc7+joUO+jlJPacsFA8bt53DSy/v7+fnWGlHSn8edimUiaGxLQ70dfX5+ZmFB7npa5JJ+ZFn0nZ4VZudu1pBeciFALgooSnTsIKsqK5VCzyjTJ1Ig04cfmUq2eLLeWvrGi0KwleLR84WeddVbpSLo8B7BHsDU5L9uh5UmT994q00aKZf3a6LusUy6TpI0wWxFqUrJbnguOtqTPhg0bzO+ptdHKI26ZcbyMmxgdHR2q2SavrRG/3EFQUaJzB0FFic4dBBWlbTa3FV1m2WVel5m17K22vhNQjICTtpFlZ2vttWzYvI3Hjx8/w/a3bChuZ1t5wK0ldj1Lvsr6re8i6+f2vjYjz6qjo6Oj8Cy0MQJAn1UF6M/M+45t2LDBXEJamxlnrXcn28j3+bOV98PbDo345Q6CihKdOwgqyorJciu6TIuSamRJIg05+d5ynWgRUxJLUmvnyckhHBl1xeuwTBNe5nVxNeKi5FiReNaz0KLGLHNDukC9JoYly8uOm56expo1a8z6tXtiRQ7KMs18ku45K8IuXGFB8AEmOncQVJTo3EFQUdpmc0t71pt4wWs3yfo0m0Tag9zmsdxd3ow18jht7a2Ojg7VpSPdR9ZMLa/bsNVYoalem1WOaXhdj1aSBGs2FqesLLe5vWMt1rpqWpssrOt63z9O/HIHQUWJzh0EFaVtstzKX21F93C3kBXpI2WRJkstiWRJWa88k8dZbj0tl1YjslzDe5ykWfNDm+Em3XqamZVSUk0TaUrJZZm0dlh51jnS1eiVx5YUt8o87Ugpud8rjfjlDoKKEp07CCpK22S5NxECoE8MkHLVK/+saDhOM6Oact+KtpO5v06fPl3b57Lcup51bW875LZ2nmVKWfV7osHKPufXs1bQ1NoE2BFfGo2YMNo9kO93M+aNd1kqwF4FNSd+uYOgokTnDoKKEp07CCrKikWocRvFmnBvTczn+143lrQjOVb0l2X3Wt+F29ncrp6cnCzY3NaytFYb+ff21sGRthv/btYyPtb3tGZEeZ8tr6+RRAhaLnFrVqJMQGnZ4Nq1ve+wVZ/lGrSWjdaIX+4gqCietcJ6iOhpItpPRC8Q0Z9ln28joqeIaISIvk9EvmUmgiBoCx5ZPgPgqpTSRLba5+NE9EMAfwLg6yml+4norwHcAuDbWiVS6lguLk3+NZK/mksyS/55I5esOjQXjtzn25OTk4V9KyGDlUtLq8OK2OPtnZ6edrt3tFVPgaLJYeWC13J4z8/Pq1GLVvIKS/ZbslxLbDE/P+9+JyzXlbVKKb932nFzc3Pu/OwadX+50yL5mqxd2b8E4CoAD2af7wFwXd2rBUHQNsjjbCeiDgDPAtgO4K8A/A8AT6aUtmflWwD8MKV0UX7O+Ph4reKRkZEWNzsIAmBxbfOcgYGBguRwjZanlBYAXEJE6wH8AMBHyg6z6njuuecK+/39/bxRhTKefphv8wXkgWIeMplnq0yS/epXv8J5551XOG4lZPnc3By6urpaIsu15Zaseev59rFjxzA0NNSULOcj/cCimZHTqCwfGhrCsWPH3N/Fyq+mzfX2yPKJiQmsXbu2KVneyEqkU1NTte2JiYna9smTJwEAO3bswEsvvVQoO3HiRKGO0dFRAMBXv/pVta0NucJSSmNE9GMAlwFYT0SdKaV5AJsBHLbOtcLpLJvbyoFdth5Y2fX4g7XyaHvDSht5kJZNxfct14aV+FCzMWUCRi1p4cmTJ90hrFbn5i/s9PS0Wgevnz+XmZmZMzpqjjdBhdxvZmac5SoF/KHR1vgE39feHfl+LIvNTURnZ7/YIKJeAL8P4EUAjwG4PjtsF4CH614tCIK24fnl3gRgT2Z3rwLwQErpESL6BYD7iei/Avg5gLuXsZ1BEDRI3c6dUnoewEdLPn8VwMe9F5LShMsMS5Zb9qwVlaaVeeRMjjcayfouZdKqs7MT8/Pz6nmWu85yB/Lj+PI++TXlcb29vTXbLceS5dKFxuGynNvc1nLHXEJPTk4WTAnrOVmJPjj82lbCB16HZd7Jfa87txWy3Jr1qBERakFQUaJzB0FFeU9MHPHKckvyWitBcvnnlXuyzVrSBblvjZZLWc7bok0wkW22orWs3GJlrrAtW7bg7bffhoZ3dUqgKNN5+61ca1yGy9F371I6zbgy5XmWN8WaJGR5O/j9kc9Ti1qU21okm2yHRvxyB0FFic4dBBUlOncQVJQVs7mtWTOaDWuFVFo0YlNxeJst+8f6LpoLbWFhwR3eaiVS1K4t7ze/d/n92LJlC44fP67eA8st5HH5AfbSvJyZmRm1TH7O67eSUlihtJz8PZIhwWVoLkt5PzRbWu7z5xk2dxAELqJzB0FFWbHlhKyILC5HtJxYgD+fNT/Pkl2W68eS3lYOds3lsmrVKtV1ZU0OaTZJAoffm3feeaehqL2yOgD9HkhJrT0XGZFlrb5qRQt6V+iUCRqAxXfNyh8v22JFoVmuTe2ZSVnuXUpLI365g6CiROcOgooSnTsIKsqKrRVmTUTX7BVruVZvojxp/3hnQXmT1Vl5xWXYJLdHrdBRrb2yLdYywJqrbXJyUs0zLtthLSWshYtaz0yizXSynq3X5rYSceSsXbsW09PT7jZaodFeV5i17XUXa8QvdxBUlOjcQVBR2ibLrcnmlpvMkuXe6KRcGnd0dJjL51hJGCz3i9Uma0lZLQmDlRBQoklIr4kxNzenJmq05Kks49+tLDFE2b7cbma2l/UsrHxommRvdjkhKze5Jam1KLT5+fnlz6EWBMH7k+jcQVBRVixCzZusgWPl85aUjT739fWZqyVa7bBGXrkk9U5msc6TMpxH6cn0v/z78POsSR/WBBPvqH0rEijIc7xpib2j9to58jiZZtsyTbyy3JLUnmQkIcuDIFCJzh0EFSU6dxBUlLbZ3FbecmmvaFFd1nGWi4vbco3Y3J7c0LJ+C69d6l0PTB7Ly6wINf69urq6Cm3xjhlYkXhliSHKjrNcg1bkoBWZqJVZEWrW0lYS7X3xJMasVyaTlFjutHrtBOKXOwgqi7tzE1EHEf2ciB7J9rcR0VNENEJE3yei7np1BEHQPhqR5bdhcQHAddn+nwP4ekrpfiL6awC3APi2drKUETzyTEpBTdbJ47y5tHh9lvT2SJ0yLAlZ5rLo6upySzVZv/yefNkgfp6U77wOfj+sZAoSywUlJbbWXks283tl5Ym3XFXeZA0a1kQU2S4r4Yi3zNq2kkG0TJYT0WYA/xrAXdk+AbgKwIPZIXsAXOepKwiC9kCev25E9CCA/w6gH8B/BHAzgCdTStuz8i0AfphSuig/Z3x8vFbxyMhIa1sdBAEAYHh4uLY9MDBQkDF1ZTkR/SGAYymlZ4nod/OPSw41/0o88MADhf1169bVtvv6+gplvb29tW0uO+WC8j09PaXHAeXSc3BwECdOnCgc1wpZziWwlX453+7t7T1j+RxtsgxQTG08OTlZKOOra2orbQLlExmuueYa7N27V43Ek3hluRVRV/bMrrjiCjzxxBOF58uPs8y2Vsry4eHhM36IvLLcms8tzQr+PCcmJmrb+bO96aabcO+995aWyfrvvltfOdtjc18J4Boi+gyAHiza3N8AsJ6IOlNK8wA2AzhsVWIlkLPsTy2hAQDVRpP72rbc90zoL8N6kPLlBhY7N++IgL0OF9+33HOWq81K1Mg7j3fGmzfE1JqdJt2cWp3yD6a1RpyGNSbjXUtOllszBZe6Fp61xHPZfhl1be6U0ldSSptTSlsB3ADgn1JKNwF4DMD12WG7ADxc92pBELSNpfi57wDwJ0T0CoCNAHR9EARB22koQi2l9GMAP862XwXwce+5fIlXoGhXS+nK973Sx5r145Xl3plDlnlgRUzlbNy4EePj46pbyMrjbs0O8kbUNYtXAls5zrQZaTMzM+q4g5Wv3lsmZblm7s3PzzeVrMF6N5t1hXmXhtaICLUgqCjRuYOgorRt4oiU5Vx2STcWlyCWW8WSPtrkk1bIcitfliWbeR1jY2NumdtMtJN3KR2ZJIHTyDJDmlkh76m2auv09LQ6Gu91u8k6eZn8jtqyRgsLC02tHNrIaLn2vshlhrx52DTilzsIKkp07iCoKNG5g6CirJjNzSO0ZFipZmdbiQosm8drc3uXpZX2jrWMj2Z7TU1NuW1a75K13mV7JZodad0riTZDzxrH4Lbu1NSUe2YZfw+sqEIrCYUWKWclhpDlzdrcWjRiI64waynqnPjlDoKKEp07CCrKislyvi9nMHGZ7s21Jss0t4qUSF7J7l1x0Uq0IGWXJoG9Mq5Z+L2Zm5tTzY9GTBjN5WeZA/J7aTK6bPJNjjVBRnOHAvbEIm+EmncpKmsikHbc7OysuuwQcGZ/KiN+uYOgokTnDoKKEp07CCpK22xuaVdzm0HaDzwLh5XMryzLSY5mN0nbyLK5vbPCLCxbtBXJ9vg+vwfWemacd999V21jI6G6zayrZo0zeBNlWM/CqkNro3wuEm18xUqm4F3CV25rZWX7ZcQvdxBUlOjcQVBR3pOynO9bstyKYtKQ8smKpvImQvAuEyRlrRbh5J1ZJuts5Nr8Wu1MbMGxZmp5Zbl1P7z1eb8zoD8nKymFJal5v5AyXCsr2y8jfrmDoKJE5w6CitI2WS5HDL2y3JsT3Fq2xmqHJec1WWetcOldTbKzs1OVg7JNvMzrFZCfWyPiVuSZVqfVRgtrlNpKv8zh+80sGSTPk8/BmwPOm57bStbApbfsEzFaHgRBKdG5g6CiROcOgorSNptb2gjclrZsbj5DzErW4LWdZTuaWXfKsu9lGbeR+bV6e3vVMQPLxrRmtVlo96e7u9u0MTlWmZaD3WqfzDHO7wdPminHWvi+lbecY42FeF2UgG5ny/fK68bi9ckEicu+nFAQBO9PXL/cRPQ6gFMAFgDMp5QuJaJBAN8HsBXA6wD+KKU0ujzNDIKgURqR5b+XUnqb7e8G8GhK6U4i2p3t36GdLGUEl0LNynJr4gjHkuWWm0nDShAgJSRvM7+WlOX8PK+Mk3VakWze1Tut6DLvSpkceT+0c3p6etQEDY08dy1JhxXJxllYWHC7uCxXlRWhpuUqb/sqnwbXAtiTbe8BcN0S6gqCoMWQJwiAiF4DMAogAfifKaW/IaKxlNJ6dsxoSmlDvj8+Pl6rWC5oHgRBaxgeHq5tDwwMFCSlV5ZfmVI6TERDAPYR0UuNNuKLX/xiYb+/v7+23dfXVyhbv772NwMDAwOl5wDFlUL5HHCgONqay7orrrgCTz/9dOE4a+RVw5LlVlkuBS+++GIcOHBAldvtkuWf/exn8dBDD6nHWZNDvHgmjtxwww24//77l1WWW16MnE9+8pP4yU9+0hJZfvr06dr2xMREoWxycrK2PTY2Vts+deoUAOCb3/wmbr311oJ5ylOBA78xXZ9//vkzvkeOq3OnlA5n/x8joh9gceneo0S0KaV0hIg2AThm1WFN7reSv/FOKx+ydzYWfynl7DTL5tZeFCsEVJbxPx68HatXr25qVptEyxfeiEtLu7Y1a8vrgrPgdXR3d6sd2noullvSmhWm/RGTCSOtpXO9Nrd857Q/1o3Y3J6c93XfKCJaQ0T9+TaAPwBwEMBeALuyw3YBeLju1YIgaBueX+5zAPwg+4vYCeBvU0r/QEQ/BfAAEd0C4NcAPrd8zQyCoFHqdu6U0qsAdpZ8/g6Aq70XkrJCW8pVHqvNoAH8rjBv3m/L9cPbaElByxXGkctz4DQjAAAJgklEQVQWWy4oS15qubMt21/KWu+MLkuKNzOjiz+z3t5eNTGHZbJ4owotV5hMwLDU/GdA8V2VElrLW27lcpN1eJ5ZRKgFQUWJzh0EFSU6dxBUlLbNCrMSzVk2Cd9uxN+p2VuWzS3r0OxlyxVm2Zjyc2sNLA3LlvauI2a57rwJF73jDrJ+LeRWhp9aWXa8trQ1jsFtZ5kw0psBxWtzW+40LcniwsKCao+X7ZcRv9xBUFGicwdBRXlPyHLpJtMmuvMZYrIOKyzTkuUyaaFWZiVqtLCiujRXmzdHulUm5bAm/7q6uszZUxwrYYUmqa2IPfn9tcQZVjJG73JC1jJM0vXlTcLgdZlZslxbUkku8exdHqpwTN0jgiB4XxKdOwgqSttkuYU1qcSbz8oaLedImSUlqtWuHEu+e0c1LZnVrOz3SmqO975JrNztVnSZJtk7OztVs6KR2Wia3LYiJGXUmTVDT4s8syaOWJLdaq8m34EYLQ+CDzTRuYOgokTnDoKK8p60uTUXhjUzRtpUmgvKcpl57RppmzdjG8m1sTierCFl12smX3hXV5eaINGy67wJK6xEC/I4LfmG1Q7LxdVMdpt6Nrc3WYNlj2tjAdaaZRGhFgRBjejcQVBR2ibLLWnpnRhgSTAph7UJFd4Eg/LaHCvnm+XGklFoWpSRN9GCrLPVstyT3LBs35vjTG5rUryR6EbNBSWfuyapZ2Zm3HJeS7ogryfbr71zUpZbOQJDlgfBB5jo3EFQUaJzB0FFaZvNbc2Iso712l5WckNr9hgvs0IxeZvk2mZa6KWEt3F2dlY9T9pT3mQQjYSt5kgXlNcV5n2e1nHyWtqzthYG8CYwtHKHS9eUd/ld71peVo5xzRUrbW6JJ298/HIHQUWJzh0EFaVtslxKRstdwuWqJT+sCB7NjSNlFr+WJfs5zebR5skmZmdnm8qhZuG9p9YSSs3Kcu/yQtpzkRFZXneXJZu5vJayXJP2s7Oz7mWCrDZ6ltgF/Ik5rFx3GvHLHQQVxdW5iWg9ET1IRC8R0YtEdDkRDRLRPiIayf7fUL+mIAjahVeW/yWAf0gpXU9E3QD6APwpgEdTSncS0W4AuwHcoV7IkOUSTZZbo8hWhBrHO0mlrM0erGg7Xv/U1FRBpnOJ7o3qkniTRngjn5pJ/gDoI+JW/TJnmHfihZTbmhS3VneVx1kRjVobG1mFU5PUjaz66llu2rPK5zoAvwPgbgBIKc2mlMYAXAtgT3bYHgDX1b1aEARtg+rFqBLRJQD+BsAvsLgg4LMAbgPwZkppPTtuNKVUk+bj4+O1ikdGRlrc7CAIAGB4eLi2PTAwUJBhHt3ZCeBjAG5NKT1FRH+JRQneEF/60pcK+5YE6e3trW339PTUtmVqY14mV80sky1f/vKX8a1vfavwmSVRtVUnZTt4e/v6+tSy/LzLLrsMTz75ZMtluRX4USYnL730UjzzzDPqnOJGcpdxGpXln/jEJ/D444+3XJafPn1aPa5Mlt9+++342te+Zq4WorXRCqaR8CCoiYmJ2nb+ve677z7ceOONZq6B/NiDBw+q1/F07kMADqWUnsr2H8Ri5z5KRJtSSkeIaBOAY1Yl3tlSgL2UDMfrfrGwls/hL1WzucS1a50+fbrwkvI/TvJ+WC4RzYb1Lt9k2boSr2vMyrettXF2dla1YS13lDVDz9tJrXzh1h9Ja7kfCy3/u3c2nfd6dW3ulNJbAN4got/OProaixJ9L4Bd2We7ADxc92pBELQN73DwrQDuzUbKXwXwb7H4h+EBIroFwK8BfG55mhgEQTO4OndK6TkAl5YUXe2+UAsi1BpJ+KBhyRsrWUMjiQs0pCuM29lcNjYiyzWsBAHS3rQm4HCatcE52rWmp6dVG9Y7OSSvp6zMK8vn5uZctq6swzsWAthLWPHPG1mRtoyIUAuCihKdOwgqSnTuIKgoK5Yg0Tvs36xbTLNJrOQBVl70VrjCpM3N28hdYd6c4BbWLDm+PTk56fZtN2Nzy3M0N5l0DVoJGazEhN4EhtoMtIWFBTPRgnYfvbnlJTJZJd/25uVX6657RBAE70vqhp82Cw8/DYJg+ZHhp/HLHQQVJTp3EFSUZZPlQRCsLPHLHQQVZdk7NxF9moh+SUSvZBlb2gIRfYeIjhHRQfZZW1NDEdEWInosS031AhHdtkLt6CGip4lof9aOP8s+30ZET2Xt+H42d2DZIaIOIvo5ET2yUu0goteJ6AARPUdEz2SftT112HKmMFvWzk1EHQD+CsC/AnAhgBuJ6MLlvCbjHgCfFp/txmJqqGEAj6KJeekNMg/g9pTSRwBcBuCPs+/f7nbMALgqpbQTwCUAPk1ElwH4cwBfz9oxCuCWZW5Hzm0AXmT7K9WO30spXZJSyudNtPu5AL9JYbYDi8lQXmxZO/KVHpbjH4DLAfwj2/8KgK8s5zXF9bcCOMj2fwlgU7a9CcAv29WW7JoPA/jUSrYDi/nvfgbgXwJ4G0Bn2bNaxutvzl7YqwA8AoBWqB2vAzhLfNbW5wJgHYDXkI19tbodyy3LzwXwBts/lH22UpyTUjoCANn/Q+26MBFtBfBRAE+tRDsyKfwcFpNq7APw/wCMpZTy0Kd2PZtvAPhPAPLwq40r1I4E4EdE9CwR/bvss3Y/l98CcBzAdzMz5S4iWtOqdix35y6Ly/zADc8T0VoAfwfgP6SUTq5EG1JKCymlS7D4y/lxAB8pO2w520BEfwjgWErpWf5xu9uRcWVK6WNYNBn/mIh+pw3XlOQpzL6dUvoogEm00BRY7s59CMAWtr8ZwOFlvqbF0SwlFDypoVoBEXVhsWPfm1J6aKXakZMWM9f+GItjAOuJKJ9f0I5ncyWAa4jodQD3Y1Gaf2MF2oGU0uHs/2MAfoDFP3jtfi5lKcw+1qp2LHfn/imA4Ww0tBvADVhMz7RStDU1FC3O+LgbwIsppb9YwXacTUTrs+1eAL+PxYGbxwBc3652pJS+klLanFLaisV34Z9SSje1ux1EtIaI+vNtAH8A4CDa/FzScqcwa8PAxWcAvIxFG+8/L/f12HXvA3AEwBwW/0LegkX77lEAI9n/g8vchk9gUWI+D+C57N9nVqAd/wLAz7N2HATwX7LPfwvA0wBeAfC/AKxu4/P5XQCPrEQ7suvtz/69kL+X7X4u2TUvAfBM9mz+N4ANrWpHRKgFQUWJCLUgqCjRuYOgokTnDoKKEp07CCpKdO4gqCjRuYOgokTnDoKKEp07CCrK/we4YHuBGfBGvgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def mean_face(faces):\n",
    "    return faces.mean(axis=1).reshape((64, 64))\n",
    "\n",
    "plt.imshow(mean_face(faces), cmap='gray');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Loops in Python are slow, and most of the time you want to utilise the fast native code provided by Numpy without explicitly using\n",
    "for loops. To put things into perspective, we can benchmark the two different implementation with the `%time` function\n",
    "in the following way:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 59.2 ms, sys: 9 µs, total: 59.2 ms\n",
      "Wall time: 59.1 ms\n",
      "CPU times: user 222 µs, sys: 0 ns, total: 222 µs\n",
      "Wall time: 191 µs\n"
     ]
    }
   ],
   "source": [
    "# We have some HUUUGE data matrix which we want to compute its mean\n",
    "X = np.random.randn(20, 1000)\n",
    "# Benchmarking time for computing mean\n",
    "%time mean_naive(X)\n",
    "%time mean(X)\n",
    "pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 50.4 ms, sys: 8.31 ms, total: 58.7 ms\n",
      "Wall time: 52.4 ms\n",
      "CPU times: user 802 µs, sys: 65 µs, total: 867 µs\n",
      "Wall time: 715 µs\n"
     ]
    }
   ],
   "source": [
    "# Benchmarking time for computing covariance\n",
    "%time cov_naive(X)\n",
    "%time cov(X)\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, using Numpy's functions makes the code much faster! Therefore, whenever you can use something that's implemented in Numpy, be sure that you take advantage of that."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Affine Transformation of Datasets\n",
    "In this week we are also going to verify a few properties about the mean and\n",
    "covariance of affine transformation of random variables.\n",
    "\n",
    "Consider a data matrix $\\boldsymbol X$ of size $(D, N)$. We would like to know\n",
    "what is the covariance when we apply affine transformation $\\boldsymbol A\\boldsymbol x_i + \\boldsymbol b$ for each datapoint $\\boldsymbol x_i$ in $\\boldsymbol X$, i.e.,\n",
    "we would like to know what happens to the mean and covariance for the new dataset if we apply affine transformation.\n",
    "\n",
    "For this assignment, you will need to implement the `affine_mean` and `affine_covariance` in the cell below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# GRADED FUNCTION: DO NOT EDIT THIS LINE\n",
    "def affine_mean(mean, A, b):\n",
    "    \"\"\"Compute the mean after affine transformation\n",
    "    Args:\n",
    "        x: ndarray, the mean vector\n",
    "        A, b: affine transformation applied to x\n",
    "    Returns:\n",
    "        mean vector after affine transformation\n",
    "    \"\"\"\n",
    "    ### Edit the code below to compute the mean vector after affine transformation\n",
    "    affine_m = np.zeros(mean.shape) # affine_m has shape (D, 1)\n",
    "    ### Update affine_m\n",
    "    affine_m = A @ mean + b\n",
    "\n",
    "    return affine_m\n",
    "\n",
    "def affine_covariance(S, A, b):\n",
    "    \"\"\"Compute the covariance matrix after affine transformation\n",
    "    Args:\n",
    "        S: ndarray, the covariance matrix\n",
    "        A, b: affine transformation applied to each element in X\n",
    "    Returns:\n",
    "        covariance matrix after the transformation\n",
    "    \"\"\"\n",
    "    ### EDIT the code below to compute the covariance matrix after affine transformation\n",
    "    affine_cov = np.zeros(S.shape) # affine_cov has shape (D, D)\n",
    "    ### Update affine_cov\n",
    "    affine_cov = A @ S @ A.T\n",
    "\n",
    "    return affine_cov"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once the two functions above are implemented, we can verify the correctness our implementation. Assuming that we have some $\\boldsymbol A$ and $\\boldsymbol b$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "random = np.random.RandomState(42)\n",
    "A = random.randn(4,4)\n",
    "b = random.randn(4,1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we can generate some random matrix $\\boldsymbol X$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = random.randn(4,100) # D = 4, N = 100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assuming that for some dataset $\\boldsymbol X$, the mean and covariance are $\\boldsymbol m$, $\\boldsymbol S$, and for the new dataset after affine transformation $\\boldsymbol X'$, the mean and covariance are $\\boldsymbol m'$ and $\\boldsymbol S'$, then we would have the following identity:\n",
    "\n",
    "$$\\boldsymbol m' = \\text{affine}_\\text{mean}(\\boldsymbol m, \\boldsymbol A, \\boldsymbol b)$$\n",
    "\n",
    "$$\\boldsymbol S' = \\text{affine}_\\text{covariance}(\\boldsymbol S, \\boldsymbol A, \\boldsymbol b)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "X1 = (A @ X) + b  # applying affine transformation to each sample in X\n",
    "X2 = (A @ X1) + b # twice"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One very useful way to compare whether arrays are equal/similar is use the helper functions\n",
    "in `numpy.testing`.\n",
    "\n",
    "Check the Numpy [documentation](https://docs.scipy.org/doc/numpy-1.13.0/reference/routines.testing.html)\n",
    "for details. The mostly used function is `np.testing.assert_almost_equal`, which raises AssertionError if the two arrays are not almost equal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.testing.assert_almost_equal(mean(X1), affine_mean(mean(X), A, b))\n",
    "np.testing.assert_almost_equal(cov(X1),  affine_covariance(cov(X), A, b))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.testing.assert_almost_equal(mean(X2), affine_mean(mean(X1), A, b))\n",
    "np.testing.assert_almost_equal(cov(X2),  affine_covariance(cov(X1), A, b))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "coursera": {
   "course_slug": "mathematics-machine-learning-pca",
   "graded_item_id": "YoDq1",
   "launcher_item_id": "vCPZ0"
  },
  "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}