07_visualization.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Visualization\n",
    "\n",
    "\n",
    "The visualization of the results of an analysis represents one of the most delicate and important steps in data analysis, and typically tells a lot about the quality of the analysis itself. A good and bad examples are:\n",
    "   * the discovery of the [Higgs boson](https://twiki.cern.ch/twiki/pub/CMSPublic/Hig13002TWiki/HZZ4l_animated.gif)\n",
    "   * the claim of [superluminal neutrinos](https://arxiv.org/pdf/1109.4897v1.pdf)\n",
    "\n",
    "Graphics (plots) should convey as much information as possible, in an intuitive way (should require little or no explanation), and without sacrificing the clarity. There are plenty of examples of bad practices around...\n",
    "\n",
    "There are essentially two kind of plots you'll ever want to make, histograms and scatter plots."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Histograms\n",
    "\n",
    "Histograms are representations of the distribution of a certain feature of the data. Histograms can be thought as scalar functions, $f(x): {\\rm I\\!R}^n \\to {\\rm I\\!R}$, where at most $n=2$. Histograms are often used to compare different distributions, for instance the data and its expectation (e.g. simulation). \n",
    "\n",
    "Histograms are **binned**, which means that the domain is divided in intervals (\"bins\", or \"buckets\"), and for the $i$th bin $x_i<x\\le x_{i+1}$ and the content of the bin (say $n_i$) is  $n_i=\\frac{\\int_{x_{i}}^{x_{i+1}} f(x)}{x_{i+1}- x_{i}}$.\n",
    "\n",
    "Bins doesn't need to have the same size, i.e. it could well be that $(x_{i+1}- x_{i}) \\neq (x_{j+1}- x_{j})$ for $i\\neq j$.\n",
    "\n",
    "If the bin size is small, $\\lim_{x_{i+1}- x_{i} \\to 0} n_i \\simeq f(\\frac{x_{i+1}+x_{i}}{2})$.\n",
    "\n",
    "\n",
    "\n",
    "Choosing the number of bins and thus the bin size requires a special care. Typically, the content of each $i$th bin, $n_i$, should be statistically significant, i.e. the corresponding Poisson uncertainty, $1/\\sqrt{n_i}$, should be small compared to $n_i$.\n",
    "\n",
    "In general, the uncertainties affecting the data (or the expectation) should be reported as error bars, especially if the distribution is normalized or comes from a simulation (and thus needs to be weighted).\n",
    "\n",
    "**Note:** the Poisson uncertainty is valid only for random counts; in all other cases, the appropriate uncertainty should be derived.\n",
    "\n",
    "#### Shape and Normalization\n",
    "\n",
    "The are essentially two main quantities that are conveyed by an histogram:\n",
    "\n",
    "   * The total amount of data, which is referred to as \"normalization\"\n",
    "   * The \"shape\" or differential distribution, which should follow the underlying p.d.f\n",
    "   \n",
    "Both quantities can be expressed separately, by quoting $N=\\sum_i n_i$ the number of entries and displaying the original bins rescaled by $n_i / N$, such that $\\sum_i n_i = 1$.\n",
    "\n",
    "In fact, the definition of normalization factor could be ambiguous; it could either be $N$ or $I=\\int f(x) dx$. The former is the most commonly used, but the latter is often the default of built-in methods in several visualization packages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import Image\n",
    "Image(\"https://matplotlib.org/stable/_images/sphx_glr_histogram_features_001_2_0x.png\", width=600)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Scatter plots\n",
    "\n",
    "Scatter plots are used to compare two features, in order to visually detect possible correlations between them. \n",
    "\n",
    "Scatter plots can be *unbinned* or *binned*.\n",
    "- Unbinned scatter plots are often prefered until the number of point isn't too large, after which the plot is not readable anymore. The display of different features is also possible by using a different marker or color for the data points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Image(\"https://matplotlib.org/stable/_images/sphx_glr_scatter_demo2_001_2_0x.png\", width=600)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Binned scatter plots are, as a matter of fact, 2D histograms. The bin content $n_i$ is represented through a color code indicating the \"3rd dimension\". Alternatively, this can be shown as a \"contour plot\", i.e. with smoothened edges and color-coded levels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Image(\"https://matplotlib.org/stable/_images/sphx_glr_contourf_log_001_2_0x.png\", width=600)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The condensation of the scatter plots into \"profiles\" is also very useful. The data are binned along the $x$ (if you had to bin on the other variable, just invert the axes), and for every bin the mean and the standard deviation of the corresponding $y$ values are calculated and reported as data points and their uncertainties, respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Image(\"images/profile.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Other kinds of plots\n",
    "\n",
    "There are an almost infinite number of other plots, but in some way or another, they can be tracked back to scatter plots and histograms.\n",
    "\n",
    "A non-exhaustive list can be seen in the [Matplotlib documentation](https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.show.html#matplotlib.pyplot.show).\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Matplotlib (pyplot)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Matplotlib is a comprehensive library for creating static, animated, and interactive visualizations in Python. It provides an intuitive set of tools to visualize the data, and is well integrated with Python libraries e.g. Numpy.\n",
    "With relatively few commands, it allows everybody to produce readable and high-quality plots.\n",
    "Link to the [documentation](https://matplotlib.org/stable/index.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "# always useful\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A plot is a container of nested Python objects. A *Figure* object is the outermost container for a matplotlib graphic object, which can contain multiple *Axes* objects. The name can be confusing, as an Axes actually translates into what we think of as an individual plot or graph (rather than the plural of “axis,” as one might expect).\n",
    "\n",
    "A Figure object is like a container of one or more Axes (a.k.a. plots). Below the Axes in the hierarchy are individual *Axis*, and smaller objects such as tick marks, individual lines, legends, and text boxes. Every element of a chart is a Python object, which can be manipulated if needed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import Image\n",
    "Image(url='https://files.realpython.com/media/fig_map.bc8c7cabd823.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, _ = plt.subplots() # create a Figure\n",
    "print(type(fig))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the example above, we created two variables with `plt.subplots()`. The first is a top-level Figure object. The second is a “throwaway” variable that we don’t need just yet, denoted with an underscore. Using attribute notation, it is easy to go through the figure hierarchy and access the sub-objects, e.g. to see the first tick of the Y axis of the first Axes object:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(fig.axes)\n",
    "one_tick = fig.axes[0].yaxis.get_major_ticks()[0]\n",
    "print(type(one_tick))\n",
    "ticks, labels = plt.xticks()\n",
    "print(ticks)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the example above, `fig` (an instance of the Figure class) has multiple Axes (a list, for which we take the first element). Each Axes has a yaxis and xaxis, each of which have a collection of “major ticks,” and we grab the first one.\n",
    "\n",
    "Matplotlib documentation presents this as a figure anatomy, rather than an explicit hierarchy:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib.ticker import AutoMinorLocator, MultipleLocator, FuncFormatter\n",
    "\n",
    "np.random.seed(19680801)\n",
    "\n",
    "X = np.linspace(0.5, 3.5, 100)\n",
    "Y1 = 3 + np.cos(X)\n",
    "Y2 = 1 + np.cos(1 + X / 0.75) / 2\n",
    "Y3 = np.random.uniform(Y1, Y2, len(X))\n",
    "\n",
    "fig = plt.figure(figsize=(8, 8))\n",
    "ax = fig.add_subplot(1, 1, 1, aspect=1)\n",
    "\n",
    "\n",
    "def minor_tick(x, pos):\n",
    "    if not x % 1.0:\n",
    "        return \"\"\n",
    "    return \"%.2f\" % x\n",
    "\n",
    "ax.xaxis.set_major_locator(MultipleLocator(1.000))\n",
    "ax.xaxis.set_minor_locator(AutoMinorLocator(4))\n",
    "ax.yaxis.set_major_locator(MultipleLocator(1.000))\n",
    "ax.yaxis.set_minor_locator(AutoMinorLocator(4))\n",
    "ax.xaxis.set_minor_formatter(FuncFormatter(minor_tick))\n",
    "\n",
    "ax.set_xlim(0, 4)\n",
    "ax.set_ylim(0, 4)\n",
    "\n",
    "ax.tick_params(which='major', width=1.0)\n",
    "ax.tick_params(which='major', length=10)\n",
    "ax.tick_params(which='minor', width=1.0, labelsize=10)\n",
    "ax.tick_params(which='minor', length=5, labelsize=10, labelcolor='0.25')\n",
    "\n",
    "ax.grid(linestyle=\"--\", linewidth=0.5, color='.25', zorder=-10)\n",
    "\n",
    "ax.plot(X, Y1, c=(0.25, 0.25, 1.00), lw=2, label=\"Blue signal\", zorder=10)\n",
    "ax.plot(X, Y2, c=(1.00, 0.25, 0.25), lw=2, label=\"Red signal\")\n",
    "ax.plot(X, Y3, linewidth=0,\n",
    "        marker='o', markerfacecolor='w', markeredgecolor='k')\n",
    "\n",
    "ax.set_title(\"Anatomy of a figure\", fontsize=20, verticalalignment='bottom')\n",
    "ax.set_xlabel(\"X axis label\")\n",
    "ax.set_ylabel(\"Y axis label\")\n",
    "\n",
    "ax.legend()\n",
    "\n",
    "\n",
    "def circle(x, y, radius=0.15):\n",
    "    from matplotlib.patches import Circle\n",
    "    from matplotlib.patheffects import withStroke\n",
    "    circle = Circle((x, y), radius, clip_on=False, zorder=10, linewidth=1,\n",
    "                    edgecolor='black', facecolor=(0, 0, 0, .0125),\n",
    "                    path_effects=[withStroke(linewidth=5, foreground='w')])\n",
    "    ax.add_artist(circle)\n",
    "\n",
    "\n",
    "def text(x, y, text):\n",
    "    ax.text(x, y, text, backgroundcolor=\"white\",\n",
    "            ha='center', va='top', weight='bold', color='blue')\n",
    "\n",
    "\n",
    "# Minor tick\n",
    "circle(0.50, -0.10)\n",
    "text(0.50, -0.32, \"Minor tick label\")\n",
    "\n",
    "# Major tick\n",
    "circle(-0.03, 4.00)\n",
    "text(0.03, 3.80, \"Major tick\")\n",
    "\n",
    "# Minor tick\n",
    "circle(0.00, 3.50)\n",
    "text(0.00, 3.30, \"Minor tick\")\n",
    "\n",
    "# Major tick label\n",
    "circle(-0.15, 3.00)\n",
    "text(-0.15, 2.80, \"Major tick label\")\n",
    "\n",
    "# X Label\n",
    "circle(1.80, -0.27)\n",
    "text(1.80, -0.45, \"X axis label\")\n",
    "\n",
    "# Y Label\n",
    "circle(-0.27, 1.80)\n",
    "text(-0.27, 1.6, \"Y axis label\")\n",
    "\n",
    "# Title\n",
    "circle(1.60, 4.13)\n",
    "text(1.60, 3.93, \"Title\")\n",
    "\n",
    "# Blue plot\n",
    "circle(1.75, 2.80)\n",
    "text(1.75, 2.60, \"Line\\n(line plot)\")\n",
    "\n",
    "# Red plot\n",
    "circle(1.20, 0.60)\n",
    "text(1.20, 0.40, \"Line\\n(line plot)\")\n",
    "\n",
    "# Scatter plot\n",
    "circle(3.20, 1.75)\n",
    "text(3.20, 1.55, \"Markers\\n(scatter plot)\")\n",
    "\n",
    "# Grid\n",
    "circle(3.00, 3.00)\n",
    "text(3.00, 2.80, \"Grid\")\n",
    "\n",
    "# Legend\n",
    "circle(3.70, 3.80)\n",
    "text(3.70, 3.60, \"Legend\")\n",
    "\n",
    "# Axes\n",
    "circle(0.5, 0.5)\n",
    "text(0.5, 0.3, \"Axes\")\n",
    "\n",
    "# Figure\n",
    "circle(-0.3, 0.65)\n",
    "text(-0.3, 0.45, \"Figure\")\n",
    "\n",
    "color = 'blue'\n",
    "ax.annotate('Spines', xy=(4.0, 0.35), xycoords='data',\n",
    "            xytext=(3.3, 0.5), textcoords='data',\n",
    "            weight='bold', color=color,\n",
    "            arrowprops=dict(arrowstyle='->',\n",
    "                            connectionstyle=\"arc3\",\n",
    "                            color=color))\n",
    "\n",
    "ax.annotate('', xy=(3.15, 0.0), xycoords='data',\n",
    "            xytext=(3.45, 0.45), textcoords='data',\n",
    "            weight='bold', color=color,\n",
    "            arrowprops=dict(arrowstyle='->',\n",
    "                            connectionstyle=\"arc3\",\n",
    "                            color=color))\n",
    "\n",
    "ax.text(4.0, -0.4, \"Made with http://matplotlib.org\",\n",
    "        fontsize=10, ha=\"right\", color='.5')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stateful Versus Stateless Approaches\n",
    "\n",
    "There is an important difference between the stateful (state-based, state-machine) and stateless (object-oriented, OO) interfaces.\n",
    "\n",
    "Above, we used `import matplotlib.pyplot as plt` to import the pyplot module from matplotlib and name it plt.\n",
    "\n",
    "Almost all functions from pyplot, such as `plt.plot()`, are implicitly either referring to an existing current Figure and current Axes, or creating them anew if they do not exist.\n",
    "\n",
    "* The **stateful interface** makes its calls with `plt.plot()` and other top-level pyplot functions. There is only one Figure or Axes that you’re manipulating at a given time (and this is the last one created or called), and you don’t need to explicitly refer to it.\n",
    "\n",
    "* Modifying the underlying objects directly is the **object-oriented** approach. We usually do this by calling methods of an Axes object, which is the object that represents a plot itself.\n",
    "\n",
    "This all boils down to the definition of `plt.plot()` as:\n",
    "``` python\n",
    "# matplotlib/pyplot.py\n",
    "def plot(*args, **kwargs):\n",
    "    \"\"\"An abridged version of plt.plot().\"\"\"\n",
    "    ax = plt.gca()\n",
    "    return ax.plot(*args, **kwargs)\n",
    "\n",
    "def gca(**kwargs):\n",
    "    \"\"\"Get the current Axes of the current Figure.\"\"\"\n",
    "    return plt.gcf().gca(**kwargs)\n",
    "```\n",
    "\n",
    "Both approaches have their own advantages and disadvantages.\n",
    "The steteful approach is convenient for simple plots, while the stateless (object-oriented) approach is typically more customizable and comes in handy as graphs become more complex."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A few simple examples\n",
    "\n",
    "1) Let's start with a **simple histogram**. In the example below:\n",
    "   - create the dataset\n",
    "   - create the fig, axes, etc.\n",
    "   - use the method of `pyplot.plt` to create and display the histogram in the most recent fig (*stateful approach*)\n",
    "   - save the relevant histogram quantities, possibly for futher analysis\n",
    "   - do cosmetic changes to the plot\n",
    "   - `tight_layout()` applies to the Figure object to clean up whitespace padding"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nrnd = np.random.normal(loc=0., scale=1., size=(100,)) # generate 100 random numbers using a Gaussian pdf\n",
    "\n",
    "fig = plt.figure(figsize=(10, 6)) # create the figure and the axes\n",
    "h, bins, _ = plt.hist(nrnd, bins=25, range=(-3, +3), density=True, cumulative=False) # create and plot the histogram\n",
    "# get the axes from the figure and pick the first (and only) element (beacuse there is only one axes)\n",
    "ax = fig.axes[0]\n",
    "ax.set_title('Gaussian distribution')\n",
    "ax.set_xlabel('x')\n",
    "ax.set_ylabel('Entries')\n",
    "fig.tight_layout()\n",
    "\n",
    "# print histogram content and bins\n",
    "print(\"Histogram content:\", h)\n",
    "print(\"Bin boundaries:\", bins)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2) A slightly more complex example: **stacked histograms**. In detail, that's what it is done:\n",
    "\n",
    "   - create three random time series\n",
    "   - define one Figure (fig) containing one Axes (a plot, ax)\n",
    "   - methods of `ax` have been called directly to create a stacked area chart and to add a legend, title, and y-axis label (*stateless approach*). By using the specific `ax` object, there is no ambiguity about which Axes will contain the plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rng = np.arange(50)\n",
    "rnd = np.random.randint(0, 10, size=(3, rng.size)) # create a 3 x 50 array with random int numbers from 0 to 10\n",
    "yrs = 1950 + rng\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(10, 6)) # create the figure and the axes\n",
    "# plt.subplots instead of plt.figure because it allows to get the axes directly\n",
    "ax.stackplot(yrs, rng + rnd, labels=['Eastasia', 'Eurasia', 'Oceania']) # create a stack plot, where the 3 components are summed\n",
    "ax.set_title('Combined debt growth over time')\n",
    "ax.legend(loc='upper left')\n",
    "ax.set_xlabel('Year')\n",
    "ax.set_ylabel('Total debt')\n",
    "ax.set_xlim(left=yrs[0], right=yrs[-1])\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3) Subplots, scatter plots, and **more than one histogram in a plot**\n",
    "\n",
    "   - as a “1x2” Figure is created, the returned result of plt.subplots(1, 2) is now a Figure object and a NumPy array of Axes objects. Alternatively:\n",
    "\n",
    "``` python\n",
    "fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(8, 4))\n",
    "# ax is now a numpy array of axes's\n",
    "```\n",
    "\n",
    "   - we need to do operations on `ax1` and `ax2` individually, which would be difficult to do with the stateful approach. The example below is a good illustration of the object hierarchy, where we are modifying the yaxis belonging to the second Axes, placing its ticks and ticklabels to the right.\n",
    "\n",
    "   - text inside dollar signs utilizes TeX markup (you can try fancier TeX stuff)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.random.randint(low=1, high=11, size=50)\n",
    "y = x + np.random.randint(1, 5, size=x.size)\n",
    "\n",
    "# create one figure with 2 subplots\n",
    "fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(8, 4))\n",
    "# unpack ax1 and ax2 inline\n",
    "\n",
    "ax1.scatter(x=x, y=y, marker='o', c='r', edgecolor='b')\n",
    "ax1.set_title('Scatter: $x$ versus $y$')\n",
    "ax1.set_xlabel('$x$')\n",
    "ax1.set_ylabel('$y$')\n",
    "\n",
    "# either this..\n",
    "data = np.column_stack((x, y)) # join the two numpy arrays\n",
    "ax2.hist(data, bins=np.arange(data.min(), data.max()), label=('x', 'y'))\n",
    "# ..or this\n",
    "#ax2.hist(x, bins=np.arange(data.min(), data.max()), label='x', alpha=0.5)\n",
    "#ax2.hist(y, bins=np.arange(data.min(), data.max()), label='y', alpha=0.5)\n",
    "\n",
    "ax2.legend(loc=(0.65, 0.8))\n",
    "ax2.set_title('Frequencies of $x$ and $y$')\n",
    "ax2.yaxis.tick_right() # move the y axis ticks to the right of the plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4) A more advanced example: **scatter plot with profiles**\n",
    "\n",
    "Let's consider some data from the housing market in California:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from io import BytesIO\n",
    "import tarfile\n",
    "from urllib.request import urlopen\n",
    "\n",
    "# in this example, you can also see how to access a remote tgz compressed file and read it, all within Python\n",
    "\n",
    "url = 'http://www.dcc.fc.up.pt/~ltorgo/Regression/cal_housing.tgz'\n",
    "b = BytesIO(urlopen(url).read())\n",
    "fpath = 'CaliforniaHousing/cal_housing.data'\n",
    "\n",
    "with tarfile.open(mode='r', fileobj=b) as archive:\n",
    "    housing = np.loadtxt(archive.extractfile(fpath), delimiter=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(housing)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's extract the area's average home value, which becomes our independent variable $z$, while `pop` and `age` are the area's population and average house age, respectively:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "z = housing[:, -1]\n",
    "pop, age = housing[:, [4, 7]].T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(pop, '\\n', age, '\\n', z)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's create a fancy layout by means of `plt.subplot2grid`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gridsize = (3, 2)\n",
    "fig = plt.figure(figsize=(12, 8))\n",
    "\n",
    "ax1 = plt.subplot2grid(gridsize, (0, 0), colspan=2, rowspan=2)\n",
    "scatter = ax1.scatter(x=age, y=pop, c=z, cmap='RdYlGn') # create the scatter plot in ax1\n",
    "plt.colorbar(scatter, ax=ax1, format='$%d') # create the color bar on the right (note that we need to specify the ax)\n",
    "ax1.set_title('Home value as a function of home age and area population', fontsize=14)\n",
    "ax1.set_xlabel('Age')\n",
    "ax1.set_ylabel('Area population')\n",
    "ax1.set_yscale('log') # use logarithmic y axis scale\n",
    "\n",
    "# create the two 1D histograms of single variables\n",
    "ax2 = plt.subplot2grid(gridsize, (2, 0))\n",
    "n2, bins2, patches2 = ax2.hist(age, bins='auto')\n",
    "\n",
    "ax3 = plt.subplot2grid(gridsize, (2, 1))\n",
    "n3, bins3, patches3 = ax3.hist(pop, bins='auto', log=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Seaborn\n",
    "\n",
    "Seaborn is a wrapper around Matplotlib which makes it better, with a better graphics and advanced functionalities, especially for what concerns statistical analysis.\n",
    "\n",
    "It allows advanced data and statistical analysis with a very compact synthax, and is very well integrated with Pandas.\n",
    "\n",
    "Unfortunately, Seaborn is not actively developed anymore, and it's quite rigid in the sense that it's very difficult (or impossible) to do things that were not implemented by the original developer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import seaborn as sns\n",
    "#sns.set(style=\"darkgrid\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Scatter plots\n",
    "\n",
    "Let's load some data and start play with it. (Seaborn is nice as it comes with plenty of example datasets you can load as Pandas DataFrames)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tips = sns.load_dataset(\"tips\")\n",
    "tips.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at some example of scatter plot in seaborn. `sns.relplot()` (link to [documentation](https://seaborn.pydata.org/generated/seaborn.relplot.html)) is most commonly used, which refers to `sns.scatterplot()`, also available.\n",
    "\n",
    "Essentialy we want to play with the \"extra dimensions\" it allows to display. Bear in mind that the more dimensions you put in there, the more difficult the plot results to read."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.relplot(x=\"total_bill\", y=\"tip\", hue=\"smoker\", data=tips)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.relplot(x=\"total_bill\", y=\"tip\", hue=\"smoker\", style=\"time\", data=tips)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Not only categorical but also quantitative \"3rd dimension\" is of course possible. With color (hue):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.relplot(x=\"total_bill\", y=\"tip\", hue=\"size\", data=tips)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And with marker size:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.relplot(x=\"total_bill\", y=\"tip\", size=\"size\", sizes=(15, 200), data=tips)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Histograms\n",
    "\n",
    "Seaborn provides advanced functionalities with histograms, not just in terms of graphics and visualization, but in terms of statistical analysis too."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.random.randn(100)\n",
    "sns.histplot(x=x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `displot()` ([documentation](https://seaborn.pydata.org/generated/seaborn.displot.html)) and `histplot()` ([documentation](https://seaborn.pydata.org/generated/seaborn.histplot.html)) methods extend the functionalities of `pyplot.hist()`, for instance can easily provide the kernel density estimate ([KDE](https://en.wikipedia.org/wiki/Kernel_density_estimation)), the individual entries at the bottom of the histogram, and much more."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# with unbinned distribution too\n",
    "sns.displot(x=x, kde=True, rug=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Kernel Density Estimate\n",
    "\n",
    "The **Kernel density estimation** (KDE) is a *non-parametric* technique that allows to estimate the PDF of a distribution (e.g. from a histogram).\n",
    "\n",
    "It consists of replacing each entry of the histogram (each small vertical line above the X axis) with a Gaussian, called a *kernel* with:\n",
    "\n",
    "  a) the mean corresponding to the entry of the histogram\n",
    "  \n",
    "  b) the width (the bandwidth) is not fixed, but it has to be chosen accordingly\n",
    "  \n",
    "  c) the integral of the Gaussian is often chosen to be 1 for entries with unitary weight.\n",
    "\n",
    "The kernels are summed to make the kernel density estimate.\n",
    "\n",
    "Note that the bandwidth of the kernels is a free parameter that exhibits a strong influence on the resulting estimate. There are some common choices, however these are generally rules of thumb."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Scatter plots and projections\n",
    "\n",
    "The `sns.jointplot()` ([documentation](https://seaborn.pydata.org/generated/seaborn.jointplot.html)) allows to visualize at the same time both the scatter plot between two features and show their projections along the axes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "mean, cov = [0, 1], [(1, .5), (.5, 1)]\n",
    "data = np.random.multivariate_normal(mean, cov, 500)\n",
    "df = pd.DataFrame(data, columns=[\"x\", \"y\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.jointplot(x=\"x\", y=\"y\", data=df) # note that you can provide the DataFrame and the names of the columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, ax = plt.subplots(figsize=(6, 6))\n",
    "cmap = sns.cubehelix_palette(as_cmap=True, dark=0, light=1, reverse=False) # define color palette\n",
    "sns.kdeplot(x=df[\"x\"], y=df[\"y\"], cmap=cmap, n_levels=60, shade=True) # alternate method to pass the DataFrame data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g = sns.jointplot(x=\"total_bill\", y=\"tip\", data=tips, kind=\"hex\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "iris = sns.load_dataset(\"iris\") # another well known dataset integrated in Seaborn\n",
    "g = sns.jointplot(x=\"sepal_width\", y=\"petal_length\", data=iris, kind=\"kde\", space=0, color=\"g\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g = sns.jointplot(x=\"total_bill\", y=\"tip\", data=tips, kind=\"reg\") # perform a regression on the data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Chasing correlations\n",
    "\n",
    "Seaborn allows you to look for correlations systematically checking pair-wise your dataset's features:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.pairplot(iris)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g = sns.PairGrid(iris)\n",
    "g.map_diag(sns.kdeplot)\n",
    "g.map_offdiag(sns.kdeplot, n_levels=6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot accessibility (Bokeh, Plotly, etc.)\n",
    "\n",
    "Often, plots need to be accessible by others from remote, and static formats like `.png` and `.pdf` are not the best option anymore. There are packages, like Bokeh, which are useful to publish plots on web servers.\n",
    "\n",
    "Other tools like [Dash](https://plot.ly/products/dash/) are extremely powerful and worth checking out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from bokeh.plotting import figure, output_file, show\n",
    "\n",
    "# prepare some data\n",
    "x = [0.1, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0]\n",
    "y0 = [i**2 for i in x]\n",
    "y1 = [10**i for i in x]\n",
    "y2 = [10**(i**2) for i in x]\n",
    "\n",
    "# output to static HTML file\n",
    "output_file(\"images/log_lines.html\")\n",
    "\n",
    "# create a new plot\n",
    "p = figure(\n",
    "   tools=\"pan,box_zoom,reset,save\",\n",
    "   y_axis_type=\"log\", y_range=[0.001, 10**11], title=\"log axis example\",\n",
    "   x_axis_label='sections', y_axis_label='particles'\n",
    ")\n",
    "\n",
    "# add some renderers\n",
    "p.line(x, x, legend_label=\"y=x\")\n",
    "p.circle(x, x, legend_label=\"y=x\", fill_color=\"white\", size=8)\n",
    "p.line(x, y0, legend_label=\"y=x^2\", line_width=3)\n",
    "p.line(x, y1, legend_label=\"y=10^x\", line_color=\"red\")\n",
    "p.circle(x, y1, legend_label=\"y=10^x\", fill_color=\"red\", line_color=\"red\", size=6)\n",
    "p.line(x, y2, legend_label=\"y=10^x^2\", line_color=\"orange\", line_dash=\"4 4\")\n",
    "\n",
    "# show the results\n",
    "show(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}