4_Hypothesis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 数据探索之假设检验\n",
    "\n",
    "\n",
    "### 单个样本均值的假设检验\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import brfss\n",
    "from scipy import stats  # 导入scipy的统计模块\n",
    "\n",
    "%config InlineBackend.figure_format = 'retina'\n",
    "\n",
    "df = brfss.ReadBrfss()  # 读取数据"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "df2 = df[['bmi', 'income']].dropna()\n",
    "bmi = df2.bmi\n",
    "bmi_rich = df2[df2.income == 8].bmi  # 富人bmi数据\n",
    "bmi_ord = df2[df2.income != 8].bmi  # 普通人bmi数据"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "28.18812531332513"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mu = np.mean(bmi)  # 计算bmi均值\n",
    "mu"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "p value is 0.0000000000\n"
     ]
    }
   ],
   "source": [
    "t_stats, p_value = stats.ttest_1samp(bmi_rich, mu)\n",
    "print(\"p value is %.10f\" % (p_value/2))  # 除以2是因为我们使用的是单边检验"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "p value is 0.2906513576\n"
     ]
    }
   ],
   "source": [
    "t_stats, p_value = stats.ttest_1samp(bmi_rich[:500], mu)  # 选择前500个数据进行检验\n",
    "print(\"p value is %.10f\" % (p_value/2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "p value is 0.0000660545\n"
     ]
    }
   ],
   "source": [
    "t_stats, p_value = stats.ttest_1samp(bmi_rich[:4000], mu)  # 选择前4000个数据进行检验\n",
    "print(\"p value is %.10f\" % (p_value/2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "#### 蒙特卡洛模拟之 bootstrap 方法\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def bootstrap_replicate_1d(data, func):   # 进行一次重新抽样，并返回检验统计量\n",
    "    return func(np.random.choice(data, size=len(data)))  \n",
    "\n",
    "def draw_bs_reps(data, func, size=1):\n",
    "    bs_replicates = np.empty(size)  # 初始一个空数组\n",
    "    for i in range(size):   # 进行多次重新抽样\n",
    "        bs_replicates[i] = bootstrap_replicate_1d(data, func)  \n",
    "    return bs_replicates  # 返回多次抽样的检验统计量数组\n",
    "\n",
    "\n",
    "def bootstrap_pvalue_1samp(data, pop_stats, func, size=1):\n",
    "    sample_stats = func(data)  # 计算原有样本的检验统计量    \n",
    "    translated_data = data - sample_stats + pop_stats  # 数据平移\n",
    "    bs_replicates = draw_bs_reps(translated_data, func, size) # 重新抽样\n",
    "    p = np.sum( bs_replicates < sample_stats) / size # 计算抽样统计量小于原有统计量的概率\n",
    "    return p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bootstrap_pvalue_1samp(bmi_rich, mu, np.mean, size=10000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.71260000000000001"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bootstrap_pvalue_1samp(bmi_rich[:500], mu, np.mean, size=10000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0001"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bootstrap_pvalue_1samp(bmi_rich[:4000], mu, np.mean, size=10000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 两个独立样本均值的假设检验\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "p value is 0.0000000000\n"
     ]
    }
   ],
   "source": [
    "t_stats, p_value = stats.ttest_ind(bmi_rich, bmi_ord)\n",
    "print(\"p value is %.10f\" % (p_value/2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "p value is 0.1219871502\n"
     ]
    }
   ],
   "source": [
    "t_stats, p_value = stats.ttest_ind(bmi_rich[:500], bmi_ord[:500])\n",
    "print(\"p value is %.10f\" % (p_value/2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "p value is 0.0000846839\n"
     ]
    }
   ],
   "source": [
    "t_stats, p_value = stats.ttest_ind(bmi_rich[:1000], bmi_ord[:1000])\n",
    "print(\"p value is %.10f\" % (p_value/2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 蒙特卡洛模拟之Permutation方法\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def diff_of_means(data_1, data_2):  \n",
    "    diff = np.mean(data_2) - np.mean(data_1)  # 计算两组数据均值的差\n",
    "    return diff\n",
    "\n",
    "\n",
    "def permutation_sample(data1, data2):  # 产生新的分组数据\n",
    "    data = np.concatenate((data1, data2))  # 合并两组数据\n",
    "    permuted_data = np.random.permutation(data)  # 对合并后的数据进行重新排列\n",
    "    perm_sample_1 = permuted_data[:len(data1)]  # 分成新的两组数据\n",
    "    perm_sample_2 = permuted_data[len(data1):]\n",
    "    return perm_sample_1, perm_sample_2\n",
    "\n",
    "\n",
    "def draw_perm_reps(data_1, data_2, func, size=1):  # 进行多次重新分组的操作\n",
    "    perm_replicates = np.empty(size)\n",
    "    for i in range(size):\n",
    "        perm_sample_1, perm_sample_2 = permutation_sample(data_1, data_2)\n",
    "        perm_replicates[i] = func(perm_sample_1, perm_sample_2)\n",
    "    return perm_replicates\n",
    "\n",
    "\n",
    "def permutation_pvalue(data_1, data_2, func, size=1):  # 计算P值\n",
    "    empirical_test_stats = func(data_1, data_2)\n",
    "    perm_replicates = draw_perm_reps(data_1, data_2, func, size)\n",
    "    p = np.sum(perm_replicates > empirical_test_stats) / len(perm_replicates)\n",
    "    return p\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "permutation_pvalue(bmi_rich, bmi_ord, diff_of_means, size=10000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.1142"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "permutation_pvalue(bmi_rich[:500], bmi_ord[:500], diff_of_means, size=10000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0001"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "permutation_pvalue(bmi_rich[:1000], bmi_ord[:1000], diff_of_means, size=10000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}