confidence_interval.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import statistics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Keywords"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Normal distribution / standard normal distribution\n",
    "- t-distribution\n",
    "- Confidence interval\n",
    "- population mean / population proportion\n",
    "- t-statistic / z-statistic"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [Normal Distribution](https://en.wikipedia.org/wiki/Normal_distribution)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Gaussian distribution with 2 parameters: mean and variance.\n",
    "- Symmetric distribution whose meam, mmedian, and mode are the same.\n",
    "- Formula:\n",
    "$$\n",
    "    f(x) = \\frac{1}{\\sigma\\sqrt{2\\pi}}e^{\\bigl({-\\frac{1}{2}}\\bigl(\\frac{x-\\mu}{\\sigma} \\bigr)^2  \\bigr)},\n",
    "$$\n",
    "where &mu; is the mean of the distribution, &sigma; is the standard deviation.  Here, &mu; can be median or mode as well.\n",
    "- The variance of the distribution is &sigma;<sup>2</sup>.\n",
    "- The normal distribution is important because of the [central limit theorem](https://en.wikipedia.org/wiki/Central_limit_theorem)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "See [numpy.random.normal](https://docs.scipy.org/doc/numpy-1.14.0/reference/generated/numpy.random.normal.html) and [scipy.stats.norm](https://docs.scipy.org/doc/scipy-0.16.1/reference/generated/scipy.stats.norm.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAEiZJREFUeJzt3W2MXdV97/HvLw5Jqqa6QHEt17ZqlLqqSKU41RSo6AtuooADV4VILQK1iRshOZVASqSoDeQN6QOSI92Em+qmSE5x49ybxrGaVFipVeoSpCgvAgyJS7BpxJQYYcvB00KeFBXJ5N8XZ7k9mBnPmfF4zhmv70c6mr3XXvuc/96y5zdrP5ydqkKS1J/XjbsASdJ4GACS1CkDQJI6ZQBIUqcMAEnqlAEgSZ0yACSpUwaAJHXKAJCkTr1+3AWczWWXXVabN28edxmStKo88cQT/1ZVaxfqN9EBsHnzZqanp8ddhiStKkmeG6Wfh4AkqVMGgCR1ygCQpE4ZAJLUKQNAkjplAEhSpwwASeqUASBJnTIAJKlTE30nsDTJNt/194vqf3TnjeepEmlpHAFIUqcMAEnqlIeApAUs9lCPtFo4ApCkThkAktQpA0CSOmUASFKnDABJ6pQBIEmd8jJQqfFyT/XGAJBWyHwB41dEaFw8BCRJnXIEoO54qEcacAQgSZ1aMACSvCnJY0n+OcnhJH/S2i9P8miSmSRfTPKG1v7GNj/Tlm8eeq+7W/t3klx/vjZKkrSwUUYALwPvqKq3AVuBbUmuBj4O3FdVvwy8BNze+t8OvNTa72v9SHIFcCvwVmAb8JdJ1iznxkiSRrdgANTAj9vsRe1VwDuAv23te4Cb2/RNbZ62/J1J0tr3VtXLVfVdYAa4clm2QpK0aCOdA0iyJskh4CRwEPhX4PtVdap1OQZsaNMbgOcB2vIfAD8/3D7HOpKkFTZSAFTVK1W1FdjI4K/2Xz1fBSXZkWQ6yfTs7Oz5+hhJ6t6iLgOtqu8neQT4TeDiJK9vf+VvBI63bseBTcCxJK8H/gfw70Ptpw2vM/wZu4BdAFNTU7W4zZFWH28Q07iMchXQ2iQXt+mfAd4FPA08AvxO67YdeLBN72/ztOVfrapq7be2q4QuB7YAjy3XhkiSFmeUEcB6YE+7Yud1wL6q+kqSI8DeJH8OfAt4oPV/APh/SWaAFxlc+UNVHU6yDzgCnALuqKpXlndzJEmjWjAAqupJ4O1ztD/LHFfxVNV/AL87z3vdC9y7+DIlScvNO4ElqVMGgCR1ygCQpE4ZAJLUKQNAkjplAEhSpwwASeqUTwTTBcsnf0ln5whAkjplAEhSpwwASeqUASBJnTIAJKlTBoAkdcoAkKROGQCS1CkDQJI6ZQBIUqcMAEnqlAEgSZ0yACSpUwaAJHVqwQBIsinJI0mOJDmc5IOt/WNJjic51F43DK1zd5KZJN9Jcv1Q+7bWNpPkrvOzSZKkUYzyPIBTwIer6ptJfg54IsnBtuy+qvrfw52TXAHcCrwV+EXgn5L8Slv8aeBdwDHg8ST7q+rIcmyIJGlxFgyAqjoBnGjTP0ryNLDhLKvcBOytqpeB7yaZAa5sy2aq6lmAJHtbXwNAksZgUU8ES7IZeDvwKHANcGeS9wHTDEYJLzEIh28MrXaM/w6M589ov2pJVUsdmO+JZkd33rjClehCNfJJ4CRvBr4EfKiqfgjcD7wF2MpghPCJ5SgoyY4k00mmZ2dnl+MtJUlzGCkAklzE4Jf/56vqywBV9UJVvVJVPwU+w38f5jkObBpafWNrm6/9VapqV1VNVdXU2rVrF7s9kqQRjXIVUIAHgKer6pND7euHur0HeKpN7wduTfLGJJcDW4DHgMeBLUkuT/IGBieK9y/PZkiSFmuUcwDXAO8Fvp3kUGv7KHBbkq1AAUeBDwBU1eEk+xic3D0F3FFVrwAkuRN4CFgD7K6qw8u4LZKkRRjlKqCvA5lj0YGzrHMvcO8c7QfOtp60WPOdKJW0MO8ElqROGQCS1CkDQJI6ZQBIUqcMAEnqlAEgSZ0yACSpUwaAJHXKAJCkThkAktQpA0CSOmUASFKnFvVEMEnj55PCtFwcAUhSpwwASeqUASBJnTIAJKlTBoAkdcqrgLQq+OhHafk5ApCkThkAktQpA0CSOmUASFKnFgyAJJuSPJLkSJLDST7Y2i9NcjDJM+3nJa09Sf4iyUySJ5P8+tB7bW/9n0my/fxtliRpIaOMAE4BH66qK4CrgTuSXAHcBTxcVVuAh9s8wLuBLe21A7gfBoEB3ANcBVwJ3HM6NCRJK2/BAKiqE1X1zTb9I+BpYANwE7CnddsD3NymbwI+VwPfAC5Osh64HjhYVS9W1UvAQWDbsm6NJGlkizoHkGQz8HbgUWBdVZ1oi74HrGvTG4Dnh1Y71trmaz/zM3YkmU4yPTs7u5jyJEmLMHIAJHkz8CXgQ1X1w+FlVVVALUdBVbWrqqaqamrt2rXL8ZaSpDmMFABJLmLwy//zVfXl1vxCO7RD+3mytR8HNg2tvrG1zdcuSRqDUa4CCvAA8HRVfXJo0X7g9JU824EHh9rf164Guhr4QTtU9BBwXZJL2snf61qbJGkMRvkuoGuA9wLfTnKotX0U2AnsS3I78BxwS1t2ALgBmAF+ArwfoKpeTPJnwOOt359W1YvLshWSfFKYFm3BAKiqrwOZZ/E75+hfwB3zvNduYPdiCpQknR/eCSxJnTIAJKlTBoAkdcoAkKROGQCS1CkDQJI6ZQBIUqcMAEnq1Ch3AksrZr67WSUtP0cAktQpA0CSOmUASFKnDABJ6pQBIEmdMgAkqVMGgCR1ygCQpE4ZAJLUKQNAkjplAEhSp/wuIOkCd7bvVzq688YVrESTxhGAJHVqwQBIsjvJySRPDbV9LMnxJIfa64ahZXcnmUnynSTXD7Vva20zSe5a/k2RJC3GKCOAzwLb5mi/r6q2ttcBgCRXALcCb23r/GWSNUnWAJ8G3g1cAdzW+kqSxmTBcwBV9bUkm0d8v5uAvVX1MvDdJDPAlW3ZTFU9C5Bkb+t7ZNEVS5KWxbmcA7gzyZPtENElrW0D8PxQn2Otbb52SdKYLDUA7gfeAmwFTgCfWK6CkuxIMp1kenZ2drneVpJ0hiVdBlpVL5yeTvIZ4Ctt9jiwaajrxtbGWdrPfO9dwC6AqampWkp9mnw++lEavyWNAJKsH5p9D3D6CqH9wK1J3pjkcmAL8BjwOLAlyeVJ3sDgRPH+pZctSTpXC44AknwBuBa4LMkx4B7g2iRbgQKOAh8AqKrDSfYxOLl7Crijql5p73Mn8BCwBthdVYeXfWskSSMb5Sqg2+ZofuAs/e8F7p2j/QBwYFHVSZLOG+8ElqROGQCS1CkDQJI6ZQBIUqcMAEnqlAEgSZ0yACSpUwaAJHXKAJCkThkAktQpA0CSOmUASFKnlvQ8AEkXhvmey3B0540rXInGwRGAJHXKAJCkThkAktQpzwHovPLZv9LkcgQgSZ0yACSpUwaAJHXKAJCkThkAktSpBQMgye4kJ5M8NdR2aZKDSZ5pPy9p7UnyF0lmkjyZ5NeH1tne+j+TZPv52RxJ0qhGuQz0s8D/BT431HYX8HBV7UxyV5v/CPBuYEt7XQXcD1yV5FLgHmAKKOCJJPur6qXl2hBJy8eviOjDgiOAqvoa8OIZzTcBe9r0HuDmofbP1cA3gIuTrAeuBw5W1Yvtl/5BYNtybIAkaWmWeg5gXVWdaNPfA9a16Q3A80P9jrW2+dolSWNyzieBq6oYHNZZFkl2JJlOMj07O7tcbytJOsNSA+CFdmiH9vNkaz8ObBrqt7G1zdf+GlW1q6qmqmpq7dq1SyxPkrSQpQbAfuD0lTzbgQeH2t/Xrga6GvhBO1T0EHBdkkvaFUPXtTZJ0pgseBVQki8A1wKXJTnG4GqencC+JLcDzwG3tO4HgBuAGeAnwPsBqurFJH8GPN76/WlVnXliWZK0ghYMgKq6bZ5F75yjbwF3zPM+u4Hdi6pOknTeeCewJHXKAJCkThkAktQpnwimZeGTv6TVxxGAJHXKAJCkThkAktQpA0CSOmUASFKnDABJ6pQBIEmd8j4ASSPzUZEXFkcAktQpA0CSOmUASFKnDABJ6pQBIEmdMgAkqVNeBqpF8WufpQuHIwBJ6pQBIEmdMgAkqVMGgCR16pxOAic5CvwIeAU4VVVTSS4FvghsBo4Ct1TVS0kCfAq4AfgJ8AdV9c1z+XxJk8HvCFqdlmME8D+ramtVTbX5u4CHq2oL8HCbB3g3sKW9dgD3L8NnS5KW6HwcAroJ2NOm9wA3D7V/rga+AVycZP15+HxJ0gjONQAK+MckTyTZ0drWVdWJNv09YF2b3gA8P7Tusdb2Kkl2JJlOMj07O3uO5UmS5nOuN4L9VlUdT/ILwMEk/zK8sKoqSS3mDatqF7ALYGpqalHrSpJGd04BUFXH28+TSf4OuBJ4Icn6qjrRDvGcbN2PA5uGVt/Y2jSBvONXuvAt+RBQkp9N8nOnp4HrgKeA/cD21m078GCb3g+8LwNXAz8YOlQkSVph5zICWAf83eDqTl4P/E1V/UOSx4F9SW4HngNuaf0PMLgEdIbBZaDvP4fPliSdoyUHQFU9C7xtjvZ/B945R3sBdyz18yRJy8tvA5V03niD2GTzqyAkqVMGgCR1ykNAnfNyT6lfjgAkqVMGgCR1ygCQpE4ZAJLUKU8CS1px3h8wGRwBSFKnDABJ6pQBIEmd8hxAJ7zhS9KZHAFIUqccAUiaGF4dtLIMgAuMh3okjcpDQJLUKQNAkjrlISBJE89zA+eHIwBJ6pQjgFXIE72SloMBIGnV8tDQuVnxAEiyDfgUsAb4q6raudI1rBb+pS/pfFrRAEiyBvg08C7gGPB4kv1VdWQl65B0YXNkMJqVHgFcCcxU1bMASfYCNwFdBIB/0Uvjdbb/gz2Gw0oHwAbg+aH5Y8BVK1zDovmLW7rwLdf/89UUJBN3EjjJDmBHm/1xku+Msx7gMuDfxlzDYq3GmmF11m3NK2dV1J2Pv2p2XDX/0iidVjoAjgObhuY3trb/UlW7gF0rWdTZJJmuqqlx17EYq7FmWJ11W/PKWY11T3rNK30j2OPAliSXJ3kDcCuwf4VrkCSxwiOAqjqV5E7gIQaXge6uqsMrWYMkaWDFzwFU1QHgwEp/7jmYmMNRi7Aaa4bVWbc1r5zVWPdE15yqGncNkqQx8MvgJKlTBsCQJLuTnEzy1FDbx5IcT3KovW4YZ41nSrIpySNJjiQ5nOSDrf3SJAeTPNN+XjLuWk87S82Tvq/flOSxJP/c6v6T1n55kkeTzCT5YrvAYSKcpebPJvnu0L7eOu5az5RkTZJvJflKm5/Y/XzaHDVP9H42AF7ts8C2Odrvq6qt7TVp5y9OAR+uqiuAq4E7klwB3AU8XFVbgIfb/KSYr2aY7H39MvCOqnobsBXYluRq4OMM6v5l4CXg9jHWeKb5agb4o6F9fWh8Jc7rg8DTQ/OTvJ9PO7NmmOD9bAAMqaqvAS+Ou47FqKoTVfXNNv0jBv/4NjD4io09rdse4ObxVPhaZ6l5otXAj9vsRe1VwDuAv23tk7av56t5oiXZCNwI/FWbDxO8n+G1Na8GBsBo7kzyZDtENDGHUs6UZDPwduBRYF1VnWiLvgesG1NZZ3VGzTDh+7oN8Q8BJ4GDwL8C36+qU63LMSYszM6suapO7+t7276+L8kbx1jiXP4P8MfAT9v8zzPh+5nX1nzaxO5nA2Bh9wNvYTB8PgF8YrzlzC3Jm4EvAR+qqh8OL6vBpV4T91ffHDVP/L6uqleqaiuDu9ivBH51zCUt6Myak/wacDeD2n8DuBT4yBhLfJUk/ws4WVVPjLuWUZ2l5ondz2AALKiqXmj/gX4KfIbBf/qJkuQiBr9IP19VX27NLyRZ35avZ/DX38SYq+bVsK9Pq6rvA48AvwlcnOT0PTWv+XqTSTFU87Z2GK6q6mXgr5msfX0N8NtJjgJ7GRz6+RSTvZ9fU3OS/z/h+9kAWMjpX6LNe4Cn5us7Du3Y6APA01X1yaFF+4HtbXo78OBK1zaf+WpeBft6bZKL2/TPMHiuxdMMfqn+Tus2aft6rpr/ZeiPgzA4lj4x+7qq7q6qjVW1mcHXxXy1qn6PCd7P89T8+5O8n2ECvw10nJJ8AbgWuCzJMeAe4Np26VYBR4EPjK3AuV0DvBf4djvOC/BRYCewL8ntwHPALWOqby7z1XzbhO/r9cCeDB5s9DpgX1V9JckRYG+SPwe+xSDcJsV8NX81yVogwCHgD8dZ5Ig+wuTu5/l8fpL3s3cCS1KnPAQkSZ0yACSpUwaAJHXKAJCkThkAktQpA0CSOmUASFKnDABJ6tR/AtGZpkoE3XovAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "obtained mean = 30.02 / exact mean = 30.00\n",
      "obtained median = 30.02 / exact median = 30.00\n",
      "obtained variance = 16.02 / exact mean = 16.00\n"
     ]
    }
   ],
   "source": [
    "mu, sigma = 30, 4 # mean and standard deviation\n",
    "num_samples=50000\n",
    "\n",
    "np.random.seed(1)\n",
    "s = np.random.normal(mu, sigma, num_samples) # generate random numbers that follow normal distribution\n",
    "plt.hist(s,bins=50)\n",
    "plt.show()\n",
    "\n",
    "print(\"obtained mean = %0.2f / exact mean = %0.2f\"%(np.mean(s),mu))\n",
    "print(\"obtained median = %0.2f / exact median = %0.2f\"%(np.median(s),mu))\n",
    "print(\"obtained variance = %0.2f / exact mean = %0.2f\"%(np.var(s),sigma**2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [Standard normal distribution](https://en.wikipedia.org/wiki/Normal_distribution#Standard_normal_distribution)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- The simplest case of the normal distribution (mean = 0, standard deviation = 1)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl0VGW67/HvQ0hCQECEKMgYMbaizDHCkSnBhCFKtAUcjh6ktZVuOPa1FaUvfdq7tD2O3dLatoDzjEza0JEhBgRRQAJEMSAYwNZE5kkgQEjy3D9SYRURSCVU5a3h+ayVRe1d7079SuGXnXfv2ltUFWOMMZGhnusAxhhj6o6VvjHGRBArfWOMiSBW+sYYE0Gs9I0xJoJY6RtjTASx0jfGmAhipW+MMRHESt8YYyJIfdcBqmrRooV26NDBdQxjjAkpq1ev3q2q8dWNC7rS79ChA7m5ua5jGGNMSBGRf/syzqZ3jDEmgljpG2NMBLHSN8aYCGKlb4wxEcSn0heRwSKyUUQKRGTCGcbdKCIqIkle6/7g2W6jiAzyR2hjjDG1U+3ZOyISBbwApAGFwCoRmaOq66uMawz8Dljpta4TcDNwOXAh8LGIXKKqZf57C8YYY3zly55+MlCgqltUtQSYBmSeYtyjwJPAUa91mcA0VT2mqluBAs/3M8YY44Avpd8a+MFrudCz7gQR6QG0VdWsmm5rTCj44YcfePHFF5k+fTqHDh1yHceYWjvrD2eJSD3gr8AdZ/E97gbuBmjXrt3ZRjLGL1SVqVOn8sorr7Bq1aoT6xs0aMCgQYMYP348V199tcOExtScL3v6RUBbr+U2nnWVGgNXAJ+IyHdAL2CO52BuddsCoKpTVTVJVZPi46v9FLExAaeq3HfffYwZM4bS0lIef/xxvvnmG5YsWcLdd9/NF198wcCBA8nKqvrLrTFBTlXP+EXFbwNbgAQgBvgSuPwM4z8BkjyPL/eMj/VsvwWIOtPr9ezZU41xqaysTO+55x4F9N5779Xy8vKfjdm9e7f27NlTo6OjddasWQ5SGnMyIFer6XNVrX5PX1VLgXHAAmADMF1V80XkEREZVs22+cB0YD0wHxirduaOCWKqyl133cWUKVOYMGECkyZNQkR+Nq558+Z8/PHHJCUlMXLkSGbMmOEgrTE1JxU/IIJHUlKS2gXXjCvvvfcet956KxMnTuTRRx89ZeF7O3ToEOnp6eTn57N+/Xpat7bzFIwbIrJaVZOqHWelb0yFPXv2cNlll9GhQweWL19OVFSUT9sVFBTQuXNnhgwZwuzZswOc0phT87X07TIMxng88MAD7N27l5deesnnwge4+OKLefjhh/nggw/44IMPApjQmLNnpW8MkJOTw+uvv8748ePp2rVrjbe///776dq1K2PHjuXAgQMBSGiMf1jpm4hXUlLCmDFjuPjii/nTn/5Uq+8RHR3NSy+9xI4dO/jjH//o54TG+I+Vvol47733HgUFBTz77LPExcXV+vtceeWV3HXXXUydOpVt27b5MaEx/mOlbyJaeXk5Tz75JF26dCEjI+Osv9+DDz5IaWkpkyZN8kM6Y/zPSt9EtLlz57JhwwYmTJhQ7emZvujYsSMjR47kxRdfZP/+/X5IaIx/WembiKWqPP744yQkJDBixAi/fd+HHnqIgwcP8uKLL/rtexrjL1b6JmItXbqUlStXMn78eOrXP+trD57QrVs3Bg8ezKRJkzhy5Ijfvq8x/mClbyLWE088wfnnn88dd9zh9+/90EMPsXPnTl5//XW/f29jzoaVvolIGzduZP78+dx7771ndcbO6fTv35/k5GQmTZpEsH3q3UQ2K30TkV577TWioqK48847A/L9RYR77rmHTZs2sXz58oC8hjG1YaVvIk5paSlvvvkmQ4YMoWXLlgF7nREjRtCwYUNee+21gL2GMTVlpW8izsKFC9m2bRujR48O6Os0btyYESNG8P7773P48OGAvpYxvrLSNxHntddeo0WLFlx77bUBf61f/epXHDx40K6+aYKGlb6JKHv27GHOnDncdtttxMTEBPz1+vbtS8eOHXn11VcD/lrG+MKn0heRwSKyUUQKRGTCKZ4fIyLrRCRPRJaJSCfP+g4icsSzPk9EJvv7DRhTE++++y4lJSUBn9qpJCLccccdfPLJJ2zZsqVOXtOYM6m29EUkCngBGAJ0Am6pLHUv76pqZ1XtBjwF/NXruc2q2s3zNcZfwY2pjddee40ePXrQpUuXOnvNUaNGISK88cYbdfaaxpyOL3v6yUCBqm5R1RJgGpDpPUBVf/JabATYickm6Kxfv561a9cG5MNYZ9K2bVuuueYa3nrrLTtn3zjnS+m3Bn7wWi70rDuJiIwVkc1U7Onf6/VUgoisFZElItL3VC8gIneLSK6I5O7atasG8Y3x3cyZMxERv15nx1cjR45k69at5OXl1flrG+PNbwdyVfUFVe0IPARU3kViG9BOVbsDvwfeFZEmp9h2qqomqWpSfHy8vyIZc5KZM2fSp0+fgJ6bfzrXX389UVFRzJw5s85f2xhvvpR+EdDWa7mNZ93pTAOuB1DVY6q6x/N4NbAZuKR2UY2pvU2bNrFu3TpuvPFGJ6/fokUL+vfvz8yZM22KxzjlS+mvAhJFJEFEYoCbgTneA0Qk0WsxA/jWsz7ecyAYEbkISATsFAZT52bNmgXAL3/5S2cZhg8fzqZNm8jPz3eWwZhqS19VS4FxwAJgAzBdVfNF5BERGeYZNk5E8kUkj4ppnFGe9f2ArzzrZwJjVHWv39+FMdWYNWsWV111FW3btq1+cIDccMMNiMiJH0DGuCDB9qtmUlKS5ubmuo5hwsjWrVu56KKLePrpp3nggQecZunXrx/79+/nq6++cprDhB8RWa2qSdWNs0/kmrBXuWftaj7f24033si6devYtGmT6ygmQlnpm7A3a9YsevToQUJCgusoJ44p2BSPccVK34S1oqIiVqxYERR7+VDxQa2rrrrKTt00zljpm7CWlZUFQGZmZjUj605mZiZr1qzhxx9/dB3FRCArfRPWsrKyaN++PZ06Vb1clDsZGRkAfPTRR46TmEhkpW/C1tGjR/n444/JyMhARFzHOaFz5860adPmxG8hxtQlK30TtpYsWUJxcfGJPetgISJkZGSQnZ3NsWPHXMcxEcZK34StrKws4uLiSElJcR3lZzIyMjh8+DBLly51HcVEGCt9E5ZUlaysLFJTU4mLi3Md52dSU1OJjY21KR5T56z0TVjauHEjW7ZsCbqpnUqNGjUiJSXFSt/UOSt9E5YqyzRYSx8qshUUFNinc02dstI3YSkrK4srrriCdu3auY5yWpU/kGxv39QlK30Tdn766Sc+/fTToN7LB0hISOCyyy6z0jd1ykrfhJ2cnBxKS0sZOnSo6yjVysjIYOnSpRw+fNh1FBMhrPRN2Fm4cCHnnHMOvXv3dh2lWunp6Rw/fpwlS5a4jmIihJW+CTvZ2dmkpKQQHR3tOkq1+vTpQ2xsLNnZ2a6jmAjhU+mLyGAR2SgiBSIy4RTPjxGRdSKSJyLLRKST13N/8Gy3UUQG+TO8MVVt3ryZzZs3k56e7jqKT+Li4ujXrx8LFy50HcVEiGpL33OP2xeAIUAn4BbvUvd4V1U7q2o34Cngr55tO1FxT93LgcHAPyrvmWtMIFTuMYdK6UNF1vXr11NYWOg6iokAvuzpJwMFqrpFVUuAacBJ16lV1Z+8FhsBlfdgzASmqeoxVd0KFHi+nzEBsXDhQtq1a0diYqLrKD5LS0sDsCkeUyd8Kf3WwA9ey4WedScRkbEispmKPf17a7jt3SKSKyK5u3bt8jW7MScpLS1l0aJFpKenB9VVNavTuXNnLrjgAit9Uyf8diBXVV9Q1Y7AQ8Afa7jtVFVNUtWk+Ph4f0UyEWbVqlUcOHAgpKZ2AOrVq0daWhrZ2dmUl5e7jmPCnC+lXwS09Vpu41l3OtOA62u5rTG1tnDhQkSE1NRU11FqLC0tjd27d5OXl+c6iglzvpT+KiBRRBJEJIaKA7NzvAeIiPcEagbwrefxHOBmEYkVkQQgEfji7GMb83PZ2dkkJSXRvHlz11FqzOb1TV2ptvRVtRQYBywANgDTVTVfRB4RkWGeYeNEJF9E8oDfA6M82+YD04H1wHxgrKqWBeB9mAh34MABVqxYEXJTO5VatWpF586d7dRNE3D1fRmkqh8BH1VZ9yevx787w7aPAY/VNqAxvvjkk08oKys7scccitLS0vj73/9OcXExDRs2dB3HhCn7RK4JCzk5OcTFxdGrVy/XUWpt4MCBlJSU8Pnnn7uOYsKYlb4JC4sWLaJv377Exsa6jlJrffv2pX79+uTk5LiOYsKYlb4Jedu3byc/Pz8kz9rx1rhxY5KTk630TUBZ6ZuQt2jRIqBieiTUDRw4kNWrV7N//37XUUyYstI3IW/RokWce+65dO/e3XWUszZw4EDKy8vtUssmYKz0TcjLyclhwIABREWF/rX8evXqRVxc3InfXozxNyt9E9K2bt3Kd999FxZTOwCxsbH06dPH5vVNwFjpm5BWWY7hUvpQ8V7y8/PZsWOH6ygmDFnpm5CWk5NDq1atuPTSS11H8ZvKs5BsiscEgpW+CVmqyqJFi0hNTQ2pSylXp0ePHjRt2tRK3wSElb4JWfn5+ezcuTPkz8+vKioqigEDBti8vgkIK30Tsir3hMOt9KFiXr/yILUx/mSlb0LW4sWLSUhIoEOHDq6j+F1KSgpQ8R6N8ScrfROSKj/AVFmO4ebyyy8nPj7eSt/4nZW+CUlffvkl+/btC9vSFxEGDBjA4sWLUVXXcUwY8an0RWSwiGwUkQIRmXCK538vIutF5CsRyRGR9l7PlYlInudrTtVtjamNyvn8cC19qDhWUVhYSEFBgesoJoxUW/oiEgW8AAwBOgG3iEinKsPWAkmq2gWYCTzl9dwRVe3m+RqGMX6wePFiLrnkElq3bu06SsDYvL4JBF/29JOBAlXdoqolVNz4PNN7gKouVtViz+IKKm6AbkxAlJaWsnTp0rDeywe45JJLaNWqlZW+8StfSr818IPXcqFn3encCczzWm4gIrkiskJErq9FRmNOsmbNGg4ePBj2pS8ipKSk2Ly+8Su/HsgVkduAJOBpr9XtVTUJuBWYJCIdT7Hd3Z4fDLm7du3yZyQThir3fAcMGOA2SB1ISUlhx44dfPPNN66jmDDhS+kXAW29ltt41p1ERK4BJgLDVPVY5XpVLfL8uQX4BPjZRc9VdaqqJqlqUnx8fI3egIk8ixYt4vLLL+eCCy5wHSXg7Do8xt98Kf1VQKKIJIhIDHAzcNJZOCLSHZhCReHv9FrfTERiPY9bAFcD6/0V3kSekpISli1bFvZTO5USEhJo166dzesbv6lf3QBVLRWRccACIAp4VVXzReQRIFdV51AxnXMOMMNz4avvPWfqXAZMEZFyKn7APKGqVvqm1latWkVxcXHElH7lvP6//vUvysvLqVfPPlpjzk61pQ+gqh8BH1VZ9yevx9ecZrvPgc5nE9AYb5V7vP3793ecpO6kpKTwxhtv8PXXX9OlSxfXcUyIs90GE1IWL15M165dad68uesodcbO1zf+ZKVvQsbRo0f57LPPwvKqmmfSrl07OnbsaAdzjV9Y6ZuQsWLFCo4dOxYx8/neUlJSWLJkCWVlZa6jmBBnpW9CxuLFi6lXrx79+vVzHaXOpaSkcODAAfLy8lxHMSHOSt+EjMWLF5+4lWCksXl94y9W+iYkFBcXs2LFioic2gFO3PzdSt+cLSt9ExI+++wzjh8/HrGlDxV7+0uXLuX48eOuo5gQZqVvQsLixYupX78+ffr0cR3FmZSUFA4dOsTq1atdRzEhzErfhITFixdz5ZVX0rhxY9dRnKm8wJxN8ZizYaVvgt7BgwdZtWpVRE/tAMTHx9O5c2crfXNWrPRN0Fu2bBllZWURX/pQMcWzbNkyjh07Vv1gY07BSt8EvUWLFhETE8N//Md/uI7iXGpqKkeOHOGLL75wHcWEKCt9E/RycnLo3bs3DRs2dB3Fuf79+1OvXj1ycnJcRzEhykrfBLU9e/aQl5fHwIEDXUcJCueeey49e/a06/CYWrPSN0Htk08+QVUj7iJrZ5KamsqKFSs4fPiw6ygmBFnpm6C2aNEiGjVqRHJysusoQSM1NZXjx4+zbNky11FMCPKp9EVksIhsFJECEZlwiud/LyLrReQrEckRkfZez40SkW89X6P8Gd6Ev0WLFtGvXz+io6NdRwkaffr0ITo62qZ4TK1UW/oiEgW8AAwBOgG3iEinKsPWAkmq2gWYCTzl2fY84GHgKiAZeFhEmvkvvglnRUVFfPPNNzafX0XDhg3p3bu3Hcw1teLLnn4yUKCqW1S1BJgGZHoPUNXFqlrsWVwBtPE8HgRkq+peVd0HZAOD/RPdhLvKDyHZfP7PDRw4kDVr1rBv3z7XUUyI8aX0WwM/eC0Xetadzp3AvFpua8wJOTk5nHfeeXTt2tV1lKCTmpqKqvLJJ5+4jmJCjF8P5IrIbUAS8HQNt7tbRHJFJHfXrl3+jGRClKqyaNEiUlJSqFfPzjeoKjk5mUaNGtm8vqkxX/41FQFtvZbbeNadRESuASYCw1T1WE22VdWpqpqkqknx8fG+ZjdhbPPmzXz//fc2n38aMTEx9O3b1+b1TY35UvqrgEQRSRCRGOBmYI73ABHpDkyhovB3ej21AEgXkWaeA7jpnnXGnFHlHqzN559eamoqGzZs4Mcff3QdxYSQaktfVUuBcVSU9QZguqrmi8gjIjLMM+xp4Bxghojkicgcz7Z7gUep+MGxCnjEs86YM8rOzqZNmzZccsklrqMErWuuuQbA9vZNjYiqus5wkqSkJM3NzXUdwzhUVlZGfHw8119/Pa+++qrrOEGrvLycli1bMnjwYN58803XcYxjIrJaVZOqG2dHyEzQqTwVMS0tzXWUoFavXj0GDhzIxx9/TLDtvJngZaVvgk52djaAHcT1QVpaGtu2bSM/P991FBMirPRN0MnOzqZbt26cf/75rqMEvcrfhip/UBpTHSt9E1QOHz7MZ599ZlM7Pmrbti2/+MUvrPSNz6z0TVBZunQpx48ft9KvgbS0NJYsWWK3UDQ+sdI3QSU7O5vY2Fj69OnjOkrISEtLo7i4mOXLl7uOYkKAlb4JKtnZ2fTt25e4uDjXUULGgAEDiIqKsike4xMrfRM0tm3bxtdff21TOzXUpEkTevXqZaVvfGKlb4LGxx9/DGClXwtpaWnk5uayd6994N2cmZW+CRoLFiwgPj7eLqVcC+np6ajqiR+cxpyOlb4JCuXl5SxYsIBBgwbZpZRrITk5mWbNmjFv3rzqB5uIZv+6TFDIzc1l9+7dDBkyxHWUkBQVFUV6ejrz58+nvLzcdRwTxKz0TVCYP38+IkJ6errrKCFryJAhbN++na+++sp1FBPErPRNUJg3bx5XXnklLVq0cB0lZA0aNAjApnjMGVnpG+f27NnDypUrbWrnLLVs2ZLu3btb6ZszstI3zmVnZ6OqVvp+MGTIED7//HMOHDjgOooJUj6VvogMFpGNIlIgIhNO8Xw/EVkjIqUiMrzKc2Weu2mduKOWMd7mzZtH8+bNSUqq9v4PphqDBw+mrKzMTt00p1Vt6YtIFPACMAToBNwiIp2qDPseuAN49xTf4oiqdvN8DTvF8yaCVZ6qmZ6eTlRUlOs4Ia937940bdqU+fPnu45iglR9H8YkAwWqugVARKYBmcD6ygGq+p3nOTtXzNRIXl4eO3bsYPDgwa6jhIX69etzzTXXMG/ePFQVEXEdyQQZX6Z3WgM/eC0Xetb5qoGI5IrIChG5vkbpTNirPOhYeeaJOXtDhgyhqKiIdevWuY5iglBdHMht77lZ763AJBHpWHWAiNzt+cGQu2vXrjqIZILF3LlzSU5O5oILLnAdJWwMHToUgH/961+Ok5hg5EvpFwFtvZbbeNb5RFWLPH9uAT4Bup9izFRVTVLVpPj4eF+/tQlx27dvZ+XKlQwbZod6/KlVq1YkJyczZ46dN2F+zpfSXwUkikiCiMQANwM+/W0SkWYiEut53AK4Gq9jASayVe6JXnfddY6ThJ/rrruOlStXsn37dtdRTJCptvRVtRQYBywANgDTVTVfRB4RkWEAInKliBQCI4ApIpLv2fwyIFdEvgQWA0+oqpW+AWDOnDm0b9+ezp07u44Sdip/e7IpHlOVqKrrDCdJSkrS3Nxc1zFMgBUXF9OiRQvuuusunnvuOddxwo6qkpCQQNeuXfnnP//pOo6pAyKy2nP89IzsE7nGiZycHI4cOWLz+QEiIgwbNozs7GyKi4tdxzFBxErfODFnzhyaNGlCv379XEcJW8OGDePIkSPk5OS4jmKCiJW+qXPl5eXMnTuXwYMHExMT4zpO2OrXrx9NmjSxs3jMSaz0TZ1btWoVO3bssKmdAIuJiWHw4MHMnTvXbqxiTrDSN3Vuzpw5REVF2VU168CwYcPYsWMHX3zxhesoJkhY6Zs6parMnDmTAQMGcN5557mOE/YyMjKIjo5m1qxZrqOYIGGlb+rUunXr2LRpEyNGjHAdJSKce+65pKWlMWPGDILt9GzjhpW+qVMzZ86kXr163HDDDa6jRIzhw4fz73//m9WrV7uOYoKAlb6pM6rKjBkz6N+/P+eff77rOBEjMzOT+vXrM2PGDNdRTBCw0jd1Jj8/n2+++Ybhw4dXP9j4zXnnncfAgQNtiscAVvqmDs2cORMR4Ze//KXrKBFnxIgRbN26lbVr17qOYhyz0jd1ZsaMGfTr14+WLVu6jhJxrr/+eqKiomyKx1jpm7qxfv161q9fb1M7jjRv3pzU1FRmzpxpUzwRzkrf1Amb2nFvxIgRFBQUkJeX5zqKcchK3wScqvLOO+/Qr18/LrzwQtdxItYNN9xA/fr1effdd11HMQ5Z6ZuAW7VqFZs2beL22293HSWitWjRgqFDh/LOO+9QVlbmOo5xxKfSF5HBIrJRRApEZMIpnu8nImtEpFREhld5bpSIfOv5GuWv4CZ0vPnmmzRo0MDm84PA7bffzrZt2+xyyxGs2tIXkSjgBWAI0Am4RUQ6VRn2PXAH8G6Vbc8DHgauApKBh0Wk2dnHNqGipKSEadOmMWzYMJo2beo6TsS79tpradq0KW+99ZbrKMYRX/b0k4ECVd2iqiXANCDTe4CqfqeqXwFVr986CMhW1b2qug/IBgb7IbcJEfPnz2fPnj02tRMkGjRowMiRI5k9ezaHDh1yHcc44EvptwZ+8Fou9KzzhU/bisjdIpIrIrm7du3y8VubUPDWW28RHx/PoEGDXEcxHv/1X/9FcXExH3zwgesoxoGgOJCrqlNVNUlVk+Lj413HMX6yf/9+5s6dyy233EJ0dLTrOMbj6quvJiEhgTfffNN1FOOAL6VfBLT1Wm7jWeeLs9nWhLgZM2Zw7Ngxm9oJMiLCbbfdRk5ODkVF9s8x0vhS+quARBFJEJEY4GbA15tuLgDSRaSZ5wBuumediQCvvvoql112GT179nQdxVRx++23o6q8/vrrrqOYOlZt6atqKTCOirLeAExX1XwReUREhgGIyJUiUgiMAKaISL5n273Ao1T84FgFPOJZZ8JcXl4eK1as4J577kFEXMcxVSQmJpKamspLL71k5+xHGAm263AkJSVpbm6u6xjmLP3mN7/h9ddf58cff6RZMztLNxjNmDGDkSNHkpWVxdChQ13HMWdJRFaralJ144LiQK4JLwcPHuTtt9/m5ptvtsIPYpmZmVxwwQVMnjzZdRRTh6z0jd+98847HDp0iDFjxriOYs4gJiaGO++8k6ysLL7//nvXcUwdsdI3fqWqTJ48mW7dupGcnOw6jqnGr3/9a1SVl19+2XUUU0es9I1frVy5ki+//JIxY8bYAdwQ0KFDB4YMGcLLL7/M8ePHXccxdcBK3/jVP/7xD8455xxuvfVW11GMj8aMGcO2bdv48MMPXUcxdcBK3/hNYWEh7733HqNHj6Zx48au4xgfDR06lI4dO/LMM8/YXbUigJW+8ZtJkyahqvz+9793HcXUQFRUFPfffz9ffPEFn376qes4JsCs9I1f7N+/nylTpnDTTTfRoUMH13FMDd1xxx3Ex8fz1FNPuY5iAsxK3/jF5MmTOXToEOPHj3cdxdRCXFwc//3f/01WVhZff/216zgmgKz0zVk7evQof/vb30hPT6dbt26u45ha+u1vf0vDhg155plnXEcxAWSlb87a22+/zfbt220vP8Q1b96cu+66i3feeYfCwkLXcUyAWOmbs1JSUsLjjz9Ojx49GDhwoOs45ixVHoR/4oknHCcxgWKlb87KK6+8wpYtW3jkkUfsw1hhoH379tx5551MnTqVrVu3uo5jAsBK39RacXExjz76KH369LGrNIaR//mf/yEqKoqHH37YdRQTAFb6ptaef/55tm3bxuOPP257+WGkdevW3Hvvvbz99tt2Jk8Y8qn0RWSwiGwUkQIRmXCK52NF5H3P8ytFpINnfQcROSIieZ4vu4ZrmNi3bx9PPPEEGRkZ9OnTx3Uc42cPPfQQTZo0YeLEia6jGD+rtvRFJAp4ARgCdAJuEZFOVYbdCexT1YuBZ4EnvZ7brKrdPF92rd0w8fTTT7N//34ee+wx11FMAJx33nmMHz+eOXPm8Pnnn7uOY/zIlz39ZKBAVbeoagkwDcisMiYTeMPzeCYwUOz3/bD17bff8pe//IX//M//pGvXrq7jmAD53e9+R6tWrbj33nvtlophxJfSbw384LVc6Fl3yjGee+oeAJp7nksQkbUiskRE+p5lXuOYqjJ27FgaNGhgH+IJc+eccw7PPvssq1ev5sUXX3Qdx/hJoA/kbgPaqWp34PfAuyLSpOogEblbRHJFJHfXrl0BjmTOxvTp08nOzuaxxx6jZcuWruOYABs5ciRpaWlMnDiRbdu2uY5j/MCX0i8C2nott/GsO+UYEakPNAX2qOoxVd0DoKqrgc3AJVVfQFWnqmqSqibFx8fX/F2YOvHTTz9x33330bNnT37zm9+4jmPqgIjwwgsvcOzYMe6//37XcYwf+FL6q4CNNWLbAAALJ0lEQVREEUkQkRjgZmBOlTFzgFGex8OBRaqqIhLvORCMiFwEJAJb/BPd1LU//vGPbN++ncmTJxMVFeU6jqkjiYmJTJgwgffee4/s7GzXccxZqrb0PXP044AFwAZguqrmi8gjIjLMM+wVoLmIFFAxjVN5Wmc/4CsRyaPiAO8YVd3r7zdhAm/hwoU8//zzjBs3jqSkJNdxTB2bMGECv/jFLxg9ejR79uxxHcecBQm2O+UkJSVpbm6u6xjGy86dO+nSpQvNmzcnNzeXuLg415GMA2vWrKFXr15kZGQwe/Zs+0BekBGR1apa7R6ZfSLXnJGqMnr0aPbv38+0adOs8CNYjx49ePzxx/nwww+ZMmWK6zimlqz0zRk999xzfPTRRzzzzDN07tzZdRzj2H333Ud6ejr33Xcf+fn5ruOYWrDSN6eVk5PDAw88wHXXXcfYsWNdxzFBoF69erzxxhs0adKEzMxMdu/e7TqSqSErfXNK69ev58Ybb+TSSy/lrbfesvlbc0LLli358MMPKSws5Prrr+fo0aOuI5kasNI3P7Njxw4yMjKIi4sjKyuLpk2buo5kgkzv3r156623+Oyzzxg9ejTl5eWuIxkfWembk+zfv5/rrruOHTt2MHfuXNq1a+c6kglSI0aM4IknnmDatGk8+OCDBNuZgObU6rsOYILH7t27SU9P5+uvv2b27Nl2Pr6p1oMPPkhhYSF/+ctfKCkpYdKkSdSrZ/uSwcxK3wAVUzrXXHMNBQUF/POf/2TIkCGuI5kQICI899xzxMTE8Ne//pWjR48yefJkK/4gZqVv2LBhA5mZmRQVFZGVlUVqaqrrSCaEiAjPPPMMcXFxPPbYY+zbt4/XXnuNc845x3U0cwr24zjCzZ49m+TkZA4cOMDChQut8E2tiAh//vOfefrpp5k9eza9e/emoKDAdSxzClb6EerYsWNMmDCBG2+8kU6dOrF69Wquvvpq17FMiHvggQeYN28eP/74I0lJSXz44YeuI5kqrPQj0PLly+nRowdPPvkkv/71r1m6dClt2rRxHcuEifT0dFavXk3Hjh254YYbuOmmm9ixY4frWMbDSj+C7Ny5k3HjxnH11Vdz8OBBsrKymDp1KrGxsa6jmTDToUMHli9fzp///Gc+/PBDLrvsMqZMmcLx48ddR4t4VvoRYO/evUycOJGLLrqIF198kbFjx5Kfn8/QoUNdRzNhLCYmhokTJ/Lll19yxRVXMGbMGC699FLefPNNu+euQ1b6YWzdunWMHTuWDh068L//+79ce+215Ofn8/zzz9O4cWPX8UyEuPTSS1myZAlz586lSZMmjBo1iksuuYSnnnoKuz1q3bPr6YeZ77//ng8++ID333+f5cuXExsby4gRIxg/fjxdunRxHc9EuPLycj744AP+9re/8emnnxITE0NmZiY33ngjQ4cOtZ2Rs+Dr9fR9Kn0RGQz8DYgCXlbVJ6o8Hwu8CfQE9gA3qep3nuf+ANwJlAH3quqCM72WlX7NHDhwgGXLlrFkyRJycnJYs2YNAFdccQWjR49m1KhRNG/e3HFKY34uPz+fKVOm8P7777Nz505iY2NJSUlhwIAB9O/fn549exIdHe06ZsjwW+l77nG7CUgDCqm4Z+4tqrrea8xvgS6qOkZEbgZuUNWbRKQT8B6QDFwIfAxcoqqnndCz0j+14uJiNm/ezLfffsvGjRvJy8sjLy+Pb7/9FlUlOjqaq666iuuuu44bbriBxMRE15GN8UlZWRmff/45s2bNYuHChWzYsAGABg0a0LlzZ7p160aXLl1ITEwkMTGR9u3b2z2aT8Gfpd8b+H+qOsiz/AcAVX3ca8wCz5jlIlIf2A7E47lXbuVY73Gne71wLP3y8nJKSkooKSnh6NGjHD16lCNHjnD48GGKi4s5dOgQBw4c4MCBA+zfv5/du3eza9cudu7cSVFREUVFRT+7bnmHDh3o3r073bt3p0+fPvTq1cvuamXCws6dO1m6dCkrVqxg7dq1rF27ln379p14PioqilatWtGmTRsuvPBC4uPjadGiBS1atKBp06Y0adKEpk2b0qhRIxo1akTDhg2Ji4ujQYMGxMbGEhMTQ3R0dNhdLtzX0vflMgytgR+8lguBq043RlVLReQA0NyzfkWVbVv78Jo1tnfvXvr06VPtuNP9kPNeX/lYVU98eS+Xl5ef9Ljyq6ys7MSfpaWllJaWcvz48RpfdjY2Npb4+Hji4+Np27YtvXv3pk2bNnTs2JHExEQuvvhiu9yxCVvnn38+w4cPZ/jw4UDFv7tt27ZRUFDAt99+y5YtW07sDH3zzTcsW7aM3bt31/jfWf369U/6ioqKol69eif+FJGf/Vn1CzjpceXyqfjyQ6Zr16689957NXofNRUU194RkbuBu4FaX8q3fv36XHHFFb6+XrXrq/4P9V72/gsQFRV1Yp33X5ro6OgTf5FiY2OJjo4mJiaGBg0anPiq3BNp1KgRTZs2pWnTppx77rk0atQo7PZCjKktEeHCCy/kwgsvpF+/fqccU15efuK35cqv4uJiDh8+zOHDhzl27NiJ37JLSko4fvw4JSUlJ3bQKnfOysrKKCsrO2mH7lQ7e1V3BCv5slN5JgkJCTX8r1NzvpR+EdDWa7mNZ92pxhR6pneaUnFA15dtUdWpwFSomN7xNby3Jk2aMH369NpsaowJcfXq1aNZs2Y0a9bMdZSg58t5+quARBFJEJEY4GZgTpUxc4BRnsfDgUVa8aNtDnCziMSKSAKQCHzhn+jGGGNqqto9fc8c/ThgARWnbL6qqvki8giQq6pzgFeAt0SkANhLxQ8GPOOmA+uBUmDsmc7cMcYYE1j24SxjjAkDvp69Y5dhMMaYCGKlb4wxEcRK3xhjIoiVvjHGRBArfWOMiSBBd/aOiOwC/u06Ry20AHZXOyq82HuODPaeQ0N7VY2vblDQlX6oEpFcX06XCif2niODvefwYtM7xhgTQaz0jTEmgljp+89U1wEcsPccGew9hxGb0zfGmAhie/rGGBNBrPQDQETuFxEVkRauswSaiDwtIt+IyFci8oGInOs6UyCIyGAR2SgiBSIywXWeQBORtiKyWETWi0i+iPzOdaa6IiJRIrJWRP7lOksgWOn7mYi0BdKB711nqSPZwBWq2gXYBPzBcR6/E5Eo4AVgCNAJuEVEOrlNFXClwP2q2gnoBYyNgPdc6XfABtchAsVK3/+eBR4EIuJgiaouVNVSz+IKKu6OFm6SgQJV3aKqJcA0INNxpoBS1W2qusbz+CAVJRiQ+1sHExFpA2QAL7vOEihW+n4kIplAkap+6TqLI78C5rkOEQCtgR+8lguJgAKsJCIdgO7ASrdJ6sQkKnbaanaX9RASFDdGDyUi8jHQ8hRPTQT+LxVTO2HlTO9ZVf/pGTORiimBd+oymwksETkHmAX8H1X9yXWeQBKRa4GdqrpaRAa4zhMoVvo1pKrXnGq9iHQGEoAvRQQqpjnWiEiyqm6vw4h+d7r3XElE7gCuBQZqeJ4DXAS09Vpu41kX1kQkmorCf0dVZ7vOUweuBoaJyFCgAdBERN5W1dsc5/IrO08/QETkOyBJVUPtok01IiKDgb8C/VV1l+s8gSAi9ak4SD2QirJfBdyqqvlOgwWQVOy5vAHsVdX/4zpPXfPs6T+gqte6zuJvNqdvztbfgcZAtojkichk14H8zXOgehywgIoDmtPDufA9rgZuB1I9/1/zPHvAJsTZnr4xxkQQ29M3xpgIYqVvjDERxErfGGMiiJW+McZEECt9Y4yJIFb6xhgTQaz0jTEmgljpG2NMBPn/5fMQ/1K/lYIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def standard_norm_pdf(x): \n",
    "    return np.exp(-0.5 * x**2 )/np.sqrt(2*np.pi)\n",
    "\n",
    "x = np.linspace(-5,5,101)\n",
    "plt.plot(x,standard_norm_pdf(x), 'k')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [t-distribution (student's t-distribution)](https://en.wikipedia.org/wiki/Student%27s_t-distribution)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Features"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- This is used to calculate the condidence intervals and do hypothesis testing.\n",
    "- This distribution has one shape parameter called [degree of freedom](https://en.wikipedia.org/wiki/Degrees_of_freedom_(statistics)).\n",
    "- As the degree of freedom increases, the t-distribution becomes closer to the standard normal distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[scipy.stats.t](https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.stats.t.html)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Probability density function (pdf) is visualized. See also [scipy.stats.t](https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.stats.t.html).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3Xl8VNXd+PHPmSX7Tlayr2QngRAQkEUQQRRwacXlqdrFWqvWtrZ1edTWlqebVeuvWuujtu48roiiouyL2SYQAmENSSALS0jIRvaZ8/tjJiErDCHJZDnv1ysvM/eee+c7Ct+5fu+53yOklCiKoijjg8bWASiKoijDRyV9RVGUcUQlfUVRlHFEJX1FUZRxRCV9RVGUcUQlfUVRlHFEJX1FUZRxRCV9RVGUcUQlfUVRlHFEZ+sAevL29pZhYWG2DkNRFGVUyc3NPSOl9LnYuBGX9MPCwjAYDLYOQ1EUZVQRQhyzZpwq7yiKoowjKukriqKMIyrpK4qijCMq6SuKoowjViV9IcRiIcQhIUShEOKRC4y7SQghhRBpXbY9ajnukBDimsEIWlEURRmYi87eEUJogReBq4EyIEcIsVZKub/HOFfgZ0BWl23xwEogAZgIbBBCxEgpjYP3ERRFURRrWXOlnw4USimLpJStwGpgeR/jfg/8GWjusm05sFpK2SKlLAYKLedTFEVRbMCaefqBQGmX12XA9K4DhBBTgGAp5TohxK96HJvZ49jAAcaqKDaTkZHBjh07cHFxYenSpYSEhNg6JEUZkMu+kSuE0ADPAr+8jHPcI4QwCCEMlZWVlxuSogyas2fPct111zFz5kx+/etfc9999xEVFcUf/vAHxtL60r/97W955plnOHjwICkpKaSmpnL06NGLHvfBBx+QkJCARqNRD1WOEtYk/XIguMvrIMu2Dq5AIrBFCFECzADWWm7mXuxYAKSUr0gp06SUaT4+F32KWFGGRW1tLXPnzuXrr7/mr3/9K1VVVRw+fJibbrqJJ554ggcffHBMJX6ANWvWcPPNN7N7924iIyMvOj4xMZGPP/6YOXPmDEN0ymCwpryTA0QLIcIxJ+yVwG0dO6WUtYB3x2shxBbgYSmlQQjRBLwrhHgW843caCB78MJXlKEhpeSOO+7gwIEDrFu3jkWLFgHg5eXFu+++S0BAAM899xyJiYn8+Mc/tnG0A7Nq1SreeOMNfH19CQ4OJi4ujpdeegmtVsvGjRvZvHnzRc8RFxc3DJEqg+miSV9K2S6EuB9YD2iB16WUBUKIpwGDlHLtBY4tEEK8D+wH2oGfqpk7ymjw2muv8fnnn/P3v/+9M+F3EELwzDPPsG/fPn7xi1+wcOFCq66K+/O7zwrYX1F3uSF3Ez/RjaeuT+h3f25uLqtXryYvL4/29namTJnC1KlTuffee3FxceHhhx8G4Morr6S+vr7X8c888wwLFy4c1JiV4WFVwzUp5RfAFz22PdnP2Hk9Xq8CVg0wPkUZdrW1tTzyyCPMnTuX+++/H6PJyDfHv8Fw0oCD1oFrwq4hySeJ119/nbi4OH75y1+yZs0aW4d9SbZv384NN9yAk5MTAMuWLet3nDK2jLgum4piax31+7/97W+cbTnLg5sfJL8yH1e9Ky3GFt7Y/wb/Ff9fPJz2MI899hiPPfYY27ZtG3Bd+0JX5LamrvTHHpX0FaWL06dP89xzz7Fy5Upik2O544s7KKsv439m/w9LI5bS1N7E33f9nbf2v0WbsY1fPPQLnn/+ef74xz+OqpuZc+bM4a677uLRRx+lvb2dzz77rM97E+pKf+xRvXcUpYt//etfNDY28uSTT/LUt09RVFvE36/6O9dHXo9GaHDWO/No+qN8L/57rD60mg0VG/jZz37GV199xZ49e2wdvtWmTJnCLbfcwuTJk1myZAnTpk0b0Hk++eQTgoKCyMjIYOnSpVxzjeq0MtKJkTblLC0tTar5voottLa2EhoaSkpKCr955Tc8sOkB7k+5nx9P7n0F3G5q5+6v7qakroQ3575Jckwyy5cv5+2337ZB5IoCQohcKWXaxcapK31FsXj//fc5efIk991/H3/K/hOR7pF8P/H7fY7VaXQ8dcVTNLQ28HbJ29x999188MEHVFVVDXPUinJpVNJXFItXX32VqKgomiKbKG8o5+FpD6PX6vsdH+UZxU0xN/Fx4ccsu20Zra2tvPPOO8MYsaJcOpX0FQUoKSlh69at3PG9O3i94HUSJyQya+Ksix73w6QfIhB8K78lLS2NV199dcw9pauMLSrpKwp01uJD54VS3lDOPcn3IIS46HH+zv4si1zG2qNrue3O29i7dy+7du0a6nAVZcBU0lfGPSklb775JnPnzmVb4zaCXIKYGzzX6uNvjb2VFmMLdlPs0Ol0vP/++0MYraJcHpX0lXEvPz+fI0eOsOiGReSczOGmmJvQCOv/akzymkSqbyqfn/icBQsX8MEHH6gSjzJiqaSvjHtr1qxBCEHzpGZ0QseKqBXdB9RVwCc/gWdi4IVU2PYMGNu6DflOzHcoaygj/Zp0iouLyc3NHcZPcPkG2lr5iSeeIDk5mZSUFBYtWkRFRcUwRKtcDpX0lXHvk08+4YqZV7C9bjtzg+fi7eh9fmflIXhlPhR8DBHzwCMUNv0e3rkZ2s4vErcgZAGOOkda41rR6/WjtsRzqa2Vf/WrX5Gfn09eXh7XXXcdTz/99DBEqVwOlfSVca24uJg9e/aQtiCNM01nWBy++PzOprPw7i0gTfCjzXDjK/C9NbDsH1C0BT7/eedQJ70T84Pns716O/PmzeOzzz4b/g9ziVatWkVMTAyzZ8/m0KFDNDY28vzzz/PPf/6T+fPnW3UONze3zt/PnTtn1c1vxbZU7x1lXOvojqlJ1ODY6sicwC79c755EmqOw91fgF/8+e1T/su8fdtfIHYpxF0HwNKIpXxR/AWps1L55rffUFRURERExMWD+PIROLl3MD8W+CfBkj/1u3swWys//vjjvPnmm7i7u1vVg1+xLXWlr4xra9asITEpkV2mXcwJmoOT3txqmLJc2PUWzPgJhMzofeDcX4NfIqx/DNpbALhi4hV42HvQFmOu969bt264PsYl69pa2c3N7YKtlfPy8nr9dO2wuWrVKkpLS7n99tv5xz/+MVwfQRkgdaWvjFs1NTXs3LmT2++7nV3Nu1gc1qW0s/kP4OwN8x7p+2CtHq7+Hbx9Exj+DTPuRa/RMydoDptLNxMdE826det44IEHLh7IBa7Ibe1SWivffvvtXHvttfzud78brvCUAbDqSl8IsVgIcUgIUSiE6PW3QAhxrxBirxAiTwixQwgRb9keJoRosmzPE0K8PNgfQFEGavPmzRiNRuwT7HHQOjAr0PIEbsVuOLoJZtwH9q79nyByAYTOhm//X+dsnquCr6K+tZ6p86ayZcsWzp07Nwyf5NLNmTOHNWvW0NTURH19fb/3IC52pX/kyJHOsZ9++imxsbHDEr8ycBe90hdCaIEXgauBMiBHCLFWSrm/y7B3pZQvW8YvA54FOi6bjkopUwY3bEW5fOvXr8fFxYUSjxLSvdNx1Dmad+z8O9i7w7QfXPgEQsDM++G9lbD/U0i6mSsmXoG91h7nZGdaWlrYuHFjv6UTW+raWtnX13fArZUfeeQRDh06hEajITQ0lJdfVtd1I5015Z10oFBKWQQghFgNLMe87i0AUsquC3w6A+rJFGVEk1Kyfv16Zlw5gxPNJ/hBoCXBnzsDBz6D9B+Dg/vFTxR9DXhFQtbLkHQzTnonZgTM4JD2EC4uLnz55ZcjMumD+Qbs448/flnn+OijjwYpGmW4WFPeCQRKu7wus2zrRgjxUyHEUeAvwINddoULIXYLIbYKIa7s6w2EEPcIIQxCCENlZeUlhK8oA3P06FFKSkrwT/UHYE6QZdbOntVgaocp37PuRBoNTL0LynLgjLnUMT94PidbT5I2M41NmzYNQfSKMnCDNntHSvmilDIS+A3w35bNJ4AQKWUq8AvgXSGEWx/HviKlTJNSpvn4+AxWSIrSr6+//hqAxshGoj2jCXAJAClh91sQNA18L6E2nfxdEBrIexeg895AQEoAhw8fpqysbNDjV5SBsibplwPBXV4HWbb1ZzWwAkBK2SKlrLL8ngscBWIGFqqiDJ6vv/6asLAwjmqPnp+bX7EbKg9C6h2XdjJXf4haCPn/ByYj/s7+RLhH0BZuvrmrrvaVkcSapJ8DRAshwoUQdsBKYG3XAUKI6C4vlwJHLNt9LDeCEUJEANFA0WAErigDZTQa2bJlC3Ez4jBi5MogS9Vx/6eg0UHcAGrwk1dCXTkczwBg5sSZFDkV4e3trZK+MqJcNOlLKduB+4H1wAHgfSllgRDiactMHYD7hRAFQog8zGWcOy3b5wD5lu0fAvdKKasH/VMoyiXIz8+ntrYWl0kuOOocSfZONpd29n8K4XPAyevSTxp9DWjt4cDngPlBrVbZSvKMZDZu3Ki6biojhlUPZ0kpvwC+6LHtyS6//6yf4z4C1O19ZUTZtm0bALVBtUzxm2JeEvHkXjhbDLMfGthJ7V0g8irzzJ/FfyTNLw2dRodnoidln5dRWFhIdHT0xc+jKENMtWFQxp2tW7cSGhZKhb6C6f7TzRv3f2q+GRt73cBPHHc91JVBxW6c9E5M8Z1CQ2gDABs3bhyEyIfOQFsrf/DBByQkJKDRaDAYDN32/fGPfyQqKopJkyaxfv36oQpduUQq6SvjipSSbdu2ETUlCoD0gHTzjgOfQ+gsc+uFgZq0BITWfLUPzAiYQZlDGf7+/uzYseNyQx8Wl9paOTExkY8//pg5c+Z0275//35Wr15NQUEBX331Fffddx9Go3GowlYugUr6yriyf/9+qqqqcI51xtXOlVjPWKgtg8oDEHPN5Z3cyQtCZ8IR83TQNP80hBDETI0ZkUl/MForx8XFMWnSpF7bP/30U1auXIm9vT3h4eFERUWRnZ092B9BGQDVcE0ZV7Zu3QpATWAN0/ymodVoodBSeolaeIEjrRS1ADb8FupOkDghEQetA+6T3Dm27hilpaUEBwf3OuTP2X/mYPXBy3/vLmK9YvlN+m/63T+YrZX7Ul5ezowZ57uTBgUFUV5+oZneynBRSV8ZV7Zt24b/RH+qnKrOl3YKvwG3QPAZhGZhUQvNSf/oJvSptzPZZzIlwSUA7Ny5k5UrV17+ewyCrq2VgQu2VlbGFpX0lXFDSsnWrVuJnhbNWXGWdP90c3fMoq2QsMLcQO1y+SWCix8UboDU25nqP5Ws8ixcXFzYsWNHn0n/QlfktjbQK/3AwEBKS893bykrKyMwsFf3FsUGVE1fGTcKCws5efIkbrFueNh7EOURZe6Z01IHUVcPzpsIYW65fHQTmIyk+aWBFmJTY0dUXX+wWiv3Z9myZaxevZqWlhaKi4s5cuQI6enpQ/FRlEukkr4ybnz77bcANAQ2kOKTYl7PtXCDecZNxNzBe6OoBdBcA+W7SPZJRq/R4xXnRX5+PjU1NYP3Ppeha2vlJUuWDLi18ieffEJQUBAZGRksXbqUa64x3wxPSEjgu9/9LvHx8SxevJgXX3wRrVY7mB9BGSBV3lHGjYyMDFzdXDnjdoZUv1TzxqObzQ3WrGmjbK3IqwABRzdhHzyNJO8kjoccR0pJRkYGS5YsGbz3ugyD0Vr5hhtu4IYbbhiy8yuDT13pK+NGZmYmMZNjEBpBqm8qNNfBiTxz64WLOHiyjue+Ocyr24uoaWy98GAnL3Ntv8R8EzTNP40zvmfQarUjqsSjjE8q6SvjQn19PXv37sVzkid2GjsSJiTA8UyQJgibfcFj/y/nONe9sIO/bzzCH9YdYOGz29hfUXfBYwi/0ny/oK2ZqX5TwQ5ikkbmfH1lfFFJXxkXDAYDJpOJ9uB2ErwTsNPama/EtXbm8k4/dhae4ZGP93JF5AR2P3E1nz8wG71W8KM3DRe+4g+bDe3NUJ5LsncyAoF/vD85OTm0tbUNwSdUFOuopK+MCxkZ5pbHZ3zPmEs7ACU7IDAN7Jz6PKap1cjDH+wh0seFl++YiqezHYmB7rx8x1RO1zfzh3UH+n/D0JmAgJLtuNi5EOUZhSZUQ1NTE/v27RvkT6co1lNJXxkXMjIyCIsKA0eY4jvlfD3/AqWd13YUcaK2mf+5IQln+/NzHiYHe3D3rHA+2lVGQUVt3wc7eoJ/kvmLBUjxSaHa19xVPCsra9A+l6JcKpX0lTFPSklmZiYT4ycCkOKbctF6flOrkdd2FLMwzpf08N799X86PwoXex0vb73AmkBhV0JpNrQ1M9lnMq0erXhO8FRJX7EplfSVMe/o0aOcOXMGXYSOSPdI3O3dL1rP/2R3OWcb2/jRlRF97nd31LNyWjBf7D1BRU1T328cfiUYW6AshxRf83MB4UnhIybpr1q1ioSEBJKTk0lJSemM6/nnn6exsXHQ3icsLIwzZ84M+PgtW7Zw3XWX0fJ6CLm4uFx0TEfbaoAnn3ySDRs29Dt2zZo17N+/f9Di64tVSV8IsVgIcUgIUSiEeKSP/fcKIfYKIfKEEDuEEPFd9j1qOe6QEOIy2xgqyqXLzMwEoMavxnyVDxes50sp+c+3xSRMdOvzKr/D964IQ0rJW5nH+h4QcoW5R3/JDkJcQ/C098Ql0oUDBw7Y/CGtjIwMPv/8c3bt2kV+fj4bNmzobAY32En/Ug1XC2ZbtHp++umnL/g084hI+pY1bl8ElgDxwK1dk7rFu1LKJCllCvAX4FnLsfGY19RNABYDL3WsmasowyUjIwNnF2eMfsbu8/P7Ke0UVNRx+FQDt08PNT+1249gLyfmT/Llk13lGE19LIfo6AG+CXA8AyEEk30m0xLYAkBOTs6gfLaBOnHiBN7e3tjb2wPg7e3NxIkTeeGFF6ioqGD+/Pmd7ZV/8pOfkJaWRkJCAk899VTnOcLCwnjqqaeYMmUKSUlJHDxo7hRaVVXFokWLSEhI4Ic//GG3pSJXrFjB1KlTSUhI4JVXXunc7uLiwi9/+UsmT55MRkYGX331FbGxsUyZMoWPP/64z8/wn//8hxtvvJHFixcTHR3Nr3/968597733HklJSSQmJvKb35zvbdTzfcLCwnj00UdJSUkhLS2NXbt2cc011xAZGcnLL78MQENDAwsWLOj8nJ9++ulF//32bFvd4a677uLDDz8E4JFHHiE+Pp7k5GQefvhhvv32W9auXcuvfvUrUlJSrFrEZiCseSI3HSiUUhYBCCFWA8uBzq8jKWXXScvOQMd/5eXAaillC1AshCi0nC9jEGJXFKtkZmYSkRiB1EiSfJLM8+elCUKv6HP8p3nl6LWCa5P8L3ruFamBbDx4mqyiKmZG9bEAS8gMyHsXjO1M9p3MRj9zG+esrCyuvtrc7+ehhx4iLy9v4B+wDykpKTz//PP97l+0aBFPP/00MTExLFy4kFtuuYW5c+fy4IMP8uyzz7J582a8vc2fZ9WqVXh5eWE0GlmwYAH5+fkkJycD5i+LXbt28dJLL/HMM8/w6quv8rvf/Y7Zs2fz5JNPsm7dOl577bXO93399dfx8vKiqamJadOmcdNNNzFhwgTOnTvH9OnT+dvf/kZzczPR0dFs2rSJqKgobrnlln4/R15eHrt378be3p5JkybxwAMPoNVq+c1vfkNubi6enp4sWrSINWvWsGLFim7v0yEkJIS8vDx+/vOfc9ddd7Fz506am5tJTEzk3nvvxcHBgU8++QQ3NzfOnDnDjBkzWLZsWb8XBP21re6qqqqKTz75hIMHDyKEoKamBg8PD5YtW8Z1113HzTfffPH/yANkTXknECjt8rrMsq0bIcRPhRBHMV/pP3iJx94jhDAIIQyVlZXWxq4oF9XY2MiePXvwjPHEVe9KmFuY+eYqwlze6cFokqzdU8HcGF88nOwuev6FcX642Ov4ZHc/veJDZkDbOTi1jxSfFLTOWkKiQmxe13dxcSE3N5dXXnkFHx8fbrnlFv7zn//0Ofb9999nypQppKamUlBQ0K38cOONNwIwdepUSkpKAHP76jvuuAOApUuX4unp2Tn+hRdeYPLkycyYMYPS0lKOHDkCgFar5aabbgLg4MGDhIeHEx0djRCi81x9WbBgAe7u7jg4OBAfH8+xY8fIyclh3rx5+Pj4oNPpuP322zvXRe76Ph062konJSUxffp0XF1d8fHxwd7enpqaGqSUPPbYYyQnJ7Nw4ULKy8s5depUvzF1bVvt5ubWZ9vqjph/8IMf8PHHH3e2uB4Og9Z7R0r5IvCiEOI24L+BOy/h2FeAVwDS0tL6+P9kRRmY3bt3YzQaMQYbSfBOQCM0UJYNfgng4NZrfHZxNafqWvjvpROtOr+jnZZFCX6sLzjJ/9yYhF7b4zoq2LIG7/FMEtLuRCd0+Mf6k5WVhZQSIcQFr8iHklarZd68ecybN4+kpCTeeOMN7rrrrm5jiouLeeaZZ8jJycHT05O77rqL5ubmzv0d5SGtVkt7e/sF32/Lli1s2LCBjIwMnJycmDdvXue5HBwcBtSQreP9rY2hr/fpOIdGo+l2Po1GQ3t7O++88w6VlZXk5uai1+sJCwvr9u9gIHQ6HdnZ2WzcuJEPP/yQf/zjH2zatOmyzmkta670y4Guy/0EWbb1ZzWwYoDHKsqgys3NBeCsz1mSvJPAZIIyQ7+zdtYXnMRBr2FhnJ/V77Eo3p+65nZySqp77/QIBrcgKM3EUefIJK9JaMO0VFZWdl4Z28KhQ4c6r7LBXCYJDQ0FwNXVtbOHfl1dHc7Ozri7u3Pq1Cm+/PLLi557zpw5vPvuuwB8+eWXnD17FoDa2lo8PT1xcnLi4MGDnTfYe4qNjaWkpKSzpv3ee+9d0mdLT09n69atnDlzBqPRyHvvvcfcuQPvolpbW4uvry96vZ7Nmzdz7Fg/N+4trGlb3dDQQG1tLddeey3PPfcce/bsAbr/ux8q1lzp5wDRQohwzAl7JXBb1wFCiGgpZcefoKVAx+9rgXeFEM8CE4FoQC2UqQwbg8GAt583Gg8Nid6JUHnQ3D+/4wq8CyklGw+eYnaUN4521l91zonxxk6n4Zv9p5gZ2VddfzocywApSfZJZp+f+YncrKwswsPDB/zZLkdDQwMPPPAANTU16HQ6oqKiOm+s3nPPPSxevJiJEyeyefNmUlNTiY2NJTg4mFmzZl303E899RS33norCQkJzJw5k5CQEAAWL17Myy+/3LmubtflFLtycHDglVdeYenSpTg5OfW7kEt/AgIC+NOf/sT8+fORUrJ06VKWL19u9fE93X777Vx//fUkJSWRlpZGbOyFV1jr2rba19e3z7bV9fX1LF++nObmZqSUPPvsswCsXLmSH/3oR7zwwgt8+OGHfPPNNwDce++9A46/FynlRX+Aa4HDwFHgccu2p4Fllt//DhQAecBmIKHLsY9bjjsELLnYe02dOlUqymCJi4uTqfNSZeJ/EuXpc6elNPxbyqfcpDxT2Gvs4ZN1MvQ3n8t3Mo9d8vvc/e9sOfvPG6XJZOq9M/Nf5vc8e0yuLVwrE15NkPYO9vKhhx4awCdSlL4BBmlFPreqpi+l/AL4ose2J7v8/rMLHLsKWGXN+yjKYKqvr+fgwYPMnDETR2dHfJx8zDdxnSaAV++HrjYePA3AVbG+l/xeC+P82HTwNIdPNTDJ37X7zpCOun4WiaFpCJ0gPH7kPKSljC/qiVxlzNq9ezdSSpoDms31fDAn/aD0PtfD3XTgNImBbvi7O1zyey2IM39RbD50uvdO3wSwc4XSTELdQnHVu+IR7dE5pU9RhpNK+sqYZTAYADjnd85cz2+shqojENx7rdaGlnZyj59lXsylX+UD+Lk5MMnPlR1H+mg3oNVBUBocz0QjNCR4J2AKNNHU1MSBAxfo1KkoQ0AlfWXMMhgM+E70ReeuM1/pl1megu0j6eceO4vRJLkicsKA329WlDfZJdU0t/XxeH/IDDhVAM21JHknUetT2xmjogwnlfSVMctgMBAQG2C+up6QYC7tCC1MTO01Nru4Cp1GkBriMeD3uzLam9Z2E4aSs713Bk8HJJTlkOSdhM5Ph5OzU+eUUkUZLirpK2NSTU0NR44cwT7Ungj3CJz0TlCaZe5xb+fca3x2cTVJQe442Q38ecX0cC/0WsGOwj5KPEFp5uZrpdkkeiciNIKQ+BB1pa8MO5X0lTFp165dAJzzP2cu7RjboXxXn6Wd5jYje0prL9hR0xrO9jpSQzzZUdhHKxF7V/NTwKVZ+Dj54O/sj0uEC3l5eWr5RGVYqaSvjEkdV9DGQKP5Ju7pAnMPnD4eytp9vIZWo4npl5n0Aa6M8qagoo7qc32snxs83fw0sMlIkncSbQFttLS0UFBQcNnvqyjWUklfGZMMBgP+wf7oXCw3cUstD4L30X4hu7gaIWBq6OUn/ZlRE5DSfI+gl+Dp0NoAp/eT6J1IY0BjZ6yKMlxU0lfGJIPBgE+MDw5aB6I8o8xJ38UfPEJ6jc0uqSLO3w13R/1lv29SoAcOeg1ZxX304ekoLZVmkeSdhJ2vHS5uLirpK8NKJX1lzKmqqqK4uBhdiI64CXHoNXpzZ83gab0eymptN5F77Oxl1/M72Ok0TAnxJLuvpO8Rav7iKc0mfkI8GqEhMDZQzeBRhpVK+sqY05FEG/wazPX8htNwtqTPev6+ilqa2wannt8hPdyL/SfqqG3qcYNWCPPVfmkWznpnIj0icQxzZM+ePbS0tAza+yvKhaikr4w5HeUSXXDPen7vmTsdV+TTBjHpTw831/Vzj/VV4plu/gKqP0WidyLNAc20tbWxb9++QXt/RbkQlfSVMcdgMBAQGoDWWWt5EjcbtHYQMLnX2OziaqJ8XfB2se/jTAOTGuKBXiv6qetb/m+jLJsk7yRMQabOmBVlOKikr4w5BoMBr2gvPO09CXQJNF/pB0wGffdGakaTJKeketDq+R0c9FomB3n0XdcPSAatfefNXL23HlcPV5X0lWGjkr4yppw6dYrS0lJEsDA/+Wpsg4rdfZZ2Dp6so765fVDr+R3Sw73YW1ZLY2uPLpo6e3MbiNJsojyjcNA5EDApQCV9ZdhYlfSFEIsn15gSAAAgAElEQVSFEIeEEIVCiEf62P8LIcR+IUS+EGKjECK0yz6jECLP8rN2MINXlJ46buI2+jeaSzsn90J7c59P4nbW88MGP+lPj5hAu0my61hN753B6VCxG73JRJxXHHZhduzbt++y111VFGtcNOkLIbTAi8ASIB64VQgR32PYbiBNSpkMfAj8pcu+JilliuWn97LwijKIDAYDQgjsQ+zNM3fKLDdx+0j6WUXVBHs5MtHDcdDjmBrqiUZc4CEtYyuc2EOidyJN/k20t7eTn58/6HEoSk/WXOmnA4VSyiIpZSvmhc+7LTgppdwspWy0vMzEvAC6ogy73Nxc/MP90TpabuKWZoF7MLhN7DZOSkl2STXpYQNvpXwhLvY6EgPdrXpISxtiXo9XlXiU4WBN0g8ESru8LrNs688PgC+7vHYQQhiEEJlCiBUDiFFRrGYwGPCI9CDYNRgPBw8ozemz9cLRygaqz7UOST2/Q3qYF7tLa2hp79Ff38UXPMPheKb5Zq6XHjcvN5X0lWExqDdyhRB3AGnAX7tsDpVSpgG3Ac8LISL7OO4eyxeDobKyjw6FimKFiooKKioqzjdZqy2HurI+H8rquAIf7Jk7XU0L96K13UR+WW3vncHToTSbIJdAPB088Y3xVUlfGRbWJP1yILjL6yDLtm6EEAuBx4FlUsrOxwullOWWfxYBW4BeK1hIKV+RUqZJKdN8fHwu6QMoSoeOm7jtge3n5+eDuf1CD9nF1fi62hM6wWnI4um4Qdzn1M3gdDh3GlFzjETvRLQhWgoKCmhsbOw9VlEGkTVJPweIFkKECyHsgJVAt1k4QohU4F+YE/7pLts9hRD2lt+9gVnA/sEKXlG6MhgMaDQaHEIcLPX8HNA5gF9St3FSSrKKzPPzRR8LpA8WL2c7on1dyCm5wENapeaHtJoDmjGZTOTl5Q1ZPIoCViR9KWU7cD+wHjgAvC+lLBBCPC2E6JiN81fABfigx9TMOMAghNgDbAb+JKVUSV8ZEgaDAb9wP+wc7Ij1ijXfxJ04BXR23caVnW3iZF3zkNbzO0wL9yK3xLz+bje+cWDnCqVZJHon4hDm0PkZFGUoWbU2nJTyC+CLHtue7PL7wn6O+xZI6mufogwmKSUGgwH3ye5Ee0bjIIETe+CKn/Yae76ePzQzd7pKD/Pi3azjHDxZR8JE9/M7NFrzEoql2SQu/G/0nnrcfdxV0leGnHoiVxkTysrKOH36NG0T28ylnRN5YGrr56GsKjyc9ET7ugx5XB2N3Pqu60+H0wV4oSPQJZAJ0RNUm2VlyKmkr4wJHVfImhCNeeZOaZZ5Rz+dNdPDvNBohq6e3yHQw5FAD8d+6vrpIE1QnkuydzIEwYEDB2hoaBjyuJTxSyV9ZUwwGAxodVocgh1I9kk2N1nzDAeX7rPBTtU1U1LVOKRTNXtKD/ciu/gsUvao6welAcJc4vFOxBhoRErJ7t27hy02ZfxRSV8ZEwwGAz7hPrg6uxLmGmpO+heYnz99GOr5HaaFeXGmoYWSqh7TMR3cwTfe/GSuTxKOYeZ2EKqurwwllfSVUa/jJq5TuBMJExLQ1pXBudP91vNd7HXEBbgOW3zp4Z6d791LyHQoyyHWIwZ7D3s8/DxU0leGlEr6yqhXUlJCdXU1rRNbLfX8/pusZRdXMzXUE512+P7oR/q44OVsR3bx2d47g6dDSx2OZ48R7RmNR5RK+srQUklfGfU6kqR9mP35Jmt2LubSSRfV51o5fKphWOv5AEII0kI9+7+ZC53N14yBRg4fPkxtbR+tGxRlEKikr4x6BoMBnV6HfaD9+TVxA6ea58J30ZF0h+OhrJ7Sw704Xt3IqboePfM9w8HZp/PJXE2w+a/krl27hj1GZXxQSV8Z9QwGA96R3vi7+eOnc4ZT+/q8iZtdXI29TkNSkHsfZxla6f3N1xfC0nzN/GSuupmrDDWV9JVRzWQykZubi12onbmeX7HLPPe9n6SfGuKBvU7bx5mGVnyAG8522v6br1UfJULvhqunK54BnirpK0NGJX1lVDt69Ci1tbW0T2wnySepy0NZU7uNq29uo6CidlhaL/RFp9Uwpd+6vvkLSluWS8KEBFwiXFTSV4aMSvrKqNaRHB3DHc/P3PGJBUfPbuNyj53FJG1Tz+8wLcyLQ6fqqW1s674jIAU0+s75+sZAI0VFRVRX9/EFoSiXSSV9ZVQzGAzo7fU4TnQkwTMOyvpeKSuruBqdRpAa4mGDKM2mhXkhJRiO9UjmegeYmNJ5M9cuzNwVVPXhUYaCSvrKqGYwGPCK8CJiQgSu9aeg6Wy/9fykIHec7KxqLDskUkM80GtF/83XKnaR5DEJx1B1M1cZOirpK6OW0Whk165daEO05tJOWd8PZTW2tpNfVsOMCNvU8zs46LUkB3mQ3d98/fZm/OpO4efth2egupmrDA2V9JVR6/DhwzQ0NCCChXl+/rEMcPSCCdHdxu06VkObUdq0nt9hWpgXe8tqaWrtsVi6pRuoKDM3X3MMd1RJXxkSViV9IcRiIcQhIUShEOKRPvb/QgixXwiRL4TYKIQI7bLvTiHEEcvPnYMZvDK+db2Ja076OyHkCtB0/2OdWVSFViNIC7N90k8P96TdJNld2qMlg1sAeIR0PplrCjRx/PhxKisrbROoMmZdNOkLIbTAi8ASIB64VQgR32PYbiBNSpkMfAj8xXKsF/AUMB1IB54SQniiKIPAYDCgd9DjGuhKjN4NzhZD6BW9xmUWVZEU6I6Lve3q+R2mhnohBOT014enNJukCUk4hpvr+upmrjLYrLnSTwcKpZRFUspWYDWwvOsAKeVmKWVH39hMIMjy+zXAN1LKainlWeAbYPHghK6MdwaDAc9IT+K849CX5pg3hszsNqap1cieEVDP7+DuqGeSnyvZJX103AyeDvUnSLBzxyFUrZmrDA1rkn4gUNrldZllW39+AHw5wGMVxSrt7e3s3r0bESzMi6YczwC9MwQkdxu36/hZcz0/wvalnQ7Tw70s9xlM3XdYbkC7ndxPpF8kHsGq46Yy+Ab1Rq4Q4g4gDfjrJR53jxDCIIQwqBqmYo39+/fT1NSEPlTPZN/JcOxbCJ4GWn23cZ31/NCRU1WcFu5FU5uRgoq67jt8E8xfXJa6vl2onUr6yqCzJumXA8FdXgdZtnUjhFgIPA4sk1K2XMqxUspXpJRpUso0Hx+fnrsVpZeuN3FTXMPhVEGv0g6Yk35ioDuuDvpe+2wl3XJDOafnfH2tztw+wtJ8TROsoby8nBMnTtggSmWssibp5wDRQohwIYQdsBJY23WAECIV+BfmhH+6y671wCIhhKflBu4iyzZFuSwGgwE7ZzuCw4LxryoBJIT2rufnldYwYwSVdgB83RwIm+BEZlE/df2T+0h2j1I3c5UhcdGkL6VsB+7HnKwPAO9LKQuEEE8LIZZZhv0VcAE+EELkCSHWWo6tBn6P+YsjB3jask1RLovBYMA5zJkUvxTzVE2N3rLQ+Hkd9fwZNmqydiEzo7zJKq7uo64/A6SRSU0NuIe7IzRClXiUQWXVHDYp5RfAFz22Pdnl94UXOPZ14PWBBqgoPbW2trJnzx7cFriR4psC374JE1NB79htXGZRFRoBaWEjp57fYXaUN+9mHSe/rIapoV3+T8TyxaUvz2Vy0GSKg4pV0lcGlXoiVxl18vPzaW1txTHckcmesVCxu8/5+duPnGFysMeIqud3uCJiAkLAjiM9SjyOHuATB6XZpPikoAnRkGPIQUppm0CVMUclfWXUycoy98x3j3InrrEBTG0QOqvbmNqmNvLLargyytsWIV6Up7MdiRPd2Vl4pvfO4HQ4nkWqz2Qcwhw4feo05eW95j8oyoCopK+MOpmZmTh4OpAyKQX98UwQml6dNTOOVmGSMDt65M4GmxXlza7jZznX0t59R/gcaKllslHbuXxidna2DSJUxiKV9JVRJzMrE/twe3M9v3gbBEw2l0W62FFYibOd1qb98y9mdpQ37SbZu9Vy2JUAeJTvIjY5Fo1e0/l/N4pyuVTSV0aV6upqCo8U4hDhQIpnrHnRlPA5vcbtOHKGGRET0GtH7h/xtDBP7HQadvQs8bj6mVf/Kt5GWmAaTqFOZGZm2iZIZcwZuX8jFKUPHWUOx0hHJre0mOv5PZJ+aXUjJVWNzI4emfX8Dg56LdPCPPuu64fPgeMZpExIxD7cnuycbNrb23uPU5RLpJK+MqpkZmaCgMjESLzL88zz80O6z9zpuHK+coQnfTDX9Q+erKeyvqX7jvA50NZIikmHY6QjzU3N7N271zZBKmOKSvrKqJKVlYVzkDNTQ6aa6/lB08DOuduYHUfO4O/mQKSPi42itN5sy+yiXlf7YbMBQdipw/jH+QOour4yKFTSV0YNKSWZWZnYhdsxxTMWTuzpVdoxmiQ7j55hdrQ3QggbRWq9xInuTHC2Y/Oh0913OHpCwGREyXamx0/Hzt1O1fWVQaGSvjJqFBYWUnO2BsdIR6a2mkCaeiX9vNKz1DS2MW/SyJ2q2ZVGI5g7yYethysxmno8gBU+B0qzSPVJxD7Cnp0ZO20TpDKmqKSvjBodV7p+sX5EnDwIOkdzeaeLjQdOo9MIrhzB8/N7mj/Jl5rGNvJ6LqEYPheMraSYdDhFOlF4uJDqatW6Srk8Kukro0ZWVhZaBy2zUmehKdlubr2gs+s2ZtPB06SFeeLuOPJaL/RnTowPWo1g08EeJZ6QGaDRkVBVikuU+f6EekhLuVwq6Sujxo6MHTiEOTBtQixUHuhV2imvaeLgyXoWxPrZKMKBcXfUMzXUk00HeywgZO8CgWk4FG1lWto0EOpmrnL5VNJXRoWmpib25e/DMcKRqU3N5o1R3Zu7dlwpXxXnO9zhXbarYn05cKKOk7XN3XdELYSKPK4ISsEh0EHV9ZXLppK+Mirs3r0bY7sRz0meTCrfCy5+4JfYbcymA6cIneBEhLdzP2cZueZPMn9R9ZrFE7UAkKS1SRwjHcnMzMRkMvU+gaJYSSV9ZVTouIk7PX06uqLN5ivgLlMyG1vb2Xm0iqtifUfFVM2eYvxcCPRw7F3XD0gBJ29STh3BJcqF+tp6jhw5YpsglTHBqqQvhFgshDgkhCgUQjzSx/45QohdQoh2IcTNPfYZLatpda6opSiXasu2Leh99MwJjIDmGssV8HnfFlbR2m4adfX8DkIIFsb5sv1IJY2tXdotaDQQtQCno1uYnDYZQM3XVy7LRZO+EEILvAgsAeKBW4UQ8T2GHQfuAt7t4xRNUsoUy8+yPvYrygVJKdm+czvO0c5Mbagzt1KOmN9tzPqCk7ja65gWPvJWybLWNYn+NLeZ2Hqoxw3dqIXQeIarJiWgcdSwfed22wSojAnWXOmnA4VSyiIpZSuwGljedYCUskRKmQ+oYqMy6AoLC6k5U4PbJDcSy/IgMA2czi8x2GY08c2BUyyM98Nep7VhpJcnPcwLL2c7vtx3svuOyKsAwfSWFpyinNi8bbNN4lPGBmuSfiBQ2uV1mWWbtRyEEAYhRKYQYsUlRacowI4dOwCYmp6MXfnuXrN2MouqqGlsY3Givy3CGzQ6rYar4/zYdPA0Le3G8zucvWFiCinl+3CJcaHoUJF6SEsZsOG4kRsqpUwDbgOeF0JE9hwghLjH8sVgqKys7H0GZVzbuHUjWmctV4eFAbJX0v9y30mc7LTMjRk9T+H2Z3GSPw0t7ew40qMBW9TVOJUZSEyLA2DnTjV1UxkYa5J+ORDc5XWQZZtVpJTlln8WAVuA1D7GvCKlTJNSpvn4jP6/uMrg2rp9K07RTsysPQ1OE2BiSuc+o0nydcFJ5sf64qAfvaWdDrMivXF10PUu8UQtBGni6vhIhFawacsm2wSojHrWJP0cIFoIES6EsANWAlbNwhFCeAoh7C2/ewOzgP0DDVYZfyorKykrKsMz1pOEogyIWQKa88k9p6SaMw2tLBnlpZ0OdjoNC+P8+Gb/KdqMXW6RBaWBkzezGqtxDHdkw9YNtgtSGdUumvSllO3A/cB64ADwvpSyQAjxtBBiGYAQYpoQogz4DvAvIUSB5fA4wCCE2ANsBv4kpVRJX7FaRxljWko0upZaiL222/4v9p7AXqfpfLhpLFic6E9tUxvfHq06v1GjhUlLSC3JxmWSCwf2HKCpqcl2QSqjllU1fSnlF1LKGCllpJRylWXbk1LKtZbfc6SUQVJKZynlBCllgmX7t1LKJCnlZMs/Xxu6j6KMRV9t+gqhEywNcTN31ewyVbPNaOKzPRUsjPfD2V5nwygH19wYH1wddKzZ3aOKGrsUp+Y6kidHYGw3qj48yoCoJ3KVEW3Ttk04Rjgy+/QhiJwPdk6d+7YdruRsYxs3pFzKZLKRz0Gv5brkAL7ad5JzLV0e1IqYB3onlkZ5gYD1m9fbKkRlFFNJXxmxGhsbObrvKBNiPYk6WwaTupd21uRV4OmkZ84YmLXT04qUQJrajHy9v8sNXb0jRF7F1XVF2Afaq6SvDIhK+sqIlZWVhcloYlpcAAIBMYs799U3t/F1wUmuS56InW7s/TGeFuZFoIcjn+yu6L4jdimxZyvwjPWgwFBAe3t73ydQlH6Mvb8typjx0VcfgYAV/kYIng4u56/o1xecoqXdxIrUsVXa6aDRCFakTmTHkUpO13VptxyzGK3QMDXen9amVvbs2WO7IJVRSSV9ZcT6ZuM3OIQ6sKi+BOK7t236KLeMEC8npoR42Ca4YXBDahAmCWv3dLnad/KC0FksD2gD4OP1H9soOmW0UklfGZEaGxsp3FNISII3fkYjxJ/v4FF85hwZRVXcMi14VLZRtlaUrwuTg9z5wFCGlF0WTU+4gevaTmDna8eXG760XYDKqKSSvjIibdi6AVO7iXkRdhAyE9zPl3FW5xxHqxF8Z2qQDSMcHivTQzh0qp5dx7ssmh6/nAAT+Cd4UJCl6vrKpVFJXxmRVq9bDRq43asWEm/s3N7abuJDQxkLYn3xdXOwYYTDY9nkibjY63gn6/j5jc7eEDGPWZH2tDa28m32tzaLTxl9VNJXRqTtW7fjEu7ETGGE+POdvDccOEXVuVZunR5iw+iGj7O9jhWpE1mXf4KaxtbzOxJv4nY/8xO5b69920bRKaORSvrKiFNXV0f5/nLiY5zRhV0JLudbLLybdZxAD0fmRI+9ufn9uS09lJZ2Ex/t6vKEbuxSFtmDY5ADGzdttF1wyqijkr4y4qz+ajXSJLk2yNittHP4VD07Cs9w2/QQtJqxewO3p/iJbqSGePBO5jFMJssNXUcP9FFXExvjREleCc3NzRc+iaJYqKSvjDgfffkRQiu4O0DXbdbO6zuKcdBruC19fJR2urprZhhFZ86x+VCXhdOTbuLaUImpzcR769+zXXDKqKKSvjLiZG/NxjvCkZCkpeBonodf1dDCx7vLuXFKEJ7OdjaOcPhdmxRAgLsDr24vPr9x0lJ+GOoMAj5Y94HtglNGFZX0lRElvyifmuIaZkfrIOW2zu3vZB2ntd3E92eF2S44G9JrNdw9K4yMoir2lddaNjoQlvZdvMMcyN6uOm4q1lFJXxlRXnr/JQC+H+fV2Ua5uc3ImxnHmBvjQ5Svqy3Ds6lbpoXgbKfltR1drvZT72B6lJ6qw9UcKjtku+CUUUMlfWVE+ear9di5arl20fc6V8j6v5xSzjS0cO/cXssrjyvujnpumRbC2j0VlFY3mjcGpHBnchCY4KUPXrJtgMqoYFXSF0IsFkIcEkIUCiEe6WP/HCHELiFEuxDi5h777hRCHLH83DlYgStjT9W5Ko7lljJ5kj2aKXcA0NJu5OWtR5kW5smMCC8bR2h798yJQCsEL24uNG8QghtX3IfeScMX66xaxVQZ5y6a9IUQWuBFYAkQD9wqhIjvMew4cBfwbo9jvYCngOlAOvCUEMLz8sNWxqLXv3oNY4ORm6dEgHc0AB/llnOitpkHrooe0312rOXv7sDK9GA+zC3rvNrXpq4kaZIDJYYyTp87fZEzKOOdNVf66UChlLJIStkKrAaWdx0gpSyRUuYDph7HXgN8I6WsllKeBb4BFqMoffjo4zdBwF23/xIwt1x4aUshKcEeXBntbePoRo6fzItEIwQvbTlq3uDszS1XJNNe285rX/zTtsEpI541ST8QKO3yusyyzRpWHSuEuEcIYRBCGCorK608tTKW1LXWsS+7mKAQB3xn3Q7Ae9nHKTvbxM8Wqqv8rgLcHbllWjAf5pZyrOocAHf+8LcAfPixasmgXNiIuJErpXxFSpkmpUzz8Rk/j9cr56359lXOHW1kyRUpoNVT39zGCxuPMCPCi3ljcDnEy3X/VVHoNBr+st48Y8cv9RqCQ5w4mFvBmUZ14aT0z5qkXw4Ed3kdZNlmjcs5VhlHXv/gf0HC93/03wD877Yiqs618uiSOHWV3wc/NwfumRPBuvwT5B4zt12+/qrZNBY2smbr32wcnTKSWZP0c4BoIUS4EMIOWAlYO01gPbBICOFpuYG7yLJNUTpVVh9hd3YFbh72pM9bwqm6Zv53ezFLkwOYHDx2V8a6XD+eG4Gvqz2r1u1HSskd338MJLzxqXo6V+nfRZO+lLIduB9zsj4AvC+lLBBCPC2EWAYghJgmhCgDvgP8SwhRYDm2Gvg95i+OHOBpyzZF6fTZpt/TsLeB65YsRKPRsGrdAYxS8utrJtk6tBHNyU7HLxfFsOt4Dev2nmD6rCtx93QgP/cMx4s22To8ZYSyqqYvpfxCShkjpYyUUq6ybHtSSrnW8nuOlDJISukspZwgpUzocuzrUsooy8+/h+ZjKKNWWxNvbN2CqcXE9+58gG+PnmHtngrunRtJ6ARnW0c34t08NZj4ADd+//l+zrUaWb58BfX59Xy47Y+2Dk0ZoUbEjVxl/Dqa9f/Ym9+Ao7M9s66cy5OfFhDs5ch988b307fW0moE/3NjEqfrW/jb14e547bvI1slbxkOYKo6auvwlBFIJX3FdoxtfLj3Tep317F4ybW8nlFK4ekGfnt9Ag56ra2jGzVSgj343oxQ3sgowSMyBWc3J0ryzpG77fe2Dk0ZgVTSV2ymNfc/vHu8gfZ6IzOvup4XNh5h2eSJLIjzs3Voo87D10zC19WeJ9Ye4PrrV1CfV88npduhtszWoSkjjEr6im20NbM5+znKdzWgt9PzVa0/Hk52/G5ZwsWPVXpxddDzhxVJHDxZjz5yJsZzRj463k79llW2Dk0ZYVTSV2wj9998JFqoz64jZuocCmtM/PHGpHG5QMpguTrej5XTgtl6LgB7BwfO7Gpg7dHPoPKwrUNTRhCV9JXh11JP6c5n2XDcRGtNG6d907h5ahBXx6uyzuV64rp4wv08cYmZQYOhgfecXJGbVG1fOU8lfWX4bf8b7+lbqc2qQ2PnQOKM+Ty9XJV1BoOzvY7nV6aij5lLa30rew82klO8HsoMtg5NGSFU0leGV3UR5zJf4iNnd+pyGnCJuYJ/3T0TJzudrSMbM1KCPXj6/tvQOLpSn3mO1V7e8MWvwNSzCa4yHqmkrwyvr59gjYsLJ/fW0t7Yyv0/upNov/G7BOJQ+eGcaBJnL6Ymt4avTYKK03tg91u2DksZAVTSV4ZP4QZMBz/n356+VG9pxtljAr/9yW0XP065ZEIInn/8AWS7kRpDHa8FRMGG30Kj6oIy3qmkrwyPlgb47Od8OSGMspoazu09w09+eDd6vd7WkY1Z8+bMJiw8gpqtbXygaaO6tQ42/s7WYSk2ppK+Mjw2/R5TbSm/tfOkbmcb0mTiBz/4ga2jGtOEENz743s4V3iK5pPneMJjCuT+B45utnVoig2ppK8MveNZyKx/8XuHK2iyP4PMaWLWrFnExsbaOrIx76677kKn06HNcmSbcy2HtUGYPr0fmutsHZpiIyrpK0OruY6WD35IBRNYO0Fgf9yRkyUn1FX+MPHz82PFihWc3H4ck+kc9zqlI+sqaPniMVuHptiISvrK0JGSU+/9FF1dKb9y/w6tdqXYZejw8PDgu9/9rq2jGzfuueceaqprCCwKpMHvAC9yDfb5b1FreN/WoSk2YFXSF0IsFkIcEkIUCiEe6WO/vRDi/yz7s4QQYZbtYUKIJiFEnuXn5cENXxnJ8j7/J37H1vK24600hB1kQvMEcr7J4Uc/+hHOzqpX/nBZsGAB0dHRnFp/iqb2cxybmUiejEb3+YMcO7LX1uEpw+yiSV8IoQVeBJYA8cCtQoj4HsN+AJyVUkYBzwF/7rLvqJQyxfJz7yDFrYxgUkre+/QzJhmeokCfhP76+RTXHcU3zxcpJT/96U9tHeK4otFoeOihh8jflU/yuWS2nf6Uxpuepw0tTe/cQc6RCluHqAwja67004FCKWWRlLIVWA0s7zFmOfCG5fcPgQVCrWY9LjW3GXn87U3M2fUzmvXuBN7zFq/t/xdxLnFsfH8jN9xwA6GhobYOc9y588478fLyovqraozSyLq6L2m97p/EUsKpt77PGzuLkFLaOkxlGFiT9AOB0i6vyyzb+hxjWVO3Fphg2RcuhNgthNgqhLjyMuNVRrAjp+r5zj+2cMORR/DVnsPj+x/wZulnnG48TfD+YKqrq/n5z39u6zDHJWdnZ+69917Wr1vPtW7XsvboWspDgmma+yTXaTKo+fJpfv5/eTS2tts6VGWIDfWN3BNAiJQyFfgF8K4Qwq3nICHEPUIIgxDCUFlZOcQhKYNNSsl72cdZ8Y+tPFT7Z6ZpDqO/8SWKndz5d8G/WRK4hNUvr2bevHnMmjXL1uGOW/fffz96vZ7yteX4OfmxKnMV+jkPIlPu4Ge6T9DuXc2yf+wkv6zG1qEqQ8iapF8OBHd5HWTZ1ucYIYQOcAeqpJQtUsoqACllLnAUiOn5BlLKV6SUaVLKNB8fn0v/FIrNnKpr5sdv5fLox/m86PoWC8iCa/6ITLiR32f+HrPt1/UAABuNSURBVCedEz77fKioqOCJJ56wdbjjWkBAAD/+8Y95+823ucPvDg6dPcR7h/5/e2ceX0WRLuynzn5OkpOdLCSQBBN2kUV2BVz5cAVlhFHRufONG4vO1auCy+edmd8ojus416swznXg0xl3QR3CpiibSBDQsBogIQlk306Ws9f9ozshQUICoicx9fBrqrq6+vTbneq33nqruuqfiKufh/RJPGNZyoWNXzD95S08v/YgvoCaoO3nSGeU/nYgUwiRLoSwALOAlSflWQncpsdvBD6VUkohRLzeEYwQIgPIBA6fG9EVoSQYlLy57SiXPfs5nx8s46PMT5jcmA0XPwjj7uG9794jpzSHeUPn8dJzLzF+/HimTJkSarF7PA8//DAmk4nP/+dzLk65mBe/fpEjDcUw601Eyij+GHyBhRlHeHH9d1z3l83sKKgOtciKc0yHSl/30c8DVgP7gLellHuEEL8TQlyrZ3sNiBVC5KG5cZqHdV4MfCOE2IXWwXuXlFLN+NTN2VVYw8xXt7Log28ZmhxOzrBPGFr4Joy9B6YsoqCugKe3P82YpDEcX3OcwsJCnnjiCVTffuhJTk7mzjvvZPny5dwafytWo5VHNz9KwGyHm99BJA7l/x5/gg+mVFDV4OWG/97Cg+/uprLeE2rRFecI0dV67EeNGiVzctSCD12RYzVNPJ29nw93HSMu3MrDl2dwQ9GTiG/fhovuh0sewyf93LbqNvLr8nlt4mtMGDaB0aNHk52dHWrxFTolJSVkZmZy6aWXcsdzd/DQxoeYd8E87hx2pzYL55s3QXEOniuf4bnq8by28QgOi5G7JvfjV+PTsVuMob4FxSkQQuyQUo7qKJ/6IlfRIcdqmnh8RS6Tn9nAqtwS5k05jw1zh3Lj3nmawr/0cW0Tgpd2vsS3Fd/y+LjHWfLcEmpra3nmmWdCfQuKViQmJvLII4+wYsUKzIfMTEufxsu7X2brsa3giIE5H0K/S7Fm/zsLbR+Qfe8ERqXF8HT2AS7+02cs25qP16/8/d0VZekr2qWwqpH//vwQ7+QUIiXMHJXCvEsy6e0tgH/cBHXH4br/gvNnApB9JJv/+OI/+EXWL5jhnMHw4cO5/fbbWbJkSYjvRHEybrebwYMHY7fb2bJ9C3NWz6HKXcXb17xNYlgiBHzw0b2w6w3ofxVMf4WcEj9Prz7AV0eqSHBa+bcJ6cwe0wenTU2P3RXorKWvlL6iDVJKNudV8vqWI6zfX4bZYOAXF6Zw16R+pETZtdWXVj0ElnCY/Q9I0crYgaoD3LrqVgbEDGDJpUuYMmkKeXl57Nu3j7i4uBDfleJUfPjhh0yfPp2nnnqKmXfNZPbHs+nr7MvrU1/HYXaAlLDtVVi9CGIy4Bd/R/YaxKa8Cl79/DCb8ioIt5r45Zg+3DymD31j1dQaoUQpfcUZUVHvYeWuY7z51VHyyuqJDbMwe3Qfbhnbl8RIGzTVwMe/hT3vQ/rFMH0JOJMAKHQVMmfVHAzCwFtXv8U/X/sn8+fPZ/ny5dxyyy0hvjPF6bjxxhv56KOP2LlzJxXOChZ8uoCxyWN56ZKXMBt0Cz5/E7zzK3DXaG68sXPBYCC3uJZXvzjMv749TiAomXheHL8c04fLBiZgMSnP8U+NUvqKDnH7AqzfV8b7Xxex4WA5gaBkWEokc8alcdX5SdjMRs3a27sCVj0IDRUwZRFM/C0YtM68iqYK5qyaQ523jr9P/Tv+Uj8jR45k4sSJrFq1So3Y6eKUlpYyePBg+vXrx+bNm1lxeAVPbH2CazKu4fcTfo9R/zvTUKG5e/Z/DH0nwjUvQFwmACW1bt7OKeSt7YUU1zQRG2bhqvOTuGZYMiP7RGMwqDLwU6CUvuKUuNw+PjtQzpo9JWw4UE69x0+i08b1w3szY0RvslovUl51GLIXwcFVkHg+XPtnSB7ecri0oZTfrP0NJQ0lLL1iKf0j+jNmzBiKi4vZvXs3vXufPFuHoivy1ltvMWvWLB566CGeeuopXt39Kn/Z9Reuzria30/4PSaDScsopebjz14EvkYYPx8ufgAsmlsnEJR88V057+QUsn5fGR5/kKRIG1efn8SVgxMZ3icao6oAfjQ6q/RNP4UwitAhpeRQeQObvitnw8FytuRV4g0EiQu3cM2wJKYNTWJ8v7i2L2NDJXzxJ9j+VzCa4Yo/wJi7wXiiuBS6CvnNmt9Q46nh5UtfZlj8MO6++26++eYbPvnkE6XwuxE33XQTn332GYsXL2b8+PHcee2dCH0kli/o448T/4jFaAEhYPgtkHkFrP1/sOk5+OYtmPQQXPBLjEYzU/r3Ykr/XtR7/KzbW8pHu4/x+pZ8lm48QpTDzKSseC4Z0ItJWfFEOSyhvvUeibL0f4aU1rn58nAlG7+rYHNeBcdr3QCkxTq4fFBC+1ZXQwVse0XrvPPWw/BbYfLCFt99M1+Xfs1vN/yWgAzwymWvMCRuCC+99BILFizgwQcfZPHixSi6F263m4kTJ5KXl8emTZsYMmQIr+e+zrM7nmV4r+G8MOUFYmwxbU8q2AprHoXiHK2jd/JCGDyjjXEAUOf2sfFgBev3l7LhQDlVDV4MAob0jmRsRixjM2K4MC2GCDUK6Aeh3Ds9BF8gyP7jLnYUVLHjaA1fF1RTXNMEQKTdzITzYpl4XjwXZcaRGuM49Y9UHYEvX4avl4PfDQOugkseg15t17CVUvLOwXd48qsn6R3emz9P+TMZURmsXLmS6dOnc8011/Dee+9hNKqPd7ojR48eZdy4cQgh2Lp1K6mpqWTnZ/PopkeJs8fx7ORnGRw7uO1JUsLBbPj0D1CaC5GpMOZOGDEHbJHfu0YgKNldVMOGA+V8ebiSXUdr8AaCGAQM7R3JiL7RXJAaxbCUKPrGOlSf0BmglP7PkCZvgAOlLvYeq2Pv8Vr2Hqtj33EXTb4AAIlOGyP7RjOibzQXpkUzODmyfR+q3wP7P4Gv/w6HN4DBDMNugvH3Qvz35sSj2l3N77b+jnVH1zEheQKLL15MpDWSVatWMX36dIYOHcqGDRvUiljdnG+++YaLLrqI5ORk1q9fT3JyMrkVudz32X1UuiuZP3w+tw++HYM4aXROMKj1/Wx9GQo2aUN6h86EC27WhvW2o7ybvAF2Hq3my8OVfHmkim+LalvKc6TdzLDUKIalRDIwyUn/xAjSYsNUv0A7KKXfjWnw+DlS0cCh8noOl2vh/hIXh8vrCep/rgibiYFJTgYnOxnRJ5qRfaNJjrKf/ocDfijYDHs/1EbkNFZCZB8YcavmyjnJjQOadb+mYA1PbnuSWm8t84fP57ZBt2E0GFm5ciUzZ85kyJAhrF27lpiYmFNcVNHd2LhxI9OmTSMhIYF169aRlpZGraeW/9z6n6wtWMuIXiN4ZOwjZEV/3zgA4NguzUW45wPwN0FsJlwwGwZeB3Hnnfba/kCQ78rq2V1Yw+6iGnYV1nKgpK6l3FtNBjITwumf4GRAYgRZiRGkx4bRO9re4ysDpfS7MFJKapt8FFU36VsjBZWNHK7QlHyzDx40Aykl2s6ARCeDkpwMStbClGh755q+7lo4shG+W6NZ9o0VYHZA1pVap1zGlJbhlydzoOoAi7cvZnvJdgbEDOAPE/5A/5j+SCl5/vnneeCBBxg1ahSrV68mOjr6XD0eRRdg27ZtTJ06FZvNxvvvv8+4ceOQUrLi0AqezXkWl9fF7AGzueP8O4i2tfO3d9dpxsWuN+HoFi0trj8MvBr6T9NGgrVT9tr8jC/Ad6X17C+p40CJiwOlLvYdd1HRahI4s1GQGuMgIy6MtNgw0uLCSI8Lo3eUnaQoG1bTz9/lqJR+iJBS4vL4KavzUOZyU+7yUFrnprhFwTdRXNNEvaftCkURVhMZ8WH0iw8nIz6MjPhw+sWH0zfWoY2X7yzeRji2U/ug5tCnULQdZEBrbmdeAYOvh/MuB0s7/n1gX+U+ln67lLUFa4m0RrJg+AJuyLwBo8FIbW0tc+fO5Y033mDGjBksW7ZMuXR+puzZs4drr72WwsJCnn/+ee655x6EENR6annx6xd59+C72E12bh54M3MGzSHKFtX+j9UUwoF/wb6PoGCLViatkZA2ETImaR/8xfUHQ+c/6qqs95BXVk9+ZQOHKxrIr2ggv6KR/MoGPCfNDRQfYSU5yk7vKBvJkXYtHm0n0Wmjl9NKbJi1239QppT+OSQQlNQ0eqlu9FHd6KWqwUtNo5fKBi9ldR7KXZqCL3N5KKvztPgkWxNhNdE72k5KtIOUaHvL1jtK249ymM+80yrgh6pDcHw3FH4FRV9BSa72QiEg+QLod4m2pYwGU/tD5HwBH58Vfsa7B99l6/GthJvD+eXAXzJn0BwirVqH3OrVq7njjjsoLi7m8ccf59FHH8VwBi+povtRVVXFzTffTHZ2NpdddhlLly4lLS0NgEM1h3hl9yuszl+NxWhhWvo0bhpw0/c7e0+msUozSI58Doc/h5oCLd0WCckjoPdIfRsBEYlnLHMwKCmpc5Nf2cCxGjfHapoorm7iWK1mcB2racLt+/6EcdEOM/ERVm0L18JeETbiI6xEh1mIdpiJdliIcpgJt5q6XCezUvon4fUHcbl9uNx+ffNRp4f1nhNp1Y0+qhu8VDV6qWn0UdXgpc7to73HFG410UsvKL2cNhIirPRyaoWlV3Pcafthk1IFfFBzVBtlU7ZX20r3QPkBCOhNXHOY9pKkjtYUfOpobcbE0+AP+tlZtpN1BevIzs+myl1FUlgSN2bdyKwBs3BatJUt9+7dy6JFi1ixYgWZmZksW7aMsWPHnv39KLoVUkpeffVVHnjgAfx+PwsWLGDhwoUtLr286jze2P8Gnxz+hCZ/E1nRWUxNm8qVaVfSx9mn4wtU52st06Icbfhn6V7dcAEcsRA/EHoN1EaT9Rqk9RGExbXbOdyZ+6lu9FFc3URpnZvyes1YK6/XWublLk9L2skthmbMRkGkvW1FEO2wEBWmhRE2E+FWE06bmXCbiQibiQibVlmEW00/Sv9Dt1b6G7dso9Hrp9EboMkXoNEboNHrp8l7It6ox5vTmnxt05oVfLNib++P1xqb2UC0w6Jt+h9Pi1uIcZj12v7E8ZgwCw7LOfi+zecG13FwlYDrGNQWaQq++ogW1hadeAkAIpK0wp8wCHoNhsQh2othPL0sUkqK6ovIKclhe8l2NhVvotpTjdVoZWLviczInMGE5AkYDUaklGzatIlnn32WFStWEBYWxmOPPcZ9992H1Wr94fes6HYUFhby2GOPsWzZMpxOJ7/+9a+ZO3cuGRkZALi8Lj469BGrjqxiV/kuAPpF9mNs8ljGJo1lVMIowi3hHV/I2wgl32huyrK9ULYfyvaB13UijzkMotP0ra8WRvXRWgYRSRAW36n+gtPR7Kotd3mobjjR0m9u9dc0eqluaE47EXo7scxks/KPsJn0SsFMuNWI3WzCYTHisBixt4QmHObWaaY2xx1mE3aLEavZeO6UvhBiKvAiYAT+KqV86qTjVmAZMBKoBG6SUubrxxYCvwYCwAIp5erTXcualCmTbnuhQ5lOXBscZu3BhFmN2PWHE24zE2Ez4dQfaIS1VW2r17xOPU9zDXzOfHoBn9aEbapqGzZW6vHqtkq+6RRL0tljICYdotPbhvEDOrTgQSuwJQ0l7K/az/7q/RyoOkBuRS6ljaUAxNhiGJM0hsv6XMbE3hNxmB0Eg0F27tzJxx9/zPLlyzl06BDR0dHMnz+f+fPnq9kyFQDs3r2bJ598kvfee49AIMCkSZOYOXMm119/PcnJyQCUNJSwJn8Nm49tZkfpDjwBDwJBRmQGg+MGMyh2EANiBpDmTCPGFtOxq0RKqCvWlH/lIa11UJ2vuYaq87VpIVojDBCeoFUC4YkQkQCOOO3dscecFEaDLeqM+hPaF1PS5Au08Si43H7dm9Da0/B9L0O9x4/bF2wxajtjqLamYPHV50bp62vcHgQuB4rQ1sydLaXc2yrPPcD5Usq7hBCzgOlSypuEEIOAfwCjgWRgHZAlpfy+01snNWuIXLRkRduazHKqGlDbt5oMZ+dbk1L7EMnv1ixtv1sbu+5v0kJf04l9bwN46jVLw1Ovfa3aErpOhJ7mfVf71zXZtMLWbJE4k/R4shY6k7V0e/udYlJKGnwNVLorqWyqpNJdSVljGUWuIorqi7TQVYQ7oI0CEgj6OvsyIGYAIxNGcmHihWREZlBdXc3OnTvZuXMnOTk5rF+/noqKCoQQTJkyhTlz5nDDDTcQHt4J60zR4yguLmbp0qW89dZb7N+/H4CsrCwmT57M2LFjGTp0KIMGDcJoNbK7bDc7SneQW5lLbkUuVe4Tq6ZGWCJId6bT19mXxLBEEhwJ9HL0IiFMC6OsUSfm/zkVUkJDueYCdZWcMKjqS/R9Pa2pGmQ7ilQYtD4FSwRYI8AaroWWcD3u1OP6MbNDe5fN9o7Ds2xxBIJS93S09nI0ezf8rbwgAZq8fuZfmnXOlP444Akp5ZX6/kIAKeWTrfKs1vNsFUKYgBIgHn2t3Oa8rfO1d71RA/vInKX/DkGf1lEZ9EHQfyIe8LU65m+bL6DnDfpPyufTFXgrJR84yzU/DWawRiCt4QQs4VpoDiNoDUNa9DR7FAGbk6A9CmmPwm+JwGsNw2t24DUY8Aa9eAP6FvTiC/ha0jwBD42+Rlw+Fw3eBup99TT4ToQur4vKpsoWhS6lhAAEfUGsQSsJ1gTiTfHEmeJwBpw4A07sXjuuGhfl5eUUFBSQn5/PkSNHqKysbLmt1NRUJk+ezBVXXMHll19OQkLC2T0fRY9DSsmePXtYvXo1GzZs4IsvvqCurg4AIQTp6en06dOHlJQUUlNTSUhIwBBmoMHUQL2xnhpRQ1mgjBJvCVX+KqRJIkwC0crvHW4OJ9IaidPiJNIaSaQ1knBzOHaTHbvJjs1k+17cZrRhMpgwG8xaKIyY/G5MngbM3gZMbhdmjwuTuw6zuwaDuw6DtwHhqcfgdWHw1iM89a2MuXrgLNzhBrOm/I0WfTPrmwUMprZpBj3dqKcbWuU1mvX8ZhBGLW4waZWKwYSYeO85m3CtN1DYar8IGNNeHimlXwhRC8Tq6V+edO5pZ+Lak1fIoOn3t3m0EgEIECClACH04+LEcdEq3iqEVnmFaPktiUE/RzsuEZT5TljoWl0oW8Wb0xuB2rZCy+/HZevEU5WTU5xzMkK/h9Zhyz8hCAaC+Dw+gsHONQNtNht9+vQhPT2dkSNHkpGRwfDhwxk+fLhy3SjOGiEEQ4YMYciQIdx///0EAgEOHTpEbm4uubm57Nu3j6KiIjZu3EhxcTF+v7/jHwVMZhMms0lT/gZoVgMSiRRSi4sTcSFES55zfX8t8db/n3QdAViFmQRzuK40glooJdAIsqFFO2lpUn/35UlpbY+J7x1rPufs6BKzbAoh7gDuAHA47AyZNBWEASEMLU+2zYM/TdrZHg/IALkVuS0KtWXT9xFg0JcUNhgMbfIBGIShTbx1mhACozBiNBgxYMBoNLbsG4URgzBgMpgwCANGYcRsMmMUxg7vxWg0Yrfbsdlsp9yioqKIjY0lJiaG2NhYHA41l4nix8doNJKVlUVWVhYzZsxocywQCFBVVUVNTQ01NTVUV1dTXV1NY2Mjbrcbj8fTEjbHg8Fgh5s/4Mfn9+H1ewnIAIFggIAMIJEEZZCgDLbEpZRt9pvTZIuRp8elXsG0Sgc9TZ6UV093mBykR6af8rmczaCZDs+RUndZSfa+/2GnfrMzSr8YSG21n6KnnSpPke7eiUTr0O3MuUgplwBLQBu98/a773dKeIVC0b0wGo3Ex8cTHx8falF+dnTWoOtMd/V2IFMIkS6EsACzgJUn5VkJ3KbHbwQ+lVoVtRKYJYSwCiHSgUzgq05JplAoFIpzToeWvu6jnwesRhuy+Tcp5R4hxO+AHCnlSuA1YLkQIg+oQqsY0PO9DewF/MDc043cUSgUCsWPS5f8OKurTcOgUCgUXZ3OfpGrJk5RKBSKHoRS+gqFQtGDUEpfoVAoehBK6SsUCkUPQil9hUKh6EF0udE7QohyoOBHvEQcUPEj/v6PjZI/tCj5Q0t3lv/Hlr2vlLLDr966nNL/sRFC5HRmWFNXRckfWpT8oaU7y99VZFfuHYVCoehBKKWvUCgUPYieqPSXhFqAH4iSP7Qo+UNLd5a/S8je43z6CoVC0ZPpiZa+QqFQ9Fh6rNIXQswXQuwXQuwRQjwdannOBiHE/UIIKYToVsteCSH+pD/7b4QQHwgh2l8QuIsghJgqhDgghMgTQjwcannOBCFEqhDiMyHEXr283xtqmc4GIYRRCLFTCPFxqGU5U4QQUUKId/Vyv09fhjYk9EilL4SYAlwHDJNSDgaeCbFIZ4wQIhW4AjgaalnOgrXAECnl+cBBYGGI5TktQggj8F/A/wEGAbOFEINCK9UZ4Qful1IOAsYCc7uZ/M3cC+wLtRBnyYtAtpRyADCMEN5Hj1T6wN3AU1JKD4CUsizE8pwNzwMP8kMWywwRUso1UsrmhVK/RFtRrSszGsiTUh6WUnqBf6IZDd0CKeVxKeXXetyFpnBOu1Z1V0MIkQJcBfw11LKcKUKISOBitHVHkFJ6pZQ1oZKnpyr9LOAiIcQ2IcTnQogLQy3QmSCEuA4ollLuDrUs54B/A1aFWogO6A0UttovopspzWaEEGnAcGBbaCU5Y15AM3KCoRbkLEgHyoH/0d1TfxVChIVKmC6xMPqPgRBiHZB4ikOPoN13DFpT90LgbSFEhuxCQ5k6kH8Rmmuny3I6+aWUK/Q8j6C5Ht74KWXrqQghwoH3gPuklHWhlqezCCGuBsqklDuEEJNDLc9ZYAJGAPOllNuEEC8CDwOPhUqYnyVSysvaOyaEuBt4X1fyXwkhgmjzYpT/VPJ1RHvyCyGGolkOu/WFkFOAr4UQo6WUJT+hiKfldM8fQAhxO3A1cGlXqmzboRhIbbWfoqd1G4QQZjSF/4aU8v1Qy3OGTACuFUJMA2yAUwjx/6WUt4RYrs5SBBRJKZtbV++iKf2Q0FPdOx8CUwCEEFmAhW4yiZOU8lspZS8pZZqUMg2tQI3oSgq/I4QQU9Ga6tdKKRtDLU8n2A5kCiHShRAWtDWgV4ZYpk4jNOvgNWCflPK5UMtzpkgpF0opU/TyPgv4tBspfPR3s1AI0V9PuhRt3fCQ8LO19Dvgb8DfhBC5gBe4rRtYmz8n/gJYgbV6a+VLKeVdoRWpfaSUfiHEPGA1YAT+JqXcE2KxzoQJwK3At0KIXXraIinlv0IoU09jPvCGbjQcBn4VKkHUF7kKhULRg+ip7h2FQqHokSilr1AoFD0IpfQVCoWiB6GUvkKhUPQglNJXKBSKHoRS+gqFQtGDUEpfoVAoehBK6SsUCkUP4n8BFKaKlBRYqiUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from scipy.stats import t\n",
    "\n",
    "df_list=[1,3,10] # list of the degree of freedom\n",
    "\n",
    "for j in range(len(df_list)):\n",
    "    df=df_list[j]\n",
    "    x = np.linspace(t.ppf(0.0001, df), t.ppf(0.9999, df), 500000)\n",
    "    plt.plot(x, t.pdf(x, df), label='df=%d'%df)\n",
    "    \n",
    "x = np.linspace(-10,10,1001)\n",
    "plt.plot(x,standard_norm_pdf(x), 'k', label='Standard normal dist.') # standard normal distribution\n",
    "\n",
    "plt.legend(frameon = False)\n",
    "plt.xlim([-7,7])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Useful commands"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Percent point function - Find the point where you can get specific probability."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can use <strong>percent point function (ppf)</strong>. If you want to find the point which gives the probability of 0.1 with the specific degree of freedom, df, you can type:\n",
    "`t.ppf(0.1, df)`\n",
    "\n",
    "See also [scipy.stats.t](https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.stats.t.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xt8zvX/x/HHi23OJUVkvsysMZRYRCdFmg6sr1T6OaQ00+T0VRERfn0rqW+FL6YphS8r5TuHCEUHxzlEaMyhjBUiItuyvX5/7OK31tg1ru1zbdfrfrtdN9f1ud6fa88Vz332eX8OoqoYY4zxDaWcDmCMMaboWOkbY4wPsdI3xhgfYqVvjDE+xErfGGN8iJW+Mcb4ELdKX0QiRCRJRJJFZMgFxnUSERWRcNfrOiJyWkQ2ux6TPRXcGGNMwfnlN0BESgMTgbuAFGC9iCSo6vZc4yoB/YG1uT5it6o28VBeY4wxl8CdLf3mQLKq7lHVDGA20DGPcWOAV4E0D+YzxhjjQflu6QM1gf05XqcALXIOEJGmQC1VXSgiz+RaP0hENgEngOGq+lXuLyAiUUAUQIUKFZrVr1+/AN+CKRRJSdl/hoY6m8MY45YNGzYcUdWq+Y1zp/QvSERKAW8Aj+XxdirwN1X9RUSaAfNEpKGqnsg5SFVjgViA8PBwTUxMvNRY5lK1bp3954oVTqYwxrhJRH5wZ5w7u3cOALVyvA50LTurEtAIWCEi+4CbgAQRCVfVdFX9BUBVNwC7gWvdCWaMMcbz3Cn99UCIiASJSADwCJBw9k1VPa6qV6lqHVWtA6wBOqhqoohUdU0EIyJ1gRBgj8e/C2OMMW7Jd/eOqp4Rkb7AEqA0ME1Vt4nIaCBRVRMusPptwGgR+QPIAqJV9agnghtjjCk4t/bpq+oiYFGuZSPOM7Z1judzgbmXkM8YY4wH2Rm5xhjjQ6z0jTHGh1jpG2OMD7HSN8YYH2Klb4wxPsRK3xhjfIiVvjHG+BArfWOM8SFW+sYY40Os9I0xxodY6RtjjA+x0jfGGB9ipW+MMT7ESt8YY3yIlb4xxvgQt0pfRCJEJElEkkVkyAXGdRIRFZHwHMuGutZLEpG7PRHaGGPMxcn3Jiqu2x1OBO4CUoD1IpKgqttzjasE9AfW5lgWRvbtFRsC1wDLRORaVc303LdgjDHGXe5s6TcHklV1j6pmALOBjnmMGwO8CqTlWNYRmO26QfpeINn1ecYYYxzgTunXBPbneJ3iWnaOiDQFaqnqwoKua4wxpuhc8kSuiJQC3gD+cQmfESUiiSKSePjw4UuNZIwx5jzcKf0DQK0crwNdy86qBDQCVojIPuAmIME1mZvfugCoaqyqhqtqeNWqVQv2HRhjjHGbO6W/HggRkSARCSB7Yjbh7JuqelxVr1LVOqpaB1gDdFDVRNe4R0SkjIgEASHAOo9/F8YYY9yS79E7qnpGRPoCS4DSwDRV3SYio4FEVU24wLrbRCQe2A6cAWLsyB1jjHFOvqUPoKqLgEW5lo04z9jWuV6/BLx0kfmMMcZ4kJ2Ra4wxPsRK3xhjfIiVvjHG+BArfWOM8SFW+sYY40Os9I0xxodY6RtjjA+x0jfGGB9ipW+MMT7ErTNyjfFVqsqmTZtYtWoVBw4cwM/Pj3r16tGmTRsCAwOdjmdMgVnpG5OHrKws3n//fV555RWSkpIA8Pf3JzMzk6ysLAAiIiIYNWoUzZvbfYFM8WG7d4zJZd++fdx666307NmTChUqEBcXR0pKCunp6aSnp7NlyxZGjRrFxo0badGiBU8//TRpaWn5f7AxXsC29I3JYeXKlURGRpKVlcV7771H9+7dEZFz7/v5+dG4cWMaN27MwIEDGT58OG+//TZr165lwYIFVKtWzcH0xuTPtvSNcVm6dCnt2rWjRo0abNq0iR49evyp8HOrVKkSb731Fp988gnfffcdt9xyC4cOHSrCxMYUnJW+McC6det44IEHqF+/Pl9//TV169Z1e93IyEiWLl1KSkoKERERnDhxohCTGnNprPSNz9u5cyf33HMPV199NYsXL6ZKlSoF/oybb76ZuXPnsnXrVjp27EhGRkYhJDXm0rlV+iISISJJIpIsIkPyeD9aRLaKyGYR+VpEwlzL64jIadfyzSIy2dPfgDGX4vTp03Tu3BkR4bPPPqNGjRoX/Vnt27fn3XffZcWKFQwbNsyDKY3xnHwnckWkNDARuAtIAdaLSIKqbs8xbJaqTnaN7wC8AUS43tutqk08G9sYzxg4cCBbtmzh008/JTg4+JI/r2vXrqxZs4Zx48bRunVr7r33Xg+kNMZz3NnSbw4kq+oeVc0AZgMdcw5Q1Zw7MSsA6rmIxhSO+Ph4pkyZwnPPPUdERET+K7hp3LhxNGnShO7du5OSkuKxzzXGE9wp/ZrA/hyvU1zL/kREYkRkNzAW6JfjrSAR2SQiK0Xk1ry+gIhEiUiiiCQePny4APGNuThHjhzhqaeeonnz5owZM8ajn122bFnmzJlDWloaMTExqNo2kPEeHpvIVdWJqhoMPAcMdy1OBf6mqjcAg4BZInJZHuvGqmq4qoZXrVrVU5GMOa/Bgwdz/Phx4uLi8Pf39/jnX3vttYwaNYqEhAQ++eQTj3++MRfLndI/ANTK8TrQtex8ZgORAKqarqq/uJ5vAHYD115cVGM8Y/ny5UyfPp1nn32WRo0aFdrXGTBgAE2aNKFv374cP3680L6OMQXhTumvB0JEJEhEAoBHgIScA0QkJMfLe4FdruVVXRPBiEhdIATY44ngxlyMjIwM+vTpQ7169Rg+fHj+K1wCPz8/YmNj+fnnn3nhhRcK9WsZ4658S19VzwB9gSXADiBeVbeJyGjXkToAfUVkm4hsJns3Tg/X8tuALa7lHwHRqnrU49+FMW6aNGkSu3bt4q233qJcuXKF/vVuvPFGoqKimDRpEjt37iz0r2dMfsTbJpnCw8M1MTHR6RimdevsP1escDKFRx07dox69erRtGlTPvvsswteYsGTDh06dO5yzLZ/3xQWEdmgquH5jbMzco3PeOmllzh27Bjjxo0rssIHqFatGkOGDGHevHmsXLmyyL6uMXmx0jc+4YcffmD8+PE89thjXH/99UX+9QcMGEBgYCDPPPOMHcJpHGWlb3zCSy+9BMCoUaMc+frly5fnxRdfZP369SxcuNCRDMaAlb7xAXv37uXdd9/lySefpFatWvmvUEi6d+9O3bp1efHFF21r3zjGSt+UeC+99BKlS5dm6NChjubw9/dn+PDhbNiwgfnz5zuaxfguK31Tou3Zs4f33nuPqKgoatb8y9VDily3bt0IDg62rX3jGCt9U6K99tprlC5dmiFD/nJFcEf4+fkxfPhwNm3axOLFi52OY3yQlb4psQ4dOnTuPrfXXHON03HOefTRR6lZsyavvfaa01GMD7LSNyXWhAkTSE9PZ/DgwU5H+ZOAgAAGDhzIF198wfr1652OY3yMlb4pkU6dOsXEiRPp0KEDoaGhTsf5iyeffJLLL7/ctvZNkbPSNyXStGnTOHr0KM8++6zTUfJ02WWX0adPH+bOncvu3budjmN8iJW+KXHOnDnD66+/TqtWrWjVqpXTcc6rX79++Pn58cYbbzgdxfgQK31T4nz44Yf88MMPXruVf1aNGjXo1q0b06ZNw+4YZ4qKlb4pUVSV119/ndDQUO6//36n4+Rr8ODBpKWlMWnSJKejGB9hpW9KlDVr1rBhwwb69etHqVLe/9e7fv36REREMHnyZDIyMpyOY3yAW/8qRCRCRJJEJFlE/nKWi4hEi8hWEdksIl+LSFiO94a61ksSkbs9Gd6Y3MaPH89ll11G9+7dnY7itqeffprU1FQ+/vhjp6MYH5Bv6btudzgRaA+EAV1ylrrLLFVtrKpNgLHAG651w8i+vWJDIAL499nbJxrjaampqXz44Yf07NmTihUrOh3HbREREQQHBzN+/Hinoxgf4M6WfnMgWVX3qGoG2Tc+75hzgKqeyPGyAnD2oiIdgdmuG6TvBZJdn2eMx8XGxnLmzBliYmKcjlIgpUqVIiYmhlWrVrFx40an45gSzp3Srwnsz/E6xbXsT0QkRkR2k72l36+A60aJSKKIJNpRDOZiZGRkMHnyZNq3b09ISIjTcQqsZ8+elC9fngkTJjgdxZRwHpvpUtWJqhoMPAcML+C6saoarqrhVatW9VQk40Pmzp3LTz/9xNNPP+10lItSuXJlunfvzqxZszhy5IjTcUwJ5k7pHwBy3nki0LXsfGYDkRe5rjEXZfz48dSrV4+77y6+xwr07duX9PR03nnnHaejmBLMndJfD4SISJCIBJA9MZuQc4CI5Px9+l5gl+t5AvCIiJQRkSAgBFh36bGN+X+bN29m9erVxMTEFIvDNM+nYcOG3HnnnUyaNInMzEyn45gSKt9/Iap6BugLLAF2APGquk1ERotIB9ewviKyTUQ2A4OAHq51twHxwHZgMRCjqva32XhUbGwsZcuWLVaHaZ5PdHQ0P/74I0uWLHE6iimhxNvu3hMeHq6JiYlOxzCtW2f/uWKFkynyderUKWrUqEFkZCTvv/++03EuWUZGBrVq1aJly5bMmzfP6TimGBGRDaoant+44vu7sDHAnDlz+O233+jdu7fTUTwiICCAnj17smDBAg4csOkv43lW+qZYi42NJSwszKuvpllQTz75JJmZmcTFxTkdxZRAVvqm2Pr2229Zu3YtUVFRiIjTcTwmODiYu+66i3feeccmdI3HWembYmvq1KmUKVOGbt26OR3F43r37s3+/fv59NNPnY5iShgrfVMs/f7773zwwQd07tyZKlWqOB3H4zp06ED16tWZMmWK01FMCWOlb4qlOXPmcOLECaKiopyOUij8/f15/PHHWbRoEfv3789/BWPcZKVviqXY2Fjq16/PLbfc4nSUQtOrVy9U1c7QNR5lpW+KnS1btrBmzZoSN4GbW1BQEO3ateOdd97hzJkzTscxJYSVvil2pk6dSkBAQIk4Azc/vXv35uDBgyxatMjpKKaEsNI3xcrZCdwHH3yQK6+80uk4he6+++6jRo0aNqFrPMZK3xQrH374IcePHy+xE7i5+fv788QTT/Dpp5/yww8/OB3HlABW+qZYiY2N5dprr+W2225zOkqR6dWrF4BN6BqPsNI3xcZ3333HqlWrSvwEbm61a9cmIiKCadOm2YSuuWRW+qbYODuB26NHD6ejFLmoqCib0DUeYaVvioXTp0/z/vvv06lTJ6666iqn4xS5e++9lxo1ahAbG+t0FFPMuVX6IhIhIkkikiwiQ/J4f5CIbBeRLSKyXERq53gvU0Q2ux4Judc1xh0fffQRv/76q89M4OZ29gzdTz/9lB9//NHpOKYYy7f0RaQ0MBFoD4QBXUQkLNewTUC4ql4HfASMzfHeaVVt4np0wJiLEBsbS0hICLfffrvTURzzxBNPoKpMmzbN6SimGHNnS785kKyqe1Q1g+wbn3fMOUBVv1DV310v15B9A3RjPGL79u18/fXXPjeBm9vZM3Tj4uJsQtdcNHdKvyaQ84pPKa5l5/MEkPN6sGVFJFFE1ohIZF4riEiUa0zi4cOH3YhkfMnUqVPx9/f3yQnc3Hr37k1KSgqLFy92Oooppjw6kSsiXYFw4LUci2u77tv4KPCmiATnXk9VY1U1XFXDq1at6slIppg7ffo006dP5+9//zv2dyP7DN3q1avbhK65aO6U/gGgVo7Xga5lfyIibYFhQAdVTT+7XFUPuP7cA6wAbriEvMbHzJ07l2PHjvnsBG5uZyd0Fy5cSEpKitNxTDHkTumvB0JEJEhEAoBHgD8dhSMiNwBTyC78QzmWXyEiZVzPrwJuBrZ7Krwp+WJjY6lXrx6tW7d2OorX6NWrF1lZWTahay5KvqWvqmeAvsASYAcQr6rbRGS0iJw9Guc1oCLwYa5DMxsAiSLyLfAF8IqqWukbt+zYsYOvvvqKXr16UaqUnVJyVs5LLts9dE1B+bkzSFUXAYtyLRuR43nb86y3Cmh8KQGN75oyZQr+/v707NnT6SheJyoqigcffJAlS5Zwzz33OB3HFCO2+WS8Us4J3GrVqjkdx+t06NCBq6++2iZ0TYFZ6RuvFB8fz6+//kp0dLTTUbzS2d+AFixYwIEDfzmuwpjzstI3XmnKlCmEhob69Bm4+enVqxeZmZm8++67TkcxxYiVvvE6W7ZsYfXq1fTu3dunz8DNT3BwMG3btrUJXVMgVvrG60yZMoUyZcrYGbhuiIqK4ocffmDp0qVORzHFhJW+8SonT57kgw8+oHPnzlSpUsXpOF6vY8eOVK1a1SZ0jdus9I1XmT17Nr/99ptN4LopICCAnj17kpCQQGpqqtNxTDFgpW+8ypQpU2jYsCGtWrVyOkqxYRO6piCs9I3X2LBhA4mJiURHR9sEbgGEhIRw5513MnXqVLKyspyOY7yclb7xGlOmTKFcuXJ07drV6SjFTlRUFPv27bMJXZMvK33jFU6cOMGsWbPo0qULlStXdjpOsRMZGUm1atWYOHGi01GMl7PSN15h5syZnDp1it69ezsdpVgqU6YMTz75JAsWLGDfvn1OxzFezErfOE5VmTRpEk2aNOHGG290Ok6x1bt3b0qVKsWkSZOcjmK8mJW+cdyXX37J1q1biYmJsQncS1CrVi06duxIXFwcp0+fdjqO8VJW+sZx48eP54orruDRRx91Okqx17dvX3755RfmzJnjdBTjpdwqfRGJEJEkEUkWkSF5vD9IRLaLyBYRWS4itXO810NEdrkedl69+ZP9+/czb948evXqRfny5Z2OU+y1bt2asLAwJkyYgKo6Hcd4oXxLX0RKAxOB9kAY0EVEwnIN2wSEq+p1wEfAWNe6VYCRQAugOTBSRK7wXHxT3E2aNAlV5amnnnI6SokgIsTExLBhwwbWrVvndBzjhdzZ0m8OJKvqHlXNAGYDHXMOUNUvVPV318s1ZN88HeBuYKmqHlXVY8BSIMIz0U1xl5aWxtSpU7n//vupU6eO03FKjG7dulGpUiUmTJjgdBTjhdwp/ZrA/hyvU1zLzucJ4NOCrCsiUSKSKCKJhw8fdiOSKQlmz57NkSNHePrpp52OUqJUqlSJHj16EB8fz6FDh5yOY7yMRydyRaQrEE72jdLdpqqxqhququFVq1b1ZCTjpVSV8ePHExYWxp133ul0nBLnqaeeIiMjg3feecfpKMbLuFP6B4BaOV4Hupb9iYi0BYYBHVQ1vSDrGt+zevVqNm7cSN++fe0wzULQoEED2rRpw7///W/++OMPp+MYL+JO6a8HQkQkSEQCgEeAhJwDROQGYArZhZ/z98klQDsRucI1gdvOtcz4uPHjx3P55ZfTrVs3p6OUWAMGDODAgQN8+OGHTkcxXiTf0lfVM0Bfsst6BxCvqttEZLSIdHANew2oCHwoIptFJMG17lFgDNk/ONYDo13LjA87ePAgH330ET179qRixYpOxymx7rnnHkJDQ3n99dft8E1zjp87g1R1EbAo17IROZ63vcC604BpFxvQlDwTJkwgMzOTvn37Oh2lRCtVqhQDBw4kOjqaL7/80m4ybwA7I9cUsZMnTzJ58mQeeOABgoODnY5T4nXv3p0rr7ySN954w+koxktY6Zsi9e6773Ls2DEGDx7sdBSfUK5cOZ566inmz5/Pzp07nY5jvICVvikymZmZ/Otf/6Jly5a0bNnS6Tg+IyYmBn9/f958802noxgvYKVviswnn3zC3r17bSu/iF199dV07dqV9957j19++cXpOMZhVvqmSKgq48aNIzg4mI4dO+a/gvGoQYMGcfr0aSZPnux0FOMwK31TJFatWsXatWsZOHAgpUuXdjqOz2nYsCERERFMmDCB9PT0/FcwJZaVvikSY8eOpUqVKjz22GNOR/FZgwYN4qeffmLGjBlORzEOstI3hW7r1q0kJCTQr18/KlSo4HQcn9W2bVuaNWvGK6+8QmZmptNxjEOs9E2he/nll6lYsaJdTdNhIsLzzz9PcnKyXZrBh1npm0KVnJzMnDlz6NOnD1WqVHE6js+LjIykQYMG/POf/yQrK8vpOMYBVvqmUL366qv4+/szaNAgp6MYsi/NMHToULZu3crChQudjmMcYKVvCs3+/fuZPn06vXr1onr16k7HMS5dunShTp06vPTSS3YhNh9kpW8Kzbhx41BVnnnmGaejmBz8/Px47rnnWLt2LZ9//rnTcUwRs9I3heKnn35i6tSpdO3aldq1azsdx+Ty2GOPUaNGDcaMGWNb+z7GSt8UipdffpmMjAyGDRvmdBSTh7JlyzJ06FBWrlxpW/s+xq3SF5EIEUkSkWQRGZLH+7eJyEYROSMiD+Z6L9N1Y5VzN1cxJVtKSgqTJ0/mscceo169ek7HMefx5JNPEhgYyAsvvGBb+z4k39IXkdLARKA9EAZ0EZGwXMN+BB4DZuXxEadVtYnr0SGP900Jc3aC8IUXXnA6irmAsmXLMnz4cFavXs3ixYudjmOKiDtb+s2BZFXdo6oZwGzgT1fMUtV9qroFsAN/fdy+ffuIi4ujV69eti+/GOjZsyd16tRhxIgRtrXvI9wp/ZrA/hyvU1zL3FVWRBJFZI2IROY1QESiXGMSDx8+XICPNt5mzJgxlCpVyvblFxMBAQGMGDGCxMREEhJs76svKIqJ3NqqGg48CrwpIn+5R56qxqpquKqGV61atQgimcKwa9cupk+fTnR0NDVrFmS7wDipW7duhISEMGLECDtL1we4U/oHgFo5Xge6lrlFVQ+4/twDrABuKEA+U4w8//zzlC1bliFD/jLXb7yYn58fo0aNYsuWLcycOdPpOKaQuVP664EQEQkSkQDgEcCt3wNF5AoRKeN6fhVwM7D9YsMa77V69Wo++ugjnnnmGTv7thh6+OGHadasGcOGDeP06dNOxzGFKN/SV9UzQF9gCbADiFfVbSIyWkQ6AIjIjSKSAnQGpojINtfqDYBEEfkW+AJ4RVWt9EsYVWXw4MFUr16df/zjH07HMRehVKlSjBs3jv379/P22287HccUIj93BqnqImBRrmUjcjxfT/Zun9zrrQIaX2JG4+U+/vhjVq1axdSpU6lYsaLTccxFat26Nffffz///Oc/efzxx7H5tZLJzsg1lyQjI4MhQ4bQsGFDevbs6XQcc4leffVVTp06xZgxY5yOYgqJlb65JFOmTCE5OZmxY8favW9LgAYNGtCrVy8mTZrErl27nI5jCoGVvrlohw8fZsSIEbRp04b27ds7Hcd4yIsvvki5cuUYOHCg01FMIbDSNxft+eef5+TJk7z99tuIiNNxjIdUr16dkSNHsnDhQubPn+90HONhVvrmoqxbt464uDj69+9PWFjuSzGZ4q5fv36EhYXRv39/O4SzhLHSNwWWmZlJTEzMuS1CU/L4+/szYcIE9u7dy9ixY52OYzzISt8U2LRp00hMTGTcuHFUqlTJ6TimkNxxxx08/PDDvPLKK+zdu9fpOMZDrPRNgRw+fJihQ4dy22230aVLF6fjmEI2btw4SpcuTb9+/ewqnCWElb4pkAEDBnDixAkmTpxok7c+IDAwkFGjRrFgwQLi4+OdjmM8wErfuG3BggXMmjWLYcOG0ahRI6fjmCLSv39/brzxRp5++mmOHDnidBxziaz0jVtOnDhBdHQ0jRo1YujQoU7HMUXIz8+PuLg4jh07ZsfulwBW+sYtzz33HKmpqcTFxREQEOB0HFPEGjduzPPPP8+MGTNYtGhR/isYr2Wlb/K1YsUKJk+ezIABA2jevLnTcYxDnn/+ecLCwujduzfHjx93Oo65SFb65oKOHTtG9+7dqVevHqNHj3Y6jnFQmTJlmDZtGqmpqcTExDgdx1wkK31zQTExMRw8eJCZM2dSoUIFp+MYh7Vo0YIXXniBmTNn8p///MfpOOYiuFX6IhIhIkkikiwif7kXnojcJiIbReSMiDyY670eIrLL9ejhqeCm8P3888/85z//4cUXX7TdOuacYcOG0bJlS/r06cOPP/7odBxTQPmWvoiUBiYC7YEwoIuI5L7Yyo/AY8CsXOtWAUYCLYDmwEgRueLSY5vClpaWxs5du7j55pvtaB3zJ35+fsyYMYPMzEy6d+9OZmam05FMAbizpd8cSFbVPaqaAcwGOuYcoKr7VHULkJVr3buBpap6VFWPAUuBCA/kNoXojz/+YPuOHQB88MEHdp188xd169Zl/PjxrFy5kpdfftnpOKYA3Cn9msD+HK9TXMvc4da6IhIlIokiknj48GE3P9oUlmeffZYTJ04Qeu21BAUFOR3HeKkePXrw6KOPMnLkSJYtW+Z0HOMmr5jIVdVYVQ1X1XC7L6ez4uPjefPNN6lZsybVqlVzOo7xYiJCbGwsDRo0oEuXLuzfvz//lYzj3Cn9A0CtHK8DXcvccSnrmiK2Y8cOnnjiCVq2bElwcLDTcUwxUKFCBebOnUtaWhqdO3cmIyPD6UgmH+6U/nogRESCRCQAeARIcPPzlwDtROQK1wRuO9cy42VOnDhBp06dKFeuHPHx8ZSyi6kZN4WGhvLuu++ydu1aBgwY4HQck498S19VzwB9yS7rHUC8qm4TkdEi0gFARG4UkRSgMzBFRLa51j0KjCH7B8d6YLRrmfEiZ86c4ZFHHmHnzp3Mnj2bwMBApyOZYubBBx9k8ODBTJo0iYkTJzodx1yAnzuDVHURsCjXshE5nq8ne9dNXutOA6ZdQkZTyAYNGsSnn37K5MmTufPOO52OY4qpV155haSkJPr160dwcDAREXagnjfyiolc45yJEycyfvx4Bg0aRO/evZ2OY4qx0qVLM2vWLK677joeeughvvvuO6cjmTxY6fuwRYsW0a9fPzp06GD3QTUeUbFiRebPn0/FihW57777SE1NdTqSycVK30d99dVXdOrUiSZNmjBz5kw7Act4TGBgIPPnz+fIkSO0a9eOo0dtGs+bWOn7oE2bNnHfffdRu3ZtFi9eTMWKFZ2OZEqYZs2a8d///pedO3dyzz33cPLkSacjGRcrfR+TlJTE3XffTeXKlVm6dCl2MpwpLG3atGHOnDkkJiYSGRlJenq605EMVvo+Zee1knB9AAAPzElEQVTOnbRp0waApUuXUqtWrXzWMObSREZGEhcXx/Lly+nUqRNpaWlOR/J5Vvo+YseOHdx+++2kp6ezbNkyrr32WqcjGR/Ro0cPJk2axMKFC+nYsSOnT592OpJPs9L3AVu3buX2229HVVmxYgXXXXed05GMj4mOjiYuLo6lS5dy7733curUKacj+Swr/RJu7dq13HHHHfj7+7Ny5UoaNmzodCTjox5//HHef/99Vq5cyd13321H9TjESr8ES0hI4I477uDyyy/nyy+/JDQ01OlIxsd17dqV2bNns379em655RZ++OEHpyP5HCv9Emry5Mk88MADNGrUiNWrV9tVM43X6Ny5M5999hmpqancdNNNbN682elIPsVKv4TJzMzk2WefpU+fPrRv354vvvjCrotvvM7tt9/O119/jb+/P7feeivz5893OpLPsNIvQY4ePUr79u157bXXeOqpp5g3bx4VKlRwOpYxeWrYsCGrV68mNDSUDh06MHr0aLKyct9x1XialX4JsXXrVm688UZWrlzJO++8w8SJE/Hzc+siqsY4pmbNmnz11Vd069aNkSNH0qlTJ06cOOF0rBLNSr+YU1WmTp1KixYtSEtLY+XKlTzxxBNOxzLGbeXKlWP69Om8+eabzJ8/n2bNmrFhwwanY5VYbpW+iESISJKIJIvIkDzeLyMic1zvrxWROq7ldUTktIhsdj0meza+bzt27BgPPfQQUVFR3HLLLWzYsIGbbrrJ6VjGFJiI0L9/f7744gvS0tJo2bIlr7/+uu3uKQT5lr6IlAYmAu2BMKCLiITlGvYEcExV6wH/Al7N8d5uVW3iekR7KLfP+/zzz2nSpAnz5s1j7NixLF68mOrVqzsdy5hLcuutt/Ltt99y3333MXjwYNq3b283XPcwd7b0mwPJqrpHVTOA2UDHXGM6AtNdzz8C2ojYTVYLw4kTJ4iOjqZNmzaUKVOGb775hmeeeYZSpWxPnSkZqlSpwty5c5k0aRJff/01DRs2ZOrUqaiq09FKBHeaoiaQ80dtimtZnmNc99Q9Dlzpei9IRDaJyEoRufUS8/q0hQsX0qhRI6ZOncrgwYPZvHkzzZs3dzqWMR4nIkRHR7N161bCw8OJioqiXbt2JCcnOx2t2CvszcNU4G+qegMwCJglIpflHiQiUSKSKCKJhw8fLuRIxc+ePXvo0KED9913HxUrVmTVqlW89tprlC9f3uloxhSqunXrsmzZMiZPnszatWtp2LAhL7zwAr///rvT0Yotd0r/AJDzGryBrmV5jhERP+By4BdVTVfVXwBUdQOwG/jL5R1VNVZVw1U13K7v/v9+++03Ro4cSVhYGJ9//jljx45l8+bNtGjRwuloxhSZUqVK0bt3b5KSknjooYf43//9Xxo0aEB8fLxN9F4Ed0p/PRAiIkEiEgA8AiTkGpMA9HA9fxD4XFVVRKq6JoIRkbpACLDHM9FLrvT0dN5++22Cg4MZPXo0f//730lKSuKZZ54hICDA6XjGOKJGjRp88MEHfPnll1SuXJmHH36Y5s2bs3TpUtvfXwD5lr5rH31fYAmwA4hX1W0iMlpEOriGxQFXikgy2btxzh7WeRuwRUQ2kz3BG62qdmm988jMzOT9998nNDSU/v3707hxY9atW8esWbOoWTP3NIoxvunWW29l48aNTJ8+/dx9eNu2bcu6deucjlY8qKpXPZo1a6a+5vTp0zplyhStV6+eAtq0aVP97LPPNCsry7lQt9+e/TDGi6Wlpelbb72lVatWVUDbtWuny5cvd/bfjkOARHWjY+04PwcdP36cV199laCgIHr37k3lypX56KOPWL9+PXfddRd21KsxF1amTBn69evH7t27efnll/n2229p06YNLVq04OOPPyYzM9PpiF7HSt8BW7ZsoU+fPgQGBjJkyBAaN27MsmXLWLduHZ06dbJj7o0poEqVKjFkyBD27dvH5MmTOXr0KJ06dSIkJIRXX30VOyrw/1m7FJHTp08zY8YMbr75Zq6//nree+89OnXqRGJiIp999hlt2rSxLXtjLlHZsmXPHekTHx/P3/72N4YMGUJgYCD/8z//w1dffeXzk75W+oUoMzOTzz//nMcff5yrr76abt26ceTIEd544w0OHDjAe++9R7NmzZyOaUyJU7p0aTp37syKFSvYtm0b0dHRLFy4kNtuu43g4GBGjBhBUlKS0zEdId72Uy88PFwTExOdjnHRMjMzWbNmDZ988gmzZ8/mwIEDVKpUiQcffJBu3brRunXr4rFF37p19p8rVjiZwhiPOXXqFHPnzmXGjBksX76crKwsbrzxRh5++GEiIyOL/d3lRGSDqobnO85K/9KlpaWxfPly5s2bR0JCAocOHcLf35+7776bbt26cf/991OuXDmnYxaMlb4pwQ4ePMjs2bOZMWMGmzZtAqBRo0ZERkYSGRlJ06ZNi8fGWQ5W+oVIVfn+++9ZtmwZy5YtY/ny5Zw6dYpKlSpxzz33EBkZSfv27bn88sudjnrxrPSNj9i7dy///e9/mTdvHl999RVZWVkEBgaeO/6/TZs2xeKWo1b6HqSq7N27l2+++eZc0R88eBCA4OBg7rrrLiIjI2ndujVlypRxOK2HWOkbH3TkyBEWLlzI/PnzWb58Ob/++isA119/PW3btuWOO+6gZcuWVKlSxeGkf2WlfwnS09PZtGkTq1at4ptvvmHVqlX89NNPAFx55ZW0adOGtm3b0rZtW4KCghzNWmis9I2Py8zMZOPGjSxbtoylS5fyzTffkJGRAUBYWBitWrXi5ptvplWrVoSEhDi+O8hK300nT57k22+/ZdOmTWzcuJFNmzaxbds2/vjjDyD7Kn+tWrU692jcuLFvHEdvpW/Mn/z++++sW7fu3Mbg6tWrOXbsGACXXXYZN9xwAzfccANNmzblhhtuoH79+kV6n2p3S99n7px96tQpvv/+e77//nt27NjBjh072LZtGzt37jx33G7VqlVp2rQp7du3Jzw8nJYtW1KjRg2HkxtjvEH58uVp3bo1rV0bRFlZWSQlJbFq1So2btzIxo0bmTJlCqdPnwayzxlo3LgxYWFhhIWF0aBBAxo0aEBQUBClS5d27PsoUVv6J0+eZN++fezdu5e9e/eyZ8+ecyX/448/nhtXunRpgoODCQsL+9NP5muuucbxX9G8hm3pG1NgmZmZJCUlndtzsHnzZnbs2EFqauq5MWXKlCE0NJTQ0FCCg4OpW7cudevWJTg4mMDAwIv+7aDE7d7JzMzk0KFDpKamcvDgQQ4cOHCu3Pfu3cu+ffv+cqp1+fLlCQ0NPfcTtkGDBtSvX5+QkBC7RHF+rPSN8Zhff/313B6GHTt2sH37dnbt2sW+ffvO7UoG8PPzo3bt2tStW5fatWsTGBhIYGAgNWvWPPf88ssvz3PjtNju3jly5Ahjxozh4MGDf3r8/PPPf7l4kr+/P7Vr1yYoKIjIyEiCgoLOPerUqUO1atVsy90Y47jKlSvTsmVLWrZs+aflmZmZpKSksGfPnnOP3bt3s3v3brZs2cLPP//8l88qX778uR8E11xzDdWrV6d69epuZ/G6LX0RUYCrrrqKa6655ryPGjVqUKNGDUf3jZVotqVvjOMyMjJITU0lJSXl3OPAgQPnnqempvLTTz+RlpYGUDy39Bs3bsz69etLzvHuxhhzkQICAqhduza1a9c+7xhV5bfffnP7ZFC3jj0UkQgRSRKRZBEZksf7ZURkjuv9tSJSJ8d7Q13Lk0Tk7vy+VkBAgBW+Mca4SUS47LLL3B6fb+m77nE7EWgPhAFdRCQs17AngGOqWg/4F/Cqa90wsu+p2xCIAP599p65xhhjip47W/rNgWRV3aOqGcBsoGOuMR2B6a7nHwFtJHsGtSMwW1XTVXUvkOz6PGOMMQ5wZ59+TWB/jtcpQIvzjVHVMyJyHLjStXxNrnX/codvEYkColwv00XkO7fSO+sq4IjTIdxwaTmL7ugn3/jvWXQsp+cUh4wAoe4M8oqJXFWNBWIBRCTRnRlop1lOz7KcnmU5Pac4ZITsnO6Mc2f3zgGgVo7Xga5leY4RET/gcuAXN9c1xhhTRNwp/fVAiIgEiUgA2ROzCbnGJAA9XM8fBD7X7BMAEoBHXEf3BAEhwDrPRDfGGFNQ+e7ece2j7wssAUoD01R1m4iMBhJVNQGIAz4QkWTgKNk/GHCNiwe2A2eAGFXNzPML/b/Yi/92ipTl9CzL6VmW03OKQ0ZwM6fXnZFrjDGm8PjAheGNMcacZaVvjDE+xKtLX0T+ISIqIlc5nSUvIjJGRLaIyGYR+UxErnE6U15E5DUR+d6V9RMRqex0pryISGcR2SYiWSLiVYfI5XcpEm8hItNE5JA3n+siIrVE5AsR2e76/93f6Ux5EZGyIrJORL515RzldKYLEZHSIrJJRBZcaJzXlr6I1ALaAT/mN9ZBr6nqdaraBFgAjHA60HksBRqp6nXATmCow3nO5zvg78CXTgfJyc1LkXiL98i+5Ik3OwP8Q1XDgJuAGC/975kO3Kmq1wNNgAgRucnhTBfSH9iR3yCvLX2yr+HzLOC1M82qeiLHywp4aVZV/UxVz7heriH7fAmvo6o7VDXJ6Rx5cOdSJF5BVb8k+wg6r6Wqqaq60fX8N7KL6i9n6jtNs510vfR3Pbzy37iIBAL3Au/kN9YrS19EOgIHVPVbp7PkR0ReEpH9wP/gvVv6OT0OfOp0iGImr0uReF1JFUeuK/LeAKx1NkneXLtMNgOHgKWq6pU5gTfJ3kjOym+gY5dhEJFlQF63exkGPE/2rh3HXSinqv5XVYcBw0RkKNAXGFmkAV3yy+kaM4zsX61nFmW2nNzJaXyDiFQE5gIDcv3W7DVc5xU1cc2DfSIijVTVq+ZLROQ+4JCqbhCR1vmNd6z0VbVtXstFpDEQBHzrutVhILBRRJqr6k9FGBE4f848zAQW4VDp55dTRB4D7gPaqIMnZxTgv6c3scuJeJiI+JNd+DNV9WOn8+RHVX8VkS/Ini/xqtIHbgY6iMg9QFngMhGZoapd8xrsdbt3VHWrqlZT1TqqWofsX6WbOlH4+RGRkBwvOwLfO5XlQkQkguxf/Tqo6u9O5ymG3LkUiXGT67LrccAOVX3D6TznIyJVzx7pJiLlgLvwwn/jqjpUVQNdffkI2ZfBybPwwQtLv5h5RUS+E5EtZO+O8spDz4AJQCVgqevw0slOB8qLiDwgIilAS2ChiCxxOhNkX4qE7F13S8iedIxX1W3OpsqbiPwHWA2EikiKiDzhdKY83Ax0A+50/X3c7NpK9TY1gC9c/77Xk71P/4KHQxYHdhkGY4zxIbalb4wxPsRK3xhjfIiVvjHG+BArfWOM8SFW+sYY40Os9I0xxodY6RtjjA/5Pz8TDqCE91eaAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "df=3\n",
    "\n",
    "x = np.linspace(-4, 4, 300)\n",
    "y = t.pdf(x, df)\n",
    "plt.plot(x, y, 'k')\n",
    "plt.plot(np.array([t.ppf(0.1,df),t.ppf(0.1,df)]),np.array([0,1]),color='red')\n",
    "\n",
    "plt.xlim([-4,4])\n",
    "plt.ylim([0,0.45])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Comparison for the percent point function for the probabilities of 0.1 and 0.9.\n",
      "0.1 probability:  -1.6377443572159065\n",
      "0.9 probability:  1.637744357215907\n",
      "These results show that the t-distribution is symmetric.\n"
     ]
    }
   ],
   "source": [
    "print(\"Comparison for the percent point function for the probabilities of 0.1 and 0.9.\")\n",
    "print(\"0.1 probability: \",t.ppf(0.1,df))\n",
    "print(\"0.9 probability: \",t.ppf(0.9,df))\n",
    "print(\"These results show that the t-distribution is symmetric.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cumulative density function - Find the probability of a specific point."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can use <strong>cumulative density function (cdf)</strong>. If you want to find the probability at the point of -1, you can type:\n",
    "`t.cdf(-0.1, df)`\n",
    "\n",
    "See also [scipy.stats.t](https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.stats.t.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.19550110947788527"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t.cdf(-1, df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The percent cumulative density function is the inverse of the cumulative density function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cdf at x = -1: 0.196\n",
      "ppf of cdf (exact value = -1): -1.000\n"
     ]
    }
   ],
   "source": [
    "print(\"cdf at x = -1: %0.3f\"%t.cdf(-1, df))\n",
    "print(\"ppf of cdf (exact value = -1): %0.3f\"%t.ppf(t.cdf(-1, df) , df))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [Confidence interval](https://en.wikipedia.org/wiki/Confidence_interval) for population mean"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example\n",
    "\n",
    "Let us assume that we want to know the average starting salary of all STEM students who graduated last year in Switzerland. Idealy, it is nice to ask all of the STEM students. However, it is almost impossible. Instead of asking all of the STEM students, we can randomly select n students to make a sample.\n",
    "\n",
    "- Population: all STEM students in Switzerland who graduated last year. Mean = &mu;, Std = &sigma;\n",
    "- Sample: random sample of n students. Mean = x_bar, Std = s.\n",
    "\n",
    "The confidence interval shows the lower and upper limit of the population mean &mu; with the chosen confidence level (%).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## z-statistic"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### General concept"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "According to the central limit theorem, the distribution of the sample mean is described by the normal distribution with N(&mu;, &sigma; (n)<sup>-1/2</sup>):\n",
    "$$\n",
    "    \\bar{x} = N\\biggl(\\mu,\\frac{\\sigma}{\\sqrt{n}}\\biggr)\n",
    "$$\n",
    "This relation can be rewritten as\n",
    "$$\n",
    "    \\frac{\\bar{x}-\\mu}{\\sigma/\\sqrt{n}} = N(0,1).\n",
    "$$\n",
    "Here,\n",
    "$$\n",
    "    \\frac{\\bar{x}-\\mu}{\\sigma/\\sqrt{n}}\n",
    "$$\n",
    "is the z-statistic and equal to the standard normal distribution. As you can see, <strong>the z-statistic contains the population standard deviation &sigma;</strong>. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A (1-&alpha;) confidence interval for the population mean &mu; is given by\n",
    "$$\n",
    "    \\bar{x}-|z_{\\alpha/2}|\\frac{\\sigma}{\\sqrt{n}} < \\mu < \\bar{x}+|z_{\\alpha/2}|\\frac{\\sigma}{\\sqrt{n}}\n",
    "$$\n",
    "Here, \n",
    "$$\n",
    "    |z_{\\alpha/2}|\\frac{\\sigma}{\\sqrt{n}}\n",
    "$$\n",
    "is the margin of error, and &alpha; is the probability outside the confidence level."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculation example\n",
    "Let us calculate the 85 % confidence interval for the population mean &mu;. Let us assume that the population standard deviation is &sigma; = 10 and that the number of samples is n = 20. The sample mean x_bar is 70."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "|Z_alpha/2| = 1.440\n",
      "Lower limit: 66.781\n",
      "Upper limit: 73.219\n",
      "\n",
      "A 85.0 % confidence interval for the population mean is: \n",
      " 66.78 < mu < 73.22\n"
     ]
    }
   ],
   "source": [
    "from scipy.stats import norm\n",
    "\n",
    "alpha=0.15\n",
    "sigma = 10\n",
    "n = 20\n",
    "\n",
    "z = norm.ppf(0.5*alpha) \n",
    "print(\"|Z_alpha/2| = %0.3f\"%np.abs(z))\n",
    "\n",
    "me= np.abs(z) * sigma / np.sqrt(n)  # margin of error\n",
    "print(\"Lower limit: %0.3f\"%(70-me))\n",
    "print(\"Upper limit: %0.3f\"%(70+me))\n",
    "print(\"\\nA %0.1f %% confidence interval for the population mean is: \\n %0.2f < mu < %0.2f\"%((1-alpha)*100,70-me,70+me))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us caluclate the confidence intervals for different confidence levels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "A 85.0 % confidence interval for the population mean is: \n",
      " 66.78 < mu < 73.22\n",
      "\n",
      "A 95.0 % confidence interval for the population mean is: \n",
      " 65.62 < mu < 74.38\n",
      "\n",
      "A 99.0 % confidence interval for the population mean is: \n",
      " 64.24 < mu < 75.76\n"
     ]
    }
   ],
   "source": [
    "sigma = 10\n",
    "n = 20\n",
    "\n",
    "alpha_list=[0.15,0.05, 0.01] # 85, 95, and 99 % confidence interval\n",
    "\n",
    "for j in range(len(alpha_list)):\n",
    "    alpha = alpha_list[j]\n",
    "    z = norm.ppf(0.5*alpha) \n",
    "    me= np.abs(z) * sigma / np.sqrt(n)  # margin of error\n",
    "    print(\"\\nA %0.1f %% confidence interval for the population mean is: \\n %0.2f < mu < %0.2f\"%((1-alpha)*100,70-me,70+me))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the lower and upper limits become infinite, a 100 % confidence interval can be achieved."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Size calculation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In some cases, it is necessary to determine the required size n to get a specific confidence interval. Let us assume the following conditions:\n",
    "- Margin of error: +/- 0.3 \n",
    "- Confidence level: 95 % (this is a commonly used setting)\n",
    "- Population standard deviation: 0.9\n",
    "\n",
    "Recall that the margin of error can be written as\n",
    "$$\n",
    "    ME = |z_{\\alpha/2}|\\frac{\\sigma}{\\sqrt{n}}.\n",
    "$$\n",
    "From this equation, we can get\n",
    "$$\n",
    "    n = \\biggl(|z_{\\alpha/2}|\\frac{\\sigma}{ME}\\biggr)^2.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 120,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "In order to get the specified confidence interval with 95 % confidence level, 34.6 samples are needed.\n"
     ]
    }
   ],
   "source": [
    "alpha=0.05\n",
    "sigma = 0.9 # population standard deviation\n",
    "me = 0.3 # margin of error \n",
    "\n",
    "z = norm.ppf(0.5*alpha) # z-statistc\n",
    "\n",
    "n = (np.abs(z) * sigma / me )**2\n",
    "print(\"In order to get the specified confidence interval with 95 %% confidence level, %0.1f samples are needed.\"%n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## t-statistic\n",
    "### Gemeral concept\n",
    "When the population standard deviation &sigma; is not available, we replace &sigma and z-statistic (z<sub>&alpha;/2</sub>) by sample standard deviation s and t-statistic (t<sub>&alpha;/2</sub>), respectively.  \n",
    "\n",
    "A (1-&alpha;) confidence interval for the population mean &mu; is given by\n",
    "$$\n",
    "    \\bar{x}-|t_{\\alpha/2}|\\frac{s}{\\sqrt{n}} < \\mu < \\bar{x}+|t_{\\alpha/2}|\\frac{s}{\\sqrt{n}}.\n",
    "$$\n",
    "\n",
    "The t-statistic |t<sub>&alpha;/2</sub>| can be calculated by percent position function (ppf):   \n",
    "`scipy.stats.t.ppf(alpha/2, df=n-1)`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculation example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "|t_alpha/2| = 1.500\n",
      "Lower limit: 66.645\n",
      "Upper limit: 73.355\n",
      "\n",
      "A 85.0 % confidence interval for the population mean is: \n",
      " 66.65 < mu < 73.35\n"
     ]
    }
   ],
   "source": [
    "from scipy.stats import t\n",
    " \n",
    "alpha=0.15 \n",
    "n = 20 # number of samples\n",
    "s=10 # sample standard deviation\n",
    "\n",
    "t_statistic = t.ppf(0.5*alpha, df=n-1) # t-statistic\n",
    "print(\"|t_alpha/2| = %0.3f\"%np.abs(t_statistic))\n",
    "\n",
    "me= np.abs(t_statistic) * s / np.sqrt(n)  # margin of error\n",
    "print(\"Lower limit: %0.3f\"%(70-me))\n",
    "print(\"Upper limit: %0.3f\"%(70+me))\n",
    "print(\"\\nA %0.1f %% confidence interval for the population mean is: \\n %0.2f < mu < %0.2f\"%((1-alpha)*100,70-me,70+me))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Confidence interval for population proportion\n",
    "## General concept\n",
    "A (1-&alpha;) confidence interval for the population proportion p is given by\n",
    "$$\n",
    "    \\hat{p}-|z_{\\alpha/2}|\\sqrt{\\frac{\\hat{p}(1-\\hat{p})}{n}} < p < \\hat{p}+|z_{\\alpha/2}|\\sqrt{\\frac{\\hat{p}(1-\\hat{p})}{n}} \n",
    "$$\n",
    "Here, p and p_hat are the population proportion and the sample proportion from our sample, respectively. Only z-statistic is defined.\n",
    "\n",
    "Let us calculate a 90 % confidence interval for an example. The number of sample data is 210. The sample proportion is 51 %."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "|Z_alpha/2| = 1.645\n",
      "Lower limit: 0.453\n",
      "Upper limit: 0.567\n",
      "\n",
      "A 90.0 % confidence interval for the population mean is: \n",
      " 0.453 < p < 0.567\n"
     ]
    }
   ],
   "source": [
    "from scipy.stats import norm\n",
    "\n",
    "alpha=0.1\n",
    "n = 210\n",
    "p_hat=0.51\n",
    "\n",
    "z = norm.ppf(0.5*alpha) \n",
    "print(\"|Z_alpha/2| = %0.3f\"%np.abs(z))\n",
    "\n",
    "me= np.abs(z) * np.sqrt(p_hat*(1-p_hat)/n)  # margin of error\n",
    "print(\"Lower limit: %0.3f\"%(p_hat-me))\n",
    "print(\"Upper limit: %0.3f\"%(p_hat+me))\n",
    "print(\"\\nA %0.1f %% confidence interval for the population mean is: \\n %0.3f < p < %0.3f\"%((1-alpha)*100,p_hat-me,p_hat+me))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Size calculation\n",
    "The margin of error is \n",
    "$$\n",
    "    ME = |z_{\\alpha/2}|\\sqrt{\\frac{\\hat{p}(1-\\hat{p})}{n}}.\n",
    "$$\n",
    "From this equation, we can get \n",
    "$$\n",
    "    n = |z_{\\alpha/2}|^2\\frac{\\hat{p}(1-\\hat{p})}{ME^2}.\n",
    "$$\n",
    "We can determine the z-statistic and margin of error. However, we cannot know p_hat until we start investigation. Because of this reason, we can make a conservative estimate. In the equation above, n is proportional to\n",
    "$$\n",
    "    \\hat{p}(1-\\hat{p}).\n",
    "$$\n",
    "This factor takes the largest value when p_hat = 0.5. In the conservative estimate, p_hat = 0.5 is used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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.1"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}