HW2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Homework 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Import needed packages\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import nltk\n",
    "import re\n",
    "import matplotlib.pyplot as plt\n",
    "from tqdm import tqdm\n",
    "from collections import Counter\n",
    "import scipy\n",
    "import pickle\n",
    "import lda\n",
    "from sklearn.metrics import roc_auc_score\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(23155, 3)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>president</th>\n",
       "      <th>speech</th>\n",
       "      <th>year</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Washington</td>\n",
       "      <td>Fellow-Citizens of the Senate and House of Rep...</td>\n",
       "      <td>1790</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>Washington</td>\n",
       "      <td>I embrace with great satisfaction the opportun...</td>\n",
       "      <td>1790</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>Washington</td>\n",
       "      <td>In resuming your consultations for the general...</td>\n",
       "      <td>1790</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>Washington</td>\n",
       "      <td>Among the many interesting objects which will ...</td>\n",
       "      <td>1790</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>Washington</td>\n",
       "      <td>A free people ought not only to be armed, but ...</td>\n",
       "      <td>1790</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    president                                             speech  year\n",
       "0  Washington  Fellow-Citizens of the Senate and House of Rep...  1790\n",
       "1  Washington  I embrace with great satisfaction the opportun...  1790\n",
       "2  Washington  In resuming your consultations for the general...  1790\n",
       "3  Washington  Among the many interesting objects which will ...  1790\n",
       "4  Washington  A free people ought not only to be armed, but ...  1790"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Read data set\n",
    "\n",
    "data = pd.read_table(\"speech_data_extend.txt\",encoding=\"utf-8\")\n",
    "print(data.shape)\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import topicmodels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "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": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "docsobj = topicmodels.RawDocs(data.speech, \"long\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "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",
      "['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"
     ]
    }
   ],
   "source": [
    "print(data.speech.values[3]) # fourth paragraph (note that Python uses 0-indexing)\n",
    "print(docsobj.tokens[3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "docsobj.token_clean(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['come', 'also', 'to', 'welcome', 'you', 'as', 'you', 'take', 'up', 'your', 'duties', 'and', 'to', 'discuss', 'with', 'you', 'the', 'manner', 'in', 'which', 'you', 'and', 'should', 'fulfill', 'our', 'obligations', 'to', 'the', 'american', 'people', 'during', 'the', 'next', 'years']\n"
     ]
    }
   ],
   "source": [
    "print(docsobj.tokens[3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['come', 'welcome', 'duties', 'discuss', 'manner', 'fulfill', 'obligations', 'american', 'people', 'next', 'years']\n"
     ]
    }
   ],
   "source": [
    "docsobj.stopword_remove(\"tokens\")\n",
    "print(docsobj.tokens[3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['come', 'welcom', 'duti', 'discuss', 'manner', 'fulfil', 'oblig', 'american', 'peopl', 'next', 'year']\n"
     ]
    }
   ],
   "source": [
    "docsobj.stem()\n",
    "print(docsobj.stems[3])\n",
    "docsobj.stopword_remove(\"stems\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "docsobj.term_rank(\"stems\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f0c4e5f6278>]"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD8CAYAAABw1c+bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmcFOW97/HPr6dngWEZlgFZhkVElBhZMhIXXEBFJSaY\n5XogiSFGLzHqyXqTaHJzPJrNnEQTc5JXIkdJTK5b4hKNUQSXiCSKDALKKsgiDNuwjwww2+/+0YUZ\noWemmemZ6qn+vl+vfnXVU09V/7peNb+ueeqpeszdERGR7BELOwAREWlfSvwiIllGiV9EJMso8YuI\nZBklfhGRLKPELyKSZZT4RUSyjBK/iEiWUeIXEcky8eYqmFkBMA/ID+o/4u63mNlQ4CGgF7AIuMrd\nq5OsfzNwDVAHfNndn23uM3v37u1Dhgw5nu8hIpLVFi1atNPdi1Op22ziBw4DE939XTPLBeab2TPA\n14Gfu/tDZvZbEsn9Nw1XNLORwFTgA0B/4DkzO9nd65r6wCFDhlBWVpZK/CIiApjZxlTrNtvU4wnv\nBrO5wcuBicAjQfl9wBVJVp8CPOTuh919PbAWGJdqcCIikn4ptfGbWY6ZLQF2AHOBt4G97l4bVNkM\nDEiy6gBgU4P5xuqJiEg7SSnxu3udu48GBpI4Yz8l3YGY2QwzKzOzsoqKinRvXkREAsfVq8fd9wIv\nAmcBRWZ25BrBQKA8ySrlQEmD+cbq4e4z3b3U3UuLi1O6PiEiIi3QbOI3s2IzKwqmOwEXAytJ/AB8\nKqg2HXgiyepPAlPNLD/oBTQceC0dgYuISMuk0qunH3CfmeWQ+KH4k7s/ZWYrgIfM7AfAYuBeADP7\nGFDq7v/h7svN7E/ACqAWuKG5Hj0iItK2LBNH4CotLXV15xQRSZ2ZLXL30lTqRurO3V8+v4aX3tKF\nYRGRpkQq8d/90tvMU+IXEWlSpBJ/p7w4VdW6hCAi0pRIJf78eIzDtUr8IiJNiVTiz80xausy72K1\niEgmiVTij+fEqK2vDzsMEZGMFq3EHzNqdMYvItKkSCX+3JwYNXU64xcRaUrEEr/a+EVEmhOpxB/X\nGb+ISLMilfhzc4zaep3xi4g0JVKJPx6LUaszfhGRJkUq8efmqFePiEhzIpX44zG18YuINCdSib9z\nXo6e1SMi0oxIJf5unXLZf7Am7DBERDJasyNwmVkJ8AegL+DATHe/y8weBkYE1YqAvcGA7EevvwGo\nBOqA2lQHCmiJ7p1yqTxcS21dPfGcSP2miYikTSpDL9YC33D3182sK7DIzOa6+78dqWBmdwD7mtjG\nBHff2cpYm9WzMA+APVU1FHfNb+uPExHpkJo9LXb3re7+ejBdSWKg9QFHlpuZAVcCD7ZVkKnq3SWR\n7LfsPRhyJCIimeu42kPMbAgwBljQoPhcYLu7r2lkNQfmmNkiM5vRxLZnmFmZmZVVVLRsFK3+RQUA\nrN5e2aL1RUSyQcqJ38y6AI8CX3X3/Q0WTaPps/3x7j4WuAy4wczOS1bJ3We6e6m7lxYXF6ca1vuM\nGlhEt4I4C9fvbtH6IiLZIKXEb2a5JJL+/e7+WIPyOPAJ4OHG1nX38uB9B/A4MK41ATclFjPOPLEX\nL6/Zibtu5BIRSabZxB+04d8LrHT3O49afBGwyt03N7JuYXBBGDMrBCYBy1oXctPOPbmYbfsPsUBn\n/SIiSaVyxn8OcBUw0cyWBK/JwbKpHNXMY2b9zezpYLYvMN/MlgKvAX9z99lpij2pj48ZQF48xh9f\n3diWHyMi0mE1253T3ecD1siyzycp2wJMDqbXAaNaF+Lx6ZIf5+qzh3D3vHXcdGkVJT07t+fHi4hk\nvEje5fTZMwcD8MO/rQw5EhGRzBPJxF/SszPXXzCM2cu38fKalnUNFRGJqkgmfoBrxg+ld5c8/u9f\nlnFQD24TEXlPZBN/ry753DV1DBt3VXH9/YuortXjmkVEIMKJH+Cck3pz02Wn8OLqCj577wLerng3\n7JBEREIX6cQPcN35w/jBFaexaut+LvvFy9zyxDLW7zwQdlgiIqGxTLzDtbS01MvKytK6zR37D/Hj\nZ1bxtze3Ul/vXHPuUCaf1o+R/buRq0c4i0gHZ2aLUn3sfdYk/iMqKg/z/adW8OTSLQD07pLHJz80\nkFEDizhnWG+6d85tk88VEWlLSvwpqKg8zD/f3slfl27hhVU7qHfIy4lx0cg+fKB/dwb26MSp/box\nvE8XEk+tEBHJXEr8x+nA4VpWbavkL4vLmbtiO9v2H3pvWVHnXD40qAcjTujKkN6FnHJCV07u25WC\n3Jx2i09EpDlK/K1UVV3L5j0HWfLOXl7bsJvF7+xh464qausT+yonZgwrLuTUft0Y3KuQsYOKOPPE\nXvoxEJHQKPG3gdq6et7ZXcXqbZWs2LqfFVv2s3LrfrbtP0S9QzxmnD6wO2cM7clFp/Zl7KAe5MTU\nRCQi7UOJvx0dqqnjn2/v5NV1u3lt/W5WbNlPdV09vQrz+Oio/nzk9H6UDu6h6wQi0qaU+ENUeaiG\nF1bt4Nnl23huxQ6q6+rp372Aiaf24exhvbnw1D7kx9UkJCLppcSfIaqqa3lq6VaeW7mdl9fs5GBN\nHV3y43xy7AAuPLUvZw3rpXsIRCQt0pr4zawE+AOJQVUcmOnud5nZfwL/Gzjy+MvvuPvTSda/FLgL\nyAHucffbmwsqKom/oZq6ev6xdicPL9zE7OXbcIcBRZ0YPaiI84cXc9awXho7QERaLN2Jvx/Qz91f\nD4ZRXARcAVwJvOvuP2ti3RzgLeBiYDOwEJjm7iua+swoJv6Gtu8/xEtvVTB72TbeLN9HReVhAEaV\nFHHF6P5ceEpfBvXSj4CIpO54En8qI3BtBbYG05VmthIYkGIs44C1wUhcmNlDwBSgycQfdX27FXBl\naQlXlpZQU1fP4nf2smTTHh5fvIVb/7qCW/+6ghOLC7nqzMFMGzdI3URFJK2Oq4HZzIYAY4AFQdGN\nZvaGmc0ysx5JVhkAbGowv5nUfzSyQm5OjHFDezLjvGE885Vzeerfx/PF806kuraeW/+6gtG3zeGb\nf17K5j1VYYcqIhGRcuI3sy7Ao8BX3X0/8BtgGDCaxH8Ed7QmEDObYWZlZlZWUZG9o2adNqA7N08+\nlZe/NYFffXoMowYW8edFmxn/kxf5wu8X8tr63Ryq0cAyItJyzTb1AJhZLomkf7+7Pwbg7tsbLP8f\n4Kkkq5YDJQ3mBwZlx3D3mcBMSLTxpxJXlJkZl5/en8tP78+a7ZU89cZW7nl5HS+s2sFpA7px+ydO\n57QB3cMOU0Q6oGbP+C1x59G9wEp3v7NBeb8G1T4OLEuy+kJguJkNNbM8YCrwZOtCzj7D+3blaxef\nzN+/OYHrzh/GsvL9XP7f8/n1i2s1rKSIHLdUevWMB14G3gSOjF/4HWAaiWYeBzYAX3T3rWbWn0S3\nzcnB+pOBX5DozjnL3X/YXFBR79XTWmt3vMsN97/O6u2VxAxuuuwUrh1/IjE9IkIka+kGrixQVV3L\ni6squGPOatbtPEDvLvlcdtoJfPWi4fTqkh92eCLSzpT4s8juA9U8sGAjZRv38PfVFeTHY3x0VH++\nP+U0OuWpG6hItlDiz1Lz3qrgZ3NW88bmfeTHY3z2zMH0617A9LOH6NEQIhGnxJ/F3J2HF25i5rx1\nbNh1gHqHE4sLmTX9DAb26ERcPwAikaTEL+/58TMrufuldUBiNLE5XzuPPl0LQo5KRNLteBK/Tv8i\n7ubLTuV3V5/BteOHsreqhnE/fJ5HF21m6aa9YYcmIiFJ6QYu6dgmjOjDhBF9yIkZd89bxzf+vBSA\n6WcN5paPfkDdQEWyjJp6ssym3VXsrarhrufX8NzK7Qzp1ZkvXTCMDw/txZDehWGHJyItpKYeaVRJ\nz858cGB3fv2ZMXz6w4PYsKuKbz/6Jp+5ZwHVtfXNb0BEOjyd8We5ykM1/PHVjfzX7NWcckJXPja6\nP6WDezJmUJG6gIp0IGl9Hr9EW9eCXK47bxgHq+v47xfWsmr2agD+z6STuXHi8JCjE5G2oFM6IRYz\nvjFpBK9990Ieu/5senfJ42dz3mJZ+b6wQxORNqCmHjnGk0u38OUHFwMwvE8XOuflMHXcIKaNGxRy\nZCLSGDX1SKt8bFR/uhXEefrNrew7WMMbm/cxc946JX6RiFDil6QuGNGHC0b0AWDW/PXc9tQKfvT0\nSq4ZP5S8nBixmNG9U27IUYpIS6ipR5pVVV3LWT9+gX0Ha95X/tNPnc7/Ki1pZC0RaU9q6pG06pwX\n57mvn8/r7+xh696DOHDrX1dw7/z1SvwiHVCzid/MSoA/AH1JjLY1093vMrOfAh8FqoG3gavd/ZgH\nwJjZBqASqANqU/1FksxS3DWfSz5wwnvz+w/W8vPn3mLjrgMM7qU7fkU6klSGXuwH9HP3182sK7AI\nuILEwOkvuHutmf0EwN2/nWT9DUCpu+9MNSg19WS+19bv5sq7XwFgRN+uxHOMeMwoyM3he5eP1EDw\nIu0srY9scPet7v56MF0JrAQGuPscd68Nqr1K4odAssQZQ3pw25QP8JHT+zGkd2dO6FZAUec8Fqzf\nzUML3wk7PBFpwnG18ZvZEGAMsOCoRV8AHm5kNQfmmJkDd7v7zEa2PQOYATBokLoNZjoz43NnDeFz\nZw15X/kp33uGLXsPhROUiKQk5Tt3zawL8CjwVXff36D8u0AtcH8jq45397HAZcANZnZeskruPtPd\nS929tLi4OOUvIJnlgpP7sHTTXvYcqGbPgWoqD9U0v5KItKuUzvjNLJdE0r/f3R9rUP554HLgQm/k\nYoG7lwfvO8zscWAcMK+VcUuG+uDA7sxevo0x35/7XtnpA7tzzfihTBk9IMTIROSIVHr1GHAvsNLd\n72xQfinwLeB8d69qZN1CIObulcH0JOC2tEQuGemzHx5Mt0651NUlHvH8xuZ9PLa4nG898gZnD+tN\ncdf8kCMUkVTO+M8BrgLeNLMlQdl3gF8C+cDcxG8Dr7r7dWbWH7jH3SeT6AL6eLA8Djzg7rPT/B0k\ng3TvnMtVZw5+X9n1E07iojtf4okl5Vx77okhRSYiRzSb+N19PpBsbL6nG6m/BZgcTK8DRrUmQOn4\nTurThV6FeSzfsr/5yiLS5nTnrrSLYX268PjicgpycwAwg2lnDOKDA9XfX6S9KfFLu/j0uEG8s6uK\n51ZuB6Ci8jAPLHiHqWeUcPU5QxlxQteQIxTJHkr80i6uGDOAK8b8q1fPT59dxRNLtvDQwk3EYsaP\nPv7BEKMTyS4agUtC8c1LTmH+tycyqqSIjbsOhB2OSFZR4pdQlfToRPmeg2GHIZJV1NQjoRpQ1Ik5\ny7ezfMs+YmbBC+I5MQb37EwslqxDmYi0hhK/hGpwr0Kq6+r5yC/nH7PsO5NPYcZ5w0KISiTalPgl\nVJ/80AD6FRVQXVuPu1PvUO/OjQ8s5i+Lt1BVXcenPjSQgT06hx2qSGQo8Uuo8uM5TAjG9m3oxVUV\nPLZ4Myu27uePr2zknzdPJD+eE0KEItGji7uSke64chTrf/wRvnHxyew6UK0LwCJppMQvGe3k4Mau\nquq6kCMRiQ4lfslonfMSzTtK/CLpozZ+yWid8xKH6B9f3cjfV+84Znlhfpxrzx2q9n+R46DELxlt\nUM/O9CrMY/ayrccsq3eoq3fGlBRx9km9Q4hOpGNS4peMVtw1n0Xfuzjpsre2VzLp5/PYdaC6naMS\n6diabeM3sxIze9HMVpjZcjP7SlDe08zmmtma4L1HI+tPD+qsMbPp6f4Ckr16dM4D4K7n1/DlBxdz\nuFbXAURSkcrF3VrgG+4+EjiTxIDpI4GbgOfdfTjwfDD/PmbWE7gF+DCJsXZvaewHQuR49SrM44rR\n/Snfc5Anl27hhvsXhx2SSIfQbOJ3963u/nowXQmsBAYAU4D7gmr3AVckWf0SYK6773b3PcBc4NJ0\nBC4Sixm/mDqGxf+RaAratDvp0M8icpTj6s5pZkOAMcACoK+7H7nito3E+LpHGwBsajC/OSgTSZuC\n3Bw+PmYAVTW1YYci0iGknPjNrAvwKPBVd3/f4Knu7oC3JhAzm2FmZWZWVlFR0ZpNSRbqlJdD1WG1\n8YukIqXEb2a5JJL+/e7+WFC83cz6Bcv7Acd2soZyoKTB/MCg7BjuPtPdS929tLi4ONX4RQAozMvh\nQHUth2rqmnzV1tWHHapI6JrtzmlmBtwLrHT3OxssehKYDtwevD+RZPVngR81uKA7Cbi5VRGLJNG1\nIJdDNfWc8r3ZTdbr3imXl789gW4Fue0UmUjmSaUf/znAVcCbZrYkKPsOiYT/JzO7BtgIXAlgZqXA\nde5+rbvvNrPvAwuD9W5z991p/QYiwNQzSsiPx6jzxlsc125/l8cWl7Nl70G6naDEL9mr2cTv7vOB\nxoZBujBJ/TLg2gbzs4BZLQ1QJBV9uhXwxfObHrTln2/v5LHF5ew5UNNOUYlkJj2kTbJGUafEDV97\nq3Snr2Q3PbJBskaPwkTzzhNLtrBhVxVnDOlB6ZCeIUcl0v6U+CVr9CrM54RuBcxevo3Zy7fRJT/O\nslsvCTsskXanxC9ZIy8e4x83TaS2vp7bn1nF7/6xgfp6JxZr7BKWSDSpjV+ySk7MyI/n0LdbAQCH\na9WvX7KPEr9kpfx44tA/VKO7fSX7KPFLVirITYzYpTN+yUZK/JKVCnITh/5BnfFLFtLFXclKvbvk\nA/DvD75OYV7qfwZXnzOUS087oa3CEmkXOuOXrDS6pIiLTu17XEl/yaa9PP3msWP/inQ0OuOXrNS1\nIJd7ppce1zqTfv4S1bomIBGgM36RFOXFY1Trsc4SAUr8IinKj+doQHeJBCV+kRTl5cTU1CORoMQv\nkqL83Jj6/UskKPGLpCgvJ8aBw7Xsevdw2KGItEqzid/MZpnZDjNb1qDsYTNbErw2NBiZ6+h1N5jZ\nm0G9snQGLtLeuhTEebviAB/6wXPMXrYt7HBEWiyV7py/B34F/OFIgbv/25FpM7sD2NfE+hPcfWdL\nAxTJFF+/+GRO69+d255awZa9B8MOR6TFmj3jd/d5QNJxcoOB2K8EHkxzXCIZZ2CPzkwdVwKgbp3S\nobW2jf9cYLu7r2lkuQNzzGyRmc1oakNmNsPMysysrKKiopVhibSN3JzEn4x690hH1trEP42mz/bH\nu/tY4DLgBjM7r7GK7j7T3UvdvbS4uLiVYYm0jXjMMIManfFLB9bixG9mceATwMON1XH38uB9B/A4\nMK6lnyeSCcxM/fmlw2vNGf9FwCp335xsoZkVmlnXI9PAJGBZsroiHUleXP35pWNLpTvng8ArwAgz\n22xm1wSLpnJUM4+Z9Tezp4PZvsB8M1sKvAb8zd1npy90kXDk65k90sE1253T3ac1Uv75JGVbgMnB\n9DpgVCvjE8k4eTkxnl22jVVb96dlexePPIEvXTAsLdsSSYUeyyxynD571mBeeXtXWra1cut+/rK4\nXIlf2pUSv8hxuv6Ck7j+gpPSsq0bH3id5VvS85+DSKr0rB6REOXF1UNI2p8Sv0iIdKFYwqDELxIi\n3RMgYVDiFwmRmnokDEr8IiHSOL4SBiV+kRDl5eRQV+/U1XvYoUgWUeIXCVFeXE/7lPanxC8SIiV+\nCYMSv0iI8nIMgMN1dSFHItlEiV8kRDrjlzAo8YuE6Ejir6nTxV1pP0r8IiHqlJsDwMQ7/s6En/1d\nvXukXSjxi4Ro/PBivnnJCM4Z1pv1Ow9wqEZt/dL2UhmIZZaZ7TCzZQ3K/tPMys1sSfCa3Mi6l5rZ\najNba2Y3pTNwkSjokh/nhgknMfGUPoDG8pX2kcoZ/++BS5OU/9zdRwevp49eaGY5wK9JDLQ+Ephm\nZiNbE6xIVMWD3j1q6pH20Gzid/d5wO4WbHscsNbd17l7NfAQMKUF2xGJvJyYEr+0n9a08d9oZm8E\nTUE9kiwfAGxqML85KBORo8SDxF+rxC/toKWJ/zfAMGA0sBW4o7WBmNkMMyszs7KKiorWbk6kQ4mZ\nzvil/bQo8bv7dnevc/d64H9INOscrRwoaTA/MChrbJsz3b3U3UuLi4tbEpZIh6U2fmlPLUr8Ztav\nwezHgWVJqi0EhpvZUDPLA6YCT7bk80SiLieW+FNUU4+0h2YHWzezB4ELgN5mthm4BbjAzEYDDmwA\nvhjU7Q/c4+6T3b3WzG4EngVygFnuvrxNvoVIBxfXxV1pR80mfneflqT43kbqbgEmN5h/Gjimq6eI\nvF/Oexd31Y9f2l6ziV9E2l5OcHH3/736Dn275Ycczb9cfnp/TurTJewwJM2U+EUywKBenSnIjfHg\na++EHcr7bN9/iB9/4vSww5A0U+IXyQAn9+3KytuS3SAfnrNvf4FaPTU0kpT4RTKEBc09mSJmhtJ+\nNOnpnCKSlBnUu1J/FCnxi0hSZqC8H01K/CKSVMwMV+aPJCV+EUkqZobuJ4smJX4RScpQG39UKfGL\nSFJq448uJX4RSSrRnVOZP4qU+EUkKTPQo4OiSYlfRJJKXNzVGX8UKfGLSFKmO3cjS4lfRJIyUD/+\niFLiF5GkYjHUjz+imk38ZjbLzHaY2bIGZT81s1Vm9oaZPW5mRY2su8HM3jSzJWZWls7ARaRt6c7d\n6ErljP/3wNHPi50LnObupwNvATc3sf4Edx/t7qUtC1FEwpC4gSvsKKQtNJv43X0esPuosjnuXhvM\nvgoMbIPYRCREpl49kZWONv4vAM80ssyBOWa2yMxmNLURM5thZmVmVlZRUZGGsESkNWKZNTyApFGr\nEr+ZfReoBe5vpMp4dx8LXAbcYGbnNbYtd5/p7qXuXlpcXNyasEQkDXTGH10tTvxm9nngcuAz3sgV\nIHcvD953AI8D41r6eSLSvmK6czeyWpT4zexS4FvAx9y9qpE6hWbW9cg0MAlYlqyuiGQe07N6IiuV\n7pwPAq8AI8xss5ldA/wK6ArMDbpq/jao29/Mng5W7QvMN7OlwGvA39x9dpt8CxFJu5ipV09UNTvY\nurtPS1J8byN1twCTg+l1wKhWRScioTEMd7X1RJHu3BWRpGIxPY8/qpT4RSQpPZ0zupT4RaRRauOP\nJiV+EUlKz+qJLiV+EUkqZqgzZ0Qp8YtIUrpzN7qU+EUkKd25G11K/CKSlIZejC4lfhFJSkMvRpcS\nv4gkpX780aXELyJJ6c7d6FLiF5GkDJ3xR5USv4gkZaYz/qhS4heRpGLq1RNZSvwikpQZauqJqJQS\nv5nNMrMdZrasQVlPM5trZmuC9x6NrDs9qLPGzKanK3ARaVvq1RNdqZ7x/x649Kiym4Dn3X048Hww\n/z5m1hO4BfgwifF2b2nsB0JEMova+KMrpcTv7vOA3UcVTwHuC6bvA65IsuolwFx33+3ue4C5HPsD\nIiIZKDECV9hRSFtoTRt/X3ffGkxvIzHG7tEGAJsazG8OykQkw8XUxh9Zabm464n7ult1hJjZDDMr\nM7OyioqKdIQlIq2QeB5/2FFIW2hN4t9uZv0AgvcdSeqUAyUN5gcGZcdw95nuXurupcXFxa0IS0TS\nIRbTGX9UtSbxPwkc6aUzHXgiSZ1ngUlm1iO4qDspKBORjGcaejGiUu3O+SDwCjDCzDab2TXA7cDF\nZrYGuCiYx8xKzeweAHffDXwfWBi8bgvKRCTDxQw0Blc0xVOp5O7TGll0YZK6ZcC1DeZnAbNaFJ2I\nhCbRjz/sKKQt6M5dEUlKd+5GlxK/iCQVM6Nep/yRlFJTj4hkHzOoPFzLxXe+FHYoWaNH5zz+dN1Z\nbf45SvwiktRHR/VnR+VhDb/YjroV5LbL5yjxi0hSYwf1YOyn9WitKFIbv4hIllHiFxHJMkr8IiJZ\nRolfRCTLKPGLiGQZJX4RkSyjxC8ikmWU+EVEsoxl4l15ZlYBbGzh6r2BnWkMJ6q0n1Kj/ZQa7afU\ntOV+GuzuKY1ilZGJvzXMrMzdS8OOI9NpP6VG+yk12k+pyZT9pKYeEZEso8QvIpJlopj4Z4YdQAeh\n/ZQa7afUaD+lJiP2U+Ta+EVEpGlRPOMXEZEmRCbxm9mlZrbazNaa2U1hx9PezKzEzF40sxVmttzM\nvhKU9zSzuWa2JnjvEZSbmf0y2F9vmNnYBtuaHtRfY2bTw/pObcnMcsxssZk9FcwPNbMFwf542Mzy\ngvL8YH5tsHxIg23cHJSvNrNLwvkmbcfMiszsETNbZWYrzewsHU/HMrOvBX9zy8zsQTMryPjjyd07\n/AvIAd4GTgTygKXAyLDjaud90A8YG0x3Bd4CRgL/BdwUlN8E/CSYngw8AxhwJrAgKO8JrAveewTT\nPcL+fm2wv74OPAA8Fcz/CZgaTP8W+FIwfT3w22B6KvBwMD0yOM7ygaHB8ZcT9vdK8z66D7g2mM4D\ninQ8HbOPBgDrgU4NjqPPZ/rxFJUz/nHAWndf5+7VwEPAlJBjalfuvtXdXw+mK4GVJA7KKST+gAne\nrwimpwB/8IRXgSIz6wdcAsx1993uvgeYC1zajl+lzZnZQOAjwD3BvAETgUeCKkfvpyP77xHgwqD+\nFOAhdz/s7uuBtSSOw0gws+7AecC9AO5e7e570fGUTBzoZGZxoDOwlQw/nqKS+AcAmxrMbw7KslLw\n7+MYYAHQ1923Bou2AX2D6cb2WTbsy18A3wLqg/lewF53rw3mG37n9/ZHsHxfUD/q+2koUAH8LmgS\nu8fMCtHx9D7uXg78DHiHRMLfBywiw4+nqCR+CZhZF+BR4Kvuvr/hMk/8T5nV3bjM7HJgh7svCjuW\nDBcHxgK/cfcxwAESTTvv0fEEwTWOKSR+KPsDhXSA/2iikvjLgZIG8wODsqxiZrkkkv797v5YULw9\n+Jeb4H1HUN7YPov6vjwH+JiZbSDRJDgRuItE00Q8qNPwO7+3P4Ll3YFdRH8/bQY2u/uCYP4REj8E\nOp7e7yJgvbtXuHsN8BiJYyyjj6eoJP6FwPDgSnoeiYsmT4YcU7sK2gnvBVa6+50NFj0JHOlJMR14\nokH554LeGGcC+4J/4Z8FJplZj+BsZlJQFgnufrO7D3T3ISSOkxfc/TPAi8CngmpH76cj++9TQX0P\nyqcGvTTXApkJAAAA6UlEQVSGAsOB19rpa7Q5d98GbDKzEUHRhcAKdDwd7R3gTDPrHPwNHtlPmX08\nhX1VPF0vEr0K3iJxNfy7YccTwvcfT+Lf7jeAJcFrMon2w+eBNcBzQM+gvgG/DvbXm0Bpg219gcTF\npbXA1WF/tzbcZxfwr149JwZ/aGuBPwP5QXlBML82WH5ig/W/G+y/1cBlYX+fNtg/o4Gy4Jj6C4le\nOTqejt1PtwKrgGXAH0n0zMno40l37oqIZJmoNPWIiEiKlPhFRLKMEr+ISJZR4hcRyTJK/CIiWUaJ\nX0Qkyyjxi4hkGSV+EZEs8/8BN0KwLvyU5b4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f0c4f981be0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot([x[1] for x in docsobj.tfidf_ranking])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "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": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "9488"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(docsobj.stems)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(docsobj.stems[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "all_stems = list(set(all_stems))\n",
    "corp_stems = [stem for doc in docsobj.stems for stem in doc]\n",
    "all_stems_dict = dict(zip(all_stems,range(len(all_stems))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Nd = len(corp_stems)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_stems_dict\n",
    "D = len(docsobj.stems)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Uncollapsed Gibbs sampler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def perp(thetas,betas,countsmat,docsobj,all_stems_dict,D,Nd):\n",
    "    s1 = 0\n",
    "    ltb = np.log(np.dot(thetas,betas))\n",
    "    for d in range(D):\n",
    "        s2 = 0\n",
    "        #counts = Counter(docsobj.stems[d])\n",
    "        for v in set(docsobj.stems[d]):\n",
    "            s2 += ((countsmat[d][v]) * ltb[d,all_stems_dict[v]])\n",
    "            #s2 += (counts[v] * ltb[d,all_stems_dict[v]])\n",
    "        s1 += s2\n",
    "    num = -1 * s1\n",
    "    den = Nd\n",
    "    return(np.exp(num/den))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 200/200 [15:11<00:00,  4.49s/it]\n"
     ]
    }
   ],
   "source": [
    "K = 10\n",
    "D = len(docsobj.stems)\n",
    "V = len(all_stems)\n",
    "\n",
    "# Initialise params\n",
    "\n",
    "eta = 200/V\n",
    "alpha = 50/K\n",
    "\n",
    "# Initialise betas and thetas\n",
    "betas = np.random.dirichlet(alpha = [eta]*V, size = K)\n",
    "thetas = np.random.dirichlet(alpha = [alpha]*K, size = D)\n",
    "\n",
    "perplexity = []\n",
    "countsmat = [Counter(doc) for doc in docsobj.stems]\n",
    "\n",
    "burnin = 200\n",
    "# Run iterations\n",
    "for l in tqdm(range(burnin)):\n",
    "    den = thetas.dot(betas)\n",
    "    Z = []\n",
    "    for d,doc in enumerate(docsobj.stems):\n",
    "        z_d = [np.argmax(np.random.multinomial(n=1, size=1, pvals=((thetas[d,:].dot(np.diag(betas[:,all_stems_dict[stem]]))/den[d,all_stems_dict[stem]]).tolist()))) for stem in doc]\n",
    "        Z.append(z_d)\n",
    "        \n",
    "    thetas = np.zeros(shape=(D,K))\n",
    "    N = np.zeros(shape=(D,K))\n",
    "    for i in range(D):\n",
    "        counts = Counter(Z[i])\n",
    "        N[i] = [counts[j] for j in range(K)]\n",
    "        theta = np.random.dirichlet(alpha=alpha + N[i], size=1)\n",
    "        thetas[i,:] = theta\n",
    "\n",
    "    betas = np.zeros(shape=(K,V))\n",
    "    z_s = [z for doc in Z for z in doc]\n",
    "    M = np.zeros(shape=(K,V))\n",
    "    for k in range(K):\n",
    "        stems = [corp_stems[i] for i in range(Nd) if z_s[i] == k]\n",
    "        counts = Counter(stems)\n",
    "        for stem in set(stems):\n",
    "            M[k,all_stems_dict[stem]] = counts[stem]\n",
    "        beta = np.random.dirichlet(alpha=eta + M[k], size=1)\n",
    "        betas[k,:] = beta\n",
    "    \n",
    "    if l%20 == 0:\n",
    "        perplexity.append(perp(thetas,betas,countsmat,docsobj,all_stems_dict,D,Nd))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1336.016028212824,\n",
       " 1305.4309505367669,\n",
       " 1285.3356394974182,\n",
       " 1261.8332868370364,\n",
       " 1238.6579964637383,\n",
       " 1212.1753259103232,\n",
       " 1191.6357751252826,\n",
       " 1177.903194602289,\n",
       " 1166.5075571946554,\n",
       " 1154.1395793416207]"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "perplexity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 500/500 [38:05<00:00,  4.46s/it]\n"
     ]
    }
   ],
   "source": [
    "m = 5\n",
    "n_samples = 100\n",
    "labels = []\n",
    "perplexity_pb = [] #perplexity post burn in\n",
    "for l in tqdm(range(m*n_samples)):\n",
    "    den = thetas.dot(betas)\n",
    "    Z = []\n",
    "    for d,doc in enumerate(docsobj.stems):\n",
    "        z_d = [np.argmax(np.random.multinomial(n=1, size=1, pvals=((thetas[d,:].dot(np.diag(betas[:,all_stems_dict[stem]]))/den[d,all_stems_dict[stem]]).tolist()))) for stem in doc]\n",
    "        Z.append(z_d)\n",
    "        \n",
    "    thetas = np.zeros(shape=(D,K))\n",
    "    N = np.zeros(shape=(D,K))\n",
    "    for i in range(D):\n",
    "        counts = Counter(Z[i])\n",
    "        N[i] = [counts[j] for j in range(K)]\n",
    "        theta = np.random.dirichlet(alpha=alpha + N[i], size=1)\n",
    "        thetas[i,:] = theta\n",
    "\n",
    "    betas = np.zeros(shape=(K,V))\n",
    "    z_s = [z for doc in Z for z in doc]\n",
    "    M = np.zeros(shape=(K,V))\n",
    "    for k in range(K):\n",
    "        stems = [corp_stems[i] for i in range(Nd) if z_s[i] == k]\n",
    "        counts = Counter(stems)\n",
    "        for stem in set(stems):\n",
    "            M[k,all_stems_dict[stem]] = counts[stem]\n",
    "        beta = np.random.dirichlet(alpha=eta + M[k], size=1)\n",
    "        betas[k,:] = beta\n",
    "    \n",
    "    if l%20 == 0:\n",
    "        perplexity_pb.append(perp(thetas,betas,countsmat,docsobj,all_stems_dict,D,Nd))\n",
    "    if l%m == 0:\n",
    "        #Z_s = [i for sublist in Z for i in sublist]\n",
    "        #j = np.int(i/m)\n",
    "        labels.append(Z)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "s1 = 0\n",
    "ltb = np.log(np.dot(theta_hat,beta_hat))\n",
    "for d in tqdm(range(D)):\n",
    "    s2 = 0\n",
    "    counts = Counter(docsobj.stems[d])\n",
    "    for v in set(docsobj.stems[d]):\n",
    "        s2 = s2 + (counts[v] * ltb[d,all_stems_dict[v]])\n",
    "    s1 = s1 + s2\n",
    "num = -1 * s1\n",
    "den = Nd\n",
    "perplex = np.exp(num/den)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 2 - Collapsed Gibbs sampling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "#ldaobj = topicmodels.LDA.LDAGibbs(docsobj.stems,10)\n",
    "#ldaobj.sample(0,5,100)\n",
    "#perp_coll = ldaobj.perplexity.tolist()\n",
    "\n",
    "# Done locally and reading from pickle\n",
    "with open(\"perp_values.txt\", \"rb\") as fp:   # Unpickling\n",
    "    perp_coll = pickle.load(fp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "perp_uncoll = perplexity + perplexity_pb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1145.9582550588393,\n",
       " 1140.0075408208832,\n",
       " 1136.1620095546934,\n",
       " 1130.8840733393151,\n",
       " 1126.1619659429737,\n",
       " 1125.3222809279741,\n",
       " 1120.1367240362799,\n",
       " 1117.5880299992764,\n",
       " 1115.6602966891207,\n",
       " 1113.9469299401314,\n",
       " 1109.8496199670797,\n",
       " 1109.6886209730935,\n",
       " 1106.6981598718305,\n",
       " 1104.5330305359876,\n",
       " 1103.0905874962293,\n",
       " 1103.9394173186133,\n",
       " 1103.8548235656399,\n",
       " 1102.4366737419161,\n",
       " 1100.5079732125766,\n",
       " 1100.1581177628598,\n",
       " 1099.2055481079776,\n",
       " 1097.2098259703744,\n",
       " 1095.9387679518231,\n",
       " 1094.9546720217372,\n",
       " 1094.6513861399242]"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "perplexity_pb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEKCAYAAADq59mMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd8lGXW8PHfyaSRkAYJNQkBhVBC6MXFgm3F3lfsva7r\nlsd31cdd2677uo/su7ZdVx7FXVexrL3giqiIriggIL2JlNACgSQEEkKS8/5x3YEhJGQCmZlkcr6f\nz3wyc93tzIzO4Sr3dYmqYowxxjRFVLgDMMYY0/pY8jDGGNNkljyMMcY0mSUPY4wxTWbJwxhjTJNZ\n8jDGGNNkQUseIjJJRApFZJFf2e9EZIGIzBeRqSLSzSsfKyIlXvl8EbnP75hxIrJcRFaJyN3BitcY\nY0zgJFj3eYjI8UAZ8IKq5nllyapa6j2/A+ivqreIyFjgTlU9q845fMAK4FSgAJgNXKqqS4IStDHG\nmIAEreahqjOA7XXKSv1eJgKNZa6RwCpVXa2qlcArwLnNGqgxxpgmiw71BUXkYeAqoAQ40W/TMSLy\nHbARVwtZDHQH1vvtUwCMCuQ66enpmpOT0ywxG2NMW/Dtt99uU9WMQPYNefJQ1XuBe0XkHuB24H5g\nLtBDVctE5AzgbaB3U88tIjcBNwFkZ2czZ86c5gvcGGMinIisDXTfcI62egm4EFxzlqqWec+nADEi\nkg5sALL8jsn0yuqlqhNVdbiqDs/ICCh5GmOMOQwhTR4i4l+bOBdY5pV3ERHxno/04irCdZD3FpGe\nIhILjAfeDWXMxhhjDha0ZisReRkYC6SLSAGueeoMEckFaoC1wC3e7hcBt4pIFVAOjFc3DKxKRG4H\nPgJ8wCSvL8QYY0wYBW2obrgNHz5crc/DmJZr7969FBQUUFFREe5Q2pz4+HgyMzOJiYk5oFxEvlXV\n4YGcI+Qd5sYYA1BQUEBSUhI5OTl4rdYmBFSVoqIiCgoK6Nmz52Gfx6YnMcaERUVFBR07drTEEWIi\nQseOHY+4xmfJwxgTNpY4wqM5PndLHn6qqmv4y2er+HzF1nCHYowxLZolDz++KGHijNVMXbw53KEY\nYyLA2LFj992snJOTw7Zt28IaQ3Oy5OFHVDkttYDyzSvCHYoxxrRoljwOoPy++G5GbH0r3IEYY0Jg\nzZo15OXl7Xs9YcIEHnjgAcaOHctdd93FyJEj6dOnD1988QUA1dXV3HnnneTl5ZGfn8+TTz4JwCef\nfMKQIUMYOHAg1113HXv27Dnkdc877zyGDRvGgAEDmDhx4r7y9u3b88tf/pIBAwZw8skns3Wra0J/\n4okn6N+/P/n5+YwfPx6AXbt2cd111zFy5EiGDBnCO++8A0B5eTnjx4+nX79+nH/++ZSXlzffB+bH\nhur6i/KxI/EoMkt+YNeeKhLj7OMxJhQefG8xSzaWNr5jE/Tvlsz9Zw847OOrqqqYNWsWU6ZM4cEH\nH2TatGlMnDiRNWvWMH/+fKKjo9m+fTsVFRVcc801fPLJJ/Tp04errrqKp59+ml/84hcNnnvSpEl0\n6NCB8vJyRowYwYUXXkjHjh3ZtWsXw4cP589//jMPPfQQDz74IE899RSPPPIIP/zwA3FxcRQXFwPw\n8MMPc9JJJzFp0iSKi4sZOXIkp5xyCs888wwJCQksXbqUBQsWMHTo0MP+DA7Fah51VHbsS9+odawp\n2hXuUIwxYXTBBRcAMGzYMNasWQPAtGnTuPnmm4mOdv+w7NChA8uXL6dnz5706dMHgKuvvpoZM2Yc\n8txPPPEEgwYNYvTo0axfv56VK1cCEBUVxSWXXALAFVdcwZdffglAfn4+l19+OS+++OK+a0+dOpVH\nHnmEwYMHM3bsWCoqKli3bh0zZszgiiuu2Hdcfn5+M34q+9k/reuI6TaQjLVvMq9gHQO6DQx3OMa0\nCUdSQzgS0dHR1NTU7Hvtf+9DXFwcAD6fj6qqqma75vTp05k2bRozZ84kISFh3w9/fWqH1H7wwQfM\nmDGD9957j4cffpiFCxeiqrzxxhvk5uY2W2xNYTWPOtJyBgNQtm5BmCMxxgRb586dKSwspKioiD17\n9vD+++8fcv9TTz2VZ555Zl8y2b59O7m5uaxZs4ZVq1YB8M9//pMTTjihwXOUlJSQlpZGQkICy5Yt\n4+uvv963raamhtdffx2AyZMnc+yxx1JTU8P69es58cQT+eMf/0hJSQllZWWcdtppPPnkk9ROMTVv\n3jwAjj/+eCZPngzAokWLWLAgOL9lljzqiOvuahtRhTb/ojGRLiYmhvvuu4+RI0dy6qmn0rdv30Pu\nf8MNN5CdnU1+fj6DBg1i8uTJxMfH8/zzz3PxxRczcOBAoqKiuOWWWxo8x7hx46iqqqJfv37cfffd\njB49et+2xMREZs2aRV5eHp9++in33Xcf1dXVXHHFFQwcOJAhQ4Zwxx13kJqaym9/+1v27t1Lfn4+\nAwYM4Le//S0At956K2VlZfTr14/77ruPYcOGNc+HVYdNjFiP4odymBM7nFPufr2ZozLG1Fq6dCn9\n+vULdxgtSvv27SkrKwvJter7/JsyMaLVPOqxLeEoulSsJlITqzHGHClLHvWo6NCXo3UdRTuDMz7a\nGGPqE6paR3Ow5FEPX9c84mUvG39YEu5QjDGmRbLkUY/U2hFXa78LcyTGGNMyWfKoR6deg6hWoWaL\njbgyxpj6WPKohy8ugQ2+biTsWB7uUIwxpkWy5NGAwnZH0bl8VbjDMMa0EA888AATJkwA4Jprrtl3\nM1+4Ygg3Sx4N2J2aS9eaLVRXtJ7RD8YYEyqWPBoQ1WUAUaJsW22d5sZEshdeeGHfHeNXXnkla9as\n4aSTTiI/P5+TTz6ZdevWHfL4hx56iBEjRpCXl8dNN9207/6wsWPH8vOf/5zBgweTl5fHrFmzAPj8\n888ZPHgwgwcPZsiQIezcuROARx99lBEjRpCfn8/999+/7/wPP/wwffr04dhjj2X58pbTlG4TIzYg\nKXswfAvFa+bRuf+YcIdjTGT78G7YvLB5z9llIJz+yCF3Wbx4Mb///e/56quvSE9PZ/v27Vx99dX7\nHpMmTeKOO+7g7bffbvAct99+O/fddx8AV155Je+//z5nn302ALt372b+/PnMmDGD6667jkWLFjFh\nwgT+8pe/MGbMGMrKyoiPj2fq1KmsXLmSWbNmoaqcc845zJgxg8TERF555RXmz59PVVUVQ4cODdp0\nI01lNY8GdOvZl10aR/WmReEOxRgTJJ9++ikXX3wx6enpgJtifebMmVx22WWASwa106I35LPPPmPU\nqFEMHDiQTz/9lMWL94/SvPTSSwE3WWFpaSnFxcWMGTOGX/3qVzzxxBMUFxcTHR3N1KlTmTp1KkOG\nDGHo0KEsW7aMlStX8sUXX3D++eeTkJBAcnIy55xzTpA+iaazmkcD0pPiWUo3oou/D3coxkS+RmoI\nLVVFRQW33XYbc+bMISsriwceeOCA6dVrp1T3f3333Xdz5plnMmXKFMaMGcNHH32EqnLPPfdw8803\nH7D/Y489FpL3cTiCWvMQkUkiUigii/zKficiC0RkvohMFZFuXrmIyBMissrbPtTvmKtFZKX3uDqY\nMftdkx0J2STvWhuKyxljwuCkk07iX//6F0VFRYCbYv1HP/oRr7zyCgAvvfQSxx13XIPH1yaK9PR0\nysrKDhqB9eqrrwLw5ZdfkpKSQkpKCt9//z0DBw7krrvuYsSIESxbtozTTjuNSZMm7ZueZMOGDRQW\nFnL88cfz9ttvU15ezs6dO3nvvfea/TM4XMGuefwdeAp4wa/sUVX9LYCI3AHcB9wCnA709h6jgKeB\nUSLSAbgfGA4o8K2IvKuqO4IcO9LhKDIKZrC7fDcJ7RKCfTljTIgNGDCAe++9lxNOOAGfz8eQIUN4\n8sknufbaa3n00UfJyMjg+eefb/D41NRUbrzxRvLy8ujSpQsjRow4YHt8fDxDhgxh7969TJo0CXC1\nic8++4yoqCgGDBjA6aefTlxcHEuXLuWYY44B3Oy6L774IkOHDuWSSy5h0KBBdOrU6aDzh1PQp2QX\nkRzgfVXNq2fbPUC2qt4qIs8A01X1ZW/bcmBs7UNVb/bKD9ivIUcyJXutJR8+Q/9vfs28cz5myNCR\nR3QuY8yBIn1K9rFjxzJhwgSGDw9ohvOQa5VTsovIwyKyHrgcV/MA6A6s99utwCtrqLy+894kInNE\nZM7WrVuPOM6s3m7t343fW6e5Mcb4C0vyUNV7VTULeAm4vRnPO1FVh6vq8IyMjCM+X1I3tzZw2cZl\nR3wuY0zbMn369BZb62gO4R6q+xJwofd8A5Dlty3TK2uoPPgSOrDLl0zUjtXU1NjCUMY0N1twLTya\n43MPefIQkd5+L88Fav9Z/y5wlTfqajRQoqqbgI+AH4tImoikAT/2ykKiIimHbtUbWb3NpikxpjnF\nx8dTVFRkCSTEVJWioiLi4+OP6DxBHW0lIi/jOrzTRaQAN2rqDBHJBWqAtbiRVgBTgDOAVcBu4FoA\nVd0uIr8DZnv7PaSq24MZt7/Yzr3J2TGdL9fu4OhOSaG6rDERLzMzk4KCApqjf9I0TXx8PJmZmUd0\njqAmD1W9tJ7i5xrYV4GfNrBtEjCpGUMLWPuuuSQtf4MFP2zikhHZ4QjBmIgUExNDz549wx2GOUzh\n7vNo8aTjUQBsWbM0zJEYY0zLYcmjMV7y8BX/wI5dlWEOxhhjWgZLHo3p4JJHL9nM3HVBv6ndGGNa\nBUsejYlPRhM70StqM9+uteRhjDFgySMg0vEo+sdtZd664nCHYowxLYIlj0B0PIosNrN4Y4mNSTfG\nGCx5BKbDUSRXFVFTUUrBjvJwR2OMMWFnySMQ3oirHrKFRRtKwhyMMcaEnyWPQHgjro7ybWHRRkse\nxhhjySMQHXoBMKz9dhZtKA1zMMYYE36WPAIRmwDJ3RkQt9U6zY0xBksegUvtQaZsY1tZJYU794Q7\nGmOMCStLHoFK6U5aVSGAdZobY9o8Sx6BSu5GbPkWoqTG+j2MMW2eJY9AJWci1ZUM6VBtI66MMW2e\nJY9AJXcDYHR6BYut2coY08ZZ8ghUSncA8pPL2FhSwXabnt0Y04ZZ8ghUskseveNdf8dia7oyxrRh\nljwClZAOvli6i1s+3TrNjTFtmSWPQEVFQVJX4so3k9WhnQ3XNca0aZY8miIlE0o2kNctxZqtjDFt\nmiWPpkjuBqUbyOuewpqi3ZRW7A13RMYYExaWPJoiuTuUbmRA1/YALLZ+D2NMG2XJoymSu0PNXgam\nuRqH9XsYY9qqoCUPEZkkIoUissiv7FERWSYiC0TkLRFJ9cpzRKRcROZ7j7/5HTNMRBaKyCoReUJE\nJFgxN8q716Nj1Va6psTbnebGmDYrmDWPvwPj6pR9DOSpaj6wArjHb9v3qjrYe9ziV/40cCPQ23vU\nPWfoeHeZU7qRvO4pLLSahzGmjQpa8lDVGcD2OmVTVbXKe/k1kHmoc4hIVyBZVb9Wt4jGC8B5wYg3\nIMleuKVuxNUP23ZRtqfq0McYY0wECmefx3XAh36ve4rIPBH5XESO88q6AwV++xR4ZeGR0BF8sVC6\ngYGZyajCko3WaW6MaXvCkjxE5F6gCnjJK9oEZKvqEOBXwGQRST6M894kInNEZM7WrVubL+BaUVGu\n6cq71wOs09wY0zaFPHmIyDXAWcDlXlMUqrpHVYu8598C3wN9gA0c2LSV6ZXVS1UnqupwVR2ekZER\nnDfgDdftlBxPp6Q4Sx7GmDYppMlDRMYBvwbOUdXdfuUZIuLznvfCdYyvVtVNQKmIjPZGWV0FvBPK\nmA+S3B1KXUtaXvcUG3FljGmTgjlU92VgJpArIgUicj3wFJAEfFxnSO7xwAIRmQ+8DtyiqrWd7bcB\nzwKrcDUS/36S0EvuBqWboKaGvO4prCosY3eldZobY9qW6GCdWFUvraf4uQb2fQN4o4Ftc4C8Zgzt\nyKRkQs1e2LWVvG7J1Cgs3bSTYT3Swh2ZMcaEjN1h3lT77vXYwMBM6zQ3xrRNljyaylsUitINdEmO\np2NirCUPY0ybY8mjqfYlj42IiN1pboxpkyx5NFWiW1GQEjfiqm/XJFZv3cXe6powB2aMMaFjyaOp\nRLwRVxsByO2cRGV1DWuLdoU5MGOMCR1LHocjOXNfzaNP5yQAlm3eGc6IjDEmpCx5HI60HlC8DoCj\nO7UnSmCFJQ9jTBtiyeNwpOXAzo2wt4L4GB856Yks32LJwxjTdljyOBxpOe6vV/vI7ZzEii1l4YvH\nGGNCzJLH4ahNHjvWAK7fY03RLsorq8MWkjHGhJIlj8OR2sP99ZJH3y5JqMKqQqt9GGPaBkseh6N9\nJ4hut7/m0cWNuLJ+D2NMW2HJ43CIuKYrL3n06JBAbHQUKyx5GGPaCEseh8sveUT7ojg6oz3Lbbiu\nMaaNCCh5iMi3IvJTEbF5x2ul5UDxWnCLIdK3S5IlD2NMmxFozeMSoBswW0ReEZHTvJX92q60HKgs\ng91FgOv32FxaQcnuveGNyxhjQiCg5KGqq1T1Xty64pOBScBaEXlQRDoEM8AWq85w3VxvmpIVhVb7\nMMZEvoD7PEQkH/gT8Chu1b+LgVLg0+CE1sLVvdejdsSVNV0ZY9qAgJahFZFvgWLcMrJ3q+oeb9M3\nIjImWMG1aKnZ7u+OHwDolhJPUly0JQ9jTJsQ6BrmF6vqav8CEempqj+o6gVBiKvli02A9p331TxE\nhD5dkuxeD2NMmxBos9XrAZa1LWk5sGPtvpe5XZJYtqmUmhoNX0zGGBMCh6x5iEhfYACQIiL+NYxk\nID6YgbUKaTmw9qt9LwdnpjL5m3Ws3raLozu1D19cxhgTZI01W+UCZwGpwNl+5TuBG4MVVKuRlgML\n/wVVlRAdy9AeqQDMXbfDkocxJqIdMnmo6jvAOyJyjKrODFFMrUdaDmgNlKyHjkfRK709yfHRzFu3\ng58Mzwp3dMYYEzSNNVv9WlX/B7hMRC6tu11V7whaZK2B/3DdjkcRFSUMyU5j7tricEZljDFB11iH\n+VLv7xzg23oeDRKRSSJSKCKL/MoeFZFlIrJARN4SkVS/bfeIyCoRWS4ip/mVj/PKVonI3U18f8FV\n514PgKHZaawo3Elphd1pboyJXIdMHqr6nvf0VVX9h/8D+KCRc/8dGFen7GMgT1XzgRXAPQAi0h8Y\nj+ucHwf8VUR8IuID/gKcDvQHLvX2bRnadwFf3IHJo0cqqvDdeqt9GGMiV6BDdWeJyOjaFyJyIfDV\nIfZHVWcA2+uUTVXVKu/l10Cm9/xc4BVV3aOqPwCrgJHeY5WqrlbVSuAVb9+WISrK3SzolzwGZ6Ui\ngjVdGWMiWqA3CV4OTBKR6bgJEjsCJx3hta8DXvWed8clk1oFXhnA+jrlo47wus0rLWffXeYASfEx\n9OmUxNx1O8IXkzHGBFmgEyMuBB4GbgFOBG5X1YLDvaiI3AtUAS8d7jkaOO9NIjJHROZs3bq1OU/d\nsE59YetyqN7fxzG0Ryrz1u2wmwWNMREr0PU8ngN+AeQD1wLvi8hPD+eCInIN7t6Ry1W19td1A+A/\ntjXTK2uovF6qOlFVh6vq8IyMjMMJr+m6DobqSihcuq9oSHYapRVVrN5ma5obYyJToH0eC4ETvbms\nPsI1HQ1t6sVEZBzwa+AcVd3tt+ldYLyIxIlIT6A3MAuYDfQWkZ4iEovrVH+3qdcNqq6D3d9N3+0r\nGprt1syyfg9jTKQKtNnqMSBeRHK91yWqev2hjhGRl4GZQK6IFIjI9cBTQBLwsYjMF5G/eedbDLwG\nLAH+DfxUVau9zvXbgY9ww4Zf8/ZtOTr0gtgk2DR/X1Gv9ERS2sVYv4cxJmIFOiX72cAEIBboKSKD\ngYdU9ZyGjlHVg24qxE3p3tD+D+P6VeqWTwGmBBJnWERFQdf8A2oe7mbBVEsexpiIFWiz1QO4YbPF\nAKo6H+gVpJhan66DYfMiqK7aVzQipwMrtpSxuaQijIEZY0xwBJo89qpqSZ2ymuYOptXqOgiqymHb\nin1Fpw3oDMCHizaFKypjjAmaQJPHYhG5DPCJSG8ReZJGbhJsU7oOcn/9+j2O7pREbuckpiy05GGM\niTyBJo+f4aYO2QO8jFu7/BfBCqrVSe8NMQkH9HsAnDGwK3PW7mBLqTVdGWMiS6CjrXar6r2qOsK7\nj+JeVbVfxFpRPugy8KDkcWZ+F1ThQ6t9GGMiTGNTsr8HNHib9KFGW7U5XQfDvBehptolE/Y3XX2w\ncBPXjOkZ5gCNMab5NDZUd0JIoogEXQfBrGeg6HvI6LOv+IyBXXnskxVsKa2gc7Kt3GuMiQyNTcn+\nee0Dd8PfDtxMuTO9MlOrW+2d5vMPKLamK2NMJAp0bqszge+BJ3B3ia8SkdODGVirk54L0fEH9Xvs\nH3W1OUyBGWNM8wt0tNWfcHNbjVXVE3Az6/45eGG1Qr5o6JwHG+cftOmMgV2ZvXY7G4rLwxCYMcY0\nv0CTx05VXeX3ejWwMwjxtG7dBruaR031AcUXDnNLk7w6a104ojLGmGYXaPKYIyJTROQaEbkaeA+Y\nLSIXiMgFQYyvdckcCZU7D5ieHSAzLYET+mTw6pz1VFXbjfnGmNYv0OQRD2wBTgDGAluBdsDZuLU5\nDEDWSPd3/TcHbbpsZDZbSvfwybLCEAdljDHNr9FZdUXEByxQVevjaExaDiR2csljxIEz1p/UtxOd\nk+N4edY6ThvQJTzxGWNMM2m05qGq1UB906ubukQge1S9NY9oXxSXDM/i8xVbWb99dz0HG2NM6xFo\ns9V/ROQpETlORIbWPoIaWWuVNQp2rIGdWw7adMnIbAR4dfb6kIdljDHNKdDkMRg3MeJDuGG7f8Lu\nPq9f1ij3t2DWQZu6p7ZjbG4nXp2znr3WcW6MacUCWklQVU8MdiARo+sg8MW6pqt+Zx+0efyILD5d\nVsjny7dySv/OYQjQGGOOXKB3mHcWkedE5EPvdX9vTXJTV3QcdBsC6w+ueQCc2LcTHRJjeXNeQYgD\nM8aY5hNos9XfgY+Abt7rFdh6Hg3LGgkb50HVnoM2xfiiOGdQN6YtKaR4d2UYgjPGmCMXaPJIV9XX\n8JaeVdUqoPrQh7RhWaOhurLeqUoALhqWSWV1De8vsMkSjTGtU6DJY5eIdMRb20NERgN11zQ3tQ5x\nsyDAgG7J5HZO4o251nRljGmdAk0evwLeBXqJyH+AF3BL05r6tO8EaT0bTB4iwgVDuzNvXTHfby0L\ncXDGGHPkAk0eS4C3gNm4aUr+F9fvYRqS5d0sqPUvxHjekO5ECbw1d0OIAzPGmCMXaPJ4AegL/AF4\nEugD/DNYQUWEo06EXVuhYHa9mzsnx3Ns7wzemreBmpoGV/o1xpgWKdDkkaeqN6jqZ97jRtxNgw0S\nkUkiUigii/zKLhaRxSJSIyLD/cpzRKRcROZ7j7/5bRsmIgtFZJWIPCEi0tQ3GRa5Z4AvDha92eAu\nFw3LZENxOVOXHHw3ujHGtGSBJo+5Xic5ACIyCpjTyDF/B8bVKVsEXADMqGf/71V1sPe4xa/8aeBG\noLf3qHvOlik+GXqfCkvehpr67yY/I68LR3dqzx//vYzKKrvj3BjTegSaPIYBX4nIGhFZg1vPfIRX\nI1hQ3wGqOgO33rl/2VJVXR5ocCLSFUhW1a9VVXHNZ+cFenzYDTgfdm6CdTPr3Rzti+K/z+jLD9t2\n8dI3a0McnDHGHL6ApichNP/a7yki84BS4Deq+gXQHfAfz1rgldVLRG4CbgLIzs4OYqgB6jMOotvB\n4jchZ0y9u5yY24ljj07n8U9WcsGQTFISYkIcpDHGNF1ANQ9VXXuoRzPEsQnIVtUhuGHBk0Ukuakn\nUdWJqjpcVYdnZGQ0Q1hHKK499DkNlrwD1VX17iIi/PcZ/Sgp38tTn60McYDGGHN4Am22CipV3aOq\nRd7zb4HvcSO6NgCZfrtmemWtR94FbtTVmi8a3KV/t2QuHpbJ379aY/d9GGNahRaRPEQkw1uxEBHp\nhesYX62qm4BSERntjbK6CngnjKE2Xe8fQ2x713R1CHf+OJfEuGh++tJcyitt5hdjTMsWtOQhIi/j\nOtZzRaRARK4XkfNFpAA4BvhARD7ydj8eWCAi84HXgVtUtbaz/TbgWWAVrkbyYbBiDoqYdpB7Oix9\nr96JEmt1So7nsUsGs3zLTn7z9iK0gZsLjTGmJQi0w7zJVLWhpWvfqmffN4A3GjjPHCCvGUMLvcGX\nwcJ/wYLXYOiVDe42NrcTPzupN098spIROWmMH9kCOv2NMaYeLaLZKuL1OhG6DISvnmjwno9aPz+5\nN8f1Tue+dxfz9eqiEAVojDFNY8kjFERgzC9g2wpY8e9D7uqLEh4fP4SstHZc+dw3vDbH1js3xrQ8\nljxCpf95kJoN/3ms0V07JMby5q1jGNWzI79+fQH/98OlNv+VMaZFseQRKr5oOOZnbqbddV83untK\nQgzPXzuCK0Zn88znq/nNO9aJboxpOSx5hNKQy6FdB/jP4wHtHuOL4nfn5nHzCb2Y/M06Js5YHeQA\njTEmMJY8Qik2EUbeBMunwJYlAR0iItx1Wl/OzO/K//1wGR/Y0rXGmBbAkkeojboZYpPg80cCPiQq\nSvjTxYMY1iONX742nzlrtjd+kDHGBJElj1BL6ACjb3HzXW1e1Pj+nvgYH/971XAyU9tx7fOzWVBQ\nHMQgjTHm0Cx5hMMxP4W45CbVPsCNwnrpxlGkJsZw5XOzWLKxNEgBGmPMoVnyCId2aTD6NjdlyaZ6\nl0NpUNeUdky+YTQJsT6ueO4b3vtuI2uLdtlQXmNMSFnyCJfRt0JcCkxvWu0DIKtDApNvHE2sL4qf\nvTyPEx6dTv6DU5nw0XIbzmuMCYmgzW1lGtEu1TVfTf8DFMyBzOGNH+OnZ3oi0//PWJZv3snSTaV8\ntryQpz5bRUq7GG48vleQgjbGGMdqHuF0zG3QvjN8+OtG57yqT3yMj0FZqYwfmc3Tlw/jzIFdeXjK\nUt77bmMQgjXGmP0seYRTXBKc+hBs+Ba+m3xEp4qKEv70k0GMyEnjv177ji9XbmumII0x5mCWPMIt\n/xLIHAlugVoQAAAVc0lEQVTTHoCKkiM6Ve1w3uyOCVzx3Df84pV5bCgub544jTHGjyWPcBOBMx6F\nXdsOq/O8rtSEWN667UfcNvYoPly0mRMnTOf37y9h/fbdzRCsMcY4Eqmjc4YPH65z5swJdxiBe+/n\nMPefcPMM6NI8a19tKC7nTx8t553vNlKjykm5nbjimB4c3zsDX5Q0yzWMMZFDRL5V1YBG71jyaCl2\nFcFfRkJyN7jhE4iObbZTbyop5+Vv1jF51jq2lVXSOTmO84dkcsmILHqmJzbbdYwxrVtTkoc1W7UU\niR3h7Mdh8wL4YkKznrprSjt+9eNc/nP3SfztiqEM7J7C/36xmtMem8GHC22iRWNM01nyaEn6nQX5\n42HGBNgwt9lPHxftY1xeV569egRf3X0SA7uncNvkubwwc02zX8sYE9ms2aqlKS+Gvx7jhvHePANi\n4oN2qYq91dw+eR7Tlm7h0pHZ9O2SRHxMFJ2S4zmhdwZR1i9iTJvSlGYru8O8pWmXCuc+BS9eAFN/\nA2c2bxOWv/gYH3+7YigPvLeYF79ed8C2UT078McL88mxPhFjTD2s5tFSfXQvzHwKfvJP6H9O0C+3\nu7KK8spqKqpqmLFiK3/4YCl7a2q488e5XHlMD+KifUGPwRgTXjbaighIHlWVMOk02P493PwFpPUI\n6eU3l1Twm7cXMm1pIV2S47nhuJ5cOjKbxDirrBoTqVpE8hCRScBZQKGq5nllFwMPAP2Akao6x2//\ne4DrgWrgDlX9yCsfBzwO+IBnVTWgO+laffIA2P4DPHM8ZOTCtR+CLyakl1dVvli5jb9OX8XXq7eT\n0i6G0wZ05rQBXRhzdDol5XtZsrGUlYU72VutxPiEGF8U+ZkpDMlKsz4TY1qZlpI8jgfKgBf8kkc/\noAZ4BrizNnmISH/gZWAk0A2YBvTxTrUCOBUoAGYDl6pqowuAR0TyAFj0Jrx+LYy8Gc74n7CFMXfd\nDv7x1Ro+XVrIzj1V+KKE6kOsIZLePpZT+nXmvCHdGdWzAyKWSIxp6VpEh7mqzhCRnDplS4H6fkjO\nBV5R1T3ADyKyCpdIAFap6mrvuFe8fRtNHhEj7wIomA1f/xU6D4BhV4cljKHZaQzNTmNPVTUzvy9i\n5uoiuibH079bCrldkoiLjqKqRtldWcXM74v4eMkW3l+wiVdmrye3cxJX/agHxx6dTlpiLElx0ZZM\njGnlWkoDdnfga7/XBV4ZwPo65aNCFVSLcervYOsy+OC/IL039PhR2EKJi/YxNrcTY3M71bu9fVw0\n5w7uzrmDu1NeWc17323k71+t4d639q/XHh0lZHVIYFBmCvmZqRzXO53enZNC9RaMMc2gpSSPZiEi\nNwE3AWRnZ4c5mmbki4aLJsGzp8CrV8KNn4a8A/1wtIv18ZMRWVw8PJPvCkpYVVjGjl2VbN9dyarC\nMmauLuLt+W7tkWOPTufaMTmcmNvJ+kqMaQVaSvLYAGT5vc70yjhE+UFUdSIwEVyfRzPHGF7t0uDS\nV+DZk+HvZ8KVb0P60eGOKiAiwuCsVAZnpR60bVNJOW/N28ALX63l+n/MoWtKPCf0yeD4PhmMOSqd\nlITQDhIwxgQmqEN1vT6P92s7zP3Kp3Ngh/kAYDL7O8w/AXoDguswPxmXNGYDl6nq4sauHTEd5nVt\n+g7+eYGbyv2KN6FrfrgjahZ7q2v4cNFmPliwka9WFbFzTxUAvTISGZyVSv+uyaQlxJLcLoaUdjEM\nzkolNtpm1zGmObWU0VYvA2OBdGALcD+wHXgSyACKgfmqepq3/73AdUAV8AtV/dArPwN4DDdUd5Kq\nPhzI9SM2eQBsWwUvnAt7SuHy1yE7srqB9lbXMH99Md+sLmL++mLmry9mW1nlAfvM++2ppCU238zD\nxpgWkjzCLaKTB0BJAfzjHCjf7vpAOvQKd0RBo6qUlO894PGjo9JtTRJjmplNyd4WpGTC5f9yzydf\n4iZUjFAiQmpCLD06Jnqjs2wxK2PCzZJHa9bxKDf31fbV7kbC6qpwR2SMaSMsebR2PY+Ds/4M338K\n794O1XvDHZExpg1oKUN1zZEYehWUboLpf4Cdm+EnL0B8crijMsZEMKt5RIqxd8E5T8GaL2DSONeh\nbowxQWLJI5IMvdIN3S1ZDxPHwg8zwh2RMSZCWfKINEedCDdMg3Yd3L0gMyZATU24ozLGRBhLHpEo\nI9fd+zHgfPj0d/Di+bB1ebijMsZEEEsekSquPVz4HJz5J9gwF57+EXx4N5TvCHdkxpgIYMkjkonA\niBvgZ3NhyBUw6xl4chgs+BdE6MwCxpjQsOTRFrTPgLMfh5s+h7QcePMGeOki2LE23JEZY1opSx5t\nSdd8uP5jGPdHWDsTnhwKk8fDknegak+4ozPGtCJ2k2BbE+WD0bdA3zNh1kRY8Bqs+BASM+CEu2DY\nNeCzNTSMMYdmNY+2KjULfvw7+OViuPwNyOgLU+6Ev46Gpe9bn4gx5pAsebR1vmjofQpc/Z5bqVCi\n4NXL4ekxrmPdJls0xtTDkodxRCD3dLh1Jpz3N9Bq17H+5FD4+m+wpyzcERpjWhBLHuZAvmgYfKlL\nIuNfhqQu8O+74M/94eP7oWxruCM0xrQAljxM/aKioO8ZcP1UuH4a9DoRvnoCHh8E0x6A3dvDHaEx\nJowseZjGZY2An/wDfjrLNW19+ZhLIp8/CpW7wh2dMSYMLHmYwKX3houeg1u/gp7Hw2e/h8cHw+xn\nYc/OcEdnjAkh0Qgdkjl8+HCdM2dOuMOIbOtnwcf3wbqZEBUNmSOh11gYeJFbItcY06qIyLeqOjyg\nfS15mCOiCmu/glUfw/efwabvXHnvH7ubEbNGQ3S860MxxrRoljyw5BE2OzfDnOdhznOwy29kVnQ7\n6Nwfcs+Avme5aeNFwhenMeYgljyw5BF2VXtg2ftQvB72lkNlmWve2vCt256QDh2Pds1bXfLh6FPc\nc0soxoRNU5KHzW1lgiM6DvIuPLi8dCMs/9A1bxV9D6umwfyX3La0HMg+BlKzISUTkrpBuzRolwrt\nO7s1SowxLULQkoeITALOAgpVNc8r6wC8CuQAa4CfqOoOERkLvAP84B3+pqo+5B0zDngc8AHPquoj\nwYrZhEByNxhx/YFl239wSWTVNFj9OezcBNSpEUsUdB/m7jfpPtQNES7fATXVrpPemsGMCamgNVuJ\nyPFAGfCCX/L4H2C7qj4iIncDaap6l5c87lTVs+qcwwesAE4FCoDZwKWquqSx61uzVStWVQk7N7r+\nk/JiqCiGolWwerpr9tJ61mTv0MvNFNz3bMgcYR30xhyGFtFspaozRCSnTvG5wFjv+T+A6cBdhzjN\nSGCVqq4GEJFXvHM0mjxMKxYd65qw0nIOLD/pN662sW0lxKe4Jq3qSljxESyf4ubg+upJ18R19ClQ\nvRfKNsOuIkhM987ZAxI6uuPjkt1U9O07u+1RvjC8WWNap1D3eXRW1U3e881AZ79tx4jId8BGXC1k\nMdAdWO+3TwEwKiSRmpapXRpkjTywbMT17lFRAis/hqXvuX6VuCQ3N1daDzfya/mUA0eA+ZMo6Ngb\nskdB1ijoOsjVZmIT3faqSihZ7zr+o+PdI6Gj9cOYNitsHeaqqiJS22Y2F+ihqmUicgbwNtC7qecU\nkZuAmwCys7ObLVbTSsSnuBsUB17U8D6Vu13tpaLEPXZvg7Itrols03duVcW5L+zfP6mb60sp3chB\n/TAAHY6CboNdraZ6rxtlJlGu079DT0jt4Wo37dLcpJOqbp+qClfzseY100qFOnlsEZGuqrpJRLoC\nhQCqWlq7g6pOEZG/ikg6sAHI8js+0yurl6pOBCaC6/MIxhswrVxsgnukdK9/e00NbFsOhUth+/du\nRBi4JJCa7RJUVYV7lG50CWf9bFj0plcjiYOaKldDOeja7d2wZa12r8XnmssSO7l4UjLdgAJk/zVq\nE01VpTtvbQKLjoO4FFe7imkHvli3AmRse1cjSujo3mdNlRtUoNUucWmNu267VJfQ4lNcsrPBBqaJ\nQp083gWuBh7x/r4DICJdgC1ebWQkbs6tIqAY6C0iPXFJYzxwWYhjNm1JVBR06uceh0sVdhe5UWQl\n69wMxLuLXE0npp1rCvPFuRrQrkIoK4TSDbD+G1dWyxfrbq6MjnOPKB/g/chX7YE9pfUnqcMi7vxR\n0RAV4xJRfDLEp7pEk5K1vyZVe/29u91Ahi2LYdsKV5PKyIX0Pi55RUW72lZVpYtzz0533vZdIKmz\ne28VJbCnxJ0zubv36Oau62vg56mm2g2kaJdq/VRhFMyhui/jOsfTRaQAuB+XNF4TkeuBtcBPvN0v\nAm4VkSqgHBivbhhYlYjcDnyEG6o7yesLMablEvFqFOluRuKmqNztagK+2MCatGqqXc2kunL/j3Rt\nstq720sG0e6cUT4XW3WVG8FWvgMqSl1tRKvduWqq3KM2OVWUuPOt+MglurpiElyiPfoUt++2FbDi\n314tqY7oeNe0V1vzakztoIi4ZK+Jz+f6nYrXQ81e955qBzy077S/ebBsC+xY6/ZV9ZJwrNuW2Ana\nZ7jPYHeRe+zZ6T6rvbtdUm/fyfWVRce576Nyl993mgFJXd2yzRm5rrnSF+PiVXXNn0Wr3D8GJMpt\ni27n9uvQ053TX/Ve94+HXYUuacd7tcnyHVC8DkoK3LUT0r0aZZpLrHHJbp8tC2HzQncOEXfNuGQ4\n/s7APuMjYHeYG2MCU7nL/XCLeLWhdu7HtG6Sq97rmudqqtxzX4z7QfTFuAS1u8j9yFZX7h/1pjWu\nGbC0wG3bvR3Kt7u/e3a6R3UlpGa52k9SF7etbDPs3OLV4La6H9T2Ge7HOiXLJc7qSpdgd293AybK\nCl387dIgoYO7fmyiS4RVe1zyKdvsknFsontoDeza5o4vr7OWTVS0S4w11VBV3vDnJ1GuaVKi3Lmr\nKryaZjP8BsckuvPUVLsk96vDG5DaIobqGmMiTGwidOrb+H6+mP3/Gq8ryuf+Zd++08HbkrsCw44o\nxJDYU+ZqWFuXu9pBbd+UiEtaHY92/WOqroZUuQu2r3bH7FgDiKsJRce7GkVSF/d51FS5mmBFyf6m\nwlSvy3dXkRvcUV7sDfYodv1bXfOhc55LgiFmycMYY5oirr2b5aD70MCPyQzoH/MN69DryI4PAhsn\naIwxpskseRhjjGkySx7GGGOazJKHMcaYJrPkYYwxpskseRhjjGkySx7GGGOazJKHMcaYJovY6UlE\nZCtu/qzDkQ5sa8ZwWoO2+J6hbb7vtvieoW2+76a+5x6qmhHIjhGbPI6EiMwJdH6XSNEW3zO0zffd\nFt8ztM33Hcz3bM1WxhhjmsyShzHGmCaz5FG/ieEOIAza4nuGtvm+2+J7hrb5voP2nq3PwxhjTJNZ\nzcMYY0yTWfLwIyLjRGS5iKwSkbvDHU+wiEiWiHwmIktEZLGI/Nwr7yAiH4vISu9vWrhjbW4i4hOR\neSLyvve6p4h8433nr4pIbLhjbG4ikioir4vIMhFZKiLHRPp3LSK/9P7bXiQiL4tIfCR+1yIySUQK\nRWSRX1m93604T3jvf4GINGFBkoNZ8vCIiA/4C3A60B+4VET6hzeqoKkC/ktV+wOjgZ967/Vu4BNV\n7Q184r2OND8Hlvq9/iPwZ1U9GtgBXB+WqILrceDfqtoXGIR7/xH7XYtId+AOYLiq5gE+YDyR+V3/\nHRhXp6yh7/Z0oLf3uAl4+kgubMljv5HAKlVdraqVwCvAuWGOKShUdZOqzvWe78T9mHTHvd9/eLv9\nAzgvPBEGh4hkAmcCz3qvBTgJeN3bJRLfcwpwPPAcgKpWqmoxEf5d41ZJbSci0UACsIkI/K5VdQZQ\nZ1H1Br/bc4EX1PkaSBWRrod7bUse+3UH1vu9LvDKIpqI5ABDgG+Azqq6ydu0GegcprCC5THg10CN\n97ojUKyqVd7rSPzOewJbgee95rpnRSSRCP6uVXUDMAFYh0saJcC3RP53Xauh77ZZf+MsebRhItIe\neAP4haqW+m9TNwwvYobiichZQKGqfhvuWEIsGhgKPK2qQ4Bd1GmiisDvOg33r+yeQDcgkYObdtqE\nYH63ljz22wBk+b3O9MoikojE4BLHS6r6ple8pbYa6/0tDFd8QTAGOEdE1uCaJE/C9QWkek0bEJnf\neQFQoKrfeK9fxyWTSP6uTwF+UNWtqroXeBP3/Uf6d12roe+2WX/jLHnsNxvo7Y3IiMV1sL0b5piC\nwmvrfw5Yqqr/z2/Tu8DV3vOrgXdCHVuwqOo9qpqpqjm47/ZTVb0c+Ay4yNstot4zgKpuBtaLSK5X\ndDKwhAj+rnHNVaNFJMH7b732PUf0d+2noe/2XeAqb9TVaKDEr3mryewmQT8icgauXdwHTFLVh8Mc\nUlCIyLHAF8BC9rf//zeu3+M1IBs3I/FPVLVuZ1yrJyJjgTtV9SwR6YWriXQA5gFXqOqecMbX3ERk\nMG6QQCywGrgW9w/HiP2uReRB4BLcyMJ5wA249v2I+q5F5GVgLG723C3A/cDb1PPdeon0KVwT3m7g\nWlWdc9jXtuRhjDGmqazZyhhjTJNZ8jDGGNNkljyMMcY0mSUPY4wxTWbJwxhjTJNZ8jCmBRGR6SLS\nptbZNq2TJQ9jjDFNZsnDmEaISKKIfCAi33nrQ1wiIveJyGzv9UTvBqzamsOfRWSOt3bGCBF501tb\n4ffePjne2hovefu8LiIJ9Vz3xyIyU0Tmisi/vLnIEJFHxK3FskBEJoT20zDGseRhTOPGARtVdZC3\nPsS/gadUdYT3uh1wlt/+lao6HPgbbmqInwJ5wDUi0tHbJxf4q6r2A0qB2/wvKCLpwG+AU1R1KDAH\n+JV3/PnAAFXNB34fnLdszKFZ8jCmcQuBU0XkjyJynKqWACd6q9ItxE2yOMBv/3f9jlvsrZ+yBzc1\nSO3EdOtV9T/e8xeBY+tcczRuUbL/iMh83BxFPXDTi1cAz4nIBbhpJowJuejGdzGmbVPVFd6SnWcA\nvxeRT3C1ieGqul5EHgDi/Q6pnS+pxu957eva/+fqzgtU97UAH6vqpXXjEZGRuMn+LgJuxyUvY0LK\nah7GNEJEugG7VfVF4FHclOYA27x+iIsaPLhh2SJyjPf8MuDLOtu/BsaIyNFeDIki0se7XoqqTgF+\niVtW1piQs5qHMY0bCDwqIjXAXuBW3NKei3Artc0+jHMux60dPwk3XfgB60mr6lYRuQZ4WUTivOLf\nADuBd0QkHlc7+dVhXNuYI2az6hoTYt7Sv+97ne3GtErWbGWMMabJrOZhjDGmyazmYYwxpskseRhj\njGkySx7GGGOazJKHMcaYJrPkYYwxpskseRhjjGmy/w+9kUVzYvTMAgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f0c4f9af940>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(perp_uncoll,label = 'uncollapsed')\n",
    "plt.plot(perp_coll,label = 'collapsed')\n",
    "plt.legend()\n",
    "plt.ylabel('perplexity')\n",
    "plt.xlabel('samples')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### We can see that the collapsed sampler stabilizes much faster as compared to the uncollapsed sampler. I take the perplexity every 20 iterations for the uncollapsed sampler"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Predictive theta and beta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "beta_hat = np.zeros(shape=(K,V))\n",
    "for i in tqdm(range(K)):\n",
    "    den = (M[i,:]+eta).sum()\n",
    "    for j in range(V):\n",
    "        beta_hat[i,j] = (M[i,j] + eta)/den"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "theta_hat = np.zeros(shape=(D,K))\n",
    "for i in tqdm(range(D)):\n",
    "    den = (N[i,:]+alpha).sum()\n",
    "    for j in range(K):\n",
    "        theta_hat[i,j] = (N[i,j] + alpha)/den"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LogisticRegressionCV\n",
    "from sklearn.model_selection import train_test_split\n",
    "#from sklearn.model_selection import GridSearchCV\n",
    "from sklearn.metrics import classification_report\n",
    "from sklearn.feature_extraction.text import CountVectorizer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Create a new variable with only presidents and years\n",
    "parties = data.loc[:,['president', 'year']]\n",
    "parties = parties.reset_index()\n",
    "# Create a new variable with 1 when presidents are Democrats, 0 when they are Republicans\n",
    "zero_len = pd.Series(np.zeros(len(parties.index)))\n",
    "parties['parties'] = zero_len\n",
    "parties['parties'] = (parties.president == \"Truman\") | (parties.president == \"Kennedy\") | (parties.president == \"Johnson\") | (parties.president == \"Carter\") | (parties.president == \"Clinton\") | (parties.president == \"Obama\")\n",
    "parties.parties = list(map(lambda x: 1 if x else 0, parties.parties))\n",
    "parties = parties.iloc[:][parties.year >= 1947]\n",
    "# Set the parties variable as y\n",
    "y = parties.parties"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Accuracy: 0.6761503336845802\n",
      "             precision    recall  f1-score   support\n",
      "\n",
      "          0       0.66      0.74      0.70      1444\n",
      "          1       0.69      0.61      0.65      1403\n",
      "\n",
      "avg / total       0.68      0.68      0.67      2847\n",
      "\n",
      "Tuned Model Parameters: [ 0.35938137]\n",
      "AUC score when using document-term matrix: 0.7399974431520899\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:lda:n_documents: 9488\n",
      "INFO:lda:vocab_size: 4742\n",
      "INFO:lda:n_words: 250000\n",
      "INFO:lda:n_topics: 10\n",
      "INFO:lda:n_iter: 1000\n",
      "WARNING:lda:all zero row in document-term matrix found\n",
      "INFO:lda:<0> log likelihood: -2661826\n",
      "INFO:lda:<10> log likelihood: -2077618\n",
      "INFO:lda:<20> log likelihood: -1987666\n",
      "INFO:lda:<30> log likelihood: -1953114\n",
      "INFO:lda:<40> log likelihood: -1931842\n",
      "INFO:lda:<50> log likelihood: -1921545\n",
      "INFO:lda:<60> log likelihood: -1915514\n",
      "INFO:lda:<70> log likelihood: -1910984\n",
      "INFO:lda:<80> log likelihood: -1906653\n",
      "INFO:lda:<90> log likelihood: -1904568\n",
      "INFO:lda:<100> log likelihood: -1903775\n",
      "INFO:lda:<110> log likelihood: -1901340\n",
      "INFO:lda:<120> log likelihood: -1900211\n",
      "INFO:lda:<130> log likelihood: -1900286\n",
      "INFO:lda:<140> log likelihood: -1898358\n",
      "INFO:lda:<150> log likelihood: -1899064\n",
      "INFO:lda:<160> log likelihood: -1898420\n",
      "INFO:lda:<170> log likelihood: -1898299\n",
      "INFO:lda:<180> log likelihood: -1898357\n",
      "INFO:lda:<190> log likelihood: -1896960\n",
      "INFO:lda:<200> log likelihood: -1896936\n",
      "INFO:lda:<210> log likelihood: -1896832\n",
      "INFO:lda:<220> log likelihood: -1896101\n",
      "INFO:lda:<230> log likelihood: -1896264\n",
      "INFO:lda:<240> log likelihood: -1896902\n",
      "INFO:lda:<250> log likelihood: -1895331\n",
      "INFO:lda:<260> log likelihood: -1896191\n",
      "INFO:lda:<270> log likelihood: -1896194\n",
      "INFO:lda:<280> log likelihood: -1895482\n",
      "INFO:lda:<290> log likelihood: -1895822\n",
      "INFO:lda:<300> log likelihood: -1895494\n",
      "INFO:lda:<310> log likelihood: -1895735\n",
      "INFO:lda:<320> log likelihood: -1895354\n",
      "INFO:lda:<330> log likelihood: -1895662\n",
      "INFO:lda:<340> log likelihood: -1895981\n",
      "INFO:lda:<350> log likelihood: -1895578\n",
      "INFO:lda:<360> log likelihood: -1895446\n",
      "INFO:lda:<370> log likelihood: -1896001\n",
      "INFO:lda:<380> log likelihood: -1895157\n",
      "INFO:lda:<390> log likelihood: -1895382\n",
      "INFO:lda:<400> log likelihood: -1895830\n",
      "INFO:lda:<410> log likelihood: -1894915\n",
      "INFO:lda:<420> log likelihood: -1895231\n",
      "INFO:lda:<430> log likelihood: -1895467\n",
      "INFO:lda:<440> log likelihood: -1895541\n",
      "INFO:lda:<450> log likelihood: -1895082\n",
      "INFO:lda:<460> log likelihood: -1893754\n",
      "INFO:lda:<470> log likelihood: -1895273\n",
      "INFO:lda:<480> log likelihood: -1894879\n",
      "INFO:lda:<490> log likelihood: -1895589\n",
      "INFO:lda:<500> log likelihood: -1894444\n",
      "INFO:lda:<510> log likelihood: -1895085\n",
      "INFO:lda:<520> log likelihood: -1894216\n",
      "INFO:lda:<530> log likelihood: -1895415\n",
      "INFO:lda:<540> log likelihood: -1894355\n",
      "INFO:lda:<550> log likelihood: -1895692\n",
      "INFO:lda:<560> log likelihood: -1895247\n",
      "INFO:lda:<570> log likelihood: -1895446\n",
      "INFO:lda:<580> log likelihood: -1895182\n",
      "INFO:lda:<590> log likelihood: -1895582\n",
      "INFO:lda:<600> log likelihood: -1895114\n",
      "INFO:lda:<610> log likelihood: -1895666\n",
      "INFO:lda:<620> log likelihood: -1895888\n",
      "INFO:lda:<630> log likelihood: -1895216\n",
      "INFO:lda:<640> log likelihood: -1894604\n",
      "INFO:lda:<650> log likelihood: -1894898\n",
      "INFO:lda:<660> log likelihood: -1896388\n",
      "INFO:lda:<670> log likelihood: -1895314\n",
      "INFO:lda:<680> log likelihood: -1896231\n",
      "INFO:lda:<690> log likelihood: -1894641\n",
      "INFO:lda:<700> log likelihood: -1895369\n",
      "INFO:lda:<710> log likelihood: -1895893\n",
      "INFO:lda:<720> log likelihood: -1895072\n",
      "INFO:lda:<730> log likelihood: -1894962\n",
      "INFO:lda:<740> log likelihood: -1895376\n",
      "INFO:lda:<750> log likelihood: -1894371\n",
      "INFO:lda:<760> log likelihood: -1895052\n",
      "INFO:lda:<770> log likelihood: -1894714\n",
      "INFO:lda:<780> log likelihood: -1895667\n",
      "INFO:lda:<790> log likelihood: -1894668\n",
      "INFO:lda:<800> log likelihood: -1895757\n",
      "INFO:lda:<810> log likelihood: -1895970\n",
      "INFO:lda:<820> log likelihood: -1895031\n",
      "INFO:lda:<830> log likelihood: -1895653\n",
      "INFO:lda:<840> log likelihood: -1895328\n",
      "INFO:lda:<850> log likelihood: -1895483\n",
      "INFO:lda:<860> log likelihood: -1894107\n",
      "INFO:lda:<870> log likelihood: -1895051\n",
      "INFO:lda:<880> log likelihood: -1894647\n",
      "INFO:lda:<890> log likelihood: -1895376\n",
      "INFO:lda:<900> log likelihood: -1894637\n",
      "INFO:lda:<910> log likelihood: -1894316\n",
      "INFO:lda:<920> log likelihood: -1894538\n",
      "INFO:lda:<930> log likelihood: -1894119\n",
      "INFO:lda:<940> log likelihood: -1892682\n",
      "INFO:lda:<950> log likelihood: -1893762\n",
      "INFO:lda:<960> log likelihood: -1893804\n",
      "INFO:lda:<970> log likelihood: -1893042\n",
      "INFO:lda:<980> log likelihood: -1892981\n",
      "INFO:lda:<990> log likelihood: -1892956\n",
      "INFO:lda:<999> log likelihood: -1893598\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Accuracy: 0.580962416578855\n",
      "             precision    recall  f1-score   support\n",
      "\n",
      "          0       0.59      0.65      0.62      1480\n",
      "          1       0.57      0.50      0.54      1367\n",
      "\n",
      "avg / total       0.58      0.58      0.58      2847\n",
      "\n",
      "Tuned Model Parameters: [ 0.35938137]\n",
      "AUC score when using topics: 0.604546847505882\n"
     ]
    }
   ],
   "source": [
    "# Select only documents which appear after 1947 (stored in data)\n",
    "X = [' '.join(doc) for doc in docsobj.stems]\n",
    "vectorizer = CountVectorizer()\n",
    "X = vectorizer.fit_transform(X)\n",
    "vocab = vectorizer.get_feature_names()\n",
    "\n",
    "# Split the X and y variables into a training and test set\n",
    "X_train, X_test, y_train, y_test = train_test_split(X.toarray(), y, test_size = 0.3, random_state = 42)\n",
    "\n",
    "# Create the logistic regression estimator with an l1 loss parameter\n",
    "log_reg = LogisticRegressionCV(penalty = \"l1\", solver=\"liblinear\")\n",
    "\n",
    "# Fit the model, predict and report the accuracy and the best parameter\n",
    "log_reg.fit(X_train, y_train)\n",
    "\n",
    "y_pred = log_reg.predict(X_test)\n",
    "\n",
    "print(\"Accuracy: {}\".format(log_reg.score(X_test, y_test)))\n",
    "\n",
    "print(classification_report(y_test, y_pred))\n",
    "\n",
    "print(\"Tuned Model Parameters: {}\".format(log_reg.C_))\n",
    "\n",
    "auc = roc_auc_score(y_test, log_reg.decision_function(X_test))\n",
    "print('AUC score when using document-term matrix: {}'.format(auc))\n",
    "\n",
    "## Run a logistic regression using the topics instead of the document-term matrix\n",
    "X = [' '.join(doc) for doc in docsobj.stems]\n",
    "vectorizer = CountVectorizer()\n",
    "X = vectorizer.fit_transform(X)\n",
    "\n",
    "K, S, alpha, eta = 10, 1000, 0.1, 0.01\n",
    "\n",
    "Col_Gibbs = lda.LDA(n_topics=K, n_iter=S, alpha=alpha, eta=eta)\n",
    "\n",
    "X = Col_Gibbs.fit_transform(X)\n",
    "\n",
    "# Split the X and y variables into a training and test set\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.3, random_state = 59)\n",
    "\n",
    "# Create the logistic regression estimator\n",
    "log_reg = LogisticRegressionCV(penalty=\"l2\", solver=\"newton-cg\")\n",
    "\n",
    "# Fit the model, predict and report the accuracy and the best parameter\n",
    "log_reg.fit(X_train, y_train)\n",
    "y_pred = log_reg.predict(X_test)\n",
    "print(\"Accuracy: {}\".format(log_reg.score(X_test, y_test)))\n",
    "print(classification_report(y_test, y_pred))\n",
    "print(\"Tuned Model Parameters: {}\".format(log_reg.C_))\n",
    "\n",
    "auc = roc_auc_score(y_test, log_reg.decision_function(X_test))\n",
    "print('AUC score when using topics: {}'.format(auc))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# We can see from the AUC scores that using document-term matrix is better for prediction as compared to using the topics from LDA"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}