{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using the $\\chi^2$ test of a distribution: Is this a fair die?\n",
    "\n",
    "In this notebook I provide two data sets from the roll of a nominally standard 6-sided die.\n",
    "The question is: Is this a fair die, or not?  The data is an array of integers chosen \n",
    "from the set {0,1,2,3,4,5}. (We could make this set be {1,2,3,4,5,6} to correspond to the \n",
    "numbers on a die, but with python it's easier to start at 0.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import stats\n",
    "from scipy import optimize\n",
    "\n",
    "import matplotlib as mpl       # As of July 2017 Bucknell computers use v. 2.x \n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Following is an Ipython magic command that puts figures in the  notebook.\n",
    "%matplotlib notebook\n",
    "\n",
    "# As of Aug. 2017 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": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "xk = np.arange(6)\n",
    "# or xk = np.linspace(0,5,6, dtype=int)\n",
    "pk = np.ones(6)/6.\n",
    "custm = stats.rv_discrete(name='custm', values=(xk, pk))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.16466667 0.16466667 0.16466667 0.16466667 0.16466667 0.17666667]\n",
      "0.9999999999999999\n"
     ]
    }
   ],
   "source": [
    "d = 0.01\n",
    "pk2 = pk - d/5\n",
    "pk2[5] += d/5 + d\n",
    "custm2 = stats.rv_discrete(name='custm', values=(xk, pk2))\n",
    "print(pk2)\n",
    "print(sum(pk2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The data is an array of integers from 0 to 5, "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = custm2.rvs(size=50000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[8228 8183 8414 8304 8121 8750]\n"
     ]
    }
   ],
   "source": [
    "o = np.bincount(data, minlength=6)\n",
    "print(o)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "e  = len(data)*np.ones(6)/6."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "31.17112\n"
     ]
    }
   ],
   "source": [
    "chi2_data = np.sum((o-e)**2/e)\n",
    "print(chi2_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.351167146075195e-05"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "1 - stats.chi2.cdf(chi2_data, 6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.bincount?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0, 1, 2, 3, 4, 5])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.arange(6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0, 1, 2, 3, 4, 5])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linspace(0,5,6, dtype=int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt('die_50000.dat',data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt?"
   ]
  },
  {
   "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.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}