math/thinkstats/hypothesis-testing.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 假设检验\n",
    "## 经典假设检验\n",
    "- “直观效应”：第一胎和其他胎的区别。\n",
    "- 经典假设检验的目的是回答一个问题：“给定一个样本和一个直观效应，这个效应是偶然出现的概率为多少？”\n",
    "- 回答步骤：\n",
    "    - 选择一个检验统计量（test statistic），对直观效应进行量化。直观效应：第一胎和其他胎的妊娠时间存在差异，因此很自然选择这两个群组的均值作为检验统计量。\n",
    "    - 定义原假设（null hypothesis）。原假设是系统的一个模型，所基于的假设是直观效应不为真（直观效应是偶然产生的）。原假设：第一胎和其他胎没有区别，即两个群组的妊娠时间具有相同的分布。\n",
    "    - 计算p值（p value）。p值是在原假设为真时，直观效应的概率。我们将计算两个群组均值的实际差异，然后计算在原假设为真的情况下，这个差异等于或大于实际值的概率。\n",
    "    - 解释结果。如果p值很低，那么我们称这个效应是统计显著（statistically significant）的，即不太可能是偶然发生的。在这种情况下，我们推断，这个效应在大规模总体中出现的可能性更大。\n",
    "    \n",
    "假设投一枚硬币250次，结果得到140次正面，110次反面。基于这个结果怀疑这个硬币质地不均匀。结果在0.06，如果硬币是质地均匀的，预期6%可能性看到正面和反面的差值达到30的情况。按照惯例，5%是统计显著的阈值，如果p值小于5%，那么我们认为该效应是显著的，否则不是。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from __future__ import print_function, division\n",
    "\n",
    "import numpy as np\n",
    "import pandas\n",
    "import random\n",
    "import math\n",
    "import thinkstats2\n",
    "import thinkplot\n",
    "import nsfg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "coin test p-value 0.059\n"
     ]
    }
   ],
   "source": [
    "class CoinTest(thinkstats2.HypothesisTest):\n",
    "    \"\"\"Tests the hypothesis that a coin is fair.\"\"\"\n",
    "\n",
    "    def TestStatistic(self, data):\n",
    "        \"\"\"Computes the test statistic.\n",
    "\n",
    "        data: data in whatever form is relevant        \n",
    "        \"\"\"\n",
    "        heads, tails = data\n",
    "        test_stat = abs(heads - tails)\n",
    "        return test_stat\n",
    "\n",
    "    def RunModel(self):\n",
    "        \"\"\"Run the model of the null hypothesis.\n",
    "\n",
    "        returns: simulated data\n",
    "        \"\"\"\n",
    "        heads, tails = self.data\n",
    "        n = heads + tails\n",
    "        sample = [random.choice('HT') for _ in range(n)]\n",
    "        hist = thinkstats2.Hist(sample)\n",
    "        data = hist['H'], hist['T']\n",
    "        return data\n",
    "\n",
    "# run the coin test\n",
    "ct = CoinTest((140, 110))\n",
    "pvalue = ct.PValue()\n",
    "print ('coin test p-value', pvalue)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 检验均值差和其他统计量\n",
    "最常检验的效应之一是两组样本的均值差。例如全国家庭增长数据，测试第一胎和其他胎妊娠时间差异。\n",
    "\n",
    "- 置换（permutation）：从两组中取值混排，分为大小为n和m两个组，以此模拟原假设。\n",
    "- 单侧（one-sided）：之前是双侧（two-sided）检验。如果我们有理由认为第一胎可能出生较晚，那么就不应该检验差值，而是使用单侧检验。\n",
    "- 检验标准差的差值\n",
    "\n",
    "三种检验p值都很大，所以这个效应不是统计显著的。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def MakeFrames():\n",
    "    \"\"\"Reads pregnancy data and partitions first babies and others.\n",
    "\n",
    "    returns: DataFrames (all live births, first babies, others)\n",
    "    \"\"\"\n",
    "    preg = nsfg.ReadFemPreg()\n",
    "\n",
    "    live = preg[preg.outcome == 1]\n",
    "    firsts = live[live.birthord == 1]\n",
    "    others = live[live.birthord != 1]\n",
    "\n",
    "    assert len(live) == 9148\n",
    "    assert len(firsts) == 4413\n",
    "    assert len(others) == 4735\n",
    "\n",
    "    return live, firsts, others\n",
    "live, firsts, others = MakeFrames()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "prglngth\n",
      "\n",
      "means permute two-sided\n",
      "p-value = 0.151\n",
      "actual = 0.0780372667775\n",
      "ts max = 0.183327219593\n",
      "\n",
      "means permute one-sided\n",
      "p-value = 0.088\n",
      "actual = 0.0780372667775\n",
      "ts max = 0.166034355154\n",
      "\n",
      "std permute one-sided\n",
      "p-value = 0.1\n",
      "actual = 0.17600895914\n",
      "ts max = 0.339693796308\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEZCAYAAABmTgnDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHNBJREFUeJzt3XuYXFWd7vHvG0gQSJDbCMr1gUZwvDAiIjAqrQ4QECfK\nGQeBkcsAD6MiKOLgjHjoOepBxBH1gJeMkTEqknGQCSJoeEhaBOQSMEElARoQQggwDCoCgjH8zh97\ndWd3paqrqqt212W/n+fph127Vu1atan022utvdZWRGBmZjaRaZ2ugJmZdT+HhZmZ1eWwMDOzuhwW\nZmZWl8PCzMzqcliYmVldDguzAkn6vaRdO10Ps1Y5LKyjJP1a0rOSnpK0RtIlkjbrdL0AJL0gabcm\nyi+R9Pf5fRExKyJ+XUDdHpD01jYc53hJP21Hnay/OSys0wJ4e0RsAewD7Auc0+xBJG3U7oqR1a3f\niXJ8TmuRw8K6gQAiYg1wDfAqAElbSPq6pEckrZL0SUlKzx0v6QZJn5f0BHBuxb7fSLpX0oGSTpD0\nkKRHJR039qYVLYH8X9mSfpLqdWdq9bxb0paSfiDpcUn/k7Zflsp/CngTcFEq/6W0f6x1kj7P/PT6\nByR9vPK9JV0g6UlJ90maXfVkSfOBnYEfpPc6K+3fX9KN6bP/XNJBudeckI75VPrv0ZL2Ar4CHJC6\ny55s6f+i9TWHhXUNSTsBhwN3pF3fBP4I7Aa8FjgYODn3kjcAI8BLgE+nffsBy4CtgcuABcDrgN2B\n95L9Mp+omysAImL0F+2rI2KLiPge2b+XbwA7kf2yfha4OJU/B/gpcFoqf3r+eMlFwCxgV2AQOE7S\nibnn9wNWANsAFwDzqlYw4jjgIeCI9F6fS6F1FfB/ImIr4CzgcknbpM/7ReDQ1II7EFgWESuBfwB+\nlrrLtp7gvFjJOSysG/xX+qv2emAJcJ6klwCHAR+OiOci4gngC8DRudetjogvR8QLEfF82vdARMyP\nbNGzBcDLgH+JiLURcS1Z+Aw0UTeNbkTEkxFxRUQ8HxHPAOcBb27k9ZKmAUcBH4uIZyPiQeBfyQJs\n1IMR8Y1U928C26fzULduwN8BP4yIH6e6XgcsJQtfgHXAqyW9KCIei4gVdeptNs7Gna6AGTAnIpbk\nd0jaBZgOrBnteUo/D+WKrapyrMdy238ASEGT3zdzMpWUtClZYB0KbJnqM1OSov6KnNuS/XvL1/9B\nYIfc40dHNyLiD6nLbSbweAPV2wX4W0nvGK1uer/FEfGspKOAjwLfkHQDcFZE3N3Acc0AtyysO6jK\nvlXAc8A2EbF1RGwVEVtGxGtyZVodmH0GyHdJbV+n/EeAPYDXR8SWrG9VjNZ/ovo8Aawl+6U+ahdg\ndcO1Ha/yvVYB89O5Gj1fsyLiswARcW1EHEL2Ge8G5jZQZ7MxDgvrShHxKLAIuFDSLGV2k1Sv26dS\ntSAatQw4UtKmkgaAkyqef5RsvGTULLKWyVOStgaGKso/VlF+TES8APwH8GlJM1PL6cPAtxr9IHXq\n9m3gHZIOkTRN0oskHSTpZZJeIumv09jFWuBp4IVcnXeUNH2S9bCScFhYp030l+1xwAzgLuBJ4HvU\n/+u/3vHzjy8k++X5KHAJ2S/cvCFgfro66W9S+c3IWgk3AVdXlP8i8O50pdQXqrzf6WSD4veTjc98\nOyIuaaLueZ8BPpHqdmZEPAzMAf4Z+G+yLq6zyP6NTwPOJGvFPEHWInpfOs5i4FfAo5Ia6e6yklKR\nNz+SNA84AnisovsgX+ZLZAOZzwAnRMSywipkZmaTUnTL4hKywcCqJB0G7B4RewCnAl8tuD5mZjYJ\nhYZFRNwA/GaCInOA+ansLcCLJW1XZJ3MzKx5nR6z2IHxlz+uZvylhGZm1gU6HRZmZtYDOj0pbzXZ\n0gmjdqTGdeeSfD24mdkkRMREl5A3ZCrCYnTmbTVXAh8AFkjaH/htRDxWoyxFXrlVNkNDQwwNDXW6\nGk0bGRkZ2x4YaGbVjmL16vmcyMLFy1lwzVKe/+PaKX3fu26+ij/f/4gpfc9OuPyL/zAl75NWQGhZ\noWEh6VKyBdO2kfQQcC7ZdfMREXMj4mpJh0saIbt09sTaRzOzInQqFBq1yYzpHHXYvsx5696drkqp\nFRoWEXFMA2VOK7IOZrZeO4OhqF/iQ0OPMjQ0NX91W+M6PWZhHTI4ONjpKvSVTp/PoloHnfirvtPn\n0qordAZ3OzW2sKf1u24ds+iEXmglWOdJ6pkBbjNrkYPBOs1hYdalmg0Ih4AVyWFh1kUcENatHBZm\nHdRMODgYrJMcFmZTZDLjDg4I6xYOC7MpsHDxcuYv/FlDZR0Q1o0cFmYFcPeS9RuHhVkbNBMOx805\nwMFgPcdhYTZJvnLJysRhYdYEdy9ZWTkszCbgcDDLOCzMqmg0JBwQVhYOC7MK9S5zdUBYGTkszJJa\nrQmHg5nDwkquXneTL3M1yzgsrHQaGY9wa8JsPIeFlYZDwmzyHBbW9+qFhAPCrD6HhfW1Wlc2OSDM\nmuOwsL4zUUvCIWE2OQ4L6ysTzZHwlU1mk+ewsL7hLiez4jgsrOfV6nZyS8KsfRwW1tNqtSYcFGbt\npYjodB0aIil6pa5WnJGREQCuu2WEnyxb7aU5zOqQRESo1eO4ZWE957pbRrhiya/YfPPNx+13a8Ks\nOA4L6xkLFy9n3veG+ePaP43b79aEWfEcFtYzFlyzdIOgcGvCbGo4LKzrVbvaacb0jR0UZlPIYWFd\nq9YlsTOmb8znP/J2BgYGOlQzs/KZ1ukKmNVSKyje/sY9O1Qjs/Jyy8K6TrUWxegg9it33nyCV5pZ\nURwW1lWqTbLbZMZ0Lr3gJGD9PAszm1ruhrKuUSsojjps3w7VyMxGFd6ykDQb+AJZMM2LiPMrnt8G\n+DbwUmAj4F8j4t+Lrpd1l2pB4audzLpHoWEhaRpwEfA24BHgNkkLI2JlrthpwLKIOEzStsDdkr4d\nEX+qckjrM14E0Kw3FN0NtR9wb0Q8GBFrgcuAORVlHgVmpe1ZwP84KMphtDXhoDDrfkV3Q+0ArMo9\nfpgsQPL+DbhO0iPATOCogutkXWCi8QkHhVn36Yarof4JWB4Rb5G0O3CtpNdExNOVBYeGhsa2BwcH\nGRwcnLJKWvt4fMKsOMPDwwwPD7f9uIUuUS5pf2AoImanxx8DIj/ILelq4NMRcWN6fB1wdkQsrTiW\nlyjvce0Yn8hfOusZ3Gb1tWuJ8qLHLG4DBiTtImkG8B7gyooyK4C/ApC0HfBy4P6C62VTzOMTZr2t\n0G6oiFgn6TRgEesvnV0h6dTs6ZgLnAdcImk5IOAfI+LJIutlU8vjE2a9z3fKs0K1e3zC3VBmzfGd\n8qyref6EWX9xWFjbVWtNgIPCrJc5LKytPD5h1p8cFtZWC64Zd8WzWxNmfcKrzlrbLFy8fNwYhYPC\nrH84LKwtKrufNpkx3UFh1kccFtayauMUvgeFWX9xWFjLPE5h1v8cFtYSj1OYlYPDwibN4xRm5eGw\nsEmr7H7yOIVZ//I8C2tataU83P1k1t8cFtaUWjO0HRRm/c1hYQ2ptTDg6FIeZtbfHBZWlxcGNDOH\nhU3ICwOaGTgsbALtvnGRmfUuXzprVTkozCzPYWEbcFCYWSWHhY3joDCzahwWNo4XBTSzahwWNsaL\nAppZLQ4LG5NvVXhWtpnlOSwM2LBV4VnZZpbnsDAvNW5mdXlSXonVWu/JrQozq+SWRYlVCwoPaptZ\nNW5ZlFTlGIXXezKziTgsSqZa19MmM6Zz6QUndbBWZtbtHBYlUmupcY9RmFk9DosSqZyd7a4nM2uU\nw6IkPDvbzFrhq6FKwPMozKxVDosSqOx+8hiFmTXL3VB9rNqVT+5+MrPJKLxlIWm2pJWS7pF0do0y\ng5J+LumXkpYUXaeyqHaJrIPCzCaj0JaFpGnARcDbgEeA2yQtjIiVuTIvBi4GDomI1ZK2LbJOZVFr\n0p2Z2WQU3Q21H3BvRDwIIOkyYA6wMlfmGODyiFgNEBFPFFynUqhcbtyT7sysFUV3Q+0ArMo9fjjt\ny3s5sLWkJZJuk/TeguvU97zcuJm1WzcMcG8M7AO8Fdgc+Jmkn0XESGer1Zt8mayZFaHosFgN7Jx7\nvGPal/cw8EREPAc8J+l6YG9gg7AYGhoa2x4cHGRwcLDN1e1t1ZbzcKvCrFyGh4cZHh5u+3EVEW0/\n6NjBpY2Au8kGuNcAtwJHR8SKXJm9gP8HzAY2AW4BjoqIuyqOFUXWtR8c89F5fX+Z7MjI+r8hBgYG\nOlgTs94giYhQq8cptGUREesknQYsIhsfmRcRKySdmj0dcyNipaQfA3cC64C5lUFhjen3oDCzzil8\nzCIifgTsWbHvaxWPPwd8rui69LOFi5ePe+ygMLN28nIffaLyUlkzs3ZyWPQBXyprZkVzWPQ4Xypr\nZlPBYdHjvKKsmU0Fh0UP8w2NzGyqOCx6WOWgtoPCzIrisOhRHtQ2s6nUDWtDWROq3dDIrQozK5pb\nFj2mMijArQozK55bFj2k1g2N3Kows6I5LHpEtfkUvqGRmU2VCbuhJP17bvv4wmtjNXk+hZl1Ur0x\ni3z/xhlFVsRq83wKM+u0emHhG0h0mJfzMLNuUG/MYkdJXwKU2x4TEacXVjMD3P1kZt2hXlh8NLe9\ntGYpK4S7n8ysW0wYFhHxzamqiG3Iy3mYWbeoOylP0vGS7pD0TPpZKum4qahcmXk5DzPrJhO2LNLl\nsh8CzgTuIBu72Ae4QFJExLeKr2I5uVVhZt2k3pjF+4B3RcSvc/sWS/pfwGWAw6LNqq395FaFmXVa\nvW6oLSqCAoC0b4siKlR2XiTQzLpRvbD4wySfs0motfaTmVmn1euGeoWkO6vsF7BbAfUpLa/9ZGbd\nrF5Y7A1sB6yq2L8T8GghNSopT74zs25WrxvqQuB3EfFg/gf4XXrO2sCT78ys29ULi+0i4heVO9O+\nXQupUQn5Mlkz63b1wmLLCZ7btJ0VKStPvjOzXlAvLJZKOqVyp6STgduLqVJ5eEVZM+sV9Qa4PwRc\nIelY1ofDvsAM4F1FVqwMPKhtZr2i3kKCjwEHSnoL8Kq0+4cRsbjwmvU5D2qbWS9p6B7cEbEEWFJw\nXUrD3U9m1mvqrjpr7efuJzPrNQ6LDnD3k5n1GofFFFu4ePm4xw4KM+sFDospVjkBz8ysFzgsppAn\n4JlZryo8LCTNlrRS0j2Szp6g3OslrZV0ZNF16hQv62FmvarQsJA0DbgIOBR4JXC0pL1qlPsM8OMi\n69NJblWYWS8rumWxH3BvWq12LdmtWOdUKfdB4D+BxwuuT8e4VWFmvazosNiB8ffCeDjtGyPpZcA7\nI+IrZDdV6jtuVZhZr+uGAe4vAPmxjL4KDM/WNrN+0NByHy1YDeyce7xj2pe3L3CZJAHbAodJWhsR\nV1YebGhoaGx7cHCQwcHBdte37Txb28ym0vDwMMPDw20/riKi7QcdO7i0EXA38DZgDXArcHRErKhR\n/hLgBxHx/SrPRZF1bbeFi5ez4Jqlnq3dZiMjI2PbAwMDHayJWW+QRES03GNTaMsiItZJOg1YRNbl\nNS8iVkg6NXs65la+pMj6TKXKoHD3k5n1sqK7oYiIHwF7Vuz7Wo2yf190faZC5YD2JjOmu/vJzHpa\n4WFRRpWXyV56wUkdrI2ZWeu64WqovuLLZM2sHzks2siXyZpZv3JYtJEvkzWzfuWwaBPfU9vM+pnD\nog3c/WRm/c5h0QbufjKzfuewaJG7n8ysDBwWLfLS42ZWBg6LFnhOhZmVhcOiBW5VmFlZOCwmya0K\nMysTh8UkuVVhZmXisJgEtyrMrGwcFpPgVoWZlY3DokluVZhZGTksmuRWhZmVkcOiCW5VmFlZOSwa\n5MUCzazMHBYN8mKBZlZmDosGeLFAMys7h0UDPKhtZmXnsKjDg9pmZg6LCXlQ28ws47CYgAe1zcwy\nDosaPKhtZraew6IGD2qbma3nsKjCg9pmZuM5LKpwq8LMbDyHRQW3KszMNuSwyPGlsmZm1Tkscnyp\nrJlZdQ6LHF8qa2ZWncMiWbh4+bjHDgozs/UcFknlFVBmZraewyLxFVBmZrUVHhaSZktaKekeSWdX\nef4YScvTzw2SXl10nSq5C8rMbGKFhoWkacBFwKHAK4GjJe1VUex+4M0RsTfwKeDfiqxTNe6CMjOb\nWNEti/2AeyPiwYhYC1wGzMkXiIibI+J36eHNwA4F12kD7oIyM5tY0WGxA7Aq9/hhJg6Dk4FrCq1R\nBXdBmZnVt3GnKzBK0luAE4E31iozNDQ0tj04OMjg4GBL71ltxraZWS8bHh5meHi47cdVRLT9oGMH\nl/YHhiJidnr8MSAi4vyKcq8BLgdmR8R9NY4V7a7rMR+d54l4PWZkZGRse2BgoIM1MesNkogItXqc\noruhbgMGJO0iaQbwHuDKfAFJO5MFxXtrBUVRHBRmZo0ptBsqItZJOg1YRBZM8yJihaRTs6djLvAJ\nYGvgy5IErI2I/YqsF3iswsysGYWPWUTEj4A9K/Z9Lbd9CnBK0fWo5MtlzcwaV8oZ3L5nhZlZc0oX\nFr5nhZlZ80oXFr5nhZlZ80oXFr4CysyseaULizwHhZlZY0oVFpWXy5qZWWNKFRa+XNbMbHJKFRa+\nXNbMbHJKExaesW1mNnmlCQt3QZmZTV4pwsIzts3MWtP3YeEZ22Zmrev7sPCMbTOz1vV1WFR2P3nG\ntpnZ5PRtWLj7ycysffoyLCqDAtz9ZGbWir4Mi8pxCnc/mZm1pu/CwuMUZmbt13dhUTn5zkFhZta6\nvgoLT74zMytGX4WFWxVmZsXoq7Bwq8LMrBh9ExZeVdbMrDh9ERbVJuCZmVn79EVYeP0nM7Ni9XxY\neF6FmVnxej4sfAWUmVnxejosPK/CzGxq9HRYuFVhZjY1ejos3KowM5saPRsWnldhZjZ1ejYsKrug\nzMysOD0ZFh7YNjObWj0ZFh7YNjObWj0XFm5VmJlNvcLDQtJsSSsl3SPp7BplviTpXknLJP1FrWNV\nWwPKrQozs+IVGhaSpgEXAYcCrwSOlrRXRZnDgN0jYg/gVOCrtY7nNaDaZ3h4uNNV6Cs+n+3jc9md\nim5Z7AfcGxEPRsRa4DJgTkWZOcB8gIi4BXixpO2qHcxrQLWP/0G2l89n+/hcdqeiw2IHYFXu8cNp\n30RlVlcpM467n8zMplbPDXCDu5/MzKaaIqK4g0v7A0MRMTs9/hgQEXF+rsxXgSURsSA9XgkcFBGP\nVRyruIqamfWxiFCrx9i4HRWZwG3AgKRdgDXAe4CjK8pcCXwAWJDC5beVQQHt+bBmZjY5hYZFRKyT\ndBqwiKzLa15ErJB0avZ0zI2IqyUdLmkEeAY4scg6mZlZ8wrthjIzs/7QFQPcrUzca+S1ZTKJc/na\n3P5fS1ou6eeSbp26WneveudT0p6SbpL0nKQzm3ltGbV4Pv39zGngXB6TztdySTdIek2jr60qIjr6\nQxZYI8AuwHRgGbBXRZnDgB+m7TcANzf62jL9tHIu0+P7ga06/Tm65afB87kt8Drgk8CZzby2bD+t\nnM/0nL+fzZ3L/YEXp+3Zrf7e7IaWRSsT9xp5bZm0OglSdElrs0vUPZ8R8URE3A78qdnXllAr5xP8\n/cxr5FzeHBG/Sw9vZv38tUl9N7vhxE9m4t5omUZeWyatToIM4FpJt0k6pbBa9o5Wvl/+bm6o1XPi\n7+d6zZ7Lk4FrJvlaoPhLZ4viy2iL8ZcRsUbSn5H9o1wRETd0ulJmib+fkyDpLWRXmb6xleN0Q8ti\nNbBz7vGOaV9lmZ2qlGnktWXSyrkkItak//43cAVZc7XMWvl++bu5oZbOib+f4zR0LtOg9lzgryPi\nN828tlI3hMXYxD1JM8gm7l1ZUeZK4DgYmxU+OnGvkdeWyaTPpaTNJM1M+zcHDgF+OXVV70rNfr/y\nLV5/Nzc06fPp7+cG6p5LSTsDlwPvjYj7mnltNR3vhooWJu7Vem2HPkrHtXIuge2AK9KyKhsD34mI\nRZ34HN2ikfOZLg5YCswCXpB0BvDnEfG0v5vjtXI+gT/D388xjZxL4BPA1sCXJQlYGxH7Tfb3pifl\nmZlZXd3QDWVmZl3OYWFmZnU5LMzMrC6HhZmZ1eWwMDOzuhwWZmZWl8PCJk3SuaPLSEv6F0lvTdtv\nlPRLSXdI2kTSBZJ+Ien8iY/YOZJ6ftkISTMk/SRdU1/E8R+QtHWDZf9D0q5F1MM6o+OT8qw/RMS5\nuYfHAv83Ii4FSIu+bRUNTuqRtFFErCugmjVFREvr5nSJY4GrGj3Pk9DMcecCZwKnF1QXm2JuWVhT\nJH1c0t2Srgf2zO2/RNKRkk4C/hb4pKRvSVoIzARul/RuSdtK+k9Jt6SfA9Lrz5U0P/2FP1/SNEmf\nTWWWja4yKukgSUskfU/SCknfytXh9ZJuTOVvlrR5reNU+Vy/r3f8ivJLJH0+rYB6V3rv76dz88lc\nuWPTe98h6Sujf/VL+rKkW1OL69xc+QckDUm6XdlNa16e9r9Z2U1/7kjPbV6lWscAC1P5iyQdkbav\nkPT1tH3iaP0mqNvBym5AtFTSAkmbjVYvPb+ppKslnaRsGY6rUt3ulPTuVHYYOLzaubMe1embePin\nd36AfYDlwCZkyzHcS7pBDXAJcGTldnr8VG77O8CBaXsn4K60fS7ZmjUz0uNTgH9O2zPSc7sABwG/\nAV5K9svrJuBAspu43Afsk14zE9io1nGqfLan0n+rHr9K+SXAeWn7dOAR4CXpPVYBWwF7ka25s1Eq\ndzHwd2l7y/TfaelYr0qPHwDen7bfB8xN21cCB6TtzYBpFfWZBqzJPT4KOD9t3wLclLa/ARxcq27A\nNsBPgE3T/n8Ezknb96f/B9cCx6Z9RwJfy73vrNz2MCW/4VM//bgbyprxJuCKiHgeeF5Sowvj5fvQ\n/wp4Ra5ffWbuL9crI+KPafsQ4NW5v1S3APYA1gK3RlqBVNIyYFfgKeCRiLgDICKeTs/XOs6DE9S3\n2vFvqlJu9PP/AvhFRDyeXnMfWRC+iSxgb0uf90XAY+k170mtnI2B7cnWPxpdGO+K9N/bgXel7RuB\nCyV9B/h+RFSuErptOgejfgp8SNIrgLuALSVtDxwAfBA4oUbd9k91uTHtn5777AL+C/hsRHw399k/\nJ+k8sjsw5sd+1qRzt7LKubMe47CwqZDv6xbwhsju0LV+Z5Ydz1SU+2BEXFtR7iDg+dyudaz/Hlcb\n2K16nDpqHb9WuRcqXjO62J2Ab0bEx8dVKBv4/Qjwuoh4StIlZL+sK4879t4Rcb6kq4C3k/0iPyQi\n7qmoz9jnj4hHJG0JHErWUtiarHvw9xHxTAqCanU7AlgUEcfW+Mw3kt2i87vpfe6VtA9Zl9OnJF0X\nEaPdcKK5cQ7rYh6zsGZcD7xT2RVOs4B3NPi6/C/xRcAZY09Ie9d4zY+B90vaOJXbI9cCqeZuYHtJ\nr0vlZ0raqMZxNq1Tx3a5DvgbZTfrQdJWypaN3gJ4Gvi9slVWD6t3IEm7RcSvIuKzZF1pe1UUeYKs\n6y3vZuDDZP/fbgDOImtxTFS3m4G/lLR72r+ZpD1yx/zfwG8lXZyefynwh8guZrgAeG2u7EuZuAVn\nPcQtC2tYRPxc0gLgTrIui1vzT9fYrnx8BnCxpOVkYwrXA++v8nZfJ+vCuCP9Ffw48M5q1Up1Wyvp\nKOCiFAbPknV5NXWcNuzP12mFpHOARZKmAX8EPhARt6burRVk4xs3VL62ig8pu+PZOuBXrL9FJum9\nXlB2ufLLcy2OnwIHR8T9kh4iG0e5voG6nQB8V9ImqT7nkI1PjX6uMyTNk/QZYDFwgaQX0jHeB5DC\neYeIcBdUn/AS5WZ9QtLxwPYR0fH5LJIOBo6IiDPqFrae4G4os/7xXeDw3MUDnXQycGGnK2Ht45aF\nmZnV5ZaFmZnV5bAwM7O6HBZmZlaXw8LMzOpyWJiZWV0OCzMzq+v/AwGMACtgETbRAAAAAElFTkSu\nQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11b407b50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "class DiffMeansPermute(thinkstats2.HypothesisTest):\n",
    "    \"\"\"Tests a difference in means by permutation.\"\"\"\n",
    "\n",
    "    def TestStatistic(self, data):\n",
    "        \"\"\"Computes the test statistic.\n",
    "\n",
    "        data: data in whatever form is relevant        \n",
    "        \"\"\"\n",
    "        group1, group2 = data\n",
    "        test_stat = abs(group1.mean() - group2.mean())\n",
    "        return test_stat\n",
    "\n",
    "    def MakeModel(self):\n",
    "        \"\"\"Build a model of the null hypothesis.\n",
    "        \"\"\"\n",
    "        group1, group2 = self.data\n",
    "        self.n, self.m = len(group1), len(group2)\n",
    "        self.pool = np.hstack((group1, group2))\n",
    "\n",
    "    def RunModel(self):\n",
    "        \"\"\"Run the model of the null hypothesis.\n",
    "\n",
    "        returns: simulated data\n",
    "        \"\"\"\n",
    "        np.random.shuffle(self.pool)\n",
    "        data = self.pool[:self.n], self.pool[self.n:]\n",
    "        return data\n",
    "\n",
    "\n",
    "class DiffMeansOneSided(DiffMeansPermute):\n",
    "    \"\"\"Tests a one-sided difference in means by permutation.\"\"\"\n",
    "\n",
    "    def TestStatistic(self, data):\n",
    "        \"\"\"Computes the test statistic.\n",
    "\n",
    "        data: data in whatever form is relevant        \n",
    "        \"\"\"\n",
    "        group1, group2 = data\n",
    "        test_stat = group1.mean() - group2.mean()\n",
    "        return test_stat\n",
    "\n",
    "\n",
    "class DiffStdPermute(DiffMeansPermute):\n",
    "    \"\"\"Tests a one-sided difference in standard deviation by permutation.\"\"\"\n",
    "\n",
    "    def TestStatistic(self, data):\n",
    "        \"\"\"Computes the test statistic.\n",
    "\n",
    "        data: data in whatever form is relevant        \n",
    "        \"\"\"\n",
    "        group1, group2 = data\n",
    "        test_stat = group1.std() - group2.std()\n",
    "        return test_stat\n",
    "\n",
    "    \n",
    "def PrintTest(p_value, ht):\n",
    "    \"\"\"Prints results from a hypothesis test.\n",
    "\n",
    "    p_value: float\n",
    "    ht: HypothesisTest\n",
    "    \"\"\"\n",
    "    print('p-value =', p_value)\n",
    "    print('actual =', ht.actual)\n",
    "    print('ts max =', ht.MaxTestStat())\n",
    "\n",
    "\n",
    "def RunTests(data, iters=1000):\n",
    "    \"\"\"Runs several tests on the given data.\n",
    "\n",
    "    data: pair of sequences\n",
    "    iters: number of iterations to run\n",
    "    \"\"\"\n",
    "\n",
    "    # test the difference in means\n",
    "    ht = DiffMeansPermute(data)\n",
    "    p_value = ht.PValue(iters=iters)\n",
    "    print('\\nmeans permute two-sided')\n",
    "    PrintTest(p_value, ht)\n",
    "\n",
    "    ht.PlotCdf()\n",
    "    thinkplot.Config(root='hypothesis1',\n",
    "                   title='Permutation test',\n",
    "                   xlabel='difference in means (weeks)',\n",
    "                   ylabel='CDF',\n",
    "                   legend=False) \n",
    "    \n",
    "    # test the difference in means one-sided\n",
    "    ht = DiffMeansOneSided(data)\n",
    "    p_value = ht.PValue(iters=iters)\n",
    "    print('\\nmeans permute one-sided')\n",
    "    PrintTest(p_value, ht)\n",
    "\n",
    "    # test the difference in std\n",
    "    ht = DiffStdPermute(data)\n",
    "    p_value = ht.PValue(iters=iters)\n",
    "    print('\\nstd permute one-sided')\n",
    "    PrintTest(p_value, ht)\n",
    "\n",
    "# compare pregnancy lengths\n",
    "print('\\nprglngth')\n",
    "data = firsts.prglngth.values, others.prglngth.values\n",
    "RunTests(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "第一组数据（置换检验）：CDF与观察到的差值（0.07）在0.82（即p值0.18的补）处相交。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "birth weight\n",
      "means permute two-sided\n",
      "p-value = 0.0\n",
      "actual = 0.124761184535\n",
      "ts max = 0.0899775334883\n"
     ]
    }
   ],
   "source": [
    "# compare birth weights\n",
    "print('\\nbirth weight')\n",
    "data = (firsts.totalwgt_lb.dropna().values,\n",
    "        others.totalwgt_lb.dropna().values)\n",
    "ht = DiffMeansPermute(data)\n",
    "p_value = ht.PValue(iters=1000)\n",
    "print('means permute two-sided')\n",
    "PrintTest(p_value, ht)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "使用新生儿体重进行分析，得到p值为0。运行1000次模拟，均值差都没有达到观察到的差值0.12磅。因此，我们将报告p<0.001，得出结论为，新生儿体重的差值是统计显著的。\n",
    "\n",
    "## 检验相关性\n",
    "- 检验的原假设是母亲年龄和新生儿体重之间没有相关性。我们可以将观察值混排进行模拟。\n",
    "- TestStatistic计算Pearson相关性的绝对值。RunModel将xs进行混排，返回模拟数据。\n",
    "- 实际数据的相关性为0.07。检验计算得到的p值为0。在1000次重复试验中，模拟得到的最大相关性0.04。因此，虽然观察到的变量相关性很小，但这种相关性是统计显著的。\n",
    "- “统计显著”不一定说明一个效应是重要的，或者在实践中是显著的。只说明一个效应不太可能是偶然产生的。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "age weight correlation 2\n",
      "p-value = 0.0\n",
      "actual = 0.0688339703541\n",
      "ts max = 0.031025636795\n"
     ]
    }
   ],
   "source": [
    "class CorrelationPermute(thinkstats2.HypothesisTest):\n",
    "    \"\"\"Tests correlations by permutation.\"\"\"\n",
    "\n",
    "    def TestStatistic(self, data):\n",
    "        \"\"\"Computes the test statistic.\n",
    "\n",
    "        data: tuple of xs and ys\n",
    "        \"\"\"\n",
    "        xs, ys = data\n",
    "        test_stat = abs(thinkstats2.Corr(xs, ys))\n",
    "        return test_stat\n",
    "\n",
    "    def RunModel(self):\n",
    "        \"\"\"Run the model of the null hypothesis.\n",
    "\n",
    "        returns: simulated data\n",
    "        \"\"\"\n",
    "        xs, ys = self.data\n",
    "        xs = np.random.permutation(xs)\n",
    "        return xs, ys\n",
    "# test correlation\n",
    "live2 = live.dropna(subset=['agepreg', 'totalwgt_lb'])\n",
    "data = live2.agepreg.values, live2.totalwgt_lb.values\n",
    "ht = CorrelationPermute(data)\n",
    "p_value = ht.PValue()\n",
    "print('\\nage weight correlation 2')\n",
    "PrintTest(p_value, ht)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 检验比例与卡方检验\n",
    "检验骰子是否有问题，将骰子掷了60次，结果如下：\n",
    "\n",
    "| 点数 | 频数 |\n",
    "| --- | --- |\n",
    "| 1 | 8 |\n",
    "| 2 | 9 |\n",
    "| 3 | 19 |\n",
    "| 4 | 5 |\n",
    "| 5 | 8 |\n",
    "| 6 | 11 |\n",
    "\n",
    "### 偏差总和统计量\n",
    "每个点数平均出现10次，在这个数据集中，3出现的次数较多，4较少。这些差异统计显著嘛？在这个例子中，我们计算出每个值的预期频数、预期频数与观察频数的差值，以及差值绝对值的和。检验统计量是差值绝对值的和。计算的到p值为0.13，我们预期检验统计量达到或超过观察值的概率为13%，因此，这个直观效应不是统计显著的。\n",
    "\n",
    "### 卡方统计量\n",
    "$x^2=\\sum_i\\frac{(O_i-E_i)^2}{E_i}$\n",
    "\n",
    "其中$O_i$是观察频数，$E_i$是预期频数。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dice test 0.1279\n",
      "dice chi test 0.0438\n"
     ]
    }
   ],
   "source": [
    "class DiceTest(thinkstats2.HypothesisTest):\n",
    "    \"\"\"Tests whether a six-sided die is fair.\"\"\"\n",
    "\n",
    "    def TestStatistic(self, data):\n",
    "        \"\"\"Computes the test statistic.\n",
    "\n",
    "        data: list of frequencies\n",
    "        \"\"\"\n",
    "        observed = data\n",
    "        n = sum(observed)\n",
    "        expected = np.ones(6) * n / 6\n",
    "        test_stat = sum(abs(observed - expected))\n",
    "        return test_stat\n",
    "\n",
    "    def RunModel(self):\n",
    "        \"\"\"Run the model of the null hypothesis.\n",
    "\n",
    "        returns: simulated data\n",
    "        \"\"\"\n",
    "        n = sum(self.data)\n",
    "        values = [1,2,3,4,5,6]\n",
    "        rolls = np.random.choice(values, n, replace=True)\n",
    "        hist = thinkstats2.Hist(rolls)\n",
    "        freqs = hist.Freqs(values)\n",
    "        return freqs\n",
    "\n",
    "class DiceChiTest(DiceTest):\n",
    "    \"\"\"Tests a six-sided die using a chi-squared statistic.\"\"\"\n",
    "\n",
    "    def TestStatistic(self, data):\n",
    "        \"\"\"Computes the test statistic.\n",
    "\n",
    "        data: list of frequencies\n",
    "        \"\"\"\n",
    "        observed = data\n",
    "        n = sum(observed)\n",
    "        expected = np.ones(6) * n / 6\n",
    "        test_stat = sum((observed - expected)**2 / expected)\n",
    "        return test_stat\n",
    "\n",
    "def RunDiceTest():\n",
    "    \"\"\"Tests whether a die is fair.\n",
    "    \"\"\"\n",
    "    data = [8, 9, 19, 5, 8, 11]\n",
    "    dt = DiceTest(data)\n",
    "    print('dice test', dt.PValue(iters=10000))\n",
    "    dt = DiceChiTest(data)\n",
    "    print('dice chi test', dt.PValue(iters=10000))\n",
    "# run the dice test\n",
    "RunDiceTest()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "使用卡方统计量计算的p值为0.04，明显小于使用偏差和的值0.13，如果我们坚持使用5%的阈值，就会认为受检测的效应是统计显著的。但是，将这两个检验放在一起考虑，结果并不明确。\n",
    "\n",
    "这个例子告诉我们：p值取决于检验统计量的选择和原假设模型，有时这些因素决定了一个效应是否统计显著。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "pregnancy length chi-squared\n",
      "p-value = 0.0\n",
      "actual = 101.501414829\n",
      "ts max = 26.3476744126\n"
     ]
    }
   ],
   "source": [
    "class PregLengthTest(thinkstats2.HypothesisTest):\n",
    "    \"\"\"Tests difference in pregnancy length using a chi-squared statistic.\"\"\"\n",
    "\n",
    "    def TestStatistic(self, data):\n",
    "        \"\"\"Computes the test statistic.\n",
    "\n",
    "        data: pair of lists of pregnancy lengths\n",
    "        \"\"\"\n",
    "        firsts, others = data\n",
    "        stat = self.ChiSquared(firsts) + self.ChiSquared(others)\n",
    "        return stat\n",
    "\n",
    "    def ChiSquared(self, lengths):\n",
    "        \"\"\"Computes the chi-squared statistic.\n",
    "        \n",
    "        lengths: sequence of lengths\n",
    "\n",
    "        returns: float\n",
    "        \"\"\"\n",
    "        hist = thinkstats2.Hist(lengths)\n",
    "        observed = np.array(hist.Freqs(self.values))\n",
    "        expected = self.expected_probs * len(lengths)\n",
    "        stat = sum((observed - expected)**2 / expected)\n",
    "        return stat\n",
    "\n",
    "    def MakeModel(self):\n",
    "        \"\"\"Build a model of the null hypothesis.\n",
    "        \"\"\"\n",
    "        firsts, others = self.data\n",
    "        self.n = len(firsts)\n",
    "        self.pool = np.hstack((firsts, others))\n",
    "\n",
    "        pmf = thinkstats2.Pmf(self.pool)\n",
    "        self.values = range(35, 44)\n",
    "        self.expected_probs = np.array(pmf.Probs(self.values))\n",
    "\n",
    "    def RunModel(self):\n",
    "        \"\"\"Run the model of the null hypothesis.\n",
    "\n",
    "        returns: simulated data\n",
    "        \"\"\"\n",
    "        np.random.shuffle(self.pool)\n",
    "        data = self.pool[:self.n], self.pool[self.n:]\n",
    "        return data\n",
    "# compare pregnancy lengths (chi-squared)\n",
    "data = firsts.prglngth.values, others.prglngth.values\n",
    "ht = PregLengthTest(data)\n",
    "p_value = ht.PValue()\n",
    "print('\\npregnancy length chi-squared')\n",
    "PrintTest(p_value, ht)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "用卡方统计量检验第一胎和其他胎的妊娠时间：全国家庭增长调查数据的总卡方统计量为102，这个数字本身没有什么意义。但在原假设下，1000次重复产生的最大验证统计量是30。对比之下，我们认为观察到的卡方统计量在原假设下不太可能出现，因此这个直观效应是统计显著的。\n",
    "\n",
    "这个例子说明卡方检验存在一个局限：卡方检验可以证明两个群组之间存在差异，但不能揭示这个差异是什么。\n",
    "\n",
    "## 误差\n",
    "- 误报率（false positive rate）：如果一个效应的确是偶然发生的，那么我们将它误判为统计显著的概率是多少？\n",
    "- 漏报率（false negative rate）：如果一个效应不是偶然发生的，那么假设检验失败的概率是多少？\n",
    "\n",
    "如果阈值为5%，那么漏报率就是5%，理由如下：\n",
    "- 如果效应不为真，那么原假设成立，因此，通过模拟原假设就可以计算出检验统计量的分布。我们将这个分布称为$CDF_T$。\n",
    "- 重复运行实验，每次得到一个来自$CDF_T$的测试统计量t。随后，我们计算出p值。p值是$CDF_T$中随机值大于t的概率，即$1-CDF_T(t)$。\n",
    "- 当$CDF_T(t)$大于95%，即t大于第95百分秩时，p值小于5%。而$CDF_T$中随机抽取的值大于第95百分秩的概率为多少？答案是5%。\n",
    "\n",
    "因此，如果你进行一个阈值为5%的假设检验，20次里会有1次漏报。\n",
    "\n",
    "## 功效\n",
    "误报率受实际效应大小的影响，而通常我们无法得知实际效应的大小，因此误报率较难计算。一个办法是计算一个假定效应大小的误报率。\n",
    "\n",
    "检验妊娠时间，结果约为70%。这个结果说明，如果妊娠时间均值的实际差异为0.78周，那么我们预期，如果使用这个规模的样本进行实验，结果有70%的可能性为误报。\n",
    "\n",
    "另一种方式描述结果：如果实际差异为0.78周，那么我们预期检验通过的可能性只有30%。这个“正确通过率”称为检验的功效(power)，有时也称为“敏感度”（sensitivity）。这个值反映了一个检验检测出指定大小效应的能力。\n",
    "\n",
    "在这个示例中，这个检验结果通过的可能性只有30%（假设实际差异为0.78周）。一般来说，我们认为80%的功效是可接受的，因此示例中的检验属于“低功效的”（underpowered）。\n",
    "\n",
    "通常，假设检验失败并不说明两个群组之间不存在差异，而是说，如果差异的确存在的话，这个差异太小，以至于无法在这种规模的样本中检测到。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "false neg rate 0.701\n"
     ]
    }
   ],
   "source": [
    "def FalseNegRate(data, num_runs=1000):\n",
    "    \"\"\"Computes the chance of a false negative based on resampling.\n",
    "\n",
    "    data: pair of sequences\n",
    "    num_runs: how many experiments to simulate\n",
    "\n",
    "    returns: float false negative rate\n",
    "    \"\"\"\n",
    "    group1, group2 = data\n",
    "    count = 0\n",
    "\n",
    "    for i in range(num_runs):\n",
    "        sample1 = thinkstats2.Resample(group1)\n",
    "        sample2 = thinkstats2.Resample(group2)\n",
    "        ht = DiffMeansPermute((sample1, sample2))\n",
    "        p_value = ht.PValue(iters=101)\n",
    "        if p_value > 0.05:\n",
    "            count += 1\n",
    "\n",
    "    return count / num_runs\n",
    "# compute the false negative rate for difference in pregnancy length\n",
    "data = firsts.prglngth.values, others.prglngth.values\n",
    "neg_rate = FalseNegRate(data)\n",
    "print('false neg rate', neg_rate)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 复现\n",
    "第一版使用的是全国家庭增长调查的第6周期数据，发布于2002年。2011年10月，疾病控制和预防中心发布了基于2006~2010年调查附加数据。nsfg2.py中包含读取和清洗这些数据的代码。在这个心数据集中：\n",
    "- 妊娠时间均值差为0.16周，p<0.001，是统计显著的（原始数据中均值差为0.078周）\n",
    "- 新生儿体重的差值为0.17磅，p<0.001（原始数据中差值为0.12磅）\n",
    "- 新生儿体重与母亲年龄之间的相关性为0.08，p<0.001（原始数据中相关性为0.07）\n",
    "- 卡方检验结果p<0.001，是统计显著的（原始数据也是如此）\n",
    "\n",
    "原始数据中所有统计显著的效应，在新数据集中得到了复现。此外，妊娠时间差值在原始数据中不显著，在新数据集中这个差值变大，由不显著变为显著。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "prglngth2\n",
      "means permute two-sided\n",
      "p-value = 0.0\n",
      "actual = 0.164708338667\n",
      "ts max = 0.14194274494\n",
      "\n",
      "birth weight 2\n",
      "means permute two-sided\n",
      "p-value = 0.0\n",
      "actual = 0.170870412946\n",
      "ts max = 0.0930818843361\n",
      "\n",
      "age weight correlation 2\n",
      "p-value = 0.0\n",
      "actual = 0.0808183040947\n",
      "ts max = 0.0305377196587\n",
      "\n",
      "pregnancy length chi-squared 2\n",
      "p-value = 0.0\n",
      "actual = 145.214547488\n",
      "ts max = 25.6167521177\n"
     ]
    }
   ],
   "source": [
    "def ReplicateTests():    \n",
    "    \"\"\"Replicates tests with the new NSFG data.\"\"\"\n",
    "    import nsfg2\n",
    "    live, firsts, others = nsfg2.MakeFrames()\n",
    "\n",
    "    # compare pregnancy lengths\n",
    "    print('\\nprglngth2')\n",
    "    data = firsts.prglngth.values, others.prglngth.values\n",
    "    ht = DiffMeansPermute(data)\n",
    "    p_value = ht.PValue(iters=1000)\n",
    "    print('means permute two-sided')\n",
    "    PrintTest(p_value, ht)\n",
    "\n",
    "    print('\\nbirth weight 2')\n",
    "    data = (firsts.totalwgt_lb.dropna().values,\n",
    "            others.totalwgt_lb.dropna().values)\n",
    "    ht = DiffMeansPermute(data)\n",
    "    p_value = ht.PValue(iters=1000)\n",
    "    print('means permute two-sided')\n",
    "    PrintTest(p_value, ht)\n",
    "\n",
    "    # test correlation\n",
    "    live2 = live.dropna(subset=['agepreg', 'totalwgt_lb'])\n",
    "    data = live2.agepreg.values, live2.totalwgt_lb.values\n",
    "    ht = CorrelationPermute(data)\n",
    "    p_value = ht.PValue()\n",
    "    print('\\nage weight correlation 2')\n",
    "    PrintTest(p_value, ht)\n",
    "\n",
    "    # compare pregnancy lengths (chi-squared)\n",
    "    data = firsts.prglngth.values, others.prglngth.values\n",
    "    ht = PregLengthTest(data)\n",
    "    p_value = ht.PValue()\n",
    "    print('\\npregnancy length chi-squared 2')\n",
    "    PrintTest(p_value, ht)\n",
    "# run the tests with new nsfg data\n",
    "ReplicateTests()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [Root]",
   "language": "python",
   "name": "Python [Root]"
  },
  "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}