simulate_bernoulli_trials.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are two sports team A and B. A match is played best 4 of 7 games (or 4/7 or 5/9). The probability that A wons a game is p1. Suppose the playing follows a markov chain process, P(A wins game k|A wins game k-1)=p11, P(A lost game k|A wins game k-1)=p01, how to simulate the probability that the 7th game is played? and how to construct the CI for that estimation. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import random\n",
    "import numpy as np\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def bernoulli(p):\n",
    "    return 1 if random.random() < p else -1\n",
    "def sim(p1, p01, p11):\n",
    "    result = [bernoulli(p1)]\n",
    "    while result.count(1) < 4 and result.count(-1) < 4:\n",
    "        if result[-1] == 1:\n",
    "            result.append(bernoulli(p11))\n",
    "        else:\n",
    "            result.append(bernoulli(p01))\n",
    "    return result\n",
    "def pGameSeven(N, p1, p01, p11):\n",
    "    res = [0]*N\n",
    "    for i in range(N):\n",
    "        trial_res = sim(p1, p01, p11)\n",
    "        if len(trial_res)==7:\n",
    "            res[i] = 1\n",
    "    return res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get a point estimate\n",
    "total_runs = 10000\n",
    "res = pGameSeven(10000,0.5,0.4,0.55)\n",
    "prob_full_game = sum(res)*1.0/len(res)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.2644"
      ]
     },
     "execution_count": 63,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "prob_full_game"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [],
   "source": [
    "# use bootstrap to build CI\n",
    "sample_means = []\n",
    "for i in range(total_runs):\n",
    "    sample = np.random.choice(res, total_runs-1)\n",
    "    sample_mean = sum(sample)*1.0/(total_runs-1)\n",
    "    sample_means.append(sample_mean)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x116bd81d0>"
      ]
     },
     "execution_count": 78,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXUAAAD4CAYAAAATpHZ6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmUW+d93vEvthnMgtkxOznD9ZIUKZJaKEoURcqWLcub\nnKSOkyZNEzfOUrVJXJ8ktaucHuU4p6c5rtKmqRPHS2Q3seNaibzIoSVrIUVKFBdxHZJzh9uQs+/7\niuX2D2AoiiY5Qwww9wJ4PufMOSAGAzy8AH548b7vfV+XZVmIiEhmcNsdQEREkkdFXUQkg6ioi4hk\nEBV1EZEMoqIuIpJBVNRFRDKIdyE3MgzjAeC/m6a52zCM1cBzgAU0AU+Zphk1DOMzwG8DYeCLpmm+\nON/99vWNOW4+ZWlpPkNDk3bHWBBlTY10yZouOUFZky0YDLhu9bt5W+qGYfwR8DXAH7/qWeBp0zR3\nAi7gScMwqoHfA3YAjwP/zTCM3MUGt4PX67E7woIpa2qkS9Z0yQnKupQW0v1yEfj56/59L7AvfnkP\n8BiwDXjTNM0Z0zRHgAvA3ckMKiIi85u3+8U0zX8yDKPxuqtcpmnOdZuMAcVAETBy3W3mrr+t0tJ8\nR34qBoMBuyMsmLKmRrpkTZecoKxLZUF96jeIXnc5AAwDo/HLN15/W07stwoGA/T1jdkdY0GUNTXS\nJWu65ARlTbbbfegkMvvluGEYu+OXnwD2A4eBnYZh+A3DKAbWExtEFRGRJZRIS/1zwFcNw8gBzgHP\nm6YZMQzjL4kVeDfwX0zTnE5iThERWYAFFXXTNFuB7fHLLcCum9zmq8BXkxlORETujE4+EhHJICrq\nIiIZREVdRCSDJDJQKpLx9p7ouOn1u7fULXESkTujlrqISAZRURcRySAq6iIiGURFXUQkg6ioi4hk\nEM1+kbRys1kpmpEi8i611EVEMoiKuohIBlFRFxHJICrqIiIZREVdRCSDqKiLiGQQFXURkQyioi4i\nkkFU1EVEMoiKuohIBlFRFxHJICrqIiIZRAt6ScYanwqx70QHZ1uHiESiRIH6igI2NJaxvrGUAr/P\n7ogiSaeiLhknHInyvdcvsvdEB6Fw9Nr1LhdcaB9h74lOvB4396+r5H331LGqrtjGtCLJpaIuGWVm\nNsKXv9/E6UsDVBT7eey+ZTy8qYZ8v5dINEpr1xhnWgc5eKaHg2e6OXimm7XLSvj4jkbWN5Ticrns\n/i+ILIqKumSMyekwf/G9E1zsGGXTynI2rizD53VzuLnn2m12b4m1zD/2UCPnrgzx8pE2Tl0c4Ev/\neAJjWQm//NgagsGAjf8LkcVRUZeMYFkWf7fnHBc7Rtm+oYpPf2Q9B0533fL2LpeLDY1lbGgso7V7\nlO/vv8ypiwM889wRPrS9kZICHz6v5hFI+lFRl4yw72Qn75h9rF1Wwm9+dANu98K7URqri/iDT26m\n6dIA33n1PHsOthLI97Hz7hoqSvJSF1okBdQUkbTX0TfOd145T4Hfy2997M4K+vU2riznmU9v41+9\nbw1jkyH2HLpK85WhJKcVSS0VdUlrUcviay+eIxSO8ukPr6esyL+o+/N63Pzbj2zgg/cvI9fn4fC5\nXo6f78eyrCQlFkktdb9IWrvYMcqVnjFW1hYxMjl7042pE1Fdns8T25fzytF2Tl8cYGY2zAMbqpJy\n3yKppKIuaWs2HOF4Sx9ej4utaysW9DcLKfqBwlhrP5Cfw4ceiBX2lrYR/DleHt1av6jMIqmm7hdJ\nW00XB5mejXDXirKUnR2al+vlsfvqKczzceriAPtPdqbkcUSSRUVd0tL4VIizrUPk+73ctaIspY+V\nl+vl/ffWk+Nz882fmJhXNXgqzqWiLmmp6dIgUcti65oKvJ7Uv4yLC3N4dGsdAH/7o7OMT4VS/pgi\niUioT90wDB/wTaARiACfAcLAc4AFNAFPmaYZvcVdiCRscjrEhfYRCvN8rKgpWrLHrSrL5+MPN/L9\n/Zf51k+a+d1PbNSyAuI4iTZxPgx4TdN8CPhT4M+AZ4GnTdPcCbiAJ5MTUeS9zlweImpZbFpZlvCc\n9ER99MFG1tYXc9Tsu+0ZqyJ2SbSotwBewzDcQBEQAu4F9sV/vwd4bPHxRN5raiZMS9swBX4vK21Y\nXdHtdvGZj91FXq6H7756gdGJ2SXPIHI7iU5pHCfW9dIMVAAfBR4xTXPuDI0xYN53XGlpPl6vJ8EI\nqZNOCzplW9aLnWNEohb3rq+ipCh1p/DPTWu8UTAYIBgM8Gsf3sBXXjjNi29f5fd/aWvKcswn257/\npZJOWW+UaFH/LPCSaZqfNwxjGfAakHPd7wPA8Hx3MjQ0meDDp04wGKCvb8zuGAuSbVlnQhGaLvWT\n6/NQX5HP2Ph0ktK9V6DQf8v7nvs/3LemnH+pLOSVI1fZti7Iahu+NWTb879U0iHr7T50Eu1+GQJG\n4pcHAR9w3DCM3fHrngD2J3jfIjf19pluZkNR1i4rXpIZL7fjcbv51Q+uBeDvXzKJRrWMgDhDou+M\nvwDuMQxjP7FW+heAp4BnDMM4SKzV/nxyIorEltZ95Wg7LhcYy0vtjgPAmvoSdmys5mrvOG+f7bY7\njgiQYPeLaZrjwC/e5Fe7FhdH5OaarwzR0T9BY02AfL9zVrf4xM6VHDrXww8OXGbb+irbv0GIOOfd\nIXIbPz3aDsD6Bntb6TdbO2bX5jpePdbOm6e72LWlzoZUIu9Ss0Icb3B0mpMX+llREyDowE0rPvJQ\nAz6vmx++2UooHLE7jmQ5FXVxvDdPd2EBj2yutTvKTZ240M+a+mKGxmb4+o/PJW35X5FEqKiLo0Ut\niwOnu8jxudm23rnrmd+1ogyP28XZ1iHNhBFbqU9dHGvviQ66BybpG55mVW0Rh8712B3plvJyvayq\nK6KlbYQrPc6e4yyZTS11cbQLHbHTIVbXL/3JPXdqQ2NsCeCzlwe1/Z3YRkVdHGs2FOFK9xiBfB+V\npc4bIL1RUUEOy6sKGRidoaVt3hOqRVJCRV0c63JXbJ2X1XXFabPE7V3x1vqeQ1dtTiLZSkVdHOtC\nxwguYJUN66okKliaR7DEz6mLA/Q6cG0jyXwq6uJI7b3jDIxMUxsscNQZpAsxt4zB3hPaz1SWnoq6\nONL+U7ENKNakwQDpjRqqCynM83HgVJdORpIlp6IujhOORDl4pht/joe6YKHdce6Yx+1m5+YaxqdC\nHGnutTuOZBkVdXGcE+f7GZ8KsbK2CM8Sb1eXLLu31OECXj+us0tlaamoi+PMdb2kw9z0WwmW5LFx\nZTkXO0a5qpORZAmpqIujDI5O03R5gFW1RZQU5todZ1Ee3RpbsXGvWuuyhFTUxVHebOrGsuDhu2vs\njrJod68qp7wol4NnepiaCdsdR7KEiro4RtSyOHCq0/GLdy2U2+1i15Y6ZkIR3mrSzkiyNFTUxTFa\nrg7TNzzN/UYlebnpNTf9VnZursXjdrH3eIfWg5ElkRnvHMkIcwOkmdD1Mqe4IId7jSCHz/Vyvn2E\ntctKfuY2N1t/fbd2UJIEqaUujjA5HeYds5fK0rybFr50Njdg+tqxdpuTSDZQURdHOHS2m9lwlJ13\n16TN4l0LtXZZCbUVBbxj9jE6OWt3HMlwKuriCPtOduJ2udixKXO6Xua4XC52ba4lErV483SX3XEk\nw6moi+1au0e52jPO5tXlaT83/VYe2lSNz+tm34lOohowlRTSQKnY7o34aoa7tjhzY+lE3Gzwc1ll\nIZc6R2m+MnRtlySRZFNLXWw1PRvm7bM9lBXlsnFFud1xUmpuAFhL8koqqaiLrY6c62V6NsLDm2pw\np+niXQsVLPFTHyzgeEsfI+MzdseRDKWiLrZ642QnLmDn3ZnT9XIrLlfsDNNI1OKABkwlRVTUxTbt\nfeNc7Bxl48pyyov9dsdZEg/eVU2OTwOmkjoq6mKbuQHSRzZnfit9Tr7fy7b1VfSPTHP28qDdcSQD\nqaiLLULhCAfPdFNUkMPm1Zk9QHqjuSUANGAqqaCiLrY4avYxMR1mx6ZqvJ7sehmuqAmwvLKQE+f7\nGRrTgKkkV3a9m8QxsrHrZY7L5WLX1rprSw2LJJOKuiy57sFJzLZh1i0voao03+44tti+oYrcHA/7\nTmrAVJJLZ5TKktsfb50+Ej+D9GZnX2a6vFwv2zdUse9EJ519E9RXFtodSTKEWuqypMKRKG+e7qbA\n7+XetUG749hqbsC0pW3Y5iSSSRJuqRuG8Xng40AO8GVgH/AcYAFNwFOmaUaTkFEyyLdfaWF0YpZ1\nDSW8meVbvDVUB2isDnCle4yJqRAFeT67I0kGSKilbhjGbuAhYAewC1gGPAs8bZrmTsAFPJmkjJJB\nzrePALCmPrM2wkjU7q11WLx7XEQWK9GW+uPAaeAFoAj4Q+AzxFrrAHuAD8Z/LwLA4Og0nX0TVBT7\nKQ1k5hK787lx/CAUjuLzuLnQPsLdq8ozfv0bSb1Ei3oF0AB8FFgB/BBwm6Y5N4w/BhTPdyelpfl4\nvZ4EI6ROMBiwO8KCpVPW45cGsYBNqysIFDp7WYClzGc0lNJ0aYDB8VlW1MbeNgt9XtPp+VfWpZFo\nUR8Amk3TnAVMwzCmiXXBzAkA847+DA1NJvjwqRMMBujrG7M7xoKkU9by8kJ+8lYrXo+L6tI8xsan\n7Y50S4FC/5Lma6wupOnSACdb+qgoin2DWcjzmk7Pv7Im1+0+dBKd/XIA+JBhGC7DMGqBAuDVeF87\nwBPA/gTvWzLQyfN9DIxO01hThM+rSVfXKyvyU1Hsp6N/gvHJkN1xJM0l9O4yTfNF4DhwGPgR8BTw\nOeAZwzAOEpsR83yyQkr6e+nQFQDW1M/bK5eV5jbQON+u6Y2yOAlPaTRN849ucvWuRWSRDDU6Ocuh\npi7qggVUZMkSu3eqsSbAkeZeLnSMsHl1hd1xJI3pe7Ck3MGmbsIRi0fursXl0uyOm/F63KyqLWJq\nJkJ737jdcSSNqahLSlmWxYFTXXg9bh7cWG13HEdbVRfrmrrc5exBOnE2FXVJqas943T0T3D/hioK\ndcbkbZUV5VJckEN77zhTM2G740iaUlGXlNl7ooP/9/oFAEoKc7Ny4a474XK5WFETIBK1ONbSZ3cc\nSVMq6pIy0ajF5a5RcnxuGmrS92SOpbSitgiAt89k97o4kjgVdUmZroFJpmcjNFYX4XHrpbYQgfwc\nKor9nL0yxMi4dkWSO6d3mqTMpc7YIlWr4q1PWZgVtUVYFhxu7rU7iqQhbZIhi3azvvJQOMrVnnEC\n+T4qSjQ3/U40Vgc42tzLobM9fOC+ZfP/gch11FKXlGjrHScStVhRU6S56XcoL9fLuuWlXOocZWDE\nuWvkiDOpqEtKtHaNArBCA6QJuX9dJQBHTXXByJ1RUZekm5mN0Nk/QWkgl+LC7Fw3fbHuMYK4XHBE\n/epyh1TUJemu9o4RtdRKX4yi/BzWN8S6YPpHpuyOI2lERV2Sbu4098ZqzXpZjPvmumCadSKSLJyK\nuiTV1EyYnoFJgiV+CvO1LMBi3LM2iNvlUheM3BEVdUmqK91jWKiVngxF+Tmsayjhctco/cPqgpGF\n0Tx1Saor3bGul4Zq9acvxtzc/6KCHAC++/oF7lpRxu4tdXbGkjSglrokzdRMmJ6hKYIleeT71V5I\nhuVVhbhc735YisxHRV2Spq0ntrlDQ3WhzUkyhz/HS3VZPv0j09q/VBZERV2S5kpPrDW5vEpdL8nU\nGO/Kmju+Irejoi5JMT0boXtwkvJivzbDSLJl8S6YVnXByAKoqEtStPeOY1nQUKWul2Sb64IZGJnW\nLBiZl4q6JIW6XlJrrgvmiNaCkXmoqMuizYYidMXXepmbgifJtawqEFsL5pyKutyeirosWnvfBFFL\nc9NTyZ/joaY8n9buMfrUBSO3oaIuizY3h3q5+tNTau5D86iWDZDbUFGXRZmeDdPZP0FxQQ4lWmY3\npZZVBvC4tRaM3J6KuizK6UuDRKKWul6WgD/Hw/qGUlq7x+hVF4zcgoq6LMpcV4C6XpbG3I5I76i1\nLregoi4Jmw1FOHVxgEC+j9KAul6Wwta1QTxuF4dV1OUWVNQlYU2XB5kJRVheFdDm0kukMM/H+sZS\nrnSP0T0wYXcccSAVdUnYO/ETYbSA19Ka64I5cLLT5iTiRCrqkpBQOMqJCwOUF+VSXuS3O05W2bom\n1gVz4GSH3VHEgVTUJSHnrgwyNRPmXqNSXS9LrDDPx4bGMi62j9A7NGl3HHEYFXVJyFEzthnyvUbQ\n5iTZaa4LRnPW5UYq6nLHwpEox1v6KC7MYVVdsd1xstLWtRV4PS6ONvfZHUUcZlF7jhmGUQm8A3wA\nCAPPARbQBDxlmmZ0sQHFecy2YSamw7zvnjrc6nqxRYHfx5a1lRw910PP0CRVpfl2RxKHSLilbhiG\nD/gKMHdq27PA06Zp7gRcwJOLjydO9M61rpdKm5Nkt51bagE42NRtcxJxksV0v3wJ+Btgbl7VvcC+\n+OU9wGOLuG9xqGjU4lhLH4V5PtYuU9eLnR7cVEuuz8NbTd1ELcvuOOIQCXW/GIbx60CfaZovGYbx\n+fjVLtM0515ZY8C87/jS0ny8Xk8iEVIqGEyfdUyWOmvTxX5GJ2Z5fHsD1VWxpzhQuLApjQu9nRM4\nNeuNz/fDW2p59UgbvWOzbFpVYVOqhdH7amkk2qf+acAyDOMxYAvwLeD67+IBYHi+Oxly4HSsYDBA\nX1967AVpR9ZXDl0B4K6GkmuPPTY+Pe/fBQr9C7qdEzg56/XPdzAY4N7VFbx6pI0f779IdZFzl2rQ\n+yq5bvehk1D3i2maj5imucs0zd3ACeDXgD2GYeyO3+QJYH8i9y3OFbUs3jF7KfB7Wbe81O44Aqxd\nXkJFsZ+jzX1Mz4btjiMOkMwpjZ8DnjEM4yCQAzyfxPsWB7jUOcrw+Cxb1wTxejQb1gncLhcPbaxm\nJhS5NoAt2W1RUxoB4q31ObsWe3/iXHPL7ObkuNl7Qqeo2+H64z7XTfTQphp++GYr+091sWNTjY3p\nxAnU3JIFseJdLz6vm5pyzYl2ksqSPNY3lNLSNkyXVm7MeirqsiCt3WMMjM6wrLIQj1svGyfZe6KD\nYElsts4//LRF36KynN6dsiBHTe1w5GTLqgL4czxc7BglEtWJ3NlMRV3mZVkW7zT3kZvjobaiwO44\nchMet4tVdUXMhCJc7Rm3O47YSEVd5tXWO07v8BSbV5Vr1ouDrakvAeB824jNScROeofKvOa6Xu7T\nWi+OVlSQQ3VZPt2Dk3QPOu/EPlkaKupyW5ZlcbS5jxyvm00ry+2OI/NYE1+P540T2uouW6moy211\n9k/QPTjJplXl5OY4b50eea/lVYXk+jwcON1FKKwB02ykoi63NbezjnY4Sg8et5tVdUWMT4U4fl5n\nmGYjFXW5pdePt7P3eAcet4uxyZDmP6eJuQHTfeqCyUoq6nJLg6MzjE6GWFZZiM+rl0q6KC7MwVhW\nwrkrQ/RowDTr6J0qt3S5axSAxpr0XVs6W+3aGtsVSa317KOiLjcVtSxau8bI8bqpC+qEo3Rz79pK\nAvk+9p/qZDYUsTuOLCEVdbmp823DTM6EWV4d0FovacjndfPI5lompsMcPtdrdxxZQnq3yk29fbYH\ngBXqeklbu7fU4XLBa8fa7Y4iS0hFXX5GKBzlaHMvebkeqsq0zG66Ki/2s2V1Ba3dY1zqHLU7jiwR\nFXX5GcfP9zExHWZFTRFul8vuOLIIj95TB8Draq1njUXvfCSZ58DpLgBW1xfbnEQSNXdOgWVZBPJ9\nHDzbQ11lAR/a1mBzMkk1tdTlPQZHpzlzaZCVtUWUFDp3d3pZGJfLhbG8hGjU4kK7Vm/MBirq8h5v\nNXVjAQ/frb0uM8XqumI8bhctbSNEo5bdcSTF1P0i11iWxYHTXeR43WxbV8Xh5h67I0kS5Pg8rKgt\n4kL7CN9+pYX6yvfuXrV7S51NySQV1FKXa1rahukdmuIeI0i+X5/3mWTd8th6MObVYZuTSKqpqMs1\nrx6LDa7t2lxrcxJJtrIiP8ESPx39E4xNztodR1JIRV2A2ADpMbOP+mAha5eV2B1HUsBYXgqotZ7p\nVNQFgNePdxC1LB67rx6X5qZnpIbqQvw5Hi50jBCOaAONTKWiLoTCEfad6KTA72X7hiq740iKeNxu\n1tQXMxuK0to1ZnccSREVdeHQ2V7Gp0I8sqWWHJ+2rMtka5aV4EJdMJlMRT3LRS2Ll49cxeWCR7dq\nalumK8zzUV9ZyMDoNP3DU3bHkRRQUc9yJ8/30943wfYNVVQU59kdR5aAEZ/e2KzWekZSUc9ilmXx\nw7dacQEffrDR7jiyRGrK8ykuyKG1a5TJ6bDdcSTJVNSzWNPlQa50j3GvEaSuQrsbZQuXy8X6xlKi\nFphXh+yOI0mmop6lLMviR2+1AvDRhxptzSJLb2VtEbk+D2bbMDOz2u4uk6ioZ6mmy4NcaB9hy+oK\nlldpd6Ns4/W4MZaXMBuK8mZTl91xJIlU1LNQNGrxvdcv4AJ+7pGVdscRmxjLS3C7Xbx8pI2opdUb\nM4VWbcpCbzV10943wY5N1SyrLLy2oYJkl7xcLyvjqzeevNDP1jVBuyNJEiRU1A3D8AHfABqBXOCL\nwFngOcACmoCnTNPUucgO88o7bXx//2U8bhfVZfkq6FluQ0MpF9pHePlwm4p6hki0++VXgQHTNHcC\nHwL+CngWeDp+nQt4MjkRJZnOtQ4xOR1mfUMpBXk+u+OIzUoCuWxcUYbZNkxrtzanzgSJFvXvAX8S\nv+wCwsC9wL74dXuAxxYXTZJtYGSa05cG8Od42LiyzO444hAf3LYMgJePtNmcRJIhoe4X0zTHAQzD\nCADPA08DXzJNc260ZQyYd9fi0tJ8vF7nrTUSDKbPbJA7yfr1Pc2EIxa7ttZSXrr089IDhf4lf8xE\npUvWZOTcfX8Dz++7xJFzvfz2z2+moiQ1ZxZn6vvKaRIeKDUMYxnwAvBl0zS/bRjGn1/36wAw7znI\nQ0OTiT58ygSDAfr60mMFuzvJeqZ1kDdPdhIs8VNbnsfY+HSK071XoNC/5I+ZqHTJmqyc/f3jvH9r\nHX+3p5nv/OQcv/T+NUlI916Z+r6yy+0+dBLqfjEMowp4Gfhj0zS/Eb/6uGEYu+OXnwD2J3Lfknyh\ncJRv/7QFF7BtQ5XWS5efsf2uakoDuew90cHohHZGSmeJ9ql/ASgF/sQwjL2GYewl1gXzjGEYB4Ec\nYt0y4gB73r5C18Akj95TR3lRenQryNLyed18eHsDs6EoLx25anccWYRE+9R/H/j9m/xq1+LiSLJ1\nDUzw4sFWSgpz+IVdqzh0rsfuSOJQj2yu4cWDrbx2rIMnHmigULOj0pLOKM1glmXxf18yCUcsfuUD\na8nL1blmcms+r4cnti1nZjaimTBpTEU9g33txbM0Xx2mvrKQ0clZnWgk89q1tY6ifB+vHG1jbFJ9\n6+lITbcMNTAyzdHmPnxeNw9sqNTgqNzSjR/2xvJSjjT38uODV1IyE0ZSSy31DGRZFs/tOUcoEuX+\ndZUU+NU3Kgu3dnkxFcV+XjvWzsCI86d2ynupqGegN052cqZ1iLqKAlbVFdkdR9KMx+3mEztXEI5Y\n/ODAZbvjyB1SUc8wAyPTfPe1C+Tletm+UXPSJTHbN1RTFyzgzaYu2vvG7Y4jd0BFPYPMdbtMz0b4\npfevVreLJMztdvHJ3auxLPjOK+extN562tBAaQaYG+hqaRu+1u0SjkTVSpeEzb2m6oIFnLsyxLde\nMmmoDrB7S53NyWQ+aqlniLHJWd6Jz3ZRt4sky/3rKnG74GhzL+GItkdIByrqGSAatdh/sotQJMq2\n9ZrtIslTVJDD+sZSJqbDnLk8aHccWQAV9Qxw/Hw//SPTrKgJsLJWs10kuTatKicv18vpi4N09k/Y\nHUfmoaKe5s5cHuTM5UEC+T6231WtbhdJuhyvhwc2VBK1LJ7b06xNqh1ORT2NDY/N8LUXz+Jywc7N\ntfi8ejolNZZXBWioDnChY4TXj2m5CSdTFUhTUcviL/7xGCMTs9yzNkhFsZbUldSKjdd4eX7vRXoc\nuMGNxKiop5m9JzrYe6KDv/l+E8eae6mtyGdDY6ndsSQL5OV6+ZUPrmUmFOErPzij2TAOpaKehnoG\nJ3mnpY+8XC87NtWoH12WzPYN1ezYWE1r9xj//MYlu+PITaiop5mJ6RD7TnQC8Pj2Bq2RLkvuX39g\nLZWlefzk0FVOXRywO47cQEU9jYTCEfYd72R6NsJ9RiV1wUK7I0kWysv18jtP3oXX4+YrP2zSNEeH\nUVFPE1HL4hv/0kz/yDQra4tY11BidyTJYo3VRfzGh9cxNRPhL58/xfhUyO5IEqeinia+v/8Sh872\nECzx8+BdWgZA7DE3UL/3RAczoQibVpbROzzFX/3TKWZCEbvjCSrqaeGNk528+NYVKkvzePSeOjwe\nPW3iDFvWVHD/ukpa2kf4q386RSiswm43VQeHO3imm2/uaaYwz8dnP7kZf44GRsU5XC4Xn/nYBras\nruBM6xD/54UmFXabqag72JHmXr724lnycr187lNbqCrLtzuSyM/wetz87ic2snFlGacuDvDsd08y\nOa0+druo2edQ+0918s09Jrk+D5/91GYaqgN2RxK5qbm11zevKmd0YhazbZinv3aIx+6rJz++Ymig\n0M/Y+LTWY18Caqk7jGVZ/OjNy/zdvzSTl+vhP31qC6tqi+2OJTIvj8fNI1tqMZaXMDw+y48PXqFv\neMruWFlHLXUHmZmN8K2Xmjl4pofyolwevruG9r5x7REpacPtcrFtfSWBPB/vmH28dKiNBzZUcs/6\narujZQ0VdYfo7J/gr7/fREf/BCtqivgPP7+Jkxf77Y4lcsdcLhcbVpRREsjljZOdHDzTw+DYLPes\nrbA7WlZQUbdZOBLlpcNX+cGBVsKRKO+/t55PvW81Xk1blDRXW1HARx5sYP/JLsyrQ3T1j7O6rpgV\nNdrIJZVU1G1kXh3iO6+e52rPOP4cDw/fXUtdsIADp7vsjiaSFIH8HB5/YDlnWoc40dLHF791lMfv\nX86TO1dZC4x6AAAH+ElEQVSQ6/PYHS8jqajboL13nH9+4xInLsS6V1bVFXGfUUlujl7kknk8bhc7\n7q6lstjPifP9/OTwVQ439/DxHSvYsakaj1vfSpNJRX2JWJZFS9swe65b2W5NfTG/+L7VtPVqIFQy\nX3V5Ps88so0fHrjMT4+289yeZn58sJWHN9XwwF3VVJbkXbvt3DTJG2lK5PxU1FNsdGKWQ+d62Hei\n89pqdsGSPDatLKMuWKCCLlkl1+fhk4+u5rH7lvHiW63sP9XFC/sv88L+y1SV5tFYU0R9sIDe4SkK\n83wU5vnw53i01tEdUFFPgdGJWY619HGkuZfmq0NYVuwr6Lb1lZQW5VJVqjNDJbuVBnL5N48b/MKu\nVdfeKxc7Rjh0todDN9zW7XZR4PdSmOfjas84q2qLWF1fTGVJnor9TaioJ0E4EqW1a4zmq0OcbR3E\nbBtmbsP1imI/jdUBVtQWaUMLyXo361bZvaWOh++uwbIseoen6Oyf4O2zPUxMhRifCjExFWZiOkTX\nwCRdA5PsPR67j6J8H6vrS1hbX8yGxtg3XxV5FfU7FolG6R2KvfA6+iY43z6M2TZMOGJdu02wxE9D\ndYCGqgAFeT4b04qkD5fLRVVpPlWl+YxMzP7M70PhKGvqi7nQPsKFjhHOt49wrKWPYy19AOTleqgp\nL6CmPJ+f27mSsqLs3IzdZVnW/LdaIMMw3MCXgc3ADPCbpmleuNXt+/rGkvfgSVJSms/5ywMMjc0w\nODbN0NgMQ6MzDI3N0DM0RffgxHsKOEBJYQ5VZflUl+VTVZa3ZCspzq2nkQ6UNfnSJSekJqtlWUxM\nh+kemKRrYIKugUmmZ99dIbK8KJdVdcUsqyxkWWUhwZI8yor8806lDAYD9PWNJTVrsgWDgVt+JUl2\n9fkE4DdN80HDMLYD/wN4MsmPAcT26gyHo0St2JMbjVpEAStqEbWs2PVRi1AkysxshJnQuz8TU2HG\np0KMT80yNhlibCrE+GQoft2tV5fzelzUBwupqyigNlhAXUUBDdVFHD/fl4r/oojchsvlojDPx+r6\nYlbXF2NZFsPjM3T2TxKNWpxvH+HwuV4On+t9z9/l5XrI8XnIve7H63Hhdsd+/Lk+ugcmcLtij+F2\nxda1aagM4PO6r/3keN34fB58Hjc5Pjc+jxufz02O13PD7Ty43bEaPFeJvR7XtcXOki3ZRf1h4CcA\npmm+bRjGfUm+fwBeO9bO37/ckrT7cwG5ObEnty5YSK7PTb7fS77fS4HfR35u7PKNo/CDYzMMjqmg\niziBy+WiNOCnNBDrdtm4sozxqRBDYzMMj88yPhUi1+tmdDJ0rYE3NjnL9GyEhXRYXGgfSWref/eR\n9ezYVJPU+4TkF/Ui4Pr/ecQwDK9pmuGb3fh2XyFu51OPr+dTj69P5E9FRDJask/lGgWuX/jbfauC\nLiIiyZfsov4m8GGAeJ/66STfv4iI3Eayu19eAD5gGMZbxLqqfyPJ9y8iIreR1CmNIiJiLy2PJiKS\nQVTURUQySMYvEzDfWa6GYfwy8AdAmNjA7r83TTNqGMYxYrN5AC6bpvkbhmGsBp4DLKAJeMo0zagD\nc24FXgTOx6//a9M0v5uMnIvM+nng40AO8GXTNL+eymOagqyOO67ArwG/Hr+JH9gCVAMVOOy43ibr\nCpx3XD3AN4FGIAJ8xjTN5lS/XpMhG1rq185yBf4zsbNcATAMIw/4IvCoaZo7gGLgo4Zh+AGXaZq7\n4z9zA77PAk+bprmT2EBwMs+WTWbOe4Fnr7s+aW+QRWTdDTwE7AB2Acvif5LKY5rsrI47rqZpPjeX\nB3gH+D3TNIdx4HG9TVbHHVdis/i8pmk+BPwp8GfxP0n1cV20bCjq7znLFbj+LNcZ4CHTNCfj//YC\n08Q+0fMNw3jZMIzX4tMzIfbi2xe/vAd4zME5P2IYxhuGYXzdMIzrzx2wK+vjxFpBLwA/ItYym8ua\nqmOaiqxOO64AxM/evss0zb+9LqvTjuvtsjrtuLYA3ngrvwiYWz8k1cd10bKhqN/0LFcA0zSjpmn2\nABiG8R+BQuCnwCTwJWJv7t8B/iH+Ny7TNOemC40R+1R3Ys7DwB+apvkIcAn4r0nMmWjWCmJvpk9e\nl9VFao9psrM68bjO+QLwzHX/duJxvVVWJx7XcWJdL83AV4G/jP9tqo/romV8nzrznOUa/yT+c2At\n8AumaVqGYbQAF+JPXothGANADXB931kAGHZozhfiX2sh1tr830nMmWjWAaDZNM1ZwDQMYxoIktpj\nmuysjjuu8etLAMM0zdev+1vHHdfbZHXccTUM47PAS6Zpft4wjGXAa4ZhbCL1x3XRsqGlPt9Zrl8h\nNmjzieu+gn2aeL+bYRi1xD7pu4Dj8f5WgCeA/Q7N+ZJhGNvit3k/sf7LZEok6wHgQ4ZhuOJZC4AB\nUntMk53ViccV4BHg1Rtu68TjequsTjyuQ7zbuh8EfMQGT1N9XBct408+um7k+27ePcv1HmJfs47G\nf/YTG80G+F/Aj4mNcC+PX//Hpmm+ZRjGWmJfxXKAc8RGxN9dwNk5Oe8h1toJAd3Ab5mmOUqSJJLV\nNM0XDMP4c+BRYo2JL5im+VIqj2kKsjr1uP4hEDJN839ed19OPa43y+q440qsC+YbxL755sTzfzvV\nxzUZMr6oi4hkk2zofhERyRoq6iIiGURFXUQkg6ioi4hkEBV1EZEMoqIuIpJBVNRFRDKIirqISAb5\n/12adzdlCKeoAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x116b4b9d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.distplot(sample_means)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(0.2557255725572557, 0.27302730273027304)\n"
     ]
    }
   ],
   "source": [
    "# 1. direct use quantile\n",
    "# For example, if we were interested in a confidence interval of 95%, \n",
    "# then alpha would be 0.95 and we would select the value at the 2.5% percentile as the lower bound \n",
    "# and the 97.5% percentile as the upper bound on the statistic of interest.\n",
    "alpha = 0.95\n",
    "ordered = sorted(sample_means)\n",
    "lower = np.percentile(ordered, (1-alpha)*100.0/2)\n",
    "upper = np.percentile(ordered, alpha*100+((1-alpha)*100.0/2))\n",
    "print(lower,upper)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(0.25572742560156447, 0.27310989813080888)\n"
     ]
    }
   ],
   "source": [
    "# 2. using std\n",
    "# according to central limit theorem, the sample mean should follow gaussian distribution. \n",
    "# and the 95% confidence CI for sample_mean is mean(t)+/-1.96*se (10000 sample sizes, treat as std known, using z*)\n",
    "# se here is the std of sample mean distribution, the std of a estimator\n",
    "se = np.std(sample_means) \n",
    "# or se = np.std(res)/np.sqrt(total_runs)\n",
    "mean_value = np.mean(sample_means)\n",
    "lower = mean_value - 1.96*se\n",
    "upper = mean_value + 1.96*se\n",
    "print(lower, upper)"
   ]
  }
 ],
 "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
}