4_Gibbs_Sampling_Tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gibbs Sampling Tutorial ###"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This tutorial walks through the intuition of running Gibbs sampling for Latent Dirichlet Allocation. The method here is not the most efficient way of performing LDA, but hopefully it provides the intution for the more efficient version, known as collapsed Gibbs sampling. A nice implementation of LDA with collapsed Gibbs sampling can be found on Stephen Hansen's github at <a>https://github.com/sekhansen/text-mining-tutorial</a>. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import re\n",
    "from itertools import chain\n",
    "import numpy as np \n",
    "from collections import Counter\n",
    "\n",
    "data = [[\"rugby\",\"football\",\"competition\",\"ball\",\"games\"],\n",
    "        [\"macro\",\"economics\",\"competition\",\"games\"],\n",
    "        [\"technology\",\"computers\",\"apple\",\"AAPL\",\"internet\"],\n",
    "        [\"football\",\"score\",\"touchdown\",\"team\"],\n",
    "        [\"keynes\",\"macro\",\"friedman\",\"policy\"],\n",
    "        [\"stocks\",\"AAPL\",\"gains\",\"analysis\"],\n",
    "        [\"playoffs\",\"games\",\"season\",\"compete\",\"ball\"],\n",
    "        [\"analysis\",\"economy\",\"economics\",\"government\"],\n",
    "        [\"apple\",\"team\",\"jobs\",\"compete\",\"computers\"]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Roadmap ###\n",
    "\n",
    "We would like to extract a document-term matrix, $\\Theta$, whose rows sum to 1 and represent the distribution over topics of a document, and a term-topic matrix $\\mathbf{\\beta}$, whose columns sum to 1 and represent the distrubtion of words of a topic k.\n",
    "\n",
    "Ultimately, we would like to find out for a word in a document, what is the probability that the word belongs to a topic k, $P(z_i=j|\\mathbf{z}_{-i})$, where $z_i$ represents the assignment of word $i$ and $z_{-i}$ is the assignment of all other words.\n",
    "\n",
    "This is where Gibbs sampling will help us. \n",
    "\n",
    "* Step 1-We'll start by randomly assigning all the words in every document to a topic, 1,2,...K. \n",
    "* Step 2-Compute two matrices: \n",
    "    * Document-Topic Counts: Each row is a document. Each columns is a topic. Element i,j is the number of words in document i assigned to topic j.\n",
    "    * Term-Topic Counts: Each row is a term in our vocabular. Each column is a topic. Element i,j is the number of times in the entire corpus (across all documents) term i is assigned to topic j.\n",
    "* Step 3-For each word in each document, compute for each k topic:\n",
    "    * $P(z_i=j|\\mathbf{z}_{-i})= \\frac{n_{-i,j}^{w_i}+\\beta}{n_{-i,j}^{corpus}+V\\beta}* \\frac{n_{-i,j}^{doc}+\\alpha}{n_{-i}^{doc}+K\\alpha}$, where\n",
    "        * $n_{-i,j}^{w_i}$: number of times word $w_i$ was assigned to topic j, excluding current word\n",
    "        * $n_{-i,j}^{corpus}$: number of words assigned to topic j across entire corpus\n",
    "        * $n_{-i,j}^{doc}$: number of words in document assigned to topic j, excluding current word\n",
    "        * $n_{-i}^{doc}$: number of words, excluding current word, in document\n",
    "        * $\\beta$ is a smoothing parameter, typically equal to 200/V\n",
    "        * $\\alpha$ is a smoothing paramets, typically 50/K.\n",
    "    * We will have K values, for this specific word in this document, of $P(z_i=j|\\mathbf{z}_{-i})$. We will normalize these probabilities so that they sum to 1, and then pick 1 out of the K topics at random. \n",
    "    * <b>Note</b> $\\frac{n_{-i,j}^{w_i}+\\beta}{n_{-i,j}^{corpus}+V\\beta}$, is how much a topic likes a particular word. So if over the corpus this word shows up in that topic quite often, it will increase the probability topic K is chosen when we pick 1 topic at random.\n",
    "    * <b>Note</b> $\\frac{n_{-i,j}^{doc}+\\alpha}{n_{-i}^{doc}+K\\alpha}$, is how much a document likes a particular topic. So if a document uses a particular topic k in all other words in that document, it will increase the probability topic k is chosen when we pick 1 topic at random.\n",
    "* Step 4- Repeat step 3 N times (this is where Gibbs sampling says we will approximate the joint distribution $P(\\mathbf{z})$, the topic assignment of all words in all documents. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Code ##\n",
    "\n",
    "Let's start off by defining the number of topics, number of unique words."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "text = pd.Series(data)\n",
    "tokens = list(set(chain(*text)))\n",
    "V = len(tokens)\n",
    "K = 3\n",
    "alpha = 50.0/K\n",
    "beta = 200.0/V"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll create the Document-Topic Count matrix and the Term-Topic Count matrix, using pandas dataframes to make it easy to see the columns/indices.\n",
    "\n",
    "As the first step suggests, we will create the Document-Topic Count matrix by initializing a matrix with random numbers, 1 through K, for each word in each document."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0    [0, 2, 2, 1, 0]\n",
      "1       [1, 0, 0, 2]\n",
      "2    [2, 0, 2, 1, 2]\n",
      "3       [1, 0, 0, 2]\n",
      "4       [2, 1, 1, 2]\n",
      "dtype: object\n",
      "   0  1  2\n",
      "0  2  1  2\n",
      "1  2  1  1\n",
      "2  1  1  3\n",
      "3  2  1  1\n",
      "4  0  2  2\n",
      "             0  1  2\n",
      "economics    2  0  0\n",
      "playoffs     0  1  0\n",
      "apple        0  0  2\n",
      "AAPL         0  1  1\n",
      "competition  1  0  1\n"
     ]
    }
   ],
   "source": [
    "doc_topic = text.apply(lambda x: np.random.randint(0,K,size=(len(x),)))\n",
    "doc_topic_counts = doc_topic.apply(lambda x: Counter(x)).apply(pd.Series)\n",
    "doc_topic_counts = doc_topic_counts.fillna(0)\n",
    "\n",
    "term_topic_count = pd.DataFrame(index=tokens,columns=range(K),\n",
    "\t\t\t\t\t\t\t\t\tdata=np.zeros((V,K)))\n",
    "for doc_ind in range(text.shape[0]):\n",
    "\tfor (topic,word) in zip(doc_topic.ix[doc_ind],text.ix[doc_ind]):\n",
    "\t\tterm_topic_count.loc[word,topic]+=1\n",
    "\n",
    "print doc_topic.head()\n",
    "print doc_topic_counts.head()\n",
    "print term_topic_count.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we will iterate through each word in each document. For each step, we will remove the word from  <i>doc_topic_counts</i>, and then from the  <i>term_topic_count</i>. \n",
    "\n",
    "Then, we will compute $P(z_i=j|\\mathbf{z}_{-i})= \\frac{n_{-i,j}^{w_i}+\\beta}{n_{-i,j}^{corpus}+V\\beta}* \\frac{n_{-i,j}^{doc}+\\alpha}{n_{-i}^{doc}+K\\alpha}$. \n",
    "\n",
    "Lastly, we will draw a topic using the multinomial $(P(z_i=Topic1|\\mathbf{z}_{-i}),P(z_i=Topic2|\\mathbf{z}_{-i}),...,P(z_i=TopicK|\\mathbf{z}_{-i}))$, and then assign that word in that document to the topic. We will update both the  <i>doc_topic_counts</i> and <i>term_topic_count</i> matrices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0    [2, 0, 1, 0, 1]\n",
      "1       [0, 0, 2, 0]\n",
      "2    [2, 2, 1, 2, 0]\n",
      "3       [2, 1, 0, 0]\n",
      "4       [2, 0, 0, 2]\n",
      "dtype: object\n",
      "   0  1  2\n",
      "0  2  2  1\n",
      "1  3  0  1\n",
      "2  1  1  3\n",
      "3  2  1  1\n",
      "4  2  0  2\n",
      "             0  1  2\n",
      "economics    2  0  0\n",
      "playoffs     1  0  0\n",
      "apple        0  2  0\n",
      "AAPL         0  0  2\n",
      "competition  0  1  1\n"
     ]
    }
   ],
   "source": [
    "\tfor doc_ind,doc in enumerate(text):\n",
    "\t\ttopics = doc_topic.ix[doc_ind]\n",
    "\n",
    "\t\tfor word_ind,word in enumerate(doc):\n",
    "\t\t\t# Remove current word from current calculations\n",
    "\t\t\tdoc_topic_counts.loc[doc_ind,topics[word_ind]]-=1\n",
    "\t\t\tterm_topic_count.loc[word,topics[word_ind]]-=1\n",
    "\n",
    "\t\t\t# Find conditional probability \n",
    "\t\t\t# Multiply how much a word likes a given topic by\n",
    "\t\t\t# \thow much a document likes that topic\n",
    "\t\t\tword_given_topic = (term_topic_count.ix[word]+beta)/\\\n",
    "\t\t\t\t\t\t\t\t(doc_topic_counts.sum()+V*beta)\n",
    "\t\t\ttopic_given_doc = (doc_topic_counts.ix[doc_ind]+alpha)/\\\n",
    "\t\t\t\t\t\t\t\t(doc_topic_counts.sum(1).ix[doc_ind]+K*alpha)\n",
    "\t\t\tweights = word_given_topic*topic_given_doc\n",
    "\t\t\tweights = weights/weights.sum()\n",
    "\n",
    "\t\t\tnew_topic = np.where(np.random.multinomial(1,weights)==1)[0][0]\n",
    "\t\t\ttopics[word_ind] = new_topic\n",
    "\n",
    "\t\t\t# Add back the removed word to appropriate topic\n",
    "\t\t\tdoc_topic_counts.loc[doc_ind,new_topic]+=1\n",
    "\t\t\tterm_topic_count.loc[word,new_topic]+=1\n",
    "\n",
    "\t\tdoc_topic.ix[doc_ind] = topics\n",
    "        \n",
    "print doc_topic.head()\n",
    "print doc_topic_counts.head()\n",
    "print term_topic_count.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can see how the topic assignments within a particular document changed. \n",
    "\n",
    "We will implement hundreds or even thousands of iterations, each time getting a better approximation of the distribution. We do this finally in the class below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#################################\n",
    "### Author: Paul Soto \t\t  ###\n",
    "### \t\tpaul.soto@upf.edu ###\n",
    "#\t\t\t\t\t\t\t\t#\n",
    "# This file is a class to #######\n",
    "# run  Gibbs                #\n",
    "# sampling for Latent Dirichlet #\n",
    "# Dirichlet Allocation  \t\t#\n",
    "#################################\n",
    "\n",
    "import pandas as pd\n",
    "from Text_Preprocessing import *\n",
    "import re\n",
    "from itertools import chain\n",
    "import numpy as np \n",
    "from collections import Counter\n",
    "\n",
    "class Gibbs():\n",
    "\t\"\"\"\n",
    "\tA class for the uncollapsed Gibbs sampling on text\n",
    "\t\"\"\"\n",
    "\tdef __init__(self, text, K=2):\n",
    "\t\t\"\"\"\n",
    "\t\ttext: Pandas series with each row a list of words in the document\n",
    "\t\tK: number of topics\n",
    "\t\t\"\"\"\n",
    "\t\tself.tokens = list(set(chain(*text.values)))\n",
    "\t\tself.V = len(self.tokens)\n",
    "\t\tself.K = K\n",
    "\t\tself.text = text\n",
    "\n",
    "\t\t### Create objects we'll need in updating parameters\n",
    "\t\tself.doc_topic = text.apply(lambda x: np.random.randint(0,self.K,size=(len(x),)))\n",
    "\t\tself.doc_topic_counts = self.doc_topic.apply(lambda x: Counter(x)).apply(pd.Series)\n",
    "\t\tself.doc_topic_counts = self.doc_topic_counts.fillna(0)\n",
    "\n",
    "\t\t# Fill missing columns (typically if K is too large)\n",
    "\t\tif list(self.doc_topic_counts.columns)!=range(self.K):\n",
    "\t\t\tneeded = [el for el in range(self.K) if el not in self.doc_topic_counts.columns]\n",
    "\t\t\tfor col in needed:\n",
    "\t\t\t\tself.doc_topic_counts[col] = 0\n",
    "\n",
    "\t\tself.term_topic_count = pd.DataFrame(index=self.tokens,columns=range(self.K),\n",
    "\t\t\t\t\t\t\t\t\t\t\tdata=np.zeros((self.V,self.K)))\n",
    "\t\tfor doc_ind in range(self.text.shape[0]):\n",
    "\t\t\tfor (topic,word) in zip(self.doc_topic.ix[doc_ind],self.text.ix[doc_ind]):\n",
    "\t\t\t\tself.term_topic_count.loc[word,topic]+=1\n",
    "\n",
    "\t\t# Set priors\n",
    "\t\tself.alpha = 50.0/self.K\n",
    "\t\tself.beta = 200.0/self.V\n",
    "\t\tself.perplexity_scores = []\n",
    "\n",
    "\tdef iterate(self,n=1000):\n",
    "\t\t\"\"\"\n",
    "\t\tRun n steps of the Gibbs sampler\n",
    "\n",
    "\t\tRelies on two calculations: \n",
    "\t\t\tword_given_topic: \"probability\" of observing a word given a topic\n",
    "\t\t\ttopic_given_doc: \"probability\" of observing topic j \n",
    "\n",
    "\t\t\tEach is calculated by removing the current word from document\n",
    "\t\t\"\"\"\n",
    "\t\tfor step in range(n):\n",
    "\t\t\tif step%25==0: \n",
    "\t\t\t\tprint \"Step %s of Gibbs Sampling Completed\" % step\n",
    "\t\t\t\tself.perplexity()\n",
    "\n",
    "\t\t\tfor doc_ind,doc in enumerate(self.text):\n",
    "\t\t\t\ttopics = self.doc_topic.ix[doc_ind]\n",
    "\n",
    "\t\t\t\tfor word_ind,word in enumerate(doc):\n",
    "\t\t\t\t\t# Remove current word from current calculations\n",
    "\t\t\t\t\tself.doc_topic_counts.loc[doc_ind,topics[word_ind]]-=1\n",
    "\t\t\t\t\tself.term_topic_count.loc[word,topics[word_ind]]-=1\n",
    "\n",
    "\t\t\t\t\t# Find conditional probability \n",
    "\t\t\t\t\t# Multiply how much a word likes a given topic by\n",
    "\t\t\t\t\t# \thow much a document likes that topic\n",
    "\t\t\t\t\tword_given_topic = (self.term_topic_count.ix[word]+self.beta)/\\\n",
    "\t\t\t\t\t\t\t\t\t\t(self.doc_topic_counts.sum()+self.V*self.beta)\n",
    "\t\t\t\t\ttopic_given_doc = (self.doc_topic_counts.ix[doc_ind]+self.alpha)/\\\n",
    "\t\t\t\t\t\t\t\t\t\t(self.doc_topic_counts.sum(1).ix[doc_ind]+self.K*self.alpha)\n",
    "\t\t\t\t\tweights = word_given_topic*topic_given_doc\n",
    "\t\t\t\t\tweights = weights/weights.sum()\n",
    "\n",
    "\t\t\t\t\tnew_topic = np.where(np.random.multinomial(1,weights)==1)[0][0]\n",
    "\t\t\t\t\ttopics[word_ind] = new_topic\n",
    "\n",
    "\t\t\t\t\t# Add back the removed word to appropriate topic\n",
    "\t\t\t\t\tself.doc_topic_counts.loc[doc_ind,new_topic]+=1\n",
    "\t\t\t\t\tself.term_topic_count.loc[word,new_topic]+=1\n",
    "\n",
    "\t\t\t\tself.doc_topic.ix[doc_ind] = topics\n",
    "\n",
    "\tdef perplexity(self):\n",
    "\t\t\"\"\"\n",
    "\t\tCompute perplexity scores of samples (currently insample)\n",
    "\t\t\"\"\"\n",
    "\t\tdt = (self.doc_topic_counts+self.alpha).apply(lambda x: x/x.sum(),1).fillna(0)\n",
    "\t\ttt = (self.term_topic_count+self.beta)/(self.term_topic_count+self.beta).sum().fillna(0)\n",
    "\n",
    "\t\tdef prob(row):\n",
    "\t\t\tword_list = row[0]\n",
    "\t\t\tindex = row['index']\n",
    "\t\t\tdoc_perp= 0\n",
    "\t\t\tfor each in word_list:\n",
    "\t\t\t\tdoc_perp+=np.log((tt.ix[each]*dt.ix[index]).sum())\n",
    "\t\t\treturn doc_perp\n",
    "\n",
    "\t\tperplexity = self.text.reset_index().apply(prob,1)\n",
    "\t\tperplexity = perplexity.sum()\n",
    "\t\tself.perplexity_scores.append(np.exp( - np.sum(perplexity) / self.text.apply(len).sum()))\n",
    "\n",
    "\tdef top_n_words(self,n=10):\n",
    "\t\t\"\"\"\n",
    "\t\tReturns the n most frequent words from each topic\n",
    "\t\t\"\"\"\n",
    "\n",
    "\t\tfor topic in range(self.K):\n",
    "\t\t\ttop_n = self.term_topic_count.sort(topic,ascending=False)[topic].head(n)\n",
    "\t\t\tprint \"Top %s terms for topic %s\" % (n,topic)\n",
    "\n",
    "\t\t\tfor word in top_n.index: print word"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Step 0 of Gibbs Sampling Completed\n",
      "Step 25 of Gibbs Sampling Completed\n",
      "Step 50 of Gibbs Sampling Completed\n",
      "Step 75 of Gibbs Sampling Completed\n",
      "Step 100 of Gibbs Sampling Completed\n",
      "Step 125 of Gibbs Sampling Completed\n",
      "Step 150 of Gibbs Sampling Completed\n",
      "Step 175 of Gibbs Sampling Completed\n",
      "Step 200 of Gibbs Sampling Completed\n",
      "Step 225 of Gibbs Sampling Completed\n",
      "Step 250 of Gibbs Sampling Completed\n",
      "Step 275 of Gibbs Sampling Completed\n",
      "Step 300 of Gibbs Sampling Completed\n",
      "Step 325 of Gibbs Sampling Completed\n",
      "Step 350 of Gibbs Sampling Completed\n",
      "Step 375 of Gibbs Sampling Completed\n",
      "Step 400 of Gibbs Sampling Completed\n",
      "Step 425 of Gibbs Sampling Completed\n",
      "Step 450 of Gibbs Sampling Completed\n",
      "Step 475 of Gibbs Sampling Completed\n"
     ]
    }
   ],
   "source": [
    "data = [[\"rugby\",\"football\",\"competition\",\"ball\",\"games\"],\n",
    "        [\"macro\",\"economics\",\"competition\",\"games\"],\n",
    "        [\"technology\",\"computers\",\"apple\",\"AAPL\",\"internet\"],\n",
    "        [\"football\",\"score\",\"touchdown\",\"team\"],\n",
    "        [\"keynes\",\"macro\",\"friedman\",\"policy\"],\n",
    "        [\"stocks\",\"AAPL\",\"gains\",\"analysis\"],\n",
    "        [\"playoffs\",\"games\",\"season\",\"compete\",\"ball\"],\n",
    "        [\"analysis\",\"economy\",\"economics\",\"government\"],\n",
    "        [\"apple\",\"team\",\"jobs\",\"compete\",\"computers\"]]\n",
    "\n",
    "\n",
    "gibbsobj = Gibbs(pd.Series(data),K=3)\n",
    "gibbsobj.iterate(n=500)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Top 5 terms for topic 0\n",
      "games\n",
      "touchdown\n",
      "technology\n",
      "policy\n",
      "playoffs\n",
      "Top 5 terms for topic 1\n",
      "ball\n",
      "team\n",
      "analysis\n",
      "computers\n",
      "rugby\n",
      "Top 5 terms for topic 2\n",
      "economics\n",
      "games\n",
      "macro\n",
      "apple\n",
      "AAPL\n",
      "None\n"
     ]
    }
   ],
   "source": [
    "print gibbsobj.top_n_words(5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see the topics aren't necessarily coherent, but with more iterations they should converge to easily understood topics. "
   ]
  }
 ],
 "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}