AB testing.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Test if 'setosa' and 'versicolor' flower dimensions taken the iris dataset are statistically different using Permutation Test."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from sklearn.datasets import load_iris\n",
    "\n",
    "data = load_iris()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check our data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['DESCR', 'data', 'feature_names', 'filename', 'target', 'target_names']\n",
      "['sepal length (cm)', 'sepal width (cm)', 'petal length (cm)', 'petal width (cm)']\n",
      "['setosa' 'versicolor' 'virginica']\n",
      "(150, 4)\n"
     ]
    }
   ],
   "source": [
    "print(dir(data))\n",
    "print(data['feature_names'])\n",
    "print(data['target_names'])\n",
    "print(data['data'].shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Select two kinds of irises of interest."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "setosa = data['target'] == 0\n",
    "versicolor = data['target'] == 1\n",
    "setosa_data = data['data'][setosa]\n",
    "versicolor_data = data['data'][versicolor]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we can select the size of a 'setosa' sample set ('A') and 'versicolor' ('B'). By setting 'replace = True' we can go beyond the original sample size. Also lets scale the features by the respective mean values in set 'A'."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "setosa_sample_size = 50\n",
    "versicolor_sample_size = 50\n",
    "\n",
    "sample = np.random.choice(len(setosa_data), setosa_sample_size, replace = False)\n",
    "A_sample = setosa_data[sample, :]\n",
    "sample = np.random.choice(len(versicolor_data), versicolor_sample_size, replace = False)\n",
    "B_sample = versicolor_data[sample, :]\n",
    "\n",
    "scale_mean = np.array([np.mean(col) for col in A_sample.transpose()])\n",
    "A_sample = A_sample/scale_mean\n",
    "B_sample = B_sample/scale_mean"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here the author has to confess his complete lack of any knowledge in botany. Therefore lets treat available features (lengths of different parts of a flower) as equally important. We can use Principal Component Analysis to determine a linear combination of features that has the largest variance. This would give us a metric to check how close are sets 'A' and 'B' (Null hypothesis - they are the same species, so the underlying distribution is the same for both sets)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.decomposition import PCA\n",
    "\n",
    "pca = PCA(n_components=1)\n",
    "pca.fit(np.concatenate((A_sample, B_sample), axis=0))\n",
    "\n",
    "A_scaled = pca.transform(A_sample)\n",
    "B_scaled = pca.transform(B_sample)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Simple function for the Permutation Test. The output is the difference of means of 'A' and 'B' and an array ('shuffled') of differences of means for permutted sets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def perm_test(A, B, n):\n",
    "    l_A, k = len(A), 0\n",
    "    shuffled = np.zeros(n)\n",
    "    delta = np.mean(A) - np.mean(B)\n",
    "    S = np.concatenate([A, B])\n",
    "    for i in range(n):\n",
    "        np.random.shuffle(S)\n",
    "        shuffled[i] = np.mean(S[:l_A]) - np.mean(S[l_A:])\n",
    "    return delta, shuffled\n",
    "\n",
    "delta, shuffled = perm_test(A_scaled, B_scaled, 1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets plot the histogram for 'shuffled'."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAngAAAGgCAYAAAA5EhrlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdd5hU5dnH8e9NlS4iHQUbKggooCJGQRGDYnkj9hY0CEokNoJIjCHGxIaUmNgVxBYjsStYUMSGCHbBqMgiERTQAFFAYPd+/3jOrLOzM7tn++7w+1zXXLtz6n36Pc95znPM3RERERGR7FGrqgMQERERkfKlBE9EREQkyyjBExEREckySvBEREREsowSPBEREZEsowRPREREJMsowavmzGyOmc2p6ji2BWY20MwWmtlGM3Mz61SKaUwzs03lH50UxcyGlmGblXrc0jCzXmb2qpn9L5pv/8qYb1Uzsxwzm1bVcUjZ6bpUMyjBq0RJF5I+GfrfZmZlbpjQzPYxs/GVdcHKBmbWDPgnYMAo4CxgdYZhO0Trd99KDFEqkZmdaWYXV8B0awMPA+2B3xL2s8XlPR/ZtpXHNcDMfhZNY/vyi6x6qajjvLqoU9UBSLGOLMU4+wB/AOYAOeUZTBbbD9geuNrdHy9m2A6E9ZsDvFfBcUnVOBPYC5hcztPdCdgNuMTdbyvnaYsklMc14GfRNKYBa1P6lea6VB1V1HFeLagEr5pz983uvrmq4ygpM2tQ1TGUUKvob+qJTKQ8lft+Zmb1opLBkoxT047PGsvMGlZ1DOWtpl6Xtjnurk8lfYChgAN9MvS/LWySAt3mAHNSup0MvA2sB9YBHwJXpcwj9TM0afzjgLeADcB/gceAvdLE87NouE3AMmAMcE40vU5Jw+UAs4DDgHnR8OOT5vU08BXwYzSd64H6KfOaFo3XHngc+B+wErg06r9HNI/vo2mNKMF6/xXwQTT9VcB0oF3KOk5dX3MyTKt/hvU7PmU5WhNuxa2P1vEdwHZppndq0rZYBzwFdImxTIk4TgPGAcujdfM40IJQOn9ttA5/AP4BNE4znYHAy9H6/iH6v2/KMB2BW4BPojjXRnF2LSKmS6P9YhPwJrBfyrCtgbuiuH8EVkT7SbcYy94HeCOa9pfAWNLslyVYvqHJ42bYHzxp+NHA68CaKIYPgV/FiHtamunmJPXfB3gyWr8bovV2VIZ1fAbw+2j5c1OXO2WcHMrv+Cx2vwbqRdP4OlrnLxJKSXKAaSnD7gw8kLQu3wPOThmmU7TMY4HzgM+j9fMSsEs0TGJ/2xgta9sSnI8PA/4WxfA9MCPd+EDvaF2tjeYzHzimiGlOjtaBp/TrB0wknIvWAfcA2wGNgNujONYBfwfqpFkPQ9PElr9uKeYaABxCqI6yjJ+OvduB7ZOmNz7DNPonHSOp16UGwHVRLD8CS4FrKLwv5UTb6ADgtWhdLgcujnk+7wk8S6hCsxH4gnBOb5Q0jAEX8tN5/1vCftYhaZg56ZYxqX/G62xN+egWbdVoZmY7pum+XXEjmtkRhIv1S8AVhJP7nsCh0SBzCSeGXwN/4af6PW9E459G2NHfB64EmhHqnL1hZr3cfWk0XA/gecKB8SdgM+Hk+kOG0HYH/kW4aN9DuPAAnAtsBW4mXBD6Ei6QOwGnp0yjFjCTcBEaA5wC3GRm66NYHyNcAM8FbjOz+e7+bjHraywh0ZlLqPO0M+HAP8TM9nP3tcCfgY9S1tk3GSa5mHDyG0+4uL0adf8gZTmej7qNAQ4krLs1hGQsEdvlhBPio4QTVGNgJPC6mfVMbItijCFsmxsJidjFwFTCBbsz4QTbHRhOuNjk1zeJ9oX7CUnPlVHc5wIvmVk/d38rGnR/wkXpUcLJuR1wPjDXzLq6+9cpMV0G1CVs87qE9f64me3u7luiYWYA3QgX1qVAS8I+vCfhRJqWmXUhJAzro2XbHC3b92mGjbt8qf4MNI+W85I0/S8h7Kf/JFwUfgHcZWa13f2OTLETLqJLKbjvfB/F2pmQNP4ITIq6nwM8bWYnuvtjKdO6Ipr3FMLFrNDypyiv47PY/TpazqGEbfwyYf95gZTzW3QOfB3YIZr/CsIPnnvNrLm7T0mZ/8lAQ+BWoAlwOWG/+gcwhLDe2hH2v1uB/ytmnSRMIZzXriYcQ6OAPaPz4eYo1n7Ac4R98xrCdjoFeDLaPo+mTPNmQiL4F8I5NtlkQnLyR6AXYTv/QLh1v4GQuB9OOBd8RslvHxZ5DSCsxx0I+8LXQA9gGOEHxsHRMI8SkvJTCPv7mqh72vqiZmbROIMI57K3CAUEvwO6Eo6RZLsQfiROI1yPTgEmmdkid38+04KZWUvCvrQGuIGwz+5E+KHSmJ+uT7dGy3Qf4cdpG8J27Zty3k97nMe4ztYMVZ1hbksfMv+ySvsLIhpnDkm/lAgnsXVA7SLmcypJv7aSutcllOh8QsFfO/sRduD7k7o9Qfjl0zGp247Ad6QvwXPguDSxNEzT7Uogj4K/pqZF0/h9UrfGhAt5HnBuUvd2RBelYtb3jtEyvEzBX8LHR/O6prh1lmG6fcj8S3pa6rSj7o8Dq5O+7wxsISpJSerelnBhuLuYGPpH81kM1Evqfk/UfQ5QK6n7s9G6tOh7o2hbTkvdXoQkZHYx23D3aN3+Lk1Mn5FUqkO40DowOPreLPo+uhTH0L8ISd1uSd1aRussf78s4fINTbNPzyKpdC3GPv0C8Flp9x1CMrSFpJJ0oCmhlGVZYlsmreMvSTqGi5lnDuV3fBa3X3eLhrs7Zbi/RN2nJXWbEHUbkNStHiE5+B5oFnXrFA23BmieNOxVUfdPSSolIlzQ84CWxayXxHb/IGX8c6Puw6PvRjhnvkTBY6oWoZR1SZppvkXSOSel3wtEx2HU/aUo3vtShl8EfJj0PbEe0p13clLWbcbzWYZtfmY0/MFJ3caSplQ86jeHgtelY6Jh/5Qy3KSo+xFp9sfkbvUJP6ofKWabJc7dvYsYpm+69UT4obuFgteYtMc5Ma6zNeGjOnhV4zeEW0epn6djjLuecPEqTSXXXoRfMre4e35JnIdSsBeBoy2oHcXzpLsvSxpuDfBghmn/x92fTO3o7hsAzKyWmSVKLucSTpo900znrqRxvwc+JiSf9yV1X0Eo0t+tmOU9gnDimOTuW5PGfwL4NzC4mPHL4paU768AO5pZk+j7CYTbqP8wsx0TH8IJ6C3CL/g4pnvBujCJUqlp7p6X0r0JIRmCsH2bAw+kzL8hYV84xMzqwk/bEEJ9IjNrQUioPiXsU6mmuntyUzGvRH93jf5uipazv5ntEHM5E0+gDgKecvclie7uvppQCpAs9vKVVNI+XdfMdoim+xKwe/Q0domkLNcnSfNZT6i2sTMhcUo2PfkYjqG8js/i9utjor+ppU6T0kzrGOBdd5+dFM/maNxGhGQ22b/c/b9J3xP7+oPu/mNKdyOUEsVxe8r40wn799HR9x6E0pv7gR2S9qUdCCW5u5pZx5Rp3pl8zklxj0dZREq8d6cM9xbFn+NKLGmbm5k1TSpJhfTHcxyJBO+mlO43RH9Tz7WfufuLSTH9SLhzsytFW5+YXxHH78mEHwjPphz7Kwg/PuOcW8tyna02dIu2arzt7vNSO5rZiTHGvZWwAz9rZisIF6tHCcmYFzlm+AUI4ddoqkWEnbkp4SLYgHAwpPo0w7S/SNfRzLoSbh/2j6aZLPXx+y3uvjKl2zrga//p1l5y9+YZYknoFP1Nt7yLKXwBKS9boiQ0WeLCtAOhPljnpDjS2ZChe6ovU76vi/4uz9C9OaHuT2L+GW+HELbPajPbjnD76kxCCWOyNYXGCiVO+dz9v+EODjtE33+Mbp3fAHxjZm8RShjvd/fU5UnWkrBv/jtNv9RusZeviP5pmdnxhNto+wKpDzc046d1HVdLwsUk03EJYV9+P6n7ksKDFqm8js/i9uuO/FSqls/dV5vZf1PG7UQ4d6VKXuZkpdnX4yhwnnP3rWa2lLAs8NO+lJqAJWtFwf2+qO1TkuVoYGb1UxLQMjGznQjb/GjCj75kpW0SpRPhPF3g4SF3X2lmaym8LZdR2H8JpWxFmUPYZ/4AXGpmrxBu9T4YFQhA2F6NyVzNprjrJJTtOlttKMGrYdz966h+3EDCr/5BwNnATDMbXIadz8o43MZCA4bSjJcJyco4wklvI+FBimkUfoo7j/RySxhLHGUZtziZliN5vollP4pwuzlVpmWOO1xx6ywx/6GE+nrpJC48Uwj1WW4m/NJfR1jGyaR/Er/Y7eXuE83sMcItl4GEhGmcmR2fXKKTYfx0+3jq9izJ8sVmZgcT6oK+BowgVHnYTLhYXkL5t0wQ+3grRkUen8lxFnVclfUcU9p9vTjF7U+J9TAWWJhhGqk/MIraPqVdjqLO7bGeojazRD3KloS6yYsJ9dZqEW5XVsRdvXTboVTbLLq+DTGzAwilhgMJdT7HmdmB7v4NYRm+JdymTqfYku8KvM5WKiV4NVB0G+MZ4Jmocuu1hArHfQkX4Ew7X070dy8Kl2zsRfgFtZ5QvL2R8ORqqnTdMjmMcCLp7+6J23SY2cASTKMscqK/e1G45HEvSt8+VHkc3Ilf+F+6+6Iih6wYifmvTr5VksEphFuCBRoENbPmpC/Bi8XDQySTgclRqcK7hItopgRvFSEZ2StNv84p30uyfGnDy9D9JMIt5iOTb0Ob2WGlmEfCasJFJ91yJbrllGH6mVTE8ZlDuEh3JulhmahyfGrpUA6Vv8zpdCbUiwPAzOoQSpwS6ySxL/2vlPtSefku+ltgPZpZfQqXrGfaf7sT1u9Qd783aRrpzuslOc/lAEea2fbJpXhm1oZQqp1TgmkVy93nE55ivsrMjiLcARhGeHBiCSExe8vd/1fcpIqYR3HX2WpPdfBqmKj+U77ol0TiSdLELYkfUr4nLCQ8NXVBcjtYSb9UnvUgl3DCOy65bklUjyH1ybqiJH6l5f8qi35BXlqCaZTFi4Sn3S6ypHbCzOxYQp2aZ0o53UzrtyRmEEru/hitkwKiC2JFmkWoZ3RldIEoav65pPyyjp5QbVeaGUf1+ArcDnT35YQELuM6jfbL54BjzSy/blIUa+p+WZLlS+cH0t+uyiVcFPK3WZTonlvM9DKKlmsWoV5RfqIa1WsbQbill/HJ4jKoiOMzcUylvh0g3dPITwP7JSfHUb2qiwiJ/JwyxFESI1L2kbMJ235m9H0h4Tbu6HR1LCvhWAUgSlZWExLzZOdTuAQv0zmq0DaPjE4zy5Kc556Oppm6nX8b/S3tubYAM2seJVvJ3on+JuL8B+H4HJ9mfLOCLVikPc5jXmerPZXg1Tx3RTvobEKdjfaEZj9W8tMvzncIt1OusPCamY2EXzNLzexSQoX0183sPn5qJmUd4TZZwh+AnwOvmdmthErx5xF+iTUn3q+71wlF5fea2c3RNE4k1I+ocO6+xszGE355vWhmjxIeqR9FWI4JpZz0Z4SSzgvM7HtC3aOP3P2jEsS21MzGENrDmhfF9h2h3s/RhArW55cyvjjz/5+ZDQceAt43swcIyX8HwgXkB8LtYwhN05xtobmajwh1z04hQ72uGDoTmip5hPAQzY+EZd6bny4ImVxF2C9fNbO/E/ap4YQ6Pfkn6hIuXzoLCLeCphC2RZ67/4NQ3+dS4IXo+NmBcFx8TXiAqbR+R/iRlViuRDMpOwMnpTwwU17K/fh09w+i9XKumTXlp2ZSBlK4tPc6wm20p8zsr4Rz2CmEJ40vcfcS30Ivg5fN7CFCyd0oQj3AaQDunmdm5xJ+XCwys3sI+1tb4CDCMdu1kuK8nfCjZSrhoYTewAAKr9u01wBCPc/PCM1PdSCcc44iHBepFkR//xKtm83AS+6+Ks2wzxB+pFxlZjsT2o/rS2iv8YlyLPn8JfDrqHrHEkK90XMIiesMAHd/NdqfLo0KL2YSfjDsQmiu5W7CvpdYxnTHeZzrbPVXUY/n6lP4Qzk0dExo72kW4YLyI+HX/TSiBj+ThvsV4UDeSsoj44R6T/MJB/1aQlMH6Ro6PjQaLtEA6hjCyc+B1knD5QCzMizTgYT6Sj8QKr3ewk9NKSTHNA3YlGb8TI+xvwfMi7nehxFKQH4knAinA+1ThondTEo0/LHRNDdH440vZjkS275TSvfBhIvgesJJ6DNCO3b7FzP//tH0To2zHIQSFU/dzoST8NOEE/0mQhMi/6BgEwZNCe22fRNtx1cIF+3UfTNtTFG/5HXUglCfbxEhOV5PONGeW9Qyp8T8JvEaOo6zfIW2DeGhh+mEBCiPgg2gnk2ou7SJ0OjupZnmnyb2oprY2YeQQK4jHJtvAkfH2e7FzDOHijs+0627eoRK/KuI19Dxg9F63kR4kCRjQ8dx1iU/NcszqJj1koj9MEK7cYmGjh8lqSH0pOG7Ehp5XkU4lyyPtteJaaZZ6ByfqR8ZmiMhVF9wCjY5tB3hAYDvonX7DOHJ03TrNu01gHD3Yla0n30Xrf/WJB2jSdP4HT81pp1/XiFzQ8fXR8NvjmLK2NBxmvUzjQxNEyUNsx+hgCKHnxqufw44NM2wZxOS4B8I55nFhHY39yzuOCfmdba6fxJtYonEEv3SGU54K0LcBwFERKoVMxtK+DF1kKdp1UCkplMdPMkotZ5UVNfkLGCukjsREZHqS3XwpCg5ZnY/oQmAdoRbnY0JbaKJiIhINaUET4ryDOGNC20I9TgWAGe4e414RFxERGRbpTp4IiIiIllGdfBEREREskzW3qLdcccdvVOnTlUdhkiV+N//whuVmjQp7bvDRUSkVBZGb7TrVfHn34ULF65x97SNbWftLdrevXv7ggULih9QJAvNmRMae+/fPzuPbxGRaivxso1KyK/MbKG7907XT7doRURERLKMEjwRERGRLKMET0RERCTLKMETERERyTJK8ERERESyjBI8ERERkSyjBE9EREQkyyjBExEREckylfomCzM7FBgN9ALaAee4+7RixukG/A04APgOuB34k2drC80iUqHWr1/PqlWr2LJlS1WHIlKt1K1bl1atWtG0adOqDkXKQWW/qqwx8BEwPfoUycyaAi8Ac4H9gT2BacAPwE0VFqWIZKX169fzzTff0L59exo0aIAlWpwX2ca5Oxs3buSrr74CUJKXBSo1wXP3Z4FnAcxsWoxRzgAaAr90943AR2a2N3CpmU1UKZ6IlMSqVato3749DRs2rOpQRKoVM6Nhw4a0b9+eFStWKMHLAtW9Dt5BwKtRcpfwHOH2bqcqiUhEaqwtW7bQoEGDqg5DpNpq0KCBqi9kieqe4LUBvknp9k1SPxGREtFtWZHMdHxkj+qe4AGk3oa1DN0xs+FmtsDMFqxevbriIxMRERGphqp7gvc1hUvqWkV/U0v2cPc73L23u/du2bJlhQcnUhadxj5T6CMilePCCy+kf//++d+HDh3KMccck/89Ly+PESNG0KJFC8yMOXPmpO0mUl1V9wTvTeAQM9suqdtAYAWQUyURiYhUkXfffZfatWtz8MEHZxxm8+bNNGvWjPfff585c+ZgZvmfhg0b0r17d6ZPL7YRg23OlClTuP/++/O/P/vss0ydOpWnnnqKlStX0rdv37TdRKqrSk3wzKyxme1rZvtG8945+r5z1P9aM5udNMqDwAZgmpntY2YnAGMBPUErItucO++8k5EjR/LRRx+xePHitMO8/PLLNG/enB49euR3+/jjj1m5ciUfffQRZ5xxBkOHDmXevHmVFXaN0KxZM7bffvv8759//jlt27alb9++tGnThnr16qXtVhp6iEEqQ2WX4PUG3o0+DYA/Rv9fHfVvC+yWGNjd1xFK7NoBC4C/E9q/m1h5IYuIVL2NGzfy4IMPct5553HiiSdy9913px3uiSee4Pjjjy/QrVWrVrRp04Zdd92Vyy+/nB122IF3330347wSJX8zZ86kV69eNGjQgEMOOYT//Oc/vPLKK/To0YPGjRtzzDHH8O233xYYd+rUqXTp0oXtttuOzp07M2nSJPLy8vL7T5w4ke7du9OoUSPat2/PsGHDWLt2bX7/adOm0bhxY2bPns0+++xDo0aNOOyww1i6dGn+MMuXL+f4449nhx12oGHDhuy111784x//yLg8ubm5jB49mubNm9O8eXMuvvhicnNzCwyTfIt26NChXHLJJXz55ZeYGZ06dUrbDUL7cTfccAO77bYbDRo0oFu3bgVKAnNycjAzHnroIQ4//HAaNGjA7bffDsAbb7xBv3798psnueCCC1i/fn3+uP3792fkyJGMGzeOHXfckVatWjF69OgC63Pz5s2MGzeOjh07Ur9+fXbddVf++te/5vdftGgRgwcPpkmTJrRq1YrTTjuNr7/+OuO6kizi7ln56dWrl4tUZx0vf7rQp7y8/DL+8suU2/SyxaJFi6o6hFKbPn26d+/e3d3dX375ZW/ZsqVv3ry5wDB5eXnevn17nz17dv5wgK9evdrd3XNzc/2xxx7z2rVr+9tvv51xXonx9t9/f587d66///773rVrV+/bt68ffvjhPm/ePH/77be9U6dOfuGFF+aPd8cdd3ibNm38kUce8S+++MKffPJJb926td988835w0yaNMlnz57tS5cu9Tlz5ni3bt38zDPPzO8/depUr1Onjg8YMMDfeustf//9933ffff1I488Mn+YY445xo844gh/7733/IsvvvCZM2f6zJkzMy7P9ddf702bNvWHH37YFy9e7BdeeKE3adLE+/Xrlz/ML3/5Sx88eLC7u69du9avuuoq79Chg69cudJXrVqVtpu7+7hx47xz584+c+ZM/+KLL/yBBx7whg0b+tNPh+N56dKlDnjHjh3z18vy5cv9gw8+8EaNGvmECRP8008/9Xnz5nmfPn18yJAh+TH169fPmzZt6r///e/93//+tz/88MNeu3Ztf/DBB/OHOfXUU719+/Y+Y8YMX7Jkib/00kt+7733urv7ihUrvEWLFj5mzBhftGiRv//++37MMcf4/vvv77m5uRnXV00+TqoFCJ9KmRULPEMeVOWJWEV9lOBJdacEr/KlvXAlTsaV/SmhQw891G+88UZ3D4lcx44dfcaMGQWGmT9/vjdv3ty3bNni7j8lao0aNfJGjRp5nTp1vHbt2n7TTTcVOa/EeLNmzcrvdvPNNzvgCxcuzO/2hz/8wbt27Zr/faeddvLp06cXmNakSZN87733zjivmTNner169fITjqlTpzrgn3zySf4w999/v9etWzd/mG7duvn48eOLXIZkbdu29WuuuSb/e25uru+xxx4ZEzx39xtvvNE7duxYYDqp3b7//nvfbrvtfO7cuQWGu+iii/yoo45y958SvAkTJhQY5qyzzvJzzz23QLd3333XAf/mm2/cPSR4ffr0KTDMEUcc4b/61a/c3f3TTz91IGNy+/vf/94PP/zwAt2+++47B/ytt95KO467ErwyqyYJXmW/qkxEREro888/5/XXX+ehhx4CQltlZ5xxBnfddRdDhgzJH+6JJ55g8ODB1KlT8NSeqJf3448/Mn/+fEaNGkWTJk0477zzipxv9+7d8/9v3bo1AN26dSvQbdWqVQCsXr2a5cuXM2LECC644IL8YbZu3RpKEyIvvfQS1157LYsXL2bdunXk5uayefNmvv76a9q1awdA/fr12XPPPfPHadeuHVu2bGHt2rXssMMOXHTRRZx//vnMmjWLAQMG8Itf/IJevXqlXYZ169axcuVKDjrooPxutWrV4sADD2T58uVFLn9xFi1axKZNmxg0aFCB9uO2bNmSfws3oXfv3gW+L1y4kM8//5yHH344v1tiPS1ZsoRWrUKDEcnbAMK6SKzzd999l1q1anHYYYeljW/hwoXMnTuXxo0bF+q3ZMkSDjjggJhLKjWREjwR2bYlJR/V1V133UVubi4777xzfrdEMrB8+XJ22mknAB5//HHGjx9faPxddtmFHXfcEYCuXbvy1ltvcfXVVxeb4NWtWzf//0QCk9otUR8s8fe2227L+HTpsmXLGDx4MOeddx5XX301LVq04J133uG0005j8+bN+cOlJqiJeSfm8atf/Yqf//znPPvss7z44ov07duXK664Iu2yV6REPE899VSBbQMF1xNAo0aNCo07bNgwLrnkkkLTbd++fcbpJK9zL2bfzcvLY/DgwUyYMKFQv0TCLtlLCZ6ISDW2detW7r33Xq699toC7bQBnHXWWUydOpWrrrqKJUuW8PnnnzNo0KBip1m7dm02bNhQrnG2bt2a9u3bs2TJEs4+++y0wyxYsIDNmzczadIkateuDcDTTz9dqvl16NCB4cOHM3z4cK6//nqmTJmSNsFr1qwZbdu2Zd68eRx++OFASIzmz59P27ZtSzXvhC5dulC/fn2WLVuWP+24evbsyccff8zuu+9e6vn37NmTvLw8Xn755bTbvWfPnvzzn/+kY8eOhRJFyX5K8EREqrFnnnmGNWvWcN5559GiRYsC/U499VRuvfVWrrzySp544gkGDBiQ9nbcqlWr2Lp1a/4t2vvuu48TTzyx3GMdP348o0aNYvvtt+foo49my5YtvPPOO3z11VdcccUV7LHHHuTl5TF58mROOOEE5s2bx+TJk0s8n4suuoijjjqKzp07s379embNmkWXLl2KHP7aa6+lc+fOdOvWjVtuuYWVK1eWOcFr0qQJo0ePZvTo0bg7hx56KN9//z3z5s2jVq1aDB8+POO4l19+OX369OH8889nxIgRNGnShE8++YSnnnoq/ynb4uyxxx6cfPLJDBs2jClTptCzZ0/+85//kJOTw1lnncWvf/1r7rzzTk455RQuv/xyWrZsyRdffME///lPbrrpJpo0aVKm5Zfqrbo3dCwisk27++67OeywwwoldwAnnXQSy5Yt48UXX+Txxx8v1DxKQteuXWnbti277747l19+OSNGjODmm28u91iHDRvGPffcw3333UePHj045JBDuOOOO9hll12AUJ9sypQpTJw4kS5dunDXXXelvX1YnLy8PEaNGkWXLl0YOHAgrVu35t577804/GWXXcY555zDsGHDOPDAA8nLy+OMM84o9XIm+49fZC8AACAASURBVNOf/sT48eOZMGECXbt2ZeDAgfzrX//KX+ZMunfvzty5c8nJyaFfv3706NGDK664osS3TqdPn87pp5/Ob37zG/baay+GDh3KunXrgFBf7/XXX6dWrVoMGjSIrl278utf/5r69etTv379Ui+z1AxW3D38mqp3796+YMGCqg5DJKN0rybLuW5wuUx7zpxQZ6l//+w8vktr8eLF7L333lUdRrlbs2YNbdu2Zfny5bRpk/p2R5GSydbjpNIkHriphPzKzBa6e+90/VSCJyJSw3333XdMnDhRyZ2I5FMdPBGRGq5z58507ty5qsMQkWpEJXgiIiIiWUYJnoiIiEiWUYInIiIikmWU4ImIiIhkGSV4IiIiIllGCZ6IiIhIllGCJyIiIpJllOCJpNFp7DOFPiLyk3322Yfx48dXyLQvvPBC+vfvXyHTLo1p06alfcdvSc2ZMwczY82aNeUQlUjRlOCJiFRzq1evZuTIkXTq1In69evTunVrBgwYwAsvvJA/TKdOnUr1XteaKCcnBzPL/zRu3Jg999yTYcOG8cEHH5T7/E455RS++OKLEo2Tbnv07duXlStXpn2vsEh505ssRESquSFDhrBhwwbuvvtudt99d1atWsUrr7zCt99+W9WhlVpeXh7uTu3atUs9jVmzZtGjRw82btzI4sWLue222+jVqxf33Xcfp556arnF2qBBAxo0aFDm6dSrV0+vk5NKoxI8EZFqbO3atbz66qtcd911DBgwgI4dO7L//vszevTo/CSmf//+LFu2jN/+9rf5pVoA3377LaeddhodOnSgQYMGdO3alalTpxaYfv/+/Rk5ciTjxo1jxx13pFWrVowePZq8vLz8YVatWsXxxx9PgwYN6NixI/fcc0+hOCdOnEj37t1p1KgR7du3Z9iwYaxduza/f+I257PPPss+++xDvXr1WLx4Mbm5uYwePZrmzZvTvHlzLr74YnJzc2OtmxYtWtCmTRt22WUXjj76aJ588klOOukkzj///ALzfuONN+jXrx8NGzakffv2XHDBBaxfvx6A22+/ndatW7N169YC0z799NM5/vjjC8SesGTJEo4//njatGlDo0aN6NmzJ08//XSBdZpue6S7Rfvoo4/SrVs36tevz0477cSf//xnPOkl9Z06deKaa65hxIgRNG3alA4dOnDjjTfGWj+ybVOCJyJSjTVu3JjGjRvz5JNPsmnTprTDPProo3To0IGrrrqKlStXsnLlSgA2bdqUn3x8/PHHXHTRRYwYMYLZs2cXGP+BBx6gTp06vPHGG/ztb39j8uTJPPzww/n9hw4dyueff86LL77I448/zvTp08nJySkwjVq1ajF58mQ+/vhjHnzwQebPn8+oUaMKDLNp0yauueYabr/9dhYtWkTHjh256aabuPPOO7n99tt58803yc3N5YEHHij1+ho9ejTr1q3jxRdfBODDDz/kyCOP5LjjjuP999/n0Ucf5b333uPcc88F4OSTT2bt2rX5wwP88MMPPPHEE5x55plp5/H9999z1FFH8cILL/D+++8zZMgQTjjhBD755BMg8/ZItXDhQk466SROOOEEPvzwQ6677jquvfZa/va3vxUYbtKkSXTr1o133nmHyy+/nDFjxvDmm2+Weh3JtkG3aEVkmzZnjlXJfPv39+IHAurUqcO0adM477zzuOOOO9hvv/04+OCDOemkkzjwwAMB2GGHHahduzZNmjQpcAuwffv2/Pa3v83/Pnz4cF566SUeeughBgwYkN+9S5cuXH311QB07tyZO++8k9mzZ3Paaafx6aefMnPmTF577TUOPvhgAO6991523XXXAnFefPHF+f936tSJG264geOPP557772XWrVCWUJubi4333wzvXr1yh928uTJjBkzhpNPPhmAKVOm8Nxzz8VaN+l06dIFIL/O3I033sgpp5zCZZddlj/Mrbfeyn777ceqVato1aoVRx99NA888ACDBg0C4LHHHqNOnToce+yxaefRo0cPevTokf/9d7/7HU899RQzZszgyiuvzLg9Uk2cOJF+/frxxz/+EQjr/rPPPuP6668vkBwfeeSRXHjhhQCMGjWKv/71r8yePZuDDjqoNKtIthEqwRMRqeaGDBnCihUreOqppzjqqKN444036NOnD3/5y1+KHC83N5c///nPdO/enRYtWtC4cWMeffRRvvzyywLDde/evcD3du3asWrVKgAWL15MrVq1OOCAA/L7d+zYkXbt2hUY56WXXmLgwIF06NCBJk2acMIJJ7B582a+/vrr/GHq1KnDvvvum/993bp1rFy5skCiUqtWrfzEtTQStzcTt0UXLlzI/fffn18S2rhx4/xEdcmSJQCceeaZPP7442zYsAEIJZonnngi2223Xdp5/PDDD4wZM4YuXbrQvHlzGjduzIIFCwqt1+IsXrw4P5aEn/3sZ3z11Vf5t5Ch6O0jkolK8ERkmxa3JK2qbbfddgwcOJCBAwdy1VVXMWzYMMaPH8/o0aOpV69e2nEmTJjATTfdxJQpU+jWrRuNGzdm3LhxhZKDunXrFvhuZvl18JLrg2WybNkyBg8ezHnnncfVV19NixYteOeddzjttNPYvHlz/nD169cv00MVcSxatAggv4QxLy+PYcOGcckllxQatn379gAcc8wx1KlThyeeeIIBAwbw4osv8vzzz2ecx+jRo5k1axYTJkxgjz32oGHDhpx99tkFljUOd89PRFMldy9q+4hkogRPRKQG6tKlC1u3bmXTpk3Uq1ePevXqFXo44bXXXuPYY4/lrLPOAkJC8emnn7L99tvHns/ee+9NXl4eb7/9Nn379gXgyy+/ZMWKFfnDLFiwgM2bNzNp0qT8BC75oYNMmjVrRtu2bZk3bx6HH354fozz58+nbdu2sWNMNmHCBJo1a8YRRxwBQM+ePfn444/ZfffdM45Tv359TjzxRB544AHWrFlDmzZt6NevX8bhX3vtNc4++2yGDBkChLqFS5YsoXPnzvnDpNseqbp06cJrr71WaNqJUlCRstAtWhGRauzbb7/l8MMP5/777+eDDz5g6dKlPPLII9xwww0MGDCApk2bAqHe26uvvspXX32V/5Rm586dmT17Nq+99hqffPIJF154IUuXLi3R/Pfcc08GDRrEiBEjePPNN3nvvfcYOnRogWZD9thjD/Ly8pg8eTJLly7loYceYvLkybGmf9FFF3HDDTcwY8YM/v3vf3PxxRdnfCgh3br5+uuvWbp0KTNnzuS4445jxowZ3HbbbTRr1gyAyy+/nPnz53P++efz7rvv8vnnn/P0008zYsSIAtM688wzee6557jttts4/fTT8+sNptO5c2cee+wx3nnnHT788EPOPPPMQg/ApNseqS677DJeeeUVxo8fz6effsoDDzzATTfdxJgxY2Itv0hRlOCJiFRjjRs3pk+fPkyZMoV+/frRtWtXxo0bx+mnn17gSderr76a5cuXs9tuu9GyZUsArrzySg444ACOOuooDj30UBo1asQZZ5xR4himTZvGLrvswuGHH86xxx7L6aefTqdOnfL7d+/enSlTpjBx4kS6dOnCXXfdFbvR5csuu4xzzjmHYcOGceCBB5KXlxc7xkGDBtG2bVu6du3KRRddRMuWLVmwYEGBNvC6d+/O3LlzycnJoV+/fvTo0YMrrriC1q1bF5jWoYceSvv27Vm0aFHGp2cTJk6cSKtWrTjkkEM46qij6NOnD4ccckiBYdJtj1Q9e/bkkUce4V//+hf77LMPY8eOZezYsfkPVIiUhcWpX1ET9e7d2xcsWFDVYUgNle7VZDnXDa4x80g8GVpT6pdVlsWLF7P33ntXdRgi1ZqOkzJK1J+shPzKzBa6e+90/VSCJyIiIpJllOCJiIiIZBkleCIiIiJZRgmeiIiISJZRgici25RsfbBMpDzo+MgeSvBEZJtRt25dNm7cWNVhiFRbGzduLPTmDKmZlOCJyDajVatWfPXVV2zYsEElFSJJ3J0NGzbw1Vdf0apVq6oOR8qBXlUmItuMxFsfVqxYwZYtW6o4GpHqpW7durRu3Tr/OJGaTQmeiGxTmjZtqguYiGQ9JXiyzUv3RgkREZGaTHXwRERERLKMEjwRERGRLKMET0RERCTLKMETERERyTJK8ERERESyjJ6ilayU7snYnOsGV0EkIiIilU8leCIiIiJZRgmeiIiISJZRgiciIiKSZZTgiYiIiGQZJXgiIiIiWUYJnoiIiEiWUYInIiIikmWU4ImIiIhkGSV4IiIiIllGCZ6IiIhIllGCJyIiIpJlKj3BM7ORZrbUzDaZ2UIzO6SY4U83s/fMbIOZfW1m95tZm8qKV0RERKSmqdQEz8xOAaYAfwH2A94AZprZzhmGPxi4D7gX6Ar8H9AFeKBSAhYRERGpgSq7BO9SYJq73+nui919FLASuCDD8AcB/3H3Se6+1N3nATcDB1ZSvCIiIiI1TqUleGZWD+gFPJ/S63mgb4bRXgfamtmxFuwInAo8W3GRioiIiNRslVmCtyNQG/gmpfs3QNo6de7+JnAa4ZbsZmA1YMAv0w1vZsPNbIGZLVi9enV5xS0iIiJSo1TFU7Se8t3SdAs9zLoAfwX+RCj9G0RIBm9PO2H3O9y9t7v3btmyZflFLCIiIlKD1KnEea0BcilcWteKwqV6CVcA8939xuj7B2b2A/Cqmf3O3ZdXTKgiIiIiNVelleC5+2ZgITAwpddAwtO06TQkJIXJEt+t/KITERERyR6VWYIHMBG4z8zmEx6gOB9oB9wGYGbTAdz97Gj4p4A7zewC4DmgLTAZeMfdv6zk2EVERERqhEpN8Nz9YTNrAVxJSNY+Ao5292XRIDunDD/NzJoAFwI3AeuAl4ExlRe1ZItOY58p9/FzrhtcpmmKiIhUhMouwcPdbwFuydCvf5puNxPavhMRERGRGPQuWhEREZEsowRPREREJMsowRMRERHJMkrwRERERLKMEjwRERGRLKMET0RERCTLKMETERERyTJK8ERERESyjBI8ERERkSyjBE9EREQkyyjBExEREckySvBEREREsowSPBEREZEsowRPREREJMvUiTOQmfUDNrn7W9H3ocAw4GPgMnf/vsIiFKlhOo19plC3nOsGV0EkIiKyrYpbgjcZaANgZnsCtwMfAAcBN1ZMaCIiIiJSGnETvN2AD6P/hwAvuPtI4Dzg2IoITERERERKJ26C50Dt6P8BwKzo/6+BFuUdlIiIiIiUXtwE723g92Z2FnAIMDPq3omQ5ImIiIhINRE3wbsE2Bf4G/Bnd18SdT8JeKMiAhMRERGR0on1FK27fwh0T9NrNLC1XCMSERERkTKJVYJnZi+Z2fZpetUDni/fkERERESkLOLeou1PSOZSbUeokyciIiIi1USRt2jNrGfS1+5m9l3S99rAz4GvKiIwERERESmd4urgLSA0keKkvxW7ERhV3kGJZJt0b7cQERGpKMUleLsABnwBHACsTuq3GVjl7rkVFJuIiIiIlEKRCZ67L4v+jVtXT0RERESqWKxmUgDMbCfCAxWtSEn43H1iOcclIiIiIqUUK8EzszOAewht3q0m1MlLcEAJnoiIiEg1EbcE72rgJuD3qnMnIiIiUr3FrVvXGrhLyZ2IiIhI9Rc3wXsWOLAiAxERERGR8hH3Fu0LwPVm1hX4ENiS3NPdHy3vwERERESkdOImeLdHf8el6eeEt1qIiIiISDUQK8Fzd7WDJyIiIlJDxG4HT0QqXkleaZZz3eAKjERERGqyuO3gXVpUfzV0LCIiIlJ9xC3BG5XyvS7QFtgIrEINHYuIiIhUG3Hr4O2S2s3MWgNTgTvLOygRERERKb1SPzzh7t8AvwNuKL9wRERERKSsyvp0bC3CWy5EREREpJqI+5DFCamdCHXwfg28Wt5BiYiIiEjpxX3IYkbKdwdWAy8Bl5VrRCIiIiJSJmroWERERCTLKHETERERyTKx32RhZoOBy4EuhFu0i4Dr3f3ZCopNpNoryZsnREREKkusEjwzGwY8BiwhJHljgaXAY2Z2bsWFJyIiIiIlFbcE73LgUnf/W1K3u81sISHZu6fcIxMRERGRUolbB29nYFaa7jOBjuUXjoiIiIiUVdwE70tgYJruRwLLyi8cERERESmruLdoJwA3m1lP4A3CQxY/A84CRlVQbCIiIiJSCnHbwbvdzFYRGjVOvNViMXCyuz9RUcGJiIiISMnFbibF3R8jPEkrIiIiItVY3GZS+plZvwzdDy3/sERERESktOI+ZDEJaJ6me9Oon4iIiIhUE3ETvD2B99N0/zDqF5uZjTSzpWa2ycwWmtkhxQxfz8yujsb50cy+NLPflGSeIiIiItuSuHXwNgLtCG+vSNYB2Bx3ZmZ2CjAFGAm8Fv2daWZd3P3LDKM9BOwEDAc+A1oDDeLOU0RERGRbEzfBew64zsyOc/f/ApjZDsBfon5xXQpMc/c7o++jzGwQcAFwRerAZnYkcASwm7uviTrnlGB+IiIiItucuLdoRwNtgBwze9XMXiWU5rUlNJ1SLDOrB/QCnk/p9TzQN8No/we8DVxqZv8xs8/M7K9m1jhm3CIiIiLbnLjt4K00sx7AGcC+gAH3Ag+6+4aY89oRqA18k9L9G0IpXTq7EhpU/hEYAmwP3Ey4XXxi6sBmNpxwK5edd945ZlgiIiIi2aUk7eBtAO4sdsAYk0r5bmm6JdSK+p3u7usAzOxC4Dkza+3uBZJFd78DuAOgd+/emaYpIiIiktXi3qItD2uAXMKt3mStKFyql7AS+CqR3EUWR39VRCciIiKSRqUleO6+GVgIDEzpNZDwftt0XgfapdS56xz9XVa+EYqIiIhkh8oswQOYCAw1s2FmtreZTSHUp7sNwMymm9n0pOEfBL4FpppZVzM7mNDMygx3X1XJsYuIiIjUCLHr4JUHd3/YzFoAVxKewP0IONrdE6VxO6cM/72ZHUF4sOJt4L/A48DYyotaREREpGap1AQPwN1vAW7J0K9/mm7/Bo6s4LBEREREskasBC9q1PjPwADCQxEFbu26e9PyD01ERERESiNuCd7dwH6EJkhWkLlZExERERGpYnETvAHAQHd/qyKDEREREZGyi/sU7Srg+4oMRERERETKR9wE73fA1XoHrIiIiEj1F/cW7ZVAJ2CVmS0DtiT3dPfu5RyXiIiIiJRS3ARvRoVGISIl1mnsM4W65Vw3uAoiERGR6iZWgufuf6zoQERERESkfFT2q8pEREREpILFSvDMrJ6Z/dHMPjWzTWaWm/yp6CBFREREJL64JXh/An4J3ATkAb8F/g58C4ysmNBEREREpDTiPmRxMnC+u88yswnAE+6+xMwWAwOB2yssQpEkerBARESkeHFL8FoDi6L/vwe2j/6fBRxZ3kGJiIiISOnFTfC+BNpF/38O/Dz6/yBgY3kHJSIiIiKlFzfBe4zwPlqAKcAfzWwpMA24qwLiEhEREZFSitsO3hVJ/88ws/8AfYFP3f3pigpOREREREquVO3gufs8d5+o5E6keuk09pkCD6KkeyhFRESyX+wEz8yOMrOnzWyRme0UdRtmZgOKG1dEREREKk/cho7PAP4JfAbsAtSNetUGxlRMaCIiIiJSGnFL8MYA57n7JcDWpO7zgH3LPSoRERERKbW4Cd4ewJtpun8PNC2/cERERESkrOImeCuAzmm6HwosKb9wRERERKSs4iZ4dwB/NbODo+87mdkvgRuAWyskMhEREREplbjt4N1gZs2AF4DtgJeBH4EJ7v73CoxPREREREooVoIH4O6/M7M/A10IJX+L3P37CotMREREREoldoIH4O4bgAUVFIuIiIiIlINYCZ6Z1QdGAocBrUipu+fuB5R/aCIiIiJSGnFL8O4EjgGeABYBXmERiYiIiEiZxE3wjgOOd/dXKjIYERERESm7uM2krALWVGQgIiIiIlI+4iZ444C/mFnzigxGRERERMou7i3a54HhwCoz+xrYktzT3Xct78BEREREpHTiJnjTga7AZOAb9JCFiIiISLUVN8EbCBzu7m9VZDAiIiIiUnZx6+B9SXg1mYiIiIhUc3ETvEuAG8xs94oMRkRERETKLu4t2keA+sC/zexHYGtyT3dvWt6BiYiIiEjpxE3wLqzQKESkwnQa+0yhbjnXDa6CSEREpLLESvDc/d6KDkREREREykfcOngiIiIiUkMowRMRERHJMkrwRERERLJM3IcsRKqtdA8RiIiIbMtUgiciIiKSZTKW4JnZPXEn4u7nlk84IiIiIlJWRd2ibZny/VAgD/gw+r4PoQRwbgXEJSIiIiKllDHBc/djE/+b2RXARuAcd/8h6tYIuJufEj4RERERqQbi1sH7DTA+kdwBRP//CRhVEYGJiIiISOnETfAaA+3SdG8LNCy/cERERESkrOImeP8CpprZqWbWKfqcSrhF+2jFhSciIiIiJRW3HbwLgJuAaUDdqNtWQoI3uvzDEhEREZHSipXguftGYKSZ/RbYDTDg8+Q6eSIiIiJSPZS0oeMG0ecTJXciIiIi1VOsEjwzawLcAwwBHNgD+MLMbgO+dvfxFRahiFSKdK98y7lucBVEIiIiZRW3BO96wlO0PQnt4SU8DfyiJDM0s5FmttTMNpnZQjM7JOZ4PzOzrWb2UUnmJyIiIrKtiZvgHQdc7O7vEUrwEhYDu8admZmdAkwB/gLsB7wBzDSznYsZrzkwHZgdd14iIiIi26q4CV5z4Ns03ZsAuSWY36XANHe/090Xu/soYCXhKd2i3A3cC7xZgnmJiIiIbJPiJnhvE0rxEhKleCMIpXDFMrN6QC/g+ZRezwN9ixhvJNAGuCZmrCIiIiLbtLjt4I0DnjOzrtE4l0b/HwAcGnMaOwK1gW9Sun8DHJFuBDPrBvwB6OPuuWYWc1YiIiIi265YJXju/gahlK0esAQYAKwADnL3d0o4T0/5bmm6YWb1gX8Ao919aZwJm9lwM1tgZgtWr15dwrBEREREskPcEjzc/UPgl2WY1xpCfb02Kd1bUbhUD8J7brsQXpE2NepWCzAz2woc7e4Fbve6+x3AHQC9e/culDSKiIiIbAtiJ3gAZtaOkJAVKPmLU4rn7pvNbCEwEHgkqddAwrtuU30FdEvpNjIa/hdATuzARURERLYhcRs63g+4H9iLcEs1mRPq1sUxEbjPzOYDrwPnE9rXuy2az3QAdz/b3bcABdq8M7NVwI/urrbwRERERDKIW4J3B7AcOI9Q965Utz/d/WEzawFcSbgF+xHhVuuyaJAi28OT7KA3JlS9dNtARESyR9wErwuwn7t/WtYZuvstwC0Z+vUvZtzxwPiyxiAiIiKSzeK2g/chhR+OEBEREZFqKG6CNw64wcyOMLPWZrZD8qciAxQRERGRkol7i/bF6O/zFKx/l2jDLu5DFiIiIiJSweImeIdVaBQiIiIiUm5iJXju/kpFByLZSU9rioiIVL7YDR2bWUNgX9I3dPxoOcclIiIiIqUUt6HjI4CHgBZpeqsOnoiIiEg1Evcp2inAM0AHd6+V8lFyJyIiIlKNxL1F2wk4zt1XVGAsIiIiIlIO4pbgvQ7sWZGBiIiIiEj5iFuCdxswwczaEd5qsSW5p7u/U96BiYiIiEjpxE3wZkR/70jTTw9ZiIiIiFQjcRO8XSo0ChEREREpN3EbOl5W0YGIiIiISPnImOCZ2QnAU+6+Jfo/IzV0LGWV7o0XOdcNroJIJJm2i4hIzVRUCd4MoA2wip/q4KWjOngiIiIi1UjGBM/da6X7X0RERESqNyVuIiIiIlkm7rto67v7j9H/7YHhQENCHb25FRifiIiIiJRQkSV4ZranmX0MbDCzd82sCzAfuJSQ5L1kZv9XCXGKiIiISEzF3aKdAKwEjgM+Ap4FngOaAc2B24GxFRmgiIiIiJRMcbdo+wAD3f09M5sLrANucfc8ADO7GZhXwTGKiIiISAkUV4LXAlgB4O7/A34Avkvq/1+gScWEJiIiIiKlEecpWi/mu4iIiIhUI3Geor3fzH6M/t8OuNPMNkTf61dMWCIiIiJSWsUlePemfL8/zTDTyykWkQLSvSZLREREildkgufu51RWICIiIiJSPvQmCxEREZEsowRPREREJMsowRMRERHJMrHeRSsShx6KEBERqR5UgiciIiKSZZTgiYiIiGQZJXgiIiIiWUYJnoiIiEiWUYInIiIikmWU4ImIiIhkGSV4IiIiIllGCZ6IiIhIllGCJyIiIpJllOCJiIiIZBkleCIiIiJZRgmeiIiISJZRgiciIiKSZZTgiYiIiGQZJXgiIiIiWaZOVQcg1Vunsc+k7Z5z3eBKjkRERETiUgmeiIiISJZRgiciIiKSZZTgiYiIiGQZJXgiIiIiWUYJnoiIiEiWUYInIiIikmWU4ImIiIhkGSV4IiIiIlmm0hM8MxtpZkvNbJOZLTSzQ4oY9gQze97MVpvZ/8zsLTM7rjLjFREREalpKvVNFmZ2CjAFGAm8Fv2daWZd3P3LNKP0A14CrgS+A84AHjOz/u7+aiWFLSJJ0r3dRG82ERGpXir7VWWXAtPc/c7o+ygzGwRcAFyROrC7X5TS6Y9mNhj4P0AJnoiIiEgalXaL1szqAb2A51N6PQ/0LcGkmgD/La+4RERERLJNZdbB2xGoDXyT0v0boE2cCZjZr4EOwH0Z+g83swVmtmD16tVliVVERESkxqqKp2g95bul6VaImQ0BbgTOcPdlaSfsfoe793b33i1btix7pCIiIiI1UGUmeGuAXAqX1rWicKleAVFydx9wtrs/WTHhiYiIiGSHSkvw3H0zsBAYmNJrIPBGpvHM7GTgfmCou8+ouAhFREREskNlP0U7EbjPzOYDrwPnA+2A2wDMbDqAu58dfT+VUHI3GphrZonSv83u/l0lxy4iIiJSI1RqgufuD5tZC0K7dm2Bj4Cjk+rUSSe3RAAADwZJREFU7ZwyyvmEGCdHn4RXgP4VG62IiIhIzVTZJXi4+y3ALRn69S/qu4iIiIgUT++iFREREckySvBEREREsowSPBEREZEsowRPREREJMsowRMRERHJMkrwRERERLJMpTeTIiLbrk5jn0nbPee6wZUciYhIdlMJnoiIiEiWUYInIiIikmWU4ImIiIhkGSV4IiIiIllGD1lIqWSqLC/bpnT7gx6cEBGpOirBExEREckySvBEREREsowSPBEREZEsowRPREREJMsowRMRERHJMnqKVkQqhJ60FhGpOirBExEREckySvBEREREsowSPBEREZEsowRPREREJMvoIYttgCq7i4iIbFtUgiciIiKSZZTgiYiIiGQZJXgiIiIiWUYJnoiIiEiW0UMWIlItpXs4KOe6wVUQiYhIzaMSPBEREZEsowRPRERE5P/bu/egu6ryjuPfXxIuBVqwpjR4SYxCIUUtl0jBQQjWaLl0htsYLoONkWQiQg2UirFCA2MREQhRKgFihaggWgdaRAYMkuoIBAJjY0yJWBPABCghifACAQlP/1jrpDsn57w5ed+Ts7N3fp+ZPclee+11nr3eS56svddeNeMEz8zMzKxmnOCZmZmZ1YwnWdSMV62wKur0+9YTL8zMOuMRPDMzM7OacYJnZmZmVjNO8MzMzMxqxgmemZmZWc04wTMzMzOrGc+irTDPmDUzM7NWPIJnZmZmVjNO8MzMzMxqxgmemZmZWc04wTMzMzOrGU+yqAhPqDAzM7NOeQTPzMzMrGac4JmZmZnVjBM8MzMzs5pxgmdmZmZWM55kMUitJj8sv+zYEiIx2z51OgFpMD+X7T7DP+tmtq3yCJ6ZmZlZzTjBMzMzM6sZJ3hmZmZmNeMEz8zMzKxmnOCZmZmZ1UzPZ9FKOgv4B2Av4JfAtIj4aT/1jwSuAvYHVgKXR8TsXsQ6UF5WzKwaPAvezOqqpyN4kiYAs4BLgQOB+4G7JI1sU3808MNc70Dgi8BXJZ3Um4jNzMzMqqfXt2jPA26MiBsi4r8j4hzgaeCTbepPBVZGxDm5/g3ATcD5PYrXzMzMrHJ6luBJ2hE4GLin6dA9wPvbnHZYi/p3A2Ml7dDdCM3MzMzqoZcjeMOBocCzTeXPAiPanDOiTf1huT0zMzMza1LGUmXRtK8WZZur36ocSVOAKXm3T9LSAUVYb8OBVWUHUQPbdD8e9aXG344rM4xO9Kwf9aXN19mSeltadyvbpr8fK8T92B3uRwBp83X610k/jmp3oJcJ3ipgPZuO1u3JpqN0Dc+0qf868Hxz5Yi4Hrh+cGHWm6SFETG27Diqzv3YHe7H7nA/dof7sTvcj90x2H7s2S3aiHgNeAQY33RoPGmWbCsPAB9qUX9hRPy+uxGamZmZ1UOvZ9FeBUyUdKakMZJmAW8BZgNImitpbqH+bOBtkq7O9c8EJgJX9DhuMzMzs8ro6TN4EXGrpDcDnye96HgxcExEPJGrjGyqv0zSMcBM0qtUVgJ/FxHf72HYdeNb2N3hfuwO92N3uB+7w/3YHe7H7hhUPyqiv/kNZmZmZlY1XovWzMzMrGac4JmZmZnVjBO87ZSkKZLuk7RWUkh6R9kxVYWksyQtk7RO0iOSPlB2TFUi6QhJ/yFpRf7em1h2TFUkabqkhyW9IOk5SXdIenfZcVWNpE9JWpT78QVJD0g6tuy4qkzS5/LP9jVlx1Ilkmbkfituzwy0PSd4269dSMvAzSg5jkqRNAGYBVwKHEh6xc9dkkb2e6IV7UaaYPVp4JWSY6myccDXSEs9fpD0ftB5kv64zKAq6LfABcBBwFjgx8Dtkt5balQVJelQYDKwqOxYKmopaRJqY3vPQBvyJIvtnKSxwMPA6IhYXnI42zxJC4BFETG5UPY48G8RMb28yKpJUh9wdkTcWHYsVSdpN+B3wPERcUfZ8VSZpNXA9Ii4ruxYqkTS7sCjpATvImBxRJxdblTVIWkGcHJEdGUk3iN4Zh2StCNwMGnks+ge0iiKWZn+kPQ7fU3ZgVSVpKGSTiGNMrd7Ab+1dz3pP7s/LjuQCntnfnxlmaTvSHrnQBsqYy1as6oaDgxl06X1nmXTFVfMem0W8HPSCkC2BSS9h9RvOwN9wAkR8Ytyo6oWSZOBvYEzyo6lwhaQFnN4jLQs6+eB+yXtHxGbLM+6OR7BqxFJX2jxgGbzNq7sOGug+bkGtSgz6xlJVwGHAydFxPqy46mgpcABwKHAtcBNnrDSOUn7kp5LPj0vS2oDEBF3RcR3I2JRRMwDjiPlaX87kPY8glcvVwPf2kydJ3sRSE2tAtYDI5rK92TTUT2znpA0EzgFOCoiflN2PFWUk5Jf592Fkt4HnAt8oryoKuUw0h2OxZIaZUOBIyRNBXaNiFfLCq6qIqJP0i+BfQZyvhO8GomIVaQkxLaCiHhN0iPAeOB7hUPjAS+fZz2X1/M+BRgXEY+VHU+NDAF2KjuICrkdWNhU9g3gcdLInkf1BkDSzsB+wH0DOd8J3nZK0gjSSNSf5aI/l7QH8GRErC4vsm3eVcA3JT0E/AyYCrwFmF1qVBWSZ3vunXeHACMlHQCsjgiPMHdI0r+Qnnc6HliTf6YB+iKir7zIqkXSZcCdwFOkiSqnkV5B43fhdSgi1gJri2WSXiL9TC8uJ6rqkXQFcAfpTtuewIXArsBNA2nPCd72ayrwT4X9O/OfHwdu7Hk0FRERt0p6M+nh171I73M7JiKeKDeyShnLxv8jvThvN5EeMLbOnJX/vLep/GL8fsstMYL0aMsI0mtmFgFHR8TdpUZl26O3AbeQbnc/BzwIHDrQf1/8HjwzMzOzmvEsWjMzM7OacYJnZmZmVjNO8MzMzMxqxgmemZmZWc04wTMzMzOrGSd4ZmZmZjXjBM/MakVSn6SJW1B/oqTt8sXAkuZLumYrtX2yJL+Hy6wkTvDMbCOSbpQUkua0OHZ5PvaDMmLrBUkzJG0vb98/EZje2JG0XNL5JcZjZl3iBM/MWnkKmCBp10aBpGGkpbG8nFhNRMTqiHix7DjMrPuc4JlZK4tIC4V/tFB2LLAOmN9cWdLHJS2RtE7SrySdK2lI4fh5khZJeknSCklz8trHjeMT863Vv5K0ONe7T9Lo/oKUtHe+zbhO0lJJx7Wo81ZJ35G0Jm93StqnTXsTSUv47Z9HKqNxu3dz19CmveWSLsqjoi9KekrSBEl75Jj6JD0u6cOFc4ZK+rqkZZJeycc/09SfwyTNLFzTTEnXSppfqDNf0tckXSpplaT/lXRFUzsbbtHmc0cBX25ce/Fr03Rd43Kd4YWyj0l6QtLLeYT3T1v0x99IeiR/vZZJ+mdJO/bXh2Y2ME7wzKydrwOTCvuTgG8AGz1XJWkycClwETAG+HvgAv5/rVSAN4BpwP6kxdwPAb7a9Hk7kW4XTgIOA/YAZrcLLicqt5F+jx2Wz5uR22nU2YW07u064Mhc72lgXj7W7FbgSmApaa3hvXJZp9fQyjTgIeAg4LukNXdvBn4IHAD8BPiWpJ1z/SHAClJyPQb4R+BzpHWiG84nrdt7JnBoPue0Fp99OvA68H7g7BzLhDZxngj8FrikcO0dkfSXpDWsr8/XdEdup1jnI8C3gWtIfTgJOJn0vWNm3RYR3rx587ZhI/1D/QPgTcArwD6khdhfBUY2jhfqPwmc0dTGNGBJP5/x17m9IXl/Iilx3LdQ53TgtUadFm18GFgPjCyUHZ7bmZj3J5FGIlWoMxR4Hvho4bP7CsdnAIs76KeNrqFNneXALYX93XJ8XymUvSOXje2nncuAeYX9p4HPFvYFPAbML5TNBx5oaudHwJymOtc0xXt+0zkb9U8uG5djHp73bwZ+1FRnTvonZsP+T4ALm+ocD/QVvz7evHnrzjYMM7MWImKNpNtISdJaUvLwpKQNdST9CfB24DpJ1xZOH0ZKOhr1PkganRsD7E5KsnYkJY4rc7VXI2JpoY2VwA6kkbzVLUIcA6yIiOIzgQtII20NBwOjgReLcQO7AO/q7/qbdXgNrSxq/CUi+iS9DPyicPzZ/Oeehc+aShqdGwX8AakfnsjHds+f+VCh3ZD0MOlr0fKzs5XFz+miMaRRu6IHgE8U9g8GDpF0QaFsCOn6RpCSVjPrEid4ZtaffyXdUuwj3YJt1njMYypwf6sGJI0C7gRuyG08T7pdeQspQWp4venUxq3gdo+SqE15c3w/B05pcaxV0tj6gzq/hlZ+37QfTWUbXaekCcDVpNuw9wMvAJ8CTmjRzua0+uwtfTTnDTbt6x2a9jv9WlwMfK/Fsee2MCYz2wwneGbWn3tJt0mHA7c3H4yIZyWtAN4VEXPbtDGWlASdGxHrAVpNhhiAJcBbJb09Ip7KZYewcQLzKHAqsCoi1nbY7muk0bmirXUNrRwOLIiIDe+nk7RhtDEififpGdK13pePC3gf8MwgP7vVtT8H7CLpjyLihVx2QFOdJaRnAYua9x8F9ouIXw8yRjPrgCdZmFlbERHAe4HREfFqm2ozgM/kmbP7Snp3nlHZeL/a46TfNdMkjZZ0KukZvcGaR3rubK6kAyQdBsxk45HAb5Nugf67pCPz5x8h6cp2M2lJz6GNknSQpOGSdtqK19DKr4CDJB0taR9JF5ImiBTNIvX5CZL2JU0M2YvORvX6sxz4QJ553JghuwB4Cfii0qzlk9h4Ag3AV4APSZqeY57MpiOOlwCnSbokf4/sp/Qy5MsHGbOZteAEz8z6FREvFkZuWh2fQ3pO7wzgv4CfAlOAZfn4IuDTwHmkkZ4zSbcfBxvXG6QkYggpCZkLfIE08aFR52XgCOA3pFuDj5FuOb8JWNOm6e+TZrjeSxq9OnVrXUMb15Fm294MPEyahHFlU50rgG+SZjU/mMtuI80WHoyLSM/x/Q/5tmlErCZNeBlPenZwCnBh8aSIeJD0vN0nSc/9nUhK/It17ia9auco0vODDwGfxe9VNNsqlP6DbmZmVSbpUeBnEXFO2bGYWfn8DJ6ZWcXkSR8fAf6T9Ht8CvAX+U8zMyd4ZmYV9AbwMeDLpFvUS4CjI2JhqVGZ2TbDt2jNzMzMasaTLMzMzMxqxgmemZmZWc04wTMzMzOrGSd4ZmZmZjXjBM/MzMysZpzgmZmZmdXM/wF4F/DAdgv+GwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 648x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "std = np.sqrt(np.var(shuffled))\n",
    "\n",
    "plt.rcParams[\"font.size\"] = 14\n",
    "fig = plt.figure(figsize=(9, 6))\n",
    "ax = fig.add_subplot(111)\n",
    "\n",
    "# the histogram of the data\n",
    "n, bins, patches = ax.hist(shuffled, bins = 50, density=1)\n",
    "\n",
    "ax.axvline(x = np.abs(delta), linewidth=2, color='r', label = 'A/B means difference')\n",
    "ax.axvline(x = std, linewidth=2, color='y', label = 'Standard Deviation')\n",
    "ax.set_xlabel('Mean delta magnitude')\n",
    "ax.set_ylabel('Binned mean counts')\n",
    "ax.set_title(r'Histogram of the means delta for random permutation sets')\n",
    "\n",
    "plt.figlegend(loc = (0.5, 0.8))\n",
    "# Tweak spacing to prevent clipping of ylabel\n",
    "fig.tight_layout()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Yellow line represents the Standard Deviation (usually referred to as 'sigma') of 'shuffled' and the red line is the difference between 'A' and 'B' means. We conclude that 'A' and 'B' are at least 8 sigma away and are extremely unlikely to represent the samples of the same kind of flower."
   ]
  },
  {
   "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}