notebook/ch07_hypothesis_and_inference.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from probability import normal_cdf, inverse_normal_cdf # 이전장에서 작업한 확률 코드를 임포트 해야함\n",
    "import math, random"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def normal_approximation_to_binomial(n, p):\n",
    "    \"\"\"finds mu and sigma corresponding to a Binomial(n, p)\"\"\"\n",
    "    mu = p * n\n",
    "    sigma = math.sqrt(p * (1 - p) * n)\n",
    "    return mu, sigma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(500.0, 15.811388300841896)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mu_0, sigma_0 = normal_approximation_to_binomial(1000,0.5)\n",
    "\n",
    "mu_0, sigma_0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "normal_probability_below = normal_cdf\n",
    "\n",
    "# it's above the threshold if it's not below the threshold\n",
    "def normal_probability_above(lo, mu=0, sigma=1):\n",
    "    return 1 - normal_cdf(lo, mu, sigma)\n",
    "    \n",
    "# it's between if it's less than hi, but not less than lo\n",
    "def normal_probability_between(lo, hi, mu=0, sigma=1):\n",
    "    return normal_cdf(hi, mu, sigma) - normal_cdf(lo, mu, sigma)\n",
    "\n",
    "# it's outside if it's not between\n",
    "def normal_probability_outside(lo, hi, mu=0, sigma=1):\n",
    "    return 1 - normal_probability_between(lo, hi, mu, sigma)\n",
    "\n",
    "def normal_upper_bound(probability, mu=0, sigma=1):\n",
    "    \"\"\"returns the z for which P(Z <= z) = probability\"\"\"\n",
    "    return inverse_normal_cdf(probability, mu, sigma)\n",
    "    \n",
    "def normal_lower_bound(probability, mu=0, sigma=1):\n",
    "    \"\"\"returns the z for which P(Z >= z) = probability\"\"\"\n",
    "    return inverse_normal_cdf(1 - probability, mu, sigma)\n",
    "\n",
    "def normal_two_sided_bounds(probability, mu=0, sigma=1):\n",
    "    \"\"\"returns the symmetric (about the mean) bounds \n",
    "    that contain the specified probability\"\"\"\n",
    "    tail_probability = (1 - probability) / 2\n",
    "\n",
    "    # upper bound should have tail_probability above it\n",
    "    upper_bound = normal_lower_bound(tail_probability, mu, sigma)\n",
    "\n",
    "    # lower bound should have tail_probability below it\n",
    "    lower_bound = normal_upper_bound(tail_probability, mu, sigma)\n",
    "\n",
    "    return lower_bound, upper_bound"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(469.01026640487555, 530.9897335951244)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "normal_two_sided_bounds(0.95, mu_0, sigma_0) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def two_sided_p_value(x, mu=0, sigma=1):\n",
    "    if x >= mu:\n",
    "        # if x is greater than the mean, the tail is above x\n",
    "        return 2 * normal_probability_above(x, mu, sigma)\n",
    "    else:\n",
    "        # if x is less than the mean, the tail is below x\n",
    "        return 2 * normal_probability_below(x, mu, sigma)   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.06207721579598857"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "two_sided_p_value(529.5, mu_0, sigma_0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def count_extreme_values():\n",
    "    extreme_value_count = 0\n",
    "    for _ in range(100000):\n",
    "        num_heads = sum(1 if random.random() < 0.5 else 0    # count # of heads\n",
    "                        for _ in range(1000))                # in 1000 flips\n",
    "        if num_heads >= 530 or num_heads <= 470:             # and count how often\n",
    "            extreme_value_count += 1                         # the # is 'extreme'\n",
    "\n",
    "    return extreme_value_count / 100000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.06121"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "count_extreme_values()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p_hat = 525 / 1000\n",
    "mu = p_hat\n",
    "sigma = math.sqrt(p_hat * (1 - p_hat) / 1000)\n",
    "\n",
    "lower, upper = normal_two_sided_bounds(0.95, mu, sigma)\n",
    "\n",
    "lower < 0.5  < upper"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p_hat = 540 / 1000\n",
    "mu = p_hat\n",
    "sigma = math.sqrt(p_hat * (1 - p_hat) / 1000)\n",
    "\n",
    "lower, upper = normal_two_sided_bounds(0.95, mu, sigma)\n",
    "\n",
    "lower < 0.5  < upper"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "upper_p_value = normal_probability_above\n",
    "lower_p_value = normal_probability_below    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def run_experiment():\n",
    "    \"\"\"flip a fair coin 1000 times, True = heads, False = tails\"\"\"\n",
    "    return [random.random() < 0.5 for _ in range(1000)]\n",
    "\n",
    "def reject_fairness(experiment):\n",
    "    \"\"\"using the 5% significance levels\"\"\"\n",
    "    num_heads = len([flip for flip in experiment if flip])\n",
    "    return num_heads < 469 or num_heads > 531"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[508, 491, 484, 483, 502, 499, 483, 511, 481, 538]"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "random.seed(0)\n",
    "experiments = [run_experiment() for _ in range(1000)]\n",
    "\n",
    "[sum(experiment) for experiment in experiments][:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "46"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "num_rejections = len([experiment for experiment in experiments if reject_fairness(experiment)])\n",
    "num_rejections"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def estimated_parameters(N, n):\n",
    "    p = n / N\n",
    "    sigma = math.sqrt(p * (1 - p) / N)\n",
    "    return p, sigma\n",
    "\n",
    "def a_b_test_statistic(N_A, n_A, N_B, n_B):\n",
    "    p_A, sigma_A = estimated_parameters(N_A, n_A)\n",
    "    p_B, sigma_B = estimated_parameters(N_B, n_B)\n",
    "    return (p_B - p_A) / math.sqrt(sigma_A ** 2 + sigma_B ** 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def B(alpha, beta):\n",
    "    \"\"\"a normalizing constant so that the total probability is 1\"\"\"\n",
    "    return math.gamma(alpha) * math.gamma(beta) / math.gamma(alpha + beta)\n",
    "\n",
    "def beta_pdf(x, alpha, beta):\n",
    "    if x < 0 or x > 1:          # no weight outside of [0, 1]    \n",
    "        return 0        \n",
    "    return x ** (alpha - 1) * (1 - x) ** (beta - 1) / B(alpha, beta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-1.1403464899034472"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "z = a_b_test_statistic(1000, 200, 1000, 180)\n",
    "z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.254141976542236"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "two_sided_p_value(z) # 두 광고 효과의 평균이 같을 확률이 높음 (p > 0.05) = 광고효과가 같음"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-2.948839123097944"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "z = a_b_test_statistic(1000, 200, 1000, 150)\n",
    "z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.003189699706216853"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "two_sided_p_value(z) # 두 광고 효과의 평균이 같을 확률이 매우 낮음 (p < 0.05) = 광고효과가 다름"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "##\n",
    "#\n",
    "# Bayesian Inference\n",
    "#\n",
    "##\n",
    "\n",
    "def B(alpha, beta):\n",
    "    \"\"\"a normalizing constant so that the total probability is 1\"\"\"\n",
    "    return math.gamma(alpha) * math.gamma(beta) / math.gamma(alpha + beta)\n",
    "\n",
    "def beta_pdf(x, alpha, beta):\n",
    "    if x < 0 or x > 1:          # no weight outside of [0, 1]    \n",
    "        return 0        \n",
    "    return x ** (alpha - 1) * (1 - x) ** (beta - 1) / B(alpha, beta)"
   ]
  }
 ],
 "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.4.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}