tutorial_notebook.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook introduces how to use the topicmodels module for implementing Latent Dirichlet Allocation using the collapsed Gibbs sampling algorithm of Griffiths and Steyvers (2004).  The module contains three classes: one for processing raw text, another for implementing LDA, and another for querying.  This tutorial will go through the main features of each, for full details see the documented source code.\n",
    "\n",
    "To illustrate LDA, the tutorial uses text data from State of the Union Addresses at the paragraph level.  These are available for download from http://www.presidency.ucsb.edu/sou.php.  They are contained in the tab-separated text file speech_data_extend.txt distributed with this tutorial.\n",
    "\n",
    "To interact with this data, we begin by importing some libraries that are not strictly speaking necessary for using topicmodels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To begin, we read in the data, specifying the encoding of the text data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "data = pd.read_table(\"speech_data_extend.txt\",encoding=\"utf-8\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The data object is called a pandas DataFrame, and is similar to a Data Frame in R.  One can see the data has three fields."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index([u'president', u'speech', u'year'], dtype='object')"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.columns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the tutorial we focus on State of the Union addresses made since the television era, which began in 1947."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "9488"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = data[data.year >= 1947]\n",
    "len(data) # The number of documents (paragraphs of State of the Union Addresses) in the dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Cleaning Raw Text Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now import topicmodels, the module used in most of the analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import topicmodels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before implementing a topic model, it is important to pre-process the data.  The first class in the topicmodels module is called RawDocs and facilitates this pre-processing.  We are going to pass this class the text data contained in the DataFrame along with a list of stopwords specified by \"long\".  Stopwords are common words in English that tend to appear in all text, and so are not helpful in describing content.  There is no definitive list of stopwords, and another option we describe below is to let the data itself reveal which words are useful for discriminating among documents.  The list of words comes from http://snowball.tartarus.org/algorithms/english/stop.txt, but one need not use all of them.  In Hansen, McMahon, and Prat (2014), for example, we use just a subset of these, which you can use by specifying \"short\" instead of \"long\".  (You can view the stopwords by typing docsobj.stopwords into the interpreter)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "docsobj = topicmodels.RawDocs(data.speech, \"long\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Rather than passing the text as a DataFrame column, one can also pass a text file, in which case each new line will be read as a separate document.\n",
    "\n",
    "docsobj is now an object with several attributes.  The most important is its tokens attribute.  This is the outcome of taking each raw document, converting each contraction into its underlying words (e.g. \"don't\" into \"do not\"), coverting it into lowercase, and breaking it into its underlying linguistic elements (words, numbers, punctuation, etc.).  To illustrate, compare the fourth paragraph in the 1947 State of the Union Address as a raw document and after tokenization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "I come also to welcome you as you take up your duties and to discuss with you the manner in which you and I should fulfill our obligations to the American people during the next 2 years. \n",
      "[u'i', u'come', u'also', u'to', u'welcome', u'you', u'as', u'you', u'take', u'up', u'your', u'duties', u'and', u'to', u'discuss', u'with', u'you', u'the', u'manner', u'in', u'which', u'you', u'and', u'i', u'should', u'fulfill', u'our', u'obligations', u'to', u'the', u'american', u'people', u'during', u'the', u'next', u'2', u'years', u'.']\n"
     ]
    }
   ],
   "source": [
    "print(data.speech.values[3]) # fourth paragraph (note that Python uses 0-indexing)\n",
    "print(docsobj.tokens[3])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At this point, all tokens are in the dataset, which can be useful for some purposes.  For example, one might want to count the number of question marks in each speech.  For implementing LDA, though, one generally wishes to focus on words.  docsobj has a method to clean tokens.  This will remove all non-alphanumeric tokens, and, by default, all numeric tokens as well (to keep numeric tokens, pass False as the second argument in parentheses).  The number passed as an argument removes all tokens whose length is less that - this can be useful if some symbols in the data like copyright signs trascribed as a single 'c' that the user would like to remove, in which case one would pass 1.  In this case, we pass 1 for illustration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "docsobj.token_clean(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[u'come', u'also', u'to', u'welcome', u'you', u'as', u'you', u'take', u'up', u'your', u'duties', u'and', u'to', u'discuss', u'with', u'you', u'the', u'manner', u'in', u'which', u'you', u'and', u'should', u'fulfill', u'our', u'obligations', u'to', u'the', u'american', u'people', u'during', u'the', u'next', u'years']\n"
     ]
    }
   ],
   "source": [
    "print(docsobj.tokens[3])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we remove the stopwords we passed when creating docsobj.  Here we pass \"tokens\" as an argument to specify that the stopwords should be removed from docsobj.tokens.  The other option would be to pass \"stems\", which we discuss below.  Notice by how much we have reduced the size of the document by removing stopwords!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[u'come', u'welcome', u'duties', u'discuss', u'manner', u'fulfill', u'obligations', u'american', u'people', u'next', u'years']\n"
     ]
    }
   ],
   "source": [
    "docsobj.stopword_remove(\"tokens\")\n",
    "print(docsobj.tokens[3])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next step is to attempt to group together words that are grammatically different but themeatically identical.  For example, the document above has the token \"obligations\" but another may have \"obligation\" and yet another \"oblige.\"  Ultimately these three words denote the same concept, and so we might want them to share the same symbol.  One way to achieve this is through stemming, a process whereby words are transformed through a deterministic algorithm to a base form.  One popular stemmer is the Porter stemmer, which docsobj applies (via its implementation in Python's Natural Language Toolkit).  This creates a new stems attribute."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[u'come', u'welcom', u'duti', u'discuss', u'manner', u'fulfil', u'oblig', u'american', u'peopl', u'next', u'year']\n"
     ]
    }
   ],
   "source": [
    "docsobj.stem()\n",
    "print(docsobj.stems[3])\n",
    "docsobj.stopword_remove(\"stems\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that the outcome of stemming need not be an English word.  These stems are the data on which we will run the topic model below.  We make an additional call to remove stopwords from stems, since the stemmed forms of tokens not in the stopword list may themselves be in the stopword list.\n",
    "\n",
    "The final step in pre-processing is to drop remaining words that are not useful for identifying content.  We have already dropped standard stopwords, but there may also be data-dependent common words.  For example, in data from Supreme Court proceedings, \"justice\" might be treated as a stopword.  Also, words that appear just once or twice in the collection are not informative of content either.  Ideally, one would like a measure of informativeness that both punishes common words in the data, and rare words.  One such option is to give each stem a tf-idf (term frequency - inverse document frequency) score.  This is standard in the language processing literature, so we omit details here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "docsobj.term_rank(\"stems\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This call produces two .csv files in the working directory, df_ranking.csv and tfidf_ranking.csv.  df_ranking.csv ranks each stem according to its document frequency, or the number of documents it appears in.  tfidf_ranking.csv ranks each stem accroding to the tf-idf measure, according to which highly informative words are those that appear frequently in the entire dataset, but in relatively few documents.  Stems with the highest scores include \"gun\", \"iraq\", and \"immigr\".\n",
    "\n",
    "At this stage, the user can decide how many stems to drop based on either the df or tf-idf scores.  The first argument to rank_remove specifies the ranking method to use (\"df\" or \"tfidf\"), the second whether to drop from \"tokens\" or \"stems\" (since we formed the rankings based on stems above, we should specify stems), and finally the value to use for the cutoff for dropping stems.  One might instead prefer to provide a number $n$ such that all stems with a tf-idf value less than or equal to the $n$th ranked stem are then dropped, which we illustrate below.  One can determine the cutoff from exploring the csv files, but here we plot the ranking in Python, which indicates a reasonable cutoff is 5000.  (When using df rather than tfidf, substitute docsobj.df_ranking)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x10f78bdd0>]"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEACAYAAAC9Gb03AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFwFJREFUeJzt3XuUXWV9//H3N0QmCQkxYkhCwyU2BQUJ4S4EyqCA8HMV\nsD/rBUHQH1384a9Y7ZL7WgmKXVKXCLW1rUhpwPJb1ksJIMhFHAFZ3Cck5RJBSLg0GS7BcJVC8vz+\neHZgGCbMzJlzzj57z/u11lmz9845e39zZuYzz3n2s58dKSUkSfU1ruwCJEmtZdBLUs0Z9JJUcwa9\nJNWcQS9JNWfQS1LNDRn0EdEVEbdHRG9ELI+IhcX2aRFxXUSsiIhrI2Jq68uVJI1UDGccfURMSim9\nHBGbAb8BTgb+N/BsSunvIuJUYFpK6bTWlitJGqlhdd2klF4uFruA8UACjgIWF9sXA0c3vTpJ0qgN\nK+gjYlxE9AJrgOtTSncCM1JKfQAppTXA1q0rU5LUqOG26DeklHYHZgP7RMQu5Fb9W57W7OIkSaM3\nfiRPTik9HxE9wOFAX0TMSCn1RcRM4KnBXhMR/gGQpAaklKIZ+xnOqJv3bhxRExETgUOBB4ArgBOK\npx0PLNnUPlJKHfVYuHBh6TVUoaZOrcuarGks1NVMw2nRzwIWR8Q48h+GH6WUro6I24D/iIgvAKuA\nTza1MklSUwwZ9Cml5cAeg2xfCxzSiqIkSc0zJq+M7e7uLruEt+nEmqAz67Km4bGm4evUupplWBdM\njeoAEanVx5CkuokIUrtOxkqSqs2gl6SaM+glqeYMekmqOYNekmrOoJekmjPoJanmDHpJqjmDXpJq\nzqCXpJoz6CWp5gx6Sao5g16Sas6gl6Saa0vQ/8//tOMokqTBtCXoH3mkHUeRJA2mLUHf09OOo0iS\nBtOWoL/hhnYcRZI0mLbcSnDatMSqVTBlSksPJUm1UblbCe6/P1x4YTuOJEkaqC0t+qVLE4ccAo89\nBhMntvRwklQLlWvR77Yb7L03fO1r7TiaJKm/trToU0qsXg077QQ33QTz57f0kJJUeZVr0QPMmgXn\nnAPHHQevvNKuo0qS2taiB0gpB/2qVXnIZVdXSw8tSZVVyRY9QARccgnMng177QU//WkOf0lS67S1\nRb/R+vVwzTVw1ll5/a//Go45BjbfvKWlSFJlNLNFX0rQb5QSXHstfOtb8NBDcPTRcPDBcPjhDsOU\nNLbVJug3SgmWLcut/Ouvh97ePCTzwx+GP/9z2Hnn3O0jSWNF7YJ+oNWrc/BffTUsWQJr18KOO+bH\nzjvnMfm77gozZ8I4Z9SXVEO1D/qB1q7NXTu//S0sXw533gn33w/PPw8f/CDMm5c/Aey8M2y7LWy3\nnV0/kqptzAX9prz0Eixdmlv/vb3w8MPw+OOwZg3suSd86ENw0EGwYAFsuWVLSpCkljDoh/DCC3DL\nLXD77fDrX+dPAPPm5ZO9Rx+du4AkqZMZ9CP0hz/AddfBz38OV10FkyfDEUfklv4RR+R1SeokbQ36\niJgNXALMADYA308pfTciFgJ/CTxVPPWMlNIvBnl96UHf3/r1uZvn+uvzvDs335zD/i/+Ag44ALbZ\npuwKJan9QT8TmJlSWhoRk4G7gaOATwEvpJTOG+L1HRX0Az33HFx0EfzqVzn0d9kF9t0XDj00j+mf\nNKnsCiWNRaV23UTE5cB3gQOAF1NK3x7i+R0d9P2tXZtH8yxZkvv1770X9tsPjjwSPvpRmDOn7Aol\njRWlBX1E7AD0AB8E/gY4AVgH3AX8TUpp3SCvqUzQD9TXlydfu/LK/Hj/++Gzn4XPfCbPxilJrVJK\n0BfdNj3A11NKSyJiOvBMSilFxDnArJTS/xnkdWnhwoVvrHd3d9Pd3d2M2ttq/Xr45S9h8WL48Y/h\nc5/L/fqHHeZVu5JGr6enh56enjfWzz777PYGfUSMB64CrkkpXTDIv28PXJlSmjfIv1W2Rb8pjz8O\nP/whXHZZnojtkEPgm9808CU1T9tb9BFxCbn1/pV+22amlNYUy18G9k4pHTPIa2sX9ButXw833ghf\n/WoewnnWWXDssWVXJakO2j3qZgFwE7AcSMXjDOAYYD55yOVK4KSUUt8gr69t0G/0yiv5wqxPfSpf\nkXvMMXDiiWVXJanKvGCqQz31FPzmN/lk7T775JE6xx2X596RpJGo7B2m6m7rreHjH899+CedBCtX\n5snWNt4+UZLKYIu+xXp74fzz8y0UTz4Zpk3LI3be976yK5PUyWzRV8juu+chmcuXw/Tp8MgjMH8+\nXHppXpakVrNFX4Kf/Qy+8Y18Fe4ZZ8App0BXV9lVSeoknoytid7ePF/+QQfBn/xJHrWz775lVyWp\nE9h1UxO7757n0znsMJgwAfbfHy64AF5/vezKJNWJLfoO8vOfw1e+km+WPnNmnkN/woSyq5JUBrtu\nauzll+F3v4MTTshTJp9xRp5MTdLYYtdNjU2aBLvuCldckUP/Ax/IF1ydfHLZlUmqKlv0He7pp2H1\n6nzh1c0357tgSaq/ZrboxzdjJ2qd6dPz4/vfhwMPzHe+mjUrt/LPPhvG+ZlM0hBs0VfIE0/ALbfk\nmTLPPjvPi7/XXmVXJakVPBkrzjwTrrkGvvc9mDEjb3vve2HKlHLrktQcBr1ICb7wBbj7bnjhBdiw\nAbbYApYtg/F2yEmVZ9DrbVLKc+h8+9v5jleSqs2TsXqbCPjYx+D003MXTldXbtm/612w1VZ26Uhj\nmS36GnniiTfnwd+wAV57Ld/9atYsuOuusquTNBJ23WjY1q2Dd787B77TKUjV4ZWxGrapU/MUCvff\nX3YlkspiH/0YsNde+a5W06fn9Rkz4BOfgKOPdoSONBbYdTMGrF4NDz6Yl1OCX/4SfvIT+OIXnUNH\n6lT20WvUbrsNjj0WHn647EokDcag16i9/nqeKfO55/KFVpI6i0Gvpth339yVs/32ef0974F//uc8\nJl9SuQx6NcWTT8Ktt765/ld/BdtuCzvtBJdeauBLZTLo1RKPPgpr1sCCBdDX9+YoHUnt5xQIaok5\nc/Jj771hxQqDXqoLL5jS28ydm1v3kurBrhu9zSmnwOWX55O048blR0SeM+cHP7DvXmoH++jVUs8+\nC729eUTOhg35sX49/Nmf5Yut5s7NtzKU1DoGvUpx0kl5Fsx77oElS+DII8uuSKovg16l+ta38iRp\nF19cdiVSfTl7pUq13Xbw8stlVyFpuAx6jdjEiXl+e0nV4Dh6jdikSXDHHfDpTw/+7wcemGfGlNQZ\nhuyjj4jZwCXADGADcGFK6e8jYhrwI2B7YCXwyZTSukFebx99zbzyClx1VR6JM9Cjj+YpkO++u/11\nSXXS1pOxETETmJlSWhoRk4G7gaOAzwPPppT+LiJOBaallE4b5PUG/RiyahXstx/cfnu+SfnEiWVX\nJFVTW0/GppTWpJSWFssvAg8As8lhv7h42mLg6GYUpGqbOROmTYM99oA994Tf/77siiSN6GRsROwA\nzAduA2aklPog/zEAtm52caqeri647748/HLlynyBlaRyDTvoi26bnwBfKlr2A/tj7J/RG6ZPh09+\nEl54oexKJA1r1E1EjCeH/KUppSXF5r6ImJFS6iv68Z/a1OsXLVr0xnJ3dzfd3d0NF6zqmDzZoJeG\nq6enh56enpbse1hXxkbEJcAzKaWv9Nt2LrA2pXSuJ2M1mNNPhy23zF8ljUy7R90sAG4ClpO7ZxJw\nBnAH8B/AtsAq8vDKt516M+jHru98B846a+h70h5zDJx/fntqkqrCuW5UCRs2wDPPvPNzbr01/0H4\n9a/bU5NUFd5hSpUwbhxsPcRYrLlz4emn21OPNFY5141KNX26QS+1mkGvUm21FaxbB3/6p3DwwfCL\nX5RdkVQ/Br1KNX483HknnHMOTJ0KV19ddkVS/dhHr9Lttlv++uCDOfQlNZctenWMSZOc515qBYNe\nHWPiRO9cJbWCQa+OMWkSvPRS2VVI9WMfvTrG9OnQ2wvHHz/814wbB1//Osye3bq6pKoz6NUx5s/P\nUyG89trwX3PeebB8uUEvvRODXh1j/Pg8781ILFkCf/hDa+qR6sI+elXahAmO1JGGYtCr0iZMsEUv\nDcWgV6UZ9NLQDHpVml030tAMelXaxImweDGceCK8/nrZ1UidyRuPqNKeeAJuuSUH/apVeTZMqQ68\nw5Q0wKxZcPfdsM02ZVciNUczg96uG9WCJ2WlTTPoVQtdXfDqq2VXIXUmg161YIte2jSDXrVgi17a\nNOe6US1MmADPPpvvP9sMU6bkmTGlOjDoVQtz5458QrRNefVVOOUU+NrXmrM/qWwOr5QG+O53YcUK\n+Id/KLsSjWUOr5RaaMIE+/tVLwa9NIAjeFQ3Br00gCN4VDcGvTSALXrVjUEvDWCLXnVj0EsD2KJX\n3Rj00gBdXQa96sWglwZweKXqxqCXBrBFr7ox6KUBbNGrbgx6aYCuLnjmmTzXzeWXl12NNHpDBn1E\nXBQRfRGxrN+2hRHxRETcUzwOb22ZUvvMmAFnngmPPQbf+EbZ1UijN+SkZhFxAPAicElKaV6xbSHw\nQkrpvCEP4KRmqqilS+H44+Hee8uuRGNRWyc1SyndAjw3WB3NKEDqVJttBuvXl12FNHqj6aP/vxGx\nNCJ+EBFTm1aR1CHGj4fXXy+7Cmn0Gg367wHvSynNB9YAQ3bhSFVji1510dAdplJKT/dbvRC48p2e\nv2jRojeWu7u76e7ubuSwUlsZ9Gqnnp4eenp6WrLvYd1hKiJ2AK5MKe1arM9MKa0plr8M7J1SGvRG\nbp6MVVWtXAkHHQSrVpVdicaiZp6MHbJFHxGXAd3AVhHxGLAQODgi5gMbgJXASc0oRuoktuhVF0MG\n/SZa6he3oBapoxj0qouG+uilsWCzzfKom04bebPZZhAObtYIGPTSJmyxRZ7zZsKEsit504YN8Ld/\nC6edVnYlqhKDXtqEyZPh+efLruKtzj0X1q4tuwpVjZOaSRUyblxu1UsjYdBLFWLQqxEGvVQhBr0a\nYdBLFWLQqxEGvVQhBr0aYdBLFWLQqxEGvVQhEeDUURopg16qEFv0aoRBL1WIQa9GGPRShRj0aoRB\nL1WIQa9GGPRShRj0aoRBL1WIQa9GGPRShUQY9Bo5g16qkHHjHEevkTPopQqx60aNMOilCjHo1QiD\nXqoQg16NMOilCjHo1QiDXqoQg16NMOilCnF4pRph0EsV4vBKNcKglyrErhs1wqCXKsSgVyMMeqlC\nDHo1wqCXKsSgVyMMeqlCDHo1wqCXKsThlWqEQS9ViMMr1QiDXqoQu27UCINeqhCDXo0w6KUKMejV\nCINeqhCDXo0YMugj4qKI6IuIZf22TYuI6yJiRURcGxFTW1umJDDo1ZjhtOgvBj46YNtpwA0ppZ2A\nG4HTm12YpLdzeKUaMWTQp5RuAZ4bsPkoYHGxvBg4usl1SRqEwyvViEb76LdOKfUBpJTWAFs3ryRJ\nm2LXjRoxvkn7ecc2xqJFi95Y7u7upru7u0mHlcYWg76+enp66Onpacm+Iw3jc2BEbA9cmVKaV6w/\nAHSnlPoiYibwq5TSBzbx2jScY0ga2rJlcOyx+avqLSJIKUUz9jXcrpsoHhtdAZxQLB8PLGlGMZLe\nmS16NWI4wysvA24FdoyIxyLi88A3gUMjYgXwkWJdUosZ9GrEkH30KaVjNvFPhzS5FklDMOjVCK+M\nlSokwuGVGjmDXqoQW/RqRLOGV0pqg64uWLUK5s0ru5Kx5V/+Bfbbr+wqGjes4ZWjOoDDK6WmWrEC\nXn217CrGljlzYMqU9h6zmcMrDXpJ6kBljKOXJFWUQS9JNWfQS1LNGfSSVHMGvSTVnEEvSTVn0EtS\nzRn0klRzBr0k1ZxBL0k1Z9BLUs0Z9JJUcwa9JNWcQS9JNWfQS1LNGfSSVHMGvSTVnEEvSTVn0EtS\nzRn0klRzBr0k1ZxBL0k1Z9BLUs0Z9JJUcwa9JNWcQS9JNWfQS1LNGfSSVHMGvSTVnEEvSTU3fjQv\njoiVwDpgA/BaSmmfZhQlSWqe0bboNwDdKaXdqxTyPT09ZZfwNp1YE3RmXdY0PNY0fJ1aV7OMNuij\nCftou078pnZiTdCZdVnT8FjT8HVqXc0y2pBOwPURcWdE/GUzCpIkNdeo+uiBBSml1RExnRz4D6SU\nbmlGYZKk5oiUUnN2FLEQeCGldN6A7c05gCSNMSmlaMZ+Gm7RR8QkYFxK6cWI2AI4DDh74POaVagk\nqTGj6bqZAfxn0WIfD/x7Sum65pQlSWqWpnXdSJI6U8uGRkbE4RHxYET8NiJObdVx+h3voojoi4hl\n/bZNi4jrImJFRFwbEVP7/dvpEfFQRDwQEYf1275HRCwr6j5/FPXMjogbI+K+iFgeESeXXVOxr66I\nuD0ieou6FnZIXeMi4p6IuKIT6in2tzIi7i3eqzs6oa6ImBoRPy6OcV9E7Fvyz/mOxftzT/F1XUSc\n3AHv05cj4r+K/f17RGxedk3F/r5U/N61NxNSSk1/kP+APAxsD7wLWAq8vxXH6nfMA4D5wLJ+284F\nTimWTwW+WSzvDPSSu5x2KGrd+OnmdmDvYvlq4KMN1jMTmF8sTwZWAO8vs6Z+tU0qvm4G3AbsU3Zd\nwJeBHwJXlP2961fTI8C0AdvKfp/+Dfh8sTwemFp2Tf1qGwf8N7BtmTUB2xTfu82L9R8Bx5f9PgG7\nAMuALvLv3nXAH7ejrlF9Y9/hP/Qh4Jp+66cBp7biWAOOuz1vDfoHgRnF8kzgwcHqAa4B9i2ec3+/\n7Z8G/qlJtV0OHNJhNU0C7gL2LrMuYDZwPdDNm0Ff+vsEPApsNWBbme/TlsDvBtle+ntV7Ocw4Oay\nayIH/SpgGjkkr+iE3z3gE8CF/dbPAr4KPNDqulrVdfNHwOP91p8otrXb1imlPoCU0hpg62L7wPqe\nLLb9EbnWjZpSd0TsQP60cRv5G1pqTUU3SS+wBrg+pXRnyXV9h/wD3/+EUenvE2+9IPDEDqhrDvBM\nRFxcdJV8P/Lot054rwA+BVxWLJdWU0rpv4FvA48V+1+XUrqhzJoK/wUcWHTVTAL+F/nTT8vrqtz0\nBaPU9jPPETEZ+AnwpZTSi4PU0PaaUkobUkq7k1vS+0TELmXVFREfA/pSSkvJU2psShmjBhaklPYg\n/0J+MSIOHKSOdtY1HtgD+MeirpfIrb7Sf6Yi4l3AkcCPN1FD22qKiHcDR5E/4W8DbBERny2zJoCU\n0oPkbprryd0tvcD6wZ7a7GO3KuifBLbrtz672NZufRExAyAiZgJPFdufJP8l3WhjfZva3pCIGE8O\n+UtTSks6oab+UkrPAz3A4SXWtQA4MiIeAf4f8OGIuBRYU/b7lFJaXXx9mtz1tg/lfv+eAB5PKd1V\nrP+UHPyd8DN1BHB3SumZYr3Mmg4BHkkprU0prQf+E9i/5JoASCldnFLaK6XUDfyefO6u5XW1Kujv\nBOZGxPYRsTm5D+mKFh2rv+CtrcIrgBOK5eOBJf22f7o4Ez8HmAvcUXxsWhcR+0REAJ/r95pG/Cu5\nL+2CTqkpIt678ax+REwEDiX3EZZSV0rpjJTSdiml95F/Tm5MKR0HXFlGPRtFxKTi0xjx5gWByynx\n+1d8vH88InYsNn0EuK/Mmvr5DPkP9UZl1vQY8KGImFDs6yPA/SXXBEDk6WKIiO2Aj5O7ulpf12hP\nwLzDiYfDyX+tHgJOa9Vx+h3vMvIZ/1fJ3+jPk0/G3FDUcR3w7n7PP518FvsB4LB+2/ck/0I/BFww\ninoWkD+WLSV/RLuneE/eU1ZNxb52LWpZSh4BcGaxvdS6iv0dxJsnY8t+n+b0+94t3/gz3AF17UZu\nSC0FfkYedVN2TZOAp4Ep/baVXdPCYv/LgMXk0X+d8DN+E7mvvpc8xXtb3isvmJKkmhtrJ2Mlacwx\n6CWp5gx6Sao5g16Sas6gl6SaM+glqeYMekmqOYNekmru/wMtMk0vzCo2twAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10cc565d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot([x[1] for x in docsobj.tfidf_ranking])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "number of unique stems = 4742\n",
      "number of total stems = 250000\n"
     ]
    }
   ],
   "source": [
    "docsobj.rank_remove(\"tfidf\",\"stems\",docsobj.tfidf_ranking[5000][1])\n",
    "all_stems = [s for d in docsobj.stems for s in d]\n",
    "print(\"number of unique stems = %d\" % len(set(all_stems)))\n",
    "print(\"number of total stems = %d\" % len(all_stems))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After pre-processing, we have 4742 unique stems, and 250000 total stems.  We now proceed to estimate a topic model on them."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Estimating a Topic Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first step in estimation is to initialize a model via topicmodels' LDA class.  We will pass docsobj.stems as the set of documents, and we also need to decide on a number of topics.  Here we choose 30 topics. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ldaobj = topicmodels.LDA.LDAGibbs(docsobj.stems,30)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are three main parameters in LDA, the number of topics, and the two hyperparameters of the Dirichlet priors.  topicmodels.LDA follows the advice of Griffiths and Steyvers (2004) and sets the hyperparameter of the Dirichlet prior on topics to $200/V$, where $V$ is the number of unique vocabulary elements, and the hyperparameter of the Dirichlet prior on document-topic distributions to $50/K$, where $K$ is the number of topics."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "30\n",
      "1.66666666667\n",
      "0.0421762969211\n"
     ]
    }
   ],
   "source": [
    "print(ldaobj.K) # number of topic, user defined.\n",
    "print(ldaobj.alpha) # hyperparameter for document-topic distribution, automatically defined\n",
    "print(ldaobj.beta) # hyperparameter for topics, automatically defined"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Should users wish to define their own priors, they can do so by calling ldaobj.set_prior(alpha,beta).  \n",
    "\n",
    "Another quantity set automatically by topicmodels.LDA is a random allocation of stems to topic assignments.  It is a 250000-dimensional vector of integers in $\\{0,\\ldots,29\\}$.  Should the user wish to define another seed, call ldaobj.set_seed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[22  8 17  7 17  8  4 27 12 15]\n",
      "(250000,)\n"
     ]
    }
   ],
   "source": [
    "print(ldaobj.topic_seed[:10])\n",
    "print(ldaobj.topic_seed.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have initialized our topic model, we are ready to sample.  To sample, we pass three parameters.  The first is the number of iterations we want the chain to burn in before beginning to sample.  The second is a thinning interval, the number of iterations to let the chain run between samples.  Allowing for a thinning interval reduces autocorrelation between samples.  The third is the number of samples to take.  So, for example, if the user passes (1000,50,20) the following will happen.  First, the chain will run for 1,000 iterations.  Then 20 samples will be taken corresponding to the $\\{1050,1100\\ldots,1950,2000\\}$ iterations for a total of 2000 iterations overall.\n",
    "\n",
    "In order not to waste time in the tutorial, we start with a relatively short chain with no burnin, a thinning interval of 50, and 10 samples, for a total of 500 iterations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration 10 of (collapsed) Gibbs sampling\n",
      "Iteration 20 of (collapsed) Gibbs sampling\n",
      "Iteration 30 of (collapsed) Gibbs sampling\n",
      "Iteration 40 of (collapsed) Gibbs sampling\n",
      "Iteration 50 of (collapsed) Gibbs sampling\n",
      "Iteration 60 of (collapsed) Gibbs sampling\n",
      "Iteration 70 of (collapsed) Gibbs sampling\n",
      "Iteration 80 of (collapsed) Gibbs sampling\n",
      "Iteration 90 of (collapsed) Gibbs sampling\n",
      "Iteration 100 of (collapsed) Gibbs sampling\n",
      "Iteration 110 of (collapsed) Gibbs sampling\n",
      "Iteration 120 of (collapsed) Gibbs sampling\n",
      "Iteration 130 of (collapsed) Gibbs sampling\n",
      "Iteration 140 of (collapsed) Gibbs sampling\n",
      "Iteration 150 of (collapsed) Gibbs sampling\n",
      "Iteration 160 of (collapsed) Gibbs sampling\n",
      "Iteration 170 of (collapsed) Gibbs sampling\n",
      "Iteration 180 of (collapsed) Gibbs sampling\n",
      "Iteration 190 of (collapsed) Gibbs sampling\n",
      "Iteration 200 of (collapsed) Gibbs sampling\n",
      "Iteration 210 of (collapsed) Gibbs sampling\n",
      "Iteration 220 of (collapsed) Gibbs sampling\n",
      "Iteration 230 of (collapsed) Gibbs sampling\n",
      "Iteration 240 of (collapsed) Gibbs sampling\n",
      "Iteration 250 of (collapsed) Gibbs sampling\n",
      "Iteration 260 of (collapsed) Gibbs sampling\n",
      "Iteration 270 of (collapsed) Gibbs sampling\n",
      "Iteration 280 of (collapsed) Gibbs sampling\n",
      "Iteration 290 of (collapsed) Gibbs sampling\n",
      "Iteration 300 of (collapsed) Gibbs sampling\n",
      "Iteration 310 of (collapsed) Gibbs sampling\n",
      "Iteration 320 of (collapsed) Gibbs sampling\n",
      "Iteration 330 of (collapsed) Gibbs sampling\n",
      "Iteration 340 of (collapsed) Gibbs sampling\n",
      "Iteration 350 of (collapsed) Gibbs sampling\n",
      "Iteration 360 of (collapsed) Gibbs sampling\n",
      "Iteration 370 of (collapsed) Gibbs sampling\n",
      "Iteration 380 of (collapsed) Gibbs sampling\n",
      "Iteration 390 of (collapsed) Gibbs sampling\n",
      "Iteration 400 of (collapsed) Gibbs sampling\n",
      "Iteration 410 of (collapsed) Gibbs sampling\n",
      "Iteration 420 of (collapsed) Gibbs sampling\n",
      "Iteration 430 of (collapsed) Gibbs sampling\n",
      "Iteration 440 of (collapsed) Gibbs sampling\n",
      "Iteration 450 of (collapsed) Gibbs sampling\n",
      "Iteration 460 of (collapsed) Gibbs sampling\n",
      "Iteration 470 of (collapsed) Gibbs sampling\n",
      "Iteration 480 of (collapsed) Gibbs sampling\n",
      "Iteration 490 of (collapsed) Gibbs sampling\n",
      "Iteration 500 of (collapsed) Gibbs sampling\n"
     ]
    }
   ],
   "source": [
    "ldaobj.sample(0,50,10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because we allowed no burn in and started sampling straight away, one would imagine the initial draws were poor in terms of describing topics.  A formalization of this idea is to compute the perplexity of each of the samples.  Perplexity is a common goodness-of-fit meausure in natural language processing and information theory literature that describes how well a probability model explains data.  Lower values indicate better goodness-of-fit.  Calling ldaobj.perplexity() returns the perplexity of each sample."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([910.00942281, 874.09497579, 862.04496391, 852.09503063,\n",
       "       847.52895528, 843.21069781, 842.30736504, 840.07617437,\n",
       "       840.15844958, 838.29737893])"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ldaobj.perplexity()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Just as we suspected, the first sample has a much higher perplexity than the last.  Moreover, it might be that if we had kept sampling the chain, we could get even lower perplexity.  Once we call ldaobj.sample the first time, all further calls extend the existing chain by default rather than start from scratch.  So let's draw another ten samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration 10 of (collapsed) Gibbs sampling\n",
      "Iteration 20 of (collapsed) Gibbs sampling\n",
      "Iteration 30 of (collapsed) Gibbs sampling\n",
      "Iteration 40 of (collapsed) Gibbs sampling\n",
      "Iteration 50 of (collapsed) Gibbs sampling\n",
      "Iteration 60 of (collapsed) Gibbs sampling\n",
      "Iteration 70 of (collapsed) Gibbs sampling\n",
      "Iteration 80 of (collapsed) Gibbs sampling\n",
      "Iteration 90 of (collapsed) Gibbs sampling\n",
      "Iteration 100 of (collapsed) Gibbs sampling\n",
      "Iteration 110 of (collapsed) Gibbs sampling\n",
      "Iteration 120 of (collapsed) Gibbs sampling\n",
      "Iteration 130 of (collapsed) Gibbs sampling\n",
      "Iteration 140 of (collapsed) Gibbs sampling\n",
      "Iteration 150 of (collapsed) Gibbs sampling\n",
      "Iteration 160 of (collapsed) Gibbs sampling\n",
      "Iteration 170 of (collapsed) Gibbs sampling\n",
      "Iteration 180 of (collapsed) Gibbs sampling\n",
      "Iteration 190 of (collapsed) Gibbs sampling\n",
      "Iteration 200 of (collapsed) Gibbs sampling\n",
      "Iteration 210 of (collapsed) Gibbs sampling\n",
      "Iteration 220 of (collapsed) Gibbs sampling\n",
      "Iteration 230 of (collapsed) Gibbs sampling\n",
      "Iteration 240 of (collapsed) Gibbs sampling\n",
      "Iteration 250 of (collapsed) Gibbs sampling\n",
      "Iteration 260 of (collapsed) Gibbs sampling\n",
      "Iteration 270 of (collapsed) Gibbs sampling\n",
      "Iteration 280 of (collapsed) Gibbs sampling\n",
      "Iteration 290 of (collapsed) Gibbs sampling\n",
      "Iteration 300 of (collapsed) Gibbs sampling\n",
      "Iteration 310 of (collapsed) Gibbs sampling\n",
      "Iteration 320 of (collapsed) Gibbs sampling\n",
      "Iteration 330 of (collapsed) Gibbs sampling\n",
      "Iteration 340 of (collapsed) Gibbs sampling\n",
      "Iteration 350 of (collapsed) Gibbs sampling\n",
      "Iteration 360 of (collapsed) Gibbs sampling\n",
      "Iteration 370 of (collapsed) Gibbs sampling\n",
      "Iteration 380 of (collapsed) Gibbs sampling\n",
      "Iteration 390 of (collapsed) Gibbs sampling\n",
      "Iteration 400 of (collapsed) Gibbs sampling\n",
      "Iteration 410 of (collapsed) Gibbs sampling\n",
      "Iteration 420 of (collapsed) Gibbs sampling\n",
      "Iteration 430 of (collapsed) Gibbs sampling\n",
      "Iteration 440 of (collapsed) Gibbs sampling\n",
      "Iteration 450 of (collapsed) Gibbs sampling\n",
      "Iteration 460 of (collapsed) Gibbs sampling\n",
      "Iteration 470 of (collapsed) Gibbs sampling\n",
      "Iteration 480 of (collapsed) Gibbs sampling\n",
      "Iteration 490 of (collapsed) Gibbs sampling\n",
      "Iteration 500 of (collapsed) Gibbs sampling\n"
     ]
    }
   ],
   "source": [
    "ldaobj.sample(0,50,10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([910.00942281, 874.09497579, 862.04496391, 852.09503063,\n",
       "       847.52895528, 843.21069781, 842.30736504, 840.07617437,\n",
       "       840.15844958, 838.29737893, 837.77127923, 838.24312563,\n",
       "       835.7131597 , 834.27288412, 835.88357615, 834.14504198,\n",
       "       834.5741682 , 834.92096362, 836.03689584, 834.48549032])"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ldaobj.perplexity()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Indeed, the topic model has improved, and may yet improve more with further sampling, but at this point we will stop to continue with the tutorial.  (In research applications, one would normally apply some convergence criterion to determine the stopping point.)  Ideally we'd like to throw away the initial samples and only keep the last ones.  ldaobj.samples_keep(n) keeps the last n samples of the chain (users can also pass a list of numbers corresponding to the indices they'd like to keep - remember that Python uses 0-indexing).  We will keep the last four samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ldaobj.samples_keep(4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So far we have only been sampling topic assignments for all the words in the dataset.  What we really care about are the topics and the distribution of topics in each document.  ldaobj has been carrying these around for us while we have been sampling.  ldaobj.tt are the esitmated topics, and ldaobj.dt are the estimated document-topic distributions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(4742, 30, 4)\n",
      "(9488, 30, 4)\n"
     ]
    }
   ],
   "source": [
    "print(ldaobj.tt.shape)\n",
    "print(ldaobj.dt.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The estimated topics are represented by $4742 \\times 30$ matrices whose columns sum to one, one for each sample, while the estimated distributions of topics within each document are represented by $9488 \\times 30$ matrices whose rows sum to one.  To get an idea of the topics that have been estimated, and whether they make sense, ldaobj.topic_content(n) produces topic_description.csv in the working directory.  Its rows contain the first n stems in each topic ranked according their probability, using the final stored sample.  It's a good idea to check the topics are \"reasonable\" before proceeding with any analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ldaobj.topic_content(20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Most economics researchers will probably be most interested initially in the distributions of topics within each document.  To generate these, one should average the matrices in ldaobj.dt.  Here we have only taken four samples for purposes of illustration, but in actual research one should ideally take as many as is computationally feasible.  ldaobj has a convenience method for doing this average, which will both return it as well as, by default, write it to dt.csv in the working directory (to disable printing, pass False to the method)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dt = ldaobj.dt_avg()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One might also be interested in the average topics themselves, in which case there is a similar convenience function available that writes tt.csv to the working directory by default.  Each unique stem in the data is associated to a number corresponding to rows of the topic matrices.  Therefore, in most cases one will probably want to print out this key too, available as ldaobj.dict_print."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "tt = ldaobj.tt_avg()\n",
    "ldaobj.dict_print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One might also want to replace the speech field in the original dataset with the estimated topics in order to have a ready-to-go dataset for regression or other econometric analysis.  The following code builds this dataset, and also writes it to file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "data = data.drop('speech',1)\n",
    "for i in range(ldaobj.K): data['T' + str(i)] = dt[:,i]\n",
    "data.to_csv(\"final_output.csv\",index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If one wishes to analyze some function of the estimated document-topic distributions, this function should be computed for each separate sample and then averaged.  Since the relevant functions are context-specific, topicmodels.LDA does not provide them, but it can be easily extended to accomodate this."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Querying Using Estimated Topics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After estimating a topic model, one is often interested in estimating the distribution of topics for documents not included in estimation.  In this case, one option is to $\\textit{query}$ those documents by holding fixed the topics estimated from LDA, and only estimating the distribution of topics for the out-of-sample documents.  The topicmodels module also provides a class for such querying, which this section introduces.\n",
    "\n",
    "We will apply querying to the corpus of entire State of the Union Addresses since 1947 (recall that we estimated topics on the level of the paragraph within each speech).  In terms of estimating topics, the paragraph level is preferable to the speech level since individual paragraphs are more likely to be based around a single theme.  But, in terms of econometric work, the entire speech is a more natural unit of analysis.  At the same time, there is no general way of \"adding up\" probability distribution at the paragraph level in order to arrive at a speech-level distribution.  Hence the need for querying, which allows us to estimate the speech-level distributions.  (Extra credit: after the tutorial, estimate LDA on the entire speech level, and judge for yourself how the topics compare to those estimated at the paragraph level).\n",
    "\n",
    "The Query class is initialized in much the same way as LDA, but takes two additional objects: a 3-D array of estimated topics (number of tokens in the estimated topics $\\times$ number of estimated topics $\\times$ number of samples from the estimation); and a dictionary that maps tokens into an index.  We can just pass these directly from ldaobj, which contains data from the above estimated LDA model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "data['speech'] = [' '.join(s) for s in docsobj.stems] # replace the speech field in the original data with its cleaned version from docsobj\n",
    "aggspeeches = data.groupby(['year','president'])['speech'].apply(lambda x: ' '.join(x)) # aggregate up to the speech level\n",
    "aggdocs = topicmodels.RawDocs(aggspeeches) # create new RawDocs object that contains entire speech stems in aggdocs.tokens\n",
    "queryobj = topicmodels.LDA.QueryGibbs(aggdocs.tokens,ldaobj.token_key,ldaobj.tt) # initialize query object with ldaobj attributes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before continuing, suppose that we instead wanted to query a document whose constitutent parts had not been included in estimation, for example a State of the Union Address from the 1930s.  How to proceed?  First, create a RawDocs object with the text to be queried (recall that RawDocs can take a basic text file, which each new line treated as a separate documents).  Second, perform the same cleaning steps as were done for the documents that went into the estimated model.  However, there is no need to do any stopword removal.  When you initialize a Query object, tokens in the documents to be queried that are not present in the estimated model are automatically stripped out.\n",
    "\n",
    "Since we don't need to estimate topics when querying, we can use far fewer iterations.  Let's start with 10."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sample 0 queried\n",
      "Sample 1 queried\n",
      "Sample 2 queried\n",
      "Sample 3 queried\n"
     ]
    }
   ],
   "source": [
    "queryobj.query(10) # query our four samples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To convince yourself that we don't need many iterations, let's look at the perplexity of the data at the entire speech level.  Notice that it is much higher than the perplexity of the data at the paragraph level.  This indicates that the topic model predicts content at the paragraph level much better."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1228.92265336, 1230.21991452, 1230.41268569, 1230.958315  ])"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "queryobj.perplexity()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's triple the number of iterations to 30 and again look at the perplexity.  (Unlike LDA's sampling, each call to query starts sampling from scratch)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sample 0 queried\n",
      "Sample 1 queried\n",
      "Sample 2 queried\n",
      "Sample 3 queried\n"
     ]
    }
   ],
   "source": [
    "queryobj.query(30) # query our four samples using more iterations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1228.8451597 , 1230.29149861, 1230.47619068, 1231.00879337])"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "queryobj.perplexity()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note these values are nearly exactly the same as for the 10-iteration querying.\n",
    "\n",
    "Finally, we follow similar steps as for LDA to output the estimated distribution of topics for entire speeches."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dt_query = queryobj.dt_avg()\n",
    "aggdata = pd.DataFrame(dt_query,index=aggspeeches.index,columns=['T' + str(i) for i in range(queryobj.K)])\n",
    "aggdata.to_csv(\"final_output_agg.csv\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Assessing Output"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At this point, you can use all of the csv files this tutorial has generated with your statistical software of choice (should this not be Python!) to analyze the topics.  Before finishing, though, we can perform an initial test of whether our output makes sense intuitively.  The following code determines each President's top topics, as measured in terms of deviations from the sample average."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/stephenhansen/anaconda/lib/python2.7/site-packages/pandas/core/computation/check.py:17: UserWarning: The installed version of numexpr 2.4.4 is not supported in pandas and will be not be used\n",
      "The minimum supported version is 2.4.6\n",
      "\n",
      "  ver=ver, min_ver=_MIN_NUMEXPR_VERSION), UserWarning)\n"
     ]
    }
   ],
   "source": [
    "def top_topics(x):\n",
    "\ttop = x.values.argsort()[-5:][::-1]\n",
    "\treturn(pd.Series(top,index=range(1,6)))\n",
    "\n",
    "temp = aggdata.reset_index()\n",
    "ranking = temp.set_index('president')\n",
    "ranking = ranking - ranking.mean()\n",
    "ranking = ranking.groupby(level='president').mean()\n",
    "ranking = ranking.sort_values('year')\n",
    "ranking = ranking.drop('year',1)\n",
    "ranking = ranking.apply(top_topics,axis=1)\n",
    "ranking.to_csv(\"president_top_topics.csv\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this particular topic model, for example, George W. Bush's top topic contains words relating to military force, and Obama's employment and economic activity.  The topic model you estimate will of course vary, so I encourage you to open president_top_policy_topics.csv and topic_description.csv to have a look for yourself.  Note too that some topics probably relate to policy, while some others relate to pure rhetoric.  Depending on the nature of the analysis you want to do with the data, it may make sense to restrict attention to some subset of the estimated topics."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's all for now, I hope you enjoyed the tutorial, and begin to use topic modelling in your own work!"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}