{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0da5ab55-df7a-48d9-a9e8-a84473eb2a72",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Paired-Sample Tests\n",
    "\n",
    "In [The Design of Experiments](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2458144/https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2458144/), Fisher describes an experiment conducted by Charles Darwin to measure the effect of cross-fertilization on plant growth. In Darwin's experiment, pairs of plants were grown from the same batch of seed in the same pot under the same conditions, except that one was self-fertilized and the other was cross-fertilized. \n",
    "\n",
    "> The evident object of these precautions is to increase the sensitiveness of the experiment, by making such differences in growth rate as were to be observed as little as possible dependent from environmental circumstancces, and as much as possible, therefore, from intrinsic differences due to their mode of origin.\n",
    "\n",
    "The `x` and `y` arrays below record the height of the cross-fertilized and self-fertilized plants, respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "77e2b91f-6d6c-43c0-b6f1-0be309885c9d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "x = np.array([23.5, 12, 21, 22, 19.125, 21.5, 22.125, 20.375, 18.25, 21.625, 23.25, 21, 22.125, 23, 12])  # (in) cross-fertilized\n",
    "y = np.array([17.375, 20.375, 20, 20, 18.375, 18.625, 18.625, 15.25, 16.5, 18, 16.25, 18, 12.75, 15.5, 18])  # (in) self-fertilized\n",
    "assert len(x) == len(y)  # elements at corresponding positions form a pair"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "24b37c06-fb34-493e-ab4b-c49a14185c69",
   "metadata": {},
   "source": [
    "The null hypothesis was that the method of fertilization would have no effect, and Darwin's alternative hypothesis was that cross-fertilized plants would be taller. Fisher recommends the t-test, which is implemented for paired (\"related\") samples by `stats.ttest_rel`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "005bb7e7-8751-4298-a13e-c796e149df9e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Ttest_relResult(statistic=2.1479874613311205, pvalue=0.024851472010900447)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy import stats\n",
    "# alternative hypothesis 'greater': the mean of the distribution underlying \n",
    "# `x` is greater than the mean of the distribution underlying `y` \n",
    "res_t = stats.ttest_rel(x, y, alternative='greater')\n",
    "res_t"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "809a24f6-4345-47f5-b7b8-9c0b64b7384f",
   "metadata": {},
   "source": [
    "However, the t-test is derived under the assumption that samples were drawn from a normal population. He reports that \n",
    "\n",
    "> There has, however, in recent years, been a tendency for theoretical statisticians, not closely in touch with the requirements of experimental data, to stress the element of normality, in the hypothesis tested, as though it were a serious limitation to the test applied.\n",
    "\n",
    "(The same is true today!) And so he proposes a different test that does not rely on the assumption of normality:\n",
    "\n",
    "> On the hypothesis that the two series of seeds are random sampled from identical populations, and that their sites have been assigned to members of each pair independently at random, the 15 differences [in height between pairs] would each have occured with equal frequency with a positive or with a negative sign... \n",
    "Since *ex hypothesi* each of these $2^{15}$ combinations will occur by chance with equal frequency, a knowledge of how many of them are equal to or greater than the the value actually observed affords a direct arithmetical test of the significance of this value.\n",
    "\n",
    "He continues to perform the following test."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "5c328b74-1774-4008-a8fc-d922fa0d94b9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Observed Frequency')"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEGCAYAAACtqQjWAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAAAjeklEQVR4nO3df5xddX3n8debIKCtbLSMNgZwAhvE4FLAEaiKP1pdE7BGfYCG2oJom8aCW9e2a1C3ujy6u6z4o1KRiBIBqwT8Uc0WLCKt0roNMEFEAkRDiBJIIZUVUFhi8L1/nO/ozc2dmTOTe2bunXk/H4/7mHvO98f93Dv33s8953vO98g2ERER3bDXdAcQEREzR5JKRER0TZJKRER0TZJKRER0TZJKRER0zd7THcB0OuCAAzw4ODjdYURE9JX169f/m+2BTmWzOqkMDg4yPDw83WFERPQVST8YrSy7vyIiomsaTSqSFkvaKGmTpJUdyiXp/FJ+q6RjyvqDJP2jpDskbZD0Jy1tni7pWknfL3+f1lJ2dulro6RXNfncIiJid40lFUlzgAuAJcAi4FRJi9qqLQEWltty4MKyfifwp7afCxwPnNnSdiVwne2FwHVlmVK+DDgCWAx8vMQQERFTpMktlWOBTbY3294BrAGWttVZClzmyjpgrqR5trfZvhnA9iPAHcD8ljaXlvuXAq9tWb/G9uO27wY2lRgiImKKNJlU5gP3tCxv5ZeJoXYdSYPA0cANZdUzbW8DKH+fMYHHQ9JyScOShrdv3z6R5xMREeNoMqmow7r22SvHrCPpV4EvAu+w/XAXHg/bF9kesj00MNDxiLiIiJikJpPKVuCgluUDgfvq1pH0JKqE8lnbX2qpc7+keaXOPOCBCTxeREQ0qMmkchOwUNICSftQDaKvbauzFjitHAV2PPCQ7W2SBFwM3GH7wx3anF7unw58pWX9Mkn7SlpANfh/Y/efVkREjKaxkx9t75R0FnANMAdYbXuDpBWlfBVwNXAi1aD6o8AZpfmLgN8HvivplrLu3bavBs4FrpT0VuCHwCmlvw2SrgRupzp67EzbTzT1/CIiYneazRfpGhoacs6oj14yuPIqtpx70m73I3qJpPW2hzqV5Yz6iIjomiSViGkwuPKqrtaL6BVJKhER0TVJKhER0TVJKhER0TVJKhFTqM4YyeDKq2rXi+g1SSoRU6BTophMUujUR5JL9JIklYiI6JoklYgeMNbWRrZEop/M6mvUR0y3bp6vkjPwoxdkSyViio2WICayRZKtl+hVSSoREdE1SSoRDZmOrYlswcR0S1KJmAGSTKJXJKlEdNFUfLkngUQvS1KJiIiuySHFEQ3KVkXMNo1uqUhaLGmjpE2SVnYol6TzS/mtko5pKVst6QFJt7W1uULSLeW2ZeRyw5IGJT3WUraqyecW0aqXk0cvxxYzT2NbKpLmABcArwS2AjdJWmv79pZqS4CF5XYccGH5C3AJ8DHgstZ+bb+x5TE+BDzUUnyX7aO6+kQiahr58p7uL/GRx8+JkDEdmtxSORbYZHuz7R3AGmBpW52lwGWurAPmSpoHYPt64MHROpck4A3A5Y1EHxERE9ZkUpkP3NOyvLWsm2id0ZwA3G/7+y3rFkj6tqRvSjqhUyNJyyUNSxrevn17zYeKiIg6mkwq6rDOk6gzmlPZdStlG3Cw7aOBdwKfk7T/bp3bF9kesj00MDBQ86Ei+s9074aL2anJpLIVOKhl+UDgvknU2Y2kvYHXA1eMrLP9uO0flfvrgbuAwyYVecQE9PKXd6+M88Ts0WRSuQlYKGmBpH2AZcDatjprgdPKUWDHAw/Z3laj71cAd9reOrJC0kA5OABJh1AN/m/uxhOJmCmSXKJpjR39ZXunpLOAa4A5wGrbGyStKOWrgKuBE4FNwKPAGSPtJV0OvAw4QNJW4H22Ly7Fy9h9gP4lwDmSdgJPACtsjzrQHxER3Se77hDGzDM0NOTh4eHpDiP6WD/+8s+hxrGnJK23PdSpLNO0RExSPyYU6N+4oz8kqURMQr6YIzpLUomYoJmQUGbCc4jelKQSERFdk6QSERFdk6QSUUN2F0XUk6QSMcslYUY3JalERETX5MqPETXNtF/0M+35RG/IlkrEOPLlG1FfkkpERHRNkkpERHRNkkpERHRNkkpERHRNkkpERHRNkkrEKGbTUV+57HB0S5JKRER0TZJKxBhm0y/32fRcozmNJhVJiyVtlLRJ0soO5ZJ0fim/VdIxLWWrJT0g6ba2Nu+XdK+kW8rtxJays0tfGyW9qsnnFhERu2ssqUiaA1wALAEWAadKWtRWbQmwsNyWAxe2lF0CLB6l+4/YPqrcri6PtwhYBhxR2n28xBAxYbP5V/vgyqtm9fOPPdPklsqxwCbbm23vANYAS9vqLAUuc2UdMFfSPADb1wMPTuDxlgJrbD9u+25gU4khIiKmyLhJRdLzJtn3fOCeluWtZd1E63RyVtldtlrS0ybSl6TlkoYlDW/fvr3GQ8Vskl/pEXumzpbKKkk3SvpjSXMn0Lc6rPMk6rS7EDgUOArYBnxoIn3Zvsj2kO2hgYGBcR4qZpMkk93lNYmJGjep2H4x8CbgIGBY0uckvbJG31tLmxEHAvdNok57PPfbfsL2z4FP8stdXBPuKyIiuqvWmIrt7wPvBd4FvBQ4X9Kdkl4/RrObgIWSFkjah2oQfW1bnbXAaeUosOOBh2xvGyuWkTGX4nXAyNFha4FlkvaVtIBq8P/GOs8vIr/IR5fXJiZi3It0SToSOAM4CbgW+B3bN0t6FvAvwJc6tbO9U9JZwDXAHGC17Q2SVpTyVcDVwIlUg+qPlscZedzLgZcBB0jaCrzP9sXAByQdRbVrawvwR6W/DZKuBG4HdgJn2n5iQq9GRETsEdljD2FIup5qN9MXbD/WVvb7tj/TYHyNGhoa8vDw8HSHET0gv8bHtuXck6Y7hOghktbbHupUVudywicCj4386pe0F7Cf7Uf7OaFERET31RlT+Trw5Jblp5R1ERERu6iTVPaz/ZORhXL/Kc2FFBER/apOUvlp25xczwceG6N+RMxAGXeKOuqMqbwD+LykkXM+5gFvbCyiiCmUL8p68jpFXXVOfrwJOBx4G/DHwHNtr286sIjoPUkuMZ46WyoALwAGS/2jJWH7ssaiioiIvlTn5MfPUM21dQswcjKhgSSViIjYRZ0tlSFgkcc7SzKiz2RXTkT31Tn66zbg15sOJCIi+l+dLZUDgNsl3Qg8PrLS9msaiyoiIvpSnaTy/qaDiIiImWHcpGL7m5KeDSy0/XVJT6GadTgiImIXdS4n/IfAF4BPlFXzgS83GFNE17UOymeAfs/lNYzR1BmoPxN4EfAw/OKCXc9oMqiIiOhPdZLK47Z3jCxI2pvxryMfERGzUJ2k8k1J7waeXK5N/3ngfzcbVkSzsvtm8vLaxVjqJJWVwHbgu1SX7r2a6nr145K0WNJGSZskrexQLknnl/Jb22ZDXi3pAUm3tbU5T9Kdpf7fSppb1g9KekzSLeW2qk6MMXsMrrwqX4hdlNcyOqkzoeTPbX/S9im2Ty73x939JWkOcAGwBFgEnCppUVu1JcDCclsOXNhSdgmwuEPX1wLPs30k8D3g7Jayu2wfVW4rxosxIvZMEku0qzP31910GEOxfcg4TY8FNtneXPpZAywFbm+psxS4rCSpdZLmSppne5vt6yUNdnjcr7UsrgNOHu85RETE1Kg799eI/YBTgKfXaDcfuKdleStwXI0684FtNfoHeAtwRcvyAknfpjpS7b22/6m9gaTlVFtFHHzwwTUfJiIi6qiz++tHLbd7bf8V8Fs1+lan7iZRp3Pn0nuAncBny6ptwMG2jwbeCXxO0v67dW5fZHvI9tDAwECdh4o+l100zctrHCPq7P46pmVxL6otl6fW6HsrcFDL8oHAfZOo0ymm04FXA789Mr5j+3HK3GS210u6CzgMGK4Ra0REdEGd3V8farm/E9gCvKFGu5uAhZIWAPcCy4DfbauzFjirjLccBzxke8xdX5IWA+8CXmr70Zb1A8CDtp+QdAjV4P/mGnHGDJZf0FNncOVVbDn3pOkOI6ZZnbm/Xj6Zjm3vlHQWcA3VXGGrbW+QtKKUr6I6PPlEYBPwKHDGSHtJlwMvAw6QtBV4n+2LgY8B+wLXSgJYV470eglwjqSdVBcTW2H7wcnEHhERk6Pxjg6W9M6xym1/uKsRTaGhoSEPD2fv2EyVrZTpka2VmU/SettDncrqHv31AqpdVQC/A1zPrkdtRURE1L5I1zG2HwGQ9H7g87b/oMnAIiKi/9SZpuVgYEfL8g5gsJFoIiKir9VJKp8BbpT0fknvA24ALms2rIjoVxnLmt3qHP313yV9FTihrDrD9rebDSsiIvpRnS0VgKcAD9v+KLC1nHsS0bPya3l65fWfvepcTvh9VCcbjswG/CTgb5oMKiIi+lOdLZXXAa8Bfgpg+z7qTdMSERGzTJ2ksqPMr2UASb/SbEgREdGv6iSVKyV9Apgr6Q+BrwOfbDasiIjoR2Me/aVqcq0rgMOprlHyHOAvbF87BbFFTEoGiSOmz5hJxbYlfdn286ku4xsRETGqOru/1kl6QeORROyhbKH0nvxPZp86c3+9HFghaQvVEWCi2og5ssnAIiYjX2K9I/+L2WnUpCLpYNs/BJZMYTwREdHHxtr99WUA2z8APmz7B623KYkuYhz5NRzRW8ZKKmq5f0jTgURERP8bK6l4lPsREREdjZVUfkPSw5IeAY4s9x+W9Iikh+t0LmmxpI2SNkla2aFcks4v5bdKOqalbLWkByTd1tbm6ZKulfT98vdpLWVnl742SnpVnRijf43s+hpceVV2g/Ww/J9ml1GTiu05tve3/VTbe5f7I8v7j9expDnABVQD/YuAUyUtaqu2BFhYbsuBC1vKLgEWd+h6JXCd7YXAdWWZ0vcy4IjS7uMlhoiImCJ1p76fjGOBTbY3294BrAGWttVZClzmyjqqqWDmAdi+HniwQ79LgUvL/UuB17asX2P7cdt3A5tKDBERMUWaTCrzgXtalreWdROt0+6ZtrcBlL/PmEhfkpZLGpY0vH379nGfRERE1NdkUlGHde0D/nXqdPPxsH2R7SHbQwMDA5N8qIiI6KTJpLIVOKhl+UDgvknUaXf/yC6y8veBPegrIqZABuhnj1GTyshRXqPdavR9E7BQ0gJJ+1ANoq9tq7MWOK0cBXY88NDIrq0xrAVOL/dPB77Ssn6ZpH3L5Y4XAjfWiDMiIrpk1GlabD8VQNI5wL8Cn6HaxfQmalz50fZOSWcB1wBzgNW2N0haUcpXAVcDJ1INqj8KnDHSXtLlwMuAAyRtBd5n+2LgXKprvLwV+CFwSulvg6QrgduBncCZtp+o/1JEv8iv3v428v/bcu5J0xxJNEHVRR3HqCDdYPu48db1o6GhIQ8PD093GDFBSSozQ5JK/5K03vZQp7I6YypPSHqTpDmS9pL0JiBbABGxR/LjYGaqk1R+F3gDcH+5nVLWRUy5fBFF9LZxr6diewu7n7QYERGxm3G3VCQdJum6kTm4JB0p6b3NhxYREf2mzu6vTwJnAz8DsH0r1eHBEVMqu75mnvxPZ546SeUpttvP99jZRDAREdHf6iSVf5N0KGXKE0knA+OdoBgREbPQuAP1wJnARcDhku4F7qY6ATIiImIXdZLKD2y/QtKvAHvZfqTpoCIioj/V2f11t6SLgOOBnzQcT8QuMpA7e+R/PTPUSSrPAb5OtRvsbkkfk/TiZsOKiIh+NG5Ssf2Y7Sttvx44Gtgf+GbjkUVERN+pdT0VSS+V9HHgZmA/qmlbIhqV3SER/afOGfV3A+8A/gl4nu032P5i04FFtEqCmbnyv51Zxjz6S9Ic4NO2z5mieCIioo+NuaVSLnL18imKJSIi+lyd81T+j6SPAVcAPx1ZafvmxqKKaJHdIzNf/sczR52B+hcCRwDnAB8qtw/W6VzSYkkbJW2StLJDuSSdX8pvlXTMeG0lXSHplnLbIumWsn5Q0mMtZavqxBi9J18wkfdA/6pzPZVJ7f4q4zEXAK8EtgI3SVpr+/aWakuAheV2HHAhcNxYbW2/seUxPgQ81NLfXbaPmky8ERGx5+oc/fVMSRdL+mpZXiTprTX6PhbYZHuz7R3AGna/2NdS4DJX1gFzJc2r01aSqA5tvrxGLBERMQXq7P66BLgGeFZZ/h7VIcbjmQ/c07K8tayrU6dO2xOA+21/v2XdAknflvRNSSd0CkrScknDkoa3b99e42nEdMjuj9kp//f+VyepHGD7SuDnALZ3Ak/UaKcO61yzTp22p7LrVso24GDbRwPvBD4naf/dOrEvsj1ke2hgYGDU4GN65EslRuS90J/qJJWfSvo1fnk9lePZdRxjNFuBg1qWDwTuq1lnzLaS9gZeT3VEGgC2H7f9o3J/PXAXcFiNOKNH5EskYNf3Qd4T/adOUnknsBY4VNK3gMuAt9dodxOwUNICSftQXYJ4bVudtcBp5Siw44GHbG+r0fYVwJ22t46skDRQBviRdAjV4P/mGnFGRESX1Dn662ZJL6WarVjARts/q9Fup6SzqMZj5gCrbW+QtKKUrwKuBk4ENgGPAmeM1bal+2XsPkD/EuAcSSO751bYfnC8OCMionvGTSqSTgH+viSE9wLHSPrLOic/2r6aKnG0rlvVct9UU+rXattS9uYO674IZE6yiIhpVGf313+1/Ui5hsqrgEupzieJ6JrsO4+YGeoklZEjvU4CLrT9FWCf5kKKiIh+VSep3CvpE1QnGl4tad+a7SIiYpapkxzeQDVgvtj2j4GnA3/eZFAREdGf6lxO+FFgC7BE0tuBeba/1nRgERGQ8bZ+U2fur7+gGpz/NeAA4NPlKLCIiCmRxNI/6lxP5VTgaNv/D0DSuVTXqv/LJgOLmW/ki2LLuSdNcyTRD/J+6Q91xlS2APu1LO9LNQVKRETELkbdUpH011TzfT0ObJB0bVl+JfDPUxNezAbZtRExc4y1+2u4/F0P/G3L+m80Fk1ERPS1UXd/2b7U9qVUMwGvp0oyV7Ssj5iQzD4b3ZD3UW8bNalI2lvSB6imob8U+BvgHkkfkPSkqQowIiL6x1gD9edRnei4wPbzy8WvDgXmAh+cgtgiIqLPjJVUXg38oe1HRlbYfhh4G9V09RGTkl0WETPXWEnFZWr69pVPsPulfSMiIsZMKrdLOq19paTfA+5sLqSIiOhXYx1SfCbwJUlvoTr6y8ALgCcDr5uC2GIGyq6v6Ia8j3rXqEnF9r3AcZJ+CziC6lLCX7V93VQFFxER/aXOLMX/YPuvbZ8/0YQiabGkjZI2SVrZoVySzi/lt0o6Zry2kt4v6V5Jt5TbiS1lZ5f6GyW9aiKxRkTEnqszoeSkSJoDXEA1rctW4CZJa23f3lJtCbCw3I6jukzxcTXafsT2Loc1S1oELKPaqnoW8HVJh5UDCyIiYgo0eQXHY4FNtjfb3gGsAZa21VkKXObKOmCupHk127ZbCqyx/bjtu4FNpZ+IiJgiTSaV+cA9Lctby7o6dcZre1bZXbZa0tMm8HhIWi5pWNLw9u3bJ/J8Yg9kYDWaMvLeynusNzSZVNRhXfv5LaPVGavthVRn9h8FbAM+NIHHw/ZFtodsDw0MDHRoEt2UD3rE7NLYmArVlsJBLcsHAvfVrLPPaG1t3z+yUtIngb+bwOPFFEpCiZh9mtxSuQlYKGmBpH2oBtHXttVZC5xWjgI7HnjI9rax2pYxlxGvA25r6WuZpH0lLaAa/L+xqScXERG7a2xLxfZOSWcB1wBzgNW2N0haUcpXAVdTzSO2CXgUOGOstqXrD0g6imrX1hbgj0qbDZKuBG4HdgJn5siv3pAtlmha3mO9Qx2m95o1hoaGPDw8PH7FmJDBlVex5dyT8kGPaTHy3su17Jsjab3toU5lTe7+ilksCSWmS9570ytJJSIiuiZJJSIiuiZJJSIiuiZJJboq+7OjVwyuvCrvx2mQpBJdkQ9w9Kq8L6dWkkpERHRNkkrssfwSjF6X9+jUSVKJiIiuSVKJiIiuSVKJiIiuSVKJCWs90iv7qqNf5D07NZJUIiKia5JUIiKia5JUYtKyGyH6Td6zzUtSiYiIrklSiYhZJ1sszWnscsIAkhYDH6W6JPCnbJ/bVq5SfiLV5YTfbPvmsdpKOg/4HWAHcBdwhu0fSxoE7gA2lu7X2V7R5PObTfIhjJmm9T2dq0R2T2NbKpLmABcAS4BFwKmSFrVVWwIsLLflwIU12l4LPM/2kcD3gLNb+rvL9lHlloQSEbXkR1P3NLn761hgk+3NtncAa4ClbXWWApe5sg6YK2neWG1tf832ztJ+HXBgg88hIiImoMmkMh+4p2V5a1lXp06dtgBvAb7asrxA0rclfVPSCZ2CkrRc0rCk4e3bt9d7JhERUUuTSUUd1rlmnXHbSnoPsBP4bFm1DTjY9tHAO4HPSdp/t07si2wP2R4aGBgY5ylErpMSs0Xe593RZFLZChzUsnwgcF/NOmO2lXQ68GrgTbYNYPtx2z8q99dTDeIf1pVnMotkKouY7fLe3zNNJpWbgIWSFkjaB1gGrG2rsxY4TZXjgYdsbxurbTkq7F3Aa2w/OtKRpIEywI+kQ6gG/zc3+PwiIqJNY4cU294p6SzgGqrDglfb3iBpRSlfBVxNdTjxJqpDis8Yq23p+mPAvsC11RHJvzh0+CXAOZJ2Ak8AK2w/2NTzmw3yiy1mm7zn91yj56nYvpoqcbSuW9Vy38CZdduW9f9+lPpfBL64J/HOdvlARcSeyhn1sZskl5jtMrY4eUkqs1yO7ororPVzkc9JfUkqERHRNUkqAWQzP2I0+WxMTJLKLJUPSkQ0IUklIqKmjK2ML0klIiK6Jklllmg/kiUiJi+fodElqcxQdZJIPhgRe679czTbP1dJKjNY3uwRzcpnbHcqk/zOSkNDQx4eHp7uMBqRN3fE9JnplyeWtN72UKeybKlERETXJKlERDRgtu4tSFKZQTIJXkRvmM2fxSSVPtP6Zu10ItZsfBNH9JOZnnCSVCIiomty9FcfHf01U3/ZRMxm/XikWI7+6lFj7bpq3bWVZBIxc431We/H3duNJhVJiyVtlLRJ0soO5ZJ0fim/VdIx47WV9HRJ10r6fvn7tJays0v9jZJe1eRz21N13xxJLBGzw0Q/6736ndBYUpE0B7gAWAIsAk6VtKit2hJgYbktBy6s0XYlcJ3thcB1ZZlSvgw4AlgMfLz0ExERU6SxMRVJvwm83/aryvLZALb/Z0udTwDfsH15Wd4IvAwYHK3tSB3b2yTNK+2f096/pGtKH/8yWozdHFMZXHnVLvtGe/VXRETMLu1jNiPfTXsyljPWmMrek+51fPOBe1qWtwLH1agzf5y2z7S9DaAklme09LWuQ1+7kLScaqsI4CclSXWF/tekmh0A/Fu3Ypgi/RZzv8UL/Rdz4m3epGIe7Xtpkt9XI549WkGTSUUd1rVvFo1Wp07byTweti8CLhqnrykjaXi0jN+r+i3mfosX+i/mxNu8fom5yYH6rcBBLcsHAvfVrDNW2/vLbi/K3wcm8HgREdGgJpPKTcBCSQsk7UM1iL62rc5a4LRyFNjxwENl19ZYbdcCp5f7pwNfaVm/TNK+khZQDf7f2NSTi4iI3TW2+8v2TklnAdcAc4DVtjdIWlHKVwFXAycCm4BHgTPGalu6Phe4UtJbgR8Cp5Q2GyRdCdwO7ATOtP1EU8+vi3pmV9wE9FvM/RYv9F/Mibd5fRHzrD6jPiIiuitn1EdERNckqURERNckqUwjSW8vU8pskPSBlvU9O92MpD+TZEkHtKzryXglnSfpzjIF0N9KmttS1qsxjzm10XSTdJCkf5R0R3nf/klZP+r0Sb1C0hxJ35b0d2W5Z2OWNFfSF8r79w5Jv9nL8bZKUpkmkl4OLAWOtH0E8MGyvmenm5F0EPBKqgMkRtb1bLzAtcDzbB8JfA84G3o35ppTG023ncCf2n4ucDxwZomx4/RJPeZPgDtalns55o8Cf2/7cOA3qOLu5Xh/IUll+rwNONf24wC2R863WQqssf247bupjow7dppibPcR4L+w60mlPRuv7a/Z3lkW11GduwS9G/OxwCbbm23vANZQxdozbG+zfXO5/wjVl918qjgvLdUuBV47LQGOQtKBwEnAp1pW92TMkvYHXgJcDGB7h+0f06PxtktSmT6HASdIukHSNyW9oKwfbeqaaSXpNcC9tr/TVtST8XbwFuCr5X6vxtyrcXUkaRA4GriBtumTgGeM0XQ6/BXVD6Kft6zr1ZgPAbYDny676z4l6Vfo3Xh30eQ0LbOepK8Dv96h6D1Ur/3TqHYhvIDq3JtDmNwUNV0xTrzvBv5jp2Yd1k3ZcepjxWz7K6XOe6h223x2pFmH+r1wbH2vxrUbSb8KfBF4h+2HpU6h9wZJrwYesL1e0sumOZw69gaOAd5u+wZJH6VHd3V1kqTSINuvGK1M0tuAL7k6UehGST+nmjBu2qabGS1eSf8BWAB8p3x5HAjcLOlYpnl6nLFeYwBJpwOvBn7bvzwpq1en9OnVuHYh6UlUCeWztr9UVt8vaV7L7OEPjN7DlHsR8BpJJwL7AftL+ht6N+atwFbbN5TlL1AllV6NdxfZ/TV9vgz8FoCkw4B9qGYg7bnpZmx/1/YzbA/aHqR60x9j+1/pwXhHSFoMvAt4je1HW4p6NeY6UxtNK1W/Ki4G7rD94Zai0aZPmna2z7Z9YHnvLgP+wfbv0aMxl8/VPZKeU1b9NtVMIT0Zb7tsqUyf1cBqSbcBO4DTyy/pvppupsenx/kYsC9wbdnCWmd7Ra/GPM70RL3iRcDvA9+VdEtZ925GmT6px/VyzG8HPlt+XGymmsJqL3o33l/INC0REdE12f0VERFdk6QSERFdk6QSERFdk6QSERFdk6QSERFdk6QSM4Kkb7TPNizpHZI+Pk6boYbjurzMkvyf29a/drKTRUo6qpzIV7f+XEl/PNF6kp4l6Qvdqh+zQ5JKzBSXU53Y1mpZWT8tJP068ELbR9r+SFvxa6lmIp6Mo6guw13XXGDcpNJez/Z9tk/uYv2YBZJUYqb4AvBqSfvCLyY7fBbwz5IulDRcrv/x3zo1lvSTlvsnS7qk3B+Q9EVJN5Xbizq03U/SpyV9t0wA+PJS9DXgGZJukXRCS/0XAq8Bzitlh5bb30taL+mfJB1e6p4i6TZJ35F0fTkZ7hzgjaXtG9tiOULSjaXsVkkLqU7yO7SsO0/Sr0q6TtLNJeaRmZDb6w2Wk3Pr9ttaf46kD5b+b5X09tr/yehvtnPLbUbcgKuApeX+SuC8cv/p5e8c4BtU17Ch3B8q93/S0s/JwCXl/ueAF5f7B1NNT9L+uH8KfLrcP5zqbOf9gEHgtlFivQQ4uWX5OmBhuX8c1VQiAN8F5pf7c8vfNwMfG6XfvwbeVO7vAzy5PQ6qmTT2L/cPoJr6Xx3q/WK5Zr+t9d9GNT/Y3q3/g9xm/i3TtMRMMrIL7Cvl71vK+jdIWk71ZTqParfTrTX7fAWwSL+chXd/SU91dS2RES+m+tLF9p2SfkB1aYOH6zyAqhl/Xwh8vuVx9i1/vwVcUqaV+VKH5u3+BXiPquuHfMn297X7DMIC/oekl1BNBT8feGYX+m31CmCVy/VsbD9YI/aYAZJUYib5MvBhSccAT7Z9c5kw8s+AF9j+v2W31n4d2rbOV9Ravhfwm7YfG+Nx93Te972AH9s+areg7BWSjqO6wNQtknar01b/c5JuKPWvkfQHVHNHtXoTMAA83/bPJG2h82sy0X5biR6dtj+alTGVmDFs/4Rql9ZqfjlAvz/wU+AhSc+kulRvJ/dLeq6kvYDXtaz/GnDWyMIoX+rXU31Rj8w4fTCwcZxwHwGeWuJ+GLhb0imlD0n6jXL/UNs32P4LqlmsD2pt207VNXk22z6falbbIzvU/3dU1xf5WRn/eXZ7TJPst9XXgBWS9i7tnz7O6xEzRJJKzDSXU13Tew2AqytVfhvYQJVsvjVKu5XA3wH/AGxrWf+fgKEy2Hw7sKJD248DcyR9F7gCeLPLZaLHsAb48zKwfyhVUnqrpO+UWEcGz88rg923USWv7wD/SLVLbreBeuCNwG2qZhA+HLjM9o+Ab5UB//OoLlY2JGm4PO6d5bVqrzfRflt9imps6dbynH53nNcjZojMUhwREV2TLZWIiOiaJJWIiOiaJJWIiOiaJJWIiOiaJJWIiOiaJJWIiOiaJJWIiOia/w/6uzjBlVsrLQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from itertools import product\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def statistic(d):\n",
    "    # d is the differences in height between paired plants\n",
    "    # the statistic is the sum of these differences\n",
    "    return np.sum(d)\n",
    "\n",
    "# compute the statistic for all possible combinations of signs on the elements of `d`\n",
    "def null_distribution(d):\n",
    "    signs = (-1, 1) \n",
    "    n = len(d)  # number of observations per sample\n",
    "    null_distribution = []\n",
    "    for dsigns in product(*[signs]*n):\n",
    "        stat = statistic(d * dsigns)\n",
    "        null_distribution.append(stat)\n",
    "    return null_distribution\n",
    "\n",
    "d = x - y\n",
    "null_dist = null_distribution(d)\n",
    "assert len(null_dist) == 2**15\n",
    "bins = np.unique(null_dist).tolist()\n",
    "bins.append(np.max(null_dist)+1)\n",
    "plt.hist(null_dist, density=True, bins=bins)\n",
    "plt.xlabel(\"Value of test statistic\")\n",
    "plt.ylabel(\"Observed Frequency\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fa460986-ab98-448f-878d-35abe9e7c82a",
   "metadata": {},
   "source": [
    "> In just 863 cases out of 32,768 the total deviation will have a positive value as great as or greater than that observed.\n",
    "\n",
    "The $p$-value is simply the ratio of the two."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "ee2ae5b0-8ed4-4f20-be2d-904365f7ff6a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "t-test p-value: 0.024851472010900447\n",
      "permutation test p-value: 0.026336669921875\n"
     ]
    }
   ],
   "source": [
    "assert np.sum(null_dist >= statistic(d)) == 863\n",
    "pvalue = np.sum(null_dist >= statistic(d)) / len(null_dist)\n",
    "print(f\"t-test p-value: {res_t.pvalue}\")\n",
    "print(f\"permutation test p-value: {pvalue}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e105edca-1c7e-4223-ab29-ed5db55fcf4f",
   "metadata": {},
   "source": [
    "The two tests agree remarkably well in this case, suggesting the applicability of the t-test even when the original data are not strictly normally distributed."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bbcb4fdb-7191-4bc6-82c0-0fdfac445b4d",
   "metadata": {},
   "source": [
    "`permutation_test` minimizes the code required to perform the same procedure. The most essential difference here compared to the independent \n",
    "sample tests above is that we pass `permutation_type='samples'` instead of `permutation_type='independent'`. Note also that we can simply pass `np.sum` as the statistic because it satisfies the required interface, even with `vectorized=True`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "dd6d90a0-e443-4ec2-a6c5-04eddccea877",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Note that the first argument of `permutation_test` is a sequence containing all the samples.\n",
    "# We only have one sample here, but it still needs to be in a sequence.\n",
    "res = stats.permutation_test((d,), np.sum, alternative='greater', permutation_type='samples', vectorized=True, n_resamples=np.inf)\n",
    "assert res.pvalue == pvalue"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ef5fc8c-00a0-4f70-b385-76f189b1f66b",
   "metadata": {},
   "source": [
    "The value of the `permutation_type`, \"samples\", comes from the fact that permuting the signs of the differences is equivalent to permuting the *sample* to which each paired observation is assigned. That is, we can obtain the same $p$-value by performing a two-sample test in which observations are exchanged between the two samples in all possible ways *without* breaking up pairs. This is what `permutation_test` does with `permutation_type='samples'` when there are two or more samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "c27182a9-ba9d-4b7f-8fb4-66b49eee11bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "def statistic(x, y, axis=0):\n",
    "    # observations will be permuted between samples `x` and `y`,\n",
    "    # changing the sign of the corresponding element of `d` \n",
    "    d = x - y\n",
    "    return np.sum(d, axis=axis)\n",
    "\n",
    "res = stats.permutation_test((x, y), statistic, alternative='greater', permutation_type='samples', n_resamples=np.inf)\n",
    "assert res.pvalue == pvalue"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe585e67-6955-43a5-a509-0b3e1743d22b",
   "metadata": {},
   "source": [
    "Other common hypothesis tests can be performed as paired-sample permutation tests. We continue with several comparisons against `scipy.stats.wilcoxon` to further illustrate the usage of `permutation_test` with `permutation_type='samples'`.\n",
    "\n",
    "#### Two-sample Test"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d6b915f5-6182-4060-806c-60aa242071ae",
   "metadata": {},
   "source": [
    "[`wilcoxon`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.wilcoxon.html) is to `ttest_rel` as `mannwhitneyu` is to `ttest_ind`; it is a nonparametric permutation test of the null hypothesis that each observation in a pair is drawn from the same distribution. Because it is a paired-sample test, it computes the $p$-value by reassigning observations between the samples in all possible ways while maintaining the pairs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "4b4c4b12-0e1d-4a41-a202-3682cae871b1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "WilcoxonResult(statistic=10.0, pvalue=0.765625)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rng = np.random.default_rng()\n",
    "x = rng.random(7)\n",
    "y = rng.random(7)\n",
    "# Each element in `x` is paired with the corresponding element of `y` \n",
    "stats.wilcoxon(x, y, alternative='greater')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d542cee8-1f20-440f-be69-9812f593099a",
   "metadata": {},
   "source": [
    "Using `wilcoxon` to compute only the statistic, we can use `permutation_test` to calculate precisely the same $p$-value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "65a8170d-26be-4c64-9be3-d89e49699cc8",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\matth\\Desktop\\scipy\\scipy\\stats\\_morestats.py:3380: UserWarning: Sample size too small for normal approximation.\n",
      "  warnings.warn(\"Sample size too small for normal approximation.\")\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.765625"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def statistic(x, y):\n",
    "    return stats.wilcoxon(x, y, alternative='greater', method='approx').statistic\n",
    "\n",
    "res = stats.permutation_test((x, y), statistic, alternative='greater', permutation_type='samples')\n",
    "res.pvalue"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8bd72655-a67b-4493-aa7b-54a4cee5ef98",
   "metadata": {},
   "source": [
    "(The warning can be safely ignored because it is referring to the $p$-value returned by the test, whereas we are using only the statistic value.)\n",
    "\n",
    "Again, the advantage of `permutation_test` is its ability to handles ties in the data; `wilcoxon` does not."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "3f9a43ec-861c-418d-b736-646075f65785",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\matth\\Desktop\\scipy\\scipy\\stats\\_morestats.py:3366: UserWarning: Exact p-value calculation does not work if there are zeros. Switching to normal approximation.\n",
      "  warnings.warn(\"Exact p-value calculation does not work if there are \"\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "WilcoxonResult(statistic=4.5, pvalue=0.8539232992870668)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = rng.integers(0, 5, size=7)\n",
    "y = rng.integers(0, 5, size=7)\n",
    "stats.wilcoxon(x, y, method='exact')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "226a7dd7-7dc9-4ef1-9161-b2a59bf324fd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(4.5, 1.0)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res = stats.permutation_test((x, y), statistic, permutation_type='samples')\n",
    "res.statistic, res.pvalue"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f10349d8-9232-482b-92ed-74e59023dd08",
   "metadata": {},
   "source": [
    "#### One-sample Test\n",
    "\n",
    "The `wilcoxon` statistic does not depend on the specific values in `x` and `y`, only on the values of `x - y`. Instead of passing `x` and `y` into `wilcoxon` separately, `wilcoxon` can accept a single sample - the differences between paired observations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "c77d4c4c-efec-429b-9af1-b900308f05a5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "WilcoxonResult(statistic=25.0, pvalue=0.0390625)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = np.array([209, 200, 177, 169, 159, 169, 187])\n",
    "y = np.array([151, 168, 147, 164, 166, 163, 176])\n",
    "stats.wilcoxon(x - y, alternative='greater', method='exact')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d4d542a-a922-47ca-aebc-9328209a1da3",
   "metadata": {},
   "source": [
    "As described above, this is relatively common in paired-sample tests (e.g. see also `ttest_rel`), so `permutation_test` also supports single samples as input and forms the null distribution by permuting the *signs* of each observation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "9667f070-749c-45c2-8e51-1e284fda6a1c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0390625"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def statistic(d):\n",
    "    # return just the Mann-Whitney U statistic\n",
    "    return stats.wilcoxon(d, alternative='greater').statistic\n",
    "\n",
    "res = stats.permutation_test((x-y,), statistic, alternative='greater', permutation_type='samples')\n",
    "res.pvalue"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ecb0eb66-4c5c-4dd5-8403-a8c1a46f0ebd",
   "metadata": {},
   "source": [
    "#### Gotchas\n",
    "\n",
    "Suppose that we wish to perform the Wilcoxon test with a two-sided alternative. We might expect that `alternative='two-sided'` should be specified everywhere the `alternative` parameter is accepted."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "353b33cc-2bfc-43ea-b870-47fb600aa879",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.078125 0.15625\n"
     ]
    }
   ],
   "source": [
    "alternative = 'two-sided'\n",
    "res1 = stats.wilcoxon(x - y, alternative=alternative)\n",
    "\n",
    "def statistic(d):\n",
    "    return stats.wilcoxon(d, alternative=alternative).statistic\n",
    "\n",
    "res2 = stats.permutation_test((x-y,), statistic, alternative=alternative, permutation_type='samples')\n",
    "\n",
    "print(res1.pvalue, res2.pvalue)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cefd2711-9f37-46fb-97c1-073362fa7041",
   "metadata": {},
   "source": [
    "This guess was not correct; the `pvalue` returned by `permutation_test` is greater by a factor of two. Note that the documentation of `wilcoxon` states:\n",
    "> If `alternative` is “two-sided”, [the `statistic` is] the sum of the ranks of the differences above or below zero, whichever is smaller.\n",
    "\n",
    "The sign information that should be carried by the statistic (\"above or below zero\") is not preserved by `statistic`. This can be corrected by passing `alternative='greater'` or `alternative='less'` in the call to `wilcoxon` so that that the statistic is always:\n",
    "> the sum of the ranks of the differences above zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "2aa4d597-f85e-4e31-a5f3-d38b8dd6fef4",
   "metadata": {},
   "outputs": [],
   "source": [
    "def statistic(d):\n",
    "    # `alternative='less'` or `alternative='greater'` will work\n",
    "    return stats.wilcoxon(d, alternative='less').statistic\n",
    "\n",
    "res2 = stats.permutation_test((x-y,), statistic, alternative=alternative, permutation_type='samples')\n",
    "\n",
    "np.testing.assert_allclose(res2.pvalue, res1.pvalue, atol=1e-15)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "31159ba3-500b-468f-90cd-00647704cd1b",
   "metadata": {},
   "source": [
    "### Other Tests\n",
    "`permutation_test` with `permutation_type='samples'` is a versatile tool for comparing paired samples. Provided only data and a statistic, it can produce the null distribution and replicate the $p$-value of similar tests in SciPy, and it may be more accurate than these existing implementations, especially for small samples and when there are ties:\n",
    "\n",
    "- [`ttest_rel`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.ttest_rel.html)\n",
    "- [`wilcoxon`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.wilcoxon.html)\n",
    "- [`page_trend_test`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.page_trend_test.html)\n",
    "\n",
    "In addition, `permutation_test` with `permutation_type='samples'` can be used to perform tests not yet implemented in SciPy.\n",
    "\n",
    "However, there are yet other types of permutation tests that assume neither that samples are independent nor paired. We conclude the study of `permutation_test` with [Correlated-Sample Tests](https://nbviewer.org/github/scipy/scipy-cookbook/blob/main/ipython/ResamplingAndMonteCarloMethods/resampling_tutorial_2c.ipynb)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}