{ "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, 2nd Edition](https://www.crcpress.com/Introduction-to-Probability-Second-Edition/Blitzstein-Hwang/p/book/9781138369917), 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+CbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJztnX2YXXV94D/fmdyEO4EyicbWDISw\n2CY1BDJlViNsnzWsLbQIjiCmCr71he3uurummBqUp7wUl+yTtdju9tku3da+iJhg6FQMbtCHWCs1\naOJMgEiiohC4oEaTQcgMyc3Md/8450zOnDm/83Jf5tx75/t5nnmSe8/b955z7+/7+31fRVUxDMMw\njK6iBTAMwzBaA1MIhmEYBmAKwTAMw/AxhWAYhmEAphAMwzAMH1MIhmEYBmAKwSgQEblORB4qWo4A\nEflVETlYtBxpiMj7ReRrCdub9jna5R7NFiKyXERUROYVLUsjMIXQQETkaREZF5GXReSoiOwQkbP9\nbTeJyFdjjnm1iJwQkfP9H7qKyMbIPs+JyJtn6WM0hbgfjqreo6q/XqRcYVT1n1V1RdFy5MW/r68L\nXjfzc+Q5d5riMloPUwiN50pVPR14LfAj4H/67/89cLGInBvZ/7eAx1X1Cf/1EeAjIvJzsyJtgxCR\n7qJlMJpLJ8yCO+EzNBNTCE1CVV8BPge83n/9HPAw8J7Iru8F/jb0+kng68CGLNcRkVeJyAMi8jMR\n+aaI3BGelYnIShH5kogcEZGDIvLO0La/EZE/91cyL4nIoyJyXo5j/7eIPCgix4B1InKFiAz7sjwr\nIreGRA1WR6P+CupN0RmkiFzsf4YX/X8vDm37ioj8sYg84sv6kIi82nFPZsxMw7NoEflNEfm2f56K\niHzYf//NIvJc6JinReTDIvKYL9NWETkttP0PReQFEXleRH43OlOPXP8r/rP5F//zP+A/u3tCz265\nv++M1ZR//O/GnDe4r/v8864Pfw4R2SQin4sc86ci8mf+/z8gIk/69+L7IvLvQ/u92V+dfkREfgh8\nKuYebRKRp/zjvy0ib/ff/2XgL4A3+XKN+u8vEJH/ISKHRORHIvIXIlKOu2cxctwoIj/27/kHQtvP\nFJG/E5HDIvKMiNwsIl3+tvf735m7ROQIcGvkvVH/c1/sv/+sf433hc6f9L3uLFTV/hr0BzwNvMX/\nfw/eQP93oe3XAd8NvV4BnACW+K/fD3wNWAOMAov9958D3uy45mf9vx485fMs8DV/20L/9QeAecCv\nAD8BVvnb/wZvRfIGf/s9wGdzHPsicAnexOI04M3Aav/1BXgrpEF//+WAAvNCsr8/JOti4CiewpwH\nvMt//Sp/+1eAp4BfAsr+682OezJ13tB7CrzO//8LwK/6/18E/Ir//zcDz0We5zeApb58TwK/72+7\nHPghsMq/938fvkaMTF8BvgecB5wJfBv4DvAW//P+HfCphHv1FeB34z5f9LrhzwGcA4wBP+e/7vY/\n/1r/9RW+TAL8W3/f8P04Cfx3YIF/36P36Fr//nQB64FjwGsTnsMngc/79/MM4AHgzpTfVSDH7UAJ\n+E1fzkX+9r8D/tE/33L/vv5OSIaTwH/273M59N4H/PtxB3AI+HP/c/468BJweuj6mb/X7fxnK4TG\nM+TPhn4G/BqwJbTtH4CfD8183wt8UVUPh0+gqiPAQ8BHki4knpnmGuAWVR1T1W8zfbXxVuBpVf2U\nqp5U1W8B24F3hPa5X1W/oaon8RTCmhzH/qOqPqKqk6r6iqp+RVUf918/BtyLN8hk4Qo8Zfn3/vXu\nBQ4AV4b2+ZSqfkdVx4FtIVnzUgVeLyI/p6pH/c/m4s9U9XlVPYI3eAXXfKcvz35VHQNuy3DdT6nq\nU6r6IvBF4ClV/bJ/7+8D+mv8PE5U9RngW8Cg/9alwJiq7va37/BlUlX9J7zv3a+GTjGJ9/067t/3\n6Pnv8+/PpKpuBb6LN8GYgYgI8HvABlU9oqovAf8Nz2yaRhW4XVWrqvog8DKwwv8NrAduUtWXVPVp\n4BNMX4k/r6r/0/9eBZ/hB/53ewLYCpztn/+4qj6EN1F7nf8Z6/letxWmEBrPoKr24s00Pgj8k4j8\nAoA/cNwHvNf/cVzH9AE8zB8B/yE41sESvFnPs6H3wv8/B3ijvywe9RXVdUD4nD8M/X8MOD3HseFr\nISJvFJFd/tL9ReD3gVizTgxLgWci7z0D9GWQNS/X4M0ynxGRfxKRNyXs67rmUtz33cWPQv8fj3ld\n6+dJ4zN4Ky6Ad/uvARCR3xCR3eKZBUfx7kv4mR1Wz/wZi4i8V0RGQt+R83E/8yV4q6m9of3/n/9+\nGj/1FWdA8CxeDcxn+ncn+r2JezbRe4+qxj6POr/XbYUphCahqhOqej8wAfyb0Ka/xZtd/hreEvcL\njuMPAPcDH024zGG8pe9ZoffODv3/WeCfVLU39He6qv6HDB8hy7HRUrmfwTMHnK2qZ+LZkMWxb5Tn\n8ZRQmGVAJYOsUY7hDTwARJWqqn5TVd8GvAYYwltt5OUF3Pe9Xo75//aE3kuaGKRxH/BmETkLeDu+\nQhCRBXirvv8B/Lw/kXmQU88MEp6biJwD/CXexOdV/vFP4H7mP8EbaFeFvlNnqheEUSs/wVs9hL87\n0e9NvSWdk77XHYUphCYhHm/Ds1E/Gdr0z3j+gbvx7PUnEk5zG56dszduo7/cvR/PUdYjIivxzFAB\nXwB+SUTeIyIl/+9f+w6/NGo59gzgiKq+IiJvwJuNBhzGMz/8K8exD/rXe7eIzBOR9Xg+kViFmcI+\nYJWIrBHPCXxrsEFE5ouX/3CmqlbxTHsTNVxjG/ABEfllEenBW9E1BN+EWAGuF5FuEfltPDu/ix/h\nvq/B+b4CfArPVBJ8H+fjrWQPAydF5Dfw7OdZWYg32B4Gz0GNt0IIy3WWiMz35ZjEUyB3ichr/GP6\nROSyHNechv8b2AZ8XETO8JXUHwCfrvWcMSR9rzsKUwiN5wEReRlvoPk48D5V3R9sVFXFc4Kd4//r\nRFV/gOesXJiw2wfxnJQ/9Pe9FzjuH/8S3g/8t/Bm4D/klIMwkRqP/Y/A7SLyEt4AOTXz9s1lHwce\n8c0FayPX+yme3+JG4KfAHwJvVdWfpMkaI/t38ByQX8azaUdj4d8DPC0iP8Nb/l9fwzW+CPwZsAvP\nWfx1f9PxvOdy8HvARrx7sQr4l4R9bwX+1r+v73Ts8xk8B/aUuch/xv8F7zkdxRvoPp9VQN9n9Qm8\nz/4jPMfrI6FdHgb2Az8UkeA5fgTvfu327/+X8YIr6uE/462qvo/3rD8D/HWd5wzj/F53GuKNT0an\nICL/HfgFVX1f6s5Gw/BXTk8ACyK2bsNoG2yF0OaIlytwgW+iegPwO3jRTEaTEZG3+yaoRXirpwdM\nGRjtjCmE9ucMPD/CMbyl7CfwYrKN5vPv8eznT+H5IbI4640YROSj4iWwRf++WLRscwkzGRmGYRiA\nrRAMwzAMn7Yq9PTqV79aly9fXrQYhmEYbcXevXt/oqqpCYBtpRCWL1/Onj17ihbDMAyjrRCRaBWA\nWMxkZBiGYQCmEAzDMAwfUwiGYRgGYArBMAzD8DGFYBiGYQCmEAzDMAyfwsJO/bLEX8WrnjkP+Jyq\n3lKUPIZhdC5DwxW27DzI86PjLO0ts/GyFQz296UfOMcoMg/hOHCpqr4sIiXgayLyxaC1n2EYRiMY\nGq5w0/2PM1712l5URse56f7HAUwpRCjMZOT3cH3Zf1ny/6ywkmEYDWXLzoNTyiBgvDrBlp0HC5Ko\ndSnUh+B3gxoBfgx8SVUfjdnnBhHZIyJ7Dh8+PPMkhmEYCTw/Op7r/blMoQrB7zu8Bq837RtE5PyY\nfe5W1QFVHViyJEsvbsMwjFMs7S3nen8u0xJRRqo6itfz9fKCRTEMo8UZGq5wyeaHOXfTDi7Z/DBD\nw5XE/TdetoJyqXvae+VSNxsvq7dzZ+dRmEIQkSUi0uv/v4zX7/VAUfIYhtH6BA7iyug4yikHcZJS\nGOzv486rV9PXW0aAvt4yd1692hzKMRQZZfRavMbg3XiKaZuqfqFAeQzDaHGSHMRJA/xgf58pgAwU\nphBU9TGgv6jrG4bRfpiDuLm0hA/BMAwjC+Ygbi6mEAzDaBvMQdxc2qpjmmEYc5vAD9BqZSg6pTSG\nKQTDMNqKVnMQd1JpDDMZGYZh1EEnlcYwhWAYhlEHnRT5ZCYjwzAS6RT7eLNY2lumEjP4t2Pkk60Q\nDMNwUktm8FyjkyKfbIVgGIaTWjODw2RdYbTrSqRVI59qwRSCYRhO6rWPZ43AaedInXZVZHGYycgw\nDCf1ZgZnjcBx7fehrSOZKpoWRaeZ1EwhGIbhpF77eNYVRtKKYzYG2bwltQM6KeQUTCEYhpFAvaWj\ns64w0lYczRxk65nld1LIKZgPwTCMFOrJDN542YppvgGIX2HE7RelWYNsPY7zTgo5BVshGIbRRLKu\nMML7uWjWIFvPLL+TQk4BRFWLliEzAwMDumfPnqLFMAyjiUQjjsAbZMOKpFGRPUPDFW7cto+JmHFw\nUU+JnvnzOiJcVkT2qupA6n6mEAzDKJrooLpu5RJ2HTgcO8hmURhZr+kyU5W6BRSqk6fGx1qu0Spk\nVQjmQzAMo1DichC27604B99GJMu5zgPQLcLC+fMYHa/WfY12w3wIhjGHqTXcspHkDd1sVGSPa/9J\nVV6MKIOAyuh42+YYZMFWCIYxR2lkdnBg8qmMjtMtwoQqfRnt6XkH+EZF9qSdJ24bUFcGdav7G2yF\nYBhzlEYlVYXj+IEpB21ldJyN9+1LnVHnzYZuVGRP0nnitgXUmhPRDlnNhSkEETlbRHaJyJMisl9E\n/mtRshjGXKRRpheXLR48p+ytn98/7b2omWrdyiW5Bvh6k+Wi5+ktl6beO63UNW2bi1pyItohq7lI\nk9FJ4EZV/ZaInAHsFZEvqeq3C5TJMOYM9ZhewqaPtDjFsHPW5UC+5qI+Z1RRHI1so3n85OTU/4+O\nVaeZhAIzWJRaciLaIau5MIWgqi8AL/j/f0lEngT6AFMIhtFkhoYrjJ04OeP9LKaXpHDNNFyz5F0H\nDvPIpktzn69e0iKWsmZaZ6EdsppbwqksIsuBfuDRmG03ADcALFu2bFblMox6aFUHomtA7y2XuPWq\nVakyJpmI4ljUc8okkzRLbuT9GhqucOvn90+tThb1lLjlypmfLW3WnqfXQZr8jVQuzaJwhSAipwPb\ngQ+p6s+i21X1buBu8BLTZlk8w6iJVq7v7xrQFy6Yl0m2PCaOUrdwy5Wrpl6fVupivDo5Yz+BzH0T\n0gbnoeEKG+/bNy2p7OhYlY2f2zfjfFlm7VnMU1medzs00ik0ykhESnjK4B5Vvb9IWQyjkbSSAzHq\nxHWFU2Yd6F0mjr7eMp9cv2aas3fLOy6cNuCF7fVhJiH1fmWN0tmy8+A0ZRBQndAZ/RUaFbGU9XkP\n9vfxyKZL+cHmK3hk06UtpQygwBWCiAjwV8CTqvonRclhGM2gVRyIcTNXgVhHcDDQ12P6SJtNx4zT\niYTvV1ITnS07D05dP+0eV0bH2bB1hD3PHOGOwdVT565n1t6s5z3bZsciTUaXAO8BHheREf+9j6rq\ngwXKZBgNoVUciHGDqMIMpRAM6M02fQRJa1kJ368sTXSCY1yroAAF7tl9iIFzFjckYqkZz7sIs2Nh\nJiNV/ZqqiqpeoKpr/D9TBkZH0CplkV2DqEJsHH+zTB+B2SqPMojer6xNdDZetoJSl6SeX6FhJrxm\nPO8izI6FO5UNoxNpFQeia+ba11uODfOs178QR54w1b7esvN+rVu5hE/vPpR4/POj41PHhKOMkvZv\nBLU+7ySTUBFmR1MIhlEHST/oRiZP1cq6lUu4Z/ehWPNQlKHhSqp/Ie6YtEEwa5iqS0kF7HjshdRz\nBHKG7/3QcIUNW0diP1eXCEPDFedzymPDdz1v1znSTEJFmB2tlpFh1Eir16YZGq6wfW9l2kAowDUX\nxQ9cW3YejB00BZwKJMvnzzKjzWJeOTqWPNsHTwFGGezv47q1y4gzIk2oOp9ZI55v0jnSTEJFmB1N\nIRhGjbRSaGkcLofyrgOHY/d3mYuUeCem6/PfuG3ftHLaaTPabhFnLaJwyGwWtu+txA7Ydwyu5q71\na+iWmWrB9czSnm+W0uG3PbDfeY4sSXGNqNmUBzMZGUaNtEpoqYs88iWZi4I+x1HTh0uBhKudfmjr\nCD2l5HnnpGrsIHfz0OMzzF1pJDWxGezvY8PWkZij4u9JWlZ1WgTQ0HDFuapJuod5k+IaiSkEw6iR\nokNL0+zbeeRLMxflyWeIMhaTmZwmz81Dj6c6kF0kKeQ89yRp36TV0YatIyztLXPs+MxaUeFztGIp\nCzMZGUaNFBlamsW+vfGyFV5v4BClbomVL8lctGHrCDdu2+fMZ6iHqH9iaLjCmtseqlkZQLJCjntm\nQrzvIen5upTOhOrU80iKcAqU92ybhNKwFYJh1EiRoaVZsnb3PHOE6kRkDh8zpU8yFwWHuPIHgnyG\n50fH6cqZdAZw3dpl00wsWcJTu0WYVKW3p8TLr5ycVqYiqpDjVlHXXNQ3zRSleL6HgXMWA9Ofp6ss\nt6ssdhZ6y6WWikQLI5rzARbJwMCA7tmzp2gxDKNwlqc4WUtdElvPB2aGdybVN0pDgLvWr4kNo0yj\n1CVsufZUraN65ADoElgwr4tXqpMs7S2zbuUStu+tzDDJLJjXFTt7X9RT4pXq5Iz942bttZYAd52v\n2YjIXlUdSNvPVgiGURBZK3dG94GZpSeiuJQBzLSx1+MED7J9wzPdLAlhgYxhB3C9zvhJZaqSamV0\nPNYhPV6dcA7icQ5gl5M6ujp0rY4W9ZTomT+vZaubRjGFYBgFkDVKJW6fBfO6ckXeRDmzXJqajXeL\n1HUumD6QB4oha4RQ+NikyKWsDuwwjbJ9uBRVWCnEOdnLpe7YHgytjCkEw6iBPBmscfumdeoCt5+g\nlk5lYUbHq1Mz+Lw2/zi6RDh3045pK5h7H30204AcdgDHRd0AiMB1b1zGrgOH6zIpBSzqKTE6Vs2s\nMJKytMPyhosGLuopoeo55MOVWFsdUwiGkZOss/u4mWOwr2tQD89G6zGhdOH1GJgNwnkHGz+3DzSb\nook6gIN7d9sD+6eZb1SZ6rsc9QnkJZi1f8iRjxBHXAQSJFSSFXj5+Mkph34rNUdKw8JODSMnWTJY\ng5BQmGm6GK9OxGbMwvTZaK35DH29Zc4Mta2cTaoTmui/6BZJDLEc7O+jZ/7MeWrQd/nOq1eToZCp\nk+CafTnurSuz21lJVpkR3dVKGexJmEIwjAhpJQnSMoCzFHObUE3NYYiLg08jOMdohro/s43gfe60\nvsRJFVcH+/sSm+z0lt2KsK+3PBUNNXbCnTQWd9048irsVslgT8IUgmGEyJLw5RoIgvez/PAXzu9O\nTUoKEpfyEMxEz0wYGIsibDbbsHWEm4cen9o2NFyh//aHEk05S3vLUzkTcfT1lhm55df55Po1TmUb\nPN8shfLC140jr8Ke7eZItWA+BMMIkcXZm1ZyIEvHrmMnvGOTyj0DNSVBPT86Tm8DTUaXnLeYbx16\nsW5ndhgFPr37EDsee4GjY9XUKKIgwzpLRdZoSOiZ5RIinoM3b/JcUuZ5cJ0bt+3LdM7K6DiXbH64\npR3MtkIwjBBZCsKllRzIOnPMalOuZSbaKJNRtwjXDizjmov6nH6Peghm6mnD6UTIQRtHcHxg6gsi\ne+5av4bjJyc56kcV5Y2qSksiG+zv4xPvvDDz82m1EulRLFPZMEK4smVdzVtc4adZ4vAF+MHmKzKd\nE7LNRINM2HpKK0QpdQtocrLbbNBbLvHSKydj70GgqqJ5AKeVunKZh8KkNewJE44qa/S5G0HWTGVb\nIRhGiDwF65L8DbsOHE6d9cbZlF3nBK9MdBLdIlPNb7L2Fc5CUuRQ49cMbkbHq4k1leKiuWpVBnmL\nFAY9prPej1Z1MJsPwTBC5ClYl+RvSPvBu7qQJZ0zzTcxoTqtSNtsjNazvWaoJWM57/nrKTGRxX8U\n7NeKFKoQROSvgbcCP1bV84uUxTACokohsPVHB4gkf0PawODqQpZ0zrvWr0ktqBaOd59R6bQDCGcD\nZyVrkl4jzDiubOswRfc8SKJok9HfAJcXLINhTCPObBMNk4Tk8NM0R3BcYtTQcIWuhIS1qDPbRWV0\nvGH+g1ZEIVdyWhZl0KhBOi7g4Pq1y1qq50ESha4QVPWrIrK8SBkMI4qrJME9uw8xcM7iqR/zupVL\nYhu5rFu5hMH+PvY8cyR2e6lrZpOaJCd0eLCKFlSrh6ASZ7spj95yKVM11ayEfS+NoNV6HOSh5X0I\nInIDcAPAsmXLCpbG6FTCkT1JjWLC+QiukgbB+67tUcfo0HAlMSIpWvagljr8cVxxwWu5Y3A1Q8MV\nNmwdmXV/QK00Ovo17Htp14G8URRtMkpFVe9W1QFVHViyJL7IlGHUQ9RElESW4nPB+67tk8q0WHRX\nslWYINrotgf2NyxBbMdjLwDejHY2lEGDgp4YHasmlqiohWbUGkorgdKKtLxCMIxmk6X2UECW4nPB\n+0mRJOEBKGsIYj1hlHEcHavSf/tDqd3XGoEAb/pXixtyrqW9ZW69alXDB69Gms6ylEBpRUwhGHOe\nrANyluJz4X3SHMuV0XHO3bTD6UieDRqpYJJQ4JGnjtR9HsG7b1t2HuTda5c5VwrlUjfXr12WK8O7\nkZnYaRVxW5VCFYKI3At8HVghIs+JyO8UKY8xN8kSE95bLjmLz8VFkAQ+ibSVRy3lFOYSpS5hkV+X\nKdpXYvveCrdetYqnN1/BJ9evmfEc7hhczTUX9WVOx5hQbZhpJ0sJlFak6CijdxV5fcOAbLHjCxfM\ni3U4xkWU1NqAvVuESb889LqVS9j6jWcLLxdRBN1+Abq+UIJYXEkRV7/jMF/Y90Iu/0ijmtm48lBa\nNSEtoOWjjAyj2WRpDu+a2Q0NV6YdF8xma3H8Tqpy1/o1bNl5MDZcda7wiXdemCsJ0NXBbs8zR2oK\nT82iaNJIq4jbqpgPwTB8XnrF3TRFYYY5YWi4wsb79k0bdI6OVWu2yyteieZ2ywtoNHF29iQHvste\nf++jz9YsQ73PIK0ibqtiKwRjzhPMMNNs+VFzwpadBxtu0mlXA1EjawxV/Fl/ePBct3LJjFyNYMa9\nwdFUpx7fjMAMGZJwVb1tdQUQxRSCMefJE3YaNifkcRA2uyjbbNPdJUyElGGjP1tY8Q4NV9i+tzLt\nGgJT2cVJpr40FvWUGPV7JYSJJiGGiQ7+61YuYfveygyTVSB/O2EmI2POkzfyI0u4aG+5NM1ccF3O\nEMhWZlFPiTMWNHcuGQ7RdJUS+cK+F+i//aGalUG51M0tV65yKrO470VcfsE9uw+1ZYhpHLZCMOY8\nWUsWh0kKFy11C7detSp2dpjWNKfVuX7tMu4YXM25s5DMFjiMXc8mjyIQoLenhCq8OF6d1njItXqL\n81u4lJNL/nbDFIIx53FFhFxzUd80U0AWFs7v5uNvn+k8jDN7tBuLekrcMbgaqE2JdufsZ9zbU5oy\nvdRDtwhP3fmbsdsu2fxwao/mMHkG+VYPMY3DFIIx5wkG79se2D8VIbRgXhcD5yxm4JzFqUXvwvT2\nzJ86X9jWnLe5eysS7tOcJXcjyoQq5VJ3pmNK3RJr268F131PWn24+lW4FGF0ldEOIaZxmEIw2g5X\nREc95+ntKfFyKOx0dLzKTfc/zp1Xr55qmuLqtxwmmEFGY+PbXRmANxhG79mCeV2ZTTdBmekdj70Q\nG5orAqreSuTlV042bDXl6j2RtPqIOwaSV5O7Dhyu+ztZNKYQjLbClYQE+SI6hoYrbPzcvqmuYnED\n1Hh1gtse2D81AJ5ZLlHqlsROZIGZIE/kUjsgwPJXlafd+6NjVcql7sz9CYIy03de7ZmdXEr9ks0P\nN6zGkmumnvR8kmb3eVqstiOibTRzGRgY0D179hQthlEDeWf1rv1ds/S87Q/7b3+opkGn1CWcfto8\njo5VZ5gJgtd9NdjX5xJpz+rcTTtyrQ7CfZAh22CddI1Prl/TMQN8gIjsVdWBtP1shWA0nbyz+qT9\nG1U0rNYZaHVSOTpWpc+PP9914PDU4B8uvGa4Sbs/eRzWcQ7jLIO56xp9fqvSuYrlIRhNJ28p4KT9\n03oQzBaV0XE+vfsQz79og38tLE9oGrPxshW5KpTWQlxpcsHLiJ7LmEIwmk7eWX3S+2k9CLJSLjXm\nq99GFteWozI6zoatI9w8NN25O9jfx3Vrl2VSCi7nbxqD/X0zSmMrsH1vpeWb2DQTMxkZTSdvKeCk\n/fM69eJ8EQAn52BZ6VZE8ZL1ov2M7xhcPRXym2Q+qmdGv+vA4Rl+hEZUOm1nbIVgNJ28s/q0/Qf7\n+9h42QqW9pZ53u+eFTeriyszsGHrCB/aOpIYKdQoekpdHVOuopkEdYOiDPb38cimSxNXAbsOHAZq\n61/crk1smokpBKPp5C0FnLZ/1n61ecoMNIP/dvUFU5/DSCZpEE4yB4b7IeTtX9wq/qhWwsJOjbYj\na+jpbDSPTyIcvpglqW0ukxaKuua2+CJ2gbKtJRQ5rrNdudQ9Y/LRCTkHWcNObYVgtB1pS/2h4Qr9\ntz80myLFcuO2fVMmjOWvip91XnLe4lmWqjj6estcct7iGc7iLEEBt161ymlGdH0fgqq0LhNSo1ai\nnYStEIy2I2mFUEuNnTBZa+3kpdP6IdSKABeft5infzqee9adN1kxTHTmn4VGJUG2ApaYZnQsGy9b\nMa3sBHjF0DZetqKukhG95RK3XrUqNbKlFkwZeCjwL08d4a6EbOCk7mNxx2SZBNQSPTQXnc6FKgQR\nuRz4U6Ab+L+qurlIeYw2Iq7FFek/VlcJ5lLX9B4G9awyjGTSupG5stQhPtw4GorcqP4EecOlO4HC\nTEYi0g18B/g14Dngm8C7VPXbrmPMZDQ3SHPkJS3lwV0aIWw2SLqGOYBnh6AGUVACJKlM+KKeEq9U\nJxMdwAGu59ctwifeeWHmlUm0NWbSNVuddjAZvQH4nqp+H0BEPgu8DXAqBKOzGRquTOtJAPF1j5KW\n8netX+Oc3Z8Wyk5OaoDeySZJ0HXaAAAcJ0lEQVSB2SIoCZ3UIS5w1H5696Gp91ylKFzVaONWGi4T\n0oSqs4ZW3Mpk+95Kx5S1zkqRCqEPeDb0+jngjdGdROQG4AaAZcuWzY5kRlMIz8DOLJcQ8ZquuGZj\nAdEfvmspH5gigh9xZXR8mjP36Fg1cUDopGY2RROeRYcH/EYTRBLFmZBu3LZvxnN0KRFX/axdBw63\nnQO5HooMO40rVTLjV6iqd6vqgKoOLFkytwtPtTPREL7R8SpH/Y5YrkblYcKz9rhM5oBgxnns+El6\nyyVnaYIk2YpSBqVu4fq17T/pCVcMDVpu1kPQc8FFXEjoYH8fk47nGLcCzONAriUrul0oUiE8B5wd\nen0W8HxBshhNJi36J20IDjvywvHjLkbHq86mLRU/uzWrbLNFdUKbOpueDYJorzC1ZGp3i0zLDYjL\nQ4gSVfZ5MpGz7tvpuQlFmoy+CfyiiJwLVIDfAt5doDxGE6nHLl/qEsZOnIw1DeRtphIQNh2Zz6CB\nKOx55sg052zP/PzzTpfz99bP70/szhZdScZlIsclwWXdN6k0eyf4FgpbIajqSeCDwE7gSWCbqu4v\nSh6judQaqlfqAoRp5qXwjKzW84Znk/WGEfaWS1bEzqc66a1ywjPo7/74WK5zLOopxQ6ug/19LFyQ\nPId1rSTTamhl3bfTcxMS766IfAh4BBj2B/CGoqoPAg82+rxG61FrBvGkChOTbsdgPZnJwY+4nnMI\nXlkFYMasOO9AaHiz8luuXOXcnjTwxs3ok6LJomTZt9NzE9JMRmfhJY6tFJHHgH/BUxBfV9UjzRbO\n6ByiyUPhKKPenpKzpaXLwRsMDMF5o+Gq4JmaEJylroMfcXCOP9g2Qt42CRo6PporYeSjWyQ1xt81\nIGc5thHkMUO1I4kKQVU/DCAi84EB4GLgt4G/FJFRVX1980U0OgXXDOySzQ87FYIrszhqGnAlm0G8\nsgh+xMExtSaiuaJfLLEtP5OqqQO6a0CerWSxvA2a2o2sTuUy8HPAmf7f88DjiUcYRkaSzADveuPZ\nsdmicTMyl8IJK4vK6DjdIoxXJ7j18/s5duJkXc1yTpyciM1wtWJ2pwgK2u3+/tHEkN4sZpdWGJDz\nmKHajcTSFSJyN7AKeAl4FNgN7FbVo7Mj3nSsdEVrUm/NeFepgXKpi8ULF0wN4hOqUxVN85w/LgO6\nkZS6ZZpSaQVl0CXkNn818zzh55alD0FAp/QjKJpGla5YBiwAvosXGvocMFq/eEankFSMLOsPN84M\nUOoSTk7qlKKYUJ1aGSQ1LwGmrQRmI8ksusIoWhlA8iCeVWF9cv2a1DDPJMLXiftepA30jfhuGflI\nLW4nIoK3SrjY/zsfOILnWL6l6RKGsBVC69GomvHRwX3sxMnYGX1w3rhZZpoT2fAG6Z753Rw7kRxR\ndf3aZdwxuDqx61yf7+CNKphyqZvTSl2Jzy8LndSPoGgaVtxOPY3xhIiMAi/6f2/FK043qwrBaD0a\nFZcdtcue6xiIKqPjnHfTg7Ez/2ojbBsdjkKiMoiadlyrifCgHLdS27B1JPb8eb4XnR7z34okJqaJ\nyH8Rkc+KyLPAV/EUwUHgamDu9P4znDSyUXm4RozEVbrymQ0zUFA6YS4RDPJhk47rTh87fnKqlg/A\nI5su5a71awDYsHWELscDVMhc/6eR3y0jG2krhOXA54ANqvpC88Ux2o11K5fMKHFcS1x21ARUdLHR\nCVV6y6Wa7eftRtwzS5qJB/clsOvveebItGiwJKWd1RfQ6TH/rYj1VDZqYmi4EutwFOA63/6chyKa\n0ixKSIibC/T1lnl+dJzenhKq8OL4qVLkQfnwrLgc+GmO/bSoMYsyagzt0CDHaFPiHLoBCuw6cDj3\nOYtI5Do6VqXUBdXJ5l2jW4RJ1ZaIPArTWy5NOec33rdvyv8SbViTFdeg7ypBHZC2WujkmP9WpMjy\n10abklYuOlpeOq1+fOC8nG2E5ioD8AbKIpXBgnldXvRViKB/NHjVQ7M64/t6y87M7G6Hz2Bpb9m5\nLSCuR4VRDKYQjNxkifIIKpJmqR+f5LxsJq02a28Gx09OsuXaC6dV8dxy7anS0ll9JILnOI7rS1Au\ndfOuN54d+/66lUsyBQFY5FBrYCYjIzeuAmNhwrM+V/14IFMdoVbI/C2S3nKJVUvP4JGnaq8nWW/c\nfpcIQ8OVxKSygXMWxzapz4JFDrUG5lQ2cpPkQwgTGApc37BSlySaKwL7+2mlLo6fnGRSvVIKKIQt\nPZ2uMJ7efEXdTvdFPSVuuXLVDHt8/+0PZXasR8tLpDl8s8o8m8Xp5ipZncpmMjJyE20mkmQ/Tpr5\npdmuA/v7eHVyqhTDpIJ0ncoR6BbhurXLeHrzFYl9d9uVoP1kvSaVo2PV2FaPt1y5ilL39OfX3RWf\ngxFe2WUxBSbJvKinlNq0xph9zGRk1EQ4+sNVrCyIF9/4uX0NLScRbpgzocr2vRUGzlnMrVetmhYx\n0+pEi+JFCd/DLGa6NKKtHoMZfnVCZxQPTMs0ztJK0iVzvWVNLPS0edgKwaibpPaDg/19LJzf3HlH\neOY6f177fKVPXzCPhfPdrTeDz3Xz0OOMnWhMw8JgQA/P8GFm8cDenvjV1tKUFUu0p3GcozlPYlmn\nN7VvNWyFYMyglhlZUrz4i7OQ7VsZHW/4SiQPQZnoqD9DgNe9ZiHf+/GxGX4OLw9CElcKteYFuAgG\ndNcM/7YHvLbmL78yUwGVuiV1xRJtXBRcq9bZfac3tW81TCEYwCklEK1e2YiSw40wd6TRVUeV00V+\npm6tZSoC53dcVq4C3z885nR6Vyd1Wt2mZjrIw7Nz1wz/6FjVmZuwcP68qe9A1rIS9SaWWYG72aWQ\n9bWIXCsi+0VkUkRSPd9Gc4maD6JDQb2JQ7NRe6Yet8HRsSovjlfpLZcSTTguAue3K94+LQ4/vLlR\nyqDUJVy/dtkMMx540T9J13EpxvBKL2om7C2XOK3UxYatI5mL12XBCtzNLkWtEJ7Aq5j6fwq6vhEi\nLfMY6puRDfb3NbVjWSNQal8htBoC05LPArKGC7uIDsLhXtbNamRjBe5ml0JWCKr6pKparnqLkGWw\nr3dGdsuVMzNc51p56dmgXOrmrvVrYgfiLIofPBNaHmdwkp2/XpICFozGYz4EI9XG34gZWdjBGNdl\nK0ynJ5rVQtCdLK3HcdJgmUXxl0vd3HKlV+coqzPYdd7K6DiXbH647nBRK3A3ezRNIYjIl4FfiNn0\nMVX9xxznuQG4AWDZsmUNks4IE9fTIBiU+/wSBFt2HmTD1pGG/LCTMliD64Vr6891wnH7adm/waw8\n7vmkKf5oKeqsz9h1XuFUFVvrh9weNM1kpKpvUdXzY/4yKwP/PHer6oCqDixZsqRZ4s5ZhoYrbN9b\nmaEMguzfjZetYPveSu448KQKp64ZZVBA7Y7B1TMcluHM1uvXLkutoNkphEM9IX2Wn/R84vICAsI5\nCHmJO2/cKs+qmrY+ZjKa48TZf8M9DWqJA09zMmaNYU8anAbOWcyHHNm0WekWYcE8YazZNbDrIBzq\nCfkKC0bvX/D6xm37ZkQ+1RPbH5dv4JLRwkVbm6LCTt8uIs8BbwJ2iMjOIuQw0uO8a4kDT3MyumaU\nldFxzrvpQZY7+iaEGezvY5EjmzYL5VI3n3jnhYw3SRlE6wPVSjSpL2mWH8b1fAb7+5xNa2oZrIOV\nYFDq4q71a3hk06VTNZiiWLhoa1NUlNE/qOpZqrpAVX9eVS8rQg4jPc67ljhw1+wweD8cOQLTzQvB\nzDWLaeqWK1fNaP4Sh0BsTH6wWmkGC+fPcw6KYdJMYHGhnndevTrVZJb0uRoV259UVqIRZSuM2ad9\nCr8YmUjrThYlabZ+yeaHWbdySe4ftmuwCr8/2N83NZN0Bc2k2ZwH+/vYcu2FiVVOwz2eH9l0KT/Y\nfAWPbLp0Kn7+2PHG1AiK8uJ4lY2XrUhcKQRO3O17K7HJa677PNjfxyfeeWGqP8BFowbrNHOihYu2\nH+ZD6CBqSRBKCgetjI6zfW+Fay7qY9eBw87wwWjtozwZu2lmirTtUV9DuARHUEpi14HD05q7BPvV\nk6SVxtLe8tT1Pnr/YzP8FEE3sTh7PnjKM2kAjT63aLXSpIG3ETWGIN2caOGi7YcphA6i1kJgSeGg\n49UJdh047CxXHNek3UWczT/NSRo1Y6QV3gv+H1WMG+/bx20P7Gd0rMrS3jLHjp9smjIQTpXrCGfz\nxnUTS2pOn6WgIJwa2LMog/CxcfvlKWyYJTjAaC9MIXQQ9RYCq+X4PE3a48a+oPZ+3BkEL0ciSG46\nrdQ1zQnsWgHFKcbqpE6Vzmh2ob2e+d1s2DrClp0HpwbU6AB8yeaHExVSlkG10SUj8p7Pykp0HuZD\n6CDqdRbWcnye+j9xZbAH+/u4bu2y2DIW8+d1sfUbz045LeMiguL8DEWFNgpedNGxExOpORuNyAxv\ndMmIvOczP0HnYSuENiDrMr7eGVtcxnIjZ3xximVouMKuA4dRZiYzHT+ZLSQ0qgAaWW67XOrizqsv\nyNRrQZlZgjvOZJfm6M86qDa6NHQt5zM/QWdhCqHFybOMr8dZ6MpYvuai5B/8op5SpiqmcYol+tlq\nrV8UVTRxiq1WxquTU59/w7aRWLNXGpXRcc7dtGOa7yCJektG1GrDz3s+a23ZeZjJqMWpZRkfDa+s\n9TrhjGUXcU3ao7giZrJW30wi7MAFuHno8YYpg4Bg4EtSBkK80zwgMCHds/tQ4mfOkrsQ0OhY/zzn\ns9aWnYkphBZntjpG1Xqdwf4+trzjwmlJZmGCjOA4xdSIz3DxeYunNY2vVRm48ttKXUxrHpQkR1yJ\n7yhJsuUdzBttw89zvmaWvDaKw0xGLc5shfbVc52wHTnOjADElkHOYutP6jcM8PRPTx2/ZefBmpRB\nudTNNRf18ZlHD80oLV2dhOpk+irmW4de5NoBz/4ffP48sqTlHbhotA0/6/mstWVnYiuEFme2SgA0\n6jpRkxWQq7xBqUumVTbd8o4LE69XGR1n+aYdLN+0IzUHwlU99c6rV3PH4OpMNYJcjFcnuHHbvmk1\nfVzmnzyrqFbFWlt2JrZCaHEalVVa1HWSTAuBwki7ZpCNWw9XXPBa7hhcnbjPsRP1+TOidZiuuahv\nRl+HYDWSlPndDlgOQmciWkvYREEMDAzonj17ihbDyMG5m3Y4k85+sPmKTOeIZkPXQrnUnWqSWb5p\nR6Zz5enotqinhKqXg9Gug78LizJqH0Rkr6oOpO1nKwSjqWT1TSQNLsG/t35+f65EuDBZSnj0lkup\n58/b0e3oWDWxz3E7YzkInYetEIymEldErtQlnH7avKm6QnEDbLhKaRjXiiMLaauStJVI+Pih4Uou\nBRVug2kYs42tEIyWIOqbOLNc4tiJk9PqCn1696EZxylwz+5DDJyzeNosNC0yKVA2cclyXSLcPPT4\nlP2+N8acs+XaC50VSKOrmqyZ1GDRN0Z7YFFGRtMJRx4tXDAvtQREgMKMuPaNl61wNsXpLZfYcu2F\nznyACVU+vfvQVMTT0bEqo+PVadFPQGyvgajDNG9SnUXfGO2ArRCMukiy/cdtyztTjt0/og9K3cKW\nd8wM23TN9F2MVye47YH99Myfx3h1IrHHQNLnKJe6LfrGaEvmxAohbxcxIxtJ5Qtc23pz9kGOzqy3\n7Dw4Y4VRndAZK4mk3sFJHB2rTpmkJlSnBvOosnHN+IO8BqsAarQjHb9CaHTN+GbTTqF8aeUL4rYt\nmNc1YwadxLHjJxkarrDnmSPc++izzhl/eMY+VXsoz4dx4IpOSorDryX6pp2eu9G5dLxCqLWLWBEU\npbxqHYxqKV8QROUE5phFPSVeHK/OKBkR3v8Pto6Q5r4NZuzNaI0Z93kamcjXbpMWo3MpRCGIyBbg\nSuAE8BTwAVUdbca12qnmShHKK20wSlIWaTkGSdFAgTnmigtey9ZvPJto3klTBmEbfZKzN8ghiItq\nSsJlHmpUHH47TVqMzqYoH8KXgPNV9QLgO8BNzbpQO9VcKUJ5JQ1GcX6ADVtHuHnIUxhJ9Y/itkUZ\nr05w76PP1pWBHLXRu+6VAI9supQ7Blc7awwt6ik1rG5UHr9VO01ajM6mEIWgqg+p6kn/5W7grGZd\na7aKwzWCIpRX0mDk6pFwz+5DDA1XEsslR7e5yBMFFEfUTJPlHsaFrpa6hFuuXBX7eYBcQQl5ewW0\n06TF6GxaIcrot4EvNuvk7dT3tQjllTQYuZRFOD8gqSFPeJtrVt4tyc110ojLU8h0D6OX9V/nqdaa\nJFOeXgHtNGkxOpumKQQR+bKIPBHz97bQPh8DTgL3JJznBhHZIyJ7Dh9O7t7lotYuYrNNEcoraTBK\nmqHmNWe4rvOuN55dV9npQI7ARLNh6winlbroLZec9zBr6Gqwb95GMHlNQO00aTE6m6Y5lVX1LUnb\nReR9wFuBf6cJBZVU9W7gbvBqGTVUyBZktguGpUXLbNg6Ehu+mdeckXSdgXMW11ziemlveYZjPK2g\nXNYBe2i44pQpSSHW0mzICsUZrUBRUUaXAx8B/q2qjhUhg3EK12A02N/HnmeOzGhLWas5I+k6g/19\nuQvXBXLkjdLJMmAHvZldJA3u1ivAaFeK8iH8L+AM4EsiMiIif1GQHEYKdwyunur+1WxzhmuQ7RaZ\namIfZwrKa6JJs9mn9WZOG9zNBGS0K4WsEFT1dUVc16iNeswZeZLeXDPrtME0r4kmzUyWluWcZXA3\nE5DRjnR8prJRHLUkvYWb1GfN/nUlm61buWSaLNHzuvoTpDnMwxFWhtFJmEIwmkZahE6csrjz6tW5\nG8nsOhAffRa8n7c0RFrPBSstYXQqrZCHYHQoeZPexqsT3LhtX+6qtGk+hEbkBURJCz01jHbEFILR\nNGpJeptQzZwAluU6kK4womUmgExZ1lZawug0TCEYTaPWpLcA1yw8OoCvW7kkMWooSWG4ykwAU8mM\nveX4Hg5WWsLoNEwhGHXjKuSWFH6ZxSwD8cli0QF8+94K11zU5wzzTFJMaeakoeEKx06cJEqpSyyv\nwOg4zKls1EWawzYpGQ1OhX52+f0RosR1TIsbwHcdOOx0RieFmW7YOhJ7TNj/ENcD+vTT5plD2eg4\nTCEYdVFPLf+wsohrbBOXAJbFHxA38LsUU1oOg+t6o2PVxM9mGO2IKYQamcstD8Of3ZXAldfhmrUD\nWdIAXkvnsbQyE7XUJTKMdsUUQg3M5ZaHWVtU1jJgZsnuTRrAa1mtpCkiq0tkzCVMIdTAXG55mNSi\nMqCZA2Y9/oCkc9aqMAyjkzCFUANzueVh0mcUmJUBs1Z/QKOvZxidhimEGpjLdmXXZ+/rLecuOdFo\nzLxjGPVheQg1MJdbHoYLxmV5fzaxstOGUR+2QqiBuWxXTiskVzRm3jGM2jGFUCNzdeCZy/4Tw+h0\nzGRk5CKtkJxhGO2LKQQjF3PZf2IYnY6ZjIxczGX/iWF0OqYQjNzMVf+JYXQ6phBajLlcI8kwjGIp\nRCGIyB8DbwMmgR8D71fV54uQpZWYyzWSDMMonqKcyltU9QJVXQN8AfijguRoKfL2/jUMw2gkhSgE\nVf1Z6OVCcFZRnlNYjL9hGEVSmA9BRD4OvBd4EViXsN8NwA0Ay5Ytmx3hCmIu10gyDKN4mrZCEJEv\ni8gTMX9vA1DVj6nq2cA9wAdd51HVu1V1QFUHliwpvl5OM7EYf8MwiqRpKwRVfUvGXT8D7ABuaZYs\n7YLF+BuGUSRFRRn9oqp+1395FXCgCDlaEYvxNwyjKIryIWwWkRV4YafPAL9fkByGYRiGTyEKQVWv\nKeK6hmEYhhsrbmcYhmEAphAMwzAMH1MIhmEYBmAKwTAMw/AxhWAYhmEAphAMwzAMH1MIhmEYBmAK\nwTAMw/AxhWAYhmEA1kKzJbE2moZhFIEphBbD2mgahlEUZjJqMayNpmEYRWEKocWwNpqGYRSFKYQW\nw9Uu09poGobRbEwhtBjWRtMwjKIwp3KLYW00DcMoClMILYi10TQMowjMZGQYhmEAphAMwzAMH1MI\nhmEYBmAKwTAMw/AxhWAYhmEAIKpatAyZEZHDwDORt18N/KQAcRqByT77tKvcYLIXQbvKDdNlP0dV\nl6Qd0FYKIQ4R2aOqA0XLUQsm++zTrnKDyV4E7So31Ca7mYwMwzAMwBSCYRiG4dMJCuHuogWoA5N9\n9mlXucFkL4J2lRtqkL3tfQiGYRhGY+iEFYJhGIbRAEwhGIZhGECHKQQR+bCIqIi8umhZsiIifywi\nj4nIiIg8JCJLi5YpCyKyRUQO+LL/g4j0Fi1TVkTkWhHZLyKTItLyIYUicrmIHBSR74nIpqLlyYOI\n/LWI/FhEnihaljyIyNkisktEnvS/K/+1aJmyIiKnicg3RGSfL/ttWY/tGIUgImcDvwYcKlqWnGxR\n1QtUdQ3wBeCPihYoI18CzlfVC4DvADcVLE8engCuBr5atCBpiEg38OfAbwCvB94lIq8vVqpc/A1w\nedFC1MBJ4EZV/WVgLfCf2ui+HwcuVdULgTXA5SKyNsuBHaMQgLuAPwTaykuuqj8LvVxIm8ivqg+p\n6kn/5W7grCLlyYOqPqmqB4uWIyNvAL6nqt9X1RPAZ4G3FSxTZlT1q8CRouXIi6q+oKrf8v//EvAk\n0BZNStTjZf9lyf/LNK50hEIQkauAiqruK1qWWhCRj4vIs8B1tM8KIcxvA18sWogOpQ94NvT6Odpk\nYOoURGQ50A88Wqwk2RGRbhEZAX4MfElVM8neNh3TROTLwC/EbPoY8FHg12dXouwkya6q/6iqHwM+\nJiI3AR8EbplVAR2kye3v8zG85fU9sylbGllkbxMk5r22WEV2AiJyOrAd+FBkNd/SqOoEsMb37f2D\niJyvqql+nLZRCKr6lrj3RWQ1cC6wT0TAM118S0TeoKo/nEURnbhkj+EzwA5aRCGkyS0i7wPeCvw7\nbbGElhz3vNV5Djg79Pos4PmCZJlTiEgJTxnco6r3Fy1PLajqqIh8Bc+Pk6oQ2t5kpKqPq+prVHW5\nqi7H+wH9SqsogzRE5BdDL68CDhQlSx5E5HLgI8BVqjpWtDwdzDeBXxSRc0VkPvBbwOcLlqnjEW92\n+VfAk6r6J0XLkwcRWRJE/YlIGXgLGceVtlcIHcBmEXlCRB7DM3u1S3jb/wLOAL7kh8z+RdECZUVE\n3i4izwFvAnaIyM6iZXLhO+4/COzEc2xuU9X9xUqVHRG5F/g6sEJEnhOR3ylapoxcArwHuNT/fo+I\nyG8WLVRGXgvs8seUb+L5EL6Q5UArXWEYhmEAtkIwDMMwfEwhGIZhGIApBMMwDMPHFIJhGIYBmEIw\nDMMwfEwhGEaN+OGrI5G/SRH5jaJlM4xasLBTw2gQInIDXj2qdao6WbQ8hpEXUwiG0QBE5JeAh4GL\nVbXdSrAbBmAmI8OoG7/mzWeAD5syMNoZWyEYRp2IyGbgtar6vqJlMYx6aJtqp4bRiojIm4FrgF8p\nWBTDqBtbIRhGjYjIIuBbwLtV9etFy2MY9WIrBMOond8HXgP8b78XR8Cdqrq1GJEMo3ZshWAYhmEA\nFmVkGIZh+JhCMAzDMABTCIZhGIaPKQTDMAwDMIVgGIZh+JhCMAzDMABTCIZhGIbP/wdyGJtg2zdR\n5gAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "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": [ "covariance matrix for Z and W:\n", "[[0.99653582 0.70792068]\n", " [0.70792068 1.00873687]]\n", "\n", "estimated covariance of Z and W is 0.7079206838367835\n" ] } ], "source": [ "# numpy.cov(x, y) returns a 2 x 2 covariance matrix\n", "# cov(x,x) cov(x,y)\n", "# cov(y,x) cov(y,y)\n", "est_cov = np.cov(Z, W)\n", "print('covariance matrix for Z and W:\\n{}'.format(est_cov))\n", "\n", "print('\\nestimated covariance of Z and W is {}'.format(est_cov[0][1]))" ] }, { "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+CbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJztnX2UHHWZ779P91SSnoB0gnNXMuYF\nWU00hszICEjuuoZlCYrgyFtk0dXdey/X3cNdk3WzDsiShMVDvLNu8Lieu3JXry9EDBCcBYM30UtY\nNWvAxJmAkWQVkYSBXaPJRMh0kp6Z5/5RVZ3q6vpV/eqtq1+ezzl9ku6qrvpV1fTz/H7PKzEzBEEQ\nBCGX9QAEQRCExkAUgiAIggBAFIIgCIJgIQpBEARBACAKQRAEQbAQhSAIgiAAEIUgtBhEdBMRba/z\nOecR0atElK/neesJEf2SiC7LehxCuohCaHCsH2LJEjhHiWgrEc21tt1KRN/z+M5riegUEb2ViD5C\nRExEa1z7vEhE76rTZaQCES2wrq3D/oyZNzHz5fUcBzMfZOYzmHky7Het5zNpPV/na04aY603loJ2\nXte49cwuUOz/BBGdcOx/wOfYy4loBxEdI6JfemxfYG0fJ6L9otCCEYXQHFzFzGcAOAfAfwD4nPX5\n1wBcQkTnuvb/AIBnmPkn1vsjAD5BRK+py2gTopVn3C5+aCkU5+ulrAeVBJaCrlwXgD8H8AsAP/b5\n2i2O7yz02e84gC8BWKPYfj+AYQBnA/gkgIeIqCv8VbQPohCaCGY+AeAhAG+x3r8I4HEAH3Lt+scA\nvuJ4/yyAHwJYrXMeIjqbiB4lot8S0Y+I6C4i+oFj+yIi+g4RHSGiA0R0g2Pbl4no89ZK5hUiepKI\nzgvx3f9FRI8R0XEAy4noSiIatsZyiIjWOYZqr47GrNnkO6wZt3Osl1jXcMz69xLHtieI6G+JaKc1\n1u1E9FrFPXmWiN7reN9BRL8more5VypE9CfW/q8Q0S+I6L/r3HePc55n3ae3We/nWOd8V9B5iOhd\n1irwr4noV0T0MhH1E9F7iOjfrOPe5th/HRE9RESbreP9mIiWKsaVI6IBInqOiH5DRA8Q0WzNy/ow\ngK9yAiUSmPkpZv4aTAXjHuObALwNwFpmLjHzFgDPALg27nlbGmaWVwO/APwSwGXW/zthCvqvOrbf\nBOBnjvcLAZwC0GW9/wiAHwDoATAGYLb1+YsA3qU45zesVydM5XMIwA+sbTOt938CoAPmj+7XABZb\n278Mc0VyobV9E4BvhPjuMQDLYE5WZgB4F4Al1vvzYa6Q+q39FwBgAB2OsX/EMdbZAI7CVJgdAG60\n3p9tbX8CwHMA3gSgYL3foLgndwDY5Hh/JYD9XuOwtp0HgAD8PoBxAG9THLcyXsX2/wZToXcC2Abg\n71xj8DyPdd8mrHEb1nEOA/g6gDMBLAZwAsAbrP3XASgDuM7a/68APA/A8Pg7XAVgF4DXA5gO4AsA\n7tf4W54PYBLAuT77PGGN89cAdkLxN+r6zmUAfun67P0AnnV99g8APpf1b7qRX5kPQF4BD8j8Ib4K\nU5hPAHgJwBLH9k4AvwVwifX+UwD+2bHdKSAfAPBp6/+eCgFA3hIMCx2f3eU4xkoA33d95wswZ2KA\nKdT/ybHtPQ7BqfPdr3rdB8f+9wDYaP1/AfwVwocAPOX6/g8BfMT6/xMAbnds+3MA/1dx3t8F8AqA\nTuv9JgB3qMbh+u4QgI8ptn3Eeq5jjtdzrn0egTm7fRrAdJ97UzkPTIVQApC33p9pjfEix/57cFq5\nrgOwy7EtB+BlAL/n+Du0FcKzAP7Ase851t+M5/U79vsbAE8E7HORNdbpMFcTrwA4L+A7XgrhQ87r\ncfw2vpz0b7SVXmIyag76mbkI80dyC4B/IaLXAQAzjwN4EMAfExHBXDF8RXGcOwD8mf1dBV0wZ9OH\nHJ85/z8fwEVENGa/rHM6j/nvjv+PAzgjxHed5wIRXWQ5Bg8T0TEAHwXgadbxYA6AF1yfvQCgW2Os\nVTDzz2EKwquIqBPA1TBn2zUQ0buJaJdllhmDqRT9xryLmYuO13mu7f8bwFthzm5PhjjPb/i0o7tk\n/fsfju0l1/VW7j0zT8GcNHg5t+cD+KbjGT4Lc+b/Oz7XCNSaMmtg5ieZ+RVmPsnMX4G5SnhPwHG9\neBWA22f2GpgKRlAgCqGJYOZJZn4Y5o/vPzs2fQXADQD+EObs6luK7+8H8DCA27y2WxyGOWN9veOz\nuY7/HwLwLy4BdgYz/5nGJeh8121b/jrMGfJcZj4LwD/CNJF47evmJZjCy8k8AKMaY/Xifphmp/cB\n+KmlJKogoukAtgD4OwC/YynyxxxjDgURnQFzVfRFAOtsW33S57GoPGciysH8G/Bybh8C8G7Xc5zB\nzMr7SkTLYCqXh0KOiRHtmvYBeAMRnen4bKn1uaBAFEITQSbvAzAL5qzM5vswTQ33wrTXn/I5zHqY\nNvyi10ZrRvkwTOHTSUSLYM7sbL4F4E1E9CEiMqzX24nozRqXEOW7ZwI4wswniOhCAH/k2HYYwBSA\nNyi++5h1vj+ynMArYfpEPBWmBt8AcDmAP4NidQBgGsyV3GEAE0T0bus7UfksgD3M/F8BbIWpENM4\nDwBcQETXWM7xVQBOwvQVuPlHAJ8iovkAQERd1t+lHx8GsIWZlTN0IioS0QoimmE9r5sAvBOm78Rr\n/xwRzYDp8yDre9MAgJn/DcAIgLXW5++H6YPaEjDOtkYUQnPwKBG9CtNX8CkAH2bmykyHTQPpV2HO\nhr/qdyBmfh5muOpMn91uAXAWTHPK12DOjE9a338FpuD5AMzZ478D+DRM4eRLxO/+OYA7iegVmCav\nBxzHG4d5P3Za5ouLXef7DYD3Avg4gN8A+GsA72XmXweNVTH+l2H6IC4BsNnnGv/CGudRmArskYBD\nv4Nq8xDebgnZK2CayQDgLwG8jYhuinieIP4Zpp/HdsRfw8xlj/0+a51ru/VcdsG0/XtiCe0b4GEu\nIqLbiOjb1lsDpr/Kdir/D5jm0gPWvr9n/Q5s3gnT7PUYzJVfCYAzKfEDAPqs69kA4DpmPux3A9od\nMmWJIKghok8DeB0zfzjrsQjpQGY47+8y8wezHouQHbJCEGogM1fgfMtEdSGA/wLgm1mPSxCEdOkI\n3kVoQ86EaSaaA+BXAD4D05wgCEILIyYjQRAEAYCYjARBEASLpjIZvfa1r+UFCxZkPQxBEISmYs+e\nPb9m5sDCfk2lEBYsWIDdu3dnPQxBEISmgojcGfueiMlIEARBACAKQRAEQbAQhSAIgiAAyFAhWPVF\nniKivUS0j4jWZzUWQRAEIVun8kkAlzLzq0RkAPgBEX2bmb2KaQmCIAgpk5lCsAqy2YWqDOslWXKC\nIAgZkWnYKZlN1PfA7Eb1eWZ+0mOfmwHcDADz5s2r7wAFQagbQ8OjGNx2AC+NlTCnWMCaFQvR39sd\n/EUhMTJ1KlsNX3pgNuK4kIje6rHPvczcx8x9XV2BeRWCIDQhQ8OjuPXhZzA6VgIDGB0r4daHn8HQ\ncNReRkIUGiLKiJnHYPa3vSLjoQiCkAGD2w6gVJ6s+qxUnsTgtgOJnmdoeBTLNjyOcwe2YtmGx0Xh\nuMgyyqiLiIrW/wswG2Xvz2o8giBkx0tjpVCfR0FWIcFkuUI4B8AOInoawI8AfIeZo7Y2FAShiZlT\nLIT6PAr1WoU0M1lGGT0NoDer8wtCM9AujtY1Kxbi1oefqRLYBSOPNSsWJnaOeqxCmp2mKm4nCO2E\nbeKwhaRt4gDQckrBvp40ld+cYgGjHsI/yVVIsyMKQRAaFD8TR6spBMBUCmleVz1WIc2OKARBaFDE\nxJEs9ViFNDuiEAShQWkGE0ez+TjSXoU0Ow2RhyAIQi1rVixEwchXfdZIJg4J42w9RCEIQoPS39uN\nu69Zgu5iAQSgu1jA3dcsaZgZroRxth5iMhKEBqaRTRzi42g9ZIUgCEIk6pFMJtQXUQiCIESi0X0c\nQnjEZCQIQiT8wjibLfpIMBGFIAhCZLx8HO2UYd1qiMlIEIREkeij5kUUgiAIiSLRR82LKARBEBJF\noo+aF/EhCIIQGS/nsRSRa15khSAIQiRUpSsANFyGtbTO1ENWCIIgRMLPebxz4NKGiSiSqCd9RCEI\nQguTZj5AsziP262vRBxEIQhCixJnZqyjSBqlPHfQWJtFcTUC4kMQhBYlaj6AblnrrEtXDA2PovfO\n7Vi1ecR3rBL1pI8oBEFoUaLOjHUVSZzy3HGdvLbSOjpertnmHmvWiquZEJORILQoUU06YRRJlPLc\nSTh5vZSWaqzSOlMfUQiC0KJEzQdI2zeQhJM3aJXjHmsj95VoJMRkJAhNiI7JJapJJ20TSxJOXj/l\nJOag6GS2QiCiuQC+CuB1AKYA3MvMn81qPILQLIQxuUSZGadtYom7AhkaHsXxkxOe24oFA+uuXiyr\ngYhkaTKaAPBxZv4xEZ0JYA8RfYeZf5rhmASh4alHXH2aJpY4pS3cytBmVqeBtVeJIohLZgqBmV8G\n8LL1/1eI6FkA3QBEIQiCD80eVx9nBaJyJndO6/DsyyCO5HA0hFOZiBYA6AXwpMe2mwHcDADz5s2r\n67gEoREJa3LxE4zNJjR1laGUq4hG5k5lIjoDwBYAq5j5t+7tzHwvM/cxc19XV1f9BygIDcbyRV0g\n12cqk4tfkpluAlrSxDmvbpKZNOmJRqYKgYgMmMpgEzM/nOVYBCFrdCKHhoZHsWXPKNjxGQG49gJv\nm7+fYMxKaMY5r1cEFMFUkk6a3ayWFVlGGRGALwJ4lpn/PqtxCEIjoGvi8BKmDGDH/sOex40iGEfH\nSlgwsBV5Ikwyo7uBiuL193Zj9wtHsGnXwYpSZABb9oyib/7syhgbpc5Ss5GlD2EZgA8BeIaIRqzP\nbmPmxzIckyBUUS8bu27kUFhhGiQYvbbZTDJX9nErJ/d9Wb6oCzv2H9a6T3GF9Y79h6tWSEDtvZIm\nPdHIMsroB0CNKVQQGoZ6OiZ1BX1YYeolGAnmtRQLBvI5wuSUW7zWUipPYtXmEQxuO4Dli7qw+alD\nKE+dVhj37TpY2TfoPsUV1jr3SspVRCNzp7IgNCr1tLHrOkvDZhE7s5UBUxnY4n+sVNZSBk5s4V8O\n+J7ffYpTFA+Q6qVp0hBhp4LQiNTTMak7a44y87WTzJZteNzXRJQ0fvcpTuLb8kVdVT4EoPZeSdhp\nNEQhCIKCejsmZxi5igDzK8EQVZjWO8ImjfukG2WVdjZ3vXxL9c4TEYUgCArq5Zj0KsdwcmIqkeM6\nhclZBQNjpdr+AWngdZ+SSJBTRVndt+sgduw/XDmnaiWUhFKs1+oji1WOKARBcOEUTmcVDMwwchgb\nL6c2Q0tjNuslTIw8wchRoP3fhgBsXNnjWTvID68wVT/hBkBb8AWFy655cK9vqEoSq5Z69WjOohe0\nKARBcOAWXGOlMgpGHhtX9qT2I0zDV+ElTMqTjFmdBjqndWB0rFTlYPbirIJR5bPQ+Q4A7By4VGs8\nTsezruBTmfFs/JSdl58hijmmXr6lLJLrJMpIEBxkkb2bRtSMSmiMjZcrkUpBgp2oWmjmiQK/kyfv\n6bmfcAsj+LyirHRxRjLVo3xGXLKIphKFIAgOVMJpdKyUWo2fNBrS+AmToPaTNkfHy1jz4N6K0LQT\n1fxQ7VPsNDw/ZwA5hRI5q1D7HXcYrS7dxYK20zmIevVozqIXtCgEQXDgN/tKq/Bb2Lh8nZpHfsIk\njMlB199gYwtq9xhP+CgglRI5fmpC2Qlu58CluGdlT801GjlCPletYIw81QjRuOUz4uRR6FKv8zgh\n1tD6jUJfXx/v3r0762EILYyqAYtNd7HgaSOvF17jKxh5XHtBd03pCOC07d9Zl2j81ASOjicfbVQw\n8rj7miUAENoRrcLrfjvNWMVOA8zAsVK5UkLDmUUNmEpi8PqlVYJUlZNhn6/ZyoIHQUR7mLkvcD9R\nCIJQzdDwKFZtHvHcRgCe33BlaucNEkIqQeZ29hKAmy6eh775s2uEc5hIozAUjBzuvub8ihJKAvf9\nHhoexZoH9yoFfs/67Z6htXkiTDFXKUsvxapSaPa2ZlUKugpBTEaC4KK/t1tpo1bZwv3QLWut4+RU\nmTTc4p0BbNp1EOse2VcbbTTFUJjtY1EqT1WuISncJrx1j+yrUWblKca6R/ZhaHhUmWcxyVx1XwEo\nzTHt3EtBFIIgeLBmxUIY+Vqp+eoJb7u2Ci9Bv2rzCHrWb686jq4QChNhwoBSQOoaBszcBe1T+pqJ\ncoSqmkq156n+1MuBqrqesVJZW2A7Q1p3DlyK5zdciZ0Dl1Zm/+3cS0EUgiB40N/bjZnTatN0ylMc\naqaoiugZK5WrVgC6QkjVICYs3cUCOgMkfbFgYPC6pfhPr0kmzHGKUXEGOyOIZnWa5xm8fimKjs9n\nhNFECCew/fZt5+J5ohAEQcExxWw0KcHjXAHoCiF35MmsTsNXcM7qNDyjjZYv6sKJsn95jLFSGas2\nj4Q2AalyEbqLhcqKyTnTd47DWbLj6Hi5xmw2S2Gym9VphBLYfv2nj5+cqPm8XXopiEIQBAVJzBSD\n9rUVhmrmPzpWqvE72KaOjSt7cKI8hZJCsBeMPK48/5wqhVEsGLj7miV4eM+LiF8tqRYjR7jxornK\nkNe4LT3XXrW4xpRn5Alrr1rseQ91TFG2j2fBwFas3jxSY5aa1Wk0tUM5DFK6QmgJ4nTwUqFb3M4v\nOsjrGE5yRDh3YCvmFAuV0FF3iQhVx7KPP7BXGcPfbd2DLXtGq859zJr1p0V5inHfroOY1WmAwBi3\nlNX0DlMp+SX+qQjb+Ma9zW9/dxiv193snNbRFsoAkLBToQUIyh0AoocNBoWCep3bFubdDoG0/tF9\ngbH/9hhVYZvFgoGRtZcHXq8dqplW/wM7pyEsBSOP6R250BVX08z90LlHaYYa1wvdsFNZIQhNj04p\nhqhVIoN6D6jKMQOnZ/Z3X7MEw3dcXqVcch5CtVSe9J31j5XKeMvffLsy61Zhm6nSiIopFgylbyWI\nUnkSM4wcCkZeO2ktbdu9zj1qB2eyjfgQhKZHV/ClISCDjqmKX1cJ/aCZd5AycArQNATZ8VMTnjWG\ndBkbL1c5xf2oR6mGoHvULs5kG1EIQtOjK/jSEJA6x3zJKoznzEdIA4IZqrl68wiWbXgcyxd1Ra4M\nqqI8aSa1RT1ujgirLR/GTRfP841IcuYGpIVfGG89FFKjIT4EoemJ4kNIqlaNzrntZKw0+xkbOQLI\nFNg2zhpHur0MdLAb5+j4RaJg5AmD1y2teR5JPrOkAxAaHfEhCG2DV+SJ3488amtCP4GkaiBjmxxW\npxjZkyfCGTM6aoRzqTyJHfsPVxyy9vjjKqY5Vinp/t5uDA2PeiqGWZ0Grjz/nMoz8PKZKPHYLal2\nkl7H2bJntO1WAipkhSC0Hb13bvec2fpFs6iqjLoFiVPoJllhtGDkcHJiCu6adPZsevXmEc/Zv1eE\njN+qJp8jTAYUvvvgxfNwV/+SyrHcxeac47LvzbkDW0OtTtzPIqg6qS6q47iL37WacmiK4nZE9CUi\n+hUR/STLcQjtw9DwqFIw+zmIdWsN9fd2V+zS9ox4dKyEV09MeNZGcpMnqhRbu2dlD3654Urcs7IH\nM4x8jTIgy0Q0uO2Asuiel4/Dznb2Mt8HKQMA+Nbelyv/H9x2wLNyqj0uv3H44X4WUfIXdI5r4y5+\nl1YzpEYna6fylwFckfEYhDbCrw6RLbSGhkfRe+d2LBjYigUDW9GzfrtS8DgFjJ3xumrziGeF0Y4c\nBXb6OnNGBzau7Kk4VIeGR7Hmob2eSsxe3KsUTlCETFTjwFipHFiDyR6XvV/Y1pc5oiqhrFIoBIQS\n3jqKqV0qm3qRqUJg5u8BOJLlGIT2wk+ArVmx0FMA+yVSOZVIUOnnUnkKx09OKOvx2OdyzlAHtx2o\nchSrKE8xZk7r0O6uFVfgBdVgslm9eQS3Dz1TWZWooorcTDJX3Yc1KxZ6hqkywl2LrmJqh8qmXmS9\nQgiEiG4mot1EtPvw4cNZD0fICJ2eAjqoBFixYFRq4asEsFsgOWfgun2Kx0rlQPORnaA2NDwaSjCN\nlcoYt/IEXhorYXDbAeV9iivwnDWY3LWCnDCA+3YdrCiFqRDLEudMvb+3W+mDCHMt7uKAKgXVTslo\nThpeITDzvczcx8x9XV1dWQ9HyADd5jE6qHoNr7t6MQB/4WKXowBMQVIqT2LdI/vQe6fapOSFczav\nwp4hh23Ic3S8jLFSuar3Qu+d22vuVZRGP05sgdnf243B65ciyD1y366D6Fm/PXTYq/N5qO5XWOHt\n7IPwmRuW1r2RfSPT8ApBEJLsYBXUuNxPuOSJKsletsN4rFSOFD10rFTGzoFLfZVCqTyZSJy/Vxnp\nOMGFXgJTw6oVuoYRUP08VMo8jPB2rjR71m/H+kfNjnL2SqEdk9GcSB6C0PDoOnR1k5ZU9YlUtfBt\nJpmxadfBRJK7GKhkE2/+0SEtP0EcSuVJrNo8gsFtB7BmxcLI9YiKBQPrrl5cU100DYw8VQl7nUqn\nfrjDbZ0KapK5olzaVRkAGSsEIrofwLsAvJaIXgSwlpm/mOWYhMZiaHhUmWHrdujGSVrSyTiGYhxR\nsZOiOnKUukJwnlOVs6CDs4GNTRx/BAHKpLWZHmWng4oN+hHk54laALGVyFQhMPONWZ5fiE9S5QRU\nrH90n1J4LV9k+pT8TEq6Y9F1CkdhVqeBzmkdniudtM7pRxzVUypPYv2j+yrPPI4vwk4qWzCw1XN7\nFBOTHzqKq12ji2zEZCREFupJlRPwO76fDX3HfjPqTKcfsdc1AqfND2GEZNiaQGPjZQzfcXnobN1G\n5ej4ab9JVB+H0/av6q+gG6Kqw9DwqFb5jHaNLrIRhdDmhBXqOjX9k1p2B9mmbYE/p1jwnH37mZTW\nPLQXYHhm2frhLhinQ44Itw89A6J4ztxmIU+E6R2kLNWdJ6py3EYtBa6L/fyDjtfO0UU2EmXU5oSJ\n4HGHf6p+YEHLbt2cgqDj2AI/KPrE6xrLkxxaGdj9iO/qXxJKcEyy2VZSdboEJ8INwRSzss8zAHzm\nhtM1jmwfkRdBWd26qMyBBFitPiW6yEZWCG2OjrnFRtfO7rfsDrMiUc38gWqBHxR9kpRd2OlQTTKy\nhjm8GaqedBq5wMY8Tuznr2oD6o5QUl13UrN1v+c/fMfliZyjVRCF0Ob4mVvcdncdE0nQsjtoReIu\nYe1uEg+Ys7q1V1WHPvpFn+iOPYhSeRK3Pfy0b5vLKHQXC1hwdgE7n4tXxSVqr+N7VvbgLx8YUa5g\nTk7oH9PInQ4V9aoOaycA6pTi3v3CkURm7EEmRT/SDppoNMRk1OaozC3LF3XVZAerlvbOCp1By26/\nqpXu823adbAmaeielT0YvuNyz+Yp7oJ0zjo4OpVGjRwh71OGATBbWCapDOx7/dTzR2Mfa5I5sC2l\nm26rt8GMDrUoCHO9Rp4qylmVAKhT9wkANu06mEjV0agJbUlmyDcLskJoc1TmFlXzeK8GMGFsr6rZ\nml0Kwn0+IDhpyC5I54zlHyuVsebBvZVrXPfIPt8wxjwRVl44F33zZyfSREYHO8nrtoefDu3P8KLb\nWlXdt+ug1v72PR0aHvU1CYVxho+XpzA0PFpRCl7PS9f0aBeuizsjj5rQlkQ4c7MhCkHw/OGqOnzZ\n9XxUP6ygJfaaFQs9TQlBAsLvh6gqSFee4sp3gjJzJ5mxZc8o+ubPxs6BS7UT1eIwsvZy3D70TCj7\nvIocAcdPTmCTpjIAgGsvMO/l6gfU3dwI5srpVIjEuY8/sBerN48oBW8YZeu1ooxixomS0BbGv9Yq\niEIQPFHN5MN0FfNyGPutSIIERZQf6OhYCecObNWKQS+VJ7F680jFnk4AZk7LY/zUJOYUC3jpWCmx\nsNHuYgG3Dz2jPZsPYorDJ3Jt2fNiYCmOS86bHdq34WwMdOvDz2D3C0cqrTTDJrK57fxp5764zx3V\n99CsiA9B8CSK3TVMVzG72qTdCMbOOvZD9UMM+oH6hch67WtbbxjA8VOTuOniedg5cCluumie1jGC\nKBh5LDi7kJgyiEqpPBUY2fTL38SbDZfKk7hv18GKHT5MIpvX31uShQ6DSKKYXrMhCkHwJKgqqBdx\nlth21rEKgncY4tDwKMZPqQvSJcH9Tx7C0PBo4Bh1ef2sGbEjiurBrE6jruaRe1b2BP691dOME+U3\n0OyIyUhQEtbuGmWJrRN+CJiz9dWOap3OaJW06wHZvQmSOs/PfnU8keMEYTeOP6tghDYnGXnCleef\ng/ufPJRoVJUfzmerol5mHLefYuPKnpZWBDayQhASI+wSWzf80MYd+qeKVtGtgVMwctpZwlkUoYsD\nwcwIfn7DlRhZGy75qlgwsPLtc7Flz6inMjBypBXGGxadsM56mHHaMdzURhSCEIhuqYmwS+yoFUZt\nm7HKTKA7o509czo23tATOna/GbBDNu1nFaYMxMjay7Fj/2Glsh28fikGr1taOWbY+2fkCcWCt3M5\nyB9QDzNOPf0UjYaYjAQlQ8OjWP/ovipHYFBURxgzk5/dNyjrNkz2tN8x+nu7sUoRYtsM5Mh0gnvd\nL+ez8gr39cIW8qpnM8VcEzHmNPvZ41DlLhABg9eZtYxU1V/tc6vCS+P0RNAhip+iVTKaRSEInvjZ\n5905AUE/BtV2ndDWZRseV9qMdYWcCtvu3J1QaYu4RKlnNMWmMxaAp2Kz+z7PnN5RyfqeZEaxYOD4\nqYma/I2jx09iaHhU+WzOKhhYtuFxTyHtRBVSe8kbZle1K/Urm1Kv8FKv84fxU2Q51qQRk5HgSZA5\nxzmL87O3+m1XhZr++tWTFfPUgrMLNSYJZ9by3dcswawITVpyOB215GWXzoJLzpsd6XurN48oEwkB\nMz/BFnB21ve6qxdj5dvn1uw7Xp7Cmof2et53I0c4fmpCy7auisj61+eOVJUUUfkDsjTbhPVTtJKJ\nSRRCm6PyD+iWntYpVqfa/q29L3se++TEVEXg7HzuSNWsmWBm2DrNFp3Twi90p2DOqJ2N1rOkYOQi\nx/wzwq0s7PuvEtrlSfa872ziSNGgAAAgAElEQVRt8zqWk6HhUeWKy/ZtAP7+gCyzhMP6KVopo1lM\nRm2M31JXt/R00I/Bb3uUYEbG6dnn0PBoYI2iIJJu0xiFHIC7rznfd5afNGGFFQOYUNRbcnems/+G\ndM6t8gdknSUcxk+R9ViTRFYIbYzf7F1lRrGbxDjtwF7Yn0fNLvbjpbGSWdDuwb0NIdDjMgVg3SP7\nYvUnDkux00Auoc48OaLKylInckzn2TdTlnAzjTUIWSG0CV6OXb/Zu26FSFWxOqd93r3dyBOOn4ye\nXTynWDAL2iVQIbRRqLdiOzZeRvySeiZ24h4QvPIgmCvRZRse943EiVqhNAuaaaxBEDdRk9e+vj7e\nvXt31sNoOrwihgpGHjOMnGdtGb8Cdqrj60YZFTsNvHpiIrIwt8ttr948klqHsaiNZhqVgpHDCY26\nRXGxQ1ZVpsa4pdOF6BDRHmbuC9xPFELrowrd9IoVT/NHOjQ8GrvbWI6Ad7xhNnb94mhqQvuDF88L\nrALaLPgp/qQhABtX9qQ6+RCioasQxIfQ4vhGfLgknts/kCS3Dz2D1ZtHfIW4TvjoFAM7nzuS6gx+\ny57RVEoz2OiW1oiL/TzHIiqDPFGokN45Vvc1rwgd1RiSjsTRzaoXvMlUIRDRFUR0gIh+TkQDWY6l\nFdGJ+HAyc3pHaiuDoBk3Abjy/HNSFcS6lMqToRrChGFWp4EbL5qLgE6diWA/z7MUZSKCmGTG2qsW\nw9AYrNNv5FXePI3gAjftXIMoKTJTCESUB/B5AO8G8BYANxLRW7IaTysStlZQWnHTg9sOBJpf7HBS\nZ42ces2k68nR8TK27Bn17WGcFLbztjwZzX1MMIX7GTP8Y09mTgs2M9YjEqeVEsSyIssoowsB/JyZ\nfwEARPQNAO8D8NMMx9RShBXwXrO1qGUpoozDjm5yfv+8Wx9rKQcvUN/KqXFKcjDM5+tncvrgxfNw\nV/+SwGPVIxKnlRLEskJLIRDR1wB8D8D3mXl/QufuBnDI8f5FABd5nPtmADcDwLx5yXSsahdUCTOz\nOg2cKE8pQ0Vtgmq06NZw0S1C56WQbrxobqzOYkaePPstC3rc+vAzKHYaSoewjjKwSbsoXSsliGWF\n7rr1/wA4B8DniOg5ItpCRB+LeW4ve0DNL5eZ72XmPmbu6+oKbrMonEa1TF971WKt1Pw4ZSmCxuHG\nyBHGT03UOAP75s+OZdcUZRCPUnkSzGiKxKtWShDLCq0VAjM/TkT/AuDtAJYD+CiAxQA+G+PcLwJw\nVtd6PYCXYhxPcBG0TA+arcUpS+E1DlWZCSMHlKe4Mgt1Nme//6lDiSVQCdEYK5Vxz8qehk+8aqUE\nsazQNRn9PwAzAfwQwPcBvJ2ZfxXz3D8C8EYiOhfAKIAPAPijmMcUXMRZpgctwVXbGUDP+u1Yd/Xi\nimlpcNsBHCuVUSwYIALGxsuYUyxg+aIubPIwCZXKk/j6kwfRQsnI2jRaYpy9lPfKF2i0PgBpm6Va\nHd3V+NMATgF4K4DzAbyViGIZ5ph5AsAtALYBeBbAA8y8L84xW42sY6q9luBO0874qQllSOJYqYw1\nD+7F7UPPVIUCjpXKOFGewsaVPdg5cCl27D+sjEBqR2Vwz8oeTDWQMgCqK5Q6kTDP1kNLITDzamZ+\nJ4D3A/gNTJ/CWNyTM/NjzPwmZj6PmT8V93itRCP82NxJRsWCAZAZOskw/wVB2Ze4PMW4b9dBXz+D\nRICcxr6PSTlBZ3UaoVpn+jFqFRR0ImGerYeWQiCiW4hoM4ARAP0AvgQzf0BIiTR/bGFWHs4ko5nT\nO2qctOVJ9myVGIStCFTCj2DW4MmKZefNTkyY6sJsRvUsX9QVu8+zkSesvWox1qxYmFiyn3tCImGe\nrYduHkIBwN8D2GOZeoSUCfqxRbXdxmn3l+QP3c6eVbXB7JyWx/FT2TWt2fWLo5hiRsHIoVSun1u7\nVJ7E/U8e0q6jZOTMpveAvzPV3Rs76ticrVMlzLP10DUZDTLzk6IM6ofqR8UAeu/cjjUP7o1kTtJZ\neahWEEn+0MdKZSwY2IqPP7AXb5t3VpVZyshTpsoAMMs2MKBUBtNTzDTWdSjP6jQweP3SGkV+/OQE\n1j+6r/L8AGD4jst9Vx26qwjnpCBOmGfW/jHBGylu16D4xe4fHS/XlI/WNSfprDxUvos0eg9Pstmu\ncfmiLqVZqhE5NTFVd5OSGzuD2P3Mxkrlip/H+fxUCn1Wp1EpGWLnpRQV9Y+cxwjbatKmEfxjgjfS\nIKdBccZU65Yf0DHpBC3z1z1S21/YVjZ22GHcEtZe3P/kIfTNnx2r1EI9mVMsZG4rZ5h2/ekdOd9y\nGPbzW76oyzPr+8rzz6kJ11T10HDP/qOEefqtUiVkNFtkhdDA2A5dXZdg1NaEdgezBQNblZ27bOHX\n39udSliks+uWFzr3oLNOTmhbMDaCrbxUntTqtjY6VsL9Tx7y3Gb3qHYSdfavgzijGxdZITQBOrWA\ndG237mxOu4NZkFBxCj/d2kRh8Zvl6qig6UYe4yk7gIsFo5Jw9/kdP0v1XEmjWtXZIaVuYZ9Wkpc4\noxsXWSE0AapZfbFgRJq9OUNJO6d1aLWzdCqbsLVh6lXEOmojmLDY9/lnvzpel/PVgyRs+LqOYqk5\n1LjICqEJiFKjRRWW6v5cZ6Y/q9OoOld/b7eyLpGNs38uO953FwtYcHah0gIzT4QbL5qLHfsPx151\n1MMVPVYqe86m06RYMPDKyQlMppi6HcWG79cr2y+cWWoONS7SU7kFUTkEr72gG1v2jIaqx6/qsTw0\nPBqp0X234sfvNeZGxe4DvGBga6rnse99mMCCOBCA5zdcqbWv7vOSnsmNgfRUbmNUURz3P3kolMD1\nM0X193ZHmpH7hRimGdufJLbzM81un857Xy9naxgbvm43PnEUNxdiMmpBVD9C3VDRYsHAyNrLffeJ\nY292myeaaXUAnBacaaVLFAtG1aw6KSe+bbbzyr4Oa8PXFfTiKG4ummNKJvjiduYVO6M1Vbc5FhBx\nZAvwODgFStjez1lSMPJYvqirkgEcBwI8q8USVStcnaCCWQHPvLtYwMaVPbhnZQ/cbn4CcO0F4SKK\ndAS9OIqbD1EITY5X1uerJyZqShGEsW7kiHwjRZIQ4E6B0izJaATgbfPOwpY9o4mMmQEMXr+0Jiv4\n6HgZax7aW7n3dk6AU+jPnNaBdVcvxvMbrsTOgUux9qrFymds2/H7e7s9nx3DOxfBj6Qj34TGQExG\nTY7XD7w8xSgWDMyc3hEqmsjGNi05I0Xsc71kKZ44OGeOtw/FW2nUEwYq0VFJ0F0soL+3G+sfrW0D\nUp5k3Pbw01UC9YTDzDNWKldF8fT3dmP3C0ewadfBqufjnqUnlRQmkUKtiSiEJkf1Qz5WKlf8AEPD\no77lJooFA8dKZeQ8OnWVypNY/+g+vHpyIpEaQ3miqpmjKnu2UUmyZIctqFVVSMfLU5VVgtfzc/ti\n7upfgr75s32FdJJJYdKdrPUQhdDkBP3AbZOSnyCzfQaqfXTKJjvzDvyYYq4SIvVsFak7Rj+Sam/p\nzu1Qsf7RfThRnlKe06t/td9xvcqNi61fsBEfQpMTlPWpY+9nxBeUut93z0TzqnZrCUIAZk7LJ2Lq\nuvGiubEbzhSMPNZetRhAcLTW0fGy7/Nz38+gbOE0axQJzY+sEFKkHg3Ig2y5jRQHnqPashc3XjTX\nswJnkjAQq78CAZX7CgCbn4pu5nIn5sXpgOee2es2PxJTj6BCFEJKxOlMFha/H3hahejCYuSAwet7\nmk4QuTNte9Zv16r95MUHL56Hu/qXVGbxQQ56I09Kv02OUDOzl7LS9ZmEtTKiEFKiUX6cqhaV9aYj\nf9qsZf9oG0FRBTE6VsKCga0oFgy8d+k5WqWmVdgrC53yIXkizJzWoTzfa2bU+iDavax0PSdhrYr4\nEFKiUX6cbpux3aKy3tjRSs68iWZirFSObdoqT7FW+RAjR/jMDUt9EwS9tqkihdolW1inPazgj6wQ\nUqKRar57dcOyl9UgoF6BPkfHy54d2doJrQglS1/7mfu8/o7aPYKoUSZhzYysEFIizZrvcRuUO/sh\n1KVmtANdk8sHL54XWI7BjxyhbquiMIFSOlFV5UnG+kf3Yc2KhZ7jNnLk+XfU7hFE7b5CSoJMFAIR\nXU9E+4hoiogCS7I2I2n9OJNuUN6oP5atT78cWFPJjylGVdP5mdM6MKvTCGwiH4RbPBeMPC55w2yt\n7xo5s/eDe6LghZ37MXjd0irFWCwYGLx+qfLvyKns7XIV7YI03olPJv0QiOjNAKYAfAHAXzGzVpMD\n6YcALNvwuKcZIWrd+WarNBoHZ2+HKNdt5AgrLzSb+TijWHQc5GSZ5uzENve/XkgvgfBIlJE3uv0Q\nMvEhMPOzAEB1SEpqNfzspGF+DM59Z4RoTt+RI0xohl0mkRmcJKXyJD7+wF6s3jyCOcUCrr2gO1Sn\ntpUXzsVd/UtqPl+9ecT3e4TTfhpb+E8yg+DvUxDbd3gkxyIeDe9DIKKbiWg3Ee0+fDhcRcZWRGXi\nKXYa2qYkt9nJXRvfD11lECczWMekEpVJ5sr9uW/XQRw/OaH93S17Rj3vZ5DZTXUfgu5Po5rzhNYl\nNYVARN8lop94vN4X5jjMfC8z9zFzX1dXV1rDbRpUdlJmaIfc1aP/wKmJqVBO4ZnT8lW+lnoxVipr\nlwZX3U+vZxKXZrB9xw1uEBqP1ExGzHxZWsduZ1SlKlRmi9GxUiUr1t63HqaI8hSD2RRsOsrnRHkK\nG1eezmRe/+g+z6J6nUYO5SlOpPKqTZgjed27oDGHwVkmo5FNH5IE1ppk4lSunJzoCYhTORFUzmY3\nRp5wxvSO2IJLBwKwcWVPlfLyG6Pb6bvmob1Vgj9HwIyOHMZDmLiSJsjRGycLu5mcyEkHNwjpoutU\nzirs9P1E9CKAdwDYSkTbshhHK6FrtihPMk6WJ1O109vMsRrA7By4FBtX9gTa650mmf7ebgxet7Qq\nlwBAjTIw8hQrXyEMOmYc+3rDhkukaSJKw7QjSWCtSSYKgZm/ycyvZ+bpzPw7zLwii3G0Eu68Bz8h\nOV6ewt3XLIkci6+DU8DZ5gWdpDSnQHHG1M+c3gEvf3Z58rRpKk1mdRqh8kjO8rm3BGDZebPrkkCW\ndN6KjSSBtSZSuqKFsEPubCGgu69t0vHqmBaGHJkJYV4lnnWd2CqB4jfzPFYqV0xTSdZIimPPV0VU\nEwEbb6hf1de0iiy2e5mMVkUUQgsSJICdKwNn3HbcJLUpPi0UnMJGV0gTavsl2ATV9bHPt+bBvTXl\nqY0cENbtENcWPqby0XB9na5pmXakp3JrIgqhAYmabanj0DRyhHVXL67a33meu69ZEitaxmv2qdt2\nkmEKmFWbR2oyeWd1GsjBTG93knfU9Vn3yD7PXgVhlUESM12VAit2GjVRX2kK0TSLLEoSWOshCqEO\nhM0gjhLOpzu7L08xBrcdwO4XjlTV5bfPc+0F3TgRM4rHPfvUNUMRTq8mnBm9gFnbJ2/bpBw4nWBh\nexXkySwzDVSHjE7viO9a8zKpGHnCqycmKudxP9s0yi6IaUcIgyiElAkr4KPafMPY6UfHSti062BN\n/H2pPIn7nzwUu4l8jghDw6MVIWfX8fFDp8zFpOfsnyPXu3dep1MJjpXKvs/IuRKzVzBuv4mXSeX4\nyYkapeWMrEojrl9MO0IYMs1DCEsz5iGEjdc+d2Crp2AkwCxXrUD1PUDfZJM0BSOHCY0ksk4jfm6B\nbgKc1/emd+Q8Vxdez8hvJebMo/DC79mqTDsS1y8kQUPnIbQTYZ16UcP5VNu7iwVMJaAMopQhLJWn\ntDKKS+WpWCGweaLIjvBSeVJpavIS0H4rsaDuXKpndFbBkLh+oSEQhZAyYQV81Jruft9TnUtXyBeM\nPG66eF5qeQt2z4KopLX68WpmEySg/bavWbEQRq72mMdPTSjzFiSuX6gnohBSJqyAj9pYx+97qjFc\nct5srQ5ed1+zBHf1L8HI2stxz8qeyjlaHS9FE3WlBpjP6IwZtW678iSDqDa5Tpy/Qr0RH0IdiBM9\nklTkifs4yxd1VUUZqfCzYevWT2pk6ulDAPz9CO66T+L8FZKioRvktBtR47WTrCjpHsOyDY8HKoOg\nGapXSKMfHlGjmdBp5FAqT1UU47f2vlyzj+ranVE7flFGKvzyAiSuX8gaUQgNTFplBwD/7GHdkg39\nvd3Y/cKRmhBWO4TUHUqqqwzSjoqaNXM6fjpwqXK2P6vTwNqrFvv2LY56/yUvQGhkxIfQwKQVeTI0\nPKr0AczqNLBxZQ8AszVkUHXMHfsP15hAGKZQjyLSC0Yen7nBrHIaBz8fx6ij3ajX6qZzWkdqM/Wo\nPiJBqAeyQvAgrN0+rcbeKvPCWYV45Q8Gtx1QCutj4+WqPgRBZiqVctKd4Rt5wsxpHThWKtdcS02m\nb45wxozgXg4FI49rL+j2TbLzM3WlHeoppiGhURGF4CKs3d5vfyBahqgzE9ZtdjFyhOOnTme8RvEr\n+Am8KQBTrtwBPzOVSmnpmH3yRFj5drNx/dDwKNY/ug+rNo9g1eYRFAsGrr2gGzv2H665f16mHvs+\nOW35ffNnKwV/qTypHKNf6WpBaGVEIbgIa7dX7b/ukX04OTEVuyaRU1wVCwaIUDNDDutXCOpc5oVK\niahs4tde0B0YxTTJjC17THPU5h8dqkpiGyuVsfmpQxi8fmnNdemWY7Dfr1K0F51khpGjmoJ4x09N\nVEpvAOmtAAWh0RAfgguVoFR9rhKUY6WydtN7J36ZsCcnppTmEl0zx9DwKMZP+Xcu80IVX6+yid/V\nvwR3X7Mk8Lh2/SSvjOYwdYp2v3DEsytYf2+30h/RXSwo8wJWWf6T24eeSaXBjCA0IrJCcKEyI6gS\nuMLOtuNkuvqZOXQyWqP2O7CjYFQzZZVNvL+3W6tpjZ9pyet+eJnp7tt1sLLdvRrzi+xZrVg92MdR\nFQG0FZWsHIRWQlYILlTCSfW5KgtY1cIyTqarPY6oGa1hKqLmiapm/ABqZsqrNo/gzX/zbd9evWtW\nLAzMavbLlva6HzrX4e7PrIrsCbrfKlVlKx1ZOQithKwQXHT7VJ30QmXPBmojWXRrEvnN4m2naZSZ\naZjomSnmquqqqkS2klWlVOUjUeUq2Nj+BrcPwWb5oq7I1+Huz+z2CazePIKzCgaMPGkV4XPiVVCv\nVJ7Exx/YWzmfIDQbohBcREkc8gsjDCu47e3rHtlXU07B2Z4yisAJY95yz5x1hLDKuX1X/xL0zZ/t\nm93bN382bnv46Zoy2Fv2jKJv/uyqY+peh9fs321uGiuVYeQIszoNpX/GHenlV2p7kjmRPgaCkAVS\ny8iDRokqCRpHlHwJHR+CV02eMHWLfunTt8EP3d4ROtdBAG66eB7u6q92bPudwy9iyh3+GuQbkT4G\nQiMhtYxikHTiUFQF4zeOKHWO3HV4vMgTeWbO2g7YoOkDWWOLcv90M7O9zHQLzi7gX587Uhkfw3t1\n4XeOsN3FskxuE4Q0EIWQMkkWqLOPpxLoOvkItpLxmmX7VesM8gXYMBC51lJRYbZxtuR0X4fNsg2P\nK6OBdMxNtnlJdzJg7/PxB/ZGjvoShEYjkygjIhokov1E9DQRfZOIilmMox74JbqFxRbifqYK3Zlp\nlJo6d/UvwUaNfghRZsdDw6N49YR3foRtl/eL4NFdXURtQORFf283PnPDUuljILQMWa0QvgPgVmae\nIKJPA7gVwCcyGkuq6AgqXZOSTrhlmJlpGNOYe4x27X6/2XYYBrcdqMkYdhK0+gma+dsk3XRemtgL\nrUQmCoGZtzve7gJwXRbjqAdBgsrLpLR68wh2v3CkxiEaNPMOMzMN49dQmb28ylNEnR3rrCqC2lPq\nRocl7SOSYnVCq9AIiWl/CuDbqo1EdDMR7Sai3YcPH67jsJIhyEThNetnAJt2HawxkfjNvMOUUXaa\nnnSSqlRmrx37DydWyllnVRHUnlLKSgtCPFILOyWi7wJ4ncemTzLzP1v7fBJAH4BrWGMgrdBCs9hp\ngBmVcs9hQhfDOoJV6IZ32vi1fXzeI8Q0SlRVUChplOsUBMEk87BTZr7MbzsRfRjAewH8gY4yaGZU\nkT1e5a2d6IRbRrFXh228o2ufB6KVD3cqy+kdORwrlWsUp9jlBSF9MvEhENEVMJ3Iv8/M41mMIQtU\n5iEVc4oFz9m2TsKT3yw9jIAHwtnnw5QPdyuPo+NlFIw8Nq7sEeEvCBmQlQ/hHwCcCeA7RDRCRP+Y\n5smGhkc9SyPXGz+nqDuMs2DksXxRV6QCakE+Aq/6QH6fh7HPh1l9JBmSKwhCfLKKMvrdep0r6cSw\nOKhm5qqCdTqzba+VQND3duz3ds6rPgf0I2nCrD7S6hktCEI0GiHKKFUaaRbqF3HU39uNnQOX4vkN\nV2LnwKXo7+0OFJiqlUBQM580BXGYxC+ViUqyfAUhG1peITTSLDRsaGSQwFQpO79mPjrHjUOYa0wy\na1gQhPi0fC0jlQkjR4RzB7ZqJWVFCaFUfSdMElOQM1el1OwmOqrvRSnxHYaw9YAky1cQGoOWL3+t\nUypZFeMeJe7f6zuqUsy641cJzKBSznFLZzdKGXBBEOKhm4fQ8goBqBZsOUVPYq+krLAJXH7fISDx\ncMqkEtXqfWxBEOqLrkJoeR8CgCqH7ZRCAXqZX6L4H1Tb7LLQSZJmuYZGcsYLglAfWt6H4CZMWGTY\nBC6/7wDpOLLTKqzWSM54QRDqQ1usEJyEiWyJEgWzZsVCZa+AZgqnlJBQQWg/2k4hhDGzRDHJ9Pd2\n46aL53lmHjdTOKWEhApC+9EWTuUsaIUInVa4BkEQJMpIEARBsMi8/LWQDTKrFwQhKqIQGpCoQr2R\nCvkJgtB8tJ1TudEJ297SieQOCIIQB1khNBhhGsy4aaXcATF9CUL9EYXQYMQR6lES6ZIkKSEupi9B\nyAYxGTUYcRLCsswdiGPqciOmL0HIBlEIDUYcoZ5mbaMgkhTirWT6EoRmQkxGDUbcHgFp1TYKIkkh\nnrXpSxDaFVEIDUhWQj0OSQrxtBv4CILgjZiMhERI0n+RpelLENoZWSEIiZB0O8xmXCUJQrMjCkFI\nDBHigtDciMlIEARBAJCRQiCivyWip4lohIi2E9GcLMYhCIIgnCarFcIgM5/PzD0AvgXgjozGEYqh\n4VEs2/A4zh3YimUbHo+UdCUIgtCoZOJDYObfOt7OhNmDvqGRcgqCILQ6mfkQiOhTRHQIwE3wWSEQ\n0c1EtJuIdh8+fLh+A3Qh5RQEQWh1UlMIRPRdIvqJx+t9AMDMn2TmuQA2AbhFdRxmvpeZ+5i5r6ur\nK63hBiLlFARBaHVSMxkx82Wau34dwFYAa9MaSxJIOQVBEFqdrKKM3uh4ezWA/VmMIwxZVhIVBEGo\nB1klpm0gooUApgC8AOCjGY1Dm6QzcQVBEBqNrKKMrs3ivHGRTFxBEFoZyVQWBEEQAIhCEARBECxE\nIQiCIAgARCEIgiAIFqIQBEEQBACiEARBEAQLYm74unIViOgwzLyFOLwWwK8TGE7WyHU0Dq1wDUBr\nXEcrXAOQ/HXMZ+bA2j9NpRCSgIh2M3Nf1uOIi1xH49AK1wC0xnW0wjUA2V2HmIwEQRAEAKIQBEEQ\nBIt2VAj3Zj2AhJDraBxa4RqA1riOVrgGIKPraDsfgiAIguBNO64QBEEQBA9EIQiCIAgA2lQhENHf\nEtHTRDRCRNuJaE7WY4oCEQ0S0X7rWr5JRMWsxxQWIrqeiPYR0RQRNV24IBFdQUQHiOjnRDSQ9Xii\nQERfIqJfEdFPsh5LVIhoLhHtIKJnrb+nj2U9prAQ0QwieoqI9lrXsL7uY2hHHwIRvYaZf2v9/y8A\nvIWZG75JjxsiuhzA48w8QUSfBgBm/kTGwwoFEb0ZZqOkLwD4K2benfGQtCGiPIB/A/CHAF4E8CMA\nNzLzTzMdWEiI6J0AXgXwVWZ+a9bjiQIRnQPgHGb+MRGdCWAPgP5mehZERABmMvOrRGQA+AGAjzHz\nrnqNoS1XCLYysJgJoCm1IjNvZ+YJ6+0uAK/PcjxRYOZnmflA1uOIyIUAfs7Mv2DmUwC+AeB9GY8p\nNMz8PQBHsh5HHJj5ZWb+sfX/VwA8C6CpulmxyavWW8N61VU2taVCAAAi+hQRHQJwE4A7sh5PAvwp\ngG9nPYg2oxvAIcf7F9FkQqgVIaIFAHoBPJntSMJDRHkiGgHwKwDfYea6XkPLKgQi+i4R/cTj9T4A\nYOZPMvNcAJsA3JLtaNUEXYe1zycBTMC8loZD5xqaFPL4rClXm60CEZ0BYAuAVS5LQFPAzJPM3ANz\ntX8hEdXVhJdJT+V6wMyXae76dQBbAaxNcTiRCboOIvowgPcC+ANuUIdQiGfRbLwIYK7j/esBvJTR\nWNoey+6+BcAmZn446/HEgZnHiOgJAFcAqJuzv2VXCH4Q0Rsdb68GsD+rscSBiK4A8AkAVzPzeNbj\naUN+BOCNRHQuEU0D8AEAj2Q8prbEcsh+EcCzzPz3WY8nCkTUZUcKElEBwGWos2xq1yijLQAWwoxu\neQHAR5l5NNtRhYeIfg5gOoDfWB/tarZoKSJ6P4DPAegCMAZghJlXZDsqfYjoPQDuAZAH8CVm/lTG\nQwoNEd0P4F0wSy7/B4C1zPzFTAcVEiL6zwC+D+AZmL9rALiNmR/LblThIKLzAXwF5t9SDsADzHxn\nXcfQjgpBEARBqKUtTUaCIAhCLaIQBEEQBACiEARBEAQLUQiCIAgCAFEIgiAIgoUoBEEQBAGAKARB\nEATBQhSCIMSAiD5q9dUYIaLniWhH1mMShKhIYpogJIBVR+dxAP+TmR/NejyCEAVZIQhCMnwWZrMi\nUQZC09Ky1U4FoV4Q0UcAzEcDl1EXBB3EZCQIMSCiC2AWJPs9Zj6a9XgEIQ5iMhKEeNwCYDaAHZZj\n+Z+yHpAgREVWCIIgCAd6kTgAAAA5SURBVAIAWSEIgiAIFqIQBEEQBACiEARBEAQLUQiCIAgCAFEI\ngiAIgoUoBEEQBAGAKARBEATB4v8DcdaG5Bhddm0AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "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": [ "covariance matrix for z and w:\n", "[[1.03751418 0.69336422]\n", " [0.69336422 0.94874494]]\n", "\n", "estimated covariance of z and w is 0.6933642208178616\n" ] } ], "source": [ "est_cov2 = np.cov(z, w)\n", "print('covariance matrix for z and w:\\n{}'.format(est_cov2))\n", "\n", "print('\\nestimated covariance of z and w is {}'.format(est_cov2[0][1]))" ] }, { "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": "iVBORw0KGgoAAAANSUhEUgAAAz4AAAF1CAYAAAApyz9eAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzs3XmYnWV9//H3d2ayQFiyEAjZWMOm\nyDYFLIoLIAEtYAUFpKCl8rNCpVKtYC1a1LoWKhapKLZoVUREiBoXiijVCiUgoqwJAZKQQAYCCZB1\nZr6/P86TyZmTM8kkmSVznvfruuaaZ7mfZ+6TEM58zn3f3ycyE0mSJElqZE2D3QFJkiRJ6m8GH0mS\nJEkNz+AjSZIkqeEZfCRJkiQ1PIOPJEmSpIZn8JEkSZLU8Aw+GnIi4omIOHaw+7GpIuJdEfHrDZz/\nSUScM5B9kiQ1pqH6Xin1J4OPtkhEnBkRsyLipYhYVPzy/prB7ldvRMTrI6Kz6PuLEfFIRLy76vzD\nEfGXda67MCJmFduviIifR8TzEfFCRNwTESduTn8y84TMvK4X/c6I2HtzfoYkaeA16ntlROxevCe9\nVHw9ExE/iojjau7xRESsqGr3UkRMHJxXpDIz+GizRcRFwL8C/wzsAkwFvgycPJj92kQLM3M7YAfg\nA8BXI2Lf4tx1wNl1rvmL4hzAD4Fbqbz+nYH3A8v6tcf9KCJaBrsPktRIGvC98sNU3isPqDo/ujh/\nEJX3xB9ExLtq7vFnmbld1dfCAem5VMXgo80SETsClwHnZ+ZNmflyZq7JzB9m5oeKNodHxG+LkZBF\nEfFvETG8OLf2U6KWqnv+MiL+qmr/PRHxUPEJ04MRcWhVFw6OiPsjYmlEfDciRhbX/DEi/qzqHsMi\n4tmIOHhDrycrZgJLgFcVh78JvCYidqu63/7F+e9ExE7AHsBXM3N18fWbzOxxOltxjy8UI0SPR8QJ\n9V5/ROwdEb8qXt+zEfHd4vgdRfPfF5+YvaPqz2pORCyJiBnVn6RFxJuKT+iWRsSXi/uu/Tnviojf\nRMQVEbEE+HhE7BURv4iI54qf/a2IGF11vyci4kPFn//LEXFtROxSfIL5YkT8d0SM2dCfgSSVQYO+\nV94MPA8cUOf805n5ReDjwGcjwt8ztVXxP0htrlcDI4EfbKBNB5VRlJ2K9scA7+vNzSPiNCr/4zyb\nyidMJwHPVTV5OzCdSvB4FfCu4vg3gLOq2p0ILMrM+zby85oi4qSir3MAMnMBcDuVEZ61zgZmZuaz\nRX/mAP8VEadExC69eGlHAI8UP+dzwLUREXXafQL4OTAGmAx8qejT0cX5g4pPzL4bEW8EPk3lz2RX\n4Eng+uJ17QTcCFwCjCt+9p/W6dNcKiNWnwKiuN9EYH9gCpW/i2pvA44D9gH+DPgJ8JHidTVRGfmS\npLJrxPfKtwKjgT9soOlNVN5T9t1AG2nAGXy0ucYBz2Zme08NMvOezLwzM9sz8wngK8Drenn/vwI+\nl5l3F58wzcnMJ6vOX5mZCzNzCZXpZms/pfov4MSI2KHY/wsqIzc9mRgRLwArqLwxXZSZv6s6f11x\nD4pPrt5ZHCMzE3gD8ATwL8CiiLgjIqZt4Oc9mZlfzcyO4j67Upn6UGsNsBswMTNXbmQU6Z3A1zPz\n3sxcRSXkvDoidqfyZvZA8UljO3Al8HTN9Qsz80vF39OK4s/61sxclZltwOWs//f2pcx8JjOfAv4H\nuCszf1f8/B8Ah2ygv5JUFo32Xvks8DHgLzLzkQ20XzuNbWzVsZuLUa0XIuLmXr4+qU8ZfLS5ngN2\nig2sCYmIfaKyyPHpiFhGZX7zTr28/xTgsQ2cr/7lfTmwHUAxZ/g3wNuK6VknAN/awH0WZuZoKp+U\nXQm8seb8TcCuEXEk8HpgW+DHa09m5oLMvCAz96ISVF6m8knaRvudmcuLze3qtPt7KiMv/xcRD0Sd\nIgtVJlIZ5Vl735eo/P1MKs7NrzqXwIKa6+dX70TEzhFxfUQ8Vfy9/Rfr/709U7W9os5+vdckSWXT\nUO+VmTk2Mw/OzOs30q9JxfclVcdOKe4xOjNP2cj1Ur8w+Ghz/RZYCWzof15XAw8D0zJzBypTodZO\n63q5+L5tVfsJVdvzgb02s2/XURnCPw34bTEqsUHFSMWHgQMj4pSq48upTBU7m8onYtdn5uoe7jEf\nuAp45Wb2u/peT2fmezJzIvD/gC9Hz5XcFlIJXQBExCgqnzI+BSyiMlVu7bmo3l/742r2P10ce1Xx\n93YW6/7eJEm911DvlZvgrcBiKtOrpa2GwUebJTOXApcCVxXrW7YtFkeeEBGfK5ptT6XC2UsRsR/w\n11XXt1H5xfysiGguRjSq/+f9NeCDEXFYVOwdVUUGNuJm4FDgQjY8+lL7mlZTmbJ2ac2p64B3UFnX\n0lVuOiLGRMQ/FX1rKtbT/CVwZ29/Zk8i4rSIWBtQnqcSRDqK/WeAPauafxt4d0QcHBEjqHxaeFcx\nZeLHFGGu+MTxfLq/adazPfAS8EJETAI+tKWvR5LKqBHfKzekKHRzAZXpcJdkZmdf3FfqKwYfbbbM\nvBy4CPgo0Eblk6cLqPzPFOCDwJnAi8BXge/W3OI9VH6pfg54BfC/Vff+HpWF9t8urr+Z7nOFN9Sv\nFcD3qSzmvGkTX9bXganV1W6AO4ClwFOZeXfV8dXA7sB/U3nT+iOwinWLR7fEnwB3RcRLwAzgwsx8\nvDj3ceC6Yp702zPzNuAfqbzmRVTeFE8HKIownEalkMJzVKrwzCr62ZN/ovJmuJRKcNrUP0NJUqFB\n3ytrvRARL1MpeHAicFpmfn0L7yn1uahM+ZcaS0RcCuyTmWdttHGJFAUaFgDvzMzbB7s/kqTB43ul\nysYRHzWciBgLnAtcM9h92RpExPERMbqYBrd27vgWT8eTJA1dvleqjAw+aigR8R4q0wh+kpl3bKx9\nSbyaStWfZ6k8c+eUYoqDNKRFxPSoPJx3TkRcXOf8eyPiDxFxX0T8OoonzUfloZAriuP3RcS/D3zv\npcHje6XKyqlukqQhJyKagUepPEh3AXA3cEZmPljVZofMXFZsnwS8LzOnF8+4+lFmbnEFRknS0NGr\nEZ+I+EDxLJE/RsR3ImJkROwREXdFxOyI+G5EDC/ajij25xTnd+/PFyBJKqXDgTmZObeoyHg9cHJ1\ng7WhpzCK9Uu3S5JKZKPBpyhn+36gtfh0rJlKxajPAldk5jQq5XbPLS45F3g+M/cGrijaSZLUlybR\n/eG7C1j30MQuEXF+RDxGpbLh+6tO7RERv4uIX0XEa/u3q5KkrUGPTxKu026biFhD5SFai6g84f7M\n4vx1VErsXk3lE7ePF8dvBP4tIiI3MKdup512yt13331T+y5J6kP33HPPs5k5frD70Uv1Hqq73vtM\nZl5F5RkqZ1IpJ3wOlfewqZn5XEQcBtwcEa+oGSEiIs4DzgMYNWrUYfvtt19fvwZJ0ibY0vepjQaf\nzHwqIr4AzANWAD8H7gFeyMz2oln1J21dn8JlZntELKXyFPlnq+9b/YYydepUZs2atbmvQZLUByLi\nycHuwyZYAEyp2p8MLNxA++upfDhHZq6ieJZVZt5TjAjtQ+UZV10y8xqKiletra3p+5QkDa4tfZ/q\nzVS3MVRGcfYAJlKZJ31CnaZrP2nr7adw12Rma2a2jh8/VD5glCRtJe4GphXrTYdTmYI9o7pBREyr\n2n0zMLs4Pr4ojkBE7AlMA+YOSK8lSYOmN1PdjgUez8w2gIi4CfhTYHREtBSjPtWftK39FG5BRLQA\nOwJL+rznkqTSKmYUXAD8jMra069n5gMRcRkwKzNnABdExLHAGiprUc8pLj8auCwi2oEO4L2Z6fuU\nJDW43gSfecCREbEtlalux1CZDnA7cCqV6QPnALcU7WcU+78tzv9iQ+t7JEnaHJk5E5hZc+zSqu0L\ne7ju+8D3+7d3kqStzUanumXmXVSKFNwL/KG45hrgw8BFETGHyhqea4tLrgXGFccvAtZ7qJwkSZIk\nDaReVXXLzI8BH6s5PJfKcxRq264ETtvyrkmSJElS3+jVA0wlSZIkaSgz+EiSJElqeAYfSZIkSQ3P\n4CNJkiSp4fWquIEk9daKFSvo6OgY7G6UVkSwzTbb0NTk51qSJFUz+EjqE48//jizH32QzjWraWlp\nHuzulFYmrOnoZOLkqRx44EE0N/t3IUkSGHwk9YF58+Yx55H7OeygA9hxxx0Guzult2rVKh58aDb3\n3jOLPzn8iMHujiRJWwXnQkjaYk/MncMr99/b0LOVGDFiBAe96gDannmKVatWDXZ3JEnaKjjiI9Xx\n8MQpm33tfgvn92FPhoYXXniOceP2H+xuqEpTUxM77jCKpUuXsvPOOw92dyRJGnSO+EjaIpkJiYvp\nt0LNTU10dnYOdjckSdoq+JuKpNIYNmoCcx57fFD7cNmnPs/Zf3n+oPZBkqQyMvhI6hd779/KpN1f\nycsvv9x17Nr//BbHTH9rr64/ZvpbufY/v7XZP39Lr5ckSY3F4COp37R3dPClL39tsLshSZJk8JHU\nf/7uwvdx+Rev5oUXltY9/7933s2Rrz2ecbtO48jXHs//3nk3AP/48U/z69/cxYUXfYTRO+/J+y+6\nBICHH5nN9Le8nZ0n78crDj6K733/lrr37el6gNtuv4P9X/Vqxk/al7/5wMWVNUqF/7ju2xx46GsZ\nP2lfTjzpdJ6cV79QxZtPPoOr/v3abscOPeKN/OCWHwPwgQ9+lD32OZSxE/bm8KPexK9/c2fd+/zq\njt+w+7RDuh3be/9WbvvFHQB0dnbyuS98iX1feQS7TNmfM/7iPSxZ8jwAK1eu5Oy/PJ9dpuzPThP3\n4cjXHs8zz7TV/TmSJMngI6kfHXboQbzutX/K5V+8er1zS5Y8z8lvO4sL/vqveGb+Q/zt37yXk992\nFs89t4RPfPwSXnPUEXzx8n/mhcVzufLyT/Pyyy9zwp+9g9Pf/lYWPvFHvvkfV/M3H7iEBx58eL17\n17t+rZk/uZXf3vFT7rnzNm68aQY/v/V2AG754U/47Beu5IZvX8uiJx/gNUcdwVnv+uu6r+v0t7+V\n737v5q79Bx96hHnzF3Di9GMBaD3sYGb99r9ZvOBhznj7Wzn9rPewcuXKTf7z+9KXv8YtP/oJt/3s\nB8ybcx+jR4/m/R+ohLhvfOsGli1bxuOP3MMz8x/iqis/xzbbjNzknyFJUlkYfCT1q4999ENc9e/X\n0tb2bLfjM3/63+y91x6cdeZptLS0cPrb38q+++zNj2b+vO59fvyTW9lttym86+wzaGlp4dBDXsVb\nT34zN938o03qz4f+7m8YPXpHpk6ZzOuPPorf3/8AAF+99hv8/Qffz/777UNLSwsXf+hCfn//A3VH\nfU456UR+f/8fu85957s3ccpJJzJixAgA3nnGqYwbN5aWlhY+cOFfs2r1ah559LFN6ifA177+TT7x\nsUuYPGkiI0aM4NKPfJDv3/wj2tvbGdbSwnNLnmfOY0/Q3NzMYYccxA47bL/JP0OSpLIw+EjqV698\nxf6cOP04PvcvX+p2fNGip9lt6uRux3abOpmFi56ue5958xbwf3ffy04T9+n6+s53v8/TzyzepP5M\n2GXdM2222WYbXiqKL8ybv4CLPvTRrnvvPHk/MpOnFq7fn+23344Tpx/LDcWozw033syZ73hb1/kr\nvng1Bx76WsbtOo2dJu7D0qXLeO65JZvUT4An5y3g1DPe3dWnAw97Lc3NzTyzuI2zzjyNNx37es56\n1/9j6l4HcfE/XMaaNWs2+WdIklQWPsBUUr/72Ec/xOFHHccH3v/ermO77jqBJ4s1MWvNm/8Ubzr2\nDQBERLdzkydP4ujXvJqf/uiGXv3M2us3ZvKkSVz8ob/lzNPftvHGwDtOeyuf+Od/4TWveTUrVqzk\n9a87CoBf/+ZOPn/5v/GzH9/IKw7Yl6amJsZP2rfbWqK1th21LctXrOja7+jooO3Z57r2p0yeyDVX\nX8FRrz68bh/+8SMf5B8/8kGeeHIeJ/35O9lnn735y3PO3JSXLUlSaTjiI6nf7b3XHpz2tpP5t6vX\nFQQ44fhjmD1nLt/57k20t7dzw40389DDj/LmE44DYJedx/P44092tX/zCccxe85j/Ne3v8eaNWtY\ns2YNd9/zOx56+NG6P7P2+o0576/O5nNfuLJrzdDSpcu48aYZPbY/4fhjmDd/Af/0ic9x2ttO7nqA\n64svvkRLSwvjdxpHe3s7n/z0v7Bs2Yt177HP3nuxcuUqZv70VtasWcM/f/YKVq1a3XX+PeeezaX/\n9OmuKXVtbc8y40c/BeCXv/o1f/jjQ3R0dLDD9tvT0jKMZh8iK0lSj3yXlDQgPnrJRbz88vKu/XHj\nxnLzjd/kiiv/nV2m7M8XrvgyN9/4TXbaaRwAF7zvPdx0848YP2lf/vaD/8D222/HzBnf5YYbb2bq\n3gcxec9X8ZGPfrJbUKhWe/3GnHLSiXzwogt45znvZeyEvTn4T17PT3/+ix7bjxgxglNOOpHbbr+D\n09+x7tlEbzruDRz/pjdywMF/yl77tTJyxAimTJ5Y9x477rgDX7ri0/y/9/0du+19MKO23ZbJk3bt\nOv/+89/DW048nhNPOp0xu+zFa97wZv7v7nsBePqZNk4/668YO2FvDjz0tRz92lfzzjNO3ejrlCSp\nrKLe9IuB1tramrNmzRrsbkhdHp44ZbOv3W9h/RLIjSoz+eHNN3Li8UcPdldUY9Y997P7tAOZMGFC\nr9pHxD2Z2drP3RqSfJ+SpMG3pe9TjvhIkiRJangGH0mSJEkNz+AjSZLUwNZ0dPLgwqU8+9Kqwe6K\nNKgsZy1JktTAbrx7Po8+Xaku+b5j9mb89iMHuUfS4HDER1K/OKj1aH51x2/65F6/uuM37D7tkD65\n16aaN38Bo3fek46ODgCOmf5Wrv3Pbw1KXyRpc6wNPQCzn3lpEHsiDS6Dj6R+8ftZd/C6o4/arGuH\njZrAnMce7+MeVVz3zet53bEn9br91CmTeWHxXJqbm/ulP5I0oLaCar7SYDH4SJIkSWp4Bh9J/WLv\n/Vu57Rd3AHDZpz7PGX/xHt71VxcwZpe9OKj1aGbde1/d697wplMAOOzINzJ65z254cabu85d8cWr\nmbjbK5iy56v4z298p+v4qlWr+PtLPs6e+x7GpN1fyfve//esWLFivXs/9PCjnH/hh7nzrlmM3nlP\ndpq4DwAzf3orra8+lrET9maPfQ7lsk99vuuaJ56cx7BRE2hvb1/vfnMee5w3Hn8K43adxoSpB3Dm\n2edtxp+UJPWf2uc1Ot6jMjP4SBoQP/zxz3nHqafw7MJHecuJx3PhRR+p2+72n1eCzj13/oIXFs/l\n7adWgtDTzyxm6bJlPDnnPq758uW8/6KP8PzzLwBwyUc/yew5c5n12//m4T/cycKFi/jkpy9f7977\n77cPV33xsxx5RCsvLJ7LswsfBWDUttvyH1/9Es8ufJRbvv9ffOWr13HLD3+y0df08cs+y3HHvJ62\npx7hiUfv5X3vPXez/mwkqb90dHaPOp0mH5WYwUfSgDjq1YdzwvRjaW5u5p1nnMr9f3hwk64fNmwY\nH73k7xg2bBgnTD+W7UZtyyOzHyMzufY//4svfPYyxo4dw/bbb8eHP3Rht5GijXnd0Udx4Cv3p6mp\niVcdeADvOO2t3PE/v93odS3DWnhy3gIWLnqakSNH8po/PWKTXpMk9bfaJT2dJh+V2EaDT0TsGxH3\nVX0ti4i/jYixEXFrRMwuvo8p2kdEXBkRcyLi/og4tP9fhqSt3S677Ny1ve2227By5cq608d6Mm7s\nGFpa1lXg33bbbXj5pZdpa3uW5ctXcMRr3sROE/dhp4n78JZTzqDt2ed6fe+77r6XY0/4c3bd7QDG\n7TqNa679Bs8+t2Sj133mk5eSmfzp0SdwUOvR/Md13+71z5SkgdBRk3xq96Uy2WjwycxHMvPgzDwY\nOAxYDvwAuBi4LTOnAbcV+wAnANOKr/OAq/uj45IEsNNO49hmm234/axf8ezCR3l24aM8t2g2Lyye\nW7d9RKx37Ox3/zVvOfF4Hn/kXp5bNJvzzj17vXnx9UyYsDNfuepfmPfY77nqys/zNx+4pN+q0UnS\n5uisXeNj8FGJbepUt2OAxzLzSeBk4Lri+HXAKcX2ycA3suJOYHRE7NonvZVUCrvsPJ7HH3+yV22b\nmpo4913v5IMfvpTFi9sAeGrhIn5+6+093vuppxayevXqrmMvvvgSY8eMZuTIkfzfrHu5/oabevWz\nb7xpBgueWgjAmNE7EhE0NzuDWNLWo3ZqmzPdVGab+g59OrC2lNIumbkIoPi+dh7LJGB+1TULimOS\n1Cv/+A8f5C/Pez87TdyH733/lo22//QnP8pee+7Ba97wZsZO2Jvpb3k7j8x+rG7bN7z+NRyw/75M\n3vNVTJh6AABf+tfP8E+f/BxjdtmLT336ck59W++e8zPrnvs46nUnMnrnPfnzt5/D5Z//BHvsvlvv\nX6gk9bPaoFNb7EAqk+jtkGdEDAcWAq/IzGci4oXMHF11/vnMHBMRPwY+nZm/Lo7fBvx9Zt5Tc7/z\nqEyFY+rUqYc9+WTvPt2VBsLDE6ds9rX7LZy/8UYNJDP54c03cuLxRw92V1Rj1j33s/u0A5kwYUKv\n2kfEPZnZ2s/dGpJaW1tz1qxZg90NaZO9sHw1X/z5o137h+85lhNeNXEQeyRtvi19n9qUEZ8TgHsz\n85li/5m1U9iK74uL4wuA6t8aJ1MJTN1k5jWZ2ZqZrePHj9/0nkvaajhnXIMhIqZHxCNFMZ2L65x/\nb0T8oSjM8+uIOKDq3CXFdY9ExPED23Np4NSO8DjiozLblOBzBuumuQHMAM4pts8Bbqk6fnZR3e1I\nYOnaKXGSGk9E0NTcTEdHx2B3RTXaOzppbm4e7G70i4hoBq6i8qHcAcAZ1cGm8O3MPLAozvM54PLi\n2gOoTN1+BTAd+HJxP6nhrFfO2tyjEutV8ImIbYHjgOoVv58BjouI2cW5zxTHZwJzgTnAV4H39Vlv\nJW2Vxo7diba23pePVv9rb29n6bKXGT169MYbD02HA3Myc25mrgaup1Jcp0tmLqvaHcW6h9afDFyf\nmasy83Eq71eHD0CfpQFXW9Wtdl8qk5aNN4HMXA6Mqzn2HJUqb7VtEzi/T3onaUjYa9q+/P6eO4kI\nxo8fR1OTlc0G07JlL/Lgw7OZvNseDBs2bLC701/qFdJZ7wmyEXE+cBEwHHhj1bV31ly7XhGemrWo\nfdJpaaDVVnVzarLKrFfBR5I2ZMKECXDYkcx59GHuue8hmprWf1aOBkYmjNx2FFOm7sG+++472N3p\nT/X+I1vvN7rMvAq4KiLOBD5KZWp2b6+9BrgGKsUNtqi30iBZ7wGmnYPUEWkrYPCR1CcmTJjAhAkT\nyEzX+wyipqamsoy49aqQTpXrWfdA7U29Vhqyatf0ONVNZWbwkdSnIoKWFv/Xon53NzAtIvYAnqJS\nrODM6gYRMS0zZxe7bwbWbs8Avh0RlwMTgWnA/w1Ir6UBtv4DTA0+Ki9/O5EkDTmZ2R4RFwA/A5qB\nr2fmAxFxGTArM2cAF0TEscAa4HmKSqRFuxuAB4F24PzMdJhSDak26FjOWmVm8JEkDUmZOZNKJdHq\nY5dWbV+4gWs/BXyq/3onbR1qB3gc8FGZlWIiuCRJUhmtV9zA5KMSM/hIkiQ1KMtZS+sYfCRJkhrU\n+g8wHaSOSFsBg48kSVKDqs05tSNAUpkYfCRJkhqUM9ukdQw+kiRJDap2TY/P8VGZGXwkSZIaVO3M\nNme6qcwMPpIkSQ0qsaqbtJbBR5IkqUH5AFNpHYOPJElSg1ov+AxON6StgsFHkiSpQa3/HB+jj8rL\n4CNJktSgnOomrWPwkSRJalC1xQwsbqAyM/hIkiQ1qNqYYzlrlZnBR5IkqUHVrumpLW8tlYnBR5Ik\nqUG5xkdax+AjSZLUoNZf4zNIHZG2AgYfSZKkBrX+Gh+Tj8rL4CNJktSgaosZGHtUZgYfSZKkBmU5\na2kdg48kSVKDsriBtI7BR5IkqUHVlq92jY/KzOAjSZLUoNZb42PuUYkZfCRJkhrUelPdLG+gEjP4\nSJIkNajaYgadnYPUEWkrYPCRJElqUPXGd6zsprIy+EiSJDWoeiHH3KOy6lXwiYjREXFjRDwcEQ9F\nxKsjYmxE3BoRs4vvY4q2ERFXRsSciLg/Ig7t35cgSZKkemqLG1SOmXxUTr0d8fki8NPM3A84CHgI\nuBi4LTOnAbcV+wAnANOKr/OAq/u0x5IkSdps5h6V1UaDT0TsABwNXAuQmasz8wXgZOC6otl1wCnF\n9snAN7LiTmB0ROza5z2XJEnSBnXWGfKxspvKqjcjPnsCbcB/RMTvIuJrETEK2CUzFwEU33cu2k8C\n5lddv6A41k1EnBcRsyJiVltb2xa9CEmSJK2vXsSpN/1NKoPeBJ8W4FDg6sw8BHiZddPa6ok6x9b7\nJ5aZ12Rma2a2jh8/vledlSRJUu/VW8/jVDeVVW+CzwJgQWbeVezfSCUIPbN2ClvxfXFV+ylV108G\nFvZNdyVJktRrdUKO5axVVhsNPpn5NDA/IvYtDh0DPAjMAM4pjp0D3FJszwDOLqq7HQksXTslTpIk\nSQOn3oiPVd1UVi29bPc3wLciYjgwF3g3ldB0Q0ScC8wDTivazgROBOYAy4u2kiRJGmB1H2A64L2Q\ntg69Cj6ZeR/QWufUMXXaJnD+FvZLkiRJW6jetDaLG6isevscH0mSJA0x9Wa1ucZHZWXwkSQNSREx\nPSIeiYg5EbFetdGIuCgiHoyI+yPitojYrepcR0TcV3zNGNieSwOn3uiOuUdl1ds1PpIkbTUiohm4\nCjiOSjXRuyNiRmY+WNXsd0BrZi6PiL8GPge8ozi3IjMPHtBOS4Og3sNKLW6gsnLER5I0FB0OzMnM\nuZm5GrgeOLm6QWbenpnLi907qTxeQSqVulPdBr4b0lbB4CNJGoomAfOr9hcUx3pyLvCTqv2RETEr\nIu6MiFP6o4PS1qDeeh7X+KisnOomSRqKos6xur/NRcRZVCqTvq7q8NTMXBgRewK/iIg/ZOZjNded\nB5wHMHXq1L7ptTTAXOMjreOJb2CDAAAgAElEQVSIjyRpKFoATKnanwwsrG0UEccC/wCclJmr1h7P\nzIXF97nAL4FDaq/NzGsyszUzW8ePH9+3vZcGSL2M4xoflZXBR5I0FN0NTIuIPYqHa58OdKvOFhGH\nAF+hEnoWVx0fExEjiu2dgKOA6qIIUsOoO63N3KOScqqbJGnIycz2iLgA+BnQDHw9Mx+IiMuAWZk5\nA/g8sB3wvYgAmJeZJwH7A1+JiE4qHwB+pqYanNQw6uUeH2CqsjL4SJKGpMycCcysOXZp1faxPVz3\nv8CB/ds7aetQb1pbvRLXUhk41U2SJKlB1V/jM+DdkLYKBh9JkqQGZTlraR2DjyRJUoOqW9vA3KOS\nMvhIkiQ1qHprfCxnrbIy+EiSJDUoR3ykdQw+kiRJDar+Y3xMPiong48kSVKDqhdyrOqmsjL4SJIk\nNah6Icepbiorg48kSVKDspy1tI7BR5IkqUHVyzhWdVNZGXwkSZIaVL2IY+xRWRl8JEmSGlS90R0H\nfFRWBh9JkqQGVf85PiYflZPBR5IkqUHVCzmWs1ZZGXwkSZIaVP01PiYflZPBR5IkqUG5xkdax+Aj\nSZLUqHyAqdTF4CNJktSg6q3n8Tk+KiuDjyRJUoOqt57H3KOyMvhIkiQ1qHojPpazVlkZfCRJkhqV\n5aylLr0KPhHxRET8ISLui4hZxbGxEXFrRMwuvo8pjkdEXBkRcyLi/og4tD9fgCRJkuqznLW0zqaM\n+LwhMw/OzNZi/2LgtsycBtxW7AOcAEwrvs4Dru6rzkqSJKn36hc3GPh+SFuDLZnqdjJwXbF9HXBK\n1fFvZMWdwOiI2HULfo4kSZI2Q731PK7xUVn1Nvgk8POIuCciziuO7ZKZiwCK7zsXxycB86uuXVAc\nkyRJ0gCql3HMPSqrll62OyozF0bEzsCtEfHwBtpGnWPr/RMrAtR5AFOnTu1lNyRJktRb9dbz+Bwf\nlVWvRnwyc2HxfTHwA+Bw4Jm1U9iK74uL5guAKVWXTwYW1rnnNZnZmpmt48eP3/xXIEmSpLo6O9c/\nZu5RWW00+ETEqIjYfu028Cbgj8AM4Jyi2TnALcX2DODsorrbkcDStVPiJEmSNDB6WstjVTeVVW+m\nuu0C/CAi1rb/dmb+NCLuBm6IiHOBecBpRfuZwInAHGA58O4+77UkSZI2qKeRHau6qaw2Gnwycy5w\nUJ3jzwHH1DmewPl90jtJkiRtlp7W8jjVTWW1JeWsJUmStJXqKeBYzlplZfCRJElqQD2t5TH2qKwM\nPpIkSQ2op7U8lrNWWRl8JEmSGpBrfKTuDD6SJEkNyDU+UncGH0mSpAbU43N8zD0qKYOPJElSA3KN\nj9SdwUeSJKkB9RRvzD0qK4OPJGlIiojpEfFIRMyJiIvrnL8oIh6MiPsj4raI2K3q3DkRMbv4Omdg\ney4NjB6nug1wP6SthcFHkjTkREQzcBVwAnAAcEZEHFDT7HdAa2a+CrgR+Fxx7VjgY8ARwOHAxyJi\nzED1XRooPU1pc6qbysrgI0kaig4H5mTm3MxcDVwPnFzdIDNvz8zlxe6dwORi+3jg1sxckpnPA7cC\n0weo39KA6WmNj7lHZWXwkSQNRZOA+VX7C4pjPTkX+MmmXBsR50XErIiY1dbWtoXdlQaB5aylbgw+\nkqShKOocq/vbXEScBbQCn9+UazPzmsxszczW8ePHb3ZHpcHS81S3Ae6ItJUw+EiShqIFwJSq/cnA\nwtpGEXEs8A/ASZm5alOulYa6noJPWt5AJWXwkSQNRXcD0yJij4gYDpwOzKhuEBGHAF+hEnoWV536\nGfCmiBhTFDV4U3FMaiyu8ZG6aRnsDkiStKkysz0iLqASWJqBr2fmAxFxGTArM2dQmdq2HfC9iACY\nl5knZeaSiPgElfAEcFlmLhmElyH1Kx9gKnVn8JEkDUmZOROYWXPs0qrtYzdw7deBr/df76TB19OU\nNnOPysqpbpIkSQ3IER+pO4OPJElSA7JstdSdwUeSJKkBOeIjdWfwkSRJakA9lrM296ikDD6SJEmN\nqCrgRNVjew0+KiuDjyRJUgOqHvFpalqXfHyAqcrK4CNJktSAquNNS1Xw6Wntj9ToDD6SJEkNqNuI\nT9VcN6u9qawMPpIkSQ2oOt80V091M/eopAw+kiRJDainER/LWausDD6SJEkNyBEfqTuDjyRJUgOq\nXsvT3OSIj2TwkSRJakDV1duafI6PZPCRJElqRNXP62lpaqo6LpWTwUeSJKkBVY/sVOUep7qptHod\nfCKiOSJ+FxE/Kvb3iIi7ImJ2RHw3IoYXx0cU+3OK87v3T9clSZLUE4sbSN1tyojPhcBDVfufBa7I\nzGnA88C5xfFzgeczc2/giqKdJEmSBpAPMJW661XwiYjJwJuBrxX7AbwRuLFoch1wSrF9crFPcf6Y\nor0kSZIGSGePVd0GozfS4OvtiM+/An8PdBb744AXMrO92F8ATCq2JwHzAYrzS4v23UTEeRExKyJm\ntbW1bWb3JUmSVE+PU90sb6CS2mjwiYi3AIsz857qw3WaZi/OrTuQeU1mtmZm6/jx43vVWUmSJPVO\nz1PdBqM30uBr6UWbo4CTIuJEYCSwA5URoNER0VKM6kwGFhbtFwBTgAUR0QLsCCzp855LkiSpRxY3\nkLrb6IhPZl6SmZMzc3fgdOAXmflO4Hbg1KLZOcAtxfaMYp/i/C/SVXSSJEkDKnsY8bGctcpqS57j\n82HgooiYQ2UNz7XF8WuBccXxi4CLt6yLkiRJ2lTV8aal2xofqZx6M9WtS2b+EvhlsT0XOLxOm5XA\naX3QN0mSJG2mzh6nuhl9VE5bMuIjSZKkrVS3qW6u8ZEMPpIkSY2op+IGrvFRWRl8JEmSGpDlrKXu\nDD6SJEkNqPsan3XbPsBUZWXwkSRJakDVa3yam9b9ytfZORi9kQafwUeSJKkBVY/rVC3xqZxzvptK\nyOAjSZLUgKrX+EQEVct8XOejUjL4SJIkNaDqcBN0L3BgZTeVkcFHkiSpAWVNVbduIz6D0B9psBl8\nJEmSGlB1VbemqEx3W8s1Piojg48kSVIDql3jU13goNPcoxIy+EiSJDWi6jU+6434DEJ/pEFm8JEk\nSWpAnbVrfKrOOdVNZWTwkSRJakCdNSM+VnVT2Rl8JElDUkRMj4hHImJORFxc5/zREXFvRLRHxKk1\n5zoi4r7ia8bA9VoaSNVrfLCqm0qvZbA7IEnSpoqIZuAq4DhgAXB3RMzIzAerms0D3gV8sM4tVmTm\nwf3eUWkQda/q1r2ctcUNVEYGH0nSUHQ4MCcz5wJExPXAyUBX8MnMJ4pznYPRQWmwZW3wqV7l41Q3\nlZBT3SRJQ9EkYH7V/oLiWG+NjIhZEXFnRJxSr0FEnFe0mdXW1rYlfZUGRe06HstZq+wMPpKkoSjq\nHNuUX+WmZmYrcCbwrxGx13o3y7wmM1szs3X8+PGb209p0GS3qm7dy1lb3EBlZPCRJA1FC4ApVfuT\ngYW9vTgzFxbf5wK/BA7py85JW4MNrfEx9qiMDD6SpKHobmBaROwREcOB04FeVWeLiDERMaLY3gk4\niqq1QVKj6DaqE/gcH5WewUeSNORkZjtwAfAz4CHghsx8ICIui4iTACLiTyJiAXAa8JWIeKC4fH9g\nVkT8Hrgd+ExNNTip4TRFdHuOj7lHZWRVN0nSkJSZM4GZNccurdq+m8oUuNrr/hc4sN87KA2y6hGf\ngJpy1iYflY8jPpIkSQ2o2xqfJkd8JIOPJElSA6pex9McQXOTVd1UbgYfSZKkBtTRWVvOet25Th/k\noxIy+EiSJDWgzs51201N3Ud8Osw9KiGDjyRJUgPqrJnq1uQDTFVyBh9JkqQG1K2qW01xA6e6qYwM\nPpIkSQ2oW1W32jU+5h6VkMFHkiSpAXUvbmBVN8ngI0mS1ICqy1k3ucZH2njwiYiREfF/EfH7iHgg\nIv6pOL5HRNwVEbMj4rsRMbw4PqLYn1Oc371/X4IkSZJqdVQHH9f4SL0a8VkFvDEzDwIOBqZHxJHA\nZ4ErMnMa8DxwbtH+XOD5zNwbuKJoJ0mSpAFUPajTHNBc9VufuUdltNHgkxUvFbvDiq8E3gjcWBy/\nDjil2D652Kc4f0xE9XI6SZIk9bfaNT7hVDeVXK/W+EREc0TcBywGbgUeA17IzPaiyQJgUrE9CZgP\nUJxfCoyrc8/zImJWRMxqa2vbslchSZKkbrqVs65Z49PhkI9KqFfBJzM7MvNgYDJwOLB/vWbF93qj\nO+v968rMazKzNTNbx48f39v+SpIkqReqs01zU6Wk9bpzBh+VzyZVdcvMF4BfAkcCoyOipTg1GVhY\nbC8ApgAU53cElvRFZyVJktQ7ndlzOWtzj8qoN1XdxkfE6GJ7G+BY4CHgduDUotk5wC3F9oxin+L8\nLzL95yVJkjSQOjewxsepbiqjlo03YVfguohophKUbsjMH0XEg8D1EfFJ4HfAtUX7a4FvRsQcKiM9\np/dDvyVJktSDzOw2qhOBDzBV6W00+GTm/cAhdY7PpbLep/b4SuC0PumdJEmSNln1iE5zU6xX3MAB\nH5XRJq3xkSRJ0tavdrQHaoobmHxUQgYfSZKkBtOR3Ud8gO7lrJ3qphIy+EiSJDWY2mf4ADR1q+pm\n8FH5GHwkSZIaTPVUtuYi+DR3e47PQPdIGnwGH0mSpAZTHWyait/2LGetsjP4SJIkNZjah5dCTTlr\ng49KyOAjSZLUYDo61w8+0W2qm8FH5WPwkSRJajDVuWbtQE+zz/FRyRl8JEmSGkxnvXLW1VPdHPFR\nCRl8JEmSGkz1Gp6uctZhOWuVm8FHkiSpwXR7gGlX8Kk63znQPZIGn8FHkiSpwVSv4Vk70ONUN5Wd\nwUeSJKnBdHuAadPaB5gafFRuBh9JkqQGU+85Pj7AVGVn8JEkSWow3YJP8dtek8/xUckZfCRJkhpM\nZ7fn+NQrZz3QPZIGn8FHkiSpwVSv8VkbfJotZ62SM/hIkiQ1mLoPMO1Wztrgo/Ix+EiSJDWYjZaz\n9jk+KiGDjyRpSIqI6RHxSETMiYiL65w/OiLujYj2iDi15tw5ETG7+Dpn4HotDYyOOlPdmixnrZIz\n+EiShpyIaAauAk4ADgDOiIgDaprNA94FfLvm2rHAx4AjgMOBj0XEmP7uszSQ6pWz9gGmKjuDjyRp\nKDocmJOZczNzNXA9cHJ1g8x8IjPvB2on9RwP3JqZSzLzeeBWYPpAdFoaKNW5pt4aH5f4qIwMPpKk\noWgSML9qf0FxrL+vlYaE6qluXWt8fICpSs7gI0kaiqLOsd7+JterayPivIiYFRGz2traNqlz0mDr\nVtWtzhofy1mrjAw+kqShaAEwpWp/MrCwL6/NzGsyszUzW8ePH7/ZHZUGQ7c1Pk3rr/HpMPiohAw+\nkqSh6G5gWkTsERHDgdOBGb289mfAmyJiTFHU4E3FMalhVE9lW7vGp6XJqW4qN4OPJGnIycx24AIq\ngeUh4IbMfCAiLouIkwAi4k8iYgFwGvCViHiguHYJ8Akq4elu4LLimNQw6gUf1/io7FoGuwOSJG2O\nzJwJzKw5dmnV9t1UprHVu/brwNf7tYPSIGqvCjZrR3pamqPueaksHPGRJElqMBud6taRFjhQ6Rh8\nJEmSGkx7x7rHV60NPhFBU9Vvfk53U9kYfCRJkhpMvRGfynZT3TZSGWw0+ETElIi4PSIeiogHIuLC\n4vjYiLg1ImYX38cUxyMiroyIORFxf0Qc2t8vQpIkSet0dFvjs+7XvWZLWqvEejPi0w78XWbuDxwJ\nnB8RBwAXA7dl5jTgtmIf4ARgWvF1HnB1n/dakiRJPWrvYcSnep1Pe4fBR+Wy0eCTmYsy895i+0Uq\nZUMnAScD1xXNrgNOKbZPBr6RFXcCoyNi1z7vuSRJkurqNuLTXD/4ONVNZbNJa3wiYnfgEOAuYJfM\nXASVcATsXDSbBMyvumxBcUySJEkDoNsan6rn9zQZfFRivQ4+EbEd8H3gbzNz2Yaa1jm23r+siDgv\nImZFxKy2trbedkOSJEkb0d65flU3qJnqVtVGKoNeBZ+IGEYl9HwrM28qDj+zdgpb8X1xcXwBMKXq\n8snAwtp7ZuY1mdmama3jx4/f3P5LkiSpRvepbut+3WuxqptKrDdV3QK4FngoMy+vOjUDOKfYPge4\nper42UV1tyOBpWunxEmSJKn/VRcu6FbOutmpbiqvll60OQr4C+APEXFfcewjwGeAGyLiXGAecFpx\nbiZwIjAHWA68u097LEmSpA3qXs66+jk+1VPdDD4ql40Gn8z8NfXX7QAcU6d9AudvYb8kSZK0mXp6\ngKlV3VRmm1TVTZIkSVu/np7j02zwUYkZfCRJkhpMb6a6GXxUNgYfSZKkBtObqW6Ws1bZGHwkSZIa\nTEe35/g01d2urvwmlYHBR5IkqcG092KqW6dT3VQyBh9JkqQG09NUN8tZq8wMPpIkSQ2kszPJItM0\nNUFTj2t8DD4qF4OPJElSA6kuWtAU3R/F2L2qm8UNVC4GH0mSpAbSbX1Pc/df9SxnrTIz+EiSJDWQ\nzh4KG9TuO9VNZWPwkSRJaiDtPRQ2gO4jQB2Ws1bJGHwkSZIayIaCj1XdVGYtg90BqdE8PHHKZl+7\n38L5fdgTSVIZtbevK1owrGaNT/V+e4fFDVQujvhIkiQ1kDUd1cGn+4hP9f5qg49KxuAjSZLUQKoD\nzfCWnkd81hh8VDIGH0mSpAaypqpoQe1Ut+ogtMbiBioZg48kSVIDWd3LNT7V7aQyMPhIkiQ1kOop\nbMNrg0+LU91UXgYfSZKkBtKtuIFrfKQuBh9JkqQGsuGpbuuquq1xqptKxuAjSZLUQLoXN+heznp4\ns8UNVF4GH0mSpAbS/Tk+3X/Va24KmopDHZ1JR6fhR+Vh8JEkSWog1VPYap/jExFWdlNpGXwkSZIa\nyIZGfGqPWeBAZWLwkSRJaiCruwWfWO+8wUdlZfCRJElqIBsf8Ym6baVGZ/CRJA1JETE9Ih6JiDkR\ncXGd8yMi4rvF+bsiYvfi+O4RsSIi7iu+/n2g+y71p+pqbbVrfGqPWdJaZdIy2B2QJGlTRUQzcBVw\nHLAAuDsiZmTmg1XNzgWez8y9I+J04LPAO4pzj2XmwQPaaWmAbOg5PrXHLGmtMnHER5I0FB0OzMnM\nuZm5GrgeOLmmzcnAdcX2jcAxEbH+ggepwXSb6lZnxKdbVTenuqlEDD6SpKFoEjC/an9Bcaxum8xs\nB5YC44pze0TE7yLiVxHx2no/ICLOi4hZETGrra2tb3sv9aNVa9aFmZF1gs/IYeuOrWrvGJA+SVsD\ng48kaSiqN3JTO2enpzaLgKmZeQhwEfDtiNhhvYaZ12Rma2a2jh8/fos7LA2EzGRlVZgZMax5vTYj\nq46tXO2Ij8rD4CNJGooWAFOq9icDC3tqExEtwI7AksxclZnPAWTmPcBjwD793mNpALR3Jh3Fup3m\n5qi7xqc6DK1c44iPymOjwScivh4RiyPij1XHxkbErRExu/g+pjgeEXFlUUHn/og4tD87L0kqrbuB\naRGxR0QMB04HZtS0mQGcU2yfCvwiMzMixhfFEYiIPYFpwNwB6rfUr1asXhdktqkz2gM1Iz5OdVOJ\n9GbE5z+B6TXHLgZuy8xpwG3FPsAJVN5ApgHnAVf3TTclSVqnWLNzAfAz4CHghsx8ICIui4iTimbX\nAuMiYg6VKW1r36uOBu6PiN9TKXrw3sxcMrCvQOof1Wt2RtRZ3wPdA9HK1QYflcdGy1ln5h1rn31Q\n5WTg9cX2dcAvgQ8Xx7+RmQncGRGjI2LXzFzUVx2WJAkgM2cCM2uOXVq1vRI4rc513we+3+8dlAZB\ndZAZ2cOIz4huxQ1c46Py2Nw1PrusDTPF952L472psgNYLUeSJKmvrawKMtsM3/hUtxWu8VGJ9HVx\ng95U2akctFqOJElSn6ouVjCipYfg0+JUN5XT5gafZyJiV4Di++LieG+q7EiSJKkfdJ/qVv/XPKe6\nqaw2N/hUV8o5B7il6vjZRXW3I4Glru+RJEkaGNVT3Xpa47ONU91UUhstbhAR36FSyGCniFgAfAz4\nDHBDRJwLzGPd4tGZwInAHGA58O5+6LMkSZLqqJ7q1nNxg3XHV63pIDOJqLdaQWosvanqdkYPp46p\n0zaB87e0U5IkSdp0y1e1d233VNyguSkYMayJVWs6yayM+mw7fKO/EkpDXl8XN5AkSdIgeakq+Gw3\noucwM6rq3MtV10iNzOAjSZLUIKpDzKgNBJ/qUPTSSoOPysHgI0mS1CC6jfiMdMRHqmbwkSRJagCZ\nyfJV64ob9HrEx+CjkjD4SJIkNYAVazro6Kw8N37EsCaGNff8a171aJAjPioLg48kSVID6O36ntrz\nrvFRWRh8JEmSGkB1gNlQRbfa8051U1kYfCRJkhrA0hVrura3Hzlsg22rzy+ruk5qZAYfSZKkBvDC\n8tVd22NGbTj4jK46//zyNVSeQS81NoOPJElSA3j+5XUjN6O3Hb7BttsMa2bEsMqvgWvaO1m+umOD\n7aVGYPCRJElqAEtXVI34bCT4RES3cFQ9WiQ1KoOPJElSA6ge8dlx2w1PdQMYXdXmheWu81HjM/hI\nkiQNcWs6OnlxZSW8RMCO2/Qm+Kwb8XnupVX91jdpa2HwkSRJGuLalq1ibX2CsaOG07KBh5euNX77\nEV3bi5cZfNT4DD6SJElD3OIXV3Zt77zDyF5dU92urep6qVEZfCRJkoa4xcuqg8+IDbRcZ+eqEZ9n\nX1pFe0dnn/dL2poYfCRJkoa4p5du+ojPiGHNXQUOOjuh7UWnu6mxGXwkSZKGsM7O5KnnV3TtTxqz\nba+vrW47f8nyPu2XtLUx+EiSJA1hTy9dyer2yjS1HbYZ1quKbmtNHWfwUXkYfCRJkoawJ557uWu7\nOsj0xpSx69o/8ezL5NrScFIDMvhIkiQNYbOffrFre8/x223StbvsMJJthjcD8NLKdha9YHU3NS6D\njyRJ0hC1fFU785asG/GZtsumBZ+mpmDaLtt37T+0aFmf9U3a2hh8JEmShqg/LFhKZ1GFevLYbdhu\nZO/X96y1367rgs/v5z9PZ6fT3dSYDD6SJElDUGdnMuuJJV37B08ds1n32WfCDmw3sgWAF1e08+gz\nL27kCmloMvhIkiQNQQ8uXMqzxbN3Rgxr4hWTdtys+zQ3RbfQ9D+PtFnkQA2pZbA7IPWnhydOGewu\nSJLU51a3d/CLhxZ37R+x5zhGDmve7Pu17jGWOx97lvaOZOELK7h/wQscNGXzRpCkrZUjPpIkSUNI\nZvKzPz7N8y+vBmDk8GaO2GvcFt1zx22GceReO3Xt//T+RV33lxqFwUeSJGmIyEz+59E27n3i+a5j\n0w+cwLbDt3wSz2um7cTobSvFEVau6eRbv32CF1eu2eL7SlsLg48kSdIQsLq9gx/9fiG3V01xe+Xk\nHXnV5NF9cv8Rw5o59U+m0NwUADz30mq+9qu5zK96QKo0lBl8JEmStmLtHZ3c88QSvvyLOd1GevYY\nP4qTD5lERPTZz5o0Zlv+vHUyTcVviMtWrOE/fv04M373FM++tKrPfo40GCxuIG1FtqQYw34L5/dh\nTyRJg2nVmg7mL1nOw4uW8dCiZSxf1dHt/Csn78hJh0yipbnvP8M+YOKODDuiiZvuWcDK1R1kwu+e\nfJ7fPfk8e4wfxX677sC0XbZn9LbD+jR0Sf2tX4JPREwHvgg0A1/LzM/0x8+RJJXXxt5rImIE8A3g\nMOA54B2Z+URx7hLgXKADeH9m/mwAuy4BlefwLF/dzksr21m6Yg1tL66i7cVVLF62kmeWraReRelt\nhjcz/cBdOXDyjv0aOqbtsj3/7/V7MfP+Rcx+et1zfR5ve5nH217mJyxi1IgWdh09kgk7jmTMtsMZ\nM6rytd2Iln4JZNKW6vPgExHNwFXAccAC4O6ImJGZD/b1z9LQYVnprZsjTRpqevlecy7wfGbuHRGn\nA58F3hERBwCnA68AJgL/HRH7ZGb3j9TV8DKTTOjMJAGqtjOL80W7zoSOzqSjM2nv7Kzazq7trnMd\nle3VHZ2sXNPJqjUdrGrvZFV78X1NBy+tauflVe10dvaur9tv08KRe47jsN3HMmILylZvitHbDufM\nI3djbttL3PXYc8x+5sVuYezlVe3MeeYl5jzz0nrXjhjWxKgRLYwa0cI2w5oZ3tLE8JYmRhTfh7c0\nM6KliZbmoKWpieam6Ppqqdqu7DfR1BQE0BQQEURAU8S6fSr7UexL9fTHiM/hwJzM/9/eucbYVZVh\n+HntFWhpZVpJU25TUwkkJlKLkij8gShtlBEVU//YRBNiAonEmFhDQhoTfqDBH0YjqZEIBK1XYmMw\n4N3EWKBgKW3KpaUYxpZWWymUwrQNnz/2mrJnevbMmXb27Zz3SXZnnW/tffY7X78571ln77NWvAgg\naSMwBHjgY0yJeHBp+oxuvGYIWJ/avwS+p+wd0RCwMSJGgD2SdqXn+0fRyQ4eGeH+v+8ZE4vCB/nw\nqR35N45FS0R2Wjyym/Uki/YZ1VHc372O7s438XMUn2/C0034fBHZ75kNWPIDmBQji5H263bAURcS\nnH/uXC5edA6XLTmXiwbOru0N/bLF81i2eB6H3zzO86+8xvOvvM7woaO8dbw4iSPH32bk+DEOHal+\nSuxs8DN2QPSu3IBIgqwna5887uQ/Y/t18nnf2VnKHZPr07j+oj7R+f/yTP6Luzm22xrqtFdXz1/w\ne3Vzgiqqu4yBz1Ig/xHwMPDhEs4zhn76xNpvcE2T6Kd6bNtrRY/Tjdec3CciTkg6DAyk+OZxxy6d\n6GQjJ95mz388s5WZfs6aPYP5c2cyf+4sBubNZtG8OSyaP4clC886owVJy2DBWbO4cnCAKwcHiAgO\nvnGMvf97k4NHRnj16DEOvXGMV48e5+ix7q9klcHo4PedIXYXI2rTF5Qx8Ok0YDul4iTdDNycHo5I\n2l6Clu6Y+vB6EfDfEjYKFRoAAAaySURBVJSUSds0t00vtE+z9U6V3n+tuLRuAVOgG68p2ue0fGr9\nje+vz6dOj7bVn/WWT9s0t00vtE9z2/SekU+VMfAZBvIfAV8A7B2/U0RsADYASNoSEStL0FIKbdML\n7dPcNr3QPs3WWz5t0yxpS90apkA3XjO6z7CkmcAC4FCXx7bap6B9mq23fNqmuW16oX2a26j3TI4v\nY8qNJ4DlkgYlzSb7AummEs5jjDGmf+nGazYBa1P7s8CfIvuCyCZgjaQ5kgaB5cDjFek2xhhTE9N+\nxSfdR30r8AjZFKP3RsSO6T6PMcaY/qXIayR9E9gSEZuAHwEPpMkLDpENjkj7/ZxsIoQTwC2e0c0Y\nY3qfUtbxiYiHgYencMiGMnSUSNv0Qvs0t00vtE+z9ZZP2zS3Sm8nr4mIO3Ltt4CbCo69E7hzCqdr\nVW4SbdNsveXTNs1t0wvt09xXejXZVJXGGGOMMcYY03a8rK4xxhhjjDGm56l84CPpJkk7JL0taeW4\nvm9I2iXpOUkfz8WvT7FdktZVrTmn42eStqbtJUlbU/wSSW/m+u6pS2MeSesl/Tuna3Wur2Ou60bS\ntyU9K2mbpIckLUzxRuYYmlOfRUi6UNKfJe1Mf3tfSfHC+mgC6W/smaRtS4qdJ+n3kl5IP99dt04A\nSZfm8rhV0muSbmtajiXdK+lAfvmAopwq47uprrdJWlGf8mqxT1VL27zKPjX92KfKxz6VyFY0rm4D\nLiObg/svwMpc/HLgaWAOMAjsJvvC6ozUXgbMTvtcXrXuDr/H3cAdqX0JsL1uTR00rge+1iHeMdd1\n603aPgbMTO27gLsanuNG1uc4jUuAFak9H3g+1UDH+mjKBrwELBoX+xawLrXXjdZHk7ZUE68AFzct\nx8A1wIr831JRToHVwO/I1ry5Cnisbv0V5sk+Va3OVnmVfaoUjfap6muiL32q8is+EbEzIp7r0DUE\nbIyIkYjYA+wCPpS2XRHxYkQcAzamfWtDkoDPAT+tU8cZUJTr2omIRyPiRHq4mWx9jSbTuPocT0Ts\ni4inUvt1YCeTrFLfYIaA+1L7PuBTNWop4lpgd0T8q24h44mIv5HNbpanKKdDwP2RsRlYKGlJNUrr\nxT7VGBrpVfap6cc+VTl961NN+o7PUuDl3OPhFCuK18nVwP6IeCEXG5T0T0l/lXR1XcI6cGu6/Hdv\n7nJrE3PaiS+SjeRHaWKO25JLILsVA7gCeCyFOtVHUwjgUUlPSro5xc6PiH2QGSXwntrUFbOGsW82\nm5xjKM5pq2q7IuxT5dFWr7JPTTP2qUroW58qZeAj6Q+StnfYJvqEQR1iMUG8FLrU/nnGFsw+4KKI\nuAL4KvATSeeWpXEKen8AvBf4QNJ49+hhHZ6qsun9usmxpNvJ1td4MIVqy/Ek1JrLqSBpHvAr4LaI\neI3i+mgKH4mIFcAq4BZJ19QtaDKULaR5A/CLFGp6jieiNbV9Otinqn0NbZtX2afqwT5VPv3uU2Wt\n43PdaRw2DFyYe3wBsDe1i+LTzmTaJc0EPg18MHfMCDCS2k9K2g28D9hSls7cubvKtaQfAr9NDyfK\ndel0keO1wCeAayPdxFlnjieh1lx2i6RZZGbyYET8GiAi9uf68/XRCCJib/p5QNJDZLdr7Je0JCL2\npcvZB2oVeSqrgKdGc9v0HCeKctqK2j5d7FPVvoa2zavsU9Vjn6qMvvapJt3qtglYI2mOpEFgOfA4\n8ASwXNJgGqWuSfvWxXXAsxExPBqQtFjSjNReRqb9xZr0nURj73O8ERidIaMo17Uj6Xrg68ANEXE0\nF29kjmlefZ6CJJGtYL8zIr6TixfVR+1IOkfS/NE22ZeJt5Pldm3abS3wm3oUFjLmU/Ym5zhHUU43\nAV9QxlXA4dFbDfoY+1QJtM2r7FPTj32qUvrbp7qdZWG6NrKkDpN9KrIfeCTXdzvZzCPPAaty8dVk\nM3zsBm6vWvM4/T8Gvjwu9hlgB9lMKU8Bn6xTY07XA8AzwLZUHEsmy3XdG9mXV18GtqbtnibnOGlr\nTH0W6Pso2aXfbbm8rp6oPureyGYfejptO0bzCgwAfwReSD/Pq1trTvPZwEFgQS7WqByTmd0+4Hh6\nHf5SUU7JbiH4fqrrZ8jNbtbrm32qcr2t8ir7VCn67FPVaO57n1I60BhjjDHGGGN6libd6maMMcYY\nY4wxpeCBjzHGGGOMMabn8cDHGGOMMcYY0/N44GOMMcYYY4zpeTzwMcYYY4wxxvQ8HvgYY4wxxhhj\neh4PfIwxxhhjjDE9jwc+xhhjjDHGmJ7n/8FqiNLQk168AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "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", "ax1.text(0.69, 0.95, \n", " 'Note the values\\nin the tails', \n", " transform=ax1.transAxes,\n", " fontsize=12,\n", " verticalalignment='top',\n", " bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.3))\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", "Joseph K. Blitzstein and Jessica Hwang, Harvard University and Stanford University, © 2019 by Taylor and Francis Group, LLC" ] } ], "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.6.3" }, "notebook_info": { "author": "Joseph K. Blitzstein, Jessica Hwang", "chapter": "7. Joint Distributions", "creator": "buruzaemon", "section": "7.7", "title": "Introduction to Probability, 1st Edition" } }, "nbformat": 4, "nbformat_minor": 2 }