example- Randomisation test on continuously variable data.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "DC_stats version:  0.3.0\n"
     ]
    }
   ],
   "source": [
    "import dcstats\n",
    "print(\"DC_stats version: \", dcstats.__version__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "RANTEST performs a randomisation test to compare two independent samples.  According to the null hypothesis of    no-difference, each outcome would have been the same regardless of which group the individual happened to be allocated to.  Therefore all N = n1 + n2 observations are pooled and, as in the actual experiment, divided at random into groups of size n1 and n2.  The fraction of randomisations that gives rise to a difference between the groups at least as large as that observed gives the P value.\n",
    "\n",
    "In the binomial case, in which the measurement is the fraction of 'successes' in each sample (say r1 out of n1, and r2 out of n2) a 'success' is given a score of 1, 'failure' scores 0.\n",
    "Randomisation test on binomial data could be done using as criterion: \n",
    "(1) number of successes in set 1 (r1) or \n",
    "(2) difference between the p=r/n values.\n",
    "Both criteria give the same one-tail P value. Use of r1 as criterion is the direct Monte Carlo equivalent of summing the the exact Fisher test probabilities for the observed 2x2 table with those for all tables that depart further from the null hypothesis in the observed direction. A 2-tail probablilty can be found by doubling the one-tail value, at least if the displayed distribution is symmetrical. Use of (p1-p2) as criterion gives both one and two-tail     probabilities directly by seeing how many random allocations of the observations to groups of size n1 and n2 produce and absolute value of (p1-p2) at least as big as that observed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from dcstats.rantest import RantestContinuous\n",
    "from dcstats.basic_stats import TTestContinuous\n",
    "from dcstats.Hedges import Hedges_d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Samples from treatment T1 and T2\n",
    "T1 = [100, 108, 119, 127, 132, 135, 136]#, 164]\n",
    "T2 = [122, 130, 138, 142, 152, 154, 176]\n",
    "are_paired = False"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate basic statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "ttc = TTestContinuous(T1, T2, are_paired)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n \t\t 7     \t7\n",
      "Mean \t\t 122.428571    \t  144.857143\n",
      "SD \t\t 14.010200     \t  17.808505\n",
      "SDM \t\t 5.295358     \t  6.730982\n",
      "95% confidence intervals:\n",
      "lower \t\t 112.049670    \t  131.664418\n",
      "upper \t\t 132.807473    \t  158.049868\n",
      "\n",
      " Mean difference (dbar) = \t -22.428571\n",
      "  s(d) = \t 8.100558 \t s(dbar) = \t 3.061723\n",
      "\n",
      " Two-sample unpaired Students t-test:\n",
      " t = \t -2.618850\n",
      " two tail P = \t 2.243E-02\n",
      "\n",
      " statistical power = 0.745\n"
     ]
    }
   ],
   "source": [
    "print(ttc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate Hedges statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "Effect size\n",
      "  Hedges unbiased d = \t 1.310483\n",
      "  approximate 95% confidence intervals \n",
      "  upper 95% CI =\t 2.394100\n",
      "  lower 95% CI =\t 0.639740\n"
     ]
    }
   ],
   "source": [
    "#calculation of hedges d and approximate 95% confidence intervals\n",
    "#not tested against known values yet AP 170518\n",
    "hedges_calculation = Hedges_d(T1, T2)\n",
    "hedges_calculation.hedges_d_unbiased()\n",
    "#paired needed for degrees of freedom\n",
    "lowerCI, upperCI = hedges_calculation.bootstrap_CI(5000)\n",
    "#option to have bootstrap calculated CIs should go here\n",
    "print(hedges_calculation)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run randomisation tests"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "nran = 50000\n",
    "rnt = RantestContinuous(T1, T2, are_paired)    \n",
    "rnt.run_rantest(nran)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Rantest:  50000 randomisations\n",
      "Two-tailed P = 2.078e-02\t(greater than or equal in absolute value to observed)\n",
      "Pequal = 2.540e-03\t(equal in absolute value to observed)\n"
     ]
    }
   ],
   "source": [
    "print(rnt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Red line: observed difference\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAADTCAYAAAAiRcCrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAGXVJREFUeJzt3X2cHWV99/HP14RHeQiBhcaEuOFuRJEXSlgBFS0KYgAl2IKNtSVQNNUbEW29b4JYQ0V6g1qp9FYBSyRQCwZQiKKGgARoy0MSAuEhpImBQprchBoIII+B3/3H/LYMm7O7Zzd75pyz+b5fr/PamWuumfnN5Owv114zc40iAjMzq8Ybmh2AmdmWxEnXzKxCTrpmZhVy0jUzq5CTrplZhZx0zcwq5KRrZlYhJ10zswo56ZqZVWhkswOo2m677RadnZ3NDuM1y5cXP/feu7lxmLWDFv59Wbx48X9FREd/9ba4pNvZ2cmiRYuaHcZrDj20+LlgQTOjMGsPLfz7Iuk/6qnn7gUzswo56ZqZVahhSVfSLEnrJN1fKhstab6kFflzlyyXpAskrZS0VNKk0jrTsv4KSdNK5QdIui/XuUCSGnUsZmZDpZEt3UuByT3KZgA3RcRE4KacBzgSmJif6cD3oUjSwEzgIOBAYGZ3os4600vr9dyXmVnLaVjSjYhbgfU9iqcAs3N6NnBsqfyyKNwBjJI0BvgwMD8i1kfEk8B8YHIu2ykibo9iQODLStsyM2tZVd+9sEdErAWIiLWSds/yscBjpXqrs6yv8tU1ymuSNJ2iVcz48eM38xBsS9E54/qG7+ORc49u+D6stbTKhbRa/bExiPKaIuLiiOiKiK6Ojn5vozMza5iqW7qPSxqTrdwxwLosXw3sWao3DliT5Yf2KF+Q5eNq1LctRBWtULNGqLqlOxfovgNhGnBdqfyEvIvhYGBDdkPMA46QtEteQDsCmJfLnpF0cN61cEJpW2ZmLathLV1JV1C0UneTtJriLoRzgTmSTgYeBY7P6r8AjgJWAs8BJwFExHpJZwMLs97XIqL74txnKe6Q2A74ZX7MzFpaw5JuRHyil0WH1agbwCm9bGcWMKtG+SJg382J0cysaq1yIc3MbIvgpGtmViEnXTOzCjnpmplVyEnXzKxCTrpmZhVy0jUzq5CTrplZhZx0zcwq5KRrZlYhJ10zswo56ZqZVchJ18ysQk66ZmYVctI1M6uQk66ZWYWcdM3MKuSka2ZWISddM7MKOemamVXISdfMrEJOumZmFWpK0pX0RUkPSLpf0hWStpU0QdKdklZI+rGkrbPuNjm/Mpd3lrZzRpYvl/ThZhyLmdlAVJ50JY0FPg90RcS+wAhgKnAecH5ETASeBE7OVU4GnoyI3wfOz3pI2ifXezswGfiepBFVHouZ2UA1q3thJLCdpJHA9sBa4IPA1bl8NnBsTk/JeXL5YZKU5VdGxIsR8TCwEjiwovjNzAal8qQbEf8JfAt4lCLZbgAWA09FxMasthoYm9Njgcdy3Y1Zf9dyeY11XkfSdEmLJC164oknhvaAzMwGoBndC7tQtFInAG8C3ggcWaNqdK/Sy7LeyjctjLg4Iroioqujo2PgQZuZDZFmdC8cDjwcEU9ExMvAT4D3AKOyuwFgHLAmp1cDewLk8p2B9eXyGuuYmbWkZiTdR4GDJW2ffbOHAQ8CNwPHZZ1pwHU5PTfnyeW/jojI8ql5d8MEYCJwV0XHYGY2KCP7rzK0IuJOSVcDdwMbgSXAxcD1wJWSvp5ll+QqlwCXS1pJ0cKdmtt5QNIcioS9ETglIl6p9GDMzAao8qQLEBEzgZk9ildR4+6DiHgBOL6X7ZwDnDPkAdqgdc64vtkhmLU0P5FmZlYhJ10zswo1pXvBzApVdcc8cu7RlezH+ueWrplZhZx0zcwq5KRrZlYhJ10zswo56ZqZVchJ18ysQk66ZmYVctI1M6uQk66ZWYWcdM3MKuSka2ZWobqSrqR9Gx2ImdmWoN6W7oWS7pL0PyWNamhEZmbDWF1JNyIOAT5J8U6yRZL+WdKHGhqZmdkwVHefbkSsAL4CnA78AXCBpIck/WGjgjMzG27q7dPdT9L5wDLgg8BHI+JtOX1+A+MzMxtW6h3E/P8CPwC+HBHPdxdGxBpJX2lIZGZmw1C9Sfco4Pnut+1KegOwbUQ8FxGXNyw6M7Nhpt4+3RuB7Urz22eZmZkNQL1Jd9uIeLZ7Jqe3H+xOJY2SdHVeiFsm6d2SRkuaL2lF/twl60rSBZJWSloqaVJpO9Oy/gpJ0wYbj5lZVepNur/rkewOAJ7vo35/vgP8KiLeCryD4gLdDOCmiJgI3JTzAEcCE/MzHfh+xjAamAkcBBwIzOxO1GZmrarePt0vAFdJWpPzY4A/HswOJe0EvB84ESAiXgJekjQFODSrzQYWUNyeNgW4LCICuCNbyWOy7vyIWJ/bnQ9MBq4YTFxmZlWoK+lGxEJJbwX2BgQ8FBEvD3KfewFPAD+U9A5gMXAasEdErM39rZW0e9YfCzxWWn91lvVWvglJ0ylayYwfP36QYZuZbb6BDHjzLmA/YH/gE5JOGOQ+RwKTgO9HxP7A73itK6EW1SiLPso3LYy4OCK6IqKro6NjoPGamQ2Zeh+OuBz4FnAIRfJ9F9A1yH2uBlZHxJ05fzVFEn48uw3In+tK9fcsrT8OWNNHuZlZy6q3T7cL2Cf7VTdLRPw/SY9J2jsilgOHAQ/mZxpwbv68LleZC3xO0pUUF802ZPfDPOBvSxfPjgDO2Nz4zMwaqd6kez/we8DaIdrvqcCPJG0NrAJOomh1z5F0MvAocHzW/QXFwxkrgeeyLhGxXtLZwMKs97Xui2pmZq2q3qS7G/CgpLuAF7sLI+KYwew0Iu6hdvfEYTXqBnBKL9uZBcwaTAxmZs1Qb9I9q5FBmJltKeq9ZewWSW8GJkbEjZK2B0Y0NjQzs+Gn3rsXPk1xl8FFWTQWuLZRQZmZDVf13qd7CvBe4Gn47wHNd+9zDTMz20S9SffFfFwXAEkj6eVBBDMz6129SfcWSV8Gtst3o10F/KxxYZmZDU/1Jt0ZFOMl3Af8BcW9s35jhJnZANV798KrFK/r+UFjwzEzG97qSrqSHqZGH25E7DXkEZmZDWMDGXuh27YUj+iOHvpwzMyGt7r6dCPit6XPf0bE31O8ft3MzAag3u6FSaXZN1C0fHdsSERmZsNYvd0Lf1ea3gg8Anx8yKOxhuqccX2zQ7AmqeLf/pFzj274PoaDeu9e+ECjAzEz2xLU273wl30tj4hvD004ZmbD20DuXngXxVscAD4K3MrrXwxpZmb9GMgg5pMi4hkASWcBV0XEpxoVmJnZcFTvY8DjgZdK8y8BnUMejZnZMFdvS/dy4C5JP6V4Mu1jwGUNi8rMbJiq9+6FcyT9EnhfFp0UEUsaF5aZ2fBUb/cCwPbA0xHxHWC1pAkNisnMbNiq93U9M4HTgTOyaCvgnxoVlJnZcFVvS/djwDHA7wAiYg1+DNjMbMDqTbovRUSQwztKeuPm7ljSCElLJP085ydIulPSCkk/lrR1lm+T8ytzeWdpG2dk+XJJH97cmMzMGq3epDtH0kXAqHwz8I1s/oDmpwHLSvPnAedHxETgSeDkLD8ZeDIifh84P+shaR9gKvB2YDLwPUl+LbyZtbR6h3b8FsUr2K8B9ga+GhH/MNidShoHHA38Y86LYqjIq7PKbODYnJ6S8+Tyw7L+FODKiHgxIh4GVgIHDjYmM7Mq9HvLWLYe50XE4cD8Idrv3wP/m9f6hXcFnoqIjTm/Ghib02PJx40jYqOkDVl/LHBHaZvldXoew3RgOsD48eOH6BDMzAau35ZuRLwCPCdp56HYoaSPAOsiYnG5uNau+1nW1zqvL4y4OCK6IqKro6NjQPGamQ2lep9IewG4T9J88g4GgIj4/CD2+V7gGElHUbz6ZyeKlu8oSSOztTsOWJP1VwN7UtwbPBLYGVhfKu9WXsfMrCXVeyHteuCvKUYWW1z6DFhEnBER4yKik+JC2K8j4pPAzcBxWW0acF1Oz815cvmv806KucDUvLthAjARuGswMZmZVaXPlq6k8RHxaETM7qveEDkduFLS14ElwCVZfglwuaSVFC3cqQAR8YCkOcCDFG+zOCW7QszMWlZ/3QvXApMAJF0TEX80lDuPiAXAgpxeRY27DyLiBYq3D9da/xzgnKGMycyskfrrXihfrNqrkYGYmW0J+ku60cu0mZkNQn/dC++Q9DRFi3e7nCbnIyJ2amh0ZmbDTJ9JNyL8WK2Z2RAayHi6Zma2mZx0zcwq5KRrZlYhJ10zswo56ZqZVchJ18ysQk66ZmYVctI1M6uQk66ZWYWcdM3MKuSka2ZWISddM7MKOemamVXISdfMrEJOumZmFXLSNTOrkJOumVmFnHTNzCpUedKVtKekmyUtk/SApNOyfLSk+ZJW5M9dslySLpC0UtJSSZNK25qW9VdImlb1sZiZDVR/L6ZshI3AX0XE3ZJ2BBZLmg+cCNwUEedKmgHMAE4HjgQm5ucg4PvAQZJGAzOBLoo3FS+WNDcinqz8iMyMzhnXN3wfV676LQfvtWvD99NIlSfdiFgLrM3pZyQtA8YCU4BDs9psYAFF0p0CXBYRAdwhaZSkMVl3fkSsB8jEPRm4orKDGUJVfGHNrPma2qcrqRPYH7gT2CMTcndi3j2rjQUeK622Ost6K6+1n+mSFkla9MQTTwzlIZiZDUjTkq6kHYBrgC9ExNN9Va1RFn2Ub1oYcXFEdEVEV0dHx8CDNTMbIk1JupK2oki4P4qIn2Tx49ltQP5cl+WrgT1Lq48D1vRRbmbWsppx94KAS4BlEfHt0qK5QPcdCNOA60rlJ+RdDAcDG7L7YR5whKRd8k6HI7LMzKxlNePuhfcCfwbcJ+meLPsycC4wR9LJwKPA8bnsF8BRwErgOeAkgIhYL+lsYGHW+1r3RTUzs1bVjLsX/oXa/bEAh9WoH8ApvWxrFjBr6KIzM2ssP5FmZlYhJ10zswo56ZqZVchJ18ysQk66ZmYVctI1M6uQk66ZWYWcdM3MKuSka2ZWISddM7MKOemamVXISdfMrEJOumZmFXLSNTOrkJOumVmFnHTNzCrUjDdHmJkN2h2rfsvUGdc3dB+PnHt0w7btpNuPzgb/41656rfFxMEN3Y2ZtQh3L5iZVchJ18ysQk66ZmYVctI1M6tQ2yddSZMlLZe0UtKMZsdjZtaXtk66kkYA3wWOBPYBPiFpn+ZGZWbWu7ZOusCBwMqIWBURLwFXAlOaHJOZWa8UEc2OYdAkHQdMjohP5fyfAQdFxOd61JsOTM/ZvYHlOb0b8F8VhTtU2jFmaM+42zFmaM+42zFmeH3cb46Ijv5WaPeHI1SjbJP/RSLiYuDiTVaWFkVEVyMCa5R2jBnaM+52jBnaM+52jBkGF3e7dy+sBvYszY8D1jQpFjOzfrV70l0ITJQ0QdLWwFRgbpNjMjPrVVt3L0TERkmfA+YBI4BZEfHAADaxSZdDG2jHmKE9427HmKE9427HmGEQcbf1hTQzs3bT7t0LZmZtxUnXzKxCW1zSlXS2pKWS7pF0g6Q3ZbkkXZCPEy+VNKnZsZZJ+qakhzK2n0oaVVp2Rsa9XNKHmxlnmaTjJT0g6VVJXT2WtWTM3drh8XJJsyStk3R/qWy0pPmSVuTPXZoZYy2S9pR0s6Rl+f04LctbNnZJ20q6S9K9GfPfZPkESXdmzD/OC/p9i4gt6gPsVJr+PHBhTh8F/JLi3t+DgTubHWuPuI8ARub0ecB5Ob0PcC+wDTAB+A0wotnxZmxvo3gYZQHQVSpv2ZgzvhEZ017A1hnrPs2Oq0ac7wcmAfeXyr4BzMjpGd3fk1b6AGOASTm9I/Dv+Z1o2dgzL+yQ01sBd2aemANMzfILgc/2t60trqUbEU+XZt/Iaw9TTAEui8IdwChJYyoPsBcRcUNEbMzZOyjuSYYi7isj4sWIeBhYSfF4dNNFxLKIWF5jUcvGnNri8fKIuBVY36N4CjA7p2cDx1YaVB0iYm1E3J3TzwDLgLG0cOyZF57N2a3yE8AHgauzvK6Yt7ikCyDpHEmPAZ8EvprFY4HHStVWZ1kr+nOKVjm0V9zdWj3mVo+vL3tExFookhuwe5Pj6ZOkTmB/ipZjS8cuaYSke4B1wHyKv4aeKjWG6vqeDMukK+lGSffX+EwBiIgzI2JP4EdA9zgNdT1S3Ej9xZ11zgQ2UsQOTY67nphrrVajrJXuXWz1+IYFSTsA1wBf6PEXaEuKiFci4p0Uf2UeSNF9tkm1/rbT1g9H9CYiDq+z6j8D1wMzaYFHivuLW9I04CPAYZGdSDQ57gGc67Kmn+t+tHp8fXlc0piIWJvdY+uaHVAtkraiSLg/ioifZHFbxB4RT0laQNGnO0rSyGzt1vU9GZYt3b5ImliaPQZ4KKfnAifkXQwHAxu6/9RpBZImA6cDx0TEc6VFc4GpkraRNAGYCNzVjBgHoNVjbufHy+cC03J6GnBdE2OpSZKAS4BlEfHt0qKWjV1SR/cdQ5K2Aw6n6Iu+GTguq9UXc7OvCjbhKuQ1wP3AUuBnwNjS1cnvUvTT3EfpansrfCguNj0G3JOfC0vLzsy4lwNHNjvWUlwfo2g1vgg8Dsxr9ZhL8R1FcVX9N8CZzY6nlxivANYCL+d5PhnYFbgJWJE/Rzc7zhpxH0LxZ/jS0vf5qFaOHdgPWJIx3w98Ncv3omgwrASuArbpb1t+DNjMrEJbXPeCmVkzOemamVXISdfMrEJOumZmFXLSNTOrkJNuG5P0So6Wdr+kn5VHHtvM7XaWR65qFEnP5s83Sbq6VH5Fjqb2RUlvzWNcIul/NDqmwZB0qYo3U9dbv1PSnzQyJmtdTrrt7fmIeGdE7Esx8MkpzQ5oMCJiTUQcByDp94D3RMR+EXE+xQAi10XE/hHxm3q2J2lEA8MdCp2Ak+4Wykl3+LidHGxD0g6SbpJ0t6T7usdByBbWMkk/yDFBb8ina5B0QI4Vejul5J3jiP4wt7NE0gey/ERJ12YL+2FJn5P0l1nnDkmjewaYT3jdLmmhpLNL5eWW9Q3A7tm6nQl8AfiUpJuz7p/muKb3SLqoO8FKelbS1yTdCbw7j+cWSYslzcvHSpG0QNJ5uY1/l/S+LB8h6Vt5nEslnVo6L5tsp4bDJd2W2/xIaZvfzONdKukvsu65wPvyGL4o6ReS9st1lkj6ak6fLelTOf2/Stv5m9K56+t8nJP/pndI2qPGv8dZkmbn9+ARSX8o6Rt5Dn6l4lHdXs+BpE9nTPdKukbS9ll+qYqxqf9N0qruvwIkjZF0q1776+x9vZzL4a3ZT3r4s1lPyTybP0dQPA0zOedHkuMGA7tRPC0jihbWRuCduWwO8Kc5vRT4g5z+JjlGK/BXwA9z+q3Ao8C2wIm53R2BDmAD8Jmsdz7FICY9450LnJDTp5Ti7yzt77+nc/4s4Es5/TaKpwi3yvnvlbYXwMdzeivg34COnP9jipeWQjG279/l9FHAjTn9WYqnFbvHLB7d13Z6HNelwK8oGjETKZ4O2xaYDnwl62wDLKIYP/hQ4Oel9Wfk+diJ4hHkeVl+M8V4xEdQvABRuY+fU4yl29/5+GhOf6M7jh5xnwX8Sx7nO4DnyKcDgZ9S/JXR17nctbStrwOnls7HVRnrPhTDZELxXTqz9J3dsdm/Q834DMsBb7Yg26kYaq4TWEwx3BwUv5x/K+n9wKsULeDuls7DEXFPTi8GOiXtDIyKiFuy/HLgyJw+BPgHgIh4SNJ/AG/JZTdHMR7qM5I2UCQAKB6j3q9GvO8F/qi0j/MGeLyHAQcACyUBbMdrg6K8QpE0oUhU+wLzs94Iisdlu3UPsLKY4txB8Sz9hZHD9EXEekn79rOdsjkR8SqwQtIqiv+gjgD202v9vTtTJOWXeqx7G8WA+g9TDMD0oWw1dkbEckmfzm0tyfo75Hb26+N8vESRnLuP80O9xP3LiHhZ0n15fL/K8vvy3PR1LveV9HVgVMY0r7Tda/N8PFhqZS8EZmUL+trS93CL4qTb3p6PiHdm0vw5RWvpAopxgjuAA/IX6hGKlhcU4yB0e4XiF1X0PiRdrWEOu5W39Wpp/lV6/25tznPnAmZHxBk1lr0QEa+U6j0QEe/uZTvdcb7Ca3HWOgf9baes57qR658aEeVkhKRDe9RdCHQBqyj+49wN+DRFsuyO4/9ExEU9tnMqvZ+PlyOblLz+OHt6ESAiXpVUXqf737Cvc3ApcGxE3CvpRIoW/Ou2W4qfiLg1GwJHA5dL+mZEXNZLXMOW+3SHgYjYQNFS+lK2InYG1mXC/QDw5n7WfwrYIOmQLPpkafGt3fOS3gKMpxikZjD+lWLErp77qNdNwHGSds94RkuqdWzLgQ5J7856W0l6ez/bvgH4jKSR3dse4HaOl/QGFXdY7JXrzgM+W+obfYukNwLPUHTLABDF2ykeAz5O8VaQ24Av5U9yO3+uYvxZJI3Nc1Dv+dgcfZ2DHYG1eXz9/ntmbOsi4gcUo4y11HsIq+KkO0xExBKKd3lNpRjgvEvSIopfhof6WjedBHxXxYW050vl3wNG5J+fPwZOjIgXa22gDqcBp0haSPEfw4BExIPAV4AbJC2laBVucmErk9hxwHmS7qUYxeo9/Wz+Hyn6q5fmOn8ywO0sB26heKPHZyLihdzmg8DdKi4UXkTRelwKbMwLUF/M9W8DHo9i2M7bKMZmvS2P5waKsZ9vz3+Hqyn6Q+s6H5ujn3Pw1xRvfJhPfd+xQ4F7JC2h6Gb6zlDG2i48ypiZWYXc0jUzq5CTrplZhZx0zcwq5KRrZlYhJ10zswo56ZqZVchJ18ysQv8f+0YWpmF9HIAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 360x216 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax  = plt.subplots(1,1, figsize=(5,3))\n",
    "ax.hist(rnt.randiff, 10)\n",
    "ax.axvline(x=rnt.dbar, color='r')\n",
    "ax.axvline(x=-rnt.dbar, color='r')\n",
    "ax.set_ylabel(\"Frequency\")\n",
    "ax.set_xlabel('Random difference between means')\n",
    "print('Red line: observed difference')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "What id measurements were paired?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "ttc = TTestContinuous(T1, T2, True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n \t\t 7     \t7\n",
      "Mean \t\t 122.428571    \t  144.857143\n",
      "SD \t\t 14.010200     \t  17.808505\n",
      "SDM \t\t 5.295358     \t  6.730982\n",
      "95% confidence intervals:\n",
      "lower \t\t 112.049670    \t  131.664418\n",
      "upper \t\t 132.807473    \t  158.049868\n",
      "\n",
      " Mean difference (dbar) = \t -22.428571\n",
      "  s(d) = \t 8.100558 \t s(dbar) = \t 3.061723\n",
      "\n",
      " Paired Students t-test:\n",
      "  t(6)= \t dbar / s(dbar) \t = \t -7.325473\n",
      "  two tail P =\t 0.000331\n"
     ]
    }
   ],
   "source": [
    "print(ttc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "nran = 50000\n",
    "rnt = RantestContinuous(T1, T2, True)    \n",
    "rnt.run_rantest(nran)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Rantest:  50000 randomisations\n",
      "Two-tailed P = 1.604e-02\t(greater than or equal in absolute value to observed)\n",
      "Pequal = 1.604e-02\t(equal in absolute value to observed)\n"
     ]
    }
   ],
   "source": [
    "print(rnt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Red line: observed difference\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVcAAADTCAYAAAA1Z1BiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAFVVJREFUeJzt3X+0VWWdx/H3R0gBTRFBcwC7WKSZU/64qaVZpqOoJdbIRMuZ0FwxraGpJluK1pJ+OUvHynJKJ8wf5DgqWilmZYi/+gUKgiKiwRIS0hRFIX+kod/5Yz9XNnedc+++P5577zl8XmvddfZ+9vPs/X0O937Z5zl7P1sRgZmZ9a5t+jsAM7Nm5ORqZpaBk6uZWQZOrmZmGTi5mpll4ORqZpaBk6uZWQZOrmZmGTi5mpllMLi/A8hh5MiR0dLSUr3BI48Ur3vtlSUeMxuAuvF3v2jRoqcjYlSVuk2ZXFtaWli4cGH1Bh/4QPF65505wjGzgagbf/eS/li1rocFzMwycHI1M8vAydXMLAMnVzOzDJxczcwyaMqrBSy/lum3ZD/G6vOOz34Ms1x85mpmloHPXM36QF+c6YPP9gcSn7mamWXg5GpmloGTq5lZBk6uZmYZOLmamWXg5GpmloGTq5lZBk6uZmYZOLmamWXg5GpmloGTq5lZBk6uZmYZZE2ukv5D0jJJD0q6RtIQSeMkLZC0QtJ1krZNdbdL6yvT9pbSfs5K5Y9IOiZnzGZmvSFbcpU0Gvgs0BoR+wKDgMnA+cCFETEeeBY4LTU5DXg2It4KXJjqIWmf1O4dwATgYkmDcsVtZtYbcg8LDAaGShoMDAOeAD4I3JC2zwJOTMsT0zpp+5GSlMqvjYiXI2IVsBI4KHPcZmY9ki25RsSfgG8Cj1Ek1Q3AIuC5iNiUqq0FRqfl0cCa1HZTqr9LubxGm9dJmippoaSF69at6/0OmZl1Qc5hgZ0pzjrHAX8HbA8cW6NqtDWps61e+ZYFETMjojUiWkeNGtW9oM3MeknOYYGjgFURsS4i/gb8BHgvMDwNEwCMAR5Py2uBsQBp+07A+nJ5jTZmZgNSzuT6GHCIpGFp7PRI4CHgDuCkVGcKcFNanpPWSdtvj4hI5ZPT1QTjgPHAPRnjNjPrsWzP0IqIBZJuAO4DNgGLgZnALcC1kr6Ryi5LTS4DrpK0kuKMdXLazzJJsykS8yZgWkS8mituM7PekPUBhRExA5jRrvhRanzbHxF/BSbV2c+5wLm9HqCZWSa+Q8vMLAMnVzOzDJxczcwycHI1M8vAydXMLAMnVzOzDJxczcwycHI1M8vAydXMLAMnVzOzDJxczcwyyDq3gFkjaJl+S3+HYE3IZ65mZhk4uZqZZeBhARuw/HHdGpnPXM3MMnByNTPLwMnVzCwDJ1czswycXM3MMnByNTPLwMnVzCwDX+fahHx9qFn/c3I1ayJ98R/r6vOOz36MZuBhATOzDLImV0nDJd0g6WFJyyW9R9IISXMlrUivO6e6knSRpJWSHpB0QGk/U1L9FZKm5IzZzKw35D5z/S7wy4jYG3gXsByYDsyLiPHAvLQOcCwwPv1MBS4BkDQCmAEcDBwEzGhLyGZmA1Wl5Cpp367uWNKOwOHAZQAR8UpEPAdMBGalarOAE9PyROBHUZgPDJe0O3AMMDci1kfEs8BcYEJX4zEz60tVz1z/R9I9kv5N0vCKbfYE1gFXSFos6YeStgd2i4gnANLrrqn+aGBNqf3aVFavfAuSpkpaKGnhunXrKoZoZpZHpeQaEYcBJwNjgYWS/k/SP3TSbDBwAHBJROwPvMDmIYBaVOvQHZS3j3FmRLRGROuoUaM6Cc3MLK/KY64RsQL4MnAm8H7govRF1UfrNFkLrI2IBWn9Bopk+2T6uE96fapUf2yp/Rjg8Q7KzcwGrKpjru+UdCHFF1IfBD4cEW9PyxfWahMRfwbWSNorFR0JPATMAdq+8Z8C3JSW5wCfSFcNHAJsSMMGtwJHS9o5fZF1dCozMxuwqt5E8D3gUuDsiHiprTAiHpf05Q7a/TtwtaRtgUeBUykS+mxJpwGPAZNS3Z8DxwErgRdTXSJivaSvA/emel+LiPUV4zYz6xdVk+txwEsR8SqApG2AIRHxYkRcVa9RRCwBWmtsOrJG3QCm1dnP5cDlFWM1M+t3VcdcbwOGltaHpTIzM6uhanIdEhHPt62k5WF5QjIza3xVk+sL7W5HPRB4qYP6ZmZbtapjrp8HrpfUdgnU7sDH8oRkZtb4KiXXiLhX0t7AXhQX9T8cEX/LGpmZWQPrynyu7wZaUpv9JRERP8oSlZlZg6uUXCVdBbwFWAK8mooDcHI1M6uh6plrK7BPuhbVzMw6UfVqgQeBN+UMxMysmVQ9cx0JPCTpHuDltsKIOCFLVGZmDa5qcv1KziDMzJpN1Uux7pL0ZmB8RNwmaRgwKG9oZmaNq+qUg5+imI/1B6loNHBjrqDMzBpd1S+0pgGHAhvh9Ymzd+2whZnZVqxqcn05Il5pW5E0mBqPWjEzs0LV5HqXpLOBoenZWdcDN+cLy8yssVVNrtMpnuS6FPhXiqcGdPQEAjOzrVrVqwVeo3jMy6V5wzEzaw5V5xZYRe3HWe/Z6xGZmTWBrswt0GYIxUMFR/R+OGZmzaHSmGtEPFP6+VNEfIfisdpmZlZD1WGBA0qr21Ccyb4xS0RmZk2g6rDAt0rLm4DVwD/1ejRmZk2i6tUCR+QOxMysmVQdFvhCR9sj4tu9E46ZWXPoytUC7wbmpPUPA3cDa3IEZWbW6KreoTUSOCAiTo+I04EDgTER8dWI+GpHDSUNkrRY0s/S+jhJCyStkHSdpG1T+XZpfWXa3lLax1mp/BFJx3Sno2Zmfalqct0DeKW0/grFk2Cr+BywvLR+PnBhRIwHngVOS+WnAc9GxFuBC1M9JO0DTAbeAUwALpbkuWTNbECrmlyvAu6R9BVJM4AFVHjyq6QxwPHAD9O6KK6PvSFVmQWcmJYnpnXS9iNT/YnAtRHxckSsAlYCB1WM28ysX1S9WuBcSb8A3peKTo2IxRWafgc4g83XxO4CPBcRm9L6WoqJt0mva9LxNknakOqPBuaX9llu8zpJU4GpAHvssUeVbpmZZVP1zBVgGLAxIr4LrJU0rqPKkj4EPBURi8rFNapGJ9s6arO5IGJmRLRGROuoUaM6Cs3MLLuql2LNoLhiYC/gCuANwP9SPJ2gnkOBEyQdRzEfwY4UZ7LDJQ1OZ69jgMdT/bXAWIrEPRjYCVhfKm9TbmNmNiBVPXP9CHAC8AJARDxOJ7e/RsRZETEmIloovpC6PSJOBu4ATkrVpgA3peU5aZ20/faIiFQ+OV1NMA4YD9xTMW4zs35R9TrXVyIiJAWApO17cMwzgWslfQNYDFyWyi8DrpK0kuKMdTJARCyTNBt4iOLW22kR8WoPjm9mll3V5Dpb0g8oPtJ/CvgkXZg4OyLuBO5My49S49v+iPgrxVSGtdqfC5xb9XhmZv2t6tUC30zPztpIMe56TkTMzRqZmVkD6zS5pgv2b42IowAnVDOzCjr9QiuNb74oaac+iMfMrClUHXP9K7BU0lzSFQMAEfHZLFGZmTW4qsn1lvRjZmYVdJhcJe0REY9FxKyO6pmZ2ZY6G3O9sW1B0o8zx2Jm1jQ6S67l+/r3zBmImVkz6Sy5Rp1lMzPrQGdfaL1L0kaKM9ihaZm0HhGxY9bozMwaVIfJNSI847+ZWTd0ZT5XMzOryMnVzCwDJ1czswycXM3MMnByNTPLwMnVzCwDJ1czswycXM3MMnByNTPLwMnVzCwDJ1czswycXM3MMnByNTPLwMnVzCyDbMlV0lhJd0haLmmZpM+l8hGS5kpakV53TuWSdJGklZIekHRAaV9TUv0VkqbkitnMrLfkPHPdBJweEW8HDgGmSdoHmA7Mi4jxwLy0DnAsMD79TAUugSIZAzOAg4GDgBltCdnMbKDKllwj4omIuC8t/wVYDowGJgJtT5OdBZyYlicCP4rCfGC4pN2BY4C5EbE+Ip4F5gITcsVtZtYb+mTMVVILsD+wANgtIp6AIgEDu6Zqo4E1pWZrU1m98vbHmCppoaSF69at6+0umJl1SWfP0OoxSTsAPwY+HxEbJdWtWqMsOijfsiBiJjAToLW1dUA+TLFl+i39HYJZj/XF7/Hq847Pfozcsp65SnoDRWK9OiJ+koqfTB/3Sa9PpfK1wNhS8zHA4x2Um5kNWDmvFhBwGbA8Ir5d2jQHaPvGfwpwU6n8E+mqgUOADWnY4FbgaEk7py+yjk5lZmYDVs5hgUOBfwGWSlqSys4GzgNmSzoNeAyYlLb9HDgOWAm8CJwKEBHrJX0duDfV+1pErM8Yt5lZj2VLrhHxG2qPlwIcWaN+ANPq7Oty4PLei87MLC/foWVmloGTq5lZBk6uZmYZOLmamWXg5GpmloGTq5lZBk6uZmYZOLmamWXg5GpmloGTq5lZBk6uZmYZOLmamWXg5GpmloGTq5lZBk6uZmYZZH+GlplZV/XFc7quffQZDtlzl2z7d3JN5j/6DJP9AEEz6yUeFjAzy8DJ1cwsAydXM7MMnFzNzDJwcjUzy8DJ1cwsAydXM7MMnFzNzDJwcjUzy6BhkqukCZIekbRS0vT+jsfMrCMNkVwlDQK+DxwL7AN8XNI+/RuVmVl9DZFcgYOAlRHxaES8AlwLTOznmMzM6mqUiVtGA2tK62uBg8sVJE0FpqbV5yU90uWjnP+h7sQ2Eni6Ow0bRLP3D5q/j83eP+hGH98DRVaRutLszVUrNkpyrdX72GIlYiYws2/C2UzSwoho7evj9pVm7x80fx+bvX8wMPvYKMMCa4GxpfUxwOP9FIuZWacaJbneC4yXNE7StsBkYE4/x2RmVldDDAtExCZJnwFuBQYBl0fEsn4Oq02fD0X0sWbvHzR/H5u9fzAA+6iI6LyWmZl1SaMMC5iZNRQnVzOzDJxcu0HSBZIelvSApJ9KGl7adla6RfcRScf0Z5w9IWmSpGWSXpPU2m5bs/Sx6W6plnS5pKckPVgqGyFprqQV6XXn/oyxJySNlXSHpOXp9/NzqXzA9dHJtXvmAvtGxDuBPwBnAaRbcicD7wAmABenW3cb0YPAR4G7y4XN0scmvqX6Sop/l7LpwLyIGA/MS+uNahNwekS8HTgEmJb+3QZcH51cuyEifhURm9LqfIrrbqG4JffaiHg5IlYBKylu3W04EbE8Imrd5dYsfWzKW6oj4m5gfbviicCstDwLOLFPg+pFEfFERNyXlv8CLKe4g3PA9dHJtec+CfwiLde6TXd0n0eUV7P0sVn6UcVuEfEEFMkJ2LWf4+kVklqA/YEFDMA+NsR1rv1B0m3Am2ps+lJE3JTqfIniY8rVbc1q1B+w17pV6WOtZjXKBmwfO9As/dgqSdoB+DHw+YjYqK7ND9AnnFzriIijOtouaQrwIeDI2HyxcEPdpttZH+toqD52oFn6UcWTknaPiCck7Q481d8B9YSkN1Ak1qsj4iepeMD10cMC3SBpAnAmcEJEvFjaNAeYLGk7SeOA8cA9/RFjRs3Sx63pluo5wJS0PAWo96lkwFNxinoZsDwivl3aNOD66Du0ukHSSmA74JlUND8iPp22fYliHHYTxUeWX9Tey8Am6SPAfwOjgOeAJRFxTNrWLH08DvgOm2+pPrefQ+oxSdcAH6CYgu9JYAZwIzAb2AN4DJgUEe2/9GoIkg4Dfg0sBV5LxWdTjLsOqD46uZqZZeBhATOzDJxczcwycHI1M8vAydXMLAMnVzOzDJxctwKSXpW0RNKDkm4uz+LVw/22lGdfanZpNrRlki7opN5qSSPT8u9qtZc0StICSYslvS937Nb3fCnWVkDS8xGxQ1qeBfyhN67pTPd2/ywi9u3pvhqBpI3AqIh4uZN6q4HWiHi6Xfnr7SVNBo6NiCm19lFnv4Mi4tVuhG79wGeuW5/fkyYokbSDpHmS7pO0VNLEVN6S5su8NJ1p/UrS0LTtQEn3S/o9MK1tp5KGSLoi7WexpCNS+SmSbkxnzKskfUbSF1Kd+ZJGtA9Q0pWSLknzdj4q6f1pntLlkq4s1Tta0u9T/Nen+82RdI6ke9OZ+sx0Vw+S7pR0vqR7JP2h1hmjChektkslfSyVzwG2Bxa0lZXa7JLeo8WSfkBp3gJJz9dofybwX8Bx6RPF0A76sjr15zfAJElvkfRLSYsk/VrS3qX37CJJv0vv2UmlGM5Ifblf0nmprN5+JqW+3y9pi+kmrYsiwj9N/gM8n14HAdcDE9L6YGDHtDySYvpAAS0Ud1/tl7bNBv45LT8AvD8tXwA8mJZPB65Iy3tT3CUzBDgl7feNFHd7bQA+nepdSHGHV/t4r6SYAlAUU8ltBP6e4mRgEbBfivduYPvU5kzgnLQ8orSvq4APp+U7gW+l5eOA22oc+x8p5usdBOyW+rF7+X2s0eai0rGPp5gAZmT7Nu2WTwG+V3rv6/VlNXBGqd08YHxaPhi4vfSeXZ/eo30oplOEYr7a3wHDyu9NB/tZCoxOy8P7+3e3kX88ccvWYaikJRRJcxFF8oAief2npMMpbiUcTZFQAFZFxJK0vAhokbQTxR/cXan8Koo/XoDDKG6XJSIelvRH4G1p2x1RzL35F0kbgJtT+VLgnXVivjkiQtJS4MmIWAogaVnqxxiKJPLbdGK6LcVZOcARks4AhgEjgGWlY7ZN9LEo7ae9w4Brovj4/aSku4B30/G8A4dTTCxORNwi6dkO6tZySAd9AbgOXp8J6r3A9do8C9R2pXo3RsRrwEOS2v4dj6L4T+/FFN/6TvbzW+BKSbPZ/F5ZNzi5bh1eioj9UnL8GcXH+YuAkynOJg+MiL+lscIhqU15XPFVYChFMq43SN/RnG/lfb1WWn+N+r+D5Trt2w9OMc2NiI9vEYQ0BLiYYsxzjaSvsLlP5f2+WufY3Z27ridfXogafSl5Ib1uAzwXEfvVqVd+n1R6bR9b3f1ExKclHUxxBr5E0n4R8Uz7etY5j7luRSJiA/BZ4Isqpm3bCXgqJdYjgDd30v45YIOKyTOgSM5t7m5bl/Q2igk0aj3JoLfMBw6V9NZ0zGHpuG2J9Ol0hnZSvR3UcTfwMUmDJI2iOCvtbNavct+PBbr6/KZ6fdlCRGwEVkmalOpJ0rs62fevgE9KGpbajOhoP5LeEhELIuIc4Gm2nJbRusDJdSsTEYuB+ymm2LsaaJW0kCI5PFxhF6cC31fxhdZLpfKLgUHpY/x1wCnRybfqPRER6yjGLa+R9ABFgto7/QdwKcWQw40UUwt2xU8pxpXvB26nGO/8cydtvgocLuk+4GiKcdrK6vWlTvWTgdMk3U8x3NHho2ki4pcUQxoL09DQFzvZzwXpy68HKf7TuL8rfbHNfCmWmVkGPnM1M8vAydXMLAMnVzOzDJxczcwycHI1M8vAydXMLAMnVzOzDP4fBGJTZM+3LC0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 360x216 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax  = plt.subplots(1,1, figsize=(5,3))\n",
    "ax.hist(rnt.randiff, 10)\n",
    "ax.axvline(x=rnt.dbar, color='r')\n",
    "ax.axvline(x=-rnt.dbar, color='r')\n",
    "ax.set_ylabel(\"Frequency\")\n",
    "ax.set_xlabel('Random mean of differences')\n",
    "print('Red line: observed difference')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}