01-Data-Science-From-Scratch/04-확률&통계(Probability&Statistics)/2. 가설과 추론(Hypothesis and Inference).ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 가설과 추론\n",
    "\n",
    "통계와 확률 이론을 통해 가설을 검정하고 통계적 추론을 할 수 있다.  \n",
    "가설(hypothesis)이란 '데이터 과학자는 MATLAB보다 R이나 Python을 선호한다.' 등과 같은 주장이며, 데이터 통계에 대한 얘기가 될 수도 있다.\n",
    "통계치들은 다양한 가정에서 특정 분포에 대한 확률변수의 관측치로 이해할 수도 있고, 가정이 얼마나 타당한지 알 수 있게 해주기도 한다.  \n",
    "고전적 가설검정에는 기본적인 가설 귀무가설(H0, null hypothesis)과 비교할 대립가설(H1, alternative hypothesis)로 구성된다.\n",
    "통계를 통해 H0을 기각할지 말지를 정할 수 있다.\n",
    "\n",
    "동전 던지기의 상황에 대해서 동전의 앞 뒷면이 나올 확률이 공평하다면 'p=0.5'가 귀무가설(H0)이 되며, p가 0.5가 아닌 경우가 대립가설이 된다.\n",
    "동전을 n번 던져서 나온 횟수 X를 세는 것으로 검정을 진행할 수 있다. 베르누이 분포를 따르면 X가 이항분포를 따르는 확률변수라는 걸 의미한다. \n",
    "또한 이항분포는 정규분포로 근사할 수 있다.\n",
    "\n",
    "### 단측검정\n",
    "\n",
    "가설을 '남자의 키가 클 것이다.'와 같이 세울 수 있다. 이처럼 가설을 세울 때 방향이 정해지는 것을 단측검정이라고 한다. 이미 가설의 방향을 알고 있거나 잠정적으로 결정한 상태에서 검증에 들어간다.\n",
    "\n",
    "### 양측검정\n",
    "\n",
    "반면에 같은 내용이라고 해도 '남녀 사이에 차이가 있을 것이다.'처럼 어느 쪽으로 방향을 정하지 않고 가설을 세울 때, 양측검정이라고 한다. \n",
    "\n",
    "## 오류\n",
    "\n",
    "### 1종 오류\n",
    "\n",
    "귀무가설H0이 참일 때 H0을 기각하는 false positive(가양성) 오류를 뜻한다.\n",
    "**유의수준(significance)**는 제 1종 오류를 어느정도 허용할 것인가를 의미한다. 유의수준은 보통 5%나 1%로 설정하는 경우가 많다.\n",
    "\n",
    "### 2종 오류\n",
    "\n",
    "2종 오류를 범하지 않을 확률을 구하면 모델의 검정력(power)을 알 수 있다.\n",
    "제 2종 오류란 H0가 거짓이지만 H0를 기각하지 않는 오류를 의미한다.\n",
    "\n",
    "유의 수준에 맞게 가설을 기각하고 오류를 줄이면 검정력이 좋아졌다고 볼 수 있다.\n",
    "\n",
    "### p-value\n",
    "\n",
    "가설을 바라보는 또 다른 관점은 **p-value**이다. 이는 확률 값을 기준으로 구간을 선택하지 않고 H0가 참이라는 전제 하에 실제 관측된 값보다 더 극단적이 나올 확률을 구하는 것이다.\n",
    "\n",
    "Wikipedia - p-value(유의 확률)은 귀무가설이 맞을 경우 대립가설 쪽의 값이 나올 확률을 나타내는 값. 확률 값이라고도 한다. 표본 평균이 귀무가설 값에서 멀수록 작아지게 된다.\n",
    "\n",
    "p-value가 유의수준보다 작을수록 귀무가설을 지지하는 정도가 약하므로 귀무가설을 기각하게 되며(=대립가설 채택), p-value가 클수록 귀무가설을 지지하는 정도가 커짐으로 귀무가설을 채택하게 된다.\n",
    "\n",
    "#### p-value의 모순?\n",
    "\n",
    "\"귀무가설이 'p=0.5'이고 유의수준 5%에서 극한값을 매겨서 p-value를 측정하면 p=0.5+-5%에서 멀어질 때마다 극한값이 \n",
    "카운트되니까 p-value 값이 높을 수록 귀무가설에서 멀어지는 게 아닐까?\"\n",
    "\n",
    "(유의수준 = H0을 기각하는 구간)  \n",
    "\\>\\> 95% 구간의 신뢰구간 바깥에 양 끝 꼬리 2.5%가 유의수준에 해당하는데, p-value가 5%보다 커야 P값이 5%의 유의수준에서 귀무가설을 기각하지 않을 때 신뢰구간에 들어갈 수 있다.\n",
    "\n",
    "\n",
    "### 연속 수정(continuity correction)\n",
    "\n",
    "\n",
    "In probability theory, a continuity correction is an adjustment that is made when a discrete distribution is approximated by a continuous distribution.  \n",
    "확률 이론에서 연속 수정은 이항분포를 연속분포로 근사했을 때 일어나는 수정을 뜻한다.  \n",
    "동전의 앞면이 530번 나올 확률을 확인하는 방법으로 529 - 531를 범위로 잡는 것보다 529.5 - 530.5를 범위로 잡는 것이 더 정확하다.\n",
    "\n",
    "### 신뢰구간\n",
    "\n",
    "사건에 대한 분포를 모를 때 관측된 값에 대한 신뢰구간을 사용해 가설을 검정할 수 있다. p의 정확한 값을 알고 있다면 중심극한정리를 이용해 베르누이 확률변수들의 평균은 대략 평균이 p이고 표준편차는 다음과 같은 정규분포로 추정할 수 있다.\n",
    "\n",
    "```python\n",
    "math.sqrt(p * (1-q) / 1000)\n",
    "```\n",
    "\n",
    "정확한 p의 값을 모른다면 추정치를 사용할 수 있다."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "46\n"
     ]
    }
   ],
   "source": [
    "import math\n",
    "\n",
    "def normal_cdf(x, mu=0,sigma=1):\n",
    "    return (1 + math.erf((x - mu) / math.sqrt(2) / sigma)) / 2\n",
    "\n",
    "def inverse_normal_cdf(p, mu=0, sigma=1, tolerance=0.00001):\n",
    "    \"\"\"find approximate inverse using binary search\"\"\"\n",
    "\n",
    "    # if not standard, compute standard and rescale\n",
    "    if mu != 0 or sigma != 1:\n",
    "        return mu + sigma * inverse_normal_cdf(p, tolerance=tolerance)\n",
    "\n",
    "    low_z, low_p = -10.0, 0            # normal_cdf(-10) is (very close to) 0\n",
    "    hi_z,  hi_p  =  10.0, 1            # normal_cdf(10)  is (very close to) 1\n",
    "    while hi_z - low_z > tolerance:\n",
    "        mid_z = (low_z + hi_z) / 2     # consider the midpoint\n",
    "        mid_p = normal_cdf(mid_z)      # and the cdf's value there\n",
    "        if mid_p < p:\n",
    "            # midpoint is still too low, search above it\n",
    "            low_z, low_p = mid_z, mid_p\n",
    "        elif mid_p > p:\n",
    "            # midpoint is still too high, search below it\n",
    "            hi_z, hi_p = mid_z, mid_p\n",
    "        else:\n",
    "            break\n",
    "\n",
    "\n",
    "#####\n",
    "#\n",
    "# probabilities a normal lies in an interval\n",
    "#\n",
    "######\n",
    "\n",
    "# the normal cdf _is_ the probability the variable is below a threshold\n",
    "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",
    "######\n",
    "#\n",
    "#  normal bounds\n",
    "#\n",
    "######\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\n",
    "\n",
    "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)\n",
    "\n",
    "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\n",
    "\n",
    "upper_p_value = normal_probability_above\n",
    "lower_p_value = normal_probability_below\n",
    "\n",
    "\n",
    "# p-value Hacking\n",
    "#\n",
    "# 동전 1000번 던지기 실험을 1000번 반복해 5%의 유의수준을 두고 가설검정을 한다.\n",
    "import random\n",
    "\n",
    "def run_experiment():\n",
    "    return [random.random() < 0.5 for _ in range(1000)]\n",
    "\n",
    "def reject_fairness(experiment):\n",
    "    num_heads = len([flip for flip in experiment if flip])\n",
    "    return num_heads < 469 or num_heads > 531\n",
    "\n",
    "random.seed(0)\n",
    "experiments = [run_experiment() for _ in range (1000)]\n",
    "num_rejections = len([experiment for experiment in experiments if reject_fairness(experiment)])\n",
    "\n",
    "print(num_rejections)\n",
    "\n",
    "# 놀랍게도 p-value가 0.05 이하로 나오며, p-value의 관점에서 p=0.5의 가설은 기각된다."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A/B 테스트 해보기\n",
    "\n",
    "Google Adsense를 관리한다고 가정해보자. 블로그의 포스트 끝 부분에 등록한 광고 A를 N_A명의 사용자가 보았고 n_A명의 사용자가 클릭했다. 또 블로그의 배너에 등록한 광고 B를 N_B명의 사용자가 보았고 n_B명의 사용자가 클릭했다.  \n",
    "각 사용자가 광고를 클릭할 확률을 pA라고 했을 때 n_A/N_A는 평균이 p_A이고 평균이 p_A이고 표준편차가\n",
    "\n",
    "$$\\sigma_A = \\sqrt {P_A(1 - P_A) / N_A}$$\n",
    "인 정규분포에 근접하다는 것을 알 수 있고, 확률을 p_B라고 했을 때 n_B/N_B는 평균이 p_B이고 평균이 p_B이고 표준편차가\n",
    "\n",
    "$$\\sigma_B = \\sqrt {P_B(1 - P_B) / N_B}$$\n",
    "인 정규분포에 근접하다는 것을 알 수 있다."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.254141976542236"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "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_statistics(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)\n",
    "\n",
    "# 두 광고를 모두 1000명이 봤을 때 A는 200명이 클릭을, B는 180명이 클릭을 했다.\n",
    "z = a_b_test_statistics(1000, 200, 1000, 180)\n",
    "two_sided_p_value(z)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "값이 크기 때문에 두 분포가 다르다고 결론지을 수 없다.  \n",
    "만약 150명이 B광고를 클릭했다면 통계는 다음처럼 나타난다."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.003189699706216853"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "z = a_b_test_statistics(1000, 200, 1000, 150)\n",
    "two_sided_p_value(z)\n",
    "# 두 광고가 동일하게 효과적일 때 50만큼의 차이가 발생활 확률은 0.003에 불과하다."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 베이지안 추론\n",
    "\n",
    "여태까지의 방법은 통계적 검정을 확률적으로 설명하는 것에 (H_0이 사실이라면 이런 통계치가 발생할 확률은 3%이다.) 중점을 두었다. 또 다른 추론의 방법으로 알려지지 않은 파라미터를 변수로 보는 방법이 있다. 분석가에게 파라미터에 대한 사전분포가 주어지고, 관측된 데이터와 베이즈 정리를 사용해 사후분포를 갱신할 수 있다. 통계적 검정에 대한 확률적으로 결론을 내는 대신에 파라미터에 대해 확률적으로 결론을 낼 수 있다.  \n",
    "동전 던지기, 주사위 굴리기 처럼 알려지지 않은 파라미터가 '확률'이라면, 보통 모든 확률 값이 0에서 1 사이로 정의되는 베타분포를 사전분포로 사용한다.\n",
    "\n",
    "이러한 과정을 통하면 가설에 대해 새로운 관점에서 확률적으로 말할 수 있게 된다. 예컨데 \"관측된 데이터와 사전분포로 볼 때 동전의 앞면이 나올 확률이 49% ~ 51%일 경우는 5%밖에 되지 않는다.\"가 그것이다. \"동전이 공평하다면 이런 편향된 데이터를 관측할 가능성은 5%밖에 되지 않는다\"와는 완전히 다른 관점이다."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def B(alpha, beta):\n",
    "    \"\"\"모든 확률 값의 합이 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:\n",
    "        return 0\n",
    "    return x ** (alpha - 1) * (1 - x) ** (beta - 1) / B(alpha, beta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Counter({'False': 4989, 'True': 5011})"
      ]
     },
     "execution_count": 75,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 대수의 법칙 놀이\n",
    "\n",
    "from collections import Counter\n",
    "big_number = []\n",
    "for _ in range(10000):\n",
    "    big_number.append(random.choice([\"True\",\"False\"]))\n",
    "\n",
    "Counter(big_number)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [default]",
   "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}