docs/exercises/08_probabilities.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Probability distributions & meteor shower gazing\n",
    "\n",
    "\n",
    "**Our goals for today:**\n",
    "- Discuss some key statistics topics: samples versus populations and empirical versus theorectical distributions\n",
    "- Simulate a head/tail coin toss and well drilling i.e. Binomial distribution\n",
    "- Simulate meteors entering Earth's atmosphere i.e. Poisson distribution\n",
    "- Simulate geomagnetic polarity reversals i.e. Gamma distribution\n",
    "- Use Gutenberg-Richter to assess earthquake probability\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup\n",
    "\n",
    "Run this cell as it is to setup your environment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from math import factorial"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Case 1: Flipping a coin\n",
    "\n",
    "Let's pretend we are flipping a coin 10 times using ```np.random.choice([0, 1])```. How many times will be get heads? 1 is heads, 0 is tails. Let's use a for loop and get Python to simulate such a coin flip scenario for us.\n",
    "\n",
    "What **for loops** do is take a chunk of code (in Python the chunk that is indented) being run multiple times. In this case, the code will get looped through 10 times -- specified by ```range(0,10)```."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# np.random.seed(0)\n",
    "for flip in range(0,10):\n",
    "    ## simulate a coin flip using np.random.choice\n",
    "    flip_result = \n",
    "    print(f\"Flip {flip} result: {flip_result}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What does `np.random.seed(0)` do? Try running the code block above without it and with it. What do you notice?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's record how many times the result was heads. We will make a list called `flip_results` and have it be blank to start. Each time we go through the code we will append the result to the list:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "flip_results = []\n",
    "\n",
    "for flip in range(0,10):\n",
    "    ## simulate a coin flip using np.random.choice\n",
    "    flip_result = \n",
    "    flip_results.append(\n",
    "    \n",
    "print(f\"{flip_results = }\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can calculate how many times were heads by taking the sum of the list:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## use sum() to calculate the number of heads\n",
    "number_heads = \n",
    "print(f\"{number_heads = }\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's flip the coin 10 times and do that 10 times. Each time we flip it, let's record how many heads resulted from the flip."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "number_heads = []\n",
    "\n",
    "## do 100 experiments\n",
    "number_experiments = 100\n",
    "number_flips = 10\n",
    "for flip_experiment in range(number_experiments):\n",
    "\n",
    "    flip_results = []\n",
    "    \n",
    "    ## each experiment, do 10 flips\n",
    "    for flip in range(0,number_flips):\n",
    "        flip_result = \n",
    "        flip_results.append(\n",
    "\n",
    "    ## count the number of heads in the current experiment using sum()\n",
    "    current_experiment = \n",
    "    \n",
    "    ## save the result\n",
    "    number_heads.append(\n",
    "\n",
    "## print the number of heads\n",
    "print(f\"{number_heads[:10] = }\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's visualize the distribution of the number of heads."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## plot the histogram of the number of heads using plt.hist\n",
    "plt.figure()\n",
    "bins = np.arange(0, number_flips+1, 1)-0.5 # make the bins centered\n",
    "plt.\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Instead of 10 experiments, let's do 1000 experiments. Plot the histogram of the result and check how the distribution changes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "number_heads = []\n",
    "\n",
    "## do 10 experiments\n",
    "number_experiments = 1000\n",
    "number_flips = 10\n",
    "for flip_experiment in range(number_experiments):\n",
    "\n",
    "    flip_results = []\n",
    "    \n",
    "    ## each experiment, do 10 flips\n",
    "    for flip in range(0,number_flips):\n",
    "        flip_result = \n",
    "        flip_results.append(\n",
    "\n",
    "    ## count the number of heads in the current experiment\n",
    "    current_experiment =\n",
    "    \n",
    "    ## save the result\n",
    "    number_heads.append(\n",
    "\n",
    "## print the number of heads\n",
    "print(f\"{number_heads[:10] = }\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## plot the histogram of the number of heads using plt.hist\n",
    "plt.figure()\n",
    "bins = np.arange(0, number_flips+1, 1) - 0.5 # make the bins centered\n",
    "plt.\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How does the distribution change as the number of experiments increases? Is the center of the distribution changing? Is the width of the distribution changing?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Binomial distribution:\n",
    "\n",
    "### Theoretical\n",
    "\n",
    "A relatively straight-forward distribution is the _binomial_ distribution which describes the probability of a particular outcome when there are only two possibilities (yes or no, heads or tails, 1 or 0).\n",
    "\n",
    "For example, in a coin toss experiment (heads or tails), if we flip the coin  $n$ times, what is the probability of getting $x$ 'heads'?  We assume that the probability $p$ of a head for any given coin toss is 50%; put another way $p$ = 0.5.  \n",
    "\n",
    "The binomial distribution can be described by an equation: \n",
    "\n",
    "$$P=f(x,p,n)= \\frac{n!}{x!(n-x)!}p^x(1-p)^{n-x}$$\n",
    "\n",
    "We can look at this kind of distribution by evaluating the probability for getting $x$ 'heads' out of $n$ attempts. We'll code the equation as a function, and calculate the probability $P$ of a particular outcome (e.g., $x$ heads in $n$ attempts). \n",
    "\n",
    "Note that for a coin toss, $p$ is 0.5, but other yes/no questions can be investigated as well (e.g., chance of finding a fossil in a sedimentary layer; whether or not a landslide occurs following an earthquake). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def binomial_probability(number_positive,positive_rate,number_attempts):\n",
    "    \"\"\"\n",
    "    This function computes the probability of getting x particular outcomes (heads) in n attempts, where p is the \n",
    "    probability of a particular outcome (head) for any given attempt (coin toss).\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    number_positive : number of a particular outcome\n",
    "    positive_rate : probability of that outcome in a given attempt\n",
    "    number_attempts : number of attempts\n",
    "    \n",
    "    Returns\n",
    "    ---------\n",
    "    prob : probability of that number of the given outcome occuring in that number of attempts\n",
    "    \"\"\"\n",
    "    x = number_positive\n",
    "    p = positive_rate\n",
    "    n = number_attempts\n",
    "\n",
    "    ## compute the binomial probability following the equation above. Hint: use the factorial() function for x! n! (n-x)!\n",
    "    prob = \n",
    "\n",
    "    return prob"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use this function to calculate the probability of getting 10 heads ($x=10$) when there are 10 coin tosses ($n=10$) given with the $p$ (probability) of 0.5.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "binomial_probability(10,0.5,10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's calculate the probability of getting [ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10] heads."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## create an array of the number of heads we want to calculate the probability for\n",
    "x = np.arange(11)\n",
    "\n",
    "prob_heads = 0.5\n",
    "number_flips = 10\n",
    "\n",
    "probabilities = []\n",
    "for number_head in x:\n",
    "    ## use the binomial_probability function to calculate the probability\n",
    "    prob = \n",
    "    probabilities.append(\n",
    "    \n",
    "    print(f\"Number of heads: {number_head}, probability: {prob:.5f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These probabilities are the theoretical probabilities of getting a particular number of heads in 10 coin flips.\n",
    "Let's verify that the histogram of the number of heads in 1000 coin flips matches the theoretical probabilities.\n",
    "\n",
    "Make a plot where you both plot the histogram from 1000 coin flips (using ```plt.hist()``` with ```density=True```) and you plot the results head_numbers probabilities (using ```plt.plot()```)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10,5))\n",
    "\n",
    "## plot the histogram from 1000 coin flips using plt.hist()\n",
    "bins = np.arange(0, number_flips+1, 1)-0.5 # make the bins centered\n",
    "plt.\n",
    "\n",
    "## plot the theoretical probabilities calculated based on the binomial distribution using plt.plot\n",
    "plt.\n",
    "\n",
    "plt.xlabel(f'Number of heads out of {number_flips} attempts')\n",
    "plt.ylabel('Fraction of times with this number of heads') \n",
    "\n",
    "plt.title(f'Coin flip results (n={number_flips})')\n",
    "\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How does the histogram of the number of heads in 1000 coin flips match the theoretical probabilities? Does this imply that the theoretical distribution i.e., the binomial distribution, is a good model for the coin flip experiment?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Empirical\n",
    "\n",
    "The type of sampling we were doing above where we were flipping coins is called a _Monte Carlo simulation_. We can use simulate data from all sorts of distributions. Let's keep focusing on the binomial distribution and look at using the ```np.random.binomial``` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "help(np.random.binomial)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`np.random.binomial( )` requires 2 parameters, $n$ and $p$, with an optional keyword argument `size` (if `size` is not specified, it returns a single trial). We could have used this function earlier to get the number of heads that were flipped, but the way we did it also worked.\n",
    "\n",
    "Let's follow the example the is given in the `np.random.binomial( )` docstring.\n",
    "\n",
    "A company drills 9 wild-cat oil exploration wells (high risk drilling in unproven areas), each with an estimated probability of success of 0.1. All nine wells fail. What is the probability of that happening? \n",
    "\n",
    "*Note that success in this context means that liquid hydocarbons came out of the well. In reality, you may not consider this a success given that the result is that more hydrocarbons will be combusted as a result, leading to higher atmospheric carbon dioxide levels and associated global warming.*\n",
    "\n",
    "If we do ```np.random.binomial(9, 0.1, 100)``` we will get a list of 100 values that represent the number of wells that yielded oil when there is a 10% (p = 0.1) chance of each individual well yielding oil."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "success_rate = 0.1\n",
    "number_wells = 9\n",
    "number_simulations = 100\n",
    "\n",
    "## simulate the number of successful wells using np.random.binomial\n",
    "number_success = \n",
    "\n",
    "print(f\"{number_success[:10] = }\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Count the number of failures\n",
    "number_failures = \n",
    "\n",
    "## Calculate the failure rate\n",
    "failure_rate = \n",
    "\n",
    "print(f\"The number of failures is {number_failures} out of {number_simulations} simulations\")\n",
    "print(f\"The failure rate is {failure_rate:.2f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can write a function that uses this process to simulate fraction of times that there no successful wells for a given number of wells, a given probability and a given number of simulations;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def wildcat_failure_rate(number_wells,success_rate,number_simulations):\n",
    "    '''\n",
    "    Simulate the number of times that there are no successful wells for a given number of wells and a given probability for each well.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    number_wells : number of wells drilled in each simulation\n",
    "    success_rate : probability that each well will be successful\n",
    "    number_simulations : number of times that drilling number_wells is simulated\n",
    "    '''\n",
    "    \n",
    "    simulations = np.random.binomial(number_wells, success_rate, number_simulations)\n",
    "    number_failures = np.sum(simulations == 0)\n",
    "    failure_rate = number_failures/number_simulations\n",
    "\n",
    "    return failure_rate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color=goldenrod>**Put the `wildcat_failure_rate` function to use**</font>\n",
    "\n",
    "Use the function to simulate the failure rate for the above scenario (9 wells drilled, 0.1 probability of success for each well) and do it for 10 simulations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "number_wells = \n",
    "success_rate = \n",
    "number_simulations = \n",
    "\n",
    "failure_rate = \n",
    "\n",
    "print(f\"The failure rate is {failure_rate:.2f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the function to simulate the failure rate for the same scenario for 1000 simulations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "number_wells = \n",
    "success_rate = \n",
    "number_simulations = \n",
    "\n",
    "failure_rate = \n",
    "\n",
    "print(f\"The failure rate is {failure_rate:.2f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the function to simulate the failure rate for 100,000 simulations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "number_wells = \n",
    "success_rate =  \n",
    "number_simulations =  \n",
    "\n",
    "failure_rate =  \n",
    "\n",
    "print(f\"The failure rate is {failure_rate:.2f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color=goldenrod>**Put the `binomial_probability` function to use**</font>\n",
    "\n",
    "Instead of the examples above we are simulating the result, we could directly use the theoretical binomial_probability distribution to calculate the probability without doing many simulations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## calculate the probability using the binomial_probability function with number_positive=0, positive_rate=0.1, number_attempts=9\n",
    "failure_rate = \n",
    "\n",
    "print(f\"The failure rate is {failure_rate:.2f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**How well does the calculated binomial_probability match the simulated wildcat_failure rates? How many times do you need to simulate the problem to get a number that matches the theoretical probability?** "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Poisson distribution:\n",
    "\n",
    "A Poisson Distribution gives the probability of a number of events in an interval generated by a Poisson process: the average time between events is known, but the exact timing of events is random. The events must be independent and may occur only one at a time.\n",
    "\n",
    "Within Earth and Planetary Science there are many processes that approximately meet this criteria.\n",
    "\n",
    "### Theoretical\n",
    "\n",
    "The Poisson distribution gives the probability that an event (with two possible outcomes) occurs $k$ number of times in an interval of time where $\\lambda$ is the expected rate of occurance. The Poisson distribution is the limit of the binomial distribution for large $n$. So if you take the limit of the binomial distribution as $n \\rightarrow \\infty$ you'll get the Poisson distribution:\n",
    "\n",
    "$$P(k) = e^{-\\lambda}\\frac{\\lambda^{k}}{k!}$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "def poisson_probability(k,lamb):\n",
    "    \"\"\"\n",
    "    This function computes the probability of getting k particular outcomes when the expected rate is lambda.\n",
    "    \"\"\"\n",
    "    \n",
    "    ## compute the poisson probability of getting k outcomes when the expected rate is lambda using the equation above\n",
    "    prob = \n",
    "    \n",
    "    return prob"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Observing meteors\n",
    "\n",
    "<img src=\"./images/AMS_TERMINOLOGY.png\" width = 600>\n",
    "\n",
    "From https://www.amsmeteors.org/meteor-showers/meteor-faq/:\n",
    "\n",
    "> **How big are most meteoroids? How fast do they travel?** The majority of visible meteors are caused by particles ranging in size from about that of a small pebble down to a grain of sand, and generally weigh less than 1-2 grams. Those of asteroid origin can be composed of dense stony or metallic material (the minority) while those of cometary origin (the majority) have low densities and are composed of a “fluffy” conglomerate of material, frequently called a “dustball.” The brilliant flash of light from a meteor is not caused so much by the meteoroid’s mass, but by its high level of kinetic energy as it collides with the atmosphere.\n",
    "\n",
    "> Meteors enter the atmosphere at speeds ranging from 11 km/sec (25,000 mph), to 72 km/sec (160,000 mph!). When the meteoroid collides with air molecules, its high level of kinetic energy rapidly ionizes and excites a long, thin column of atmospheric atoms along the meteoroid’s path, creating a flash of light visible from the ground below. This column, or meteor trail, is usually less than 1 meter in diameter, but will be tens of kilometers long.\n",
    "\n",
    "> The wide range in meteoroid speeds is caused partly by the fact that the Earth itself is traveling at about 30 km/sec (67,000 mph) as it revolves around the sun. On the evening side, or trailing edge of the Earth, meteoroids must catch up to the earth’s atmosphere to cause a meteor, and tend to be slow. On the morning side, or leading edge of the earth, meteoroids can collide head-on with the atmosphere and tend to be fast.\n",
    "\n",
    "> **What is a meteor shower? Does a shower occur “all at once” or over a period of time?** Most meteor showers have their origins with comets. Each time a comet swings by the sun, it produces copious amounts of meteoroid sized particles which will eventually spread out along the entire orbit of the comet to form a meteoroid “stream.” If the Earth’s orbit and the comet’s orbit intersect at some point, then the Earth will pass through this stream for a few days at roughly the same time each year, encountering a meteor shower. The only major shower clearly shown to be non-cometary is the Geminid shower, which share an orbit with the asteroid (3200 Phaethon): one that comes unusually close to the sun as well as passing through the earth’s orbit. Most shower meteoroids appear to be “fluffy”, but the Geminids are much more durable as might be expected from asteroid fragments.\n",
    "\n",
    "## Observing the Southern Taurids meteor shower\n",
    "\n",
    "Let's say you are planning a camping trip to go out try to see shooting stars next Fall in a rural location. You are looking at a date in October and that there is an active shower:\n",
    "\n",
    "> **Southern Taurids**\n",
    "\n",
    "> *Active from September 28th to December 2, 2021. The peak is November 4-5, 2021*\n",
    "\n",
    "> The Southern Taurids are a long-lasting shower that reaches a barely noticeable maximum on October 9 or 10. The shower is active for more than two months but rarely produces more than five shower members per hour, even at maximum activity. The Taurids (both branches) are rich in fireballs and are often responsible for increased number of fireball reports from September through November. https://www.amsmeteors.org/meteor-showers/meteor-shower-calendar/\n",
    "\n",
    "At a rate of 5 observed meteors per hour, what is the probability of observing 6 meteors in an hour?\n",
    "\n",
    "We can use the Poisson probability function to answer this question:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lamb = 5\n",
    "k = 6\n",
    "## use the poisson_probability function to calculate the probability\n",
    "prob = \n",
    "print(f\"The probability of getting {k} outcomes when the expected rate is {lamb} is {prob*100:.3f}%\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So that result tells us that there is a 14.6% chance of observing exactly 6, but it would be much more helpful to be able to visualize the probability distribution. So let's go through and calculate the probability of seeing any number between 0 and 10. First, we can make an array between 0 and 11:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Let's calcualte the theoretical probability of seeing different numbers of meteors when the expected rate is 5 (Southern Taurids)\n",
    "taurid_meteor_rate = 5\n",
    "number_meteors_seen = np.arange(0,11)\n",
    "taurid_meteor_sighting_probability = []\n",
    "\n",
    "for n in number_meteors_seen:\n",
    "    ## use the poisson_probability function to calculate the probability\n",
    "    prob = \n",
    "    taurid_meteor_sighting_probability.append(\n",
    "    print(f\"The probability to see {n} meteors is {prob:.3f}\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Based on the probability of seeing 0 meteors, what is the theoretical probability of seeing at least 1 meteor?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## use the poisson_probability function to calculate the probability of seeing at least 1 meteor. Hint: use 1 - the probability of seeing 0 meteors\n",
    "prob_w_meteor_sighting = \n",
    "print(f\"The probability of seeing at least 1 meteor is {prob_w_meteor_sighting*100:.3f}% at the Southern Taurids rate of {taurid_meteor_rate} meteors per hour\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Plot the probability distribution using plt.plot\n",
    "plt.figure()\n",
    "plt.\n",
    "plt.xlabel('Number of meteors seen')\n",
    "plt.ylabel('Probability')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When there is not an active shower the background meteor rate is about 2 an hour (although it is variable depending on time of night and season; see more here: https://www.amsmeteors.org/meteor-showers/meteor-faq/).\n",
    "\n",
    "Let's update our code to calculate the probability of seeing different numbers of meteors when the background rate is 2 an hour (lambda = 2).\n",
    "\n",
    "- **Calculate the probability of seeing different numbers of meteors when the background rate is 2 an hour (lambda = 2).**\n",
    "- **Plot these probabilities alongside the probability of seeing those same numbers during the Southern Taurids shower.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Let's calcualte the theoretical probability of seeing different numbers of meteors when the background rate is 2 an hour\n",
    "\n",
    "background_meteor_rate = 2\n",
    "number_meteors_seen = np.arange(0,11)\n",
    "background_sighting_probability = []\n",
    "\n",
    "for n in number_meteors_seen:\n",
    "    ## use the poisson_probability function to calculate the probability\n",
    "    prob = \n",
    "    background_sighting_probability.append(\n",
    "\n",
    "    print(f\"The probability to see {n} meteors is {prob:.3f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Based on the probability of seeing 0 meteors, what is the theoretical probability of seeing at least 1 meteor?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## use the poisson_probability function to calculate the probability of seeing at least 1 meteor. Hint: use 1 - the probability of seeing 0 meteors\n",
    "prob_w_meteor_sighting = \n",
    "print(f\"The probability of seeing at least 1 meteor is {prob_w_meteor_sighting*100:.3f}% at the background rate of {background_meteor_rate} meteors per hour\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Plot the probability distribution using plt.plot\n",
    "plt.figure()\n",
    "## plot the Southern Taurids probability distribution\n",
    "plt.\n",
    "## plot the background probability distribution\n",
    "plt.\n",
    "\n",
    "plt.ylabel('Probability')\n",
    "plt.xlabel('Number of meteors seen')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulate meteor observing\n",
    "\n",
    "There are many cases where it can be useful to simulate data sets. In this case, one could simulate what your experience could be in terms of the number of hours you could spend looking at the night sky and seeing 1 meteor or more on a normal night vs. a night with the Southern Taurids shower ongoing.\n",
    "\n",
    "We can use the `np.random.poisson` function to simulate 'realistic' data. \n",
    "\n",
    "Recall the Poisson distribution represents the probability of observing $k$ events when the expected rate is $\\lambda$. \n",
    "\n",
    "The Poisson distribution assumes that the events are independent and occur at a constant rate. Each call to `np.random.poisson( )` is doing an experiment following the Poisson distribution to see how many events occur.\n",
    "\n",
    "Let's try it with $\\lambda = 2$ (the background rate) and watch the sky for 1000 hours."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lamb = 2\n",
    "number_hours = 1000\n",
    "number_meteors = []\n",
    "\n",
    "for hour in np.arange(number_hours):\n",
    "    ## use the np.random.poisson function to simulate the number of meteors per hour\n",
    "    number_meteor_per_hour = \n",
    "    number_meteors.append(\n",
    "\n",
    "## Calculate the fraction of hours watched with at least 1 meteor sighting\n",
    "number_meteors = np.array(number_meteors)\n",
    "## use np.sum to calculate the number of hours with at least 1 meteor sighting\n",
    "number_hours_w_meteor_sighting = \n",
    "prob_w_meteor_sighting = \n",
    "\n",
    "print(f\"The probability of seeing at least 1 meteor per hour is {prob_w_meteor_sighting*100:.3f}% from {number_hours} hours of watching\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- **Do the same meteor watching simulation with $\\lambda = 5$ (the Southern Taurids rate).**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lamb = 5\n",
    "number_hours = 1000\n",
    "number_meteors = []\n",
    "\n",
    "for hour in np.arange(number_hours):\n",
    "    ## use the np.random.poisson function to simulate the number of meteors per hour\n",
    "    number_meteor_per_hour = \n",
    "    number_meteors.append(\n",
    "    \n",
    "## Calculate the fraction of hours watched with at least 1 meteor sighting\n",
    "number_meteors = np.array(number_meteors)\n",
    "## use np.sum to calculate the number of hours with at least 1 meteor sighting\n",
    "number_hours_w_meteor_sighting = \n",
    "prob_w_meteor_sighting = \n",
    "\n",
    "print(f\"The probability of seeing at least 1 meteor per hour is {prob_w_meteor_sighting*100:.3f}% from {number_hours} hours of watching\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above estimate is based on Monte Carlo simulation of hundreds and thousands of hours. How does it compare to the theoretical probability?\n",
    "\n",
    "Hint: If your results are not similar, try increasing the number of hours simulated. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extension: Earthquake Probability\n",
    "\n",
    "The occurrence of earthquakes is also a Poisson process, events occur randomly in time, and the average recurrence can be determined from Gutenberg-Richter. \n",
    "\n",
    "The Gutenberg-Richter statistic which gives the annual rate of earthquakes above a given magnitude.\n",
    "\n",
    "Assumeing the Gutenberg-Richter relationshipfor the San Franciso Bay Area: log10(N)= 3.266 - 0.797M, for each year.\n",
    "\n",
    "Let's apply the Poisson distribution to this problem. \n",
    "\n",
    "So $\\lambda = N * {\\Delta}time$, where N is the annual rate. It is common to consider ${\\Delta}time=30 yrs$.\n",
    "\n",
    "Use the Poisson's distribution to find the probability of 0 to 9 M5 events in a 30 year period."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Calculate the expected number of M5 earthquakes in a year\n",
    "M = 5\n",
    "logN = \n",
    "N = \n",
    "print(f\"The expected number of M5 earthquakes in 1 year is {N:.2f}\")\n",
    "\n",
    "## Calculate the expected number of M5 earthquakes in 30 years\n",
    "lamb = \n",
    "print(f\"The expected number of M5 earthquakes in 30 years is {lamb:.2f}\")\n",
    "\n",
    "## Calculate the probability of observing 0-9 M5 earthquake in a year\n",
    "number_earthquakes = np.arange(0,10)\n",
    "earthquake_probability = []\n",
    "for n in number_earthquakes:\n",
    "    ## use the poisson_probability function to calculate the probability\n",
    "    prob = \n",
    "    earthquake_probability.append(\n",
    "    print(f\"The probability of observing {n} M5 earthquakes in a year is {prob*100:.3f}%\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Plot the probability distribution using plt.plot\n",
    "plt.figure()\n",
    "plt.\n",
    "\n",
    "plt.ylabel('Probability')\n",
    "plt.xlabel('Number of M5 earthquakes in 30 years')\n",
    "plt.legend()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How does the probability change with the number of events? How does that reconcile with the rate of M5 earthquakes in 30 years?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Poisson Probability of 1 or more earthquakes\n",
    "\n",
    "The Poisson probability of zero events has an interesting use in characterizing earthquake hazard.\n",
    "\n",
    "$P(k=0)=e^{-\\lambda}$\n",
    "\n",
    "The complement of the zero event probability is the probability of 1 or more earthquakes occuring in the period of time. It is this probability that is used in earthquake forecast reports. The probability of one or more events is written as; \n",
    "\n",
    "$P(k >= 1) = 1 - e^{-\\lambda}$\n",
    "\n",
    "Determine the probability of 1 or more M4, M5, M6 and M7 in a 30 year period.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Calculate the probability of observing at least one M4, M5, M6 and M7 earthquake in 30 years\n",
    "magnitudes = [4, 5, 6, 7]\n",
    "number_years = 30\n",
    "for M in magnitudes:\n",
    "    logN = \n",
    "    N = \n",
    "    lamb = \n",
    "    prob0 = \n",
    "\n",
    "    ## calculate the probability of observing at least one event, hint: use 1 - prob0\n",
    "    prob = \n",
    "    print(f\"The probability of observing at least one M{M} earthquake in {number_years} years is {prob*100:.3f}%\") "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How do the probabilities change if a 10 year period is considered?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Calculate the probability of observing at least one M4, M5, M6 and M7 earthquake in 10 years\n",
    "magnitudes = [4, 5, 6, 7]\n",
    "number_years = 10\n",
    "for M in magnitudes:\n",
    "    logN = \n",
    "    N = \n",
    "    lamb = \n",
    "    prob0 = \n",
    "\n",
    "    ## calculate the probability of observing at least one event, hint: use 1 - prob0\n",
    "    prob = \n",
    "    print(f\"The probability of observing at least one M{M} earthquake in {number_years} years is {prob*100:.3f}%\") \n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}