{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Monte Carlo simulation of a simple experiment using `scipy.stats` tools\n",
    "\n",
    "Question: I flip a fair coin 100 times.  What is the probability that I get 60 or more \"heads\"?\n",
    "\n",
    "One way to answer this question is to use the known properties of the binomial distribution (see the end of this notebook), but I want to get the answer by using tools from the `scipy.stats` submodule."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import stats\n",
    "\n",
    "import matplotlib as mpl       # As of July 2017 Bucknell computers use v. 2.x \n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib notebook\n",
    "# As of Aug. 2017 M.L. reverting to 1.x defaults.\n",
    "# In 2.x text.ustex requires dvipng, texlive-latex-extra, and texlive-fonts-recommended, \n",
    "# which don't seem to be universal\n",
    "# See https://stackoverflow.com/questions/38906356/error-running-matplotlib-in-latex-type1cm?\n",
    "mpl.style.use('classic')\n",
    "        \n",
    "# M.L. modifications of matplotlib defaults using syntax of v.2.0 \n",
    "# More info at http://matplotlib.org/2.0.0/users/deflt_style_changes.html\n",
    "# Changes can also be put in matplotlibrc file, or effected using mpl.rcParams[]\n",
    "plt.rc('figure', figsize = (6, 4.5))            # Reduces overall size of figures\n",
    "plt.rc('axes', labelsize=16, titlesize=14)\n",
    "plt.rc('figure', autolayout = True)             # Adjusts supblot parameters for new size"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Simulating spin flips\n",
    "\n",
    "I am going to generate 100 random 0's and 1's, and let 1 represent a \"heads\" and 0 represent a \"tails.\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1 1 1 1 0 1 0 0 0 1 1 0 0 0 0 1 0 0 1 1 1 0 0 1 1 0 1 1 1 1 0 0 1 0 0 0 0\n",
      " 0 1 0 0 1 0 1 1 1 0 1 1 0 1 0 1 0 0 0 0 0 1 0 1 0 0 0 1 1 0 1 0 0 0 0 1 1\n",
      " 0 1 1 1 1 0 1 1 1 0 0 0 0 0 1 1 1 0 1 1 0 0 1 0 0 1]\n"
     ]
    }
   ],
   "source": [
    "n_flips = 100 # Number of spin flips\n",
    "data = stats.randint.rvs(0, 2, size=n_flips)\n",
    "print(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I can count the \"heads\" by summing the array `data`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "48\n"
     ]
    }
   ],
   "source": [
    "n_heads = np.sum(data)\n",
    "print(n_heads)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Simulating many trials of 100 spin flips, to see how many times I get 60 or more heads\n",
    "\n",
    "Use a variable `counter` to keep track of the number of times I get 60 or more heads."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2840 0.0284\n"
     ]
    }
   ],
   "source": [
    "n_expts = 100000               # Number of experiments to simulat\n",
    "\n",
    "counter = 0\n",
    "for i in range(n_expts):\n",
    "    if np.sum(stats.randint.rvs(0, 2, size=n_flips)) > 59:\n",
    "        counter += 1\n",
    "\n",
    "print(counter, counter/n_expts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### I get 60 or more heads in about 2.8% of the trials"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Find average number of \"heads\"; should be close to 50"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using properties of binomial distribution to analyze the same experiment\n",
    "\n",
    "The probablity of obtaining 60 or more heads is \n",
    "\n",
    "$$1 - \\mbox{probability of obtaining 59 or fewer  heads} = 1 - C_{DF}(59)$$  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "probability = 0.02844396682049044\n"
     ]
    }
   ],
   "source": [
    "p = 0.5    # probability of getting a \"head\" in a single trial\n",
    "n = 100    # number of trials\n",
    "\n",
    "print(\"probability =\", 1 - stats.binom.cdf(59, n, p))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Agrees with result of Monte Carlo simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}