Tutorials/PythonVisualize.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Visualizing data in Python\n",
    "**by [Richard W. Evans](https://sites.google.com/site/rickecon/), June 2018**\n",
    "\n",
    "A picture is worth a thousand words. In fact, most journal articles can be summarized by one good figure or one good table. This conjecture makes for a good rule for any researcher. \"What is the best picture or table that can summarize the story I want to tell with my results.\" This section will focus on using data visualization as preliminary a way for understanding and interpreting your data. This chapter is about using visualization for exploring your data. But the principles can be applied to the later steps of presenting your results. \n",
    "\n",
    "The standard visualizations of any data are the scatterplot, the line plot, the bar graph, the histogram (a particular type of bar graph), There are 3-D versions of all of these. Then there are ways to allow each of these methods to display more than three dimensions of the data. Finally, recent developments in D3 (data-driven documents) JavaScript visualizations for web browsers have made interactive graphics available to the common internet user via flyover display, toggle switches, slider bars, and other methods.\n",
    "\n",
    "We will focus on Python's `matplotlib` visualization package. But many other libraries exist in Python (e.g.scikitlearn, chaco, mayavi, bokeh), and many other programming platforms also have powerful visualization packages (e.g., R, Stata, SAS, MATLAB, Mathematica)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 0. Some matplotlib preliminaries\n",
    "Plots in matplotlib are stored as `Figure` objects. The command to create a new figure object is,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "# This next command is specifically for Jupyter Notebook\n",
    "%matplotlib notebook\n",
    "\n",
    "fig = plt.figure()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can reference and change options in the active figure using the `plt.gcf()` function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 0.1. Subplot commands\n",
    "When you declare a figure using the `plt.figure()` command, you must then declare how many subplots will be in the figure. This is done using the `add_subplot()` command. You have to create one or more subplots."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax1 = fig.add_subplot(2, 2, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax2 = fig.add_subplot(2, 2, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax3 = fig.add_subplot(2, 2, 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's fill in the third subplot with a two-way line plot of some white noise. If we don't specify the axes of the specific subplot that we want, matplotlib will use the most recently created subplot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "plt.plot(np.random.randn(50).cumsum(), 'k--')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `k--` is a style option instructing matplotlib to plot the data with a black dashed line. The objects returned by `fig.add_subplot()` avobe are `AxesSubPlot` objects on which you can directly plot on the other empty subplots by calling each one's instance methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = ax1.hist(np.random.randn(100), bins=20, color='y', alpha=0.3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax2.scatter(np.arange(30), np.arange(30) + 3*np.random.randn(30))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since creating a figure with multiple subplots according to a particular layout is such a common task, there is a convenience method, `plt.subplots()`, that creates a new figure and returns a NumPy array containing the created subplot objects. This is very useful as the `axes` array can be easily indexed like a two-dimensional array--for example, `axes[0, 1]`. You can also indicate that subplots should have the same $x$ and $y$ axis using `sharex` and `sharey` options, respectively. This is especially useful when comparing data on the same scale; otherwise, matplotlib auto-scales plot limits independently"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(2, 3)\n",
    "axes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 0.2. Colors, markers, and line styles\n",
    "Matplotlib can accept a single string argument that can encapsulate a subset of the options available in the `color`, `marker`, and `linestyle` options. For example, the following two commands do the same thing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xx = np.arange(45)\n",
    "yy = 0.5 * xx ** 2 - 3 * xx * np.random.randn(45)\n",
    "\n",
    "plt.plot(xx, yy, 'go--')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(xx, yy, linestyle='--', color='g', marker='o')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 0.3. Ticks, labels, and legends\n",
    "The following commands show how to add a title to a plot, add axis labels, and set the minimum and maximum to be plotted on both axes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(xx, yy, linestyle='--', color='g', marker='o')\n",
    "plt.title('This is the best plot ever', fontsize=20)\n",
    "plt.xlabel(r'$X$')\n",
    "plt.ylabel(r'$Y$')\n",
    "plt.xlim([10, 30])  # This gives the xmin and xmax to be plotted\n",
    "plt.ylim([200, 700])  # This gives the ymin and ymax to be plotted"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Adding a legend to a plot is essential if you are plotting more than one series. In the example below, I plot two time series. Both series are with periods 1 and 50. Note that I put grid lines on the plot by using the `MultipleLocator` function from the `matplotlib.ticker` package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib.ticker import MultipleLocator\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "my_age = np.arange(3, 45)\n",
    "knees_hurt = 0.03 * np.arange(50) ** 2\n",
    "all_ages = np.arange(50)\n",
    "plt.plot(all_ages[4:46], my_age, marker='D', label='My age')\n",
    "plt.plot(all_ages, knees_hurt, marker='o',\n",
    "         label='How much my knees hurt')\n",
    "# for the minor ticks, use no labels; default NullFormatter\n",
    "minorLocator = MultipleLocator(1)\n",
    "ax.xaxis.set_minor_locator(minorLocator)\n",
    "plt.grid(b=True, which='major', color='0.65', linestyle='-')\n",
    "plt.title('My knee pain and my age', fontsize=20)\n",
    "plt.xlabel(r'Age')\n",
    "plt.ylabel(r'Pain units')\n",
    "plt.xlim([-1.0, 50])\n",
    "plt.ylim([-2.0, 80])\n",
    "plt.legend(loc='upper left')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now suppose I wanted to change the tick labels on the x axis to highlight some key years (8, 28, and 41) in addition to the 0, 10, 20,... currently listed. And I want those three extra years to be labeled with the strings ('bt', 'mt', 'cr'). And assume I wanted the labels to be vertically oriented (rather than the current horizontal orientation)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xtickvals = [0, 8, 10, 20, 28, 30, 40, 41, 50]\n",
    "xticklabs = [0, 'bt', 10, 20, 'mt', 30, 40, 'cr', 50]\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "my_age = np.arange(3, 45)\n",
    "knees_hurt = 0.03 * np.arange(50) ** 2\n",
    "all_ages = np.arange(50)\n",
    "plt.plot(all_ages[4:46], my_age, marker='D', label='My age')\n",
    "plt.plot(all_ages, knees_hurt, marker='o',\n",
    "         label='How much my knees hurt')\n",
    "# for the minor ticks, use no labels; default NullFormatter\n",
    "minorLocator = MultipleLocator(1)\n",
    "ax.xaxis.set_minor_locator(minorLocator)\n",
    "plt.grid(b=True, which='major', color='0.65', linestyle='-')\n",
    "plt.title('My knee pain and my age', fontsize=20)\n",
    "plt.xlabel(r'Age')\n",
    "plt.ylabel(r'Pain units')\n",
    "plt.xticks(xtickvals, xticklabs, rotation='vertical')\n",
    "plt.xlim([-1.0, 50])\n",
    "plt.ylim([-2.0, 80])\n",
    "plt.legend(loc='upper left')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 0.4. Annotations and drawing on a Subplot\n",
    "Annotations and text can be added using the `text`, `arrow`, and `annotate` functions. `text` draws text at given coordinates `(x, y)` on the plot with optional custom styling.\n",
    "```python\n",
    "ax.text(x, y, 'Hello world!', family='monospace', fontsize=10)\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "from datetime import datetime\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "data = pd.read_csv('DataFiles/spx.csv', index_col=0, parse_dates=True)\n",
    "spx = data['SPX']\n",
    "spx.plot(ax=ax, style='k-')\n",
    "\n",
    "crisis_data = [\n",
    "    (datetime(2007, 10, 11), 'Peak of bull market'),\n",
    "    (datetime(2008, 3, 12), 'Bear Stearns Fails'),\n",
    "    (datetime(2008, 9, 15), 'Lehman Bankruptcy')\n",
    "]\n",
    "\n",
    "for date, label in crisis_data:\n",
    "    ax.annotate(label, xy=(date, spx.asof(date) + 50),\n",
    "                xytext=(date, spx.asof(date) + 200),\n",
    "                arrowprops=dict(facecolor='black'),\n",
    "                horizontalalignment='left', verticalalignment='top')\n",
    "\n",
    "# Zoom in on 2007-2010\n",
    "ax.set_xlim(['1/1/2007', '1/1/2011'])\n",
    "ax.set_ylim([600, 1800])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 0.4. Saving figures\n",
    "Good housekeeping in your file directories is an essential skill for any significant coding project. The directory of the main Python scripts can get cluttered with datafiles, scripts, modules, figure files, and other output files. It is good practice to save files to an \"images\" folder. I prefer to have that folder in the same directory as my calling Python scripts or modules. But in the following format, one can put the images folder anywhere to which the user's terminal can write. This code won't work in this Jupyter notebook.\n",
    "```python\n",
    "from matplotlib.ticker import MultipleLocator\n",
    "import os\n",
    "\n",
    "plot = True\n",
    "\n",
    "if plot:\n",
    "    '''\n",
    "    --------------------------------------------------------------------\n",
    "    cur_path    = string, path name of current directory\n",
    "    output_fldr = string, folder in current path to save files\n",
    "    output_dir  = string, total path of images folder\n",
    "    output_path = string, path of file name of figure to be saved\n",
    "    xx          = (45,) vector, values of xx\n",
    "    yy          = (45,) vector, values of yy\n",
    "    --------------------------------------------------------------------\n",
    "    '''\n",
    "    # Create directory if images directory does not already exist\n",
    "    cur_path = os.path.split(os.path.abspath(__file__))[0]\n",
    "    output_fldr = 'images'\n",
    "    output_dir = os.path.join(cur_path, output_fldr)\n",
    "    if not os.access(output_dir, os.F_OK):\n",
    "        os.makedirs(output_dir)\n",
    "\n",
    "    # Plot line plot of xx and yy\n",
    "    xx = np.arange(45)\n",
    "    yy = 0.5 * xx ** 2 - 3 * xx * np.random.randn(45)\n",
    "    plt.plot(xx, yy, 'go--')\n",
    "    # for the minor ticks, use no labels; default NullFormatter\n",
    "    minorLocator = MultipleLocator(1)\n",
    "    ax.xaxis.set_minor_locator(minorLocator)\n",
    "    plt.grid(b=True, which='major', color='0.65', linestyle='-')\n",
    "    plt.title('Line plot of X and Y', fontsize=20)\n",
    "    plt.xlabel(r'$X$')\n",
    "    plt.ylabel(r'$Y$')\n",
    "    output_path = os.path.join(output_dir, 'xyplot')\n",
    "    plt.savefig(output_path)\n",
    "    # plt.show()\n",
    "    plt.close()\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This code won't work in this Jupyter notebook because of the line in which `cur_path` is set. This script is not being run from a traditional script and therefore cannot use that `os.path.split()` function. But lets go through some of the main principles in this code.\n",
    "* It is often wise to set a Boolean in your script (`plot` in this example) to allow yourself to easily turn the plotting on and off within your code. Sometimes your plotting is just a diagnostic that you don't need to run every time with your analysis.\n",
    "* Have your script intelligently choose the folder from which the plotting functionality is being called. This is the `cur_path = os.path.split()` line. Those lines check what the current directory is, check if a subdirectory named \"images/\" already exists, and makes that directory if it does not already exist. Then the `plt.savefig()` command saves the figure to that folder.\n",
    "* At the end of a script that creates a figure, make sure to include the plt.close() command. If you don't do that, the matplotlib objects can pile up in your computer's memory--regardless of whether the plots are open and visible on your screen--and crash your computer."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. The Almighty Histogram and its Friends"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.1. The cumulative distribution function\n",
    "All the standard univariate descriptive statistics (e.g., mean, standard deviation, quantiles, autocorrelation) are simply statistics derived from the cumulative distribution function (*cdf*). The CDF $F$ of a random variable $X$ is defined as fraction of the values of $X$ that are greater than some realization of the random variable $x$.\n",
    "$$F(x) \\equiv Pr(X\\leq x)$$\n",
    "The estimator for the CDF is the empirical cumulative distribution function (*ecdf*) defined as:\n",
    "$$F_n = \\frac{\\#\\{x_i\\leq x\\}}{n} = \\frac{\\#x_i\\in(-\\infty, x]}{n} = \\frac{1}{n}\\sum_{i=1}^n I_{(-\\infty, x]}(x_i)$$,\n",
    "where $\\{x_1,x_2,...x_n\\}$ is a random sample of $n$ draws from $F$ and ,\n",
    "$$I_A(x) = \\begin{cases}1\\quad\\text{if}\\quad x\\in A \\\\ 0 \\quad\\text{if}\\quad x\\notin A.\\end{cases}$$\n",
    "\n",
    "Because the data for the *ecdf* are discrete, the *ecdf* is a stepwise linear (staircase) function. The Figure below shows the *ecdf* of Old Faithful eruption durations lasting longer than 3 seconds. The smooth red dotted line is a fitted continuous *cdf*.\n",
    "![title](Images/OldFaithfulCDF.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2. The univariate probability density function and the histogram\n",
    "The most visually instructive transformation of the *cdf* is the probability density function (*pdf*). The value of the *pdf* at a given point does not have an intuitive interpretation. The *pdf* is the derivative or slope or rate of change of the *cdf*. But the *pdf* is good for describing what percent of the distribution lies between two values. The *pdf* $f$ is a function of particular values $x$ from the support of the random variable $X$ defined in the following way,\n",
    "$$f(x)\\equiv \\frac{\\partial F(x)}{\\partial x}.$$\n",
    "\n",
    "Scott (2015, p. 39) states that \"the empirical probability density function *epdf* is always a discrete Uniform density over the data, that is, the probability mass $n^{-1}$ at each data point.\" This means that the *epdf* is not continuous, nor is it continuously differentiable. It is simply a one-dimensional scatter diagram with each marker lining up along $n^{-1}$ on the $y$ axis and at each data point on the $x$-axis.\n",
    "\n",
    "A much more instructive representation of the empirical content of a probability density function is the *histogram*. It is more analogous to the discrete *ecdf* from the previous section. A histogram is an empirical estimator for the probability density function. A histogram $h(x)$ for a random variable $X$ is a set of nonoverlapping, contiguous intervals that include the entire support of $X$ that gives the count or frequency of observations in each bin. As such, a histogram $h(x)$ with $B$ bins is completely characterized by bin cutoffs $\\{c_1,c_2,...c_B, c_{B+1}\\}$. Let $bin_b$ be characterized by $x\\in(c_b, c_{b+1}]$ The specification of the count version of the histogram is the following.\n",
    "$$h(x) = \\sum_{i=1}^N I_{bin_b}(x_i)$$\n",
    "The bins need not have equal width, but that is usually the default for most histogram programming and graphical functions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the penny thickness dataset `penny.csv`. Read it in to pandas. Reformat it to a Series with 90 elements such that each observation is a penny thickness. Create a histogram with 10 equally spaced bins."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from pandas import DataFrame, Series\n",
    "\n",
    "pennies = pd.read_csv('DataFiles/penny.csv', names =['year', 'penny1', 'penny2'],\n",
    "                      index_col='year', skiprows=[0, 1, 2, 3, 4, 5])\n",
    "pennies_all = pennies.stack()\n",
    "pennies_all"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(pennies_all, 10, edgecolor='k')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that `plt.hist()` outputs three objects. The first is a list of counts in each bin (`n`). The second is a list of bin endpoints (`bin_cuts`), and the third is a list of patch objects. I have no idea what the patch objects are, but the `n` and `bin_cuts` are straightforward."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now try a histogram with four equally spaced bins. Let's make the bin color green."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n, bin_cuts, patches = plt.hist(pennies_all, 4, facecolor='green', edgecolor='k')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now try a histogram with 90 equally spaced bins. Note that there are only 90 observations in the `pennies_all` data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n, bin_cuts, patches = plt.hist(pennies_all, 90, edgecolor='k')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The cardinal sin of visualizations is not labeling your plot and your axes. We have so far committed that cardinal sin. We can label the plot and the axes by using the following additional matplotlib options."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n, bin_cuts, patches = plt.hist(pennies_all, 10, edgecolor='k')\n",
    "plt.title('Histogram of penny thickness (mils): 1945-1989', fontsize=17)\n",
    "plt.xlabel(r'Thickness (mils)')\n",
    "plt.ylabel(r'Number of observations')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now what if we wanted our histogram to be in frequency form rather than count form. The frequency form of a histogram is more analogous to the probability density function. It is simply a normalized version of the count histogram in which each observation is weighted by $n^{-1}$. We can implement this by creating a weight vector of $N$ elements, each of which is $n^{-1}$, we can then feed this directly into the histogram function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_bins = 10\n",
    "weights = (1 / pennies_all.shape[0]) * np.ones_like(pennies_all)\n",
    "n, bin_cuts, patches = plt.hist(pennies_all, num_bins, weights=weights, edgecolor='k')\n",
    "plt.title('Histogram of penny thickness (mils): 1945-1989', fontsize=17)\n",
    "plt.xlabel(r'Thickness (mils)')\n",
    "plt.ylabel(r'Percent of observations in bin')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can verify that the percents in each bin sum to one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(n.sum())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It would probably be more instructive if the labels on the $x$ axis were closer to the boundaries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n, bin_cuts, patches = plt.hist(pennies_all, num_bins, weights=weights, edgecolor='k')\n",
    "plt.xticks(np.round_(bin_cuts, 1))\n",
    "plt.title('Histogram of penny thickness (mils): 1945-1989', fontsize=17)\n",
    "plt.xlabel(r'Thickness (mils)')\n",
    "plt.ylabel(r'Percent of observations in bin')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is no universal rule for what is the optimal number of bins and what are the optimal bin cutoff values. Some of the principles for making these decisions are based on minimizing error and asymptotic bias between the true probability density function $f(x)$ and the estimator of the *pdf* $h(x)$. A great reference for many optional rules is Scott (2015, chap. 3).\n",
    "\n",
    "One of the most common of these rules is Sturges' number of bins rule for equal width bins. Let $k$ be the number of equal-width bins,\n",
    "$$\\textbf{Sturges' rule}:\\quad k = \\log_2 N$$\n",
    "rounded to the nearest integer. This rule is based on the property that the observation counts in each bin of histogram $h(x)$ are themselvess binomial random variables. Sturges' rule gives a frequency histogram of normal density with mean $(k-1)/2$ and variance $(k-1)/4$. This rule is good for distributions that should be normal. In our penny thickness case, $k$ should equal 8."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_bins = int(np.ceil(1 + np.log2(90)))\n",
    "n, bin_cuts, patches = plt.hist(pennies_all, num_bins, weights=weights, edgecolor='k')\n",
    "plt.xticks(np.round_(bin_cuts, 1))\n",
    "plt.title('Histogram of penny thickness (mils): 1945-1989', fontsize=17)\n",
    "plt.xlabel(r'Thickness (mils)')\n",
    "plt.ylabel(r'Percent of observations in bin')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, from previous versions of this histogram, it is clear that the distribution is not normal. In this case, Sturges' rule oversmoothes the distribution and we don't see the richness between thicknesses 52.0 and 55.0. A histogram with the correct thickness takes a little exploration."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.3. The 3D histogram\n",
    "A 3D histogram can be used as an estimator for a bivariate probability density function. Load in Old Faithful geyser duration data file `geyser.csv`. Make a 3D frequency histogram with geyser durations in the current period `duration(t)` on the $x$ axis and geyser durations in the next period `duration(t+1)` on the $y$ axis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geyser_dur = pd.read_csv('DataFiles/geyser.csv', names =['duration'])\n",
    "geyser_dur"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dur_t = geyser_dur['duration'][:-1]\n",
    "dur_tp1 = geyser_dur['duration'][1:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To make a 3D histogram using matplotlib, you have to use the NumPy function `histogram2d()` and then port the objects into a 3D bar chart using matplotlib's `bar3D()` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "\n",
    "'''\n",
    "--------------------------------------------------------------------\n",
    "bin_num  = integer > 2, number of bins along each axis\n",
    "hist     = (bin_num, bin_num) matrix, bin percentages\n",
    "xedges   = (bin_num+1,) vector, bin edge values in x-dimension\n",
    "yedges   = (bin_num+1,) vector, bin edge values in y-dimension\n",
    "x_midp   = (bin_num,) vector, midpoints of bins in x-dimension\n",
    "y_midp   = (bin_num,) vector, midpoints of bins in y-dimension\n",
    "elements = integer, total number of 3D histogram bins\n",
    "xpos     = (bin_num * bin_num) vector, x-coordinates of each bin\n",
    "ypos     = (bin_num * bin_num) vector, y-coordinates of each bin\n",
    "zpos     = (bin_num * bin_num) vector, zeros or z-coordinates of\n",
    "            origin of each bin\n",
    "dx       = (bin_num,) vector, x-width of each bin\n",
    "dy       = (bin_num,) vector, y-width of each bin\n",
    "dz       = (bin_num * bin_num) vector, height of each bin\n",
    "--------------------------------------------------------------------\n",
    "'''\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111, projection ='3d')\n",
    "bin_num = int(9)\n",
    "hist, xedges, yedges = np.histogram2d(dur_t, dur_tp1, bins=bin_num)\n",
    "hist = hist / hist.sum()\n",
    "x_midp = xedges[:-1] + 0.5 * (xedges[1] - xedges[0])\n",
    "y_midp = yedges[:-1] + 0.5 * (yedges[1] - yedges[0])\n",
    "elements = (len(xedges) - 1) * (len(yedges) - 1)\n",
    "ypos, xpos = np.meshgrid(y_midp, x_midp)\n",
    "xpos = xpos.flatten()\n",
    "ypos = ypos.flatten()\n",
    "zpos = np.zeros(elements)\n",
    "dx = (xedges[1] - xedges[0]) * np.ones_like(bin_num)\n",
    "dy = (yedges[1] - yedges[0]) * np.ones_like(bin_num)\n",
    "dz = hist.flatten()\n",
    "ax.bar3d(xpos, ypos, zpos, dx, dy, dz, color='g', edgecolor='k', zsort='average')\n",
    "ax.set_xlabel('duration(t) (mins.)')\n",
    "ax.set_ylabel('duration(t+1) (mins.)')\n",
    "ax.set_zlabel('Percent of obs.')\n",
    "plt.title('Histogram by Old Faithful duration(t) and duration(t-1)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.4. The countour plot\n",
    "A contour plot is a 2D version of the 3D histogram. It is like a topographical map. The contour plot fits smooth level curves to the density data. Contour plots can also include color scale for monotonically increasing levels of the variables. [TODO: Add material here.]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. The scatter plot\n",
    "The scatter plot is usually first in line among visualization tools for 2D data. Let the $i$th observation of the two dimensional data be represented by $(x_i, y_i)$. A scatter plot places markers at each point $(x_i, y_i)$ in a two-dimensional plot."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The famous Fisher-Anderson Iris data is a canonical dataset of three species of Iris with five variables on each of the observations (sepal length, sepal width, petal length, petal width, species). We can creat a scatter plot of two of these variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "from pandas import Series, DataFrame\n",
    "\n",
    "iris_data = pd.read_csv(\n",
    "    'DataFiles/iris.csv',\n",
    "    names=['length_sep', 'width_sep', 'length_pet', 'width_pet', 'species'], skiprows=[0])\n",
    "length_sep = iris_data['length_sep']\n",
    "width_sep = iris_data['width_sep']\n",
    "fig, axes = plt.subplots(1, 1)\n",
    "plt.scatter(length_sep, width_sep)\n",
    "# plt.xticks(np.round_(bin_cuts, 1))\n",
    "plt.title('Sepal length vs. sepal width', fontsize=17)\n",
    "plt.xlabel(r'sepal length (mm)')\n",
    "plt.ylabel(r'sepal width (mm)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We could color the markers associated with different species with different colors. Note below that, with the colors, we are plotting three dimensions of data (sepal length, width, and species)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "length_sep_set = iris_data[iris_data['species']=='setosa']['length_sep']\n",
    "width_sep_set = iris_data[iris_data['species']=='setosa']['width_sep']\n",
    "length_sep_ver = iris_data[iris_data['species']=='versicolor']['length_sep']\n",
    "width_sep_ver = iris_data[iris_data['species']=='versicolor']['width_sep']\n",
    "length_sep_vrg = iris_data[iris_data['species']=='virginica']['length_sep']\n",
    "width_sep_vrg = iris_data[iris_data['species']=='virginica']['width_sep']\n",
    "\n",
    "fig, axes = plt.subplots(1, 1)\n",
    "plt.scatter(length_sep_set, width_sep_set, color='k', label='setosa')\n",
    "plt.scatter(length_sep_ver, width_sep_ver, color='r', label='versicolor')\n",
    "plt.scatter(length_sep_vrg, width_sep_vrg, color='y', label='virginica')\n",
    "\n",
    "# plt.xticks(np.round_(bin_cuts, 1))\n",
    "plt.title('Sepal length vs. sepal width by species', fontsize=15)\n",
    "plt.xlabel(r'sepal length (mm)')\n",
    "plt.ylabel(r'sepal width (mm)')\n",
    "plt.legend(loc='upper right')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We could plot four dimensions by doing a three dimensional scatter plot (sepal length, sepal width, and petal length) with colors for each species."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "\n",
    "length_sep_set = iris_data[iris_data['species']=='setosa']['length_sep']\n",
    "width_sep_set = iris_data[iris_data['species']=='setosa']['width_sep']\n",
    "length_pet_set = iris_data[iris_data['species']=='setosa']['length_pet']\n",
    "length_sep_ver = iris_data[iris_data['species']=='versicolor']['length_sep']\n",
    "width_sep_ver = iris_data[iris_data['species']=='versicolor']['width_sep']\n",
    "length_pet_ver = iris_data[iris_data['species']=='versicolor']['length_pet']\n",
    "length_sep_vrg = iris_data[iris_data['species']=='virginica']['length_sep']\n",
    "width_sep_vrg = iris_data[iris_data['species']=='virginica']['width_sep']\n",
    "length_pet_vrg = iris_data[iris_data['species']=='virginica']['length_pet']\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111, projection ='3d')\n",
    "ax.scatter(length_sep_set, width_sep_set, length_pet_set,\n",
    "           color='k', label='setosa')\n",
    "ax.scatter(length_sep_ver, width_sep_ver, length_pet_ver,\n",
    "           color='r', label='versicolor')\n",
    "ax.scatter(length_sep_vrg, width_sep_vrg, length_pet_vrg,\n",
    "           color='y', label='virginica')\n",
    "\n",
    "# plt.xticks(np.round_(bin_cuts, 1))\n",
    "plt.title('Sepal length, sepal width, and petal length by species', fontsize=15)\n",
    "ax.set_xlabel(r'sepal length (mm)')\n",
    "ax.set_ylabel(r'sepal width (mm)')\n",
    "ax.set_zlabel(r'petal length (mm)')\n",
    "plt.legend(loc=6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One way to visualize all five dimensions of the Iris data is to use pairwise scatter plots. This method shows the scatterplot of each possible pair of variables. Pandas has a nice prepackaged graphical function that does this. Note in the figure below that there are only  6 unique scatter plots (either the lower triangle or the upper triangle)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "axes = pd.tools.plotting.scatter_matrix(iris_data, alpha=0.2)\n",
    "plt.tight_layout()  # adjust the spacing between the plots"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we wanted to customize our figure, we would simply use matplotlib's `add_subplot()` command. This matrix of scatterplots allows us, in a restricted way, to look at more than three dimensions at a time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "\n",
    "width_pet_set = iris_data[iris_data['species']=='setosa']['width_pet']\n",
    "width_pet_ver = iris_data[iris_data['species']=='versicolor']['width_pet']\n",
    "width_pet_vrg = iris_data[iris_data['species']=='virginica']['width_pet']\n",
    "\n",
    "ax1 = fig.add_subplot(3, 3, 1)\n",
    "ax4 = fig.add_subplot(3, 3, 4)\n",
    "ax5 = fig.add_subplot(3, 3, 5)\n",
    "ax7 = fig.add_subplot(3, 3, 7)\n",
    "ax8 = fig.add_subplot(3, 3, 8)\n",
    "ax9 = fig.add_subplot(3, 3, 9)\n",
    "\n",
    "ax1.scatter(length_sep_set, width_sep_set, color='k', label='setosa')\n",
    "ax1.scatter(length_sep_ver, width_sep_ver, color='r', label='versicolor')\n",
    "ax1.scatter(length_sep_vrg, width_sep_vrg, color='y', label='virginica')\n",
    "ax1.set_xlabel('sepal length (mm)')\n",
    "ax1.set_ylabel('sepal width (mm)')\n",
    "ax1.set_xticks([4, 5, 6, 7, 8])\n",
    "\n",
    "ax4.scatter(length_sep_set, length_pet_set, color='k', label='setosa')\n",
    "ax4.scatter(length_sep_ver, length_pet_ver, color='r', label='versicolor')\n",
    "ax4.scatter(length_sep_vrg, length_pet_vrg, color='y', label='virginica')\n",
    "ax4.set_xlabel('sepal length (mm)')\n",
    "ax4.set_ylabel('petal length (mm)')\n",
    "ax4.set_xticks([4, 5, 6, 7, 8])\n",
    "\n",
    "ax5.scatter(width_sep_set, length_pet_set, color='k', label='setosa')\n",
    "ax5.scatter(width_sep_ver, length_pet_ver, color='r', label='versicolor')\n",
    "ax5.scatter(width_sep_vrg, length_pet_vrg, color='y', label='virginica')\n",
    "ax5.set_xlabel('sepal width (mm)')\n",
    "ax5.set_ylabel('petal length (mm)')\n",
    "ax5.set_xticks([1, 2, 3, 4, 5])\n",
    "\n",
    "ax7.scatter(length_sep_set, width_pet_set, color='k', label='setosa')\n",
    "ax7.scatter(length_sep_ver, width_pet_ver, color='r', label='versicolor')\n",
    "ax7.scatter(length_sep_vrg, width_pet_vrg, color='y', label='virginica')\n",
    "ax7.set_xlabel('sepal length (mm)')\n",
    "ax7.set_ylabel('petal width (mm)')\n",
    "ax7.set_xticks([4, 5, 6, 7, 8])\n",
    "\n",
    "ax8.scatter(width_sep_set, width_pet_set, color='k', label='setosa')\n",
    "ax8.scatter(width_sep_ver, width_pet_ver, color='r', label='versicolor')\n",
    "ax8.scatter(width_sep_vrg, width_pet_vrg, color='y', label='virginica')\n",
    "ax8.set_xlabel('sepal width (mm)')\n",
    "ax8.set_ylabel('petal width (mm)')\n",
    "ax8.set_xticks([1, 2, 3, 4, 5])\n",
    "\n",
    "ax9.scatter(length_pet_set, width_pet_set, color='k', label='setosa')\n",
    "ax9.scatter(length_pet_ver, width_pet_ver, color='r', label='versicolor')\n",
    "ax9.scatter(length_pet_vrg, width_pet_vrg, color='y', label='virginica')\n",
    "ax9.set_xlabel('petal length (mm)')\n",
    "ax9.set_ylabel('petal width (mm)')\n",
    "\n",
    "plt.tight_layout()  # Adjust the spacing between the plots"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These pairwise scatter plots allow you to start to see how the multiple dimensions (more than three) might work together. The plot above shows 5 dimensions of the data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Linked brushing* or *painting* is the use of a pointing device such as a mouse to highlight certain markers in one plot and have them automatically also be highlighted in the other plots.\n",
    "* A python library that only works with Python 2.7 to 3.4 is the `mpld3` package. An example of a linked brushing pairwise system of scatter plots using the Iris data is [here](https://mpld3.github.io/examples/linked_brush.html).\n",
    "* A page showing Bokeh's capability for linked brushing is [here](http://bokeh.pydata.org/en/0.10.0/docs/user_guide/interaction.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. The line plot\n",
    "The two-way line plot is the staple for visualizing time series data. It is also an alternative to a scatterplot when there is a clear order to the points in the data. The line plot connects the points."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1. Line plot\n",
    "We did some examples of the line plot in the Section 0 on Python preliminaries. The `matplotlib.pyplot` package has a function `plot` for making two-way line plots. And you can stack them on top of each other in the same figure by simply calling one after another. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "my_age = np.arange(3, 45)\n",
    "knees_hurt = 0.03 * np.arange(50) ** 2\n",
    "all_ages = np.arange(50)\n",
    "plt.plot(all_ages[4:46], my_age, marker='D', label='My age')\n",
    "plt.plot(all_ages, knees_hurt, marker='o',\n",
    "         label='How much my knees hurt')\n",
    "# for the minor ticks, use no labels; default NullFormatter\n",
    "minorLocator = MultipleLocator(1)\n",
    "ax.xaxis.set_minor_locator(minorLocator)\n",
    "plt.grid(b=True, which='major', color='0.65', linestyle='-')\n",
    "plt.title('My knee pain and my age', fontsize=20)\n",
    "plt.xlabel(r'Age')\n",
    "plt.ylabel(r'Pain units')\n",
    "plt.xlim([-1.0, 50])\n",
    "plt.ylim([-2.0, 80])\n",
    "plt.legend(loc='upper left')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2. 3D surface plots vs. 2D level plots\n",
    "You can productively plot 3D surface visualizations in matplotlib. If you are working from the terminal on your computer or creating a D3 web manipulatable graphic, a 3D image is a productive summary of data. However, in print, a 3D image is often not ideal. If the image must remain static, level curves are often more instructive."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we did with the 3D scatter plot, we will use matplotlib's `mpl_toolkits.mplot3d` toolkit. Read in as a NumPy array the `emat.csv` file. This 80 x 7 NumPy array contains relative estimates of hourly compensation by age (21 to 100) and income percentile (0-24%, 25-49%, 50-69%, 70-79%, 80-89%, 90-98%, 99-100%)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "emat = np.loadtxt('DataFiles/emat.csv', delimiter=',')\n",
    "emat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib\n",
    "from matplotlib import cm\n",
    "cmap1 = matplotlib.cm.get_cmap('summer')\n",
    "\n",
    "ages = np.arange(21, 101)\n",
    "abil_midp = np.array([12.5, 37.5, 60, 75, 85, 94, 99.5])\n",
    "abil_mesh, age_mesh = np.meshgrid(abil_midp, ages)\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.gca(projection='3d')\n",
    "ax.plot_surface(age_mesh, abil_mesh, emat, rstride=8,\n",
    "                cstride=1, cmap=cmap1)\n",
    "ax.set_xlabel(r'age-$s$')\n",
    "ax.set_ylabel(r'ability type -$j$')\n",
    "ax.set_zlabel(r'ability $e_{j,s}$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we can't manipulate this graphic, then we might want to just plot the level curves."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "abil_pcts = np.array([0.0, 0.25, 0.50, 0.70, 0.80, 0.90, 0.99, 1.0])\n",
    "J = 7\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = plt.subplot(111)\n",
    "linestyles = np.array([\"-\", \"--\", \"-.\", \":\",])\n",
    "markers = np.array([\"x\", \"v\", \"o\", \"d\", \">\", \"|\"])\n",
    "pct_lb = 0\n",
    "for j in range(J):\n",
    "    this_label = (str(int(np.rint(pct_lb))) + \" - \" +\n",
    "        str(int(np.rint(pct_lb + 100*abil_pcts[j]))) + \"%\")\n",
    "    pct_lb += 100*abil_pcts[j]\n",
    "    if j <= 3:\n",
    "        ax.plot(ages, np.log(emat[:, j]), label=this_label,\n",
    "            linestyle=linestyles[j], color='black')\n",
    "    elif j > 3:\n",
    "        ax.plot(ages, np.log(emat[:, j]), label=this_label,\n",
    "            marker=markers[j-4], color='black')\n",
    "ax.axvline(x=80, color='black', linestyle='--')\n",
    "box = ax.get_position()\n",
    "ax.set_position([box.x0, box.y0, box.width * 0.8, box.height])\n",
    "ax.legend(loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "ax.set_xlabel(r'age-$s$')\n",
    "ax.set_ylabel(r'log ability $log(e_{j,s})$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. More than 3 dimensions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.1. Plotting a time component with dynamic graphics\n",
    "[Here](https://sites.google.com/site/evanjmagnusson/) is an example of a gif that a researcher made of the time series of a three-dimensional object. Python's `PIL` package has the functionality to do this (e.g., `from PIL import Image, ImageDraw`). See examples [here](http://stackoverflow.com/questions/13041893/generating-an-animated-gif-in-python) and [here](http://stackoverflow.com/questions/34126016/trying-to-create-a-gif-from-a-set-of-images)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.2. Brushing\n",
    "*Linked brushing* or *painting* is the use of a pointing device such as a mouse to highlight certain markers in one plot and have them automatically also be highlighted in the other plots.\n",
    "* A python library that only works with Python 2.7 to 3.4 is the `mpld3` package. An example of a linked brushing pairwise system of scatter plots using the Iris data is [here](https://mpld3.github.io/examples/linked_brush.html).\n",
    "* A page showing Bokeh's capability for linked brushing is [here](http://bokeh.pydata.org/en/0.10.0/docs/user_guide/interaction.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.3. Using glyphs\n",
    "Glyphs, or images with many characteristics, can be used as scatter plot markers in order to incorporate more than two dimensions of data. You can use marker color, marker length, marker shape, marker size to encapsulate more dimensions of the data. Below is a rich example.\n",
    "![title](Images/glyphexamp.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.4. Chernoff faces\n",
    "The Chernoff face is an interesting way to compare data visually in a way that simultaneously presents as many dimensions of the data as you want. Chernoff (1973) proposed using an image that associates different attributes of the image to variables in the data. For example, the size of the nose might represent the variable \"wind speed\", the color of the skin might represent \"temperature\", the openness of the eyes might represents \"time of day\", and the length of the hair might represent \"age\". The only limit to the dimensionality is how many characteristics you can put on the face. Also note that it is valuable that there be some intuitive association between the variable and the characteristic of the glyph. And, of course, the image could be anything besides a face."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![title](Images/ChernoffFaces.png)\n",
    "The above image represents data for lawyers' ratings of twelve judges. Each facial feature represents a dimension of the ratings. Because the features of the faces vary in perceived importance, the way in which variables are mapped to the features should be carefully chosen (e.g. eye size and eyebrow-slant have been found to carry significant weight)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.5. D3 javascript visualizations\n",
    "D3.js or D3 or Data Driven Documents is a Javascript library for producing dynamic, interactive data visualizations in web browsers. D3 graphics are the future of visualization. These techniques have not yet become mainstream among researchers. The features available in D3 visualizations include slider bars, toggle buttons, flyover info, and dynamic time advancement.\n",
    "\n",
    "The *New York Times* graphics and data department have been the best at producing D3 graphical representations. My favorit example is [this graphic](http://www.nytimes.com/interactive/2015/03/19/upshot/3d-yield-curve-economic-growth.html?_r=0) that shows the user the ins and outs of U.S., German, and Japanese yield curve data (interest rates across different maturities at different points in time).\n",
    "\n",
    "D3 visualizations allow us to use objects like slider bars to envision more than three dimensions of data at a time. The [Bokeh package](http://bokeh.pydata.org/en/latest/) is Python's library for D3 visualizations. Below is a screenshot of an example. This [*Wall Street Journal* article](http://blogs.wsj.com/economics/2016/08/18/why-its-so-hard-to-get-rid-of-tax-loopholes/) has an example that my group at the Open Source Policy Center built using Bokeh.\n",
    "![title](Images/BokehSlide.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "* Chernoff, H., \"The Use of Faces to Represent Points in $k$-dimensional Space Graphically,\" _Journal of the American Statistical Association_, 94, pp. 807-823 (1973).\n",
    "* McKinney, Wes, Python for Data Analysis, O'Reilly Media, Inc. (2013).\n",
    "* [Python labs](http://www.acme.byu.edu/?page_id=2067), Applied and Computational Mathematics Emphasis (ACME), Brigham Young University.\n",
    "* Scott, D.W., Gotto, A.M., Cole, J.S., and Gorry, G.A., \"Plasma Lipids as Collateral Risk Factors in Coronary Artery Disease: A Study of 371 Males with Chest Pain,\" *Journal of Chronic Diseases*, 31, pp. 337-345.\n",
    "* Scott, David W., *Multivariate Density Estimation: Theory, Practice, and Visualization*, 2nd edition, John Wiley & Sons, Inc. (2015)."
   ]
  }
 ],
 "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}