{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Testing Distributions \n",
"\n",
"Knowing what distribution your data follow is important for knowing which statistical tests to apply to it. If you are using statistical tests that assume data to have a particular distribution, you need to test whether this is indeed the case.\n",
"\n",
"As we saw previously, different statistical distributions have different 'shapes'. But just looking at the shape isn't enough for saying definitively what distribution the data follows. Here, we will explore some statistical tests for checking the distribution of data. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"
\n",
"By 'testing distributions' we mean statistical tests that evaluate whether observed data follow a particular distribution.\n",
"
"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"# Imports\n",
"%matplotlib inline\n",
"\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"import scipy.stats as stats\n",
"from scipy.stats import normaltest"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# Set up a helper function for checking p-values against an alpha level, and printing result\n",
"def check_p_val(p_val, alpha):\n",
"\n",
" if p_val < alpha:\n",
" print('We have evidence to reject the null hypothesis.')\n",
" else:\n",
" print('We do not have evidence to reject the null hypothesis.')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Evaluating a Normal Distribution"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here, we will focus on the most common case: testing whether a dataset is normally distributed. "
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"# Create a dataset of normally distributed data\n",
"d1 = stats.norm.rvs(size=100000)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXxU5dnA/d81kwQQIzsqm4kWRBRkiaCIsogIqGAtKtpWrFoe3kpX+1hsVXixWh/11Yraqq1LaxW0qIgILqyCgBBIANl3CGtIWBOyzMz1/jHDZBJCMglJzizX9/OZD+c+5z4zV4bkmjP3uRdRVYwxxsQul9MBGGOMqV2W6I0xJsZZojfGmBhnid4YY2KcJXpjjIlxCU4HUFbz5s01JSXF6TCMMSaqrFix4pCqtijvWMQl+pSUFNLT050OwxhjooqI7DzTMWu6McaYGGeJ3hhjYpwlemOMiXER10ZvjIlcxcXFZGVlUVBQ4HQocat+/fq0adOGxMTEsM+xRG+MCVtWVhbJycmkpKQgIk6HE3dUlZycHLKyskhNTQ37PGu6McaEraCggGbNmlmSd4iI0KxZsyp/o7JEb4ypEkvyzqrO+x9WoheRwSKyUUS2iMi4CuqNEBEVkbSQfY8GztsoIjdVOUJjjDFnpdJELyJu4FVgCNAJuFtEOpVTLxn4FfBdyL5OwEjgcmAw8LfA8xljTLWICA8//HCw/PzzzzNhwoQ6jeG+++5j6tSp5e5PTU3lyiuvpEOHDtx7773s2bMneHzo0KEcOXLkjM/717/+lfz8/BqPN5wr+p7AFlXdpqpFwBRgeDn1ngSeBUIbj4YDU1S1UFW3A1sCz2dMZJvQqPyHcVy9evX4+OOPOXToULXO93g8NRxRac899xyrVq1i48aNdOvWjf79+1NUVATAzJkzady48RnPra1EH06vm9bA7pByFtArtIKIdAPaquoMEfl9mXOXljm3ddkXEJHRwGiAdu3ahRe5MbUg50QhmbuPkFF8Bxu0LRdKLre5v6W7bMaapiNDQkICo0eP5sUXX+Spp54qdWznzp3cf//9ZGdn06JFC95++23atWvHfffdR9OmTcnIyKB79+4kJyezfft29u3bx6ZNm3jhhRdYunQps2bNonXr1nz22WckJiYyceJEPvvsM06ePEnv3r15/fXXw24jFxF++9vf8sknnzBr1iyGDx8enOKlQYMG3HnnnWRlZeH1enn88cc5cOAAe/fupX///jRv3px58+bV3HsWTrzl7AuuPygiLuBF4L6qnhvcofoG8AZAWlqarW1o6lzOiUIembqaORsOBvb8MHjsXe8gUmUfI9zfcF+hh4b1rFcyQMq4z2vtuXc8c3OFxx966CG6dOnCI488Umr/2LFjuffeexk1ahRvvfUWv/rVr5g2bRoAmzZtYvbs2bjdbiZMmMDWrVuZN28e69at45prruGjjz7i2Wef5Yc//CGff/45t912G2PHjuWJJ54A4Kc//SkzZszg1ltvrdLP0r17dzZs2MDw4SUNIV988QWtWrXi88/97+HRo0dp1KgRL7zwAvPmzaN58+ZVeo3KhNN0kwW0DSm3AfaGlJOBK4D5IrIDuBqYHrghW9m5xjhuTdZRbn15UUiSP912vZDnPHcx8o2l5OYV1WF0pjznnXce9957L5MmTSq1f8mSJdxzzz2APzEvWrQoeOyOO+7A7S65RThkyBASExPp3LkzXq+XwYMHA9C5c2d27NgBwLx58+jVqxedO3dm7ty5rF27tsqxlrcud+fOnZk9ezZ/+MMfWLhwIY0a1W6zYDiJfjnQXkRSRSQJ/83V6acOqupRVW2uqimqmoK/qWaYqqYH6o0UkXoikgq0B5bV+E9hTDX9N303P3ptMXuP+m8tiUD3do15wD2T5xP/zkj3XJIpaTNds+cod72+hAPHbGSo037zm9/w5ptvkpeXd8Y6oc0sDRs2LHWsXr16ALhcLhITE4N1XS4XHo+HgoICfvGLXzB16lTWrFnDz3/+82qNCM7IyOCyyy4rta9Dhw6sWLGCzp078+ijjzJx4sQqP29VVPodVFU9IjIW+BJwA2+p6loRmQikq+r0Cs5dKyIfAusAD/CQqnprKHZjzsq7S3bw+KclV2jJ5PFSwqsMOJgJgdHlI9wLmZDwL/7jHchTnh+juNh88AR3vLaE9x7sRdum5zgTfASorHmltjVt2pQ777yTN998k/vvvx+A3r17M2XKFH7605/y3nvv0adPn2o//6mk3rx5c06cOMHUqVMZMWJE2OerKi+//DL79u0Lfls4Ze/evTRt2pSf/OQnnHvuubzzzjsAJCcnc/z4cUeablDVmaraQVUvUdWnAvueKC/Jq2q/wNX8qfJTgfMuVdVZNRe6MdX3/Z6jPDljfbDcXrKYnvQ4A9yZp9WtL8U8mDCLlxJfJcHlv+rblZvPXa8v4Ui+NeM46eGHHy7V+2bSpEm8/fbbdOnShXfffZeXXnqp2s/duHFjfv7zn9O5c2duu+02rrrqqrDO+9///d9g98rly5czb948kpKSStVZs2YNPXv2pGvXrjz11FM89thjAIwePZohQ4bQv39/AB588MEaWZ9Dyms/clJaWprawiOmNp0o9HDLpIXsyPE3yVwu2/kg6UnOlcq/ls/xduP/Kf41Rfj/cG9zLeKvSX+DCUdrNeZIsX79+tOaIUzdK+//QURWqGpaefVtCgQTV1SVP368JpjkG3KSVxJfDivJA9zgzuDlxFeC5Wm+PnzpLfdvy5iIYYnexJUP03czfVVJx6+nE/9Jqmt/lZ7jJnc6t7sWBst/Kn7AeuKYiGaJ3sSN7OOFTPxsXbB8V1pbhruXVOu5xif+m5YcBuAQjRg/verd7oypK5boTdx4Ze5m8or8nb4ubtGQCcMur/ZzNZI8nkn8R7D82aq9fLm2at8MjKkrluhNXNiVk8/7y3YFy38aehkNks5ufr0B7kxGuBcEy89+sQGvL7I6NxgDluhNnPj/vt5IsdefhK9KacKAji1r5HkfS/hPcEDV1uw8Zqy2gd8m8tikHSbmrd17lE8zSxLwuCEda2zxjMaSx8/cXzDJezsAk+Zs5pYurXC74mQGtJqe0TPMbqqffPIJt99+O+vXr6djx44V1n3nnXcYNGgQrVq1qlZI8+fP5/nnn2fGjBmn7R8+fDgXX3wx+fn5nH/++TzyyCPccsstlT5fUlISvXv3rlY81WFX9CbmPfvFxuD2QFc6Pd5OrdFphx9ImEVyYKIzu6qvG5MnT6ZPnz5MmTKl0rrvvPMOe/fWzv/JddddR0ZGBhs3bmTSpEmMHTuWOXPmVHjO/PnzWbx4ca3EcyaW6E1MW7othwWbsgFw4eORhA9q/DUaSR4/61OyUPOkOZutrb4WnThxgm+//ZY333zztET/7LPP0rlzZ6688krGjRvH1KlTSU9P58c//jFdu3bl5MmTpKSkBEfTpqen069fPwCWLVtG79696datG71792bjxo1lX7pCXbt25YknnuCVV/zjLD777DN69epFt27dGDhwIAcOHGDHjh289tprvPjii3Tt2pWFCxeWW6+mWaI3Me0f32wLbt/uXkgH154KalffA9em2lV9HZk2bRqDBw+mQ4cONG3alJUrVwIwa9Yspk2bxnfffceqVat45JFHGDFiBGlpabz33ntkZmbSoEGDMz5vx44d+eabb8jIyGDixIn88Y9/rHJsp6YkBujTpw9Lly4lIyODkSNH8uyzz5KSksKYMWP47W9/S2ZmJtddd1259WqatdGbmLUzJ4+5G0umHv6F+9Nae61G5yTysz6pTJqzGYCX5mzm1i6tcMVLW30dmjx5Mr/5zW8AGDlyJJMnT6Z79+7Mnj2bn/3sZ5xzjn+iuaZNm1bpeY8ePcqoUaPYvHkzIkJxcXGVYwudUiYrK4u77rqLffv2UVRURGpqarnnhFvvbNgVvYlZ/16yk1N/d/1cmVxcxRGwVTKhEQ8s6ksy/ilzt2XnsXD89bX3enEqJyeHuXPn8uCDD5KSksJzzz3HBx98gKqiqmHdZE9ISMDn8wGUmnb48ccfp3///nz//fd89tlnZz0l8S9/+UvGjh3LmjVreP3118/4fOHWOxuW6E1Myiv08OHykhUw73N/Ueuv2UjyuTOkX/1/vANr/TXjzdSpU7n33nvZuXMnO3bsYPfu3aSmprJo0SIGDRrEW2+9FVxzNTc3FyiZ+veUlJQUVqxYAcBHH30U3H/06FFat/avdHpq2uCqWL16NU8++SQPPfTQac/3r3/9K1ivbDxnqleTrOnGxKSPM/ZwvNC/CPTFzRty/fE1dfK697jn8KZ3KABzfN3Ze+QkrRqfuV046tXxrJ2TJ09m3Lhxpfb96Ec/4v333+fvf/87mZmZpKWlkZSUxNChQ3n66ae57777GDNmDA0aNGDJkiWMHz+eBx54gKeffppevUqWv37kkUcYNWoUL7zwAgMGDAgrnoULF9KtWzfy8/Np2bIlkyZN4oYbbgBgwoQJ3HHHHbRu3Zqrr76a7du3A3DrrbcyYsQIPv30U15++eUz1qtJNk2xiTmqysAXFrA129+MMuHWTtz3ddc6e/17iv7IYt8VAPzqhvb87sYOdfbatc2mKY4MNk2xiXuLthwKJvlz6yXwox5t6vT1f+KeHdyesmwXxV5fnb6+MWWFlehFZLCIbBSRLSIyrpzjY0RkjYhkisgiEekU2J8iIicD+zNF5LWa/gGMKeudb3cEt0f0aENy/cQ6ff0bXStoEZjZ8uDxQmavq/l+0cZURaWJXkTcwKvAEKATcPepRB7ifVXtrKpdgWeBF0KObVXVroHHmJoK3JjyHDxWwLyQLpWjeqfUeQyJ4mWke36w/J/vdtZ5DLUp0pp740113v9wruh7AltUdZuqFgFTgOFlXvhYSLEhYL8JxhGfZu7l1KDUqy9uSmrzho7EMTJhLi78TTbfbslhW/YJR+KoafXr1ycnJ8eSvUNUlZycHOrXr1+l88LpddMa2B1SzgJ6la0kIg8BvwOSgNBb1qkikgEcAx5T1YXlnDsaGA3Qrl27sIM3pqyPVmYFt2/vXrdt86FaSw4DXBnM9vUA4IP03Tw6JPpvYrZp04asrCyys7OdDiVu1a9fnzZtqva7HU6iL28Ewmkf56r6KvCqiNwDPAaMAvYB7VQ1R0R6ANNE5PIy3wBQ1TeAN8Df66ZKP4ExAev2HmPDfn//5HoJLoZccYGj8Yx0zwsm+s8y9/KHmzpG/UjZxMTEWhm5aWpXOIk+C2gbUm4DVDSRxxTg7wCqWggUBrZXiMhWoANg/SdNjfsko+Rq/ibfQpKfGelgNHC9axWNz0nkSH4xe48WsHxHLr0ubuZoTCY+hdNGvxxoLyKpIpIEjASmh1YQkfYhxZuBzYH9LQI3cxGRi4H2wDaMqWEer49pIXPO3+4+rYWwziWJl5s7Xxgsh8ZnTF2qNNGrqgcYC3wJrAc+VNW1IjJRRIYFqo0VkbUikom/nX5UYP/1wGoRWQVMBcaoam6N/xQm7i3acojs44UAtOAwfVzfOxyR323dWge3Z67ZR5HH+tSbuhfWFAiqOhOYWWbfEyHbvz7DeR8BH5V3zJia9PHKkumHh7sXkyCRkVB7tGtC68YN2HPkJEdPFrNgUzY3djrf6bBMnLGRsSbqHS8o5qt1JTNTRkKzzSkulzCsa8kSdp9m1s58+MZUxBK9iXpfrj1AQbH/Cr6j7KSTa5fDEYWY0IjbFv8oWJy9egcnxtsVvalbluhN1Ju1Zl9we7i7btfiDMelriw6iv/Dp4B6fOm7yuGITLyxRG+i2vGCYhZuPhQsD3V952A0Zzbc/W1we5r3WgcjMfHIEr2JanM3HKQoMDvk5a3O4yLXwUrOcMawkG8ai32XczivyMFoTLyxRG+i2syQZhunR8JWpLXk0E3868l6cTN7vc1oaeqOJXoTtfIKPczfWDLnypCQwUmRaLB7eXD7y7W1uH6tMWVYojdRa/7GbAoDA5AuPT+ZS1qc63BEFbvJVZLov9l8iLzAUofG1DZL9CZqzfw+pNmmc+Q225yS4jpAR/HPTV/k8ZX6NmJMbbJEb6JSQbGXeRtKbrwOjfBmm1NCm2++sOYbU0cs0ZvoM6ERCyYOIr/IC8DFspf2f2sDExo5HFjlQptv5q4/QEGx18FoTLywRG+i0ixvz+D2UNd3SJRM895RdnOR+K/k84q8LN56qJIzjDl7luhN1ClWN3N83YLlwe5lDkZTNSIwOOSq/ovvrfnG1D5L9CbqLPddynH8a8G2JpvLJboW374ppJ1+9vqDeLyRMdOmiV2W6E3Ume3rHty+wZ0RNc02p3SVrZx/Xj0AcvOKWLbDlmgwtcsSvYkqqsqckEQ/0LXCwWiqxyXKoE4l3UHnrI/MaRtM7LBEb6LKloMn2Kn+JNmQk/RyrXc4ouoZGLL4yJz1B1BVB6MxsS6sRC8ig0Vko4hsEZFx5RwfIyJrRCRTRBaJSKeQY48GztsoIjfVZPAm/swOufrt61pNPYnO0aW9UptyTpIbgB05+Ww7lOdwRCaWVZroA4t7vwoMAToBd4cm8oD3VbWzqnYFngVeCJzbCf9i4pcDg4G/nVos3JjqCJ0M7Ab3SgcjOTv1E930+UHzYHmuNd+YWhTOFX1PYIuqblPVImAKMDy0gqoeCyk2BE59Dx0OTFHVQlXdDmwJPJ8xVZZzopCVuw4D4MJHf1emwxGdhQmNGLhpYrA4e9ZHUTHgy0SncBJ9a2B3SDkrsK8UEXlIRLbiv6L/VRXPHS0i6SKSnp1t83+Y8s3dcJBTTdk9ZBNN5bizAZ2lfu6SD6p0vZSj2tDBaEwsCyfRl9d57bQ7R6r6qqpeAvwBeKyK576hqmmqmtaiRYswQjLxKLR3ysAobrY5paUc5UrZCvjnqF/g6+JwRCZWhZPos4C2IeU2wN4K6k8BbqvmucaUq6DYyzebS77t3eCK/kQPMCDkA2uOt3sFNY2pvnAS/XKgvYikikgS/pur00MriEj7kOLNwObA9nRgpIjUE5FUoD0QPePVTcT4bntucBKzVNnHJRIb1wuhH1jzfVfaKFlTKxIqq6CqHhEZC3wJuIG3VHWtiEwE0lV1OjBWRAYCxcBhYFTg3LUi8iGwDvAAD6mqTddnqix0SuL+rugbDXsml8tOzieXAzTlKOeyctcReqY2dTosE2MqTfQAqjoTmFlm3xMh27+u4NyngKeqG6AxAPM3liT6AdHc26YMERjgzmCy9wbAP3jKEr2paTYy1kS87Yfy2JGTD8A5SW6ucm1wOKKaNTCk+WbOButPb2qeJXoT8UKv5ntf0jxqR8OeSW/XWpIoAvxTPGQdznc4IhNrLNGbiDcvZG3V/h1jr/ttAyni6pA5e2wtWVPTLNGbiJZf5GHptpxgud+lLR2Mpvb0c60KbluiNzXNEr2JaEu25lDk8Xc57HD+ubRu3MDhiGpHv5AbzIu3HqLQY53TTM2xRG8i14RGzPvPX4LF/ocmx+x8MKmyn3bin7Atv8hL+o7DDkdkYoklehOxVP2DiE7pF0PdKssSKdt8Y71vTM2xRG8i1lZtRZb62+TPJZ801yaHI6pdoR9k1k5vapIlehOx5vm6Brf7uL4nUWK73foa1zqSEvx/kpsPnmDPkZMOR2RihSV6E7HmhyT6qJ57PkwNpIheIaNirfnG1BRL9CYi5RV6WO67NFju615VQe3Y0T+k+6g135iaYoneRKSl23IoIhGAjrKLCyQ+eqH0+2pwcHvxuh0UjW8asz2NTN2xRG8i0oJNJVezfV3xcTUPpbtZ5tGA9JBvNcZUlyV6E5HiNdGX7Wa5IKR7qTHVZYneRJwdh/LYeWq2SgroEePdKsvqWyrR2/KC5uxZojcRJ/RqvrdrbczNVlmZq13rSMT/M2/QizigjR2OyES7sBK9iAwWkY0iskVExpVz/Hcisk5EVovIHBG5KOSYV0QyA4/pZc81pqxv4rTZ5pSGUlhqzv0FXmu+MWen0kQvIm7gVWAI0Am4W0Q6lamWAaSpahdgKvBsyLGTqto18BhWQ3GbGFXo8bJ4a8lslX1dqx2MxjmhP7c135izFc4VfU9gi6puU9UiYAowPLSCqs5T1VOrJSwF2tRsmCZepO84zMli/wjYFNlPO1d8DhoK/SazyNcZr08djMZEu3ASfWtgd0g5K7DvTB4AZoWU64tIuogsFZHbyjtBREYH6qRnZ9sgkXgWr71tyrpUdnM+uQAc5VxWZR1xOCITzcJJ9FLOvnIvL0TkJ0Aa8FzI7naqmgbcA/xVRC457clU31DVNFVNa9Ei9lYQMuFbsNESPfi7WV7vDmm+sVGy5iyEk+izgLYh5TbA3rKVRGQg8CdgmKoWntqvqnsD/24D5gPdziJeE8P2HT3JxgPHAUhyu0otrxePQj/ovtlsid5UXziJfjnQXkRSRSQJGAmU6j0jIt2A1/En+YMh+5uISL3AdnPgWmBdTQVvYktob5urUptwjhRWUDv29XF9jwv/6lqrdh/hcF6RwxGZaFVpoldVDzAW+BJYD3yoqmtFZKKInOpF8xxwLvDfMt0oLwPSRWQVMA94RlUt0ZtyfbPpUHC7X4fYXBu2KhpLHlfKVgB8Cou2HKrkDGPKlxBOJVWdCcwss++JkO2BZzhvMdD5bAI08cHj9bEwpHmi76UtYK6DAUWIvu5VZHjaA/4b1bde2crhiEw0spGxxnkTGrHq/+3JsQL/aNALyKH936yHLpRpp9+Ujap1szRVZ4neRITQ0Z993auR8vp6xaEuso3G+G9QHzxeyIb9xx2OyEQjS/QmIoSO/oznbpVluUW5zrUmWA4dZ2BMuCzRG8flajKr9WIAXPi41rXW4Ygiy/Uu609vzo4leuO4hb4r0MCvYjfZTCPJcziiyNI3ZOBU+s5cThTG12ye5uxZojeO+8Yb0mzjjs9JzCrSUo5w2YXnAVDsVZaETPpmTDgs0RtH+XxaahWl6+N0tsrK9O1QMjXIgk3xOdGbqT5L9MZR6/cf4xD+hTWacJzOss3hiCJT6URv3SxN1ViiN44KHQ17nWs1brEEVp4eFzWhYZIbgN25J9mRk1/JGcaUsERvHDV/Y0kzxPXWPn9GSQkurrmkebC8YKM135jwWaI3jjleUMyKnYeDZes/X4EJjei7+S/B4oLP34MJjRwMyEQTS/TGMd9uOYQnsHLSFbKdFnLM4YgiW+jygkt8nSjQRAejMdHEEr1xzHxbZKRK2rkOcrH4l4IooB7LfB0djshEC0v0xhGqWmo4fz+3JfpwhH4gzg/plmpMRSzRG0dsOnCCfUcLADiPPLrJZocjig79SiX6rg5GYqKJJXrjiNDeNte51pAgPgejiR69XOupj3/lrW3ail3WzdKEwRK9cYS1z1dPfSmmd8ikb/NtlKwJQ1iJXkQGi8hGEdkiIuPKOf47EVknIqtFZI6IXBRybJSIbA48RtVk8CY6nSj0kL4zN1jua+3zVVKq+cZmszRhqDTRi4gbeBUYAnQC7haRTmWqZQBpqtoFmAo8Gzi3KTAe6AX0BMaLSJOaC99Eo8VbDlHs9XervOzC8zhfjjgcUXQJTfSLtx6ioNjrYDQmGoRzRd8T2KKq21S1CJgCDA+toKrzVPVUY+FS4NQ6cDcBX6tqrqoeBr4GBtdM6CZazQ/tbXNpiwpqmvKU6mZZ7OO77bmVnGHiXTiJvjWwO6ScFdh3Jg8As6pyroiMFpF0EUnPzravojFrQiN0fCMWfLciuKvvYmvNq45+rszg9nybDsFUIpxEX97qneXOPCUiPwHSgOeqcq6qvqGqaaqa1qKFXeHFss3amj34/4+TyaeHdausFmunN1URTqLPAtqGlNsAe8tWEpGBwJ+AYapaWJVzTfyY6+sW3L7OtYZEsfbl6ujp2kAD/OMQth/KY2eOrcplziycRL8caC8iqSKSBIwEpodWEJFuwOv4k3zo98gvgUEi0iRwE3ZQYJ+JU3O9JYm+vyvDwUiim7+b5bpg2a7qTUUqTfSq6gHG4k/Q64EPVXWtiEwUkWGBas8B5wL/FZFMEZkeODcXeBL/h8VyYGJgn4lDR7UhK7RDsNzPnVlBbVOZ0Hb6ORusnd6cWUI4lVR1JjCzzL4nQrYHVnDuW8Bb1Q3QxI4Fvi548S+ecaVssdkqz1J/dyYE1glfujWHvEIPDeuF9Sdt4oyNjDV1JrTZZoDbmm3OVhs5RMcLkgEo8vr4dsuhSs4w8coSvakT3jKLgA9wWbNNTRjQsWVwe64135gzsERv6kTm7sMcxn/12YLDXC47nA0oRtxwWelE7/PZmrvmdJboTZ0Ivdoc4M7EZYuA14iubZvQ5Bz/SlMHjxeydq/d9zCns0Rv6sSc9SWJ3rpV1hy3S+h/aclV/ZwNBxyMxkQqS/Sm1u09cpIN+48DkIiHPq7vHY4otgy4zNrpTcUs0ZtaNy9kLpZervWcKwUORhN7ru/QggSXf7aR1VlHOXDM3l9TmiV6U+u+XlfSnDDAmm1q3Hn1E+mZ2jRYnmdX9aYMS/SmVp0o9LB4S06wfKMr3cFoYldoN0sbJWvKskRvatU3m7Ip8vrXg+0oO2nrskE9NWpCI5jQiBu+HhLctWjdTluMxJRiid7Uqq/W7g9uD3KtqKCmORuprv3BxUhOUt9GyZpSLNGbWlPs9ZXqBXKj2xJ9bQptFvtqrXWzNCUs0Ztas3x7LscK/LNuXdioPlfIdocjim03uUsS/ez1B/DaKFkTYIne1JqvQnrb3NjpfKS89cZMjekqW2nJYQBy8opYsfOwwxGZSGGJ3tQKVS3VrfLGTuc7GE18cImWah77MuT+iIlvluhNrVi/7zh7jpwEILleAr1SmzkcUXwYFNpOv24/qtZ8YyzRm1oSejXfr2NLkhLsV60uXONaSzL5AOzOLZl6wsS3sP76RGSwiGwUkS0iMq6c49eLyEoR8YjIiDLHvIHlBYNLDJrY9/X6kmYDa7apO0niLTVpnDXfGAgj0YuIG3gVGAJ0Au4WkU5lqu0C7gPeL+cpTqpq18BjWDnHTYzZ/cQP+H6Pf7rcRDz0+7i7f2CPqROD3NbN0hhm5esAABVDSURBVJQWzhV9T2CLqm5T1SJgCjA8tIKq7lDV1YCvFmI0UWamr1dw+zrXas6Tkw5GE3/6uVaR5Pb/aa/bd4zdufkOR2ScFk6ibw3sDilnBfaFq76IpIvIUhG5rbwKIjI6UCc9Ozu7Ck9tItFMb8/g9hDXMgcjiU/nSgHX/qDk5rc135hwEn15vZ+rciu/naqmAfcAfxWRS057MtU3VDVNVdNatGhRhac2kSbrcD6r9AcAJOCx0bAOuenyC4LbM9fsczASEwnCSfRZQNuQchtgb7gvoKp7A/9uA+YD3aoQn4kyX3xfcvV4ret7Gkueg9HEr5suvwB3YI76lbuOBLu6mvgUTqJfDrQXkVQRSQJGAmH1nhGRJiJSL7DdHLgWWFfdYE3k+zzk6nGoNds4pknDJK79QfNgeZZd1ce1ShO9qnqAscCXwHrgQ1VdKyITRWQYgIhcJSJZwB3A6yKyNnD6ZUC6iKwC5gHPqKol+hi198hJMnYdAcCNt1TvD1PHJjTilm1/DhZnzJwenNLYxJ+EcCqp6kxgZpl9T4RsL8ffpFP2vMVA57OM0USJ0Gab3q61NJETDkZjBrnT+aPnATwkkKnt2e1rbusBxCkbrmhqzMxSzTbfORiJAWgseVznWhMszwrp9mriiyV6UyP2Hy0gPTBbojXbRI6b3UuD2597r3YwEuMkS/SmRoRezV/tWk8zsTlWIsGNrhUk4l8TYJVewm6fdV+OR5boTY2YlrknuH2La4mDkZhQjSSf612rg+UZPruqj0eW6M1Z25p9gtVZRwFIcrsY6rZulZEktPlmhjXfxCVL9OasfZpRcjXfv2MLGtkgqYgy0LWCJIoAWKupbDlozWrxxhK9OSuqyrTMkoHSt3WtyjRIpi6cJycZGDJ18ccr91RQ28QiS/TmrKzcdYRdgdkRk+sn0L9jS4cjMuW53b0wuP1Jxh58tnB4XLFEb87KpyE3YYdecSH1E90ORmPOpK9rFU3xrxGw72gBS7flOByRqUuW6E21FY9vyowlJQNybls12obYR6hE8TLMvThY/jjDmm/iiSV6U20LfZ3J5TwALiSHXq4NDkdkKvJD96Lg9qw1+8gv8jgYjalLluhNtX3i7RPcHub+FpdYu28k6yLbuET8V/J5RV5bZjCOWKI31XIkv4gvfWnB8m3ubx2MxoRDpPRNWWu+iR+W6E21TMvYQxFJAHSWbVzm2l3JGSYShH4gL9qczYFjBQ5GY+qKJXpTZarKlOUliX2ke66D0ZiqaC05XHOxfz1Zn8JHK7McjsjUBUv0psoydx9hw37/6MoGFDDMbXPbRJMRPUqWjpiybLf1qY8DYSV6ERksIhtFZIuIjCvn+PUislJEPCIyosyxUSKyOfAYVVOBG+d8EHI1f4t7Kcli65FGk5u7XEijBokA7MrNZ9EWW4wk1lWa6EXEDbwKDAE6AXeLSKcy1XYB9wHvlzm3KTAe6AX0BMaLSJOzD9s45UShh+mrSqY8GOme52A0pjrqP9WUHxV9Giy/984rNv4hxoVzRd8T2KKq21S1CJgCDA+toKo7VHU14Ctz7k3A16qaq6qHga+BwTUQt3HIjFV7yS/yAtBesugumx2OyFTHPSH3VWb7enBAGzsYjalt4ST61kBol4qswL5whHWuiIwWkXQRSc/Ozg7zqY0TJoc029zlnoeIg8GYavuBay89ZT0AXtx86O3nbECmVoWT6Mv7Uw737k1Y56rqG6qapqppLVrYCjiRat3eY6zafQTwzzt/e8hISxN9fpwwO7g92TMAr92UjVnhJPosoG1IuQ2w9wx1a/JcE2He/nZ7cPumKy6gqS0XGNUGu5YHJzrbS3MWbDrocESmtoST6JcD7UUkVUSSgJHA9DCf/0tgkIg0CdyEHRTYZ6LMoROFfBoy7/x9vVOcC8bUiHriYYT7m2D5vaW7HIzG1KZKE72qeoCx+BP0euBDVV0rIhNFZBiAiFwlIlnAHcDrIrI2cG4u8CT+D4vlwMTAPhNl3lu6iyKv/177lW0b072d3byLBXeH3JSdu/EgW7NPOBiNqS2iGlntcmlpaZqenu50GCZE4fhmXFs4iUP4k/tLiS8z3AZJxYwHin7PHF93AO7p1Y6nf9jZ4YhMdYjIClVNK++YjYw1lZrhuyaY5M8nl6EuW/w7ljzo/jy4/dGKLHJOFDoYjakNluhNhVSVtzwlQx/uTfiKRPE6GJGpaVe71tNZtgFQ6PHx7tKdDkdkapolelOhZdtzWaupANSjqNRAGxMbRODBhJKr+neX7KSg2D7MY4klelOhfyzcFty+3b2QJmI362LRUNcyWjWqD0BOXhEfr7S56mOJJXpzRt/vOcrs9SV9q+93f+FgNKY2JYqX+/ukBsv/XLTNZrWMIZbozRm9PLdkHpuhru9o77KrvFh211VtSa6XAMC27Dy+Wrff4YhMTbFEb8q1ft8xvgxZU/SXCR87GI2pC8nPNOceT8n/84vvTcM33sZLxAJL9KZcoVfzN7mW2VKBcWJ0wuecg395wY3ajs99vRyOyNQES/TmNBv3H2fmmpKv7b9M+MTBaExdaibH+VnIvZi/en5kk53FAEv05jSvzNsS3B542flc4bJ+1fHk5wmfk0w+AFu1NZ9m2r2ZaGeJ3pSY0Ih1T3RmxqqSBaN/tfVBBwMyTmgsedzvnhUsvzRnMx5v2TWFTDSxRG+CVOHPnp+ggV+LAa6VdHFtr+QsE4seSJhJI/xjJnbm5Fu/+ihnid4EzfV1Y7HvCgDceHk0YbLDERmnnCcnGR0yWvb5rzZyotDjYETmbFiiNwAUe3085flxsHyPe471m49zo9xf0oLDABw8XsirIfduTHSxRG8AeG/pTrZpKwCSyeM3CR85HJFx2rlSwKOJJd/q3ly4ne2H8hyMyFSXJXrD0fxi/jqnpN/8LxOm0cyWCTTAba5vg4vMFHl9/HnGOocjMtVhid7w3FcbOJJfDEA7OcAot632aPxcokzY/xCCv9fNnA0Hmfd4X5jQyOHITFWElehFZLCIbBSRLSIyrpzj9UTkg8Dx70QkJbA/RUROikhm4PFazYZvztbSbTn8J2St0EcT3qee2E03U6KLazt3uecHyxM9P6VAE50LyFRZpYleRNzAq8AQoBNwt4h0KlPtAeCwqv4AeBH4v5BjW1W1a+AxpobiNjXgZJGXcR+tDpYHulYw2LXcwYhMpPp9wock42+f364X8qJnhMMRmaoI54q+J7BFVbepahEwBRheps5w4F+B7anADSIiNRemqQ0vzt7Ejhz/CMjkegn8OfEt7H/NlKe5HOMPCVOC5Te8N5O+I9fBiExVhJPoWwOhM1plBfaVW0dVPcBRoFngWKqIZIjIAhG5rrwXEJHRIpIuIunZ2dlV+gFM9WTuPsI/QxYVeeyWy7hADjsYkYl0P3bP4TqX/xug4uLh/64iz/rWR4VwEn1513hlZzk6U519QDtV7Qb8DnhfRM47raLqG6qapqppLVq0CCMkczbyx7fk93/7gFNzVfVxreHOGZ2dDcpEPBH4v8Q3gk04O3PyeWbWBoejMuEIJ9FnAW1Dym2AvWeqIyIJQCMgV1ULVTUHQFVXAFuBDmcbtKk+VeWx4vvZom0AaEABf0n4pzXZmLC0klwmJP47WH536U4WbLJv4ZEunES/HGgvIqkikgSMBKaXqTMdGBXYHgHMVVUVkRaBm7mIyMVAe2AbxjGTl+3mY19JC9qEhH/T1mV/qCZ8t7sWcqMrPVj+9ZQMdufmOxiRqUyliT7Q5j4W+BJYD3yoqmtFZKKIDAtUexNoJiJb8DfRnOqCeT2wWkRW4b9JO0ZV7Q6OQ9ZkHWXC9LXB8h3u+dyVMN+5gExUEoGnE/9Jy+R6ABzJL+Z/3l3BySKvw5GZMxHVyFpUIC0tTdPT0yuvaKrkaH4xt7yykN25JwHoKDv5JGk8DaTI4chMtFrxsx2MfGMJxV5/Dhl2ZSteGtkV63DnDBFZoapp5R2zkbFxoKDYy4P/Xh5M8ueSz98TX7Ikb85Kj7dTmChvBMvTV+3lH4//uIIzjFMs0cc4j9fH2PdXsnxHSdfJ5xJfJ9W1v4KzjAnP3QnzuMc9O1j+i+duPsnIquAM4wRL9DFMVfnjJ2uYvf5gcN9jN1/GELeNfjU1Z0LCv0iTjUCgf/2Hq5i1Zp/DUZlQluhjlKry9GO/4MP0kqurMe7pPDinm4NRmViUJF7+mfQ8HcW/trBP4VdTMpi38WAlZ5q6Yok+Bnm8PsZ9tIZ/eG8J7hvhXlBqCLsxNamx5PFu0l+4WPxDbIq9yph3VzB3wwGHIzNgiT7mFBR7+cV7K/kgvWTWihtd6TyT8A8bFGVqVQs5xntJT9OmSQMACj0+HvxXOu8u3elwZMYSfQzJOVHIqLeW8dW6kquo213f8LfEl0gQn4ORmXhxoeTyft7/0Eb8zTY+hcenfc9f/vQ/+HyR1ZU7nliijxHLd+Ry86RFfLe9ZDzaz90zeD7xdRLFBrKYutPOdZCPk8bTRbYG973uHcb//GcFh/OsS68TLNFHOZ9PeW3BVka+sZT9xwqC+8cN6cifEt/HJXYVZepeSznKlKQ/M9C1Irjv63UHGPLSQhZvOeRgZPHJEn0U23LwOPf8cynPzNqAN/C1uDHHeTvxWcbM6+5wdCbenSOFvJ74Ave7Zwb37T9WwI/f/I6/zFpvUybUIZsCIQqdLPLy8tzN/GPhtuDwc4DusomXk16mteQ4GJ0xp/va250/JP2R3JCmm1aN6vPo0Mu4pcuFNm1CDahoCgRL9FGkoNjLf9N387f5W9l3tKSZxo2Xn7s/5+GE/1p7vIlYB7UxDxePYaGvS6n9PVOa8sjgS0lLaepQZLHBEn2UO1Ho4b9P3sNrnls5QOk/hh6ykT8nvsVlrt1nONuYyOFT4UNvX57z3EUOjUodS7uoCWP6XsKAji1xuewKv6os0UchVSVj9xGmLNvFjNX7yC/Tntmco/w+4QPudC+wG64m6hzTBkzy3M473pvwkFDqWGrzhvyoe2t+2L0NrRs3cCjC6GOJPkp4fcrKXYf5et0Bvlq7P7hwd6gWHGZMwmfc455rs0+aqLfNdwF/7/wh0zL3lLrfdMrVFzdlUKcL6N+xJanNGzoQYfSwRB+hijw+Nuw/xrLtuSzbnsvyHbkczi8ut257yeKn7q+50z2f+lJ+HWOi1T5tyluewUz2DuAE55RbJ6XZOVxzSXN6XNSEHhc1IaXZOXYTN4QleoedLPKyMzePHYfy2ZGTx6b9x1m37xhbs0+UexVzyjkUcLN7KSPd8+gum20KAxPzTmoSX/nS+Mh7HYv0SioaTNv4nEQ6nJ/Mpecn0+GCZFKbNaRt0wa0atyARHf89Rw/60QvIoOBlwA38E9VfabM8XrAv4EeQA5wl6ruCBx7FHgA8AK/UtUvK3qtSE30Xp9S6PFysshLfpGXk8Ve8go9HC/wcKLQw/GCYg7nF5ObV0RuXhHZxws5cKyA/ccKOHKGq/TytOAwA90rGeRawTWutXb1buLWAW3MHG935vq68q3vCk5SP6zzXAIXnFefFufVp2VyPVok16PpOUk0PieRJuckcV6DRBrWc5Ncz/9vgyQ3DRLd1E90Uy/BFbXfEs4q0QcW994E3Ahk4V8s/G5VXRdS5xdAF1UdIyIjgR+q6l0i0gmYDPQEWgGzgQ6qesY+gNVN9OM+Ws2eIydRBcX/M6mCTzWwDwiUvar41D+q1KeK11fy8PgUj9dHkVfx+HwUe3wUenx4ammejrZykB6yiatcG+jl2sAlsteu3I0po1ATWOHrwEptzwpfBzJ8P+AIybXyWkluF0kJ/keCS0h0u0hwC26X4JbAv4GHiOAWcIngEgHxf9AIgoh/fd1T28BpHyKnSqd2/3JAe3pc1KRacVeU6BPK21lGT2CLqm4LPNkUYDiwLqTOcGBCYHsq8Ir4f6LhwBRVLQS2BxYP7wksqc4PUpH0nYfZcvBETT9tjUjAQ1vJ5iI5QIrsJ1X2c5lrJx1lF+fJSafDMybi1RMPvd3r6B1IO6qwh+Zs8rVho7Zls681u7Ulu7UF+2l2Vq9V5PVR5PVBYU1EXjX39GxXK88bTqJvDYR20s4Cep2pjqp6ROQo0Cywf2mZc1uXfQERGQ2MDhRPiASWq6ldzYE6m3Rja+VVnFan70eUsPfkdBH0nhwDtjkdBNTge3LT/53V6Red6UA4ib68hoSy7RhnqhPOuajqG8Ab5dStNSKSfqavOfHI3o/T2XtyOntPThcN70k4t6azgLYh5TbA3jPVEZEEoBGQG+a5xhhjalE4iX450F5EUkUkCRgJTC9TZzowKrA9Apir/ru804GRIlJPRFKB9sCymgndGGNMOCptugm0uY8FvsTfvfItVV0rIhOBdFWdDrwJvBu42ZqL/8OAQL0P8d+49QAPVdTjpo7VaVNRFLD343T2npzO3pPTRfx7EnEDpowxxtSs+Bs+ZowxccYSvTHGxLi4T/Qi8nsRURFp7nQsThOR50Rkg4isFpFPRKSx0zE5RUQGi8hGEdkiIuOcjsdpItJWROaJyHoRWSsiv3Y6pkggIm4RyRCRGU7HUpG4TvQi0hb/1A67nI4lQnwNXKGqXfBPe/Gow/E4IjDtx6vAEKATcHdgOo945gEeVtXLgKuBh+w9AeDXwHqng6hMXCd64EXgEcoZxBWPVPUrVfUEikvxj3uIR8FpP1S1CDg17UfcUtV9qroysH0cf3I7bZR7PBGRNsDNwD+djqUycZvoRWQYsEdVVzkdS4S6H5jldBAOKW/aj7hOaqFEJAXoBnznbCSO+yv+C0Wf04FUJpwpEKKWiMwGLijn0J+APwKD6jYi51X0nqjqp4E6f8L/Vf29uowtgoQ1dUc8EpFzgY+A36jqMafjcYqI3AIcVNUVItLP6XgqE9OJXlUHlrdfRDoDqcCqwLShbYCVItJTVffXYYh17kzvySkiMgq4BbhB43eQhU3dUQ4RScSf5N9T1Y+djsdh1wLDRGQoUB84T0T+o6o/cTiuctmAKUBEdgBpqhohs/I5I7DAzAtAX1XNdjoepwTma9oE3ADswT8NyD2qutbRwBwUmHb8X0Cuqv7G6XgiSeCK/veqeovTsZxJ3LbRm3K9AiQDX4tIpoi85nRATgjckD417cd64MN4TvIB1wI/BQYEfjcyA1ezJgrYFb0xxsQ4u6I3xpgYZ4neGGNinCV6Y4yJcZbojTEmxlmiN8aYGGeJ3hhjYpwlemOMiXH/PxecMgfwNVbIAAAAAElFTkSuQmCC\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"# Plot a histogram of the observed data\n",
"# Included is expected distribution, if the data is normally distributed, with the same mean and std of the data. \n",
"xs = np.arange(d1.min(), d1.max(), 0.1)\n",
"fit = stats.norm.pdf(xs, np.mean(d1), np.std(d1))\n",
"plt.plot(xs, fit, label='Normal Dist.', lw=3)\n",
"plt.hist(d1, 50, density=True, label='Actual Data');\n",
"plt.legend();"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The plot above shows a histogram (in orange) of 'observed' data, and the expected distribution of the data, under the assumption of it being normally distributed. There seems to be a pretty good alignment! Now, let's asses this statistically. \n",
"\n",
"To do so, we will use the `normaltest` function, from scipy, which tests whether a sample of data differs from a what would be expected under a normal distribution. The null hypothesis is that the data come from a normal distribution. We can use `normaltest` to check this null, and to decide whether we should reject the null (to claim the data are not normal), or accept it. To assess this, the `normaltest` function compares the skew and kurtosis of the observed data to what would be expected if it's normally distributed. "
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"# Check the documentation for normaltest\n",
"normaltest?"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"# Run normal test on the data\n",
"stat, p_val = normaltest(d1)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"Normaltest p-value is: 0.49 \n",
"\n",
"We do not have evidence to reject the null hypothesis.\n"
]
}
],
"source": [
"# Check the p-value of the normaltest\n",
"print('\\nNormaltest p-value is: {:1.2f} \\n'.format(p_val))\n",
"\n",
"# With alpha value of 0.05, how should we proceed\n",
"check_p_val(p_val, alpha=0.05)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this case, we do not reject the null, meaning we do not reject that this data is normally distributed. \n",
"\n",
"This statistical test is therefore consistent with the data being normally distributed. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Testing a non-Normal distribution\n",
"\n",
"Let's test another dataset for normality.\n",
"\n",
"For this example, we will explore case when visual inspection can be misleading, to emphasize the point that although visually inspecting datasets can give you a quick sense about whether they are normally distributed or not, this is not enough. Visual inspection can be somewhat misleading, since non-normally data can 'look normal'.\n",
"\n",
"This is why, as well as visually checking data when it looks pretty normal, it is important to perform other checks.\n",
"\n",
"For example, under some parameters, the Beta distributed data can look like it is normally distributed. "
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"# Generate some data from a beta distribution\n",
"d2 = stats.beta.rvs(7, 10, size=100000)"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD6CAYAAACxrrxPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAVOklEQVR4nO3df5CV1Z3n8fd3BAddmWSRTsYApsmOrklMAqRXHcZsGeNYDGUkG3ElkxnF0qFMqRlrk3KNs1FiqtyJGq0xpDRELczUBExhomDJToVRyh8VNA0iBhknTCRlb6jYgkFZRCX57h/3Ytr2dvdz4Xb35fB+Vd3y+XH6uV+fxg/Hc89zbmQmkqSD3x+MdgGSpNYw0CWpEAa6JBXCQJekQhjoklQIA12SCjFkoEfEuIh4MiKejohNEfG1Bm3mR0RvRGyovy4ennIlSQMZU6HN68DpmbkrIsYCj0XEqsxc26/dPZl5WdU3njhxYnZ2djZRqiRp3bp1L2VmR6NzQwZ61p482lXfHVt/HfDTSJ2dnXR3dx/oZSTpkBIRvxzoXKUx9Ig4LCI2AC8CP87MJxo0OyciNkbE8oiYsp+1SpL2U6VAz8zfZuY0YDJwUkSc2K/JSqAzMz8KrAbubnSdiFgQEd0R0d3b23sgdUuS+mlqlktm/gZYA8zqd3x7Zr5e3/0u8PEBfn5xZnZlZldHR8MhIEnSfhpyDD0iOoA3M/M3EXEEcAbwjX5tjsnMbfXds4HNLa9U0qh488036enpYc+ePaNdyiFl3LhxTJ48mbFjx1b+mSqzXI4B7o6Iw6j16H+QmQ9ExHVAd2auAL4YEWcDe4EdwPymq5fUlnp6ehg/fjydnZ1ExGiXc0jITLZv305PTw9Tp06t/HNVZrlsBKY3OH5Nn+2vAF+p/K6SDhp79uwxzEdYRHD00UfT7GeNPikqaUiG+cjbn3tuoEtSIaqMoUvS7y18V4uvt7NSsx/96Ed89rOfZfPmzZxwwgmDtl2yZAlnnnkm73vf+/arpDVr1nDTTTfxwAMPvOP4nDlz+MAHPsDu3bt573vfy5VXXslZZ5015PUOP/xwZs6cuV/1VGWgq30NFBwVA0BlWbp0KaeeeirLli1j4cKFg7ZdsmQJJ5544n4H+mA+8YlPvBX0GzZs4DOf+QxHHHEEn/rUpwb8mTVr1nDUUUcNe6A75CKp7e3atYvHH3+cO++8k2XLlr3t3A033MBHPvIRPvaxj3HVVVexfPlyuru7+fznP8+0adN47bXX6Ozs5KWXXgKgu7ub0047DYAnn3ySmTNnMn36dGbOnMlzzz3XVF3Tpk3jmmuuYdGiRQCsXLmSk08+menTp3PGGWfw61//mq1bt3L77bdzyy23MG3aNB599NGG7VrBHroOPvbcDzn33Xcfs2bN4vjjj2fChAmsX7+eGTNmsGrVKu677z6eeOIJjjzySHbs2MGECRNYtGgRN910E11dXYNe94QTTuCRRx5hzJgxrF69mquvvpp77723qdpmzJjBjTfeCMCpp57K2rVriQjuuOMObrjhBr75zW9yySWXcNRRR/HlL38ZgJdffrlhuwNloEtqe0uXLuWKK64AYN68eSxdupQZM2awevVqLrzwQo488kgAJkyY0NR1d+7cyQUXXMDPf/5zIoI333yz6dpq6xfW9PT0cN5557Ft2zbeeOONAeeQV23XLIdcJLW17du389BDD3HxxRfT2dnJjTfeyD333ENmkpmVpveNGTOG3/3udwBve+L1q1/9Kp/85Cf52c9+xsqVK/fradinnnqKD37wgwBcfvnlXHbZZTzzzDN85zvfGfB6Vds1y0CX1NaWL1/O+eefzy9/+Uu2bt3KCy+8wNSpU3nsscc488wzueuuu9i9ezcAO3bsAGD8+PG8+uqrb12js7OTdevWAbxtSGXnzp1MmjQJqH2Q2qyNGzfy9a9/nUsvvfQd17v77t+vUdi/noHaHSiHXCQ1Z4Q/q1i6dClXXXXV246dc845fP/73+e2225jw4YNdHV1cfjhhzN79myuv/565s+fzyWXXMIRRxzBT37yE6699louuugirr/+ek4++eS3rnPllVdywQUXcPPNN3P66adXqufRRx9l+vTp7N69m/e85z3ceuutb81wWbhwIeeeey6TJk3ilFNO4fnnnwfg05/+NHPnzuX+++/nW9/61oDtDlT0Hf8ZSV1dXekXXGhQzc539kPRYbF58+a3hhQ0shrd+4hYl5kNP+21h65yOPtFhzjH0CWpEAa6pCGN1tDsoWx/7rmBLmlQ48aNY/v27Yb6CNq3Hvq4ceOa+jnH0FW+wT5cdXx9SJMnT6anp6fptbl1YPZ9Y1EzDHSNvlav3qeWGjt2bMueZNTwcshFkgphoEtSIQx0SSqEY+g6tPkwkgpiD12SCjFkoEfEuIh4MiKejohNEfG1Bm3+MCLuiYgtEfFERHQOR7GSpIFVGXJ5HTg9M3dFxFjgsYhYlZlr+7S5CHg5M/8kIuYB3wDOG4Z6dTBzeqI0rIbsoWfNrvru2Pqr/yNjc4B9i/ouBz4VVVadlyS1TKUx9Ig4LCI2AC8CP87MJ/o1mQS8AJCZe4GdwNENrrMgIrojotunziSptSoFemb+NjOnAZOBkyLixH5NGvXG37HwQ2YuzsyuzOzq6OhovlpJ0oCamuWSmb8B1gCz+p3qAaYARMQY4F3AjhbUJ0mqqMosl46IeHd9+wjgDOBf+zVbAVxQ354LPJQuzSZJI6rKLJdjgLsj4jBqfwH8IDMfiIjrgO7MXAHcCfxjRGyh1jOfN2wVS5IaGjLQM3MjML3B8Wv6bO8Bzm1tadIo8glSHYR8UlSSCuFaLmo9HyCSRoU9dEkqhIEuSYUw0CWpEAa6JBXCQJekQhjoklQIA12SCmGgS1IhDHRJKoRPikrNcI0XtTEDXfvPR/yltuKQiyQVwkCXpEIY6JJUCANdkgphoEtSIQx0SSqEgS5JhTDQJakQBrokFWLIQI+IKRHxcERsjohNEfG3DdqcFhE7I2JD/XXN8JQrSRpIlUf/9wJfysz1ETEeWBcRP87MZ/u1ezQzz2p9idJBwDVe1AaG7KFn5rbMXF/ffhXYDEwa7sIkSc1pagw9IjqB6cATDU7/aUQ8HRGrIuLDA/z8gojojoju3t7epouVJA2scqBHxFHAvcAVmflKv9Prgfdn5seAbwH3NbpGZi7OzK7M7Oro6NjfmiVJDVQK9IgYSy3M/ykzf9j/fGa+kpm76tsPAmMjYmJLK5UkDarKLJcA7gQ2Z+bNA7T543o7IuKk+nW3t7JQSdLgqsxy+TPgr4FnImJD/djVwLEAmXk7MBf4QkTsBV4D5mVmDkO9Gg1+kYV0UBgy0DPzMSCGaLMIWNSqoiRJzfNJUUkqhIEuSYUw0CWpEAa6JBXCQJekQhjoklSIKvPQJe0vV2HUCLKHLkmFMNAlqRAGuiQVwkCXpEIY6JJUCANdkgrhtEX9nsvkSgc1e+iSVAgDXZIKYaBLUiEMdEkqhIEuSYUw0CWpEAa6JBViyECPiCkR8XBEbI6ITRHxtw3aRETcGhFbImJjRMwYnnIlSQOp8mDRXuBLmbk+IsYD6yLix5n5bJ82fwEcV3+dDNxW/6ckaYQMGeiZuQ3YVt9+NSI2A5OAvoE+B/heZiawNiLeHRHH1H9WUn9+8YWGQVNj6BHRCUwHnuh3ahLwQp/9nvoxSdIIqRzoEXEUcC9wRWa+0v90gx/JBtdYEBHdEdHd29vbXKWSpEFVCvSIGEstzP8pM3/YoEkPMKXP/mTgV/0bZebizOzKzK6Ojo79qVeSNIAqs1wCuBPYnJk3D9BsBXB+fbbLKcBOx88laWRVmeXyZ8BfA89ExIb6sauBYwEy83bgQWA2sAXYDVzY+lIlSYOpMsvlMRqPkfdtk8ClrSpKktQ8v+DiUOQXWUhF8tF/SSqEgS5JhTDQJakQBrokFcJAl6RCGOiSVAgDXZIK4Tx0qZ24rK4OgD10SSqEgS5JhTDQJakQBrokFcJAl6RCGOiSVAgDXZIKYaBLUiEMdEkqhIEuSYUw0CWpEAa6JBXCQJekQgy52mJE3AWcBbyYmSc2OH8acD/wfP3QDzPzulYWqf0w0Kp9kopVZfncJcAi4HuDtHk0M89qSUWSpP0yZKBn5iMR0Tn8pUgakOukq4JWjaH/aUQ8HRGrIuLDLbqmJKkJrfjGovXA+zNzV0TMBu4DjmvUMCIWAAsAjj322Ba8tSRpnwPuoWfmK5m5q779IDA2IiYO0HZxZnZlZldHR8eBvrUkqY8DDvSI+OOIiPr2SfVrbj/Q60qSmlNl2uJS4DRgYkT0ANcCYwEy83ZgLvCFiNgLvAbMy8wctoolSQ1VmeXyuSHOL6I2rVGSNIp8UlSSCmGgS1IhDHRJKoSBLkmFMNAlqRAGuiQVwkCXpEK0Yi0XjSbXPT+0uQqj+rCHLkmFMNAlqRAGuiQVwkCXpEIY6JJUCANdkgphoEtSIQx0SSqEgS5JhfBJ0YOFT4RKGoI9dEkqhIEuSYUw0CWpEI6hSyVyFcZD0pA99Ii4KyJejIifDXA+IuLWiNgSERsjYkbry5QkDaXKkMsSYNYg5/8COK7+WgDcduBlSZKaNWSgZ+YjwI5BmswBvpc1a4F3R8QxrSpQklRNKz4UnQS80Ge/p37sHSJiQUR0R0R3b29vC95akrRPKwI9GhzLRg0zc3FmdmVmV0dHRwveWpK0TysCvQeY0md/MvCrFlxXktSEVgT6CuD8+myXU4CdmbmtBdeVJDVhyHnoEbEUOA2YGBE9wLXAWIDMvB14EJgNbAF2AxcOV7GSpIENGeiZ+bkhzidwacsqkiTtFx/9l6RC+Oi/dChxSYCi2UOXpELYQ283fpGFpP1kD12SCmGgS1IhDHRJKoSBLkmFMNAlqRAGuiQVwkCXpEI4D13S4M8/+BTpQcMeuiQVwkCXpEIY6JJUCANdkgphoEtSIZzlMlpcVVFSi9lDl6RCGOiSVAgDXZIK4Ri6pMH5PaQHjUo99IiYFRHPRcSWiLiqwfn5EdEbERvqr4tbX6okaTBD9tAj4jDg28CfAz3ATyNiRWY+26/pPZl52TDUKEmqoEoP/SRgS2b+IjPfAJYBc4a3LElSs6oE+iTghT77PfVj/Z0TERsjYnlETGl0oYhYEBHdEdHd29u7H+VKkgZSJdCjwbHst78S6MzMjwKrgbsbXSgzF2dmV2Z2dXR0NFepJGlQVQK9B+jb454M/Kpvg8zcnpmv13e/C3y8NeVJkqqqMm3xp8BxETEV+L/APOAv+zaIiGMyc1t992xgc0urPJj5iL9K5XTGtjNkoGfm3oi4DPhn4DDgrszcFBHXAd2ZuQL4YkScDewFdgDzh7FmSVIDlR4syswHgQf7Hbumz/ZXgK+0tjRJUjN89F+SCmGgS1IhDHRJKoSBLkmFcLVFSa3ldMZRY6C3ivPNJY0yh1wkqRAGuiQVwkCXpEI4hi5pZPhh6bCzhy5JhTDQJakQDrk0y+mJktqUgS5pdDm23jIOuUhSIQx0SSqEgS5JhXAMfSB++CmNLsfWm2YPXZIKYQ9d0sHFnvuADHSHViQVolKgR8Qs4B+Aw4A7MvPv+53/Q+B7wMeB7cB5mbm1taVK0iDsuQ8d6BFxGPBt4M+BHuCnEbEiM5/t0+wi4OXM/JOImAd8AzhvOAqWpKYcQkFfpYd+ErAlM38BEBHLgDlA30CfAyysby8HFkVEZGa2sNYD49CKpMJVCfRJwAt99nuAkwdqk5l7I2IncDTwUiuKbIrBLamKVmZFm/T2qwR6NDjWv+ddpQ0RsQBYUN/dFRHPVXj/4TCR0fjLZmjtWhe0b23tWhe0b23tWhe0b22D1/W1RhE4bN4/0Ikqgd4DTOmzPxn41QBteiJiDPAuYEf/C2XmYmBxhfccVhHRnZldo11Hf+1aF7Rvbe1aF7Rvbe1aF7Rvbe1aV39VHiz6KXBcREyNiMOBecCKfm1WABfUt+cCD7XV+LkkHQKG7KHXx8QvA/6Z2rTFuzJzU0RcB3Rn5grgTuAfI2ILtZ75vOEsWpL0TpXmoWfmg8CD/Y5d02d7D3Bua0sbVqM+7DOAdq0L2re2dq0L2re2dq0L2re2dq3rbcKREUkqg4tzSVIhig30iJgVEc9FxJaIuKrB+f8aEesjYm9EzG2z2v5HRDwbERsj4l8iYsBpSiNc1yUR8UxEbIiIxyLiQyNRV5Xa+rSbGxEZESMyI6HCPZsfEb31e7YhIi4eibqq1FZv89/rf9Y2RcT326GuiLilz/36t4j4zUjUVbG2YyPi4Yh4qv7f5+yRqq2SzCzuRe3D238HPgAcDjwNfKhfm07go9TWoJnbZrV9Ejiyvv0F4J42qeuP+myfDfyfdrln9XbjgUeAtUBXO9QFzAcWjdSfryZrOw54CviP9f33tENd/dpfTm0iRrvcs8XAF+rbHwK2jvTvdrBXqT30t5YryMw3gH3LFbwlM7dm5kbgd21Y28OZubu+u5ba3P92qOuVPrv/gQYPj41WbXVfB24A9rRZXaOhSm1/A3w7M18GyMwX26Suvj4HLB2BuqBabQn8UX37XbzzmZxRVWqgN1quYNIo1dJfs7VdBKwa1opqKtUVEZdGxL9TC84vjkBdlWqLiOnAlMx8YIRqqlRX3Tn1/z1fHhFTGpwfDlVqOx44PiIej4i19VVV26EuAOpDjVOBh0agLqhW20LgryKih9rMv8tHprRqSg30SksRjJLKtUXEXwFdwI3DWlH97Roce0ddmfntzPxPwP8E/tewV1UzaG0R8QfALcCXRqiet966wbH+92wl0JmZHwVWA3cPe1U1VWobQ23Y5TRqPeE7IuLdbVDXPvOA5Zn522Gsp68qtX0OWJKZk4HZ1J6/aZscbZtCWqzKcgWjpVJtEXEG8HfA2Zn5ervU1ccy4DPDWtHvDVXbeOBEYE1EbAVOAVaMwAejQ96zzNze5/f3XWrfGTASqi7ZcX9mvpmZzwPPUQv40a5rn3mM3HALVKvtIuAHAJn5E2ActXVe2sNoD+IP04cbY4BfUPvftX0fbnx4gLZLGNkPRYesDZhO7cOZ49qsruP6bH+a2pPCbVFbv/ZrGJkPRavcs2P6bP83YG273DNgFnB3fXsiteGGo0e7rnq7/wxspf6sTBvds1XA/Pr2B6kF/ojVOOS/w2gXMIy/nNnAv9WD8e/qx66j1uMF+C/U/kb+f9S+ZWlTG9W2Gvg1sKH+WtEmdf0DsKle08ODhepI19av7YgEesV79r/r9+zp+j07oV3uGbUhhpupfbfBM8C8dqirvr8Q+PuRuldN3LMPAY/Xf58bgDNHusbBXj4pKkmFKHUMXZIOOQa6JBXCQJekQhjoklQIA12SCmGgS1IhDHRJKoSBLkmF+P9HSf5nOKSJ6AAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"# Plot a histogram of the observed data\n",
"plt.hist(d2, 50, density=True, color='C1', label='Actual Data');\n",
"plt.legend();"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This data, as plotted above, we might think looks quite normal, based on the visualization. \n",
"\n",
"However, there are more specific checks that we can do. "
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXhTZdrH8e+dpKUFStm3shQFRGQX2QQVRUVGBBQUFYVxYRh1fHWc13F0VAZfHbfBGcUZNxBUBJUdEXRAUHYo+77IImUtBVpaStskz/tHQpqWlqaQ9qTJ/bmuXnOe5CT9mXZuTp/zLGKMQSmlVPlnszqAUkqp4NCCrpRSYUILulJKhQkt6EopFSa0oCulVJhwWPWNa9asaRITE6369kopVS6tWbPmuDGmVmHPWVbQExMTSUpKsurbK6VUuSQi+4t6TrtclFIqTBRb0EUkRkRWicgGEdkiIn8r5JxhIpIiIuu9X4+UTlyllFJFCaTLJRu40RiTISJRwBIRmWuMWVHgvK+MMU8EP6JSSqlAFFvQjWdtgAxvM8r7pesFKBXGcnNzSU5O5uzZs1ZHiVgxMTE0aNCAqKiogF8T0E1REbEDa4CmwPvGmJWFnHaXiFwH7ASeNsYcCDiFUiqkJCcnExcXR2JiIiJidZyIY4whNTWV5ORkmjRpEvDrAropaoxxGWPaAQ2ATiLSqsAps4FEY0wbYD4wobD3EZHhIpIkIkkpKSkBh1RKla2zZ89So0YNLeYWERFq1KhR4r+QSjTKxRhzClgE9C7weKoxJtvb/Bi4uojXf2SM6WiM6VirVqHDKJVSIUKLubUu5vMvtstFRGoBucaYUyISC/QC3ihwTj1jzGFv8w5gW4mTKHURcl1uth5KZ/W+EyTtO8np7Fya1KzE5bUq07R2ZVonxFO1YrTVMZUqE4H0odcDJnj70W3A18aYb0VkFJBkjJkFPCkidwBO4AQwrLQCKwVw7PRZXp+7nXmbj3Amx5XvuaW7U33H0XYbv702kcd6NiU+NvCbS8p6IsIf//hH/vGPfwDw9ttvk5GRwciRI8ssw7Bhw7j99tsZOHDgeY//9NNPVKlShaysLLp06cLf//53EhISAOjTpw9ffvklVatWLfR9//nPfzJ8+HAqVqwY1LzFdrkYYzYaY9obY9oYY1oZY0Z5H3/JW8wxxvzFGHOVMaatMaanMWZ7UFMq5eV2Gyau3M9N//iJaWsPnlfMC8pxufnw5z1c/9ZCxi3ZS47TXUZJ1aWqUKEC06ZN4/jx4xf1eqfTGeRE+b311lts2LCBHTt20L59e3r27ElOTg4A3333XZHFHDwF/cyZM0HPZNnUf6VKavexDJ6bupGk/SfzPZ5QNZZrEqtxTZPq1ImLYc/xDHYfy2DDgTR2HD0NwKkzuYz6ditfrvqVT4ddQ8Pqwb0yUsHncDgYPnw477zzDq+++mq+5/bv389DDz1ESkoKtWrV4tNPP6VRo0YMGzaM6tWrs27dOjp06EBcXBx79+7l8OHD7Ny5k9GjR7NixQrmzp1LQkICs2fPJioqilGjRjF79myysrLo1q0bH374YcB92CLC008/zfTp05k7dy79+vXzLW0SGxvL3XffTXJyMi6XixdffJGjR49y6NAhevbsSc2aNVm4cGHwPrOgvZNSpWjDgVMMGbuS02fzrroSa1Tk1QGtubZpzQJn1wE8Q7/mbDrMG/O2c+BEFuD5R2HgB8v44uHONKsTV1bxy7XE5+aU2nvve/03F3z+8ccfp02bNjz77LP5Hn/iiSd48MEHGTp0KOPGjePJJ59kxowZAOzcuZP58+djt9sZOXIkv/zyCwsXLmTr1q107dqVqVOn8uabbzJgwADmzJlD//79eeKJJ3jppZcAeOCBB/j222/p27dvif5bOnTowPbt2+nXr5/vsXnz5lG/fn3mzPF8hmlpacTHxzN69GgWLlxIzZoFf3cvja7lokLexuT8xdxhE57o2ZR5T11XSDHPIyLc3qY+8/94PX/9zZVEOzy/7kfTs7n7w+VsOHCqTPKri1elShUefPBB3n333XyPL1++nPvuuw/wFOAlS5b4nhs0aBB2u93Xvu2224iKiqJ169a4XC569/YM0mvdujX79u0DYOHChXTu3JnWrVvz448/smXLlhJnLWx/5tatWzN//nz+/Oc/s3jxYuLj40v8viWhBV2FrpHxbHqpHUPG/OAr5tU4zYzHr+VPt15BTJS9mDfwqOCw80iPy5jw205Uiva85uSZXO77eAXLfrm4/llVdp566inGjh1LZmZmkef4d49UqlQp33MVKlQAwGazERUV5TvXZrPhdDo5e/Ysjz32GFOmTGHTpk08+uijFzVDdt26dVx55ZX5HmvevDlr1qyhdevW/OUvf2HUqFElft+S0C4XFbI2uxszJOcvpOP5P2hVTjMx+lVaJgy+qPfrenkNvny0C0M/XcWpM7lk5rh4aPxqZjx+LS3qVglm9LBSXLdIaatevTp33303Y8eO5aGHHgKgW7duTJ48mQceeICJEyfSvXv3i37/c8W7Zs2aZGRkMGXKlPNGtVyIMYb33nuPw4cP+67+zzl06BDVq1dnyJAhVK5cmfHjxwMQFxfH6dOntctFRYbUjGwezvlf0qgMQDwZTIx+jZa2X2FkfOFfAWjbsCpf/64rteM8V21nc938/ou1nD6bW2r/LerSPfPMM/lGu7z77rt8+umntGnThs8//5x//etfF/3eVatW5dFHH6V169b079+fa665JqDX/e///i9t27alefPmrF69moULFxIdnX/Ow6ZNm+jUqRPt2rXj1Vdf5a9//SsAw4cP57bbbqNnz54APPLII0HZH0IK6/cpCx07djS6wYUqjDGGhyck8eP2YwBUIZMvo1+llW3fhV84Mi3g77Hz6Gn6jVlKVq5n2GPvq+rynyEddHak17Zt287rPlBlr7Cfg4isMcZ0LOx8vUJXIWfc0n2+Yg7wz6j3iy/mJdS8Thyv39Xa15635Qhjl+wN6vdQqqxpQVchZVNyGq/PzVs54hH7HG60ry+V79WvXQJDuzb2tf8+dzur9p4ole+lVFnQgq5CRka2kz9MWkuuy9MN2Fr28Kxjcql+zxd+05J2DT0z+lxuw5OT1pGRXbozDJUqLVrQVcj4v2+3si/VMx26cgUH70W9R7RceGr/pYp22Pj3/R2oVtGzzsuR9LO889+dpfo9lSotOmxRhYR1v55k8uq8PVFeHdCKxBlHS/YmRY10KeZmaf2qsbzc9yqe+srTtTN+2T7u6tCAlvV1KKMqX/QKXVnO7Ta8PCtvZt7NLevQr11CmWbo164+3S6vAXi6Xv46YxNut+60qMoXvUJXlvs66QAbkz1X0dEOGy/d3rLMM4gIo/q14rZ//Uyuy7D211N8lXSAezs1KvMsISnAcf6Bv19gQ0ynT5/OnXfeybZt22jRosUFzx0/fjy33HIL9evXv6hIixYt4u233+bbb7897/F+/fpx2WWXcebMGerUqcOzzz7L7bffXuz7RUdH061bt4vKczH0Cl1ZKu1MLm9+v8PXHnH95ZathNi0dmV+d93lvvbrc7eTmpF9gVeo0jZp0iS6d+/O5MnF3xwfP348hw4dKpUcPXr0YN26dezYsYN3332XJ554ggULFlzwNYsWLWLZsmWlkqcoWtCVpUb/dwcnMj1rSCeQwu+XXFuimZ8BKWpmaSHf44kbm9KweiwAaVm5vDFPl/a3SkZGBkuXLmXs2LHnFfQ333yT1q1b07ZtW5577jmmTJlCUlIS999/P+3atSMrK4vExETf7NKkpCRuuOEGAFatWkW3bt1o37493bp1Y8eOHQW/9QW1a9eOl156iTFjxgAwe/ZsOnfuTPv27enVqxdHjx5l3759fPDBB7zzzju0a9eOxYsXF3pesGlBV5bZdjidz1fs97VfjPqCWMmxMBHERNkZdUfeHuhT1iSz+9hpCxNFrhkzZtC7d2+aN29O9erVWbt2LQBz585lxowZrFy5kg0bNvDss88ycOBAOnbsyMSJE1m/fj2xsbFFvm+LFi34+eefWbduHaNGjeL5558vcbZzS+UCdO/enRUrVrBu3ToGDx7Mm2++SWJiIiNGjODpp59m/fr19OjRo9Dzgk370JVl3pi3nXP3HXvYNnKrbbW1gbx6tqjNdc1r8fPOFNwGRv93J/++v9B9z1UpmjRpEk899RQAgwcPZtKkSXTo0IH58+fz29/+1rd9W/Xq1Uv0vmlpaQwdOpRdu3YhIuTmlnwdH/8lU5KTk7nnnns4fPgwOTk5NGnSpNDXBHrepdArdGWJtb+eZNGOFABsuHnJ8TmhtIzKn25p7jv+btMRNh8MfJ0YdelSU1P58ccfeeSRR0hMTOStt97iq6++whiDMSagNXccDgdut2fLQf/lcF988UV69uzJ5s2bmT179iUvlfuHP/yBJ554gk2bNvHhhx8W+X6BnncptKArS/xz/i7f8R22ZTSzHbQwzfnaNKjKrVfV8bVH62SjMjVlyhQefPBB9u/fz759+zhw4ABNmjRhyZIl3HLLLYwbN863J+eJE57lGs4tSXtOYmIia9asAWDq1Km+x9PS0nybOZ9bzrYkNm7cyCuvvMLjjz9+3vtNmDDBd17BPEWdF0za5aLK3Jr9J/h5p/fqXOBJxzTrwlxgMtIzt1zBD1uPYgz8uP0Ya/af5OrG1co2X6gowUqWwTBp0iSee+65fI/dddddfPnll/znP/9h/fr1dOzYkejoaPr06cNrr73GsGHDGDFiBLGxsSxfvpyXX36Zhx9+mNdee43OnTv73ufZZ59l6NChjB49mhtvvDGgPIsXL6Z9+/acOXOG2rVr8+6773LTTTcBMHLkSAYNGkRCQgJdunRh717PIm99+/Zl4MCBzJw5k/fee6/I84JJl89VZW7IJytZstsz+uDO9gmM3na9xYkK4S1gT01ex4z1nqFwXS+rwaThXaxMVWZ0+dzQEPTlc0UkRkRWicgGEdkiIn8r5JwKIvKViOwWkZUikniR+VU4GxnP6pc6+Yq5HRd/2HKPxaEu7KlezbHbPP21y/eksnS3blmnQlcgfejZwI3GmLZAO6C3iBS8THkYOGmMaQq8A7wR3JgqXLzjvMt33N+2hCa2IxamKV5izUoMurqBr/2vBbsucLZS1iq2oBuPDG8zyvtVsJ+mH3Cul38KcJPo1i+qgFXuK1jm9ozxtuPiScd0ixMF5g83NcPhvUpftfcE6w+csjhR2bCqO1Z5XMznH9AoFxGxi8h64BjwX2PMygKnJAAHvCGcQBpQo5D3GS4iSSKSlJKSUuKwqnz70NnXd3ynfTGNbccucHboSKgaS9+2eeuDfPTzLxamKRsxMTGkpqZqUbeIMYbU1FRiYmJK9LqARrkYY1xAOxGpCkwXkVbGmM1+pxR2NX7eb4Ix5iPgI/DcFC1RUlWu7T6WwQJ3B197hH22hWlK7tEelzF9nWdo5bzNR9ifmknjGpUsTlV6GjRoQHJyMnrhZZ2YmBgaNGhQ/Il+SjRs0RhzSkQWAb0B/4KeDDQEkkXEAcQDupeX8hm3NG+IVi/bGi63HbYwTcm1rF+FHs1qsnjXcdwGPlm8l1f6tyr+heVUVFRUqcxkVKUrkFEutbxX5ohILNALKLhi0SxgqPd4IPCj0b/VlFdqRjZT1yT72o845liY5uL5r8T4zZoDvkXFlAoVgVyh1wMmiIgdzz8AXxtjvhWRUUCSMWYWMBb4XER247kyH1xqiVW5M3Hlr2Q7PVOwW8leOks5WMGwkAlH1xpoWe9bth5O52yum8+W7+OpXs3Pf61SFglklMtGY0x7Y0wbY0wrY8wo7+MveYs5xpizxphBxpimxphOxpg9pR1clQ9nc118tnyfr/2oY05IrdlSEiLwu+sv87U/W76frJzS3fNUqZLQtVxU8PmtNz7rbwM4nuHpmqhHKn1sBQdIlS99WtcjoapnadYTmTlMXZtczCuUKjta0FWpMQY+cfXxtYc5vidKyvcVbZTdxkPd824WfrZ8nw7tUyFDC7oqNUvcrdhpGgJQiSwG23+0OFFwDOrYgIrRdgB2Hs1g1V4d0KVCgxZ0VWo+d93sOx5k/4l4OWNhmuCpEhPFgPYJvvZnfrsuKWUlLeiqVBw21ZnvztvlZ4h9voVpgu+Bro19x99vPsKx9OBvVqBUSWlBV6VisrMnbu+vV1fbFpraSmc3dqu0qFuFTomerc+cbsOkVQcsTqSUFnRVCnKNncmunr52uF2dnzPE7yr9y1X7yXW5LUyjlBZ0VQoWuDtwFM/Va01OcbMtPDcy6X1VXWpWrgDA0fRs5m89anEiFem0oKug+8LVy3c82L6Q6HI+VLEo0Q4b93Zq6Gt/rjdHlcW0oKug2ns8kyXu1gDYcHOvIzyGKhblvs6NfDsaLfslld3HThfzCqVKj24SrYJqot9V6o22dSRIqoVpSkGBNV7qAb14iu/pBMDkVQf46+0tLQimlBZ0dSkKFLezJoop2WOAOADuD9OboQXdZ/+R792egj5t3UGe7d2CaIf+8avKnv7WqaD53n0Np7zFvKEc43rbRosTlY3utk3Uj/fsLHMiM4f52/TmqLKGFnQVNF+7bvAd32NfiE0iY40TuxgGdcy7OTp5tY5JV9bQgq6C4oC7Jku9G0DbcHOXfbHFicrWoI4NfMsCL96VwsFTWdYGUhFJC7oKim9c1/uOr7NtpJ5E1oJVDapVpHvTmoBnlclvkvQqXZU9LejqkrmMMMWvoN9tX2RdGAvdc01et8s3Scm43ZHR5aRChxZ0dcmWultxCM/VaXXS6WVbY3Eia9zcsg7VKkYBcPBUFkt/OW5xIhVptKCrS+Z/M7S/fUnYzgwtTgWHnQHtG/jaenNUlTUt6OqSnDKV+MHd0de+2/6ThWms59/t8t8tRzmRmWNhGhVptKCrSzLDdS05eLoZ2sovtLBF9lXpFXXjaNuwKgA5Ljez1h+0OJGKJMUWdBFpKCILRWSbiGwRkf8p5JwbRCRNRNZ7v14qnbgq1Ph3twyK0JuhBQ26Oq/bZepaLeiq7AQy9d8JPGOMWSsiccAaEfmvMWZrgfMWG2NuD35EFaq2uhux1SQCUIEc+tqXWxvISn7LIPQ1lRjF++QQzaaDaew8eprmdeIsDKciRbFX6MaYw8aYtd7j08A2IOHCr1KRYJqrh+/4VtvqsNkz9FLFSyY329b62lPXJFuYRkWSEvWhi0gi0B5YWcjTXUVkg4jMFZGrinj9cBFJEpGklJSUEodVocPpcjPDda2vfWeEzQwtzl32n33H09cdxKm7GakyEHBBF5HKwFTgKWNMeoGn1wKNjTFtgfeAGYW9hzHmI2NMR2NMx1q1al1sZhUCFu86znE8N/9qcZLuts0WJwot19k2UpNTABw7nc2S3TomXZW+gAq6iEThKeYTjTHTCj5vjEk3xmR4j78DokSkZlCTqpAydW1eN8IA+1Icoleg/hzipr99qa89TW+OqjIQyCgXAcYC24wxo4s4p673PESkk/d9w2xnA3VOWlYuP/jtn6ndLYXz73b5fssR0s/mWphGRYJARrlcCzwAbBKR9d7HngcaARhjPgAGAr8XESeQBQw2xuhCFuGiwEYWc503kOMcDkBL2RfxY8+LcqXtAC3rVWHr4XSynW6+23iYwZ0aWR1LhbFiC7oxZgkgxZwzBhgTrFAqtE11Xec79r8KVee7s0MCW+d4bjlNXZusBV2VKp0pqkpkv7s2q00LAOy4uMO+zOJEoa1fuwTfJtKr951kf2qmxYlUONOCrkpkuru77/h62wZqScEBT8pfrbgK3NA8b0TXzPWHLEyjwp0WdBUwY2C6K6+gR9quRBerf/u8eXjT1x1Eby+p0qIFXQVsrWnGflMXgDgyuclvNqQq2s0t61C5gud21d7jmWxITrM4kQpXWtBVwPxnhv7GvpIY0WF4gYiJsnNbq7q+9ox1OiZdlQ4t6CogOcbOt66uvnZ/+xIL05QjI+NhZDwDNo7wPTR72QZydSkAVQq0oKuA/Oxuy0k8KwbW5zidZIfFicqXzrZt1PXOtUslnsW7dC0jFXxa0FVApvt1t/SzL8UmemOvJOxi6Oe3FMD0dTraRQWfFnRVrHQTy3z31b72AO1uuSgD/Ar6D1uOcFqXAlBBpgVdFWueqxPZRAOeqf7NbXpT72K0sB2ghewHINvpZt7mIxYnUuFGC7oq1gx3XneLXp1fmgH5ul30H0YVXFrQ1QUdTstiubslAIJbp/pfon72pQieES7L96RyJO2sxYlUOAlktUUVKQqsqggwy3k7hvsAuNa2hTpyqqxThZW6cpKutq0sc7fCGJi94RCPXneZ1bFUmNArdHVB+Ua32JZe4EwVqP427XZRpUMLuirSDncDtpvGAFQgh9721RYnCg+97auIdnj+r7f1cDq7jp62OJEKF1rQVZH8p/r3sq0hTrIsTBM+qkgWva6s7WvPWK9X6So4tKCrQrmNMNPlP7pFu1uCqV+7vBUYZ6w7hNutE7XUpdOCrgq12lzBITz7fFflNNfZNlicKLzccEUt4mOjADh4Kos1v560OJEKB1rQVaEKrqwYLS4L04SfCg47fVrX87V1BUYVDFrQ1XmyjYPvXJ197f7a3VIq+rer7zues+kwOU5dgVFdGi3o6jyL3O1IozIADeQYV8tOixOFp2sSq5NQNRaAU2dy+WmnrsCoLk2xBV1EGorIQhHZJiJbROR/CjlHRORdEdktIhtFpEPpxFVlYWa+sefLdGXFUmKzCXf4XaVrt4u6VIFcoTuBZ4wxVwJdgMdFpGWBc24Dmnm/hgP/CWpKVWY8Kyu297W1u6V09fcb7TJ/21FdgVFdkmILujHmsDFmrff4NLANSChwWj/gM+OxAqgqIvVQ5c48Vydy/FZWbKYrK5YO705GV3zQgCtlH6ArMKpLV6I+dBFJBNoDKws8lQAc8Gsnc37RR0SGi0iSiCSlpGh/YSjyX1lRr87Lhv8Yf51kpC5FwAVdRCoDU4GnjDHpBZ8u5CXndbwaYz4yxnQ0xnSsVatWyZKqUnfEVNOVFS1wh32ZbwXGZb/oCozq4gVU0EUkCk8xn2iMmVbIKclAQ792A0D32CpnZrm6Yby/Et1sW6krOtmlLJxbgRHwrcCo1MUIZJSLAGOBbcaY0UWcNgt40DvapQuQZow5HMScqgz4r6zY36YbWZQlXYFRBUMgV+jXAg8AN4rIeu9XHxEZISIjvOd8B+wBdgMfA4+VTlxVWrYfSWebSQR0ZUUr9LavIpocwLMC405dgVFdhGI3uDDGLKHwPnL/cwzweLBCqVJWyEYWM3IHA3cAcLOurFjmqkgWN9vWMsfdBfCMSX+2dwuLU6nyRmeKqkJWVtTuFiv08xvtMnO9rsCoSk4LumKluwWHqQFANU5znW2jxYki0w229VStmLcCY9J+vSmtSkYLumKGu7vvuK99OVG6sqIlosXFb/xWYNSbo6qktKBHuLMmqsDKitrdYqUB7fPm483ZeIizufqPqwqcFvQI96O7PaepCEBjOUJ72W1xosh2deNqNKzuWYEx/ayTH7cfsziRKk+0oEe4aa687pZ+tqXIBcczqdImIgxo38DXnrZWu11U4LSgR7BUE8cidztf+07tbgkJ/t0ui3Yc40RmjoVpVHmiBT2CzXZ1xemdinC17CDRdtTiRIqR8TQZU5/2sgsAp9sw++/3WxxKlRda0CPYNFcP37GOPQ8td9oX+479u8WUuhAt6BFqt7s+G83lAESTy+32FRYnUv5ut68gCicAG0xTfknJsDiRKg+0oEco/6vzm2xrqSqZFqZRBVWTDHra1vna0/XmqAqAFvQI5DbCdL8/4/3/vFehw//nMn3dQV0KQBVLC3oEWuG+Mt9U/+ttGyxOpArT07aeeDxdLQdPZbFq3wmLE6lQpwU9Ak1z53W33GFfRrRO9Q9JFcSZ797GtLXJFqZR5YEW9AhzJsfJXFcnX1u7W0LbXfaffcdzNh7mTI7TwjQq1GlBjzDzNh8hE8/U8svlIG1kj8WJ1IW0l91cJp4t6TJzXMzddMTiRCqUaUGPMN8k5f3Zfqd9sU71D3EiMMj+k689ZY12u6iiaUGPIAdOnGH5nlQAbLi5S7tbyoUB9iXYvP/wLt+TyoETZ6wNpEKWFvQI4n9118O2kbqiGyiUB3XlJD2a1fK1dcEuVRQt6BHC7Tb5Crr/n/Eq9A28Om8FxilrD+iYdFUoLegRYsWeVA6e8mz8HE8GvWxrLU6kSuLmlnWoEuNZSO3ACR2TrgrnKO4EERkH3A4cM8a0KuT5G4CZwF7vQ9OMMaOCGVJdhJHx+Zrf5Pwe8Iw/72dfRozkWhBKXayYKDt3tKvPFyt+BTzdZ10uq2FxKhVqArlCHw/0LuacxcaYdt4vLeYhJt3EMtedN/Zcu1vKp4FXN/Qdf7fpMJnZOiZd5VdsQTfG/Azo33fl2BxXF85SAYAW8iutZG8xr1AhZ2Q8bT9pTFPx3Ac5k+Nizit3WhxKhZpg9aF3FZENIjJXRK4q6iQRGS4iSSKSlJKSEqRvrYrzjet63/FA+0869ryc8oxJz5s5+pXrBuvCqJAUjIK+FmhsjGkLvAfMKOpEY8xHxpiOxpiOtWrVKuo0FUS73fVZa5oD4MCpG1mUc3fZf8bhXSd9jbmCXUdPW5xIhZJLLujGmHRjTIb3+DsgSkRqXnIyFRSTXDf6jm+yraWGaAEoz2pKOjfb1vjak1cfsDCNCjWXXNBFpK6I5494Eenkfc/US31fdenOmqh8G1nca//RwjQqWAbbF/qOp61NJtupq2Uqj2ILuohMApYDV4hIsog8LCIjRGSE95SBwGYR2QC8Cww2xuishxDwvbsjJ4kDIIEUetg2WZxIBUN32yYS8NyDOnkmlx+26ObeyqPYcejGmHuLeX4MMCZoiVTQTPbrbrnHsQi76L+z4cAuhrsdi3jHOQiAr1YfoG/b+hanUqFAZ4qGqb3uuix3ewYc2XDr2PMwM8j+EzbcACzZfZxfU3XBLqUFPWxN9hvSdqNtHfVEpxKEk/pyIt/WgV8n6c1RpQU9LOU43a7EgNQAABcaSURBVEz1G3vufxNNhY97/H6u36w5gNPltjCNCgVa0MPQgm1HOY5nLZe6pHKDbb3FiVRpuMm2jpqVPTOAj6Zns2D7MYsTKatpQQ9DX6761Xd8t/0nHKJXbuEoSlzcc03esrpfrNhvYRoVCrSgh5l9xzNZvOs4AIKbux2LrA2kStW9nRr5djNavOs4e1IyrA2kLKUFPcx87neV1tO2ngZy3MI0qrQ1qFaRG1vU8bUnrvz1AmercFfsOHQV4vzWPc8y0XyT/T5QCYAH7T9YFEqVmZHxPOBqw3yeA+CbJZv50+rrif2bLn4XifQKPYzMdHUj3VvME+UI1+nM0IjQw7aJxnIEgHQqMdvV1eJEyipa0MOEMTDBdYuvPcT+X2w6MzQi2MQwxD7f1/7MdQu6+kZk0oIeJtaY5mwziQDEkJ1v3WwV/gbaf6YCOQBsNk3YkJxmcSJlBS3oYWKCM+/qvL99KfGSaWEaVdaqSQZ97ct97c+X6xDGSKQFPQwcM/HM89sz9AH7fy1Mo6zi/3OfvfEQxzOyLUyjrKAFPQxMdt1IrnfAUkfZwVU2vTqLRG1te2grvwCe5R8mrtAhjJFGC3o5l20cfO7s5Ws/6NChipHsIcd3vuPPV+zjbK5ufhFJtKCXc7Nc3UihGgB1OEFv2yqLEykr9bGtop53w7DjGTnM2nDI4kSqLOnEovLCbwLROcbAWNfrvvZQx/dEi16RRbIocTHU8T2vO+8DYNySvQy6ugHeXSJVmNMr9HJssbs1200jACpylvt1z1AF3GtfSGyUHYDtR06z7Bfd4jdSaEEvxz52/cZ3fLd9kQ5VVADESyZ3d8xbhfGTxXssTKPKkhb0cmq7uyGL3W0AzxZzD9nnWZxIhZLfXtuEc70sC3eksPuYrsIYCbSgl1OfuPr4jm+1raaRTTc3UHkSx9TnJknytT/914uF3odR4aXYgi4i40TkmIhsLuJ5EZF3RWS3iGwUkQ7Bj6n8HTNVmem61td+xG+omlLn+P9eTHFdR4qpYmEaVRYCuUIfD/S+wPO3Ac28X8OB/1x6LHUhE5y3+CYStZddXG3bZXEiFYo6y3bfRKNsohnr7FPMK1R5V2xBN8b8DFxoy/h+wGfGYwVQVUTqBSugyi/NVOQzv1UVH3XMsTCNCmUi8Jhjpq/9hasXaWdyLUykSlsw+tATgAN+7WTvY+cRkeEikiQiSSkpugD/xZjgupXTVATgMjnErbbVFidSoexm2xqaSTIAGVRkwvJ9luZRpSsYBb2wGQuFLsZsjPnIGNPRGNOxVq1aQfjWkSXDxDDOmdf79YRjBnZd81xdgE1Mvqv0cUv3kpnttDCRKk3BKOjJQEO/dgNA5xuXgomuXpwiDoCGcow7bMssTqTKg7625TQUzyioU2dymbRKF+0KV8Eo6LOAB72jXboAacaYw0F4X+Uny0TzsTNvItFj9pk4xG1hIlVeOMTNCPssX/vjxXvIduoSEeEokGGLk4DlwBUikiwiD4vICBEZ4T3lO2APsBv4GHis1NJGsMmunhzHM464HqncpTsSqRK4y76Y2pwE4Gh6NlPWJFucSJWGYhfnMsbcW8zzBng8aInUebKdLj503u5rj3DM1kW4VInESC6POubwqnMIAP9e+AsDr25ABYfd4mQqmHSmaDnwdVIyR6gBQE1OcY99ocWJVHl0n30BNSpFA3DwVBaTVmpferjRgh7isnJcvLsgb+LQ7xzfEiM6lliVXCXJ5rHsT3ztMbOXcebl2rokQBjRgh7iPl22l5TTnr0h63CCIfb5FidS5dn99gV5G2BQlfF+k9RU+acbXIQav6ulNFOJD7L/CVQC4H8c04iVHIuCqXAQI7n8wTGd552PAPChsy/32xeg1+jhQa/QQ9gHzttJ9xbzRDnCIPtPFidS4WCQ/ScayxEA0qjMJ37DYVX5pgU9RB01VfnUlTcr9BnH10TpyBYVBFHi4mnHVF97rOs2jmdkW5hIBYsW9BD1rvNOzlIBgKtkL7+xrbQ4kQonfW3LuEI8o1zOEMOYH3dbnEgFgxb0ELTXXZevXDf42s86vsKma7aoILKL4RnHN7725yv2s+voaQsTqWDQgh6C/s85BKf3fnUX2xaus220OJEKRzfb1tDFtgUAl9vwf3O2WZxIXSot6CFmoasdC9yeTZ8EN887vvTtDalUMInAS47PseFZE+innSks3K5bGZZnWtBDSI7TzSjnA772PfZFtLHttTCRCnctbb/mm3n8ypyt5Lp00bfySgt6CPl06V72Gs9mT3Fk8ifH1xYnUpHgGcc3xHEGgD0pmXz28r2e+RDnvlS5oQU9RBxLP5tviv/TjqnUlHQLE6lIUVPSedIxzdf+p/MuUk2chYnUxdKCHiJen7edzBzPOPNmkswD9v9anEhFkqH272kinm0MTlOJ13LvsziRuhha0EPA8l9Smbb2oK/9suMznUSkylS0uHjJ8bmvPdV9PYtdrSxMpC6GFnSLZeW4eG5a3rDE3rZVdLdvtjCRilQ97eu53bbc137e+QhZJtrCRKqktKBb7B8/7GB/queGVFyMg79Fjbc0j4psL0dNIJ4MAA6Y2rzjHGhxIlUSWtCtMjKedS91YNySX3wPveh8nzpyysJQKtLVknRecEz0tT9x9WFTcpqFiVRJaEG3SLZx8Gzu73B7fwQ9bBt1NUUVEgbZf6KbzdPt58bGn6du1LHp5YQWdIu87+zPLtMAgIqc5TXHJzojVIUEEXjNMZYKeNbe33o4Pd+QWhW6tKBbYM3+E7zv6udr/9kxmYa24xYmUiq/RNvRfIt3jVm4m5V7Ui1MpAIRUEEXkd4iskNEdovIc4U8P0xEUkRkvffrkeBHDQ9pZ3J5ctJ6XHh2W+8k23TMuQpJj9i/83W9GANPf/QtaS/X1xmkIazYgi4iduB94DagJXCviLQs5NSvjDHtvF+fFPJ8xDPG8OepGzl4KguAeDJ4J/rfujSuCkk2MYyO+g9V8Syre4iaPJ/7MEZ/XUNWIFfonYDdxpg9xpgcYDLQr5jXqEJ8sWI/87Yc8bXfjPqIBNE/Y1XoqisneSPqY197jrsL37iutzCRupBACnoCcMCvnex9rKC7RGSjiEwRkYZBSRdGth5K5xW/9aaH2r/nVnuShYmUCsyt9iTus8/3tV92DmWru5GFiVRRAinohY29KPhH12wg0RjTBpgPTCj0jUSGi0iSiCSlpKSULGk5dupMDo9NXEOO0zP068p6VfiL40uLUykVuBcdX9BUkgHIIobhuX/kRGaOxalUQYEU9GTA/4q7AXDI/wRjTKox5twusx8DVxf2RsaYj4wxHY0xHWvVqnUxecudXJeb33+xln3e2aCxUXbG3NeeGMm1OJlSgYuVHP4T9S8qe5fZTTa1eWziGh2fHmICKeirgWYi0kREooHBwCz/E0Sknl/zDkD3ssJzE/Tll55hud9wr3d4m8vfL6zHSqnQ1sx2kH9G/Rvx7nC0Ys8J/u/brRanUv4cxZ1gjHGKyBPA94AdGGeM2SIio4AkY8ws4EkRuQNwAieAYaWYudyYsGwfX7p6+dp/cnxFb/tqCxMpdWl62dfyJ/MNbznvAWDC8v20WP1X7nXk7XrESF0qwCrFFnQAY8x3wHcFHnvJ7/gvwF+CG618W7jjGKP8rl762ZbyuH2mhYmUCo7H7DPZ6m7MHHcXAP7qfIiaksbN9rUWJ1M6U7QUrNiTyu+/WIPbe+u4nezijaiPdGq/Cgsi8FbUh7SUfQC4sPN47pMscxU2PUWVJS3oQbb215M8PH41Z3M9/YwJpPBR9Gi9CarCSkXJZnz0GzQWz7yKHKJ5NPcZNrqbWJwssmlBD6Ith9IYNm6Vbyu52nEVmBj9GrVF+xRV+KktaXwR9XfqcAKATGIZmvMcu46etjhZ5NKCHiTbj6TzwNhVpJ91AlC9UjQTH+lMou2oxcmUKj0NbSl8Hv133/IAJ4njno9WsPmgXsRYQQt6EKzck8qgD5b7JlpUiXHw2UOdaFZHd05X4a+57SATot+gEp41ik5k5jD4vR9Y+VLnvIW8dDGvMqEF/RLN23yEB8at4rT3yrwyZxjvfp5WHzfSX2IVMdra9vBZ9OtUIROADCryYM5zLHC1tzhZZNGCfgkmrtyfb0p/LU7yVfQrdLDttjiZUmXvatsuvo4eRS1OApBNNMNz/8gkZ0+Lk0UOLegXIcfp5uWZm3lh+mbf0MREOcK06JFcZdtvbTilLNTCdoAp0X+joRwDPEMa/+J8lBdyH/Jd+KjSowW9hA6dyuKej5YzYXle4W7TIJ4p0SNpaIucBceUKkpj2zGmRI/kSu84dYCJrl7c+/EKjqWftS5YBBBj0Wr1HTt2NElJ5Wv52CW7jvPk5HX5Vpm7rVVd3h7Ulkp/r2FhMqVCzxlTgT/nPspsdzffY7U5yeio/9DdvjnvRF0qoEREZI0xpmNhz+kVegAys528PHMzQ8au9BVzu014oc+V/Pv+DlSqENAKCkpFlIqSzbtRY3jB8QU274Jex6jGkNzneTF3GJmmgsUJw49WomIs2XWc56ZtJPlklu+xWnEVGHNvezpfplflSl2ICDzq+I4r5VeezH2CE1QB4HPXLfzkbsvbUR/QyeKM4US7XIpwNP0sb3+/g2/WJOd7vKdtHW9EfUxtOWVRMqXKpxRThRdyH+YH9zX5Hh/QPoE/925B3fgYi5KVLxfqctGCXkBmtpMPf97Dxz/vISvX5Xs8PjaKl/u2ZMDMq3SRLaUukjEw3d2dl3OHcppKvsdjOcsIx2yG2+cQK957VNq3XqgLFXTtcvE6k+Pkq9UHeH/hLxzPyM733K1X1eGV/q2oHRdTYGsPpVRJiMCd9iV0tW1lZO6DfO/2dLhkEcM7zkFMdPZiuGMO99kXUNHirOVRxF+hp2ZkM2H5fj5bvo9TZ/KviNiibhwv/OZKejTz2y5PZ38qFTTLXC15xTmEbSYx3+PVOM1vb76aB7s2pmrFaGvChSjtcinAGMPKvSf4JimZOZsO+Za6PadulRieuaU5d3ZogN1WoH9FC7pSQeUywteuG3jHOZBjVMv3XLTDxm9a12PwNQ3p1KQ6ov2dWtDP2X0sg7mbDjNlbTL7vZs2+2tYPZZHe1zGoKsbEvta9TLNplSkO2uimOrqwQeuOzhgap/3/GVyiP72pfR58j2a1o7che8itqA7XW42JKcxf9tRvt9yhD0pmYWe1yqhCr+77nJua1UXh907NF+vxJWyhNPYmOXuxqfO3mwylxV6TtPalbmtVV2ub16Ltg2rEmWPnCk1EVPQs50uth0+zeq9J1j2y3FW7ztJRraz0HPjYhz0a1efuzs2pHVC/Pl/ymlBV8pym92JTHLdyExXNzKKuE1aiSy62LbR9bb7ad+oGlfVr0JMlL2Mk5adsCzoJzNz2Hn0NLuOZbD1cDqbktPYfiSdXFfR/z2xUXaub16L21rX5dar6np+6Fq4lQp5maYCC9wdmOfqxEJ3W7Ioesx6lF24sl4VWifE06JuHM3rxHFF3biwubkaNsMWdx09zYszN7P7WAbHM3KKfwFQPz6Ga5vW5Jar6tK9aU1io8P3X26lwlUlyeYO+3LusC8ny0Tzk7sti9xtWexqzUFq5Ts312XYmJzGxuT849hrVIqmUY2KNK5ekUY1KtGgaix14mOoFx9DnSoxVIlxlPubrgEVdBHpDfwLsAOfGGNeL/B8BeAz4GogFbjHGLMvuFEhJsrOij0nLnhOYo2KtG1Yla6X1aDr5TVoVL1iuf8hKaXyxEoOve2r6W1fjXHAr6Y2i92tWeNuznrTlL2mXqGvS83MITUzh3W/Fj7LO9puo0blaGpUjqZ6pQrEx0YRH+sgPjaKuJgoKlVwULmCnUrRDipVcBATZSc2yk5stJ0KDhsVHDaiHTYqOOxE2cWSulNsQRcRO/A+cDOQDKwWkVnGmK1+pz0MnDTGNBWRwcAbwD3BDptQNZaK0XbO5LiIibLRtHZlmteOo2mdyrRJqErrhHjiK0YV/mLtWlEq7IhAYzlGY9sChrAAgFOmEuvdl7PdNGKnuyE7TAN2mQRyuHCXS47LzeG0sxxOC84Svw6bEGW34bB7/tduE6Jsgt0uRNttLHjmhqB8n3zfM4BzOgG7jTF7AERkMtAP8C/o/YCR3uMpwBgRERPkDnqbTRj/207Ui48hoWostoJjxEELt1IRrqpkcoN9Izew0feYywhHqM5+dx1+NbXZb+pwxFTnCNU9/2uqXbBf/mI43Qan2wW55z8X7SidUTmBFPQE4IBfOxnoXNQ5xhiniKQBNYDj/ieJyHBguLeZISI7LiZ0ENSkQLYQEaq5IHSzhWouCN1soZoLSjVbGrD3Yl8c9Fzy6kW/tHFRTwRS0AvrCCp45R3IORhjPgI+CuB7lioRSSrqLrGVQjUXhG62UM0FoZstVHNB6GYL1VwFBXLdnww09Gs3AA4VdY6IOIB44MJ3L5VSSgVVIAV9NdBMRJqISDQwmPPXHJwFDPUeDwR+DHb/uVJKqQsrtsvF2yf+BPA9nmGL44wxW0RkFJBkjJkFjAU+F5HdeK7MB5dm6CCwvNunCKGaC0I3W6jmgtDNFqq5IHSzhWqufCybKaqUUiq4ImdFG6WUCnNa0JVSKkyEbUEXkd4iskNEdovIc4U8f52IrBURp4gMDLFsfxSRrSKyUUQWiEiR407LONcIEdkkIutFZImItCyLXIFk8ztvoIgYESmTIWYBfGbDRCTF+5mtF5FHyiJXINm859zt/V3bIiJfhkIuEXnH7/PaKVJ2O7IHkK2RiCwUkXXe/3/2KatsATHGhN0Xnpu3vwCXAdHABqBlgXMSgTZ41qAZGGLZegIVvce/B74KkVxV/I7vAOaFymfmPS8O+BlYAXQMhVzAMGBMWf1+lTBbM2AdUM3brh0KuQqc/wc8AzFC5TP7CPi997glsK+sf7YX+grXK3TfcgXGmBzg3HIFPsaYfcaYjYC7sDewONtCY8y5LZVW4Bn7Hwq50v2alShk8phV2bxeAd4EgrMYR/ByWSGQbI8C7xtjTgIYY46FSC5/9wKTyiAXBJbNAFW8x/GcPyfHUuFa0AtbriDBoiwFlTTbw8DcUk3kEVAuEXlcRH7BUzifLINcAWUTkfZAQ2PMt2WUKaBcXnd5/zyfIiINC3m+NASSrTnQXESWisgK76qqoZALAG9XYxPgxzLIBYFlGwkMEZFk4Ds8f0GEjHAt6AEtRWCRgLOJyBCgI/BWqSbyfrtCHits+Yb3jTGXA38G/lrqqTwumE1EbMA7wDNllMf3rQt5rOBnNhtINMa0AeYDE0o9lUcg2Rx4ul1uwHMl/ImIVA2BXOcMBqYYY1ylmMdfINnuBcYbYxoAffDMvwmZOhoyQYIskOUKrBJQNhHpBbwA3GGMyQ6VXH4mA/1LNVGe4rLFAa2ARSKyD+gCzCqDG6PFfmbGmFS/n9/HePYMKAuBLtkx0xiTa4zZC+zAU+CtznXOYMquuwUCy/Yw8DWAMWY5EINn4a7QYHUnfind3HAAe/D8uXbu5sZVRZw7nrK9KVpsNqA9npszzUIsVzO/4754ZgqHRLYC5y+ibG6KBvKZ1fM7HgCsCJXPDOgNTPAe18TT3VDD6lze864A9uGd/BhCn9lcYJj3+Eo8Bb/MMhb732B1gFL84fQBdnoL4wvex0bhueIFuAbPv8iZeHZZ2hJC2eYDR4H13q9ZIZLrX8AWb6aFFyqqZZ2twLllUtAD/Mz+7v3MNng/sxah8pnh6WIYjWdvg03A4FDI5W2PBF4vq8+qBJ9ZS2Cp9+e5HrilrDNe6Eun/iulVJgI1z50pZSKOFrQlVIqTGhBV0qpMKEFXSmlwoQWdKWUChNa0JVSKkxoQVdKqTDx/ypBveekvzbRAAAAAElFTkSuQmCC\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"# Plot the data with with expected distribution, under the hypothesis that it is normally distributed. \n",
"# The blue line is the expected data density, with the same mean and standard deviation, if the data are normal. \n",
"xs = np.arange(d2.min(), d2.max(), 0.01)\n",
"fit = stats.norm.pdf(xs, np.mean(d2), np.std(d2))\n",
"plt.plot(xs, fit, label='Normal Dist.', lw=3)\n",
"plt.hist(d2, 50, density=True, label='Actual Data');\n",
"plt.legend();"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the above plot, we can see _some_ difference from the expected normal distribution. \n",
"\n",
"Overall, though, we might think it looks pretty normal, especially if examined without adding the probability density function. \n",
"\n",
"Let's check this statistically. "
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"# Run normal test on the data\n",
"stat, p_val = normaltest(d2)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"Normaltest p-value is: 1.05e-170 \n",
"\n",
"We have evidence to reject the null hypothesis.\n"
]
}
],
"source": [
"# Check the p-value of the normaltest\n",
"print('\\nNormaltest p-value is: {:1.2e} \\n'.format(p_val))\n",
"\n",
"# With alpha value of 0.05, how should we proceed\n",
"check_p_val(p_val, alpha=0.05)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this example, the plot is ambiguous, but 'normaltest' suggests these data are actually very unlikely to come from a normal distribution. \n",
"\n",
"We happen to know that this is indeed true, as the 'ground truth' of the data is that they were generated from a beta distribution.\n",
"\n",
"Therefore, using this data in statistical tests that expect normally distributed inputs is invalid, since we have violated the assumptions upon which these tests are based. We will have to use different methods to perform statistical comparisons with these data."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Evaluating Different Distributions\n",
"\n",
"So far we have examined checking whether a dataset is normally distributed.\n",
"\n",
"More generally, you can use the Kolmogorov-Smirnov test to check if a set of data has some other distribution (that you specify). \n",
"\n",
"This test is implemented in `scipy`, and you can explore using it. "
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"from scipy.stats import kstest"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"# Check out the documentation for kstest\n",
"kstest?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The `kstest` has a similar form to the `normaltest`, in that we can compare a sample of data to a specified distribution. \n",
"\n",
"The null hypothesis is that the data comes from the specified distribution. \n",
"\n",
"If we reject the null, we can conclude that the data do not conform to the tested distribution."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Statistic: \t0.56 \n",
"P-Value: \t0.00e+00\n",
"\n",
"We have evidence to reject the null hypothesis.\n"
]
}
],
"source": [
"# Let's continue using our data from above, from the beta distribution. \n",
"# We can confirm kstest considers it not normally distributed \n",
"stat, p_val = kstest(d2, 'norm')\n",
"print('Statistic: \\t{:1.2f} \\nP-Value: \\t{:1.2e}\\n'.format(stat, p_val))\n",
"check_p_val(p_val, alpha=0.05)"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Statistic: \t0.00 \n",
"P-Value: \t6.00e-01\n",
"\n",
"We do not have evidence to reject the null hypothesis.\n"
]
}
],
"source": [
"# Now compare the data to the beta distribution.\n",
"# Note that in this case, we have to specify some parameters for the beta \n",
"# distribution we are testing against, so we will use the simulation parameters\n",
"stat, p_val = kstest(d2, 'beta', [7, 10])\n",
"print('Statistic: \\t{:1.2f} \\nP-Value: \\t{:1.2e}\\n'.format(stat, p_val))\n",
"check_p_val(p_val, alpha=0.05)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this case, the data sample is consistent with being from a beta distribution!"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.4"
}
},
"nbformat": 4,
"nbformat_minor": 2
}