exploration- Effect Size and Statistical Power.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from pylab import*\n",
    "import math\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def my_pdf(x, mu, sd):\n",
    "    \"\"\"Normal distribution density \"\"\"\n",
    "    return 1 / (sd * math.sqrt(2 * math.pi)) * math.exp( - (x - mu)**2 / (2 * sd**2)) \n",
    "\n",
    "def my_cdf(x, mu, sd):\n",
    "    \"\"\"Cumulative density function for normal distribution \"\"\"\n",
    "    return (1.0 + math.erf((x - mu) / (sd * math.sqrt(2.0)))) / 2.0\n",
    "\n",
    "def phi(x):\n",
    "    \"\"\"Cumulative density function for the standard normal distribution \"\"\"\n",
    "    return (1.0 + math.erf(x / math.sqrt(2.0))) / 2.0\n",
    "\n",
    "def beta_two_tail(m1, s1, n1, m2, s2, n2):\n",
    "    \"\"\" beta-error two-tailed for alpha = 0.05\"\"\"\n",
    "    z025 = 1.96\n",
    "    z = (m1 - m2) / math.sqrt(s1**2  / n1 + s2**2 / n2)\n",
    "    return phi(z025 - z) - phi(-z025 - z)\n",
    "\n",
    "def power(m1, s1, n1, m2, s2, n2):\n",
    "    \"\"\"Statistical power (1 - beta error)\"\"\"\n",
    "    return 1 - beta_two_tail(m1, s1, n1, m2, s2, n2)\n",
    "\n",
    "def sd_pooled(m1, s1, n1, m2, s2, n2):\n",
    "    \"\"\" Pooled standard deviation\"\"\"\n",
    "    return math.sqrt(((n2 - 1) * s2 ** 2 + (n1 - 1) * s1 ** 2) / (n1 + n2 - 2))\n",
    "\n",
    "def hedges_biased(m1, s1, n1, m2, s2, n2):\n",
    "    \"\"\"Biased effect size\"\"\"\n",
    "    return (m2 - m1) / sd_pooled(m1, s1, n1, m2, s2, n2)\n",
    "\n",
    "def hedges_unbiased(m1, s1, n1, m2, s2, n2):\n",
    "    \"\"\"\"\"\"\n",
    "    correction = 1.0 - 3.0 / (4 * (n1 + n2 - 2) - 1)\n",
    "    return correction * hedges_biased(m1, s1, n1, m2, s2, n2)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu1, sd1 = 0.0, 1.0 # mean and SD of first sample\n",
    "mu2, sd2 = 1.0, 1.0 # mean and SD of second sample\n",
    "ns1, ns2 = 16, 16    #number of obs per sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "xmin = mu1 - 4 * sd1\n",
    "xmax = mu1 + 4 * sd1\n",
    "increase = (xmax - xmin) / 100\n",
    "x = np.arange(xmin, xmax, increase)\n",
    "pdf1 = [my_pdf(x1, mu1, sd1) for x1 in x]\n",
    "cdf1 = [my_cdf(x1, mu1, sd1) for x1 in x]\n",
    "pdf2 = [my_pdf(x1, mu2, sd2) for x1 in x]\n",
    "cdf2 = [my_cdf(x1, mu2, sd2) for x1 in x]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztnXd4VOXyxz+T0KU3UYogTUFFEFTEggqCvfeu196u/SKKCOjPjgUVGxbs14peFSsiFgSUjihdmvROIGV+f8wmWWLKSdjNOSHzeZ73yZ7+3c3umfPOO++MqCqO4ziOA5AStgDHcRwnOrhRcBzHcXJwo+A4juPk4EbBcRzHycGNguM4jpODGwXHcRwnBzcKjuM4Tg5uFBzHcZwc3Cg4juM4OVQIW0BxqV+/vjZv3jxsGY7jOGWKCRMmrFDVBkXtV+aMQvPmzRk/fnzYMhzHccoUIjI/yH7uPnIcx3FycKPgOI7j5OBGwXEcx8nBjYLjOI6TgxsFx3EcJ4ekGQURGSYiy0RkagHbRUSeEJFZIjJZRDolS4vjOI4TjGT2FF4Gehey/WigdaxdDjyTRC2O4zhOAJI2T0FVR4tI80J2ORF4Va0e6M8iUltEdlHVJcnS5DiOE09WFmzeDJs22d+qVaFBA1v/00+wdau19HTIyIBWrWCvvWDLFnjzTcjMzG1ZWbD//tClC6xbB88/D6q2PivLXvfoYduXL4ehQ20d5P49/njo1AkW/qW8NHQLFbZuokJGGhUy0mjZvgonXLlr0j+TMCevNQb+ilteGFv3D6MgIpdjvQmaNWtWKuIcxyk7bN4MixbZzbdNG1v34IOweDGsWmVt3Tro3h0GDLDtderAmjXbnufKK+GZZ+wmffDB/7zOrbfaeTdvhosv/uf2AQPspr92LdxyC6SQSUOW0YilNGIpu6xvSZcubVg9dRGN+vWnNmtyWk3WsWjJQBh6Bmu/GMtd93Xd5tw/NjsLrnxz+z+sIijSKIhIqqpmJuHaks86zW9HVX0OeA6gc+fO+e7jOM6Oi6rd9Fetgn32sXVXXQU//wzz58Pq1bauVy/4/HN7/cwztn+dOtZq17aeQDbXXQepqVCtmq2vWhXat7dtqakwciRUqgQVK+b+3Xln216zJsyZAxUqQGraRirNnkHFuX9QuW1z4CAaV1pOZrMuyOJFSEZGzjUzd7oP6EPr3TNp3eiTXGG1GkDN3WlzWj0A2h+zG9x7r4mrVg2qVOGgli2T9vnGE6SnMEtE3gVeUtXpCbz2QqBp3HITYHECz+84Thnmo4/gyy9h8mRra9fCnnvC9NhdKC0Ndt0VunaFxo2t7bFH7vGzZtnNvSCyewwFcdRR+azcuhWoRIpm0uLm02HiRJg3L9f/c/nlcMRBpNSrA4ceAk2bQpMmsMsu0KgRqa1bAyC7NYMlhXjKd9kF7rijcIFJIohR2Ac4C3hBRFKAYcBbqrpuO689ArhWRN4CDgDW+niC45Q/MjLgt9/gm29g/Hh45x0QgQ8/hPfes57BOefA3ntDu3a5x730UuHnLcwgBGbpUrNM330H339vvqmPP7aTp6WZr+iii6yLsccesPvudlyFCjB8eAIElD6iGtwbIyKHAm8CtYF3gYGqOquAfd8EugP1gb+Bu4GKAKo6VEQEGIJFKG0CLlbVIjPdde7cWT0hnuOUfX75BR5/HD77LNf9064djB4N9erB+vWw006QEtZsqgsvhFdftdd16tggQ+/ecPXVIQnaPkRkgqp2Lmq/QGMKwLHAxUBz4BHgdeAQ4FOgTX7HqerZhZ03FnV0TVHXdxxnx2DTJvjkE+jc2R6oly6FL76AE06we2337tCoUe7+NWqUorjJk80AfPYZjBtnfvwjjrCn/969oUOHEK1T6RLEffQn8C3wkKr+GLf+3VjPwXEcp0CmTIFnnzVvyrp1MGgQ9O0LxxxjhiEhbp6SsH69GYJhw+DXX20kuVcvWLECmjWznkI5JIhRuEBVx8SvEJFuqvqDql6fJF2O45RxMjOhZ0/49luoXBlOP93COA87zLZXCCsgPivLnvpnzoRrr4WOHeGJJ2zgol69kERFhyD/lieAvCkonsxnneM45ZyMDPjqK/O4pKZaZNBxx9lDd+j325kzYeBAiz19/nnzY02dmhuH6gCFGAUR6QocBDQQkZviNtUEwurwOY4TQVThgw/MLfT77+ai33tvC7UPnblzoV8/eOMNqFIFboq7nblB+AeFjZxUAqpjhqNGXFsHnJZ8aY7jlAXGjrUewamnWijpBx9YKohI8M47Nrnhvffg5pvNQAwcGLaqSFNgT0FVvwO+E5GXVTVQbU/HccoXGzaYq6haNXjxRbjgghDHCrJRNWE1api1OuccMwSNG4csrGxQmPvoMVX9NzBERP4xmUFVT0iqMsdxIss338Dhh0P16jBiBOy7bymHkBbEwoWW/2LjRvj6a5tRPGxY2KrKFIXZ9OzpeA+XhhDHcaLPqlU2d+vtt62dcQYcckjYqrDewSuvwPXXW9jToEG2TvJLseYURmHuowmxv99lrxOROkBTVZ1cCtocx4kQ33wD551naZ8HDYJTTglbUYx16yy96ZtvWrzrsGG56SacYlPkFD0RGSUiNUWkLjAJeElEHk2+NMdxosKTT9qcg1q1LD1F374RGDvIRtWSJg0aZC4jNwjbRZB/ay1VXSci/8Iypd4tIt5TcJxyxJ57wtlnW2GY6tXDVhPj00/hyCPNUk2ebOGmznYTJJlHBRHZBTgD+CTJehzHiQhz5uSO0fboAa+9FhGDkJkJt90Gxx5rXRhwg5BAgvQUBgAjgTGqOk5EdsfyITmOs4Py449w4on2+uSTLUloJNi0ybosI0bYiPcNN4StaIejyJ6Cqv5XVfdR1atjy3NU9dTkS3McJwzee88ShNauDT/8ECGDsGqVVb75+GMYMgSeesqS2DkJJUjq7AbAZVja7Jz9VfWS5MlyHCcMXnkFLrkEDjjAHsbr1w9bURyLF8Off9os5dM8qUKyCOI++gj4HvgKSEatZsdxIsLatdZL+PBDK3ATCdats6LIe+1lAx2REbZjUmTlNRGZqKr7lpKeIvHKa46TeP7+O7cofWZmiDUO8jJ/vlXfufZay13klJigldeCRB99IiLHJECT4zgR5JlnoHVri+qECBqE1atzizA4SSeIUbgBMwxpIrJORNaLyLpkC3McJ/m88ooF8XTvbnMRIsPSpTYHYfVqK9DQucgHXCdBFDmmoKpRSHPlOE6CefddG1Tu0cPGbiMTyJOebrU6lyyxGcpuEEqVINFHApwLtFDVgSLSFNhFVX9JujrHcZLCxIlw7rlw4IE2qBypuV8VK9r8g513NoFOqRLEffQ00BU4J7a8AXgqaYocx0k67dvbpOARIyIUzJOVBdOn2+sLL7RCDU6pE8QoHKCq1wBpAKq6GqvK5jhOGWPZMmsVK1rdmdDrJsfTrx907AgzZoStpFwTZJ5CuoikAgo5k9mykqrKcZyEs2WLpa7YsMHcR5GJMgJ4/XUr6Pyvf8Eee4StplwTxCg8AXwANBSRe7H6zHcmVZXjOAlF1aKMfv7ZBpgjZRB++gkuvdRCoJ56ygvjhEyQ6KPXRWQCcCQgwEmq6v07xylDPPOMZTy98044NUqZy/7+26r1NGli1qqSe6bDprAazXXjFpcBb8ZvU9VVyRTmOE5iGDPGgnmOOw7uuSdsNXmoV89cRmedFbEBjvJLYT2FCdg4ggDNgNWx17WBBUCLpKtzHGe7adsWLrgAHn0UUoKElpQWGzda6NPAgWErceIo8Cuiqi1UdXeslsLxqlpfVesBxwHvl5ZAx3FKRmamtQYN4MUXrUBZZHj3XcutMXNm2EqcPAR5buiiqp9mL6jqZ4AnInGciPN//2eZIjZuDFtJHubMsanUzZpBC3c4RI0gRmGFiNwpIs1FZDcR6QusDHJyEektIjNFZJaI/Cef7c1E5FsR+U1EJnviPcdJDKNHw9132/httWphq4kjI8OmUqekwNtv+8ByBAliFM4GGmBhqR/EXp9d1EGxuQ1PAUcD7YCzRaRdnt3uBN5R1Y7AWdjsacdxtoNVq+Ccc6BlS4s6ilSE54ABFhc7dCjstlvYapx8CBKSugrLlFpc9gdmqeocABF5CzgRmB5/eqBm7HUtYHEJruM4ThzXXGORnj//DDWilM4yKwsmTYKLLrJoIyeSBJm8VlIaA3/FLS8EDsizT3/gCxG5DtgJ6JFEPY6zw7NihRmDu++G/fYLW00eUlIs+97WrWErcQohmQFq+XVa85Z5Oxt4WVWbAMcAw0XkH5pE5HIRGS8i45cvX54EqY6zY1C/vhXL+c8/RvBC5rHHYN4882VVrhy2GqcQkmkUFgJN45ab8E/30KXAOwCq+hNQBfhHqXBVfU5VO6tq5wYNGiRJruOUXVTNTZ+WZi6jCsn0ARSXzz6DG2+EF14IW4kTgCKNgoi0EZGvRWRqbHkfEQmS+2gc0FpEWohIJWwgeUSefRZg6TMQkT0xo+BdAccpJi+8AFddZeH/kWL1apux3L493HVX2GqcAATpKTwP9AHSAVR1MnaDLxRVzQCuxSa/zcCijKaJyAAROSG2283AZSIyCUujcZGq5nUxOY5TCIsWwS23wOGHW7RnpLjxRhv1fvlldxuVEYJ0Mqup6i+ybVxbRpCTxya9fZpnXb+419OBbkHO5TjOP1G1HkJ6Ojz3XMTCT7/4wopA9+3rJTXLEEGMwgoRaUluPYXTgCVJVeU4TiDefhs+/hgefhhatQpbTR4OOAD694/gqLdTGEGMwjXAc8AeIrIImIvVbHYcJ2T23huuuAL+/e+wleRB1ZIt3X132EqcYhJkTGG+qvbAZjLvoaoHq+r8JOtyHCcA7dtb1FGkiub89puV1fz997CVOCUgiFGYKyLPAQcCG5Ksx3GcAPz4I5x3nqW0iBSZmXD55bB0Key8c9hqnBIQxH3UFjgecyO9KCKfAG+p6pikKnMcJ1/S0+HKK2HNmgjmkxs6FMaPhzffhDp1wlbjlIAiewqqullV31HVU4COWK6i75KuzHGcfHniCZgyxf5Wrx62mjiWL7d6n0ccAWeeGbYap4QEmtEsIoeJyNPAr9gEszOSqspxnHxZtMjGbo87Dk48MWw1eXjiCdiwAYYMiVhsrFMcinQfichcYCKWjuJWVY1ayQ7HKTf062clCZ54IoL33X79oFcv2HPPsJU420GQMYUOqrou6UocxymSe++1XkKkCpZlZloPoVYtOPjgsNU420mBRkFEblPVB4F7ReQfqSdU9fqkKnMcJ4esLOsZNGoEJ58ctpo8DBtmYwm//OKFc3YACuspzIj9HV8aQhzHKZgXX7SMESNGQN26YauJY+1aS2PRtq3VXHbKPAUaBVX9OPZyk6r+N36biJyeVFWO4+SwZg3ccYe56iMX5XnffVbZ57PPIjjI4ZSEINFHfQKucxwnCdx7L6xcCY8/HrH77pw5VjznwgsjWObNKSmFjSkcjVVDaywiT8RtqknALKmO42wfc+dapNFFF1nmiEjx1ltQsaJZLWeHobCewmJsPCENmBDXRgC9ki/NcZyHHrIqaoMGha0kH+64w2bR7bpr2EqcBCJF1bQRkYqqml5Keoqkc+fOOn68j3075YPNmy1rxCGHhK0kjqwsWLIEGjcOW4lTDERkgqoWWdgiyJhCcxF5V0Smi8ic7JYAjY7jFIAqbNkCVatGzCAA/Pe/0LIlTJgQthInCQQxCi8Bz2DjCIcDrwLDkynKcco7771n0UZz54atJA9bt5rbqG1b2HffsNU4SSCIUaiqql9jrqb5qtofOCK5shyn/JKebvfdatUiGPo/dKhFHT34YMSKODiJIkiaizQRSQH+FJFrgUVAw+TKcpzyy7Bh8Oef8NFHEbvvrl0LAwZAjx5w1FFhq3GSRJCewr+BasD1wH7A+cCFyRTlOOWVjRutrHG3bnD88WGrycPXX5theOCBiE2YcBJJkT0FVR0Xe7kBuDi5chynfPP661a07N13I3jfPeUUWLAAdtklbCVOEils8trHQIHxqqp6QlIUOU455l//gj32sJ5CpFi4EJo0cYNQDiisp/BwqalwHIf0dJsgfOihYSvJw9y5Zqkef9zqgDo7NIUlxPOSm45TSixebOmDnn0WTohaH7x/f0hJieAgh5MMglZey6+ewu5JUeQ45ZB777Vko3vtFbaSPEybBsOHwy23+AzmckKQkNT4adFVgNOBKGV0d5wyzbx58PzzcOmlsHvUHrX69YMaNeD228NW4pQSRYakqurKuLZIVR/DJ685TsK45x7zztx5Z9hK8rB8OXz1Fdx0E9SrF7Yap5QI4j7qFLeYgvUcaiRNkeOUIxYtMu/MDTdYcE+kaNDAZi9XqhS2EqcUCeI+eiTudQYwDzgjKWocp5zRuDH8+CO0aBG2kjysXg21a3sPoRwSZPLa4aUhxHHKG6o2QW3//cNWkg9nnGHxsZ9+GrYSp5QpckxBRGqLyPUi8qiIPJHdgpxcRHqLyEwRmSUi/ylgnzNiabmnicgbxX0DjlNWOfdcuPHGsFXkw+jRNpbQo0fYSpwQCJL76FOgOTCFbSuwFYqIpAJPAUcD7YCzRaRdnn1aY/Weu6lqeyzPkuPs8EyeDG++aYE9kaNfP5u5fNVVYStxQiDImEIVVb2pBOfeH5ilqnMAROQt4ERgetw+lwFPqepqAFVdVoLrOE6Zo39/qFUrgj2FUaPgu++sMHTVqmGrcUIgSE9huIhcJiK7iEjd7BbguMbAX3HLC2Pr4mkDtBGRH0TkZxHpnd+JRORyERkvIuOXL18e4NKOE11++w0++MAMQp06YavJw9ChVnP5ssvCVuKERJCewlbgIaAvuTObFShqmk1+OR7zzoyuALQGugNNgO9FZC9VXbPNQarPAc+B1WgOoNlxIsugQRbY8+8oOktfecWKOVSpErYSJySCGIWbgFaquqKY514INI1bbgIszmefn1U1HZgrIjMxIzEOx9lBefRRyx5Rq1bYSuJQtYx8lStHMNeGU5oEcR9NAzaV4NzjgNYi0kJEKgFnASPy7PMhVvcZEamPuZPmlOBajlNm2G03OOaYsFXk4dtvoWVLmDIlbCVOyAQxCpnARBF5tjghqaqaAVwLjARmAO+o6jQRGSAi2XkgRwIrRWQ68C1wq6quLNlbcZxoM3Ei9OplmagjhaqNfGdlQevWYatxQiaI++jDWCs2qvopFtIav65f3GvF3FMliW5ynDLFgAEwdmwEB5e/+w6+/94ijnwsodwTZEbzK6UhxHF2ZCZNsoiju++2QeZIcc890KiRlX1zyj1eT8FxSoGBA6FmTUt8FynGjbO5CYMH+7wEB/B6Co6TdKZOhffeg7vuiqDrqHNn+N//oHv3sJU4EcHrKThOkmna1CqrRW5eQnZGvmOOgWrVwlbjRASvp+A4SaZWLbjjjrBV5MMZZ0CXLnDbbWErcSJESeopzMXrKThOIPr3hw4d4OSTw1aSh19+gXffNfeR48Th9RQcJ0nMmmUDzDfeGEGjMGgQ1K0LV18dthInYgSpp3CfiNSOW64jIoOSK8txyj733WeVLG+5JWwleZg4ET7+2KxVJHN3O2ESZEbz0fEJ6mJprqM2Sd9xIsW8eVZ7+fLLbQpApBg0yAY6rrsubCVOBAkyppAqIpVVdQuAiFQFKidXluOUbe6/H1JSIjqG27evDTJHKiOfExWCGIXXgK9F5CVsEtslgM9ydpxCOOQQC0VtnLeCSBTo2NGa4+RDkIHmB0VkMtADq5EwUFVHJl2Z45Rhzj03bAX5MHOmuY7uvz+i1sqJAkHGFFDVz1X1FlW92Q2C4xTM0qVWL2HjxrCV5MN999nU6ooVw1biRJhARsFxnGA88gjceissWRK2kjzMmQOvvw5XXgkNG4atxokwbhQcJ0GsWAHPPANnnw2tWoWtJg//939QoUIE42OdqFGgURCR50TkZBHxQGbHCcDgwbBpkwX3RIoFC6z28qWXwq67hq3GiTiF9RSGAR2AT0XkaxG5XUQ6lJIuxylTrF4NTz4Jp50Ge+4Ztpo8VKkCV10Ft98ethKnDFBg9JGq/gz8DPQXkXrAUcDNIrI38Bvwuaq+UzoyHSfarFhhUZ533hm2knxo2BAefzxsFU4ZIcg8BWJ1k9+MNURkP6B3EnU5TpmidWurahk5nnsO2rTxeglOYEo00KyqE1T13kSLcZyyyNdfRzDaCGDZMivi8NJLYStxyhAefeQ428HGjXDWWXDFFWEryYdHH4W0tIgWc3CiihsFx9kOhg618YQ+fcJWkoeVK2HIEDjzTGjbNmw1ThkiSOrs8SJyjYhErbqs44TKpk3w0EPQowd07Rq2mjw89ph1YyIXH+tEnSA9hbOAXYFxIvKWiPQSEUmyLseJPM8/D3//DXfdFbaSfNhlF/Np7bVX2EqcMkaQhHizgL4ichdwHDZ/IUtEhgGPq+qqJGt0nEgybx4cfjgcemjYSvLBK6o5JSTQmIKI7IPVan4IeA84DVgHfJM8aY4TbQYPhpFRSw+5dq3lOMrICFuJU0YJMqYwARgMjAP2UdXrVXWsqj4CzEm2QMeJGlu2wNSp9jpyCUeHDIHzzssV6DjFJEhP4XRVPVJV34irvtYCQFVPSao6x4kgL70Ee+9tpY4jxfr1FoZ63HGw775hq3HKKEGMwrsB1znODs/WrZZw9MADoUPUMoE99RSsWgX9+oWtxCnDFDjQLCJ7AO2BWiIS3yOoCVRJtjDHiSIvv2xJR597DiIVg7dhgxVzOPpo6NIlbDVOGaawnkJbLNqoNnB8XOsEXBbk5CLSW0RmisgsEflPIfudJiIqIp2DS3ec0mXrVitedsABcNRRYavJw4IF0KCB9xKc7aawLKkfAR+JSFdV/am4JxaRVOApoCewEJvnMEJVp+fZrwZwPTC2uNdwnNJk4kRYvtxmMUeqlwDQrp0NLqd4kgJn+yjMfXSbqj4InCMiZ+fdrqrXF3Hu/YFZqjondr63gBOB6Xn2Gwg8CHhJKCfS7L+/PZDXrRu2kjz89JNNUqvh9bCc7aewyWszYn/Hl/DcjYG/4pYXAgfE7yAiHYGmqvqJiLhRcCLL4sU2SbhevW3XL1gAo0fD7NnWi1i+HFTNk9OgAbRsCYccAs2bJ0nY+vUWbdSzJ7z1Vu76Vatg0iSYPh1mzIB163K31a1rPYt27Wy03I2JE0dh7qOPY39fyV4nIilAdVVdV9BxceTXwdY85xoMXFTkiUQuBy4HaNasWYBLO07i2LLFoo2OPRaefhrGjYMXX7SJa/Pn5+5Xp44ZAhEzDqvi5vo3a2bjEJdeamMSCXM/DRliF7rpJsvf/f778N57VtwhK8v2qVkz15qpWkrtTZtsuVIlMyinngonnhjBbpBT6qhqoQ14A4s42gn4HVgC3BrguK7AyLjlPkCfuOVawApgXqylAYuBzoWdd7/99lPHKU2efloVVG+8UbVjR3tdrZrqKaeoPv646sSJqlu3/vO49HTVyZNVn3xS9bTTVHfayY7dZx/VZ5/N/5hisXatat26qocdpnrRRaoVKtgF9thDtW9f1ZEjVRcuVM3K2va4zEzVuXNVP/nE3tRuu9lxVauq3nCD6oIF2ynMiSLAeC3ivq32TSjy5j4x9vdc4FGgIjA5wHEVsBnPLYBKwCSgfSH7jyrKIKgbBaeU2bRJtV49u1+CaocOZiTWrCn+udauVR06NNewtGmj+v77/7xnB+bmm+1E2Tf0669XnTat+OfJylIdN071wgtVU1NVK1ZUveIK1RUrSijMiSKJNArTYobgv8BhsXWTAp0cjgH+AGYDfWPrBgAn5LOvGwUnUsyapdq6tf1KmjRR/eij7biBx5GVZQ/p7drZuQ85RHXmzGKcID1d9ZFH7Aaemqrap4/q339vvzBV60FcfbX1OurXV33llcS8aSd0EmkUrgcWAZ9i4wS7Ad8HOXkymhsFJ9lkZZl7p1o1u+e2bq26ZUvir5Oebj2H2rXtWk8/HeD+O2NGblfj2GNtORlMnqzatatd54gjzA3llGkSZhTyPQgqlOS4RDQ3Ck4yWblS9bjjcu+Ff/yhOn9+cq+5cKHqUUfZNXv3Vl2+vIAdX3/dBibq1g1oQbaTzEzVZ56xazZooPrFF8m9npNUghqFIFlSK4vIOSJyh4j0E5F+gBd9dXY4pk+3uQgjR8KDD8L//getW1vkUDJp3Bg+/9xSF337rWnYJsnpli1w1VVw7rnQsSOccQbcfLPVAU0mKSlw5ZUWbtWwIfTqBf3750Y1OTskQaY/foRNOssANsY1x9lh+OQTCzvdsAFGjYLNm2HPPa08QWkgYnVxRo+GtDQr7/nRR8CaNdC7t02jvvVWeOMNeOUVCyFt0KB0xO25J4wdC+efD/fcA+ecYyKdHZOiuhLA1CBdjtJq7j5yEs0zz6iKqHbqZNGYK1eq1qypetJJ4ehZuFC1SxfVpszXFbu0t2ig4cNtY/Yg8J9/lr6wrCzVBx/UnNHxlStLX4NTYkiU+wj4UUT2Tq5pcpzSR9XSYF91FRxzDHz/PTRtaq6j9eth4MBwdDVuDKOfmcbEygeSumQhr547Ej33PJgzx9KzXnoptGpV+sJErLfy5pvWc+jWDRYuLH0dTlIJYhQOBibEsp1OFpEpIjI52cIcJ5mowu23wx13mDfkgw+gWjWbFPzEE7YutJr3kyZRpXd36tQTHjx+DBe+fDg33QRZ34yCKlXgrrtCEhbjrLPgyy8t98dhh207rdsp8xSW+yibo5OuwnFKEVXLCvHYY+bHf/LJ3OSiL70E6ek2nhoKv/5qaSeqVUO+/ZZBu7di442mdevVlzBk/slI3TohiYvj0EPNMPTqZYbh22+hRYuwVTkJoMiegqrOB5oCR8RebwpynONEEVXo08dusjfcYKmD4rNN9+kDv/wSjneGX3+FI4+E6tUtd1GrVqSkmNb7L/mDp5+GmwfVQbXoU5UK++8PX39tyfYOOwzmzg1bkZMAgoSk3g3cjuUuApvd/FoyRTlOsujfHx54wMYRBg/eNjHdhg223LFjCMJ+/92eumvVMoMdOwStAAAgAElEQVSw++45m+SnH7l9WFuG9X6HwYPN5RUZw9CpE3zzjX14PXua/80p0wR54j8ZOIFYGKqqLgY8165T5nj4YRgwwMZphwzZ1iCMG2cDvN99F4Kw+fPthpqaai6Z+Dzbqja426gRF71zDFdcAfffbwPkkWHffeGzz2DpUjNs8elhnTJHEKOwNRbOpAAislNyJTlO4hk+3O6tZ5wBzz67rcsoe9C5cuUQegl//w09etiT9hdf2Gy5eD78EH78EQYMQGpU5+mnbRC8b18YNqyUtRbGAQfYxIqZMy3H+IYNYStySkpRMatYRbRnsYynlwE/AdcFiXdNRvN5Ck5x+ewzC+0/4gjVtLR/bv/0Uwu9f/zxUha2fr3qfvtZ4qMffvjn9q1bLZXqnntaoqQYW7ZYWozUVNURI0pRbxA++EA1JUX1mGO20eyED4nMfYTVWX4IeBjoGeSYZDU3Ck5xGDfO7rkdO1rq6rykp1u20pYtk5P0rkDS0y2hXUqK6scf57/Pb7/ZLLp87vzr16t27qxapYrqjz8mWWtxGTrUbi2XX+4ZViNEUKMQJCQVVf0S+DLh3RTHSSLz51ulyoYN4dNPrQBZXkaNspxH775rRchKBVW47jpLrvTMMyYyP/bdF+bNg9q1/7GpenV7T127WsG0n3/eZmw6XK64wj78//s/Gx/p06fIQ5zoUOCYgoisF5F1BbXSFOk4xWXtWrvXpqXZzbNRo/z369EDJk6EU04pRXEPPWS5jP7zH0s4lx+TJlniuTp1Cqzd2aCBvbeMDHPjr16dRM3FZdAgG/y44w7L1+SUGQo0CqpaQ1VrAo8B/wEaA02w8NRBpSPPcYpPeroNKP/+u5Ur3nPP/PfLDpLp0CGBNZOL4sMPzRiceSbce2/++8ybZwO399xT5OnatLFTzp4Np50GW7cmVm6JSUmxkfBDD4VLLoGffgpbkROUovxLwNgg60qr+ZiCUxRXX20u7RdeKHifBQtUq1dXffHF0tOlv/5qAxwHHGB1PgvizDOtvGYxaiW/8oq958sui5gbf/ly1d13V23YUHXevLDVlGtIYEK8TBE5V0RSRSRFRM4FMpNopxynxDzzDDz9NNxyi81HKIjbbjO3yxFHlJKwJUvghBOgXj17tK9aNf/9fvgB3n7bBDZtGvj0F1xgYarPP29zMCJD/fqWl3zLFnv/69eHrcgpiqKsBtAcq6mwAlgOfAg0D2JxktG8p+AUxDffWJjmsceqZmQUvN/o0fZU3a9fKQnbvFl1//2tgtnEiQXvl5lpIUWNG6tu2FDsy2RmWrrvlJQIFkkbOdL+OSeeaEKdUodkluMMs7lRcPJj1iyrUtmuXf6hp9lkZKjuu69q06aqGzeWgrCsLNULLrCf2vvvF77vggWqLVrk1k4oAevXq+69t9V9njmzxKdJDo89Zp/DXXeFraRcEtQoBApJdZwos369eSYARozIP/Q0m4kTYcYMePVVS5WddAYPtov17w8nn1z4vk2bwrRpNrW6hFSvbp9Bly72mYwda+mUIsH111tU1cCBsPfecPrpYSty8sGznTplmqws86fPnAnvvAMtWxa+/377WaROqdyPvvjCcmucemrRNRA+/RQ2bbKxhpTt+1k2b27zLmbPhvPOi1BJZREb9OnaFS66yAyEEzmCZElNLQ0hjlMSBg60cdtHHrGs04UxaZLNG2vcuBRCUGfNsrDT9u3h5ZcLv9FPmgTHH19wiGoJOOwwePxxG+Pt1y9hp91+KleG99+3+RcnnQQrVoStyMlDkEeSWSLykIi0S7oaxykGH35oXpkLLzTPRGH88oslu3v++VIQtmGD3fBSUkxk9eoF76sK114LdetayFQCueoq+Ne/zNb8978JPfX20aiRlbpbssQMZ0ZG2IqceIoadMDSZF8G/Aj8DFwO1AwyYJGM5gPNjqrqtGk2z6BLFwvuKYz0dMt9tOuuhQ9CJ4SsLNVTTrEQoC+/LHr/7AkGhU2q2A7S0lS7drXpEZMmJeUSJSf7vf/732ErKReQjOgj4FBgEVZb4RWgVXGOT0Rzo+CsXq3aqpXqzjur/vVX0fs//rh90995J/nadNAgu9gjjxS97/LlqvXq2V07iWGaixer7rKLBTatXJm0y5SMG26wz+uVV8JWssOTMKMApGJFdj4AfgNuAnYGTgP+CHKRRDY3CuWbjAzVo49WrVhRdcyYovdftEi1Rg3VXr1KYabvxx+riqiee26wi82erXrYYapTpiRZmOpPP6lWqqTas2fEMlpv3ap6+OGqlStbSlsnaSTSKMwBXgQOymfbE0EuksjmRqF806ePfWuHDg22//ffqzZrpvrnn8nVpTNmWJrrTp0KT2ERIi++aJ/dzTeHrSQPy5ap7rabapMmqkuXhq1mhyWRRuHgfNZ1C3LyZDQ3CuWXt9/WnPw+xWHr1uToyWHNGtW2bVUbNFCdP7/o/TdtMj96CDfAa66xz/C110r90oXz66+W7+ngg0u5sEX5IahRCBJ99EQ+654sxli242w3kybBxRfDQQcFy+2zciU88YQFtlSsmERhmZk2GWD2bJsc0KxZ0cfccw889phNVCtlBg+2xKX/+hdMmFDqly+Yjh3hxRdhzBj497/DVlOuKXBGs4h0BQ4CGojITXGbamLjDI5TKqxYYRGedepYKuwgxXBuusnS+PfoAe2SGUzdr59NBhgyxO62RTF2rNVTuPTSUszGl0vFihae2qWLTbAeP96KEEWCs8+2KecPPmj5zK+4ImxF5ZLCegqVgOqY4agR19Zhg8xFIiK9RWSmiMwSkf/ks/0mEZkuIpNF5GsR2a34b8HZkcmujbBkiYX8F1QsJ57PPrPMEn36JNkgvP023HcfXHYZXH110ftv3mwzeRs3hkcfTaKwwmnY0D7LFSsiVoMB7PPs3dvmbnz/fdhqyidF+ZeA3YL4ofI5LhWYDeyOGZhJQLs8+xwOVIu9vgp4u6jz+phC+SK7NsKrrwbbf8UKC79s395i9JPGhAnmA+/WLbgP/M477c1EJIXpm2/mjtFEqgbD6tWqbdrYGI3XYEgYbG9CPBF5TFX/DQwREc3HmJxQhL3ZH5ilqnNi53sLOBGYHneOb+P2/xk4r4hzOuWIoUOtNsKtt8L55wc75tpr7Qn400+3K69c4Sxdav6s+vWD+7OyxTVrBj17JklY8TjrLJgyxR7O997bykZHgtq1Lavf/vtbAeoxYwqfFe4klMKypA6P/X24hOduDPwVt7wQOKCQ/S8FPstvg4hcjs2kplmQgTynzPPtt3aTOuYYq/8elOuusxxI++6bJGFpaeaMX7nS3Bs771z0MZs2meHYeWdzNUWIgQNtvPvGG2GPPSJjr6BtW3PPHXusPRG89952Jwp0AhKkO1GSBpwOvBC3fD7wZAH7nof1FCoXdV53H+34/PGH1UbYc0+L9gxC0sNOVc3Hcu655nN5773gx51/vmr37oVX/gmRdetU99rLajDMmBG2mjwMHmyf9x13hK2kzMP2hqSKyJTYAHC+LYC9WQjE1xNsAizO5zo9gL7ACaq6JZAlc3ZYVq2C446zLKYffxysFkBGBhx1FNxxR5LF3X8/vP46DBoEp5wS7JjXX4fhw6F7d0iNZtBejRr2WVesaJ99pBKX3nCD9a7uuw9eey1sNeWDgqwFsFthrShrg7mm5gAtyB1obp9nn47YYHTrIBZMvaewQ7Nli2U8qFjRSmYG5a67ijcYXSKyZ86dc07wUdnZsy3HRrduEcstkT8//mjZJg45JMmD9MVlyxbraVWqVLwvhrMNRKEcJ3AM8Efsxt83tm4A1isA+Ar4G5gYayOKOqcbhR2TrCzVSy8t/s39yy8t3dDFFydPm/7wg90tu3UrOiVrNmlpVpe5Zk3VuXOTKC6xZEckXXBBxCKSVq60WeN160awzmjZYLuNAjAm9nc9Njdhm79BTp6M5kZhx2TgQPs23nln8GOWLLFMqe3alajOfTD+/NMymbZqZVlNgzJ3rh1TnLGHiDBggP0v+vULW0keZs1SrV9ftWVLy5fkFItI9BSS0dwo7Hi89FLJnk4//9weHKdOTZKw5ctVW7c2o/DHH8U/PlI+mODE99qefz5sNXnI9nF17aq6cWPYasoUQY1CoBgvEekkIteLyHUi0rHYAxeOUwBffGHjiD16WFW04pTJ7NUL5s61ipcJZ+NGC4f86y+b/tu6dbDjpk61MnBpaUmcKJFcsksp9+4NV15pcz4iQ9euNnj/88820cKrtiWcIDWa+2EFdeoB9YGXReTOZAtzdnzGj7ea9u3aFW8O2PDhlsYCoGbNJAhLT7f8D+PHw1tvwcEHBztu5Uqbw/DOO7BmTRKElR7ZOZI6dIDTT7d7cGQ49VTLNfXxx2a19B9za53toaiuBDADqBK3XBWYEaQbkozm7qMdg+nTzSvTvLkVwgnKL7+Y9+CII5JUrCwz0+YVgOpzzwU/Lj5CJkj1nzLCkiXmwq9Tp1RqARWP7LCzvn3DVlImIIHuo3lAlbjlylg0keOUiPnzbV5BhQrw5Zew667BjzvhBEuK9/bbSZjgqmpx8cOHw4ABwWcfq1oKi1Gj4IUXoFu3BAsLj0aN7H9Utar9z+bMCVtRHPfcY/+je++Fh0uaeMHJS2G5j54EFNgCTBORL2PLPYExpSPP2dFYutRuLhs2wHffQatWwY5bs8ZSXmzeDF99ZWmHEoqqpVUdMsTybt9ZDA/pvHnmZurTJ3iSpjJEixY29nPooZYGY/RoS/QaOtmDH+vWWYKsatWCZat1CqWw3EfjY38nYPWZsxmVNDXODs2yZVZCYNEiu8nss0/wY0eMgD//hM8/T9LA8r33wgMPmI/64YeLN+LdooXVAWjePAnCokH79jbg3LOn/Q9HjYJddglbFTZLfPhwyy91zTWw005w4YVhqyrbBPExRan5mELZZNkyS2ddrZrqd9+V7ByzZiVWUw733We+6fPPL95AxRdfWG6ecsSYMao77aS6xx4RK6e8ebNqjx6qKSlJntpediFRYwoi0lpE3o0Vw5mT3UrBXjk7CMuXW8jpnDlWpCxIgTIwj86tt+ZGvrRsmQRxAwda0qRzzoFhw4IPVIwbZ5FGw4ZZ+Gk5oVs36zEsWGA9hr//DltRjCpVLHS4e3frKbz8ctiKyixBfgEvAc8AGVhRnFfJTavtOIWyaJEZgT//NBfQ4YcHO04VbrnFPDmff54EYapw991WTvP88y3GtUJh3tQ4Jk2Co4+GBg1g5Ei7IZUjDj0U/vc/G0o55BAzEJFgp50sTLVHD7jkEhv0d4pPUV0JYELs75S4dd8H6YYko7n7qOwwe7ZqixaWE644LqOsrNwiZdddl4QcPJmZqjfeaBe4+OLipbSeONFiaZs0SaI/q2zwww+qtWqpNmtm2UAiw+bNqr172//3kUfCVhMZSGBIapqIpAB/isi1InIyEJVS305EmTLFniLXroWvvy6ey+i22yw79aWXwmOPFW/Mt0jS0829MHiwzTx+4YXipbT+9Vd7Ih01Kkn+rLLDQQdZMaRNm+x/PXFi2IpiZLuSTjsNbr4Zbr/dJ7gVh6KsBtAFqI7VQ3gJeB84MIjFSUbznkL0+eorSw66666qkycX79iMDNVTTlG95pokTE7bsEH16KPtCfLee4vXBVm3btvzODlMn67atKlq9eqWjyoyZGSoXnllbo+wVCoxRRcSnRAPqAnUCLp/spobhWjzyiuqFSpYpNGCBcGP27LFZs+q2m834S6jv/5S3Xdfi04pzkxlVdVPPrHMez/+mGBROw4LF6rus49qaqrqiy+GrSaOrCzVu++2W12PHqqrVoWtKDSCGoUg0UedRWQKMBmYIiKTRGS/5PRbnLJKZqbN97rwQnMljBkDTZsWfRzA6tWWfK1nT9i61fLuJNRlNG4cdOkCs2db+FNx6iS//LIVj2/Roty7iwqjcWMrWX3EEeb2u/12+06Ejgj07w8vvWSzJbt2tagHp2CKshqYMTgkbvlgYHIQi5OM5j2F6LF6teqxx9rD2KWX2lN/UGbPtpj3ihVVhw9PgrjXXlOtUsWSLBUneU9Ghuptt9mb6tlzW/eRUyBbt6pecYV9bL16RezB/LvvLEigTh3VkSPDVlPqkCj3EfBDkHWl1dwoRIspU1TbtDGX0dNPF8/t8+WXVjOlTh3VUaMSLCwtLdeffMghqn//Xbzjhw2zY6+8snhWzlFV1WefNUPfsqUFbEWGWbPMtymi2r9/krIqRpPtNgpAp1gbDDwLdAcOA54G7g1y8mQ0NwrRICtLdehQewjfeefil87NzFTt0sV+nwmvrjhnjmrnzvb1vvXW4tVHzi6Mk5Gh+vHHCRZWvvjhB9VddrGstk8+GaHynhs25GbC7dWr3FRxS4RR+LaQ9k2QkyejuVEIn1WrVE891b49Rx1VvHQHS5dauV1VG5xMaCBPVpY94VevbuFPH3xQvGOHDLGg++wRb2e7WbZM9Zhj7Lty4omqK1aErShGVpZ1ZypXVm3YUHXEiLAVJZ2EuY+i1twohMuHH6o2amTuogcfLF7v+/33zV103nlJEPb336onnWRf6e7dVefNC37s0qWqJ5xgxx5zTLl5ciwtsrIsRVTFinb/fffdsBXFMXmyaocO9r+/7LIdeuwokWMKtYBHsayp44FHgFpBTp6M5kYhHJYtUz3rLPvGdOig+uuvwY9dsSK3t96pk+q0aQkUlpmp+sILFjJaqZLNYA1qqbKyLIa2Th07dvDgCPk4djwmTbL/P6iefnqEEuqlpanefruNMzRpYk8+OyCJNArvAfcAu8fa3cD7QU6ejOZGoXRJT1d94gnV2rXtSW/AgOLNAfr2W7tfp6ZaoayEzh+aOtUGkbMHk6dPL/45TjtNtVs31RkzEijMKYj0dEtKW6mSefgGD47QnLKfflLde2/N8XUVp7dZBkikUZgYZF1pNTcKpUNWlkUHZf9GevQo3j1340b7u2SJ6nHHJbiU45IlqpdfbhPR6tSx2VJBewerV6vedFOuoHXrylUESlT4/ffc9ETt2ql+9llEOmlbt6o+8IBq1ao23nD77apr1oStKiEk0ij8BBwct9wN+CnIyZPR3Cgknx9/VD38cPt2NG9uYwFBf7DTptlD1kEHJeFHvmqVar9+NpBcoYLq9derLl8e7NiNG821VL++uQmeeCLB4pzikpWl+tFHqrvvbt+1Qw9V/f77sFXFWLAg1+dZv759d7KfdMooiTQKHYBJWK3mecBvwD5BTp6M5kYhOWRlqX7zTe7TW8OGdt/MjtAsil9/VT3zTHt4r1lTddCgBLoFli9XveMOS7cKlhzpjz+CH//00/aGsiei/fZbgoQ5iSAtzQK/GjXK/Rd99VVEeg4TJqgecUTuj+KBB8rsYHRCjAJWb+GM2OuaQM0gJ01mc6OQWNLSVN94Izesv2FD8/muXx/8HB98YMfWrGm97YSFHU6apPqvf9lkCBEbnZw0Kdix8+fnuoX+8x+704wZkyBhTjLYuNEi2rKNQ8eONst98+awlal1YY46KveLfuONNh2/DJHInsLoICcqreZGITH8/rvqLbdYzxhUW7e2sO0gP8C//rIko9lpKTZuVH344QS5XtevV33ppdwB5KpVbfwgSMhSWprFOx57rHVZ/vc/W1+ceglO6KSlWUBZ27b2FahXz4aBShJHkHDGjlU9+2xzX4pY1t133omI5SqcRBqFu4BbgKZA3ewW5OTJaG4USs78+aoPPZQbFpiaap6YkSOLHmtdtsyMxmGH2W8BLMdNQti0ycIAzzvPCgBnW6kHHsid6VYYa9eamLp17dhdd1Xt08dmxzlllsxMC3Y47TS7B4Mlun3ggQgEBi1caJWgGjc2YXXq2Hfwiy8iFE61LUGNgti+BSMic/NZraq6e6EHJonOnTvr+PHjw7h0mSMjwxKEfvqpVSmcNMnWd+4MZ59tbZdd8j82Kwt+/x3atbPlo4+2spht28K559qxrVpth7i//rJSlp9/bm3jRqhTxwqjXHSRZbMsKFXqsmXwxReWUvWSS0zsHntAp052bM+exSuc40SepUvh7bfhjTfgl19s3T77wHHHwbHHWhLcihVDEJaZCd98Y9l0P/oo93t87LG5qX8bRqMmmYhMUNXORe5XlFGIGm4UCmbTJpgwwQrdjxplqYzXr7da9N26wfHHw0knQevW/zw2LQ1++w3GjoXRoy3L8KpVVpi9YUMzLhUrQocOJUhrrQpz58IPP1hO7dGjzeIA7LqrCTv1VCu6XtAv+/PPzbKNGQOTJ9u6jh2tEhqYYUgJUkjQKevMng0ffGBZ0MeMsfty9eqWsr17dzjwQNhvPyuQV6ps3gxffgnvvWdPYitW2Pp99zVxBx9s5eoaN05wbvhgJMwoiEgV4GosZbYC3wNDVTUtgIjewONAKvCCqt6fZ3tl4FVgP2AlcKaqzivsnG4U7B67eLHdV6dMsR7AxIn2OjuHfdu2cPjh1nr0gLp1c49dtAhmzoSpU61UQPPmMGyY5cEHWz78cMuNf9JJ9oMLzKZNJmzaNLvAr79aW7XKtteqZT+MI4+EXr2gffvcH0j8sdOmwYwZ9utPSbEaCG++accedpgd26mTG4JyzurVVu7122/tgT37WSM1Ffbay+7HHTrA3ntbZ7LU7sdZWfa9HznShP38s32/AXbe2axWx44msn17aNMGKldOqqREGoV3gPXAa7FVZwN1VPX0Io5LBf4AegILgXHA2ao6PW6fq7Hw1itF5CzgZFU9s7DzlgejkJ5uHpLFi+0GvnAhzJ8P8+bBnDnwxx+wYUPu/jvvbF3pLl3se9aihXlWFi2y79wee8D06eb2mTVr22Nffx3OOceuNXYsHHCAPbzni6r9Cv/+2w746y9YsMDEzZ5tbeHC3P0rVrRfY6dOJmzPPaF2bXtzCxfasdddB/XrwyOPwC23bHtsu3b25NWggV23Rg2oUCGRH7Wzg7F8uX2Px441N9PkyeZ6ymannez+u/vu9vDTvLkZisaN7XvfsCFUqpQEYenp9uT2009mLCZMsB9lVpZtT0mBZs2skFPLlrDbblalqmlT8/E2agQ1a26XRUukUZikqh2KWpfPcV2B/qraK7bcB0BV/y9un5GxfX4SkQrAUqCBFiIqCkZB1f6XGRnWtm61//nWrda2bLGWlmY9ys2b7SFhwwZz56xda23NGltev97ueUuWwMqVtpyXSpXsIbt2bbs3Vq1q983TT4err7Zj27SBDRuUCmRQkXQqks6gfulcd2U6f/+1lTtv3ULLJlto1SSN5rtsYfdGm6hbZZOJW78e1q3LFbd6da7ANWvsRr5ihb3hvDRqZF/cjAyoVs0Kp6ek2Ifx+OP2VPTaa3D++dseJ2K/3i5d7Enqq6/MgrVvbwMWoTiJnR2NZcusFz1zprU//rAHrHnz7Deal9q17TmkXj3rYdepY+tq1rRWvbq1nXayr3u1avZ7rFLFHvYrV7bfa6VK9hXObhUq2M8i576elmZisnvFs2fbU9ucObmup3iqVIEhQ3K79MUkqFEI8tj1m4gcqKo/x058APBDgOMaA3/FLS8EDihoH1XNEJG1QD0gn09k+/juwmHs8+rNVGMT8RYnnUosqWB1IxtmLKYKm7c5bguVWZLSBFVorAupzJZttqdTlaXYaG1TFlCRbW+am6jG3+yMoDRlAalkIXEKNlGNDZXqkZIK9TYvRMgi/lkgvUJVKtevQQpZZC5dDsvJOV5Qsr6vADen0igjg3UZGfzjOWKAtZ2B54v5mQF2c27Xzp5iPvnkn9uHD4fzzrMucs+eZrGyf0W1a+cakQMPhCefNAPSqJE9Ae26a+6N/8ADrTlOgmnY0LyVRx657XrV3E7vokXWli2z3sayZbljajNm5D4rJaLEaIUK5t5KTa1Cauo+pKTsQ0qKrUtJAUmFnRptYteshTTO+ouGWUtpmLWUBllL2W1he3pvv4TC9QXY5wDgAhFZEFtuBsyI1W1WVd2ngOPy6+fk7QEE2QcRuRy4HKBZs2YBJP+TarvVZ0WVJuyUvmabC29JrcLaXfcCgQpLoGr6to/pWyrsRFoz2756QRaVMzduu71yTVJb7YWkwNrfM6iYuQUkdn6B9Op1aNBhT1IrwOafMkjJSkcERMT+1m9AkwP2ABE2jPzBjEL2thQhpfGupHaw7fL1KCQFUlKE2ElIadkS2rdHsrLM1ZKSkttSU82vtN9+dnP+4INtH18qV7aBg0MOsaf6t9+2R57sfSpVMv99u3bWc/juu9xHoqpV7REp+//Rvbtdo6DubatWcO21JfrfOU4yEMl9RunUqej9Va3Xv2GDBRlt2GCd7GyPQFparocgr/cgM9N+Hunpua8zM83jkJlpLdsDkZUFqtXIympDVlYb1iusU/hToc2hpfC5BHAf7VbYdlWdX8BxO6z7yHEcp6yRMPdRQTf9AIwDWotIC2ARcBZwTp59RgAXYkn3TsMqupWtGFnHcZwdiKSFcsTGCK4FRmIhqcNUdZqIDMBm1o0AXgSGi8gsYBVmOBzHcZyQSGp8n6p+CnyaZ12/uNdpQKGhrY7jOE7p4TN/HMdxnBzcKDiO4zg5uFFwHMdxcnCj4DiO4+TgRsFxHMfJocylzhaR5UBJ507UJwkpNBJEVLVFVRdEV5vrKj5R1RZVXVB8bbupaoOidipzRmF7EJHxQWb0hUFUtUVVF0RXm+sqPlHVFlVdkDxt7j5yHMdxcnCj4DiO4+RQ3ozCc2ELKISoaouqLoiuNtdVfKKqLaq6IEnaytWYguM4jlM45a2n4DiO4xRCuTUKInKLiKiI1A9bSzYiMlBEJovIRBH5QkQKqpZcqojIQyLye0zbByJSO2xNACJyuohME5EsEYlEhIiI9BaRmSIyS0T+E7YeABEZJiLLRGRq2FriEZGmIvKtiAsMOCAAAAfxSURBVMyI/R9vCFtTNiJSRUR+EZFJMW33hK0pHhFJFZHfRCSfcojbR7k0CiLSFOgJLChq31LmIVXdR1X3BT4B+hV1QCnxJbBXrMreH0CfkPVkMxU4BRgdthCwHyrwFHA00A44W0TahasKgJch6VUcS0IGcLOq7gkcCFwTkc8LYAtwRKwW/b5AbxGJUr3YG4AZyThxuTQKwGDgNvIp/RkmqroubnEnIqJPVb9Q1ezi0z8DTcLUk42qzlDVmWHriGN/YJaqzlHVrcBbwIkha0JVR2P1SiKFqi5R1V9jr9djN7nG4aoy1NgQW6wYa5H4PYpIE+BY4IVknL/cGQUROQFYpKqTwtaSHyJyr4j8BZxLdHoK8VwCfBa2iIjSGPgrbnkhEbnJRR0RaQ50BMaGqySXmItmIrAM+FJVo6LtMeyhNisZJ09qkZ2wEJGvgEb5bOoL3AEcVbqKcilMm6p+pKp9gb6xmtbXAndHQVdsn75Yl//10tAUVFeEkHzWReLpMsqISHXgPeDfeXrLoaKqmcC+sTG0D0RkL1UNdVxGRI4DlqnqBBHpnoxr7JBGQVV75LdeRPYGWgCTRATMDfKriOyvqkvD1JYPbwD/o5SMQlG6RORC4DjgyNKso12MzysKLASaxi03ARaHpKVMICIVMYPwuqq+H7ae/FDVNSIyChuXCXuwvhtwgogcA1QBaorIa6p6XqIuUK7cR6o6RVUbqmpzVW2O/Yg7lZZBKAoRaR23eALwe1ha4hGR3sDtwAmquilsPRFmHNBaRFqISCWs5viIkDVFFrEnsxeBGar6aNh64hGRBtlRdiJSFehBBH6PqtpHVZvE7l9nAd8k0iBAOTMKZYD7RWSqiEzGXFxRCdEbAtQAvoyFyw4NWxCAiJwsIguBrsD/RGRkmHpig/HXAiOxQdN3VHVamJoARORN4CegrYgsFJFLw9YUoxtwPnBE7Hs1MfYEHAV2Ab6N/RbHYWMKCQ//jCI+o9lxHMfJwXsKjuM4Tg5uFBzHcZwc3Cg4juM4ObhRcBzHcXJwo+A4juPk4EbB2S5EpL+I3BJ7PUBECpxsJiInFZbwTESuFJELYq9HFSfzqYjUFpGr45Z3FZF3gx6fKGLx7WNjGSwPCXjMvLCz9Ubl83PCx42CkzBUtZ+qflXILidh2UP/gYhUUNWhqvpqCS9fG8i5qanqYlU9rYTn2h6OBH5X1Y6q+n0I18/O1lpcovL5OSHjRsEpNiLSN1Yz4Cugbdz6l0XktNjr+0VkeqwGw8MichA2S/uh2CSllrHewH0i8h1wQ3yvI8Z5IvJjbELf/rHzbrNPbFtz4H6gZezcD4lI8+z6AbHc+C+JyJTYE/zhsfUXicj7IvK5iPwpIg/G1qfG3svU2DE35vMZ7CYiX8fe39ci0kxE9gUeBI6J6aia55gjY9efIlbjoHLc5lvF8vf/IiKtYvufHtMwSURGx2l7SETGxa59RWx9d7HaBG8AU0TkgTxP/v1F5GYRqR7T+2tMR3YW11L9/JwIo6revAVuwH7AFKAaUBOYBdwS2/YycBpQF5hJ7uTI2vHb4841Cng6brl/3LlGAc/HXh8KTM27T2x5KtA81qbGrW8ed8zNwEux13tgdTSqABcBc4BaseX5WO6i/bAZrMTrz/M5fAxcGHt9CfBh7PVFwJB89q+CZVBtE1t+FUsABzAPS/AHcAHwSez1FKBxns/wcuDO2OvKwHgsn1d3YCPQIratI/Bd3PWnA82wfGc1Y+vqx/5/Utqfn7foNu8pOMXlEOADVd2kltEyv9w+64A04AUROQUoLF/S24VsexNy6gHUlJJXfDsYGB471+/YzatNbNvXqrpWVdOwG+du2I1udxF5UizvU36ZO7tiSQuJnfvgIjS0Beaq6h+x5VcwY5fNm3F/u8Ze/wC8LCKXAdkuoaOAC8RSOo8F6gHZObN+UdW5sff5G9AwNjbQAVitqgswA3CfWPqGr7DU3jsXoT0Zn58TUdwoOCWh0NwoajmA9seyX54EfF7I7huLcR3FUnfHf2+rFKYlRn4prbPZEvc6E6igqquBDlhv5RqCFTMpKl9MYRryHq8AqnolcCf29D1RROrFznOdqu4bay1U9YvYcXk/y3exntuZWMEfsDodDYD91Cr8/U3Rn2FpfH5ORHCj4BSX0cDJIlJVRGoAx+fdQSw/fi1V/RT4N1bOEGA9llgvKGfGzncwsFZV12Kulk6x9Z0w10lR5x6N3QwRkTaYG6XAim2xSKAUVX0PuCv7enn4EctSSezcY4p4L78DzbPHC7BEcN/FbT8z7u9PMR0tVXWsqvYDVmDGYSRwlVjKaUSkjYjsVMA134ppPA0zEGCunmWqmh4bG9gttr60Pz8nouyQ9RSc5KGqv4rI28BEzI2QX4RNDeAjEamCPWVmDzS+BTwvItdjN6qiWC0iP2JjF5fE1r1HrvtkHFYzGlVdKSI/xAZHP8NqJWfzNDBURKZgPY2LVHWLSIEPwI2Bl0Qk+6Epv5rU1wPDRORWYDlwcWFvRFXTRORi4L8iUiGmPT7bbGURGYs9qJ0dW/eQWDp1Ab4GJgGTMX//r2JvYDnWG8vvmtNihnuRqi6JrX4d+FhExmP/w99j+5b25+dEFM+S6jiO4+Tg7iPHcRwnBzcKjuM4Tg5uFBzHcZwc3Cg4juM4ObhRcBzHcXJwo+A4juPk4EbBcRzHycGNguM4jpPD/wP+1sLb5LRdVgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot(x, pdf1, 'b-')\n",
    "plot(x, cdf1, 'b--')\n",
    "plot(x, pdf2, 'r-')\n",
    "plot(x, cdf2, 'r--')\n",
    "xlabel('distributions of observations')\n",
    "ylabel('probability density / cumulative density');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "effect size (Hedges biased) = 1.0\n"
     ]
    }
   ],
   "source": [
    "print('effect size (Hedges biased) =', hedges_biased(mu1, sd1, ns1, mu2, sd2, ns2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "effect size (Hedges un-biased) = 0.9747899159663865\n"
     ]
    }
   ],
   "source": [
    "print('effect size (Hedges un-biased) =', hedges_unbiased(mu1, sd1, ns1, mu2, sd2, ns2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "statistical power = 0.8074205649484903\n"
     ]
    }
   ],
   "source": [
    "print('statistical power =', power(mu1, sd1, ns1, mu2, sd2, ns2))"
   ]
  },
  {
   "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
}