{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Chapter 7: Joint Distributions\n", " \n", "This Jupyter notebook is the Python equivalent of the R code in section 7.7 R, pp. 318 - 320, [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": [ "## Multinomial\n", "\n", "The functions for the Multinomial distribution represented by [`scipy.stats.multinomial`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.multinomial.html) are `pmf` (which is the joint PMF of the Multinomial distribution) and `rvs` (which generates realizations of Multinomial random vectors). The joint CDF of the Multinomial is a pain to work with, so like R it is not built in `multinomial`.\n", "\n", "To use `pmf`, we have to input the value at which to evaluate the joint PMF, as well as the parameters of the distribution. For example," ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "from scipy.stats import multinomial\n", "\n", "# to learn more about scipy.stats.multinomial, un-comment out the following line\n", "#print(multinomial.__doc__)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "multinomial.pmf(x, n, p) = 0.041152263374485576\n" ] } ], "source": [ "x = [2, 0, 3]\n", "n = 5\n", "p = [1/3, 1/3, 1/3]\n", "\n", "ans = multinomial.pmf(x, n, p)\n", "print('multinomial.pmf(x, n, p) = {}'.format(ans))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "returns the probability $P(X_1 = 2, \\, X_2 = 0, \\, X_3 = 3)$, where\n", "\n", "\\begin{align}\n", " X = (X_1, \\, X_2, \\, X_3) \\sim Mult_3\\left(5, \\, (\\frac{1}{3}, \\frac{1}{3}, \\frac{1}{3})\\right)\n", "\\end{align}\n", "\n", "Of course, `n` has to equal `numpy.sum(x)`; if we attempted to do `multinomial.pmf(x, 7, p)`, the return value would simply be 0.0." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For `rvs`, the named function parameter `size` is the number of Multinomial random vectors to generate, and the other inputs are the same. When we ran `rvs(n, p, size=10)` with `n` and `p` as above, `multinomial` gave us the following matrix:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "matrix of Multinomial random vectors has shape (10, 3)\n", "\n", "[[2 0 3]\n", " [2 2 1]\n", " [1 1 3]\n", " [3 1 1]\n", " [1 3 1]\n", " [4 0 1]\n", " [0 3 2]\n", " [2 3 0]\n", " [1 2 2]\n", " [3 1 1]]\n" ] } ], "source": [ "# seed the random number generator\n", "np.random.seed(46368)\n", "\n", "rv_vector = multinomial.rvs(n, p, size=10)\n", "\n", "print('matrix of Multinomial random vectors has shape {}\\n'.format(rv_vector.shape))\n", "\n", "print(rv_vector)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Each row of the matrix corresponds to a draw from the $Mult_3\\left(5, \\, (1/3, 1/3, 1/3)\\right)$ distribution. In particular, the sum of each row is 5." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Multivariate Normal" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Functions for the Multivariate Normal distribution are located in the [`scipy.stats.multivariate_normal`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.multivariate_normal.html) module. As expected, `pdf` can be used for calculating the joint PDF, and `rvs` can be used for generating random vectors. For example, suppose that we want to generate 1000 independent Bivariate Normal pairs $(Z, W)$, with correlation $\\rho = 0.7$ and $N(0, 1)$ marginals. To do this, we can execute the following:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "from scipy.stats import multivariate_normal\n", "\n", "# to learn more about scipy.stats.multivariate_normal, un-comment out the following line\n", "#print(multivariate_normal.__doc__)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "matrix r has shape: (1000, 2)\n" ] } ], "source": [ "np.random.seed(75025)\n", "\n", "meanvector = [0, 0]\n", "rho = 0.7\n", "covmatrix = np.array([[1, rho],\n", " [rho, 1]])\n", "\n", "r = multivariate_normal.rvs(mean=meanvector, cov=covmatrix, size=10**3)\n", "print('matrix r has shape: {}'.format(r.shape))\n", "\n", "# take the 1st column of matrix r as Z\n", "Z = r[:, 0]\n", "\n", "# take the 2nd column of matrix r as W\n", "W = r[:, 1]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The covariance matrix here is\n", "\n", "\\begin{align}\n", " \\begin{pmatrix}\n", " 1 & \\rho \\\\\n", " \\rho & 1\n", " \\end{pmatrix}\n", "\\end{align}\n", "\n", "because\n", "\n", "* $Cov(Z, Z) = Var(Z) = 1$ (this is the upper left entry)\n", "* $Cov(W, W) = Var(W) = 1$ (this is the lower right entry)\n", "* $Cov(Z, W) = Corr(Z, W) \\, SD(Z) \\, SD(W) = \\rho$ (this is the other two entries).\n", "\n", "Now `r` is a 1000 $\\times$ 2 matrix, with each row a BVN random vector. To see these as points in the plane, we can use [`matplotlib.pyplot.scatter(Z, W)`](https://matplotlib.org/api/_as_gen/matplotlib.pyplot.scatter.html) to make a scatter plot, from which the strong positive correlation should be clear. To estimate the covariance of `Z` and `W`, we can use [`numpy.cov(Z, W)`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.cov.html), which computes the true covariance matrix." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEWCAYAAABmE+CbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJztnX2YXXV94D/fmdyEO4EyicbWDISw2CY1BDJlViNsnzWsLbQIjiCmCr71he3uurummBqUp7wUl+yTtdju9tku3da+iJhg6FQMbtCHWCs1aOJMgEiiohC4oEaTQcgMyc3Md/8450zOnDm/83Jf5tx75/t5nnmSe8/b955z7+/7+31fRVUxDMMwjK6iBTAMwzBaA1MIhmEYBmAKwTAMw/AxhWAYhmEAphAMwzAMH1MIhmEYBmAKwSgQEblORB4qWo4AEflVETlYtBxpiMj7ReRrCdub9jna5R7NFiKyXERUROYVLUsjMIXQQETkaREZF5GXReSoiOwQkbP9bTeJyFdjjnm1iJwQkfP9H7qKyMbIPs+JyJtn6WM0hbgfjqreo6q/XqRcYVT1n1V1RdFy5MW/r68LXjfzc+Q5d5riMloPUwiN50pVPR14LfAj4H/67/89cLGInBvZ/7eAx1X1Cf/1EeAjIvJzsyJtgxCR7qJlMJpLJ8yCO+EzNBNTCE1CVV8BPge83n/9HPAw8J7Iru8F/jb0+kng68CGLNcRkVeJyAMi8jMR+aaI3BGelYnIShH5kogcEZGDIvLO0La/EZE/91cyL4nIoyJyXo5j/7eIPCgix4B1InKFiAz7sjwrIreGRA1WR6P+CupN0RmkiFzsf4YX/X8vDm37ioj8sYg84sv6kIi82nFPZsxMw7NoEflNEfm2f56KiHzYf//NIvJc6JinReTDIvKYL9NWETkttP0PReQFEXleRH43OlOPXP8r/rP5F//zP+A/u3tCz265v++M1ZR//O/GnDe4r/v8864Pfw4R2SQin4sc86ci8mf+/z8gIk/69+L7IvLvQ/u92V+dfkREfgh8KuYebRKRp/zjvy0ib/ff/2XgL4A3+XKN+u8vEJH/ISKHRORHIvIXIlKOu2cxctwoIj/27/kHQtvPFJG/E5HDIvKMiNwsIl3+tvf735m7ROQIcGvkvVH/c1/sv/+sf433hc6f9L3uLFTV/hr0BzwNvMX/fw/eQP93oe3XAd8NvV4BnACW+K/fD3wNWAOMAov9958D3uy45mf9vx485fMs8DV/20L/9QeAecCvAD8BVvnb/wZvRfIGf/s9wGdzHPsicAnexOI04M3Aav/1BXgrpEF//+WAAvNCsr8/JOti4CiewpwHvMt//Sp/+1eAp4BfAsr+682OezJ13tB7CrzO//8LwK/6/18E/Ir//zcDz0We5zeApb58TwK/72+7HPghsMq/938fvkaMTF8BvgecB5wJfBv4DvAW//P+HfCphHv1FeB34z5f9LrhzwGcA4wBP+e/7vY//1r/9RW+TAL8W3/f8P04Cfx3YIF/36P36Fr//nQB64FjwGsTnsMngc/79/MM4AHgzpTfVSDH7UAJ+E1fzkX+9r8D/tE/33L/vv5OSIaTwH/273M59N4H/PtxB3AI+HP/c/468BJweuj6mb/X7fxnK4TGM+TPhn4G/BqwJbTtH4CfD8183wt8UVUPh0+gqiPAQ8BHki4knpnmGuAWVR1T1W8zfbXxVuBpVf2Uqp5U1W8B24F3hPa5X1W/oaon8RTCmhzH/qOqPqKqk6r6iqp+RVUf918/BtyLN8hk4Qo8Zfn3/vXuBQ4AV4b2+ZSqfkdVx4FtIVnzUgVeLyI/p6pH/c/m4s9U9XlVPYI3eAXXfKcvz35VHQNuy3DdT6nqU6r6IvBF4ClV/bJ/7+8D+mv8PE5U9RngW8Cg/9alwJiq7va37/BlUlX9J7zv3a+GTjGJ9/067t/36Pnv8+/PpKpuBb6LN8GYgYgI8HvABlU9oqovAf8Nz2yaRhW4XVWrqvog8DKwwv8NrAduUtWXVPVp4BNMX4k/r6r/0/9eBZ/hB/53ewLYCpztn/+4qj6EN1F7nf8Z6/letxWmEBrPoKr24s00Pgj8k4j8AoA/cNwHvNf/cVzH9AE8zB8B/yE41sESvFnPs6H3wv8/B3ijvywe9RXVdUD4nD8M/X8MOD3HseFrISJvFJFd/tL9ReD3gVizTgxLgWci7z0D9GWQNS/X4M0ynxGRfxKRNyXs67rmUtz33cWPQv8fj3ld6+dJ4zN4Ky6Ad/uvARCR3xCR3eKZBUfx7kv4mR1Wz/wZi4i8V0RGQt+R83E/8yV4q6m9of3/n/9+Gj/1FWdA8CxeDcxn+ncn+r2JezbRe4+qxj6POr/XbYUphCahqhOqej8wAfyb0Ka/xZtd/hreEvcLjuMPAPcDH024zGG8pe9ZoffODv3/WeCfVLU39He6qv6HDB8hy7HRUrmfwTMHnK2qZ+LZkMWxb5Tn8ZRQmGVAJYOsUY7hDTwARJWqqn5TVd8GvAYYwltt5OUF3Pe9Xo75//aE3kuaGKRxH/BmETkLeDu+QhCRBXirvv8B/Lw/kXmQU88MEp6biJwD/CXexOdV/vFP4H7mP8EbaFeFvlNnqheEUSs/wVs9hL870e9NvSWdk77XHYUphCYhHm/Ds1E/Gdr0z3j+gbvx7PUnEk5zG56dszduo7/cvR/PUdYjIivxzFABXwB+SUTeIyIl/+9f+w6/NGo59gzgiKq+IiJvwJuNBhzGMz/8K8exD/rXe7eIzBOR9Xg+kViFmcI+YJWIrBHPCXxrsEFE5ouX/3CmqlbxTHsTNVxjG/ABEfllEenBW9E1BN+EWAGuF5FuEfltPDu/ix/hvq/B+b4CfArPVBJ8H+fjrWQPAydF5Dfw7OdZWYg32B4Gz0GNt0IIy3WWiMz35ZjEUyB3ichr/GP6ROSyHNechv8b2AZ8XETO8JXUHwCfrvWcMSR9rzsKUwiN5wEReRlvoPk48D5V3R9sVFXFc4Kd4//rRFV/gOesXJiw2wfxnJQ/9Pe9FzjuH/8S3g/8t/Bm4D/klIMwkRqP/Y/A7SLyEt4AOTXz9s1lHwce8c0FayPX+yme3+JG4KfAHwJvVdWfpMkaI/t38ByQX8azaUdj4d8DPC0iP8Nb/l9fwzW+CPwZsAvPWfx1f9PxvOdy8HvARrx7sQr4l4R9bwX+1r+v73Ts8xk8B/aUuch/xv8F7zkdxRvoPp9VQN9n9Qm8z/4jPMfrI6FdHgb2Az8UkeA5fgTvfu327/+X8YIr6uE/462qvo/3rD8D/HWd5wzj/F53GuKNT0anICL/HfgFVX1f6s5Gw/BXTk8ACyK2bsNoG2yF0OaIlytwgW+iegPwO3jRTEaTEZG3+yaoRXirpwdMGRjtjCmE9ucMPD/CMbyl7CfwYrKN5vPv8eznT+H5IbI4640YROSj4iWwRf++WLRscwkzGRmGYRiArRAMwzAMn7Yq9PTqV79aly9fXrQYhmEYbcXevXt/oqqpCYBtpRCWL1/Onj17ihbDMAyjrRCRaBWAWMxkZBiGYQCmEAzDMAwfUwiGYRgGYArBMAzD8DGFYBiGYQCmEAzDMAyfwsJO/bLEX8WrnjkP+Jyq3lKUPIZhdC5DwxW27DzI86PjLO0ts/GyFQz296UfOMcoMg/hOHCpqr4sIiXgayLyxaC1n2EYRiMYGq5w0/2PM1712l5URse56f7HAUwpRCjMZOT3cH3Zf1ny/6ywkmEYDWXLzoNTyiBgvDrBlp0HC5KodSnUh+B3gxoBfgx8SVUfjdnnBhHZIyJ7Dh8+PPMkhmEYCTw/Op7r/blMoQrB7zu8Bq837RtE5PyYfe5W1QFVHViyJEsvbsMwjFMs7S3nen8u0xJRRqo6itfz9fKCRTEMo8UZGq5wyeaHOXfTDi7Z/DBDw5XE/TdetoJyqXvae+VSNxsvq7dzZ+dRmEIQkSUi0uv/v4zX7/VAUfIYhtH6BA7iyug4yikHcZJSGOzv486rV9PXW0aAvt4yd1692hzKMRQZZfRavMbg3XiKaZuqfqFAeQzDaHGSHMRJA/xgf58pgAwUphBU9TGgv6jrG4bRfpiDuLm0hA/BMAwjC+Ygbi6mEAzDaBvMQdxc2qpjmmEYc5vAD9BqZSg6pTSGKQTDMNqKVnMQd1JpDDMZGYZh1EEnlcYwhWAYhlEHnRT5ZCYjwzAS6RT7eLNY2lumEjP4t2Pkk60QDMNwUktm8FyjkyKfbIVgGIaTWjODw2RdYbTrSqRVI59qwRSCYRhO6rWPZ43AaedInXZVZHGYycgwDCf1ZgZnjcBx7fehrSOZKpoWRaeZ1EwhGIbhpF77eNYVRtKKYzYG2bwltQM6KeQUTCEYhpFAvaWjs64w0lYczRxk65nld1LIKZgPwTCMFOrJDN542YppvgGIX2HE7RelWYNsPY7zTgo5BVshGIbRRLKuMML7uWjWIFvPLL+TQk4BRFWLliEzAwMDumfPnqLFMAyjiUQjjsAbZMOKpFGRPUPDFW7cto+JmHFwUU+JnvnzOiJcVkT2qupA6n6mEAzDKJrooLpu5RJ2HTgcO8hmURhZr+kyU5W6BRSqk6fGx1qu0SpkVQjmQzAMo1DichC27604B99GJMu5zgPQLcLC+fMYHa/WfY12w3wIhjGHqTXcspHkDd1sVGSPa/9JVV6MKIOAyuh42+YYZMFWCIYxR2lkdnBg8qmMjtMtwoQqfRnt6XkH+EZF9qSdJ24bUFcGdav7G2yFYBhzlEYlVYXj+IEpB21ldJyN9+1LnVHnzYZuVGRP0nnitgXUmhPRDlnNhSkEETlbRHaJyJMisl9E/mtRshjGXKRRpheXLR48p+ytn98/7b2omWrdyiW5Bvh6k+Wi5+ktl6beO63UNW2bi1pyItohq7lIk9FJ4EZV/ZaInAHsFZEvqeq3C5TJMOYM9ZhewqaPtDjFsHPW5UC+5qI+Z1RRHI1so3n85OTU/4+OVaeZhAIzWJRaciLaIau5MIWgqi8AL/j/f0lEngT6AFMIhtFkhoYrjJ04OeP9LKaXpHDNNFyz5F0HDvPIpktzn69e0iKWsmZaZ6EdsppbwqksIsuBfuDRmG03ADcALFu2bFblMox6aFUHomtA7y2XuPWqVakyJpmI4ljUc8okkzRLbuT9GhqucOvn90+tThb1lLjlypmfLW3WnqfXQZr8jVQuzaJwhSAipwPbgQ+p6s+i21X1buBu8BLTZlk8w6iJVq7v7xrQFy6Yl0m2PCaOUrdwy5Wrpl6fVupivDo5Yz+BzH0T0gbnoeEKG+/bNy2p7OhYlY2f2zfjfFlm7VnMU1medzs00ik0ykhESnjK4B5Vvb9IWQyjkbSSAzHqxHWFU2Yd6F0mjr7eMp9cv2aas3fLOy6cNuCF7fVhJiH1fmWN0tmy8+A0ZRBQndAZ/RUaFbGU9XkP9vfxyKZL+cHmK3hk06UtpQygwBWCiAjwV8CTqvonRclhGM2gVRyIcTNXgVhHcDDQ12P6SJtNx4zTiYTvV1ITnS07D05dP+0eV0bH2bB1hD3PHOGOwdVT565n1t6s5z3bZsciTUaXAO8BHheREf+9j6rqgwXKZBgNoVUciHGDqMIMpRAM6M02fQRJa1kJ368sTXSCY1yroAAF7tl9iIFzFjckYqkZz7sIs2NhJiNV/ZqqiqpeoKpr/D9TBkZH0CplkV2DqEJsHH+zTB+B2SqPMojer6xNdDZetoJSl6SeX6FhJrxmPO8izI6FO5UNoxNpFQeia+ba11uODfOs178QR54w1b7esvN+rVu5hE/vPpR4/POj41PHhKOMkvZvBLU+7ySTUBFmR1MIhlEHST/oRiZP1cq6lUu4Z/ehWPNQlKHhSqp/Ie6YtEEwa5iqS0kF7HjshdRzBHKG7/3QcIUNW0diP1eXCEPDFedzymPDdz1v1znSTEJFmB2tlpFh1Eir16YZGq6wfW9l2kAowDUXxQ9cW3YejB00BZwKJMvnzzKjzWJeOTqWPNsHTwFGGezv47q1y4gzIk2oOp9ZI55v0jnSTEJFmB1NIRhGjbRSaGkcLofyrgOHY/d3mYuUeCem6/PfuG3ftHLaaTPabhFnLaJwyGwWtu+txA7Ydwyu5q71a+iWmWrB9czSnm+W0uG3PbDfeY4sSXGNqNmUBzMZGUaNtEpoqYs88iWZi4I+x1HTh0uBhKudfmjrCD2l5HnnpGrsIHfz0OMzzF1pJDWxGezvY8PWkZij4u9JWlZ1WgTQ0HDFuapJuod5k+IaiSkEw6iRokNL0+zbeeRLMxflyWeIMhaTmZwmz81Dj6c6kF0kKeQ89yRp36TV0YatIyztLXPs+MxaUeFztGIpCzMZGUaNFBlamsW+vfGyFV5v4BClbomVL8lctGHrCDdu2+fMZ6iHqH9iaLjCmtseqlkZQLJCjntmQrzvIen5upTOhOrU80iKcAqU92ybhNKwFYJh1EiRoaVZsnb3PHOE6kRkDh8zpU8yFwWHuPIHgnyG50fH6cqZdAZw3dpl00wsWcJTu0WYVKW3p8TLr5ycVqYiqpDjVlHXXNQ3zRSleL6HgXMWA9Ofp6sst6ssdhZ6y6WWikQLI5rzARbJwMCA7tmzp2gxDKNwlqc4WUtdElvPB2aGdybVN0pDgLvWr4kNo0yj1CVsufZUraN65ADoElgwr4tXqpMs7S2zbuUStu+tzDDJLJjXFTt7X9RT4pXq5Iz942bttZYAd52v2YjIXlUdSNvPVgiGURBZK3dG94GZpSeiuJQBzLSx1+MED7J9wzPdLAlhgYxhB3C9zvhJZaqSamV0PNYhPV6dcA7icQ5gl5M6ujp0rY4W9ZTomT+vZaubRjGFYBgFkDVKJW6fBfO6ckXeRDmzXJqajXeL1HUumD6QB4oha4RQ+NikyKWsDuwwjbJ9uBRVWCnEOdnLpe7YHgytjCkEw6iBPBmscfumdeoCt5+glk5lYUbHq1Mz+Lw2/zi6RDh3045pK5h7H30204AcdgDHRd0AiMB1b1zGrgOH6zIpBSzqKTE6Vs2sMJKytMPyhosGLuopoeo55MOVWFsdUwiGkZOss/u4mWOwr2tQD89G6zGhdOH1GJgNwnkHGz+3DzSbook6gIN7d9sD+6eZb1SZ6rsc9QnkJZi1f8iRjxBHXAQSJFSSFXj5+Mkph34rNUdKw8JODSMnWTJYg5BQmGm6GK9OxGbMwvTZaK35DH29Zc4Mta2cTaoTmui/6BZJDLEc7O+jZ/7MeWrQd/nOq1eToZCpk+CafTnurSuz21lJVpkR3dVKGexJmEIwjAhpJQnSMoCzFHObUE3NYYiLg08jOMdohro/s43gfe60vsRJFVcH+/sSm+z0lt2KsK+3PBUNNXbCnTQWd9048irsVslgT8IUgmGEyJLw5RoIgvez/PAXzu9OTUoKEpfyEMxEz0wYGIsibDbbsHWEm4cen9o2NFyh//aHEk05S3vLUzkTcfT1lhm55df55Po1TmUbPN8shfLC140jr8Ke7eZItWA+BMMIkcXZm1ZyIEvHrmMnvGOTyj0DNSVBPT86Tm8DTUaXnLeYbx16sW5ndhgFPr37EDsee4GjY9XUKKIgwzpLRdZoSOiZ5RIinoM3b/JcUuZ5cJ0bt+3LdM7K6DiXbH64pR3MtkIwjBBZCsKllRzIOnPMalOuZSbaKJNRtwjXDizjmov6nH6Peghm6mnD6UTIQRtHcHxg6gsie+5av4bjJyc56kcV5Y2qSksiG+zv4xPvvDDz82m1EulRLFPZMEK4smVdzVtc4adZ4vAF+MHmKzKdE7LNRINM2HpKK0QpdQtocrLbbNBbLvHSKydj70GgqqJ5AKeVunKZh8KkNewJE44qa/S5G0HWTGVbIRhGiDwF65L8DbsOHE6d9cbZlF3nBK9MdBLdIlPNb7L2Fc5CUuRQ49cMbkbHq4k1leKiuWpVBnmLFAY9prPej1Z1MJsPwTBC5ClYl+RvSPvBu7qQJZ0zzTcxoTqtSNtsjNazvWaoJWM57/nrKTGRxX8U7NeKFKoQROSvgbcCP1bV84uUxTACokohsPVHB4gkf0PawODqQpZ0zrvWr0ktqBaOd59R6bQDCGcDZyVrkl4jzDiubOswRfc8SKJok9HfAJcXLINhTCPObBMNk4Tk8NM0R3BcYtTQcIWuhIS1qDPbRWV0vGH+g1ZEIVdyWhZl0KhBOi7g4Pq1y1qq50ESha4QVPWrIrK8SBkMI4qrJME9uw8xcM7iqR/zupVLYhu5rFu5hMH+PvY8cyR2e6lrZpOaJCd0eLCKFlSrh6ASZ7spj95yKVM11ayEfS+NoNV6HOSh5X0IInIDcAPAsmXLCpbG6FTCkT1JjWLC+QiukgbB+67tUcfo0HAlMSIpWvagljr8cVxxwWu5Y3A1Q8MVNmwdmXV/QK00Ovo17Htp14G8URRtMkpFVe9W1QFVHViyJL7IlGHUQ9RElESW4nPB+67tk8q0WHRXslWYINrotgf2NyxBbMdjLwDejHY2lEGDgp4YHasmlqiohWbUGkorgdKKtLxCMIxmk6X2UECW4nPB+0mRJOEBKGsIYj1hlHEcHavSf/tDqd3XGoEAb/pXixtyrqW9ZW69alXDB69Gms6ylEBpRUwhGHOerANyluJz4X3SHMuV0XHO3bTD6UieDRqpYJJQ4JGnjtR9HsG7b1t2HuTda5c5VwrlUjfXr12WK8O7kZnYaRVxW5VCFYKI3At8HVghIs+JyO8UKY8xN8kSE95bLjmLz8VFkAQ+ibSVRy3lFOYSpS5hkV+XKdpXYvveCrdetYqnN1/BJ9evmfEc7hhczTUX9WVOx5hQbZhpJ0sJlFak6CijdxV5fcOAbLHjCxfMi3U4xkWU1NqAvVuESb889LqVS9j6jWcLLxdRBN1+Abq+UIJYXEkRV7/jMF/Y90Iu/0ijmtm48lBaNSEtoOWjjAyj2WRpDu+a2Q0NV6YdF8xma3H8Tqpy1/o1bNl5MDZcda7wiXdemCsJ0NXBbs8zR2oKT82iaNJIq4jbqpgPwTB8XnrF3TRFYYY5YWi4wsb79k0bdI6OVWu2yyteieZ2ywtoNHF29iQHvstef++jz9YsQ73PIK0ibqtiKwRjzhPMMNNs+VFzwpadBxtu0mlXA1EjawxV/Fl/ePBct3LJjFyNYMa9wdFUpx7fjMAMGZJwVb1tdQUQxRSCMefJE3YaNifkcRA2uyjbbNPdJUyElGGjP1tY8Q4NV9i+tzLtGgJT2cVJpr40FvWUGPV7JYSJJiGGiQ7+61YuYfveygyTVSB/O2EmI2POkzfyI0u4aG+5NM1ccF3OEMhWZlFPiTMWNHcuGQ7RdJUS+cK+F+i//aGalUG51M0tV65yKrO470VcfsE9uw+1ZYhpHLZCMOY8WUsWh0kKFy11C7detSp2dpjWNKfVuX7tMu4YXM25s5DMFjiMXc8mjyIQoLenhCq8OF6d1njItXqL81u4lJNL/nbDFIIx53FFhFxzUd80U0AWFs7v5uNvn+k8jDN7tBuLekrcMbgaqE2JdufsZ9zbU5oyvdRDtwhP3fmbsdsu2fxwao/mMHkG+VYPMY3DFIIx5wkG79se2D8VIbRgXhcD5yxm4JzFqUXvwvT2zJ86X9jWnLe5eysS7tOcJXcjyoQq5VJ3pmNK3RJr268F131PWn24+lW4FGF0ldEOIaZxmEIw2g5XREc95+ntKfFyKOx0dLzKTfc/zp1Xr55qmuLqtxwmmEFGY+PbXRmANxhG79mCeV2ZTTdBmekdj70QG5orAqreSuTlV042bDXl6j2RtPqIOwaSV5O7Dhyu+ztZNKYQjLbClYQE+SI6hoYrbPzcvqmuYnED1Hh1gtse2D81AJ5ZLlHqlsROZIGZIE/kUjsgwPJXlafd+6NjVcql7sz9CYIy03de7ZmdXEr9ks0PN6zGkmumnvR8kmb3eVqstiOibTRzGRgY0D179hQthlEDeWf1rv1ds/S87Q/7b3+opkGn1CWcfto8jo5VZ5gJgtd9NdjX5xJpz+rcTTtyrQ7CfZAh22CddI1Prl/TMQN8gIjsVdWBtP1shWA0nbyz+qT9G1U0rNYZaHVSOTpWpc+PP9914PDU4B8uvGa4Sbs/eRzWcQ7jLIO56xp9fqvSuYrlIRhNJ28p4KT903oQzBaV0XE+vfsQz79og38tLE9oGrPxshW5KpTWQlxpcsHLiJ7LmEIwmk7eWX3S+2k9CLJSLjXmq99GFteWozI6zoatI9w8NN25O9jfx3Vrl2VSCi7nbxqD/X0zSmMrsH1vpeWb2DQTMxkZTSdvKeCk/fM69eJ8EQAn52BZ6VZE8ZL1ov2M7xhcPRXym2Q+qmdGv+vA4Rl+hEZUOm1nbIVgNJ28s/q0/Qf7+9h42QqW9pZ53u+eFTeriyszsGHrCB/aOpIYKdQoekpdHVOuopkEdYOiDPb38cimSxNXAbsOHAZq61/crk1smokpBKPp5C0FnLZ/1n61ecoMNIP/dvUFU5/DSCZpEE4yB4b7IeTtX9wq/qhWwsJOjbYja+jpbDSPTyIcvpglqW0ukxaKuua2+CJ2gbKtJRQ5rrNdudQ9Y/LRCTkHWcNObYVgtB1pS/2h4Qr9tz80myLFcuO2fVMmjOWvip91XnLe4lmWqjj6estcct7iGc7iLEEBt161ymlGdH0fgqq0LhNSo1ainYStEIy2I2mFUEuNnTBZa+3kpdP6IdSKABeft5infzqee9adN1kxTHTmn4VGJUG2ApaYZnQsGy9bMa3sBHjF0DZetqKukhG95RK3XrUqNbKlFkwZeCjwL08d4a6EbOCk7mNxx2SZBNQSPTQXnc6FKgQRuRz4U6Ab+L+qurlIeYw2Iq7FFek/VlcJ5lLX9B4G9awyjGTSupG5stQhPtw4GorcqP4EecOlO4HCTEYi0g18B/g14Dngm8C7VPXbrmPMZDQ3SHPkJS3lwV0aIWw2SLqGOYBnh6AGUVACJKlM+KKeEq9UJxMdwAGu59ctwifeeWHmlUm0NWbSNVuddjAZvQH4nqp+H0BEPgu8DXAqBKOzGRquTOtJAPF1j5KW8netX+Oc3Z8Wyk5OaoDeySZJ0HXaAAAcJ0lEQVSB2SIoCZ3UIS5w1H5696Gp91ylKFzVaONWGi4T0oSqs4ZW3Mpk+95Kx5S1zkqRCqEPeDb0+jngjdGdROQG4AaAZcuWzY5kRlMIz8DOLJcQ8ZquuGZjAdEfvmspH5gigh9xZXR8mjP36Fg1cUDopGY2RROeRYcH/EYTRBLFmZBu3LZvxnN0KRFX/axdBw63nQO5HooMO40rVTLjV6iqd6vqgKoOLFkytwtPtTPREL7R8SpH/Y5YrkblYcKz9rhM5oBgxnns+El6yyVnaYIk2YpSBqVu4fq17T/pCVcMDVpu1kPQc8FFXEjoYH8fk47nGLcCzONAriUrul0oUiE8B5wden0W8HxBshhNJi36J20IDjvywvHjLkbHq86mLRU/uzWrbLNFdUKbOpueDYJorzC1ZGp3i0zLDYjLQ4gSVfZ5MpGz7tvpuQlFmoy+CfyiiJwLVIDfAt5doDxGE6nHLl/qEsZOnIw1DeRtphIQNh2Zz6CBKOx55sg052zP/PzzTpfz99bP70/szhZdScZlIsclwWXdN6k0eyf4FgpbIajqSeCDwE7gSWCbqu4vSh6judQaqlfqAoRp5qXwjKzW84Znk/WGEfaWS1bEzqc66a1ywjPo7/74WK5zLOopxQ6ug/19LFyQPId1rSTTamhl3bfTcxMS766IfAh4BBj2B/CGoqoPAg82+rxG61FrBvGkChOTbsdgPZnJwY+4nnMIXlkFYMasOO9AaHiz8luuXOXcnjTwxs3ok6LJomTZt9NzE9JMRmfhJY6tFJHHgH/BUxBfV9UjzRbO6ByiyUPhKKPenpKzpaXLwRsMDMF5o+Gq4JmaEJylroMfcXCOP9g2Qt42CRo6PporYeSjWyQ1xt81IGc5thHkMUO1I4kKQVU/DCAi84EB4GLgt4G/FJFRVX1980U0OgXXDOySzQ87FYIrszhqGnAlm0G8sgh+xMExtSaiuaJfLLEtP5OqqQO6a0CerWSxvA2a2o2sTuUy8HPAmf7f88DjiUcYRkaSzADveuPZsdmicTMyl8IJK4vK6DjdIoxXJ7j18/s5duJkXc1yTpyciM1wtWJ2pwgK2u3+/tHEkN4sZpdWGJDzmKHajcTSFSJyN7AKeAl4FNgN7FbVo7Mj3nSsdEVrUm/NeFepgXKpi8ULF0wN4hOqUxVN85w/LgO6kZS6ZZpSaQVl0CXkNn818zzh55alD0FAp/QjKJpGla5YBiwAvosXGvocMFq/eEankFSMLOsPN84MUOoSTk7qlKKYUJ1aGSQ1LwGmrQRmI8ksusIoWhlA8iCeVWF9cv2a1DDPJMLXiftepA30jfhuGflILW4nIoK3SrjY/zsfOILnWL6l6RKGsBVC69GomvHRwX3sxMnYGX1w3rhZZpoT2fAG6Z753Rw7kRxRdf3aZdwxuDqx61yf7+CNKphyqZvTSl2Jzy8LndSPoGgaVtxOPY3xhIiMAi/6f2/FK043qwrBaD0aFZcdtcue6xiIKqPjnHfTg7Ez/2ojbBsdjkKiMoiadlyrifCgHLdS27B1JPb8eb4XnR7z34okJqaJyH8Rkc+KyLPAV/EUwUHgamDu9P4znDSyUXm4RozEVbrymQ0zUFA6YS4RDPJhk47rTh87fnKqlg/AI5su5a71awDYsHWELscDVMhc/6eR3y0jG2krhOXA54ANqvpC88Ux2o11K5fMKHFcS1x21ARUdLHRCVV6y6Wa7eftRtwzS5qJB/clsOvveebItGiwJKWd1RfQ6TH/rYj1VDZqYmi4EutwFOA63/6chyKa0ixKSIibC/T1lnl+dJzenhKq8OL4qVLkQfnwrLgc+GmO/bSoMYsyagzt0CDHaFPiHLoBCuw6cDj3OYtI5Do6VqXUBdXJ5l2jW4RJ1ZaIPArTWy5NOec33rdvyv8SbViTFdeg7ypBHZC2WujkmP9WpMjy10abklYuOlpeOq1+fOC8nG2E5ioD8AbKIpXBgnldXvRViKB/NHjVQ7M64/t6y87M7G6Hz2Bpb9m5LSCuR4VRDKYQjNxkifIIKpJmqR+f5LxsJq02a28Gx09OsuXaC6dV8dxy7anS0ll9JILnOI7rS1AudfOuN54d+/66lUsyBQFY5FBrYCYjIzeuAmNhwrM+V/14IFMdoVbI/C2S3nKJVUvP4JGnaq8nWW/cfpcIQ8OVxKSygXMWxzapz4JFDrUG5lQ2cpPkQwgTGApc37BSlySaKwL7+2mlLo6fnGRSvVIKKIQtPZ2uMJ7efEXdTvdFPSVuuXLVDHt8/+0PZXasR8tLpDl8s8o8m8Xp5ipZncpmMjJyE20mkmQ/Tpr5pdmuA/v7eHVyqhTDpIJ0ncoR6BbhurXLeHrzFYl9d9uVoP1kvSaVo2PV2FaPt1y5ilL39OfX3RWfgxFe2WUxBSbJvKinlNq0xph9zGRk1EQ4+sNVrCyIF9/4uX0NLScRbpgzocr2vRUGzlnMrVetmhYx0+pEi+JFCd/DLGa6NKKtHoMZfnVCZxQPTMs0ztJK0iVzvWVNLPS0edgKwaibpPaDg/19LJzf3HlHeOY6f177fKVPXzCPhfPdrTeDz3Xz0OOMnWhMw8JgQA/P8GFm8cDenvjV1tKUFUu0p3GcozlPYlmnN7VvNWyFYMyglhlZUrz4i7OQ7VsZHW/4SiQPQZnoqD9DgNe9ZiHf+/GxGX4OLw9CElcKteYFuAgGdNcM/7YHvLbmL78yUwGVuiV1xRJtXBRcq9bZfac3tW81TCEYwCklEK1e2YiSw40wd6TRVUeV00V+pm6tZSoC53dcVq4C3z885nR6Vyd1Wt2mZjrIw7Nz1wz/6FjVmZuwcP68qe9A1rIS9SaWWYG72aWQ9bWIXCsi+0VkUkRSPd9Gc4maD6JDQb2JQ7NRe6Yet8HRsSovjlfpLZcSTTguAue3K94+LQ4/vLlRyqDUJVy/dtkMMx540T9J13EpxvBKL2om7C2XOK3UxYatI5mL12XBCtzNLkWtEJ7Aq5j6fwq6vhEiLfMY6puRDfb3NbVjWSNQal8htBoC05LPArKGC7uIDsLhXtbNamRjBe5ml0JWCKr6pKparnqLkGWwr3dGdsuVMzNc51p56dmgXOrmrvVrYgfiLIofPBNaHmdwkp2/XpICFozGYz4EI9XG34gZWdjBGNdlK0ynJ5rVQtCdLK3HcdJgmUXxl0vd3HKlV+coqzPYdd7K6DiXbH647nBRK3A3ezRNIYjIl4FfiNn0MVX9xxznuQG4AWDZsmUNks4IE9fTIBiU+/wSBFt2HmTD1pGG/LCTMliD64Vr6891wnH7adm/waw87vmkKf5oKeqsz9h1XuFUFVvrh9weNM1kpKpvUdXzY/4yKwP/PHer6oCqDixZsqRZ4s5ZhoYrbN9bmaEMguzfjZetYPveSu448KQKp64ZZVBA7Y7B1TMcluHM1uvXLkutoNkphEM9IX2Wn/R84vICAsI5CHmJO2/cKs+qmrY+ZjKa48TZf8M9DWqJA09zMmaNYU8anAbOWcyHHNm0WekWYcE8YazZNbDrIBzqCfkKC0bvX/D6xm37ZkQ+1RPbH5dv4JLRwkVbm6LCTt8uIs8BbwJ2iMjOIuQw0uO8a4kDT3MyumaUldFxzrvpQZY7+iaEGezvY5EjmzYL5VI3n3jnhYw3SRlE6wPVSjSpL2mWH8b1fAb7+5xNa2oZrIOVYFDq4q71a3hk06VTNZiiWLhoa1NUlNE/qOpZqrpAVX9eVS8rQg4jPc67ljhw1+wweD8cOQLTzQvBzDWLaeqWK1fNaP4Sh0BsTH6wWmkGC+fPcw6KYdJMYHGhnndevTrVZJb0uRoV259UVqIRZSuM2ad9Cr8YmUjrThYlabZ+yeaHWbdySe4ftmuwCr8/2N83NZN0Bc2k2ZwH+/vYcu2FiVVOwz2eH9l0KT/YfAWPbLp0Kn7+2PHG1AiK8uJ4lY2XrUhcKQRO3O17K7HJa677PNjfxyfeeWGqP8BFowbrNHOihYu2H+ZD6CBqSRBKCgetjI6zfW+Fay7qY9eBw87wwWjtozwZu2lmirTtUV9DuARHUEpi14HD05q7BPvVk6SVxtLe8tT1Pnr/YzP8FEE3sTh7PnjKM2kAjT63aLXSpIG3ETWGIN2caOGi7YcphA6i1kJgSeGg49UJdh047CxXHNek3UWczT/NSRo1Y6QV3gv+H1WMG+/bx20P7Gd0rMrS3jLHjp9smjIQTpXrCGfzxnUTS2pOn6WgIJwa2LMog/CxcfvlKWyYJTjAaC9MIXQQ9RYCq+X4PE3a48a+oPZ+3BkEL0ciSG46rdQ1zQnsWgHFKcbqpE6Vzmh2ob2e+d1s2DrClp0HpwbU6AB8yeaHExVSlkG10SUj8p7Pykp0HuZD6CDqdRbWcnye+j9xZbAH+/u4bu2y2DIW8+d1sfUbz045LeMiguL8DEWFNgpedNGxExOpORuNyAxvdMmIvOczP0HnYSuENiDrMr7eGVtcxnIjZ3xximVouMKuA4dRZiYzHT+ZLSQ0qgAaWW67XOrizqsvyNRrQZlZgjvOZJfm6M86qDa6NHQt5zM/QWdhCqHFybOMr8dZ6MpYvuai5B/8op5SpiqmcYol+tlqrV8UVTRxiq1WxquTU59/w7aRWLNXGpXRcc7dtGOa7yCJektG1GrDz3s+a23ZeZjJqMWpZRkfDa+s9TrhjGUXcU3ao7giZrJW30wi7MAFuHno8YYpg4Bg4EtSBkK80zwgMCHds/tQ4mfOkrsQ0OhY/zzns9aWnYkphBZntjpG1Xqdwf4+trzjwmlJZmGCjOA4xdSIz3DxeYunNY2vVRm48ttKXUxrHpQkR1yJ7yhJsuUdzBttw89zvmaWvDaKw0xGLc5shfbVc52wHTnOjADElkHOYutP6jcM8PRPTx2/ZefBmpRBudTNNRf18ZlHD80oLV2dhOpk+irmW4de5NoBz/4ffP48sqTlHbhotA0/6/mstWVnYiuEFme2SgA06jpRkxWQq7xBqUumVTbd8o4LE69XGR1n+aYdLN+0IzUHwlU99c6rV3PH4OpMNYJcjFcnuHHbvmk1fVzmnzyrqFbFWlt2JrZCaHEalVVa1HWSTAuBwki7ZpCNWw9XXPBa7hhcnbjPsRP1+TOidZiuuahvRl+HYDWSlPndDlgOQmciWkvYREEMDAzonj17ihbDyMG5m3Y4k85+sPmKTOeIZkPXQrnUnWqSWb5pR6Zz5enotqinhKqXg9Gug78LizJqH0Rkr6oOpO1nKwSjqWT1TSQNLsG/t35+f65EuDBZSnj0lkup58/b0e3oWDWxz3E7YzkInYetEIymEldErtQlnH7avKm6QnEDbLhKaRjXiiMLaauStJVI+Pih4UouBRVug2kYs42tEIyWIOqbOLNc4tiJk9PqCn1696EZxylwz+5DDJyzeNosNC0yKVA2cclyXSLcPPT4lP2+N8acs+XaC50VSKOrmqyZ1GDRN0Z7YFFGRtMJRx4tXDAvtQREgMKMuPaNl61wNsXpLZfYcu2FznyACVU+vfvQVMTT0bEqo+PVadFPQGyvgajDNG9SnUXfGO2ArRCMukiy/cdtyztTjt0/og9K3cKWd8wM23TN9F2MVye47YH99Myfx3h1IrHHQNLnKJe6LfrGaEvmxAohbxcxIxtJ5Qtc23pz9kGOzqy37Dw4Y4VRndAZK4mk3sFJHB2rTpmkJlSnBvOosnHN+IO8BqsAarQjHb9CaHTN+GbTTqF8aeUL4rYtmNc1YwadxLHjJxkarrDnmSPc++izzhl/eMY+VXsoz4dx4IpOSorDryX6pp2eu9G5dLxCqLWLWBEUpbxqHYxqKV8QROUE5phFPSVeHK/OKBkR3v8Pto6Q5r4NZuzNaI0Z93kamcjXbpMWo3MpRCGIyBbgSuAE8BTwAVUdbca12qnmShHKK20wSlIWaTkGSdFAgTnmigtey9ZvPJto3klTBmEbfZKzN8ghiItqSsJlHmpUHH47TVqMzqYoH8KXgPNV9QLgO8BNzbpQO9VcKUJ5JQ1GcX6ADVtHuHnIUxhJ9Y/itkUZr05w76PP1pWBHLXRu+6VAI9supQ7Blc7awwt6ik1rG5UHr9VO01ajM6mEIWgqg+p6kn/5W7grGZda7aKwzWCIpRX0mDk6pFwz+5DDA1XEsslR7e5yBMFFEfUTJPlHsaFrpa6hFuuXBX7eYBcQQl5ewW006TF6GxaIcrot4EvNuvk7dT3tQjllTQYuZRFOD8gqSFPeJtrVt4tyc110ojLU8h0D6OX9V/nqdaaJFOeXgHtNGkxOpumKQQR+bKIPBHz97bQPh8DTgL3JJznBhHZIyJ7Dh9O7t7lotYuYrNNEcoraTBKmqHmNWe4rvOuN55dV9npQI7ARLNh6winlbroLZec9zBr6Gqwb95GMHlNQO00aTE6m6Y5lVX1LUnbReR9wFuBf6cJBZVU9W7gbvBqGTVUyBZktguGpUXLbNg6Ehu+mdeckXSdgXMW11ziemlveYZjPK2gXNYBe2i44pQpSSHW0mzICsUZrUBRUUaXAx8B/q2qjhUhg3EK12A02N/HnmeOzGhLWas5I+k6g/19uQvXBXLkjdLJMmAHvZldJA3u1ivAaFeK8iH8L+AM4EsiMiIif1GQHEYKdwyunur+1WxzhmuQ7RaZamIfZwrKa6JJs9mn9WZOG9zNBGS0K4WsEFT1dUVc16iNeswZeZLeXDPrtME0r4kmzUyWluWcZXA3E5DRjnR8prJRHLUkvYWb1GfN/nUlm61buWSaLNHzuvoTpDnMwxFWhtFJmEIwmkZahE6csrjz6tW5G8nsOhAffRa8n7c0RFrPBSstYXQqrZCHYHQoeZPexqsT3LhtX+6qtGk+hEbkBURJCz01jHbEFILRNGpJeptQzZwAluU6kK4womUmgExZ1lZawug0TCEYTaPWpLcA1yw8OoCvW7kkMWooSWG4ykwAU8mMveX4Hg5WWsLoNEwhGHXjKuSWFH6ZxSwD8cli0QF8+94K11zU5wzzTFJMaeakoeEKx06cJEqpSyyvwOg4zKls1EWawzYpGQ1OhX52+f0RosR1TIsbwHcdOOx0RieFmW7YOhJ7TNj/ENcD+vTT5plD2eg4TCEYdVFPLf+wsohrbBOXAJbFHxA38LsUU1oOg+t6o2PVxM9mGO2IKYQamcstD8Of3ZXAldfhmrUDWdIAXkvnsbQyE7XUJTKMdsUUQg3M5ZaHWVtU1jJgZsnuTRrAa1mtpCkiq0tkzCVMIdTAXG55mNSiMqCZA2Y9/oCkc9aqMAyjkzCFUANzueVh0mcUmJUBs1Z/QKOvZxidhimEGpjLdmXXZ+/rLecuOdFozLxjGPVheQg1MJdbHoYLxmV5fzaxstOGUR+2QqiBuWxXTiskVzRm3jGM2jGFUCNzdeCZy/4Tw+h0zGRk5CKtkJxhGO2LKQQjF3PZf2IYnY6ZjIxczGX/iWF0OqYQjNzMVf+JYXQ6phBajLlcI8kwjGIpRCGIyB8DbwMmgR8D71fV54uQpZWYyzWSDMMonqKcyltU9QJVXQN8AfijguRoKfL2/jUMw2gkhSgEVf1Z6OVCcFZRnlNYjL9hGEVSmA9BRD4OvBd4EViXsN8NwA0Ay5Ytmx3hCmIu10gyDKN4mrZCEJEvi8gTMX9vA1DVj6nq2cA9wAdd51HVu1V1QFUHliwpvl5OM7EYf8MwiqRpKwRVfUvGXT8D7ABuaZYs7YLF+BuGUSRFRRn9oqp+1395FXCgCDlaEYvxNwyjKIryIWwWkRV4YafPAL9fkByGYRiGTyEKQVWvKeK6hmEYhhsrbmcYhmEAphAMwzAMH1MIhmEYBmAKwTAMw/AxhWAYhmEAphAMwzAMH1MIhmEYBmAKwTAMw/AxhWAYhmEA1kKzJbE2moZhFIEphBbD2mgahlEUZjJqMayNpmEYRWEKocWwNpqGYRSFKYQWw9Uu09poGobRbEwhtBjWRtMwjKIwp3KLYW00DcMoClMILYi10TQMowjMZGQYhmEAphAMwzAMH1MIhmEYBmAKwTAMw/AxhWAYhmEAIKpatAyZEZHDwDORt18N/KQAcRqByT77tKvcYLIXQbvKDdNlP0dVl6Qd0FYKIQ4R2aOqA0XLUQsm++zTrnKDyV4E7So31Ca7mYwMwzAMwBSCYRiG4dMJCuHuogWoA5N99mlXucFkL4J2lRtqkL3tfQiGYRhGY+iEFYJhGIbRAEwhGIZhGECHKQQR+bCIqIi8umhZsiIifywij4nIiIg8JCJLi5YpCyKyRUQO+LL/g4j0Fi1TVkTkWhHZLyKTItLyIYUicrmIHBSR74nIpqLlyYOI/LWI/FhEnihaljyIyNkisktEnvS/K/+1aJmyIiKnicg3RGSfL/ttWY/tGIUgImcDvwYcKlqWnGxR1QtUdQ3wBeCPihYoI18CzlfVC4DvADcVLE8engCuBr5atCBpiEg38OfAbwCvB94lIq8vVqpc/A1wedFC1MBJ4EZV/WVgLfCf2ui+HwcuVdULgTXA5SKyNsuBHaMQgLuAPwTaykuuqj8LvVxIm8ivqg+p6kn/5W7grCLlyYOqPqmqB4uWIyNvAL6nqt9X1RPAZ4G3FSxTZlT1q8CRouXIi6q+oKrf8v//EvAk0BZNStTjZf9lyf/LNK50hEIQkauAiqruK1qWWhCRj4vIs8B1tM8KIcxvA18sWogOpQ94NvT6OdpkYOoURGQ50A88Wqwk2RGRbhEZAX4MfElVM8neNh3TROTLwC/EbPoY8FHg12dXouwkya6q/6iqHwM+JiI3AR8EbplVAR2kye3v8zG85fU9sylbGllkbxMk5r22WEV2AiJyOrAd+FBkNd/SqOoEsMb37f2DiJyvqql+nLZRCKr6lrj3RWQ1cC6wT0TAM118S0TeoKo/nEURnbhkj+EzwA5aRCGkyS0i7wPeCvw7bbGElhz3vNV5Djg79Pos4PmCZJlTiEgJTxnco6r3Fy1PLajqqIh8Bc+Pk6oQ2t5kpKqPq+prVHW5qi7H+wH9SqsogzRE5BdDL68CDhQlSx5E5HLgI8BVqjpWtDwdzDeBXxSRc0VkPvBbwOcLlqnjEW92+VfAk6r6J0XLkwcRWRJE/YlIGXgLGceVtlcIHcBmEXlCRB7DM3u1S3jb/wLOAL7kh8z+RdECZUVE3i4izwFvAnaIyM6iZXLhO+4/COzEc2xuU9X9xUqVHRG5F/g6sEJEnhOR3ylapoxcArwHuNT/fo+IyG8WLVRGXgvs8seUb+L5EL6Q5UArXWEYhmEAtkIwDMMwfEwhGIZhGIApBMMwDMPHFIJhGIYBmEIwDMMwfEwhGEaN+OGrI5G/SRH5jaJlM4xasLBTw2gQInIDXj2qdao6WbQ8hpEXUwiG0QBE5JeAh4GLVbXdSrAbBmAmI8OoG7/mzWeAD5syMNoZWyEYRp2IyGbgtar6vqJlMYx6aJtqp4bRiojIm4FrgF8pWBTDqBtbIRhGjYjIIuBbwLtV9etFy2MY9WIrBMOond8HXgP8b78XR8Cdqrq1GJEMo3ZshWAYhmEAFmVkGIZh+JhCMAzDMABTCIZhGIaPKQTDMAwDMIVgGIZh+JhCMAzDMABTCIZhGIbP/wdyGJtg2zdR5gAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.scatter(Z, W)\n", "\n", "plt.xlabel('Z')\n", "plt.ylabel('W')\n", "plt.title('BVN generation using multivariate_normal')\n", "\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "estimated covariance of Z and W:\n", "[[0.99653582 0.70792068]\n", " [0.70792068 1.00873687]]\n" ] } ], "source": [ "est_cov = np.cov(Z, W)\n", "print('estimated covariance of Z and W:\\n{}'.format(est_cov))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Example 7.5.10 gives another approach to the BVN generation problem:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "np.random.seed(121393)\n", "\n", "from scipy.stats import norm\n", "\n", "rho = 0.7\n", "tau = np.sqrt(1 - rho**2)\n", "x = norm.rvs(size=10**3)\n", "y = norm.rvs(size=10**3)\n", "z = x\n", "w = rho*x + tau*y" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This gives the $Z$-coordinates in an array `z` and the $W$-coordinates in an array `w`. If we want to put them into one 1000 $\\times$ 2 matrix as we had above, we can use [`numpy.stack([z, w], axis=1)`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.stack.html) to bind the arrays together as columns." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "matrix r2 has shape: (1000, 2)\n" ] } ], "source": [ "# bind z and w into a 1000 x 2 matrix\n", "r2 = np.stack([z, w], axis=1)\n", "print('matrix r2 has shape: {}'.format(r2.shape))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's create another scatter plot now with `z` and `w`, and also check their estimated covariance." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEWCAYAAABmE+CbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJztnX2UHHWZ779P91SSnoB0gnNXMuYFWU00hszICEjuuoZlCYrgyFtk0dXdey/X3cNdk3WzDsiShMVDvLNu8Lieu3JXry9EDBCcBYM30UtYNWvAxJmAkWQVkYSBXaPJRMh0kp6Z5/5RVZ3q6vpV/eqtq1+ezzl9ku6qrvpV1fTz/H7PKzEzBEEQBCGX9QAEQRCExkAUgiAIggBAFIIgCIJgIQpBEARBACAKQRAEQbAQhSAIgiAAEIUgtBhEdBMRba/zOecR0atElK/neesJEf2SiC7LehxCuohCaHCsH2LJEjhHiWgrEc21tt1KRN/z+M5riegUEb2ViD5CRExEa1z7vEhE76rTZaQCES2wrq3D/oyZNzHz5fUcBzMfZOYzmHky7Het5zNpPV/na04aY603loJ2Xte49cwuUOz/BBGdcOx/wOfYy4loBxEdI6JfemxfYG0fJ6L9otCCEYXQHFzFzGcAOAfAfwD4nPX51wBcQkTnuvb/AIBnmPkn1vsjAD5BRK+py2gTopVn3C5+aCkU5+ulrAeVBJaCrlwXgD8H8AsAP/b52i2O7yz02e84gC8BWKPYfj+AYQBnA/gkgIeIqCv8VbQPohCaCGY+AeAhAG+x3r8I4HEAH3Lt+scAvuJ4/yyAHwJYrXMeIjqbiB4lot8S0Y+I6C4i+oFj+yIi+g4RHSGiA0R0g2Pbl4no89ZK5hUiepKIzgvx3f9FRI8R0XEAy4noSiIatsZyiIjWOYZqr47GrNnkO6wZt3Osl1jXcMz69xLHtieI6G+JaKc11u1E9FrFPXmWiN7reN9BRL8more5VypE9CfW/q8Q0S+I6L/r3HePc55n3ae3We/nWOd8V9B5iOhd1irwr4noV0T0MhH1E9F7iOjfrOPe5th/HRE9RESbreP9mIiWKsaVI6IBInqOiH5DRA8Q0WzNy/owgK9yAiUSmPkpZv4aTAXjHuObALwNwFpmLjHzFgDPALg27nlbGmaWVwO/APwSwGXW/zthCvqvOrbfBOBnjvcLAZwC0GW9/wiAHwDoATAGYLb1+YsA3qU45zesVydM5XMIwA+sbTOt938CoAPmj+7XABZb278Mc0VyobV9E4BvhPjuMQDLYE5WZgB4F4Al1vvzYa6Q+q39FwBgAB2OsX/EMdbZAI7CVJgdAG603p9tbX8CwHMA3gSgYL3foLgndwDY5Hh/JYD9XuOwtp0HgAD8PoBxAG9THLcyXsX2/wZToXcC2Abg71xj8DyPdd8mrHEb1nEOA/g6gDMBLAZwAsAbrP3XASgDuM7a/68APA/A8Pg7XAVgF4DXA5gO4AsA7tf4W54PYBLAuT77PGGN89cAdkLxN+r6zmUAfun67P0AnnV99g8APpf1b7qRX5kPQF4BD8j8Ib4KU5hPAHgJwBLH9k4AvwVwifX+UwD+2bHdKSAfAPBp6/+eCgFA3hIMCx2f3eU4xkoA33d95wswZ2KAKdT/ybHtPQ7BqfPdr3rdB8f+9wDYaP1/AfwVwocAPOX6/g8BfMT6/xMAbnds+3MA/1dx3t8F8AqATuv9JgB3qMbh+u4QgI8ptn3Eeq5jjtdzrn0egTm7fRrAdJ97UzkPTIVQApC33p9pjfEix/57cFq5rgOwy7EtB+BlAL/n+Du0FcKzAP7Ase851t+M5/U79vsbAE8E7HORNdbpMFcTrwA4L+A7XgrhQ87rcfw2vpz0b7SVXmIyag76mbkI80dyC4B/IaLXAQAzjwN4EMAfExHBXDF8RXGcOwD8mf1dBV0wZ9OHHJ85/z8fwEVENGa/rHM6j/nvjv+PAzgjxHed5wIRXWQ5Bg8T0TEAHwXgadbxYA6AF1yfvQCgW2OsVTDzz2EKwquIqBPA1TBn2zUQ0buJaJdllhmDqRT9xryLmYuO13mu7f8bwFthzm5PhjjPb/i0o7tk/fsfju0l1/VW7j0zT8GcNHg5t+cD+KbjGT4Lc+b/Oz7XCNSaMmtg5ieZ+RVmPsnMX4G5SnhPwHG9eBWA22f2GpgKRlAgCqGJYOZJZn4Y5o/vPzs2fQXADQD+EObs6luK7+8H8DCA27y2WxyGOWN9veOzuY7/HwLwLy4BdgYz/5nGJeh8121b/jrMGfJcZj4LwD/CNJF47evmJZjCy8k8AKMaY/Xifphmp/cB+KmlJKogoukAtgD4OwC/YynyxxxjDgURnQFzVfRFAOtsW33S57GoPGciysH8G/Bybh8C8G7Xc5zBzMr7SkTLYCqXh0KOiRHtmvYBeAMRnen4bKn1uaBAFEITQSbvAzAL5qzM5vswTQ33wrTXn/I5zHqYNvyi10ZrRvkwTOHTSUSLYM7sbL4F4E1E9CEiMqzX24nozRqXEOW7ZwI4wswniOhCAH/k2HYYwBSANyi++5h1vj+ynMArYfpEPBWmBt8AcDmAP4NidQBgGsyV3GEAE0T0bus7UfksgD3M/F8BbIWpENM4DwBcQETXWM7xVQBOwvQVuPlHAJ8iovkAQERd1t+lHx8GsIWZlTN0IioS0QoimmE9r5sAvBOm78Rr/xwRzYDp8yDre9MAgJn/DcAIgLXW5++H6YPaEjDOtkYUQnPwKBG9CtNX8CkAH2bmykyHTQPpV2HOhr/qdyBmfh5muOpMn91uAXAWTHPK12DOjE9a338FpuD5AMzZ478D+DRM4eRLxO/+OYA7iegVmCavBxzHG4d5P3Za5ouLXef7DYD3Avg4gN8A+GsA72XmXweNVTH+l2H6IC4BsNnnGv/CGudRmArskYBDv4Nq8xDebgnZK2CayQDgLwG8jYhuinieIP4Zpp/HdsRfw8xlj/0+a51ru/VcdsG0/XtiCe0b4GEuIqLbiOjb1lsDpr/Kdir/D5jm0gPWvr9n/Q5s3gnT7PUYzJVfCYAzKfEDAPqs69kA4DpmPux3A9odMmWJIKghok8DeB0zfzjrsQjpQGY47+8y8wezHouQHbJCEGogM1fgfMtEdSGA/wLgm1mPSxCEdOkI3kVoQ86EaSaaA+BXAD4D05wgCEILIyYjQRAEAYCYjARBEASLpjIZvfa1r+UFCxZkPQxBEISmYs+ePb9m5sDCfk2lEBYsWIDdu3dnPQxBEISmgojcGfueiMlIEARBACAKQRAEQbAQhSAIgiAAyFAhWPVFniKivUS0j4jWZzUWQRAEIVun8kkAlzLzq0RkAPgBEX2bmb2KaQmCIAgpk5lCsAqy2YWqDOslWXKCIAgZkWnYKZlN1PfA7Eb1eWZ+0mOfmwHcDADz5s2r7wAFQagbQ8OjGNx2AC+NlTCnWMCaFQvR39sd/EUhMTJ1KlsNX3pgNuK4kIje6rHPvczcx8x9XV2BeRWCIDQhQ8OjuPXhZzA6VgIDGB0r4daHn8HQcNReRkIUGiLKiJnHYPa3vSLjoQiCkAGD2w6gVJ6s+qxUnsTgtgOJnmdoeBTLNjyOcwe2YtmGx0XhuMgyyqiLiIrW/wswG2Xvz2o8giBkx0tjpVCfR0FWIcFkuUI4B8AOInoawI8AfIeZo7Y2FAShiZlTLIT6PAr1WoU0M1lGGT0NoDer8wtCM9AujtY1Kxbi1oefqRLYBSOPNSsWJnaOeqxCmp2mKm4nCO2EbeKwhaRt4gDQckrBvp40ld+cYgGjHsI/yVVIsyMKQRAaFD8TR6spBMBUCmleVz1WIc2OKARBaFDExJEs9ViFNDuiEAShQWkGE0ez+TjSXoU0Ow2RhyAIQi1rVixEwchXfdZIJg4J42w9RCEIQoPS39uNu69Zgu5iAQSgu1jA3dcsaZgZroRxth5iMhKEBqaRTRzi42g9ZIUgCEIk6pFMJtQXUQiCIESi0X0cQnjEZCQIQiT8wjibLfpIMBGFIAhCZLx8HO2UYd1qiMlIEIREkeij5kUUgiAIiSLRR82LKARBEBJFoo+aF/EhCIIQGS/nsRSRa15khSAIQiRUpSsANFyGtbTO1ENWCIIgRMLPebxz4NKGiSiSqCd9RCEIQguTZj5AsziP262vRBxEIQhCixJnZqyjSBqlPHfQWJtFcTUC4kMQhBYlaj6AblnrrEtXDA2PovfO7Vi1ecR3rBL1pI8oBEFoUaLOjHUVSZzy3HGdvLbSOjpertnmHmvWiquZEJORILQoUU06YRRJlPLcSTh5vZSWaqzSOlMfUQiC0KJEzQdI2zeQhJM3aJXjHmsj95VoJMRkJAhNiI7JJapJJ20TSxJOXj/lJOag6GS2QiCiuQC+CuB1AKYA3MvMn81qPILQLIQxuUSZGadtYom7AhkaHsXxkxOe24oFA+uuXiyrgYhkaTKaAPBxZv4xEZ0JYA8RfYeZf5rhmASh4alHXH2aJpY4pS3cytBmVqeBtVeJIohLZgqBmV8G8LL1/1eI6FkA3QBEIQiCD80eVx9nBaJyJndO6/DsyyCO5HA0hFOZiBYA6AXwpMe2mwHcDADz5s2r67gEoREJa3LxE4zNJjR1laGUq4hG5k5lIjoDwBYAq5j5t+7tzHwvM/cxc19XV1f9BygIDcbyRV0g12cqk4tfkpluAlrSxDmvbpKZNOmJRqYKgYgMmMpgEzM/nOVYBCFrdCKHhoZHsWXPKNjxGQG49gJvm7+fYMxKaMY5r1cEFMFUkk6a3ayWFVlGGRGALwJ4lpn/PqtxCEIjoGvi8BKmDGDH/sOex40iGEfHSlgwsBV5Ikwyo7uBiuL193Zj9wtHsGnXwYpSZABb9oyib/7syhgbpc5Ss5GlD2EZgA8BeIaIRqzPbmPmxzIckyBUUS8bu27kUFhhGiQYvbbZTDJX9nErJ/d9Wb6oCzv2H9a6T3GF9Y79h6tWSEDtvZImPdHIMsroB0CNKVQQGoZ6OiZ1BX1YYeolGAnmtRQLBvI5wuSUW7zWUipPYtXmEQxuO4Dli7qw+alDKE+dVhj37TpY2TfoPsUV1jr3SspVRCNzp7IgNCr1tLHrOkvDZhE7s5UBUxnY4n+sVNZSBk5s4V8O+J7ffYpTFA+Q6qVp0hBhp4LQiNTTMak7a44y87WTzJZteNzXRJQ0fvcpTuLb8kVdVT4EoPZeSdhpNEQhCIKCejsmZxi5igDzK8EQVZjWO8ImjfukG2WVdjZ3vXxL9c4TEYUgCArq5Zj0KsdwcmIqkeM6hclZBQNjpdr+AWngdZ+SSJBTRVndt+sgduw/XDmnaiWUhFKs1+oji1WOKARBcOEUTmcVDMwwchgbL6c2Q0tjNuslTIw8wchRoP3fhgBsXNnjWTvID68wVT/hBkBb8AWFy655cK9vqEoSq5Z69WjOohe0KARBcOAWXGOlMgpGHhtX9qT2I0zDV+ElTMqTjFmdBjqndWB0rFTlYPbirIJR5bPQ+Q4A7By4VGs8TsezruBTmfFs/JSdl58hijmmXr6lLJLrJMpIEBxkkb2bRtSMSmiMjZcrkUpBgp2oWmjmiQK/kyfv6bmfcAsj+LyirHRxRjLVo3xGXLKIphKFIAgOVMJpdKyUWo2fNBrS+AmToPaTNkfHy1jz4N6K0LQT1fxQ7VPsNDw/ZwA5hRI5q1D7HXcYrS7dxYK20zmIevVozqIXtCgEQXDgN/tKq/Bb2Lh8nZpHfsIkjMlB199gYwtq9xhP+CgglRI5fmpC2Qlu58CluGdlT801GjlCPletYIw81QjRuOUz4uRR6FKv8zgh1tD6jUJfXx/v3r0762EILYyqAYtNd7HgaSOvF17jKxh5XHtBd03pCOC07d9Zl2j81ASOjicfbVQw8rj7miUAENoRrcLrfjvNWMVOA8zAsVK5UkLDmUUNmEpi8PqlVYJUlZNhn6/ZyoIHQUR7mLkvcD9RCIJQzdDwKFZtHvHcRgCe33BlaucNEkIqQeZ29hKAmy6eh775s2uEc5hIozAUjBzuvub8ihJKAvf9HhoexZoH9yoFfs/67Z6htXkiTDFXKUsvxapSaPa2ZlUKugpBTEaC4KK/t1tpo1bZwv3QLWut4+RUmTTc4p0BbNp1EOse2VcbbTTFUJjtY1EqT1WuISncJrx1j+yrUWblKca6R/ZhaHhUmWcxyVx1XwEozTHt3EtBFIIgeLBmxUIY+Vqp+eoJb7u2Ci9Bv2rzCHrWb686jq4QChNhwoBSQOoaBszcBe1T+pqJcoSqmkq156n+1MuBqrqesVJZW2A7Q1p3DlyK5zdciZ0Dl1Zm/+3cS0EUgiB40N/bjZnTatN0ylMcaqaoiugZK5WrVgC6QkjVICYs3cUCOgMkfbFgYPC6pfhPr0kmzHGKUXEGOyOIZnWa5xm8fimKjs9nhNFECCew/fZt5+J5ohAEQcExxWw0KcHjXAHoCiF35MmsTsNXcM7qNDyjjZYv6sKJsn95jLFSGas2j4Q2AalyEbqLhcqKyTnTd47DWbLj6Hi5xmw2S2Gym9VphBLYfv2nj5+cqPm8XXopiEIQBAVJzBSD9rUVhmrmPzpWqvE72KaOjSt7cKI8hZJCsBeMPK48/5wqhVEsGLj7miV4eM+LiF8tqRYjR7jxornKkNe4LT3XXrW4xpRn5Alrr1rseQ91TFG2j2fBwFas3jxSY5aa1Wk0tUM5DFK6QmgJ4nTwUqFb3M4vOsjrGE5yRDh3YCvmFAuV0FF3iQhVx7KPP7BXGcPfbd2DLXtGq859zJr1p0V5inHfroOY1WmAwBi3lNX0DlMp+SX+qQjb+Ma9zW9/dxiv193snNbRFsoAkLBToQUIyh0AoocNBoWCep3bFubdDoG0/tF9gbH/9hhVYZvFgoGRtZcHXq8dqplW/wM7pyEsBSOP6R250BVX08z90LlHaYYa1wvdsFNZIQhNj04phqhVIoN6D6jKMQOnZ/Z3X7MEw3dcXqVcch5CtVSe9J31j5XKeMvffLsy61Zhm6nSiIopFgylbyWIUnkSM4wcCkZeO2ktbdu9zj1qB2eyjfgQhKZHV/ClISCDjqmKX1cJ/aCZd5AycArQNATZ8VMTnjWGdBkbL1c5xf2oR6mGoHvULs5kG1EIQtOjK/jSEJA6x3zJKoznzEdIA4IZqrl68wiWbXgcyxd1Ra4MqqI8aSa1RT1ujgirLR/GTRfP841IcuYGpIVfGG89FFKjIT4EoemJ4kNIqlaNzrntZKw0+xkbOQLIFNg2zhpHur0MdLAb5+j4RaJg5AmD1y2teR5JPrOkAxAaHfEhCG2DV+SJ3488amtCP4GkaiBjmxxWpxjZkyfCGTM6aoRzqTyJHfsPVxyy9vjjKqY5Vinp/t5uDA2PeiqGWZ0Grjz/nMoz8PKZKPHYLal2kl7H2bJntO1WAipkhSC0Hb13bvec2fpFs6iqjLoFiVPoJllhtGDkcHJiCu6adPZsevXmEc/Zv1eEjN+qJp8jTAYUvvvgxfNwV/+SyrHcxeac47LvzbkDW0OtTtzPIqg6qS6q47iL37WacmiK4nZE9CUi+hUR/STLcQjtw9DwqFIw+zmIdWsN9fd2V+zS9ox4dKyEV09MeNZGcpMnqhRbu2dlD3654Urcs7IHM4x8jTIgy0Q0uO2Asuiel4/Dznb2Mt8HKQMA+Nbelyv/H9x2wLNyqj0uv3H44X4WUfIXdI5r4y5+l1YzpEYna6fylwFckfEYhDbCrw6RLbSGhkfRe+d2LBjYigUDW9GzfrtS8DgFjJ3xumrziGeF0Y4cBXb6OnNGBzau7Kk4VIeGR7Hmob2eSsxe3KsUTlCETFTjwFipHFiDyR6XvV/Y1pc5oiqhrFIoBIQS3jqKqV0qm3qRqUJg5u8BOJLlGIT2wk+ArVmx0FMA+yVSOZVIUOnnUnkKx09OKOvx2OdyzlAHtx2ochSrKE8xZk7r0O6uFVfgBdVgslm9eQS3Dz1TWZWooorcTDJX3Yc1KxZ6hqkywl2LrmJqh8qmXmS9QgiEiG4mot1EtPvw4cNZD0fICJ2eAjqoBFixYFRq4asEsFsgOWfgun2Kx0rlQPORnaA2NDwaSjCNlcoYt/IEXhorYXDbAeV9iivwnDWY3LWCnDCA+3YdrCiFqRDLEudMvb+3W+mDCHMt7uKAKgXVTsloThpeITDzvczcx8x9XV1dWQ9HyADd5jE6qHoNr7t6MQB/4WKXowBMQVIqT2LdI/vQe6fapOSFczavwp4hh23Ic3S8jLFSuar3Qu+d22vuVZRGP05sgdnf243B65ciyD1y366D6Fm/PXTYq/N5qO5XWOHt7IPwmRuW1r2RfSPT8ApBEJLsYBXUuNxPuOSJKsletsN4rFSOFD10rFTGzoFLfZVCqTyZSJy/VxnpOMGFXgJTw6oVuoYRUP08VMo8jPB2rjR71m/H+kfNjnL2SqEdk9GcSB6C0PDoOnR1k5ZU9YlUtfBtJpmxadfBRJK7GKhkE2/+0SEtP0EcSuVJrNo8gsFtB7BmxcLI9YiKBQPrrl5cU100DYw8VQl7nUqnfrjDbZ0KapK5olzaVRkAGSsEIrofwLsAvJaIXgSwlpm/mOWYhMZiaHhUmWHrdujGSVrSyTiGYhxRsZOiOnKUukJwnlOVs6CDs4GNTRx/BAHKpLWZHmWng4oN+hHk54laALGVyFQhMPONWZ5fiE9S5QRUrH90n1J4LV9k+pT8TEq6Y9F1CkdhVqeBzmkdniudtM7pRxzVUypPYv2j+yrPPI4vwk4qWzCw1XN7FBOTHzqKq12ji2zEZCREFupJlRPwO76fDX3HfjPqTKcfsdc1AqfND2GEZNiaQGPjZQzfcXnobN1G5ej4ab9JVB+H0/av6q+gG6Kqw9DwqFb5jHaNLrIRhdDmhBXqOjX9k1p2B9mmbYE/p1jwnH37mZTWPLQXYHhm2frhLhinQ44Itw89A6J4ztxmIU+E6R2kLNWdJ6py3EYtBa6L/fyDjtfO0UU2EmXU5oSJ4HGHf6p+YEHLbt2cgqDj2AI/KPrE6xrLkxxaGdj9iO/qXxJKcEyy2VZSdboEJ8INwRSzss8zAHzmhtM1jmwfkRdBWd26qMyBBFitPiW6yEZWCG2OjrnFRtfO7rfsDrMiUc38gWqBHxR9kpRd2OlQTTKyhjm8GaqedBq5wMY8Tuznr2oD6o5QUl13UrN1v+c/fMfliZyjVRCF0Ob4mVvcdncdE0nQsjtoReIuYe1uEg+Ys7q1V1WHPvpFn+iOPYhSeRK3Pfy0b5vLKHQXC1hwdgE7n4tXxSVqr+N7VvbgLx8YUa5gTk7oH9PInQ4V9aoOaycA6pTi3v3CkURm7EEmRT/SDppoNMRk1OaozC3LF3XVZAerlvbOCp1By26/qpXu823adbAmaeielT0YvuNyz+Yp7oJ0zjo4OpVGjRwh71OGATBbWCapDOx7/dTzR2Mfa5I5sC2lm26rt8GMDrUoCHO9Rp4qylmVAKhT9wkANu06mEjV0agJbUlmyDcLskJoc1TmFlXzeK8GMGFsr6rZml0Kwn0+IDhpyC5I54zlHyuVsebBvZVrXPfIPt8wxjwRVl44F33zZyfSREYHO8nrtoefDu3P8KLbWlXdt+ug1v72PR0aHvU1CYVxho+XpzA0PFpRCl7PS9f0aBeuizsjj5rQlkQ4c7MhCkHw/OGqOnzZ9XxUP6ygJfaaFQs9TQlBAsLvh6gqSFee4sp3gjJzJ5mxZc8o+ubPxs6BS7UT1eIwsvZy3D70TCj7vIocAcdPTmCTpjIAgGsvMO/l6gfU3dwI5srpVIjEuY8/sBerN48oBW8YZeu1ooxixomS0BbGv9YqiEIQPFHN5MN0FfNyGPutSIIERZQf6OhYCecObNWKQS+VJ7F680jFnk4AZk7LY/zUJOYUC3jpWCmxsNHuYgG3Dz2jPZsPYorDJ3Jt2fNiYCmOS86bHdq34WwMdOvDz2D3C0cqrTTDJrK57fxp5764zx3V99CsiA9B8CSK3TVMVzG72qTdCMbOOvZD9UMM+oH6hch67WtbbxjA8VOTuOniedg5cCluumie1jGCKBh5LDi7kJgyiEqpPBUY2fTL38SbDZfKk7hv18GKHT5MIpvX31uShQ6DSKKYXrMhCkHwJKgqqBdxlth21rEKgncY4tDwKMZPqQvSJcH9Tx7C0PBo4Bh1ef2sGbEjiurBrE6jruaRe1b2BP691dOME+U30OyIyUhQEtbuGmWJrRN+CJiz9dWOap3OaJW06wHZvQmSOs/PfnU8keMEYTeOP6tghDYnGXnCleefg/ufPJRoVJUfzmerol5mHLefYuPKnpZWBDayQhASI+wSWzf80MYd+qeKVtGtgVMwctpZwlkUoYsDwcwIfn7DlRhZGy75qlgwsPLtc7Flz6inMjBypBXGGxadsM56mHHaMdzURhSCEIhuqYmwS+yoFUZtm7HKTKA7o509czo23tATOna/GbBDNu1nFaYMxMjay7Fj/2Glsh28fikGr1taOWbY+2fkCcWCt3M5yB9QDzNOPf0UjYaYjAQlQ8OjWP/ovipHYFBURxgzk5/dNyjrNkz2tN8x+nu7sUoRYtsM5Mh0gnvdL+ez8gr39cIW8qpnM8VcEzHmNPvZ41DlLhABg9eZtYxU1V/tc6vCS+P0RNAhip+iVTKaRSEInvjZ5905AUE/BtV2ndDWZRseV9qMdYWcCtvu3J1QaYu4RKlnNMWmMxaAp2Kz+z7PnN5RyfqeZEaxYOD4qYma/I2jx09iaHhU+WzOKhhYtuFxTyHtRBVSe8kbZle1K/Urm1Kv8FKv84fxU2Q51qQRk5HgSZA5xzmL87O3+m1XhZr++tWTFfPUgrMLNSYJZ9by3dcswawITVpyOB215GWXzoJLzpsd6XurN48oEwkBMz/BFnB21ve6qxdj5dvn1uw7Xp7Cmof2et53I0c4fmpCy7auisj61+eOVJUUUfkDsjTbhPVTtJKJSRRCm6PyD+iWntYpVqfa/q29L3se++TEVEXg7HzuSNWsmWBm2DrNFp3Twi90p2DOqJ2N1rOkYOQix/wzwq0s7PuvEtrlSfa872ziSNGgAAAgAElEQVRt8zqWk6HhUeWKy/ZtAP7+gCyzhMP6KVopo1lMRm2M31JXt/R00I/Bb3uUYEbG6dnn0PBoYI2iIJJu0xiFHIC7rznfd5afNGGFFQOYUNRbcnems/+GdM6t8gdknSUcxk+R9ViTRFYIbYzf7F1lRrGbxDjtwF7Yn0fNLvbjpbGSWdDuwb0NIdDjMgVg3SP7YvUnDkux00Auoc48OaLKylInckzn2TdTlnAzjTUIWSG0CV6OXb/Zu26FSFWxOqd93r3dyBOOn4yeXTynWDAL2iVQIbRRqLdiOzZeRvySeiZ24h4QvPIgmCvRZRse943EiVqhNAuaaaxBEDdRk9e+vj7evXt31sNoOrwihgpGHjOMnGdtGb8Cdqrj60YZFTsNvHpiIrIwt8ttr948klqHsaiNZhqVgpHDCY26RXGxQ1ZVpsa4pdOF6BDRHmbuC9xPFELrowrd9IoVT/NHOjQ8GrvbWI6Ad7xhNnb94mhqQvuDF88LrALaLPgp/qQhABtX9qQ6+RCioasQxIfQ4vhGfLgknts/kCS3Dz2D1ZtHfIW4TvjoFAM7nzuS6gx+y57RVEoz2OiW1oiL/TzHIiqDPFGokN45Vvc1rwgd1RiSjsTRzaoXvMlUIRDRFUR0gIh+TkQDWY6lFdGJ+HAyc3pHaiuDoBk3Abjy/HNSFcS6lMqToRrChGFWp4EbL5qLgE6diWA/z7MUZSKCmGTG2qsWw9AYrNNv5FXePI3gAjftXIMoKTJTCESUB/B5AO8G8BYANxLRW7IaTysStlZQWnHTg9sOBJpf7HBSZ42ces2k68nR8TK27Bn17WGcFLbztjwZzX1MMIX7GTP8Y09mTgs2M9YjEqeVEsSyIssoowsB/JyZfwEARPQNAO8D8NMMx9RShBXwXrO1qGUpoozDjm5yfv+8Wx9rKQcvUN/KqXFKcjDM5+tncvrgxfNwV/+SwGPVIxKnlRLEskJLIRDR1wB8D8D3mXl/QufuBnDI8f5FABd5nPtmADcDwLx5yXSsahdUCTOzOg2cKE8pQ0Vtgmq06NZw0S1C56WQbrxobqzOYkaePPstC3rc+vAzKHYaSoewjjKwSbsoXSsliGWF7rr1/wA4B8DniOg5ItpCRB+LeW4ve0DNL5eZ72XmPmbu6+oKbrMonEa1TF971WKt1Pw4ZSmCxuHGyBHGT03UOAP75s+OZdcUZRCPUnkSzGiKxKtWShDLCq0VAjM/TkT/AuDtAJYD+CiAxQA+G+PcLwJwVtd6PYCXYhxPcBG0TA+arcUpS+E1DlWZCSMHlKe4Mgt1Nme//6lDiSVQCdEYK5Vxz8qehk+8aqUEsazQNRn9PwAzAfwQwPcBvJ2ZfxXz3D8C8EYiOhfAKIAPAPijmMcUXMRZpgctwVXbGUDP+u1Yd/XiimlpcNsBHCuVUSwYIALGxsuYUyxg+aIubPIwCZXKk/j6kwfRQsnI2jRaYpy9lPfKF2i0PgBpm6VaHd3V+NMATgF4K4DzAbyViGIZ5ph5AsAtALYBeBbAA8y8L84xW42sY6q9luBO0874qQllSOJYqYw1D+7F7UPPVIUCjpXKOFGewsaVPdg5cCl27D+sjEBqR2Vwz8oeTDWQMgCqK5Q6kTDP1kNLITDzamZ+J4D3A/gNTJ/CWNyTM/NjzPwmZj6PmT8V93itRCP82NxJRsWCAZAZOskw/wVB2Ze4PMW4b9dBXz+DRICcxr6PSTlBZ3UaoVpn+jFqFRR0ImGerYeWQiCiW4hoM4ARAP0AvgQzf0BIiTR/bGFWHs4ko5nTO2qctOVJ9myVGIStCFTCj2DW4MmKZefNTkyY6sJsRvUsX9QVu8+zkSesvWox1qxYmFiyn3tCImGerYduHkIBwN8D2GOZeoSUCfqxRbXdxmn3l+QP3c6eVbXB7JyWx/FT2TWt2fWLo5hiRsHIoVSun1u7VJ7E/U8e0q6jZOTMpveAvzPV3Rs76ticrVMlzLP10DUZDTLzk6IM6ofqR8UAeu/cjjUP7o1kTtJZeahWEEn+0MdKZSwY2IqPP7AXb5t3VpVZyshTpsoAMMs2MKBUBtNTzDTWdSjP6jQweP3SGkV+/OQE1j+6r/L8AGD4jst9Vx26qwjnpCBOmGfW/jHBGylu16D4xe4fHS/XlI/WNSfprDxUvos0eg9PstmucfmiLqVZqhE5NTFVd5OSGzuD2P3Mxkrlip/H+fxUCn1Wp1EpGWLnpRQV9Y+cxwjbatKmEfxjgjfSIKdBccZU65Yf0DHpBC3z1z1S21/YVjZ22GHcEtZe3P/kIfTNnx2r1EI9mVMsZG4rZ5h2/ekdOd9yGPbzW76oyzPr+8rzz6kJ11T10HDP/qOEefqtUiVkNFtkhdDA2A5dXZdg1NaEdgezBQNblZ27bOHX39udSliks+uWFzr3oLNOTmhbMDaCrbxUntTqtjY6VsL9Tx7y3Gb3qHYSdfavgzijGxdZITQBOrWAdG237mxOu4NZkFBxCj/d2kRh8Zvl6qig6UYe4yk7gIsFo5Jw9/kdP0v1XEmjWtXZIaVuYZ9Wkpc4oxsXWSE0AapZfbFgRJq9OUNJO6d1aLWzdCqbsLVh6lXEOmojmLDY9/lnvzpel/PVgyRs+LqOYqk51LjICqEJiFKjRRWW6v5cZ6Y/q9OoOld/b7eyLpGNs38uO953FwtYcHah0gIzT4QbL5qLHfsPx1511MMVPVYqe86m06RYMPDKyQlMppi6HcWG79cr2y+cWWoONS7SU7kFUTkEr72gG1v2jIaqx6/qsTw0PBqp0X234sfvNeZGxe4DvGBga6rnse99mMCCOBCA5zdcqbWv7vOSnsmNgfRUbmNUURz3P3kolMD1M0X193ZHmpH7hRimGdufJLbzM81un857Xy9naxgbvm43PnEUNxdiMmpBVD9C3VDRYsHAyNrLffeJY292myeaaXUAnBacaaVLFAtG1aw6KSe+bbbzyr4Oa8PXFfTiKG4ummNKJvjiduYVO6M1Vbc5FhBxZAvwODgFStjez1lSMPJYvqirkgEcBwI8q8USVStcnaCCWQHPvLtYwMaVPbhnZQ/cbn4CcO0F4SKKdAS9OIqbD1EITY5X1uerJyZqShGEsW7kiHwjRZIQ4E6B0izJaATgbfPOwpY9o4mMmQEMXr+0Jiv46HgZax7aW7n3dk6AU+jPnNaBdVcvxvMbrsTOgUux9qrFymds2/H7e7s9nx3DOxfBj6Qj34TGQExGTY7XD7w8xSgWDMyc3hEqmsjGNi05I0Xsc71kKZ44OGeOtw/FW2nUEwYq0VFJ0F0soL+3G+sfrW0DUp5k3Pbw01UC9YTDzDNWKldF8fT3dmP3C0ewadfBqufjnqUnlRQmkUKtiSiEJkf1Qz5WKlf8AEPDo77lJooFA8dKZeQ8OnWVypNY/+g+vHpyIpEaQ3miqpmjKnu2UUmyZIctqFVVSMfLU5VVgtfzc/ti7upfgr75s32FdJJJYdKdrPUQhdDkBP3AbZOSnyCzfQaqfXTKJjvzDvyYYq4SIvVsFak7Rj+Sam/pzu1Qsf7RfThRnlKe06t/td9xvcqNi61fsBEfQpMTlPWpY+9nxBeUut93z0TzqnZrCUIAZk7LJ2LquvGiubEbzhSMPNZetRhAcLTW0fGy7/Nz38+gbOE0axQJzY+sEFKkHg3Ig2y5jRQHnqPashc3XjTXswJnkjAQq78CAZX7CgCbn4pu5nIn5sXpgOee2es2PxJTj6BCFEJKxOlMFha/H3hahejCYuSAwet7mk4QuTNte9Zv16r95MUHL56Hu/qXVGbxQQ56I09Kv02OUDOzl7LS9ZmEtTKiEFKiUX6cqhaV9aYjf9qsZf9oG0FRBTE6VsKCga0oFgy8d+k5WqWmVdgrC53yIXkizJzWoTzfa2bU+iDavax0PSdhrYr4EFKiUX6cbpux3aKy3tjRSs68iWZirFSObdoqT7FW+RAjR/jMDUt9EwS9tqkihdolW1inPazgj6wQUqKRar57dcOyl9UgoF6BPkfHy54d2doJrQglS1/7mfu8/o7aPYKoUSZhzYysEFIizZrvcRuUO/sh1KVmtANdk8sHL54XWI7BjxyhbquiMIFSOlFV5UnG+kf3Yc2KhZ7jNnLk+XfU7hFE7b5CSoJMFAIRXU9E+4hoiogCS7I2I2n9OJNuUN6oP5atT78cWFPJjylGVdP5mdM6MKvTCGwiH4RbPBeMPC55w2yt7xo5s/eDe6LghZ37MXjd0irFWCwYGLx+qfLvyKns7XIV7YI03olPJv0QiOjNAKYAfAHAXzGzVpMD6YcALNvwuKcZIWrd+WarNBoHZ2+HKNdt5AgrLzSb+TijWHQc5GSZ5uzENve/XkgvgfBIlJE3uv0QMvEhMPOzAEB1SEpqNfzspGF+DM59Z4RoTt+RI0xohl0mkRmcJKXyJD7+wF6s3jyCOcUCrr2gO1SntpUXzsVd/UtqPl+9ecT3e4TTfhpb+E8yg+DvUxDbd3gkxyIeDe9DIKKbiWg3Ee0+fDhcRcZWRGXiKXYa2qYkt9nJXRvfD11lECczWMekEpVJ5sr9uW/XQRw/OaH93S17Rj3vZ5DZTXUfgu5Po5rzhNYlNYVARN8lop94vN4X5jjMfC8z9zFzX1dXV1rDbRpUdlJmaIfc1aP/wKmJqVBO4ZnT8lW+lnoxViprlwZX3U+vZxKXZrB9xw1uEBqP1ExGzHxZWsduZ1SlKlRmi9GxUiUr1t63HqaI8hSD2RRsOsrnRHkKG1eezmRe/+g+z6J6nUYO5SlOpPKqTZgjed27oDGHwVkmo5FNH5IE1ppk4lSunJzoCYhTORFUzmY3Rp5wxvSO2IJLBwKwcWVPlfLyG6Pb6bvmob1Vgj9HwIyOHMZDmLiSJsjRGycLu5mcyEkHNwjpoutUzirs9P1E9CKAdwDYSkTbshhHK6FrtihPMk6WJ1O109vMsRrA7By4FBtX9gTa650mmf7ebgxet7QqlwBAjTIw8hQrXyEMOmYc+3rDhkukaSJKw7QjSWCtSSYKgZm/ycyvZ+bpzPw7zLwii3G0Eu68Bz8hOV6ewt3XLIkci6+DU8DZ5gWdpDSnQHHG1M+c3gEvf3Z58rRpKk1mdRqh8kjO8rm3BGDZebPrkkCWdN6KjSSBtSZSuqKFsEPubCGgu69t0vHqmBaGHJkJYV4lnnWd2CqB4jfzPFYqV0xTSdZIimPPV0VUEwEbb6hf1de0iiy2e5mMVkUUQgsSJICdKwNn3HbcJLUpPi0UnMJGV0gTavsl2ATV9bHPt+bBvTXlqY0cENbtENcWPqby0XB9na5pmXakp3JrIgqhAYmabanj0DRyhHVXL67a33meu69ZEitaxmv2qdt2kmEKmFWbR2oyeWd1GsjBTG93knfU9Vn3yD7PXgVhlUESM12VAit2GjVRX2kK0TSLLEoSWOshCqEOhM0gjhLOpzu7L08xBrcdwO4XjlTV5bfPc+0F3TgRM4rHPfvUNUMRTq8mnBm9gFnbJ2/bpBw4nWBhexXkySwzDVSHjE7viO9a8zKpGHnCqycmKudxP9s0yi6IaUcIgyiElAkr4KPafMPY6UfHSti062BN/H2pPIn7nzwUu4l8jghDw6MVIWfX8fFDp8zFpOfsnyPXu3dep1MJjpXKvs/IuRKzVzBuv4mXSeX4yYkapeWMrEojrl9MO0IYMs1DCEsz5iGEjdc+d2Crp2AkwCxXrUD1PUDfZJM0BSOHCY0ksk4jfm6BbgKc1/emd+Q8Vxdez8hvJebMo/DC79mqTDsS1y8kQUPnIbQTYZ16UcP5VNu7iwVMJaAMopQhLJWntDKKS+WpWCGweaLIjvBSeVJpavIS0H4rsaDuXKpndFbBkLh+oSEQhZAyYQV81Jruft9TnUtXyBeMPG66eF5qeQt2z4KopLX68WpmEySg/bavWbEQRq72mMdPTSjzFiSuX6gnohBSJqyAj9pYx+97qjFcct5srQ5ed1+zBHf1L8HI2stxz8qeyjlaHS9FE3WlBpjP6IwZtW678iSDqDa5Tpy/Qr0RH0IdiBM9klTkifs4yxd1VUUZqfCzYevWT2pk6ulDAPz9CO66T+L8FZKioRvktBtR47WTrCjpHsOyDY8HKoOgGapXSKMfHlGjmdBp5FAqT1UU47f2vlyzj+ranVE7flFGKvzyAiSuX8gaUQgNTFplBwD/7GHdkg39vd3Y/cKRmhBWO4TUHUqqqwzSjoqaNXM6fjpwqXK2P6vTwNqrFvv2LY56/yUvQGhkxIfQwKQVeTI0PKr0AczqNLBxZQ8AszVkUHXMHfsP15hAGKZQjyLSC0Yen7nBrHIaBz8fx6ij3ajX6qZzWkdqM/WoPiJBqAeyQvAgrN0+rcbeKvPCWYV45Q8Gtx1QCutj4+WqPgRBZiqVctKd4Rt5wsxpHThWKtdcS02mb45wxozgXg4FI49rL+j2TbLzM3WlHeoppiGhURGF4CKs3d5vfyBahqgzE9ZtdjFyhOOnTme8RvEr+Am8KQBTrtwBPzOVSmnpmH3yRFj5drNx/dDwKNY/ug+rNo9g1eYRFAsGrr2gGzv2H665f16mHvs+OW35ffNnKwV/qTypHKNf6WpBaGVEIbgIa7dX7b/ukX04OTEVuyaRU1wVCwaIUDNDDutXCOpc5oVKiahs4tde0B0YxTTJjC17THPU5h8dqkpiGyuVsfmpQxi8fmnNdemWY7Dfr1K0F51khpGjmoJ4x09NVEpvAOmtAAWh0RAfgguVoFR9rhKUY6WydtN7J36ZsCcnppTmEl0zx9DwKMZP+Xcu80IVX6+yid/VvwR3X7Mk8Lh2/SSvjOYwdYp2v3DEsytYf2+30h/RXSwo8wJWWf6T24eeSaXBjCA0IrJCcKEyI6gSuMLOtuNkuvqZOXQyWqP2O7CjYFQzZZVNvL+3W6tpjZ9pyet+eJnp7tt1sLLdvRrzi+xZrVg92MdRFQG0FZWsHIRWQlYILlTCSfW5KgtY1cIyTqarPY6oGa1hKqLmiapm/ABqZsqrNo/gzX/zbd9evWtWLAzMavbLlva6HzrX4e7PrIrsCbrfKlVlKx1ZOQithKwQXHT7VJ30QmXPBmojWXRrEvnN4m2naZSZaZjomSnmquqqqkS2klWlVOUjUeUq2Nj+BrcPwWb5oq7I1+Huz+z2CazePIKzCgaMPGkV4XPiVVCvVJ7Exx/YWzmfIDQbohBcREkc8gsjDCu47e3rHtlXU07B2Z4yisAJY95yz5x1hLDKuX1X/xL0zZ/tm93bN382bnv46Zoy2Fv2jKJv/uyqY+peh9fs321uGiuVYeQIszoNpX/GHenlV2p7kjmRPgaCkAVSy8iDRokqCRpHlHwJHR+CV02eMHWLfunTt8EP3d4ROtdBAG66eB7u6q92bPudwy9iyh3+GuQbkT4GQiMhtYxikHTiUFQF4zeOKHWO3HV4vMgTeWbO2g7YoOkDWWOLcv90M7O9zHQLzi7gX587Uhkfw3t14XeOsN3FskxuE4Q0EIWQMkkWqLOPpxLoOvkItpLxmmX7VesM8gXYMBC51lJRYbZxtuR0X4fNsg2PK6OBdMxNtnlJdzJg7/PxB/ZGjvoShEYjkygjIhokov1E9DQRfZOIilmMox74JbqFxRbifqYK3ZlplJo6d/UvwUaNfghRZsdDw6N49YR3foRtl/eL4NFdXURtQORFf283PnPDUuljILQMWa0QvgPgVmaeIKJPA7gVwCcyGkuq6AgqXZOSTrhlmJlpGNOYe4x27X6/2XYYBrcdqMkYdhK0+gma+dsk3XRemtgLrUQmCoGZtzve7gJwXRbjqAdBgsrLpLR68wh2v3CkxiEaNPMOMzMN49dQmb28ylNEnR3rrCqC2lPqRocl7SOSYnVCq9AIiWl/CuDbqo1EdDMR7Sai3YcPH67jsJIhyEThNetnAJt2HawxkfjNvMOUUXaannSSqlRmrx37DydWyllnVRHUnlLKSgtCPFILOyWi7wJ4ncemTzLzP1v7fBJAH4BrWGMgrdBCs9hpgBmVcs9hQhfDOoJV6IZ32vi1fXzeI8Q0SlRVUChplOsUBMEk87BTZr7MbzsRfRjAewH8gY4yaGZUkT1e5a2d6IRbRrFXh228o2ufB6KVD3cqy+kdORwrlWsUp9jlBSF9MvEhENEVMJ3Iv8/M41mMIQtU5iEVc4oFz9m2TsKT3yw9jIAHwtnnw5QPdyuPo+NlFIw8Nq7sEeEvCBmQlQ/hHwCcCeA7RDRCRP+Y5smGhkc9SyPXGz+nqDuMs2DksXxRV6QCakE+Aq/6QH6fh7HPh1l9JBmSKwhCfLKKMvrdep0r6cSwOKhm5qqCdTqzba+VQND3duz3ds6rPgf0I2nCrD7S6hktCEI0GiHKKFUaaRbqF3HU39uNnQOX4vkNV2LnwKXo7+0OFJiqlUBQM580BXGYxC+ViUqyfAUhG1peITTSLDRsaGSQwFQpO79mPjrHjUOYa0wya1gQhPi0fC0jlQkjR4RzB7ZqJWVFCaFUfSdMElOQM1el1OwmOqrvRSnxHYaw9YAky1cQGoOWL3+tUypZFeMeJe7f6zuqUsy641cJzKBSznFLZzdKGXBBEOKhm4fQ8goBqBZsOUVPYq+krLAJXH7fISDxcMqkEtXqfWxBEOqLrkJoeR8CgCqH7ZRCAXqZX6L4H1Tb7LLQSZJmuYZGcsYLglAfWt6H4CZMWGTYBC6/7wDpOLLTKqzWSM54QRDqQ1usEJyEiWyJEgWzZsVCZa+AZgqnlJBQQWg/2k4hhDGzRDHJ9Pd246aL53lmHjdTOKWEhApC+9EWTuUsaIUInVa4BkEQJMpIEARBsMi8/LWQDTKrFwQhKqIQGpCoQr2RCvkJgtB8tJ1TudEJ297SieQOCIIQB1khNBhhGsy4aaXcATF9CUL9EYXQYMQR6lES6ZIkKSEupi9ByAYxGTUYcRLCsswdiGPqciOmL0HIBlEIDUYcoZ5mbaMgkhTirWT6EoRmQkxGDUbcHgFp1TYKIkkhnrXpSxDaFVEIDUhWQj0OSQrxtBv4CILgjZiMhERI0n+RpelLENoZWSEIiZB0O8xmXCUJQrMjCkFIDBHigtDciMlIEARBAJCRQiCivyWip4lohIi2E9GcLMYhCIIgnCarFcIgM5/PzD0AvgXgjozGEYqh4VEs2/A4zh3YimUbHo+UdCUIgtCoZOJDYObfOt7OhNmDvqGRcgqCILQ6mfkQiOhTRHQIwE3wWSEQ0c1EtJuIdh8+fLh+A3Qh5RQEQWh1UlMIRPRdIvqJx+t9AMDMn2TmuQA2AbhFdRxmvpeZ+5i5r6urK63hBiLlFARBaHVSMxkx82Wau34dwFYAa9MaSxJIOQVBEFqdrKKM3uh4ezWA/VmMIwxZVhIVBEGoB1klpm0gooUApgC8AOCjGY1Dm6QzcQVBEBqNrKKMrs3ivHGRTFxBEFoZyVQWBEEQAIhCEARBECxEIQiCIAgARCEIgiAIFqIQBEEQBACiEARBEAQLYm74unIViOgwzLyFOLwWwK8TGE7WyHU0Dq1wDUBrXEcrXAOQ/HXMZ+bA2j9NpRCSgIh2M3Nf1uOIi1xH49AK1wC0xnW0wjUA2V2HmIwEQRAEAKIQBEEQBIt2VAj3Zj2AhJDraBxa4RqA1riOVrgGIKPraDsfgiAIguBNO64QBEEQBA9EIQiCIAgA2lQhENHfEtHTRDRCRNuJaE7WY4oCEQ0S0X7rWr5JRMWsxxQWIrqeiPYR0RQRNV24IBFdQUQHiOjnRDSQ9XiiQERfIqJfEdFPsh5LVIhoLhHtIKJnrb+nj2U9prAQ0QwieoqI9lrXsL7uY2hHHwIRvYaZf2v9/y8AvIWZG75JjxsiuhzA48w8QUSfBgBm/kTGwwoFEb0ZZqOkLwD4K2benfGQtCGiPIB/A/CHAF4E8CMANzLzTzMdWEiI6J0AXgXwVWZ+a9bjiQIRnQPgHGb+MRGdCWAPgP5mehZERABmMvOrRGQA+AGAjzHzrnqNoS1XCLYysJgJoCm1IjNvZ+YJ6+0uAK/PcjxRYOZnmflA1uOIyIUAfs7Mv2DmUwC+AeB9GY8pNMz8PQBHsh5HHJj5ZWb+sfX/VwA8C6CpulmxyavWW8N61VU2taVCAAAi+hQRHQJwE4A7sh5PAvwpgG9nPYg2oxvAIcf7F9FkQqgVIaIFAHoBPJntSMJDRHkiGgHwKwDfYea6XkPLKgQi+i4R/cTj9T4AYOZPMvNcAJsA3JLtaNUEXYe1zycBTMC8loZD5xqaFPL4rClXm60CEZ0BYAuAVS5LQFPAzJPM3ANztX8hEdXVhJdJT+V6wMyXae76dQBbAaxNcTiRCboOIvowgPcC+ANuUIdQiGfRbLwIYK7j/esBvJTRWNoey+6+BcAmZn446/HEgZnHiOgJAFcAqJuzv2VXCH4Q0Rsdb68GsD+rscSBiK4A8AkAVzPzeNbjaUN+BOCNRHQuEU0D8AEAj2Q8prbEcsh+EcCzzPz3WY8nCkTUZUcKElEBwGWos2xq1yijLQAWwoxueQHAR5l5NNtRhYeIfg5gOoDfWB/tarZoKSJ6P4DPAegCMAZghJlXZDsqfYjoPQDuAZAH8CVm/lTGQwoNEd0P4F0wSy7/B4C1zPzFTAcVEiL6zwC+D+AZmL9rALiNmR/LblThIKLzAXwF5t9SDsADzHxnXcfQjgpBEARBqKUtTUaCIAhCLaIQBEEQBACiEARBEAQLUQiCIAgCAFEIgiAIgoUoBEEQBAGAKARBEATBQhSCIMSAiD5q9dUYIaLniWhH1mMShKhIYpogJIBVR+dxAP+TmR/NejyCEAVZIQhCMnwWZrMiUQZC09Ky1U4FoV4Q0UcAzEcDl1EXBB3EZCQIMSCiC2AWJPs9Zj6a9XgEIQ5iMhKEeNwCYDaAHZZj+Z+yHpAgREVWCIIgCAd6kTgAAAA5SURBVAIAWSEIgiAIFqIQBEEQBACiEARBEAQLUQiCIAgCAFEIgiAIgoUoBEEQBAGAKARBEATB4v8DcdaG5Bhddm0AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.scatter(z, w)\n", "\n", "plt.xlabel('z')\n", "plt.ylabel('w')\n", "plt.title('BVN generation via Example 7.5.10')\n", "\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "estimated covariance of z and w:\n", "[[1.03751418 0.69336422]\n", " [0.69336422 0.94874494]]\n" ] } ], "source": [ "est_cov2 = np.cov(z, w)\n", "print('estimated covariance of z and w:\\n{}'.format(est_cov2))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Cauchy\n", "\n", "We can work with the Cauchy distribution introduced in Example 7.1.24 using the three functions `pdf`, `cdf`, and `rvs` in [`scipy.stats.cauchy`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.cauchy.html). To shift and/or scale the distribution use the `loc` and `scale` parameters. Specifically, `cauchy.pdf(x, loc, scale)` is identically equivalent to `cauchy.pdf(y) / scale` with `y = (x - loc) / scale`." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "from scipy.stats import cauchy\n", "\n", "# to learn more about scipy.stats.cauchy, un-comment out the following line\n", "#print(cauchy.__doc__)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For an amusing demonstration of the very heavy tails of the Cauchy distribution, try creating a histogram of 1000 simulated values of the Cauchy distribution:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAz4AAAF1CAYAAAApyz9eAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3XmYnWV9//H3d2ayQFiyEAjZWMOmyDYFLIoLIAEtYAUFpKCl8rNCpVKtYC1a1LoWKhapKLZoVUREiBoXiijVCiUgoqwJAZKQQAYCCZB1Zr6/P86TyZmTM8kkmSVznvfruuaaZ7mfZ+6TEM58zn3f3ycyE0mSJElqZE2D3QFJkiRJ6m8GH0mSJEkNz+AjSZIkqeEZfCRJkiQ1PIOPJEmSpIZn8JEkSZLU8Aw+GnIi4omIOHaw+7GpIuJdEfHrDZz/SUScM5B9kiQ1pqH6Xin1J4OPtkhEnBkRsyLipYhYVPzy/prB7ldvRMTrI6Kz6PuLEfFIRLy76vzDEfGXda67MCJmFduviIifR8TzEfFCRNwTESduTn8y84TMvK4X/c6I2HtzfoYkaeA16ntlROxevCe9VHw9ExE/iojjau7xRESsqGr3UkRMHJxXpDIz+GizRcRFwL8C/wzsAkwFvgycPJj92kQLM3M7YAfgA8BXI2Lf4tx1wNl1rvmL4hzAD4Fbqbz+nYH3A8v6tcf9KCJaBrsPktRIGvC98sNU3isPqDo/ujh/EJX3xB9ExLtq7vFnmbld1dfCAem5VMXgo80SETsClwHnZ+ZNmflyZq7JzB9m5oeKNodHxG+LkZBFEfFvETG8OLf2U6KWqnv+MiL+qmr/PRHxUPEJ04MRcWhVFw6OiPsjYmlEfDciRhbX/DEi/qzqHsMi4tmIOHhDrycrZgJLgFcVh78JvCYidqu63/7F+e9ExE7AHsBXM3N18fWbzOxxOltxjy8UI0SPR8QJ9V5/ROwdEb8qXt+zEfHd4vgdRfPfF5+YvaPqz2pORCyJiBnVn6RFxJuKT+iWRsSXi/uu/TnviojfRMQVEbEE+HhE7BURv4iI54qf/a2IGF11vyci4kPFn//LEXFtROxSfIL5YkT8d0SM2dCfgSSVQYO+V94MPA8cUOf805n5ReDjwGcjwt8ztVXxP0htrlcDI4EfbKBNB5VRlJ2K9scA7+vNzSPiNCr/4zybyidMJwHPVTV5OzCdSvB4FfCu4vg3gLOq2p0ILMrM+zby85oi4qSir3MAMnMBcDuVEZ61zgZmZuazRX/mAP8VEadExC69eGlHAI8UP+dzwLUREXXafQL4OTAGmAx8qejT0cX5g4pPzL4bEW8EPk3lz2RX4Eng+uJ17QTcCFwCjCt+9p/W6dNcKiNWnwKiuN9EYH9gCpW/i2pvA44D9gH+DPgJ8JHidTVRGfmSpLJrxPfKtwKjgT9soOlNVN5T9t1AG2nAGXy0ucYBz2Zme08NMvOezLwzM9sz8wngK8Drenn/vwI+l5l3F58wzcnMJ6vOX5mZCzNzCZXpZms/pfov4MSI2KHY/wsqIzc9mRgRLwArqLwxXZSZv6s6f11xD4pPrt5ZHCMzE3gD8ATwL8CiiLgjIqZt4Oc9mZlfzcyO4j67Upn6UGsNsBswMTNXbmQU6Z3A1zPz3sxcRSXkvDoidqfyZvZA8UljO3Al8HTN9Qsz80vF39OK4s/61sxclZltwOWs//f2pcx8JjOfAv4HuCszf1f8/B8Ah2ygv5JUFo32Xvks8DHgLzLzkQ20XzuNbWzVsZuLUa0XIuLmXr4+qU8ZfLS5ngN2ig2sCYmIfaKyyPHpiFhGZX7zTr28/xTgsQ2cr/7lfTmwHUAxZ/g3wNuK6VknAN/awH0WZuZoKp+UXQm8seb8TcCuEXEk8HpgW+DHa09m5oLMvCAz96ISVF6m8knaRvudmcuLze3qtPt7KiMv/xcRD0SdIgtVJlIZ5Vl735eo/P1MKs7NrzqXwIKa6+dX70TEzhFxfUQ8Vfy9/Rfr/709U7W9os5+vdckSWXTUO+VmTk2Mw/OzOs30q9JxfclVcdOKe4xOjNP2cj1Ur8w+Ghz/RZYCWzof15XAw8D0zJzBypTodZO63q5+L5tVfsJVdvzgb02s2/XURnCPw34bTEqsUHFSMWHgQMj4pSq48upTBU7m8onYtdn5uoe7jEfuAp45Wb2u/peT2fmezJzIvD/gC9Hz5XcFlIJXQBExCgqnzI+BSyiMlVu7bmo3l/742r2P10ce1Xx93YW6/7eJEm911DvlZvgrcBiKtOrpa2GwUebJTOXApcCVxXrW7YtFkeeEBGfK5ptT6XC2UsRsR/w11XXt1H5xfysiGguRjSq/+f9NeCDEXFYVOwdVUUGNuJm4FDgQjY8+lL7mlZTmbJ2ac2p64B3UFnX0lVuOiLGRMQ/FX1rKtbT/CVwZ29/Zk8i4rSIWBtQnqcSRDqK/WeAPauafxt4d0QcHBEjqHxaeFcxZeLHFGGu+MTxfLq/adazPfAS8EJETAI+tKWvR5LKqBHfKzekKHRzAZXpcJdkZmdf3FfqKwYfbbbMvBy4CPgo0Eblk6cLqPzPFOCDwJnAi8BXge/W3OI9VH6pfg54BfC/Vff+HpWF9t8urr+Z7nOFN9SvFcD3qSzmvGkTX9bXganV1W6AO4ClwFOZeXfV8dXA7sB/U3nT+iOwinWLR7fEnwB3RcRLwAzgwsx8vDj3ceC6Yp702zPzNuAfqbzmRVTeFE8HKIownEalkMJzVKrwzCr62ZN/ovJmuJRKcNrUP0NJUqFB3ytrvRARL1MpeHAicFpmfn0L7yn1uahM+ZcaS0RcCuyTmWdttHGJFAUaFgDvzMzbB7s/kqTB43ulysYRHzWciBgLnAtcM9h92RpExPERMbqYBrd27vgWT8eTJA1dvleqjAw+aigR8R4q0wh+kpl3bKx9SbyaStWfZ6k8c+eUYoqDNKRFxPSoPJx3TkRcXOf8eyPiDxFxX0T8OoonzUfloZAriuP3RcS/D3zvpcHje6XKyqlukqQhJyKagUepPEh3AXA3cEZmPljVZofMXFZsnwS8LzOnF8+4+lFmbnEFRknS0NGrEZ+I+EDxLJE/RsR3ImJkROwREXdFxOyI+G5EDC/ajij25xTnd+/PFyBJKqXDgTmZObeoyHg9cHJ1g7WhpzCK9Uu3S5JKZKPBpyhn+36gtfh0rJlKxajPAldk5jQq5XbPLS45F3g+M/cGrijaSZLUlybR/eG7C1j30MQuEXF+RDxGpbLh+6tO7RERv4uIX0XEa/u3q5KkrUGPTxKu026biFhD5SFai6g84f7M4vx1VErsXk3lE7ePF8dvBP4tIiI3MKdup512yt13331T+y5J6kP33HPPs5k5frD70Uv1Hqq73vtMZl5F5RkqZ1IpJ3wOlfewqZn5XEQcBtwcEa+oGSEiIs4DzgMYNWrUYfvtt19fvwZJ0ibY0vepjQafzHwqIr4AzANWAD8H7gFeyMz2oln1J21dn8JlZntELKXyFPlnq+9b/YYydepUZs2atbmvQZLUByLiycHuwyZYAEyp2p8MLNxA++upfDhHZq6ieJZVZt5TjAjtQ+UZV10y8xqKiletra3p+5QkDa4tfZ/qzVS3MVRGcfYAJlKZJ31CnaZrP2nr7adw12Rma2a2jh8/VD5glCRtJe4GphXrTYdTmYI9o7pBREyr2n0zMLs4Pr4ojkBE7AlMA+YOSK8lSYOmN1PdjgUez8w2gIi4CfhTYHREtBSjPtWftK39FG5BRLQAOwJL+rznkqTSKmYUXAD8jMra069n5gMRcRkwKzNnABdExLHAGiprUc8pLj8auCwi2oEO4L2Z6fuUJDW43gSfecCREbEtlalux1CZDnA7cCqV6QPnALcU7WcU+78tzv9iQ+t7JEnaHJk5E5hZc+zSqu0Le7ju+8D3+7d3kqStzUanumXmXVSKFNwL/KG45hrgw8BFETGHyhqea4tLrgXGFccvAtZ7qJwkSZIkDaReVXXLzI8BH6s5PJfKcxRq264ETtvyrkmSJElS3+jVA0wlSZIkaSgz+EiSJElqeAYfSZIkSQ3P4CNJkiSp4fWquIEk9daKFSvo6OgY7G6UVkSwzTbb0NTk51qSJFUz+EjqE48//jizH32QzjWraWlpHuzulFYmrOnoZOLkqRx44EE0N/t3IUkSGHwk9YF58+Yx55H7OeygA9hxxx0Guzult2rVKh58aDb33jOLPzn8iMHujiRJWwXnQkjaYk/MncMr99/b0LOVGDFiBAe96gDannmKVatWDXZ3JEnaKjjiI9Xx8MQpm33tfgvn92FPhoYXXniOceP2H+xuqEpTUxM77jCKpUuXsvPOOw92dyRJGnSO+EjaIpkJiYvpt0LNTU10dnYOdjckSdoq+JuKpNIYNmoCcx57fFD7cNmnPs/Zf3n+oPZBkqQyMvhI6hd779/KpN1fycsvv9x17Nr//BbHTH9rr64/ZvpbufY/v7XZP39Lr5ckSY3F4COp37R3dPClL39tsLshSZJk8JHUf/7uwvdx+Rev5oUXltY9/7933s2Rrz2ecbtO48jXHs//3nk3AP/48U/z69/cxYUXfYTRO+/J+y+6BICHH5nN9Le8nZ0n78crDj6K733/lrr37el6gNtuv4P9X/Vqxk/al7/5wMWVNUqF/7ju2xx46GsZP2lfTjzpdJ6cV79QxZtPPoOr/v3abscOPeKN/OCWHwPwgQ9+lD32OZSxE/bm8KPexK9/c2fd+/zqjt+w+7RDuh3be/9WbvvFHQB0dnbyuS98iX1feQS7TNmfM/7iPSxZ8jwAK1eu5Oy/PJ9dpuzPThP34cjXHs8zz7TV/TmSJMngI6kfHXboQbzutX/K5V+8er1zS5Y8z8lvO4sL/vqveGb+Q/zt37yXk992Fs89t4RPfPwSXnPUEXzx8n/mhcVzufLyT/Pyyy9zwp+9g9Pf/lYWPvFHvvkfV/M3H7iEBx58eL1717t+rZk/uZXf3vFT7rnzNm68aQY/v/V2AG754U/47Beu5IZvX8uiJx/gNUcdwVnv+uu6r+v0t7+V737v5q79Bx96hHnzF3Di9GMBaD3sYGb99r9ZvOBhznj7Wzn9rPewcuXKTf7z+9KXv8YtP/oJt/3sB8ybcx+jR4/m/R+ohLhvfOsGli1bxuOP3MMz8x/iqis/xzbbjNzknyFJUlkYfCT1q4999ENc9e/X0tb2bLfjM3/63+y91x6cdeZptLS0cPrb38q+++zNj2b+vO59fvyTW9lttym86+wzaGlp4dBDXsVbT34zN938o03qz4f+7m8YPXpHpk6ZzOuPPorf3/8AAF+99hv8/Qffz/777UNLSwsXf+hCfn//A3VHfU456UR+f/8fu85957s3ccpJJzJixAgA3nnGqYwbN5aWlhY+cOFfs2r1ah559LFN6ifA177+TT7xsUuYPGkiI0aM4NKPfJDv3/wj2tvbGdbSwnNLnmfOY0/Q3NzMYYccxA47bL/JP0OSpLIw+EjqV698xf6cOP04PvcvX+p2fNGip9lt6uRux3abOpmFi56ue5958xbwf3ffy04T9+n6+s53v8/TzyzepP5M2GXdM2222WYbXiqKL8ybv4CLPvTRrnvvPHk/MpOnFq7fn+23344Tpx/LDcWozw033syZ73hb1/krvng1Bx76WsbtOo2dJu7D0qXLeO65JZvUT4An5y3g1DPe3dWnAw97Lc3NzTyzuI2zzjyNNx37es561/9j6l4HcfE/XMaaNWs2+WdIklQWPsBUUr/72Ec/xOFHHccH3v/ermO77jqBJ4s1MWvNm/8Ubzr2DQBERLdzkydP4ujXvJqf/uiGXv3M2us3ZvKkSVz8ob/lzNPftvHGwDtOeyuf+Od/4TWveTUrVqzk9a87CoBf/+ZOPn/5v/GzH9/IKw7Yl6amJsZP2rfbWqK1th21LctXrOja7+jooO3Z57r2p0yeyDVXX8FRrz68bh/+8SMf5B8/8kGeeHIeJ/35O9lnn735y3PO3JSXLUlSaTjiI6nf7b3XHpz2tpP5t6vXFQQ44fhjmD1nLt/57k20t7dzw40389DDj/LmE44DYJedx/P44092tX/zCccxe85j/Ne3v8eaNWtYs2YNd9/zOx56+NG6P7P2+o0576/O5nNfuLJrzdDSpcu48aYZPbY/4fhjmDd/Af/0ic9x2ttO7nqA64svvkRLSwvjdxpHe3s7n/z0v7Bs2Yt177HP3nuxcuUqZv70VtasWcM/f/YKVq1a3XX+PeeezaX/9OmuKXVtbc8y40c/BeCXv/o1f/jjQ3R0dLDD9tvT0jKMZh8iK0lSj3yXlDQgPnrJRbz88vKu/XHjxnLzjd/kiiv/nV2m7M8XrvgyN9/4TXbaaRwAF7zvPdx0848YP2lf/vaD/8D222/HzBnf5YYbb2bq3gcxec9X8ZGPfrJbUKhWe/3GnHLSiXzwogt45znvZeyEvTn4T17PT3/+ix7bjxgxglNOOpHbbr+D09+x7tlEbzruDRz/pjdywMF/yl77tTJyxAimTJ5Y9x477rgDX7ri0/y/9/0du+19MKO23ZbJk3btOv/+89/DW048nhNPOp0xu+zFa97wZv7v7nsBePqZNk4/668YO2FvDjz0tRz92lfzzjNO3ejrlCSprKLe9IuB1tramrNmzRrsbkhdHp44ZbOv3W9h/RLIjSoz+eHNN3Li8UcPdldUY9Y997P7tAOZMGFCr9pHxD2Z2drP3RqSfJ+SpMG3pe9TjvhIkiRJangGH0mSJEkNz+AjSZLUwNZ0dPLgwqU8+9Kqwe6KNKgsZy1JktTAbrx7Po8+Xaku+b5j9mb89iMHuUfS4HDER1K/OKj1aH51x2/65F6/uuM37D7tkD6516aaN38Bo3fek46ODgCOmf5Wrv3Pbw1KXyRpc6wNPQCzn3lpEHsiDS6Dj6R+8ftZd/C6o4/arGuHjZrAnMce7+MeVVz3zet53bEn9br91CmTeWHxXJqbm/ulP5I0oLaCar7SYDH4SJIkSWp4Bh9J/WLv/Vu57Rd3AHDZpz7PGX/xHt71VxcwZpe9OKj1aGbde1/d697wplMAOOzINzJ65z254cabu85d8cWrmbjbK5iy56v4z298p+v4qlWr+PtLPs6e+x7GpN1fyfve//esWLFivXs/9PCjnH/hh7nzrlmM3nlPdpq4DwAzf3orra8+lrET9maPfQ7lsk99vuuaJ56cx7BRE2hvb1/vfnMee5w3Hn8K43adxoSpB3Dm2edtxp+UJPWf2uc1Ot6jMjP4SBoQP/zxz3nHqafw7MJHecuJx3PhRR+p2+72n1eCzj13/oIXFs/l7adWgtDTzyxm6bJlPDnnPq758uW8/6KP8PzzLwBwyUc/yew5c5n12//m4T/cycKFi/jkpy9f797777cPV33xsxx5RCsvLJ7LswsfBWDUttvyH1/9Es8ufJRbvv9ffOWr13HLD3+y0df08cs+y3HHvJ62px7hiUfv5X3vPXez/mwkqb90dHaPOp0mH5WYwUfSgDjq1YdzwvRjaW5u5p1nnMr9f3hwk64fNmwYH73k7xg2bBgnTD+W7UZtyyOzHyMzufY//4svfPYyxo4dw/bbb8eHP3Rht5GijXnd0Udx4Cv3p6mpiVcdeADvOO2t3PE/v93odS3DWnhy3gIWLnqakSNH8po/PWKTXpMk9bfaJT2dJh+V2EaDT0TsGxH3VX0ti4i/jYixEXFrRMwuvo8p2kdEXBkRcyLi/og4tP9fhqSt3S677Ny1ve2227By5cq608d6Mm7sGFpa1lXg33bbbXj5pZdpa3uW5ctXcMRr3sROE/dhp4n78JZTzqDt2ed6fe+77r6XY0/4c3bd7QDG7TqNa679Bs8+t2Sj133mk5eSmfzp0SdwUOvR/Md13+71z5SkgdBRk3xq96Uy2WjwycxHMvPgzDwYOAxYDvwAuBi4LTOnAbcV+wAnANOKr/OAq/uj45IEsNNO49hmm234/axf8ezCR3l24aM8t2g2LyyeW7d9RKx37Ox3/zVvOfF4Hn/kXp5bNJvzzj17vXnx9UyYsDNfuepfmPfY77nqys/zNx+4pN+q0UnS5uisXeNj8FGJbepUt2OAxzLzSeBk4Lri+HXAKcX2ycA3suJOYHRE7NonvZVUCrvsPJ7HH3+yV22bmpo4913v5IMfvpTFi9sAeGrhIn5+6+093vuppxayevXqrmMvvvgSY8eMZuTIkfzfrHu5/oabevWzb7xpBgueWgjAmNE7EhE0NzuDWNLWo3ZqmzPdVGab+g59OrC2lNIumbkIoPi+dh7LJGB+1TULimOS1Cv/+A8f5C/Pez87TdyH733/lo22//QnP8pee+7Ba97wZsZO2Jvpb3k7j8x+rG7bN7z+NRyw/75M3vNVTJh6AABf+tfP8E+f/BxjdtmLT336ck59W++e8zPrnvs46nUnMnrnPfnzt5/D5Z//BHvsvlvvX6gk9bPaoFNb7EAqk+jtkGdEDAcWAq/IzGci4oXMHF11/vnMHBMRPwY+nZm/Lo7fBvx9Zt5Tc7/zqEyFY+rUqYc9+WTvPt2VBsLDE6ds9rX7LZy/8UYNJDP54c03cuLxRw92V1Rj1j33s/u0A5kwYUKv2kfEPZnZ2s/dGpJaW1tz1qxZg90NaZO9sHw1X/z5o137h+85lhNeNXEQeyRtvi19n9qUEZ8TgHsz85li/5m1U9iK74uL4wuA6t8aJ1MJTN1k5jWZ2ZqZrePHj9/0nkvaajhnXIMhIqZHxCNFMZ2L65x/b0T8oSjM8+uIOKDq3CXFdY9ExPED23Np4NSO8DjiozLblOBzBuumuQHMAM4pts8Bbqk6fnZR3e1IYOnaKXGSGk9E0NTcTEdHx2B3RTXaOzppbm4e7G70i4hoBq6i8qHcAcAZ1cGm8O3MPLAozvM54PLi2gOoTN1+BTAd+HJxP6nhrFfO2tyjEutV8ImIbYHjgOoVv58BjouI2cW5zxTHZwJzgTnAV4H39VlvJW2Vxo7diba23pePVv9rb29n6bKXGT169MYbD02HA3Myc25mrgaup1Jcp0tmLqvaHcW6h9afDFyfmasy83Eq71eHD0CfpQFXW9Wtdl8qk5aNN4HMXA6Mqzn2HJUqb7VtEzi/T3onaUjYa9q+/P6eO4kIxo8fR1OTlc0G07JlL/Lgw7OZvNseDBs2bLC701/qFdJZ7wmyEXE+cBEwHHhj1bV31ly7XhGemrWofdJpaaDVVnVzarLKrFfBR5I2ZMKECXDYkcx59GHuue8hmprWf1aOBkYmjNx2FFOm7sG+++472N3pT/X+I1vvN7rMvAq4KiLOBD5KZWp2b6+9BrgGKsUNtqi30iBZ7wGmnYPUEWkrYPCR1CcmTJjAhAkTyEzX+wyipqamsoy49aqQTpXrWfdA7U29Vhqyatf0ONVNZWbwkdSnIoKWFv/Xon53NzAtIvYAnqJSrODM6gYRMS0zZxe7bwbWbs8Avh0RlwMTgWnA/w1Ir6UBtv4DTA0+Ki9/O5EkDTmZ2R4RFwA/A5qBr2fmAxFxGTArM2cAF0TEscAa4HmKSqRFuxuAB4F24PzMdJhSDak26FjOWmVm8JEkDUmZOZNKJdHqY5dWbV+4gWs/BXyq/3onbR1qB3gc8FGZlWIiuCRJUhmtV9zA5KMSM/hIkiQ1KMtZS+sYfCRJkhrU+g8wHaSOSFsBg48kSVKDqs05tSNAUpkYfCRJkhqUM9ukdQw+kiRJDap2TY/P8VGZGXwkSZIaVO3MNme6qcwMPpIkSQ0qsaqbtJbBR5IkqUH5AFNpHYOPJElSg1ov+AxON6StgsFHkiSpQa3/HB+jj8rL4CNJktSgnOomrWPwkSRJalC1xQwsbqAyM/hIkiQ1qNqYYzlrlZnBR5IkqUHVrumpLW8tlYnBR5IkqUG5xkdax+AjSZLUoNZf4zNIHZG2AgYfSZKkBrX+Gh+Tj8rL4CNJktSgaosZGHtUZgYfSZKkBmU5a2kdg48kSVKDsriBtI7BR5IkqUHVlq92jY/KzOAjSZLUoNZb42PuUYkZfCRJkhrUelPdLG+gEjP4SJIkNajaYgadnYPUEWkrYPCRJElqUPXGd6zsprIy+EiSJDWoeiHH3KOy6lXwiYjREXFjRDwcEQ9FxKsjYmxE3BoRs4vvY4q2ERFXRsSciLg/Ig7t35cgSZKkemqLG1SOmXxUTr0d8fki8NPM3A84CHgIuBi4LTOnAbcV+wAnANOKr/OAq/u0x5IkSdps5h6V1UaDT0TsABwNXAuQmasz8wXgZOC6otl1wCnF9snAN7LiTmB0ROza5z2XJEnSBnXWGfKxspvKqjcjPnsCbcB/RMTvIuJrETEK2CUzFwEU33cu2k8C5lddv6A41k1EnBcRsyJiVltb2xa9CEmSJK2vXsSpN/1NKoPeBJ8W4FDg6sw8BHiZddPa6ok6x9b7J5aZ12Rma2a2jh8/vledlSRJUu/VW8/jVDeVVW+CzwJgQWbeVezfSCUIPbN2ClvxfXFV+ylV108GFvZNdyVJktRrdUKO5axVVhsNPpn5NDA/IvYtDh0DPAjMAM4pjp0D3FJszwDOLqq7HQksXTslTpIkSQOn3oiPVd1UVi29bPc3wLciYjgwF3g3ldB0Q0ScC8wDTivazgROBOYAy4u2kiRJGmB1H2A64L2Qtg69Cj6ZeR/QWufUMXXaJnD+FvZLkiRJW6jetDaLG6isevscH0mSJA0x9Wa1ucZHZWXwkSQNSRExPSIeiYg5EbFetdGIuCgiHoyI+yPitojYrepcR0TcV3zNGNieSwOn3uiOuUdl1ds1PpIkbTUiohm4CjiOSjXRuyNiRmY+WNXsd0BrZi6PiL8GPge8ozi3IjMPHtBOS4Og3sNKLW6gsnLER5I0FB0OzMnMuZm5GrgeOLm6QWbenpnLi907qTxeQSqVulPdBr4b0lbB4CNJGoomAfOr9hcUx3pyLvCTqv2RETErIu6MiFP6o4PS1qDeeh7X+KisnOomSRqKos6xur/NRcRZVCqTvq7q8NTMXBgRewK/iIg/ZOZjNdedB5wHMHXq1L7ptTTAXOMjreOJb2CDAAAgAElEQVSIjyRpKFoATKnanwwsrG0UEccC/wCclJmr1h7PzIXF97nAL4FDaq/NzGsyszUzW8ePH9+3vZcGSL2M4xoflZXBR5I0FN0NTIuIPYqHa58OdKvOFhGHAF+hEnoWVx0fExEjiu2dgKOA6qIIUsOoO63N3KOScqqbJGnIycz2iLgA+BnQDHw9Mx+IiMuAWZk5A/g8sB3wvYgAmJeZJwH7A1+JiE4qHwB+pqYanNQw6uUeH2CqsjL4SJKGpMycCcysOXZp1faxPVz3v8CB/ds7aetQb1pbvRLXUhk41U2SJKlB1V/jM+DdkLYKBh9JkqQGZTlraR2DjyRJUoOqW9vA3KOSMvhIkiQ1qHprfCxnrbIy+EiSJDUoR3ykdQw+kiRJDar+Y3xMPiong48kSVKDqhdyrOqmsjL4SJIkNah6Icepbiorg48kSVKDspy1tI7BR5IkqUHVyzhWdVNZGXwkSZIaVL2IY+xRWRl8JEmSGlS90R0HfFRWBh9JkqQGVf85PiYflZPBR5IkqUHVCzmWs1ZZGXwkSZIaVP01PiYflZPBR5IkqUG5xkdax+AjSZLUqHyAqdTF4CNJktSg6q3n8Tk+KiuDjyRJUoOqt57H3KOyMvhIkiQ1qHojPpazVlkZfCRJkhqV5aylLr0KPhHxRET8ISLui4hZxbGxEXFrRMwuvo8pjkdEXBkRcyLi/og4tD9fgCRJkuqznLW0zqaM+LwhMw/OzNZi/2LgtsycBtxW7AOcAEwrvs4Dru6rzkqSJKn36hc3GPh+SFuDLZnqdjJwXbF9HXBK1fFvZMWdwOiI2HULfo4kSZI2Q731PK7xUVn1Nvgk8POIuCciziuO7ZKZiwCK7zsXxycB86uuXVAckyRJ0gCql3HMPSqrll62OyozF0bEzsCtEfHwBtpGnWPr/RMrAtR5AFOnTu1lNyRJktRb9dbz+BwflVWvRnwyc2HxfTHwA+Bw4Jm1U9iK74uL5guAKVWXTwYW1rnnNZnZmpmt48eP3/xXIEmSpLo6O9c/Zu5RWW00+ETEqIjYfu028Cbgj8AM4Jyi2TnALcX2DODsorrbkcDStVPiJEmSNDB6WstjVTeVVW+muu0C/CAi1rb/dmb+NCLuBm6IiHOBecBpRfuZwInAHGA58O4+77UkSZI2qKeRHau6qaw2Gnwycy5wUJ3jzwHH1DmewPl90jtJkiRtlp7W8jjVTWW1JeWsJUmStJXqKeBYzlplZfCRJElqQD2t5TH2qKwMPpIkSQ2op7U8lrNWWRl8JEmSGpBrfKTuDD6SJEkNyDU+UncGH0mSpAbU43N8zD0qKYOPJElSA3KNj9SdwUeSJKkB9RRvzD0qK4OPJGlIiojpEfFIRMyJiIvrnL8oIh6MiPsj4raI2K3q3DkRMbv4Omdgey4NjB6nug1wP6SthcFHkjTkREQzcBVwAnAAcEZEHFDT7HdAa2a+CrgR+Fxx7VjgY8ARwOHAxyJizED1XRooPU1pc6qbysrgI0kaig4H5mTm3MxcDVwPnFzdIDNvz8zlxe6dwORi+3jg1sxckpnPA7cC0weo39KA6WmNj7lHZWXwkSQNRZOA+VX7C4pjPTkX+MmmXBsR50XErIiY1dbWtoXdlQaB5aylbgw+kqShKOocq/vbXEScBbQCn9+UazPzmsxszczW8ePHb3ZHpcHS81S3Ae6ItJUw+EiShqIFwJSq/cnAwtpGEXEs8A/ASZm5alOulYa6noJPWt5AJWXwkSQNRXcD0yJij4gYDpwOzKhuEBGHAF+hEnoWV536GfCmiBhTFDV4U3FMaiyu8ZG6aRnsDkiStKkysz0iLqASWJqBr2fmAxFxGTArM2dQmdq2HfC9iACYl5knZeaSiPgElfAEcFlmLhmElyH1Kx9gKnVn8JEkDUmZOROYWXPs0qrtYzdw7deBr/df76TB19OUNnOPysqpbpIkSQ3IER+pO4OPJElSA7JstdSdwUeSJKkBOeIjdWfwkSRJakA9lrM296ikDD6SJEmNqCrgRNVjew0+KiuDjyRJUgOqHvFpalqXfHyAqcrK4CNJktSAquNNS1Xw6Wntj9ToDD6SJEkNqNuIT9VcN6u9qawMPpIkSQ2oOt80V091M/eopAw+kiRJDainER/LWausDD6SJEkNyBEfqTuDjyRJUgOqXsvT3OSIj2TwkSRJakDV1duafI6PZPCRJElqRNXP62lpaqo6LpWTwUeSJKkBVY/sVOUep7qptHodfCKiOSJ+FxE/Kvb3iIi7ImJ2RHw3IoYXx0cU+3OK87v3T9clSZLUE4sbSN1tyojPhcBDVfufBa7IzGnA88C5xfFzgeczc2/giqKdJEmSBpAPMJW661XwiYjJwJuBrxX7AbwRuLFoch1wSrF9crFPcf6Yor0kSZIGSGePVd0GozfS4OvtiM+/An8PdBb744AXMrO92F8ATCq2JwHzAYrzS4v23UTEeRExKyJmtbW1bWb3JUmSVE+PU90sb6CS2mjwiYi3AIsz857qw3WaZi/OrTuQeU1mtmZm6/jx43vVWUmSJPVOz1PdBqM30uBr6UWbo4CTIuJEYCSwA5URoNER0VKM6kwGFhbtFwBTgAUR0QLsCCzp855LkiSpRxY3kLrb6IhPZl6SmZMzc3fgdOAXmflO4Hbg1KLZOcAtxfaMYp/i/C/SVXSSJEkDKnsY8bGctcpqS57j82HgooiYQ2UNz7XF8WuBccXxi4CLt6yLkiRJ2lTV8aal2xofqZx6M9WtS2b+EvhlsT0XOLxOm5XAaX3QN0mSJG2mzh6nuhl9VE5bMuIjSZKkrVS3qW6u8ZEMPpIkSY2op+IGrvFRWRl8JEmSGpDlrKXuDD6SJEkNqPsan3XbPsBUZWXwkSRJakDVa3yam9b9ytfZORi9kQafwUeSJKkBVY/rVC3xqZxzvptKyOAjSZLUgKrX+EQEVct8XOejUjL4SJIkNaDqcBN0L3BgZTeVkcFHkiSpAWVNVbduIz6D0B9psBl8JEmSGlB1VbemqEx3W8s1Piojg48kSVIDql3jU13goNPcoxIy+EiSJDWi6jU+6434DEJ/pEFm8JEkSWpAnbVrfKrOOdVNZWTwkSRJakCdNSM+VnVT2Rl8JElDUkRMj4hHImJORFxc5/zREXFvRLRHxKk15zoi4r7ia8bA9VoaSNVrfLCqm0qvZbA7IEnSpoqIZuAq4DhgAXB3RMzIzAerms0D3gV8sM4tVmTmwf3eUWkQda/q1r2ctcUNVEYGH0nSUHQ4MCcz5wJExPXAyUBX8MnMJ4pznYPRQWmwZW3wqV7l41Q3lZBT3SRJQ9EkYH7V/oLiWG+NjIhZEXFnRJxSr0FEnFe0mdXW1rYlfZUGRe06HstZq+wMPpKkoSjqHNuUX+WmZmYrcCbwrxGx13o3y7wmM1szs3X8+PGb209p0GS3qm7dy1lb3EBlZPCRJA1FC4ApVfuTgYW9vTgzFxbf5wK/BA7py85JW4MNrfEx9qiMDD6SpKHobmBaROwREcOB04FeVWeLiDERMaLY3gk4iqq1QVKj6DaqE/gcH5WewUeSNORkZjtwAfAz4CHghsx8ICIui4iTACLiTyJiAXAa8JWIeKC4fH9gVkT8Hrgd+ExNNTip4TRFdHuOj7lHZWRVN0nSkJSZM4GZNccurdq+m8oUuNrr/hc4sN87KA2y6hGfgJpy1iYflY8jPpIkSQ2o2xqfJkd8JIOPJElSA6pex9McQXOTVd1UbgYfSZKkBtTRWVvOet25Th/koxIy+EiSJDWgzs51201N3Ud8Osw9KiGDjyRJUgPqrJnq1uQDTFVyBh9JkqQG1K2qW01xA6e6qYwMPpIkSQ2oW1W32jU+5h6VkMFHkiSpAXUvbmBVN8ngI0mS1ICqy1k3ucZH2njwiYiREfF/EfH7iHggIv6pOL5HRNwVEbMj4rsRMbw4PqLYn1Oc371/X4IkSZJqdVQHH9f4SL0a8VkFvDEzDwIOBqZHxJHAZ4ErMnMa8DxwbtH+XOD5zNwbuKJoJ0mSpAFUPajTHNBc9VufuUdltNHgkxUvFbvDiq8E3gjcWBy/Djil2D652Kc4f0xE9XI6SZIk9bfaNT7hVDeVXK/W+EREc0TcBywGbgUeA17IzPaiyQJgUrE9CZgPUJxfCoyrc8/zImJWRMxqa2vbslchSZKkbrqVs65Z49PhkI9KqFfBJzM7MvNgYDJwOLB/vWbF93qjO+v968rMazKzNTNbx48f39v+SpIkqReqs01zU6Wk9bpzBh+VzyZVdcvMF4BfAkcCoyOipTg1GVhYbC8ApgAU53cElvRFZyVJktQ7ndlzOWtzj8qoN1XdxkfE6GJ7G+BY4CHgduDUotk5wC3F9oxin+L8LzL95yVJkjSQOjewxsepbiqjlo03YVfguohophKUbsjMH0XEg8D1EfFJ4HfAtUX7a4FvRsQcKiM9p/dDvyVJktSDzOw2qhOBDzBV6W00+GTm/cAhdY7PpbLep/b4SuC0PumdJEmSNln1iE5zU6xX3MABH5XRJq3xkSRJ0tavdrQHaoobmHxUQgYfSZKkBtOR3Ud8gO7lrJ3qphIy+EiSJDWY2mf4ADR1q+pm8FH5GHwkSZIaTPVUtuYi+DR3e47PQPdIGnwGH0mSpAZTHWyait/2LGetsjP4SJIkNZjah5dCTTlrg49KyOAjSZLUYDo61w8+0W2qm8FH5WPwkSRJajDVuWbtQE+zz/FRyRl8JEmSGkxnvXLW1VPdHPFRCRl8JEmSGkz1Gp6uctZhOWuVm8FHkiSpwXR7gGlX8Kk63znQPZIGn8FHkiSpwVSv4Vk70ONUN5WdwUeSJKnBdHuAadPaB5gafFRuBh9JkqQGU+85Pj7AVGVn8JEkSWow3YJP8dtek8/xUckZfCRJkhpMZ7fn+NQrZz3QPZIGn8FHkiSpwVSv8VkbfJotZ62SM/hIkiQ1mLoPMO1Wztrgo/Ix+EiSJDWYjZaz9jk+KiGDjyRpSIqI6RHxSETMiYiL65w/OiLujYj2iDi15tw5ETG7+Dpn4HotDYyOOlPdmixnrZIz+EiShpyIaAauAk4ADgDOiIgDaprNA94FfLvm2rHAx4AjgMOBj0XEmP7uszSQ6pWz9gGmKjuDjyRpKDocmJOZczNzNXA9cHJ1g8x8IjPvB2on9RwP3JqZSzLzeeBWYPpAdFoaKNW5pt4aH5f4qIwMPpKkoWgSML9qf0FxrL+vlYaE6qluXWt8fICpSs7gI0kaiqLOsd7+JterayPivIiYFRGz2traNqlz0mDrVtWtzhofy1mrjAw+kqShaAEwpWp/MrCwL6/NzGsyszUzW8ePH7/ZHZUGQ7c1Pk3rr/HpMPiohAw+kqSh6G5gWkTsERHDgdOBGb289mfAmyJiTFHU4E3FMalhVE9lW7vGp6XJqW4qN4OPJGnIycx24AIqgeUh4IbMfCAiLouIkwAi4k8iYgFwGvCViHiguHYJ8Akq4elu4LLimNQw6gUf1/io7FoGuwOSJG2OzJwJzKw5dmnV9t1UprHVu/brwNf7tYPSIGqvCjZrR3pamqPueaksHPGRJElqMBud6taRFjhQ6Rh8JEmSGkx7x7rHV60NPhFBU9Vvfk53U9kYfCRJkhpMvRGfynZT3TZSGWw0+ETElIi4PSIeiogHIuLC4vjYiLg1ImYX38cUxyMiroyIORFxf0Qc2t8vQpIkSet0dFvjs+7XvWZLWqvEejPi0w78XWbuDxwJnB8RBwAXA7dl5jTgtmIf4ARgWvF1HnB1n/dakiRJPWrvYcSnep1Pe4fBR+Wy0eCTmYsy895i+0UqZUMnAScD1xXNrgNOKbZPBr6RFXcCoyNi1z7vuSRJkurqNuLTXD/4ONVNZbNJa3wiYnfgEOAuYJfMXASVcATsXDSbBMyvumxBcUySJEkDoNsan6rn9zQZfFRivQ4+EbEd8H3gbzNz2Yaa1jm23r+siDgvImZFxKy2trbedkOSJEkb0d65flU3qJnqVtVGKoNeBZ+IGEYl9HwrM28qDj+zdgpb8X1xcXwBMKXq8snAwtp7ZuY1mdmama3jx4/f3P5LkiSpRvepbut+3WuxqptKrDdV3QK4FngoMy+vOjUDOKfYPge4per42UV1tyOBpWunxEmSJKn/VRcu6FbOutmpbiqvll60OQr4C+APEXFfcewjwGeAGyLiXGAecFpxbiZwIjAHWA68u097LEmSpA3qXs66+jk+1VPdDD4ql40Gn8z8NfXX7QAcU6d9AudvYb8kSZK0mXp6gKlV3VRmm1TVTZIkSVu/np7j02zwUYkZfCRJkhpMb6a6GXxUNgYfSZKkBtObqW6Ws1bZGHwkSZIaTEe35/g01d2urvwmlYHBR5IkqcG092KqW6dT3VQyBh9JkqQG09NUN8tZq8wMPpIkSQ2kszPJItM0NUFTj2t8DD4qF4OPJElSA6kuWtAU3R/F2L2qm8UNVC4GH0mSpAbSbX1Pc/df9SxnrTIz+EiSJDWQzh4KG9TuO9VNZWPwkSRJaiDtPRQ2gO4jQB2Ws1bJGHwkSZIayIaCj1XdVGYtg90BqdE8PHHKZl+738L5fdgTSVIZtbevK1owrGaNT/V+e4fFDVQujvhIkiQ1kDUd1cGn+4hP9f5qg49KxuAjSZLUQKoDzfCWnkd81hh8VDIGH0mSpAaypqpoQe1Ut+ogtMbiBioZg48kSVIDWd3LNT7V7aQyMPhIkiQ1kOopbMNrg0+LU91UXgYfSZKkBtKtuIFrfKQuBh9JkqQGsuGpbuuquq1xqptKxuAjSZLUQLoXN+heznp4s8UNVF4GH0mSpAbS/Tk+3X/Va24KmopDHZ1JR6fhR+Vh8JEkSWog1VPYap/jExFWdlNpGXwkSZIayIZGfGqPWeBAZWLwkSRJaiCruwWfWO+8wUdlZfCRJElqIBsf8Ym6baVGZ/CRJA1JETE9Ih6JiDkRcXGd8yMi4rvF+bsiYvfi+O4RsSIi7iu+/n2g+y71p+pqbbVrfGqPWdJaZdIy2B2QJGlTRUQzcBVwHLAAuDsiZmTmg1XNzgWez8y9I+J04LPAO4pzj2XmwQPaaWmAbOg5PrXHLGmtMnHER5I0FB0OzMnMuZm5GrgeOLmmzcnAdcX2jcAxEbH+ggepwXSb6lZnxKdbVTenuqlEDD6SpKFoEjC/an9Bcaxum8xsB5YC44pze0TE7yLiVxHx2no/ICLOi4hZETGrra2tb3sv9aNVa9aFmZF1gs/IYeuOrWrvGJA+SVsDg48kaSiqN3JTO2enpzaLgKmZeQhwEfDtiNhhvYaZ12Rma2a2jh8/fos7LA2EzGRlVZgZMax5vTYjq46tXO2Ij8rD4CNJGooWAFOq9icDC3tqExEtwI7AksxclZnPAWTmPcBjwD793mNpALR3Jh3Fup3m5qi7xqc6DK1c44iPymOjwScivh4RiyPij1XHxkbErRExu/g+pjgeEXFlUUHn/og4tD87L0kqrbuBaRGxR0QMB04HZtS0mQGcU2yfCvwiMzMixhfFEYiIPYFpwNwB6rfUr1asXhdktqkz2gM1Iz5OdVOJ9GbE5z+B6TXHLgZuy8xpwG3FPsAJVN5ApgHnAVf3TTclSVqnWLNzAfAz4CHghsx8ICIui4iTimbXAuMiYg6VKW1r36uOBu6PiN9TKXrw3sxcMrCvQOof1Wt2RtRZ3wPdA9HK1QYflcdGy1ln5h1rn31Q5WTg9cX2dcAvgQ8Xx7+RmQncGRGjI2LXzFzUVx2WJAkgM2cCM2uOXVq1vRI4rc513we+3+8dlAZBdZAZ2cOIz4huxQ1c46Py2Nw1PrusDTPF952L472psgNYLUeSJKmvrawKMtsM3/hUtxWu8VGJ9HVxg95U2akctFqOJElSn6ouVjCipYfg0+JUN5XT5gafZyJiV4Di++LieG+q7EiSJKkfdJ/qVv/XPKe6qaw2N/hUV8o5B7il6vjZRXW3I4Glru+RJEkaGNVT3Xpa47ONU91UUhstbhAR36FSyGCniFgAfAz4DHBDRJwLzGPd4tGZwInAHGA58O5+6LMkSZLqqJ7q1nNxg3XHV63pIDOJqLdaQWosvanqdkYPp46p0zaB87e0U5IkSdp0y1e1d233VNyguSkYMayJVWs6yayM+mw7fKO/EkpDXl8XN5AkSdIgeakq+Gw3oucwM6rq3MtV10iNzOAjSZLUIKpDzKgNBJ/qUPTSSoOPysHgI0mS1CC6jfiMdMRHqmbwkSRJagCZyfJV64ob9HrEx+CjkjD4SJIkNYAVazro6Kw8N37EsCaGNff8a171aJAjPioLg48kSVID6O36ntrzrvFRWRh8JEmSGkB1gNlQRbfa8051U1kYfCRJkhrA0hVrura3Hzlsg22rzy+ruk5qZAYfSZKkBvDC8tVd22NGbTj4jK46//zyNVSeQS81NoOPJElSA3j+5XUjN6O3Hb7BttsMa2bEsMqvgWvaO1m+umOD7aVGYPCRJElqAEtXVI34bCT4RES3cFQ9WiQ1KoOPJElSA6ge8dlx2w1PdQMYXdXmheWu81HjM/hIkiQNcWs6OnlxZSW8RMCO2/Qm+Kwb8XnupVX91jdpa2HwkSRJGuLalq1ibX2CsaOG07KBh5euNX77EV3bi5cZfNT4DD6SJElD3OIXV3Zt77zDyF5dU92urep6qVEZfCRJkoa4xcuqg8+IDbRcZ+eqEZ9nX1pFe0dnn/dL2poYfCRJkoa4p5du+ojPiGHNXQUOOjuh7UWnu6mxGXwkSZKGsM7O5KnnV3TtTxqzba+vrW47f8nyPu2XtLUx+EiSJA1hTy9dyer2yjS1HbYZ1quKbmtNHWfwUXkYfCRJkoawJ557uWu7Osj0xpSx69o/8ezL5NrScFIDMvhIkiQNYbOffrFre8/x223StbvsMJJthjcD8NLKdha9YHU3NS6DjyRJ0hC1fFU785asG/GZtsumBZ+mpmDaLtt37T+0aFmf9U3a2hh8JEmShqg/LFhKZ1GFevLYbdhuZO/X96y1367rgs/v5z9PZ6fT3dSYDD6SJElDUGdnMuuJJV37B08ds1n32WfCDmw3sgWAF1e08+gzL27kCmloMvhIkiQNQQ8uXMqzxbN3Rgxr4hWTdtys+zQ3RbfQ9D+PtFnkQA2pZbA7IPWnhydOGewuSJLU51a3d/CLhxZ37R+x5zhGDmve7Pu17jGWOx97lvaOZOELK7h/wQscNGXzRpCkrZUjPpIkSUNIZvKzPz7N8y+vBmDk8GaO2GvcFt1zx22GceReO3Xt//T+RV33lxqFwUeSJGmIyEz+59E27n3i+a5j0w+cwLbDt3wSz2um7cTobSvFEVau6eRbv32CF1eu2eL7SlsLg48kSdIQsLq9gx/9fiG3V01xe+XkHXnV5NF9cv8Rw5o59U+m0NwUADz30mq+9qu5zK96QKo0lBl8JEmStmLtHZ3c88QSvvyLOd1GevYYP4qTD5lERPTZz5o0Zlv+vHUyTcVviMtWrOE/fv04M373FM++tKrPfo40GCxuIG1FtqQYw34L5/dhTyRJg2nVmg7mL1nOw4uW8dCiZSxf1dHt/Csn78hJh0yipbnvP8M+YOKODDuiiZvuWcDK1R1kwu+efJ7fPfk8e4wfxX677sC0XbZn9LbD+jR0Sf2tX4JPREwHvgg0A1/LzM/0x8+RJJXXxt5rImIE8A3gMOA54B2Z+URx7hLgXKADeH9m/mwAuy4BlefwLF/dzksr21m6Yg1tL66i7cVVLF62kmeWraReRelthjcz/cBdOXDyjv0aOqbtsj3/7/V7MfP+Rcx+et1zfR5ve5nH217mJyxi1IgWdh09kgk7jmTMtsMZM6rytd2Iln4JZNKW6vPgExHNwFXAccAC4O6ImJGZD/b1z9LQYVnprZsjTRpqevlecy7wfGbuHRGnA58F3hERBwCnA68AJgL/HRH7ZGb3j9TV8DKTTOjMJAGqtjOL80W7zoSOzqSjM2nv7Kzazq7trnMdle3VHZ2sXNPJqjUdrGrvZFV78X1NBy+tauflVe10dvaur9tv08KRe47jsN3HMmILylZvitHbDufMI3djbttL3PXYc8x+5sVuYezlVe3MeeYl5jzz0nrXjhjWxKgRLYwa0cI2w5oZ3tLE8JYmRhTfh7c0M6KliZbmoKWpieam6Ppqqdqu7DfR1BQE0BQQEURAU8S6fSr7UexL9fTHiM/hwJzM/9/eucbYVZVh+HntFWhpZVpJU25TUwkkJlKLkij8gShtlBEVU//YRBNiAonEmFhDQhoTfqDBH0YjqZEIBK1XYmMw4N3EWKBgKW3KpaUYxpZWWymUwrQNnz/2mrJnevbMmXb27Zz3SXZnnW/tffY7X78571ln77NWvAggaSMwBHjgY0yJeHBp+oxuvGYIWJ/avwS+p+wd0RCwMSJGgD2SdqXn+0fRyQ4eGeH+v+8ZE4vCB/nwqR35N45FS0R2Wjyym/Uki/YZ1VHc372O7s438XMUn2/C0034fBHZ75kNWPIDmBQji5H263bAURcSnH/uXC5edA6XLTmXiwbOru0N/bLF81i2eB6H3zzO86+8xvOvvM7woaO8dbw4iSPH32bk+DEOHal+Suxs8DN2QPSu3IBIgqwna5887uQ/Y/t18nnf2VnKHZPr07j+oj7R+f/yTP6Luzm22xrqtFdXz1/we3Vzgiqqu4yBz1Ig/xHwMPDhEs4zhn76xNpvcE2T6Kd6bNtrRY/Tjdec3CciTkg6DAyk+OZxxy6d6GQjJ95mz388s5WZfs6aPYP5c2cyf+4sBubNZtG8OSyaP4clC886owVJy2DBWbO4cnCAKwcHiAgOvnGMvf97k4NHRnj16DEOvXGMV48e5+ix7q9klcHo4PedIXYXI2rTF5Qx8Ok0YDul4iTdDNycHo5I2l6Clu6Y+vB6EfDfEjYKFRoAAAaySURBVJSUSds0t00vtE+z9U6V3n+tuLRuAVOgG68p2ue0fGr9je+vz6dOj7bVn/WWT9s0t00vtE9z2/SekU+VMfAZBvIfAV8A7B2/U0RsADYASNoSEStL0FIKbdML7dPcNr3QPs3WWz5t0yxpS90apkA3XjO6z7CkmcAC4FCXx7bap6B9mq23fNqmuW16oX2a26j3TI4vY8qNJ4DlkgYlzSb7AummEs5jjDGmf+nGazYBa1P7s8CfIvuCyCZgjaQ5kgaB5cDjFek2xhhTE9N+xSfdR30r8AjZFKP3RsSO6T6PMcaY/qXIayR9E9gSEZuAHwEPpMkLDpENjkj7/ZxsIoQTwC2e0c0YY3qfUtbxiYiHgYencMiGMnSUSNv0Qvs0t00vtE+z9ZZP2zS3Sm8nr4mIO3Ltt4CbCo69E7hzCqdrVW4SbdNsveXTNs1t0wvt09xXejXZVJXGGGOMMcYY03a8rK4xxhhjjDGm56l84CPpJkk7JL0taeW4vm9I2iXpOUkfz8WvT7FdktZVrTmn42eStqbtJUlbU/wSSW/m+u6pS2MeSesl/Tuna3Wur2Ou60bStyU9K2mbpIckLUzxRuYYmlOfRUi6UNKfJe1Mf3tfSfHC+mgC6W/smaRtS4qdJ+n3kl5IP99dt04ASZfm8rhV0muSbmtajiXdK+lAfvmAopwq47uprrdJWlGf8mqxT1VL27zKPjX92KfKxz6VyFY0rm4DLiObg/svwMpc/HLgaWAOMAjsJvvC6ozUXgbMTvtcXrXuDr/H3cAdqX0JsL1uTR00rge+1iHeMdd1603aPgbMTO27gLsanuNG1uc4jUuAFak9H3g+1UDH+mjKBrwELBoX+xawLrXXjdZHk7ZUE68AFzctx8A1wIr831JRToHVwO/I1ry5Cnisbv0V5sk+Va3OVnmVfaoUjfap6muiL32q8is+EbEzIp7r0DUEbIyIkYjYA+wCPpS2XRHxYkQcAzamfWtDkoDPAT+tU8cZUJTr2omIRyPiRHq4mWx9jSbTuPocT0Tsi4inUvt1YCeTrFLfYIaA+1L7PuBTNWop4lpgd0T8q24h44mIv5HNbpanKKdDwP2RsRlYKGlJNUrrxT7VGBrpVfap6cc+VTl961NN+o7PUuDl3OPhFCuK18nVwP6IeCEXG5T0T0l/lXR1XcI6cGu6/Hdv7nJrE3PaiS+SjeRHaWKO25JLILsVA7gCeCyFOtVHUwjgUUlPSro5xc6PiH2QGSXwntrUFbOGsW82m5xjKM5pq2q7IuxT5dFWr7JPTTP2qUroW58qZeAj6Q+StnfYJvqEQR1iMUG8FLrU/nnGFsw+4KKIuAL4KvATSeeWpXEKen8AvBf4QNJ49+hhHZ6qsun9usmxpNvJ1td4MIVqy/Ek1JrLqSBpHvAr4LaIeI3i+mgKH4mIFcAq4BZJ19QtaDKULaR5A/CLFGp6jieiNbV9Otinqn0NbZtX2afqwT5VPv3uU2Wt43PdaRw2DFyYe3wBsDe1i+LTzmTaJc0EPg18MHfMCDCS2k9K2g28D9hSls7cubvKtaQfAr9NDyfKdel0keO1wCeAayPdxFlnjieh1lx2i6RZZGbyYET8GiAi9uf68/XRCCJib/p5QNJDZLdr7Je0JCL2pcvZB2oVeSqrgKdGc9v0HCeKctqK2j5d7FPVvoa2zavsU9Vjn6qMvvapJt3qtglYI2mOpEFgOfA48ASwXNJgGqWuSfvWxXXAsxExPBqQtFjSjNReRqb9xZr0nURj73O8ERidIaMo17Uj6Xrg68ANEXE0F29kjmlefZ6CJJGtYL8zIr6TixfVR+1IOkfS/NE22ZeJt5Pldm3abS3wm3oUFjLmU/Ym5zhHUU43AV9QxlXA4dFbDfoY+1QJtM2r7FPTj32qUvrbp7qdZWG6NrKkDpN9KrIfeCTXdzvZzCPPAaty8dVkM3zsBm6vWvM4/T8Gvjwu9hlgB9lMKU8Bn6xTY07XA8AzwLZUHEsmy3XdG9mXV18GtqbtnibnOGlrTH0W6Pso2aXfbbm8rp6oPureyGYfejptO0bzCgwAfwReSD/Pq1trTvPZwEFgQS7WqByTmd0+4Hh6Hf5SUU7JbiH4fqrrZ8jNbtbrm32qcr2t8ir7VCn67FPVaO57n1I60BhjjDHGGGN6libd6maMMcYYY4wxpeCBjzHGGGOMMabn8cDHGGOMMcYY0/N44GOMMcYYY4zpeTzwMcYYY4wxxvQ8HvgYY4wxxhhjeh4PfIwxxhhjjDE9jwc+xhhjjDHGmJ7n/8FqiNLQk168AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "np.random.seed(196418)\n", "\n", "fig = plt.figure(figsize=(14, 6))\n", "\n", "# create frozen instance of a cauchy distribution\n", "cauch = cauchy()\n", "\n", "# generate 1000 random simulated values from cauchy\n", "rv = cauch.rvs(size=1000)\n", "\n", "ax1 = fig.add_subplot(121)\n", "ax1.hist(rv, bins=50, color='#d7191c')\n", "ax1.set_xlim([-100.0, 100.0])\n", "ax1.set_title('Cauchy RVS histogram')\n", "txt = 'Note the values\\nin the tails'\n", "props = dict(boxstyle='round', facecolor='wheat', alpha=0.3)\n", "ax1.text(0.69, 0.95, txt, transform=ax1.transAxes, fontsize=12,\n", " verticalalignment='top', bbox=props)\n", "\n", "# create a sequqnce of 1000 values from -100 to 100 for x\n", "x = np.linspace(-100, 100, 1000)\n", "# obtain corresponding cauchy PDF values for y\n", "y = cauch.pdf(x)\n", "\n", "ax2 = fig.add_subplot(122)\n", "ax2.plot(x, y, lw=3, alpha=0.6, color='#2c7bb6', label='cauchy pdf')\n", "ax2.set_xlim([-100.0, 100.0])\n", "ax2.set_ylim([0.0, 0.35])\n", "ax2.set_title('Cauchy PDF')\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Due to extreme values in the tails of the distribution, this histogram looks nothing like the PDF of the distribution from which it was generated." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "----\n", "\n", "© Blitzstein, Joseph K.; Hwang, Jessica. Introduction to Probability (Chapman & Hall/CRC Texts in Statistical Science)." ] } ], "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" } }, "nbformat": 4, "nbformat_minor": 2 }