{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Chapter 8: Transformations\n", " \n", "This Jupyter notebook is the Python equivalent of the R code in section 8.8 R, pp. 373 - 375, [Introduction to Probability, 1st Edition](https://www.crcpress.com/Introduction-to-Probability/Blitzstein-Hwang/p/book/9781466575578), Blitzstein & Hwang.\n", "\n", "----" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Beta and Gamma distributions" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "from scipy.stats import beta, gamma\n", "\n", "# to learn more about scipy.stats.beta, un-comment ouf the following line\n", "#print(beta.__doc__)\n", "\n", "# to learn more about scipy.stats.gamma, un-comment ouf the following line\n", "#rint(gamma.__doc__)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The Beta and Gamma distributions are implemented in [scipy.stats.beta](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.beta.html) and [scipy.stats.gamma](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.gamma.html), respectively.\n", "\n", "* To evaluate the $Beta(a, b)$ PDF or CDF at $x$, we use beta.pdf(x, a, b) and beta.cdf(x, a, b). To generate n realizations from the $Beta(a, b)$ distribution, we use beta.rvs(a, b, size=n). \n", "* To evaluate the $Gamma(a, \\lambda)$ PDF or CDF at $x$, we usegamma.pdf(x, a, scale=1/lambd) or gamma.cdf(x, a, scale=1/lambd). To generate $n$ realizations from the $Gamma(a, \\lambda)$ distribution, we use gamma.rvs(a, scale=1/lambd, size=n). \n", " * The $\\lambda$ parameter in $Gamma(a, \\lambda)$ corresponds in gamma to using scale = $\\frac{1}{\\lambda}$ and default value of loc = 0.\n", "\n", "For example, we can check that the $Gamma(3, 2)$ distribution has mean $\\frac{3}{2}$ and variance $\\frac{3}{4}$. To do this, we generate a large number of $Gamma(3, 2)$ random variables using gamma.rvs, then compute their mean and var using their corresponding methods in numpy.array (you could of course also use numpy.mean and numpy.var, passing in the array of r.v.):" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "mean of Gamma(3, 2) = 1.5012880719327166\n", "variance of Gamma(3, 2) = 0.7512319256875035\n" ] } ], "source": [ "# seed the random number generator\n", "np.random.seed(317811)\n", "\n", "alpha = 3.0\n", "lambd = 2.0\n", "\n", "y = gamma.rvs(alpha, scale=1/lambd, size=10**5)\n", "\n", "mean = y.mean()\n", "#mean = np.mean(y)\n", "print('mean of Gamma(3, 2) = {}'.format(mean))\n", "\n", "var = y.var()\n", "#var = np.var(y)\n", "print('variance of Gamma(3, 2) = {}'.format(var))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**NOTE**: [lambda](https://docs.python.org/3/tutorial/controlflow.html#lambda-expressions) in Python is a reserved keyword for declaring small, anonymous functions, and so we cannot used the name lambda for variable or functions.\n", "\n", "Try changing the numpy.random.seed input value in the code block above, and hit SHIFT+ENTER to re-execute the code block. Did you get values that are close to 1.5 and 0.75, respectively?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Convolution of Uniforms\n", "\n", "Using [scipy.stats.uniform](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.uniform.html) and [matplotlib.pyplot.hist](https://matplotlib.org/api/_as_gen/matplotlib.pyplot.hist.html), we can quickly verify that for $X, Y \\stackrel{i.i.d.}{\\sim} Unif(0, 1)$, the distribution of $T = X + Y$ is triangular in shape:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEbCAYAAAAS4RmTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAG+xJREFUeJzt3X+8ZXVd7/HXWxAUUAFBowEcsAkDk6RJSe2KYgKSDFbe8GE5eimysLIfCmo3uP4ovI8S45oWCQVkAhIq/spGkOwX4PBDfooMMMEAwRgIggmBn/vH+h7cczjnzD7nrH32HOb1fDz246z1Xd+11mev2XM++/vjrJWqQpKkPjxh3AFIkh4/TCqSpN6YVCRJvTGpSJJ6Y1KRJPXGpCJJ6o1JRZLUG5OKJKk3JhUtCkmuSXLAuOMYlyR7Jbk8ybeT/Oa445GmY1LR2CVZm+QVk8remOSfJ9arap+qunC2x3kceTtwYVU9papOGtyQZGWS+9vru0keGVj/VpKt53rSJM9O8kCSXQbKXp/k9iS7zeO4O7QE+dJJ5WckOTdJ5npsjZdJRRpCki3HHMKzgGum2lBVp1XVdlW1HfCHwGcn1qtq+6p6cK4nraobgc8CbwVI8pPAh4DDq+rW6fZLcnyS42c47j3AycBvD+zzv4G9gV8s7x+1aJlUtCgMtkKSHJPktvZN9/okByY5A9gd+Ez7hv72VvdHklzYvrFfk+SwgWPuN9Cl9IkkZyV576RzHpPkSuCBJFsmOTbJjW2fa5O8ZlL9tyW5sn27PyXJM5N8odX/UpIdZniPU8aa5ALgZcCH2nv74Rku1Y8BX5vTRZ7e+4FfTfJc4FzgzVV1SQ/H/QBwUGsNvRY4Cnh1VX2nh2NrTEwqWlSS7AW8BfiJqnoKcBCwtqp+CbiF7pfSdlX1f5M8EfgM8A/AM4DfAD7Wxie2Aj4J/DWwI/Bx4DWPOSG8DjgU2L6qHgZuBH4KeBrwf4C/GewaAn4O+Gngh4FXA18A3gnsRPf/bcrxkJliraqXA/8EvKW9t2/McIl+DLhiuo1JPtuS1lSvz061T1VdBlwCXAx8pKrOmuH8Q6uq2+iu+58BHwZWVNXtfRxb4zPuJr004VNJHh5Y3wq4bIp6jwBbA3snWV9Va2c45v7AdsAJVfU94IL2i/N1wAV0n/+TWlfLuUmm+vZ90mA3T1V9YmDbWUneAbwA+HQr+39VdSdAkn8C7qqqy9v6J4ED5xDr8TO8x0cleSqwlBmSSlX9zDDHmnTcJ9Bd9+/RtVr69AHgKuAXWvLSImdLRZuKw1v///ZVtT3w61NVqqo1dP37xwN3JTkzyQ9Oc8wfBG5tv6Qn/DuwpG27bVLf/VRjBBuUJXlDkismvt0Dz6VrhUy4c2D5v6ZY324OsQ5rX+DbwM2z2GcYfwJsD9wAvH66SoOtIOBY4NiNtYLovjw8SNetpscBk4oWnar626p6Cd3gdfH9b8+TB3dvB3Zr37Qn7A7cBtwBLJk0y2iq2UyPHjPJs4C/pOt+e3pLflcDfcxUminWYf0YcOVMg9xtfOf+aV5fmKL+r9J1Cx5Od53fNt3MrKr6mYEvBSfQtbomvihM10LaF7i6dS3qccCkokWljYe8vE2T/S7dt/9H2uY7gT0Hql8MPAC8PckT0/2dy6uBM4F/a/u9pQ3Ar6DrxprJtnRJZn2L5U10LZU+zBTrsGYcTwGoqkMGZoZNfh0yWLdNjPhDunGqO4Fz6FoWK2YR07xj1uJiUtFiszXdt+BvAv9BN6j9zrbtj4Dfb90tv1dVDwGHAYe0+h8G3lBVX2/bfhY4EvgW8It0U2ennX5bVdfSdQX9G10C+1HgX/p4UzPFOovD7EtPv6CTPIcuof1SVV3VYnyEbgzkmD7O0fQWszYNcTq41ElyMfDnVfVX445FWqxsqWizleSlSX6gdX+tBJ4H/P2445IWM6cUa3O2F3A23YysG4Gfr6o7xhuStLiNrKWS5NQkdyW5eqBsxySrktzQfu7QypPkpCRr2l8j7zewz8pW/4b2bXKi/MeTXNX2OWm6GSnSdKrq5Kp6ZlVtW1XPq6rPjTsmabEbZffXXwMHTyo7Fji/qpYB57d16AYnl7XXUcBHoEtCwHHAC+lm5hw3cJuLj7S6E/tNPpckaYGNrPurqr6SZOmk4hXAAW35NOBCupkkK4DT2/z6i5Js3259cQCwqqruBkiyCjg4yYXAU6vq31r56XTz6B8zz36ynXbaqZYunRyWJGk6l1566Teraudh6i70mMozJ/qsq+qOJM9o5UvY8C+X17WymcrXTVE+pSRH0bVq2H333Vm9evU834YkbT6S/PuwdTeV2V9TjYfUHMqn1PrOl1fV8p13HirZSpLmYKGTyp0Td3RtP+9q5evY8BYZu9LdtmKm8l2nKJckjdFCJ5XzgIkZXCv5/p1dzwPe0GaB7Q/c27rJvgi8Mt1T4nYAXgl8sW37dpL926yvNwwcS5I0JiMbU0nycbqB9p2SrKObxXUCcHaSI+meffHaVv3zwKuANcB3gDcBVNXdSd4DfLXVe/fEoD3wa3QzzJ5MN0C/0UF6SdJobXa3aVm+fHk5UC9Jw0tyaVUtH6bupjJQL0l6HDCpSJJ6Y1KRJPXGpCJJ6o13KZZ6tPTYud+Tcu0Jh/YYiTQetlQkSb2xpSJtImzl6PHAlookqTcmFUlSb0wqkqTemFQkSb0xqUiSemNSkST1xqQiSeqNSUWS1BuTiiSpNyYVSVJvTCqSpN6YVCRJvTGpSJJ6Y1KRJPXGpCJJ6o1JRZLUG5OKJKk3JhVJUm9MKpKk3phUJEm92XLcAUibmqXHfm7cIUiLli0VSVJvTCqSpN7Y/SU9Dsyny27tCYf2GIk2d7ZUJEm9MalIknpjUpEk9WYsSSXJbye5JsnVST6e5ElJ9khycZIbkpyVZKtWd+u2vqZtXzpwnHe08uuTHDSO9yJJ+r4FTypJlgC/CSyvqucCWwBHAO8HTqyqZcA9wJFtlyOBe6rqh4ATWz2S7N322wc4GPhwki0W8r1IkjY0ru6vLYEnJ9kS2Aa4A3g5cE7bfhpweFte0dZp2w9MklZ+ZlU9WFU3A2uAFyxQ/JKkKSx4Uqmq24A/Bm6hSyb3ApcC36qqh1u1dcCStrwEuLXt+3Cr//TB8in22UCSo5KsTrJ6/fr1/b4hSdKjxtH9tQNdK2MP4AeBbYFDpqhaE7tMs2268scWVp1cVcuravnOO+88+6AlSUMZR/fXK4Cbq2p9Vf03cC7wImD71h0GsCtwe1teB+wG0LY/Dbh7sHyKfSRJYzCOpHILsH+SbdrYyIHAtcCXgZ9vdVYCn27L57V12vYLqqpa+RFtdtgewDLgkgV6D5KkKSz4bVqq6uIk5wCXAQ8DlwMnA58Dzkzy3lZ2StvlFOCMJGvoWihHtONck+RsuoT0MHB0VT2yoG9GkrSBsdz7q6qOA46bVHwTU8zeqqrvAq+d5jjvA97Xe4CSpDnxL+olSb0xqUiSemNSkST1xqQiSeqNSUWS1BuTiiSpNz5OWI8783m0rqT5MalIm7n5JmGfca9Bdn9JknpjUpEk9cakIknqjUlFktQbk4okqTcmFUlSb0wqkqTemFQkSb0xqUiSemNSkST1xqQiSeqNSUWS1BuTiiSpNyYVSVJvTCqSpN6YVCRJvTGpSJJ6Y1KRJPXGxwlrk+Rz5qXFyZaKJKk3JhVJUm/s/pI0L/Ppqlx7wqE9RqJNgS0VSVJvTCqSpN6YVCRJvTGpSJJ6Y1KRJPVmqKSS5Ll9njTJ9knOSfL1JNcl+ckkOyZZleSG9nOHVjdJTkqyJsmVSfYbOM7KVv+GJCv7jFGSNHvDtlT+PMklSX49yfY9nPdPgb+vqucA+wLXAccC51fVMuD8tg5wCLCsvY4CPgKQZEfgOOCFwAuA4yYSkSRpPIZKKlX1EuD1wG7A6iR/m+Sn53LCJE8F/gdwSjv2Q1X1LWAFcFqrdhpweFteAZxenYuA7ZPsAhwErKqqu6vqHmAVcPBcYpIk9WPoMZWqugH4feAY4KXASa376mdnec49gfXAXyW5PMlHk2wLPLOq7mjnugN4Rqu/BLh1YP91rWy68sdIclSS1UlWr1+/fpbhSpKGNeyYyvOSnEjXTfVy4NVV9SNt+cRZnnNLYD/gI1X1fOABvt/VNeXppyirGcofW1h1clUtr6rlO++88yzDlSQNa9iWyoeAy4B9q+roqroMoKpup2u9zMY6YF1VXdzWz6FLMne2bi3az7sG6u82sP+uwO0zlEuSxmTYpPIq4G+r6r8AkjwhyTYAVXXGbE5YVf8B3Jpkr1Z0IHAtcB4wMYNrJfDptnwe8IY2C2x/4N7WPfZF4JVJdmgD9K9sZZKkMRn2hpJfAl4B3N/WtwH+AXjRHM/7G8DHkmwF3AS8iS7BnZ3kSOAW4LWt7ufpktoa4DutLlV1d5L3AF9t9d5dVXfPMR5JUg+GTSpPqqqJhEJV3T/RUpmLqroCWD7FpgOnqFvA0dMc51Tg1LnGIUnq17DdXw9M+qPDHwf+azQhSZIWq2FbKm8FPpFkYiB8F+AXRhOSJGmxGiqpVNVXkzwH2ItuKu/Xq+q/RxqZFj2fMy9tfmbz5MefAJa2fZ6fhKo6fSRRSZIWpaGSSpIzgGcDVwCPtOICTCqSpEcN21JZDuzdZmJJkjSlYZPK1cAPAHeMMBZJm5n5jLutPeHQHiNRX4ZNKjsB1ya5BHhworCqDhtJVJKkRWnYpHL8KIOQJD0+DDul+B+TPAtYVlVfan9Nv8VoQ5MkLTbD3vr+V+juJvwXrWgJ8KlRBSVJWpyGvU3L0cCLgfvg0Qd2PWPGPSRJm51hk8qDVfXQxEqSLZnmgViSpM3XsEnlH5O8E3hyezb9J4DPjC4sSdJiNGxSOZbuufJXAb9K94yT2T7xUZL0ODfs7K/vAX/ZXpIkTWnYe3/dzBRjKFW1Z+8RSZIWrdnc+2vCk+ge9btj/+FIkhazocZUquo/B163VdUHgZePODZJ0iIzbPfXfgOrT6BruTxlJBFJkhatYbu//mRg+WFgLfA/e49GkrSoDTv762WjDkSStPgN2/31OzNtr6oP9BOOJGkxm83sr58Azmvrrwa+Atw6iqAkSYvTbB7StV9VfRsgyfHAJ6rql0cVmCRp8Rk2qewOPDSw/hCwtPdotMmZz+NeJW1+hk0qZwCXJPkk3V/WvwY4fWRRSdJG+Hz7TdOws7/el+QLwE+1ojdV1eWjC0uStBgNe5digG2A+6rqT4F1SfYYUUySpEVq2McJHwccA7yjFT0R+JtRBSVJWpyGbam8BjgMeACgqm7H27RIkiYZNqk8VFVFu/19km1HF5IkabEaNqmcneQvgO2T/ArwJXxglyRpkmFnf/1xezb9fcBewB9U1aqRRiZJWnQ22lJJskWSL1XVqqp6W1X9Xh8JpR338iSfbet7JLk4yQ1JzkqyVSvfuq2vaduXDhzjHa38+iQHzTcmSdL8bDSpVNUjwHeSPK3nc/8WcN3A+vuBE6tqGXAPcGQrPxK4p6p+CDix1SPJ3sARwD7AwcCHk2zRc4ySpFkYdkzlu8BVSU5JctLEa64nTbIrcCjw0bYeuidJntOqnAYc3pZXtHXa9gNb/RXAmVX1YFXdDKwBXjDXmCRJ8zfsbVo+1159+SDwdr4/LfnpwLeq6uG2vg5Y0paX0O6GXFUPJ7m31V8CXDRwzMF9NpDkKOAogN13372/dyFJ2sCMSSXJ7lV1S1WdNlO92UjyM8BdVXVpkgMmiqeoWhvZNtM+GxZWnQycDLB8+fIp60iS5m9j3V+fmlhI8nc9nfPFwGFJ1gJn0nV7fZBuuvJEktsVuL0trwN2azFsCTwNuHuwfIp9JEljsLGkMtga2LOPE1bVO6pq16paSjfQfkFVvR74MvDzrdpK4NNt+by2Ttt+QftDzPOAI9rssD2AZcAlfcQoSZqbjY2p1DTLo3AMcGaS9wKXA6e08lOAM5KsoWuhHAFQVdckORu4FngYOLrNVJMkjcnGksq+Se6ja7E8uS3T1quqnjqfk1fVhcCFbfkmppi9VVXfBV47zf7vA943nxgkSf2ZMalUlX/3IUka2myepyJJ0oxMKpKk3gz7x4+S9Lgxn+fbg8+4n4ktFUlSb2ypPM7N9xuZJM2GLRVJUm9MKpKk3phUJEm9MalIknpjUpEk9cakIknqjUlFktQbk4okqTcmFUlSb0wqkqTemFQkSb0xqUiSemNSkST1xqQiSeqNSUWS1BuTiiSpNyYVSVJvTCqSpN74OOFFwEcCS1osbKlIknpjUpEk9cbuL0mapfl0Sa894dAeI9n02FKRJPXGpCJJ6o1JRZLUG5OKJKk3JhVJUm9MKpKk3ix4UkmyW5IvJ7kuyTVJfquV75hkVZIb2s8dWnmSnJRkTZIrk+w3cKyVrf4NSVYu9HuRJG1oHC2Vh4HfraofAfYHjk6yN3AscH5VLQPOb+sAhwDL2uso4CPQJSHgOOCFwAuA4yYSkSRpPBY8qVTVHVV1WVv+NnAdsARYAZzWqp0GHN6WVwCnV+ciYPskuwAHAauq6u6qugdYBRy8gG9FkjTJWMdUkiwFng9cDDyzqu6ALvEAz2jVlgC3Duy2rpVNVz7VeY5KsjrJ6vXr1/f5FiRJA8aWVJJsB/wd8Naqum+mqlOU1Qzljy2sOrmqllfV8p133nn2wUqShjKWpJLkiXQJ5WNVdW4rvrN1a9F+3tXK1wG7Dey+K3D7DOWSpDEZx+yvAKcA11XVBwY2nQdMzOBaCXx6oPwNbRbY/sC9rXvsi8Ark+zQBuhf2cokSWMyjrsUvxj4JeCqJFe0sncCJwBnJzkSuAV4bdv2eeBVwBrgO8CbAKrq7iTvAb7a6r27qu5emLcgSZrKgieVqvpnph4PAThwivoFHD3NsU4FTu0vOknSfPgX9ZKk3phUJEm9MalIknpjUpEk9cakIknqzTimFEvSZmvpsZ+b875rTzi0x0hGw6SyQObzQZKkxcLuL0lSb0wqkqTemFQkSb0xqUiSemNSkST1xqQiSeqNSUWS1BuTiiSpNyYVSVJvTCqSpN6YVCRJvTGpSJJ6Y1KRJPXGpCJJ6o1JRZLUG5OKJKk3JhVJUm9MKpKk3phUJEm98Rn1krRILD32c3Ped+0Jh/YYyfRMKrMwn39QSdoc2P0lSeqNSUWS1BuTiiSpNyYVSVJvTCqSpN6YVCRJvTGpSJJ6s+iTSpKDk1yfZE2SY8cdjyRtzhZ1UkmyBfBnwCHA3sDrkuw93qgkafO1qJMK8AJgTVXdVFUPAWcCK8YckyRtthb7bVqWALcOrK8DXji5UpKjgKPa6v1Jrp/j+XYCvjnHfUfJuGbHuGbHuGZnk4wr7wfmHtuzhq242JNKpiirxxRUnQycPO+TJauravl8j9M345od45od45qdTTUuWJjYFnv31zpgt4H1XYHbxxSLJG32FntS+SqwLMkeSbYCjgDOG3NMkrTZWtTdX1X1cJK3AF8EtgBOraprRnjKeXehjYhxzY5xzY5xzc6mGhcsQGypeswQhCRJc7LYu78kSZsQk4okqTcmFTZ+q5ckWyc5q22/OMnSgW3vaOXXJzlogeP6nSTXJrkyyflJnjWw7ZEkV7RXr5MXhojrjUnWD5z/lwe2rUxyQ3utXOC4ThyI6RtJvjWwbZTX69QkdyW5eprtSXJSi/vKJPsNbBvl9dpYXK9v8VyZ5F+T7DuwbW2Sq9r1Wr3AcR2Q5N6Bf68/GNg2sts2DRHX2wZiurp9pnZs20Z5vXZL8uUk1yW5JslvTVFn4T5jVbVZv+gG+G8E9gS2Ar4G7D2pzq8Df96WjwDOast7t/pbA3u042yxgHG9DNimLf/aRFxt/f4xXq83Ah+aYt8dgZvazx3a8g4LFdek+r9BN7FjpNerHft/APsBV0+z/VXAF+j+7mp/4OJRX68h43rRxPnoboV08cC2tcBOY7peBwCfne9noO+4JtV9NXDBAl2vXYD92vJTgG9M8X9ywT5jtlSGu9XLCuC0tnwOcGCStPIzq+rBqroZWNOOtyBxVdWXq+o7bfUiur/TGbX53BrnIGBVVd1dVfcAq4CDxxTX64CP93TuGVXVV4C7Z6iyAji9OhcB2yfZhdFer43GVVX/2s4LC/f5GuZ6TWekt22aZVwL+fm6o6oua8vfBq6ju9vIoAX7jJlUpr7Vy+R/kEfrVNXDwL3A04fcd5RxDTqS7pvIhCclWZ3koiSH9xTTbOL6udbMPifJxB+obhLXq3UT7gFcMFA8qus1jOliH+X1mq3Jn68C/iHJpelug7TQfjLJ15J8Ick+rWyTuF5JtqH7xfx3A8ULcr3Sdc0/H7h40qYF+4wt6r9T6ckwt3qZrs5Qt4mZo6GPneQXgeXASweKd6+q25PsCVyQ5KqqunGB4voM8PGqejDJm+laeS8fct9RxjXhCOCcqnpkoGxU12sY4/h8DS3Jy+iSyksGil/crtczgFVJvt6+yS+Ey4BnVdX9SV4FfApYxiZyvei6vv6lqgZbNSO/Xkm2o0tkb62q+yZvnmKXkXzGbKkMd6uXR+sk2RJ4Gl0zeJS3iRnq2EleAbwLOKyqHpwor6rb28+bgAvpvr0sSFxV9Z8Dsfwl8OPD7jvKuAYcwaSuiRFer2FMF/vYb0OU5HnAR4EVVfWfE+UD1+su4JP01+27UVV1X1Xd35Y/DzwxyU5sAtermenzNZLrleSJdAnlY1V17hRVFu4zNoqBo8X0omut3UTXHTIxuLfPpDpHs+FA/dlteR82HKi/if4G6oeJ6/l0A5PLJpXvAGzdlncCbqCnAcsh49plYPk1wEVteUfg5hbfDm15x4WKq9Xbi27QNAtxvQbOsZTpB54PZcNB1EtGfb2GjGt3unHCF00q3xZ4ysDyvwIHL2BcPzDx70f3y/mWdu2G+gyMKq62feIL57YLdb3aez8d+OAMdRbsM9bbxV7ML7qZEd+g+wX9rlb2brpv/wBPAj7R/oNdAuw5sO+72n7XA4cscFxfAu4Ermiv81r5i4Cr2n+qq4AjFziuPwKuaef/MvCcgX3/V7uOa4A3LWRcbf144IRJ+436en0cuAP4b7pvhkcCbwbe3LaH7mFzN7bzL1+g67WxuD4K3DPw+Vrdyvds1+pr7d/5XQsc11sGPl8XMZD0pvoMLFRcrc4b6SbvDO436uv1ErouqysH/q1eNa7PmLdpkST1xjEVSVJvTCqSpN6YVCRJvTGpSJJ6Y1KRJPXGpCJJ6o1JRZLUG5OKtAlozwj563HHIc2XSUUaofbwpJ9uy+9NctIcjvGjSf5lYH2/JBfMtI80Lt6lWBqt44B3t7vTPh84bA7HuAZ4dpItqruz8p8Av9tjjFJvTCrSCFXVV9oD3X4HOKA2vN0+SS6muyHpdsCOSa5om46pqi+2Y3wvyTXAPkmWAbdUeyiTtKkxqUgjlORH6R73+s3qnsq3gap6Yat3APDGqnrjNIe6CHgx3aOte3v6o9Q3x1SkEWmPa/0Y3aNcH0hy0DwOdxHwXuCTVXVbH/FJo2BSkUagPVL2XOB3q+o64D10t92fq68DDwLvn3900uh463tpEUjyIeCrVXXauGORZmJLRdqEJXl2kq8DTzahaDGwpSJJ6o0tFUlSb0wqkqTemFQkSb0xqUiSemNSkST1xqQiSeqNSUWS1Jv/D8LjHQICvipzAAAAAElFTkSuQmCC\n", "text/plain": [ "