day0/ProbabilityAndData.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a8a42a5-a8c3-4b89-b37b-814be6e02995",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.optimize import minimize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de306753-60d4-421f-85cc-0beb54157d90",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "%matplotlib widget"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11dec9f5-7820-4e57-b104-43f30faeaaa1",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "# Probability and Data\n",
    "\n",
    "##### Version 0.1\n",
    "\n",
    "***\n",
    "\n",
    "By AA Miller (Northwestern/CIERA)  \n",
    "08 September 2023"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da475079-efce-4666-befa-38305fac2359",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "\"There is a 40% chance of rain between 2 and 4 pm?\"\n",
    "\n",
    "What does this mean?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4924a2f8-84cf-4823-9086-71e5aaf4ddf4",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "*take a minute to discuss this with your partner*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f0d6c34a-4f3c-4496-b19a-0e09d22c72e0",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "\"There is a 40% chance of rain between 2 and 4 pm?\"\n",
    "\n",
    "What does this mean?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba77ebb8-8901-49f3-aa39-e4ee47dd8e1f",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "source": [
    "$~~~~$Looking at the probability distribution for rain, 40% is the mean outcome"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2b8f069-f45e-42ea-bb15-d71af7598493",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "source": [
    "$~~~~$40% is the most likely outcome"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "41a55edf-73a5-490d-b6dd-fff58de99d11",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "source": [
    "$~~~~$Given 1000s of days with identical conditions, 40% of those days have rain"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f9bc260-a78f-491d-b2a3-896978ffaa25",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "source": [
    "$~~~~$40% indicates the belief of the weather center that it would rain"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d03a2fb8-01de-4c11-a3eb-60891b54d9cf",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## Probability"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9342ab01-870b-4a11-bcfa-7f46b6450f77",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "An **experiment** is a real occurance that can be observed.\n",
    "\n",
    "Experiments have **outcomes**."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7913bfbd-9d4e-4e31-a89c-aa7b3291af08",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "all possible outcomes $\\equiv$ sample space, $\\mathcal{S}$\n",
    "\n",
    "a single outcome $\\equiv$ sample point\n",
    "\n",
    "subset of $\\mathcal{S} \\equiv$ event"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2e881c45-3f8f-4124-866b-da285225abd2",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "source": [
    "<img style=\"display: block; margin-left: auto; margin-right: auto\" src=\"images/sampleSpace.png\" alt=\"sample space\" width=\"650\" align=\"middle\">"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21a67f91-36c2-481b-a593-7304160d2bd7",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "All measurements have uncertainty$^\\dagger$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "86cdb950-1f56-4364-aa77-dcfba69c3e5c",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "$^\\dagger$Fundamentally this is set by the different uncertainty principles in quantum mechanics, but this is also seen at the macro level (e.g., repeated measurements of the same experiment with a single detector will lead to different outcomes)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a5d2b7d-8638-47c5-93dd-07288c7e7e7e",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "A system characteristic is **random** when it is not know with complete certainty.\n",
    "\n",
    "Randomness can be quantified via **probability**."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "82d6eabf-68c7-4aa5-901a-cb078964f326",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## Definitions of Probability"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "471fc8da-5f5d-4f11-b3e1-90c9aa47d98d",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "$\\mathcal{S} \\equiv$ sample space\n",
    "\n",
    "$A \\equiv$ subset of $\\mathcal{S}$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2387002c-dbf0-45bd-ae0f-2af9b8152de0",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "$P(A) \\equiv$ probability"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9e050263-6c3c-48a1-b952-0236254cc9ad",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "source": [
    "$~~~~$(i) for all $A$, $P(A) \\ge 0$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8fe384d5-618d-462f-93f6-04819d380f1f",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "source": [
    "$~~~~$(ii) if $A$ and $B$ are disjoint (i.e., $A \\cap B = \\emptyset$);  \n",
    "$~~~~~~~~~~~~$then $P(A \\cup B) = P(A) + P(B)$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c3ac7f22-16c7-47b2-9399-68ff0cc9ebe2",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "source": [
    "$~~~~$(iii) P(S) = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c4f8d81-d3cd-4243-ac65-85a68171be91",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "*Note* –  \n",
    "$\\cap$ means \"intersection\"  \n",
    "$\\cup$ means \"union\"  \n",
    "$\\subset$ means \"subset\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1467406f-bb0e-4aa2-8ba1-e099ac3c1bfc",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "From the three axioms of probability:  \n",
    "$~~~~P(\\bar{A}) = 1 - P(A)$  ($\\bar{A}$ is the complement of A  \n",
    "$~~~~P(A \\cup \\bar{A}) = 1$  \n",
    "$~~~~0 \\le P(A) \\le 1$  \n",
    "$~~~~P(\\emptyset) = 0$  \n",
    "$~~~~$if $A \\subset B$ then $P(A) \\le P(B)$  \n",
    "$~~~~P(A \\cup B) = P(A) + P(B) - P(A \\cap B)$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d0aaa30-c4a3-45df-9513-0837f742c7af",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "source": [
    "<img style=\"display: block; margin-left: auto; margin-right: auto\" src=\"images/conditionalProbability.png\" width=\"650\" align=\"middle\">"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a460be6a-b594-4ff4-9f1c-07816a4647a6",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "### Conditional Probability"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9bc351ed-11f6-4f2c-bb64-eab41a40d92d",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "$$P(A|B) = \\frac{P(A \\cap B)}{P(B)}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bb9e661b-467c-4744-9c07-d1e83c8e27a4",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "source": [
    "<img style=\"display: block; margin-left: auto; margin-right: auto\" src=\"images/conditionalProbability.png\" width=\"650\" align=\"middle\">"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "253eedfd-e14e-4f5a-9784-6664993de50b",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "If $P(A \\cap B) = P(A)P(B) \\Rightarrow A \\& B$ are independent  \n",
    "$~~~~ \\Rightarrow P(A|B) = P(A)$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d94910d9-63d7-485b-8d1a-744062d4733d",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "Fair 6-sided die as an example: $\\mathcal{S} = \\{1,2,3,4,5,6\\}$  \n",
    "$~~~~$(i) $P(\\mathrm{roll} > 3)$  \n",
    "$~~~~~~~~~~~~ A = \\{4, 5, 6\\}$  \n",
    "$~~~~~~~~~~~~ P(A) = P(4) + P(5) + P(6) = 3 \\times 1/6 = 1/2$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f9db699e-2665-4f8c-bffd-3b045836cc58",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "$~~~~$(ii) $P(\\mathrm{roll} > 3|\\mathrm{odd\\;number\\;rolled)}$  \n",
    "$~~~~~~~~~~~~ B = \\{1, 3, 5\\}$  \n",
    "$~~~~~~~~~~~~ P(B) = 1/2$  \n",
    "$~~~~~~~~~~~~ P(A \\cap B) = P(5) = 1/6$  \n",
    "$~~~~~~~~~~~~ P(A|B) = \\frac{P(A \\cap B)}{P(B)} = \\frac{1/6}{1/2} = 1/3$  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9fc0fd02-592d-40e5-8866-9230dccb7e9b",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "$A \\cap B$ is the same as $B \\cap A$  \n",
    "$P(B \\cap A) = P(A | B)P(B) = P(B | A)P(A)$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a6f0ae3-2ae0-4078-8381-1b7ad878a527",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "source": [
    "$~~~~~$\n",
    "$$ P(A | B) = \\frac{P(B | A)P(A)}{P(B)}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c1c4b9f3-2406-422e-960c-b64b6cd95341",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "(I generally find it important to emphasize that Bayes' Law *is not wrong* - it is a natural consequence of conditional probability. There are some, within our field and elsewhere, that do not like a Bayesian interpretation of probability. There is a great deal of nuance in arguments about Bayes/subjective statistics and Bryan will touch on this more in the following lecture.)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c63e2fab-f691-4e8b-b4dc-6557cda832a6",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## Measurements and models"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "47698ade-69ba-42f1-952c-9149b3b055a9",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "Physics is just a collection of models to understand the Universe."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed83e706-d81d-4f55-ab10-11b86939e36e",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "We obtain outcomes from experiments to learn about the behavior of the Universe. \n",
    "\n",
    "These measurements have associated uncertainties, which can be described via probabilty. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d1b50ee-c1f0-4a87-b87c-de20bbc02ba8",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "If out sample space $\\mathcal{S}$ has events $x_i$ for $i = 1, 2, 3, ..., N$,  \n",
    "then the probability to observe $x_i = P(x_i) \\equiv f_i$, and\n",
    "\n",
    "$$\\sum_{i=1}^\\infty f_i = 1$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2a3161f-9abc-40a3-ae4e-9656f798dbb6",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "For an experiment with single continuous variable $X$ that can take on values $x$,  \n",
    "what is $P(x)$?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b1c9479-8df2-4bb5-90b8-b5828f484998",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "source": [
    "$$P(x) = 0$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18124289-2ee3-4f2a-8a16-2ab06352cd0d",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "We need to consider probability within a small finite range: $[x, x+dx]$,  \n",
    "which is known as the **probability density function** (p.d.f.).\n",
    "\n",
    "$$P([x, x+dx]) = f(x)dx$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c7af12e-a88a-4130-a9d8-a2ce316dc1b4",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "By definition all p.d.f.s are normalized:  \n",
    "$$\\int_\\mathcal{S} f(x)dx = 1$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "88de641d-3772-491a-8797-ce90d6098dc3",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "For known uncertainties we can ask: given some model $\\mathcal{M}$, which depends upon model parameters $\\theta$, what is the probability of our observation/measurement, $x$?\n",
    "\n",
    "$$p(x|\\mathcal{M}(\\theta))$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b5bbe6e-c87b-40ac-a909-fb84bebdae63",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "In the special, and oft-assumed, case where the probability distribution for $x$ is Gaussian, then: \n",
    "\n",
    "$$p(x|\\mathcal{M}(\\theta)) = \\frac{1}{\\sqrt{(2\\pi\\sigma^2)}} \\exp\\left(- \\frac{(x - \\mu)^2}{2\\sigma^2}\\right),$$\n",
    "\n",
    "where $\\theta = (\\mu, \\sigma^2$) and $\\mu$ is the mean and $\\sigma^2$ is the variance of the Gaussian."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15ac1938-d80f-4dbe-94cc-74a911a50343",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "If individual observations are independent, (i.e., there is no correlation in successive observations) then we can define the likelihood $\\mathcal{L}$ of all observations $x_i$ as the product of the individual observation probabilities: \n",
    "\n",
    "$$\\mathcal{L} \\equiv \\prod_{i=1}^{N} p(x_i|\\mathcal{M}(\\theta))$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ced57c41-834b-453c-ad3b-dcfd7416afaa",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "Note that $\\mathcal{L}$ is not a true p.d.f., it is the product of several probabilities drawn from a p.d.f. $\\mathcal{L}$ is not normalized and often is a very small number."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e5aac7b0-a8fb-434e-99a5-be2827a27254",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "Again, for the special case where each $x_i$ follows a Gaussian p.d.f., then\n",
    "\n",
    "$$\\mathcal{L} = p(x_i|\\mathcal{M}(\\theta)) = \\prod_{i=1}^N \\frac{1}{\\sqrt{(2\\pi\\sigma_i^2)}} \\exp\\left(- \\frac{(x_i - \\mu)^2}{2\\sigma_i^2}\\right),$$\n",
    "\n",
    "where $\\sigma_i$ represents heteroskedastic uncertainties on the individual observations $x_i$. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9e27ca71-ead6-4e1d-adf8-af3524d914f2",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "Because $\\mathcal{L}$ is typically very small, numerically the log of the likelihood is more computationally stable:\n",
    "\n",
    "$$\\ln \\mathcal{L} =  -\\frac{N}{2} \\ln\\left(\\frac{1}{2\\pi}\\right) - \\sum_{i=1}^{N} \\ln \\sigma_i - \\sum_{i=1}^N \\left(\\frac{(x_i - \\mu)^2}{2\\sigma_i^2}\\right).$$\n",
    "\n",
    "The first two terms in this equation are constant, meaning: \n",
    "\n",
    "$$\\ln \\mathcal{L} \\propto - \\sum_{i=1}^N \\left(\\frac{(x_i - \\mu)^2}{2\\sigma_i^2}\\right).$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e18a117c-a4b0-4bec-88b4-96ca8c8730c5",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "Thus, as we try to determine which model (among many) \"best\" describes the Universe, the answer is that the best model **maximizes $\\mathcal{L}$**."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "34504b25-6c22-46cb-931e-7aff88837960",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "(maximizing the likelihood and the log of the likelihood is equivalent)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "07ba97c4-7074-4709-a16e-a7f564001f65",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "For a collection of outcomes ${x_i}$ (*note* – following an experiment the outcome/data  are fixed and do not change), the likelihood is a function of the model parameters $\\theta$.\n",
    "\n",
    "For any model $\\mathcal{M}$, $\\mathcal{L}$ is maximized via the choice of $\\theta$ that \"best\" explain the data. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03c3a939-a2dd-46fa-92f7-38e778f9259d",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## Frequentist Interpretation of Models and Probability"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f1549306-cc1a-4f61-bc06-0eb06cc4ca64",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "There are two fundamental approaches to statistical inference: frequentist and Bayesian statistics. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d7b5b7e-2e74-4055-9f89-4d339db094d3",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "Philosophically these approaches differ, but I want to emphasize that neither is wrong. \n",
    "\n",
    "As is true for most everything in data science, the particulars of the problem determine which approach is best. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2e0230ad-ee25-4e06-a172-978c33d07c6b",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "There are three tenets for frequentist statistics:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2ef2e34-129a-40dc-80c5-1de36d5ef723",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "source": [
    "$~~~~~~$ Probabilities are relative frequencies of events."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2017efa-f393-43cb-b032-6009ba25ec9e",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "source": [
    "$~~~~~~$ Parameters are fixed, unknown constants. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1db3e94a-b87e-4f69-b2dc-2a1b838c3a09",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "source": [
    "$~~~~~~$ Statistial procedures should have well-defined long run frequency properties."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8bcc5b70-972b-4eb5-94df-18fd04218ffb",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "(Bryan is going to discuss Bayesian interpretations in more detail during the next lecture)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "edc97efa-cccf-4ce9-977f-b819f4713f91",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## Maximum Likelihood Estimation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b5278b5b-d3c9-4958-a98f-8d592908056b",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "For simple $\\mathcal{M}$ one can differentiate the likelihood and determine its maximum analytically (but $\\mathcal{M}$ is rarely this simple). \n",
    "\n",
    "Numerical methods can be used as an alternative to analytic solutions."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d36f41bc-e616-4530-8d3b-d15ae08dae67",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "Maximizing the likelihood is equivalent to minimizing the negative log likelihood. \n",
    "\n",
    "The [`minimize`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html?highlight=minimize) function inside the `scipy.optimize` module provides a \"standard\" for determining the maximum likelihood."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "75f7e13b-a088-4b00-8dfb-bf874350fcb0",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "(the best way to minimize a function is by taking derivatives, but sometimes that is not possible. There are multiple numerical methods to address this, many of them options for the `minimize` function. The [BFGS](https://en.wikipedia.org/wiki/Broyden%E2%80%93Fletcher%E2%80%93Goldfarb%E2%80%93Shanno_algorithm) algorithm is the default.)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a3c4727-5ffd-41b8-9369-5e18970018d9",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "A caveat about numerical MLE methods: there are no guarantees of convergence to a *global* minimum. \n",
    "\n",
    "There is always a risk that the algorithm only found a local minimum."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7a4b3ad-748d-4fa3-85a2-c3389ef15aab",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "One nice aspect of the Gaussian likelihood, is that we can recast $\\ln \\mathcal{L}$ in terms of the $z$-score, $z_i = (x_i - \\mu)/\\sigma_i$: \n",
    "\n",
    "$$\\ln \\mathcal{L} = \\mathrm{constant} - \\frac{1}{2} \\sum z_i^2$$\n",
    "\n",
    "This sum is often referred to as the \"$\\chi^2$\" ($\\chi^2 = \\sum z_i^2$). "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa3d14f8-7868-4b5e-8a25-7744e5bfbc76",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "(if you have ever heard that the \"best\" model minimizes the $\\chi^2$, this is where it comes from: maximizing a Gaussian likelihood)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f9223bf2-2556-405a-8650-e2abab10515e",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "This means the distribution of $\\ln \\mathcal{L}$ can be determined from the $\\chi^2$ distribution with $N - k$ degrees of freedom, where $k$ is the number of model parameters. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9202ebb2-4988-42d2-8439-941c9247b6fa",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "(briefly, the sum of the squares of the $z$-score, $Q = \\sum z_i^2$, follows a $\\chi^2$ distribution with $k = N$ degrees of freedom, \n",
    "\n",
    "$$p(Q|k) \\equiv \\chi^2(Q|k) = \\frac{1}{2^{k/2}\\Gamma(k/2)}Q^{k/2 -1} \\exp(-Q/2),$$\n",
    "\n",
    "where $Q$ must be greater than zero and $\\Gamma$ is the [gamma function](https://en.wikipedia.org/wiki/Gamma_function).)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6cba4594-2f7b-4c19-a892-00dcf4bdb593",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "Comparing the $\\chi^2$ score to a $\\chi^2$ distribution with $N - k$ degrees of freedom provides a numerical measure of the \"goodness\" the model."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b8532ea-d6da-421d-8529-aaa7525e0b5e",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "In the frequentist interpretation of probability MLE is used to determine optimal estimates of the model parameters $\\theta$. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a7143ad-b01a-4a30-876d-eccdddb78b7e",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "(you will now use `scipy` to determine optimal model parameters $\\theta$)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a15f796-cfe6-4f0e-b801-2ac9ae603642",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## Problem 0) Helper Functions\n",
    "\n",
    "We will need to make lots of plots of data and models throughout this notebook."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6fe50a74-e97d-4dc9-beed-a7ce9ebc71b5",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "**Problem 0a**\n",
    "\n",
    "Write a function `data_plot()` that creates and returns `matplotlib` figure and axes instances using [`plt.subplots`](https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.subplots.html).\n",
    "\n",
    "The function should label the abcissa `x` and the ordinate `y`, and take optional arguments `x_obs`, `y_obs`, `y_obs_unc` that are plotted when provided by the user. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ae99f62-518c-44ed-b57b-abdd90d24f65",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def data_plot(x_obs=[], y_obs=[], y_obs_unc=[]):\n",
    "    \n",
    "    fig, ax = plt.subplots(figsize=(10,5.5))\n",
    "    \n",
    "    if len(x_obs) > 0 and len(y_obs) > 0:\n",
    "        if len(y_obs_unc) > 0:\n",
    "            ax.errorbar(x_obs, y_obs, y_obs_unc, \n",
    "                        fmt='o', mec='RebeccaPurple', mfc='white', \n",
    "                        ecolor='RebeccaPurple')\n",
    "        else:\n",
    "            ax.errorbar(x_obs, y_obs, \n",
    "                        fmt='o', mec='RebeccaPurple', mfc='white')\n",
    "    \n",
    "    ax.set_xlabel('x', fontsize=16)\n",
    "    ax.set_ylabel('y', fontsize=16)\n",
    "    ax.tick_params('both', labelsize=13)\n",
    "    return fig, ax"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d874de2-daee-4321-b54d-26cd1ca1b9fc",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "**Problem 0b**\n",
    "\n",
    "Using `np.random` simulate 15 observations that are drawn from a linear relation defined by $f(x) = 3.14\\,x + 6.626$. The observations should be collected over the range [0,10]. \n",
    "\n",
    "Assume that the observations are noisy, and the scatter is described by a Gaussian with variance = 9. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ba514c7e-e068-42da-acd4-2a5c4abe5ddc",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "n_obs = 15\n",
    "rng = np.random.default_rng(seed=2009)\n",
    "\n",
    "x_obs = rng.uniform(0, 10, n_obs)\n",
    "y_obs = rng.normal(3.14*x_obs + 6.626, 3)\n",
    "y_obs_unc = np.ones_like(y_obs)*3"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "242988bd-2717-4376-9663-b664c723b46b",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "**Problem 0c**\n",
    "\n",
    "Confirm your results from the previous two problems by using `data_plot()` to display the observations generated in **0b**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2c5e2b29-0bad-4327-90dd-1fdaad967c97",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, ax = data_plot( # complete\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8836b58b-09ab-4cd5-bd7e-35ca3b1adb95",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## Problem 1) The likelihood\n",
    "\n",
    "We just showed that the likelihood is related to the observations conditioned on the model:\n",
    "\n",
    "$$\\mathcal{L} \\equiv \\prod p(x_i|\\mathcal{M}(\\theta))$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5932f9f8-812d-4ffc-b7c3-d036259636c7",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "Our aim is to identify model parameters $\\theta$ that *maximize* the likelihood. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "000eca63-a880-441d-a8e8-0edc3c2be96e",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "**Probelm 1a**\n",
    "\n",
    "Write a function `model` with two input parameters, `theta` and `x`, where `theta` is a tuple with values $\\theta_0$ and $\\theta_1$, and the function returns $\\theta_0 + \\theta_1 x$.\n",
    "\n",
    "*Hint* – this is far more formal than necessary, but it will simplify other problems later in the notebook. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "39939d3b-b0af-4eeb-8bd9-51776a5b9282",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def model(theta, x):\n",
    "    '''\n",
    "    Return dependent variable values for f(x) = theta_0 + theta_1 x\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    theta : tuple (size=2)\n",
    "        theta[0] is the intercept and theta[1] is the slope of the line\n",
    "    \n",
    "    x : array-like\n",
    "        values of the independent variable where f(x) should be evaluated\n",
    "    '''\n",
    "    \n",
    "    # complete\n",
    "    return # complete"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4bfd430b-3dec-452e-ad22-2e198f46dc20",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "**Problem 1b**\n",
    "\n",
    "Write a function `prob` that calculates and returns the probability of observations `x` assuming that $p(x)$ follows a normal distribution with mean `mu` and standard deviation `sigma`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f38f5f71-fa74-47cf-9a68-56e0fa9ce29a",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def prob(x, mu, sigma):\n",
    "    p = # complete\n",
    "    return p"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ffbfb606-b643-47f3-98e7-3211e51213b4",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "**Problem 1c**\n",
    "\n",
    "Calculate the likelihood for the observations that were simulated in **0b**. Use the model parameters that generated the data. \n",
    "\n",
    "*Hint* – think carefully about your variable names. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5dbc607d-da90-4175-99cf-0b3bf910ac70",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "lkhd = # complete\n",
    "\n",
    "print(f'The likelihood for this data set is {lkhd}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd89ce67-3ee0-453b-b89f-2b57474d4bac",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "That is a very small number! "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "974cba22-0c2a-49a0-a221-a23a332a565d",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "**Problem 1d**\n",
    "\n",
    "Using the same model parameters as before, generate data sets of 100, 200, and 300 observations, and calculate the likelihood for each.\n",
    "\n",
    "*Hint* – do not use the same variable names (e.g., `x_obs`, `y_obs`, etc) that we used before; we do not want to over-write those variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8ce64efd-9fcd-480f-81fd-fc290183e5cd",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "for n_obs in # complete\n",
    "\n",
    "    rng = np.random.default_rng(seed=2009)\n",
    "\n",
    "    x_sim = # complete\n",
    "    y_sim = # complete\n",
    "    y_sim_unc = # complete\n",
    "\n",
    "    lkhd = # complete\n",
    "\n",
    "    print(f'The likelihood for {n_obs} observations is {lkhd}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "88040a6d-0662-4166-bd7b-73ea105c50cb",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "300 observations is a small number! Gaia has observed more than 1 billion stars, and yet, with only 300 observations in an outrageously simple dataset, the likelihood is so small it is equivalent to 0 at machine precision. \n",
    "\n",
    "\n",
    "This is why it is always a good idea to work with the log of the likelihood, not only does this turn a product into a sum, more importantly, the calculations become far more stable on your machine. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c16d8209-92d5-445c-a307-3fdd5a73fa0f",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "**Problem 1e**\n",
    "\n",
    "Write a function `lnl` that calculates the log likelihood for some observations, their uncertainties, and the model to which the observations are being compared. Assume that the likelihood is Gaussian.\n",
    "\n",
    "*Hint* – the `model` function that was created earlier should be inside the `lnl` function, this means that `theta` should be the first arguement for the `lnl` function.\n",
    "\n",
    "*Note* – likelihoods are calculated for comparison purposes, their absolute value does not have much meaning, so you can ignore constant terms for this function. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0b0ef526-b8b8-41c0-a3d3-dc8510fa403d",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def lnl(theta, y_obs, y_obs_unc, x_obs):\n",
    "    # complete\n",
    "    return # complete"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a22ad979-586e-47b8-abba-597fbeece90a",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "**Problem 1f**\n",
    "\n",
    "Using the same model parameters as before, generate data sets of 100, 200, and 300 observations, and calculate the *log likelihood* for each."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e9c395f5-c4da-4da7-bab8-80e9323e57ab",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "for n_obs in [100, 200, 300]:\n",
    "\n",
    "    rng = np.random.default_rng(seed=2009)\n",
    "\n",
    "    x_sim = rng.uniform(0, 10, n_obs)\n",
    "    y_sim = rng.normal(3.14*x_sim + 6.626, 3)\n",
    "    y_sim_unc = np.ones_like(x_sim)*3\n",
    "\n",
    "    lnlike = # complete\n",
    "\n",
    "    print(f'The ln likelihood for {n_obs} observations is {lnlike}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "99a963cb-16df-4c4e-9d30-79a58305aac6",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## Problem 2) Maximizing the Likelihood\n",
    "\n",
    "It is all well and good to calculate the likelihood, but what we truly want is to maximize the likelihood. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d93fe98b-2e0e-41d9-bb64-e9df96288c6d",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "Most algorithms are designed to minimize, rather than optimize, a function. Fortunately, minimizing the negative log likelihood is the exact same as maximizing the log likelihood.\n",
    "\n",
    "**Problem 2a**\n",
    "\n",
    "Write a function `nll` to calculate the negative log likelihood. \n",
    "\n",
    "*Hint* – this is really simple, don't overthink it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "085f8a24-aa3d-45b5-9933-fb58ed71fb42",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def nll(theta, y_obs, y_obs_unc, x_obs):\n",
    "    return # complete"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "72a4e798-c538-4c25-94c2-7e0bde2d0986",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "**Problem 2b**\n",
    "\n",
    "Using `minimize` from `scipy.optimize` determine the maximum likelihood estimation for the intercept and slope of the line that was used to generate the synthetic observations. \n",
    "\n",
    "*Hint* – for arguments `minimize` needs (1) a function, (2) an initial guess for the model parameters, which is why we've been using a tuple `theta`, and (3) a tuple containing the remaining arguments for the function to be minimized (i.e., the data/observations)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a515e1a-9233-4715-ab24-8d5a0970ab55",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "res = # complete\n",
    "\n",
    "print(f'b = {res.x[0]:.4f} and m = {res.x[1]:.4f} for the MLE')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6c55453e-1b27-4db0-aa9b-864cf445a5e3",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "**Problem 2c**\n",
    "\n",
    "Overplot the line determined by the MLE on top of the synthetic data. \n",
    "\n",
    "How does the line compare to the true line used to generate the data?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c92bfb98-cf0e-4aae-827b-b33c258084cc",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, ax = data_plot(x_obs, y_obs, y_obs_unc)\n",
    "ax.plot( # complete\n",
    "ax.plot( # complete\n",
    "ax.legend()\n",
    "fig.tight_layout()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}