{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lab 5, Exercise 1: Airline Meal Part 2\n",
    "\n",
    "We now have access to the full details of the survey carried out by the airline, with the number of customers who chose a 'beef' type meal in a single 219-seat flight, $n_{b1}$ (for simplicity, we've created the data such that the number of customers in each flight is a constant). We will use this new data to assess the quality of the model inferred in part 1.\n",
    "\n",
    "Use the posterior probability distribution from part 1 to make a posterior predictive test, using the standard deviation in $n_{b1}$ as the test quantity:\n",
    "\n",
    "$T(n_{b1}) \\equiv \\sqrt{\\dfrac{1}{n_{flight}}\\sum_{i=1}^{n_{flight}} (n_{b1,i} - \\bar{n_{b1}})^2}$\n",
    "\n",
    "where $n_{b1,i}$ is the number of customers on the $i$-th flight of the survey who chose beef, and $\\bar{n_{b1}}$ is the average $n_{b1}$ in the survey.\n",
    "\n",
    "Under the binomial distribution model inferred in part 1, what is the probability of obtaining a more extreme value of $T(n_{b1})$?\n",
    "\n",
    "Discuss the accuracy of the model based on the answer to the question above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solution\n",
    "\n",
    "Let's first calculate the test quantity for the observed data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "\n",
    "# load the survey data\n",
    "f = open('meal_survey.txt', 'r')\n",
    "nb1_obs, nc1_obs = np.loadtxt(f, dtype=int, unpack=True)\n",
    "f.close()\n",
    "\n",
    "nflight = len(nb1_obs)\n",
    "\n",
    "nb = nb1_obs.sum()\n",
    "nc = nc1_obs.sum()\n",
    "\n",
    "# calculates the standard deviation in nb1_obs\n",
    "Tobs = nb1_obs.std()\n",
    "print 'Tobs = %2.1f'%Tobs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given $p_{beef}$, the parameter of our binomial model describing the probability of a random customer to choose beef, the posterior probability distribution for $p_{beef}$ found in Part 1 is proportional to a Beta distribution:\n",
    "\n",
    "$P(p_{beef}|n_b,n_c) \\propto p_{beef}^{n_b}(1-p_{beef})^{n_c} \\propto \\rm{Beta}(n_b+1,n_c+1)$\n",
    "\n",
    "Let's simulate new survey data from our model as follows:\n",
    "- Draw a value of $p_{beef}$ from the posterior\n",
    "- For $n_{flight}$ flights, generate 219 customer choices drawn from a binomial distribution with probability $p_{beef}$. Record the value of $n_{b1}$ for each simulated flight, then calculate $T(n_{b1})$.\n",
    "- Repeat for a large number of iterations\n",
    "- Count the fraction of times for which the simulated $T$ is more extreme than the observed one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pylab\n",
    "\n",
    "\n",
    "nsim = 1000\n",
    "\n",
    "nseat = 219\n",
    "\n",
    "Tsim = np.zeros(nsim)\n",
    "for i in range(nsim):\n",
    "    # draw a value of p_beef from the posterior\n",
    "    p_beef_here = np.random.beta(nb+1, nc+1, 1)\n",
    "    \n",
    "    nb1_mock = np.zeros(nflight, dtype=int)\n",
    "    # loop over 20 flights\n",
    "    for j in range(nflight):\n",
    "        x = np.random.rand(nseat)\n",
    "        nb1_mock[j] = (x < p_beef_here).sum()\n",
    "        \n",
    "    Tsim[i] = nb1_mock.std()\n",
    "\n",
    "print (Tsim > Tobs).sum()/float(nsim)\n",
    "\n",
    "# plot histogram of simulated Tsim vs. observed one\n",
    "pylab.hist(Tsim)\n",
    "pylab.axvline(Tobs, linestyle='--', color='k')\n",
    "pylab.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "The standard deviation in $n_{b1}$ observed in the survey is around 12. It's very unlikely to obtain such a large value, based on the posterior prediction from our model. This suggests that the model is not entirely accurate.\n",
    "\n",
    "(The survey data has been generated by assuming two different values of $p_{beef}$ depending on the flight. This could be the case if people's preference varies with the time of the day, or with place of departure)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}