{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Homework 5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "from IPython.display import display, Latex, Markdown\n",
    "\n",
    "!pip install -U okpy\n",
    "from client.api.notebook import Notebook\n",
    "ok = Notebook('hw5.ok')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## This assignment\n",
    "\n",
    "In this assignment we will study a fundamental problem in crawling the web: *deciding when to re-index changing pages*.  \n",
    "\n",
    "Search engines like those deployed at Google and Microsoft periodically visit, or *crawl*, web pages to check their current status and record any changes to the page.  Some pages, such as news websites, are updated frequently while other pages, [like this one](http://youvegotmail.warnerbros.com/cmp/0frameset.html), are never updated. \n",
    "\n",
    "**Why not crawl the entire web all the time?**\n",
    "\n",
    "The web is **big** and even giant search engines don't have the resources to re-read the entire web continuously.  Fortunately, different sites change at different rates, so an intuitive strategy would be to visit sites more if they are changed more frequently.  But this requires us to know, or at least estimate, how often each web site is changed.  This is the main question you'll answer in this assignment.\n",
    "\n",
    "\n",
    "Along the way you will:\n",
    "\n",
    "* do some EDA to understand the site dataset.\n",
    "* build several *models* for the site changes.\n",
    "* use simulation to get a concrete understanding of your models.\n",
    "* compute maximum-likelihood estimates of model parameters.\n",
    "* use simulation and visualization to understand the accuracy of your estimates of site change frequencies."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Examining the data\n",
    "Our dataset consists of observations for about 1000 pages over a 30-day period.  Every hour, each page was downloaded, and it was determined whether the page had been changed since the previous hour.  For the first hour, we can't tell whether it has changed, so there are 719 such \"checks\" for each page.\n",
    "\n",
    "**Important fact about the data:** Not every check succeeded.  Sometimes the page failed to load, or it was otherwise impossible to tell whether it had changed since the previous hour.  These hours are **omitted** from the dataset.  A field in the dataset indicates the number of successful checks for each page.\n",
    "\n",
    "The dataset is in JSON format in the file `crawl.json`.  For each page, we have: the URL, called `url`; the number of successful visits to the page, `number of checks`; and the checks when a change was detected, called `positive checks`.  Examine the `crawl.json` file.  You might find it convenient to load it into python using `json.load`.  The next cell two cells are provided for that."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use this cell to examine the dataset, if you like.\n",
    "!head -n 25 crawl.json"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# This cell loads the JSON data into Python.\n",
    "import json\n",
    "with open(\"crawl.json\", \"r\") as f:\n",
    "    crawl_json = json.load(f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "---\n",
    "\n",
    "## Question 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fill in the function `json_description` to determine:\n",
    "\n",
    "1. the number of records, and\n",
    "2. the set of possible top level *keys* (field names) in each record."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "for_assignment_type": "solution"
   },
   "outputs": [],
   "source": [
    "def json_description(crawl_json_records):\n",
    "    \"\"\"Produces information about a JSON object containing crawl data.\n",
    "    \n",
    "    Args:\n",
    "      crawl_json_records (list): A list of JSON objects such as the\n",
    "                                 crawl_json variable above.\n",
    "    \n",
    "    Returns:\n",
    "      2-tuple: An (int, set) pair.  The integer is the number of records\n",
    "               in the given list.  The set is a set (constructed with\n",
    "               Python's set() function) of strings.  It contains all the\n",
    "               field top names that appear at the top level of any record\n",
    "               in crawl_json_records.\n",
    "    \"\"\"\n",
    "    n = len(crawl_json_records)\n",
    "    keys = set(k for f in crawl_json for k in f.keys())\n",
    "    return (n, keys)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = ok.grade('q01')\n",
    "_ = ok.backup()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Display results\n",
    "n, keys = json_description(crawl_json)\n",
    "print(\"Number of records:\", n)\n",
    "print(\"Keys in the records:\", keys)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## Question 2\n",
    "What is the granularity of the dataset as represented in `crawl.json`?  Write a one to two sentence description in the cell below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "for_assignment_type": "solution"
   },
   "outputs": [],
   "source": [
    "# Use this cell for your explorations.\n",
    "q2_answer = r\"\"\"\n",
    "\n",
    "**SOLUTION:** Each top-level JSON object represents all the crawls for one site.\n",
    "\n",
    "\"\"\"\n",
    "\n",
    "display(Markdown(q2_answer))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "---\n",
    "\n",
    "## Question 3\n",
    "\n",
    "It will be more convenient to work with the data in Pandas DataFrame with a rectangular format.  Fill in the function `make_crawls_dataframe` in the cell below.  **Then** run the cell below that to create the table `crawls`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "for_assignment_type": "solution"
   },
   "outputs": [],
   "source": [
    "def make_crawls_dataframe(crawl_json):\n",
    "    \"\"\"Creates a Pandas DataFrame from the given list of JSON records.\n",
    "    \n",
    "    The DataFrame corresponds to the following relation:\n",
    "\n",
    "        crawls(primary key (url, hour), updated)\n",
    "\n",
    "    `updated` is a boolean value indicating whether the check for that\n",
    "    hour found a change.\n",
    "\n",
    "    The result is sorted by URL in ascending order and **further**\n",
    "    sorted by hour in ascending order among the rows for each URL.\n",
    "    \n",
    "    Args:\n",
    "      crawl_json_records (list): A list of JSON objects such as the\n",
    "                                 crawl_json variable above.\n",
    "    \n",
    "    Returns:\n",
    "      DataFrame: A table whose schema (and sort order) is described\n",
    "                 above.\n",
    "    \"\"\"\n",
    "    pages = pd.DataFrame.from_records(crawl_json)\n",
    "    pages = pages[['url', 'number of checks']]\n",
    "    \n",
    "    positive_crawls = [{'url': page['url'], 'hour': hour-1, 'updated': True} \n",
    "                   for page in crawl_json \n",
    "                   for hour in page['positive checks']]\n",
    "\n",
    "    negative_crawls = [{'url': page['url'], 'hour': hour-1, 'updated': False} \n",
    "                   for page in crawl_json \n",
    "                   for hour in (set(range(1, page['number of checks']+1)) - \n",
    "                               set(page['positive checks']))]\n",
    "\n",
    "    crawls = pd.DataFrame.from_records(positive_crawls + negative_crawls)\n",
    "    crawls.set_index(['url', 'hour'], inplace=True)\n",
    "    crawls.sort_index(inplace=True)\n",
    "    return crawls"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run this cell before you continue.\n",
    "crawls = make_crawls_dataframe(crawl_json)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = ok.grade('q03')\n",
    "_ = ok.backup()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 4\n",
    "There are other reasonable ways to represent the data in a relational database (or in DataFrames).  One alternate schema uses 2 relations to represent the data.  The first relation is:\n",
    "\n",
    "    crawls_that_found_changes(primary key (url, hour))\n",
    "\n",
    "The primary key for the second relation is `url`.  Define a schema for that 2nd relation.  Your schema should ensure that the two resulting tables could be used to answer any question that could be answered using `crawls`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "for_assignment_type": "solution"
   },
   "outputs": [],
   "source": [
    "# Use this cell for your explorations.\n",
    "q4_answer = r\"\"\"\n",
    "\n",
    "**SOLUTION:** `pages(primary key (url), num_crawls)`\n",
    "\n",
    "\"\"\"\n",
    "\n",
    "display(Markdown(q4_answer))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "---\n",
    "\n",
    "## Question 5\n",
    "\n",
    "In the following we will construct visualizations to understand several key quantities of interest.  The following cell constructs some key summary statistics that we will be using in the next few questions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "crawl_stats = (\n",
    "    crawls['updated']\n",
    "        .groupby(crawls.index.get_level_values('url'))\n",
    "        .agg({\n",
    "            'number of crawls': 'count', \n",
    "            'proportion of updates': 'mean', \n",
    "            'number of updates': 'sum'\n",
    "        })\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part 1:\n",
    "\n",
    "What was the distribution of the number of crawls for each page?  Did most get crawled all 719 times?  (For this and the other parts of this question, create a *visualization* to answer the question.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "for_assignment_type": "solution"
   },
   "outputs": [],
   "source": [
    "crawl_stats['number of crawls'].hist(bins=20)\n",
    "plt.xlabel(\"Number of Crawls\")\n",
    "plt.ylabel(\"Number of Pages\")\n",
    "\n",
    "# Leave this for grading purposes\n",
    "q5_p1_plot = plt.gcf()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part 2\n",
    "\n",
    "What was the distribution of the number of positive checks for each page?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "for_assignment_type": "solution"
   },
   "outputs": [],
   "source": [
    "crawl_stats['number of updates'].hist(bins=20)\n",
    "plt.xlabel(\"Number of Updates\")\n",
    "plt.ylabel(\"Number of Pages\")\n",
    "\n",
    "# Leave this for grading purposes\n",
    "q5_p2_plot = plt.gcf()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part 3\n",
    "\n",
    "What is the relationship between the number of crawls for each page and the number of positive checks? Construct a scatter plot relating these two quantities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "for_assignment_type": "solution"
   },
   "outputs": [],
   "source": [
    "crawl_stats.plot.scatter('number of crawls', 'number of updates')\n",
    "\n",
    "# Leave this for grading purposes\n",
    "q5_p3_plot = plt.gcf()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part 4\n",
    "\n",
    "In 2 or 3 sentences, describe what you discovered in your initial explorations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "for_assignment_type": "solution"
   },
   "outputs": [],
   "source": [
    "q5_p4_answer = r\"\"\"\n",
    "\n",
    "**SOLUTION:** Most sites were crawled every time, but a substantial number weren't.  \n",
    "Most sites had only a few positive checks, but some sites were changed very frequently.  \n",
    "It appears that the number of positive checks increases with the number of crawls.  \n",
    "Of course, the number of positive checks is always less than the number of crawls!\n",
    "\n",
    "\"\"\"\n",
    "\n",
    "display(Markdown(q5_p4_answer))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Making a timeline from one site"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It will be useful to be able to look at timelines of positive checks or changes for sites.  The function `display_points`, defined below, will help."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def display_points(points, xlim, title):\n",
    "    \"\"\"Displays a timeline with points on it.\n",
    "    \n",
    "    Args:\n",
    "      points (ndarray): A list of floats in the range [xlim[0], xlim[1]],\n",
    "                        each one a point to display in the timeline.\n",
    "      xlim (list-like): A list/tuple/array with 2 elements giving the\n",
    "                        start and end of the timeline.\n",
    "      title (str): The title to display on the plot.\n",
    "      \n",
    "    Example:\n",
    "        >>> # plot the points at [1,3,30,50]\n",
    "        >>> display_points([1,4,30,50], [0, 75], \"Example\")\n",
    "    \"\"\"\n",
    "    fig, ax = plt.subplots(1)\n",
    "    fig.set_size_inches(8, 1)\n",
    "    ax.scatter(points, np.repeat(0, len(points)), alpha=0.5)\n",
    "    ax.axhline(0, color=\"grey\", zorder=-1, lw=.5)\n",
    "    ax.yaxis.set_visible(False)\n",
    "    ax.xaxis.set_visible(True)\n",
    "    ax.set_xlim(xlim)\n",
    "    ax.set_title(\"{} ({:d} total points)\".format(title, len(points)))\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Example Usage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "display_points([1,4,30,50], [0, 75], \"Example\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "---\n",
    "\n",
    "## Question 6\n",
    "\n",
    "We want to understand the behavior of page changes in order to determine how often a page should be visited. To do this, we need more than summary statistics of the 1000 sites.  We also need to examine the patterns of positive checks for sites. Let's examine when the positive checks occurred for a handful of sites.  The `display_points` function is helpful here.\n",
    "\n",
    "We selected a variety of sites - some changed frequently, others rarely, and some had only a few successful crawls. Their URLs are available in `pages_to_display`.  Use `display_points` to examine each page.  You should notice that the visualization is not very informative in the last 3 cases."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This cell identifies a few categories of pages and \n",
    "# associates a label in the 'Desc' column of crawl_stats\n",
    "\n",
    "crawl_stats['Desc'] = \"\" # default blank description\n",
    "\n",
    "# Normal pages have a moderate number of updates and \n",
    "# were successfully crawled most times.\n",
    "crawl_stats.loc[\n",
    "        (crawl_stats['proportion of updates'] > 0.1)\n",
    "        & (crawl_stats['proportion of updates'] < 0.9)\n",
    "        & (crawl_stats['number of crawls'] >= 700), \n",
    "    'Desc'] = 'Normal'\n",
    "\n",
    "crawl_stats.loc[\n",
    "        (crawl_stats['proportion of updates'] < .1)\n",
    "        & (crawl_stats['number of crawls'] >= 700), \n",
    "    'Desc'] = 'Rare Update'\n",
    "\n",
    "crawl_stats.loc[\n",
    "        (crawl_stats['proportion of updates'] > .9)\n",
    "        & (crawl_stats['number of crawls'] >= 700), \n",
    "    'Desc'] = 'Frequent Update'\n",
    "\n",
    "crawl_stats.loc[\n",
    "        crawl_stats['number of crawls'] < 50, \n",
    "    'Desc'] = 'Few Crawls'\n",
    "\n",
    "\n",
    "# Build a dataframe with a few examples from each type of webpage\n",
    "num_of_each = 3\n",
    "pages_to_display = pd.concat([\n",
    "    crawl_stats[crawl_stats['Desc'] == \"Normal\"].head(num_of_each),\n",
    "    crawl_stats[crawl_stats['Desc'] == 'Rare Update'].head(num_of_each),\n",
    "    crawl_stats[crawl_stats['Desc'] == 'Frequent Update'].head(num_of_each),\n",
    "    crawl_stats[crawl_stats['Desc'] == 'Few Crawls'].head(num_of_each),\n",
    "])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Construct the timeline point diagrams for each example\n",
    "for (url, desc) in pages_to_display['Desc'].iteritems():\n",
    "    crawls_for_site = crawls.loc[url]\n",
    "    display_points(\n",
    "        crawls_for_site[crawls_for_site['updated']].index,\n",
    "        [0, len(crawls_for_site)], desc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Examining the Times of Positive Checks\n",
    "\n",
    "Consider what you learned about the positive checks in the previous question. Do these occur at regular intervals? Or do they appear to occur more *randomly* in time?  Write a 2 to 3 sentence explanation of what you see:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "for_assignment_type": "solution"
   },
   "outputs": [],
   "source": [
    "q6_answer = r\"\"\"\n",
    "\n",
    "**SOLUTION:** They do seem to occur randomly in time.  However, it's hard to tell even with the\n",
    "visualization - our brains are often not very good at distinguishing random noise\n",
    "from patterns!\n",
    "\n",
    "\"\"\"\n",
    "\n",
    "display(Markdown(q6_answer))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Understanding this distribution of positive checks will help us determine the distribution of changes. However, determining this distribution is difficult to do if we do not crawl the page sufficiently often. (We will soon see that it is also hard to do if the site changed too many times.) For these reasons, we will focus on \"normal\" pages, i.e., pages that were successfully crawled at least 700 times and that had between 50 and 200 updates.\n",
    "\n",
    "Do the positive checks occur uniformly at random? Let's compare the observed distribution of the positive checks to the uniform distribution for a few of our pages. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "---\n",
    "\n",
    "## Question 7\n",
    "\n",
    "To compare the data distribution to a hypothetical probability distribution, we can compare the histogram of the data to the probability density function.  However, a more effective comparison is to compare the quantiles of the observed distribution to those of the  probability distribution. \n",
    "\n",
    "The $q^{th}$ quantile of the data, for $0 < q < 1$ is that point $x_q$ such that at least $q$ of the observations are at or below $x_q$ and at least $1-q$ of the observations are at or above $x_q$. For example, the median is the 0.5 quantile and the lower and upper quartiles are the 0.25 and 0.75 quantiles.\n",
    "\n",
    "The $q^{th}$ quantile of a continuous probability distribution, for $0 < q < 1$ is that point $x_q$ such that \n",
    "$P(X \\leq x_q) = q$ and $P(X \\geq x_q) = 1-q$. The quantiles of a uniform distribution are easy to compute. For example, for a Uniform$(0, N)$ distribution, we have $x_q = qN$.\n",
    "\n",
    "If the data come from the same probability distribution, then their quantiles should roughly match.  If we make a scatter plot of pairs of the empirical and theoretical quantiles, like these 3 points:\n",
    "\n",
    "    (lower quartile of the data, lower quartile of the pdf),\n",
    "    (median of data, median of the pdf),\n",
    "    (upper quartile of the data, upper quartile of the pdf)\n",
    "\n",
    "...the points should fall roughly on a line.\n",
    "\n",
    "Such a plot is called a quantile-quantile or Q-Q plot.  For the websites with `Desc = 'Normal'` websites in `pages_to_display`, make a Q-Q plot for that site's updates, as follows:\n",
    "\n",
    "* Let $N$ be the number of crawls for this page and $n$ be the number of positive checks.\n",
    "* Order the positive check values for a site from smallest to largest. We use these as the quantiles. That is, the $k^{th}$ smallest observation is the $k/(n+1)$-quantile.\n",
    "* Compute the $n$ corresponding quantiles for the uniform distribution on $(0, N)$. These are simply\n",
    "$$k \\times \\frac{N}{n+1},$$\n",
    "for $k = 1, \\ldots, n$.  (Why?)\n",
    "* Plot the $n$ pairs of quantiles in a scatter plot."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "for_assignment_type": "solution"
   },
   "outputs": [],
   "source": [
    "def compute_q_q_pairs(url):\n",
    "    \"\"\"Computes lists of uniform-distribution-quantiles and data-quantiles for a page.\n",
    "    \n",
    "    Args:\n",
    "      url (str): The URL of a page.\n",
    "    \n",
    "    Returns:\n",
    "      2-tuple: A pair (AKA a 2-tuple).  Both components are lists or arrays of\n",
    "               length n, where n is the number of positive checks for the page.\n",
    "               The first component contains quantiles of a uniform distribution,\n",
    "               and the second contains quantiles of the distribution of positive\n",
    "               check times for the page.  The kth element of a component is the\n",
    "               (k+1)/(n+1)-th quantile of its distribution or dataset.\n",
    "               \n",
    "               The first component of the pair contains quantiles of the uniform\n",
    "               distribution on [0, N], where N is the number of checks performed\n",
    "               for the page.\n",
    "    \"\"\"\n",
    "    updates = crawls.loc[url]['updated']\n",
    "    obs_q = np.array(updates[updates].index)\n",
    "    N = len(updates)\n",
    "    n = len(obs_q)\n",
    "    pdf_q = np.linspace(N/(n+1), N, n)\n",
    "    return (pdf_q, obs_q)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = ok.grade('q07')\n",
    "_ = ok.backup()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following code constructs the final Q-Q plot for the `\"Normal\"` sites defined above"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for url in pages_to_display[pages_to_display['Desc'] == \"Normal\"].index.values:\n",
    "    print(url)\n",
    "    (pdf_q, obs_q) = compute_q_q_pairs(url)\n",
    "    sns.regplot(pdf_q, np.array(obs_q))\n",
    "    plt.xlabel(\"True uniform quantile\")\n",
    "    plt.ylabel(\"Actual quantile\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part 2\n",
    "\n",
    "Describe your findings. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "for_assignment_type": "solution"
   },
   "outputs": [],
   "source": [
    "# SOLUTION\n",
    "q7_p2_answer = r\"\"\"\n",
    "\n",
    "**SOLUTION:** The page we examined on this site appear to have roughly uniform updates.\n",
    "There are deviations from a straight line, but it's hard to tell whether could have\n",
    "happened by chance.\n",
    "\n",
    "\"\"\"\n",
    "\n",
    "display(Markdown(q7_p2_answer))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "---\n",
    "\n",
    "## Question 8\n",
    "\n",
    "Even if the updates were distributed uniformly, the Q-Q plot will not appear as exactly a straight line.  To get a sense of what a uniform-quantile plot might look like if the data were truly distributed according to the uniform distribution, simulate $n$ observations from the uniform distribution on the interval (0, 719) and make a uniform-quantile plot for these simulated data.  $n$ is defined in the next cell; write your code in the cell below that."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "url = 'http://a.hatena.ne.jp/Syako/simple'\n",
    "n = np.count_nonzero(crawls.loc[url]['updated'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "for_assignment_type": "solution"
   },
   "outputs": [],
   "source": [
    "# SOLUTION\n",
    "obs_q = np.round(np.random.uniform(0, 719, n))\n",
    "obs_q.sort()\n",
    "\n",
    "N = 719\n",
    "pdf_q = np.arange(N/(n+1), N, N/(n+1))\n",
    "\n",
    "sns.regplot(pdf_q, np.array(obs_q))\n",
    "plt.xlabel(\"True uniform quantile\")\n",
    "plt.ylabel(\"Actual quantile\")\n",
    "\n",
    "# Leave this for grading purposes\n",
    "q8_plot = plt.gcf()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How do the empirical quantile plots from the previous quesion compare to your simulated quantile plot?  *Optionally,* we suggest looking at a few other sites and a few other simulated sets of data to see how well they match."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "for_assignment_type": "solution"
   },
   "outputs": [],
   "source": [
    "q8_answer = r\"\"\"\n",
    "\n",
    "**SOLUTION:** Our data drawn from a uniform distribution produced a Q-Q plot\n",
    "that's a *little* more \"line-like\" than the Q-Q plot for the actual positive\n",
    "checks for the first page in the previous question.  The difference was not so\n",
    "big; we would need to run a more formal test to see if the positive checks for\n",
    "the first page were really abnormal.  The other pages were even more\n",
    "uniform than our simulated data.\n",
    "\n",
    "\"\"\"\n",
    "\n",
    "display(Markdown(q8_answer))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "---\n",
    "\n",
    "# Estimating the Update Rate: A Simple Approach\n",
    "How would you estimate the change rate for a page? \n",
    "\n",
    "For example, imagine that in 720 hourly visits to a page, we observed a change in the page on 36 visits. One estimate for the rate of changes is:\n",
    "\n",
    "$$\\frac {36\\text{ changes}}{720\\text{ hours}} = \\frac{1}{20} \\text{ changes per hour}$$\n",
    "\n",
    "There is a small problem with our estimate. What if a page changes twice in an hour? We would see only one positive check.  We do not observe the true number of changes that happened to a page in 30 days, but the number of hours that had a change.  Think about how this could affect our interpretation of how often websites change.\n",
    "\n",
    "To help answer this question, we use a probability model for the change behavior of a site and examine the impact of the incomplete observations on our simple estimate of a site's rate of change. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A Model for Page Updates\n",
    "What model should we use?\n",
    "\n",
    "In our earlier Q-Q plot analysis we found that the updates appear to occur uniformly at random.  \n",
    "\n",
    "However, the number of positive checks (and therefore the number of page changes) is not the same from one page to the next. That is, both the *number* of changes and the hours of the changes appear to be random.\n",
    "\n",
    "When events can happen independently at any point in a span of time, and they're equally likely to happen at any of the times, a good model for the *number* of events that happen is the Poisson distribution.\n",
    "\n",
    "The Poisson distribution has a simple probability mass function, \n",
    "\n",
    "$$P(k) = \\frac{\\lambda^k}{k!}  e^{-\\lambda},$$\n",
    "\n",
    "for $k = 0, 1, 2, \\ldots$. The parameter $\\lambda$ is called the *rate*, and this is a Poisson$(\\lambda)$ distribution.\n",
    "\n",
    "For example, if $\\lambda$ is the hourly rate of changes to a page, then the chance of $k$ changes in one hour is:\n",
    "\n",
    "$$P(k) = \\frac{\\lambda^k}{k!}  e^{-\\lambda}.$$\n",
    "\n",
    "If we count the number of changes over $N$ hours, then the number of changes in this time period has \n",
    "a Poisson$(N \\lambda)$ distribution. That is,  \n",
    "\n",
    "$$P(k \\texttt{ updates in }N \\texttt{ hours}) = \\frac{(N\\lambda)^k}{k!}  e^{-(N\\lambda)}.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "---\n",
    "\n",
    "## Question 9\n",
    "Suppose that we observe $n$ changes for a site that was visited for $N$ hours. If the number of changes follows the Poisson distribution, show that $\\frac{n}{N}$ is the MLE for $\\lambda$.  In fancy notation, we could call this $\\lambda^{\\text{MLE}}$.  (Note that we are currently ignoring the problem of not being able to detect multiple changes in an hour.)\n",
    "\n",
    "*Note:* It's okay to write your answer in plain text.  If you know $\\LaTeX$ or would like to learn it, now is a good time to try it out.  (If you double-click on this cell or some of the cells above, you can see some $\\LaTeX$ examples.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "for_assignment_type": "solution"
   },
   "outputs": [],
   "source": [
    "q9_answer = r\"\"\"\n",
    "\n",
    "The probability of the data given $\\lambda$ is:\n",
    "\n",
    "$$L(\\lambda) = e^{-(\\lambda N)} \\frac{(\\lambda N)^{n}}{n!}$$\n",
    "\n",
    "...or in log form:\n",
    "\n",
    "$$\\log L(\\lambda) = -(\\lambda N) + n \\log(\\lambda) + n \\log(N) - \\log(n!)$$\n",
    "\n",
    "Setting the derivative to 0, we find:\n",
    "\n",
    "$$0 = \\left[\\frac{d}{d \\lambda} \\log L \\right](\\lambda^{\\text{MLE}}) = \n",
    "-N + \\frac{n}{\\lambda^{\\text{MLE}}}$$\n",
    "\n",
    "(Notice that the terms in $\\log L$ that didn't involve $\\lambda$, like $\\log(n!)$,\n",
    "conveniently disappeared when we took the derivative of the log-likelihood.)\n",
    "\n",
    "Therefore,\n",
    "\n",
    "$$\\lambda^{\\text{MLE}} = \\frac{n}{N}$$\n",
    "\"\"\"\n",
    "\n",
    "display(Markdown(q9_answer))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "---\n",
    "\n",
    "## Question 10\n",
    "\n",
    "1. Add a `'simple mle'` column to the `crawl_stats` table containing the `'simple mle'` estimator we derived earlier for each website.\n",
    "\n",
    "2. **Then,** make a plot that displays the distribution of these MLEs *for the sites with at least 700 crawls*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "for_assignment_type": "solution",
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "crawl_stats['simple mle'] = crawl_stats['proportion of updates']\n",
    "\n",
    "crawl_stats[crawl_stats['number of crawls'] >= 700]['simple mle'].hist(bins=20)\n",
    "plt.title(\"The distribution of simplified MLEs for the change rate for each site\");\n",
    "\n",
    "q10_plot = plt.gcf()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = ok.grade('q10')\n",
    "_ = ok.backup()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Impact of Hourly Observations\n",
    "\n",
    "The histogram that you made for previous problem has a small mode at 1. Why is this? It is not possible for our estimate of the rate of changes to be greater than once an hour because we only observe the page once an hour. So if the rate of changes is large, we are likely to underestimate it. Let's try to assess the impact of this.\n",
    "\n",
    "We will carry out a Monte Carlo (repeated, randomized simulation) study of the process that generates our observations. To perform such a study, we need a model for how the checks came out. For that, we need a model for the page changes themselves.\n",
    "\n",
    "We previously said that the number of changes follows a Poisson distribution with a certain rate for each page.  We also saw that the positive checks are distributed roughly uniformly.  It seems reasonable to assume that the changes themselves are also distributed uniformly.\n",
    "\n",
    "When the *number* of events follows the Poisson distribution and the *location* of the events (along the time interval) follows the uniform distribution, then we have what is called a *Poisson process*.  This is a random object that's a *set* of numbers rather than just one *number*.\n",
    "\n",
    "To perform our Monte Carlo study, for various values of $\\lambda$ we will generate the simulated changes according to the Poisson process distribution and then *reduce* the changes to the hour of the change. More specifically, we can simulate the data as follows:\n",
    "\n",
    "* Generate $M$, the number of updates, by drawing from the Poisson(720$\\lambda$) distribution.\n",
    "* Place each of the $M$ updates uniformly at random on the interval (0, 719). \n",
    "* \"Snap\" each update to the next hour.  For example, a 3:15 update time becomes 4:00, which is the time when we would have observed it. (What do you do when more than one update occurs within an hour?)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## Question 11\n",
    "Following the above description of the Poisson process, complete the function `sample_poisson_process`.\n",
    "\n",
    "*Hint:* The function `np.random.poisson` will be useful."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "for_assignment_type": "solution"
   },
   "outputs": [],
   "source": [
    "def sample_poisson_process(rate, length):\n",
    "    \"\"\"Generates n points from Poisson(rate*length) and locates them\n",
    "    uniformly at random on [0, length].\n",
    "    \n",
    "    Args:\n",
    "      rate (float): The average number of points per unit length.\n",
    "      length (float): The length of the line segment.\n",
    "    \n",
    "    Returns:\n",
    "      ndarray: An array of points scattered randomly on [0, length].\n",
    "               The number of points has a Poisson(rate*length)\n",
    "               distribution.\n",
    "    \"\"\"\n",
    "    total = np.random.poisson(rate*length)\n",
    "    def random_point():\n",
    "        return length*np.random.random()\n",
    "    return np.array(sorted([random_point() for _ in range(total)]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = ok.grade('q11')\n",
    "_ = ok.backup()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `snap_times` function in the next cell will help you simulate the hourly observations from the crawler."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def snap_times(update_times, window_length, process_length):\n",
    "    \"\"\"Given a list of change times, produces a list of the windows that detected a\n",
    "    change (that is, where at least one change happened inside the window).\n",
    "    \n",
    "    This has the effect of 'snapping' each change to the next hour (or whatever the\n",
    "    window_length is).  For periods where more than one change happened, the output\n",
    "    will still only list the period once.\n",
    "    \n",
    "    In other words, it produces a list of the positive checks given a list of true\n",
    "    change times.\n",
    "    \n",
    "    Args:\n",
    "      update_times (ndarray): A list of times when changes happened for a page.\n",
    "                              All times are between 0 and process_length.\n",
    "      window_length (float): The width of each window.  (For example, if time\n",
    "                             is being measured in hours, and an observation\n",
    "                             happens once per hour, then this should be 1.)\n",
    "      process_length (float): The last time time any change could have happened.\n",
    "    \n",
    "    Returns:\n",
    "      ndarray: A list of numbers, each one the right endpoint of a window where at\n",
    "               least one change happened.\"\"\"\n",
    "    window_ends = np.arange(0, process_length, window_length) + window_length\n",
    "    num_windows = len(window_ends)\n",
    "    event_windows = np.floor(np.array(update_times) / window_length).astype(int)\n",
    "    events_by_window = np.bincount(event_windows, minlength=num_windows)\n",
    "    event_happened = events_by_window > 0\n",
    "    return window_ends[event_happened]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "---\n",
    "\n",
    "## Question 12\n",
    "\n",
    "Use the functions `sample_poisson_process` and `snap_times` to examine what happens when we visit hourly.  Look at examples where\n",
    "* the `length` of time is 24 hours,\n",
    "* the `rate` is 1/8, 1/4, 1/2, 1, and 2 changes per hour, and\n",
    "* the `window_length` is 1.\n",
    "\n",
    "For each value of `rate`, simulate one set of change on the interval (0 hours, 24 hours), and plot the resulting points on a timeline using `display_points`.  Then snap these data to the next hour with `snap_times` and plot the resulting positive check times.\n",
    "\n",
    "**Then,** compare the actual change times to the censored times for the various rates. What happens as the rate changes? Do you think hourly visits to the site are a problem if the rate is 1/8? How about 2? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "for_assignment_type": "solution"
   },
   "outputs": [],
   "source": [
    "for rate in [1/8, 1/4, 1/2, 1, 2]:\n",
    "    draws = sample_poisson_process(rate, 24)\n",
    "    windows = snap_times(draws, 1, 24)\n",
    "    display_points(draws, (0, 24), \"Actual update times (rate={:f})\".format(rate))\n",
    "    display_points(windows, (0, 24), \"Crawler checks that found updates\")\n",
    "q12_answer = r\"\"\"\n",
    "\n",
    "**SOLUTION:** It seems problematic for rates 1/2 or higher, but possibly okay if the rate is 1/8 or 1/4.\n",
    "\n",
    "\"\"\"\n",
    "\n",
    "display(Markdown(q12_answer))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulation study of the rate estimate \n",
    "Now let us extend our Monte Carlo study to compute the MLE from the observed positive checks.  Since we know the true change rate in the simulations, we can compare the MLE to the true change rate and see how accurate it is as an estimator.  (We couldn't have done that with our real data, since we don't know the true change rates.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "---\n",
    "\n",
    "## Question 13\n",
    "\n",
    "Suppose a website is visited every hour for 30 days, so it has 719 visits.  For a particular rate $\\lambda$, \n",
    "\n",
    "* Generate the observed positive checks with our earlier functions (`sample_poisson_process` and `snap_times`).\n",
    "* Calculate the MLE from these simulated data.\n",
    "* Repeat 1000 times.\n",
    "\n",
    "This average estimate will change as we increase the true rate of changes, so vary $\\lambda$ between 0 and 4 and make a plot showing the estimate as a function of $\\lambda$.  For comparison, show the true change rates in your plot as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "for_assignment_type": "solution"
   },
   "outputs": [],
   "source": [
    "def simulate_censored_rate_estimate(true_rate, num_visits):\n",
    "    \"\"\"Simulates updates to a page and visits by a web crawler, and\n",
    "    returns the proportion of visits in which an update was observed.\n",
    "    \n",
    "    Args:\n",
    "      true_rate (float): The average number of updates per unit length\n",
    "                         of time.  (The units are irrelevant for the\n",
    "                         purposes of this function, but you can imagine\n",
    "                         that they are hours.)\n",
    "      num_visits (float): The number of visits made to the site.  One\n",
    "                          visit is made per unit length of time, so\n",
    "                          this is also equal to the duration of time\n",
    "                          simulated.\n",
    "    \n",
    "    Returns:\n",
    "      float: the MLE for true_rate lambda, based on the number of \n",
    "      observed positive checks.\n",
    "    \"\"\"\n",
    "    draws = sample_poisson_process(true_rate, num_visits)\n",
    "    windows = snap_times(draws, 1, num_visits)\n",
    "    n = len(windows)\n",
    "    return n / num_visits\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following helper function can be used to simulate many rate estimates using the function you just completed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def simulate_many_rate_estimates(true_rate, num_visits):\n",
    "    return np.mean([simulate_censored_rate_estimate(true_rate, num_visits) \n",
    "                    for _ in range(1000)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following code will simulate rate estimates for various values of $\\lambda$ (this cell may take up to 1 minute to run):"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, the following code will plot the estimated values for $\\lambda$ against their true values.  You will need this plot to answer the next question."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = ok.grade('q13')\n",
    "_ = ok.backup()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## Question 14\n",
    "\n",
    "### Part 1\n",
    "Explain why the estimated rates seem to level off at 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "for_assignment_type": "solution"
   },
   "outputs": [],
   "source": [
    "q14_p1_answer = r\"\"\"\n",
    "\n",
    "**SOLUTION:**\n",
    "\n",
    "The proportion of hours in which an update was observed can never be greater \n",
    "than 1 because we can't observe updates more than once an hour.  \n",
    "The estimate for $\\lambda$ is closest to the true value for small $\\lambda$. \n",
    "For $\\lambda$ of 0.25, the estimate is off by less than 5%. \n",
    "\n",
    "\n",
    "\"\"\"\n",
    "\n",
    "display(Markdown(q14_p1_answer))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part 2\n",
    "How far off is the estimate from the truth for $\\lambda$ less than 0.25?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "for_assignment_type": "solution"
   },
   "outputs": [],
   "source": [
    "q14_p2_answer = r\"\"\"\n",
    "\n",
    "**SOLUTION:** It's not so far off: at most ~20%.\n",
    "\n",
    "\"\"\"\n",
    "\n",
    "display(Markdown(q14_p2_answer))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Estimating the update rate: An improved approach\n",
    "Our simulation study has shown us that our MLE estimate for the rate is biased. It systematically underestimates the quantity we're trying to estimate, the rate.  This bias is small for small $\\lambda$, but can we eliminate it?\n",
    "\n",
    "We can recast the problem slightly by reconsidering our observations. Although the possible times of a page changes may be continuous, we only observe the data hourly and we only observe whether there was at least one update in that hour.  Does this remind you of a distribution that you have worked with?\n",
    "\n",
    "What if we consider the data for a page as 719 0-1 values, where 0 indicates no change and 1 indicates at least 1 change in the hour?  Do you see why these 0-1 values are equivalent to our representation of the data as positive changes?  For example, if `N` is 10 and we observe changes at hours 3, 7, and 9 then this information is equivalent to the sequence 0,0,1,0,0,0,1,0,1,0.  \n",
    "\n",
    "The distribution of the 0-1 values is Bernoulli(`p`). Under our assumption that the changes come from a Poisson process, their values are independent.\n",
    "\n",
    "This distribution incorporates the censoring of the data. How can we estimate the rate? We can use the Poisson distribution to compute the chance of a 0 or 1, which depends on the rate. Then we use the Bernoulli distribution to find the MLE of $\\lambda$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## Question 15\n",
    "What is the chance that a Poisson($\\lambda$) random variable is equal to 0?  What is the chance that it's greater than or equal to 1?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "for_assignment_type": "solution"
   },
   "outputs": [],
   "source": [
    "# SOLUTION\n",
    "q15_answer = r\"\"\"\n",
    "\n",
    "**SOLUTION:** $e^{-\\lambda}$ and $1 - e^{-\\lambda}$, respectively.\n",
    "\n",
    "\"\"\"\n",
    "\n",
    "display(Markdown(q15_answer))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To check your answers, fill in the numerical values of each probability, for $\\lambda=0.5$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# The probability that a Poisson(0.5) random variable is equal to 0\n",
    "Prob_pois_half_equals_0 = ...\n",
    "# The probability that a Poisson(0.5) random variable is greater than or equal to 1\n",
    "Prob_pois_half_gte_1 = ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = ok.grade('q15')\n",
    "_ = ok.backup()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## Question 16\n",
    "With this modified model, the MLE of $\\lambda$, given $n$ updates in $N$ visits for a site, is:\n",
    "\n",
    "$$\\lambda^* = \\log \\left(\\frac{n}{N - n} + 1 \\right)$$\n",
    "\n",
    "Show that that's true.\n",
    "\n",
    "*Hint:* If you find an expression like $\\frac{e^{-\\lambda}}{1 - e^{-\\lambda}}$, it often helps to convert such an expression to $\\frac{1}{e^{\\lambda} - 1}$ by multiplying its numerator and denominator by $e^{\\lambda}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "for_assignment_type": "solution"
   },
   "outputs": [],
   "source": [
    "q16_answer = r\"\"\"\n",
    "\n",
    "$$L(\\lambda) = (e^{-\\lambda})^{N - n} (1 - e^{-\\lambda})^{n}$$\n",
    "\n",
    "$$\\log L(\\lambda) = (N - n) (-\\lambda) + n \\log(1 - e^{-\\lambda})$$\n",
    "\n",
    "$$0 = \\left[\\frac{d}{d \\lambda} \\log L \\right](\\lambda^*) = -(N - n) + n \\frac{e^{-\\lambda^*}}{1 - e^{-\\lambda^*}}$$\n",
    "\n",
    "$$e^{\\lambda^*} = \\frac{n}{N - n} + 1$$\n",
    "\n",
    "$$\\lambda^* = \\log \\left(\\frac{n}{N - n} + 1 \\right)$$\n",
    "\n",
    "\"\"\"\n",
    "\n",
    "display(Markdown(q16_answer))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What happens when we observe a change every hour?  Then, our MLE $\\lambda^*$ takes the value $\\infty$.  \n",
    "\n",
    "Why is that true? Intuitively, if there is a change on each visit, we can always increase the likelihood of our data (which is the likelihood that every hour saw at least one change) by increasing $\\lambda$.\n",
    "\n",
    "An MLE of $\\infty$ is clearly problematic. We can use the idea of Laplace smoothing to modify the estimator a bit.  We'll use the following \"modified MLE\" estimator, which adds to every site a single imagined visit where \"half\" saw a change and the other \"half\" didn't:\n",
    "\n",
    "$$\\lambda^+ = \\log \\left(\\frac{n + .5}{N - n + .5} + 1 \\right)$$\n",
    "\n",
    "You can imagine that this estimate incorporates a prior conservative intuition that a site might not change all the time even if it happened to change every hour in the particular 30-day period we observed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "---\n",
    "\n",
    "\n",
    "## Question 17\n",
    "\n",
    "\n",
    "1. Add a `'modified mle'` column to the `crawl_stats` table containing the `'modified mle'` estimator $\\lambda^+$ for each website.\n",
    "\n",
    "2. **Then,** make a histogram that displays the distribution of the modified MLEs *for the sites with at least 700 crawls*.\n",
    "\n",
    "*Hint:* You may find your work in question 10 helpful."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "for_assignment_type": "solution"
   },
   "outputs": [],
   "source": [
    "crawl_stats[\"modified mle\"] = np.log(\n",
    "    (crawl_stats['number of updates'] + .5)\n",
    "    / (crawl_stats['number of crawls'] - crawl_stats['number of updates'] + .5)\n",
    "    + 1)\n",
    "crawl_stats[\"modified mle\"].hist(bins=20);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = ok.grade('q17')\n",
    "_ = ok.backup()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that this distribution is quite different than the earlier one we found in question 10. Now there is not a pile up of estimates at 1. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How accurate are our estimates?\n",
    "\n",
    "We don't know the true update rates of the pages, so we can't know exactly how accurate our estimates are.  But this is often the case in data science.  Let's try to figure out a reasonable way to *estimate* the accuracy of our estimates.\n",
    "\n",
    "The strategy is similar to the bootstrap.  In pseudocode:\n",
    "\n",
    "`for each page ` $i$ `:`\n",
    "\n",
    "`... Assume the true change rate ` $r_i$ ` equals our estimate ` $\\lambda^+_i$.\n",
    "\n",
    "`... for _ in range(num_simulations):`\n",
    "\n",
    "`...     Simulate changes and positive checks.`\n",
    "\n",
    "`...     Compute a new maximum likelihood estimate of ` $r_i$ ` based on the simulated data.`\n",
    "\n",
    "`... Plot these estimates in a histogram or summarize them in some other way.`\n",
    "\n",
    "This is called a *parametric* bootstrap method, because instead of resampling directly from our data, we use our data to estimate its distribution's parameters and then sample from that distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## Question 18\n",
    "Complete the definitions of `simulate_change_rate_estimate` and `modified_mle` in the cell below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "for_assignment_type": "solution"
   },
   "outputs": [],
   "source": [
    "def simulate_change_rate_estimate(change_rate, num_observations, estimator):\n",
    "    \"\"\"Simulates hourly change observations for a website and produces an\n",
    "    estimate of the change rate.\n",
    "    \n",
    "    Args:\n",
    "      change_rate (float): The hourly change rate of the website to be used in\n",
    "                           the simulation.\n",
    "      num_observations (int): The number of observations (equivalently, the\n",
    "                              number of hours) to simulate.\n",
    "      estimator (func): A function to apply to the simulated observations to\n",
    "                        produce an estimate of the change rate.  Has the same\n",
    "                        signature as the function modified_mle below.\n",
    "    \"\"\"\n",
    "    draws = sample_poisson_process(change_rate, num_observations)\n",
    "    windows = snap_times(draws, 1, num_observations)\n",
    "    return estimator(windows, num_observations)\n",
    "\n",
    "def modified_mle(positive_checks, num_observations):\n",
    "    \"\"\"Produces the modified MLE for a dataset.\n",
    "    \n",
    "    Args:\n",
    "      positive_checks (ndarray): A list or array of the hours when a\n",
    "                                 change was observed.\n",
    "      num_observations (int): The number of hours in which we could\n",
    "                              have observed a change.\n",
    "    \n",
    "    Returns:\n",
    "      (float): The modified MLE.\n",
    "    \"\"\"\n",
    "    num_updates = len(positive_checks)\n",
    "    return np.log((num_updates + .5) / (num_observations - num_updates + .5) + 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = ok.grade('q18')\n",
    "_ = ok.backup()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## Question 19\n",
    "\n",
    "Complete the definition of `plot_bootstrap_rate_estimates` in the cell below.  **Read the docstring** carefully so you know what to do. **Then** run the provided cell below that to plot the bootstrap rate estimate distributions for a few pages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "for_assignment_type": "solution"
   },
   "outputs": [],
   "source": [
    "def plot_bootstrap_rate_estimates(page_url, num_simulations):\n",
    "    \"\"\"Simulates positive check observations for a website many times.  For\n",
    "    each simulation, the page's change rate is estimated based on the simulated\n",
    "    positive checks.  The estimation method is the modified_mle function\n",
    "    implemented above.  Then this function produces a histogram that shows\n",
    "    the distribution of these estimated change rates.\n",
    "    \n",
    "    When conducting each simulation, the change rate is taken to be the\n",
    "    estimated update rate for the page.  *That* estimate is the modified\n",
    "    MLE computed from the actual positive check observations for that page.\n",
    "    \n",
    "    The modified MLE computed from the actual observations for the page\n",
    "    is also displayed as a vertical red line at an appropriate horizontal\n",
    "    position on the histogram.\n",
    "    \n",
    "    Here is a diagram that might be helpful:\n",
    "    \n",
    "                                             (This happens 1 time)\n",
    "                                                     |\n",
    "                                                     v\n",
    "    Positive check observations for this page ---modified_mle---> Estimate of change rate for this page\n",
    "    \n",
    "                                         (once per simulation)\n",
    "    Estimate of change rate for this page ---simulation---> Simulated positive checks\n",
    "    \n",
    "                              (once per simulation)\n",
    "    Simulated positive checks ---modified_mle---> Bootstrap estimate of change rate for this page\n",
    "    \n",
    "    1000 bootstrap estimates of change rate for this page ---> Histogram\n",
    "                                                               ^\n",
    "                                                               |\n",
    "    Estimate of change rate for this page --------------------/\n",
    "                                              (displayed somewhere on the histogram)\n",
    "    \n",
    "    Args:\n",
    "      page_url (str): The URL of the page to simulate.\n",
    "      num_simulations (int): The number of simulations to run.\n",
    "    \n",
    "    Returns:\n",
    "      None: This function doesn't return anything; it just makes a histogram\n",
    "      appear as described above.\n",
    "    \"\"\"\n",
    "    update_rate_estimate = crawl_stats.loc[page_url][\"modified mle\"]\n",
    "    num_observations = crawl_stats.loc[page_url]['number of crawls']\n",
    "    bootstrap_estimates = np.array([\n",
    "        simulate_change_rate_estimate(update_rate_estimate, num_observations, modified_mle) \n",
    "        for _ in range(num_simulations)])\n",
    "    plt.hist(bootstrap_estimates, bins=20)\n",
    "    plt.axvline(crawl_stats.loc[page_url][\"modified mle\"], color='red')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run this cell to make plots for several pages using your\n",
    "# function.  Since no automatic tests are provided for this\n",
    "# question, we suggest examining the plots to make sure they\n",
    "# make sense to you.\n",
    "many_updates_url = crawl_stats[np.logical_and(\n",
    "    crawl_stats['number of crawls'] >= 700,\n",
    "    np.logical_and(\n",
    "        .2 <= crawl_stats['proportion of updates'],\n",
    "        crawl_stats['proportion of updates'] >= .8))].index[0]\n",
    "plot_bootstrap_rate_estimates(many_updates_url, 1000)\n",
    "\n",
    "few_updates_url = crawl_stats[np.logical_and(\n",
    "    crawl_stats['number of crawls'] >= 700,\n",
    "    np.logical_and(\n",
    "        .05 <= crawl_stats['proportion of updates'],\n",
    "        crawl_stats['proportion of updates'] >= .15))].index[0]\n",
    "plot_bootstrap_rate_estimates(few_updates_url, 1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looking at the error distribution is good, but it's also useful to characterize error with a single number.  (This facilitates comparisons among different estimators, for example.)\n",
    "\n",
    "A common way to summarize error in estimation is the root mean squared error, or RMSE.  The RMSE is the square root of the MSE, or mean square error. We calculate the MSE as follows:\n",
    "\n",
    "* the average (across many simulations, where each simulation results in 1 estimate) of\n",
    "* the squared difference between the estimate and its actual value.\n",
    "\n",
    "Working with the RMSE, rather than the MSE, is useful when we want to understand the magnitude of error relative to the estimates.\n",
    "\n",
    "We don't know the true change rate for each website, so we can't compute the true MSE.  However, if we find the MSE using the data making up the histogram you produced in the previous question, using the original estimate (the red line) as the actual change rate, that is a *bootstrap estimate* of the RMSE."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 20\n",
    "\n",
    "---\n",
    "\n",
    "Complete the definition of `estimate_rmse` in the cell below.  **Then**, use it to estimate the RMSE for each page's modified MLE.  Add these RMSEs as a column named `'rmse'` in `crawl_stats`.\n",
    "\n",
    "**Note**: There are 1000 pages, so this will take some time.  We recommend using only 100 simulations per page, which should be sufficient for reasonable estimates.  Try your code on a single page before running it on every page."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "for_assignment_type": "solution"
   },
   "outputs": [],
   "source": [
    "def estimate_rmse(page_url, num_simulations):\n",
    "    \"\"\"Simulates update observations for a website many times.  For each\n",
    "    simulation, the page's change rate is estimated based on the simulated\n",
    "    observations.  (The estimation method is the modified MLE.)  Then this\n",
    "    function produces an estimate of the RMSE of estimates of the change\n",
    "    rate for this page.\n",
    "    \n",
    "    When conducting each simulation, the change rate is taken to be the\n",
    "    estimated change rate for the page.  *That* estimate is the modified\n",
    "    MLE computed from the actual observations for the page.\n",
    "    \n",
    "    We compute the modified MLE for each set of simulated observations.  That\n",
    "    constitutes num_simulations estimates of the change rate.\n",
    "    \n",
    "    Then we compute the RMSE of those estimates.  The \"true\" change rate in\n",
    "    that calculation is taken to be the modified MLE computed from the\n",
    "    actual observations for the page.\n",
    "    \n",
    "    Args:\n",
    "      page_url (str): The URL of the page to simulate.\n",
    "      num_simulations (int): The number of simulations to run.\n",
    "    \n",
    "    Returns:\n",
    "      float: The estimated RMSE of the modified MLE for the given page,\n",
    "             based on num_simulations simulations.\n",
    "    \"\"\"\n",
    "    change_rate_estimate = crawl_stats.loc[page_url][\"modified mle\"]\n",
    "    num_observations = crawl_stats.loc[page_url]['number of crawls']\n",
    "    bootstrap_estimates = np.array(\n",
    "        [simulate_change_rate_estimate(change_rate_estimate, num_observations, modified_mle) \n",
    "         for _ in range(num_simulations)])\n",
    "    return np.mean((bootstrap_estimates - change_rate_estimate)**2)**0.5\n",
    "\n",
    "print(estimate_rmse(many_updates_url, 1000))\n",
    "print(estimate_rmse(few_updates_url, 1000))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "for_assignment_type": "solution"
   },
   "outputs": [],
   "source": [
    "crawl_stats['rmse'] = [estimate_rmse(i, 100) for i in crawl_stats.index]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = ok.grade('q20')\n",
    "_ = ok.backup()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## Question 21\n",
    "Create a visualization to display the RMSEs you computed.  **Then,** create another visualization to see the relationship (across the 1000 pages) between RMSE and the modified MLE."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "for_assignment_type": "solution",
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "crawl_stats[\"rmse\"].hist(bins=20);\n",
    "crawl_stats[crawl_stats['number of crawls'] == 719].plot.scatter(\"modified mle\", \"rmse\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Is the model reasonable?\n",
    "All the foregoing work has assumed a particular probability model for website updates: the number of changes follows a Poisson($\\lambda$) distribution and the locations of the changes, given the number of changes, follows a uniform distribution.  Like most models, this is certainly imperfect.  However, it still may produce reasonable estimates. Let's check some of the conditions of the model to see if they are reasonable.\n",
    "\n",
    "Embedded in our model is the assumption that the rate of changes for a page is constant over time.  If we find a strong pattern in the rates of updates, that is evidence against this assumption.  For example, you might guess that changes for certain pages might happen more often at certain times of the day, or on certain days of the week.\n",
    "\n",
    "Before you group by day of week, here are some difficulties presented by the dataset:\n",
    "\n",
    "1. The day of the week is not given.\n",
    "2. The hour of the day is not given.\n",
    "3. We don't know when the checks started.\n",
    "4. When a site had a failed check, the check was just omitted from the dataset, so we don't know when the failed checks happened, only their total count for each page.\n",
    "5. Not all sites are in the same or nearby time zones.\n",
    "\n",
    "Here is what we *do* know:\n",
    "\n",
    "1. Checks for all pages started at the same time, at midnight PST (US Pacific Coast time).\n",
    "2. Checks were attempted every hour for 30 days, and then they stopped.\n",
    "3. We know each page's URL.  A URL may tell us something about its geography and therefore its rough time zone.  For example, the top-level domains `.edu`, `.com`,`.gov`, and `.net` include only US sites."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## Question 22\n",
    "Propose a ***plan*** to answer one of these questions:\n",
    "\n",
    "1. How much did the rate of changes vary by hour of the day?\n",
    "2. How much did the rate of changes vary by day of the week?\n",
    "\n",
    "**Or,** you may choose your own pattern to test.\n",
    "\n",
    "Your plan doesn't need to produce results that hold for *all* pages, only a reasonably-large subset.  Be sure to make clear exactly what question your analysis would answer.\n",
    "\n",
    "**Note**: You don't need to actually implement your plan!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "for_assignment_type": "solution"
   },
   "outputs": [],
   "source": [
    "q22_answer = r\"\"\"\n",
    "\n",
    "TODO: Write a solution.\n",
    "\n",
    "\"\"\"\n",
    "\n",
    "display(Markdown(q22_answer))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Submitting your assignment\n",
    "Congratulations, you're done with this homework!\n",
    "\n",
    "Run the next cell to run all the tests at once."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "_ = ok.grade_all()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Finally**, run the next cell to submit the assignment to OkPy so that the staff will know to grade it. You can submit as many times as you want, and you can choose which submission you want us to grade by going to https://okpy.org/cal/data100/sp17/.  After you've done that, make sure you've pushed your changes to Github as well!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "_ = ok.submit()"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
