labs/lab_02/lab02_solution.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Information Retrieval in High Dimensional Data\n",
    "## Lab 2\n",
    "\n",
    "|     |     |\n",
    "| --- | --- |\n",
    "| **Name:** | Uzair Akbar |\n",
    "| **Matriculation Number:** | 03697290 |\n",
    "| **E-mail:** | [uzair.akbar@tum.de](mailto:uzair.akbar@tum.de) |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Statistical Decision Making\n",
    "### Task 1\n",
    "Consider the two-dimensional, discrete random variable $X = [X_1\\ X_2]^\\top$ subjected to the joint probability density $p_X$ as described in the following table.\n",
    "<div style=\"text-align: center\">$\\begin{array}{c|cc} p_X(X_1, X_2)  & X_2 = 0 & X_2 = 1 \\\\ \\hline\\hline X_1 = 0 & 0.4 & 0.3 \\\\ X_1 = 1 & 0.2 & 0.1\\end{array}$</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**a)** Compute the marginal probability densities $p_{X1}, p_{X2}$ and the conditional probability $P(X_2 = 0|X_1 = 0)$ as well as the expected value $\\mathbb{E}[X]$ and the covariance matrix $\\mathbb{E}[(X - \\mathbb{E}[X])(X - \\mathbb{E}[X])^\\top]$ by hand."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Marginal probability: [p(x1=0), p(x1=1)] = [0.7 0.3]\n",
      "Marginal probability: [p(x2=0), p(x2=1)] = [0.6 0.4]\n",
      "p(x2=0 | x1=0) = 0.5714285714285715\n",
      "Expected value: [0.3 0.4]\n",
      "Covariance matrix:\n",
      " [[ 0.21 -0.02]\n",
      " [-0.02  0.24]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "# joint probability table\n",
    "p_table = np.array([[0.4, 0.3], [0.2, 0.1]])\n",
    "\n",
    "# each column of X contains a possible event. 1st row corresponds to x1,\n",
    "# 2nd row corresponds to x2\n",
    "#X = [[0 1 0 1],\n",
    "#     [0 0 1 1]]\n",
    "X = np.array([[0, 1, 0, 1], \n",
    "              [0, 0, 1, 1]])\n",
    "\n",
    "# p_table.ravel('F') = [0.4  0.2  0.3  0.1] are the joint probability values that\n",
    "# correspond to these columns\n",
    "\n",
    "# marginal probabilities: sum accross rows\n",
    "p_x1 = np.sum(p_table, axis=1)\n",
    "p_x2 = np.sum(p_table, axis=0)\n",
    "\n",
    "print('Marginal probability: [p(x1=0), p(x1=1)] = {}'.format(p_x1))\n",
    "print('Marginal probability: [p(x2=0), p(x2=1)] = {}'.format(p_x2))\n",
    "\n",
    "# conditional  p(x2 = 0 | x1 = 0)\n",
    "# p(A|B) = p(A \\intersect B) / p(B)\n",
    "p_x2equals0condx1equals0 = p_table[0,0] / p_x1[0]\n",
    "print('p(x2=0 | x1=0) = {}'.format(p_x2equals0condx1equals0))\n",
    "\n",
    "# expected value\n",
    "E_X = np.dot(X, p_table.ravel('F')) # ravel('F') while keeping column order\n",
    "print('Expected value: {}'.format(E_X))\n",
    "\n",
    "# covariance matrix\n",
    "X_centered = X - np.expand_dims(E_X, axis=1) # expand_dims is needed so that numpy is able to subtract the vector from the matrix\n",
    "CovX = np.dot(np.dot(X_centered, np.diag(p_table.ravel('F'))), X_centered.T)\n",
    "print('Covariance matrix:\\n {}'.format(CovX))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**b)** Write a PYTHON function `toyrnd` that expects the positive integer parameter `n` as its input and returns a matrix `X` of size `[2,n]`, containing `n` samples drawn independently from the distribution $p_X$, as its output."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def toyrnd(n):\n",
    "    X_out = np.zeros((2,n))\n",
    "    Q = np.zeros((n,))\n",
    "    T = np.random.rand(n)\n",
    "     \n",
    "    # Interpreting [x1, x2] as binary number and Q as its decimal representation\n",
    "    Q[T>0.4] = 1\n",
    "    Q[T>0.7] = 2\n",
    "    Q[T>0.9] = 3\n",
    "    \n",
    "    X_out[0] = Q // 2 # floor division\n",
    "    X_out[1] = Q % 2 # modulus division\n",
    "    \n",
    "    return X_out"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**c)** Verify your results in a) by generating `10000` samples with `toyrnd` and computing the respective empirical values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Empirical marginal probability: [p(x1=0), p(x1=1)] = [0.69623 0.30377]\n",
      "Empirical marginal probability: [p(x2=0), p(x2=1)] = [0.5997 0.4003]\n",
      "Empirical conditional probability P(x2=0|x1=0): 0.5700989615500625\n",
      "Empirical expected value: [0.30377 0.4003 ]\n",
      "Empirical covariance matrix:\n",
      " [[ 0.21149379 -0.02060913]\n",
      " [-0.02060913  0.24005991]]\n"
     ]
    }
   ],
   "source": [
    "N = 100000\n",
    "X_observed = toyrnd(N)\n",
    "# marginal probabilities\n",
    "p_X1equals0_empirical = np.float(np.sum(X_observed[0,:]==0))/N\n",
    "p_X1equals1_empirical = np.float(np.sum(X_observed[0,:]==1))/N\n",
    "p_X2equals0_empirical = np.float(np.sum(X_observed[1,:]==0))/N\n",
    "p_X2equals1_empirical = np.float(np.sum(X_observed[1,:]==1))/N\n",
    "\n",
    "p_x1_empirical = np.array([p_X1equals0_empirical, p_X1equals1_empirical])\n",
    "p_x2_empirical = np.array([p_X2equals0_empirical, p_X2equals1_empirical])\n",
    "\n",
    "print('Empirical marginal probability: [p(x1=0), p(x1=1)] = {}'.format(p_x1_empirical))\n",
    "print('Empirical marginal probability: [p(x2=0), p(x2=1)] = {}'.format(p_x2_empirical))\n",
    "\n",
    "# conditional probability\n",
    "X2condX1 = X_observed[1, X_observed[0,:]==0]\n",
    "P_X2equals0condX1eqzals0_empirical = np.float(np.sum(X2condX1 == 0)) / len(X2condX1)\n",
    "print('Empirical conditional probability P(x2=0|x1=0):', P_X2equals0condX1eqzals0_empirical)\n",
    "# expected value\n",
    "E_X_empirical = np.sum(X_observed, axis=1)/N\n",
    "print('Empirical expected value: {}'.format(E_X_empirical))\n",
    "\n",
    "# covariance matrix\n",
    "CovX_empirical = np.dot(X_observed - np.expand_dims(E_X_empirical, axis=1), (X_observed - np.expand_dims(E_X_empirical, axis=1)).T) / N\n",
    "print('Empirical covariance matrix:\\n {}'.format(CovX_empirical))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task 2\n",
    "The MNIST training set consists of handwritten digits from 0 to 9, stored as PNG files of size $28 \\times 28$ and indexed by label. Download the provided ZIP file from Moodle and make yourself familiar with the directory structure."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**a)** Grayscale images are typically described as matrices of `uint8` values. For numerical calculations, it is more sensible to work with floating point numbers. Load two (abitrary) images from the database and convert them to matrices `I1` and `I2` of `float64` values in the interval $[0, 1]$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "First image min/max value: 0.0/255.0\n"
     ]
    }
   ],
   "source": [
    "import imageio\n",
    "\n",
    "# define to image paths which to import\n",
    "data_folder = './mnist/'\n",
    "\n",
    "impath1 = data_folder + 'd2/d2_0075.png'\n",
    "impath2 = data_folder + 'd3/d3_0013.png'\n",
    "\n",
    "# import and convert to numpy array\n",
    "I1 = np.array(imageio.imread(impath1)).astype(np.float64)\n",
    "I2 = np.array(imageio.imread(impath2)).astype(np.float64)\n",
    "\n",
    "# check values\n",
    "print('First image min/max value: {}/{}'.format(np.min(I1), np.max(I1)))\n",
    "\n",
    "# normalize values to [0,1]\n",
    "I1 = I1 / 255.\n",
    "I2 = I2 / 255."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**b)** The matrix equivalent of the euclidean norm $\\|\\cdot\\|_2$ is the Frobenius norm. For any matrix $\\mathbf{A} \\in \\mathbb{R}^{m \\times n}$, it is defined as\n",
    "\t\t\t\\begin{equation}\n",
    "\t\t\t    \\|\\mathbf{A}\\|_F = \\sqrt{\\mathrm{tr}(\\mathbf{A}^\\top \\mathbf{A})},\n",
    "\t\t\t\\end{equation}\n",
    "\t\t\twhere $\\mathrm{tr}$ denotes the trace of a matrix. Compute the distance $\\|\\mathbf{I}_1 - \\mathbf{I}_2\\|_F$ between the images `I1` and `I2` by using three different procedures in PYTHON:\t\t\t\n",
    "-  Running the `numpy.linalg.norm` function with the `'fro'` parameter\n",
    "-  Directly applying the above equation\n",
    "-  Computing the euclidean norm between the vectorized images"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Numpy Frobenius norm: 11.71589679029833\n",
      "Implemented Frobenius norm: 11.715896790298329\n",
      "Euclidean norm of vectorized images: 11.71589679029833\n"
     ]
    }
   ],
   "source": [
    "# using frobenius norm\n",
    "frob1 = np.linalg.norm(I1 - I2, 'fro')\n",
    "\n",
    "print('Numpy Frobenius norm: {}'.format(frob1))\n",
    "\n",
    "# using formula\n",
    "frob2 = np.sqrt(np.trace(np.dot(I1 - I2, (I1 - I2).T)))\n",
    "\n",
    "print('Implemented Frobenius norm: {}'.format(frob2))\n",
    "\n",
    "# using euclidean norm of vectorized images\n",
    "frob3 = np.sqrt(np.dot((I1 - I2).ravel(), (I1 - I2).ravel()))\n",
    "print('Euclidean norm of vectorized images: {}'.format(frob3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**c)** In the following, we want to solve a simple classification problem by applying *$k$-Nearest Neighbours*. To this end, choose two digit classes, e.g. $0$ and $1$, and load `n_train = 500` images from each class to the workspace. Convert them according to subtask a) and store them in vectorized form in the matrix `X_train` of size `[784, 2*n_train]`. Provide an indicator vector `Y_train` of length `2*n_train` that assigns the respective digit class label to each column of `X_train`.\n",
    "\n",
    "From each of the two classes, choose another set of `n_test=10` images and create the according matrices `X_test` and `Y_test`. Now, for each sample in the test set, determine the `k = 20` training samples with the smallest Frobenius distance to it and store their indices in the `2*n_test, k` matrix `NN`. Generate a vector `Y_kNN` containing the respective estimated class labels by performing a majority vote on `NN`. Compare the result with `Y_test`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Loading training data...\n",
      "Loading test data...\n",
      "Computing Frobenius distances...\n",
      "Determining nearest neighbors...\n",
      "Ground truth label data:\n",
      "[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]\n",
      "Labels determined by kNN:\n",
      "[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]\n"
     ]
    }
   ],
   "source": [
    "# chose which numbers to load:\n",
    "d_id1 = 2\n",
    "d_id2 = 3\n",
    "d = [d_id1, d_id2]\n",
    "\n",
    "# define # of training and testing samples\n",
    "n_train = 500\n",
    "n_test = 10\n",
    "k = 20\n",
    "\n",
    "# initialize data matrices\n",
    "X_train = np.zeros((784, 2*n_train))\n",
    "X_test = np.zeros((784, 2*n_test))\n",
    "\n",
    "print('Loading training data...')\n",
    "for j in range(0,n_train):\n",
    "    impath = data_folder + 'd' + str(d[0]) + '/d' + str(d[0]) + '_' + str(j+1).zfill(4) + '.png'\n",
    "    X_train[:,0*n_train+j] = np.array(imageio.imread(impath)).astype(np.float).ravel()/255\n",
    "    impath = data_folder + 'd' + str(d[1]) + '/d' + str(d[1]) + '_' + str(j+1).zfill(4) + '.png'\n",
    "    X_train[:,1*n_train+j] = np.array(imageio.imread(impath)).astype(np.float).ravel()/255\n",
    "        \n",
    "Y_train = np.concatenate((np.zeros(n_train), np.ones(n_train)))\n",
    "\n",
    "\n",
    "print('Loading test data...')\n",
    "for j in range(n_test):\n",
    "    impath = data_folder + 'd' + str(d[0]) + '/d' + str(d[0]) + '_' + str(n_train+j+1).zfill(4) + '.png'\n",
    "    X_test[:,0*n_test+j] = np.array(imageio.imread(impath)).astype(np.float).ravel()/255\n",
    "    impath = data_folder + 'd' + str(d[1]) + '/d' + str(d[1]) + '_' + str(n_train+j+1).zfill(4) + '.png'\n",
    "    X_test[:,1*n_test+j] = np.array(imageio.imread(impath)).astype(np.float).ravel()/255\n",
    "    \n",
    "Y_test = np.concatenate((np.zeros(n_test), np.ones(n_test)))\n",
    "\n",
    "print('Computing Frobenius distances...')\n",
    "\n",
    "D = np.zeros((2*n_test, 2*n_train))\n",
    "for i in range(2*n_test):\n",
    "    # compute norm of distance of test sample i to all training samples\n",
    "    D[i,:] = np.sqrt(np.sum((np.expand_dims(X_test[:,i], axis=1) - X_train) ** 2, axis=0))\n",
    "    \n",
    "print('Determining nearest neighbors...')\n",
    "# np.argsort outputs indices required for sorting\n",
    "NN = np.argsort(D, axis = 1)\n",
    "# we only need the k closest neighbors, hence, we cut off after k columns\n",
    "NN = NN[:,0:k]\n",
    "\n",
    "print('Ground truth label data:')\n",
    "print(Y_test.astype(np.float))\n",
    "\n",
    "# compute nearest neighbor labelling\n",
    "# sum over labels of k nearest training examples\n",
    "# and divide by k\n",
    "# if the resulting number is smaller than 0.5, we assign label 0\n",
    "# if the resulting number is greater, we assign label 1\n",
    "kNN_mask = np.sum(Y_train[NN], axis=1)/k >= 0.5\n",
    "# convert boolean to integer\n",
    "Y_kNN= kNN_mask.astype(np.float)\n",
    "\n",
    "print('Labels determined by kNN:')\n",
    "print(Y_kNN)"
   ]
  }
 ],
 "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.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}