2-Optional-Probability.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# 2 (Optional) Probability, Distributions, and Measurement"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook is optional, but it's recommended if you want learn more about Python and/or probability theory."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## More on Python\n",
    "\n",
    "If you're new to programming in general, you should be working your way through all the chapters of the [Python for Data Science](https://www.datacamp.com/courses/intro-to-python-for-data-science) DataCamp tutorial.\n",
    "\n",
    "If you're a more accomplished programmer but new(ish) to Python, you can get a more detailed primer on the same material by going through chapters 1-4 of the [Python Data Science Handbook](https://jakevdp.github.io/PythonDataScienceHandbook/).\n",
    "\n",
    "Assuming you've completed these activities, the sections below will be a review of some key concepts."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Objects and packages\n",
    "\n",
    "Python is firmly grounded in the tradition of *object-oriented programming* (OOP). In this approach, the logic for dealing with different kinds of data is *encapsulated* or *attached* to the data. A data object can\n",
    "\n",
    "- have *attributes* that tell us about that object instance\n",
    "- have *methods* that are functions an object can perform\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### What is an instance?\n",
    "\n",
    "Every value in Python is an *instance* of some type. Each time you create a new value, the interpreter is *instantiating* an object of a defined type, allocating memory and setting initial values.\n",
    "\n",
    "Thus, any variable we define is simply pointing the variable's name (in the current namespace) to a chunk of memory containing the instance.\n",
    "\n",
    "Python keeps track of all object instances and removes them from memory when they are no longer needed.\n",
    "\n",
    "You can determine the type of a value or variable using the `type` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "x = 42\n",
    "type(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The attributes and methods of an object are accessed by using the `.` operator. You can *inspect* the attributes and methods of any object in the notebook by typing `tab` after entering the variable name and `.`, as below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# let's explore x:\n",
    "# (place the cursor after x. and press the tab key)\n",
    "x."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Let's try with a different type:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "y = 'The answer to the ultimate question.'\n",
    "type(y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# let's explore the string type:\n",
    "# (press tab after y.)\n",
    "y."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note: you can access documentation for functions and methods by using `Shift-Tab`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# place the cursor between the parentheses and type `Shift-Tab`\n",
    "y.lower()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question:** In your own words, what does the `lower` method do?\n",
    "\n",
    "(double-click the next cell and enter your answer)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### What is a package?\n",
    "\n",
    "Programming is all about abstraction, and code is meant to be reused. Python has tons of useful packages that extend its basic functionality. This means you don't have to (and shouldn't) reinvent basic algorithms. \n",
    "\n",
    "We'll talk more about how to find and evaluate packages later.\n",
    "\n",
    "You can access external packages (or modules) by **importing** them into your workspace. When you import a package, you should give it a short name that you'll use to access the functions, etc, in the package.\n",
    "\n",
    "You'll see (and start to use) lots of statements like the following:\n",
    "\n",
    "```python\n",
    "import numpy as np\n",
    "```\n",
    "\n",
    "After executing this statement, you can access the contents of the module as if they were attributes and methods of an object called `np`, like so: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "np.arange(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Arrays\n",
    "\n",
    "The `numpy` package provides Python with an important data structure, the **array**. Arrays are so important to scientific programming that you will almost inevitably import `numpy` at the beginning of every script and notebook you use.\n",
    "\n",
    "As we discussed previously, **scalar** data types like `int` and `float` can only represent a single number. But time series and point processes are inherently ordered collections of data, and arrays give us the ability to store and manipulate these large aggregates. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Array terminology\n",
    "\n",
    "- An **element** of an array is one item. It occupies a specific \"slot\" in the sequence.\n",
    "- The **length** or **size** of an array is the number of elements.\n",
    "- The **index** of an element is the numerical position of the element in the array. Python uses *zero-based indexing*, which means that the first element has the index `0`.\n",
    "- The **data type** of the array is the type of the elements in the array. In arrays, all the elements have the same type. In `numpy`, the data type is stored in the `dtype` attribute."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Array operations\n",
    "\n",
    "- **indexing** allows you to access specific elements of an array\n",
    "- **slicing** allows you to access specific subsets of an array\n",
    "- **iteration** allows you to process each element of the array in sequence\n",
    "\n",
    "In numpy, both indexing and slicing use square brackets (`[` and `]`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# first, let's initialize the random seed so that we all get the same answers\n",
    "np.random.seed(1)\n",
    "# create an array with 100 random numbers\n",
    "my_data = np.random.randn(100)\n",
    "# use indexing to get the first element:\n",
    "my_data[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Slicing also uses square brackets, but instead of a single index, you use two indices separated by `:` to specify a range:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this will retrieve the first 10 elements of my_array\n",
    "my_data[0:10]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice how indexing returns a scalar, but slicing returns a new array (technically, a view of a part of the original array). Also notice that the last index is *exclusive*; that is, element 10 is NOT returned."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Some simple problems\n",
    "\n",
    "Write an expression for the 88th value of `my_data` in the cell below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Write an expression for the last 10 elements of `my_data`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Basics of Probability\n",
    "\n",
    "Let's consider two random variables. Following standard notation, we'll use capital letters $A$ and $B$ to denote the random variable and lowercase letters $a$ and $b$ to denote specific outcomes of those random variables.\n",
    "\n",
    "**The probability of $a$:**\n",
    "\n",
    "$$Pr(A = a) \\in [0, 1]$$\n",
    "\n",
    "**The probability of not $a$:**\n",
    "\n",
    "$$Pr(A \\neq a) = 1 - Pr(A = a)$$\n",
    "\n",
    "**The probability of $a$ and $b$ (the joint probability):**\n",
    "\n",
    "$$p(a,b) = p(a|b) p(b) = p(b|a) p(a)$$\n",
    "\n",
    "Note that we're using $p(a)$ as a shorthand for $Pr(A = a)$. If $A$ and $B$ are independent,\n",
    "\n",
    "$$p(a, b) = p(a) p(b)$$\n",
    "\n",
    "**The probability of $a$ given $b$ (the conditional probability):**\n",
    "\n",
    "$$p(a \\mid b) = \\frac{p(a, b)}{p(b)} = \\frac{p(b \\mid a) p(a)}{p(b)}$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Simple exercises\n",
    "\n",
    "Using the equations above, let's try and figure out the probabilities of the following:\n",
    "\n",
    "1. Rolling a 5 on a 6-sided die.\n",
    "2. Not rolling a 3 on a 6-sided die.\n",
    "3. Rolling a 4 or a 5 on a 6-sided die.\n",
    "4. Rolling less than 4 or an even number on a 6-sided die.\n",
    "5. Rolling two 3's in a row on a 6-sided die.\n",
    "\n",
    "Enter your answers as text in the cell below:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Harder question\n",
    "\n",
    "#### The case of Tim Tebow\n",
    "\n",
    "What is the probability of becoming a Major League Baseball (MLB) player if you hit a home run (HR) in your first at-bat in the minor leagues? We have the following important information:\n",
    "\n",
    "a) 5% of future MLB players hit a home run in their first minor league at-bat.\n",
    "\n",
    "b) 1% of minor league players make it to MLB.\n",
    "\n",
    "c) Only 0.1% of players hit a homerun in their first minor league at-bat.\n",
    "\n",
    "Hint: You are trying to solve $Pr(MLB = 1 \\mid HR = 1)$.\n",
    "\n",
    "**Bonus**: What percentage of players who don't make it to MLB hit a home run in the first at bat?\n",
    "\n",
    "Write code that *prints the answer* to the question(s):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## Probability Distributions\n",
    "\n",
    "Some key terms:\n",
    "\n",
    "- A **random variable** can take on different values depending on the outcome of a random phenomenon.\n",
    "  - Random variables are usually denoted by capital letters (e.g. $X$)\n",
    "  - A concrete value of a random variable is denoted by a lowercase letter (e.g. $x$)\n",
    "- The **support** of a random variable is set of of values it can have\n",
    "  - A **discrete** support includes only a finite, or countably infinite, set of values. \n",
    "  - A **continuous** support includes all the real numbers within some range. The number of possible values is *uncountably* infinite.\n",
    "- A **distribution** describes the probability of a random variable having some value.\n",
    "  - Distributions can be mathematical functions of some number of **parameters**.\n",
    "      - Parameters are NOT observable\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Discrete vs continuous distributions\n",
    "\n",
    "A **probability mass function** (*pmf*) describes the distribution of a *discrete* variable. For example, a fair N-sided die has a probability of coming up as 1 equal to $1/N$:\n",
    "\n",
    "$$Pr(X = 1) = 1/N$$\n",
    "\n",
    "A **probability density function** (*pdf*) describes the distribution of a *continuous* variable. Density functions can be used to calculate the probability of obtaining a value in a range, using integration:\n",
    "\n",
    "$$Pr(a \\leq X \\leq b) = \\int_a^b p(x) dx$$ \n",
    "\n",
    "We will $p$ to denote probability density and mass functions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Examples\n",
    "\n",
    "Let's look at some common discrete and continuous probability distributions. Run the following two cells to import some functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# load matplotlib inline mode\n",
    "%matplotlib inline\n",
    "\n",
    "# import some useful libraries\n",
    "import numpy as np                # numerical analysis linear algebra\n",
    "import matplotlib.pyplot as plt   # plotting\n",
    "\n",
    "# import some distributions\n",
    "import sys\n",
    "sys.path.insert(0,\"/project/psyc5270-cdm8j/comp-neurosci\")\n",
    "from comp_neurosci_uva.dists import uniform, normal, beta, gamma, invgamma, exp, poisson, laplace, students_t, noncentral_t, halfcauchy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this cell defines a function that will let us plot probability distributions easily\n",
    "\n",
    "def plot_prob(dist, support=[-5, 5], npoints=100, **kwargs):\n",
    "    \"\"\"Plot a probability distribution over a support interval.\n",
    "    \n",
    "    dist - a scipy distribution function, discrete or continuous\n",
    "    support - the range of values over which to evaluate the distribution\n",
    "    npoints - the number of points within the support to evaluate (for density functions only)\n",
    "    \"\"\"\n",
    "    try:\n",
    "        x = np.linspace(support[0], support[1], npoints)        \n",
    "        prob = dist.pdf(x)\n",
    "    except AttributeError:\n",
    "        x = np.arange(support[0], support[1])\n",
    "        prob = dist.pmf(x)\n",
    "    plt.plot(x, prob, lw=3, **kwargs)\n",
    "    plt.xlabel('Value')\n",
    "    plt.ylabel('Probability')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### [Uniform](https://en.wikipedia.org/wiki/Uniform_distribution_(continuous))\n",
    "\n",
    "A continuous probability distribution assigning equal probability over a (continuous) range. The parameters of the uniform distribution are the upper and lower bounds. The density of a uniform distribution is proportional to the difference between the bounds\n",
    "\n",
    "$$U(l,u) = \\frac{1}{u - l}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the PDF\n",
    "plot_prob(uniform(lower=-2, upper=2))\n",
    "plot_prob(uniform(lower=-1, upper=1))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question**: Is this a density or a mass function? Why is the height of the distribution lower when the range is larger?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### [Normal](https://en.wikipedia.org/wiki/Normal_distribution)\n",
    "\n",
    "Normal distributions are found throughout math and nature, due to the Central Limit Theorem. Also known as Gaussian distributions, after its discoverer, Carl Friedrich Gauss. The normal distribution also has two parameters, the mean ($\\mu$) and standard deviation ($\\sigma)$. If a random variable $X$ is drawn from a normal distribution, it can be denoted as: \n",
    "\n",
    "$$X \\sim N(\\mu,\\sigma^2)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_prob(normal(mean=0, std=1))\n",
    "plot_prob(normal(mean=1, std=2), support=[-5, 8])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### [Beta](https://en.wikipedia.org/wiki/Beta_distribution)\n",
    "\n",
    "Only has support between 0 and 1. Useful to help determine the probability of a probability. Also has two parameters, $\\beta$ and $\\alpha$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_prob(beta(alpha=0.5, beta=0.5), support=[0,1])\n",
    "plot_prob(beta(alpha=2, beta=5), support=[0,1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### [Poisson](https://en.wikipedia.org/wiki/Poisson_distribution)\n",
    "\n",
    "The Poisson distribution reflects the probability of observing some number of events within a window of time (for example, the number of buses that go by some point on the road). The Poisson distribution has only a single parameter, $\\lambda$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_prob(poisson(1), support=[0,10], marker='o')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question**: Is the poisson distribution a density or a mass function? Why?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### And more...\n",
    "\n",
    "Take some time now to explore:\n",
    "\n",
    "- Gamma\n",
    "- Inverse Gamma\n",
    "- Exponential\n",
    "- Student's t\n",
    "- Half Cauchy\n",
    "- Poisson\n",
    "\n",
    "**Assignment**: choose two distributions. Try to replicate the illustrative plots on their respective pages on Wikipedia. Insert code cells below to generate the plots."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Observation and Inference\n",
    "\n",
    "Why is probability theory important to computational neuroscience? \n",
    "\n",
    "Fundamentally, the problem we face is that we can't directly determine any physical quantity. We can only take a measurement of it, and that measurement will have errors. Thus, each time we make a measurement, we are going to get a value that comes from a **distribution**.\n",
    "\n",
    "Even more troubling, we rarely are able to directly measure the actual quantities we care about. Think back to the last exercise (or run it, if you don't remember!). It was clear that a limited set of sounds were activating the neuron (i.e., generating synaptic excitation), but we didn't have a direct measurement of how strong that excitation was. We could only observe that the neuron produced action potentials at higher rates during certain intervals.\n",
    "\n",
    "The process of using observations to gain information about unobservable quantities is called **inference** or **estimation**, and probability theory gives us the tools we need to make this connection."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Sampling and Measurement\n",
    "\n",
    "Let's come up with a **observational model** that formalizes what we think is going on when we make a measurement.\n",
    "\n",
    "Observational models are fundamental to all statistical analyses and will be a part of many of the more complex models we develop.\n",
    "\n",
    "Let's say that we have a bar of length $\\mu$. We can measure the bar's length as many times as we like, and each time we do so, we'll get a value $y_i$:\n",
    "\n",
    "$$y_i = \\mu + \\varepsilon_i$$\n",
    "\n",
    "Notice the subscript $i$. This is a numerical **index** for the measurement. Let's stick with the Python convention and have $y_0$ indicate the value of the first measurement.\n",
    "\n",
    "This model simply says that any given measurement $y_i$ will be the true length of the bar plus some (hopefully small) error $\\varepsilon_i$. That is, we're assuming **additive error**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### The normal error model\n",
    "\n",
    "There are many potential sources of error in our measurement.\n",
    "\n",
    "The bar might be fluctuating slightly in length due to changes in temperature. This is an example of an **intrinsic error**.\n",
    "\n",
    "The instrument making the measurement might have limited precision due to how it's constructed or because of electrical interference. This is an example of **extrinsic error**.\n",
    "\n",
    "Thanks to the [Central Limit Theorem](https://en.wikipedia.org/wiki/Central_limit_theorem), a reasonable assumption is that the sum of all these sources of error will have a normal (or Gaussian) distribution.\n",
    "\n",
    "We can also safely assume that each measurement of the bar doesn't depend on previous measurements of the bar (this assumption would not hold if the measurement actually affected the bar's length), but only on the true length on the bar and the random error.\n",
    "\n",
    "We can formalize these assumptions by saying that $\\varepsilon_i$ is **drawn** or **sampled** from a normal distribution. This relationship is often signified with $\\sim$.\n",
    "\n",
    "$$\\varepsilon_i \\sim N(0, \\sigma^2)$$\n",
    "\n",
    "$N$ represents the normal distribution, and as we saw above, this distribution has two parameters: mean and variance.\n",
    "\n",
    "**Question**: Why is the mean for the error distribution zero? What does the variance correspond to?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Simulating measurements\n",
    "\n",
    "Python can generate numbers from a distribution. The numbers are usually not truly random, because computers are (generally) deterministic, but the values will occur with the probability specified by the underlying PDF. Let's generate 10 measurements of our bar:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from comp_neurosci_uva import data\n",
    "# true length of bar:\n",
    "mu = 12.1\n",
    "# number of measurements\n",
    "N = 10\n",
    "# errors\n",
    "epsilon = normal(mean=0.0, std=data.e2_std).rvs(N)\n",
    "# measurements\n",
    "y = mu + epsilon\n",
    "\n",
    "# inspect the variables\n",
    "print(\"errors:\", epsilon)\n",
    "print(\"measurements:\", y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An aside: numpy arrays support **broadcasting**, which is what allows us to simply add a scalar `mu` to an array `epsilon` and get a new array in which `mu` has been added to every element of `epsilon`. If you're used to lower-level langauges like C or Java, make sure you take advantage of this feature, as it's MUCH faster than iterating through the array."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Histograms\n",
    "\n",
    "It's useful to look at the raw data, but distributions are usually better visualized as **histograms**. A histogram is a plot that divides a range into a set of intervals or **bins**, and then counts the number of values in each bin.\n",
    "\n",
    "Matplotlib has a histogram function, so no need to reinvent the wheel here:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "freq, bins, _ = plt.hist(y, range=(10, 14), bins=10, density=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Oops! All the observations fall in two bins. **Assignment**: Copy the code cell above, paste it below this cell, and edit to change the `range` and/or `bins` parameters so that the plot gives you more useful information."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Fitting the observational model\n",
    "\n",
    "Once you've got a nice histogram, try fitting it to a normal PDF. Copy the plot statement from the code cell above into the cell below, then adjust the mean and std parameters in the `plot_prob` function until you get what looks like a good fit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### COPY histogram command in the line below\n",
    "\n",
    "### EDIT this line to adjust the PDF\n",
    "plot_prob(normal(mean=???, std=???), support=(11, 13))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `mean` and `std` parameters you settle on are called your `estimate`.\n",
    "\n",
    "It may not be possible to achieve a good fit with only two observations. What effect does this have on your ability to come up with good parameter estimates?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Explore whether increasing the sample size helps. Copy and paste ONLY the relevant lines of code from above and paste them into the cell below. (Normally we avoid copying and pasting like the plague, but this will help you to think about what each statement is doing so that you can choose ONLY the relevant ones)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## paste your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Summary statistics\n",
    "\n",
    "Obviously, manually fitting a PDF to the histogram is both tedious and error-prone.\n",
    "\n",
    "Because we're using a normal error model, we can estimate the parameters directly by computing the *sample mean* and *sample standard deviation*.\n",
    "\n",
    "Consult the numpy documentation (see link under the `Help` menu) and find the array methods that will compute these summary statistics, then edit the code cell below so that it prints out the values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"The mean is\", y.<insert-method-call-here>)\n",
    "print(\"The standard deviation is\", y.<insert-method-call-here>)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question**: How close were your estimates to the summary statistics? How close are the summary statistics to the values of $\\mu$ and $\\sigma$ that were used to generate the data?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Constraining models with data\n",
    "\n",
    "### Readings and other resources\n",
    "\n",
    "- If you'd like to explore the concept of maximum likelihood more, check out this [interactive visualization](https://rpsychologist.com/d3/likelihood/) from [@krstoffr](https://twitter.com/krstoffr)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Data as a random variables\n",
    "\n",
    "Most statistics operates under the assumption that any observed data are actual samples drawn from some (to be learned) model of the world.\n",
    "\n",
    "In the standard approach, these models are simply probability distributions, with parameters that govern the behavior of the model (i.e., what we can expect they will produce)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Uncertainty\n",
    "\n",
    "If the observed data samples are random draws from a probability distribution with unknown parameters, then the level of uncertainty about those parameters will decrease as we gain more data samples.\n",
    "\n",
    "Statistical inference involves figuring out what model (e.g., probability distribution, but we will be building more complicated models later in the course) and parameters generated the data. \n",
    "\n",
    "Let's spend some time trying to perform this task by hand..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# load matplotlib inline mode\n",
    "%matplotlib inline\n",
    "\n",
    "# import some useful libraries\n",
    "import numpy as np                # numerical analysis linear algebra\n",
    "import pandas as pd               # efficient tables\n",
    "import matplotlib.pyplot as plt   # plotting\n",
    "import ipywidgets as widgets      # interactive widgets\n",
    "from ipywidgets.widgets.interaction import show_inline_matplotlib_plots\n",
    "from IPython.display import display, clear_output\n",
    "import pickle\n",
    "\n",
    "# import the distributions wrapped from scipy\n",
    "from comp_neurosci_uva import dists"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# load in the data\n",
    "with open('data/random_data.pickle', 'rb') as f:\n",
    "    data = pickle.load(f)\n",
    "\n",
    "# tell us something about the data\n",
    "for i,d in enumerate(data):\n",
    "    print('Dataset %d has %d samples' % (i, len(d)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# mapper between dist name and object\n",
    "dist_dict = {'Beta': dists.beta,\n",
    "             'Exponential': dists.exp,\n",
    "             'Gamma': dists.gamma,\n",
    "             'Normal': dists.normal,\n",
    "             'Uniform': dists.uniform}\n",
    "\n",
    "# Dropdown widget to pick datasets\n",
    "ds_ind = widgets.Dropdown(options=range(len(data)),\n",
    "                          description=\"Dataset\")\n",
    "\n",
    "# Checkbox for whether to show log likelihood\n",
    "like_check = widgets.Checkbox(description='Show Log Likelihood', \n",
    "                              value=False)\n",
    "\n",
    "# Checkbox for whether to show stem plot\n",
    "stem_check = widgets.Checkbox(description='Show Stem Plot', \n",
    "                              value=False)\n",
    "\n",
    "# set up the distributions tab\n",
    "dist_tab = widgets.Tab(allow_none=False)\n",
    "children = []\n",
    "titles = []\n",
    "\n",
    "# beta params\n",
    "beta_alpha = widgets.FloatText(value=.5,\n",
    "                               description='alpha')\n",
    "beta_beta = widgets.FloatText(value=.5,\n",
    "                              description='beta')\n",
    "children.append(widgets.HBox([beta_alpha, beta_beta]))\n",
    "titles.append('Beta')\n",
    "\n",
    "# exp params\n",
    "exp_lam = widgets.FloatText(value=5.0,\n",
    "                            description='lam')\n",
    "children.append(widgets.HBox([exp_lam]))\n",
    "titles.append('Exponential')\n",
    "\n",
    "# gamma params\n",
    "gamma_alpha = widgets.FloatText(value=.5,\n",
    "                                description='alpha')\n",
    "gamma_beta = widgets.FloatText(value=.5,\n",
    "                              description='beta')\n",
    "children.append(widgets.HBox([gamma_alpha, gamma_beta]))\n",
    "titles.append('Gamma')\n",
    "\n",
    "# normal params\n",
    "normal_mean = widgets.FloatText(value=0.0,\n",
    "                                description='mean')\n",
    "normal_std = widgets.FloatText(value=1.0,\n",
    "                               description='std')\n",
    "children.append(widgets.HBox([normal_mean, normal_std]))\n",
    "titles.append('Normal')\n",
    "\n",
    "# uniform params\n",
    "uniform_lower = widgets.FloatText(value=0.0,\n",
    "                                  description='lower')\n",
    "uniform_upper = widgets.FloatText(value=1.0,\n",
    "                                  description='upper')\n",
    "children.append(widgets.HBox([uniform_lower, uniform_upper]))\n",
    "titles.append('Uniform')\n",
    "\n",
    "# add all the children and set the tab titles\n",
    "dist_tab.children = children\n",
    "for i in range(len(titles)):\n",
    "    dist_tab.set_title(i, titles[i])\n",
    "\n",
    "# set the full user interface\n",
    "ui = widgets.VBox([widgets.HBox([ds_ind, like_check, stem_check]), dist_tab])\n",
    "\n",
    "# define plotting function\n",
    "def plot_data_and_dist(*vals, **kwargs):\n",
    "    # first plot the data\n",
    "    dat = data[ds_ind.value]\n",
    "    plt.hist(dat, bins='auto', density=True, alpha=.5);\n",
    "    \n",
    "    # now plot the pdf of the dist\n",
    "    npoints = 100\n",
    "    \n",
    "    # add support for 10% of the data range on either side\n",
    "    support = (dat.min() - np.ptp(dat)*.1,\n",
    "               dat.max() + np.ptp(dat)*.1)\n",
    "    x = np.linspace(support[0], support[1], npoints)\n",
    "    \n",
    "    # get the selected dist and params\n",
    "    params = {c.description: c.value for c in \n",
    "              dist_tab.children[dist_tab.selected_index].children}\n",
    "    dist = dist_dict[dist_tab.get_title(dist_tab.selected_index)](**params)\n",
    "    \n",
    "    # calculate the pdf\n",
    "    pdf = dist.pdf(x)\n",
    "            \n",
    "    # plot the pdf and add labels\n",
    "    plt.plot(x, pdf, lw=3)\n",
    "    plt.xlabel('Value')\n",
    "    plt.ylabel('Probability')\n",
    "    if like_check.value:\n",
    "        # calculate the log like\n",
    "        log_like = np.log(dist.pdf(dat)).sum()\n",
    "\n",
    "        # add it to the plot with some formatting\n",
    "        plt.title('Log Like: {:3.4f}'.format(log_like))\n",
    "    \n",
    "    if stem_check.value:\n",
    "        # include the stem plot\n",
    "        plt.stem(dat, dist.pdf(dat), 'g')\n",
    "\n",
    "        \n",
    "# set up triggers for updating the plot\n",
    "out = widgets.interactive_output(plot_data_and_dist, \n",
    "                                 {'ds_ind': ds_ind,\n",
    "                                  'like_check': like_check,\n",
    "                                  'stem_check': stem_check,\n",
    "                                  'beta_alpha': beta_alpha,\n",
    "                                  'beta_beta': beta_beta,\n",
    "                                  'exp_lam': exp_lam,\n",
    "                                  'gamma_alpha': gamma_alpha,\n",
    "                                  'gamma_beta': gamma_beta,\n",
    "                                  'normal_mean': normal_mean,\n",
    "                                  'normal_std': normal_std,\n",
    "                                  'uniform_lower': uniform_lower,\n",
    "                                  'uniform_upper': uniform_upper,\n",
    "                                 })\n",
    "\n",
    "# wrapper for tab change\n",
    "# required b/c selecting tabs can't trigger plots like other widgets\n",
    "def tab_change(*args, **kwargs):\n",
    "    with out:\n",
    "        clear_output(wait=True)\n",
    "        plot_data_and_dist()\n",
    "        show_inline_matplotlib_plots()\n",
    "dist_tab.observe(tab_change, 'selected_index')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interactive data fitting (part 2)\n",
    "\n",
    "The widget below (h/t Per Sederberg) shows the histogram of a randomly generated distribution in blue, overlaid with the PDF from one of five different distribution types in orange. Select one of the datasets and try to figure out what distribution it was generated from. You'll likely need to adjust the parameters of the distribution.\n",
    "\n",
    "Refer back to Notebook 2 and/or Wikipedia if you want to remind yourself of what the different distributions look like and how they change shape with their parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# show everything\n",
    "display(ui, out)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q**: What dataset did you pick? What distribution and parameters best fit the observations by eye?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Likelihood calculation\n",
    "\n",
    "We don't have to do this by eye. Because our model has a probability density function (PDF), we can calculate the likelihood of observing the data given the model and parameters.\n",
    "\n",
    "For any given model and parameters, you can determine the probability of having observed any individual data point by evaluating the PDF at the value of that data point (*Turn on the Stem Plot.*)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Maximum Likelihood Estimation\n",
    "\n",
    "The goal then becomes to maximize the likelihood of observing the data $y = {y_0, y_1, \\ldots y_{n -1}$ given a model with a one or more parameters $\\theta$:\n",
    "\n",
    "$$p(y \\mid \\theta)$$\n",
    "\n",
    "As long as the data points are all independent, the likelihood of observing all of them is the product of all the probabilities.\n",
    "\n",
    "$$p(y \\mid \\theta) = \\prod_i p(y_i \\mid \\theta)$$\n",
    "\n",
    "It is more efficient and computationally tractable to maximize the log of $p(y|\\theta)$, so we typically convert the likelihood into a sum of log likelihoods.\n",
    "\n",
    "$$ \\log p(y \\mid \\theta) = \\sum_i \\log p(y_i \\mid \\theta)$$\n",
    "\n",
    "**Q** Turn on the Log Likelihood checkbox above and see if you can do better. Note that you'll get a numerical error if any of the data lie outside the support of the distribution, but this will go away once you fix that issue. In the cell below, enter your best estimates of the parameters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Automated optimization\n",
    "\n",
    "Many approaches have been developed for searching parameter spaces to find the parameters that generate the maximum or minimum value of a function.\n",
    "\n",
    "### Nelder--Mead Simplex\n",
    "\n",
    "One very popular algorithm is the Nelder--Mead simplex.\n",
    "\n",
    "It involves growing and shrinking a simplex (a generalization of a triangle to multiple dimensions) to search the parameter space efficiently to minimize a function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "![](https://upload.wikimedia.org/wikipedia/commons/thumb/d/de/Nelder-Mead_Himmelblau.gif/640px-Nelder-Mead_Himmelblau.gif)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# scipy includes lots of optimization methods\n",
    "import scipy.optimize as opt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define a likelihood function\n",
    "def like_fun(params, *args):\n",
    "    # pull the model and dat out of the args\n",
    "    model = args[0]\n",
    "    dat = args[1]\n",
    "    \n",
    "    # instantiate the model with the params\n",
    "    dist = model(*params)\n",
    "    \n",
    "    # calc the log like\n",
    "    log_like = np.log(dist.pdf(dat)).sum()\n",
    "    if np.isnan(log_like):\n",
    "        log_like = -np.inf\n",
    "    \n",
    "    # return the negative of it to minimize\n",
    "    return -log_like"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate the best-fitting params\n",
    "\n",
    "Execute the following code block to run Nelder-Mead on your dataset, using your chosen model from the widget above. Hopefully, in some small amount of time, the algorithm will converge on a solution, and you'll get a printout. See if you can figure out what the best parameter estimates are from the output.\n",
    "\n",
    "Notice that each time you run the code below, a new starting value for the optimization is generated. Try running the optimization a few times and see if you get similar answers. Why is it important to do this?\n",
    "\n",
    "If you're not sure if you've chosen the right PDF, go back to the widget, pick a different one, and rerun the cell below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# grab the distribution and data from the UI above\n",
    "dist_name = dist_tab.get_title(dist_tab.selected_index)\n",
    "model = dist_dict[dist_name]\n",
    "dat = data[ds_ind.value]\n",
    "\n",
    "# set the bounds for the distribution\n",
    "bound_dict = {'Beta': [(0, 10), (0, 10)],\n",
    "              'Exponential': [(0, 20)],\n",
    "              'Gamma': [(0, 10), (0, 10)],\n",
    "              'Normal': [(-10, 10), (0, 10)],\n",
    "              'Uniform': [(-10, 10), (-10, 10)]}\n",
    "bounds = bound_dict[dist_name]\n",
    "\n",
    "# generate a random starting point based on the bounds\n",
    "# NB: it's possible to generate invalid starting points\n",
    "x0 = [dists.uniform(*b).rvs() for b in bounds]\n",
    "\n",
    "# print some information about the distribution and starting values\n",
    "print('Dataset:', ds_ind.value)\n",
    "print('Distribution:', dist_name)\n",
    "print('Starting value:', x0)\n",
    "print()\n",
    "\n",
    "# run the optimizer\n",
    "# NOTE, not all methods make use of the bounds method\n",
    "res = opt.minimize(like_fun, x0, args=(model, dat), \n",
    "                   #bounds=bounds,\n",
    "                   #method='L-BFGS-B',\n",
    "                   #method='BFGS',\n",
    "                   method='Nelder-Mead',\n",
    "                   #method='TNC'\n",
    "                  )\n",
    "print(res)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Visualize the fit\n",
    "\n",
    "Let's visualize the likelihood as a 2D function of the parameters and see how our fit looks. If you're far away from where you think you should be, try picking some different values in the widget and/or re-running the optimization.\n",
    "\n",
    "Note that this will only work for PDFs that have two parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set the grid of points to evaluate\n",
    "x = np.linspace(bounds[0][0], bounds[0][1], 100)\n",
    "y = np.linspace(bounds[1][0], bounds[1][1], 100)\n",
    "xx, yy = np.meshgrid(x, y, sparse=True)\n",
    "\n",
    "# evaluate the likelihood\n",
    "z = np.log(model(xx, yy).pdf(dat[:, np.newaxis, np.newaxis])).sum(0)\n",
    "z = np.exp(z)\n",
    "\n",
    "# plot the contour and the best-fit value\n",
    "plt.contourf(x, y, z, 100)\n",
    "plt.plot(res.x[0], res.x[1], 'x', markersize=12, color='red')\n",
    "plt.colorbar()\n",
    "#plt.xlabel()\n",
    "#plt.ylabel()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q**: Is your estimate at the peak of the likelihood function? Are some parameters constrained more or less than others?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Assessing model fit\n",
    "\n",
    "An important question you need to address when using a statistical model is how well your model fits the data.\n",
    "\n",
    "Your fit can be bad for any number of reasons, but one rather common one is that your model is wrong.\n",
    "\n",
    "If your model is really wrong, then the parameter estimates are not worth much.\n",
    "\n",
    "However, because all models are wrong (see Box's Dictum), you can only compare models to see which are less wrong.\n",
    "\n",
    "How do we do this?\n",
    "\n",
    "We could simply compare the maximum likelihoods between the models, but that doesn't take into account the complexity of each model and amount of data. \n",
    "\n",
    "One, more principled, approach is Bayesian Information Criterion (BIC):\n",
    "\n",
    "$$BIC = \\text{ln}(n)k - 2\\text{ln}(\\hat{L}),$$\n",
    "\n",
    "where $n$ is the number of data points, $k$ is the number of parameters, and $\\hat{L}$ is the maximum likelihood value of the model $M$.\n",
    "\n",
    "***SMALLER BIC values are better!!!***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the model fit with BIC\n",
    "# recall that the optimization returns the negative log likelihood\n",
    "n = len(dat)\n",
    "k = len(res.x)\n",
    "L = -res.fun\n",
    "bic = np.log(n)*k - 2*(L)\n",
    "print('BIC:', bic)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model Comparison\n",
    "\n",
    "But how do we decide if one model is better than the other? We can compare BIC values between models, turning them into a Bayes Factor!\n",
    "\n",
    "$$BF = exp((BIC_0 - BIC_1)/2)$$\n",
    "\n",
    "This is interpreted with the help of the following guidelines:\n",
    "\n",
    "| Bayes Factor | Evidence |\n",
    "|--------------|----------|\n",
    "| 1--3         | Weak     |\n",
    "| 3--20        | Positive |\n",
    "| 20--150      | Strong   |\n",
    "| >150         | Very Strong | \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q** Copy your first BIC value into the cell below, assigning it to `bic_0`. Now try fitting your data to another  distribution functions using the widget above. Assign the second BIC value to `bic_1` below and the execute the cell to calculate the Bayes Factor. Then, in the text cell below, describe which two PDFs you considered, the maximum likelihood *parameter estimates*, and which model fit the best."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calc Bayes Factor\n",
    "# (enter numbers from the BIC assessments above)\n",
    "bic_0 = -33.93\n",
    "bic_1 = -31.07\n",
    "np.exp((bic_0 - bic_1)/2.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## t-test example\n",
    "\n",
    "We now have all the tools necessary to perform statistical inference (though we will improve on all these approaches in the coming weeks). You can:\n",
    "\n",
    "- Use optimization techniques to identify the parameters that give rise to the maximum likelihood of observing the data given the model\n",
    "- Assess model fit\n",
    "- Compare models to guide model selection\n",
    "\n",
    "Let's try a simple example of performing a t-test via model comparison approaches!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "You're hopefully quite familiar with how you can use a t-test to determine if a the mean of a population is significantly different from zero, or to determine if two populations have different means. But what's going on under the hood?\n",
    "\n",
    "Let's consider the one-sample t-test. You run a t-test and get a $p$ value, which tells you the probability that you'd get your observations if the data were from a normal distribution with mean $\\mu = 0$ (the null hypothesis).\n",
    "\n",
    "The model comparison standpoint is similar, but what you're testing is whether there's more evidence for a model where $\\mu \\neq 0$ than for one in which $\\mu = 0$. We can make this comparison using Bayes Factors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# generate some data that may or may not be significantly different from zero.\n",
    "# you can play around with the sample size, mean, and std deviation later\n",
    "A = dists.normal(.3, .5).rvs(10)\n",
    "\n",
    "# plot it\n",
    "plt.hist(A, bins='auto', density=True);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Standard t-test\n",
    "\n",
    "First we'll perform a standard one-sample t-test on our data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# do a one-sample t-test\n",
    "import scipy.stats as stats\n",
    "\n",
    "stats.ttest_1samp(A, 0.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Fit a Student's t model\n",
    "\n",
    "Next we fit the full model to the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define a likelihood function\n",
    "def students_like(params, *args):\n",
    "    # pull the model and dat out of the args\n",
    "    dat = args[0]\n",
    "    df = len(dat) - 1\n",
    "    \n",
    "    # instantiate the model with the params, \n",
    "    # the df is determined from the data\n",
    "    dist = dists.students_t(params[0], params[1], df)\n",
    "    \n",
    "    # calc the log like\n",
    "    log_like = np.log(dist.pdf(dat)).sum()\n",
    "    if np.isnan(log_like):\n",
    "        log_like = -np.inf\n",
    "    \n",
    "    # return the negative of it to minimize\n",
    "    return -log_like\n",
    "\n",
    "# pick an central starting point\n",
    "x0 = [0.0, 1.0]\n",
    "\n",
    "# run the optimization\n",
    "res = opt.minimize(students_like, x0, args=(A,), \n",
    "                   #bounds=bounds,\n",
    "                   #method='L-BFGS-B',\n",
    "                   #method='BFGS',\n",
    "                   method='Nelder-Mead',\n",
    "                   #method='TNC'\n",
    "                  )\n",
    "print(res)\n",
    "\n",
    "# calculate the BIC for this model and save it\n",
    "n = len(A)\n",
    "k = len(res.x)\n",
    "L = -res.fun\n",
    "bic_1 = np.log(n)*k - 2*(L)\n",
    "print('BIC:', bic_1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Fit a null hypothesis model\n",
    "\n",
    "Now we fit a model representing the null hypothesis that the mean of the data is actually 0.0. \n",
    "\n",
    "Note how we simply fix the mean of the Student's t distribution to zero, but still fit the standard deviation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define a likelihood function\n",
    "def students_null_like(params, *args):\n",
    "    # pull the model and dat out of the args\n",
    "    dat = args[0]\n",
    "    df = len(dat) - 1\n",
    "    \n",
    "    # instantiate the model with the params\n",
    "    # mean is fixed at zero and the df is determined from the data\n",
    "    dist = dists.students_t(0.0, params[0], df)\n",
    "    \n",
    "    # calc the log like\n",
    "    log_like = np.log(dist.pdf(dat)).sum()\n",
    "    if np.isnan(log_like):\n",
    "        log_like = -np.inf\n",
    "    \n",
    "    # return the negative of it to minimize\n",
    "    return -log_like\n",
    "\n",
    "# start at same point (though mean is fixed at zero)\n",
    "x0 = [1.0]\n",
    "\n",
    "# run the optimization\n",
    "res = opt.minimize(students_null_like, x0, args=(A,), \n",
    "                   #bounds=bounds,\n",
    "                   #method='L-BFGS-B',\n",
    "                   #method='BFGS',\n",
    "                   method='Nelder-Mead',\n",
    "                   #method='TNC'\n",
    "                  )\n",
    "print(res)\n",
    "\n",
    "# calculate and print the BIC\n",
    "n = len(A)\n",
    "k = len(res.x)\n",
    "L = -res.fun\n",
    "bic_0 = np.log(n)*k - 2*(L)\n",
    "print('BIC:', bic_0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Model Comparison\n",
    "\n",
    "Now that we have the BIC values for each model, we can use the Bayes Factor to determine whether the full model is preferred to the null model. \n",
    "\n",
    "We want a big number here. If it's less than 1.0 then there is no evidence that the alternative/full model should be preferred to the null model (i.e., the mean of the distribution is not different from 0.0)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate Bayes Factor\n",
    "bf = np.exp((bic_0 - bic_1)/2.)\n",
    "print('Bayes Factor:', bf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q**: In the cell that generates the observations, try three different values for the mean, standard deviation OR sample size. Then run the model-fitting code for each case. In the cell below, compare and contrast how changing these values affect the Bayes Factor and the p-value in the standard t-test."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  }
 ],
 "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.4"
  },
  "rise": {
   "scroll": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}