time_slopes.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "64fcf37c",
   "metadata": {},
   "source": [
    "# The Age of the Universe\n",
    "\n",
    "This exercise is from the Summer 2017 [Berkeley course\n",
    "materials](https://github.com/data-8/data8assets) with thanks.\n",
    "\n",
    "Sometimes, the primary purpose of regression analysis is to learn something about the slope or intercept of the best-fitting line.  When we use a sample of data to estimate the slope or intercept, our estimate is subject to random error, just as in the simpler case of the mean of a random sample.\n",
    "\n",
    "In this lab, we'll use regression to get an accurate estimate for the age of the universe, using pictures of exploding stars.  Our estimate will come from a sample of all exploding stars. We'll compute a confidence interval to quantify the error caused by sampling."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "94f502c5",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "# Run this cell, but please don't change it.\n",
    "\n",
    "# These lines import the Numpy and Pandas modules.\n",
    "import numpy as np\n",
    "# Make a random number generator.\n",
    "rng = np.random.default_rng()\n",
    "\n",
    "import pandas as pd\n",
    "# Safe setting for Pandas.  Needs Pandas version >= 1.5.\n",
    "pd.set_option('mode.copy_on_write', True)\n",
    "\n",
    "# Stardard plotting setup.\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "plt.style.use('fivethirtyeight')\n",
    "\n",
    "# Fancy stuff for some interactive plots\n",
    "import warnings\n",
    "# warnings.simplefilter('ignore', FutureWarning)\n",
    "from matplotlib import patches\n",
    "from ipywidgets import interact, interactive, fixed\n",
    "import ipywidgets as widgets\n",
    "\n",
    "# These lines load the tests.\n",
    "from client.api.notebook import Notebook\n",
    "ok = Notebook('time_slopes.ok')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "afebb8ef",
   "metadata": {},
   "source": [
    "\n",
    "## The Actual Big Bang Theory\n",
    "\n",
    "In the early 20th century, the most popular cosmological theory suggested that\n",
    "the universe had always existed at a fixed size.  Today, the Big Bang theory\n",
    "prevails: Our universe started out very small and is still expanding.\n",
    "\n",
    "A consequence of this is Hubble's Law, which says that the expansion of the\n",
    "universe creates the appearance that every celestial object that's reasonably\n",
    "far away from Earth (for example, another galaxy) is moving away from us at\n",
    "a constant speed.  If we extrapolate that motion backwards to the time when\n",
    "everything in the universe was in the same place, that time is (roughly) the\n",
    "beginning of the universe!\n",
    "\n",
    "Scientists have used this fact, along with measurements of the current\n",
    "*location* and *movement speed* of other celestial objects, to estimate when\n",
    "the universe started.\n",
    "\n",
    "The cell below simulates a universe in which our sun is the center and every\n",
    "other star is moving away from us.  Each star starts at the same place as the\n",
    "sun, then moves away from it over time.  Different stars have different\n",
    "directions *and speeds*; the arrows indicate the direction and speed of travel.\n",
    "\n",
    "Run the cell, then move the slider to see how things change over time.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7bf6e4ae",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Just run this cell.  (The simulation is actually not\n",
    "# that complicated; it just takes a lot of code to draw\n",
    "# everything.  So you don't need to read this unless you\n",
    "# have time and are curious about more advanced plotting.)\n",
    "\n",
    "num_locations = 15\n",
    "eg_x = rng.normal(size=num_locations)\n",
    "eg_y = rng.normal(size=num_locations)\n",
    "start_of_time = -2\n",
    "\n",
    "def scatter_after_time(t, start_of_time, end_of_time, x, y, center_name,\n",
    "                       other_point_name, make_title):\n",
    "    mx_abs_xy = np.max(np.abs(np.concatenate([x, y])))\n",
    "    max_location = (\n",
    "        1.1 * (end_of_time - start_of_time) * mx_abs_xy\n",
    "    )\n",
    "    new_x = (t - start_of_time) * x\n",
    "    new_y = (t - start_of_time) * y\n",
    "    n = len(x)\n",
    "    plt.scatter([0], [0], label=center_name, s=100, c=\"yellow\")\n",
    "    plt.scatter(new_x, new_y, label=other_point_name)\n",
    "    for i in np.arange(n):\n",
    "        plt.arrow(\n",
    "            new_x[i],\n",
    "            new_y[i],\n",
    "            x[i],\n",
    "            y[i],\n",
    "            fc='black',\n",
    "            ec='black',\n",
    "            head_width=0.025 * max_location,\n",
    "            lw=.15)\n",
    "    plt.xlim(-max_location, max_location)\n",
    "    plt.ylim(-max_location, max_location)\n",
    "    plt.gca().set_aspect('equal', adjustable='box')\n",
    "    plt.gca().set_position([0, 0, 1, 1])\n",
    "    plt.legend(bbox_to_anchor=(1.6, .7))\n",
    "    plt.title(make_title(t))\n",
    "    plt.show()\n",
    "\n",
    "interact(\n",
    "    scatter_after_time,\n",
    "    t=widgets.FloatSlider(min=start_of_time, max=5, step=.05, value=0,\n",
    "                          msg_throttle=1),\n",
    "    start_of_time=fixed(start_of_time),\n",
    "    end_of_time=fixed(5),\n",
    "    x=fixed(eg_x),\n",
    "    y=fixed(eg_y),\n",
    "    center_name=fixed(\"our sun\"),\n",
    "    other_point_name=fixed(\"other star\"),\n",
    "    make_title=fixed(lambda t: \"The world {:01g} year{} in the {}\".format(abs(t), \"\" if abs(t) == 1 else \"s\", \"past\" if t < 0 else \"future\"))\n",
    ");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f39b32c6",
   "metadata": {},
   "source": [
    "### Question 1\n",
    "\n",
    "When did the universe start, in this example?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27d6848b",
   "metadata": {
    "manual_grade": true,
    "manual_problem_id": "universe_start"
   },
   "source": [
    "*Write your answer here, replacing this text.*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e4d3105",
   "metadata": {},
   "source": [
    "### Question 2\n",
    "\n",
    "After 5 years (with the slider all the way to the right), stars with longer\n",
    "arrows are further away from the Sun.  Why?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cdf1fba3",
   "metadata": {
    "manual_grade": true,
    "manual_problem_id": "why_further"
   },
   "source": [
    "*Write your answer here, replacing this text.*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd298f0b",
   "metadata": {},
   "source": [
    "### Analogy: driving\n",
    "\n",
    "Here's an analogy to illustrate how scientists use information about stars to\n",
    "estimate the age of the universe.\n",
    "\n",
    "Suppose that at some point in the past, our friend Mei started driving in a car\n",
    "going at a steady speed of 60 miles per hour straight east.  We're still\n",
    "standing where she started."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b461375",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run this cell to see a picture of Mei's locations over time.\n",
    "mei_x = np.array([60])\n",
    "mei_y = np.array([0])\n",
    "\n",
    "interact(\n",
    "    scatter_after_time,\n",
    "    t=widgets.FloatSlider(min=-2, max=1, step=.05, value=0, msg_throttle=1),\n",
    "    start_of_time=fixed(-2),\n",
    "    end_of_time=fixed(1),\n",
    "    x=fixed(mei_x),\n",
    "    y=fixed(mei_y),\n",
    "    center_name=fixed(\"Us\"),\n",
    "    other_point_name=fixed(\"Mei\"),\n",
    "    make_title=fixed(lambda t: \"Mei's position {:01g} hour{} in the {}\".format(abs(t), \"\" if abs(t) == 1 else \"s\", \"past\" if t < 0 else \"future\")));"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "40bccb10",
   "metadata": {},
   "source": [
    "We want to know how long she's been driving, but we forgot to record the time\n",
    "when she left.  If we find out that she's 120 miles away, and she's been going\n",
    "60 miles per hour the whole time, we can infer that she left 2 hours ago.\n",
    "\n",
    "One way we can compute that number is by fitting a line to a scatter plot of\n",
    "our locations and speeds.  It turns out that the *slope* of that line is the\n",
    "amount of time that has passed.  Run the next cell to see a picture:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b7cb3322",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Just run this cell.\n",
    "mei_speed = 60\n",
    "mei_distance = 120\n",
    "# Add origin and Mei points\n",
    "plt.plot([0, mei_speed], [0, mei_distance], 'o')\n",
    "# Draw a line between them\n",
    "plt.plot([0, mei_speed], [0, mei_distance], ':')\n",
    "# Label points\n",
    "plt.text(0, 10, 'Us', fontsize=20)\n",
    "plt.text(mei_speed + 2, mei_distance + 2, 'Mei', fontsize=20)\n",
    "plt.text(mei_speed / 2 + 2, mei_distance / 2 - 4, 'Slope = 2 hours', color='blue', fontsize=20)\n",
    "# Label axes\n",
    "plt.xlabel(\"Speed moving away from us (miles per hour)\")\n",
    "plt.ylabel(\"Current distance from us (miles)\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef965ea2",
   "metadata": {},
   "source": [
    "The slope of the line is 2 hours.  (The units are vertical-axis units divided\n",
    "by horizontal-axis units, which are $\\frac{\\texttt{miles}}{\\texttt{miles} /\n",
    "\\texttt{hour}}$, or hours.)  So that's our answer.\n",
    "\n",
    "Imagine that you don't know Mei's *exact* distance or speed, only rough\n",
    "estimates.  Then if you drew this line, you'd get a slightly bad estimate of\n",
    "the time since she left.  But if you measured the distance and speed of\n",
    "hundreds of people who left you at the same time going different speeds, and\n",
    "drew a line through them, the slope of that line would be a pretty good\n",
    "estimate of the time they left, even if the individual measurements weren't\n",
    "exactly right.\n",
    "\n",
    "The `drivers.csv` dataset contains the speeds and distances-from-start of 100\n",
    "drivers.  They all left the same starting location at the same time, driving at\n",
    "a fixed speed on a straight line away from the start.  The measurements aren't\n",
    "exact, so they don't fit exactly on a line.  We've created a scatter plot and\n",
    "drawn a line through the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d73b0ea7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Just run this cell.\n",
    "drivers = pd.read_csv('drivers.csv')\n",
    "drivers.plot.scatter(x=0, y=1)  # Plot second column against first.\n",
    "# Draw a well-fitting line through the points.\n",
    "plt.plot([5, 75], [8.3333, 125], color='black')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e9858e8b",
   "metadata": {},
   "source": [
    "### Question 3\n",
    "\n",
    "By *looking at the fit line*, estimate how long ago (in hours) this cohort of\n",
    "drivers left."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cfc78243",
   "metadata": {},
   "outputs": [],
   "source": [
    "#- Fill in the start time you infer from the above line.\n",
    "driving_start_time_hours = ...\n",
    "driving_start_time_hours"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "467b44d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = ok.grade('q3')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef0abe68",
   "metadata": {},
   "source": [
    "### Back to cosmology\n",
    "\n",
    "To do the same thing for the universe, we need to know the distance-from-Earth\n",
    "and speed-away-from-Earth of many celestial objects.  Using pictures taken by\n",
    "very accurate telescopes and a lot of physics, astronomers have been able to\n",
    "estimate both.  It turns out that *nearby supernovae* -- stars that have\n",
    "recently died and exploded -- are among the best sources of this data, because\n",
    "they are very easy to see.  This picture taken by the Hubble telescope shows an\n",
    "entire galaxy, with a single supernova - as bright by itself as billions of\n",
    "stars - at the bottom left.\n",
    "\n",
    "<img src=\"supernova.jpg\">\n",
    "\n",
    "Our astronomical data for today will come from the [Supernova Cosmology\n",
    "Project](http://supernova.lbl.gov/union/) at Lawrence Berkeley Lab.  The\n",
    "original dataset is\n",
    "[here](http://supernova.lbl.gov/union/figures/SCPUnion2.1_mu_vs_z.txt), with\n",
    "(brief) documentation\n",
    "[here](http://supernova.lbl.gov/union/descriptions.html#Magvsz).  Each row in\n",
    "the table corresponds to a supernova near Earth that was observed by\n",
    "astronomers.  From pictures like the one above, the astronomers deduced how far\n",
    "away each supernova was from Earth and how fast it was moving away from Earth.\n",
    "Their deductions were good, but not perfect.\n",
    "\n",
    "Run the cell below to load the data into a table called `close_novas` and make a scatter plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "34a56fc9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Just run this cell.\n",
    "close_novas = pd.read_csv(\"close_novas.csv\")\n",
    "close_novas.plot.scatter(0, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4cd58c15",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show the first five rows.\n",
    "close_novas.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fab2ac47",
   "metadata": {},
   "source": [
    "### Question 4\n",
    "\n",
    "Looking this plot, make a guess at the age of the universe.\n",
    "\n",
    "**Note**: Make sure you get the units right!  In case you need to know what a\n",
    "parsec is, it's a big unit of distance, equivalent to 30.86 trillion\n",
    "kilometers.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8c4bcaf9",
   "metadata": {},
   "outputs": [],
   "source": [
    "#- Fill this in manually by examining the line above.\n",
    "first_guess_universe_age_years = ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f501ae0e",
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = ok.grade('q4')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "187fcceb",
   "metadata": {},
   "source": [
    "### Fitting the line yourself\n",
    "\n",
    "We would like to be able to calculate the best slope rather than estimating by\n",
    "eye.  Recall that the least-squares regression line for our supernova data is:\n",
    "\n",
    "* the line\n",
    "* with the smallest average (over all the supernovae we observe)\n",
    "* error,\n",
    "* squared,\n",
    "* where the error is\n",
    "\n",
    "$$\\text{the supernova's actual distance from Earth} - \\text{the height of the\n",
    "line at that supernova's speed.}$$\n",
    "\n",
    "We will assume the line goes through the origin, meaning that the intercept of\n",
    "the line is 0.  This corresponds to the assumption that a speed of 0 means 0\n",
    "distance."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d1877d1b",
   "metadata": {},
   "source": [
    "### Question 5\n",
    "\n",
    "Define a function called `errors`.  It should take two arguments:\n",
    "\n",
    "1.  a table like `close_novas` (with the same column names and meanings, but\n",
    "    not necessarily the same data)\n",
    "2.  the slope of a line (a number)\n",
    "\n",
    "It should return an array of the errors made when a line with that slope is\n",
    "used to predict distance from speed for each supernova in the given table.\n",
    "(The error is the actual distance minus the predicted distance.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "069785b6",
   "metadata": {
    "for_assignment_type": "#R_CODE#student"
   },
   "outputs": [],
   "source": [
    "def errors(tbl, slope):\n",
    "    ...\n",
    "    return ..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c310c82a",
   "metadata": {},
   "source": [
    "### Question 6\n",
    "\n",
    "Using `errors`, compute the errors for the line with slope `16000` on the\n",
    "`close_novas` dataset.  Name that array `example_errors`. Then make a scatter\n",
    "plot of the errors.\n",
    "\n",
    "**Hint:** To make a scatter plot of the errors, plot the error for each\n",
    "supernova in the dataset.  Put the actual speed on the horizontal axis and the\n",
    "error on the vertical axis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7670e5f1",
   "metadata": {},
   "outputs": [],
   "source": [
    "example_errors = ...\n",
    "..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5247736a",
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = ok.grade('q6')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b6b39ad",
   "metadata": {},
   "source": [
    "You should find that the errors are almost all negative.  That means our line\n",
    "is a little bit too steep.  Let's find a better one."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1c6284e7",
   "metadata": {},
   "source": [
    "### Question 7\n",
    "\n",
    "We want to use `errors` to find the slope of the line that minimizes the\n",
    "squares of the errors.\n",
    "\n",
    "Define a new function called `sos_errors`.  It should take one argument; the\n",
    "slope of a line (a number).\n",
    "\n",
    "It should return the sum of squared errors made when a line with that slope is\n",
    "used to predict distance from speed for each supernova in `close_novas`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d558e56a",
   "metadata": {
    "for_assignment_type": "#R_CODE#student"
   },
   "outputs": [],
   "source": [
    "def sos_errors(slope):\n",
    "    ...\n",
    "    return ...\n",
    "\n",
    "# Test run of your function.\n",
    "sos_errors(16000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e6339259",
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = ok.grade('q7')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21e3dc81",
   "metadata": {},
   "source": [
    "### Question 8\n",
    "\n",
    "Use your `sos_errors` function to calculate the estimated slope of the\n",
    "least-squares regression line predicting distance from speed for the table\n",
    "`close_novas`.\n",
    "\n",
    "**Note**: If you haven't tried to use the `minimize` function yet, now is a\n",
    "great time to practice. Here's an [example from the textbook](https://lisds.github.io/textbook/mean-slopes/mean_and_slopes).  Or you can use the brute-force try-everything\n",
    "approach in the same page.\n",
    "\n",
    "When you have calculated the slope, set `new_errors` equal to the errors that\n",
    "we get calling `errors` with our new line. The following line will graph the\n",
    "corresponding residual plot with a best fit line.\n",
    "\n",
    "Make sure that the residual plot makes sense.\n",
    "\n",
    "*Hint*: what qualities should the best fit line of a residual plot have?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "06608edb",
   "metadata": {},
   "outputs": [],
   "source": [
    "best_fit_slope = ...\n",
    "\n",
    "new_errors = ...\n",
    "\n",
    "# This code displays the residual plot, given your value for the\n",
    "# best_fit_slope.\n",
    "error_df = pd.DataFrame()\n",
    "error_df[\"Speed (parsecs/year)\"] = close_novas['Speed (parsecs/year)'].copy()\n",
    "error_df[\"Distance errors (million parsecs)\"] = new_errors\n",
    "error_df.plot.scatter(0, 1)\n",
    "\n",
    "# This just shows your answer as a nice string, in billions of years.\n",
    "bf_slope = float(best_fit_slope)  # Make sure slope is a simple number.\n",
    "f\"Slope: {bf_slope:g} (corresponding to an estimated age of {bf_slope/1000:.4f} billion years)\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bbab73f7",
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = ok.grade('q8')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a54ee645",
   "metadata": {},
   "source": [
    "That slope (multiplied by 1 million) is an estimate of the age of the\n",
    "universe. Nice work, data astronomer! You can compare your result to the\n",
    "[Planck project 2015 results](https://arxiv.org/pdf/1502.01589.pdf), which\n",
    "estimated the age of the universe to be 13.799±0.021 billion years.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6e79ac3d",
   "metadata": {},
   "source": [
    "## Done\n",
    "\n",
    "You're finished with the assignment!  Be sure to...\n",
    "\n",
    "- **run all the tests** (the next cell has a shortcut for that),\n",
    "- **Save and Checkpoint** from the \"File\" menu.\n",
    "- Finally, **restart** the kernel for this notebook, and **run all the cells**,\n",
    "  to check that the notebook still works without errors.  Use the\n",
    "  \"Kernel\" menu, and choose \"Restart and run all\".  If you find any\n",
    "  problems, go back and fix them, save the notebook, and restart / run\n",
    "  all again, before submitting.  When you do this, you make sure that\n",
    "  we, your humble markers, will be able to mark your notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "26955163",
   "metadata": {},
   "outputs": [],
   "source": [
    "# For your convenience, you can run this cell to run all the tests at once!\n",
    "import os\n",
    "_ = [ok.grade(q[:-3]) for q in os.listdir(\"tests\") if q.startswith('q')]"
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "cell_metadata_json": true,
   "split_at_heading": true,
   "text_representation": {
    "extension": ".Rmd",
    "format_name": "rmarkdown",
    "format_version": "1.1",
    "jupytext_version": "1.2.4"
   }
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}