2019-07-10_statistical_thinking_1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "- Notebook Author: [Trenton McKinney](https://trenton3983.github.io/)\n",
    "- Course: **[DataCamp: Statistical Thinking in Python (Part 1)](https://learn.datacamp.com/courses/statistical-thinking-in-python-part-1)**\n",
    "- This [notebook](https://github.com/trenton3983/DataCamp/blob/master/2019-07-10_statistical_thinking_1.ipynb) was created as a reproducible reference.\n",
    "- The material is from the course\n",
    "- I completed the exercises\n",
    "- If you find the content beneficial, consider a [DataCamp Subscription](https://www.datacamp.com/pricing?period=yearly)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Course Description\n",
    "\n",
    "Learn the principles of statistical inference to make clear, succinct conclusions from data. This course builds the foundation needed to think statistically and understand your data using Python-based tools. By the end of the course, you will be able to analyze data and draw meaningful conclusions efficiently."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": [
    "### Synopsis\n",
    "\n",
    "#### Graphical Exploratory Data Analysis\n",
    "- **Objectives:**\n",
    "  - Understand the importance of visualizing data.\n",
    "  - Learn various graphical techniques to explore data patterns, including histograms, box plots, scatter plots, swarm plots, and ECDFs.\n",
    "- **Methods:**\n",
    "  - Use histograms to display distributions.\n",
    "  - Create box plots to visualize data spread and detect outliers.\n",
    "  - Generate scatter plots to observe relationships between variables.\n",
    "  - Utilize swarm plots for detailed distribution views.\n",
    "  - Apply ECDFs (Empirical Cumulative Distribution Functions) to understand data distributions.\n",
    "- **Findings:**\n",
    "  - Graphical methods reveal underlying patterns, distributions, and potential outliers.\n",
    "  - Visualization aids in understanding relationships between variables and distributions of data.\n",
    "\n",
    "#### Quantitative Exploratory Data Analysis\n",
    "- **Objectives:**\n",
    "  - Perform numerical summaries of data.\n",
    "  - Understand measures of central tendency and variability.\n",
    "- **Methods:**\n",
    "  - Calculate mean, median, variance, and standard deviation.\n",
    "  - Use pandas for descriptive statistics and data summarization.\n",
    "- **Findings:**\n",
    "  - Quantitative summaries provide a numerical snapshot of data.\n",
    "  - Measures of central tendency and variability help in understanding data distribution.\n",
    "\n",
    "#### Thinking Probabilistically: Discrete Variables\n",
    "- **Objectives:**\n",
    "  - Develop a probabilistic mindset for analyzing discrete data.\n",
    "  - Understand probability mass functions and cumulative distribution functions for discrete variables.\n",
    "- **Methods:**\n",
    "  - Calculate probabilities for discrete events.\n",
    "  - Use numpy for generating random variables and computing distributions.\n",
    "- **Findings:**\n",
    "  - Probabilistic thinking aids in predicting the likelihood of discrete outcomes.\n",
    "  - Understanding distributions is key to analyzing discrete data.\n",
    "\n",
    "#### Thinking Probabilistically: Continuous Variables\n",
    "- **Objectives:**\n",
    "  - Apply probabilistic thinking to continuous data.\n",
    "  - Understand probability density functions and cumulative distribution functions for continuous variables.\n",
    "- **Methods:**\n",
    "  - Use numpy and scipy to work with continuous random variables.\n",
    "  - Compute and visualize continuous distributions.\n",
    "- **Findings:**\n",
    "  - Continuous distributions provide insights into the likelihood of continuous outcomes.\n",
    "  - Analyzing continuous data requires understanding of probability density functions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.patches import Rectangle\n",
    "import numpy as np\n",
    "from pprint import pprint as pp\n",
    "import csv\n",
    "from pathlib import Path\n",
    "import seaborn as sns\n",
    "from scipy.stats import binom\n",
    "\n",
    "from sklearn.datasets import load_iris"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Pandas Configuration Options"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.set_option('display.max_columns', 200)\n",
    "pd.set_option('display.max_rows', 300)\n",
    "pd.set_option('display.expand_frame_repr', True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data Files Location\n",
    "\n",
    "- Most data files for the exercises can be found on the [course site](https://www.datacamp.com/courses/statistical-thinking-in-python-part-1)\n",
    "  - [2008 election results (all states)](https://assets.datacamp.com/production/repositories/469/datasets/8fb59b9a99957c3b9b1c82b623aea54d8ccbcd9f/2008_all_states.csv)\n",
    "  - [2008 election results (swing states)](https://assets.datacamp.com/production/repositories/469/datasets/e079fddb581197780e1a7b7af2aeeff7242535f0/2008_swing_states.csv)\n",
    "  - [Belmont Stakes](https://assets.datacamp.com/production/repositories/469/datasets/7507bfed990379f246b4f166ea8a57ecf31c6c9d/belmont.csv)\n",
    "  - [Speed of light](https://assets.datacamp.com/production/repositories/469/datasets/df23780d215774ff90be0ea93e53f4fb5ebbade8/michelson_speed_of_light.csv)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Data File Objects**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = Path.cwd() / 'data' / '2019-07-10_statistical_thinking_1'\n",
    "elections_all_file = data / '2008_all_states.csv'\n",
    "elections_swing_file = data / '2008_swing_states.csv'\n",
    "belmont_file = data / 'belmont.csv'\n",
    "sol_file = data / 'michelson_speed_of_light.csv'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Iris Data Set**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "iris = load_iris()\n",
    "iris_df = pd.DataFrame(data=np.c_[iris['data'], iris['target']], columns=iris['feature_names'] + ['target'])\n",
    "\n",
    "def iris_typing(x):\n",
    "    types = {0.0: 'setosa',\n",
    "             1.0: 'versicolour',\n",
    "             2.0: 'virginica'}\n",
    "    return types[x]\n",
    "\n",
    "iris_df['species'] = iris_df.target.apply(iris_typing)\n",
    "iris_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Graphical exploratory data analysis\n",
    "\n",
    "Look before you leap! A very important proverb, indeed. Prior to diving in headlong into sophisticated statistical inference techniques, you should first explore your data by plotting them and computing simple summary statistics. This process, called exploratory data analysis, is a crucial first step in statistical analysis of data. So it is a fitting subject for the first chapter of Statistical Thinking in Python."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Introduction to exploratory data analysis\n",
    "\n",
    "* Exploring the data is a crucial step of the analysis.\n",
    "    * Organizing\n",
    "    * Plotting\n",
    "    * Computing numerical summaries\n",
    "* This idea is known as exploratory data analysis (EDA)\n",
    "* \"Exploratory data analysis can never be the whole story, but nothing else can serve as the foundation stone.\" - [John Tukey](https://en.wikipedia.org/wiki/John_Tukey)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "swing = pd.read_csv(elections_swing_file)\n",
    "swing.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* The raw data isn't particularly informative\n",
    "* We could start computing parameters and their confidence intervals and do hypothesis test...\n",
    "* ...however, we should graphically explore the data first"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Tukey's comments on EDA\n",
    "\n",
    "Even though you probably have not read Tukey's book, I suspect you already have a good idea about his viewpoint from the video introducing you to exploratory data analysis. Which of the following quotes is not directly from Tukey?\n",
    "\n",
    "1. Exploratory data analysis is detective work.\n",
    "1. There is no excuse for failing to plot and look.\n",
    "1. The greatest value of a picture is that it forces us to notice what we never expected to see.\n",
    "1. It is important to understand what you can do before you learn how to measure how well you seem to have done it.\n",
    "1. ~~**Often times EDA is too time consuming, so it is better to jump right in and do your hypothesis tests.**~~"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Advantages of graphical EDA\n",
    "\n",
    "Which of the following is not true of graphical EDA?\n",
    "\n",
    "1. It often involves converting tabular data into graphical form.\n",
    "1. If done well, graphical representations can allow for more rapid interpretation of data.\n",
    "1. ~~**A nice looking plot is always the end goal of a statistical analysis.**~~\n",
    "1. There is no excuse for neglecting to do graphical EDA."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting a histogram\n",
    "\n",
    "* always label the axes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bin_edges = [x for x in range(0, 110, 10)]\n",
    "plt.hist(x=swing.dem_share, bins=bin_edges, edgecolor='black')\n",
    "plt.xticks(bin_edges)\n",
    "plt.yticks(bin_edges[:-1])\n",
    "plt.xlabel('Percent of vote for Obama')\n",
    "plt.ylabel('Number of Counties')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Seaborn**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import seaborn as sns\n",
    "sns.set()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(x=swing.dem_share)\n",
    "plt.xlabel('Percent of vote for Obama')\n",
    "plt.ylabel('Number of Counties')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plotting a histogram of iris data\n",
    "\n",
    "For the exercises in this section, you will use a classic data set collected by botanist Edward Anderson and made famous by Ronald Fisher, one of the most prolific statisticians in history. Anderson carefully measured the anatomical properties of samples of three different species of iris, *Iris setosa*, *Iris versicolor*, and *Iris virginica*. The full data set is [available as part of scikit-learn](https://scikit-learn.org/stable/modules/generated/sklearn.datasets.load_iris.html). Here, you will work with his measurements of petal length.\n",
    "\n",
    "Plot a histogram of the petal lengths of his 50 samples of Iris versicolor using matplotlib/seaborn's default settings. Recall that to specify the default seaborn style, you can use `sns.set()`, where `sns` is the alias that `seaborn` is imported as.\n",
    "\n",
    "The subset of the data set containing the Iris versicolor petal lengths in units of centimeters (cm) is stored in the NumPy array `versicolor_petal_length`.\n",
    "\n",
    "In the video, Justin plotted the histograms by using the `pandas` library and indexing the DataFrame to extract the desired column. Here, however, you only need to use the provided NumPy array. Also, Justin assigned his plotting statements (except for `plt.show()`) to the dummy variable `_`. This is to prevent unnecessary output from being displayed. It is not required for your solutions to these exercises, however it is good practice to use it. Alternatively, if you are working in an interactive environment such as a Jupyter notebook, you could use a `;` after your plotting statements to achieve the same effect. Justin prefers using `_`. Therefore, you will see it used in the solution code.\n",
    "\n",
    "**Instructions**\n",
    "\n",
    "* Import `matplotlib.pyplot` and `seaborn` as their usual aliases (`plt` and `sns`).\n",
    "* Use `seaborn` to set the plotting defaults.\n",
    "* Plot a histogram of the Iris versicolor petal lengths using `plt.hist()` and the provided NumPy array `versicolor_petal_length`.\n",
    "* Show the histogram using `plt.show()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "versicolor_petal_length = iris_df['petal length (cm)'][iris_df.species == 'versicolour']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(versicolor_petal_length)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Axis labels!\n",
    "\n",
    "In the last exercise, you made a nice histogram of petal lengths of Iris versicolor, but **you didn't label the axes!** That's ok; it's not your fault since we didn't ask you to. Now, add axis labels to the plot using `plt.xlabel()` and `plt.ylabel()`. Don't forget to add units and assign both statements to `_`. The packages `matplotlib.pyplot` and `seaborn` are already imported with their standard aliases. This will be the case in what follows, unless specified otherwise.\n",
    "\n",
    "**Instructions**\n",
    "\n",
    "* Label the axes. Don't forget that you should always include units in your axis labels. Your y-axis label is just `'count'`. Your x-axis label is `'petal length (cm)'`. The units are essential!\n",
    "* Display the plot constructed in the above steps using `plt.show()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(versicolor_petal_length)\n",
    "plt.xlabel('petal length (cm)')\n",
    "plt.ylabel('count')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Adjusting the number of bins in a histogram\n",
    "\n",
    "The histogram you just made had ten bins. This is the default of matplotlib. The \"square root rule\" is a commonly-used rule of thumb for choosing number of bins: choose the number of bins to be the square root of the number of samples. Plot the histogram of Iris versicolor petal lengths again, this time using the square root rule for the number of bins. You specify the number of bins using the `bins` keyword argument of `plt.hist()`.\n",
    "\n",
    "The plotting utilities are already imported and the seaborn defaults already set. The variable you defined in the last exercise, `versicolor_petal_length`, is already in your namespace.\n",
    "\n",
    "**Instructions**\n",
    "\n",
    "* Import `numpy` as `np`. This gives access to the square root function, `np.sqrt()`.\n",
    "* Determine how many data points you have using `len()`.\n",
    "* Compute the number of bins using the square root rule.\n",
    "* Convert the number of bins to an integer using the built in `int()` function.\n",
    "* Generate the histogram and make sure to use the `bins` keyword argument.\n",
    "* Hit 'Submit Answer' to plot the figure and see the fruit of your labors!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute number of data points: n_data\n",
    "n_data = len(versicolor_petal_length)\n",
    "\n",
    "# Number of bins is the square root of number of data points: n_bins\n",
    "n_bins = np.sqrt(n_data)\n",
    "\n",
    "# Convert number of bins to integer: n_bins\n",
    "n_bins = int(n_bins)\n",
    "\n",
    "# Plot the histogram\n",
    "_ = plt.hist(versicolor_petal_length, bins=n_bins)\n",
    "\n",
    "# Label axes\n",
    "_ = plt.xlabel('petal length (cm)')\n",
    "_ = plt.ylabel('count')\n",
    "\n",
    "# Show histogram\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting all of your data: Bee swarm plots\n",
    "\n",
    "* Binning Bias: The same data may be interpreted differently depending on choice of bins\n",
    "* Additionally, all of the data isn't being plotted; the precision of the actual data is lost in the bins\n",
    "* These issues can be resolved with swarm plots\n",
    "    * Point position along the y-axis is the quantitative information\n",
    "    * The data are spread in x to make them visible, but their precise location along the x-axis is unimportant\n",
    "    * No binning bias and all the data are displayed.\n",
    "    * Seaborn & Pandas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.swarmplot(x='state', y='dem_share', data=swing)\n",
    "plt.xlabel('state')\n",
    "plt.ylabel('percent of vote for Obama')\n",
    "plt.title('% of Vote per Swing State County')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Bee swarm plot\n",
    "\n",
    "Make a bee swarm plot of the iris petal lengths. Your x-axis should contain each of the three species, and the y-axis the petal lengths. A data frame containing the data is in your namespace as `df`.\n",
    "\n",
    "For your reference, the code Justin used to create the bee swarm plot in the video is provided below:\n",
    "\n",
    "```python\n",
    "_ = sns.swarmplot(x='state', y='dem_share', data=df_swing)\n",
    "_ = plt.xlabel('state')\n",
    "_ = plt.ylabel('percent of vote for Obama')\n",
    "plt.show()\n",
    "```\n",
    "\n",
    "In the IPython Shell, you can use `sns.swarmplot`? or `help(sns.swarmplot)` for more details on how to make bee swarm plots using seaborn.\n",
    "\n",
    "**Instructions**\n",
    "\n",
    "* In the IPython Shell, inspect the DataFrame `df` using `df.head()`. This will let you identify which column names you need to pass as the `x` and `y` keyword arguments in your call to `sns.swarmplot()`.\n",
    "* Use `sns.swarmplot()` to make a bee swarm plot from the DataFrame containing the Fisher iris data set, `df`. The x-axis should contain each of the three species, and the y-axis should contain the petal lengths.\n",
    "* Label the axes.\n",
    "* Show your plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.swarmplot(x='species', y='petal length (cm)', data=iris_df, size=3)\n",
    "plt.xlabel('species')\n",
    "plt.ylabel('petal length (cm)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Interpreting a bee swarm plot\n",
    "\n",
    "Which of the following conclusions could you draw from the bee swarm plot of iris petal lengths you generated in the previous exercise? For your convenience, the bee swarm plot is regenerated and shown to the right.\n",
    "\n",
    "**Instructions**\n",
    "\n",
    "Possible Answers\n",
    "1. All I. versicolor petals are shorter than I. virginica petals.\n",
    "1. I. setosa petals have a broader range of lengths than the other two species.\n",
    "1. __**I. virginica petals tend to be the longest, and I. setosa petals tend to be the shortest of the three species.**__\n",
    "1. I. versicolor is a hybrid of I. virginica and I. setosa."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting all of your data: Empirical cumulative distribution functions\n",
    "\n",
    "* [Empirical Distribution Function](https://en.wikipedia.org/wiki/Empirical_distribution_function)\n",
    "* [Empirical Distribution Function / Empirical CDF](https://www.statisticshowto.datasciencecentral.com/empirical-distribution-function/)\n",
    "    * An **empirical cumulative distribution function** (also called the empirical distribution function, **ECDF**, or just **EDF**) and a [cumulative distribution function](https://www.statisticshowto.datasciencecentral.com/cumulative-distribution-function/) are basically the same thing; they are both probability models for data. While ***a CDF is a hypothetical model of a distribution***, **the ECDF models empirical (i.e. observed) data**. To put this another way, __**the ECDF is the [probability distribution](https://www.statisticshowto.datasciencecentral.com/probability-distribution/) you would get if you sampled from your [sample](https://www.statisticshowto.datasciencecentral.com/sample/), instead of the [population](https://www.statisticshowto.datasciencecentral.com/what-is-a-population/)**__. Lets say you have a set of experimental (observed) data $x_{1},x_{2},\\,\\ldots\\,x_{n}$. The EDF will give you the fraction of sample observations less than or equal to a particular value of $x$.\n",
    "    * More formally, if you have a set of [order statistics](https://www.statisticshowto.datasciencecentral.com/order-statistics/) ($y_{1}<y_{2}<\\ldots<y_{n}$) from an observed random sample, then the ***empirical distribution function*** is defined as a sum of [iid random variables](https://www.statisticshowto.datasciencecentral.com/random-variable/):\n",
    "    * \\$$\\hat{F}_{n}(t)=\\frac{\\text{number of elements in the sample \\(\\leq t\\)}}{n}=\\frac{1}{n}\\sum_{i=1}^n 1_{X_{i}\\leq t} $$\n",
    "    * Where $1_{A}$ = the [indicator](https://en.wikipedia.org/wiki/Indicator_function) of [event](https://en.wikipedia.org/wiki/Event_(probability_theory)) $A$.\n",
    "* x-value of an ECDF is the quantity being measured\n",
    "* y-value is the fraction of data points that have a value smaller than the corresponding x-value\n",
    "* Shows all the data and gives a complete picture of how the data are distributed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.sort(swing['dem_share'])\n",
    "y = np.arange(1, len(x)+1) / len(x)\n",
    "\n",
    "plt.plot(x, y, marker='.', linestyle='none')\n",
    "plt.xlabel('percent of vote for Obama')\n",
    "plt.ylabel('ECDF')\n",
    "plt.margins(0.02)  # keep data off plot edges\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### add annotations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(10, 5))\n",
    "ax.margins(0.05)           # Default margin is 0.05, value 0 means fit\n",
    "\n",
    "x = np.sort(swing['dem_share'])\n",
    "y = np.arange(1, len(x)+1) / len(x)\n",
    "\n",
    "ax.plot(x, y, marker='.', linestyle='none')\n",
    "plt.xlabel('percent of vote for Obama')\n",
    "plt.ylabel('ECDF')\n",
    "\n",
    "ax.annotate('20% of counties had <= 36% vote for Obama', xy=(36, .2),\n",
    "            xytext=(40, 0.1), fontsize=10, arrowprops=dict(arrowstyle=\"->\", color='b'))\n",
    "\n",
    "ax.annotate('75% of counties had < 0.5 vote for Obama', xy=(50, .75),\n",
    "            xytext=(55, 0.6), fontsize=10, arrowprops=dict(arrowstyle=\"->\", color='b'))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### plot multiple ECDFs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(10, 5))\n",
    "ax.margins(0.05)           # Default margin is 0.05, value 0 means fit\n",
    "\n",
    "for state in swing.state.unique():\n",
    "    x = np.sort(swing['dem_share'][swing.state == state])\n",
    "    y = np.arange(1, len(x)+1) / len(x)\n",
    "    ax.plot(x, y, marker='.', linestyle='none', label=state)\n",
    "\n",
    "plt.xlabel('percent of vote for Obama')\n",
    "plt.ylabel('ECDF')\n",
    "plt.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Computing the ECDF\n",
    "\n",
    "In this exercise, you will write a function that takes as input a 1D array of data and then returns the `x` and `y` values of the ECDF. You will use this function over and over again throughout this course and its sequel. ECDFs are among the most important plots in statistical analysis. You can write your own function, `foo(x,y)` according to the following skeleton:\n",
    "\n",
    "```python\n",
    "def foo(a,b):\n",
    "    \"\"\"State what function does here\"\"\"\n",
    "    # Computation performed here\n",
    "    return x, y\n",
    "```\n",
    "    \n",
    "The function `foo()` above takes two arguments `a` and `b` and returns two values `x` and `y`. The function header `def foo(a,b):` contains the function signature `foo(a,b)`, which consists of the function name, along with its parameters. For more on writing your own functions, see [DataCamp's course Python Data Science Toolbox (Part 1)](https://www.datacamp.com/courses/python-data-science-toolbox-part-1)!\n",
    "\n",
    "**Instructions**\n",
    "\n",
    "* Define a function with the signature `ecdf(data)`. Within the function definition,\n",
    "    * Compute the number of data points, `n`, using the `len()` function.\n",
    "    * The **x**-values are the sorted data. Use the `np.sort()` function to perform the sorting.\n",
    "    * The **y** data of the ECDF go from `1/n` to `1` in equally spaced increments. You can construct this using `np.arange()`. Remember, however, that the end value in `np.arange()` is not inclusive.  Therefore, `np.arange()` will need to go from `1` to `n+1`. Be sure to divide this by `n`.\n",
    "    * The function returns the values `x` and `y`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### def ecdf()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ecdf(data):\n",
    "    \"\"\"Compute ECDF for a one-dimensional array of measurements.\"\"\"\n",
    "    # Number of data points: n\n",
    "    n = len(data)\n",
    "\n",
    "    # x-data for the ECDF: x\n",
    "    x = np.sort(data)\n",
    "\n",
    "    # y-data for the ECDF: y\n",
    "    y = np.arange(1, n+1) / n\n",
    "\n",
    "    return x, y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plotting the ECDF\n",
    "\n",
    "You will now use your `ecdf()` function to compute the ECDF for the petal lengths of Anderson's *Iris versicolor* flowers. You will then plot the ECDF. Recall that your `ecdf()` function returns two arrays so you will need to unpack them. An example of such unpacking is `x, y = foo(data)`, for some function `foo()`.\n",
    "\n",
    "**Instructions**\n",
    "\n",
    "* Use `ecdf()` to compute the ECDF of `versicolor_petal_length`. Unpack the output into `x_vers` and `y_vers`.\n",
    "* Plot the ECDF as dots. Remember to include `marker = '.'` and `linestyle = 'none'` in addition to `x_vers` and `y_vers` as arguments inside `plt.plot()`.\n",
    "* Label the axes. You can label the y-axis `'ECDF'`.\n",
    "* Show your plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute ECDF for versicolor data: x_vers, y_vers\n",
    "x, y = ecdf(versicolor_petal_length)\n",
    "\n",
    "# Generate plot\n",
    "plt.plot(x, y, marker='.', linestyle='none')\n",
    "\n",
    "# Label the axes\n",
    "plt.xlabel('Versicolor Petal Length (cm)')\n",
    "plt.ylabel('ECDF')\n",
    "\n",
    "# Display the plot\n",
    "plt.margins(0.02)  # keep data off plot edges\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Comparison of ECDFs\n",
    "\n",
    "ECDFs also allow you to compare two or more distributions (though plots get cluttered if you have too many). Here, you will plot ECDFs for the petal lengths of all three iris species. You already wrote a function to generate ECDFs so you can put it to good use!\n",
    "\n",
    "To overlay all three ECDFs on the same plot, you can use `plt.plot()` three times, once for each ECDF. Remember to include `marker='.'` and `linestyle='none'` as arguments inside `plt.plot()`.\n",
    "\n",
    "**Instructions**\n",
    "\n",
    "* Compute ECDFs for each of the three species using your `ecdf()` function. The variables `setosa_petal_length`, `versicolor_petal_length`, and `virginica_petal_length` are all in your namespace. Unpack the ECDFs into `x_set`, `y_set`, `x_vers`, `y_vers` and `x_virg`, `y_virg`, respectively.\n",
    "* Plot all three ECDFs on the same plot as dots. To do this, you will need three `plt.plot()` commands. Assign the result of each to `_`.\n",
    "* A legend and axis labels have been added for you, so hit 'Submit Answer' to see all the ECDFs!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "virginica_petal_length = iris_df['petal length (cm)'][iris_df.species == 'virginica']\n",
    "setosa_petal_length = iris_df['petal length (cm)'][iris_df.species == 'setosa']\n",
    "\n",
    "# Compute ECDFs\n",
    "x_set, y_set = ecdf(setosa_petal_length)\n",
    "x_vers, y_vers = ecdf(versicolor_petal_length)\n",
    "x_virg, y_virg = ecdf(virginica_petal_length)\n",
    "\n",
    "# Plot all ECDFs on the same plot\n",
    "plt.plot(x_set, y_set, marker='.', linestyle='none')\n",
    "plt.plot(x_vers, y_vers, marker='.', linestyle='none')\n",
    "plt.plot(x_virg, y_virg, marker='.', linestyle='none')\n",
    "\n",
    "# Annotate the plot\n",
    "plt.legend(('setosa', 'versicolor', 'virginica'), loc='lower right')\n",
    "_ = plt.xlabel('petal length (cm)')\n",
    "_ = plt.ylabel('ECDF')\n",
    "\n",
    "# Display the plot\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Onward toward the whole story\n",
    "\n",
    "* Start with graphical eda!\n",
    "\n",
    "**Coming up...**\n",
    "\n",
    "* Thinking probabilistically\n",
    "* Discrete and continuous distributions\n",
    "* The power of hacker statistics using np.random()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quantitative exploratory data analysis\n",
    "\n",
    "In the last chapter, you learned how to graphically explore data. In this chapter, you will compute useful summary statistics, which serve to concisely describe salient features of a data set with a few numbers."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Introduction to summary statistics: The sample mean and median\n",
    "\n",
    "* mean - average\n",
    "    * heavily influenced by outliers\n",
    "    * `np.mean()`\n",
    "* median - middle value of the sorted dataset\n",
    "    * immune to outlier influence\n",
    "    * `np.median()`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Means and medians\n",
    "\n",
    "Which one of the following statements is true about means and medians?\n",
    "\n",
    "**Possible Answers**\n",
    "\n",
    "* ~~An outlier can significantly affect the value of both the mean and the median.~~\n",
    "* **An outlier can significantly affect the value of the mean, but not the median.**\n",
    "* ~~Means and medians are in general both robust to single outliers.~~\n",
    "* ~~The mean and median are equal if there is an odd number of data points.~~"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Computing means\n",
    "\n",
    "The mean of all measurements gives an indication of the typical magnitude of a measurement. It is computed using `np.mean()`.\n",
    "\n",
    "**Instructions**\n",
    "\n",
    "* Compute the mean petal length of Iris versicolor from Anderson's classic data set. The variable `versicolor_petal_length` is provided in your namespace. Assign the mean to `mean_length_vers`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute the mean: mean_length_vers\n",
    "mean_length_vers = np.mean(versicolor_petal_length)\n",
    "\n",
    "# Print the result with some nice formatting\n",
    "print('I. versicolor:', mean_length_vers, 'cm')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### with pandas.DataFrame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "iris_df.groupby(['species']).mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Percentiles, outliers and box plots\n",
    "\n",
    "* The median is a special name for the 50th percentile\n",
    "    * 50% of the data are less than the median\n",
    "* The 25th percentile is the value of the data point that is greater than 25% of the sorted data\n",
    "* percentiles are useful summary statistics and can be computed using `np.percentile()`\n",
    "\n",
    "**Computing Percentiles**\n",
    "\n",
    "```python\n",
    "np.percentile(df_swing['dem_share'], [25, 50, 75])\n",
    "```\n",
    "\n",
    "![](https://raw.githubusercontent.com/trenton3983/DataCamp/master/Images/2019-07-10_statistical_thinking_1/box_plot.JPG)\n",
    "\n",
    "* Box plots are a graphical methode for displying summary statistics\n",
    "    * median is the middle line: 50th percentile\n",
    "    * bottom and top line of the box represent the 25th & 75th percentile, repectively\n",
    "    * the space between the 25th and 75th percentile is the interquartile range (IQR)\n",
    "    * Whiskers extent a distance of 1.5 time the IQR, or the extent of the data, whichever is less extreme\n",
    "    * Any points outside the whiskers are plotted as individual points, which we demarcate as outliers\n",
    "        * There is no single definition for an outlier, however, being more than 2 IQRs away from the median is a common criterion.\n",
    "        * An outlier is not necessarily erroneous\n",
    "    * Box plots are a great alternative to bee swarm plots, becasue bee swarm plots become too cluttered with large data sets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_states = pd.read_csv(elections_all_file)\n",
    "all_states.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.boxplot(x='east_west', y='dem_share', data=all_states)\n",
    "plt.xlabel('region')\n",
    "plt.ylabel('percent of vote for Obama')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Computing percentiles\n",
    "\n",
    "In this exercise, you will compute the percentiles of petal length of *Iris versicolor*.\n",
    "\n",
    "**Instructions**\n",
    "\n",
    "* Create `percentiles`, a NumPy array of percentiles you want to compute. These are the 2.5th, 25th, 50th, 75th, and 97.5th. You can do so by creating a list containing these ints/floats and convert the list to a NumPy array using `np.array()`. For example, `np.array([30, 50])` would create an array consisting of the 30th and 50th percentiles.\n",
    "* Use `np.percentile()` to compute the percentiles of the petal lengths from the Iris versicolor samples. The variable `versicolor_petal_length` is in your namespace."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Specify array of percentiles: percentiles\n",
    "percentiles = np.array([2.5, 25, 50, 75, 97.5])\n",
    "\n",
    "# Compute percentiles: ptiles_vers\n",
    "ptiles_vers = np.percentile(versicolor_petal_length, percentiles)\n",
    "\n",
    "# Print the result\n",
    "ptiles_vers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Comparing percentiles to ECDF\n",
    "\n",
    "To see how the percentiles relate to the ECDF, you will plot the percentiles of Iris versicolor petal lengths you calculated in the last exercise on the ECDF plot you generated in chapter 1. The percentile variables from the previous exercise are available in the workspace as `ptiles_vers` and `percentiles`.\n",
    "\n",
    "Note that to ensure the Y-axis of the ECDF plot remains between 0 and 1, you will need to rescale the `percentiles` array accordingly - in this case, dividing it by 100.\n",
    "\n",
    "**Instructions**\n",
    "\n",
    "* Plot the percentiles as red diamonds on the ECDF. Pass the x and y co-ordinates - `ptiles_vers` and `percentiles/100` - as positional arguments and specify the `marker='D'`, `color='red'` and `linestyle='none'` keyword arguments. The argument for the y-axis - `percentiles/100` has been specified for you."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the ECDF\n",
    "_ = plt.plot(x_vers, y_vers, '.')\n",
    "_ = plt.xlabel('petal length (cm)')\n",
    "_ = plt.ylabel('ECDF')\n",
    "\n",
    "# Overlay percentiles as red diamonds.\n",
    "_ = plt.plot(ptiles_vers, percentiles/100, marker='D', color='red', linestyle='none')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Box-and-whisker plot\n",
    "\n",
    "Making a box plot for the petal lengths is unnecessary because the iris data set is not too large and the bee swarm plot works fine. However, it is always good to get some practice. Make a box plot of the iris petal lengths. You have a pandas DataFrame, `df`, which contains the petal length data, in your namespace. Inspect the data frame `df` in the IPython shell using `df.head()` to make sure you know what the pertinent columns are.\n",
    "\n",
    "For your reference, the code used to produce the box plot in the video is provided below:\n",
    "\n",
    "```python\n",
    "_ = sns.boxplot(x='east_west', y='dem_share', data=df_all_states)\n",
    "\n",
    "_ = plt.xlabel('region')\n",
    "\n",
    "_ = plt.ylabel('percent of vote for Obama')\n",
    "```\n",
    "\n",
    "In the IPython Shell, you can use `sns.boxplot?` or `help(sns.boxplot)` for more details on how to make box plots using seaborn.\n",
    "\n",
    "**Instructions**\n",
    "\n",
    "* The set-up is exactly the same as for the bee swarm plot; you just call `sns.boxplot()` with the same keyword arguments as you would `sns.swarmplot()`. The x-axis is `'species'` and y-axis is `'petal length (cm)'`.\n",
    "* Don't forget to label your axes!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(10, 7))\n",
    "# Create box plot with Seaborn's default settings\n",
    "_ = sns.boxplot(x='species', y='petal length (cm)', data=iris_df)\n",
    "\n",
    "# Label the axes\n",
    "_ = plt.ylabel('petal length (cm)')\n",
    "_ = plt.xlabel('species')\n",
    "\n",
    "# Show the plot\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Variance and standard deviation\n",
    "\n",
    "* measures of spread\n",
    "* variance:\n",
    "    * The mean squared distance of the data from the mean\n",
    "    * $$variance = \\frac{1}{n}\\sum_{i=1}^{n}(x_{i} - \\overline{x})^2$$\n",
    "    * because of the squared quantity, variance doesn't have the same units as the measurement\n",
    "* standard deviation:\n",
    "    * $$\\sqrt{variance}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Variance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dem_share_fl = all_states.dem_share[all_states.state == 'FL']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.var(dem_share_fl)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_states_var = all_states[['state', 'total_votes', 'dem_votes', 'rep_votes', 'other_votes', 'dem_share']].groupby(['state']).var(ddof=0)\n",
    "all_states_var.dem_share.loc['FL']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_states_var.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Standard Deviation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.std(dem_share_fl)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.sqrt(np.var(dem_share_fl))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_states_std = all_states[['state', 'total_votes', 'dem_votes', 'rep_votes', 'other_votes', 'dem_share']].groupby(['state']).std(ddof=0)\n",
    "all_states_std.dem_share.loc['FL']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_states_std.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Computing the variance\n",
    "\n",
    "It is important to have some understanding of what commonly-used functions are doing under the hood. Though you may already know how to compute variances, this is a beginner course that does not assume so. In this exercise, we will explicitly compute the variance of the petal length of Iris veriscolor using the equations discussed in the videos. We will then use `np.var()` to compute it.\n",
    "\n",
    "**Instructions**\n",
    "\n",
    "* Create an array called differences that is the `difference` between the petal lengths (`versicolor_petal_length`) and the mean petal length. The variable `versicolor_petal_length` is already in your namespace as a NumPy array so you can take advantage of NumPy's vectorized operations.\n",
    "* Square each element in this array. For example, `x**2` squares each element in the array `x`. Store the result as `diff_sq`.\n",
    "* Compute the mean of the elements in `diff_sq` using `np.mean()`. Store the result as `variance_explicit`.\n",
    "* Compute the variance of `versicolor_petal_length` using `np.var()`. Store the result as `variance_np`.\n",
    "* Print both `variance_explicit` and `variance_np` in one `print` call to make sure they are consistent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Array of differences to mean: differences\n",
    "differences = versicolor_petal_length - np.mean(versicolor_petal_length)\n",
    "\n",
    "# Square the differences: diff_sq\n",
    "diff_sq = differences**2\n",
    "\n",
    "# Compute the mean square difference: variance_explicit\n",
    "variance_explicit = np.mean(diff_sq)\n",
    "\n",
    "# Compute the variance using NumPy: variance_np\n",
    "variance_np = np.var(versicolor_petal_length)\n",
    "\n",
    "# Print the results\n",
    "print(variance_explicit, variance_np)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### The standard deviation and the variance\n",
    "\n",
    "As mentioned in the video, the standard deviation is the square root of the variance. You will see this for yourself by computing the standard deviation using `np.std()` and comparing it to what you get by computing the variance with `np.var()` and then computing the square root.\n",
    "\n",
    "**Instructions**\n",
    "\n",
    "* Compute the variance of the data in the `versicolor_petal_length` array using `np.var()` and store it in a variable called `variance`.\n",
    "* Print the square root of this value.\n",
    "* Print the standard deviation of the data in the `versicolor_petal_length` array using `np.std()`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute the variance: variance\n",
    "variance = np.var(versicolor_petal_length)\n",
    "\n",
    "# Print the square root of the variance\n",
    "std_explicit = np.sqrt(variance)\n",
    "\n",
    "# Print the standard deviation\n",
    "std_np = np.std(versicolor_petal_length)\n",
    "\n",
    "print(std_explicit, std_np)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Covariance and Pearson correlation coefficient\n",
    "\n",
    "![covariance](https://raw.githubusercontent.com/trenton3983/DataCamp/master/Images/2019-07-10_statistical_thinking_1/covariance.JPG)\n",
    "\n",
    "* Covariance\n",
    "    * $$covariance = \\frac{1}{n}\\sum_{i=1}^{n}(x_{i} - \\overline{x})(y_{i} - \\overline{y})$$\n",
    "    * The data point differs from the mean vote share and the mean total votes for Obama\n",
    "    * The differences for each data point can be computed\n",
    "    * The covariance is the mean of the product of these differences\n",
    "    * If both x and y tend to be above or below their respective means together, as they are in this data set, the covariance is positive.\n",
    "        * This means they are positively correlated:\n",
    "            * When x is high, so is y\n",
    "            * When the county is populous, it has more votes for Obama\n",
    "    * If x is high while y is low, the covariance is negative\n",
    "        * This means they are negatively correlated (anticorrelated) - not the case for this data set.\n",
    "        \n",
    "![pearson](https://raw.githubusercontent.com/trenton3983/DataCamp/master/Images/2019-07-10_statistical_thinking_1/pearson.JPG)\n",
    "\n",
    "* Pearson correlation\n",
    "    * A more generally applicable measure of how two variables depend on each other, should be dimensionless (not units).\n",
    "    * $$\\rho = Pearson\\space correlation = \\frac{covariance}{(std\\space of\\space x)(std\\space of\\space y)}$$\n",
    "    * $$\\rho = \\frac{variability\\space due\\space to\\space codependence}{independent\\space variability}$$\n",
    "    * Comparison of the variability in the data due to codependence (the covariance) to the variability inherent to each variable independently (their standard deviations).\n",
    "    * It's dimensionless and ranges from -1 (for complete anticorrelation) to 1 (for complete correlation).\n",
    "    * A value of zero means there is no correlation between the data, as shown in the upper left plot.\n",
    "    * Good metric for correlation between two variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10, 8))\n",
    "sns.scatterplot(x='total_votes', y='dem_share', data=swing, hue='state')\n",
    "plt.xlabel('total votes')\n",
    "plt.ylabel('% of vote for Obama')\n",
    "plt.xticks([x for x in range(0, 1000000, 100000)], rotation=40)\n",
    "plt.yticks([x for x in range(0, 100, 10)])\n",
    "\n",
    "# Create a Rectangle patch\n",
    "plt.gca().add_patch(Rectangle((400000, 52), 500000, 34, linewidth=1, edgecolor='b', facecolor='none'))\n",
    "\n",
    "plt.gca().add_patch(Rectangle((0, 5), 50000, 45, linewidth=1, edgecolor='r', facecolor='none'))\n",
    "\n",
    "# Annotate\n",
    "plt.annotate('12 largest counties; most vote for Obama', xy=(650000, 52), weight='bold',\n",
    "             xytext=(400000, 35), fontsize=10, arrowprops=dict(arrowstyle=\"->\", color='b'))\n",
    "\n",
    "plt.annotate('small counties; most vote for McCain', xy=(50000, 20), weight='bold',\n",
    "             xytext=(150000, 7), fontsize=10, arrowprops=dict(arrowstyle=\"->\", color='r'))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Scatter plots\n",
    "\n",
    "When you made bee swarm plots, box plots, and ECDF plots in previous exercises, you compared the petal lengths of different species of iris. But what if you want to compare two properties of a single species? This is exactly what we will do in this exercise. We will make a **scatter plot** of the petal length and width measurements of Anderson's Iris versicolor flowers. If the flower scales (that is, it preserves its proportion as it grows), we would expect the length and width to be correlated.\n",
    "\n",
    "For your reference, the code used to produce the scatter plot in the video is provided below:\n",
    "\n",
    "```python\n",
    "_ = plt.plot(total_votes/1000, dem_share, marker='.', linestyle='none')\n",
    "_ = plt.xlabel('total votes (thousands)')\n",
    "_ = plt.ylabel('percent of vote for Obama')\n",
    "```\n",
    "\n",
    "**Instructions**\n",
    "\n",
    "* Use `plt.plot()` with the appropriate keyword arguments to make a scatter plot of versicolor petal length (x-axis) versus petal width (y-axis). The variables `versicolor_petal_length` and `versicolor_petal_width` are already in your namespace. Do not forget to use the `marker='.'` and `linestyle='none'` keyword arguments.\n",
    "* Label the axes.\n",
    "* Display the plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "versicolor_petal_width = iris_df['petal width (cm)'][iris_df.species == 'versicolour']\n",
    "\n",
    "# Make a scatter plot\n",
    "_ = plt.plot(versicolor_petal_length, versicolor_petal_width, marker='.', linestyle='none')\n",
    "\n",
    "# Label the axes\n",
    "_ = plt.xlabel('petal length (cm)')\n",
    "_ = plt.ylabel('petal width (cm)')\n",
    "\n",
    "# Show the result\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Variance and covariance by looking\n",
    "\n",
    "![](https://raw.githubusercontent.com/trenton3983/DataCamp/master/Images/2019-07-10_statistical_thinking_1/var_covar_by_looking.JPG)\n",
    "\n",
    "Consider four scatter plots of x-y data, appearing to the right. Which has, respectively,\n",
    "\n",
    "* the highest variance in the variable x,\n",
    "* the highest covariance,\n",
    "* negative covariance?\n",
    "\n",
    "**Instructions**\n",
    "\n",
    "Possible Answers\n",
    "* ~~a, c, b~~\n",
    "* ~~d, c, a~~\n",
    "* __**d, c, b**__\n",
    "* ~~d, d, b~~"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Computing the covariance\n",
    "\n",
    "The covariance may be computed using the Numpy function `np.cov()`. For example, we have two sets of data `x` and `y`, `np.cov(x, y)` returns a 2D array where entries `[0,1]` and `[1,0]` are the covariances. Entry `[0,0]` is the variance of the data in x, and entry `[1,1]` is the variance of the data in `y`. This 2D output array is called the covariance matrix, since it organizes the self- and covariance.\n",
    "\n",
    "To remind you how the *I. versicolor* petal length and width are related, we include the scatter plot you generated in a previous exercise.\n",
    "\n",
    "**Instructions**\n",
    "\n",
    "* Use `np.cov()` to compute the covariance matrix for the petal length (`versicolor_petal_length`) and width (`versicolor_petal_width`) of *I. versicolor*.\n",
    "* Print the covariance matrix.\n",
    "* Extract the covariance from entry `[0,1]` of the covariance matrix. Note that by symmetry, entry `[1,0]` is the same as entry `[0,1]`.\n",
    "* Print the covariance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "iris_df[['petal length (cm)', 'petal width (cm)']][iris_df.species == 'versicolour'].cov()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute the covariance matrix: covariance_matrix\n",
    "covariance_matrix = np.cov(versicolor_petal_length, versicolor_petal_width)\n",
    "\n",
    "# Print covariance matrix\n",
    "covariance_matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract covariance of length and width of petals: petal_cov\n",
    "petal_cov = covariance_matrix[0, 1]\n",
    "\n",
    "# Print the length/width covariance\n",
    "petal_cov"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Computing the Pearson correlation coefficient\n",
    "\n",
    "As mentioned in the video, the Pearson correlation coefficient, also called the Pearson r, is often easier to interpret than the covariance. It is computed using the `np.corrcoef()` function. Like `np.cov()`, it takes two arrays as arguments and returns a 2D array. Entries `[0,0]` and `[1,1]` are necessarily equal to 1 (can you think about why?), and the value we are after is entry `[0,1]`.\n",
    "\n",
    "In this exercise, you will write a function, `pearson_r(x, y)` that takes in two arrays and returns the Pearson correlation coefficient. You will then use this function to compute it for the petal lengths and widths of *I. versicolor*.\n",
    "\n",
    "Again, we include the scatter plot you generated in a previous exercise to remind you how the petal width and length are related.\n",
    "\n",
    "**Instructions**\n",
    "\n",
    "* Define a function with signature `pearson_r(x, y)`.\n",
    "    * Use `np.corrcoef()` to compute the correlation matrix of `x` and `y` (pass them to `np.corrcoef()` in that order).\n",
    "    * The function returns entry `[0,1]` of the correlation matrix.\n",
    "* Compute the Pearson correlation between the data in the arrays `versicolor_petal_length` and `versicolor_petal_width`. Assign the result to `r`.\n",
    "* Print the result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "iris_df[['petal length (cm)', 'petal width (cm)']][iris_df.species == 'versicolour'].corr()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pearson_r(x, y):\n",
    "    \"\"\"Compute Pearson correlation coefficient between two arrays.\"\"\"\n",
    "    # Compute correlation matrix: corr_mat\n",
    "    corr_mat = np.corrcoef(x, y)\n",
    "\n",
    "    # Return entry [0,1]\n",
    "    return corr_mat[0,1]\n",
    "\n",
    "# Compute Pearson correlation coefficient for I. versicolor: r\n",
    "r = pearson_r(versicolor_petal_length, versicolor_petal_width)\n",
    "\n",
    "# Print the result\n",
    "print(r)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Thinking probabilistically: Discrete variables\n",
    "\n",
    "Statistical inference rests upon probability. Because we can very rarely say anything meaningful with absolute certainty from data, we use probabilistic language to make quantitative statements about data. In this chapter, you will learn how to think probabilistically about discrete quantities, those that can only take certain values, like integers. It is an important first step in building the probabilistic language necessary to think statistically."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Probabilistic logic and statistical inference\n",
    "\n",
    "* Probabilistic reasoning allows us to describe uncertainty\n",
    "* Given a set of data, you describe probabilistically what you might expect if those data were acquired repeatedly\n",
    "* This is the heart of statistical inference\n",
    "* It's the process by which we go from measured data to probabilistic conclusions about what we might expect if we collected the same data again."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### What is the goal of statistical inference?\n",
    "\n",
    "Why do we do statistical inference?\n",
    "\n",
    "**Possible Answers**\n",
    "\n",
    "* To draw probabilistic conclusions about what we might expect if we collected the same data again.\n",
    "* To draw actionable conclusions from data.\n",
    "* To draw more general conclusions from relatively few data or observations.\n",
    "* __**All of these.**__"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Why do we use the language of probablility?\n",
    "\n",
    "Which of the following is not a reason why we use probabilistic language in statistical inference?\n",
    "\n",
    "**Possible Answers**\n",
    "\n",
    "* Probability provides a measure of uncertainty.\n",
    "* __**Probabilistic language is not very precise.**__\n",
    "* Data are almost never exactly the same when acquired again, and probability allows us to say how much we expect them to vary."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Random number generators and hacker statistics\n",
    "\n",
    "* Instead o repeating data acquisition over and over, repeated measurements can be simulated\n",
    "* The concepts of probabilities originated from games of chance\n",
    "* What's the probability of getting 4 heads with 4 flips of a coin?\n",
    "* This type of data can be generated using `np.random.random`\n",
    "    * drawn a number between 0 and 1\n",
    "    * $<0.5\\longrightarrow\\text{heads}$\n",
    "    * $\\geq0.5\\longrightarrow\\text{tails}$\n",
    "    * The pseudo random number generator works by starting with an integer, called a seed, and then generates random numbers in succession\n",
    "    * The same seed gives the same sequence of random numbers\n",
    "    * Manually seed the random number generator for reproducible results\n",
    "    * Specified using `np.random.seed()`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Bernoulli Trial\n",
    "\n",
    "* An experiment that has two options, \"success\" (True) and \"failure\" (False)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Hacker stats probabilities\n",
    "\n",
    "* Determine how to simulate data\n",
    "* Simulated it repeatedly\n",
    "* Compute the fraction of trials that had the outcome of interest\n",
    "    * Probability is approximately the fraction of trials with the outcome of interest"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Simulated coin flips"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(42)\n",
    "random_numbers = np.random.random(size=4)\n",
    "\n",
    "random_numbers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "heads = random_numbers < 0.5\n",
    "\n",
    "heads"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.sum(heads)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* The number of heads can be computed by summing the array of Booleans, because in numerical contexts, Python treats True as 1 and False as 0.\n",
    "* We want to know the probability of getting four heads if we were to repeatedly flip the 4 coins"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* without `list comprehension`\n",
    "\n",
    "```python\n",
    "n_all_heads = 0  # initialize number of 4-heads trials\n",
    "\n",
    "for _ in range(10000):\n",
    "    heads = np.random.random(size=4) < 0.5\n",
    "    n_heads = np.sum(heads)\n",
    "    if n_heads == 4:\n",
    "        n_all_heads += 1\n",
    "```\n",
    "\n",
    "* with `list comprehension`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_all_heads = sum([1 for _ in range(10000) if sum(np.random.random(size=4) < 0.5) == 4])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_all_heads"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_all_heads/10000"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Generating random numbers using the `np.random` module\n",
    "\n",
    "We will be hammering the `np.random` module for the rest of this course and its sequel. Actually, you will probably call functions from this module more than any other while wearing your hacker statistician hat. Let's start by taking its simplest function, `np.random.random()` for a test spin. The function returns a random number between zero and one. Call `np.random.random()` a few times in the IPython shell. You should see numbers jumping around between zero and one.\n",
    "\n",
    "In this exercise, we'll generate lots of random numbers between zero and one, and then plot a histogram of the results. If the numbers are truly random, all bars in the histogram should be of (close to) equal height.\n",
    "\n",
    "You may have noticed that, in the video, Justin generated 4 random numbers by passing the keyword argument `size=4` to `np.random.random()`. Such an approach is more efficient than a `for` loop: in this exercise, however, you will write a `for` loop to experience hacker statistics as the practice of repeating an experiment over and over again.\n",
    "\n",
    "**Instructions**\n",
    "\n",
    "* Seed the random number generator using the seed `42`.\n",
    "* Initialize an empty array, `random_numbers`, of 100,000 entries to store the random numbers. Make sure you use `np.empty(100000)` to do this.\n",
    "* Write a `for` loop to draw 100,000 random numbers using `np.random.random()`, storing them in the `random_numbers` array. To do so, loop over `range(100000)`.\n",
    "* Plot a histogram of `random_numbers`. It is not necessary to label the axes in this case because we are just checking the random number generator. Hit 'Submit Answer' to show your plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Seed the random number generator\n",
    "np.random.seed(42)\n",
    "\n",
    "# Initialize random numbers: random_numbers\n",
    "random_numbers = np.empty(100000)\n",
    "\n",
    "# Generate random numbers by looping over range(100000)\n",
    "for i in range(100000):\n",
    "    random_numbers[i] = np.random.random()\n",
    "\n",
    "# Plot a histogram\n",
    "_ = plt.hist(random_numbers)\n",
    "\n",
    "# Show the plot\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.histplot(random_numbers, kde=True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**The histogram is nearly flat across the top, indicating there is equal chance a randomly-generated number is in any of the histogram bins.**\n",
    "\n",
    "* [Generating Random Numbers With `NumPy`](https://chrisalbon.com/python/basics/generating_random_numbers_with_numpy/)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Using `np.random.rand`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rand_num = np.random.rand(100000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.histplot(rand_num, kde=True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### The np.random module and Bernoulli trials\n",
    "\n",
    "You can think of a Bernoulli trial as a flip of a possibly biased coin. Specifically, each coin flip has a probability ***p*** of landing heads (success) and probability ***1−p*** of landing tails (failure). In this exercise, you will write a function to perform `n` [Bernoulli trials](https://en.wikipedia.org/wiki/Bernoulli_trial), `perform_bernoulli_trials(n, p)`, which returns the number of successes out of `n` Bernoulli trials, each of which has probability `p` of success. To perform each Bernoulli trial, use the `np.random.random()` function, which returns a random number between zero and one.\n",
    "\n",
    "**Instructions**\n",
    "\n",
    "* Define a function with signature `perform_bernoulli_trials(n, p)`.\n",
    "    * Initialize to zero a variable `n_success` the counter of `True` occurrences, which are Bernoulli trial successes.\n",
    "    * Write a `for` loop where you perform a Bernoulli trial in each iteration and increment the number of success if the result is `True`. Perform `n` iterations by looping over `range(n)`.\n",
    "        * To perform a Bernoulli trial, choose a random number between zero and one using `np.random.random()`. If the number you chose is less than `p`, increment n_success (use the `+= 1` operator to achieve this).\n",
    "* The function returns the number of successes `n_success`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### def perform_bernoulli_trials()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def perform_bernoulli_trials(n: int=100000, p: float=0.5) -> int:\n",
    "    \"\"\"\n",
    "    Perform n Bernoulli trials with success probability p\n",
    "    and return number of successes.\n",
    "    n: number of iterations\n",
    "    p: target number between 0 and 1, inclusive\n",
    "    \"\"\"\n",
    "    # Initialize number of successes: n_success\n",
    "    n_success = 0\n",
    "\n",
    "    # Perform trials\n",
    "    for i in range(n):\n",
    "        # Choose random number between zero and one: random_number\n",
    "        random_number = np.random.random()\n",
    "\n",
    "        # If less than p, it's a success so add one to n_success\n",
    "        if random_number < p:\n",
    "            n_success += 1\n",
    "\n",
    "    return n_success"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### With `list comprehension`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def perform_bernoulli_trials(n: int=100000, p: float=0.5) -> int:\n",
    "    \"\"\"\n",
    "    Perform n Bernoulli trials with success probability p\n",
    "    and return number of successes.\n",
    "    n: number of iterations\n",
    "    p: target number between 0 and 1, inclusive\n",
    "    \"\"\"\n",
    "\n",
    "    return sum([1 for _ in range(n) if np.random.random() < p])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### How many defaults might we expect?\n",
    "\n",
    "Let's say a bank made 100 mortgage loans. It is possible that anywhere between 0 and 100 of the loans will be defaulted upon. You would like to know the probability of getting a given number of defaults, given that the probability of a default is `p = 0.05`. To investigate this, you will do a simulation. You will perform 100 Bernoulli trials using the `perform_bernoulli_trials()` function you wrote in the previous exercise and record how many defaults we get. Here, a success is a default. (Remember that the word \"success\" just means that the Bernoulli trial evaluates to `True`, i.e., did the loan recipient default?) You will do this for another 100 Bernoulli trials. And again and again until we have tried it 1000 times. Then, you will plot a histogram describing the probability of the number of defaults.\n",
    "\n",
    "**Instructions**\n",
    "\n",
    "* Seed the random number generator to 42.\n",
    "* Initialize `n_defaults`, an empty array, using `np.empty()`. It should contain 1000 entries, since we are doing 1000 simulations.\n",
    "* Write a `for` loop with `1000` iterations to compute the number of defaults per 100 loans using the `perform_bernoulli_trials()` function. It accepts two arguments: the number of trials `n` - in this case 100 - and the probability of success `p` - in this case the probability of a default, which is `0.05`. On each iteration of the loop store the result in an entry of `n_defaults`.\n",
    "* Plot a histogram of `n_defaults`. Include the `normed=True` keyword argument so that the height of the bars of the histogram indicate the probability."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Seed random number generator\n",
    "np.random.seed(42)\n",
    "\n",
    "# Initialize the number of defaults: n_defaults\n",
    "n_defaults = np.empty(1000)\n",
    "\n",
    "# Compute the number of defaults\n",
    "for i in range(1000):\n",
    "    n_defaults[i] = perform_bernoulli_trials(100, 0.05)\n",
    "\n",
    "\n",
    "# Plot the histogram with default number of bins; label your axes\n",
    "_ = plt.hist(n_defaults, density=True)\n",
    "_ = plt.xlabel('number of defaults out of 100 loans')\n",
    "_ = plt.ylabel('probability')\n",
    "\n",
    "# Show the plot\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**This is not an optimal way to plot a histogram when the results are known to be integers. This will be revisited in forthcoming exercises.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### With `list comprehension`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(42)\n",
    "n_defaults = np.asarray([perform_bernoulli_trials(100, 0.05) for _ in range(1000)])\n",
    "\n",
    "plt.hist(n_defaults, density=True)\n",
    "plt.xlabel('number of defaults out of 100 loans')\n",
    "plt.ylabel('probability')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Will the bank fail?\n",
    "\n",
    "Using `def ecdf()` from the first section, plot the number of `n_defaults` from the previous exercise, as a CDF.\n",
    "\n",
    "If interest rates are such that the bank will lose money if 10 or more of its loans are defaulted upon, what is the probability that the bank will lose money?\n",
    "\n",
    "**Instructions**\n",
    "\n",
    "* Compute the `x` and `y` values for the ECDF of `n_defaults`.\n",
    "* Plot the ECDF, making sure to label the axes. Remember to include `marker='.'` and `linestyle='none'` in addition to `x` and `y` in your call `plt.plot()`.\n",
    "* Show the plot.\n",
    "* Compute the total number of entries in your `n_defaults` array that were greater than or equal to 10. To do so, compute a boolean array that tells you whether a given entry of `n_defaults` is `>= 10`. Then sum all the entries in this array using `np.sum()`. For example, `np.sum(n_defaults <= 5)` would compute the number of defaults with 5 or *fewer* defaults.\n",
    "* The probability that the bank loses money is the fraction of `n_defaults` that are greater than or equal to 10."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute ECDF: x, y\n",
    "x, y = ecdf(n_defaults)\n",
    "\n",
    "# Plot the ECDF with labeled axes\n",
    "plt.plot(x, y, marker='.', linestyle='none')\n",
    "plt.xlabel('Number of Defaults out of 100')\n",
    "plt.ylabel('CDF')\n",
    "\n",
    "# Show the plot\n",
    "plt.show()\n",
    "\n",
    "# Compute the number of 100-loan simulations with 10 or more defaults: n_lose_money\n",
    "n_lose_money = sum(n_defaults >= 10)\n",
    "\n",
    "# Compute and print probability of losing money\n",
    "print('Probability of losing money =', n_lose_money / len(n_defaults))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**As might be expected, about 5/100 defaults occur. There's about a 2% chance of getting 10 or more defaults out of 100 loans.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Probability distributions and stories: The Binomial distribution\n",
    "\n",
    "* [Probability Distributions in Python](https://www.datacamp.com/community/tutorials/probability-distributions-python)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Probability Mass Function (PMF)\n",
    "\n",
    "* [Probability mass function](https://en.wikipedia.org/wiki/Probability_mass_function)\n",
    "* The set of probabilities of discrete outcomes\n",
    "* PMF is a property of a discrete probability distribution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Discrete Uniform PMF\n",
    "\n",
    "* ![](https://raw.githubusercontent.com/trenton3983/DataCamp/master/Images/2019-07-10_statistical_thinking_1/discrete_uniform_pmf.JPG)\n",
    "* The outcomes are discrete because only certain values may be attained; there is not option for 3.7\n",
    "* Each result has a uniform probability of 1/6"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Probability Distribution\n",
    "\n",
    "* [Probability distribution](https://en.wikipedia.org/wiki/Probability_distribution)\n",
    "* A mathematical description of outcomes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Discrete Uniform Distribution\n",
    "\n",
    "* [Discrete uniform distribution](https://en.wikipedia.org/wiki/Discrete_uniform_distribution)\n",
    "* The outcome of rolling a single fair die, is Discrete Uniformly distributed"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Binomial Distribution\n",
    "\n",
    "* [Binomial distribution](https://en.wikipedia.org/wiki/Binomial_distribution)\n",
    "* The number ***r*** of successes in ***n*** Bernoulli trials with probability ***p*** of success, is Binomially distributed\n",
    "* The number ***r*** of heads in 4 coin flips with probability ***p = 0.5*** of heads, is Binomially distributed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.binomial(4, 0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.binomial(4, 0.5, size=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Binomial PMF\n",
    "\n",
    "* To plot the Binomial PMF, take 10000 samples from a Binomial distribution of 60 Bernoulli trials with a probability of success of 0.1\n",
    "* The most likely number of successes is 6 out of 60, but it's possible to get as many as 11 or as few as 1\n",
    "* [`scipy.stats.binom`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.binom.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(42)\n",
    "samples = np.random.binomial(60, 0.1, size=10_000)\n",
    "samples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](https://raw.githubusercontent.com/trenton3983/DataCamp/master/Images/2019-07-10_statistical_thinking_1/binomial_pmf.JPG)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n, p = 60, 0.1\n",
    "x = [x for x in range(17)]\n",
    "\n",
    "fig, ax = plt.subplots(1, 1)\n",
    "ax.plot(x, binom.pmf(x, n, p), 'bo', ms=5, label='binom pmf')\n",
    "ax.vlines(x, 0, binom.pmf(x, n, p), colors='b', lw=3, alpha=0.5)\n",
    "plt.xticks(x)\n",
    "plt.ylabel('probability')\n",
    "plt.xlabel('number of successes')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set()\n",
    "x, y = ecdf(samples)\n",
    "\n",
    "plt.plot(x, y, marker='.', linestyle='none')\n",
    "plt.margins(0.02)\n",
    "plt.xlabel('Number of Successes')\n",
    "plt.ylabel('CDF')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Sampling out of the Binomial distribution\n",
    "\n",
    "Compute the probability mass function for the number of defaults we would expect for 100 loans as in the last section, but instead of simulating all of the Bernoulli trials, perform the sampling using `np.random.binomial()`. This is identical to the calculation you did in the last set of exercises using your custom-written `perform_bernoulli_trials()` function, but far more computationally efficient. Given this extra efficiency, we will take 10,000 samples instead of 1000. After taking the samples, plot the CDF as last time. This CDF that you are plotting is that of the Binomial distribution.\n",
    "\n",
    "Note: For this exercise and all going forward, the random number generator is pre-seeded for you (with `np.random.seed(42)`) to save you typing that each time.\n",
    "\n",
    "**Instructions**\n",
    "\n",
    "* Draw samples out of the Binomial distribution using `np.random.binomial()`. You should use parameters `n = 100` and `p = 0.05`, and set the `size = 10000`.\n",
    "* Compute the CDF using your previously-written `ecdf()` function.\n",
    "* Plot the CDF with axis labels. The x-axis here is the ***number of defaults out of 100 loans***, while the y-axis is the ***CDF***."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Take 10,000 samples out of the binomial distribution: n_defaults\n",
    "np.random.seed(42)\n",
    "n_defaults = np.random.binomial(100, 0.05, size=10_000)\n",
    "\n",
    "# Compute CDF: x, y\n",
    "x, y = ecdf(n_defaults)\n",
    "\n",
    "# Plot the CDF with axis labels\n",
    "plt.plot(x, y, marker='.', linestyle='none')\n",
    "plt.xlabel('number of defaults out of 100 loans')\n",
    "plt.ylabel('CDF')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plotting the Binomial PMF\n",
    "\n",
    "As mentioned in the video, plotting a nice looking PMF requires a bit of matplotlib trickery that we will not go into here. Instead, we will plot the PMF of the Binomial distribution as a histogram with skills you have already learned. The trick is setting up the edges of the bins to pass to `plt.hist()` via the `bins` keyword argument. We want the bins centered on the integers. So, the edges of the bins should be `-0.5, 0.5, 1.5, 2.5, ...` up to `max(n_defaults) + 1.5`. You can generate an array like this using `np.arange()` and then subtracting `0.5` from the array.\n",
    "\n",
    "You have already sampled out of the Binomial distribution during your exercises on loan defaults, and the resulting samples are in the NumPy array `n_defaults`.\n",
    "\n",
    "**Instructions**\n",
    "\n",
    "* Using `np.arange()`, compute the bin edges such that the bins are centered on the integers. Store the resulting array in the variable `bins`.\n",
    "* Use `plt.hist()` to plot the histogram of `n_defaults` with the `normed=True` and `bins=bins` keyword arguments."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute bin edges: bins\n",
    "bins = np.arange(0, max(n_defaults) + 1.5) - 0.5\n",
    "\n",
    "# Generate histogram\n",
    "plt.hist(n_defaults, density=True, bins=bins)\n",
    "\n",
    "# Label axes\n",
    "plt.xlabel('number of defaults out of 100 loans')\n",
    "plt.ylabel('PMF')\n",
    "\n",
    "# Show the plot\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Poisson processes and the Poisson distribution\n",
    "\n",
    "* [Poisson distribution](https://en.wikipedia.org/wiki/Poisson_distribution)\n",
    "* The timing of the next event is completely independent of when the previous event occurred\n",
    "* Examples of Poisson processes:\n",
    "    * Natural births in a given hospital\n",
    "        * There is a well-defined average number of natural births per year, and the timing of one birth is independent of the timing of the previous one\n",
    "    * Hits on a website during a given hour\n",
    "        * The timing of successive hits is independent of the timing of the previous hit\n",
    "    * Meteor strikes\n",
    "    * Molecular collisions in a gas\n",
    "    * Aviation incidents\n",
    "* The number of arrivals of a Poisson process in a given amount of time is Poisson distributed\n",
    "* The number of arrivals ***r*** of a Poisson process in a given time interval with average rate of arrivals $\\lambda$ per interval is Poisson distributed\n",
    "    * The Poisson distribution has one parameter, the average number of arrivals in a given length of time\n",
    "* The number of hits ***r*** on a website in one hour with an average hit rate of 6 hits per hour is Poisson distributed"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Poisson PMF\n",
    "\n",
    "* ![poisson pmf](https://raw.githubusercontent.com/trenton3983/DataCamp/master/Images/2019-07-10_statistical_thinking_1/poisson_pmf.JPG)\n",
    "* For the preceding plot, for a given hour, the site is likely to get 6 hits, which is the average, but it's possible to also get 10, or none\n",
    "* This looks like the Binomial PMF from **3.3.0.5.1. Binomial PMF**\n",
    "    The Poisson distribution is a limit of the Binomial distribution for low probability of success and large number of trials, i.e. for rare events\n",
    "* To sample from the Poisson distribution, use `np.random.poisson`.\n",
    "    * It also has the size keyword argument to allow multiple samples\n",
    "    * The Poisson CDF resembles the Binomial CDF"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Poisson CDF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "samples = np.random.poisson(6, size=10_000)\n",
    "x, y = ecdf(samples)\n",
    "\n",
    "plt.plot(x, y, marker='.', linestyle='none')\n",
    "plt.margins(0.02)\n",
    "plt.xlabel('number of successes')\n",
    "plt.ylabel('CDF')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Relationship between Binomial and Poisson distribution\n",
    "\n",
    "You just heard that the Poisson distribution is a limit of the Binomial distribution for rare events. This makes sense if you think about the stories. Say we do a Bernoulli trial every minute for an hour, each with a success probability of 0.1. We would do 60 trials, and the number of successes is Binomially distributed, and we would expect to get about 6 successes. This is just like the Poisson story we discussed in the video, where we get on average 6 hits on a website per hour. So, the Poisson distribution with arrival rate equal to ***np*** approximates a Binomial distribution for n Bernoulli trials with probability ***p*** of success (with ***n*** large and ***p*** small). Importantly, the Poisson distribution is often simpler to work with because it has only one parameter instead of two for the Binomial distribution.\n",
    "\n",
    "Let's explore these two distributions computationally. You will compute the mean and standard deviation of samples from a Poisson distribution with an arrival rate of 10. Then, you will compute the mean and standard deviation of samples from a Binomial distribution with parameters n and p such that ***np=10***.\n",
    "\n",
    "**Instructions**\n",
    "\n",
    "* Using the `np.random.poisson()` function, draw `10000` samples from a Poisson distribution with a mean of `10`.\n",
    "* Make a list of the `n` and `p` values to consider for the Binomial distribution. Choose `n = [20, 100, 1000]` and `p = [0.5, 0.1, 0.01]` so that ***np*** is always 10.\n",
    "* Using `np.random.binomial()` inside the provided `for` loop, draw `10000` samples from a Binomial distribution with each `n, p` pair and print the mean and standard deviation of the samples. There are 3 `n, p` pairs: `20, 0.5`, `100, 0.1`, and `1000, 0.01`. These can be accessed inside the loop as `n[i], p[i]`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Draw 10,000 samples out of Poisson distribution: samples_poisson\n",
    "samples_poisson = np.random.poisson(10, size=10_000)\n",
    "\n",
    "# Print the mean and standard deviation\n",
    "print(f'Poisson: Mean = {np.mean(samples_poisson)} Std = {np.std(samples_poisson):0.03f}')\n",
    "\n",
    "# Specify values of n and p to consider for Binomial: n, p\n",
    "n = [20, 100, 1_000, 10_000]\n",
    "p = [0.5, 0.1, 0.01, 0.001]\n",
    "\n",
    "# Draw 10,000 samples for each n,p pair: samples_binomial\n",
    "for i in range(4):\n",
    "    samples_binomial = np.random.binomial(n[i], p[i], size=10_000)\n",
    "\n",
    "    # Print results\n",
    "    print(f'n = {n[i]} Binom: Mean = {np.mean(samples_binomial)} Std = {np.std(samples_binomial):0.03f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**The means are all about the same. The standard deviation of the Binomial distribution gets closer and closer to that of the Poisson distribution as the probability p gets lower and lower.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### How many no-hitters in a season?\n",
    "\n",
    "In baseball, a no-hitter is a game in which a pitcher does not allow the other team to get a hit. This is a rare event, and since the beginning of the so-called modern era of baseball (starting in 1901), there have only been 251 of them through the 2015 season in over 200,000 games. The ECDF of the number of no-hitters in a season is shown to the right. Which probability distribution would be appropriate to describe the number of no-hitters we would expect in a given season?\n",
    "\n",
    "* ![no hitters ecdf](https://raw.githubusercontent.com/trenton3983/DataCamp/master/Images/2019-07-10_statistical_thinking_1/no_hitter_ecdf.JPG)\n",
    "\n",
    "***Note***: The no-hitter data set was scraped and calculated from the data sets available at [retrosheet.org](https://www.retrosheet.org/) ([license](https://www.retrosheet.org/notice.txt)).\n",
    "\n",
    "**Possible Answers**\n",
    "\n",
    "* ~~Discrete uniform~~\n",
    "* ~~Binomial~~\n",
    "* ~~Poisson~~\n",
    "* __**Both Binomial and Poisson, though Poisson is easier to model and compute.**__\n",
    "* ~~Both Binomial and Poisson, though Binomial is easier to model and compute.~~\n",
    "\n",
    "**With rare events (low p, high n), the Binomial distribution is Poisson. This has a single parameter, the mean number of successes per time interval, in this case the mean number of no-hitters per season.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Was 2015 anomalous?\n",
    "\n",
    "1990 and 2015 featured the most no-hitters of any season of baseball (there were seven). Given that there are on average 251/115 no-hitters per season, what is the probability of having seven or more in a season?\n",
    "\n",
    "**Instructions**\n",
    "\n",
    "* Draw `10000` samples from a Poisson distribution with a mean of `251/115` and assign to `n_nohitters`.\n",
    "* Determine how many of your samples had a result greater than or equal to `7` and assign to `n_large`.\n",
    "* Compute the probability, `p_large`, of having `7` or more no-hitters by dividing `n_large` by the total number of samples (`10000`).\n",
    "* Hit 'Submit Answer' to print the probability that you calculated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(seed=398)\n",
    "\n",
    "# Draw 10,000 samples out of Poisson distribution: n_nohitters\n",
    "n_nohitters = np.random.poisson(251/115, size=10_000)\n",
    "\n",
    "# Compute number of samples that are seven or greater: n_large\n",
    "n_large = len(n_nohitters[n_nohitters >= 7])\n",
    "\n",
    "# Compute probability of getting seven or more: p_large\n",
    "p_large = n_large/len(n_nohitters)\n",
    "\n",
    "# Print the result\n",
    "print(f'Probability of seven or more no-hitters: {p_large}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**The result is about 0.007. This means that it is not that improbable to see a 7-or-more no-hitter season in a century. There have been two in a century and a half, so it is not unreasonable.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Thinking probabilistically: Continuous variables\n",
    "\n",
    "Probability distributions of discrete variables have been covered so far. This final section will cover continuous variables, such as those that can take on any fractional value. Many of the principles are the same, but there are some subtleties. At the end of this chapter, you will be speaking the probabilistic language required to launch into the inference techniques covered in [Statistical Thinking in Python (Part 2)](https://www.datacamp.com/courses/statistical-thinking-in-python-part-2)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Probability density functions\n",
    "\n",
    "We have talked about probabilities of discrete quantities, such as die rolls and number of bus arrivals, but what about continuous quantities?  A continuous quantity can take on any value, not just discrete ones.  For example, the speed of a train can be 45.76 km/h.  Continuous variables also have probability distributions.  Let's consider an example.  In 1879, Albert Michelson performed 100 measurements of the speed of light in air.  Each measurement has some error in it; conditions, such as temperature, humidity, alignment of his optics, etc., change from measurement to measurement.  As a result, any fractional value of the measured speed of light is possible, so it's apt to describe the results with a continuous probability distribution.  In looking at Michelson's numbers, show here in units of 1000 km/s, we see this is the case.  **What probability distribution describes these data?**  I posit, these data follow the **Normal Distribution**.  To understand what the normal distribution is, lets consider its **probability density function** (PDF).  This is the continuous analog to the **probability mass function** (PMF).  ___It describes the chances of observing a  value of a continuous variable___.  The probability of observing a single value of the speed of light, does not make sense, because there is an infinity of numbers, between 299,600 and 300,100 km/s.  Instead, areas under the PDF, give probabilities.  The probability of measuring the speed of light is greater the 300,000 km/s is an area under the normal curve.  Parameterizing the PDF based on Michelson's experiments, this is about a 3% chance, since the pink region is about 3% of the total area under the PDF.  To do this calculation, we were really just looking at the cumulative distribution function (CDF), of the Normal distribution.  Here's the CDF of the Normal distribution.  Remember, the CDF gives the probability, the measured speed of light will be less than the value on the x-axis.  Reading off the value at 300,000 km/s, there is a 97% chance, the speed of light measurement, is less than that.  There's about a 3% chance it's greater.\n",
    "\n",
    "We will study the Normal distribution in more depth in the coming exercises, but for now, let's review some of the concepts we've learned about continuous distribution functions.\n",
    "\n",
    "#### Continuous Variables\n",
    "- Quantities that can take any value, not just discrete values\n",
    "\n",
    "#### Probability Density Function (PDF)\n",
    "- Continuous analog to the PMF\n",
    "- Mathematical description of the relative likelihood of observing a value of a continuous variable\n",
    "\n",
    "![area under pdf](https://raw.githubusercontent.com/trenton3983/DataCamp/master/Images/2019-07-10_statistical_thinking_1/area_under_pdf.JPG)\n",
    "\n",
    "#### Normal Cumulative Distribution Function (CDF)\n",
    "\n",
    "![cdf](https://raw.githubusercontent.com/trenton3983/DataCamp/master/Images/2019-07-10_statistical_thinking_1/cdf.JPG)\n",
    "![cdf_97](https://raw.githubusercontent.com/trenton3983/DataCamp/master/Images/2019-07-10_statistical_thinking_1/cdf_97.JPG)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv(sol_file)\n",
    "df.drop(columns=['Unnamed: 0'], inplace=True)\n",
    "df.columns = df.columns.str.strip()\n",
    "df.head(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.histplot(df['velocity of light in air (km/s)'], bins=9, kde=True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Interpreting PDFs\n",
    "\n",
    "Consider the PDF shown here. Which of the following is true?\n",
    "\n",
    "![interpret_pdf](https://raw.githubusercontent.com/trenton3983/DataCamp/master/Images/2019-07-10_statistical_thinking_1/interpret_pdf.JPG)\n",
    "\n",
    "**Instructions**\n",
    "\n",
    "* __**x is more likely than not less than 10.**__\n",
    "* ~~___x___ is more likely than not greater than 10.~~\n",
    "* ~~We cannot tell from the PDF if ___x___ is more likely to be greater than or less than 10.~~\n",
    "* ~~This is not a valid PDF because it has two peaks.~~"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Interpreting CDFs\n",
    "\n",
    "At right is the CDF corresponding to the PDF you considered in the last exercise. Using the CDF, what is the probability that **x** is greater than 10?\n",
    "\n",
    "![interpret_cdf](https://raw.githubusercontent.com/trenton3983/DataCamp/master/Images/2019-07-10_statistical_thinking_1/interpret_cdf.JPG)\n",
    "\n",
    "**Instructions**\n",
    "\n",
    "* __**0.25: Correct! The value of the CDF at x = 10 is 0.75, so the probability that x < 10 is 0.75. Thus, the probability that x > 10 is 0.25.**__\n",
    "* ~~0.75~~\n",
    "* ~~3.75~~\n",
    "* ~~15~~"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Introduction to the Normal distribution\n",
    "\n",
    "The Normal distribution is famous, and we just used it as an example to learn about continuous distributions.  It describes a continuous variable whose PDF is symmetric and has a single peak.  The Normal distribution is parameterized by two parameters.  The **mean** determines where the center of the peak is.  The **standard deviation** is a measure of how wide the peak is, or how spread out the data are.  Note, the mean and standard deviation are the names of the parameters of the Normal distribution.  **Don't confuse these with the mean and standard deviation we computed directly from the data when doing exploratory data analysis**.  The nomenclature is confusing, but it's important to keep straight.  Adding a histogram of the Michelson measurements, shows the measured speed of light in air, looks Normally distributed.  Comparing the histogram to the PDF, suffers from binning bias, therefore, it's better to compare the ECDF of the data, to the theoretical CDF of the Normal distribution.  To compute the theoretical CDF, use `np.random.normal`, to draw samples, then compute the CDF.  As was the case with sampling out of the binomial distribution, we need to provide parameters, in this case, the mean and standard deviation, to parameterize the Normal distribution we are sampling out of.  The mean and standard deviations computed from the data are good estimates, so we'll compute and pass them into `np.random.normal` to take our samples.  We then use the `def ecdf` function we already wrote, to compute the ECDFs of the data and of the Normally distributed theoretical samples we just drew.  Finally, plot the theoretical and empirical CDFs on the same plot.  With the absence of binning bias, it's much clearer that the Michelson data are approximately Normally distributed.  Now that you can sample out of the Normal distribution, let's practice using it.\n",
    "\n",
    "#### Normal distribution\n",
    "\n",
    "- Describes a continuous variable whose PDF has a single symmetric peak\n",
    "- `mean of Normal distribution ≠ mean computed from data`\n",
    "- `standard deviation of a Normal distribution ≠ standard deviation computed from data`\n",
    "\n",
    "![normal_distribution](https://raw.githubusercontent.com/trenton3983/DataCamp/master/Images/2019-07-10_statistical_thinking_1/normal_distribution.JPG)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean = np.mean(df['velocity of light in air (km/s)'])\n",
    "std = np.std(df['velocity of light in air (km/s)'])\n",
    "\n",
    "samples = np.random.normal(mean, std, size=10000)\n",
    "x, y = ecdf(df['velocity of light in air (km/s)'])\n",
    "x_theor, y_theor = ecdf(samples)\n",
    "\n",
    "sns.set()\n",
    "plt.plot(x_theor, y_theor, label='theoretical')\n",
    "plt.plot(x, y, marker='.', linestyle='none', label='measured', color='purple')\n",
    "plt.legend()\n",
    "plt.xlabel('speed of light (km/s)')\n",
    "plt.ylabel('CDF')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### The Normal PDF\n",
    "\n",
    "In this exercise, you will explore the Normal PDF and also learn a way to plot a PDF of a known distribution using hacker statistics. Specifically, you will plot a Normal PDF for various values of the variance.\n",
    "\n",
    "**Instructions**\n",
    "\n",
    "* Draw 100,000 samples from a Normal distribution that has a mean of `20` and a standard deviation of `1`. Do the same for Normal distributions with standard deviations of `3` and `10`, each still with a mean of `20`. Assign the results to `samples_std1`, `samples_std3` and `samples_std10`, respectively.\n",
    "* Plot a histograms of each of the samples; for each, use 100 bins, also using the keyword arguments `normed=True` and `histtype='step'`. The latter keyword argument makes the plot look much like the smooth theoretical PDF. You will need to make 3 `plt.hist()` calls.\n",
    "* Hit 'Submit Answer' to make a legend, showing which standard deviations you used, and show your plot! There is no need to label the axes because we have not defined what is being described by the Normal distribution; we are just looking at shapes of PDFs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Draw 100000 samples from Normal distribution with stds of interest: samples_std1, samples_std3, samples_std10\n",
    "samples_std1 = np.random.normal(20, 1, size=100000)\n",
    "samples_std3 = np.random.normal(20, 3, size=100000)\n",
    "samples_std10 = np.random.normal(20, 10, size=100000)\n",
    "\n",
    "# Make histograms\n",
    "plt.figure(figsize=(7, 7))\n",
    "for data in [samples_std1, samples_std3, samples_std10]:\n",
    "    plt.hist(data, density=True, bins=100, histtype='step')\n",
    "\n",
    "# Make a legend, set limits and show plot\n",
    "plt.legend(('std = 1', 'std = 3', 'std = 10'))\n",
    "plt.ylim(-0.01, 0.42)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**You can see how the different standard deviations result in PDFs of different widths. The peaks are all centered at the mean of 20.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### The Normal CDF\n",
    "\n",
    "Now that you have a feel for how the Normal PDF looks, let's consider its CDF. Using the samples you generated in the last exercise (in your namespace as `samples_std1`, `samples_std3`, and `samples_std10`), generate and plot the CDFs.\n",
    "\n",
    "**Instructions**\n",
    "\n",
    "* Use your `ecdf()` function to generate x and y values for CDFs: `x_std1`, `y_std1`, `x_std3`, `y_std3` and `x_std10`, `y_std10`, respectively.\n",
    "* Plot all three CDFs as dots (do not forget the `marker` and `linestyle` keyword arguments!).\n",
    "* Hit submit to make a legend, showing which standard deviations you used, and to show your plot. There is no need to label the axes because we have not defined what is being described by the Normal distribution; we are just looking at shapes of CDFs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate and Plot CDFs\n",
    "for data in [samples_std1, samples_std3, samples_std10]:\n",
    "    x_theor, y_theor = ecdf(data)\n",
    "    plt.plot(x_theor, y_theor, marker='.', linestyle='none')\n",
    "\n",
    "\n",
    "# Make a legend and show the plot\n",
    "_ = plt.legend(('std = 1', 'std = 3', 'std = 10'), loc='lower right')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**The CDFs all pass through the mean at the 50th percentile; the mean and median of a Normal distribution are equal. The width of the CDF varies with the standard deviation.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Normal distribution: Properties and warnings\n",
    "\n",
    "- The Normal distribution is very important and widely used.\n",
    "- In practice, it's used to describe most symmetric peaked data\n",
    "- For many of the statistical procedures you've heard of, Normality assumptions about the data are present.\n",
    " - It's a very powerful distribution that seems to be ubiquitous in nature, not just in the field of statistics.\n",
    "- **There are important caveats about the distribution and we need to be careful when using it.**\n",
    " 1. A dataset may not be Normally distributed, when you think it is\n",
    " 1. Another consideration, is the _lightness_ of the distribution tails\n",
    "  - If we look at the Normal distribution, the probability of being more that 4 standard deviations from the mean, is very small.\n",
    "  - ![light tails](https://raw.githubusercontent.com/trenton3983/DataCamp/master/Images/2019-07-10_statistical_thinking_1/light_tails.JPG)\n",
    "  - When you're modeling data as Normally distributed, outliers and extremely unlikely.\n",
    "  - Real datasets often have extreme values, and when this happens, the Normal distribution might not be the best description of the data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Gauss and the 10 Deutschmark banknote\n",
    "\n",
    "What are the mean and standard deviation, respectively, of the Normal distribution that was on the 10 Deutschmark banknote?\n",
    "\n",
    "![deutschmark](https://raw.githubusercontent.com/trenton3983/DataCamp/master/Images/2019-07-10_statistical_thinking_1/deutschmark.JPG)\n",
    "\n",
    "**Instructions**\n",
    "\n",
    "* **mean = 3, std = 1**\n",
    "* ~~mean = 3, std = 2~~\n",
    "* ~~mean = 0.4, std = 1~~\n",
    "* ~~mean = 0.6, std = 6~~"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Are the Belmont Stakes results Normally distributed?\n",
    "\n",
    "Since 1926, the Belmont Stakes is a 1.5 mile-long race of 3-year old thoroughbred horses. [Secretariat](https://en.wikipedia.org/wiki/Secretariat_(horse)) ran the fastest Belmont Stakes in history in 1973. While that was the fastest year, 1970 was the slowest because of unusually wet and sloppy conditions. With these two outliers removed from the data set, compute the mean and standard deviation of the Belmont winners' times. Sample out of a Normal distribution with this mean and standard deviation using the `np.random.normal()` function and plot a CDF. Overlay the ECDF from the winning Belmont times. Are these close to Normally distributed?\n",
    "\n",
    "Note: Justin scraped the data concerning the Belmont Stakes from the [Belmont Wikipedia page](https://en.wikipedia.org/wiki/Belmont_Stakes).\n",
    "\n",
    "**Instructions**\n",
    "\n",
    "* Compute mean and standard deviation of Belmont winners' times with the two outliers removed. The NumPy array `belmont_no_outliers` has these data.\n",
    "* Take 10,000 samples out of a normal distribution with this mean and standard deviation using `np.random.normal()`.\n",
    "* Compute the CDF of the theoretical samples and the ECDF of the Belmont winners' data, assigning the results to `x_theor, y_theor` and `x, y`, respectively.\n",
    "* Hit submit to plot the CDF of your samples with the ECDF, label your axes and show the plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def time_to_sec(x):\n",
    "    \"\"\"Convert time in the form 2:28.51 to seconds\"\"\"\n",
    "    time_list = x.split(':')\n",
    "    return float(time_list[0]) * 60 + float(time_list[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv(belmont_file)\n",
    "df['Time_sec'] = df['Time'].apply(time_to_sec)\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# gets the data in the same format as that used in the exercise\n",
    "d_std = df['Time_sec'].std()\n",
    "d_avg = df['Time_sec'].mean()\n",
    "\n",
    "data = df['Time_sec'][(df['Time_sec'] >= d_avg - (2.5 * d_std)) & (df['Time_sec'] <= d_avg + (2.5 * d_std))]\n",
    "belmont_no_outliers = np.array(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute mean and standard deviation: mu, sigma\n",
    "mu = np.mean(belmont_no_outliers)\n",
    "sigma = np.std(belmont_no_outliers)\n",
    "\n",
    "print(f'Mean: {mu:0.02f}\\nStandard Deviation: {sigma:0.02f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Use `np.random.normal` with **mean** and **std** to get synthetic data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sample out of a normal distribution with this mu and sigma: samples\n",
    "samples = np.random.normal(mu, sigma, size=10_000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the CDF of the samples and of the data\n",
    "x_theor, y_theor = ecdf(samples)\n",
    "x, y = ecdf(belmont_no_outliers)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the CDFs and show the plot\n",
    "_ = plt.plot(x_theor, y_theor)\n",
    "_ = plt.plot(x, y, marker='.', linestyle='none', color='g')\n",
    "_ = plt.xlabel('Belmont winning time (sec.)')\n",
    "_ = plt.ylabel('CDF')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**The theoretical CDF and the ECDF of the data suggest that the winning Belmont times are, indeed, Normally distributed. This also suggests that in the last 100 years or so, there have not been major technological or training advances that have significantly affected the speed at which horses can run this race.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### What are the chances of a horse matching or beating Secretariat's record?\n",
    "\n",
    "Assume that the Belmont winners' times are Normally distributed (with the 1970 and 1973 years removed), what is the probability that the winner of a given Belmont Stakes will run it as fast or faster than Secretariat?\n",
    "\n",
    "**Instructions**\n",
    "\n",
    "* Take 1,000,000 samples from the normal distribution using the `np.random.normal()` function. The mean `mu` and standard deviation `sigma` are already loaded into the namespace of your IPython instance.\n",
    "* Compute the fraction of samples that have a time less than or equal to Secretariat's time of 144 seconds.\n",
    "* Print the result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(seed=398)\n",
    "\n",
    "# Take a million samples out of the Normal distribution: samples\n",
    "samples = np.random.normal(mu, sigma, size=1_000_000)\n",
    "\n",
    "# Compute the fraction that are faster than 144 seconds: prob\n",
    "prob = (len(samples[samples < 144]) / len(samples))\n",
    "\n",
    "\n",
    "# Print the result\n",
    "print(f'Probability of besting Secretariat: {prob}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Great work! We had to take a million samples because the probability of a fast time is very low and we had to be sure to sample enough. We get that there is only a 0.06% chance of a horse running the Belmont as fast as Secretariat.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Exponential distribution\n",
    "\n",
    "- Just as there are many named discrete distributions, there are many named continuous distributions.\n",
    "- For example, at the bus stop in Poissonville, we know the number of buses arriving per hour, are Poisson distributed.\n",
    " -  The amount of time between arrivals of buses is Exponentially distributed.\n",
    "- The waiting time between arrivals of a Poisson process is Exponentially distributed.\n",
    "- The single parameter is the mean waiting time.\n",
    "- The distribution is not peaked\n",
    "- ![exponential PDF](https://raw.githubusercontent.com/trenton3983/DataCamp/master/Images/2019-07-10_statistical_thinking_1/exponential_pdf.JPG)\n",
    "\n",
    "#### Possible Poisson Process\n",
    "- Nuclear incidents since 1974:\n",
    " - Timing of one is, independent of all others, so the time between incidents should be Exponentially distributed.\n",
    "- We can compute and plot the CDF we would expect based on the mean time between incidents and overlay that with the ECDF from the real data\n",
    "- Take the usual approach where we draw many samples out of the Exponential distribution, using the mean inter-incident time as the parameter\n",
    "\n",
    "```python\n",
    "mean = np.mean(inter_times)\n",
    "samples = np.random.exponential(mean, size=10000)\n",
    "x, y = ecdf(inter_times)\n",
    "x_theor, y_theor = ecdf(samples)\n",
    "_ = plt.plot(x_theor, y_theor)\n",
    "_ = plt.plot(x, y, marker='.', linestyle='none')\n",
    "_ = plt.xlabel('time (days)')\n",
    "_ = plt.ylabel('CDF')\n",
    "plt.show()\n",
    "```\n",
    "- ![Exponential inter-incident times](https://raw.githubusercontent.com/trenton3983/DataCamp/master/Images/2019-07-10_statistical_thinking_1/exponential_inter-incident_times.JPG)\n",
    "- It's close to being Exponentially distributed, indicating nuclear incidents can be modeled as a Poisson process.\n",
    "- The Exponential and Normal are just two of many examples of continuous distributions.\n",
    "- In many cases, just simulated the story to get the CDF.\n",
    " - If you can simulate a story, you can get the distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Matching a story and a distribution\n",
    "\n",
    "How might we expect the time between Major League no-hitters to be distributed? Be careful here: a few exercises ago, we considered the probability distribution for the number of no-hitters in a season. Now, we are looking at the probability distribution of the time between no hitters.\n",
    "\n",
    "**Answer the question**\n",
    "\n",
    "* ~~Normal~~\n",
    "* __**Exponential**__\n",
    "* ~~Poisson~~\n",
    "* ~~Uniform~~"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Waiting for the next Secretariat\n",
    "\n",
    "Unfortunately, Justin was not alive when Secretariat ran the Belmont in 1973. Do you think he will get to see a performance like that? To answer this, you are interested in how many years you would expect to wait until you see another performance like Secretariat's. How is the waiting time until the next performance as good or better than Secretariat's distributed? Choose the best answer.\n",
    "\n",
    "**Answer the question**\n",
    "\n",
    "* ~~Normal, because the distribution of Belmont winning times are Normally distributed.~~\n",
    "* ~~Normal, because there is a most-expected waiting time, so there should be a single peak to the distribution.~~\n",
    "* ~~Exponential: It is very unlikely for a horse to be faster than Secretariat, so the distribution should decay away to zero for high waiting time.~~\n",
    "* __**Exponential: A horse as fast as Secretariat is a rare event, which can be modeled as a Poisson process, and the waiting time between arrivals of a Poisson process is Exponentially distributed.**__"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Correct! The Exponential distribution describes the waiting times between rare events, and Secretariat is rare!**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### If you have a story, you can simulate it!\n",
    "\n",
    "Sometimes, the story describing our probability distribution does not have a named distribution to go along with it. In these cases, fear not! You can always simulate it. We'll do that in this and the next exercise.\n",
    "\n",
    "In earlier exercises, we looked at the rare event of no-hitters in Major League Baseball. Hitting the cycle is another rare baseball event. When a batter hits the cycle, he gets all four kinds of hits, a single, double, triple, and home run, in a single game. Like no-hitters, this can be modeled as a Poisson process, so the time between hits of the cycle are also Exponentially distributed.\n",
    "\n",
    "How long must we wait to see both a no-hitter and then a batter hit the cycle? The idea is that we have to wait some time for the no-hitter, and then after the no-hitter, we have to wait for hitting the cycle. Stated another way, what is the total waiting time for the arrival of two different Poisson processes? The total waiting time is the time waited for the no-hitter, plus the time waited for the hitting the cycle.\n",
    "\n",
    "Now, you will write a function to sample out of the distribution described by this story.\n",
    "\n",
    "**Instructions**\n",
    "\n",
    "* Define a function with call signature `successive_poisson(tau1, tau2, size=1)` that samples the waiting time for a no-hitter and a hit of the cycle.\n",
    " * Draw waiting times `tau1` (size number of samples) for the no-hitter out of an exponential distribution and assign to `t1`.\n",
    " * Draw waiting times `tau2` (size number of samples) for hitting the cycle out of an exponential distribution and assign to `t2`.\n",
    " * The function returns the sum of the waiting times for the two events."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### def successive_poisson"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def successive_poisson(tau1, tau2, size=1):\n",
    "    \"\"\"Compute time for arrival of 2 successive Poisson processes.\"\"\"\n",
    "    # Draw samples out of first exponential distribution: t1\n",
    "    t1 = np.random.exponential(tau1, size)\n",
    "\n",
    "    # Draw samples out of second exponential distribution: t2\n",
    "    t2 = np.random.exponential(tau2, size)\n",
    "\n",
    "    return t1 + t2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Distribution of no-hitters and cycles\n",
    "\n",
    "Now, you'll use your sampling function to compute the waiting time to observe a no-hitter and hitting of the cycle. The mean waiting time for a no-hitter is 764 games, and the mean waiting time for hitting the cycle is 715 games.\n",
    "\n",
    "**Instructions**\n",
    "\n",
    "* Use your `successive_poisson()` function to draw 100,000 out of the distribution of waiting times for observing a no-hitter and a hitting of the cycle.\n",
    "* Plot the PDF of the waiting times using the step histogram technique of a previous exercise. Don't forget the necessary keyword arguments. You should use `bins=100`, `normed=True`, and `histtype='step'`.\n",
    "* Label the axes.\n",
    "* Show your plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Draw samples of waiting times: waiting_times\n",
    "waiting_times = successive_poisson(764, 715, 100_000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make the histogram\n",
    "plt.hist(waiting_times, density=True, bins=100, histtype='step')\n",
    "plt.xlabel('waiting time')\n",
    "plt.ylabel('PDF: probability of occurence')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Notice that the PDF is peaked, unlike the waiting time for a single Poisson process. For fun (and enlightenment), I encourage you to also plot the CDF.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### CDF for observing a no-hitter and a hitting of the cycle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_theor, y_theor = ecdf(waiting_times)\n",
    "plt.plot(x_theor, y_theor, marker='.', linestyle='none')\n",
    "plt.ylabel('CDF')\n",
    "plt.xlabel('x')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Certificate\n",
    "\n",
    "![certificate](https://raw.githubusercontent.com/trenton3983/DataCamp/master/Images/2019-07-10_statistical_thinking_1/2019-11-21_statistical_thinking_in_python_1_certificate.JPG)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "py312",
   "language": "python",
   "name": "py312"
  },
  "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.12.3"
  },
  "toc-autonumbering": true
 },
 "nbformat": 4,
 "nbformat_minor": 4
}