{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Pearson's Product-Moment Correlation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this tutorial, we explore\n", "\n", "- The theory behind the Pearson test statistic and p-value\n", "- The features of the implementation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Theory" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following description is adapted from [[1]](https://arxiv.org/abs/1907.02088):\n", "\n", "Pearson's product-moment correlation is a measure of the linear correlation between two univariate random variables [[2]](https://bayes.wustl.edu/Manual/Pearson_1895.pdf). Given sample data $\\mathbf{x}$ and $\\mathbf{y}$, the sample Pearson correlation is\n", "\n", "$$\\mathrm{Pearson}_n (\\mathbf{x}, \\mathbf{y}) = \\frac{\\hat{\\mathrm{cov}} (\\mathbf{x}, \\mathbf{y})}{\\hat{\\sigma}_{\\mathbf{x}} \\hat{\\sigma}_{\\mathbf{y}}},$$\n", "\n", "where $\\hat{\\mathrm{cov}} \\left( \\mathbf{x}, \\mathbf{y} \\right)$ is the sample covariance, $\\hat{\\sigma}_{\\mathbf{x}}$ and $\\hat{\\sigma}_{\\mathbf{y}}$ are the sample standard deviations of $\\mathbf{x}$ and $\\mathbf{y}$ respectively.\n", "\n", "This implementation wraps scipy.stats.pearsonr [[3]](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.pearsonr.html) to conform to the mgcpy API." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Using Pearson's" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Before delving straight into function calls, let's first import some useful functions, to ensure consistency in these examples, we set the seed:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "import numpy as np\n", "import matplotlib.pyplot as plt; plt.style.use('classic')\n", "import seaborn as sns; sns.set(style=\"white\")\n", "\n", "from mgcpy.independence_tests.rv_corr import RVCorr\n", "from mgcpy.benchmarks import simulations as sims\n", "\n", "np.random.seed(12345678)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To start, let's simulate some linear data:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhEAAAIlCAYAAABrdaqpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deZhkdXn3//fMwAzOzrArJLjE+0FQgaBBlGV8CD/RbBqXRAOiQYxRRAGNRhQwxLgAsgdRgkGIQfQxRmNcogyKCBgXRIRbRVBEBphpppm1Z+n+/XGqhpqa6u7q01VdVV3v13X1Vd2nTp2+qyimPv1dZ4yMjCBJkjRRMztdgCRJ6k2GCEmSVIohQpIklWKIkCRJpRgiJElSKYYISQIiYkana5B6zQ6dLkCa7iJiGbBTZh46xjlnAWcCT8jMDVNUWmkRsTPwd8CfAfsCQ8BdwDXA5Zm5uXLeUcANwLGZ+ZUpqOuTwIsyc88JPGYG8F5gA/DhyrGz6KH/HlKn2BIhdYdPAM+j+DDuahHxBODbwCuBi4AXA68GbgY+ShEkqn5A8by+O8VlTsQc4Gxgbs2xnvnvIXWSLRFSF8jM3wC/6XQdTXo5sD/w7Mz8cc3x/4qIVcDZEfHBzPxRZj4G3NKRKiehx/57SB1jiJC6QH3zeaVZfl/gCuAM4KnAfcAHMvNfax43u3L/ccATgV8Bl2XmBXXXPwF4E8WH/47ALyvnXVy5f1/gXuBU4PXAPsA5mXlug3KrXQWNWjKvBDYCqyvXPYqa7oxKHVcBh1K0WhwEPAScBXwFuBg4FngMuCoz31NX35sy8/LRXrf6YiJiJvB24Hjg94AZwN3ABzPzuprrApwZEWdm5oxG142Ig4B/AJ5L0WrxXeC9mXlLXY2vAv4UeAkwC/hv4K2ZubzB6yX1NLszpO51IEUz+z8BfwT8GvhkRDyz5pzrgdOAj1XO+SxwXkT8U/WEiHgj8C/A14E/pmhJuA+4qPIhX+uDlWsdD3xplLq+AmwGvhoR/xgRR1S6OMjMBzLzg5l5zzjP7XPApyr13EfRfbAM+BnwMuAbwN9HxMvHuc54/rHydTVFt8tfUYScayPiqcCDwJGVc6+k6MLYTkQcCdwKLAT+huL1mQvcGBGH153+MeBRitf5XZXneOkkn4fUlWyJkLrXIuDwzLwDICKSoqXhT4A7IuKFle9fn5lXVR7z9YjYALwvIi6tNMs/DbgwM8+oXjgivgOsBJZSfHhXfTEzLxmrqMy8IyL+HLgc+PvK16aIuBW4DrgiMzeO89zOy8x/rtSyqVLD92paHm4AXgE8nyIYlbUPcGZmnlc9EBH3At8HjsjMqyLitspdv6m2KjTwIeB+4Ojqc4uI/6Jo1TgX+IOac7+ZmW+pfP8/EfH7wF9FxIzMdLMiTSuGCKl7ra4GiIpqH/28yu0fVm7/MyJq/1/+AkULxv8F/jUz3wEQEQuBp1OEit+vnDu77nfeQRMy8z8j4r+BIyq/53CKD9IXAG+MiKMyc+UYl/hOzfcPVW63foBn5qaIGAR2bqaeMer8K4CI2IXHn/vSyt31z72hiJhH0YXx4dpwlJlDEfEZ4B0RMb/mId+pu8RvKLqQdqRoBZGmDUOE1L3W1f6QmcMRAY93Q+5auV0xyuOfBBARTwb+GTgG2ELRZVD9oKtfG+EhmpSZmyi6Hb5R+T0LgXdTNOH/HfDOMR7+WINja+t+nvRf7RFxMHAJj8+0uIvHg1Kz60IsrpzbaEzDg5X7FtYcW1d3znDl1u5jTTu+qaXetYriL9vnAs9p8PXJyhoIXwJ+l6KVYF5m7g+8rewvjYibI+K6+uOZ+VhmvpsipOxf9vqjqAaKWXXHF4z2gIhYQDF+Y5hifMm8zDyIYtzHRKyq/P5Ga088sXLfWK0u0rRliJB61zKKJvm5mfm/1S+Kv5zPAfYCdgOeAXwyM2+uaY5/ceW2zL8B9wB/GhHPqL8jIhYDuwM/3u5Rk1Ntudin7nj9oMZa+1E8/4sz8/bM3FI5Xv/ct2z3yBqZuRa4DXhFROxYPR4RcyjWyrg1M11PQn3J7gxpauwZEY3++n8wM7f7q75J/00xffL6iPhH4HaKD85zKJrZf1Lpt78X+JvK7SMUH7x/R/EX9LyGVx7b3wNHATdHxKUUC0+tq/zutwEDwHmjPrqEzHw0Im6ieB53AQ9QTEV98hgPuxsYBN4VEesrNb4YqA56nFe59qaIWAMcFhFHVJ5PvXcDXwO+EREfrRw7laLL6ITJPDepl9kSIU2N36VYF6H+6+1lL5iZwxTTOq+qXOdrFB/w1wNLa/46/hOKdSE+QTG18iXAiRQh5IgSv/d+ivUdPk6xHsL1FNNHTwe+DDw3M0cbpzEZr6UYy3EZ8O8UYeXvxqjzMYrnPgR8uvJ1MMVr9lMen9oJxfoPz6F4TepbO8jMG4AXUkxtvQb4JMUYjiMyc9mknpXUw2aMjDjjSJIkTZwtEZIkqRRDhCRJKsUQIUmSSjFESJKkUgwRkiSpFEOEJEkqxRAhSZJKMURIkqRSDBGSJKkUQ4QkSSrFECFJkkoxREiSpFIMEZIkqRRDhCRJKsUQIUmSSjFESJKkUgwRkiSpFEOEJEkqxRAhSZJKMURIkqRSDBGSJKkUQ4QkSSrFECFJkkoxREiSpFIMEZIkqRRDhCRJKmWHThcwnoj4GDArM08c45zrgZfXHf5GZh7d1uIkSepjXRsiImIGcDZwEnDlOKcfALwL+NeaY0NtKk2SJNGlISIinkIRHA4Afj3OubOBpwG3ZebyKShPkiTRvWMingf8EngmcO845+5HEYbuandRkiTpcV3ZEpGZ1wLXAkTEeKcfAGwEzo6IY4H1wPXAOZm5oZ11SpLUz7oyREzQ/sAMIIFLKFovzgf2AV7bzAUiYg7wHOBBYEt7ypQkqWvMAvYCvpeZpccQTocQcQZwbmYOVH6+IyK2AP8eEadm5somrvEc4Nttq1CSpO50OHBT2Qf3fIjIzGFgoO7wHZXbfYBmQsSDANdeey177rlnC6uTJKn7LF++nNe85jVQ+fwrq+dDRER8BtgxM19ac/gQiimev2jyMlsA9txzT/bee+8WVyhJUteaVBd+z4WIypTOJcBAZm4EPkul6wL4AnAQcC5FF8eazlUqSdL01q1TPMdyGEXzy2EAmfkZ4ATgdcBPgPOAC4H3dag+SZL6Qte3RGTmUXU/L6OYjVF77Grg6qmrSpIk9WJLhCRJ6gKGCEmSVIohQpIklWKIkCRJpRgiJElSKYYISZJUiiFCkiSVYoiQJEmlGCIkSVIphghJklSKIUKSJJViiJAkSaUYIiRJUimGCEmSVIohQpIklWKIkCRJpRgiJElSKYYISZJUiiFCkiSVYoiQJEmlGCIkSVIphghJklSKIUKSJJViiJAkSaUYIiRJUimGCEmSVIohQpIklWKIkCRJpRgiJElSKYYISZJUiiFCkiSVYoiQJEmlGCIkSVIphghJklSKIUKSJJViiJAkSaUYIiRJUimGCEmSVIohQpIklWKIkCRJpRgiJElSKYYISZJUiiFCkiSVYoiQJEmlGCIkSVIphghJklSKIUKSJJViiJAkSaUYIiRJUimGCEmSVIohQpIklWKIkCRJpRgiJElSKYYISZJUiiFCkiSVYoiQJEmlGCIkSVIphghJklSKIUKSJJViiJAkSaUYIiRJUimGCEmSVIohQpIklbJDpwsYT0R8DJiVmSeOcc4hwIXAQcADwD9k5tVTVKIkSX2pa1siImJGRLwfOGmc83YDvgr8ADgYuAi4MiKOaX+VkiT1r65siYiIpwBXAgcAvx7n9BOBQeCUzBwG7o6Ig4HTga+1tVBJkvpYt7ZEPA/4JfBM4N5xzj0c+FYlQFQtA54fEd36/CRJ6nld2RKRmdcC1wJExHin7w38sO7Yb4G5wBJgRavrkyRJ3dsSMRFzgQ11x4YqtztNcS2SJPWN6RAi1gNz6o5Vf147xbVIktQ3pkOIuB/Yq+7YE4E1FAMuJUlSG0yHEHETcEREzKg5thT4Tt1gS0nSNLBleISBwQ3cde9KBgY3MDw80umS+lZXDqwcS0TMphgwOZCZGymmgr4TuDwiLgCOBl4NvKhzVUqS2mVw9RAnn3cDj63dyMJ5s7n4tKUsWeQQuE7oxZaIw4AHK7dk5kMUgeEgilkabwGOz8xvdqxCSVLbPDSwlsfWbgTgsbUbeWjA4W+d0vUtEZl5VN3Py4AZdcduAZ47dVVJkjpljyXzWDhv9taWiD2WzOt0SX2r60OEJEm1Fi+Yw8WnLeWhgbXssWQeixfUT9DTVDFESJJ6ysyZM1iyaCfHQXSBXhwTIUnqQs6a6D+2REiSWsJZE/3HlghJUtPGam1w1kT/sSVCktS0sVobnDXRfwwRkqSmNWptqIYIZ030H0OEJKlpY7U2tHLWxJbhEQZXD20TSGbOnDH+AzWlDBGSpKZNVWuDgzR7gyFCktS0qVqjYaxuE3UPZ2dIkrpOtdsEcJBmF7MlQpLUdRyk2RsMEZKkruPS1r3B7gxJklSKIUKS1JB7YWg8dmdIkhpymqXGY0uEJKkh98LQeAwRkqSGnGap8didIUlqyGmWGo8hQpLUkNMsNR67MyRJUimGCEmSVIohQpIklWKIkCRNSLOLULlY1fTnwEpJElB86A+uHtpmNsbMmTO2O6/ZRahcrGr6M0RIkoDmP/QbLUI1mfPUu+zOkCQBza9Q2ewiVC5WNf3ZEiFJAmD3JXNZOG/21paIPZbMbXhes4tQuVjV9GeIkCQBMGvmDM5+w/N4+NF17L7z3IbjIaD5RahcrGr6M0RIkgB4cMVa3nXpTcyfO5s16zbywTe/gMULDAAanSFCkvpQo5kYeyyZx/y5td0ZjmHQ2AwRktSHGs3EcAyDJsoQIUl9aLTpl45h0EQ4xVOS+pDTL9UKtkRIUh/q5q6LZlfOVOcZIiSpD3Xz9MvB1UO85/Lv8NtH1jB/rstldzNDhCT1gH7563zL8Aibtgxz3LH7sdviJ3D55+9wuewuZoiQpB7QL5tZDa4e4tQLbtz6PM9+w/NYsnD6Pc/pwoGVktQDmt3XotfVP8+hTVu6aryGtmWIkKQe0OuzKbYMjzAwuIG77l3JwOAGhodHGp5X/zz32mXetOy2mS7szpCkHtDNsyma0Wx3TK8/z35jiJCkHtDNsymaMdriVvV6/Xn2G7szJElt1+vdMWrMlghJUtvZTTE9GSIkSW1nN8X0ZHeGJEkqxRAhSZJKMURIkqRSDBGSpKY0u2CU+ocDKyVJTemX/TvUPFsiJElN6Zf9O9Q8Q4QkdVC3dBE0U4cLRqme3RmS1EHd0kXQTB0uGKV6hghJ6qBm95TohjpcMEr17M6QpA5qZRfBZLpG7KpQGbZESFIHtaKLYMvwCIOrh9i0ZZhTL7ixVNeIXRUqwxAhSR3Uii6CwdVDvOfy73Dcsftt7ZJYs24jK1ata/q6dlWoDEOEJHWBamtCbUvAzJkzmnrsQwNr+e0ja9ht8RNYOG82e+4yjze97FkMbdrCwOCGCV1LmghDhCR1gcnM0thjyTzmz53N5Z+/g7Pf8DwWzJtdultDmghDhCR1gcnM0qgdz7DLop1YvrI7Znxo+jNESFIXqM6OqLYeTGR2RP14hpERSl9LmghDhCR1gWZnRzQzdsKZFpoqXRkiImIWcA5wArAA+Arw5sx8aJTzrwdeXnf4G5l5dDvrlKRWaXZ2RDNjJ5xpoanSrYtNnQW8FjgeOALYG/jcGOcfALwL2Kvm6xXtLVGSttfuvTDcBEvdpOtaIiJiNnAK8NbM/Hrl2F8A90bEYZl5c4PznwbclpnLp7xgSarR7r0wJjN2Qmq1rgsRwIEUXRjLqgcy876IuA84HLi57vz9KJ7HXVNTniSNrt17YTjeQd2kG0PE3pXbB+qO/xbYp8H5BwAbgbMj4lhgPXA9cE5mbmhblZLUQLtbChzvoG7SjSFiLjCcmZvqjg8Bjf6v2R+YASRwCfBM4HyKwPHaNtYpSdsZr6VgMitTSt2mG0PEemBmROyQmZtrjs8BGo0gOgM4NzMHKj/fERFbgH+PiFMzc2Wb65WkrcZrKWj3mAlpKnVjiLi/crtXzfcAT2T7Lg4ycxgYqDt8R+V2H8AQIanjtgyP8NiaIR54ZLWrSWra6MYpnrcDq4EjqwciYl9gX+Bb9SdHxGci4vN1hw+h6P74RduqlKQJqLZA7DR7BxbOmw3Q0jET7Z5aKjXSdS0RmTkUEZcB50bECuBh4DLgxsy8pTKlcwkwkJkbgc9S6boAvgAcBJxL0cWxpjPPQpK29dDAWgbXbOTyz9/BWSceyoaNm3nSbgsazq4oM27CbhJ1Qje2REAxzuFa4BrgBuBXPL4i5WHAg5VbMvMzFCtbvg74CXAecCHwvimtWJLGUJ218bNfP8rZV97C3rsvYMminRqGg2ogeOclN3HyeTewavXQuNd3ESp1Qte1RABUBlSeVvmqv28ZxWyM2mNXA1dPSXGSVEL9rI1F80df36HMWhMuQqVO6MoQIUm9arSuiIms71AmELgIlTrBECFJLdSKsQllAoGLUKkTDBGS1EKtWPZ6vEDgglXqFoYISWqhqRib4EwMdQtDhCS10GTHJjTTytDuTb6kZhkiJKmFJjs2oZlWBmdiqFsYIiSphRq1JIxA02MYxmtl2DI8wsyZcP7bjmTFqvXstctcZ2KoYwwRktRCjVoSgG2OXXTaUcxgRsNQMV4rQ/X6a9ZtZK9d5/OBNz3fQZXqGEOEJLVQfUvCY2uHWLth8zbHlq9cx8Wf+REPPLJmuy6L8cZU1F7/gUfWOB5CHdWty15LUk+qtiQALJo/m/lzZzNnx1nbbLq16+In8OCKYmuf+iWqq2Mq9nvyLg2Xxa69vuMh1Gm2REhSC9W2JOy5yzyWr1zLJ/7zTs468VAeWbWepzxpEUObNjN/brmBka5MqW5iiJCkFqqfnTEyAstXruX0i77FXrvO57RXH8zHPn8H57/tSAYG129tWRgY3NDUwEtXplQ3MURIUhtVWw4eXLmWOTvO4p//349ZvnItO86ayX5P3gWAgcENLh6lnmSIkKQarV5SutpysHjBHFatHuLEP9l/u24IF49SrzJESFKNdi0pPVY3hItHqVcZIiSpRidaBRwsqV5liJCkGp1oFXCwpHqVIUKSaky0VcBtudXPDBGSVGOirQJuy61+5oqVkvrGluERBgY3cNe9KxkY3MDw8Mikr9loDIXUL2yJkNQ32tFq4MwK9TNDhKS+0Y6ZF86sUD8zREjqG+1oNSg7s8IBmZoODBGS+kY3tRo4IFPTgSFCUt/o5HoM9S0Pq9ZscKlr9TxDhCRNgfqWhwtPPYpF82czuMYBmepdhghJmgL1gzofeXQdF5+2lOUrO9+1IpVliJCkKdBoUOfOC3di54V2Yah3GSIkaQp006BOqVUMEZI0BdxkS9PRqMteR8Sbp7IQSZLUW8baO+PCiPifiNhnyqqRJEk9Y6wQ8QJgD+AnEfHXU1SPJEnqEaOGiMy8BTgI+AhwSUR8OSKeNGWVSZKkrjbmwMrM3AycExGfBi4E7oiIy4B1ded9oH0lSlL3ce8LqfnZGb8BvgccA7weGKq5bwQwREjqK+59ITURIiLixcClwC7A2zPz0rZXJUldrh3biku9ZtQQERG7AxcBrwD+BzgpM381VYVJUjdrx7biUq8ZqyXi7srtiZl51VQUI0m9whUopbFDxI3AmzJz+VQVI0m9whUopTFCRGa+dCoLkaRu4cwLqTnunSFJdZx5ITVnrBUrJakvNZp5IWl7hghJqlOdeTFzBuy9+3xnXkijsDtDkuosXjCHS09fytDmYVasWg+MMDw84rgIqc5Y60R8s9mLZOYLW1OOJHXezJkzGB6BUy+40XER0hjG6s64p+brQeAoYB5wB/B9YAZwOPDT9pYoSVPPcRHS+Maa4vmG6vcR8SngQ5n57tpzIuIM4LntK0+SOsMVKaXxNTsm4qUU24LXuw74+9aVI0ndwRUppfE1GyIepui6+Hnd8WOB+1takSR1AVeklMbXbIg4D7gsIg4FfkAxHuIw4FXACe0pTZIkdbOmQkRmXhoRjwFvBl4NjAA/Al6RmV9oY32SJKlLNb1ORGZ+CvhUG2uRJEk9pOkQERFPA04D/g/wV8CfAj/NzGXtKU2S2sdNtqTJa2rZ64j4A+B24KkUYyHmAPsDX4+IP2pfeZLUHtVNtt55yU2cfN4NrFo91OmSpJ7T7N4ZHwQ+mJnHABsBMvPNwIeBs9pTmiS1j4tJSZPXbIg4GPh0g+NXAvu1rhxJmhrVxaQAF5OSSmp2TMRaYHfgF3XHnw4MtrQiSZoCLiYlTV6zIeLTwEcj4gSK6Z07RcQLgUuAz7apNklqGxeTkiav2e6MdwP3AXcC8yk24fo68L3KfZLUFbYMjzAwuIG77l3JwOAGhodHOl2SNG01u9jURuBVEfEe4ECKwZV3ZuY97SxOkiaqOuvCLbyl9msqRETEL4FDMvMX1IyLiIi9gNszc/c21SdJ2xlrjYdGsy4MEVJ7jBoiIuLFwCGVH/cF3hURa+pOe/pY1ygrImYB51Dsy7EA+Arw5sx8aJTzDwEupNhp9AHgHzLz6lbXJak9RgsFtcf32mUewyMjPDSwjt12nsuF1/2AH/1sxXatDW7hLU2dsQLAvcAFFJttAbwc2FJz/wiwGji5DXWdBbwWOB5YCVwGfA54Qf2JEbEb8FXg34C/Bv4QuDIilmfm19pQm6QWG60Lonp8zbqNnHfKkZz58e9uPeesEw/lxz//1natDc66kKbOqCEiM++iaGkgIm4AXpaZj7a7oIiYDZwCvDUzv1459hfAvRFxWGbeXPeQEymmmZ6SmcPA3RFxMHA6YIiQesBoXRDV44vmz+bhR9dtc86KwfXMn1us81Db2uCsC2nqNDU7IzOXNgoQETE7Ip7f4poOpOjCWFbz+++jmB1yeIPzDwe+VQkQVcuA50dEs7NPJHXQaAs/VY+vXruR3Xeeu805v7fPzpzxuudy8WlLbW2QOqTZgZW/D3wceCaNg8esFta0d+X2gbrjvwX2GeX8HzY4dy6wBFjRwtoktcFoXRC1x3ddvBMXn3YUDw2s23rOrouf0OHKpf7W7KDIC4H1wEkU4xNOAZ5cuT2+xTXNBYYzc1Pd8SGgUfvkXGBDg3MZ5XxJLdSK3TBH64JodHzJom2Dg7txSp3TbIg4CDgiM78fEScBmZlXRMRvgTfR2lUr1wMzI2KHzNxcc3wOxfLbjc6vb8us/uyOOlKbdXpdhk7/fqmfNTtmYAbwSOX7n1N0awB8EXh2i2u6v3K7V93xJ7J9F0f1/EbnrsF9PaS26/RumJ3+/VI/azZE/AR4ceX7nwLVwZR70NrxEAC3U0wdPbJ6ICL2pVir4lsNzr8JOCIiatsvlwLfqRtsKakNOr0bZqd/v9TPmu3O+BBwXURsodiM68yI+A+KVogbWllQZg5FxGXAuRGxAniYYhzGjZl5S2UK6BJgoLIc95XAO4HLI+IC4Gjg1cCLWlmXpMbasS7DRMY5uC6E1DnNTvH8HHAocFtm/oqiVWIj8GXgDW2o6wzgWuAaipDyK4rFrgAOAx6s3FJZxfJFFOM2fgi8BTg+M7/Zhrok1akOftzvybuwZNFOLRnUWB3n8M5LbuLk825g1eqhUc9tx++X1Jyml6zOzP+t+f4GWtwCUfe7NgOnVb7q71vG46toVo/dAjy3XfVImlrufyH1hrH2zmh6tcfMPKY15UiaCt0+LdL9L6TeMFZLRKOZEJKmgW6fFuk4B6k3jLV3xuumshBJU6fbuwvc/0LqDc0ue/3qse7PzH9rTTmSpoLdBZJaodmBldeMcnwD8BuKbbgl9Qi7CyS1QlMhIjO3mQoaEbMotgn/Z+BjbahLUht1orug2wdzSpq4pqd41srMLcBdEXEq8BmKBagkaVTdPphT0sQ1u+z1aDZT7FMhSWNyjwtp+pnMwMqFFFuD39rSiiRNSw7mlKafyQys3AR8F/jb1pUjabpyMKc0/ZQaWCmpe3XrAEbXfpCmn1IDKyV1LwcwSpoqzY6J+H3gEuAAYLs2yMyc3eK6JJXU7atRSpo+mm2J+ATF1t/vANa3rxxJk+UARklTpdkQEcBzMvPOdhYjafIcwChpqjQbIn4A/A5giJC63GgDGLt1wKWk3tVsiDgJ+HxEPAf4JTBce6cbcEndzwGXklqt2RDxMuD3gLMa3DeCG3BJXc8Bl5JardkQ8TbgDOCCzFzXxnoktYkDLiW1WrMhYhbwaQOE1LsccCmp1ZpdifITwN+0sxBJ7VUdcLnfk3dhyaKdHFQpadKabYlYBBwfEX8J3EOxb8ZWmXlMqwuTJEndrdkQsSPw6XYWIkmSekuzG3C9rt2FSJKk3jJqiIiIVwOfzcyNle9HM5KZtlJIbeZiUZK6zVgtEdcA/wM8XPl+NCPY1SG1XSsWizKISGqlUUNEZs5s9L2kzmjFYlGPrRnireffwOAaV62UNHkTDgcRsUNEHBwRT2pHQZIaqy4WBUx4sagtwyMMDG7gNw+v5sy/PpSn/87OW4OIJJU15sDKiDgOOAV4WWb+OiKeAXwZ2AcYiYhPAm/MzC1tr1Tqc5NZLKq+K+SsEw/l7CtvcdVKSZMyaktERLwS+CTwE6D658qngIXAi4DDgEMplsSW1GaTWSyqvitkw8bNXHzaUletlDQpY3VnvBU4IzNPyMyVEfFs4CDg4sz8embeBrwXcPqn1AWqXRZ33buSgcENDA+PbL2vvivkSbstYOeFrlopaXLG6s54FnBizc9HU8zE+GLNsR8DT21DXZImaKzZG7VdIXvtMo/hkRHuunelMzQkTcpYIWImsLHm5yOAQeD7NceeAGxoQ12SJmis2RvVrpAli3ZiYHADJ5+3bFJTRSUJxu7OuBN4PkBELAT+L/C1zBypOefPKcZMSOqwZmdvNAobklTGWC0RlwIXR8SzKMLEE4ALACJid+DVwLuAN7a7SEnja3b2RjVsVFsinKEhqayxFpu6OiJ2Ak4CtgCvysxbKnefSTFe4sOZeXX7y5Q0nuMPfXAAABh4SURBVNoui7FMZqqoJNUac52IzLwCuKLBXf8EvC8zV7alKkltUxs2tgyPsMplsCWV1OxW4NvIzN+0uhBJU68V+3FI6l+lQoSk6aF+kOWq1Ru2HrdlQtJ4DBFSH6sdZLlo/mwWzp9jy4SkphkipD5WO8hyz13msXzl5HcKldQ/DBFSH6uf0TEygtM/JTXNECFpK6d/SpoIQ4SkrZpda0KSYOxlryVNwFi7aHbTNSWpVWyJkFqkHWsuuI6DpG5mS4TUIu3Y2MrNsiR1M0OE1CLN7qLZ6WtKUqvYnSG1SDtmNjhbQlI3M0RILdKOmQ3OlpDUzezOkCbJGRSS+pUtEdIkOYNCUr+yJUKqUaZVwRkUkvqVLRFSjTKtCrU7YTqDQlI/MURINRq1KowXIpxBIalfGSKkGmVaFVo9g2LL8AiDq4e2CSUzZ85oybUlqZUMEVKNVrcqlAkEDtSU1CsMEVKNVrcqjBYIxgoXZbpUJKkTDBFSG40WCMZqbXCgpqReYYiQ2mi0QDBWa4MDNSX1CkOE1EYL583mwlOP4sEVa9hr1/lbN9Maq7XBpa4l9QpDhNRGj63dyNs+uoyREXjq3gt56ysP5pFH17HHkrlc+o6lPLjC1gZJvcsQIbXRQwNrGVyzkZkz4LgXPYNTzl+2zTiI/Z68S6dLlKTSui5ERMTuwCXAMcBG4CrgPZm5eYzHPAzsVnf4vZl5TtsKlZpQ7baYMQNWDK531oWkaaXrQgTwOWAEOBJ4EvBJYDPwnkYnR8QeFAHiCODnNXetbmuVUhNqB0nutvPcbcZBLFn0BAYGN7iYlKSe1VUhIiKeB7wAeEpm3gvcHhHvAC6OiPdn5lCDhx1AETJuzcyNU1iuBIy9oFTtIMnh4REuPm0pD65cy5wdZ/HhT/0vy1eudTEpST2rq0IEcDjwq0qAqFoGLAAOBG5t8JgDgHsMEOqUVas38Nbzasc6HMWSRU/Y7rxqoFixah2nXXgj1Q1C7daQ1Ku6bSvwvYEH6o79tnK7zyiPOQDYHBFfiojlEfH9iDiubRWqrzXaKnz5ynXbjHV4cOW6Ma+x6+K5zJ9bTPV0MSlJvWxKWyIiYl/g3lHuHgKuATbUHszMTRExAoz2p9r+wC7AeynGTRwLXBURO2TmVa2oW6pqtNLkroufsM1Yh10Xb98KUcvFpCRNF1PdnfEAsN8o9w0DJwPb/IsaETsCM4C1ozxuKTA7M6sDKW+PiN8FTqWY2SG1TKOVJp+42zzOfsPzePjRdey+81zm7Dh2A5+LSUmaLqY0RGTmJuDu0e6PiPuBF9cdfmLltr6bo3rNIYpWjFp3AH9ZskxpVI1Wmlwwdw5btsDw8DBLFu7Ewnm2LEjqD902sPIm4EMRsU9m3l85tpRiuuaP6k+OiB0oukfOz8yP1tx1CHBnu4tV/2m0jLUtC5L6VbeFiO8CtwDXRcRbgD2AD1GEhI0AETEfmJ+ZyzNzc0R8ETgjIu4Bfgr8GXAc8JKOPANNa7XLWM+YARed6vRMSf2rq2ZnZOYI8FLgIeDbFGMargTeX3Pa6cCDNT+/HbgcuIii9eE44JWZ+bWpqFn9pbqM9WNrNzK4phgTIUn9qttaIsjM5RRBYrT7zwLOqvl5iGJWRsMVLaVWGmv3TUnqN10XIqRuVp2euWLVOnZdPNfpmZL6Wld1Z0jdrrLIJFuqy01OQKOFqiSpl9kSIU1Ao8Wmmh1YOZnHSlI3siVCmoBGi01NxWMlqRsZIqQJqA6shInvezGZx0pSN7I7Q5qAyex74Z4ZkqYbQ4Q0AZNZndKVLSVNN3ZnSJKkUgwR6gutnl7pdE1JsjtDfaLV0yudrilJtkSoT7R6eqXTNSXJEKE+0erplU7XlCS7M9QnWj290umakmSIUJ9o9fRKp2tKkt0Z0pichSFJo7MlQhpD7SyMA5++K6e86mAeeXTd1i6MmTNndLpESeoYQ4Q0huosjJkz4Phjn8Ep5y9zWqckVdidIY2hOgtjwbzZrBhc77ROSaphS4Q0htpZGLvtPJeF82ZvbYlwWqekfmeIkMZQOwtjeHjEaZ2SVMMQITXJaZ2StC3HREiSpFIMEeoY12CQpN5md4Y6xp0wJam32RKhjnEnTEnqbYYIdYw7YUpSb7M7Qx3jTpiS1NsMEeoYp0xKUm+zO0OSJJViiJDG4DRUSRqd3RnqWVuGRxhcPbTNmIpWb83tNFRJGp0hQj1rKj7gG01DNURIUsHuDPWkLcMjPLiy/etMOA1VkkZnS4S60nhdFYOrh5iz46y2b83tNFRJGp0hQl1pvK6KhwbW8on/vJOzTjyUR1at56l7L27LB7zTUCVpdHZnqCuNtyT2HkvmsXzlWk6/6Ftc/eW72HHWzJYPqpQkjc2WCJXS7pkR1bEIo3VV2M0gSZ1niFAprZoZMVoYGS8k2M0gSZ1niFAprZr6OFoYMSRIUvdzTIRKadXUx/owsnxgLY8+5gqRktQLbIlQKa0ak1A/9mH3nedy8nk3MLjGFSIlqdsZIlTKZLobasdB7LXrPC4+7SgeGljHHkvmsWrNBgbXuEKkJPUCQ4SmXKNxEPs9eZet97d7ASlJUmsYIjTlxhqU6dRNSeodhghNubHWgHBWhiT1DkOEppytDZI0PRgiNOVsbZCk6cF1IiRJUimGCEmSVIohQh2zZXiEgUFXp5SkXuWYCHVMqzbxkiR1hi0R6phG60VIknqHIUId06pNvCRJnWF3hjrG9SIkqbcZItQxrhchSb3N7gxJklSKIUKSJJViiJAkSaUYIiRJUimGCEmSVIohQpIklWKIkCRJpXTtOhERMQe4DfhIZl4zzrmvAd4H/A5wO3ByZn6v/VVKktS/urIlIiIWAJ8HntXEuUcD/wKcBxwM3AF8LSJ2a2uRkiT1ua4LEZVQ8CNgjyYf8g7g05l5RWbeBbwRGADe0KYSJUkSXRgigBdTtCwcNt6JETETeD6wrHosM4eBbwGHt6k+NbBleISBwQ387FcDDAxuYHh4pNMlSZLarOvGRGTmqdXvI2K80xcD84AH6o7/FnhOaytT1ZbhEQZXD22zcdbqNUMMPLaBFYMbmDlzJjNnwuIF7okhSdPZlIaIiNgXuHeUu4cyc6KfOnMrtxvqrwX4CdYmg6uHOPm8G3hs7UYWzpvNxactZfOWYc78+He3Hjv/bUd2ukxJUptNdUvEA8B+o9w3XOJ66yu39XtIzwHWlriemvDQwFoeW7sRgD13mcfm4eFtjj22diMrVq1njyVzx7qMJKnHTWmIyMxNwN0tvOQARVjYq+74E9m+i0MtsseSeSycN5t16zfylpc/m9MuvJH3vf5QFs6bvbUlYq9dDBCSNN113ZiIicjMkYi4GTgS+BRsHWx5BPDxTtY2nS1eMIdLT1/Kxs1FC8T7Xn8o/3XzvZx14qFs2LiZJ+22gMUL6huHJEnTTc+FiIiYD8zPzOWVQ+cDX4yIHwLfBE4FFgGf6FCJ097MmTMYHoG3X3Dj1paHs048lPf/yy1cdOpSdl7ocBRJ6gfdOMVzPKcDD1Z/yMyvACcBpwE/AJ4BHJOZKzpTXn+oHwOxYeNmLjp1KYvm2wIhSf2iq1siMnNGg2NnAWfVHbsKuGpqqhI8Pi6i2hLxpN0W2AIhSX2mq0OEutfiBXO4+LSl26wVIUnqL4YIlTJz5gyWLNqJJYtsfZCkftWLYyIkSVIXMERIkqRSDBGSJKkUQ4QkSSrFECFJkkoxREiSpFIMEZIkqRRDhCRJKsUQIUmSSjFESJKkUgwRkiSpFEOEJEkqxRAhSZJKMURIkqRSDBGSJKkUQ4QkSSrFECFJkkoxREiSpFIMEZIkqRRDhCRJKsUQIUmSSjFESJKkUgwRkiSpFEOEJEkqxRAhSZJKMURIkqRSDBGSJKkUQ4QkSSrFECFJkkoxREiSpFIMEZIkqRRDhCRJKsUQIUmSSjFESJKkUgwRkiSpFEOEJEkqxRAhSZJKMURIkqRSDBGSJKkUQ4QkSSrFECFJkkoxREiSpFIMEZIkqRRDhCRJKsUQIUmSSjFESJKkUgwRkiSpFEOEJEkqxRAhSZJKMURIkqRSDBGSJKkUQ4QkSSrFECFJkkoxREiSpFIMEZIkqRRDhCRJKsUQIUmSSjFESJKkUgwRkiSpFEOEJEkqxRAhSZJK2aHTBYwmIuYAtwEfycxrxjn3YWC3usPvzcxz2lWfJEn9ritDREQsAK4DntXEuXtQBIgjgJ/X3LW6PdVJkiTowhAREUcDHwNWNfmQA4DNwK2ZubFthUmSpG1045iIFwP/AhzW5PkHAPcYICRJmlpd1xKRmadWv4+IZh5yALA5Ir4EHAI8AFyQmZ9qT4WSJAmmOERExL7AvaPcPZSZO5W47P7ALsB7gfcAxwJXRcQOmXlVk9eYBbB8+fISv16SpN5S83k3azLXmeqWiAeA/Ua5b7jkNZcCszOzOpDy9oj4XeBUoNkQsRfAa17zmpIlSJLUk/YC7in74CkNEZm5Cbi7xdccAobqDt8B/OUELvM94HDgQWBLi0qTJKlbzaIIEN+bzEW6bkzERETEDhTdI+dn5kdr7joEuLPZ61SCyE0tLk+SpG5WugWiqudCRETMB+Zn5vLM3BwRXwTOiIh7gJ8CfwYcB7ykk3VKkjTd9VyIAE4HzgRmVH5+O/AocBFF08zdwCsz82udKU+SpP4wY2RkpNM1SJKkHtSNi01JkqQeYIiQJEmlGCIkSVIphghJklSKIUKSJJXSi1M8WyYi5gC3AR/JzGvGOfdhYLe6w+/NzHPaVV+3mODr9BrgfcDvALcDJ2fmpFZE6wURsTtwCXAMsJFiyfX3ZObmMR4z7d9TETELOAc4AVgAfAV4c2Y+NMr5hwAXAgdRLJP/D5l59dRU21klXqvrgZfXHf5GZh7dzjq7SUR8DJiVmSeOcU7fvqeqmnydSr2f+rYlIiIWAJ8HntXEuXtQ/GN/BMVaFNWvj471uOlggq/T0RTbuJ8HHEyx/PjXIqL+g3I6+hywJ3AkxYfA64CzRzu5j95TZwGvBY6neK57U7xW26m8T74K/IDi/XMRcGVEHDMllXbeWTT5WlUcALyLbd8/r2hvid0hImZExPuBk8Y5r6/fU82+ThWl3k992RJR+bD7GLCqyYccAGwGbs3MjW0rrMuUeJ3eAXw6M6+oPP6NwAuBNwAfaEuRXSAinge8AHhKZt5LsQncO4CLI+L9lWXV603791REzAZOAd6amV+vHPsL4N6IOCwzb657yInAIHBKZg4Dd0fEwRQLzE3rxeMm+lpVzn8acFtm9tX2wxHxFOBKiv+Hfj3O6f38nmr6dZrM+6lfWyJeTPEX82FNnn8AcM90/cd+DE2/ThExE3g+sKx6rPI/7bcoNjebzg4HflUJEFXLKJqkDxzlMf3wnjqQ4jVYVj2QmfcB99H4PXE48K3K+6ZqGfD8yvtrOpvoa7UfxR+Bd7W/tK7zPOCXwDMp9k4aSz+/pybyOpV+P/VlS0Rmnlr9PiKaecgBwOaI+BLF5l4PABdk5qfaU2F3mODrtBiYR/Ha1Pot8JzWVtZ19qbx8wbYB7i1wWP64T21d+W20Wuzzyjn/7DBuXOBJcCKllbXXSb6Wh1AMfbm7Ig4FlgPXA+ck5kb2lZlF8jMa4Froal/l/r2PTXB16n0+2nahYiI2JfRU9dQZu5U4rL7A7sA7wXeAxwLXBURO2TmVaUK7bA2vE5zK7f1b7ghoMxr3jXGe62Aa6h73pm5KSJGGP25T7v3VANzgeHM3FR3fLT3xFwav38Y5fzpZKKv1f4U+wclxYDeZwLnUwSO17axzl7Tz++piSj9fpp2IYIiye83yn3Doxwfz1Jgdmaurvx8e0T8LnAqxSj8XtTq12l95XZO3fE5wNoS1+sm471WJ1P3vCNiR4r/KUd77tPxPVVvPTCzEoxqZ6mM9p5YT+P3D6OcP51M9LU6Azg3MwcqP98REVuAf4+IUzNzZZvr7RX9/J6aiNLvp2kXIipJ/u4WX3OIx9Nr1R3AX7by90ylNrxOAxT/U+5Vd/yJbN9E21PGe60i4n6K8SO1nli5bfjcp+N7qoH7K7d71XwPo78n7qfx+2cNxeC46WxCr1Wlj3+g7vAdldt9AENEoZ/fU02bzPtpug8smbSI2CEi7o+It9fddQhwZydq6kaZOQLcTDHFEdg62PIIisGV09lNwFMiorbveimwGvhR/cl99J66neI1qH1P7AvsS+P3xE3AERExo+bYUuA7dQPjpqMJvVYR8ZmI+Hzd4UMogukv2lZl7+nn91TTJvN+mnYtEa0QEfOB+Zm5PDM3R8QXgTMi4h7gp8CfAccBL+lknZ1W+zpVDp0PfDEifgh8k6JpfhHwiQ6VOFW+C9wCXBcRbwH2AD4EnF+dfdGP76nMHIqIy4BzI2IF8DBwGXBjZt5SmVa2BBiovE5XAu8ELo+IC4CjgVcDL+rMM5g6JV6rz1Jpaga+QLGQ0rkUTdJrOvMsOs/3VHNa+X6yJaKx04EHa35+O3A5xUIld1L8Y//KzJzW84ybsM3rlJlfoVjU5DSKxV2eARyTmdN2BDRsbYV5KfAQ8G2KMQ1XAu+vOa1f31NnUIwQvwa4AfgVj6+KdxjFa3IYQGVlxhdR/AP2Q+AtwPGZ+c0prrlTJvJafYbHFzX7CcUCbxdSrBbbz3xPNadl76cZIyMjbatSkiRNX7ZESJKkUgwRkiSpFEOEJEkqxRAhSZJKMURIkqRSDBGSJKkUQ4QkSSrFFSulPhERx1MstrM/xcZhPwYuyszras4ZAY7LzGvaVMMngb0z8+gmz38G8OTM/K9J/M5PAE/LzKMa3PdUiiWnr8nMv6m7748oVu/7/zLzf8r+fmk6syVC6gMRcRLFFr+XAc8G/gD4L+DTEVG71e9eFEvgdosvAM9p18Uz8x6K5dnfGBFbN1Gr7IPyr8CHDRDS6GyJkPrDG4GPZ+Yna479NCICOIXiA5OafVC6xYzxT5mczLwiIv4YuDIingmsAv4d+Dnw3nb/fqmXGSKk/rAFeH5ELMrM2i2QTwfmVX+o7c6odD1sotgy+cTKNS4A/h9wBXAwkMCJmfm/9Y9vdM36oiLiz4F3AQcAIxT7G7wtM78XEcuApwJnRsQJmblvRMwBPkCxidK8yvl/l5m31FzzzZXntUel1mZaXE+k2Pr4YuCuSj0HZubmJh4r9S27M6T+8BHgucBvI+I/I+L0iDgwMx/JzPvGeNzxldvfBz5KsanYfwD/VLneRuDSMgVFxHOAzwCfBPaj2AZ7BvDxyikvA+6j2Ayo2qVxNcX28q+k2Kr4m8ANEfH0yjWPo9hN9gMUmy79BvjL8WqpbNR0EvAXFK0Pb8zMe8s8L6mfGCKkPpCZ1wMvoBgHcQRFqPhhRPwgIvYf46GPAO+ojB34aOXYv2XmlzLzDoodSw8oWdYm4G8z89LMvC8zv0cRIJ5ZqXmAovVjTWY+EhFPowgPJ2TmtzPzZ5l5NnATxc6xUAwcvSYzP56FdwHfa7KeG4FHgc3AzSWfk9RX7M6Q+kRm3gzcHBGzKFoW/hg4GfjviHhaZm5s8LB7Kludk5lriyEU3FNz/3pgTsl6fhQRqyLi3RTbxv8ecCCj/3FzUOX21kodVXNqajiAyviOGrcAz2qipI8DKyn+Xbw6Il6YmcNNPE7qW4YIaZqrzDR4N/D+zFyemVuA24DbIuLbwFcpPmT/t8HDNzU41vQHa0SM+m9MRCwF/ptiBsZ3gH8Bng5cPspDqiHneRThpdZQ5XaE7QdjNgpH9bW8EXgp8EKKfxe/TtG68ZHxHiv1M7szpOlvPcXAwVc3uG8VxQfvwy36XZuAhTU//94Y5/4t8LXMfFVmXpSZNwD7AkRENQiM1Jx/Z+V2j8z8RfULeDvwp5X7fgQcVvd7Dhmr4IjYj6Kr5tzMvDEzv0ExwPKciHj2WI+V+p0tEdI0l5krIuLDwD9FxELgcxTB4pnAOcC/ZuavW/TrvgucFBHfAWZRfDgPjXLuI8BLIuJQ4CHgj4C3Ve6bA2wAVgNPj4gnZuYvIuI64IrKDIyfAa8H/gY4pvK4c4HPRsRtFK0cr6QYC3JTowIiYieK6ZzJttM531W55jURcUhmjvYcpL5mS4TUBzLzDIq1Iv6Q4gP1TuAfKcYPnNTCX/Umiimht1IsWnUFxQyJRt5HMUXzq8D3KWZjVBe+qs7GOB84FvhxRMykaFH5MsWAzp9U7ntZpfWAzPwP4ATgzRQrch5aqWE0H6HoQnlN7ZiQzFwPHAf8H4qZKJIamDEyMjL+WZIkSXVsiZAkSaUYIiRJUimGCEmSVIohQpIklWKIkCRJpRgiJElSKYYISZJUiiFCkiSV8v8D5+zlv4QXwTsAAAAASUVORK5CYII=\n", "text/plain": [ "