CLU04-NB02.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fitting a diagonal covariance Gaussian mixture model to text data\n",
    "\n",
    "In a previous assignment, we explored k-means clustering for a high-dimensional Wikipedia dataset. We can also model this data with a mixture of Gaussians, though with increasing dimension we run into two important issues associated with using a full covariance matrix for each component.\n",
    " * Computational cost becomes prohibitive in high dimensions: score calculations have complexity cubic in the number of dimensions M if the Gaussian has a full covariance matrix.\n",
    " * A model with many parameters require more data: observe that a full covariance matrix for an M-dimensional Gaussian will have M(M+1)/2 parameters to fit. With the number of parameters growing roughly as the square of the dimension, it may quickly become impossible to find a sufficient amount of data to make good inferences.\n",
    "\n",
    "Both of these issues are avoided if we require the covariance matrix of each component to be diagonal, as then it has only M parameters to fit and the score computation decomposes into M univariate score calculations. Recall from the lecture that the M-step for the full covariance is:\n",
    "\n",
    "\\begin{align*}\n",
    "\\hat{\\Sigma}_k &= \\frac{1}{N_k^{soft}} \\sum_{i=1}^N r_{ik} (x_i-\\hat{\\mu}_k)(x_i - \\hat{\\mu}_k)^T\n",
    "\\end{align*}\n",
    "\n",
    "Note that this is a square matrix with M rows and M columns, and the above equation implies that the (v, w) element is computed by\n",
    "\n",
    "\\begin{align*}\n",
    "\\hat{\\Sigma}_{k, v, w} &= \\frac{1}{N_k^{soft}} \\sum_{i=1}^N r_{ik} (x_{iv}-\\hat{\\mu}_{kv})(x_{iw} - \\hat{\\mu}_{kw})\n",
    "\\end{align*}\n",
    "\n",
    "When we assume that this is a diagonal matrix, then non-diagonal elements are assumed to be zero and we only need to compute each of the M elements along the diagonal independently using the following equation. \n",
    "\n",
    "\\begin{align*}\n",
    "\\hat{\\sigma}^2_{k, v} &= \\hat{\\Sigma}_{k, v, v}  \\\\\n",
    "&= \\frac{1}{N_k^{soft}} \\sum_{i=1}^N r_{ik} (x_{iv}-\\hat{\\mu}_{kv})^2\n",
    "\\end{align*}\n",
    "\n",
    "In this section, we will use an EM implementation to fit a Gaussian mixture model with **diagonal** covariances to a subset of the Wikipedia dataset. The implementation uses the above equation to compute each variance term. \n",
    "\n",
    "We'll begin by importing the dataset and coming up with a useful representation for each article. After running our algorithm on the data, we will explore the output to see whether we can give a meaningful interpretation to the fitted parameters in our model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note to Amazon EC2 users**: To conserve memory, make sure to stop all the other notebooks before running this notebook."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import necessary packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import print_function # to conform python 2.x print to python 3.x\n",
    "import turicreate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also have a Python file containing implementations for several functions that will be used during the course of this assignment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from em_utilities import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load Wikipedia data and extract TF-IDF features"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load Wikipedia data and transform each of the first 5000 document into a TF-IDF representation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "wiki = turicreate.SFrame('people_wiki.sframe/').head(5000)\n",
    "#wiki['tf_idf'] = turicreate.text_analytics.tf_idf(wiki['text'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "wiki.export_csv(\"people_wiki.csv\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using a utility we provide, we will create a sparse matrix representation of the documents. This is the same utility function you used during the previous assignment on k-means with text data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "wiki = wiki.add_row_number()\n",
    "tf_idf, map_word_to_index = sframe_to_scipy(wiki, 'tf_idf')\n",
    "map_index_to_word = dict([[map_word_to_index[i], i] for i in map_word_to_index.keys()])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As in the previous assignment, we will normalize each document's TF-IDF vector to be a unit vector. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 6.47 ms, sys: 3.21 ms, total: 9.68 ms\n",
      "Wall time: 8.46 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "tf_idf = normalize(tf_idf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can check that the length (Euclidean norm) of each row is now 1.0, as expected."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0\n",
      "1.0\n",
      "0.9999999999999997\n",
      "1.0000000000000004\n",
      "0.9999999999999999\n"
     ]
    }
   ],
   "source": [
    "for i in range(5):\n",
    "    doc = tf_idf[i]\n",
    "    print(np.linalg.norm(doc.todense()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## EM in high dimensions\n",
    "\n",
    "EM for high-dimensional data requires some special treatment:\n",
    " * E step and M step must be vectorized as much as possible, as explicit loops are dreadfully slow in Python.\n",
    " * All operations must be cast in terms of sparse matrix operations, to take advantage of computational savings enabled by sparsity of data.\n",
    " * Initially, some words may be entirely absent from a cluster, causing the M step to produce zero mean and variance for those words.  This means any data point with one of those words will have 0 probability of being assigned to that cluster since the cluster allows for no variability (0 variance) around that count being 0 (0 mean). Since there is a small chance for those words to later appear in the cluster, we instead assign a small positive variance (~1e-10). Doing so also prevents numerical overflow.\n",
    " \n",
    "We provide the complete implementation for you in the file `em_utilities.py`. For those who are interested, you can read through the code to see how the sparse matrix implementation differs from the previous assignment. \n",
    "\n",
    "You are expected to answer some quiz questions using the results of clustering."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Initializing mean parameters using k-means**\n",
    "\n",
    "Recall from the lectures that EM for Gaussian mixtures is very sensitive to the choice of initial means. With a bad initial set of means, EM may produce clusters that span a large area and are mostly overlapping. To eliminate such bad outcomes, we first produce a suitable set of initial means by using the cluster centers from running k-means.  That is, we first run k-means and then take the final set of means from the converged solution as the initial means in our EM algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/robertjackson/opt/anaconda3/lib/python3.8/site-packages/sklearn/cluster/_kmeans.py:938: FutureWarning: 'n_jobs' was deprecated in version 0.23 and will be removed in 0.25.\n",
      "  warnings.warn(\"'n_jobs' was deprecated in version 0.23 and will be\"\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 58.5 s, sys: 910 ms, total: 59.4 s\n",
      "Wall time: 8.09 s\n"
     ]
    }
   ],
   "source": [
    "%%time \n",
    "\n",
    "from sklearn.cluster import KMeans\n",
    "\n",
    "np.random.seed(5)\n",
    "num_clusters = 25\n",
    "\n",
    "# Use scikit-learn's k-means to simplify workflow\n",
    "#kmeans_model = KMeans(n_clusters=num_clusters, n_init=5, max_iter=400, random_state=1, n_jobs=-1) # uncomment to use parallelism -- may break on your installation\n",
    "kmeans_model = KMeans(n_clusters=num_clusters, n_init=5, max_iter=400, random_state=1, n_jobs=1)\n",
    "kmeans_model.fit(tf_idf)\n",
    "centroids, cluster_assignment = kmeans_model.cluster_centers_, kmeans_model.labels_\n",
    "\n",
    "means = [centroid for centroid in centroids]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Initializing cluster weights**\n",
    "\n",
    "We will initialize each cluster weight to be the proportion of documents assigned to that cluster by k-means above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 899 µs, sys: 959 µs, total: 1.86 ms\n",
      "Wall time: 1.03 ms\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time \n",
    "\n",
    "num_docs = tf_idf.shape[0]\n",
    "weights = []\n",
    "for i in range(num_clusters):\n",
    "    # Compute the number of data points assigned to cluster i:\n",
    "    num_assigned = len(cluster_assignment[cluster_assignment == i])\n",
    "    w = float(num_assigned) / num_docs\n",
    "    weights.append(w)\n",
    "    \n",
    "np.sum(weights)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Initializing covariances**\n",
    "\n",
    "To initialize our covariance parameters, we compute $\\hat{\\sigma}_{k, j}^2 = \\sum_{i=1}^{N}(x_{i,j} - \\hat{\\mu}_{k, j})^2$ for each feature $j$.  For features with really tiny variances, we assign 1e-8 instead to prevent numerical instability. We do this computation in a vectorized fashion in the following code block."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [],
   "source": [
    "covs = []\n",
    "for i in range(num_clusters):\n",
    "    member_rows = tf_idf[cluster_assignment==i]\n",
    "    cov = (member_rows.multiply(member_rows) - 2*member_rows.dot(diag(means[i]))).sum(axis=0).A1 / member_rows.shape[0] \\\n",
    "          + means[i]**2\n",
    "    cov[cov < 1e-8] = 1e-8\n",
    "    covs.append(cov)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Running EM**\n",
    "\n",
    "Now that we have initialized all of our parameters, run EM."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "out = EM_for_high_dimension(tf_idf, means, covs, weights, cov_smoothing=1e-10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[3855847476.7012835, 4844053202.46348, 4844053202.46348]"
      ]
     },
     "execution_count": 78,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "out['loglik']\n",
    "#out"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interpret clustering results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In contrast to k-means, EM is able to explicitly model clusters of varying sizes and proportions. The relative magnitude of variances in the word dimensions tell us much about the nature of the clusters.\n",
    "\n",
    "Write yourself a cluster visualizer as follows.  Examining each cluster's mean vector, list the 5 words with the largest mean values (5 most common words in the cluster). For each word, also include the associated variance parameter (diagonal element of the covariance matrix). \n",
    "\n",
    "A sample output may be:\n",
    "```\n",
    "==========================================================\n",
    "Cluster 0: Largest mean parameters in cluster \n",
    "\n",
    "Word        Mean        Variance    \n",
    "football    1.08e-01    8.64e-03\n",
    "season      5.80e-02    2.93e-03\n",
    "club        4.48e-02    1.99e-03\n",
    "league      3.94e-02    1.08e-03\n",
    "played      3.83e-02    8.45e-04\n",
    "...\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fill in the blanks\n",
    "def visualize_EM_clusters(tf_idf, means, covs, map_index_to_word):\n",
    "    print('')\n",
    "    print('==========================================================')\n",
    "    \n",
    "    num_clusters = len(means)\n",
    "    for c in range(num_clusters):\n",
    "        print('Cluster {0:d}: Largest mean parameters in cluster '.format(c))\n",
    "        print('\\n{0: <12}{1: <12}{2: <12}'.format('Word', 'Mean', 'Variance'))\n",
    "        \n",
    "        # The k'th element of sorted_word_ids should be the index of the word \n",
    "        # that has the k'th-largest value in the cluster mean. Hint: Use np.argsort().\n",
    "        result = np.argsort(means[c])[::-1][:len(means[c])] # reverse order: https://stackoverflow.com/a/16486305/1253093\n",
    "        print(\"DEBUG\", \"result\", result)\n",
    "        sorted_word_ids = result\n",
    "\n",
    "        for i in sorted_word_ids[:5]:\n",
    "            print('{0: <12}{1:<10.2e}{2:10.2e}'.format(map_index_to_word[i], \n",
    "                                                       means[c][i],\n",
    "                                                       covs[c][i]))\n",
    "        print('\\n==========================================================')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "==========================================================\n",
      "Cluster 0: Largest mean parameters in cluster \n",
      "\n",
      "Word        Mean        Variance    \n",
      "DEBUG result [60179 32023 68378 ... 64971 64972     0]\n",
      "minister    7.57e-02    7.42e-03\n",
      "election    5.89e-02    3.21e-03\n",
      "party       5.89e-02    2.61e-03\n",
      "liberal     2.93e-02    4.55e-03\n",
      "elected     2.91e-02    8.95e-04\n",
      "\n",
      "==========================================================\n",
      "Cluster 1: Largest mean parameters in cluster \n",
      "\n",
      "Word        Mean        Variance    \n",
      "DEBUG result [36016 36054 35680 ... 65313 65314     0]\n",
      "film        1.76e-01    6.07e-03\n",
      "films       5.50e-02    2.97e-03\n",
      "festival    4.66e-02    3.60e-03\n",
      "feature     3.69e-02    1.81e-03\n",
      "directed    3.39e-02    2.22e-03\n",
      "\n",
      "==========================================================\n",
      "Cluster 2: Largest mean parameters in cluster \n",
      "\n",
      "Word        Mean        Variance    \n",
      "DEBUG result [10866 62192 38439 ... 64924 64925     0]\n",
      "art         1.26e-01    6.83e-03\n",
      "museum      5.62e-02    7.27e-03\n",
      "gallery     3.65e-02    3.40e-03\n",
      "artist      3.61e-02    1.44e-03\n",
      "design      3.20e-02    4.59e-03\n",
      "\n",
      "==========================================================\n",
      "Cluster 3: Largest mean parameters in cluster \n",
      "\n",
      "Word        Mean        Variance    \n",
      "DEBUG result [13611 63113 70946 ... 66180 66181 50140]\n",
      "basketball  1.86e-01    7.78e-03\n",
      "nba         1.01e-01    1.22e-02\n",
      "points      6.25e-02    5.92e-03\n",
      "coach       5.57e-02    5.91e-03\n",
      "team        4.68e-02    1.30e-03\n",
      "\n",
      "==========================================================\n",
      "Cluster 4: Largest mean parameters in cluster \n",
      "\n",
      "Word        Mean        Variance    \n",
      "DEBUG result [44208 63881 45536 ... 66384 66385     0]\n",
      "hockey      2.45e-01    1.64e-02\n",
      "nhl         1.56e-01    1.64e-02\n",
      "ice         6.40e-02    2.97e-03\n",
      "season      5.05e-02    2.52e-03\n",
      "league      4.31e-02    1.53e-03\n",
      "\n",
      "==========================================================\n",
      "Cluster 5: Largest mean parameters in cluster \n",
      "\n",
      "Word        Mean        Variance    \n",
      "DEBUG result [76386 81222 44811 ... 65846 65847     0]\n",
      "republican  7.93e-02    5.20e-03\n",
      "senate      5.41e-02    6.28e-03\n",
      "house       4.64e-02    2.41e-03\n",
      "district    4.60e-02    2.37e-03\n",
      "democratic  4.46e-02    3.02e-03\n",
      "\n",
      "==========================================================\n",
      "Cluster 6: Largest mean parameters in cluster \n",
      "\n",
      "Word        Mean        Variance    \n",
      "DEBUG result [82004 43414 60398 ... 63423 63424     0]\n",
      "she         1.60e-01    4.65e-03\n",
      "her         1.00e-01    3.14e-03\n",
      "miss        2.22e-02    7.76e-03\n",
      "women       1.43e-02    1.36e-03\n",
      "womens      1.21e-02    1.46e-03\n",
      "\n",
      "==========================================================\n",
      "Cluster 7: Largest mean parameters in cluster \n",
      "\n",
      "Word        Mean        Variance    \n",
      "DEBUG result [20600 56014 66032 ... 65805 65806     0]\n",
      "championships7.78e-02    5.17e-03\n",
      "m           4.70e-02    7.58e-03\n",
      "olympics    4.69e-02    2.59e-03\n",
      "medal       4.28e-02    2.44e-03\n",
      "she         4.18e-02    5.99e-03\n",
      "\n",
      "==========================================================\n",
      "Cluster 8: Largest mean parameters in cluster \n",
      "\n",
      "Word        Mean        Variance    \n",
      "DEBUG result [16291 73459 89498 ... 46060 46059 50140]\n",
      "book        1.45e-02    9.38e-04\n",
      "published   1.23e-02    6.16e-04\n",
      "that        1.10e-02    1.73e-04\n",
      "novel       1.07e-02    1.43e-03\n",
      "he          1.04e-02    6.05e-05\n",
      "\n",
      "==========================================================\n",
      "Cluster 9: Largest mean parameters in cluster \n",
      "\n",
      "Word        Mean        Variance    \n",
      "DEBUG result [82004 43414  6953 ... 65911 65912     0]\n",
      "she         1.37e-01    4.25e-03\n",
      "her         8.99e-02    2.74e-03\n",
      "actress     7.65e-02    4.29e-03\n",
      "film        5.98e-02    3.44e-03\n",
      "drama       5.03e-02    6.40e-03\n",
      "\n",
      "==========================================================\n",
      "Cluster 10: Largest mean parameters in cluster \n",
      "\n",
      "Word        Mean        Variance    \n",
      "DEBUG result [83970 21103 88717 ... 65952 65953 50140]\n",
      "soccer      1.15e-01    2.86e-02\n",
      "chess       4.52e-02    1.66e-02\n",
      "team        4.13e-02    2.15e-03\n",
      "coach       3.09e-02    4.45e-03\n",
      "league      3.07e-02    2.01e-03\n",
      "\n",
      "==========================================================\n",
      "Cluster 11: Largest mean parameters in cluster \n",
      "\n",
      "Word        Mean        Variance    \n",
      "DEBUG result [72208 20483 23318 ... 64757 64759 50140]\n",
      "president   2.52e-02    1.29e-03\n",
      "chairman    2.44e-02    1.97e-03\n",
      "committee   2.34e-02    2.38e-03\n",
      "served      2.24e-02    6.99e-04\n",
      "executive   2.15e-02    1.23e-03\n",
      "\n",
      "==========================================================\n",
      "Cluster 12: Largest mean parameters in cluster \n",
      "\n",
      "Word        Mean        Variance    \n",
      "DEBUG result [62213 48716 44487 ... 65195 65196     0]\n",
      "music       7.26e-02    3.48e-03\n",
      "jazz        6.07e-02    1.14e-02\n",
      "hong        3.78e-02    9.92e-03\n",
      "kong        3.50e-02    8.64e-03\n",
      "chinese     3.12e-02    5.33e-03\n",
      "\n",
      "==========================================================\n",
      "Cluster 13: Largest mean parameters in cluster \n",
      "\n",
      "Word        Mean        Variance    \n",
      "DEBUG result [93598 44054 69786 ... 64805 64806     0]\n",
      "university  3.47e-02    8.89e-04\n",
      "history     3.38e-02    2.81e-03\n",
      "philosophy  2.86e-02    5.35e-03\n",
      "professor   2.74e-02    1.08e-03\n",
      "studies     2.41e-02    1.95e-03\n",
      "\n",
      "==========================================================\n",
      "Cluster 14: Largest mean parameters in cluster \n",
      "\n",
      "Word        Mean        Variance    \n",
      "DEBUG result [89530  6924 89011 ... 64477 64478 50140]\n",
      "theatre     4.93e-02    6.17e-03\n",
      "actor       3.56e-02    2.91e-03\n",
      "television  3.21e-02    1.67e-03\n",
      "film        2.93e-02    1.16e-03\n",
      "comedy      2.86e-02    3.91e-03\n",
      "\n",
      "==========================================================\n",
      "Cluster 15: Largest mean parameters in cluster \n",
      "\n",
      "Word        Mean        Variance    \n",
      "DEBUG result [ 8199 13088 62213 ... 63419 63420     0]\n",
      "album       6.76e-02    4.78e-03\n",
      "band        5.35e-02    4.21e-03\n",
      "music       4.18e-02    1.96e-03\n",
      "released    3.13e-02    1.11e-03\n",
      "song        2.50e-02    1.81e-03\n",
      "\n",
      "==========================================================\n",
      "Cluster 16: Largest mean parameters in cluster \n",
      "\n",
      "Word        Mean        Variance    \n",
      "DEBUG result [91150 69633 74243 ... 66023 66024     0]\n",
      "tour        1.14e-01    1.92e-02\n",
      "pga         1.08e-01    2.65e-02\n",
      "racing      8.45e-02    8.26e-03\n",
      "championship6.27e-02    4.54e-03\n",
      "formula     6.06e-02    1.31e-02\n",
      "\n",
      "==========================================================\n",
      "Cluster 17: Largest mean parameters in cluster \n",
      "\n",
      "Word        Mean        Variance    \n",
      "DEBUG result [63715 74301 82443 ... 65459 65461     0]\n",
      "news        5.76e-02    8.06e-03\n",
      "radio       5.18e-02    4.62e-03\n",
      "show        3.75e-02    2.56e-03\n",
      "bbc         3.63e-02    7.41e-03\n",
      "chef        3.27e-02    1.18e-02\n",
      "\n",
      "==========================================================\n",
      "Cluster 18: Largest mean parameters in cluster \n",
      "\n",
      "Word        Mean        Variance    \n",
      "DEBUG result [36933 99102 63828 ... 66084 66085     0]\n",
      "football    1.11e-01    5.60e-03\n",
      "yards       7.37e-02    1.72e-02\n",
      "nfl         6.98e-02    9.15e-03\n",
      "coach       6.74e-02    7.85e-03\n",
      "quarterback 4.02e-02    7.16e-03\n",
      "\n",
      "==========================================================\n",
      "Cluster 19: Largest mean parameters in cluster \n",
      "\n",
      "Word        Mean        Variance    \n",
      "DEBUG result [53366 22366 80581 ... 65298 65299 50140]\n",
      "league      5.21e-02    3.13e-03\n",
      "club        5.04e-02    2.64e-03\n",
      "season      4.77e-02    2.30e-03\n",
      "rugby       4.35e-02    8.18e-03\n",
      "cup         4.22e-02    2.46e-03\n",
      "\n",
      "==========================================================\n",
      "Cluster 20: Largest mean parameters in cluster \n",
      "\n",
      "Word        Mean        Variance    \n",
      "DEBUG result [66538 62213 87812 ... 65758 65759     0]\n",
      "orchestra   1.31e-01    1.06e-02\n",
      "music       1.23e-01    6.15e-03\n",
      "symphony    8.70e-02    1.08e-02\n",
      "conductor   8.16e-02    1.01e-02\n",
      "philharmonic4.96e-02    3.27e-03\n",
      "\n",
      "==========================================================\n",
      "Cluster 21: Largest mean parameters in cluster \n",
      "\n",
      "Word        Mean        Variance    \n",
      "DEBUG result [53197 25336 49555 ... 65302 65303     0]\n",
      "law         9.52e-02    8.35e-03\n",
      "court       6.84e-02    5.24e-03\n",
      "judge       4.59e-02    4.44e-03\n",
      "attorney    3.74e-02    4.30e-03\n",
      "district    3.72e-02    4.20e-03\n",
      "\n",
      "==========================================================\n",
      "Cluster 22: Largest mean parameters in cluster \n",
      "\n",
      "Word        Mean        Variance    \n",
      "DEBUG result [36933  7481 11999 ... 66230 66231 50140]\n",
      "football    1.21e-01    6.14e-03\n",
      "afl         9.58e-02    1.31e-02\n",
      "australian  7.91e-02    1.58e-03\n",
      "club        5.93e-02    1.76e-03\n",
      "season      5.58e-02    1.83e-03\n",
      "\n",
      "==========================================================\n",
      "Cluster 23: Largest mean parameters in cluster \n",
      "\n",
      "Word        Mean        Variance    \n",
      "DEBUG result [76451 80220 93598 ... 64311 64312 50140]\n",
      "research    5.70e-02    2.68e-03\n",
      "science     3.50e-02    2.95e-03\n",
      "university  3.34e-02    7.14e-04\n",
      "professor   3.20e-02    1.26e-03\n",
      "physics     2.61e-02    5.43e-03\n",
      "\n",
      "==========================================================\n",
      "Cluster 24: Largest mean parameters in cluster \n",
      "\n",
      "Word        Mean        Variance    \n",
      "DEBUG result [13542 53366 56580 ... 65733 65734 50140]\n",
      "baseball    1.16e-01    5.57e-03\n",
      "league      1.03e-01    3.63e-03\n",
      "major       5.09e-02    1.19e-03\n",
      "games       4.66e-02    1.93e-03\n",
      "sox         4.55e-02    6.28e-03\n",
      "\n",
      "==========================================================\n"
     ]
    }
   ],
   "source": [
    "'''By EM'''\n",
    "visualize_EM_clusters(tf_idf, out['means'], out['covs'], map_index_to_word)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Quiz Question**. Select all the topics that have a cluster in the model created above. [multiple choice]\n",
    "\n",
    "- Baseball  XXX\n",
    "- Basketball XXX\n",
    "- Soccer/Football XXX\n",
    "- Music XXX\n",
    "- Politics XXX\n",
    "- Law XXX\n",
    "- Finance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparing to random initialization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create variables for randomly initializing the EM algorithm. Complete the following code block."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "np.random.seed(5) # See the note below to see why we set seed=5.\n",
    "num_clusters = len(means)\n",
    "num_docs, num_words = tf_idf.shape\n",
    "\n",
    "random_means = []\n",
    "random_covs = []\n",
    "random_weights = []\n",
    "\n",
    "for k in range(num_clusters):\n",
    "    \n",
    "    # Create a numpy array of length num_words with random normally distributed values.\n",
    "    # Use the standard univariate normal distribution (mean 0, variance 1).\n",
    "    # YOUR CODE HERE\n",
    "    mean = np.random.normal(0, 1, num_words)\n",
    "    \n",
    "    # Create a numpy array of length num_words with random values uniformly distributed between 1 and 5.\n",
    "    # YOUR CODE HERE\n",
    "    cov = np.random.uniform(1, 5, num_words)\n",
    "\n",
    "    # Initially give each cluster equal weight.\n",
    "    # YOUR CODE HERE\n",
    "    weight = 1/num_clusters\n",
    "    \n",
    "    random_means.append(mean)\n",
    "    random_covs.append(cov)\n",
    "    random_weights.append(weight)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Quiz Question**: Try fitting EM with the random initial parameters you created above. (Use `cov_smoothing=1e-5`.) Store the result to `out_random_init`. What is the final loglikelihood that the algorithm converges to?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [],
   "source": [
    "out_random_init = EM_for_high_dimension(tf_idf, random_means, random_covs, random_weights, cov_smoothing=1e-5)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2362457265.2187605"
      ]
     },
     "execution_count": 76,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "out_random_init[\"loglik\"][-1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Quiz Question:** Is the final loglikelihood larger or smaller than the final loglikelihood we obtained above when initializing EM with the results from running k-means?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2481595937.2447195"
      ]
     },
     "execution_count": 74,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# withKMeans - randomInit, with KMeans did a lot better\n",
    "4844053202.46348 - 2362457265.2187605"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Quiz Question**: For the above model, `out_random_init`, use the `visualize_EM_clusters` method you created above. Are the clusters more or less interpretable than the ones found after initializing using k-means?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "==========================================================\n",
      "Cluster 0: Largest mean parameters in cluster \n",
      "\n",
      "Word        Mean        Variance    \n",
      "DEBUG result [82004 12276 62213 ... 65460 65461     0]\n",
      "she         2.24e-02    3.54e-03\n",
      "award       1.53e-02    1.18e-03\n",
      "music       1.45e-02    1.29e-03\n",
      "university  1.43e-02    6.32e-04\n",
      "law         1.27e-02    2.59e-03\n",
      "\n",
      "==========================================================\n",
      "Cluster 1: Largest mean parameters in cluster \n",
      "\n",
      "Word        Mean        Variance    \n",
      "DEBUG result [82004 53366 80581 ... 63425 63426 50140]\n",
      "she         1.99e-02    2.80e-03\n",
      "league      1.91e-02    2.14e-03\n",
      "season      1.59e-02    1.17e-03\n",
      "football    1.53e-02    2.15e-03\n",
      "he          1.36e-02    1.13e-04\n",
      "\n",
      "==========================================================\n",
      "Cluster 2: Largest mean parameters in cluster \n",
      "\n",
      "Word        Mean        Variance    \n",
      "DEBUG result [82004 43414 62213 ... 64727 64728     0]\n",
      "she         5.65e-02    6.39e-03\n",
      "her         5.21e-02    5.62e-03\n",
      "music       1.13e-02    9.50e-04\n",
      "de          1.08e-02    1.92e-03\n",
      "opera       1.03e-02    3.26e-03\n",
      "\n",
      "==========================================================\n",
      "Cluster 3: Largest mean parameters in cluster \n",
      "\n",
      "Word        Mean        Variance    \n",
      "DEBUG result [36016 42936 44208 ... 65506 65507     0]\n",
      "film        2.98e-02    5.79e-03\n",
      "he          1.38e-02    1.07e-04\n",
      "hockey      1.33e-02    5.33e-03\n",
      "she         1.31e-02    1.69e-03\n",
      "her         1.24e-02    1.20e-03\n",
      "\n",
      "==========================================================\n",
      "Cluster 4: Largest mean parameters in cluster \n",
      "\n",
      "Word        Mean        Variance    \n",
      "DEBUG result [82004 36016 42936 ... 64134 64136     0]\n",
      "she         1.96e-02    3.21e-03\n",
      "film        1.38e-02    1.72e-03\n",
      "he          1.36e-02    1.14e-04\n",
      "her         1.35e-02    1.09e-03\n",
      "university  1.31e-02    5.36e-04\n",
      "\n",
      "==========================================================\n",
      "Cluster 5: Largest mean parameters in cluster \n",
      "\n",
      "Word        Mean        Variance    \n",
      "DEBUG result [82004 43414 42936 ... 65093 65095     0]\n",
      "she         2.03e-02    2.65e-03\n",
      "her         1.37e-02    1.09e-03\n",
      "he          1.31e-02    1.21e-04\n",
      "law         9.93e-03    1.96e-03\n",
      "court       9.54e-03    1.47e-03\n",
      "\n",
      "==========================================================\n",
      "Cluster 6: Largest mean parameters in cluster \n",
      "\n",
      "Word        Mean        Variance    \n",
      "DEBUG result [82004 43414 42936 ... 64701 64703 50140]\n",
      "she         2.70e-02    3.62e-03\n",
      "her         1.44e-02    1.18e-03\n",
      "he          1.18e-02    1.06e-04\n",
      "served      1.11e-02    3.79e-04\n",
      "state       1.03e-02    5.08e-04\n",
      "\n",
      "==========================================================\n",
      "Cluster 7: Largest mean parameters in cluster \n",
      "\n",
      "Word        Mean        Variance    \n",
      "DEBUG result [82004 43414 62213 ... 64554 64556 50140]\n",
      "she         3.41e-02    4.33e-03\n",
      "her         2.17e-02    1.89e-03\n",
      "music       1.73e-02    2.16e-03\n",
      "album       1.52e-02    2.38e-03\n",
      "marathon    1.25e-02    5.52e-03\n",
      "\n",
      "==========================================================\n",
      "Cluster 8: Largest mean parameters in cluster \n",
      "\n",
      "Word        Mean        Variance    \n",
      "DEBUG result [53366 82004 42936 ... 64774 64776 50140]\n",
      "league      1.87e-02    2.01e-03\n",
      "she         1.73e-02    2.67e-03\n",
      "he          1.46e-02    1.17e-04\n",
      "season      1.38e-02    8.80e-04\n",
      "played      1.35e-02    6.32e-04\n",
      "\n",
      "==========================================================\n",
      "Cluster 9: Largest mean parameters in cluster \n",
      "\n",
      "Word        Mean        Variance    \n",
      "DEBUG result [82004 43414 88717 ... 64607 64608     0]\n",
      "she         3.38e-02    4.81e-03\n",
      "her         1.52e-02    1.23e-03\n",
      "team        1.46e-02    8.86e-04\n",
      "played      1.45e-02    7.99e-04\n",
      "cup         1.39e-02    1.44e-03\n",
      "\n",
      "==========================================================\n",
      "Cluster 10: Largest mean parameters in cluster \n",
      "\n",
      "Word        Mean        Variance    \n",
      "DEBUG result [62213 99506 42936 ... 65215 65216 50140]\n",
      "music       1.54e-02    1.51e-03\n",
      "york        1.49e-02    8.13e-04\n",
      "he          1.31e-02    1.01e-04\n",
      "university  1.23e-02    5.72e-04\n",
      "she         1.22e-02    2.06e-03\n",
      "\n",
      "==========================================================\n",
      "Cluster 11: Largest mean parameters in cluster \n",
      "\n",
      "Word        Mean        Variance    \n",
      "DEBUG result [82004 43414 36016 ... 64427 64428 50140]\n",
      "she         2.07e-02    2.65e-03\n",
      "her         1.61e-02    1.67e-03\n",
      "film        1.55e-02    2.19e-03\n",
      "music       1.36e-02    1.71e-03\n",
      "university  1.36e-02    3.43e-04\n",
      "\n",
      "==========================================================\n",
      "Cluster 12: Largest mean parameters in cluster \n",
      "\n",
      "Word        Mean        Variance    \n",
      "DEBUG result [82004 91150 42936 ... 63921 63922     0]\n",
      "she         2.06e-02    2.98e-03\n",
      "tour        1.46e-02    3.64e-03\n",
      "he          1.42e-02    1.22e-04\n",
      "her         1.26e-02    1.11e-03\n",
      "music       1.25e-02    1.40e-03\n",
      "\n",
      "==========================================================\n",
      "Cluster 13: Largest mean parameters in cluster \n",
      "\n",
      "Word        Mean        Variance    \n",
      "DEBUG result [82004 43414 10866 ... 65642 65643     0]\n",
      "she         2.12e-02    2.25e-03\n",
      "her         1.65e-02    1.35e-03\n",
      "art         1.60e-02    2.76e-03\n",
      "nixon       1.50e-02    9.78e-03\n",
      "music       1.44e-02    1.41e-03\n",
      "\n",
      "==========================================================\n",
      "Cluster 14: Largest mean parameters in cluster \n",
      "\n",
      "Word        Mean        Variance    \n",
      "DEBUG result [62213 36016 42936 ... 65083 65084     0]\n",
      "music       1.85e-02    2.31e-03\n",
      "film        1.69e-02    2.75e-03\n",
      "he          1.31e-02    9.45e-05\n",
      "research    1.08e-02    1.06e-03\n",
      "university  1.04e-02    2.96e-04\n",
      "\n",
      "==========================================================\n",
      "Cluster 15: Largest mean parameters in cluster \n",
      "\n",
      "Word        Mean        Variance    \n",
      "DEBUG result [82004 42936 53366 ... 63913 63914 50140]\n",
      "she         2.14e-02    3.05e-03\n",
      "he          1.35e-02    1.11e-04\n",
      "league      1.32e-02    1.34e-03\n",
      "her         1.30e-02    1.32e-03\n",
      "season      1.22e-02    1.01e-03\n",
      "\n",
      "==========================================================\n",
      "Cluster 16: Largest mean parameters in cluster \n",
      "\n",
      "Word        Mean        Variance    \n",
      "DEBUG result [82004 43414 36016 ... 59431 59432 50140]\n",
      "she         2.12e-02    2.73e-03\n",
      "her         1.65e-02    1.29e-03\n",
      "film        1.36e-02    1.26e-03\n",
      "he          1.17e-02    1.01e-04\n",
      "show        1.14e-02    9.72e-04\n",
      "\n",
      "==========================================================\n",
      "Cluster 17: Largest mean parameters in cluster \n",
      "\n",
      "Word        Mean        Variance    \n",
      "DEBUG result [36016 62213 42936 ... 64960 64961     0]\n",
      "film        1.66e-02    3.18e-03\n",
      "music       1.28e-02    1.23e-03\n",
      "he          1.18e-02    9.96e-05\n",
      "she         1.11e-02    1.28e-03\n",
      "her         1.05e-02    1.17e-03\n",
      "\n",
      "==========================================================\n",
      "Cluster 18: Largest mean parameters in cluster \n",
      "\n",
      "Word        Mean        Variance    \n",
      "DEBUG result [82004 43414 13088 ... 65151 65152     0]\n",
      "she         4.03e-02    5.17e-03\n",
      "her         2.28e-02    1.75e-03\n",
      "band        1.64e-02    2.43e-03\n",
      "music       1.24e-02    9.29e-04\n",
      "university  1.19e-02    3.84e-04\n",
      "\n",
      "==========================================================\n",
      "Cluster 19: Largest mean parameters in cluster \n",
      "\n",
      "Word        Mean        Variance    \n",
      "DEBUG result [53366 80581 42936 ... 64776 64777     0]\n",
      "league      1.51e-02    1.26e-03\n",
      "season      1.43e-02    1.04e-03\n",
      "he          1.33e-02    9.20e-05\n",
      "she         1.33e-02    1.71e-03\n",
      "club        1.11e-02    8.05e-04\n",
      "\n",
      "==========================================================\n",
      "Cluster 20: Largest mean parameters in cluster \n",
      "\n",
      "Word        Mean        Variance    \n",
      "DEBUG result [82004 36016 43414 ... 64987 64988     0]\n",
      "she         2.63e-02    3.91e-03\n",
      "film        1.82e-02    1.93e-03\n",
      "her         1.45e-02    1.43e-03\n",
      "he          1.26e-02    1.05e-04\n",
      "law         1.24e-02    2.35e-03\n",
      "\n",
      "==========================================================\n",
      "Cluster 21: Largest mean parameters in cluster \n",
      "\n",
      "Word        Mean        Variance    \n",
      "DEBUG result [82004 42936 53366 ... 64295 64297     0]\n",
      "she         1.86e-02    2.57e-03\n",
      "he          1.39e-02    1.30e-04\n",
      "league      1.37e-02    1.27e-03\n",
      "season      1.22e-02    1.09e-03\n",
      "her         1.21e-02    1.09e-03\n",
      "\n",
      "==========================================================\n",
      "Cluster 22: Largest mean parameters in cluster \n",
      "\n",
      "Word        Mean        Variance    \n",
      "DEBUG result [82004 62213 42936 ... 65083 65084     0]\n",
      "she         1.50e-02    2.21e-03\n",
      "music       1.49e-02    1.63e-03\n",
      "he          1.29e-02    1.08e-04\n",
      "party       1.12e-02    8.99e-04\n",
      "her         1.00e-02    1.08e-03\n",
      "\n",
      "==========================================================\n",
      "Cluster 23: Largest mean parameters in cluster \n",
      "\n",
      "Word        Mean        Variance    \n",
      "DEBUG result [82004 43414 62213 ... 65029 65030 50140]\n",
      "she         2.27e-02    2.67e-03\n",
      "her         1.79e-02    1.67e-03\n",
      "music       1.62e-02    1.32e-03\n",
      "band        1.20e-02    1.10e-03\n",
      "york        1.18e-02    6.68e-04\n",
      "\n",
      "==========================================================\n",
      "Cluster 24: Largest mean parameters in cluster \n",
      "\n",
      "Word        Mean        Variance    \n",
      "DEBUG result [82004 43414  8199 ... 63980 63981 50140]\n",
      "she         2.62e-02    3.89e-03\n",
      "her         1.77e-02    1.67e-03\n",
      "album       1.50e-02    2.47e-03\n",
      "he          1.28e-02    1.26e-04\n",
      "soccer      1.16e-02    4.19e-03\n",
      "\n",
      "==========================================================\n"
     ]
    }
   ],
   "source": [
    "# YOUR CODE HERE. Use visualize_EM_clusters, which will require you to pass in tf_idf and map_index_to_word.\n",
    "'''By EM'''\n",
    "visualize_EM_clusters(tf_idf, out_random_init['means'], out_random_init['covs'], map_index_to_word)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note**: Random initialization may sometimes produce a superior fit than k-means initialization. We do not claim that random initialization is always worse. However, this section does illustrate that random initialization often produces much worse clustering than k-means counterpart. This is the reason why we provide the particular random seed (`np.random.seed(5)`)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Takeaway\n",
    "\n",
    "In this assignment we were able to apply the EM algorithm to a mixture of Gaussians model of text data. This was made possible by modifying the model to assume a diagonal covariance for each cluster, and by modifying the implementation to use a sparse matrix representation. In the second part you explored the role of k-means initialization on the convergence of the model as well as the interpretability of the clusters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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": 4
}