{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import csv\n",
    "import dateutil.parser # for parsing plain-text dates\n",
    "import demjson # json (loose formatting)\n",
    "import json\n",
    "import math\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy\n",
    "import sklearn\n",
    "from collections import defaultdict\n",
    "from sklearn import linear_model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Data is available at http://cseweb.ucsd.edu/~jmcauley/pml/data/. Download and save to your own directory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataDir = \"/home/jmcauley/pml_data/\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Regression (lines of best fit)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read in a small dataset of reviews of fantasy novels (from Goodreads)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "path = dataDir + \"fantasy_100.json\"\n",
    "f = open(path)\n",
    "\n",
    "data = []\n",
    "\n",
    "for l in f:\n",
    "    d = json.loads(l)\n",
    "    data.append(d)\n",
    "    \n",
    "f.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What does the data look like? (json format)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'book_id': '18245960',\n",
       " 'date_added': 'Sun Jul 30 07:44:10 -0700 2017',\n",
       " 'date_updated': 'Wed Aug 30 00:00:26 -0700 2017',\n",
       " 'n_comments': 1,\n",
       " 'n_votes': 28,\n",
       " 'rating': 5,\n",
       " 'read_at': 'Sat Aug 26 12:05:52 -0700 2017',\n",
       " 'review_id': 'dfdbb7b0eb5a7e4c26d59a937e2e5feb',\n",
       " 'review_text': 'This is a special book. It started slow for about the first third, then in the middle third it started to get interesting, then the last third blew my mind. This is what I love about good science fiction - it pushes your thinking about where things can go. \\n It is a 2015 Hugo winner, and translated from its original Chinese, which made it interesting in just a different way from most things I\\'ve read. For instance the intermixing of Chinese revolutionary history - how they kept accusing people of being \"reactionaries\", etc. \\n It is a book about science, and aliens. The science described in the book is impressive - its a book grounded in physics and pretty accurate as far as I could tell. Though when it got to folding protons into 8 dimensions I think he was just making stuff up - interesting to think about though. \\n But what would happen if our SETI stations received a message - if we found someone was out there - and the person monitoring and answering the signal on our side was disillusioned? That part of the book was a bit dark - I would like to think human reaction to discovering alien civilization that is hostile would be more like Enders Game where we would band together. \\n I did like how the book unveiled the Trisolaran culture through the game. It was a smart way to build empathy with them and also understand what they\\'ve gone through across so many centuries. And who know a 3 body problem was an unsolvable math problem? But I still don\\'t get who made the game - maybe that will come in the next book. \\n I loved this quote: \\n \"In the long history of scientific progress, how many protons have been smashed apart in accelerators by physicists? How many neutrons and electrons? Probably no fewer than a hundred million. Every collision was probably the end of the civilizations and intelligences in a microcosmos. In fact, even in nature, the destruction of universes must be happening at every second--for example, through the decay of neutrons. Also, a high-energy cosmic ray entering the atmosphere may destroy thousands of such miniature universes....\"',\n",
       " 'started_at': 'Tue Aug 15 13:23:18 -0700 2017',\n",
       " 'user_id': '8842281e1d1347389f2ab93d60773d4d'}"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To predict ratings from review length, let's start by assembling vectors of both"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "ratings = [d['rating'] for d in data]\n",
    "lengths = [len(d['review_text']) for d in data]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And generating a quick scatter plot of the data to see overall trends..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAEZCAYAAACD/A7qAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X+cXHV97/HXOzPuBgbIIgSIwZBEqxaDF3Pb1D6wstELIrZOuBV/VB9gvddr1dtSvf0l1ospetX2VsVra/vQ1lIr/obsttUi1Sz4CxKJgeWHSAsUjDFEYUIYIGEmn/vH+e5ksjuzO/tjsrMn7+fjsY+c+Z7v+X4/3zOTz3znzJlzFBGYmVl+LJrvAMzMbG45sZuZ5YwTu5lZzjixm5nljBO7mVnOOLGbmeWME7vNKUm/Ielf5juO+SDphZLu7HIfmyW9sZt9TNL3AUmr56Nvmx4n9iOcpPskPSbpEUk/lvQpSUd3uO1p6T9743UUEVdFxHndi7h3RcS3IuLn5zuOudDmDcQ/elkgnNgtgJdHxHHAmcDzgXd2uK3S9upSbIedpNyMpQu8bxYIJ3aD9B82Ih4EriVL8NkK6XxJ2yTtkfQfki5r2u769G8lzfh/SdLFkr7ZtP0BSW+W9ENJD0n6WNO6RZL+XNJuSf8u6W3jPwE01f0DSV8cV3aFpI+k5TekNh5J/762o4FnM9P3SvqWpCqwStJxkv4mfYJ5QNLlyvRJeljS6U3bn5g+8Zwo6WxJDzStWybpS5IeTDH9dirvT9s8NT1+l6QnJR2THv+JpA91GP8bJd0h6WeSvippxWz3vaT3Ar8CfCztz482dXlOq/asx0SE/47gP+Be4MVp+VTgVuBDTetfBDw3La8BdgKvSI9PA+qAmupfDNzQ9PgAMAwcCzwdeBA4N637LeA2YBmwBLgutbeoRZwrgEeBUnq8CPgx8IvA0cAe4Jlp3cnAz3c4/s3AfcBzUptF4BrgL4HFwInAjcCbUv1PApc3bf9W4Ctp+Wzg/rQs4HvAu4ACsBL4N+CctH4EuCAtXwvcDbw0Pb4eKE8S7xvTchn4IfCsFPulwLfnYt8399NJe/7rrT/P2A1gk6RHgPuBXcB7xlZExA0RcXtavg34HFkCazbVR/T3R8TeiHiALGGMfSK4ELgiInZGxB7gA+0aiIj7gW3ABanoJUA1Iramx3XgDEmLI2JXREznS8y/i4gfRMQB4KnAy4C3R8QTEfFT4CPA2CeAzzYtA/wGcFWLNtcBJ0bE+yKiHhH3kb0pvCatvwE4W1IBeB7w0fS4n+zN6oYO4n4z2b79YYr9A8CZkp7eVGfW+36cdu1ZD3FiN8hmh8eRJeznkM1SAZC0TtI30uGEClkyObFNO+3salp+DDgmLT8NeKBpXfNyK81J9bWkhBoRjwGvBt4C7JT0j5KePY34mvs9DXhKauchSQ8Df8XBMW8GjpL0i5JOA/4T2Qx/vBXA8tTGWDvvBE5K668H1gNryT4lXQcMAi8A7o6IhzuI+zTgirE+gJ+RfeexvKnOXO37qdqzHuLEbnDwGPs3gSuBP29adxWwCVgeEQPAX3Nwhj7bsyR2kh3+GbOiXcXki8CgpOVkM/fGTDkirouIc4FTgLuAT0wjjuZxPAA8AZwQEU+NiOMjYiAinpf6OQB8gWym/lrgnyKi2qLNB4B7Uhtj7SyJiF9L678DPDuN4/qI+EEa//kc/O5iKg8Abx7XxzERcWMH2061730GzALmxG7jfYTsC7Iz0uNjgIcj4klJ68gS2pjdZMddnzHDvr4AXCLpaZIGgD+YrHI6LHI98CmypHkXgKSTJL0inab5JNmx+PpMAoqInwBfAz4s6dj0pelqSS9qqvZZsk8I7Q7DAGwB9qYvfRdLKkh6rqRfSP08DtwMvI2Difw7ZMe+O03sfwVcOvZlrqQlkl7Z4bZT7ftdgM9ZX6Cc2O2QmVlKnlcC/zsVvQ24XNIe4I+BzzfVfRx4H/DtdDhg3VTtj3v8CbIkeitZkvtnoJZmxe1cRXZ8/TNNZYuAdwA7gJ+SfeH7Fmj8aOiRSdprNTO9COgD7gAeIvukcEpjg4gtQJXsi8evtmw0G8Ovkh2Dvpfsi8ZPAMc1Vbue7IvVLU2Pj2Hy4+uNeCNiE9mx8c+lw2S3Aue1qtvi8VT7/grgwnS2zUc6aM96iCL83FhvkHQe8PGIWDXfsRxpvO/zxTN2mzfpEMXL0mGK5cBlwNXzHdeRwPs+3zxjt3kj6Siyww/PBh4H/gn43Yh4dF4DOwJ43+ebE7uZWc74UIyZWc4U5zsAAEn+2GBmNk0R0fJX3z0zY5/vaysczr/LLrts3mPwmD1mj3lhj3kyPZPYzcxsbjixm5nljBP7PBgcHJzvEA47j/nI4DH3hp443VFS9EIcZmYLhSSi1788NTOzueHEbmaWM07sZmY548RuZpYzTuxmZjnjxG5mljNO7GZmOePEbmaWM07sZmY548RuZpYzTuxmZjnjxG5mljNdv4OSpPuAPcAB4MmIWNftPs3MjmSH49Z4B4DBiHh4sko333wze/fuZceOHZxwwgmcfvrpFAoF+vr62LNnD0888USj7uLFi+nr6+Ohhx5i+fLlHH300VQqFQYGBgAay6VSaUI/u3fv5pZbbmHv3r2sWLGCFStWNNpfvHgxy5Ytm7BdtVrl3nvvpVqtsnr1apYuXTqh3Wq12rLf5vKpYutUu746sXv3bnbs2MHy5ctbjmM2fY3V7evrY//+/VNu06rt8WWTtTmb/TBXeiGGPOrktWHtdf2yvZLuBX4hIn42SZ14z3veM6t++vv72b9/P4sWLaJYLFKv1ymXy6xZs6ZR5ytf+Qpbt26dtJ1CocCGDRsa242OjnL11VcfUmfdunW87GUvazweHR1leHiYQqFwSL/N5VPF1ql2fXVi/PjHj2M2fY3VBajVahQKBSS13aZV2xFxSNnatWvZtm1byzbH153p/pyN2TwX1l4nrw3v6/m/bG8A10raKulN3epk3759RAT1ep19+/ZRq9UYGhqiWq0C2Ux1qqQOUK/XG9tVq1U2bdo0oc6WLVvYvXs3kM0ihoeHqdVqh/S7e/fuQ8oni61T7frqpJ1W428ex2z6aq5bq9WAbD+226ZV25s2bZpQtmXLlpZttqo7k/05G7N5Lqy9Tl8b3teTOxyHYs6KiJ2SlgLXSbozIr41vtLmzZsbyytXrmTVqlWz7rhQKFCpVCiVSuzYsaPj7SRRqVQay63s2LGDpUuXUqlUKBQKjeQz1u+OHTsmlLeLrVPt+uqknXbjHxvHbPpqVXeybVrVX7So8zlGq7oz2Z+zMZvnwtrr9LVxJO7rkZERRkZGOqrb9cQeETvTv7slXQOsAyYk9vXr18953/V6vXFse/ny5R1vFxGN7dodqhprb2BggHq9PqHf5cuXTyhvF1un2vXVSTvtxt+ufDp9tao72Tat6h84cKDtm+h4rerOZH/OxmyeC2uv09fGkbivBwcHD7kN38aNG9vW7eqhGElHSzomLZeAc4HbutFXf38/kigUCvT391MsFimXy4139KVLl7Ju3dQn5BQKhcZ2pVKJCy64YEKddevWNWa5pVKJcrlMsVg8pN+lS5ceUj5ZbJ1q11cn7bQaf/M4ZtNXc91iMZsrFAqFttu0anvDhg0TytatW9eyzVZ1Z7I/Z2M2z4W11+lrw/t6cl398lTSKuAasuPsReAzEfGBFvXie9/7ns+K6ZDPiumdsyR6IYY88lkxU5vsy1PfzNrMbAGa77NizMzsMHJiNzPLGSd2M7OccWI3M8sZJ3Yzs5xxYjczyxkndjOznHFiNzPLGSd2M7OccWI3M8sZJ3Yzs5xxYjczyxkndjOznHFiNzPLGSd2M7OccWI3M8sZJ3Yzs5xxYjczyxkndjOznHFiNzPLGSd2M7OccWI3M8sZJ3Yzs5xxYjczyxkndjOznHFiNzPLGSd2M7OccWI3M8sZJ3Yzs5xxYjczyxkndjOznHFiNzPLGSd2M7OccWI3M8sZJ3Yzs5w5LIld0iJJ2yQNH47+zMyOZIdrxn4JcMdh6svM7IhW7HYHkk4FzgfeB7yjXb3bbruNarXKsmXL2L9/PwB9fX3s3LmTUqnEySefzP79+xkYGKBUKs0olmq1SqVSOaSN5jKASqVCX19fy77G6rZb366vsXZb1Z9JvanGP52609m2k/F3q+9umY8+rbv8nB6GxA58GPh9YMlklb785S9P2VChUEAS5XKZNWvWTCuI0dFRhoeHKRQK1Ot1yuUyEdEoq9VqHDhwoLE8vq+x7YGW69v1NdZuX19fo9+x+s319u/fz6JFiygWi5PWG7+uk3F2uq8m27aT8Xer726Zjz6tu/ycZrp6KEbSy4FdEbEdUPqbsXq9Tq1WY2hoiGq12vF21WqV4eFharUa+/bta7QxNDTUKKvX60QEtVptQl+7d+9ubN9qfXMs4/saa7e532q1OqFeRFCv16esN9n4p1N3Ots2r2s3/m713S3z0ad1l5/Tg7o9Yz8LeIWk84GjgGMl/X1EXDS+4ubNmxvLK1euZNWqVW0bLRQKVCqVjj9mVSqVxux5jCQkUa/XJ922UCiwY8eOCdu3i6VVX63qjy3PpF678bfqu9N9Ndm2k8U6kzin03e3PkrPR5/WXXl/TkdGRhgZGemoblcTe0RcClwKIOls4H+1SuoA69ev77jder3eOCbdiYGBgQkJPCKIiI76Wr58eds3gPGxtOqrXf2Z1ms3/lZ9d7qvptq2k/F3q+9umI8+rbvy/pwODg4yODjYeLxx48a2dRfUeeyFQoFisUi5XJ7WO3CpVKJcLlMsFunv72+0sWHDhkbZ2DHjYrE4oa+lS5c2tm+1vjmW8X2Ntdvcb6lUmlBPEoVCYcp6k41/OnWns23zunbj71bf3TIffVp3+Tk9SJ3MWrsehBSjo6M+K6bDej4rZu74DIr8OVKeU0lERMvvLXsmsfdCHGZmC8VkiX1BHYoxM7OpObGbmeWME7uZWc44sZuZ5YwTu5lZzjixm5nljBO7mVnOOLGbmeWME7uZWc44sZuZ5YwTu5lZzjixm5nljBO7mVnOOLGbmeWME7uZWc44sZuZ5YwTu5lZzjixm5nljBO7mVnOOLGbmeWME7uZWc44sZuZ5YwTu5lZzjixm5nljBO7mVnOOLGbmeVMsZNKkj7aongP8L2IGJrbkMzMbDY6nbEvBs4E7k5/zwNOBf6bpI90KTYzM5sBRcTUlaQbgbMiop4eF4FvAi8ERiPi9FkFIUUncZiZWUYSEaFW6zqdsR8PHNP0uAQ8NSX6fbOMz8zM5lBHx9iBPwW2SxoBBLwI+D+SSsC/dik2MzObgY4OxQBIWgasSw+3RsSP5ywIH4oxM5uWuTgUM1Z3N/Aw8ExJL5qL4MzMbG51errjB4FXA7cDB1JxADd0KS4zM5uhTs+KuQt4XkR05YtSH4oxM5ueuTgUcw/wlBl03C/pJknflzQq6bLptmFmZtPT6Vkxj5GdFfN1mk5vjIjfmWyjiNgnaX1EPCapAHxb0lcjYku7barVKvfeey/VapXVq1ezdOlSqtUqlUqFgYEBSqVSo974sl61kGKdTF7GsRB4X/eWhfZ8dJrYh9PftEXEY2mxP/XX9pjL6OgoV1999SFlq1ev5v7776dQKFCv1ymXy0QEw8PDh5StWbNmJuF13ejo6IKJdTJ5GcdC4H3dWxbi89FRYo+IK2fagaRFwM3AM4C/iIitrepVq1U2bdo0ofyee+4BoFarATA0NEREUK/XDylbtWpVz72TVqtVhoeHqdVqPR/rZPIyjoXA+7q3LNTnY9LELukLEfEqSaO0mGlHxPOm6iAiDgDPl3QcsEnS6RFxx/h67373u9m+fTsHDhxg5cqVrFq1ql1MSKJerzfKCoUClUql53Z0pVKhUCg0XhDQu7FOJi/jWAi8r3tLLz0fIyMjjIyMdFR3qhn7JenfX51NQAAR8YikzcB5wITEfvnll/OhD32IAwcOTNz40HYYfwZNvV5nYGBgtiHOuYGBgUPegKB3Y51MXsaxEHhf95Zeej4GBwcZHBxsPN64cWPbupOeFRMRO9PiWyPiP5r/gLdOFYikEyUtSctHAecAP2hVt1QqccEFF0woX716NcVikf7+forFIuVymQ0bNkwo68XZTKlUolwuL4hYJ5OXcSwE3te9ZaE+H52ex74tItaOK7t1qkMxks4AriR7A1kEfD4i3teiXuM8dp8V07vyMo6FwPu6t/Ti8zHZeeyTJnZJbyGbma8G/r1p1bHAtyPi9XMUoH+gZGY2DbNJ7EvILtn7fuCPmlbtjYiH5jBAJ3Yzs2mYcWJv0dBJZHdTAiAi7p99eE7sZmbTNetLCkj6NUl3A/cC1wP3AV+dswjNzGzOdHqtmPcCLwB+GBGrgJcAN3YtKjMzm7FOE/uTEfEzYJGkRRGxGfiFLsZlZmYz1Om1YiqSjiG7/vpnJD0IVLsXlpmZzVSn57GXgMfJZvivA5YAn0mz+NkH4S9PzcymZc7OimlqcBHw2oj4zGyDS+05sZuZTcOMz4qRdJykd0r6mKRzlfmfZDfeeFU3gjUzs9mZ6gdKQ2Q3r/4u2ZkwJwECLomI7XMWhGfsZmbTMptfno5GxBlpuQDsBFZExBNzHKATu5nZNMzmB0pPji1ERB340VwndTMzm1tTzdjrHDytUcBRZPc/FRARcdycBOEZu5nZtEw2Y5/0PPaIKHQnJDMz65ZOf3lqZmYLhBO7mVnOOLGbmeWME7uZWc44sZuZ5YwTu5lZzjixm5nljBO7mVnOOLGbmeWME7uZWc44sZuZ5YwTu5lZzjixm5nljBO7mVnOOLGbmeWME7uZWc44sZuZ5YwTu5lZzjixm5nljBO7mVnOOLGbmeVMVxO7pFMlfUPS7ZJGJf1ON/szMzNQRHSvcekU4JSI2C7pGOBmoBwRPxhXL7oZh5lZ3kgiItRqXbGbHUfET4CfpOVHJd0JLAd+ML5utVpl586dPPHEEyxevJhly5ZRKpUa5QDLli0DoFKpMDAwQKlUmnFsze0uWbKE/fv3z0mbcxFbN/RybGY2t7qa2JtJWgmcCdzUav2HP/xh6vV643GhUGDt2rVs27atUS5lb059fX3U63XK5TJr1qyZdiyjo6MMDQ0d0l+xmO2K2bQ5PDxMoVCYVWzd0MuxmdncOyxfnqbDMF8CLomIR1vVaU6yY4+3bt16SHlEEBHs27ePWq3G0NAQ1Wp1WrFUq1WGh4cn9Fer1WbdZq1Wm1Vs3dDLsZlZd3R9xi6pSJbUPx0RQ+3qbd68ubG8cuVKVq1aNWXbhUKBSqUyrUMLlUqlMfOfyzYLhQK1Wm1W7XRDL8dmZp0bGRlhZGSko7qH41DM3wJ3RMQVk1Vav379tBuu1+sMDAxMa5uBgQEm+6J2pm22+sQx3Xa6oZdjM7PODQ4OMjg42Hi8cePGtnW7fbrjWcDrgBdL+r6kbZLOa1W3UChMeLxu3bpDyiUhif7+forFIuVyedqzzlKpRLlcntBfsVicdZvFYnFWsXVDL8dmZt3R1dMdOw5CikcffdRnxXRRL8dmZtM32emOPZPYeyEOM7OFYrLE7ksKmJnljBO7mVnOOLGbmeWME7uZWc44sZuZ5YwTu5lZzjixm5nljBO7mVnOOLGbmeWME7uZWc44sZuZ5YwTu5lZzjixm5nljBO7mVnOOLGbmeWME7uZWc44sZuZ5YwTu5lZzjixm5nljBO7mVnOOLGbmeWME7uZWc44sZuZ5YwTu5lZzjixm5nljBO7mVnOOLGbmeWME7uZWc44sZuZ5YwTu5lZzjixm5nljBO7mVnOOLGbmeWME7uZWc44sZuZ5UxXE7ukv5G0S9Kt3ezHzMwO6vaM/VPASzupuGPHDqrVatv11Wp1yjqz0e32zcwOl2I3G4+Ib0k6rZO6n/70p6nX65TLZdasWXPIutHRUYaHhykUCm3rzEa32zczO5x65hj7vn37qNVqDA0NHTJrrlarDA8PU6vV2taZjW63b2Z2uHV1xj4dmzdvBqBYLPKsZz2LCy+8EIBKpUKhUKBWqzXqFgoFKpUKpVJp1v12u30zs7kwMjLCyMhIR3V7JrGvX78eyBL7+eef3ygfGBigXq8fUrderzMwMDAn/Xa7fTOzuTA4OMjg4GDj8caNG9vWPRyHYpT+JtXf30+xWKRcLh8yUy6VSpTLZYrFYts6s9Ht9s3MDjdFRPcal64CBoETgF3AZRHxqRb14kc/+hEDAwNtE2q1WqVSqUxaZza63b6Z2VySRES0nDR3NbF3SlL0QhxmZgvFZIm9Z86KMTOzueHEbmaWM07sZmY548RuZpYzTuxmZjnjxG5mljNO7GZmOePEbmaWM07sZmY548RuZpYzTuxmZjnjxG5mljNO7GZmOePEbmaWM07s86DT21vlicd8ZPCYe4MT+zzoxRdCt3nMRwaPuTc4sZuZ5YwTu5lZzvTMrfHmOwYzs4Wmp+95amZmc8eHYszMcsaJ3cwsZ+Y1sUs6T9IPJP1Q0h/OZyyzJelvJO2SdGtT2fGSvibpLknXSlrStO6jku6WtF3SmU3lF6f9cZekiw73OKZD0qmSviHpdkmjkn4nledy3JL6Jd0k6ftpvJel8pWSbkzxf1ZSMZX3SfpcGu93Ja1oauudqfxOSefO15g6JWmRpG2ShtPjXI9Z0n2SbknP9ZZUtnBe1xExL39kbyr/BpwGPAXYDjxnvuKZg/G8EDgTuLWp7IPAH6TlPwQ+kJZfBvxzWv4l4Ma0fDzw78ASYGBseb7HNsmYTwHOTMvHAHcBz8nzuIGj078F4MY0js8DF6byjwNvTstvAf4yLb8a+FxaPh34PlAEVqb/B5rvsU0x7rcD/wAMp8e5HjNwD3D8uLIF87qezxn7OuDuiPiPiHgS+BxQnsd4ZiUivgU8PK64DFyZlq/k4PjKwN+n7W4Clkg6GXgp8LWI2BMRFeBrwHndjn2mIuInEbE9LT8K3AmcSo7HHRGPpcV+siQVwHrgy6n8SmBDWm7eD18CXpyWX0GW8GoRcR9wN9n/h54k6VTgfOCTTcUvJsdjBsTEIxoL5nU9n4l9OfBA0+MfpbI8OSkidkGWBIGTU3m7sY8v38EC2SeSVpJ9YrkRODmv406HJL4P/AS4jmwWVomIA6lK8+u4Ma6IqAN7JD2VBTTe5MPA75O9iSHpBODhnI85gGslbZX031PZgnldFw9HJ9bQ7tzSlueiLhSSjiGbnV0SEY+2+F1CbsadktnzJR0HXEN26KlTC268kl4O7IqI7ZIGm1d12sTcR3VYnBUROyUtBb4m6S4mvo579nU9nzP2HcCKpsenprI82ZU+kiHpFODBVL4DeHpTvbGxL7h9kr40+xLw6YgYSsW5H3dEPAKMAL8MDEga+7/UHHtjvJIKwHER8RDt90MvOgt4haR7gM+SHVq5guxwQ17HTETsTP/uBjaRHTZaMK/r+UzsW4FnSjpNUh/wGmB4HuOZC+LQd+th4A1p+Q3AUFP5RQCSXkD2UX4XcC1wjqQlko4HzkllvexvgTsi4oqmslyOW9KJY2dCSDqKLM47gM3AhanaxRw63ovT8oXAN5rKX5POIFkFPBPY0v0RTF9EXBoRKyJiNdn/0W9ExOvJ8ZglHZ0+hSKpBJwLjLKQXtfz/M3zeWRnUtwN/NF8xjIHY7kK+DGwD7gf+E2yb8X/NY3xa8BAU/2PkZ0ZcAuwtqn8DWl//BC4aL7HNcWYzwLqZGc0fR/Ylp7Tp+Zx3MAZaYzbgVuBd6XyVcBNKfbPA09J5f3AF9K4bgRWNrX1zrQf7gTOne+xdTj+szl4Vkxux5zGNvaaHh3LTQvpde1LCpiZ5Yx/eWpmljNO7GZmOePEbmaWM07sZmY548RuZpYzTuxmZjnjxG4TSKqnS7SOShpKP5+fSTvLJH1hjmO7N117ZC7bPE3Sa5seXyzp/3W47RfTdXLarb9M0jtmH2Vnxo9lDtr7rKRnzFV7dng4sVsr1YhYGxFnkF2x8m0zaSQidkbEq+Y2tLbX55iNVcBvTLcfSacDiyK7WmFXNP1sv1OtxjJVH4VJVn+c7BK1toA4sdtUvkvTFekk/Z6kLemGAmM3mni/pLc21blM0jvS7HE0lS2S9KfKblSxXdKbUvnHJP1qWr5G0ifT8m9KurxFPGrq53WpvW2SPi5JqXyvpPemfr6TLuSEpNXKbv5wi6TLJe1NTb0feGFq55JUtlzSV9MNEj7YZt+8joM/Kx+7cczNqd/rmuo9V9JmSf8m6beb6l+Trh442nQFwbH4/2+6iuQLJL077fNbJf1VU71nSLou9fc9SavHj2WS/X62pBskDQG3p5/R/5OyG0vcKmnscgHfBP7LDN5gbD7N9893/dd7f8De9G+B7Ofh56bH5wB/nZYF/CMHbzAy0rT97cDTyG6icmsqexNwaVruI7tW0GlkN2P4YCq/CfhOWv5b4JwWsd1L9tPu55Bdo6OQyv8CeH1aPgCcn5Y/2NTvPwKvSstvBh5Jy42fyqfHF5P9PPwYsp/I3wcsbxHLCPDctHwi2aUkVqTHA+nfy4BvkV1J9QTgp00xj9VZTPbT9eOb4v/1pn6af7r+98DL0/KNwCua9uniFmNpt9/PBvY2xftfx57b9PjYpuVrgefP9+vSf53/+V3YWjlK0jZgJ3AS2XXHIbsY0jlp3Tbg2cDPRXazjaWSTpH0POChiPjxuDbPBS5Ks9CbyJLzz5HNCF8k6efJLqi1S9mV834Z+E6L2MYOkbwEWAtsTW2+mOwwBMD+iPhKWr6Z7I49pDa/lJavmmIffD0iHo2IfSmu01rUWQbsTssvAK6PiPsBIruxwph/juwGEz8DdnHwOt6/K2k7WYI+lWx/ANSAq5u2f4my29DdSnZTj+emi1Q9LSKGU3/7I+KJFjG22+8AW8biJXtjOSd9+nphROxtamM32Ru1LRC+Hru18lhErJW0mGy29jayixwJeH9EfKLFNl8ku5rfKWQXhRpPwG9HxHUTVkgDZHebuZ4s8byK7FNDdZIYBVwZEe9qsW5/03Kdg6/z5uPmU10ze1+bNpoe4JP/AAACTElEQVQ9RjZLnqrN5rYOAEVJZ5O9Gf1SROyTtLmprSciYuymFv1kn0bWRsSP0+GvsXqdXPe75X5P/Tf2b0TcLWkt2Z2S3ivp6xExdihsMfB4B31Zj/CM3VoRQJoBXgL8XjrGei3wRmWXMkXS08aOX5MdsnkN8OtkSX68a4G36uBNj39O2aVvIZuxvh24geywxe+RzeTbxgZ8HXhl0/Hz4yU9fVyd8W4EXpmWX9NUvhc4ts02k7mT7PKzY23/iqTTxuKZYtslZHch2ifpOWQz/jHN8S8me0P6WZqlvxIatyJ8QFI59deX9uf4sbTa70ePD0bSMuDxiLgK+DPg+U2rnwXcNsV4rId4xm6tNGa2kd055xbgtRHxmXTI5Lvpe8q9wOuB3RFxh6RjgR9Fun3YOJ8kOySyLX3J+SAH75P5TbLj6fdIup/scsc3TBZbRNwp6Y/J7m6ziGyW/jayW5G1O6Pl7cA/SLqULOHtSeW3AgfS4Yq/Y+K9a9u19xWyQyPfiIifSvofwDVN43tpu/iBfwF+S9LtZJeB/W6r/iJij6RPkH1vsZNDr2F+EfDXkv6EbPwXjh9LRFyh7HTMVvu92RnAn0k6kNp6C4Ckk8g+wT3YYhvrUb5srx0xJB0VEY+n5VcDr4mIC2bR3mKyG0mcFTn9jyTpd4E9EfGp+Y7FOucZux1J/rOkse8KHgbeOJvGIuKJdMx7OdkNjPPoYeDT8x2ETY9n7GZmOeMvT83McsaJ3cwsZ5zYzcxyxondzCxnnNjNzHLGid3MLGf+Pxga6gxxuDTfAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fed881013c8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(lengths, ratings, color='grey')\n",
    "plt.xlim(0, 5500)\n",
    "plt.ylim(0.5, 5.5)\n",
    "plt.xlabel(\"Review length (characters)\")\n",
    "plt.ylabel(\"Rating\")\n",
    "plt.title(\"Rating vs. review length\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To perform regression, convert to features (X) and labels (y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = numpy.matrix([[1,l] for l in lengths]) # Note the inclusion of the constant term\n",
    "y = numpy.matrix(ratings).T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fit a linear regression model to the data (sklearn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LinearRegression(copy_X=True, fit_intercept=False, n_jobs=None, normalize=False)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model = sklearn.linear_model.LinearRegression(fit_intercept=False)\n",
    "model.fit(X, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Extract the model coefficients (theta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[3.98394783e+00, 1.19363599e-04]])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "theta = model.coef_\n",
    "theta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Same thing using numpy (lstsq)..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[3.98394783e+00],\n",
       "        [1.19363599e-04]])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "theta,residuals,rank,s = numpy.linalg.lstsq(X, y, rcond=None)\n",
    "theta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Same thing using by computing the pseudoinverse directly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[3.98394783e+00],\n",
       "        [1.19363599e-04]])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "numpy.linalg.inv(X.T*X)*X.T*y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the line of best fit..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "xplot = numpy.arange(0,5501,10)\n",
    "yplot = [(theta[0] + theta[1]*x).item() for x in xplot]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAEZCAYAAACD/A7qAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt8VPWd//HXhwkJNiBBSBBRIFjxBipYadUqURcvaA1eqtVWrVbXVrd1161d29pVWltbf9aq28t2XWut9dIuVUnVFisQbhYFEY2ClxbxEilBJSAREjL5/P44J+MkM5NMLkMmh/fz8ciDc77zPd/v53tm+MyZ75w5x9wdERGJjgF9HYCIiPQuJXYRkYhRYhcRiRgldhGRiFFiFxGJGCV2EZGIUWKXXmVm55vZn/s6jr5gZp82szU57mOBmV2Syz466LvFzMb3Rd/SNUrsuzgzW2dmH5rZFjN7x8zuNrOPZbnt2PA/e+J15O73u/vJuYs4f7n7Enc/sK/j6A0Z3kD0o5d+QoldHDjV3XcHDgMmA9/MclsLt7ccxbbTmVlkxpID2jf9hBK7QPgf1t3rgLkECT54wGyGma00s81m9oaZXZ+03cLw3/rwiP+TZnaRmS1O2r7FzC43s1fN7H0z+2nSYwPM7MdmttHM/m5mV7b/BJBU9xtm9n/tym43s9vC5S+GbWwJ/z0vq4EHR6Y3mtkSM2sAys1sdzO7K/wE85aZfc8ChWa2ycwOStp+RPiJZ4SZTTOzt5IeG2Vms82sLozpq2F5UbjNHuH6t81sh5kNDte/a2a3Zhn/JWa22szeM7M/mdmYnu57M7sROAb4abg/70jqcnq69iTPuLv+duE/4HXg+HB5b+AF4Nakx48FDg6XJwLrgdPD9bFAHLCk+hcBi5LWW4AqYAiwD1AHnBg+9mXgRWAUMBT4S9jegDRxjgG2AsXh+gDgHeAI4GPAZuDj4WMjgQOzHP8CYB1wQNhmAfAw8HNgEDACWAZcFtb/X+B7SdtfATweLk8D3gyXDVgBfBuIAeOAvwHTw8ergTPC5bnAa8BJ4fpCoLKDeC8JlyuBV4EJYezfApb2xr5P7ieb9vSXX386YheAR8xsC/AmsAG4ofUBd1/k7i+Fyy8CDxIksGSdfUS/yd0/cPe3CBJG6yeCzwK3u/t6d98M/DBTA+7+JrASOCMsOgFocPfl4XocmGRmg9x9g7t35UvMX7v7y+7eAuwBnAL8m7tvd/d3gduA1k8ADyQtA5wP3J+mzanACHf/vrvH3X0dwZvC58LHFwHTzCwGHALcEa4XEbxZLcoi7ssJ9u2rYew/BA4zs32S6vR437eTqT3JI0rsAsHR4e4ECfsAgqNUAMxsqpnND6cT6gmSyYgM7WSyIWn5Q2BwuLwX8FbSY8nL6SQn1fMIE6q7fwicC3wFWG9mfzSz/bsQX3K/Y4GBYTvvm9km4L/5aMwLgN3M7AgzGwscSnCE394YYHTYRms73wTKwscXAscBUwg+Jf0FqAA+Bbzm7puyiHsscHtrH8B7BN95jE6q01v7vrP2JI8osQt8NMe+GLgH+HHSY/cDjwCj3b0E+CUfHaH39CyJ9QTTP63GZKoY+j+gwsxGExy5J46U3f0v7n4isCfwCnBnF+JIHsdbwHZguLvv4e7D3L3E3Q8J+2kBfk9wpH4e8Ki7N6Rp8y1gbdhGaztD3f0z4eNPAfuH41jo7i+H45/BR99ddOYt4PJ2fQx292VZbNvZvtcZMP2YEru0dxvBF2STwvXBwCZ332FmUwkSWquNBPOu+3azr98DV5nZXmZWAnyjo8rhtMhC4G6CpPkKgJmVmdnp4WmaOwjm4uPdCcjd/wE8AfzEzIaEX5qON7Njk6o9QPAJIdM0DMAzwAfhl76DzCxmZgeb2SfCfrYBzwJX8lEif4pg7jvbxP7fwLdav8w1s6FmdnaW23a27zcAOme9n1JilzZHZmHyvAf4z7DoSuB7ZrYZuA74XVLdbcD3gaXhdMDUztpvt34nQRJ9gSDJPQY0h0fFmdxPML9+X1LZAOBqoBZ4l+AL369A4kdDWzpoL92R6YVAIbAaeJ/gk8KeiQ3cnwEaCL54/FPaRoMxnEYwB/06wReNdwK7J1VbSPDF6jNJ64PpeH49Ea+7P0IwN/5gOE32AnByurpp1jvb97cDnw3Ptrkti/Ykj5i7nhvJD2Z2MvALdy/v61h2Ndr30aIjdukz4RTFKeE0xWjgeuChvo5rV6B9H206Ypc+Y2a7EUw/7A9sAx4F/tXdt/ZpYLsA7ftoU2IXEYkYTcWIiERMQV8HAGBm+tggItJF7p72V995c8Te19dW2Jl/119/fZ/HoDFrzBpz/x5zR/ImsYuISO9QYhcRiRgl9j5QUVHR1yHsdBrzrkFjzg95cbqjmXk+xCEi0l+YGZ7vX56KiEjvUGIXEYkYJXYRkYhRYhcRiRgldhGRiFFiFxGJGCV2EZGIUWIXEYkYJXYRkYhRYhcRiRgldhGRiFFiFxGJmJzfQcnM1gGbgRZgh7tPzXWfIiK7sp1xa7wWoMLdN3VU6dlnn+WDDz6gtraW4cOHc9BBBxGLxSgsLGTz5s1s3749UXfQoEEUFhby/vvvM3r0aD72sY9RX19PSUkJQGK5uLg4pZ+NGzfy/PPP88EHHzBmzBjGjBmTaH/QoEGMGjUqZbuGhgZef/11GhoaGD9+PKWlpSntNjQ0pO03ubyz2LKVqa9sbNy4kdraWkaPHp12HD3pq7VuYWEhTU1NnW6Tru32ZR212ZP90FvyIYYoyua1IZnl/LK9ZvY68Al3f6+DOn7DDTf0qJ+ioiKampoYMGAABQUFxONxKisrmThxYqLO448/zvLlyztsJxaLMXPmzMR2NTU1PPTQQ23qTJ06lVNOOSWxXlNTQ1VVFbFYrE2/yeWdxZatTH1lo/3424+jJ3211gVobm4mFothZhm3Sde2u7cpmzJlCitXrkzbZvu63d2fPdGT50Iyy+a1oX3d95ftdWCumS03s8ty1UljYyPuTjwep7GxkebmZubMmUNDQwMQHKl2ltQB4vF4YruGhgYeeeSRlDrPPPMMGzduBIKjiKqqKpqbm9v0u3HjxjblHcWWrUx9ZdNOuvEnj6MnfSXXbW5uBoL9mGmbdG0/8sgjKWXPPPNM2jbT1e3O/uyJnjwXklm2rw3t647tjKmYo919vZmVAn8xszXuvqR9pQULFiSWx40bR3l5eY87jsVi1NfXU1xcTG1tbdbbmRn19fWJ5XRqa2spLS2lvr6eWCyWSD6t/dbW1qaUZ4otW5n6yqadTONvHUdP+kpXt6Nt0tUfMCD7Y4x0dbuzP3uiJ8+FZJbta2NX3NfV1dVUV1dnVTfnid3d14f/bjSzh4GpQEpiP+6443q973g8npjbHj16dNbbuXtiu0xTVa3tlZSUEI/HU/odPXp0Snmm2LKVqa9s2sk0/kzlXekrXd2OtklXv6WlJeObaHvp6nZnf/ZET54LySzb18auuK8rKira3IZv1qxZGevmdCrGzD5mZoPD5WLgRODFXPRVVFSEmRGLxSgqKqKgoIDKysrEO3ppaSlTp3Z+Qk4sFktsV1xczBlnnJFSZ+rUqYmj3OLiYiorKykoKGjTb2lpaZvyjmLLVqa+smkn3fiTx9GTvpLrFhQExwqxWCzjNunanjlzZkrZ1KlT07aZrm539mdP9OS5kMyyfW1oX3csp1+emlk58DDBPHsBcJ+7/zBNPV+xYoXOismSzorJn7Mk8iGGKNJZMZ3r6MtT3cxaRKQf6uuzYkREZCdSYhcRiRgldhGRiFFiFxGJGCV2EZGIUWIXEYkYJXYRkYhRYhcRiRgldhGRiFFiFxGJGCV2EZGI2RnXYxcRkU488MADLFmyhLq6OjZs2EBdXR11dXXcddddaa8y2xEldhGRHHj00UdZsmRJmyRdV1fHzTffzLnnnptSv6ioiAMOOIBp06ZRVlbGyJEjKSsrY9iwYV3uW4ldRCQLCxcuZPHixYkE3Zqwr732Wi644IKU+o2NjQwdOpT99tsvkaTLysoYNWpU2vbPPPPMXotVl+0VkV3SihUrWLp0aUqivvzyy7n44otT6s+ePZuVK1e2SdIjR45k7NixDBkyZKfHr+uxi0jkrVmzhr/+9a8pifrcc8/l0ksvTak/e/ZsFi5cmJKo99tvP0aMGNEHI+gaJXYR6XfWrVvH8uXLUxL1SSedxOWXX55S/6GHHqKqqiolUR944IHss88+fTCC3FJiF5E+949//INVq1alJOojjzySL3/5yyn1H3nkEe69996URD1p0iQmTJjQByPIL0rsItLrNm3axOrVq1MS9aRJk9IeUT/66KPccccdKYn60EMPZfLkyX0wgv5NiV1EOvXhhx/y97//PZGgW/8dN25c2kT95z//mVmzZqUk6sMOO4xjjjmmD0awa1FiF9kFNTU1UVtbm/KDl9LSUi677LKU+vPnz+drX/ta2kQ9Y8aMPhiBdKSjxK7z2EX6iZaWFt5///2UI+rBgwdzySWXpNRfvnw5559/fpsfu4wcOZLx48enbf/444/nxRdfzPUwZCdQYhfpQ9u3b0+Zo47FYlx44YUpdZ977jlOPPHElEQ9adKktG0fffTRvPHGG7keguQhJXaRXuTu1NfXt0nUzc3NfO5zn0up+/LLL3PIIYekJOqDDjoobduHH3447733Xq6HIBGgOXaRTuzYsYONGzcmEnVDQ0Pan3/X1tZSXl7Obrvt1iZRT5gwgZtuuimlfutr3iztNKlIh/TlqUgSd2fr1q2JRF1fX5/2y8H6+nr23XdftmzZwogRIxKJety4cfzP//xPSv2WlhaampoYNGjQzhiG7OKU2CXy4vE47777LnV1dbz77rscd9xxKXWampqYMGECGzZsIBaLJRL1XnvtxR/+8IeU+u7Oe++9xx577MGAAbp1geQXJXbplz788MPEF4tHHHFEypSFuzN58mTeeecdNm3aRElJSWIK5IknnqCgIPUrpNdff52ysjKKi4t31jBEckKJXfJCS0sLmzZtoq6ujv322y9t4p0+fTpr166lrq6OHTt2JBL1ggULGDx4cEr9l156ieHDhzNixIi07YlElRK75ExjYyN1dXXsueeeDBw4MOXx8847jzVr1iSmSAYPHkxZWRmLFy+mtLQ0pf7y5csTR95DhgzRF4siGSixS9bcnc2bN1NcXJw2UV955ZU8//zzidP5tm3bRllZGUuWLGHcuHEp9RcvXkxxcTEjR46ktLSUwsLCnTAKkehTYt/F7dixAzNLO1Vx3XXX8eyzz7a5ddegQYNYunQpEydOTKk/b948Bg4cmJgiKSkp0VG1SB9QYo+Yjs5//vGPf8zTTz/d5teMW7ZsYdGiRRx55JEp9f/0pz/R0tKSSNSlpaXstttuOR+DiPSMEns/dtddd/HUU0+l/Oz84Ycf5qSTTkqpX1VVxYcfftjmBzI6XU8kepTY88js2bNZsmRJSqL+5S9/ycyZM1PqP/TQQ7z33nspV9zT6Xoiu7Y+T+xmNgBYAbzt7qenebzfJva5c+emTdQ/+MEPOOecc1Lqz549mzfffDMlUet0PRHpiny4bO9VwGpg953UX7ctXbqUp556KuXSqF//+tf5/Oc/n1J/8+bNFBQUMHny5DaJeu+9907b/tlnn53rIYjILi7nR+xmtjdwN/B94OpMR+w1NTU0NDQwatQompqaACgsLGT9+vWJ0+WampooKSnp0jTEqlWrWLZsGRs2bKC2tpba2lrq6+u57LLL+OIXv0hDQwP19fWUlJQA8Jvf/Iaamhr22GMPxowZwz777JO4hvXAgQOpr6+nsLCw01jat9u63L5+d+p1Nv6u1O3Ktq2PdTT+XPWdK33Rp+TWrvKc9vUR+0+Aa4ChHVVKd62O9mKxGO+99x5lZWUUFha2OaI+66yz+NKXvpSyzauvvspzzz2HmbF+/XqGDh3KqFGjGDduHDU1NVRVVRGLxWhubqalpYVYLMbIkSOJxWJs3LiRo446iokTJybqAjQ3NxOLxTAzKisrU04LTNduYWEh8Xi8Tf3kek1NTQwYMICCgoIO67V/rL2u1O3KttmMP1d950pf9Cm5pec0kNPEbmanAhvcfZWZVQAZT3hesWIFDQ0Nib/x48dz+OGHt6kTj8epra1l/vz5HHPMMey9997sv//+jBw5koMPPjhtu+eccw6nnnoqt912G6NGjUqUL126FHcnHo/T3NycKG9djsfjAMyZM4eRI0dSVVXVpl7y4+Xl5Ykjg4aGhkTd5PqNjY1t6gMp9eLxeEq76eq177NVur4z1W2vo22TY8g0/q7E2ZW+c3XE1Rd9Sm7pOf1Iro/YjwZON7MZwG7AEDP7jbun3B5mxYoVFBYWMnDgQEaPHs0+++yTtsGDDz6YKVOmcMEFFzB69Oisgqivr08cPbcyM8wskaAyicVi1NbWpmyf/Hh9fX3ihZOur3T1W5e7U699nx2NM1Pd9jratqNYuxNnV/rO1X/IvuhTcivqz2l1dTXV1dVZ1c1pYnf3bwHfAjCzacC/p0vqAF/+8pezbjcejyfmpLNRUlKSksDdnWy+X4jH44wePTrjG0D7WNL1lal+d+tlGn+6vrPdV51tm834c9V3LvRFn5JbUX9OKyoqqKioSKzPmjUrY91+9auVWCxGQUEBlZWVXXoHLi4uprKykoKCAoqKihJtzJw5M1HWOmfcesphcl+lpaWJ7dM9nhxL+75a203ut7i4OKWemRGLxTqt19H4u1K3K9smP5Zp/LnqO1f6ok/JLT2nH8mbHyjl6qyYZOm+LU93Vkqmsz6yOSskXV+t7eqsmO73nSu7yhkUu5Jd5Tnt8x8odaY//0BJRKQvdJTY+9VUjIiIdE6JXUQkYpTYRUQiRoldRCRilNhFRCJGiV1EJGKU2EVEIkaJXUQkYpTYRUQiRoldRCRilNhFRCJGiV1EJGKU2EVEIkaJXUQkYpTYRUQiRoldRCRilNhFRCJGiV1EJGKU2EVEIkaJXUQkYpTYRUQiRoldRCRilNhFRCJGiV1EJGKU2EVEIkaJXUQkYgqyqWRmd6Qp3gyscPc5vRuSiIj0RLZH7IOAw4DXwr9DgL2BL5nZbTmKTUREusHcvfNKZsuAo909Hq4XAIuBTwM17n5Qj4Iw82ziEBGRgJnh7pbusWyP2IcBg5PWi4E9wkTf2MP4RESkF2U1xw7cDKwys2rAgGOBH5hZMfBkjmITEZFuyGoqBsDMRgFTw9Xl7v5OrwWhqRgRkS7pjamY1robgU3Ax83s2N4ITkREele2pzv+CDgXeAloCYsdWJSjuEREpJuyPSvmFeAQd8/JF6WaihER6ZremIpZCwzsRsdFZva0mT1nZjVmdn1X2xARka7J9qyYDwnOiplH0umN7v61jjZy90YzO87dPzSzGLDUzP7k7s9k2qahoYHXX3+dhoYGxo8fT2lpKQ0NDdTX11NSUkJxcXGiXvuyfNWfYu1IVMbRH2hf55f+9nxkm9irwr8uc/cPw8WisL+Mcy41NTU89NBDbcrGjx/Pm2++SSwWIx6PU1lZibtTVVXVpmzixIndCS/nampq+k2sHYnKOPoD7ev80h+fj6wSu7vf090OzGwA8CywL/Azd1+erl5DQwOPPPJISvnatWsBaG5uBmDOnDm4O/F4vE1ZeXl53r2TNjQ0UFVVRXNzc97H2pGojKM/0L7OL/31+egwsZvZ7939HDOrIc2Rtrsf0lkH7t4CTDaz3YFHzOwgd1/dvt53vvMdVq1aRUtLC+PGjaO8vDxTTJgZ8Xg8URaLxaivr8+7HV1fX08sFku8ICB/Y+1IVMbRH2hf55d8ej6qq6uprq7Oqm5nR+xXhf+e1pOAANx9i5ktAE4GUhL79773PW699VZaWlpSN27bDu3PoInH45SUlPQ0xF5XUlLS5g0I8jfWjkRlHP2B9nV+yafno6KigoqKisT6rFmzMtbt8KwYd18fLl7h7m8k/wFXdBaImY0ws6Hh8m7AdODldHWLi4s544wzUsrHjx9PQUEBRUVFFBQUUFlZycyZM1PK8vFopri4mMrKyn4Ra0eiMo7+QPs6v/TX5yPb89hXuvuUdmUvdDYVY2aTgHsI3kAGAL9z9++nqZc4j11nxeSvqIyjP9C+zi/5+Hx0dB57h4ndzL5CcGQ+Hvh70kNDgKXu/oVeClA/UBIR6YKeJPahBJfsvQm4NumhD9z9/V4MUIldRKQLup3Y0zRURnA3JQDc/c2eh6fELiLSVT2+pICZfcbMXgNeBxYC64A/9VqEIiLSa7K9VsyNwKeAV929HDgBWJazqEREpNuyTew73P09YICZDXD3BcAnchiXiIh0U7bXiqk3s8EE11+/z8zqgIbchSUiIt2V7XnsxcA2giP8zwNDgfvCo/ieB6EvT0VEuqTXzopJanAAcJ6739fT4ML2lNhFRLqg22fFmNnuZvZNM/upmZ1ogX8huPHGObkIVkREeqazHyjNIbh59V8JzoQpAwy4yt1X9VoQOmIXEemSnvzytMbdJ4XLMWA9MMbdt/dygErsIiJd0JMfKO1oXXD3OPB2byd1ERHpXZ0dscf56LRGA3YjuP+pAe7uu/dKEDpiFxHpko6O2Ds8j93dY7kJSUREciXbX56KiEg/ocQuIhIxSuwiIhGjxC4iEjFK7CIiEaPELiISMUrsIiIRo8QuIhIxSuwiIhGjxC4iEjFK7CIiEaPELiISMUrsIiIRo8QuIhIxSuwiIhGjxC4iEjFK7CIiEaPELiISMUrsIiIRo8QuIhIxSuwiIhGT08RuZnub2Xwze8nMaszsa7nsT0REwNw9d42b7Qns6e6rzGww8CxQ6e4vt6vnuYxDRCRqzAx3t3SPFeSyY3f/B/CPcHmrma0BRgMvt6/b0NDA+vXr2b59O4MGDWLUqFEUFxcnygFGjRoFQH19PSUlJRQXF3c7tuR2hw4dSlNTU6+02Rux5UI+xyYivSuniT2ZmY0DDgOeTvf4T37yE+LxeGI9FosxZcoUVq5cmSg3C96cCgsLicfjVFZWMnHixC7HUlNTw5w5c9r0V1AQ7IqetFlVVUUsFutRbLmQz7GJSO/bKV+ehtMws4Gr3H1rujrJSbZ1ffny5W3K3R13p7GxkebmZubMmUNDQ0OXYmloaKCqqiqlv+bm5h632dzc3KPYciGfYxOR3Mj5EbuZFRAk9XvdfU6megsWLEgsjxs3jvLy8k7bjsVi1NfXd2lqob6+PnHk35ttxmIxmpube9ROLuRzbCKSverqaqqrq7OquzOmYn4FrHb32zuqdNxxx3W54Xg8TklJSZe2KSkpoaMvarvbZrpPHF1tJxfyOTYRyV5FRQUVFRWJ9VmzZmWsm+vTHY8GPg8cb2bPmdlKMzs5Xd1YLJayPnXq1DblZoaZUVRUREFBAZWVlV0+6iwuLqaysjKlv4KCgh63WVBQ0KPYciGfYxOR3Mjp6Y5ZB2HmW7du1VkxOZTPsYlI13V0umPeJPZ8iENEpL/oKLHrkgIiIhGjxC4iEjFK7CIiEaPELiISMUrsIiIRo8QuIhIxSuwiIhGjxC4iEjFK7CIiEaPELiISMUrsIiIRo8QuIhIxSuwiIhGjxC4iEjFK7CIiEaPELiISMUrsIiIRo8QuIhIxSuwiIhGjxC4iEjFK7CIiEaPELiISMUrsIiIRo8QuIhIxSuwiIhGjxC4iEjFK7CIiEaPELiISMUrsIiIRU9DXAYh01bhx43jjjTf6OgyRnWLs2LGsW7euS9uYu+cmmq4EYeb5EIf0D2aGXi+yq8j0eg/LLd02mooREYkYJXYRkYhRYhcRiRgldhGRiMlpYjezu8xsg5m9kMt+RETkI7k+Yr8bOCmbirW1tTQ0NGR8vKGhodM6PZHr9kVEdpacnsfu7kvMbGw2de+9917i8TiVlZVMnDixzWM1NTVUVVURi8Uy1umJXLcvIvnn3//93/nxj3/c12HkRN78QKmxsRGAOXPmUF5eTnFxMRAcSVdVVdHc3Exzc3PaOj2R6/Zl1/THP/6RLVu28Pe//50RI0ZwxRVXZKx79913U1tbS2FhIRMmTGDmzJmJNt5++20aGxsZM2YMZ555Zpfbzsbzzz/Pvffeyy233NLt+h2NIdexzpkzh5deeolYLMZee+3FBRdc0Gk7a9eu5fnnn+9RLPksbxL7ggULACgoKGDChAl89rOfBaC+vp5YLJZIugCxWIz6+vpeSby5bl92PZs3b+acc86hvr6ewsJCRowYwamnnsrYsakfXl988UXuvvtuFi1aBMD06dM55ZRT2LhxI6+88gpf//rXAbj00ks56aSTaG5uzrptgNmzZ3P22WdnjPXWW29lyZIllJSUZDW2dPUzjWH79u1dirU7fW/ZsoXvfve7PPvsswAceeSRzJgxg+HDh3fY1rp16xgzZkybstmzZ3PFFVfw1ltvUVRU1K0Ys/Xqq6/yz//8z1xwwQV86Utfymqb6upqqqurs6qbN4n9uOOOA4LEPmPGjER5SUkJ8Xi8Td14PJ71C7EzuW5fdj1Dhw5lxYoVieQQj8cz/lL2z3/+M+Xl5Yn1srIyli5dyrBhw3jyySe56qqrGDhwIIMHD6awsJDi4uKs2wZ46aWXOkzsV199NcOHD2fhwoVZjS1d/UxjOP7447sUa3f6XrRoEQcffHBi/dBDD2XBggWcffbZvPvuuyxcuDDxy00zY/jw4RQVFTF16lTuvffeNu0fffTRHHTQQTlL6ieccAJz585NHLwWFRVRUVGR9fYVFRVt6s+aNStj3Z1xuqOFfx0qKiqioKCAysrKNkfKxcXFVFZWUlBQkLFOT+S6fdk1tSabxYsXc+yxxzJu3Li09QYPHsyOHTsS69u3b2fNmjVMnjyZlpYWPvGJT/Czn/2ME088kYEDB3apbQh+dp5rmcYAncdaV1fH9ddf36Zs7ty5VFVVZdX322+/3eYgrKSkhNdeew2AESNGcNZZZ3HmmWcm/p02bRpvvPEG8+bN480332xzBPzkk09ywgknZD3urqitrQWCA1eApqYm1q1bx7777puT/nJ6xG5m9wMVwHAzexNtnDLNAAAQDElEQVS43t3vTlf3ggsuoKSkJG1CnThxIuXl5dTX12es0xO5bl92TQ888AAPP/wwt956a8Y6Z555JnffHfyX2Lp1K6+88gpHHHEEANdeey0//OEPueaaa7jtttu63DawU66p09EYoONYy8rKmDlzJjfeeCPXXXcd1dXVrF27lq985StZ9b1p0yYGDRqUWC8sLGTr1q0dbnPuueeydu1ampqa2LZtW6J83rx5XH755Yn11atXc8899zBt2jRWrFjBf/7nf7Jy5UqWLVvGO++8wxFHHEFzczOPPfYYv/rVrwBoaWnhpptu4sADD2TDhg0sW7aMCy64gDvvvJM999yT3/72t3zhC1/gqaeeYtiwYcydO5eXX36ZgoICrrzyyqzGnI1cnxVzfrZ1R48e3eHjxcXFOU24uW5fdq4bbrgh7UfV66+/nhtuuKHT+pnqdcV5553HaaedxuTJk5k3b17aueWysjLuvvvuxH/8SZMmUVZWxmuvvcbChQt54oknePLJJ7nkkkuYNGkSRx55ZIdtr1mzht/85jeJ6YelS5fS2NiYmIo45phjOOWUU3o0rmzHkO1+mDx5Mk1NTfzLv/wLEyZM4Gtf+1rWfQ8ZMoT3338/sb5t2zb23HPPTrcbP348ixcvblP2zDPPJBL0xo0bmTFjBsuXL6e0tJSlS5cmyg844ACefPJJbrzxRgCuueaaRBvXXXcd+++/P2eeeSb33Xcf48eP55/+6Z/41a9+xdVXX83hhx8OBG8iZ511FieddBIlJSXccsstvZrYcfc+/wvCEMlOvr9eHnvsMT/qqKMS60cddZTfcsstWW37iU98wl977TW/5ZZbvKamJlE+f/58v/nmm7vc9qxZszrt89e//rVffPHFWcWXTf3WMXQl1mXLlvlXv/pVv+OOO7rU9+OPP+4XXXRRYv3iiy/2+++/P8uRfOTll1/2008/PbH+s5/9zI899lh/4okn/L777vN33nkn8di1117rf/jDH9zd/amnnvIzzjjD3d2bm5t92LBhvn37dnd3v/TSS33ZsmXu7n7AAQe06e+oo47yuro6d3e/8cYb/fbbb88YW6bXe1ieNqfqkgIivWzAgAGJkwHcnbfeeotDDjkECE6z86TpkTfeeINDDz0UCI62x44dy8c//nHKy8upqalJ1Nu+fTuf/OQn07Y9adKkjLEk99WR9vXax9lR/UxjyDbW5557jnnz5nHHHXdwyCGH8Itf/CLrWKdNm8bKlSsT6ytXruzWPPmCBQs4/vjjefTRRwEYNGgQM2bMYPr06Zx//vnU1dXR1NQEwPz58xN93HPPPVx44YU8+uijNDQ0sPfee1NUVERTUxPPPvssU6dOZfXq1Rx00EEAPPjgg2zZsoWmpiZKS0sBuP/++zn//PN5/PHHuxx3JnlzVoxIVJx88sm8/vrr/Nd//RdvvPEG3/72t5k+fToAZ599NnfddReTJ08GYK+99uKMM87g5z//OX/729+48847gWDe+vbbb+emm26iuLiYkpKSxBRK+7ZPPPHEjLF09uXpT3/6U37/+9/z1ltvMWvWLK6++mqGDBmSEmdH9TONId1+aB9rXV0djz76KN/5zneAIFE3NTXxxz/+kc985jNZxfqNb3yDG2+8EXfnmmuuaTMNlK3y8nIWLVrEgQceCATTRz/4wQ947LHHaGxsZMiQIRx66KFs27aNYcOGMXToUCD44ri+vp6DDz6Y3XffncrKSmbPns3q1avZf//9MTP22GMPhg4dyoMPPkhFRQUrV67ktNNOS/S97777Mn/+fI499tgux52JbrQh/Y5utJG93/3ud5x77rl9HcYuYcOGDZSUlFBUVMTNN9/MuHHjOOecc3rcbndutKHELv2OErvko0svvZQpU6aw++67U1tby3/8x3/0SrtK7LJLUGKXXYlujSciIkrsIiJRo8QuIhIxSuwiIhGjxC4iEjFK7CIiEaNfnkq/M3bs2J1yOVqRfNCdG5PoPHYRkX5I57HnmWxvbxUlGvOuQWPOD0rsfSAfXwi5pjHvGjTm/KDELiISMUrsIiIRkzdfnvZ1DCIi/U1eX91RRER6j6ZiREQiRoldRCRi+jSxm9nJZvaymb1qZr1zu5E+YmZ3mdkGM3shqWyYmT1hZq+Y2VwzG5r02B1m9pqZrTKzw5LKLwr3xytmduHOHkdXmNneZjbfzF4ysxoz+1pYHslxm1mRmT1tZs+F470+LB9nZsvC+B8ws4KwvNDMHgzH+1czG5PU1jfD8jVmlvmmpXnCzAaY2UozqwrXIz1mM1tnZs+Hz/UzYVn/eV27e5/8Ebyp/A0YCwwEVgEH9FU8vTCeTwOHAS8klf0I+Ea4/B/AD8PlU4DHwuVPAsvC5WHA34GhQEnrcl+PrYMx7wkcFi4PBl4BDojyuIGPhf/GgGXhOH4HfDYs/wVwebj8FeDn4fK5wIPh8kHAcwSX9BgX/j+wvh5bJ+P+N+C3QFW4HukxA2uBYe3K+s3rui+P2KcCr7n7G+6+A3gQqOzDeHrE3ZcAm9oVVwL3hMv38NH4KoHfhNs9DQw1s5HAScAT7r7Z3euBJ4CTcx17d7n7P9x9Vbi8FVgD7E2Ex+3uH4aLRQRJyoHjgD+E5fcAM8Pl5P0wGzg+XD6dIOE1u/s64DWC/w95ycz2BmYA/5tUfDwRHjNgpM5o9JvXdV8m9tHAW0nrb4dlUVLm7hsgSILAyLA809jbl9fST/aJmY0j+MSyDBgZ1XGHUxLPAf8A/kJwFFbv7i1hleTXcWJc7h4HNpvZHvSj8YZ+AlxD8CaGmQ0HNkV8zA7MNbPlZnZpWNZvXte6uuPOlenc0n59qUIzG0xwdHaVu29N87uEyIw7TGaTzWx34GGCqads9bvxmtmpwAZ3X2VmFckPZdtE70e1Uxzt7uvNrBR4wsxeIfV1nLev6748Yq8FxiSt7x2WRcmG8CMZZrYnUBeW1wL7JNVrHXu/2yfhl2azgXvdfU5YHPlxu/sWoBo4Eigxs9b/S8mxJ8ZrZjFgd3d/n8z7IR8dDZxuZmuBBwimVm4nmG6I6phx9/XhvxuBRwimjfrN67ovE/ty4ONmNtbMCoHPAVV9GE9vMNq+W1cBXwyXvwjMSSq/EMDMPkXwUX4DMBeYbmZDzWwYMD0sy2e/Ala7++1JZZEct5mNaD0Twsx2I4hzNbAA+GxY7SLajveicPmzwPyk8s+FZ5CUAx8Hnsn9CLrO3b/l7mPcfTzB/9H57v4FIjxmM/tY+CkUMysGTgRq6E+v6z7+5vlkgjMpXgOu7ctYemEs9wPvAI3Am8DFBN+KPxmO8QmgJKn+TwnODHgemJJU/sVwf7wKXNjX4+pkzEcDcYIzmp4DVobP6R5RHDcwKRzjKuAF4NtheTnwdBj774CBYXkR8PtwXMuAcUltfTPcD2uAE/t6bFmOfxofnRUT2TGHY2t9Tde05qb+9LrWJQVERCJGvzwVEYkYJXYRkYhRYhcRiRgldhGRiFFiFxGJGCV2EZGIUWKXFGYWDy/RWmNmc8Kfz3ennVFm9vteju318NojvdnmWDM7L2n9IjP7ryy3/b/wOjmZHr/ezK7ueZTZaT+WXmjvATPbt7fak51DiV3SaXD3Ke4+ieCKlVd2pxF3X+/u5/RuaBmvz9ET5cD5Xe3HzA4CBnhwtcKcSPrZfrbSjaWzPmIdPPwLgkvUSj+ixC6d+StJV6Qzs6+b2TPhDQVabzRxk5ldkVTnejO7Ojx6rAnLBpjZzRbcqGKVmV0Wlv/UzE4Llx82s/8Nly82s++liceS+vl82N5KM/uFmVlY/oGZ3Rj281R4ISfMbLwFN3943sy+Z2YfhE3dBHw6bOeqsGy0mf0pvEHCjzLsm8/z0c/KW28c82zY71+S6h1sZgvM7G9m9tWk+g+HVw+sSbqCYGv8t4RXkfyUmX0n3OcvmNl/J9Xb18z+Eva3wszGtx9LB/t9mpktMrM5wEvhz+gfteDGEi+YWevlAhYD/9SNNxjpS33981395d8f8EH4b4zg5+EnhuvTgV+Gywb8kY9uMFKdtP1LwF4EN1F5ISy7DPhWuFxIcK2gsQQ3Y/hRWP408FS4/CtgeprYXif4afcBBNfoiIXlPwO+EC63ADPC5R8l9ftH4Jxw+XJgS7ic+Kl8uH4Rwc/DBxP8RH4dMDpNLNXAweHyCIJLSYwJ10vCf68HlhBcSXU48G5SzK11BhH8dH1YUvxnJfWT/NP13wCnhsvLgNOT9umgNGPJtN+nAR8kxXtm63Mbrg9JWp4LTO7r16X+sv/Tu7Cks5uZrQTWA2UE1x2H4GJI08PHVgL7A/t5cLONUjPb08wOAd5393fatXkicGF4FPo0QXLej+CI8FgzO5DgglobLLhy3pHAU2lia50iOQGYAiwP2zyeYBoCoMndHw+XnyW4Yw9hm7PD5fs72Qfz3H2ruzeGcY1NU2cUsDFc/hSw0N3fBPDgxgqtHvPgBhPvARv46Dre/2pmqwgS9N4E+wOgGXgoafsTLLgN3QsEN/U4OLxI1V7uXhX21+Tu29PEmGm/AzzTGi/BG8v08NPXp939g6Q2NhK8UUs/oeuxSzofuvsUMxtEcLR2JcFFjgy4yd3vTLPN/xFczW9PgotCtWfAV939LykPmJUQ3G1mIUHiOYfgU0NDBzEacI+7fzvNY01Jy3E+ep0nz5t3ds3sxgxtJPuQ4Ci5szaT22oBCsxsGsGb0SfdvdHMFiS1td3dW29qUUTwaWSKu78TTn+11svmut9p93vYf2L/uvtrZjaF4E5JN5rZPHdvnQobBGzLoi/JEzpil3QMIDwCvAr4ejjHOhe4xIJLmWJme7XOXxNM2XwOOIsgybc3F7jCPrrp8X4WXPoWgiPWfwMWEUxbfJ3gSD5jbMA84Oyk+fNhZrZPuzrtLQPODpc/l1T+ATAkwzYdWUNw+dnWto8xs7Gt8XSy7VCCuxA1mtkBBEf8rZLjH0TwhvReeJR+NiRuRfiWmVWG/RWG+7P9WNLt94+1D8bMRgHb3P1+4P8Bk5MengC82Ml4JI/oiF3SSRzZenDnnOeB89z9vnDK5K/h95QfAF8ANrr7ajMbArzt4e3D2vlfgimRleGXnHV8dJ/MxQTz6WvN7E2Cyx0v6ig2d19jZtcR3N1mAMFR+pUEtyLLdEbLvwG/NbNvESS8zWH5C0BLOF3xa1LvXZupvccJpkbmu/u7ZvbPwMNJ4zspU/zAn4Evm9lLBJeB/Wu6/tx9s5ndSfC9xXraXsP8QuCXZvZdgvF/tv1Y3P12C07HTLffk00C/p+ZtYRtfQXAzMoIPsHVpdlG8pQu2yu7DDPbzd23hcvnAp9z9zN60N4gghtJHO0R/Y9kZv8KbHb3u/s6FsmejthlV3K4mbV+V7AJuKQnjbn79nDOezTBDYyjaBNwb18HIV2jI3YRkYjRl6ciIhGjxC4iEjFK7CIiEaPELiISMUrsIiIRo8QuIhIx/x9GS6mphyO5RgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fed88117be0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(lengths, ratings, color='grey')\n",
    "plt.plot(numpy.array(xplot), yplot, color = 'k', linestyle = '--',\\\n",
    "         label = r\"$3.983 + 1.193 \\times 10^{-4} \\mathit{length}$\")\n",
    "plt.xlim(0, 5500)\n",
    "plt.ylim(0.5, 5.5)\n",
    "plt.xlabel(\"Review length (characters)\")\n",
    "plt.ylabel(\"Rating\")\n",
    "plt.title(\"Rating vs. review length\")\n",
    "plt.legend(loc='lower right')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fit a slightly more complex model with two features (review length and number of comments)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def feature(d):\n",
    "    feat = [1] # Constant feature\n",
    "    feat.append(len(d['review_text'])) # Length of review\n",
    "    feat.append(d['n_comments']) # Number of comments\n",
    "    return feat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = numpy.matrix([feature(d) for d in data])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fit the model and extract the coefficients"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[3.95432090e+00, 7.24295694e-05, 1.08071028e-01]])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model = sklearn.linear_model.LinearRegression(fit_intercept=False)\n",
    "model.fit(X, y)\n",
    "theta = model.coef_\n",
    "theta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Mean Squared Error and R^2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Extract model predictions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_pred = model.predict(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sum of squared errors (SSE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "sse = sum([x**2 for x in (y - y_pred)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Mean squared error (MSE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[1.22988869]])"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mse = sse / len(y)\n",
    "mse"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### R^2 and fraction of variance unexplained (FVU)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[0.04057361]])"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fvu = mse / numpy.var(y)\n",
    "r2 = 1 - fvu\n",
    "r2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Can also get the R^2 coefficient directly from the library..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.040573608198668176"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.score(X, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simple feature transformations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Polynomial (quadratic) function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "def feature(d):\n",
    "    feat = [1]\n",
    "    feat.append(len(d['review_text']))\n",
    "    feat.append((len(d['review_text']))**2) # Quadratic term\n",
    "    return feat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = numpy.matrix([feature(d) for d in data])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 3.86513512e+00,  4.34379750e-04, -8.37700479e-08]])"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model = sklearn.linear_model.LinearRegression(fit_intercept=False)\n",
    "model.fit(X, y)\n",
    "theta = model.coef_\n",
    "theta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the R^2 coefficient (using the library)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.026559772903043233"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.score(X, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cubic function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "def feature(d):\n",
    "    feat = [1]\n",
    "    feat.append(len(d['review_text'])/1000)\n",
    "    feat.append((len(d['review_text'])/1000)**2) # Quadratic term\n",
    "    feat.append((len(d['review_text'])/1000)**3) # Cubic term\n",
    "    return feat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = numpy.matrix([feature(d) for d in data])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 3.92656525,  0.11531522,  0.14023331, -0.03494859]])"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model = sklearn.linear_model.LinearRegression(fit_intercept=False)\n",
    "model.fit(X, y)\n",
    "theta = model.coef_\n",
    "theta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.02926072167925875"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.score(X, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bar plot of daily ratings trends"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use a (slightly larger) dataset of fantasy reviews"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "path = dataDir + \"fantasy_10000.json\"\n",
    "f = open(path)\n",
    "\n",
    "data = []\n",
    "\n",
    "for l in f:\n",
    "    d = json.loads(l)\n",
    "    data.append(d)\n",
    "    \n",
    "f.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Extract averages for each day"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "weekAverages = defaultdict(list)\n",
    "\n",
    "for d in data:\n",
    "    timeString = d['date_added']\n",
    "    t = dateutil.parser.parse(timeString)\n",
    "    weekAverages[t.weekday()].append(d['rating'])\n",
    "\n",
    "for k in weekAverages:\n",
    "    weekAverages[k] = sum(weekAverages[k]) / len(weekAverages[k])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "defaultdict(list,\n",
       "            {0: 3.6558265582655824,\n",
       "             1: 3.7186505410566517,\n",
       "             2: 3.795968234575443,\n",
       "             3: 3.7673956262425445,\n",
       "             4: 3.7895545314900154,\n",
       "             5: 3.693290734824281,\n",
       "             6: 3.755786113328013})"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "weekAverages"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the averages for each day"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "xplot = list(weekAverages.keys())\n",
    "xplot.sort()\n",
    "yplot = [weekAverages[k] for k in xplot]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEZCAYAAACervI0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAGSZJREFUeJzt3Xm4JHV97/H3ZyAM66BIiOxuKEElgAYwEjyKIAgCalwwBBONmqhx1xv1mpm5MZHodY1GjRouCLjFiEpAIcJhEARRBh1AIyqLbLKDMkqA+d4/qg40zTkzPQN1lqn363n6meqqX3d9u+HUp36/WjpVhSSpf+bNdAGSpJlhAEhSTxkAktRTBoAk9ZQBIEk9ZQBIUk8ZAJoVkrwkyTdmuo41lWRhks/OwHrfneT6JFdP97rb9R+V5P/MxLr1wBkAWiNJLkuyPMltSa5uNwQbjvja7ZOsSHLP/39VdXxV7d9dxdNiWi+qSbIt8CZgx6raajrXrbWDAaA1VcCBVbUA2AXYFXj7iK9N+/p0VFtfbA/cUFU3znQhmpsMAD0QAaiq64Bv0gRBsyB5dpLzk9ya5PIkCwded0b77y1tD2KPJC9NcubA61ckeVWSnyS5KclHB5bNS/L+dujjZ0leM9yjGGj7tiRfGpr34SQfaqf/vH2P29p/DxvpgyePSDLefr5vApsPLf9ikmuS3Ny226md/+Qk1ybJQNvnJblgivUsSHJMkuuSXJrkne38fYBTgK3a2v9tkteOJ3luO/3U9js6oH3+jCRLB9q+LMnFSW5McnKS7QaW7ZjklHbZj5K8YIpaN0ly2sR3q9nPANADlmQb4ADgkoHZvwb+rKo2BQ4E/irJwe2yvdt/F1TVgqo6t30+PIRyIPAk4A+AFybZr53/SuBZwM7AbsChk7x2wueBA5Js1NY6D3gBcFw7ZPVh4FltT+aPgEk3xJM4HjiPZsP/buClQ8tPAh4NbAGcDxwHUFXfA24A9htoezjw/6ZYz0eBTYBHAGPAEUn+oqq+RfOdX91+hy+b5LVntK+B5jv/Gfd+908DxgGSHAL8Lc33+LvAmcDn2mUb0gTNse1nfTHwL0l2HFxRks2A/wLOrKo3TPFZNNtUlQ8fq/0ALgVuax8rgFNpNuhTtf8g8P52envgbmDewPKXAksGnq8AnjLw/AvA29rpbwGvGFi2z/D7Da17CXB4O70vcEk7vSFwE/BcYP3V+OzbAv8DbDAw7zjgmCnaP6T9PJu0z98GHNtObwbcDvzeJK+bB9wBPG5g3iuB09rppwFXrKTOZwAXtNMnAy8Dzm6fjwOHttMnAX8xtN7b28/5QuCMoff9BPCudvoo4DPAMuBNM/3/pY/Ve9gD0ANxSDV7zk8DdmRgGCTJ7u1wwHVJbgFexdAwyQh+OTC9HNi4nd4K+MXAssHpyXwOmBjaOYxm752qWg68CPhr4JokX0/yuBHq2gq4uap+MzDv8omJdojqyCQ/bT/7pTQ9lInPfyxwUJINaDawS6pq8LNO2BxYF7hiaD1bj1AjwHeAxybZgqYXdQywbZKHAbtz71Dc9sCH26G2m4Ab23q3bpftObEsyc3AS4DfG1jPgcD6wCdHrEuzhAGgB2LiGMCZwNHA+weWHQ+cAGxdVQ+h2ThMjHs/0LNlrgG2GXi+3VQNW18CxpJsTbO3f/zEgqo6tar2Ax4O/DfwqRHX/9B2Az5ZDX8KPAd4RvvZH0Hz2Se+r6tpNs7Ppxn+mer00RuAO2k2whO2B64aoUbagPo+8Hrgwqq6q13vm4CfVtXNbdMrgFdV1Wbt46FVtXFVnUMTruNDyxZU1WsHVvWvwDeAk4e+E81yBoAeLB8C9k3yxPb5xjR7yXcm2Z1mr3HC9TRDIo9ew3V9EXh9kq2SPIRmSGVKVXUDzd7uUcDPq+q/AZJskeTgdpz7TprjFnevauVVdQXwPWBxkt9JshfNBn/CxjRDNze3xx7ew/1D77Nt3U8A/mOK9axoP+s/JNk4yfbAG5k6MCazBHgt9+7tjw89hyac3zFwoHrTJH/SLjuRphdxeJJ128/75OGeUlX9DU2Anphk/dWoTzPIANCaus8Grd3IHg38XTvrNcDfJ7kV+N80Y/gTbX8D/ANwVjussPuq3n/o+adoDkz+kGYP9z+Bu9oN5lSOpzlWcNzAvHk0e8NX0ext700zHESSvZLctpL3ewmwJ81wybtoPvuEY2j2qq8CLgTOnuT1X6HZm/+PqvrtStbzOprhr5/TbMyPraqjVtJ+2Bk0gbRk6Pk9AVBVJwBHAp9vh6x+COzfLvs1zQHrFwNXt48jgfmTrOuVND2GE5Kstxo1aoakqvtrV9ozL74HXFlVBw8tW4/mD+ZJNH+EL2r3sKSRJNkf+HhVPXKma1kdSX4KvLKqTpvpWtRP09UDeD1w8RTLXg7cVFU70AwjvHeaatIclWT9JAckWacd11/IFMMos1WS5wMr3PhrJnUeAO054s8GPj1Fk0O4t/v87zTddGllAiymOYXz+8BFNCEwJyQ5HfgY8OqZrkX9tu40rOODwFuBTadYvjXtaXxVdXeSW5JsVlU3TUNtmoPaYwiTHTeYE6rq6TNdgwQd9wCSHAj8sqouYOA0uFW9rMuaJEmNrnsATwUOTvJsYANgkyTHVNURA22upLni8Ook69BcTXq/vf8k03qnRUlaW1TVpDvWnfYAquodVbVdVT2K5jSy04Y2/gBf5977qLwAmPKg2HReIr1w4cIZv0y7T3XP5drnat1zufa5WvdM1L4yM3IdQJLFSQ5qn34G2DzJJcAbaG5KJUnq2HQcBAagqs6gvfikqhYOzL+D5n4okqRp5JXAUxgbG5vpEtbIXK0b5m7tc7VumLu1z9W6YXbVPi1XAj8YktRcqVWSZosk1EwcBJYkzV4GgCT1lAEgST1lAEhSTxkAktRTBoAk9ZQBIEk9ZQBIUk8ZAJLUUwaAJPWUASBJPWUASFJPGQCS1FMGgCT1lAEgST1lAEhSTxkAktRTBoAk9ZQBIEk9ZQBIUk8ZAJLUUwaAJPWUASBJPWUASFJPGQCS1FMGgCT1lAEgST1lAEhSTxkAktRTBoAk9ZQBIEk9ZQBIUk8ZAJLUUwaAJPVUpwGQZH6Sc5MsTbIsycJJ2myb5LQk5ye5IMkBXdYkSWqkqrpdQbJhVS1Psg5wFvC6qvruwPJPAudX1SeT/D5wUlU9cpL3qa5rlaS1TRKqKpMt63wIqKqWt5PzgXWB4a34CmBBO/0Q4Kqua5IkTU8PYB7wfeDRwMeq6u1Dyx8OnAI8FNgQeGZVLZ3kfewBSNJqWlkPYN2uV15VK4BdkywATkiyU1VdPNDkMOCoqvpgkj2BY4HHT/ZeixYtumd6bGyMsbGxzuqWpLlofHyc8fHxkdp23gO4z8qSdwG3V9UHBuZdCDyrqq5qn/8M2KOqbhh6rT2AWWjx4sUzXcKkFi683/kGUi/N2DGAJJsn2bSd3gDYF/jxULPLgWe2bX4fmD+88ZckPfi6HgLaEji6PQ4wD/hCVZ2UZDFwXlWdCLwF+FSSN9IcEH5pxzVJc5q9Lj1YOg2AqloG7DbJ/IUD0z8C9uqyDknS/XklsCT1lAEgST1lAEhST3V+HYA0W83Gg6keSNV0sgcgST1lD0CSVmE29hbhgfcY7QFIUk8ZAJLUUwaAJPWUASBJPWUASFJPGQCS1FMGgCT1lNcBzBKz8Txjr0qV1m72ACSppwwASeopA0CSesoAkKSeMgAkqac8C0jStPFst9nFHoAk9ZQBIEk9ZQBIUk8ZAJLUUwaAJPWUASBJPWUASFJPrVXXAczGc4yh3+cZS5q97AFIUk8ZAJLUUwaAJPWUASBJPWUASFJPGQCS1FMGgCT1VKcBkGR+knOTLE2yLMmkJ8QneWGSi9o2x3ZZkySp0emFYFV1R5KnV9XyJOsAZyU5uaq+O9EmyWOA/wU8papuS7J5lzVJkhqdDwFV1fJ2cj5N4NRQk1cAH6uq29r2N3RdkyRpGgIgybwkS4FrgVOr6ryhJo8FHpfk20nOTvKsrmuSJE1PD2BFVe0KbAPskWSnoSbrAo8B9gZeAnwqyYKu65Kkvpu2m8G14/unA/sDFw8suhI4p6pWAJcl+QmwA/D94fdYtGjRPdNjY2OMjY11WbIkzTnj4+OMj4+P1LbTAGgP6N5ZVbcm2QDYFzhyqNkJwGHA0W37HYCfT/Z+gwEgSbq/4Z3jld0luesewJY0G/Z5NMNNX6iqk5IsBs6rqhOr6ptJ9ktyEXAX8JaqurnjuiSp97o+DXQZsNsk8xcOPX8z8OYua5Ek3ZdXAktSTxkAktRTBoAk9ZQBIEk9ZQBIUk+NdBZQko9MMvtW4HtV9dUHtyRJ0nQYtQewPrALcEn72Jnm1g4vT/KhjmqTJHVo1OsAdgaeWlV3AyT5OHAmsBewrKPaJEkdGrUH8FBg44HnGwGbtYFwx4NelSSpc6P2AN4LXJBkHAjNnTv/MclGwH91VJskqUMjBUBVfSbJScDu7ax3VNXV7fRbO6lMktSp1TkNdB5wPXAz8Jgke3dTkiRpOox6Gug/AS8CLgJWtLMLWNJRXZKkjo16DOBQ4HFV5QFfSVpLjDoE9HPgd7osRJI0vUbtASynOQvoWwyc9llVr+ukKklS50YNgK+1D0nSWmLU00CP7roQSdL0WmkAJPliVb0wyTKas37uo6p27qwySVKnVtUDeH3770FdFyJJml4rPQuoqq5pJ19dVZcPPoBXd1+eJKkro54Guu8k8w54MAuRJE2vVR0D+GuaPf1HJfnhwKJNgLO6LEyS1K1VHQM4HjgZeA/wtwPzf1VVN3VWlSSpcysNgKq6leanHw8DSLIFza+DbZxk46q6ovsSJUldGOkYQJLnJLkEuBQ4A7iMpmcgSZqjRj0I/G5gT+AnVfVIYB/gnM6qkiR1btQAuLOqbgTmJZlXVacDT+6wLklSx0a9F9AtSTamuf//cUmuA27vrixJUtdG7QEcQnNH0DcC3wB+Bjynq6IkSd0b9WZwE3v7K4Cjk8yjOTPouK4KkyR1a6U9gCQLkrw9yUeT7JfGa2l+IOaF01OiJKkLq+oBfJbmR+C/A/wl8A4gwKFVdUHHtUmSOrSqAHhUVT0RIMmngWuA7arqt51XJknq1KoOAt85MVFVdwNXuvGXpLXDqnoAf5DktnY6wAbt8wBVVQs6rU6S1JlV/R7AOlW1oH1sUlXrDkyvcuOfZH6Sc5MsTbIsycKVtH1+khVJdluTDyJJWj2jXgi2RqrqjiRPr6rlSdYBzkpyclV9d7Bde5HZ6/D2EpI0bUa9EGyNVdXydnI+TeDc77eFgb8HjgTu6LoeSVKj8wBIMi/JUuBa4NSqOm9o+a7ANlXl3UUlaRpNRw9gRVXtCmwD7JFkp4llSQJ8AHjzwEvSdU2SpI6PAQyqqtuSnA7sD1zczt4EeDww3obBw4GvJjm4qs4ffo9FixbdMz02NsbY2FjXZUvSnDI+Ps74+PhIbTsNgCSb09xK+tYkG9D8uPyRE8ur6jZgi4H2pwNvqqqlk73fYABIku5veOd48eLFU7bteghoS+D0JBcA5wLfrKqTkixOctAk7QuHgCRpWnR9Gugy4H7n9VfVpNcDVNUzuqxHknSvzg8CS5JmJwNAknrKAJCknjIAJKmnDABJ6ikDQJJ6ygCQpJ4yACSppwwASeopA0CSesoAkKSeMgAkqacMAEnqKQNAknrKAJCknjIAJKmnDABJ6ikDQJJ6ygCQpJ4yACSppwwASeopA0CSesoAkKSeMgAkqacMAEnqKQNAknrKAJCknjIAJKmnDABJ6ikDQJJ6ygCQpJ4yACSppwwASeopA0CSesoAkKSeMgAkqac6DYAk85Ocm2RpkmVJFk7S5o1JLkpyQZJTk2zbZU2SpEanAVBVdwBPr6pdgV2AA5LsPtTsfOBJVbUL8GXgfV3WJElqdD4EVFXL28n5wLpADS0/o6p+2z49B9i665okSdMQAEnmJVkKXAucWlXnraT5y4GTu65JktTskXeqqlYAuyZZAJyQZKequni4XZLDgScBT5vqvRYtWnTP9NjYGGNjYw96vZI0l42PjzM+Pj5S284DYEJV3ZbkdGB/4D4BkOSZwNuBvavqzqneYzAAJEn3N7xzvHjx4inbdn0W0OZJNm2nNwD2BX481GZX4BPAwVV1Y5f1SJLu1XUPYEvg6CTzaMLmC1V1UpLFwHlVdSLwXmAj4EtJAlxeVYd2XJck9V6nAVBVy4DdJpm/cGB63y5rkCRNziuBJamnDABJ6ikDQJJ6ygCQpJ4yACSppwwASeopA0CSesoAkKSeMgAkqacMAEnqKQNAknrKAJCknjIAJKmnDABJ6ikDQJJ6ygCQpJ4yACSppwwASeopA0CSesoAkKSeMgAkqacMAEnqKQNAknrKAJCknjIAJKmnDABJ6ikDQJJ6ygCQpJ4yACSppwwASeopA0CSesoAkKSeMgAkqacMAEnqKQNAknqq0wBIMj/JuUmWJlmWZOEkbdZL8vkklyT5TpLtuqxJktToNACq6g7g6VW1K7ALcECS3YeavRy4qap2AD4EvLfLmkZ16aWXznQJa2Su1g1zt/a5WjfM3drnat0wu2rvfAioqpa3k/OBdYEaanIIcHQ7/e/APl3XNIrLLrtspktYI3O1bpi7tc/VumHu1j5X64bZVXvnAZBkXpKlwLXAqVV13lCTrYFfAFTV3cAtSTbrui5J6rvp6AGsaIeAtgH2SLLTKl6SrmuSJEGqhkdkOlxZ8i7g9qr6wMC8k4FFVXVuknWAa6pqi0leO32FStJapKom3bFet8uVJtkcuLOqbk2yAbAvcORQs68DLwXOBV4AnDbZe031ASRJa6bTAAC2BI5OMo9muOkLVXVSksXAeVV1IvAZ4LNJLgFuBF7ccU2SJKZ5CEiSNHt4JfCQJO9McmGSHyQ5P8kfznRNq5JkRZJjBp6vk+T6JF+bybpGkWSz9kLB85Nck+TKgedd91DXSJIPJHndwPNvJPnXgef/N8kbZqa6VZuL3/mgJHe3tU7UPCcuHp2N25ZZ/x97OiXZE3g2sEtV3dWejrreDJc1ituBJySZ3158ty/tqbWzXVXdBOwKkOTvgF8PniQwS51Fc7zqI0kCbA5sMrD8j4BZGwBz9DsfdHtV7TbTRayO2bptsQdwX1sCN1TVXdD8oVTVtTNc06hOAg5spw8DPjeDtaypuXKg/2yajTzA44ELgV8l2TTJesCOwPkzVdxqmivf+aC5WPOs3LYYAPd1CrBdkh8n+ViSvWe6oBEV8HngsCTzgZ1pzqpSB6rqGuDOJNvQBMHZNN/3U4AnA8sm/tDViQ0GhoC+PNPFjGhWblsMgAFVdTuwG/BK4Hrg80mOmNmqRlNVFwKPoNn7/0/m5l7SXHI28FSaAPgOcM7A87NmsK4+WF5Vu1XVrlX1/JkuZhSzddviMYAh1ZwWtQRYkmQZcARwzMpfNWt8DXgfMEYzLq3uTAwDPYFmCOhK4M3ArcBRM1iXZqnZuG2xBzAgyWOTPGZg1i7A5TNVz2qY2Nv/N2BxVV00k8X0xNnAQTR3sq2quhl4CM0w0NkzWtnab871bmfrtsUewH1tDPxzkk2Bu4Cf0nTZZrsCqKqrgI/OcC19sQx4GHDs0LwN27Ns1J25ePHSrNy2eCGYJPWUQ0CS1FMGgCT1lAEgST1lAEhSTxkAktRTBoAk9ZQBoLXKwK2CL2zvFfOm9o6dXa7zfUmWJfmnjtezfXsFqfSg8EIwrW3uuVVw+5OknwMWAIs6XOcrgIfW9FxU44U7etDYA9Baq6puoLna8rVwzx70kiTfax97tvOPTnLwxOuSHJvkOcPvN7Cn/4MkL2jnfZXmKs/vT8wbaP/DJAva6RuSHD6wvn2SzEvy3iTnJrkgySsGXvuWJN9t5y+cpJZHtT2dJz3wb0p9ZQ9Aa7WqurTd0P4u8EvgmVX1P+19WT4H/CHN71K/Efhau8F+Cs2Nuu6R5HnAzlX1xCRbAOclWVJVhyS5bYofKPk28NQkVwA/A/6Y5tYRTwH+Cng5cEtV7dH+jsBZSU4BHgvsUFW7t8NXX0uyF+2P/CR5LM3tv49o7wIrrREDQH0wcQxgPeCjSXYB7gZ2AKiqJe092h8G/Anw5apaMfQee9H+yE5VXZdknCY8TmTqm5N9G3gazU2/PgG8IslWNDeQ+02S/YAnDvQcFrQ17Qfsm+T89r03auf/AtgCOAF4XlX9eE2/EAkMAK3lkjwKuKuqrm+HUq6tqp2TrAP8ZqDpMcCfAS8G/nyUtx6YnmpcfgnwGmBb4J3Ac2kC5syB9/ibqjp1qOb9gfdU1aeG5m9Pc7vpK2h6EwaAHhCPAWhtc8+GuR32+Tjwz+2sTYFr2ukjgHUGXnc0ze/41hR71mcCLxoYTvpj7v3VtUl7AFV1Jc3vMuxQVZfR9AjeQhMMAN8EXj3xQ+xJdkiyYTv/ZUk2audv1a4T4A6aIDkiyWGr+C6klbIHoLXN+u3QyXrAncAxVfXBdtm/AF9uf4npG8DtEy9qh3V+BHxlsjetqq+0B41/AKwA3lpV108sXkk953DvjtaZwD/SBAHAp2l+xe38dqz/OuDQqjo1yY7Ad9ozWH8FHN6ul3b46CDglCS/qqoTR/hepPvxdtAS0O55/wDYrap+NdP1SNPBISD1XpJ9gIuBj7jxV5/YA5CknrIHIEk9ZQBIUk8ZAJLUUwaAJPWUASBJPWUASFJP/X+mO8ulqHsw8QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fed851d8748>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.bar(xplot,yplot,color='grey',lw=0)\n",
    "plt.xticks(numpy.arange(0.4,7.4,1), [\"S\", \"M\", \"T\", \"W\", \"T\", \"F\", \"S\"])\n",
    "plt.xlim(-0.2,7.0)\n",
    "plt.ylim(3,4)\n",
    "plt.ylabel(\"Rating\")\n",
    "plt.xlabel(\"Day of week\")\n",
    "plt.title(r\"Rating vs. day of week\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Binary (one-hot) features"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read a small dataset of beer reviews with gender attributes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "path = dataDir + \"beer_500.json\"\n",
    "f = open(path)\n",
    "\n",
    "data = []\n",
    "\n",
    "for l in f:\n",
    "    d = demjson.decode(l) # Handles the looser json format in this file (fairly slow)\n",
    "    data.append(d)\n",
    "    \n",
    "f.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What does the data look like?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'beer/ABV': 5.0,\n",
       " 'beer/beerId': '47986',\n",
       " 'beer/brewerId': '10325',\n",
       " 'beer/name': 'Sausa Weizen',\n",
       " 'beer/style': 'Hefeweizen',\n",
       " 'review/appearance': 2.5,\n",
       " 'review/aroma': 2.0,\n",
       " 'review/overall': 1.5,\n",
       " 'review/palate': 1.5,\n",
       " 'review/taste': 1.5,\n",
       " 'review/text': 'A lot of foam. But a lot.\\tIn the smell some banana, and then lactic and tart. Not a good start.\\tQuite dark orange in color, with a lively carbonation (now visible, under the foam).\\tAgain tending to lactic sourness.\\tSame for the taste. With some yeast and banana.',\n",
       " 'review/timeStruct': {'hour': 20,\n",
       "  'isdst': 0,\n",
       "  'mday': 16,\n",
       "  'min': 57,\n",
       "  'mon': 2,\n",
       "  'sec': 3,\n",
       "  'wday': 0,\n",
       "  'yday': 47,\n",
       "  'year': 2009},\n",
       " 'review/timeUnix': 1234817823,\n",
       " 'user/profileName': 'stcules'}"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Filter the dataset to include only those users who specified gender"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = [d for d in data if 'user/gender' in d]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How many users have specified gender?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "215"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Binary representation of gender attribute"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "def feat(d):\n",
    "    return [1, 1.0 * (d['user/gender'] == 'Female')]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = [feat(d) for d in data]\n",
    "y = [len(d['review/text'].split()) for d in data]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fit a model to the binary feature"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([127.07177033,   8.761563  ])"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model = sklearn.linear_model.LinearRegression(fit_intercept=False)\n",
    "model.fit(X, y)\n",
    "theta = model.coef_\n",
    "theta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The model indicates that females tend to leave slightly longer reviews than males. Plot data from the two groups."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAEZCAYAAAB1mUk3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAHhFJREFUeJzt3XuclmW97/HPlwFEPFCAYimSElqyKU8ppNmUFaPWci0qyrTUXGvX3mWtWrUC1oGZXOnWti9d7XYHC1tQS00LU3sVsDyMW3OAPKGJKZ2UDFFAUjFlgN/+475GnhnncD/znOaZ5/t+vZ7Xfb7v3zMw92+u676v61JEYGZmjW1ErQMwM7PaczIwMzMnAzMzczIwMzOcDMzMDCcDMzPDycCGOUm7JB1a6ziKJWlKir3qv6OSzpZ0R7Wva7XlZGDDXT03pKl47P0knXr+udkgOBnYcKdBHSQ1lTuQIUpkN/5B/Zxs+HAysKqTdLSkeyX9WdK1kq6R9OWC7e+VdJ+kZyTdKWlGwbbfS/oHSWvS9qsljS7Y/kVJf5L0R0nnUvAXrqTRkv63pMckbZD0DUl7pG1vl7Re0j9K2gBc2SPm0el6RxSsmyjphTSdIOmmtM9mSbeX+We2r6Tvpu+2XtIFkpS2nS3pDklflbRF0m8ltRQc+zpJt6ef9wpJX5e0JG3uinOrpGclHb/7sN7PZ8OTk4FVlaRRwFKym+144Grgbwq2HwUsAv4ubf82cGM6rssHgfcAhwBvBs5Jx7YAnwdOBqYB7+px+YuB1wNvStMDgX8t2H4A8CrgYOC/Fx4YEduBHwNnFKyeC7RHxCbgH4D1wARgf2BBrh9IfouB7cChwFHAu4G/Ldh+HPBwuv5XyX6GXa4CVqZtbcBHC7adlKb7RsS+EbEqLR/fz/lsOIoIf/yp2gd4G7C+x7o7gC+n+W8AbT22/xp4W5r/PXBGwbaLgW+k+UXAhQXbpgG7gEPT8vPAIQXbZwG/S/NvB14ERvUT+8nAbwqW7wTOTPNtwPXA1Bw/g3OB+4ENZDf5k8iS0EeBz6R9pgA7yf5gm5Ri26PgHB8Gbk3zZwOPFmzbM33v/YHJZElkTMH27wNLel6nYHtv59sJ7F/r/z/+VO4zMn/aMCuL1wJP9Fi3vmB+CvAxSeenZQGj0nFdNhbMvwC8puDcdxdse6xrRtJ+wFjgnlS7AtmNtrCu/OmI6Own9tuAPSW9BXiKrFTyk7TtErKEsEJSAN+JiIv7OM8pZCWbF4CPAJeRlUpWAF/oZf+DyX4GG7pqhtLn8YJ9nuyaiYi/pP32BvYDtkTEiwX7rgcO6ud79nY+pfM9NcBxVqecDKzaNpBVzxSaDPwmza8HvhIRFw3y3JMLlqew+5nBJrKb7/SI2NDH8f2+QRMRuyRdS3YD3wj8NCK2pW3byG7kX0jPFW6TtDoibuvlPHMLFq9In/6sJysZTIiIYt/y2QCMlzSmICFMZvd39VtDBviZgVVfB7BT0qckNUk6nay+u8t3gE9KOg5A0l6STpW0V45zXwucI+mNksZS8Dwg3US/A1yeSglIOlDSe4qM/2rgQ2QJ4aqulZJOkzQ1LT4H7CCrqimFUuxPkpUaLpO0jzKHSjqp/8MhIh4nKy21SholaRbwvoJdnk5xTu3teGscTgZWVakaZg7Zw89nyG6qNwEvpe33kD08/rqkLcCjZHXYL5+in3MvAy4Hbk3H3dJjly+RlUBWStpKdoM9rMj4VwPbyKqmfl6waRpws6TngF8A/zcibgeQ9DNJ84q5TtflCuY/BowG1gJbgOvIqpbyHHsm8Fay0tGXgWvY/fP+C/AV4BfpzaHjep6ol/PZMKTiS51m5SVpJfDNiFhc61gagaRrgIcjoq3WsdjQ4ZKBVZ2kkyRNStVEZwMzgGW1jmu4knRsqlZSev32r9j94NsM8ANkq43Dyer3xwK/A94fERv7P8RKcABZ247xwB+BT0bEmtqGZEONq4nMzMzVRGZmVqfVRKlRj5mZFSkieu2UsC6TAWTdaDSi1tZWWltbax2GWd1q5N+hgtb3r+BqIjMzczIwMzMng7rT3Nxc6xDM6pp/h3pXl6+WShpEf11mZo1NUp8PkF0yMDMzJwMzM3MyMDMznAzMzAwnAzMzw8nAzMxwMjAzM5wMzMwMJwMzM6MKyUDSHyStkXSfpNVp3aslrZD0iKTlksYV7P81Sesk3S/pyErHZ2Zm1SkZ7AKaI+KoiDgurZsH3BwRhwO3AvMBJJ0CTI2IacAngG9VIT4zs4ZXjWSgXq5zOrA4zS9Oy13rlwBExCpgnKRJVYjRzKyhVSMZBLBc0i8l/W1aN6lrAPSIeBLouuEfCKwvOPaJtM7MzCqoGiOdnRARGyTtB6yQ9AhZgijkLkjNzGqo4skgIjak6dOSfgIcB2yUNCkiNko6AHgq7f4EMLng8IPSulcoHLauubnZfZSbmfXQ3t5Oe3t7rn0rOp6BpLHAiIh4XtJewAqgDTgZ2BIRF0uaB7wqIuZJOhX4VEScJmkmcHlEzOzlvB7PwMysSP2NZ1DpksEk4HpJka71nxGxQtLdwLWSPg48BswFiIifSTpV0m+AbcC5FY7PzMzwSGdmZg3DI52ZmVm/nAzMzMzJwMzMnAzMzAwnAzMzw8nAzMxwMjAzM5wMzMwMJwMzM8PJwMzMcDIwMzOcDMzMDCcDMzPDycDMzHAyMDMznAzMzAwnAzMzw8nAzMxwMjAzM5wMzMwMJwMzM8PJwMzMcDIwMzOcDMzMDCcDMzPDycDMzHAyMDMznAzMzAwnAzMzw8nAzMxwMjAzM5wMzMyMKiUDSSMk3SvpxrT8OkkrJT0q6WpJI9P60ZKukbROUoekg6sRn5lZo6tWyeCzwNqC5YuBSyPiMGArcF5afx6wJSKmAZcDl1QpvrrQ0dFBS0sLHR0dtQ7FzIaZAZNB+qv+KEmnSXqnpP2LuYCkg4BTge8WrH4n8OM0vxj46zR/eloG+BFwcjHXGu7a2tpYvnw5bW1ttQ7FzIaZkX1tkDQV+BLwLmAd8DQwBjhM0gvAt4HFEbFrgGtcBnwRGJfOOwF4puC4PwIHpvkDgfUAEbFT0lZJ4yNiy2C+3HCzcOHCblMzs3LpMxkA/wZ8E/hEREThhlQ6+AjwUXb/Jf8Kkk4DNkbE/ZKaCzfljK/P/VpbW1+eb25uprm5ua9dh41Zs2axbNmyWodhZnWivb2d9vb2XPuqx32+rCRdCJwF7AD2BPYBfgK8BzggInZJmgksjIhTJC1L86skNQEbIuIV1VKSeuYnMzMbgCQiotc/svM8M/igpH3S/L9IWirp6DwXjogFEXFwRBwKfBi4NSLOAm4DPph2Oxu4Ic3fmJZJ22/Ncx0zMytNnreJ/iUinpN0ItkD3UVk1UelmAd8XtKjwPh0TtJ0oqR1wN+n/czMrMIGrCaSdF9EHCXpIuDBiLiqa111Quw1JlcTmZkVqaRqIuAJSd8GPgT8TNIeOY8zM7M6kadkMBZoISsVrJP0GmBGRKyoRoB9xOSSgZlZkforGfSZDCSN7++ktXz338nAzKx4/SWD/toZ3AME2bv+BwPPpPlXAY8Dh5Q5TjMzq5E+6/4j4pD0SujNwPsiYmJETADeC9SsisjMzMovzzODByNixkDrqsnVRGZmxRtsNVGXP0n6Z+AHaflM4E/lCs7MzGovzyuiZwD7AdcDS9P8GZUMyszMqqvfkkHqH2hBRHy2SvGYmVkN9FsyiIidwIlVisXMzGokzzOD+9JwldcB27pWRsTSikVlZmZVlScZjAE2k41O1iXInh+YmdkwUNHxDCrFr5aamRWv1PEMDpJ0vaSn0ufHaVxjMzMbJvK8Wvo9skFnXps+N6V1ZmY2TORpgXx/RBw50LpqcjWRmVnxSh3PYLOksyQ1pc9ZZA+UzcxsmMiTDD4OzAWeBDYAHwDOrWRQZmZWXXmqicZExItViicXVxOZmRWv1I7qfiVpI3BH+twZEX8uZ4BmZlZbudoZSDoYeBtwAnAqsNUPkM3M6ktJJYPUpuAEsmTwZuAh4M6yRmhmZjWV55nBLuCXwIURcUNVohqASwZmZsXrr2SQJxm8mazn0pPIxkJeB9weEYvKHWheTgZmZsUrKRmkE+xNlhDeBpwFEBFTyhlkMZwMzMyKV+ozg7uBPYC7yN4mOikiHitviGZmVkt5qon2i4inqxRPLi4ZmJkVb1DdUaQuKEb0lQgkTZXkUdDMzIaB/qqJJpCNcnYPcA/wNNlAN68H3g5sAuZVPEIzM6u4fquJJDWRjXB2AvAa4C/Aw8DPI+LxqkTYe1yuJjIzK1LJbxMNNU4GZmbFK7ULazMzG+Yqmgwk7SFplaT7JD0oaWFa/zpJKyU9KulqSSPT+tGSrpG0TlJH6hPJzMwqrKLJICJeAt4REUcBRwKnSDoeuBi4NCIOA7YC56VDzgO2RMQ04HLgkkrGZ2ZmmTztDPYA3g+8joK3jyLiy0VdSBoL/D/gfwI/BQ6IiF2SZgILI+IUScvS/Kr08PrJiNivl3P5mYGZWZFKfWZwA3A6sAPYVvDJe/ERku4jGyntv4DfknWBvSvt8kfgwDR/ILAeICJ2Alsljc97LTMzG5w8g9scFBEtg71AuukfJWlf4HrgDUUc3msGA2htbX15vrm5mebm5kFGaGY2PLW3t9Pe3p5r3zzVRFcA/yciHiw1MEn/QtZW4R/JV020ISL27+U8riYyMyvSYLujeFDSA2S9ld4r6RFJDxSsz3PhiZLGpfk9gXcDa4HbgA+m3c4mq4oCuDEtk7bfmuc6ZmZWmj5LBpL67aI6T8+lkmYAi8mSzgjghxHxFUmHANcArwbuA86KiM70sPr7wFHAZuDDEfGHXs7rkoGZWZFKHdzm+xHx0YHWVZOTgZlZ8Up9m2h6j5M1AceUIzAzMxsa+ntmMF/Sc8CbJD2bPs8BT7G7jt+qqKOjg5aWFjo6OmodipkNM3mqiS6KiPlViieXRq0mamlpYfny5cyePZtly5bVOhwzqzMlDXsJXCfp6B7r/gw8FhE7So7Oclu4cGG3qZlZueQpGawEjgYeIGsENgP4FTAO+B8RsaLSQfYSU0OWDMzMSlHqA+Q/AUdFxLERcQxZh3O/I2sz4I7kqsjPDMysUvJUEx0WEQ91LUTEWklviIjfSX32FmEV8LnPfY5Vq1axdetWVq5cWetwzGwYyZMMHpL0TbJGYgAfAtamBmKdFYvMzMyqJk810TnAb4C/T5/fpXWdwDsqFZi90mWXXcbs2bO57LLLah2KmQ0zA5YMIuIvwKXp09PzZY/IzMyqLs/bRCcArcAUug9uc2hFI+s/poZ8m2j69OmsXbuWI444goceemjgA8zMCpTaN9Gvgc8B9wA7u9ZHxOZyBlmMRk0Ge++9N9u2bWOvvfbi+eddKDOz4pT6aumfI+LnEfFURGzu+pQ5Rsth//337zY1MyuXPG8T3Sbpq8BS4KWulRFxb8Wisl51lYYasVRkZpWVJxkcn6bHFqwL4J3lD8f688wzz3SbmpmVS563ifz66BBxySWXsGDBAi688MJah2Jmw8yAzwwkTZK0SNLP0/IRks6rfGjW04wZMzj22GOZMWNGrUMxs2EmzwPk/wCWA69Ny4+SNT6zKmtra2P58uW0tbXVOhQzG2byJIOJEXEtsAsgdVu9s/9DrBLmzJnDhAkTmDNnTq1DMbNhJk8y2CZpAtlDYyTNJBvPwKps6dKlbN68maVLl9Y6FDMbZvIkg88DNwJTJf0CWAKcX9GorFcTJ05EEhMnTqx1KGZ1yd3A923AFsgAkkYCh5MNbvNIRNS0t9JGbYE8evRoOjs7GTVqFNu3b691OGZ1p9GHjh3UsJeS+qqYPiyd0HUVVTZ37lyuuuoq5s6dW+tQzOqSh47tW3/tDN7Xz7Yga5FsVbRp0yYigk2bNtU6FLO6NGvWrIYsEeTRZzKIiHOrGYgNbM6cOdx9991+m8jMyi7PA2QbIq688ko2b97MlVdeWetQzGyYcTKoI88991y3qZlZuTgZ1JF99tmn29TMrFzy9E10p6SvSGqR5LtQDXkMZDOrlDwjnR0CvC19ZpKNaXBHRHyu8uH1GVNDtjMwMyvFoNoZdImI30t6EdiePu8A3ljeEM3MrJbyVBP9FvgJMAlYBPy3iGjJc3JJB0m6VdJDkh6U9Jm0/tWSVkh6RNJySeMKjvmapHWS7pd05OC+1vB0xRVXMHHiRK644opah2Jmw0yeB8hfAx4HzgA+A5wtaWrO8+8APh8R04FZwKckvQGYB9wcEYcDtwLzASSdAkyNiGnAJ4BvFfNlhrsFCxawefNmFixYUOtQzOqS+ybq24DJICL+PSI+CLwLuAdoJRvTYEAR8WRE3J/mnwceBg4CTgcWp90Wp2XSdEnafxUwTtKkvF9muDvmmGO6Tc2sOB4TpG8DPjOQdClwIrA3cBfwr8AdxV5I0uuAI4GVwKSI2AhZwii44R8IrC847Im0bmOx1xuObrnllm5TMyuO+ybq24DJAOgALum6eQ+GpL2BHwGfjYjnJfV8FcivBuUwefJk/vCHPzB58uRah2JWl9w3Ud/yJIOlwEckHRIRF0g6GDggIlbnuUDq/vpHwPcj4oa0eqOkSRGxUdIBwFNp/RNA4Z3uoLTuFVpbW1+eb25uprm5OU84dW3+/PksWLCA+fPn1zoUM6sD7e3ttLe359o3TzuDb5INefnOiHijpFcDKyLiLbkuIC0BNkXE5wvWXQxsiYiLJc0DXhUR8ySdCnwqIk5LI6pdHhEzezlnQ7YzmDlzJqtWreL4449n5cqVtQ7HzOpMSe0MgOMj4mhJ9wFExDOSRue88AnAmcCD6fgAFgAXA9dK+jjwGDA3nftnkk6V9BtgG+CeU83MqiBPMuiU1MTuMZD3IyspDCgifgE09bH5XX0c8+k8525EXaWhRiwVmVll5akmOhP4EHA02WugHwD+OSKuq3x4fcbUkNVE0u7SXSN+fzMrTX/VRHnHQH4DcDLZGMi3RMTD5Q2xOI2aDEaOHMnOnTtpampix44dtQ7HzOrMoJKBpH0j4llJ43vbHhFbyhhjURo1GUyYMIEtW7Ywfvx4Nm/eXOtwzKzODDYZ/DQi3ivp93RvByAgIuLQ8oeaT6MmA5cMzKwUg3qbKCLem6aHVCowK84xxxzD6tWr3R2FmZVdngfINwJXAzdExAtViWoALhm4ZGBmxeuvZJCn19JLyQa2eVjSjyR9QNKYskZouYwbN67b1MysXPIMbnM7cHtqa/BO4O+AK4F9Kxyb9bBz585uUzOzcslTMkDSnsD7gU8Cb2F399NWRYcffni3qZlZueR5ZnAtcBywDPghcHtE5GqBXCmN+szAjc7MrBSl9k20CDgjIlw3UWP77rsvzz77LPvu6xo6MyuvPNVEdwDzJV0BIGmapPdWNizrzYsvvthtamZWLnmSwfeA7cBb0/ITwL9VLCLr0/bt27tNzczKJU8ymBoRlwCdAKmtQa91TlZZTU1N3aZmZuWSJxlsT28TdXVhPRV4qaJRWa/8aqmZVUqeB8gLyd4kmizpP4ETgHMqGZT1Lr0J0O2tIjOzcuj31VJld52DgBeAmWTVQysjYlN1wuszroZ8tXTUqFHs2LGDkSNH0tnZWetwzKzOlDSegaQHI2JGRSIbpEZNBu6byMxKUWrfRPdKekuZY7JB8DMDM6uUPCWDXwOvJxu4fhu7xzN4U+XD6zMmlwxcMjCzIpXaAnl2meOxQRozZgzbtm1jzBh3GmtWaRHBjh076OzsZMeOHcNivj95ei19rGw/XSvJ9OnTWb16NdOnT691KNagum6Q/d10KnlDq9aNs7Ozk127djFy5MiXP6NGjXrFfG/rSp0fPXo0Y8eOrcj5p0yZ0ue/7YDVRENRo1YTuaO6oSki2LlzZ0VvWkPlBttVTTnYm1Elbp7Fnrsr/oH2b2pqGnavcZdaTWRWtMIb5FC9IZbrfDt27KCpqakmN8QxY8ZU9QY7HG+QlnEyqLKef0EWM7/HHnvw0ksvMXr0aK6//vohdUPs7XwjRoyoyV+FY8aMqdpfmV1/QY4YkWtoELMhq26TwV133VWxm14lb6zAoG8+L72U9QKyfft2lixZUtRNazB1kKXcpEeOHOkbpA05HR0dtLW1sXDhQmbNmlXrcIaUun1mMGvWrKrd3Mq5byk3SD8zMCtNS0sLy5cvZ/bs2SxbtqzW4VTdsHxmcNddd9U6BDOrMwsXLuw2td3qtmRQj3GXyn0TmVkpSu2OwoaIrgTYiInQzCrLyaCOuG8iM6sUJwMzM6tsMpC0SNJGSQ8UrHu1pBWSHpG0XNK4gm1fk7RO0v2SjqxkbGZmtlulSwbf45Ud3c0Dbo6Iw4FbgfkAkk4hG295GvAJ4FsVjs3MzJKKJoOIuBN4psfq04HFaX5xWu5avyQdtwoYJ2lSJeOrN11tFNyYy8zKrRZ3lf0jYiNARDwJdN3wDwTWF+z3RFpnya5du7pNzczKZSg0OhvUe5Ktra0vzzc3N9Pc3FymcMzMhof29nba29tz7VvxRmeSpgA3dY2MJulhoDkiNko6ALgtIt4o6Vtp/odpv18Db+8qRfQ4Z0M2OnN3FGZWilo3OlP6dLkROCfNnwPcULD+YwCSZgJbe0sEjaypqanb1MyK09HRQUtLCx0dHbUOZcipaMlA0lVAMzAB2AgsBH4CXAdMJhtXeW5EbE37fx1oIRtr+dyIuLeP87pk0IDf36xU7qiuRh3VRcRH+tj0rj72/3QFwzGzBueO6vrmdxTNzGxIvE1kZlYVbW1tLF++HKAhq4n642RgZg3D1UR983gGdcQPkM2sFLV+tdTMzIY4JwMzM3MyMDMzJwMzayBugdw3v01kZg3Dr5b2zcnAzBqGXy3tm18trSMjRowgIpDkMQ3MrGh+tXSY6EqAjZgIzayynAzqSFejs8LGZ2Zm5eBkUEdcMjCzSnEyMDMzJwMzM3MyqCt77rlnt6mZWbk4GdSR888/n1GjRnH++efXOhQzG2acDOrImjVr6OzsZM2aNbUOxcyGGbdAriNuPWlmleIWyGZmDcItkM3MrF9OBmZm5mRgZmZOBmZmhpOBmZnhZGBmZjgZmJkZTgZmZoaTgZmZ4WRgZmYMwWQgqUXSryU9KulLtY7HzKwRDKlkIGkE8HVgNjAdOEPSG2ob1dDS3t5e6xDM6pp/h3o3pJIBcBywLiIei4hO4Brg9BrHNKT4P7JZafw71LuhlgwOBNYXLP8xrTMzswoaasnAzMxqYEiNZyBpJtAaES1peR4QEXFxj/2GTtBmZnWkr/EMhloyaAIeAU4GNgCrgTMi4uGaBmZmNswNqWEvI2KnpE8DK8iqsBY5EZiZVd6QKhmYmVlt+AFynXBjPLPSSFokaaOkB2ody1DkZFAH3BjPrCy+R/Y7ZL1wMqgPboxnVqKIuBN4ptZxDFVOBvXBjfHMrKKcDMzMzMmgTjwBHFywfFBaZ2ZWFk4G9eGXwOslTZE0GvgwcGONYzKrR0of68HJoA5ExE6gqzHeQ8A1boxnVhxJVwF3AYdJelzSubWOaShxozMzM3PJwMzMnAzMzAwnAzMzw8nAzMxwMjAzM5wMzMwMJwNrMJLuLGLfEyX9StK9kvaoUDxvl3TTAPv8taQ1ktZK+nYl4jBzMrCGEhEnFrH7mcCFEXF0RLxUqZiAgRr7bALeGhFHANMlvbWCsViDcjKwhiLpuTQ9QNLt6a/+BySd0GO/84C5wAWSvp/WfUHSakn3S1qY1k2R9LCk70l6RNIPJJ0s6c60fGza7y2S7pJ0T9o2rZfYxqYBWFam/d4HWdfLEbEt7TYGeLFSPx9rXENqDGSzKuj6K/wjwLKIuEiSgLHddopYJOlE4KaIWCrp3cC0iDgu7X9j2r4emAq8PyLWSrobOCMiTpT0V8A/AX8DPAycGBG7JJ0MXAR8oEds/wTcEhHnSRoHrJZ0c0T8BUDSBcBvI+Le8v9YrNE5GVij+iWwSNIo4IaIWDPA/u8B3i3pXrKOzvYCppElg99HxNq030PALWn+QWBKmn8VsCSVCILef/feA7xP0hfT8miy3mofkfQmsgGNji3ua5rl42RgDSki7pB0EnAa8B+SLo2IH/RziICLIuI73VZKU4DC5wm7CpZ3sft37ALg1oiYk465rY9rvD8i1vWybQZwe0RsH+i7mQ2GnxlYoxGApIOBpyJiEfBd4OgBjlsOfFzSXun410rar/CcAxjH7jEo+uotcznwmZcDlY4s2PYLYFGO65gNipOBNZquZwbNwJpU7TMX+Pd+9iUi/gu4CuiQ9ABwHbB3z/3o+82gS4D/Jeke+v69uwAYlR5oPwh8uWDbDLJqJLOKcBfWZmbmkoGZmTkZmJkZTgZmZoaTgZmZ4WRgZmY4GZiZGU4GZmaGk4GZmQH/H/e/9bF/6RZtAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fed85eeb2b0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "xplot = numpy.arange(-1,3,1)\n",
    "yplot = [theta[0] + theta[1]*x for x in xplot]\n",
    "plt.plot(xplot, yplot, color='k')\n",
    "plt.scatter([[x[1] for x in X]], y, color='k', s=2)\n",
    "plt.xlim(-0.3, 1.3)\n",
    "plt.ylim(0, 550)\n",
    "plt.xticks([0,1])\n",
    "plt.xlabel(\"is female?\")\n",
    "plt.ylabel(r\"review length (words)\")\n",
    "plt.title(r\"gender vs.~length\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Transformation of output variables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dataset of reddit submissions (and resubmissions)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['#image_id',\n",
       " 'unixtime',\n",
       " 'rawtime',\n",
       " 'title',\n",
       " 'total_votes',\n",
       " 'reddit_id',\n",
       " 'number_of_upvotes',\n",
       " 'subreddit',\n",
       " 'number_of_downvotes',\n",
       " 'localtime',\n",
       " 'score',\n",
       " 'number_of_comments',\n",
       " 'username']"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f = open(dataDir + \"redditSubmissions.csv\")\n",
    "cs = csv.reader(f)\n",
    "header = next(cs)\n",
    "header"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute popularity as a function of submission number (i.e., how many times has this identical image been submitted)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "pops = defaultdict(list)\n",
    "imId = None\n",
    "count = 0\n",
    "for l in cs:\n",
    "    d = dict(zip(header, l))\n",
    "    if d['#image_id'] != imId:\n",
    "        count = 0\n",
    "    imId = d['#image_id']\n",
    "    count += 1\n",
    "    try:\n",
    "        pops[count].append(int(d['number_of_upvotes']))\n",
    "    except Exception as e:\n",
    "        continue\n",
    "\n",
    "for x in pops:\n",
    "    pops[x] = sum(pops[x]) / len(pops[x])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Extract a single feature which is just the submission number (0 for first submission, etc., 1 for first resubmission, etc.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = list(pops.keys())\n",
    "x.sort()\n",
    "X = [[1,a] for a in x]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The label to predict is the number of upvotes as a function of submission number"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = [pops[a] for a in x]\n",
    "ylog = [math.log2(a) for a in y]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fit two models: one with the original output variable, one with a log-transformed output variable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "mod = linear_model.LinearRegression(fit_intercept = False)\n",
    "modlog = linear_model.LinearRegression(fit_intercept = False)\n",
    "\n",
    "mod.fit(X,y)\n",
    "modlog.fit(X,ylog)\n",
    "\n",
    "theta = mod.coef_\n",
    "thetalog = modlog.coef_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot data and fitted values for the two models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEZCAYAAABvpam5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VOW9x/HPDwj7KotJcMGyW2UREVtpBQVlU6ulgBZF2ltrr9dr6723Fbtgba1626t2tWoVEbRsVrFSNXWJNWrVIqsgEFAEJIAKIglLSH73j3MSh5iEmWROJpl836/XvGbmmXPOPGdOMt95nucs5u6IiIgkokmqKyAiIg2PwkNERBKm8BARkYQpPEREJGEKDxERSZjCQ0REEqbwkEiZ2Ttmdk6SljXDzO6txfyXmdnTyaiLSGNnOs5DomRm7wDfdPfnU10XSS0zKwV6ufumVNdFak8tDxGpK/qlmkYUHhIXM/uBmW01s71mttbMRobls8zs5pjpzjazLRVmP8PM3jKzD83sfjNrHjutmf2Pme0ws21mdpGZjTWzdWb2gZnNiFn2TDObEz5uYWZzwml2m9lrZtY1fO1KM9sY1nWjmV0alk8zs5dilvdFM3s9Zv4vxLz2gpndbGZ54XKeNrNjqvhsytbj+pj1uDLm9fZm9pCZ7Qy78X4Y52c+yczeqFD2PTN7PHw8Lvxc95a9f5zLfcHMvhHzvOLnUmpm14af3U4z+9+wvHn4WZ0cM20XMysysy7h82+Z2YZwuzxuZplh+YuAASvD+n4tLJ9gZsvC5eaZ2akxy670b07qB4WHHJWZ9QGuAYa4e3vgfODdamap+AvzMmA00BPoC/wo5rVMoDmQDcwE7gO+DgwGvgz82MxOrGTZ04D2QHfgGOBqYL+ZtQZ+DZwf1vWLwPKK85tZJ+BJ4C6gM3AnsCQsL3Np+D5dgRbAf1ezzplAu3A9/g34vZl1CF/7XfhaD2AEcIWZTa9mWWX+CvQxs54V6vRw+PhPwLfC9TwFqE3XYMVt9hXgtPB2kZl9w90PAY+GdSgzCch19w/Csa1fABOBLOA9YD6Au58dTn+qu7d394VmNhi4H/gWwTa8B3jCzDJq8DcndUzhIfEoIfiCP8XMmrn7e+7+TgLz/9bd33f3PcAtHPnlcwj4hbuXAPOALsBd7l7k7muANcDASpZZTPCl38cDy9x9X0x9TzWzlu6+w93XVjL/eGC9uz/i7qXuPg94G7ggZppZ7r7R3Q8CC4BB1azjIeBn7l7i7k8B+4C+ZtYEmAzcEK7TZuD/gMurWRYA7r4fWEz4eZlZb4LwfSLmPT9vZu3c/WN3X175kmrktnCZWwkCtmyb/Zkjt99lfBpmlwH3u/sKdy8GZgBfMLMTYqa3mMffAv7o7v8Kt+Ec4CBwJrX/m5OIKTzkqNx9I/Bd4CZgh5k9UtYdEaetMY83E/w6L/Ohf7rXxv7wfmfM6/uBtpUscw7wDDAv7Nq4zcyaunsRwZf1d4DtZvZXM+tbyfzZYV1ibSZoyZQpiHlcVEU9YtejtJLpuwDNCH6FV/U+1Yn9sr4MeDwMM4CvEoTg5rAr6sw4lxmPqrbZC0ArMxsatggHAo+Frx3xmbp7IfAhVa/ricB/mdlH4W03cByQXcXfXFZyVk2SQeEhcXH3ee7+JYJ/eIDbw/tCoHXMpJX9gx8f8/hE4P0k1Oewu//M3T9P0DV1AXBF+Nrf3f08gq6kdUBlu/e+T9CNFOsEYFtt61bBBwStpNiutxMTeJ+/A13NbCAwBXik7AV3X+ruXyHoVltM0DqKR8VtVtkPgdhtdgLhNgsDcgFBkF0KPBkGNuE05etpZm0IWoexQRRrC3CLux8T3jq5e1t3L+vqqvg3d1uc6yd1QOEhR2VmfcxspAUD3YcIWgNlv7KXA+PMrFPYGrmukkVcY2bdwwHnGwm6p2pbpxFmdkrYLbSP4Au61My6mdmF4dhHcfhaaSWL+BvQ28ymmFlTM5sM9CcYZ0ia8Mt2IXCLmbUNf61/j6DlFM/8h8P5fwl0IggTwnGBy8ysfdjl9wlBV088lgOXmFkrM+sFfLOSaf7HzDqa2fEE2zR2m/2ZoHV3GTFhFpZPN7MBZtaCYPzjn+5etgNFAfC5mOnvA642szPCdWoT7gTQ5ih/c1IPKDwkHi0IfvXtIvh12ZWgPxuCL8GVBIOZT/PZYHCCL5gcIB/YQDDuUZWKA7dV7d6ZCSwCPgbeIuhOmUPwN309wS/7DwgG3b/zmTdx/wiYQDAI/kF4P97ddx/lfeMVO/+1BN1Ym4B/AHPdfRaAmQ03s71HWdafgXOBBRW6xi4H3jGzPcBVBF/mmNnx4R5Kx1WxvDsJgrUAmAXMrWSaxcBS4E2CQH2gfMXcXydovWQBT8WUPwf8GPgLwed/EkFrqcxNwENhF9VEd19KMO7xOzP7CFhPsIMCVP83J/VApAcJhr8+/kEw8NUMWOTuPzWzHgRfMscQ/IFe7u6Hw18ZDwFDCP6hJ7v7e+GyZgDfAA4D17l7TmQVF2nETAfzSRwibXmEA3sj3X0wwZ4qY81sGEF/+f+5ex9gD582m78JfOTuvQn28Cjbv/xkgl0C+wNjgT+YmSEiIikRebdVzGBaC4LWhwMjCfYXB5hNsE85wEXhcwi6JMrOiXQhMC8cJH2XoOvjjGhrLtJo6UhwOarIw8PMmpjZMoL+1b8DG4E9MX23W/l0V77uBHtgEA4CfhwOspaXh7YR/66OIpIAd2+qLis5mrpoeZSG3VbHEbQW+iUwu7qmRETqoWZ19UbuvtfMcoEvAB3NrEnY+jiOT/d530awf/n7ZtYUaO/uH5lZWXmZ2HnKmZma2yIiNeDuCf1Yj7TlEZ40rUP4uBXB+Y3WEOxW+bVwsmkEuwVCcNqFsl31vsan5+p5AphiwYnZTgJ6Aa9X9p75+fm4e4O7zZw5M+V10DpqPbWejW8d3Wv2mzvqlkcWMDs8kKsJMN/d/2ZmawlOK/EzYBnBydEI7+eY2QaC0xpMAXD3NWa2gCB4ioF/9yrWePv27fTs2bOyl0REJEkiDQ93X0VwVs6K5e8AwyopP0iwS25ly7oVuPVo7/n++7U+84WIiBxF2h1hvn379lRXoUZGjBiR6ipErjGsI2g900ljWMeaSqvL0JqZ/+AHP+C223T+NBGReJkZXp8GzFNB3VYiItFLu/BoqN1WIiINSdqFh1oeIiLRS7vwUMtDRCR6aRcehYWFHDhwINXVEBFJa2kXHpmZmWp9iIhELO3CIzs7W+EhIhKxtAuPrKwsDZqLiEQsLcNDLQ8RkWilXXio20pEJHppFx7qthIRiV7ahYdaHiIi0Uu78FDLQ0QkemkZHmp5iIhEK+1OyV5SUkLLli355JNPaNGiRaqrJCJS7+mU7ECTJk049thjKSgoSHVVRETSVtqFB2jQXEQkamkZHho0FxGJVtqGh1oeIiLRScvwULeViEi00jI81G0lIhKttAwPtTxERKKVluGhloeISLTSNjzU8hARiU7aHWHu7pSUlNCqVSsKCwvJyMhIdbVEROo1HWEeatq0KV27dtVR5iIiEUnL8AANmouIRCltw0OD5iIi0Unr8FDLQ0QkGmkbHuq2EhGJTtqGh7qtRESik9bhoZaHiEg00jY8srOz1fIQEYlI2oaHWh4iItFJyyPMAQ4fPkzr1q0pKiqiWbNmKa6ZiEj9pSPMYzRr1ozOnTuzY8eOVFdFRCTtpG14gLquRESiktbhoUFzEZFopHV4qOUhIhKNtA4PHWUuIhKNtA4PHWUuIhKNtA8PtTxERJIvrcNDA+YiItFI6/BQy0NEJBppe4Q5QHFxMW3atGH//v00bdo0hTUTEam/dIR5BRkZGXTq1ImdO3emuioiImklrcMD1HUlIhKFtA8PDZqLiCRf2oeHWh4iIsmX9uGho8xFRJIv7cNDR5mLiCRfowgPtTxERJIr7cNDA+YiIsmX9uGhloeISPKl9RHmAIcOHaJt27YcOHCAJk3SPitFRBKmI8wr0bx5czp06MCuXbtSXRURkbSR9uEB6roSEUm2RhEeGjQXEUmuRhEeanmIiCSXwkNERBLWKMJD3VYiIsnVKMJDLQ8RkeRqFOGhloeISHI1ivBQy0NEJLkiDQ8zO87Mnjezt8xslZn9Z1jeycxyzGydmT1jZh1i5vmNmW0ws+VmNiimfJqZrQ/nuSKRemRmZlJQUEBpaWnyVk5EpBGLuuVxGLje3T8PfAG4xsz6ATcAz7p7X+B5YAaAmY0Ferp7b+DbwB/D8k7AT4ChwDBgZmzgHE3Lli1p164dH374YfLWTESkEYs0PNy9wN2Xh4/3AWuB44CLgNnhZLPD54T3D4XTvwZ0MLNjgfOBHHf/2N33ADnAmETqoq4rEZHkqbMxDzPrAQwC/gkc6+47IAgY4Nhwsu7AlpjZtoZlFcu3hWVx06C5iEjyNKuLNzGztsAi4Dp332dmFU/lW9WpfRM6yyPATTfdVP54xIgRjBgxAlDLQ0SkTG5uLrm5ubVaRuThYWbNCIJjjrsvDot3mNmx7r7DzDKBnWH5NuD4mNmPC8u2ASMqlL9Q2fvFhkcshYeISCD2hzXAT3/604SXURfdVg8Aa9z91zFlTwBXho+vBBbHlF8BYGZnAnvC7q1ngNFm1iEcPB8dlsVN3VYiIskTacvDzM4Cvg6sMrNlBN1TNwK3AwvM7BvAZmASgLv/zczGmVk+UAhMD8t3m9nPgH+Fy/hpOHAet6ysLF588cUkrZmISOMWaXi4+8tA0ypeHlXFPP9RRfmDwIM1rUtWVpZaHiIiSdIojjCHoNtKYx4iIsmR9tcwL7N//346duzIgQMHMEt4Jy4RkbSla5hXo1WrVrRu3ZqPPvoo1VUREWnwGk14gHbXFRFJlkYXHho0FxGpvUYVHho0FxFJjkYVHuq2EhFJjkYVHjrKXEQkORpVeKjlISKSHAoPERFJWKMKD3VbiYgkR6M5whygsLCQLl26UFRUpKPMRURCOsL8KNq0aUPz5s3ZsyehE/KKiEgFjSo8QOMeIiLJoPAQEZGEJRweZtbJzAZEUZm6oEFzEZHaiys8zCzXzNqb2THAm8B9ZnZHtFWLhloeIiK1F2/Lo4O77wUuAR5y92FUcSXA+k4tDxGR2os3PJqZWRbBtcafjLA+kVPLQ0Sk9uINj5uBZ4CN7v6GmX0O2BBdtaKj8BARqb1m8Uzk7guBhTHPNwFfjapSUVK3lYhI7cU7YN7HzJ4zs9Xh8wFm9qNoqxaNspZHOh1ZLyJS1+LttroPmAEUA7j7SmBKVJWKUrt27TAz9u7dm+qqiIg0WPGGR2t3f71C2eFkV6au6IqCIiK1E294fGBmPQEHMLOJQIP99tWguYhI7cQ1YA5cA9wL9DOzbcA7wNcjq1XENGguIlI78YaHu/soM2sDNHH3T8zspCgrFiW1PEREaifebqtHAdy90N0/CcsWRVOl6GVlZanlISJSC9W2PMysH/B5oIOZXRLzUnugZZQVi1J2djZvvvlmqqshItJgHa3bqi8wAegIXBBT/gnwragqFTV1W4mI1E614eHui4HFZvYFd3+1juoUOQ2Yi4jUTrxjHlvM7DEz2xneHjWz4yKtWYTU8hARqZ14w2MW8ASQHd7+GpY1SO3bt6ekpIRPPvnk6BOLiMhnxBse3dx9lrsfDm8PAl0jrFekzExHmYuI1EIiR5hPNbOm4W0q8GGUFYuauq5ERGou3vD4BsGFoAoITksyEZgeVaXqQnZ2Nlu2bEl1NUREGqR4w6PI3S90967u3s3dv+Lu70Vas4h98YtfJDc3N9XVEBFpkCye61qY2XrgXWA+8Ki774m4XjViZh7vdTry8/P58pe/zNatW2nSJN4MFRFJP2aGu1si88T1renufYAfERxt/qaZPRmOezRYvXr1ol27dixbtizVVRERaXDi/snt7q+7+/XAGcBHwOzIalVHxo8fz5IlS1JdDRGRBifey9C2N7NpZvYU8ArBoPkZkdasDig8RERqJt4xj3eAx4EF9fk0JYmMeQAcOnSIbt26sX79erp16xZhzURE6q/IxjyAz7n794C3zKxd4lWrn5o3b86oUaN46qmnUl0VEZEGJd7wGGJmq4CVwGozW2FmQyKsV51R15WISOLi7bZaCVzj7i+Fz4cDf3D3ARHXLyGJdlsBFBQU0L9/f3bu3ElGRkZENRMRqb+i7LYqKQsOAHfPAw4n8kb1VWZmJr169SIvLy/VVRERaTDiDY8XzeweMxthZmeb2R+AXDM7zcxOi7KCdUFdVyIiiYm32+qF8GHZxBY+NsDd/ZxoqpeYmnRbAbzxxhtcccUVrF27NoJaiYjUbzXptjraZWjL5FZ47gDufnMib1ZfDRkyhN27d7Np0yY+97nPpbo6IiL1XrzdVvtiboeBsUCPiOpU55o0acK4cePUdSUiEqe4uq0+M5NZC+AZdx+R9BrVQk27rQAeffRR7rvvPp5++ukk10pEpH6rSbdVTcOjE/CGu/dKeOYI1SY89u7dS/fu3dm+fTtt27ZNcs1EROqvyHbVNbNVZrYyvL0FrAPuqkkl66v27dtzxhln8Nxzz6W6KiIi9V68A+YTYh4fBna4e1oc5xGrbJfdiy66KNVVERGp12rUbVVf1abbCmDdunWce+65bNmyBbOEWnAiIg1WlEeYNwp9+vShZcuWrFixItVVERGp1xQeMcxMR5uLiMRB4VHBhAkTFB4iIkehMY8KDh48SLdu3cjPz6dr165JqpmISP2lMY8kaNGiBeecc44OFhQRqYbCoxIa9xARqZ66rSrx/vvvc8opp7Bz506aNYv3UBgRkYZJ3VZJkp2dTY8ePXjllVdSXRURkXpJ4VEFdV2JiFRN4VGF8ePH8+STT6a6GiIi9VKk4WFm95vZDjNbGVPWycxyzGydmT1jZh1iXvuNmW0ws+VmNiimfJqZrQ/nuSLKOpcZOnQou3bt4t13362LtxMRaVCibnnMAs6vUHYD8Ky79wWeB2YAmNlYoKe79wa+DfwxLO8E/AQYCgwDZsYGTlSaNm3K2LFj1XUlIlKJSMPD3fOA3RWKLwJmh49nh8/Lyh8K53sN6GBmxxKET467f+zue4AcYEyU9S6jcQ8RkcqlYsyjm7vvAHD3AuDYsLw7sCVmuq1hWcXybWFZ5M477zzy8vIoKiqqi7cTEWkw6sNBDFUdmFGjc6LfdNNN5Y9HjBjBiBEjarIYADp27Mhpp53G888/z4QJE44+g4hIA5Cbm0tubm6tlhH5QYJmdiLwV3cfED5fC4xw9x1mlgm84O79zeyP4eP54XRvA2cDI8Pprw7Lj5iuwnsl5SDBWL/85S/ZtGkTd999d1KXKyJSX9TXgwSNI1sRTwBXho+vBBbHlF8BYGZnAnvC7q1ngNFm1iEcPB8dltWJsl120+lIfBGR2op6V91HgFeAPmb2nplNB24jCIN1wDnhc9z9b8A7ZpYP3AP8e1i+G/gZ8C/gNeCn4cB5nejfvz8ZGRmsWrWqrt5SRKTe07mt4nDttdeSnZ3NjBkzkr5sEZFUq6/dVg3e+PHjWbBgARs2bFD3lYgI9WNvq3pv5MiRDBs2jJEjR5KRkcHo0aM577zzOOecczjmmGNSXT0RkTqnbqsEuDtr164lJyeHv//977z00kv079+f8847j9GjR3PmmWfSvHnzyN5fRCQKNem2UnjUwsGDB3n11VfJyckhJyeHDRs2cPbZZzN69GgmTpxIVlZWndVFRKSmFB51HB4VffDBBzz33HM888wzLF68mOnTp3PDDTfQpUuXlNVJRORoNGCeYl26dGHy5Mk88MADrFq1iv3799O3b19+8pOfsGdPne1dLCISOYVHRLKzs/n973/P0qVL2bp1K7179+bWW2+lsLAw1VUTEak1hUfEevTowQMPPEBeXh4rVqygZ8+e3HXXXRw4cCDVVRMRqTGFRx3p27cv8+bNIycnhxdeeIHevXtz7733UlxcnOqqiYgkTOFRxwYMGMDixYtZtGgRCxcupF+/fsyePVstERFpULS3VYrl5uZyyy23sGzZMiZPnsz06dMZMmQIZjU6I72ISMK0q24DDI8ymzdv5qGHHuLBBx+kVatWXHnllUydOpXMzMxUV01E0pzCowGHRxl3Jy8vj1mzZvHYY49x1llnMX36dCZMmECLFi1SXT0RSUMKjzQIj1j79u3jL3/5C7NmzWL16tVMmTKF6dOnM3jwYHVriUjSKDzSLDxivfPOO+XdWieddBKzZs3ixBNPTHW1RCQN6AjzNHbSSScxc+ZM8vPzOf/88zn99NN58MEHdYp4EUkJtTwaqJUrVzJ16lR69uzJPffcQ7du3VJdJRFpoNTyaEQGDBjAG2+8Qd++fRk4cCCLFy8++kwiIkmilkcayMvLY9q0aZx99tncddddtG/fPtVVEpEGRC2PRmr48OEsX76cZs2aMXDgQF588cVUV0lE0pxaHmlmyZIlXHXVVVx66aX8/Oc/p2XLlqmukojUc2p5COPHj2fFihVs3ryZIUOGsGzZslRXSUTSkFoeacrdeeSRR7j22mvZvXs3EPy6KDu4sOxxxefNmjWjTZs25be2bdse8Ty2vH379lx44YX069cvZespIrWngwQVHp9R9nm4e/kt9nnFx8XFxRQWFh5x27dv32fKCgsL2bVrF/PmzWPo0KFcf/31jBw5Uke+izRACg+FR53bv38/Dz/8MHfeeScZGRl873vfY8qUKToPl0gDovBQeKSMu5OTk8Mdd9zBqlWruOaaa/j2t79Nly5dElpGQUEBy5YtY9WqVfTv35/zzz9fQSQSMYWHwqNeWL16NXfeeSePPfYYkyZN4rvf/e5nxkVKSkrIz89n+fLlLFu2rPz+8OHDDB48mFNOOYU333yT1atXc+GFFzJ58mRGjRpFRkZGitZKJH0pPBQe9cqOHTu4++67ufvuuxk6dChjxoxh7dq15S2Lrl27MmjQIAYPHlx+37179yPGTbZt28bChQuZP38+GzZs4OKLL2by5MmMGDGCZs2apXDtRNKHwsPM58yZQ9OmTWnatClNmjQ54j72Fk9ZbebTwPGnysZFXn31VU499VQGDx7MwIED6dixY0LL2bx5c3mQbN68ma9+9atMnjyZL33pSzRt2jSi2oukP4WHma9fv56SkhJKS0spKSkpv1V8frSyeKavaprS0tIqgyvRQIo65Bpq6G3cuJEFCxawYMECduzYwSWXXMK4ceMYMWIErVu3TnX1RBoUhUc96bZy91qFT23miyr0ahJIUYVcxbINGzbw+OOP89RTT/Hmm29y1llnMXbsWMaOHUvv3r0bTCCKpIrCo56ER0MVb+hFEYS1mb9i6B0+fJiDBw+yf/9+SktLadOmDe3ataNjx45kZGTUm9ZeQ2rpSXpTeCg8GiV3rzJs3n77bf7xj3/w8ssvs27dOgYMGMCwYcPo168f3bt3p2XLlrUOy5oEacXQa0itPYVe+lF4KDykGnv37uW5557j6aefZunSpaxdu5ZOnTrRv39/Tj755CNunTt3rnZZhw4dYtu2bbz33nts2bKl/H7Lli3079+f66+/nqysrGqXUVXoRdn6q23olZSU4O41CqlEQi+KIFToVU3hofCQBJSWlvLee++xZs0a1q5dy5o1a8pvLVu2LA+SE044gYKCgvJweO+99/jwww/Jzs7m+OOP5/jjj+eEE07g+OOP57jjjuP5559nzpw5TJkyhe9///v06NEj1auaVNW19GobaDUJvXiXWdPQS7R1l+zWXl2EnsJD4SFJ4O68//775UGyZcsWsrKyjgiKzMxMmjatevfgnTt3ctddd3HPPfcwYcIEZsyYoRNIpliyQy+e+aIKvWR3aX75y19WeKTT+kjDt2fPHn73u9/xm9/8hrPPPpsbb7yRwYMHxz3/1q1beemll8jLyyMvL48WLVpw2WWXMWXKFDIzMyOsudQXiYReTYNw9OjRCo90Wh9JH4WFhdx777386le/YuDAgfzwhz/krLPOOmKa0tJS1qxZUx4UeXl5FBUVMXz48PLbxx9/zNy5c1m8eDFnnnkml19+OV/5yldo06ZNitZM0oG6rRQeUs8dPHiQBx98kNtvv50TTjiBq6++mnfffZe8vDxeeeUVOnfufERY9OnTp9L+7sLCQp544gnmzp3Lyy+/zAUXXMDUqVM599xz6+VpW9xdg9X1mMJD4SENxOHDh5k3bx5z586lf//+DB8+nLPOOqtGXVE7d+5kwYIFzJkzh82bN3PppZcydepUTjvttJR+YW/atIn58+czf/58Nm3axOjRoxk/fjzjxo1Tl1s9o/BQeEgjt379eh5++GHmzp1L8+bNmTRpEhMnTuSUU06pkyDZsmULCxYsYN68eUecf6xv377k5OSwZMkScnJy6N27N+PHj2f8+PEMGTKEJk10RexY7s6uXbvYuHEj+fn55Ofnc+KJJ3L55ZdHcmZphYfCQwQIvnxee+01Fi1axKJFi2jRogUTJ05k4sSJDBo0KKlBsn379vITVr799tvlZz4eOXJkpV1oxcXF5OXlsWTJEpYsWcLu3bsZN24c48ePZ/To0bRv3/6o71lSUkJRURFFRUU0adKELl26NLhusdLSUrZv314eDrFBkZ+fT0ZGBr169aJXr1707NmTl19+mS1btnDzzTczadKkpAauwkPhIfIZ7s7SpUvLg8Tdy4Pk9NNPT/hLt+xL769//Svz589n+fLlXHDBBUyZMoVRo0bRvHnzhJa3cePG8iB55ZVXGDp0KJ07dz7iksdFRUVHPD948CCtWrWiTZs2FBcXk5GRwcCBA4+49e/fP2XXfzl48CDbt29n27ZtR9zef//9Ix63b9++PCBibz179qRTp06fWe6zzz7LjTfeSHFxMbfccgtjx45NSmgqPBQeItVyd1asWMGiRYtYuHAhBw4cOKJFsmPHDgoKCti+fXv5fezjgoICdu7cSYcOHRg1ahSTJ09mzJgxtGzZMin127dvH7m5uRQVFdGmTRtat25NmzZtjri1bt2aVq1alf/ydne2bdvGihUrjrht3ryZfv36fSZUjnb2gJrIz886flbkAAAKC0lEQVTnjjvu4JVXXmHbtm3s3buXzMxMsrOz6d69e/mt4vOa7CXn7jz++OP86Ec/4phjjuHWW29l+PDhtaq/wkPhIRI3d+ett94qb5Fs2LCBzMxMMjMzycrKqvL+2GOPTbh1kQpFRUWsXr36iEBZuXIlJ598MldffTWTJk2iVatWtXqPZcuWcfvtt/Pss8/yne98h4svvpju3bvTtWvXyMdxSkpKmDt3LjNnzuTkk0/mF7/4BYMGDarRsmoSHrh72tyC1RGRmigtLU11FSJXXFzsTzzxhI8bN847d+7s1113na9ZsyahZZSWlnpubq6PGTPGs7Oz/Ve/+pXv3bs3ohof3YEDB/y3v/2tZ2Zm+uTJk339+vUJLyP87kzo+1YtDxFplN59913+9Kc/cf/999OnTx+uvvpqLrnkElq0aFHp9KWlpTz55JPcdtttfPDBB3z/+9/n8ssvr3L6ulZYWMivf/1r7rzzTsaPH89JJ51U9qP6qPe33HKLuq3SaX1EJHrFxcUsXryYe+65h5UrVzJt2jSuuuoqevXqVf76vHnzuP3222nRogUzZszg4osvrvbcZqn00UcfMWvWLPbu3QtQPqBe3f2Pf/xjhUc6rY+I1K0NGzZw7733Mnv2bAYNGsSXvvQl7r//fnr27MkNN9zAqFGjGtwuwfHQgLnCQ0SS4MCBA/zlL3/hpZde4sorr2TYsGGprlKkFB4KDxGRhNUkPHROABERSZjCQ0REEqbwEBGRhCk8REQkYQoPERFJmMJDREQSpvAQEZGEKTxERCRhCg8REUmYwkNERBKm8BARkYQpPEREJGEKDxERSZjCQ0REEtagwsPMxpjZ22a23sx+kOr6iIg0Vg0mPMysCfA74Hzg88ClZtYvtbVKntzc3FRXIXKNYR1B65lOGsM61lSDCQ/gDGCDu29292JgHnBRiuuUNI3hj7QxrCNoPdNJY1jHmmpI4dEd2BLzfGtYJiIidawhhYeIiNQTDeYa5mZ2JnCTu48Jn98AuLvfHjNNw1gZEZF6JtFrmDek8GgKrAPOBbYDrwOXuvvalFZMRKQRapbqCsTL3UvM7D+AHILutvsVHCIiqdFgWh4iIlJ/aMA8xczsXTNbYWbLzOz1VNcnWczsfjPbYWYrY8o6mVmOma0zs2fMrEMq65gMVaznTDPbamZvhrcxqaxjbZnZcWb2vJm9ZWarzOw/w/K02p6VrOe1YXm6bc8WZvZa+J2zysxmhuU9zOyf4UHYfzazanum1PJIMTPbBAxx992prksymdlwYB/wkLsPCMtuBz509/8NzxDQyd1vSGU9a6uK9ZwJfOLud6S0ckliZplAprsvN7O2wFKCY6ymk0bbs5r1nEwabU8AM2vt7kXhWPLLwHXA9cAid19oZncDy939nqqWoZZH6hlpuB3cPQ+oGIgXAbPDx7OBr9RppSJQxXpCsF3TgrsXuPvy8PE+YC1wHGm2PatYz7JjydJmewK4e1H4sAXB2LcDI4FHw/LZwMXVLSPtvrQaIAeeMbM3zOxbqa5MxLq5+w4I/lGBbimuT5SuMbPlZvanht6dE8vMegCDgH8Cx6br9oxZz9fCorTanmbWxMyWAQXA34GNwB53Lw0n2QpkV7cMhUfqneXupwPjCP5Ah6e6QnUoXftM/wD0dPdBBP+cadHdEXblLAKuC3+ZV9x+abE9K1nPtNue7l7q7oMJWpBnAAmfJ1DhkWLuvj283wU8RrAh09UOMzsWyvuXd6a4PpFw913+6WDifcDQVNYnGcLB00XAHHdfHBan3fasbD3TcXuWcfe9QC7wBaBjeAJaCEJlW3XzKjxSyMxah79yMLM2wHnA6tTWKqmMI/uKnwCuDB9PAxZXnKGBOmI9wy/SMpeQHtv0AWCNu/86piwdt+dn1jPdtqeZdSnrejOzVsBoYA3wAvC1cLKjbk/tbZVCZnYSQWvDCQatHnb321Jbq+Qws0eAEUBnYAcwE3gcWAgcD2wGJrn7nlTVMRmqWM+RBP3lpcC7wLfLxgYaIjM7C/gHsIrgb9WBGwnO8rCANNme1aznZaTX9jyVYEC8SXib7+63hN9H84BOwDJgangG88qXo/AQEZFEqdtKREQSpvAQEZGEKTxERCRhCg8REUmYwkNERBKm8BARkYQpPCQthafRvr4W83/bzKYmOE+WmS2o6XvWhdp+LiJlGsyVBEXqUnWnoq5mnu3ApAiqU2+YmbkODhPU8pAGIjyVy5PhBWxWmtnXwvJ3zOyY8PEQM3shZrZBZvZKeLGifwunOdvMcs3scTPLN7Nbzeyy8OI4K8KjbI/4hW5m/xleIGh5eER52XKWhRcHWmpmbczsRDNbFb7ewsweCOu61MxGhOXTzOxRM3sqrNftVazvO2Z2UzjvCjPrU7Fe4fNVZnZC+N5rzWxWuNy5ZnaumeWFz0+v7nMJl/XfZvZ6uJ5lFwg60czeNrPZ4bodV5vtKOlDLQ9pKMYA29x9AoCZtQvLqzuz66nAMKAdsMzMngzLBxCcRXQPsAm4z92HWXCFvGsJLooT6wdAD3cvNrP2Ydl/Af/u7q+aWWvgQIX3vwYodfcBZtYXyDGz3uFrAwlOd1EMrDOz37h7ZSeh2+nuQ8zsO8B/A1dVMk3s+vYEvurua8zsX8Cl7j7czC4Efsin12eo7HM5Fejt7meYmQFPhGd43gL0Ai539zcqeX9ppNTykIZiFTA6bCkMd/dPwvLqLtKz2N0PufuHwPN8esbiN9x9p7sfIriOQU7Me/SoZDkrgEfM7OtASVj2MnCnBZcq7RRzHYQyw4G5AO6+juCcSH3C155z933ufpDghHQnVlH/x8L7pVXUC45c/3fcfU34+C3guZj1in2Pyj6X8wg+3zeBN4G+QFnYbVZwSEUKD2kQ3H0DcBrBF+HPzexH4UuH+fTvuGXF2WIeW8zzgzHlpTHPS6m8NT4e+F34/m+YWRN3vx34JtAKeLmsW6kasV/yse9fUsV7xk4XO03s+sKR6xzvelX1udzq7qe5+2B37+Pus8LywirqJ42YwkMaBDPLAva7+yPALwm+yAHeAYaEj79aYbaLzKy5mXUGzgZq+uv5BHd/EbgBaA+0NbPPuftb7v6/4XIrXkznJeDrYd37EJx5dl0N3z/Wu4TrbmanASfFvBbvpVIr+1xygG+ElwbAzLLNrGuCy5VGRGMe0lCcCvzSzEqBQ8B3wvKbgfvN7GOCi9rEWhmWdQZudveCcPwhVrV7DllwcaC54ViHAb92971m9nMzG0nQKngLeIojL9v5B+BuM1tJMLYxLRwzqfgWVb1/VeWPAleEg9evcWQgeRWPK/rM5wIUmFk/4NWwjp8AUwlaLdq7Sj5Dp2QXEZGEqdtKREQSpvAQEZGEKTxERCRhCg8REUmYwkNERBKm8BARkYQpPEREJGEKDxERSdj/A3azVl6h8CzEAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fed85ecdd30>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(x,y,color='k')\n",
    "plt.xlim(1,30)\n",
    "plt.ylim(0,3000)\n",
    "x1 = [0,30]\n",
    "y1 = [theta[0] + theta[1]*a for a in x1]\n",
    "plt.plot(x1,y1,color='grey')\n",
    "plt.xlabel(\"submission number\")\n",
    "plt.ylabel(\"upvotes\")\n",
    "plt.yticks([0,1000,2000,3000])\n",
    "plt.title(\"submission no. vs. upvotes\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAEZCAYAAAB1mUk3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmcFPWd//HXB4HhUgPIIQwMghEwAoZDQPFnR66EiIhjPAgalXVdN+sV1wRcN6Ab4/XYhPhL3GS9orKoQQwkBgQ8RtTlkOEUEVSO4T6DwIAMM/PZP6pm7Blnhu6Z6emZnvfz8ejHVFdXV32qC+rd9f1WV5m7IyIi9VuDZBcgIiLJpzAQERGFgYiIKAxERASFgYiIoDAQEREUBhInM9tkZpdW07wmmdl/V+H948zsjeqoRaS+M/3OQOJhZpuACe7+drJrkeQys0LgbHffmOxapOp0ZCAilaVvkilEYVBPmdnPzGybmR0ys3Vm9p1w/HNm9mDUdJeY2dZSb7/AzNaa2X4ze8bMGkdPa2b3mtluM9tuZmPM7Htmtt7M9pnZpKh5TzazF8PhNDN7MZzm72a2xMzahK/daGafh7V+bmbXheN/ZGbvRc3vQjNbGvX+wVGvvWNmD5rZ++F83jCzVuV8NkXr8ZOo9bgx6vXTzOwFM9sTNpv9W4yf+dVm9mGpcXeb2axweFT4uR4qWn6M833HzG6Oel76cyk0s9vDz26PmT0Wjm8cflbnRk17hpkdNbMzwue3mNmn4XaZZWbtw/HvAgasDuv9QTj+MjNbEc73fTPrFTXvMv/NSe2gMKiHzOwc4MdAP3c/DRgJbK7gLaW/AY4DhgPdgO7A/VGvtQcaAx2AycBTwA+BbwP/D/h3M8soY94/Ak4DOgKtgH8CjplZM+A3wMiw1guBlaXfb2YtgdeBqUBr4NfA38LxRa4Ll9MGSAP+tYJ1bg+cGq7HPwC/M7PTw9d+G77WBYgAN5jZTRXMq8hfgXPMrFupmv4nHH4auCVcz/OAqjTFld5mVwB9w8cYM7vZ3fOAmWENRa4Gstx9X9g39EvgKuBMIAd4BcDdLwmn7+Xup7n7DDP7NvAMcAvBNvwD8Bcza1SJf3NSwxQG9VMBwQ77PDNr6O457r4pjvf/f3ff4e4HgYcouTPJA37p7gXAy8AZwFR3P+ruHwMfA33KmOcJgp34OR5Y4e5HourtZWZN3H23u68r4/3fBza4+3R3L3T3l4FPgNFR0zzn7p+7+3HgT8D5FaxjHvAf7l7g7nOBI0B3M2sAXANMDNdpC/CfwPUVzAsAdz8GzCb8vMzsmwRh+peoZX7LzE519y/cfWXZc6qUR8J5biMIzKJt9hIlt984vgqnccAz7r7K3U8Ak4DBZtY5anqLGr4F+L27Lwu34YvAcWAQVf83JwmmMKiH3P1z4C5gCrDbzKYXHf7HaFvU8BaCb89F9vtXZyUcC//uiXr9GNCijHm+CMwDXg6bEh4xs1Pc/SjBzvc2YKeZ/dXMupfx/g5hLdG2EBxpFNkVNXy0nDqi16OwjOnPABoSfEsubzkVid75jgNmheEEkEkQalvCpp9BMc4zFuVts3eApmY2IDxi6wP8OXytxGfq7rnAfspf1wzgHjM7ED7+DqQDHcr5N3dm9ayaVAeFQT3l7i+7+8UE/4EBHg3/5gLNoiYt6z9sp6jhDGBHNdST7+7/4e7fImgKGg3cEL62wN1HEDTdrAfKOh11B0GzTbTOwPaq1lbKPoKjmOimrow4lrMAaGNmfYBrgelFL7h7trtfQdCMNZvg6CUWpbdZWcEevc06E26zMPD+RBBM1wGvhwFMOE3xeppZc4Kjt+hgibYVeMjdW4WPlu7ewt2LmpZK/5t7JMb1kxqgMKiHzOwcM/uOBR2/eQTf1ou+Ba8ERplZy/Bo4c4yZvFjM+sYdsDeR9AcVNWaImZ2XtgMc4Rgh1toZm3N7PKw7+BE+FphGbOYA3zTzK41s1PM7BqgJ0E7fbUJd54zgIfMrEX4bfpugiObWN6fH77/caAlQTgQtquPM7PTwia2wwRNK7FYCVxpZk3N7GxgQhnT3Gtm3zCzTgTbNHqbvURw9DWOqHAKx99kZr3NLI2g/2CxuxedULAL6Bo1/VPAP5nZBeE6NQ87xZuf5N+c1AIKg/opjeBb2V6Cb39tCNqDIdiprSbo3HuDr+/onWCHMR/4DPiUoN+gPKU7Mss7HbE98CrwBbCWoPniRYJ/oz8h+Oa9j6AT+ravLcT9AHAZQafwvvDv99397ydZbqyi3387QbPRRmAhMM3dnwMwsyFmdugk83oJGAr8qVRT1PXAJjM7CPwjwc4ZM+sUnoGTXs78fk0QlLuA54BpZUwzG8gGlhME5LPFK+a+lODo4kxgbtT4t4B/B14j+PzPIjiaKTIFeCFsErrK3bMJ+g1+a2YHgA0EHfZQ8b85qQX0ozORFGf6cZjEQEcGIiKiMBCpB3T4LyelZiIREdGRgYiIBD+eqbXMTIctIiKV4O528qm+UuuPDNy9zj0mT56c9Bq0nlpHrWf9Xc/KqPVhICIiiacwEBERhUEiRCKRZJdQI+rDetaHdQStp9TyU0vNzGtzfSIitZGZ4anWgSwiIomnMBARkdofBgUFsV7FV0REKqvWh8HTTz+d7BJERFJere9Abtu2LevWraNVq1bJLkdEpE5IyQ7kzMxMfv7znye7DBGRlFbrjwz27dvHueeey/z58+nTp0+ySxIRqfVS8sigdevWPPDAA9xxxx2VvuaGiIhUrNaHAcAtt9zCoUOHeOWVV5JdiohISqr1zURF9b3//vtcd911rFu3jhYtWiS5MhGR2islm4mKDBkyhEsuuYSHH3442aWIiKScOnNkALB9+3b69OnD4sWLOfvss5NYmYhI7ZXSRwYAHTt25N577+Xuu+9OdikiIimlToUBwF133cX69euZM2dOsksREUkZdS4M0tLSmDp1KnfddRfHjx9PdjkiIimhzoUBwKhRo+jevTtTp05NdikiIimhTnUgR/vss88YNGgQq1atomPHjjVcmYhI7VWZDuQ6GwYA9913Hzk5OUybNq0GqxIRqd3qXRgcOXKEnj178tJLLzFkyJAarExEpPZK+VNLS2vRogWPPfYYt99+u26CIyJSBXU6DACuvfZaTjvtNJ566qlklyIiUmfV6WaiIqtXr2bYsGG88cYb9O3btwYqExGpvepdM1GR3r1784c//IFRo0axdu3aZJcjIlLnNEx2AdVl7Nix5ObmMnLkSN599126deuW7JJEROqMlAkDgPHjx5Obm8uwYcNYuHAhnTp1SnZJIiJ1QkqFAcCtt97KkSNHigOhXbt2yS5JRKTWS2ifgZk9Y2a7zWx11LiWZjbfzNab2TwzO726l3vPPfdw3XXXMXz4cA4cOFDdsxcRSTmJ7kB+DhhZatxE4E137w68DUxKxIInT57M8OHD+d73vsehQ4cSsQgRkZSR8FNLzSwD+Ku79w6ffwJc4u67zaw9kOXuPcp5b0ynlpbH3bnttttYt24dc+fOpVmzZpWel4hIXVFXTi1t6+67Adx9F9A2UQsyM5588kk6depEZmamLnktIlKO2vA7g4QemjRo0IA//vGPNGnShHHjxpGfn5/IxYmI1EnJOJtot5m1i2om2lPRxFOmTCkejkQiRCKRuBfYsGFDXn75ZS6//HJuvvlm/vjHP9KgQW3IQRGRqsvKyiIrK6tK86iJPoMuBH0GvcLnjwIH3P1RM/sZ0NLdJ5bz3ir1GZR29OhRvvvd7/Ktb32LJ598ErO4mtREROqEWncJazObDkSA1sBuYDIwC5gBdAK2AFe7+8Fy3l+tYQBw6NAhhg0bxq5du+jcuTPp6emkp6fTsWPH4uH09HTOPPNMGjZMuZ9hiEg9UOvCoKoSEQYA+fn55OTksG3bthKP7du3Fw/v3buXM844g/T0dDIyMhgwYACDBw+mf//+NG3atNprEhGpLgqDapSfn8+uXbvYtm0bGzduZMmSJSxatIi1a9dy7rnnMnjw4OJHRkaGmpxEpNZIyTD48ssvady4ca3Z2R47dozs7GwWLVrE4sWLWbRoEe7OoEGDisOhb9++NG/ePNmlikg9lZJh8NBDD5Gfn0+jRo1IS0ujcePGpKWlFT+Knpc3vqzh6gwWdycnJ4dFixYVPz766CM6dOhAr169SjzOPvts9UOISMKlZBi4O4WFheTl5ZGXl8fx48eLH+U9jx5fepqTBUtFoRJrsOTn5/Ppp5+yZs2aEo+dO3fSo0ePEgHRp08f2rdvn4RPV0RSVcqGQXWKJVhiGa4oWMo7WnF3du/eTU5ODhs3bmT9+vWsXr2as846i2uvvZZrrrmG1q1bV+v6ikj9ozCoYWUFS2VC5ujRoxQUFHDixAnMjObNm9OqVSuaNm0adzNYIprCRKRuURjUYYWFhezbt4/Zs2cze/ZsNm7cyNChQ4lEInTt2rXC0Ck9nJeXV2GT18mCRMEiUrcpDFLI1q1bmT59Oi+++CJHjhxh/PjxjB8/nh49yrzAawmFhYWcOHEipv6V8vpaFCwidZfCIAW5O6tWrWLatGlMnz6d/Pz8EtdVKtq5Ru9ko8c1atSoxI679E48ely7du0YOnQoQ4YMKf5hnbvH3PQVS7AU9bHEctZXRdM1btxY15cSKYfCIMUVFBSwd+9eINhJFylvuLCwkPz8/BI744qGN2/ezJtvvsnq1asZPHgww4cPZ/jw4fTu3btadrylg6UqZ4edOHGChg0bxh0s5Z1BpmCRVKIwkGrxxRdf8M4777BgwQIWLFjAwYMHGTZsWHE4pKenl/u+nJwccnJy2LJlS4m/x44do3///gwcOJBBgwbRs2fPKu2AKxMs5U0Ta7DE0gymYJHaQGEgCbFlyxYWLFjA/Pnzeeutt2jXrh2XXnpp8Q/uinb4+fn5ZGRkkJGRQefOnYv/du7cmcaNG7Ns2TIWL17M4sWL2bNnDwMGDGDQoEEMHDiQgQMH0rZtwu5zVKGKgiWeZrCygqUqP5JUsEhlKQwk4QoKClixYgVZWVmkpaWV2Om3bNky5g7iffv2sXTp0uJwWLp0Ka1bty4+chgwYAB9+vSpc7cqjQ6WyjSDlR4f3cdSlU58BUv9ojCQOquwsJD169cXh8OyZctYt24d3bp1o3///vTr149+/frVyYCorKoGS+nh8oIlnmYwBUvdoDCQlHL8+HE++ugjli1bRnZ2NtnZ2fU+ICor1mCJdZqqXCtMwZJ4CgNJeUUBkZ2dXRwSa9eupWXLlmRkZNClS5cy/1Z0Fdm8vDz27NlT/Ni9e3fx8JdffsnYsWOJRCLaaYXKC5Z4m8GqI1hKD2sbBRQGUi8VFBSwa9cuNm/ezJYtW9iyZUvx8ObNm8nJyaF58+bFwQCU2PEfOXKENm3a0K5dO9q2bVv8aNeuHYWFhUyfPp3c3FxuvvlmbrzxRjp27JjkNU4dRcFSHRehjO68r+qPJOt6sCgMRMrg7uzZs6c4HMysxA7/G9/4RoX/8d2dZcuW8fTTTzNjxgwGDx7MhAkTGD16NI0aNYq7nsOHD7Ny5Uo2bNjAhRdeSI8ePfTL7GpQE8ESbzNYsoJFYSCSYLm5ubz66qs888wzbNiwgeuvv54JEyaUe5mQL774ghUrVpCdnc3y5cvJzs5m69at9OrVi27durFw4UKaNm3KmDFjuOKKKxg0aBCnnHJKDa+VlBZLsMRzOnL0lQAq2wwWT7AoDERq0IYNG3j22Wd5/vnn6dq1KxMmTCAjI6PEjn/nzp307t27uLO7b9++9OzZs/iIwt1Zvnx58QUKd+3axejRoxkzZgzDhg3T/bZTgLtX6lph5Q3HEiwjRoxQGIjUtBMnTjBnzhyeffZZ9u/fX7zj79evHz169Ijrm/6mTZuYPXs2s2bNYsWKFQwdOpQxY8Zw2WWXJfVeF8ePH+fNN99k5syZzJkzh86dOzNixAhGjhzJoEGDKtVcJvGLNVguvvhihYFIqti/fz9/+9vfmDVrFm+99RZ9+/YlMzOTsWPH1kgndm5uLnPnzuW1115j7ty59OrViyuvvJLRo0ezdetW5s2bx/z58/n888+JRCKMHDmSkSNH0rVr14TXVpfk5+cX39Cq6LF3715uvfVWLrjggoQsU81EIinq2LFjzJ8/n5kzZ/L666/To0cPMjMzyczMpEuXLtW2nIMHD/L6668zc+ZM3nrrLQYNGkRmZiZjxowp9/ase/bsKb5cyfz582nevDkjR45kxIgRXHrppZx66qlfe8+JEyfIzc3l6NGjX/vbsWNHzjnnnDrTd+LuHDhwgC1btrBx40Y+//zzEjv+bdu20b59e7p27Vr8aNy4MVOnTuXiiy/ml7/8ZbUHqMJApB7Iy8vj7bffZubMmcyePZtOnToVB0P37t1jnk9BQQF79uxh586dLF++nJkzZ/LBBx8QiUTIzMxk9OjRtGrVKq7a3J01a9Ywb9485s2bx5IlS8jIyCAvL4/c3NziHX5hYSHNmzenWbNmNG/evPjRpEkTtm7dyo4dOzjvvPM4//zzix+9evWiRYsW8X5cVZKXl8eOHTvYvn17uY8dO3bQpEkTOnfuTLdu3Urs9Lt160bnzp1JS0v72rxzc3P51a9+xdSpU7nhhhu4//77q60pUGEgUs/k5+fz3nvvMXPmTF577TVatWpFZmYml19+OQA7duxg586dZf7du3cvLVu2pEOHDvTo0YOxY8cyatSoMr/JV9bRo0fZsGEDTZs2Ld7hN2vW7KQ3Ojp06BBr1qxh5cqVxY+1a9eSnp5eIiD69+9f7Rc4PHbsGE888QRPPPEEe/fupX379nTs2LHCR1V+Ab97924eeOABZsyYwU9/+lNuv/12mjRpUqV1UBiI1GOFhYUsXryYmTNnMnfuXNLS0ujQoQNnnnlm8d/o4Xbt2tG4ceNklx2z/Px81q9fXyIgsrOzGT16NBMnTqRnz55Vmn9BQQEvvPACP//5zxkwYAAPPPAA5557bo01V33yySdMnDiRFStW8NBDDzFu3LhK/z5BYSAi9crBgwf53e9+xxNPPMFFF13EpEmTGDBgQFzzcHfmzJnDxIkTOf3003nssce48MILE1TxyS1cuJB7772XEydO8PjjjzN06NC456EwEJF6KTc3l2eeeYbHH3+cnj17MmnSJCKRyEl/2b106VJ++tOfsnv3bh555BEuv/zyWvFrcHdnxowZTJo0ie7du3PjjTdyyimnFN/J0N1LDEf/bdSoEVdddZXCQETqr7y8PKZNm8ajjz5Kq1atmDRpEpdddtnXmls+++wz7rvvPj744AOmTJnCTTfdRMOGDZNUdfmOHz/O73//e959913M7Gv3PC/r74EDB1iwYEHdCQMzuxuYABQCa4Cb3D2v1DQKAxGJW0FBAa+99hoPP/wweXl5TJo0iWuuuYYDBw7w4IMP8vLLL3P33Xdz1113VXhF27qqzjQTmVkH4H2gh7vnmdkrwN/c/YVS0ykMRKTS3J158+bx8MMPk5OTw+HDh/nhD3/I/fffT5s2bZJdXsLUtTBYBJwPHAb+DPzG3d8sNZ3CQESqxYcffsgZZ5zBWWedlexSEq7OhAGAmd0BPAQcBea7+/VlTKMwEBGJU2XCICk9Jmb2DWAMkAF8AbxqZuPcfXrpaadMmVI8HIlEiEQiNVSliEjdkJWVRVZWVpXmkaxmoquAke5+S/j8emCgu/9Lqel0ZCAiEqfKHBkk675uOcAgM2tiwTlRQ4F1SapFRKTeS0oYuPtS4FVgBbAKMOC/k1GLiIjoR2ciIimnLjUTiYhILaIwEBERhYGIiCgMREQEhYGIiBBDGJjZRWbWPBweb2a/MrOMxJcmIiI1JZYjg/8CjppZH+Ae4HPghYrfIiIidUksYZAfnuw/Bvitu/8OqL47ZouISNLFcqG6w2Y2CbgeuNjMGgCNEluWiIjUpFiODK4BjgM3u/suIB14PKFViYhIjYrpchRhh/E33f1NM2sGnOLuhxNenC5HISISt4RcjsLMbiG4qNwfwlEdgVnxlyciIrVVLM1EPwYuAg4BuPunQNtEFiUiIjUrljA47u55RU/MrCGgthsRkRQSSxi8a2b3AU3NbDgwA/hrYssSEZGadNIO5PBU0gnACIKb0Mxz96dqoDZ1IIuIVEJlOpBjCYM73f03JxuXCAoDEZH4JermNj8qY9yN8SxERERqt3J/gWxm1wHjgLPM7C9RL50GHEh0YSIiUnMquhzF/wI7gTOA/4wafxhYnciiRESkZsX6C+R2wIDw6VJ335PQqr5arvoMRETilKhfIP8AWAr8ALgaWGJmV1WuRBERqY1iOZtoFTC86GjAzNoAb7p7n4QXpyMDEZG4JepsogalmoX2x/g+ERGpI2K5n8EbZjYPeCl8fg0wJ3EliYhITYu1A/lKYEj49D13/3NCq/pquWomEhGJU6J+gfwT4BV3316V4ipDYSAiEr9E9RmcCsw3s/fM7F/C00xFRCSFxNRMBGBmvQn6CzKBbe4+rEoLNjsdeBo4DygkuK3mklLT6MhARCROlTkyiKUDucgeYBfB2UTVcXOb3wBz3P0H4T0SmlXDPEVEpBJi6TP4Z4Ifm7UhuJfBn9z94yot1Ow0YIW7dzvJdDoyEBGJU6KODDoBd7n7ysqVVaazgH1m9hzQB1gG3Onux6pxGSIiEqOThoG7TzKzvmZ2B8HtLj9w9+XVsNy+wI/dfZmZTQUmApNLTzhlypTi4UgkQiQSqeKiRURSS1ZWFllZWVWaRyzNRP9O0Ez0WjjqCmCGu/+i0gsNzkha5O5dw+dDgJ+5++hS06mZSEQkTon6ncF6oI+7fxk+bwqsdPfula40mM+7wC3uvsHMJgPN3P1npaZRGIiIxClRfQY7gCbAl+HzNKA6foB2B/A/ZtYI2AjcVA3zFBGRSojlyGAWwb0MFhD0GQwnuKT1NgB3vyNhxenIQEQkbolqJirrHsjF3P35eBYYD4WBiEj8EhIGyaQwEBGJX0L6DMxsE0HzUAlFZwKJiEjdF0sHcv+o4SYEt79slZhyREQkGSrVTGRm2e7eLwH1lF6OmolEROKUqGaivlFPGxAcKcRzgTsREanlYtmp/2fUcD6wieAXySIikiJ0NpGISIpJ1J3OREQkxSkMREREYSAiIrGdTXRlGaO/ANa4+57qL0lERGpaLNcm+hswGHgnHBUBsgnuVvagu7+YsOLUgSwiErdEXcK6IdDT3XeHC2kHvAAMBBYCCQsDERGpGbH0GXQqCoLQnnDcAeBEYsoSEZGaFMuRQZaZvQ7MCJ9fFY5rDhxMWGUiIlJjYukzMOBKYEg46gNgZk005qvPQEQkfgnpM3B3N7P3gTyCS1kv1R5aRCS1nLTPwMyuJrjN5VUE1yRaYmZXJbowERGpObE0E60Chhf9psDM2gBvunufhBenZiIRkbgl6tpEDUr9uGx/jO8TEZE6Ipazid4ws3nAS+Hza4A5iStJRERqWkyXsDazTOCi8Ol77v7nhFb11XLVTCQiEqfKNBPpfgYiIimmWk8tNbPDBKeSfu0lgjNOT4uzPhERqaXKDQN3P7UmCxERkeTRWUEiIqIwEBGRJIeBmTUws+Vm9pdk1iEiUt8l+8jgTuDjJNcgIlLvJS0MzCwdGAU8nawaREQkkMwjg18D91L26asiIlKDYrkcRbUzs+8Du919pZlFCH67UKYpU6YUD0ciESKRSKLLExGpU7KyssjKyqrSPJLyC2Qz+yUwHsgHmgKnAq+5+w2lptMvkEVE4lQnL0dhZpcA97j75WW8pjAQEYlToi5hLSIiKS7pRwYV0ZGBiEj8dGQgIiKVojAQERGFgYiIKAxERASFgYiIoDAQEREUBiIigsJARERQGIiICAoDERFBYSAiIigMREQEhYGIiKAwEBERFAYiIoLCQEREUBiIiAgKAxERQWEgIiIoDEREBIWBiIigMBARERQGIiKCwkBERFAYiIgICgMREUFhICIiKAxERIQkhYGZpZvZ22a21szWmNkdyahDREQC5u41v1Cz9kB7d19pZi2AbGCMu39SajpPRn0iInWZmeHuFs97knJk4O673H1lOHwEWAd0TEYtIiJSC/oMzKwLcD6wJLmViIjUX0kNg7CJ6FXgzvAIQUREkqBhshZsZg0JguBFd59d3nRTpkwpHo5EIkQikYTXJiJSl2RlZZGVlVWleSSlAxnAzF4A9rn7TyqYRh3IIiJxqkwHcrLOJroIWAisATx83Ofub5SaTmEgIhKnOhMGsVIYiIjEr86cWioiIrWLwkBERBQGIiKiMBARERQGIiKCwkBERFAYiIgICgMREUFhICIiKAxERASFgYiIoDAQEREUBiIigsJARERQGIiICAoDERFBYSAiIigMREQEhYGIiKAwEBERFAYiIoLCQEREUBiIiAgKAxERQWEgIiIoDEREBIWBiIigMBARERQGIiKCwiAhsrKykl1CjagP61kf1hG0nqIwSIj68g+uPqxnfVhH0HqKwkBERFAYiIgIYO6e7BrKZWa1tzgRkVrM3S2e6Wt1GIiISM1QM5GIiCgMREREYVDtzGyzma0ysxVmtjTZ9VQXM3vGzHab2eqocS3NbL6ZrTezeWZ2ejJrrKpy1nGymW0zs+Xh47vJrLE6mFm6mb1tZmvNbI2Z3RGOT5ntWcY63h6OT6ntaWZpZrYk3N+sMbPJ4fguZrbYzDaY2Utm1vCk81KfQfUys41AP3f/e7JrqU5mNgQ4Arzg7r3DcY8C+939MTP7GdDS3Scms86qKGcdJwOH3f1XSS2uGplZe6C9u680sxZANjAGuIkU2Z4VrOM1pN72bObuR83sFOAD4E7gJ8Cr7j7DzP4LWOnuf6hoPjoyqH5GCn6u7v4+UDrgxgDPh8PPA1fUaFHVrJx1hGCbpgx33+XuK8PhI8A6IJ0U2p7lrGPH8OVU255Hw8E0oCHgwHeAmeH454GxJ5tPyu20agEH5pnZh2Z2S7KLSbC27r4bgv98QNsk15MoPzazlWb2dF1uOimLmXUBzgcWA+1ScXtGreOScFRKbU8za2BmK4BdwALgc+CguxeGk2wDOpzkPfo9AAAFOUlEQVRsPgqD6neRu/cHRhH8oxuS7IJqUCq2OT4JdHP38wn+s6VS80IL4FXgzvDbc+ntV+e3ZxnrmHLb090L3f3bBEd3FwA9KjMfhUE1c/ed4d+9wJ8JNk6q2m1m7aC4jXZPkuupdu6+17/qWHsKGJDMeqpL2KH4KvCiu88OR6fU9ixrHVN1ewK4+yEgCxgMfMPMivbv6cD2k71fYVCNzKxZ+E0EM2sOjAA+Sm5V1coo2d76F+DGcPhHwOzSb6iDSqxjuFMsciWpsz2fBT52999EjUu17fm1dUy17WlmZxQ1dZlZU2A48DHwDvCDcLKYtqXOJqpGZnYWwdGAE3Tk/I+7P5LcqqqHmU0HIkBrYDcwGZgFzAA6AVuAq939YLJqrKpy1vE7BO3NhcBm4NaidvW6yswuAhYCawj+rTpwH7AU+BMpsD0rWMdxpND2NLNeBB3EDcLHK+7+ULgvehloCawAxrv7iQrnpTAQERE1E4mIiMJAREQUBiIigsJARERQGIiICAoDERFBYSB1QHjZ4Z9U4f23mtn4ON9zppn9qbLLrAlV/VxEop30Gtcidd3JLt1bznt2AlcnoJxaw8zM9UMjCenIQGpceNmO18Mbcqw2sx+E4zeZWatwuJ+ZvRP1tvPN7H/DG6/8QzjNJWaWZWazzOwzM3vYzMaFN/tYFf4Ks8Q3aDO7I7zhycrwF8dF81kR3uwk28yam1mGma0JX08zs2fDWrPNLBKO/5GZzTSzuWFdj5azvpvMbEr43lVmdk7pusLna8ysc7jsdWb2XDjfaWY21MzeD5/3r+hzCef1r2a2NFzPohueZJjZJ2b2fLhu6VXZjpJadGQgyfBdYLu7XwZgZqeG4yu6amYvYCBwKrDCzF4Px/cmuErjQWAj8JS7D7Tg7l23E9zkI9rPgC7ufsLMTgvH3QP8s7svMrNmwJellv9joNDde5tZd2C+mX0zfK0PweUNTgDrzewJdy/romB73L2fmd0G/Cvwj2VME72+3YBMd//YzJYB17n7EDO7HPg3vro+fVmfSy/gm+5+gZkZ8Jfw6rlbgbOB6939wzKWL/WYjgwkGdYAw8Nv8kPc/XA4vqKbjsx29zx33w+8zVdXg/3Q3fe4ex7BddznRy2jSxnzWQVMN7MfAgXhuA+AX1twa8SWUdeBLzIEmAbg7usJrmlzTvjaW+5+xN2PE1wgLKOc+v8c/s0upy4ouf6b3P3jcHgt8FbUekUvo6zPZQTB57scWA50B4rCa4uCQMqiMJAa5+6fAn0Jdmy/MLP7w5fy+erfZJPSb4satqjnx6PGF0Y9L6TsI9/vA78Nl/+hmTVw90eBCUBT4IOiZpwKRO+0o5dfUM4yo6eLniZ6faHkOse6XuV9Lg+7e193/7a7n+Puz4Xjc8upT+o5hYHUODM7Ezjm7tOBxwl2zACbgH7hcGapt40xs8Zm1hq4BKjst9vO7v4uMBE4DWhhZl3dfa27PxbOt/TNQd4DfhjWfg7BVT3XV3L50TYTrruZ9QXOinot1lszlvW5zAduDi+jjpl1MLM2cc5X6hn1GUgy9AIeN7NCIA+4LRz/IPCMmX1BcJOOaKvDca2BB919V9h+H63CM2MsuNnJtLCvwIDfuPshM/uFmX2H4Fv7WmAuJW8T+CTwX2a2mqBv4Edhn0PpRZS3/PLGzwRuCDtzl1AyYLyc4dK+9rkAu8ysB7AorPEwMJ7gqEJnD0mZdAlrERFRM5GIiCgMREQEhYGIiKAwEBERFAYiIoLCQEREUBiIiAgKAxERAf4P6tkVhRZtU48AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fed852b5b00>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(x,ylog, color='k')\n",
    "plt.xlim(1,30)\n",
    "x1 = [0,30]\n",
    "y1 = [thetalog[0] + thetalog[1]*a for a in x1]\n",
    "plt.plot(x1,y1, color='grey')\n",
    "plt.xlabel(\"submission number\")\n",
    "plt.ylabel(\"log upvotes\")\n",
    "plt.yticks([2,4,6,8,10])\n",
    "plt.title(\"submission no. vs. upvotes\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Simple regression question, same form as the examples above"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "path = dataDir + \"fantasy_100.json\"\n",
    "f = open(path)\n",
    "data = []\n",
    "for l in f:\n",
    "    d = json.loads(l)\n",
    "    data.append(d)\n",
    "f.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "ratings = [d['rating'] for d in data]\n",
    "lengths = [len(d['review_text']) for d in data]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = numpy.matrix([[1,l] for l in lengths]) # Note the inclusion of the constant term\n",
    "y = numpy.matrix(ratings).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LinearRegression(copy_X=True, fit_intercept=False, n_jobs=None, normalize=False)"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model = sklearn.linear_model.LinearRegression(fit_intercept=False)\n",
    "model.fit(X, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Extract theta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[3.98394783e+00, 1.19363599e-04]])"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "theta = model.coef_\n",
    "theta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_pred = model.predict(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [],
   "source": [
    "sse = sum([r**2 for r in y - y_pred])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "MSE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[1.26575168]])"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mse = sse / len(y)\n",
    "mse"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fit a model with an additional variable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [],
   "source": [
    "def feat(d):\n",
    "    f = [1]\n",
    "    f.append(len(d['review_text']))\n",
    "    f.append(d['n_comments'])\n",
    "    return f"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = [feat(d) for d in data]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LinearRegression(copy_X=True, fit_intercept=False, n_jobs=None, normalize=False)"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model = sklearn.linear_model.LinearRegression(fit_intercept=False)\n",
    "model.fit(X, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Extract theta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[3.95432090e+00, 7.24295694e-05, 1.08071028e-01]])"
      ]
     },
     "execution_count": 61,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "theta = model.coef_\n",
    "theta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the MSE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[1.22988869]])"
      ]
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y_pred = model.predict(X)\n",
    "sse = sum([r**2 for r in y - y_pred])\n",
    "mse = sse / len(y)\n",
    "mse"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(explanation: the coefficient of length becomes smaller, as the variability in ratings is already largely explained by the number of comments)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sketch proof:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\frac{\\partial}{\\partial \\theta_0} \\sum_i (\\theta_0 - y_i)^2 = 2\\sum_i (\\theta - y_i)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "equals zero when $\\sum_i \\theta = \\sum_i y_i$, i.e., $\\theta = \\frac{1}{N} \\sum_i y_i$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sketch:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\frac{\\partial}{\\partial \\theta_0} \\sum_i | \\theta_0 - y_i | = \\sum_i \\delta(y_i > \\theta_0) - \\sum_i \\delta(y_i < \\theta_0)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "zero when $\\theta_0$ is the median value of $y$ (meaning that the two terms balance out)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sketch:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\max_\\theta \\prod_i \\frac{1}{2b}\\exp(-\\frac{|x_i\\cdot \\theta - y_i|}{b})$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$= \\max_\\theta \\sum_i -| x_i \\cdot \\theta - y_i |$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$= \\min_\\theta \\sum_i | x_i \\cdot \\theta - y_i|$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.6"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sketch:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "rewrite $\\sum_i (x_i \\cdot \\theta - y_i)^2$ as $(X\\theta - y)^T(X\\theta - y)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$= (\\theta^TX^T - y^T)(X\\theta - y)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$= \\theta^TX^TX\\theta -2y^TX\\theta -y^Ty$ (all terms are scalar)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\frac{\\partial}{\\partial \\theta} (\\theta^TX^TX\\theta -2y^TX\\theta +y^Ty) = 2(\\theta^T X^TX - y^TX)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$= 0$ when $\\theta^TX^TX = y^TX$; or when $X^T X \\theta = X^Ty$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "i.e., when $\\theta = (X^TX)^{-1}X^Ty$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.7"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similar to 2.3: when solving $\\theta$ by computing $\\frac{\\partial}{\\partial \\theta} \\sum_i (x_i \\cdot \\theta - y_i)^2 = 2\\sum_i (x_i \\cdot \\theta - y_i)$, the expression will be minimized when $\\sum_i (x_i \\cdot \\theta - y_i) = 0$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}