Skip to content
Permalink
master
Switch branches/tags

Name already in use

A tag already exists with the provided branch name. Many Git commands accept both tag and branch names, so creating this branch may cause unexpected behavior. Are you sure you want to create this branch?
Go to file
 
 
Cannot retrieve contributors at this time
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Problem 3.10.5. 5. Rounded data: it is a common problem for measurements to be observed in rounded form (for a review, see Heitjan, 1989). For a simple example, suppose we weigh an object five times and measure weights, rounded to the nearest pound, of 10, 10, 12, 11, 9. Assume the unrounded measurements are normally distributed with a noninformative prior distribution on the mean µ and variance $\\sigma^2$. \n",
"\n",
"(a) Give the posterior distribution for $(\\mu,\\sigma^2)$ obtained by pretending that the observations are exact unrounded measurements. \n",
"\n",
"(b) Give the correct posterior distribution for $(\\mu,\\sigma^2)$ treating the measurements as rounded. \n",
"\n",
"(c) How do the incorrect and correct posterior distributions differ? Compare means, variances, and contour plots. \n",
"\n",
"(d) Let z = (z1 , . . . , z5 ) be the original, unrounded measurements corresponding to the five observations above. Draw simulations from the posterior distribution of z. Compute the posterior mean of $(z_1-z_2)^2$.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"from scipy.stats import norm,chi2,invgamma,t,uniform\n",
"import matplotlib.pyplot as plt\n",
"from math import log,sqrt\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Ignoring rounding of samples: sample mean= 10.4 sample variance= 1.3\n"
]
}
],
"source": [
"data=np.array([10.0,10.0,12.0,11.0,9.0])\n",
"n=5.0\n",
"xbar=np.mean(data)\n",
"v=data.var(ddof=1)\n",
"print(\"Ignoring rounding of samples: \"\"sample mean=\",xbar,\"sample variance=\",v)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The inverse chi-square distribution which is the conjugate prior for the variance is equivalent to the inverse gamma distribution with the correct parameters: in the notation of scipy, it appears that\n",
"Inv-chi^2(df-1,s^2)=invgamma((df-1)/2,((df-1)/2)s^2). Note the weird thing that the variance of this doesn't exist in the case where df=5."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"2.6356753062363643\n",
"[0.95549603 1.54734748 2.68231072]\n"
]
},
{
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0x7f7b91da8390>]"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3X+cXHV97/HX55yZ2dnZ3WSTzebnBhMhgpEggUgSUKSAFMWC7dUKqBetyuPeC+qlthVue71caGusXpS21IpgS/1RSpHWKFH0olV+h/BDICAaICabn5sf+3Nmdn59+8fMZjebTXaSzM6ZM/N+Ph7zmDlnzpz57GTznu9+z/d8jznnEBGR+uUFXYCIiEwtBb2ISJ1T0IuI1DkFvYhInVPQi4jUOQW9iEidU9CLiNQ5Bb2ISJ1T0IuI1LlIUG88a9Yst2jRoqDeXkQklJ566qk9zrnOo3lNYEG/aNEiNmzYENTbi4iEkpn95mhfo64bEZE6p6AXEalzCnoRkTqnoBcRqXMKehGROqegFxGpcwp6EZE61zBB75yj4ApBlyEiUnUNEfTrXl3Hefecx6pvr+JLT32JXCEXdEkiIlUT2Jmx1bJ+x3puePgGTp11KvNa5vH1F77OUHaIP1v1Z0GXJiJSFXUd9PlCns+t/xxdrV187R1fIxFNMDcxl7tevIvfWvhbnLPgnKBLFBGZcnXddfPz7p+zqXcT1y6/lkQ0AcAnz/gki6Yt4i+e+Auy+WzAFYqITL2ygt7MLjazl81sk5ldP8HzHzazHjN7tnT7WOVLPXr3bbqPjngHF77uwgPrYn6MP1rxR2wd2Mr9r90fYHUiItUxadCbmQ/cBrwTWApcYWZLJ9j0X5xzp5dud1S4zqOWzCZ5ZNsjvOv17yLqRQ967tyuczll5inc+fydGokjInWvnBb9WcAm59yrzrkMcDdw2dSWdfzW71xPtpDl3K5zD3nOzLjqTVexuX8zT+x4IoDqRESqp5yDsQuArWOWu4GVE2z3X8zsXOBXwHXOua0TbFM1D297mOZIM5f/9S5wE3TRWI7WkxJ85Dt/yy8/sbr6BYqIVEk5LXqbYJ0bt/w9YJFz7jTg/wN3Tbgjs6vNbIOZbejp6Tm6So/So9sfZeW8leAO813momT7lhNp28i+9L4prUVEJEjltOi7gYVjlruA7WM3cM7tHbP4NeDzE+3IOXc7cDvAihUrxn9ZVMy+9D62Dmzl9zc/x9/E/+mw221KRfldm8cP/m4ZH+gfPPjJG/umqjwRkaoqp0X/JLDEzBabWQy4HFg7dgMzmzdm8VLgpcqVePRe2PMCAKdmMkfc7qRsljcMZ/hRS6IaZYmIBGLSoHfO5YBrgQcoBvg9zrmNZnaTmV1a2uyTZrbRzH4BfBL48FQVXI4X9ryA5xxLh48c9AC/PZTk6XicXb5fhcpERKqvrDNjnXPrgHXj1n12zOMbgBsqW9qxe37P85yYzZJwk/cOXTSU5G9mtvPjlgQf7B+oQnUiItVVl2fGvrj3Rd5URmseYFEux8nDGR5Q942I1Km6C/p96X3sS+9jSab86Q0uGkrybLyJHr/uPg4RkfoL+ld7XwXgxGz5QX9eMgXAQ83NU1KTiEiQ6i/o+0pBfxQt+iXZLPNyOX6WUNCLSP2pu6Df1LuJRCTBnHy+7NcYcG4yxWPNcYYnOj1MRCTE6i7oX+19lRPbT5zwdN4jeXsyRcrz2BCPT0ldIiJBqb+g73uVxdMXH/XrzkqnaS4U+Jn66UWkztRV0KdyKXpSPZzQdsJRv7bJwcpUmp8lmg+ZyEdEJMzqKui3Dxan4FnYtnCSLSd2birF9miE16J1fYVFEWkwdRX03QPdAHS1dR3T61en0gA81qx+ehGpH/UV9IPHF/RduTwnZLM8pn56Eakj9RX0A90kIglmNM045n2sTqVZH2/ShcNFpG7UXdB3tXVhduyD4c9OpUl5Hr/o+UUFKxMRCU5dHXXsHuw+phE3Y70llcZ3jiu/9U9kenYdWL95zSXHW56ISCDqpkXvnGPb4Dbmt84/rv20Ocdpw8NEWn5docpERIJVN0E/mB0klUsxJzHnuPe1OpXGi28DL1mBykREglU3Qb87uRuAOS3HH/Rnp9KYOSItm457XyIiQauboN+VLPanz07MPu59vWk4g8vH8dV9IyJ1oH6CfqhyQR8BckMnlvrpNSGCiIRbaEfdLLr+/oOWYx0P0TQb3vaXT4F7js3HeXJrfmgJ0WkbsdgeXKbz+HYmIhKgumnRW7SfQq4FXLQi+8sNLQEgknilIvsTEQlK3QS9F+nD5aZVbH8uO5NCZob66UUk9Oom6C3ah8tOr+QeyQ2dRKTlFaBQwf2KiFRX/QR9pJ9CrpJBD/mhkzC/NKZeRCSk6iPoLYcXGapo1w1APnkSgMbTi0io1UXQW6QfgEJFu27A5VvIp+ern15EQq0ugt4rBX2lW/RQ7L7xm39DMqvpEEQknOoi6C0yAIDLtVZ837mhkzAvzzO7n6n4vkVEqqE+gt4fAsDlKx/0+eRiXCHCY9sfq/i+RUSqoT6C/kCLvqXyO3dR8qnX8fiOxyu/bxGRKqiToB+ikEsA/pTsPz90Ei/vf5m9qb1Tsn8RkalUVtCb2cVm9rKZbTKz64+w3XvNzJnZisqVWEZ9/sCUdNuMyA0Vh1mu37l+yt5DRGSqTBr0ZuYDtwHvBJYCV5jZ0gm2awM+CTxR6SInY5GhKTkQO6KQXsC02DT104tIKJXToj8L2OSce9U5lwHuBi6bYLubgb8C0hWsryyePzilQQ8eK+et5LEdj+Gcpi0WkXApJ+gXAFvHLHeX1h1gZsuBhc6571ewtrJZZHBKu24AVs1bxc6hnWwZ2DKl7yMiUmnlBL1NsO5As9bMPOBLwKcn3ZHZ1Wa2wcw29PT0lF/lEXeaw/z01Iy4GWPVvFUAPL5do29EJFzKCfpuYOGY5S5g+5jlNuBU4D/MbDOwClg70QFZ59ztzrkVzrkVnZ2VuZiH+YPFfefbKrK/w1nYtpD5LfN5bIf66UUkXMoJ+ieBJWa22MxiwOXA2pEnnXN9zrlZzrlFzrlFwOPApc65DVNS8TgWKQX9FLfozYzV81ezfsd68oX8lL6XiEglTRr0zrkccC3wAPAScI9zbqOZ3WRml051gZMZadEXclPboodi981AdoAX97445e8lIlIpZV0z1jm3Dlg3bt1nD7PtecdfVvkOtOjzU9uiBzhr3lkAPL7jcZZ1Lpvy9xMRqYTQnxnrHei6mfoW/cz4TE6ZeYr66UUkVEIf9OYP4gpRcLGqvN/qeat5dvezmrZYREKjDoI+WZVumxGr5q0iW8hq2mIRCY06CfpE1d5v+ZzlRL2oZrMUkdAo62BsLZuqoN8cv/LgFTcW75qB5XNn89jTt8P3bzr8Dm7sq3hNIiLHQi36Y7A6leblphh7vdB/fCLSAMKfVH6q6kG/KlWct219c7yq7ysicixCHvSFUou+uarv+sZMhrZ8gccU9CISAuEOei+Nmat6i94HVqbTPNYcR5MWi0itC3XQm58CqHrQQ7H7ZmckwpZI6I9ni0idC3nQF09aCiLoV5f66R9X942I1DgF/TFamMsxP5tTP72I1DwF/bG+N7AqnWZ9PI4mLRaRWlYXQU8AQQ/FfvoB3+PFWHXm2RERORahD3rnrOrDK0esTKUx53gkoe4bEaldoQ96CnGC+jFmFgq8KZPhoeZgvmhERMoR+qAPon9+rHOTKZ5virFf0yGISI0KdTrVQtC/LZnGmfGwRt+ISI0KedBXf56b8ZZmMnTk8jyUUPeNiNSmkAd99ee5Gc8D3ppK8UhznFyglYiITCzkQT8UeIse4G3JFP2+z3NNTUGXIiJyiBAHfR7zh2si6M9OpfGd4yENsxSRGhTaoB+d0Kx614s9nDbnWJ4e5ucaZikiNSjEQT8y/UFthOu5qRS/aoqx0/eDLkVE5CB1EPTBd91AcZgloO4bEak5oQ16DnTd1EaL/sRslvnZnLpvRKTmhDbozSu2oF2hNlrQBrw9meLx5jgps6DLERE5ILxB7xeDnhpp0QOcn0yS9jzNUS8iNSW8QV9jLXqAM9PDtOUL/ERnyYpIDQlv0PspXCECrnau2RoF3p5K8bNEM7mCzpMVkdoQ2qDHS+MKtddyPn8oSa/v88zuZ4IuRUQEgLKaw2Z2MXAr4AN3OOfWjHv+vwHXAHlgELjaOfdihWs9uCY/jcvXTrfNiHNSaWIFxwf/+U6Gd+0+5PnNay4JoCoRaWSTtujNzAduA94JLAWuMLOl4zb7tnNumXPudOCvgFsqXun4urxU6aIjtSXhHKvTaSKtGwEXdDkiImV13ZwFbHLOveqcywB3A5eN3cA51z9msYUqJFyttuih2H3jxXrxmnYEXYqISFldNwuArWOWu4GV4zcys2uAPwRiwPkVqe5IvDQuO2PK3+ZYnJdM4ZwRadtIZnh+0OWISIMrp0U/0dk/h7TYnXO3OedOBD4D/NmEOzK72sw2mNmGnp6eo6t0/L5quEU/s1Agn3odkbaNQZciIlJW0HcDC8csdwHbj7D93cB7JnrCOXe7c26Fc25FZ2dn+VVOoFb76Efk+pfhx3fixQ49ICsiUk3lBP2TwBIzW2xmMeByYO3YDcxsyZjFS4BfV67EQ2XyGczL1cw8NxPJDSwrdt9Mey7oUkSkwU0a9M65HHAt8ADwEnCPc26jmd1kZpeWNrvWzDaa2bMU++mvmrKKgYHMQLG2Gm7Ru9w08slFRNqeD7oUEWlwZY2jd86tA9aNW/fZMY8/VeG6juhA0NdoH/2I3MBpxOd+Fy+2i0JmTtDliEiDCuWZsYPZQaC2W/QAuf5T1X0jIoELZdD3Z0rD9mu4jx7A5dvIJxeXgl4nT4lIMEIZ9GHoox+R6z8Nv6kHr2ln0KWISIMKZdAPZkpdNzXeRw+QG1D3jYgEK5RBP9qir+2uGwCXbyWfPJHotF+g7hsRCUI4gz47gHMGhVjQpZQl23c6XmwfXvOWoEsRkQYUzqDPDJTOig3HtVlzA8twhSjR6U8HXYqINKDauTzTURjMDNb0WbEAm+NXHrT8mWQHD017jMyN0yn775Ab+ypel4g0ntC26F2hKegyjsqlg0MM+B4/0/VkRaTKQhn0/Zn+mm/Rj7cylaYzl+N7rS1BlyIiDSaUQT+YHQzFGPqxIsC7hpI8lGhmvxfKj11EQiqUiTOQGYAQjKEf73cGhsiZ8YOWRNCliEgDCWXQD2bC16IHODmb5eThDGvb1H0jItUTuqAvuEKx6yZkffQj3jM4xMamJn4ZiwZdiog0iNAF/WB2EIcLZYse4HcGh4gVHPe2tQZdiog0iPAFfYjmuZnI9EKBi5JJ7m9tIWnhOOFLRMItdEE/Ms8NIZjn5nDe2z/IoOfxgA7KikgVhDbow9qiBzhjeJjFmay6b0SkKsIb9CHto4fiDD3vHRjkuXgTv4rqoKyITK3wBX12pEUf3q4bKE6JEHWO76hVLyJTLHxBf6CPPrwteoD2QoF3DCVZ29bCkA7KisgUqt3ZK2+cPuHqgfZpMKOdX0U/Tth7Pa7sH2BdawvfbW3hyoHBoMsRkToVuhb9oOfRXCgQ8owH4M3DGZalh/n29DYKQRcjInUrdEE/4Hm0FeonFj/QP8BvolEebg53V5SI1C4FfcAuGkoyO5fjW9Pagi5FROpUCIPeaC3Uz0W2o8D7+wd5NNHMK9HaPWQiIuEVwqCvrxY9FMfUxwpOrXoRmRKhC/rBOgz6mYUC7x4aYm1rC3t1URIRqbDQpUo9tugBPtzXT8aMb05Xq15EKitUQe8oBn1rHQb94myOC5Mp7p7WxoBOoBKRCgpV0A+bkTWjrY4Oxo71sd4+Bj2Pf1FfvYhUUFlBb2YXm9nLZrbJzK6f4Pk/NLMXzew5M3vQzF5X+VKLI24AptVhix5gaSbLOckU35jeRlqtehGpkEmD3sx84DbgncBS4AozWzpus2eAFc6504B7gb+qdKFQ7LYB6rLrZsRH+/rZ5/v8e6uuKysilVFOi/4sYJNz7lXnXAa4G7hs7AbOuZ8655KlxceBrsqWWTQS9PV4MHbEivQwb04P8/X2aWTymaDLEZE6UE7QLwC2jlnuLq07nI8CPzieog5nsAGC3oD/0dvHjkiE+359X9DliEgdKCfoJ+osnvBoqJl9EFgBfOEwz19tZhvMbENPT0/5VZY0QoseYHUqzRnpNLc/dzupXCrockQk5MoJ+m5g4ZjlLmD7+I3M7ELgT4FLnXPDE+3IOXe7c26Fc25FZ2fnURc7cjC2XkfdjDDgE/v76En1cM/L9wRdjoiEXDlB/ySwxMwWm1kMuBxYO3YDM1sOfJViyO+ufJlFjXAwdsSK9DBnzz+bO5+/k6HsUNDliEiITTqLlnMuZ2bXAg8APvB159xGM7sJ2OCcW0uxq6YV+FcrDgvc4py7tNLFDngeEedodvXdoh/x40eW07L4UZZ/+X+T2XvBIc9vXnNJAFWJSNiUNV2ic24dsG7cus+OeXxhheua0MhZsY0ywryQXkh2YCmxjp+T7V2Jy+v6siJy9EJ1Zmy9znNzJJndF4OXJdb546BLEZGQClXQD3peXc1FX45CZjbZ/SuJtq/Ha9oZdDkiEkKhCvoBz+p2+oMjGd5zIRSaaJp9P4cZ2SoiclghC/rG67oBIN/C8J4LiLT+Gr/l5aCrEZGQCV3QN8LQyolk962mMDyLpjn3A7mgyxGREFHQh0aE9K534zf1EOt4KOhiRCREQhP0OSDpeQ3ZRz8iP3QK2f5Tic16EIvuDbocEQmJssbR14Kh0lmxjRT0m+NXHrJuZ6/PZa3zeOv8m3A3fuLI5xTc2DdltYlIeISmRd/fIPPcTGZuPs+1+/t4JNHMjxLNQZcjIiEQmqBvlJkry3FF/wBvHM6wpmMmfV5o/glFJCChSQkF/agIcOOevfT6Hms6ZgRdjojUuPAFfV5BD8Xry368t5/vt7bwoLpwROQIQhf0jXQwdjIf7+3jjcMZbpo1k/3qwhGRwwhNOvSr6+YQUeDPe/bS73ncPGumJkcQkQmFJugHPA/PORINMhd9ud6QzXLN/j5+3JLg31pbgi5HRGpQqIK+tVAIT8FV9JG+flam0nyuYwabotGgyxGRGhOa3CxOaKbW/ER8YE3PHhLO8cezO0hbo1yaRUTKEaKgN/XPH8GsfIHP9exlUyzGmpnFIZeLrr8/4KpEpBaEZgqEfr9Bpyg+Cmen0ny0t48726dz2vAwpIOuSERqQYha9Ar6cly7v4/VqRR/PmsmXnxL0OWISA1Q0NeZCPCF3XuZk8vR3PUNdid3B12SiARMQV+HphcK/PWuPZg/zHX/cR2ZfCbokkQkQKEI+hzFaYp1Vmz5lmSzpLf/Ps/1PMefPvynFJw+O5FGFYqgHzpwVqyGVx6N3MCpXHfmdfxw8w+5ZcMtQZcjIgEJxaib0bno1So9Wh9500fYNbSLu168izktc/jQ0g8FXZKIVFkogl5TFB87M+NP3vIn7E7u5gtPfoHO5k4uXnxx0GWJSBWFoutGQX98fM/nc2/7HMtnL+eGh27gJ1t+EnRJIlJFoQp6HYw9dvFInNsuuI2lHUv59M8+zc+7fx50SSJSJTUZ9ONP3VeLvjJaY6185R1fYUn7Eq776XU8uv3RoEsSkSoIRR99v64udUw2x6+EGw9eNw34mufxB3Nn88kHPs4tu/dwbuowcyXc2DfVJYpIFdRki368Qc/DnKNFc9FXxPRCgTt27ubEbJZPzenkBy2JoEsSkSlUVtCb2cVm9rKZbTKz6yd4/lwze9rMcmb23koXWZyL3oXjWykkZhQK3LljN29OD/OZzg7uaWsNuiQRmSKTZqeZ+cBtwDuBpcAVZrZ03GZbgA8D3650gVDsutGB2MprdY6/39XD21Jpbp41k6+0T9PlCEXqUDl99GcBm5xzrwKY2d3AZcCLIxs45zaXnpuSNO71PaYX8lOx64YXd44v7+rh/86ayd/NaGdzNMpNe/bS5A4+KL55zSUBVikix6OcoF8AbB2z3A2snJpyJtbneUxXi37KRIGb9+xjUTbHrTPb2R6JcOuunqDLEpEKKafbe6Lr0h3TX/hmdrWZbTCzDT095QdJn+/RrhE3U8qAj/X18/929fBSLMqV8+fiNW0LuiwRqYBygr4bWDhmuQvYfixv5py73Tm3wjm3orOzs+zX9aqPvmouSqb4hx27yRkkFn2F6PQngy5JRI5TOUH/JLDEzBabWQy4HFg7tWWNKlA8GKsWffUsy2S4Z9tO8slFxOd/h/i8fyWVSwVdlogco0mD3jmXA64FHgBeAu5xzm00s5vM7FIAM3uLmXUD7wO+amYbK1XggOfhzNRHX2UzCwVSW/+A4Z7zibY/xZX3X8nL+14OuiwROQZlnRnrnFsHrBu37rNjHj9JsUun4vpKZ8W2a9RNADwyey4in1pEb/NaLr//cj6x/BNctfQqfM8PujgRKVPNn4PU6xdLnK6um8Dkh97AfZfex3ld5/Glp77ER3/0UbYN6kCtSFjU/Fw3Iy16dd0Ea/mNjwIXEJk+gw35tfz2v17KcM87yO47Gxht3Wu8vUjtqfmg7/XUoq8dRq7vTIaGTiQ+99+Jz7mf6LRnSe/8XQrpKem5E5EKqPmum35/pI9eQV8rXK6dVPdVpLo/gEX6SSy6jaY5a8FLBl2aiEwgBC16H3NOc9HXHCM3sIzc0BKaOn9IdMZjRKc/wzdfHOb9p7yfqBcNukARKan5oO/zPNoKBTTGo/o2x68sb8M++FUyyhc72vn8k5/nXx75cz69r5e3p1KY5rQXCVzNB/1+31O3TQi8IZvlqzt7eKg5zhdmzuATczs5LT3M4zd/mfzQSYydSUMHbEWqq+b76Pf6Ph15jaEPAwPOTaW5b9sOPrtnL7sjPokT7qT5dV/FT7wSdHkiDavmg36f79GhETehEgXeNzDE/Vu3k97xHrzoPhKv+xqJ1/0dkbYXyOvkN5GqqvmgV4s+vGJAtncVQ6/8Memdl2KRAZq7vsll372Me16+h3TuMNeqFZGKqumgzwK9vs9MtejDzUXJ7j+boVf+iFT3lbRGW7n58Zu56N6L+NJTX2LrwNbJ9yEix6ymD8bu94tjbdSiD69DRu5kwW2ADfEmvjEtyT8+fyf/8PydnJ1K8/sDg5ybTB38S6lROyLHraaDfl/pZCkFfX0x4C3pYd6SHman7/NvbS3c29bKp+Z0MjuX492DQ/zOYJKTstmgSxWpCzXddbNXLfq6Nzef57/39vPA1u3cuquHUzJZ7po+jd/tmsf75s/lro130ZPUZQ1FjkdNt+hHg1599PUuApyfTHF+MsVez+OHrQnub2nhixu+yC1P3cKKOSs4/4TzueCEC5jbMjfockVCpcZb9MXyZqpF31A6CgU+0D/It3fs4nvv+R5Xn3Y1+9L7WLN+De+49x1cef+V3PH8HbzW91rQpYqEQk236Pd5Pk2FAi3umK5FLnVg0fRFXHP6NVxz+jW81vcaD255kAd/8yC3Pn0rtz59KwvbFnLO/HN464K38pa5byERTQRdskjNqemg3xPx6cgXxpw8L41m0fX3j1szD/ggFunlLz/o8ci2R/juK9/l7pfvJupFOWPOGZwz/xzOmncWp8w4RVfCEqHGg36n7zM3nwu6DKlBLtfOFadcwhWnXMFwfpindz3No9sf5eFtD3PLU7cA0BptZfns5ayYu4IVc1bwxo43alZNaUg1HfQ7IhFOGx4OugwJ0BFn0LyxeNcErC7dPg3s8n02xJt4Kj7A3YPNPLTtIQCaI82c3nk6yzqXsWxW8dbR3DG1P4BIDajZoC8AuyI+84bUopejMyef55KhJJcMJfl6+tOYP4CfeI1M4jUeHtzMo9ufwKw4kquQaSefPoEbzr+IZbOWcfLMk2mJtgT8E4hUVs0G/T7fI2fG3JxG3MixO/AXQRboK96SZrwUi/FCU4znmoZ4oXkPX9zw3IHXLMxmOTmT5Q2ZDCf/3l28YcYbWNC6ADMdLZJwqtmg3+EXS5uroZVSYQnnOHN4mDPHdAvu8T1eiDXxy6Yov4rF+FUsyoOJZtxPPwWAyzdx5rw38fr21/P66a9n8fTFLJ6+mHkt8/Cspkcpi9Ru0G+LFkubn1PXjUy9WfkC56VSnJdKHViXNONU9zn8+A68+A7WD+/gqdgvscjotXFdIUIh00lhuJNCppMv/95FdLV10dXWxYymGforQGpCzQb9b0pBvzCroJdgJJzjVa6HNMVbyX7P47VolNeiEV6LRXktOsBriS1smxbhMw89OPr6QoGubI6uky5mQdsCulq7DnwJzG+ZTzwSr/4PJQ2pdoM+EmVOLkdCJ0tJjZlRKDBjeJgzxo0IGzbYGonSHYnQHY2wLeLTHYmwZWALj25/lHT+4Pn325vamZOYw9yWucxJzGFOy5wD93MTc5mdmK0TwKQiajboN0cjLFJrXkKkycFJ2Wxx1s3RHiAWbfkY4DB/EIvt468/dALbB7ezO7mbnUM72ZXcxXM9z7F/eP8h+5wWm0Zncyezmmcxs3kmHfEOZjXPoqO5dB/voKO5g5nxmUS8mv3vLAGr0d8Mx+ZolHcNDQVdiMhxO+hcAAf808Tbpc3Y7fvsivjs9H12RSLsjAyw19/JHt/nBd9jr++T9A49+GvO0V4o0JHPMyNfoL1QYHo+T3uhQPuFNzO9aTrtTe20N7UfeDwtNk1nDjeImgx6i+5lwPc4OZMJuhSRqok7xwm5HCdMMgAhacbeUuiP3jz2lB7v9z1eiUbpjTfR73nkNnxxwv04Z1CI4/IJTp+/gLZYG22xNlpjrbRFi/et0dbR9aXHrbHifUukRV8UIVGTQe83dwNw2rCCXmS8hHMkcnkWlnGOiQMGzej1Pfo8n17fo9fz6PM9ekvLfZ5H79Zu+jyPbeYx4BVvGW/yEUMt0RZao8UvhEQ0QSKSoDnaTCKSOLB84H7M48Nt0+Q3aaTSFKjRoN9Cc6HAiRldYUjkeBjQ5hxtuTwLObpzUjLAgOcxWLr1e3bg8ciXwaDXz4DnMeR5pMxIekaveSSGUyW1AAAFuElEQVQ9I1m6T03Q1XTYep0jXro1OUe8ULxv6lpBs99MU6SJJr+JuB8nHokXH5fumyPNNPlNB9bF/fiB7WNejJgfI+pHDzweuy5ikbr+gikr6M3sYuBWwAfucM6tGfd8E8WexzOBvcD7nXObj6Ug5xyR1hc5Mz1cm99CIg0iRvHaAB2F47vwT4Hi8Yex4T96byQ9b/TeM4ateEuPuU9veYJhM/rNO+j5dGn73HGGtGFEvWjxC8CPjT4e82UQ9aI88Uo/OB/nIuAiYx774CI45/GpC04hYhEi3ugt6kVHl23i9YfbZuzzxzop36RZamY+cBvwDqAbeNLM1jrnXhyz2UeB/c65k8zscuDzwPuPpaBne57Fi/Xy233JyTcWkZrnUepuyjuKsV95OTgk/Md+UWTNyJRuxccUHzO6vvgc45aNLKPrzogc/PrMmH3ngJwZf/+L/5iSn/F4lNNoPgvY5Jx7FcDM7gYuA8YG/WUcmEuQe4G/NTNzrvxB8M459qb3smb9Ggq5Vi4c2lruS0WkwUWAiHPFixQFfOVRBwdC/8C9QY7iXx7Z0uOsjWxTen7M9tnDvDYHfOQYaion6BcAY1O3G1h5uG2cczkz6wM6gD2H2+lL+17izG+cScEVyLs8juJ3QsSLMLzzcloP+oNBRCQcDIgC0ZF2boVP+pyqoJ+o82t85eVsg5ldDVxdWhx++r8+/cLEb/lso11VahZH+FJsMPosRumzGKXPYtTJR/uCcoK+G1g4ZrkL2H6YbbrNLAJMB/aN35Fz7nbgdgAz2+CcW3G0BdcjfRaj9FmM0mcxSp/FKDPbcLSvKWfc05PAEjNbbGYx4HJg7bht1gJXlR6/F/jJ0fTPi4jI1Jm0RV/qc78WeIDi8MqvO+c2mtlNwAbn3FrgTuAbZraJYkv+8qksWkREylfWUHXn3Dpg3bh1nx3zOA287yjf+/aj3L6e6bMYpc9ilD6LUfosRh31Z2HqYRERqW+6BpqISJ0LJOjN7GIze9nMNpnZ9UHUUAvMbKGZ/dTMXjKzjWb2qaBrCpKZ+Wb2jJl9P+hagmZm7WZ2r5n9svT7sTromoJgZteV/m+8YGb/bGYNdVkuM/u6me02sxfGrJtpZj82s1+X7mdMtp+qB/2YKRXeCSwFrjCzpdWuo0bkgE87594IrAKuaeDPAuBTwEtBF1EjbgV+6Jw7BXgzDfi5mNkC4JPACufcqRQHgzTaQI9/BC4et+564EHn3BLgwdLyEQXRoj8wpYJzLgOMTKnQcJxzO5xzT5ceD1D8z7wg2KqCYWZdwCXAHUHXEjQzmwacS3E0G865jHOuN9iqAhMBmkvn5yQ49Byeuuac+zmHnpN0GXBX6fFdwHsm208QQT/RlAoNGW5jmdkiYDnwRLCVBObLwJ8Q+EwlNeH1QA/wD6WurDvMrCXooqrNObcN+CKwBdgB9DnnfhRsVTVhjnNuBxQbi8DsyV4QRNCXNV1CIzGzVuA7wP90zvUHXU+1mdm7gd3OuaeCrqVGRIAzgK8455YDQ5Tx53m9KfU9XwYsBuYDLWb2wWCrCqcggr6cKRUahplFKYb8t5xz9wVdT0DOAS41s80Uu/LON7NvBltSoLqBbufcyF9391IM/kZzIfCac67HOZcF7gPODrimWrDLzOYBlO53T/aCIIK+nCkVGoIVL2lzJ/CSc+6WoOsJinPuBudcl3NuEcXfh5845xq25eac2wlsNbORyasu4OBpwRvFFmCVmSVK/1cuoAEPSk9g7JQzVwHfnewFVb+I0+GmVKh2HTXiHOBDwPNm9mxp3f8qnYksje0TwLdKjaFXObbZaUPNOfeEmd0LPE1xhNozNNgZsmb2z8B5wCwz6wb+D7AGuMfMPkrxy3DSWQl0ZqyISJ3TmbEiInVOQS8iUucU9CIidU5BLyJS5xT0IiJ1TkEvIlLnFPQiInVOQS8iUuf+E783uUBgSy3MAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"sample_v=invgamma((n-1)/2.0,scale=(n-1)/2.0*v)\n",
"samps=sample_v.rvs(10000)\n",
"# as a check, the mean is supposed to be 2s^2\n",
"print(samps.mean())\n",
"print(np.percentile(samps,[25,50,75]))\n",
"plt.figure(1)\n",
"plt.xlim(0,10)\n",
"junk=plt.hist(samps,bins=1000,density=True)\n",
"##### Gelman proposes actually sampling from chisquare and taking reciprocals -- gives the same\n",
"##### result, so we got the normalization right.\n",
"gelman=chi2(4)\n",
"samps_g=4*v/gelman.rvs(10000)\n",
"samps_g.mean()\n",
"np.percentile(samps_g,[25,50,75])\n",
"junk2=plt.hist(samps_g,bins=1000,density=True)\n",
"plt.plot(np.linspace(0,10,1000),sample_v.pdf(np.linspace(0,10,1000)))\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For the mean: for each variance in the distribution, we get a potential mean; the distribution in question is N(xbar,sigma^2/5) so we can choose from N(0,1) and adjust by z=(sqrt(5)(x-xbar)/sigma) so x=sigma/sqrt(5)z+xbar.\n",
"\n",
"If all goes well, we will end up with a very complicated way to construct the t-distribution with 4 degrees of freedom and a scale factor of sqrt(v/5)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"##### "
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XuYVNWZ7/HvS3NRNCBIYxRQ0AcVBBRtiYmamEQT9CQ4ObkMJjkxz8ThmZwhmYm5mTGHM6M5k0Qn48lknCSaOBqPaNRzMjIRA4YxinijReQioi2CNjrQKBcRBJp+zx9rV7ooqrur6apau3b9Ps+zn7WraveuV+z69aq1917b3B0REcmufrELEBGRylLQi4hknIJeRCTjFPQiIhmnoBcRyTgFvYhIxinoRUQyTkEvIpJxCnoRkYzrH+uNR4wY4WPHjo319iIiNenpp5/e4u6NvfmZaEE/duxYmpubY729iEhNMrMNvf0ZDd2IiGScgl5EJOMU9CIiGaegFxHJOAW9iEjGKehFRDJOQS8iknEKehGRjFPQi4hknIJe6stcC4tIHVHQi4hknIJeRCTjFPQiIhmnoJf6tns3PPIIvPhi7EpEKkZBL/VprsG1BqMGwwc+ACefDLNmQXt77MpEyi7afPQiUW0DrgN2AuPGwWuvwc03w7BhcPp1ndt91iMVKFI+6tFLffo/hJCfBKxdCwsWQEMDXH89rI9bmki5Keil/mwAHgcGAlcAAwaE4ZuvfAXc4d645YmUm4Je6s8DSfshIP/Om3/zNzBoECwH/rP6ZYlUioJe6stbwGOAAR8teK2xET77WXDgwapXJlIxCnqpL08B+4HJwMgir3/5y6F9AuioWlUiFVVS0JvZdDNba2YtZnZVkddvMLPlyfKCmW0rf6kiZfB40r4377n8+W+amsIfgG3A89UtTaRSejy90swagBuBi4BWYKmZzXP353LbuPvX8rb/CjC1ArWK9M2mTSG8BwBNRV7Phf05wDxCr14kA0rp0U8DWtx9nbvvBe4CLu1m+8uAO8tRnEhZLVwYxt8nAIO72S73R+BZwlk4IjWulKAfBbya97g1ee4gZnYCMA74jy5en2VmzWbW3NbW1ttaRfpmwYLQTulhu3HAkcAW4IUXKluTSBWUEvTFJu/uqpszE7jX3fcXe9Hdb3L3JndvamxsLLaJSGV0dHQG/ek9bNuPcLAWOn9GpIaVEvStwJi8x6OB17rYdiYatpE0WrYMtmyBEcCxJWyfC/qFCytYlEh1lBL0S4HxZjbOzAYSwnxe4UZmdgowjM7zGkTSY9Gi0E6m+HfUQrnhnYcegn37KlSUSHX0GPTu3g7MBhYAa4C73X21mV1jZjPyNr0MuMtdR68khR59NLQTStx+GKHnv2sXPPNMhYoSqY6SZq909/nA/ILn5hQ8/tvylSVSRh0dsGRJWD+5Fz93MvA64WenTatAYSLVoStjJfvWrIGtW2H06DBGX6pTkjb3bUCkRinoJfsWLw7teeeVNj6fk+v9L1mi8+mlpinoJftyPfLzzuvdz72bMNHZpk3w0ktlL0ukWhT0kn1PPhna972vdz9nwLnnhvXHdTKZ1C4FvWTbtm3Q0hLmmZ80qfc/f/bZoX366fLWJVJFCnrJtmXLQjt6D9wzsPc/35RMfNPcXL6aRKpMQS/ZlgvocYf482edFdpnnoH29rKUJFJtCnrJtlzQn3iIP3/00TB2bLhw6nlNUC+1SUEv2dbXHj10Dt9onF5qlIJesuuNN+Dll+Hww7uYWLtEueEbjdNLjVLQS3bleuBTp0JDH/ajA7JS40qa60akJuXOuDnrLOCxQ9vHXIO3k/Xly8MB2f762EhtUY9esmvFitCe3tOdRnpwBHDCCfDOO+GcfJEao6CX7MoF/ZSe7h1Ygtw+cvsUqSEKesmeuQa3GaxZDWZw2ml936eCXmqYgl6yaSPQAYwfD4MH931/k5N7CyropQbpqJJk06tJW45hm/z9rFwZ2rl58x1/VlMYS7qpRy/Z9ErSlivox48PE6OtXw/bt5dnnyJVUlLQm9l0M1trZi1mdlUX23zGzJ4zs9VmNre8ZYr0UrmDvn//zrH+VavKs0+RKukx6M2sAbgRuBiYCFxmZhMLthkPfAc4191PA/66ArWKlK7cQzf5+9I4vdSYUnr004AWd1/n7nuBu4BLC7b5c+BGd98K4O6by1umSC9sT5bDgCUnHjie3heF4/QiNaKUoB9FZ/8IoJWDZw45GTjZzJaY2RNmNr1cBYr0Wu63dQzlOwo116DtyrCuHr3UmFLOuinWHSo8zaA/MB64ABgNLDazSe6+7YAdmc0CZgEcf/zxvS5WpCT5QV9Ouf2tWBE+AWX6oiBSaaX0d1o58CMzGnityDb3ufs+d38ZWEsI/gO4+03u3uTuTY2NjYdas0j3WpN2dJn3OxQYORLeegveLPO+RSqolKBfCow3s3FmNhCYCcwr2ObfgA8CmNkIwlDOunIWKlKyXDekL1MTd2Vich7CxgrsW6RCegx6d28HZgMLgDXA3e6+2syuMbMZyWYLgDfM7DngIeCb7v5GpYoW6ZJ7ZwhXIugnTAitgl5qSElXxrr7fGB+wXNz8tYduDJZROLZvDlMKzwYOKoC+8/16Fu730wkTTQFgmTHXIPnkvXjqMzB0te/Elr16KWGaAoEyZZKDtvk73cjB597JpJSCnrJlkoeiIVw5s0RwC5gWw/biqSEgl6yJdejP65C+zcO7NWL1AAFvWRLpYdu8vddeDWJSEop6CU73iYMpwwERlTwfdSjlxqjoJfsyPWwj6Wyv9kKeqkxCnrJjmoM2+TvX0EvNUJBL9lR6QOxOcMJUyDvSBaRlFPQS3ZU+tTKHJ15IzVGQS/ZUa2hm/z3UNBLDVDQSzbs2gVbgAbgmCq8n4JeaoiCXrJh7dowJcExVGcGJwW91BAFvWTDc8lsZtUYtsl/HwW91AAFvWTDmjWhrfQZNzkjCBdmbQO2bq3Sm4ocGgW9ZEMu6KvVo+9H5x+V3HuLpJSCXrKh2kM3oKCXmqGgl9q3dy+0tITz24+t4vsq6KVGlBT0ZjbdzNaaWYuZXVXk9S+aWZuZLU+WK8pfqkgX/vcgaG8P4+aDqvi+o5M2921CJKV6PBHNzBqAG4GLCHfKXGpm89y98Lf71+4+uwI1inSvWlfEFlKPXmpEKT36aUCLu69z973AXcCllS1LpBeqeUVsvmMIF2ht2ABvv13lNxcpXSlBPwp4Ne9xK8U/Up80sxVmdq+ZjSlLdSKlqNZkZoX6E8LePVywJZJSpQS9FXmu8LbI/w6MdfcpwO+B24ruyGyWmTWbWXNbW1vvKhXpSqyhm/z31PCNpFgpQd8K5PfQR1NwEzV3f8Pd9yQPbwbOKrYjd7/J3ZvcvamxsfFQ6hU50P79nb+N1e7RQ2fQ64CspFgpQb8UGG9m48xsIDATmJe/gZnln9Q2A1D3RqpjwwbYBxwFHBHh/XVAVmpAj2fduHu7mc0GFhAOPd3i7qvN7Bqg2d3nAV81sxlAO/Am8MUK1izSqdpXxBbS0I3UgJLm+XP3+cD8gufm5K1/B/hOeUsTKUGMK2LzHQeYwYsvhgu3Bg6MVIhI13RlrNS2ak9mVmggMHZsOFbQ0hKpCJHuKeiltsUOeoCJE0OrA7KSUgp6qV3uneE6uvtNK2rChNBqnF5SSkEvteu112DHDjgSGBKxjlyPXkEvKaWgl9qVf8ZNscv6qiXXo9fQjaSUgl5qVy5YY47PQ2fQr10bDsqKpIyCXmpX7FMrc4YOheOOg3fegfXrIxcjcjAFvdSu2BdL5dMBWUkxBb3UrrT06EEHZCXVFPRSm9raYMsWOAwYHrsYOnv0v/1W3DpEilDQS21Kyxk3Obke/cbuNxOJQUEvtSktZ9zk5Hr0GwkXcomkiIJealOaxucBGhvDhVvvEC7kEkkRBb3UpjSdcQNhBktNWSwppaCX2pS2Hj10DiPpCllJGQW91J7t28PwyAAgTXekVI9eUkpBL7Unf2riNP0G6/6xklJp+piIlCaNwzag+8dKapV0K0GRVEnbqZVzkxP5jyZcwJW7mGvEiJhVifxRST16M5tuZmvNrMXMrupmu0+ZmZtZU/lKFCmQ6zHHvNlIMYZ69ZJKPQa9mTUANwIXAxOBy8xsYpHt3gV8FXiy3EWKHCCtQzegA7KSSqX06KcBLe6+zt33AncBlxbZ7lrgOsIlIyKV8fbbsGED9O8PI2MXU4R69JJCpQT9KODVvMetFPSlzGwqMMbdf9vdjsxslpk1m1lzW1tbr4sV4YYjwxQDx7Sn8wiTzqWXFCol6ItNGfXHyTzMrB9wA/D1nnbk7je5e5O7NzU2pukEaKkZuS7HmKhVdC133EA9ekmRUoK+lQM/VqOB/Mk83gVMAv5gZuuBc4B5OiArFZEL+rQdiM1pBAYOhFdfDTcuF0mBUoJ+KTDezMaZ2UBgJjAv96K7b3f3Ee4+1t3HAk8AM9y9uSIVS31Le4++gc6ZLFevjlqKSE6PQe/u7cBsYAGwBrjb3Veb2TVmNqPSBYocoDVpj49aRfemTAntypVx6xBJlHQ4y93nA/MLnpvTxbYX9L0skSK2bIFtwCAgzdciTZ4c2hUr4tYhktAUCFI7cj3k0aT7NzcX9OrRS0qk+eMicqBccKZ1fD4nN3SzYoXuNiWpoKCX2rFqVWjTHvTHHgvDh8O2bbBRN5GV+BT0UjtqpUdvpgOykioKeqkNHR2106MHHZCVVEnjReQiB9uwAXbuhKHAkNjF9GCuwe5kXT16SQH16KU21MqwTU6uTvXoJQUU9FIbai3oc1M0PP887NsXtRQRBb3UhloL+sOAk04KIb92bexqpM4p6KU25F8sVSt0QFZSQkEv6bdrVxgCaWionR496BRLSQ0FvaTfypXh9MoJE2Bg7GJ6Yds1of3dD+LWIXVPQS/p98wzoZ06NW4dvXVC0q6PWYSIgl7Sbq7BPV8O67UW9COBwwkzbr7+euRipJ4p6CX91iftG1fGrKL3+gFjk/XctxKRCBT0km7tdN5V6oTuNkypsUm7bFnMKqTOKegl3V4D9hGGQY6IXMuhGJu0CnqJSEEv6bYhaWuxNw8KekmFkoLezKab2VozazGzq4q8/hdmttLMlpvZo2Y2sfylSl1an7RjI9bQF8cRTgndsAHeeCN2NVKnegx6M2sAbgQuBiYClxUJ8rnuPtndzwCuA/6x7JVKfcr16MfGLKIP+tH5bUQHZCWSUnr004AWd1/n7nuBu4BL8zdw9x15D48AdP806buOjs4efa0O3YCGbyS6UuajH0XneQ8ArcB7Cjcys78EriR8Uf1QWaqT+vbyy2Fe96HAsNjF9MHYpFWPXiIppUdvRZ47qMfu7je6+0nAt4HvFt2R2Swzazaz5ra2tt5VKvVn6dLQjotbRp+NTVr16CWSUoK+lQOnkhpNOOmtK3cBf1LsBXe/yd2b3L2psbGx9CqlPj31VGhPiltGn40GBg6EF16AHTt63Fyk3EoJ+qXAeDMbZ2YDgZnAvPwNzGx83sP/ArxYvhKlbmUl6PsDZ5wR1nPfUkSqqMegd/d2YDawAFgD3O3uq83sGjObkWw228xWm9lywjj95RWrWOrDrwyWLgnrJ8YtpSzekxzWevLJuHVIXSrp5uDuPh+YX/DcnLz1vypzXVLvWoG9wDHAuyLXUg7nnAM/+Qk88UTsSqQO6cpYSad1SVvrwzY5+T1619nHUl0Kekmnl5I2K0F/4okwYgRs3gzr18euRuqMgl7SqSVpsxL0Zhqnl2gU9JI+b70FG4EGavuK2ELnnBNajdNLlSnoJX2am8MleWOorXvE9kQ9eomkpLNuRKpibnIR9m+SxydHq6Qypk0LQzjLlsGePTBoUOyKpE6oRy/p80LSnhK1ivIbOhROPRX27oXly2NXI3VEQS/p0kFn0GepRz/XwpIbp3/88bj1SF1R0Eu6vAK8AzQCwyPXUgn9/zW0ixfHrUPqioJe0iWLvfl8pybtI4/owimpGgW9pEtWx+dzjgGOArZsgTVrYlcjdUJBL+nhwPPJelaD3oAJyfojj8SsROqIgl7SYwuwlXAzyuMi11JJ+cM3IlWgoJf0eC5pTyHbv5kap5cqy/LHSWrN6qQ9LWoVlTeKMMHZxo2wbl2Pm4v0lYJe0sG9s0ef9aA34Pzzw/rDD0ctReqDgl7SYe3aMD4/hHCP1ay74ILQLloUtQypDwp6SYdc4J1G6PFm3Uc+EtoHH4SOjri1SOaVFPRmNt3M1ppZi5ldVeT1K83sOTNbYWaLzCxLk8tKNeQHfT045RQYPRra2mDFitjVSMb1GPRm1gDcCFwMTAQuM7OJBZs9AzS5+xTgXuC6chcqGbZ/P/zhD2G9XoLerLNXv3Bh3Fok80rp0U8DWtx9nbvvBe4CLs3fwN0fcvddycMnqI9RVimX7/eHrVthBDAydjFVMtdg8C1h/cEH49YimVdK0I8CXs173Jo815UvAQ/0pSipM88k7ZSoVVTfJELPfvFi2L07djWSYaUEfbFDY0Wv8jCzzwNNwPVdvD7LzJrNrLmtra30KiXbnk3aqVGrqL53AVOnhpuQaDZLqaBSgr6VcFO3nNHAa4UbmdmFwNXADHffU2xH7n6Tuze5e1NjY+Oh1CtZs3kzrAMGEI4A1ZvcOP38+XHrkEwrJeiXAuPNbJyZDQRmAvPyNzCzqcDPCSG/ufxlSmY98ED4fjgBOCx2MRF8/OOhnTdP0yFIxfQY9O7eDswGFgBrgLvdfbWZXWNmM5LNrgeOBO4xs+VmNq+L3Ykc6P77Q3tG3DKiWXduuEjs5Zdh1arY1UhGlXRzcHefD8wveG5O3vqFZa5L6sG+fbBgQVivt/H5nH6E//aHCb36yZMjFyRZpCtjJZ7Fi2HHjjAlcb2cVlnMWUl7331Ry5DsUtBLPPfcE9qz45YR3SRgILB0Kbx20HkOIn2moJc49u+HuT8L69PilhLdIELYQxi+ESkzBb1U31yD/9kfdhDuoaqZkTq/1fz611HLkGxS0EscTybtNOpjtsqeNAGDBoX56TV8I2WmoJfq6yBcnQHwnpiFpMhg4JJLwrn0d98duxrJGAW9VN8aYDvQCIyNW0qqjP5NaO+8M24dkjkKeqm+3LQu56Fhm3xTCVcHP/WU7iUrZaWgl+p66y14Klk/P2ol6TOIznPq77gjZiWSMQp6qa577oE9wCmEM27kQOcl7U/m6BaDUjYKeqmuW28N7QeiVpFekwjHLtqA3/8+cjGSFQp6qZ6WljDtwSB0kVRX+tH5R/Dmm2NWIhmioJfq+Zd/Ce17gMOjVpJuHyAcpL7vvjBfv0gfKeilOnbuhF/+Mqx/JG4pqTeccAbOvn3wVR3IkL5T0Et1fOVdYabKk4FxsYupAbmJvxcCe/fGrEQyQEEvleceAgvgo1ErqR1TCDft3Irmv5E+U9BL5T3wQLjL8DDCnC7SMwMuTtZ/9CPdZlD6REEvleUO114b1qdT4j3NBIBzgaHAs8/CokWxq5EapqCXylq0CJ54ItxRWDec7J0BdB64nn0R3KH5IuTQlBT0ZjbdzNaaWYuZXVXk9feb2TIzazezT5W/TKlZud78JYR5XKR3PkL4I7kW0L3D5RD1GPRm1gDcSBgxnAhcZmYTCzZ7BfgiMLfcBUoNW7QIHnkEjjoKLopdTI0aDHwsWb8HjdXLISmlRz8NaHH3de6+F7gLuDR/A3df7+4rCDONi4RbBX7962H9m98MgSWH5iJgCPASutWgHJJSgn4U8Gre49bkuV4zs1lm1mxmzW1tbYeyC6kVv/pVOIg4Zgx87Wuxq6lth9HZtfrGN2DPnpjVSA0qJeiLHQE6pO+P7n6Tuze5e1NjY+Oh7EJqwc6dcPXVYf3734fDNd9Bn10IHEeYL+iGG2JXIzWmlKBvBcbkPR5NOCtapLg5c+D116GpCS67LHY12dAfuDxZ/9734J8t3GRdpASlBP1SYLyZjTOzgcBMQAOFUtxTT8GPfwz9+sHPfgZ3NSiQymUScDbw9ttwK4f4vVrqUY9B7+7twGxgAeFun3e7+2ozu8bMZgCY2dlm1gp8Gvi5ma2uZNGSUnv3whVXhBtmXNwBa3UZbNn9N8LMn08Dj0auRWpGSdcpuvt8YH7Bc3Py1pcShnSknn33u7ByJZx0EnzypdjVZNPRwBeAnwO/Aq5uhdH66En3dGWs9N1cg28ZXH89NDSEM24GxS4qw84HzgR2EY6B7NsXuSBJOwW99N0bwM+S9e99D973vpjVZJ8BVxAmiXv0Ufj2tyMXJGmnoJfem5t3xsfOnfAPwFvAZGD0d3TwtRqGAl8FGginW87VRenSNQW9HLr9+8PQwSvAsYRD9vqNqp6Tgc8n65d/Dq7WKZdSnD6Wcmg6gC99CX772zDp1jcIrVTXRYTpn9uBG4D1UauRlFLQS+91ALcAt90WDrpeCbw7bkl1y4DPEW64vhv4e2Dp0qglSfoo6KV39uyBnwIPEeZL/zpwStyS6l4/4MvAWcDbwIUXwpIlcWuSVFHQS3Fzi4z3btsG06fDY4SJtr4BnBahNjnYAMLB2fcQbsL+wfPgjjsO3Cb3/1Tj+HVHQS+lWbYMTh0Gf/gDHAX8D8Il+ZIe/QkHxD8E7AM+//lw6mV7e9y6JDoFvXSvowNuvBHe+17YBBwP/C0wNmpV0pV+wJ8RbgPU0ADXXQfnnw8v6Urleqagl669BnzwgzB7dpjH5sPA3wGaYTrdjHA2zqJFMGpUuGfv6afD74D9kWuTKEqa60bqzNat4aaQC4F9j8DIkaFXv/fTsSuT3vjAB8LcQ5cMhyfehtsJB9G/ELswqTb16KXT9u3wwx+GScnuJ4zzvh+4drNCvhbNNXhgeBi3vxIYSbi7xN8TzsxZvDhqeVI96tELvPwy/NM/wS9+EaY0gOQ28MCJMQuTsjDCqZeTCX/A7ycM6yxaFE6NvRD40TswSDPRZZWCvl5t3w733gu33w4PP9z5/ETg44RQ0Fl42TIQ+ATwEcJ4/e+Atcly9xj4whfgM5+Bs88G0//8LDH3OLepaWpq8ubm5ijvXbfWrYP58+H+++GhhzpvMn3YYfDpT4ebeK85M26NUj27gSXAIsJ8RTknnACf+hR89KNw3nm652/KmNnT7t6ru/oo6LOqowOefx4eeyxcJblkCbz44oHbTADOA6YBVyS/B7qYpv440AI8Dqw8Dl7LuyX0oEFw7rlwwQUwbVro7Q8fHqdOART09ckdNm2CF14IZ1jklmceCz22fEOGhF7a8HvgdMJUtyL5Zu6Hxx+HefPgnuvCJGmFEXHSSXDmmTBhAkycGNqTTw7fDKXiDiXoSxqjN7PpwI8Js1//wt1/UPD6IMKNzc4i3IbiT919fW8KkSI6OuDNN0OQ55bXXw8HT/OX3YWJnhhGmMo2t3x7CwwYoF67dK1fv9CDP/dcOP26cJ+B54AXgR3nwtNPh4uvCi/AMmA4cOq5Yejn+OM725Ejw9LYqGGgSHrs0ZtZA/AC4RKMVmApcJm7P5e3zX8Hprj7X5jZTOAT7v6n3e03sz1693Bx0e7dYdm1q3M9t+zYEQ6GFi7btnW2mzZBW1tpl68PGwbDtsIY4NJ/hMmT4cWL1GOX8msHpjwDK1bAmjWw8AfhwrpNHNzzL+bIIztD/+ijYejQsAwZcnA7ZAgMHhy+KRx++MFLv/o8O7xSPfppQIu7r0ve5C7gUsLf+ZxLCRfGA9wL/LOZmXf3V2TbtnDWR0dHWNw71wsfV+q1/fvD/Tbb20Nb6nrh4z17Dgzycg6HDSYE9lBgymfgmGPgzZ+Eq1Mbgb/cFj4Yf+ylXwmbUchLZfQHnpsa2snJAuEPwBvAlmR5I699C9iRLDt3hmXdur7XMnBgCPzcH4IBA7pf+vfv/rV+/Q5eGhpKe66nbc0OXODg50pdDkEpQT8KeDXvcSthjryi27h7u5ltJ9yvfkuXe33ppXCmRxY1EE5lyy0DCh4fTgjwXJu/HA4cAQxJlgH5O7774Pe6/6iK/CeI9Ep/4Jhk6YoTbmj+FrAd2Ek4jrQ7eX53kcd7k2VfsfW9Ydm+vQL/QdlSStAX+xNS2GUtZRvMbBYwK3m4x2BVCe9fTSPo7o9TqfbT+Qvbd+WpqfzSWJdqKo1qKl0a6+r1HSBKCfpWwuhvzmjCqFyxbVrNrD9h4ODNwh25+03ATQBm1tzbcaZKU02lS2Ndqqk0qql0aazLzHp9cLOUoxlLgfFmNs7MBgIzgXkF28wDLk/WPwX8R7fj8yIiUjU99uiTMffZwALC6PMt7r7azK4Bmt19HvBL4HYzayH05GdWsmgRESldSefRu/t8YH7Bc3Py1t8Bentk9aZebl8Nqql0aaxLNZVGNZUujXX1uqZoV8aKiEh11OcVByIidSRK0JvZUWZ2r5k9b2ZrzOy9MerIq+cUM1uet+wws7+OWVNS19fMbLWZrTKzO80s+mQiZvZXST2rY/4bmdktZrbZzFblPTfczB40sxeTdlgKavp08m/VYWZVP3uji5quTz57K8zsN2ZW1Ysxuqjp2qSe5Wa20MyOi11T3mvfMDM3sxHVrKmruszsb81sY15eXdLTfmL16H8M/M7dTyVMr7UmUh0AuPtadz/D3c8gzNezC/hNzJrMbBTwVaDJ3ScRDoRHPchtZpOAPydcLX068DEzGx+pnFuB6QXPXQUscvfxhMl3r0pBTauA/wo8UuVacm7l4JoeBCa5+xTC9CbfSUFN17v7lOQz+FtgzkE/Vf2aMLMxhOlfXil8rUpupUhdwA25zEqOoXar6kFvZkMIN6j7JYC773X3bdWuoxsfBl5y9w2xCyEcLD88uTZhMAdfv1BtE4An3H2Xu7cDDxNuZVF17v4IB1+rcSlwW7J+G/AnsWty9zXuvraadRS8f7GaFib//wCeIFwbE7umHXkPj6C0mXMqWlPiBuBb1a4np5u6eiVGj/5EoA34VzN7xsx+YWZHRKijKzOBO2MX4e4bgX8g9CReB7a7+8K4VbEKeL+ZHW1mg4FLOPBiutiOcffXAZJ2ZOR6asGfAQ/ELgLAzP6Xmb0KfI4X06A/AAACGElEQVTq9+iL1TMD2Ojuz8aupYjZyVDXLaUMUcYI+v7AmcBP3X0q8DbV/4pdVHJB2AzgnhTUMozQQx0HHAccYWafj1mTu68Bfkj46v874FnCdFZSg8zsasL/vzti1wLg7le7+xhCPbNj1pJ0ZK4mBX9wivgpcBJwBqET+KOefiBG0LcCre7+ZPL4XkLwp8HFwDJ33xS7EMItm1929zZ33wf8P+B9kWvC3X/p7me6+/sJXylf7OlnqmiTmR0LkLSbI9eTWmZ2OfAx4HMpvIp9LvDJyDWcROhkPWtm6wnDW8vM7N1RqwLcfZO773f3DuBmwjGzblU96N39P4FXzSw3Mc+HOXDK45guIwXDNolXgHPMbLCZGeHfKepBawAzG5m0xxMOMqbl3wsOnIrjcuC+iLWkVnIjoW8DM9x9V+x6AAoO6s8Ano9VC4C7r3T3ke4+1t3HEjqoZyb5FVWuM5P4BKVMDunuVV8IXzmagRXAvwHDYtRRUNNgwuzZQ2PXklfT3xF+4VcBtwODUlDTYsIf5meBD0es407C19Z9hA/hlwhTYy8ifMtYBAxPQU2fSNb3EG7PsSAFNbUQphVfniw/S0FN/zf5PV8B/DswKnZNBa+vB0ZUs6Zu/q1uB1Ym/1bzgGN72o+ujBURyThdGSsiknEKehGRjFPQi4hknIJeRCTjFPQiIhmnoBcRyTgFvYhIxinoRUQy7v8DNjuH3fe8s6YAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# it works!\n",
"std_norm=norm(0,1)\n",
"mus=[np.sqrt(y/5)*std_norm.rvs()+xbar for y in samps]\n",
"fig,ax=plt.subplots(1,1)\n",
"x=np.linspace(6,15,1000)\n",
"plt.xlim(6,15)\n",
"junk=ax.plot(x,t.pdf(x,4,loc=xbar,scale=np.sqrt(v/5)),lw=2,color='red')\n",
"junk=ax.hist([x for x in mus if x>6 or x<15],bins=500,density=True,color='orange')\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If we treat the measurements as rounded, then a result of \"10\" means that the true result lies between 9.5 and 10.5.\n",
"To compute $p(\\mu,\\sigma^2|y)$ we need the likelihood $p(y|\\mu,\\sigma^2)$ and this comes from normal hypothesis.\n",
"For example\n",
"$$\n",
"p_{10}(\\mu,\\sigma)=p(10|\\mu,\\sigma)=\\Phi(10.5,\\mu,\\sigma)-\\Phi(9.5,\\mu,\\sigma)\n",
"$$\n",
"where $\\Phi$ is the cumulative normal distribution function. Similarly we have $p11$, $p9$, and $p12$ and\n",
"the likelihood is\n",
"$$\n",
"L(\\mu,\\sigma)=(p_{10}^2p_{11}p_{12}p_{9})\\mathrm{some\\ kind\\ of\\ prior}\n",
"$$\n",
"by independence. The noninformative prior that arises in this case is proportional to $1/\\sigma^2$.\n",
"\n",
"Another way to write this is that\n",
"$$\n",
"p(\\mu,\\sigma^2|y)\\propto \\sigma^{-2}\\int_{9.5}^{10.5}\\int_{9.5}^{10.5}\\int_{10.5}^{11.5}\\int_{11.5}^{12.5}\\int_{8.5}^{9.5} \\sigma^{-5}\\exp\\left(\\frac{-1}{2\\sigma^{2}}\\sum_{i=1}^{5}(y_i-\\mu)^2\\right) dy_1 dy_2\\cdots dy_5\n",
"$$\n"
]
},
{
"cell_type": "code",
<<<<<<< HEAD
"execution_count": 5,
=======
"execution_count": 9,
>>>>>>> f8632e30fe8b89896ed0fbf258653883205028de
"metadata": {},
"outputs": [],
"source": [
"# modified to follow Gelman\n",
"def L1(mu,sigma,dx):\n",
" l=0\n",
" for x in dx:\n",
" try:\n",
" l=l+np.log(norm.cdf(x+.5,loc=mu,scale=sigma)-norm.cdf(x-.5,loc=mu,scale=sigma))\n",
" except:\n",
" l=l-10000\n",
" return l\n",
"\n",
"def L2(mu,sigma,dx):\n",
" l=0\n",
" for x in dx:\n",
" try:\n",
" l=l+np.log(norm.pdf(x,loc=mu,scale=sigma))\n",
" except:\n",
" l=l-10000\n",
" return l"
]
},
{
"cell_type": "code",
<<<<<<< HEAD
"execution_count": 6,
=======
"execution_count": 10,
>>>>>>> f8632e30fe8b89896ed0fbf258653883205028de
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/home/jet08013/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:6: RuntimeWarning: divide by zero encountered in log\n",
" \n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD8CAYAAACSCdTiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xd8U1UbwPHfbZLuvUsXhdKy996yRYYoDhQEBFFcgKCiLFGGggIiCjhQUJaAAiJDhgzZe8+W0kFnukf2ef9IKfCyoVAo5/v5hKTNzb0noXnuuWc8RxFCIEmSJJVONiVdAEmSJOn+kUFekiSpFJNBXpIkqRSTQV6SJKkUk0FekiSpFJNBXpIkqRQrtiCvKIpKUZSDiqKsKq59SpIkSfemOGvyg4CTxbg/SZIk6R4VS5BXFCUIeAr4sTj2J0mSJBUPdTHtZxrwAeByow0URRkADABwcnKqU7FixWI6tCRJ0uNh//79aUIInzt5zT0HeUVROgEpQoj9iqK0vNF2Qojvge8B6tatK/bt23evh5YkSXqsKIpy4U5fUxzNNU2ALoqixACLgFaKovxWDPuVJEmS7tE9B3khxEdCiCAhRFngRWCTEKLnPZdMkiRJumdynLwkSVIpVlwdrwAIITYDm4tzn5IkSdLdkzV5SZKkUkwGeUmSpFJMBnlJkqRSTAZ5SZKkUkwGeUmSpFJMBnlJkqRSTAZ5SZKkUkwGeUmSpFJMBnlJkqRSTAZ5SZKkUkwGeUmSpFJMBnlJkqRSTAZ5SZKkUkwGeUmSpFJMBnlJkqRSTAZ5SZKkUkwGeUmSpFJMBnlJkqRSTAZ5SZKkUkwGeUmSpFJMBnlJkqRSTAZ5SZKkUkwGeUmSpFJMBnlJkqRSTAZ5SZKkUkwGeUmSpFJMBnlJkqRS7J6DvKIo9oqi7FEU5bCiKMcVRRlbHAWTJEmS7p26GPahB1oJIXIVRdEA/ymKskYIsasY9i1JkiTdg3sO8kIIAeQW/qgpvIl73a8kSZJ074qlTV5RFJWiKIeAFGC9EGJ3cexXkiRJujfFEuSFEGYhRE0gCKivKErV/99GUZQBiqLsUxRlX2pqanEcVpIkSbqFYh1dI4TIBDYDHa7z3PdCiLpCiLo+Pj7FeVhJkiTpBopjdI2PoijuhY8dgDbAqXvdryRJknTvimN0TQAwV1EUFdaTxu9CiFXFsF9JkiTpHhXH6JojQK1iKIskSZJUzOSMV0mSpFJMBnlJkqRSTAZ5SZKkUkwGeUmSpFJMBnlJkqRSTAZ5SZKkUkwGeUmSpFJMBnlJkqRSTAZ5SZKkUkwGeUmSpFJMBnlJkqRSTAZ5SZKkUkwGeUmSpFJMBnlJkqRSTAZ5SZKkUkwGeUmSpFJMBnlJkqRSrDiW/5MeY0IIjHoj+gIDBp0Rg856b9QZMeiNGPVGjHoTRr3R+vsrfmcymDDojZiNZowGIyaD9XcmoxmT0YzZaMJkMmM2mTGbLIX3lx9bzJbLN4soeiws1p+FRWCxWECAxWJBCEAIhBBXlP/670tRQLH+U/RYsVGK7m1sbLBR2WBjo1jvVVf/rFKrCm822FzxWK1Wo9aosFGr0NiqUduqi+7Vttbn1LZqbO1t0dipsbXToLHToLFTF95rsLW/fG+92WLrYIutvQYbG1lvk64mg/xjxKA3UpBTQH52AfmX7rPzycsusP4+R0dBTgEFuToKcnXo8qy3glwdunwD+nx94c2ALl+PoTCwixtFyjukuSLQqTSF92oVKs0VQVKjLgyilwOrWqNCsbk60F4KxFcGZuWKgI2iFB33iocAV50MhLCeyBACi+XyvbAIzJdOKIUnF7PZgtFgtj5XeDKymC2YjKbCezNmoxmT0VT02Fh4YivOz9DWwRY7Rzvsneywd7TDzvHyzw7O9tg72uPgbI+Diz0Ozg44utjj4OKAo6sDTq6OOLo6WH92scfR1RE7B1vrZyY9kmSQf8SYjCZy0nPJ1uaQlWa95WhzyE7PJSc9l9yMXHIy88jLzCM3M5+8rHzyMvPIy8rHaDDd1jHsHe2wd7YGAnunwmDhZI+btwt2DrbYORQGjsJgYmtvfWzrYK192jnYFtU6r6qB2qnR2Nta7wt/p9JYa7QqteqxDiRmsxmj3oTZaMJQeOVz6YrHoDNgMpisj/XWq6RLV0aXrpz0BQaMOiP6AutJWF9gKHyspyBPjy5PT7Y2B13hY13hifx2Ti42Khuc3BxxcnPE2d2p8N4RZ3dnnD2ccPF0xtXT2Xrv7Yqbtwtu3i64erlga2/7AD496WZkkH8I6PL1aC+mo72YQUZSJumFt8yULDJTs8hMySY7LZustBzysvJvuB+NrRpnDyecPZxxdnfEzduFMuH+OLlav6COrg44ujjg5OaIg4sDTpdqbIW/d3Sxx97ZHpVK9QDfvQSgUqlQOaoAO5we0DGFEOjy9UVXcXlZ+UVXennZ+RTk6Iqu9PKy8snPsVYacjPySDiXRG5GHrkZeejy9Tc8hr2jHa7eLrj5uOLu44qbjysevm54+Lnj4e+Op787XmU88CrjiZOb42N9or9fZJC/zwpyC0iKSSXlQiopsWmkxGlJS9CSlmAN6tqL6eRnF1zzOpVaZf1i+Fq/HAHlfHHzdsXVy6Xo5ubtgmthjcnF0xl7Rzv5JZFum6IoODjZ4+Bkj6f/3e/HoDcWXV1ma3PISs0mW3v552xtDpmp1orKhRPxZKZkYdAZr9mPrb0GzwAPvMp44BPkhXegFz5BXviEeOMX6o1/WV9cPJ3l3/gdUoqrLfBO1K1bV+zbt++BH/d+EEKQkZxJ3OmLJJxN4uK5RBLPp5B0PoXkmBSy0nKu2t5GZYN3oCfegZ54BXri5e9R9IftVcb62NPfHRdPZ9mJJpVKQgjycwqKrlovVXjSEzNIT8pAezGDtIR00uK16AsMV73W0cUBv7I++JX1ISDMjzLl/Qms4E9QRBl8Q71L/VWooij7hRB17+Q1siZ/m4QQpMZrOX80lphjcVw4GUfcqYvEnoy/qiau1qjwK+uLf5gvFWqF4R/mi19ZX/xCvfEN8cbD373U/yFK0s0oimJtQnR1JCiizA23E0KQk55LSmwaSTEpJMekWitPsdb7QxuPXdVUpLHTEBQRQHDFQEIqBlK2aghh1UIIrOD/WH/nZE3+OoQQpMSmcWrPOc7si+LsgWjOHTxPTnpu0TaeAR6EVLL+MQVXDCQ4sgyBFQLwCfZ6rP+gJOlBuXQVnXA2ifgzF4k7lUDsqQTiTiWQGJ1S1Kls72hHWPUQwmuVI6JOOSo2qEBIpcBH8kr5bmryMshj/WOJPZXAoU3HOLrtBMf+O4X2YgZg7cwsWzWY8FrlKF+zLOWqh1K2ajAuHs4lXGpJkm5EX6An9mQC54/GEnUohrMHo4k6GEN+jvWq29HVgSqNI6nWrDI1WlYmsl44KvXDXzmTQf4O6PL1HFh/hJ0r97J33aGioO4T5EWVphWp0jiSSg0qEFY9FFs7TYmWVZKke2exWIg/k8jpPec4vuM0x7efIuZ4HGAN+rVaV6PhU3Vo0KkOHr5uJVza65NB/hbMZjP71h1mw29b2bVyH7p8PU5ujtRuW526bWtQq3U1/MN8Ze+9JD0mMlOzOLLlBPv/OczedYdIjdNiY6NQ44mqtHqpGS2ea4iDs0NJF7NIiQR5RVGCgXmAP2ABvhdCfH2z1zzoIJ+bmcdfM//hr1nrSI3T4urlQrNnG9K8e0Oqt6iMWiP7nyXpcSeEIOpQDNuW7WLLkp0knE3E0cWBdr1b8szgpwgo51fSRSyxIB8ABAghDiiK4gLsB54WQpy40WseVJA36Az88fVqFn+xnNzMPGq1rkan19vSqEtdNLayCUaSpOsTQnBi5xn+mrWOLYt3IAR0fK0NvUZ3x8PPvcTK9VA01yiKsgKYIYRYf6NtHkSQP3foPJ/3nM6FE/E0eKo2fT57kfCaYff1mJIklT5pF9NZMP4PVv+wASc3RwbPfp1mzzQokbKUeJBXFKUssBWoKoTI/r/nBgADAEJCQupcuHCh2I77/3avPsCn3b/E2cOZYT8NpF6HWvftWJIkPR4unIzni1e+4ez+aPp89iIvj3j2gZfhboJ8sQ0UVRTFGVgGDP7/AA8ghPheCFFXCFHXx8enuA57jbMHohn7zGRCKwcx+9BkGeAlSSoWoZWCmL5jPK17NuOXUYtYNfuGjRUPlWIJ8oqiaLAG+PlCiD+KY593QwjBlNdm4ertwufrRuHu83AOg5Ik6dGk1qj54Je3qdOuBrOHziUjJauki3RL9zysRLGON/wJOCmEmHLvRbp7Z/ZbZ6YOmf06rl4uJVkU6Q4JITDoDNZsh9kFRTnuL+e21xfls78yl72+wGBNu1u0+MilhUdMmI3mq/K6X8r5fmlREQrzxV9JUaz/WHPSX70oyJU57S+lSFZfsfCHrZ3GmkrZXlOUetmaz/1yLndrHvdLudutGUHtnezlsN1HiI2NDQOn9qF/lSFsmr+NZ4d0Kuki3VRxjB1sAvQCjiqKcqjwdx8LIVYXw77vyPmjsQDUeKLKgz60dAWLxUK2NoeMpEwyUrLJTMkiOy2HrLRssrU55GTkWtPUZlpT1eZkWPPf326+e6AoiF5aJelS7nq15vLCI/ZOdqg0alRXrtxUeCta7UlR4FJ8FZfLLyzWRUOuPDkUrU5lNGMsXIDFaDAV/WzUmwpXyDJgKDBYTyS34VK+dpcr0kQ7e1yRo93LBTfvwoykvm54+Lnh7usmh/6WoNBKQXgGeBTFnIfZPf+VCCH+4/LXpESpNdZpySajuYRLUnqZTWZS47UkF6ZOTotPJy1BizYxA+2lbIJJmZhN1/4fKIqCs7sjLp7O1mDm4YRPsDcu7k44ezjh5OaEs7sjjq6Xc99fs3iJox0aO81DX/MVQmAymqyraOVZV9bS5erIL7w6KVqVKyuf3MJFXa48+aXGaa3pe9NzsZgt1+xfURTcvF3wLEzL6+nvYc1uGuSFb2FqXr9QH7lox30ihMBsNBXFnIdZqaoKVKhtHSJ5aNMxQisFlXBpHl1mk5nE8ynEn75I/JmLXIxK5mJUEhfPJZF8IfWaoOPs7mRNlRzoSXClQLwCPPEK8LDWOP3cihaLcPF0fmyStymKgsZWg8ZWg7P73S8DIoQgLyufzNRsslKtV0UZyVmkJ1rXItAmWtPyntkXRWZK9jUrPXmV8SCgnB+B4QGUCfcnMNyf4IqBBFbwx87B7l7f5mPr/NFYstJyCK9drqSLckulKsiHVg4mvFYYf3z9Nx1fay0nPN2GrLRszh08T/ThC0QdieH80VjiTiZc1XTi5OZIYIUAIuuH0/KFxviH+eFf1gefYC+8g7xwcLIvwXdQulmvfpxwdnciqELATbc1GoxoL2ZYF6eJTSPpfAqJ55O5eC6JvWsPkp6UedV+A8r7EVYthPLVy1KuRijla5bFL9Tnob9Kehgs/PxP7B3taPZsyYyXvxOlKsgD9B3XgxFPTWDeJ0voN+Glki7OQ8VsMnPuUAzHt5/ixM7TnNp9juQLqUXPe5XxoFz1UOq0qU5I5WBCKpYhKKKM7MR+RGhsNfiX9cW/rO91ny/I05FwNpH40xeJPZlAzPFYzh+NZcfyvUVXAK5eLkTWD6dKo0gqN46gYoMK8iT+f3avPsDmRdt5ecSzuHm7lnRxbqlUJij7qv9M1s7ZxOilw0psZtrDQAhB/JmL7F17iP3rD3N060kKcnUA+AR7UalhBJF1y1OhTjnK1ygrg/ljSpevJ+ZYLGcPnOfM3nOc2nOuKDujSq2iYoNwarWqRr0ONYmsH/7YNLldT9zpBN5tNAKfYC9m7J74wPs8SnzG6+2630HeoDMw9IlPiDoUw9jlH1Cvfc37dqyHjRCCcwfPs3nxDv77czcXzyUBEFghgNqtq1G9RWWqNKmIT5BXCZdUepjlZuZxYucZjm49wcFNRzm7PxqLReDu40rDTnVo/lwjarep/kjkYC8uSTEpDG05BkOBgW92T7zhFdP9JIP8FbK1OXzQ9lNiTyYwYuFgmjxd/74er6Rla3P4Z+5m1vy0kdiTCajUKmq2qkqjznWp37EWAWEln0FPenRla3PY989hdq3ax+6/D5CfXYC7rxvtXmlBxwFtCAy/eX/Bo+7CiTiGtx+HLk/P5I1jCK9VMnmwZJD/P9naHEZ0msiZvecYMPkVnhn8VKnrVEqJS2PJlytZ8+NG9AUGKjeKoF3vljTr3hBXT9n8IhU/g97I3jUHWT9vMzv/2o+wCBp3rUuPj58lsm75ki5esTuw4QjjXpiCxk7DhDUjKF+jbImVRQb56yjI0zGp9wz++2M3bXo1Z9DMAdg7PvpDxwrydCyc8AdLp6zCYrbQumcznh3ciXLVQ0u6aNJjRJuYwV8z17Hy27XkZOTR4vlGvP5l71LRHCiEYOmUVfw4/DdCKwcxdvkHJX5FLIP8DVgsFuZ/toxfP11C2arBfLxgMGWrBD+w4xe303vPMb7HNBKjk2ndsxl9P+uBX+j9S/omSbeSl53PsimrWDxpOSq1ine/e402PZuXdLHuWrY2hymvzWT78r00faYBw+a8iZOrY0kXSwb5W9m79iCT+nxLQU4Br03qRZc32z9yzTdbft/B572m4xngwYfz3qF688olXSRJKpIYnczkV7/l6NaTdH+vMwMm93rkvmMHNh5lcp8ZZKZk0W/iyzw7pNND8x5KNNXwo6Beh1rMPjSZqs0qMeOdn/i443hSYlNv/cKHxPble5jw0jQqNYxg5oFJMsBLD52Acn5M3jCGLm+2Z+mUv/jxw99Kuki3rSBPx7fvzuHDtp/i4GzP1zvG0/29zg9NgL9bj1VN/hIhBH/N/IcfPvgVG5UNAyb34sn+rbGxeXjPeckXUhlQfSghlQKZtHGMnKAiPdSEEEx/8wdWzV7Ppys+pFHnO6p8PnAHNh5l6oBZJJ1P4el3nqTfxJcfyr47WZO/TYqi0OXN9nx/9Csi6pZj2hvfM+yJT7hwMr6ki3ZD88b+jsloYsSiITLASw89RVEYOK0voZWD+OGDXzGbH86kgVlp2Ux+9Vs+bPspNiobvto8lre+fvWhDPB367EM8pcEhPkxacMYhv44kPNHY3mj5jB+GbUIfYG+pIt2lYI8HZsXbad9nydKZAKGJN0NWzsNLw7vRtzpi5zceaaki3MVi8XCul/+5dVKg9n42zZe/PBpvj/8ZalsAi11uWvulKIodHi1FQ061WH2sLnMH7+MjfO38tb0fjTsVKekiwdA1MHzGHRGGjxVu6SL8tCzWCwU5Oisi48UpvW9tMiIvsBgXVDEYMJkNFtzxRcuHnLJVQuEaFSFC4FYFwGxc7AtWvjDyc0RJzfHx2rG59249B06uu0UVZtWKuHSWEUfucD0t37g+PbTVG4cyeBZAwirGlLSxbpvHvsgf4mHrxvD573Lk6+2ZvpbPzCqy+fU71iLN6f1LfHZfBkp1iVzvQI9S7QcJUmXryc1TktqvNZ6n6BFm5BBRkoWmcmZZKblkKPNITcz/5p0u/eTo4sDLl7OuHm5FC3o4VXGA+8ynvgEe+FTmN/9XtINP8qc3Z1wdHEgIznz1hvfZ7mZefwyahF/zVyHs4czQ38cSLs+LR/qvrjiIIP8/6nRsgqzDk5m+Tdr+XXs7/SvMoRnBnfi5ZHP4ujiUCJlsneytg8W5OhK5PgPihCCtIR0Yo7HEXM8jtiTCcSfSSQxOvmqNLmXuHq54OnvjoefGxF1fIpWUnJyc8LJ1brE3qUFR2ztrTVxjZ0alca6cpSNyqZomb9LxxcWUbhkoHUFKIPehFFnQK8zFC5BaLAu9pFdcHllq/RcMtOySU/K5NyhGDKTM69ZFcrZ3YmAcn4ERfgTHBlIaOUgwqoGE1DeH5Wq9AYZs8mMvsBQ9DdcImUwm1k351/mjFhATnouT73ejj6fvfDYzAiXQf46NLYanhvamVYvNWXOiAX8PnkF6+dt5tUJL9Oud4sHfuYPqRgIQNShGKo1ezgueYtDdnoux3ec5sTOM5zZF8W5QzHkpOcWPe/u60ZwZBnqdahJQDk/fEO88S2sHXsFemJr93CuF2A2mclIziq68kiJSyMxOoXE6GRO7DrL5sU7i6427BxsCasWQmTd8lRsUIGqTSLxDfF+5IftXRJ95AJmk5nQyiUz+fDI1hN8N/jnou/Om9P6lljemZLyWA6hvFMnd59l5pCfObnrLOG1wnhjSm9qtHhw68gKIegT+S5+oT5MWj/6gR23uBl0Bg5tPs7+f45wcNOxonS2ao2KsGqhhNcqS/nqoZStGkzZKsGlNvVxQZ6OuJMJnD8WR/TRC5w7GMPZA9Ho8qwd/j7BXtRsWYW67WpQt30NXDycS7jEd++XUYtYMOEPFif+gIev2wM7bmJ0Mj98+Cvblu3GN8Sb/p/3pOULjR/5k6ec8XofWSwW/l24nZ8+nk9qnJYm3erz2hc9H1h7/W+fLWXumMX8eGxKidWK7oZBZ2D36oNsWbKTPWsOosvTY2uvoWrTilRvXpmqTSsSWbc8dg6P91qkZrOFmGNxHPvvJEf/O8Whf4+Trc3BRmVDjRaVafZMA5o/1whXz0cn4Bt0BnqVf5vyNUKZsHrEAzlmXlYeC8b/wZ/TV6NSq3jhw6fpPrRzqRkSKYP8A6DL17NsyioWffEnJoOJrm8/ycsjn73vta1sbQ69yr1FzVZVGfvnB/f1WMXhYlQSy2esZeP8beRk5OHu60bTp+vRsHMdarSo8tgH9VuxWCyc3hvFzr/28d/yvcSfvohao6JJt/p0GdiOqk0qPvS10mVTVzFr6Fy+WD+a2q2r3ffjZaVl07/KELLScmjzSnNeHdcD78BHP1HalWSQf4C0iRn8MnIh637ZjLOHE71GP0fnge1Qa+5fN8fCiX8yZ8QCJqwZ8dAuhBJ7KoFfP1vK1iW7UKltaPpMA9r1bkGtJ6rK4Yb3IOpQDP/M28L6X7eSm5lHxfrhvDLmOeq0rf5QBvvUeC39Kg+matOKjP/743sqY1qCFu9AL4QQt9zP/HHLqN+xFhUegQW274YM8iUg6nAMs4bO5dCmYwRFBND/85407lrvvnzxDDoDA+t8QEGOjlkHJz9UbdZ52fnMGbGIv79fj52jHV3ebM/Tb3fAK8CjpItWqhTk6dj42zYWTVpBSmwaddpW5+3prxIY7l/SRStisVgY2WkiR7ee5PsjXxFQ7u7S88Ycj+OTZyZj0BmYd27Gfa1APSruJshbh4094FudOnVEaWKxWMSuVftE30qDRBuluxjSYpQ4tefsfTnW6X3nxJN2L4qRXSYKs9l8X45xp45tPyV6lB0o2tu+KGYM+llkpGSVdJFKPYPeKJZN+1t09ewjOrn0En//sEFYLJaSLpYQQoiFE/8QbZTuYsW3a+96H0aDUfw5fbVY+d1a8UHbseL3L1cKIcRD8x5LCrBP3GG8lUG+GJmMJrFy5jrR3fdV0UbpLsa/NFUknk8u9uP88fXfoo3SXcwds7jY932n1v+2VTzp8JLoHfmuOLHrTEkX57GTGq8VH3YYJ9qqXxAz3p1T4if+vesOiXaq58S4HlPvKiCf2HVGZGmzhRBCpMSlCSGEOLn7jHgl/C2hTUwv1rI+imSQf0jkZuWJn0cuFE85viSetHtRzBo6V2Sn5xTb/i0Wi/iizzfiWZ++IjO15GrNGxdsE+00L4phbT4VWdrie3/SnTGZzGL2+/OsgX7QzyVW2406HCO6uPYSr1V/T+TnFtzRaw9vOS5eCHxNfPzUePFh+8/Emf1RVz3/RZ9vxNQBs4QQosRPZCVJBvmHTGp8mpjc91vR1uY50c2rj1g2bZUw6A3Fsm+D3iCSL6QUy77uxrlDMaKjU08xtNUnQpevL7FySFYWi0XMHDpXtFW/IFb/tPGBHz8pJkW8EPiaeDFoQFEN/HZZLBbxy+hFYtPC/4QQ1uaebwfNEbv+3l+0TXpypuhZ7k1x/lisEMJakbr02sfJ3QT5Yul4VRRlDtAJSBFCVL3V9qWp4/V2RB+5wOz353Fg/RECKwQwYFIvGnWp+1COirgdFouFQU1GkRqvZdaBSbj7uJZ0kW5IX2AgPTmL7PQ8cjLzyM3MJ/9S0jKdsShRGYBio6DWqNHYqbF3tMXR2R4nV0dcPBxx9XTGw9cVR2f7h/b/zWy28NGT4zmzL5qfT017YJOPtIkZDGk2imxtDlO3fkpYtVuvM5yfU4DJYMLF0xlFUfi443hqt6lO9/c6k5Gcydalu4g9Gc+rE14qWnZv1ez1bFq4jeCIMji4ODBgcq9Sn3fm/91Nx2txdVf/AswA5hXT/kqVctVD+XztSPauPcTsYXMZ020SNVpW4fUvX7nvQ72iDsdgMpqJrFu+2Pa5d+0hTu+LYtiPbzwUAd5kNBN3NonzJxKIPX2R+OgUki6kkRynJTcz/5avvxS0b6fCY+9oh0+QBwEh3pQp50tIRABlK5UhrFJgieZnAVCpbHh7+qsMqDGM5d+soe9nL973Y6YnZfBBm7FkJGcyacOY2wrwy2esYf64ZdRoWRlPfw/enNaXjq+1YfvyPRTk6fDwcyeyXnkSo5M5uvVkUSbL3Mw8jm07Rblqobw64aXHLsDfrWIJ8kKIrYqilC2OfZVWiqJQ/8la1G5TjdU/bGTeJ4t5q95w2vZuQd9xPfAuU7wZJmOOxzFnxAJMRjNJ0ck0ebo+/Sa+XCz73rjgP9x8XHmiR9Ni2d+dys/VcXTHWY5sP82JvdFEHYvDqDcBoFLbEFDWh4CyPlSsE4Z3gDsevm64e7vg4uGIk6sjTq4O2DvaYmuvsSYquyJBmdlkwaA3WhOR5erIzcwnNyufLG0O6SnZaBMzSYlPJ+lCGoe3n0FfYACs/78hEf5UrBNG9SaR1GweiecDnMZ/SUjFQOq2r8nGBf/R59MX7utVhzYxg/dbfUJqnJbxqz+mUoMKt3xNUkwKu//ez/eHv0Rtq2ZMt0ls+G0rfqE+ePq58+/C7XTs35rgioHocnVYLNarrFN7zpKZksWPx6cW5XKSbs8DG3iqKMoAYABASEjpzd18K2qNmi5vtqfVS01ZOME6/Xrr7ztNxlOoAAAgAElEQVR5/oOuPDesS7FMvz615ywbf9tGcGQgr33RE4DhHcaReD6ZgLC7G7N8pWP/naJ262pobB/cuOXMtBy2rdzPzjWHObrzLCajGVt7DRE1Q+nSryXlqgRTrkogQeH+qDV3N+lKURTUGhVqjQpHZ3s8fG9+lWKxWEiJSyf6eDzRx+M5fTCG7X8fZN2CHQCUrxZMoydr0KJrHYIe4Dj2Bh1rsWfNQZLOp9z1GPVbSYpJYXj7cWgvpjNhzYibJs7LSsvGzdv6WTq6OpCbkYe+wICHnzvd3+vMvrWHcPF0JqJeOFuW7KByowjKVglGY68hJTYNgMh64VSsf+uTiHStB/YtFUJ8D3wP1jb5B3Xch5WzuxOvTepFpzfa8cPw35j3ye+s/mEDr054idYvN7vrS9GstGy2Ld2Fk7sjvUY/B8DBTUfJz84vlqYVg95IWkI6wRXL3PO+bkUIwaGtp1j1y1Z2rTuCxWwhOMKfp19vTe0WlahSvzy29iWXidLGxgb/UG/8Q71p3NE6A9lsthB9LI4Dm0+yd+Mx5k/+m98mraJS3XI8+UpTWnard99PjmWrWnMbJUQl3ZcgH3c6geHtx5GfXcDn60ZRpXHkdbczGozMHjqPQ/8eo+kzDajevDI1W1Ulsl44sScT8C/rS+Mu9Ti27STxpy/S9pUWxJ1KYMprMwmpGMTxHado06sFwEPbD/IokFPISlhAOT9G/z6UY/+dZOZ7c5nUewYrZqxh4NS+N/zy3MymBf8RezqBD355G5VaRVqClphjcTTt1gAHZwdrb/s9fGH0+dZMiY7O9y+3vhCCvRuP8+sXKzl3JA5XL2e6vd6K1s81JKzyw32prlLZUKFGKBVqhPLCoA6kJWay+c+9/LNgB1PencfciSvpMeRJOrzc5L6lebj0f3Mpq2VxOr0vihEdx6PY2DB545ibpu1dP3cLqfFavto8lj2rDzLltVn8dGIq7r5unD0QTUilQGug71qPbwf9zLNDOtFzVHdqt6lG1KEY3pjS+7FdbKU4ySD/kKjatBLf7JrAxt+28dPH8xncdCQtX2jMa5N64RvsfVv7MJvMnNx9hueHdcXFwxltYgYHNx7j/NFYWr7YBLi6RnQ3Ad++cBHx/NyCO3rd7UqO1TLjwwXs23QC/1BvBk/txRPP1run3PFCCDJSc0hNyiIjLYfszHxyswvQ5RswGEyYTdZ2XxsbBY2tGnsHWxxd7HF1d8Tdyxkffze8fF3vqhnIO8Cd7m+25dmBbdj/7wkWTVvDjA8WsuLHfxkytReV6hZ/x3tBrnVxmeLOvHhg41HGPjMZV28XPl838oYZWK/8u6rUoAJu3q60faUFe9cdZOGEP2n1cjOWTfmLQ5uO0eHVVlRtWgk3H1fiTicQHBlI5UaRVG505xUc6fqKJcgrirIQaAl4K4oSD4wRQvxUHPt+nNjY2ND2lRY0fbYBv09awe+TV7Bz5T56fPQMzw3rjK39zTM32qhscHCy599F2/EO8mTVrPUYCgzUaVv9qiyAZ/ZHkXQ+hSNbTtCmV/M7auvU2Krx8HMj+XzqXb/PG/lv1QGmDv4VYREM+LQ7nV9teceBNT9Pz+nDcZw9Hk/0qUTiolOJP5+KobBj9v/ZqGxQqWxQFLBYBCaj+YbbBQR5EFzel7IV/AmvEkilmiF4+txe/iBFUajbqgp1nqjMzjWHmT1qCcM6f8lLQzvS472OxTpS5GJUEgB+oT7Fts9/F21nUu9vCIosw8Q1I67J7mg2mYuuTC4F+LzsAnLSc4uC/htf9WZQ4xF0fK01bXu35I9pqzi15xwZyZlobNXFWl7pMpmg7CGWfCGV2cPmsm3ZbgLK+TFwah8adqpz09q30WBkwktfY9AZqFi/Ag2eqk1EncvDJ3es2Mtfs9ZRrnpZXDycWDlzHcN/ffeOVqn/qOME0hMzmX1w0j29vyv9OXsj349eSmTtsnw0uz9+IbeXIlYIQczZZHasP87+/85w+khc0dJ7vmXcCQ33I7icD/5BnviUccfD2xl3T2ecXR2wd9Bc02RisVjQFxjJzdGRk5lPRloOaUlZJManE38+jbioFOJj0orG1oeU96V2kwo0eKIS1eqWve0mmLycAmZ+tJiNS3bTsltdhs3oU2zNN1Pf+J4tS3ayLOWne15aUAjB4kkr+Omj+VRrXolPl394VROK2WTmx+HzMRlNNOpcl9ptqhc9l5ag5aMnx/PhvHcIr2lt1pk9bB65GbkM/elNMpIz2bx4BxazhWeHdLqncj4uSnKcvHQf+IX6MHrJMA5sPMp3g+YwuusXt1xcXGOrYczSYRh0BvQFBhyc7YueO7b9FHNGLKDfxJep+UQVHJwd0OXpyc++s6aX6i0q8/PIRWgvpuNVDEM/183fzvejl9KkUy0++LbvbXWm5mYX8M8f+1i3dB+xUSkoikJEtUCef60lVeuWJaJaMC5ud95vYGNjg4OTHQ5Odvj4uwHXfs4GvZGoExc5fuACB3ee4+9Fu1k+bzuePi606lKLTj0a4hd48+ybTi4ODP2mN8EV/PllwgpUGhVDp/e+5w5GIQT71x+hRssq9xzgzSYz3w3+mZXfraPli014/+e3rmo2E0Lw7aCfyc3MpUHHOiyetJy40xd5sn9rbO00eAd60eyZhvw+eQVvfNUbT38P6nWoydFtJzGbzXj4udPt3Y73VEbp1mRN/hFhMppY/s0afh27BKPeyPPvd+XFj7rdtN11x8q9ZCRl8tSAtmQkZzKx53Ra9WhKh1dbFW3zfutPeGpAW1q+0OS2yxJ7MoH+1Yfy+uRePDv4qXt6X2cPX+C9pyZTo2kkY+YNvOXIk9zsApb8uIW/FuyiIE9PxRrBtO5am8Ztqly36UQIQbo2lwsxaVyMTyc5KQttWi6ZGXnk5erRFc56BWunqZ29BkdHW9zcHfHydsHXz5WAQA+CQ70IKOOBjc21QbggT8/+7WfZuOIAe7acBiFo2r4avd5pQ1DYrZsg5n+5it8m/83bk3rwVO/mt/nJXd/pvVG803gE733/Oh36PnHX+9Hl65n48tfsWLGX54Z2pv8XPa9pUsrLzufjJ8czYc0InFwd2bvuEHtWHyCyXjhtelrfhxCCr/rNRGOnpkKd8vw1cx1tejaXNfe7JPPJPwa0iRn88MGvbJy/Db9QH978ui+Nu9S74fYZyZl4+LlzZOsJNi/aTp9xLxatUv/lq9+RlZbNZyuHY9AZ2Lp0F8IiaPtKi1uW453GI8jLyufHo1/ddXuyxWJhcIcvSE/OYubmUbh43HgkhRCCDcsP8OPkNeRk5tOsQzWe69+C8MrXDuWMu6Bl3+4oDh+4wMnjCaRrLy8OrlLZ4OntjLuHEy4u9tjZWSdEKYp1+KNeZyQ3V092Zj7atBz0V7TlOzjaElkxgOq1Q6lTvxyRlcpcU1tOTcxk5fydrFq4C6PBRLfeTen1TpubdhxbLBZGvjiDU/vPM2fX2Hsa6jpt4A9snL+NhbEz73pkSkZyJiM6TeTcgfO8Nf1Vur7V4YbbTuz5NRXrV6Dbux0pyC1g69JdnN4bxUsjnima4JeelMGZfdFsXLCNeu1r0q53y7sqlySD/GPl8JbjfPPWj1w4EU+jLnV56+tXb9pxtWnBNvZvOML7c94C4Mfhv3HuUAwfLxjE7r8PcGTLCTKSM1GpVRTk6m65YPimhf/x+Ssz+HT5+zR8qs5dvYc9648ypud3DP2mN22eb3jD7Qry9Hw9+g+2rD5ClTpleePjztcE99SUbNb9fZhN644RF6sFIKCMO1WqBxNRMYDQcj4EBXvi5e1SFJgtFoFOb8RkulyTt7fTFD0vhCArM5+L8RlciEnl7OkkThyNJ/pcMkKAu7sjLdtUoX2nGoRHXD3ZKSMth1+m/sM/f+yjbIQ/o77pSZmb9DPERyXzetOxPPtmW14d1e3OP0ysFYDeEe/S+qVmDJk94K72cf7oBUZ2/pys1GxGLn6vKKXAjWxZspN9aw/SZ1wPvAI8OHsgmo2/baV93ycIqxZKzPE4giIC5IIfxUQG+ceMyWjiz69XM++T3xFC0HP0c3R/r9N1v1DnDp1nwktf07RbfQw6I3vWHGTKlrHsWnWAM3vPUatNdZo90wCAz16YwsCpfW6aasFkNNE7chC+wd5M2fzJXbUlj+//A8d3nWPewQk3HEWTk5nPqNd/5uyxBHq925bnX2tx1ZVD3AUtC+b+x6b1x7CYBTVqh9K0ZUUaNqmAf4A7FosgNj6dU2cSiTqfQmx8Oskp2WjTc8nN0xd10l7J0cEWDw9H/H3dCAxwp1yYL5HhflQo74dGoyI7K5/9e86zbfMpdm0/g9FgpnqtEF7u05Ta9a4eErl362kmf7AYlVrFhJ9eJSzyxgu/f9pnFif3RjP/6Od3dXU0a9g8ls9Yy0/HptzVSlH71x/m0+e+wsHZnk9XfHhVh/2NaBMz+GPqKpw9nOnxkfXkNLjZSAZO6YM2MYOc9Fxav9wMlVolJzQVA9nx+phRa9Q8N6wLLZ5vxLeDfuanj+azcf5Whnz/BpUbRly1bXjNMMYsHcq/C7dToXY5nhvWhYIcHaf3nKVOuxrUaVcDsM6OPb3nXFHQNRlN1z1pqDVqenz4NNPf/oldq/bTqPOdrUhmNlvY/+8JWnare8MAb9AbGTNwLlEnExk5vSeNWl8eAVSQb+DXOVv5Y/EeNBoVT3evx9Pd6xEQ6IFOZ2TXvmhmz93K/kMXyCrsWLa1VRMc6EFgGQ9qVA3G1dUeRwdbNGoVKAomkxmdzkh2jo70jFySU7LZtO0UK9ccBsDOTk2NqsE0bhBOy6YRPNG2CjnZBaxddYg/f9/Dh4MW0KBxOO8OexJff2vemnrNI5mycCDD+/7Ix/3m8PXvb+Fbxv2677dpp9rsXHOY6GPxhFe/s9QfaQnp/P39Blq92OSuAvzaOZuY+vpsQisHMW7VR7c9N8MrwIPGT9fnp4/mUybcn8h65dHYqlHbqmnU+dHNtFqayJp8KbJjxV5mvPMTaQnpdHqjHf0mXk7Tej2rf9xIzLFYXhz+NJ7+HiRGJ/Prp0uo2rQSHfu3RpuYwYLxy3iiR1OqNql4zetNRhMD6w7HaDDxw+Ev72i6fkJ0Cv0bjWHItF6069H4utvMnriK5fO28/G0l2jW/vI4/7gLWsYM/524C1o6dK5JvzeewN3DicSkTJYs38+6jcfIzdPj7uZIg7ph1KwWQoVwX/RmM1HxacQlZZKszSYju4DcfD0Go7XdXa2ywcnBDjcXB3w9nQnycycs0AtvVycuxGo5ciye3fuiib+YgUplQ9OG4fTo3oBKkQEY9CaWL93Lr3O2olLZ8MHILjRufnlCT1x0CkNenElQmDdfLRh43ZEvt/OZ3MikPt+yZclOfjw2hYAw39t+ncVi4eeRi1j0+Z/UbV+DkYvfu+nfzI3sWXOQrUt2cmLnabq82YGn33nyjvch3ZqsyT/mGnetR81WVfll1CJWzFjDzpV7eefb/tftmLVYLJzYcZqGnevg6e9BRnImy6auIiDMj+otrDVmi9lCcMVAPmz7KZM2jLkmzYJao+a1z19mZJcvWPHtWrrfwYgJbVImAD6B128Sij6VyIpfd9CpR8OrAvzJ4wl8NGQBarWKSdNfplbdMHJydUybuYGVqw+hKNCiSSSd2lfHx9+VrfujWLXnJMcXbEJfGMw1ahV+Xi54uDri6eaIrUaNAhhMZvILDMQkaNlzNIZ8nREARYEKIT40qhHGyI87Y6vYsG7jcVb/c5Qt28/QqH553n29Nc+/3IhmLSsybtQffPLREgZ/0JGOXWsDEFzOl7fHdOWLYYtZt3QvHV9ocM179iljHXaZnpx1258jWBPGbZi/jR7Dn76jAK8v0PNlv5lsXmTN/PjOt/1v2nZ+/ugFAsr7X3dE16UMq4qi3Ld0DdLdkUG+lHF0ceDNaX1p9VJTpg6YzZinJ9Hi+Ua8Nb3fVYtI2NjYEFDejyVf/YWHnzs/Dv+Nqk0r0fjpegRVsLYb+wR54RPkhW+INx5+10+bW69DTRp0rMXcT5bQpGu9206IZTRYA+6NxsQv+XELDk629B7cruh3sTFpfDRkAW5ujkya/jJ+Ae4cOHyB8V/+TXpGHp061KDXCw2Juqjll9V72XPsAmAN0F1bV8fX2wWdxUxGXj5xqVlos/NIzs5Edyn4q1S4ONrhGexKlRpB+Lo546jRkJut48jpBH5btZe5K/dQPsibHh3rsOD5/vy19gi/Lt5J74FzeOf1VnR5siZTZr7C2I+XMm3Sapxc7GnRynrSbNGxBn/N38Wi2f/Svnu9a2rzGjvr1/FGs26vx2Q0MWPQz/gEe/Hi8Kdv+3WZqVmM7voFp3afo9+El3jhw6dv2LRi0BtZMH4Ziz5fzksfP8Mrnzx/3e1k5+rDSf6vlFIV61fg272fs/iLFcwft5SDG4/x1td9eaJH06Iv88sjnsWoM7Lrr33U61CLHh91u+qLvvOvffwyehGDZ79OmfL+WCyWazoEFUXhnW/6MaDW+0x/+ycm/P3RbbXDXqoNXi+JVl6uju3rj/Pk8/VwdrVOaDIazXw2Yhmawhq8X4A7azccY9K0NQSW8WD8qG6o7NWMmrmaw6cT8PFwpt8zjfDzd2X3mThWHT5FZp41p4uDnYYgbzfcXRzwdHdCpVJAAWEWGIxmsnN1nIpLJbMwP4+dRkXdiGAG9WuFRWdm9bYTjPt+HQE+rgzr3Zp5s/sz+eu1fPXNP0RFpzJoYBs+mfAc77/7G19NWEWFCH/KBHmiKArd+jRl/KD5HNkdTa3G4Ve9b11h8jc7h5unr7jSws+XE33kAmOWDsXByf7WL8CaRXJkp4mkJaQzasnQog736zm5+yxf9fuOCyfiad2zmWyGeQTJIF+KqTVqXh75LE2fqc9X/Wcysed0ti7dyTvfvoZXgLVpoM9nL143UdnxHaf5ZfQieo7sTo0WVQBr7V+Xr0d7Mf2qGbe+Id70GfsC3w35hX/mbqF9n5a3LJtX4ZVBakL6Nc8d3Xseo8FEk7aXV5JcvmQPMedTGTf5BfwC3Nm05SSfT11N7RqhfDqiKys2H2Pm4m24ONkxtHcrLLYK89bvJzE9Gzcne+pWDMbF1Z5MvY5YbSZnU9IxZKddt2wqG4UgHzcqVw7Ay9EBs97MoTMX2X48Bmd7W55tXp1enesx549dDP3yTzo2q8LYj7syd8EOFi3bg1ptwzuvt2bkZ8/Q76VZzP5mA2O/sNZ+6zWPRK1RcWD72WuCvDbR2oTleYOrpv93/lgcCyb8SaseTWnS9cZzJa50YudpRnf9AkVRmLzpk2s66C8x6AzM++R3lny5Eq8ynoxb9RENOta+rWNIDxcZ5B8DoZWDmbrtM/6Y+jc/j1rEa9Xe491v+9Pi+cYoilIU4P9dtB03H1ci65Zj1tC5tO/zBC2ev9wBuHDinyRGJ5NwLhEnN0fGLBuGSmVtf+3yZjv++3M33w7+merNK92y2cYnyBM7B1tiTl685rnowt9FVA0CwGAwsWTBLurUL0eDJhWIidXyxbS1VK0cxITR3Zi+cCt/bDhMi7rhPN2uBl8u3cL5pHSqhfnTukEFjielsO7sOSxC4OpoR6ivB7UiAsEGTEJgtlhAAQUFjY0NNgL0OhMnL6aQnGWdSBXu70WPerVJSs5m3vp9uDk5MKhbU5ISsvh5+S7ikzP4alg3zGYzS5bvp2JEAG2fqMzzLzdi7g9bOB+VQlh5X+zsNZSt4Mf5M4nXvO8Lp62/Cyp/6yYvg97IpD7f4uzhxBtfvXLL7QG2LdvFxJ7T8QnyZMKaETdMjXF67zkm9ZlB7MkEnuzXmte/euWuOmOlh4NcJPExoVKpeG5YF2YdnExguD/je0xjwkvTyNbmFG0TUM6Xr/p9R9+Kg2j+bEOeGXQ5ZcG8T35n04JtdB7Yjskbx2Brr2HhhD+LnrexseGDX97CRmXDl/1nYS5M4HXj8tgQUSuUY7vPXfNcWnI27l7O2Dtamy0O7oshIz2Pbs/VQwjBtO/WY2urYuxHXfh5xW7+2HCYnp3qUbd2WQZ9t4ICvZGB3RqTpTIy57/9xKZl0qx6GDUrB5JrZ+KANon/4mLZERvLqfRUYvOyiM3NIiornd0J8WyNvcCelItcFHmEh/vQomZ5VCqFudsPsDfpIn27NCDEx42xv64n1VDA2Lc6ciIqieHTVtK/d3OqVCrDN7M3kpOjo3O3OqhUNmxYe7To/fkEuKNNzr7mfZ/cF43aVk3ZSrdekOWXUYuIOhzDe9+/fssZspeSjH363FeE1yrLN7smXjfAGw1Gfhm1iHcbfUxBjo4Ja0bw3g9vyAD/iJNB/jETUjGQaf+No++4HmxbtpvXqr3H3nWHAGs7/hfrR6PWqK+aJLR2ziZWzf6H8X9/RIXa5bCxsaFSg4hrauu+wd68ObUPR7ed5LfPlt6yLDWbVST6WDxZV6QdADDojNhd0SG7f080dnZqatUL48ixeA4eiaXPy004l5DG3JV76PpENTz9nJm4aBMNKofQsE5Zvt6wg3yDkXb1Isi2NbIhOpoL2ZlUDPElIswbB29b9C6CVJsC4s05xJtzSCKPAkczKg81oaEeVC3nT4EwsT7qHKdz02lZpzyBXq7M2rQbjYctL7aqyR//HWXj8SiG92vL/hNx/LpqL4PfbEtWdgHL/z6Im7sj1WqGsH9PdNH7sXPQoNcbr/k89v97ksp1y92yTf7Q5uMsm7aaTgPa0OgWM1LNZmuSsR+H/0bLFxrz5aZPcPW6NsdPzPE43mn4MfPHL6N1r+b8cPQr6rWvedN9S48GGeQfQyq1ipc+foYZeybi4unMx0+O55u3f0RfoCeoQgCzDk7GRmVDtjaH+DMXWTlzHZ/88T6+Ida0CUaDka1Ld2J7nWDUtldz2vZqzsKJf3J8++mblqNOy8oIIdi36fhVv9fYqotG3wBEnU2ifAV/bG3V/P3PEZyd7WjfuipT5m4iNMCDRvXK8eWSLbSoUR7FWcXvu4/SrlYEFheFVadOE+zjRniYNwkil71pCaRbCgjwcaJ8oBv+/o64+9ji4WOLj589ZQNdCQ1wwagycyg9kbO6dAKD3aleLoB/z0VzLj+DpxtVYf/5BDZHx9DvqQZsOHCW+Jwc2jWuyNwVu3F2tadOzVBW/3MUIQQVK5chJjq1aHESo8F8zUiU+HNJxJxMoGH76txMelImE16eTmAFf16b1POm2xr0Rsb3mMbyb9bw7OCn+Gj+oGvWJLBYLCyd8hdv1v0QbUI6Y//8gA9+fhsnN7kiU2khg/xjLLxmGN/u/ZxnBj3Fyu/W8WbdD4k6HIOrlwvPDe2Mq5cLRr2Ryg0jrlqp54M2nxIcGXjdURmKovDW133xDfVhQq/pZKTceMx3hZoheJfxYOvyqyfGuXk6kZWRh7kwp0xqcjZ+AW4IIdi97zyN6pVnz/ELXEjMoO8zjfh80b+U8/fEwd2WTcejeLpRFTbERKE3m6gVGciBzCRS9DlUDvFC4w4J5gyiClLIsynA0Ung5abC002Fi7OCSa0nVq/lvCENs7OJ8kHuCLWFnSlxBId44O/hwpJjx2lXP5K03Dw2n4umQ71Ifly9m46FHdSL1x6gRZMILiZlknAxA19/N8xmCxkZeQBkanNx+79kbFtXHgCgWdcb18wtFgtf9p9JXlY+oxa/d9PRNLmZeQxv/xnblu5iwORXeGNKn2tGRqXGa/mw3WfMHjaPuu1r8P3RKTS+zQ5c6dEhg/xjzs7BjoFT+/D5upHkZubzToOPWDrlLy7NhE67mMHxHafJycgl8Xwy416cgr2THcPmvAnA9WZMO7o4MHrxELJSs1k/b8sNj21jY0Pr7vXZ9++Jq0bZBIR4YTZZSIrPAEBvMOHgYEtqWg6ZWflUqRTIPztO4ePhTHJuHimZubRrXJGVB07SqV4llp84Sai3O0ZHwaGURKqF+pBmk0t0fhphPk54eSjgoEMrtCQZU8gSGWSLDFJNaaSY0zDZ5uPkZibc15lMcy4xRi3lgtxJ0+cRpc+gYYUQ/jx6gnZ1IziVmIqtqy0ujnb8/t8RmtcJ558dJ4mMsLZ5n4lKwb6w6elSE83FC1oCQi5PAhNCsGnJbqo2qoB3wPVTHoB1uOS+dYd5fXIvwgoX676e5AupDGk2ipM7z/DR/EE8N7TzNdvsWrWf12sO49Tuswz5/g3G/vnBVfMopNJDBnkJgDptazD70GTqd6zF7GHzGPHUBDKSM6nXvib1n6zFey1Gs/iLFTi7OzNxzUjg5mvEhtcKY9b+L64bYK70ZK9mIAQr51w+GZQv7Hg8cyweAFuNCoPBRHKKtbMywM+V/cdjaVA9lD+3H6Vm+TIs23eM8v6ebIu7gK+bM8nkY8SMj68DxzITifR3xWSXT7IhjSA3DUHu4Oakx9FBh40mG0WTjb19Pq5OenxdTYR72JNlySTLJpMIfxeS9VkU2Onxc3NiR0ostcLKsPTIcVpXr8Afe4/Rtl4EW45EUS2yDBnZBRgshVchadlFWS41ahXpKdlkpOUQVuFyfpmDW0+REJ1C+5dunMrgyNYT/Dp2Ca16NKXzG21vuF3c6QQGNx1JaryW8atH0KpH06ueNxlNzHrvF0Z1+Ry/UB9mHphMx/6tZY6ZUkwGeamIu48bY5a9zzsz+nN483FerzmM/esP03dcD0Yufo+3pvflnW/7Adamg1sFhqCIMrfcxi/Ei0ZP1uCfBdvRFxgACIv0x8nFniO7rZ2Vbu5OZKTnUVCYZsBgNpOTr8fXx5XYlEzKBnuRkJFNUBkPtPn52LtryDcZUZwFmaY8fL3UxBakUNnbCVv7PLLNqQQ7Q4S7IMA5Bz+nXPyccvF3zqG8m5FwNxUmJQOhzqa8py3pxiyw1+HtbE+8OZMgD1dO5abh4exAbH4Wig3kWn2PJocAACAASURBVIwIATlG63tITLOekAoKjGRm5APg6u7I0X0xAFSuHVr0GSz5Zh0evq40v0FTTVpCOhN6fkNAeT/e/bbfDT/T03vP8V7z0ZiMZqZu/fSqdX0BUmJTGdJ8NMum/U3Xtzrw9fZxRbObpdJLBnnpKoqi0OXN9szY8zmuXi4Mbz+OX0YvIigiAI2tpmhcfHEuPN2l3xNkp+ex4fddwP/au+/wKKq3jePfM1vTK6H3XhQVRFRApQkqIIhKUymKdEF67x1piiAgRSwgiBVRECygIh0B6TWQhCSkbi9z3j82RvkBwqtIIJ7PdeVK2dnZZ2eTeydnZp4TOL3yzlpl2P3TMaSUFC4aybmzF3NnZUrP6Spp8wQC9YLdRkSwhZ/PxHNHyUIcuphCkYIhZPucBIXoCOGlQJhOli+NSpEGwq2ZCJFM0SA3VcOhUpiLSmEuqobrlA7RMWupmIwXKRfhx2x0YTRnUyTUQrJ+kZgQM5kGBw6vh+KFIzmYlMwdJQuz6/R5woIsJKQFwj013Zb7XM7HpxEVHYLVamLH94cJiwiiXNWiABzfH8/eLUdo1a3BFScW8fv8TGo/B2e2k1GrXiU47MpTGu7ZvJ9BjcZhDbUy64dxlL6j5CW3/7JuF91rDObsoXOMWPUqvV7vcs2J4ZX8QYW8ckWlq5Xgje1TeLTjI7w34SNGPDGZrLTsKy574UwKH8364orj89ejTLWiVK5ZhpWzv8rdm69ZryLJCRmcPnaBsuULcSEpMzfkM7ICe8ZZzkAbgBPJaZQoFIXd4yHF56BIVBhHM1MoHG3GjwdhthFjlURYM7Easqgc7iXSnIpFS6CQxUHpYD+lg/0UtboINlwg1JhM2VA7cVYPaClUiLSSoadSPCwIt9FBptdJ2SJR7E1JwmTQ0Cwa59OzKBgTSnJGINyzc9o1hIVZOX40idJl4/B4fPy8+Tdq16+S27dm9dyvCQq10qTDpcMqv1syYiUHfjxCn3kvUvqOK7cf3vrxLwxrOpG44rG89t1YilX44zx7v8/P0hEfMKLZFAoUj2He9ik89PT9f+t1Um5PKuSVq7IGW+j/dnf6LujKvm8P0qvWUI7vOXXZcusWbmRB/+W83nMxuv7XF0H9LjvdxtZPttOj1lDmdF+ESfeSmpDOV+9uBeD++pXRDBrfr9tHpZxZoDJTAwGakZGzJ+90YzZpJKRnoRvAaNQ4m5WJZpXEhVlJcGYQEyqJtWr4RSoVwowEGZOJM3spF2wj1pSClVNEaslEaikEi7OEGxIoFZRO2WANjQSqhguyfUmUj7CQLdPwSTelY8I4bU/Dq/spUiCcZEfgrBmD0UB2zhuPyxV4swoPtXLqRDJV7ijG7h+P4bC5qdskMIxyfH88P3y2i+ZdHiY04vILjr778CdWv/Y5T7zciAbtrvwmsPGd7xn/zEzK1yjDrC3jL+kDn56cydCmE3l/0loe7fhIYHimwrUvtFLyFxXyyl8SQvB410bM+G4sXreXvnVGsPmDrZcs02lCW54Z2ILPF2xgeqd5uac+Xs3540l8Nn8DP368nebdGzP8/b74PV4KF4vkw9c34HJ4iIoNo2ad8mz8eBflKhTCZDJw5GACUZHBpKYE/qNwuX3oEnQpsXu9REcEI4Uk3p6J2QJFwyzYfNmYDJmUDzXi0ROoHAImzlDKaqCi1UYxcwrh2mnCtVMUNCZR3pJGxSALQo+nXIgHSKNSuMDuTyHCLCgWYSXRk47T7yU02ARGQYotEPJeXc+9iMzpCBw/yEy1oeuSGrXKsPmzPYRHBnNX7cCMS0vGf0x4TCite15+IPXs4fPM7raIKrXL02PWC1fcjuvf3sT0TvOo/nBVpm4Yecmcrr/9fITu9wzkwNbD9F/cnQFLemAJuvqk70r+pUJeuS5ValfgzV3TKF+jDJPbz2F+v2X4clr0CiF4cUp7Oo5rwzcrfvjLoNd1nS/e2khGchZth7SkSadHgMCVuPc1rErahUw+enMDAI+3rU1aSja7thzljrtKsO3HY1SqUJgTxy9gNhlwu714cy4wcvq8aAaBZgCEJNmdidA8lAozE2Tw4JUXqBiqIUihYpAfo36SYuYIKlpDKG/KorwpkwoWKG0thEmeo7Q1g0iDm6IWGx49mZLBBqKtPi560vFLL2YThAaZcfq92DweJOBwe9DIGVK6aCM8zMrRg+cJC7NSsEAYP31zkAYt7sZkNrL/52Ps+f4QT/dsdNlevD3LwaiW0zEHmRn23itXbOG7+rXPmfnSAmo0vpMJnw8hKPSPsfpP3ljPqw+Nxmw18/q2STTpXP8fvfbK7U2FvHLdouIimL5pNE/2bsraOesY1HBc7sVOQgjaj3iKzhPbsem9LYxtPQPPFS7dn//qO5w9fJ6OY5+hROWieNxetq3bRXaGndavPEa95jX48PUNpCZmULNuBQoWjeLLVb9Q75HKnI9Po0ThKM4nZlCiUBRZOQdgITCdoE/qWMwGhKajaRKb34FbZlIkCIpYBD5/EqXMdoKFm3JmH2b9AGGGWApbH6CwtS7RpvKY/EcpbrxInCmCKC2BIM1F2WAzmkjD6c/CpPmxmiDMakIXEpf/jytz07OcCCkRAk4eT6Fq5aJs+/E49z1Yns/f34bUJU+0ux+/X+etkauJKRzJE50evmT7SCmZ9fJCkk4lM2pVP+JKxF52+3sTP2LhwHd46Jn7Gffp4NwDqD6vj3l9ljCvzxLubXoXb+6cStnqpW7Qq6/crlTIK/8vRpORnnM6M/TdPhzZcZxetYZwfO8f4/Rth7ak9xsv8vNnOxn71HTczj/6xft9fjxOD73mdiYkIpiLCWlsX7+H3d/sp8p95YmMi6Dj8Bbofp0VUz9H0zSeaFubX7efpHjRSMxmI2lJgTeVcKuF8wkZCCDUYkZK8Pr9SEBoYNB0zAY/Ag8uPY0QQxalrcEYcVDIkIbVEEVpSzVi9WME+/YQ7NtFpP8gpczFiDZXIUweIsYYRFGzC5+ehFG4CTF6CTdDqNmI0MAn/ehSYjYYEIDL5cVh91AsLpLk5CwKRIRgy3bxYN0KrP/wFx5oVJUiJWL4+r0fObE/nq5jW+c2Yfvd+5M/5oc12+gyqS3V6lw+5eL7k9aybORKGnSoy9B3X8FkDpyRY89yMLL5FD55I9DCYMzagZcM3yj/XSrklb+lfru6zNoyHqlL+tUdyc4N+3Jva97jUfq99TI71u9lTKvpeD2BPXqD0YDRbGTeK0tZt+gbPnnjK/ZsOkCRsoV4oltjDEYDhUsVoPmLj7Dhg5848Mtxmj5Ti5AwK2uXbqXuI5XY+eNxShaPwZbmwOvxYzUZCTaa8Hr9ODxe3H4fINEEGIWO1eAlSPOCdCBkPMXMFkzCSGHhQdMzERFTEHHbEXF7EVFLEFo0UfopIk1VidbS0PWLhGoeIk2SSLOBIKOGJiR+dPy6ROoQbDJh1DSQkJiYQbjFgqYJ4o9eoEBcGCnxF7FluWjx3IPYs528O/0LqtxbhrrNL+3PvmfzflaMXc0jbR68bCpFKSVLR3yQG/CDlvXKnWYv8eQFet83lN3f7Kffwm50m9kx91RXRbkhIS+EaCKEOCKEOC6EGHIj1qnc+irUKMvrv0ymSNlCDH98Ep8v2JB722MvNaTfwpfZ+fU+RrWYmrtH3/v1zlSuXZ7TB+Ipe1cpHn72fp7s1eSSScBr1a9MVEwIM3suwRps5qnO9fjl20PceWdx7DY3JQpFcuZEKmajgZjQYNwuL1k2F263D6kDCAJBLzEKHbPmI1jzYMKLpp+nkKkMAomI+RAR1AqhRSK0YISlDiLmAzCWJ1rYkNJGhMFIlMlEiCHwpgE6fqnj8+t4fH58Pj9CD/xnYTUZ8Hl0Us5nUq1SEQ7ui6dps7tZs/gHqtUoRbUapVg+6TMyUrLpOv7pSy5qSjiRxIS2cyheqSh95790yW1SSt4Z8yHvT1pL0y4NGLi0Z+51Crs37afXfUPJSM5k2jejeOzFBv/mS67chv5xyAshDMA8oClQBWgrhKjyT9er3B5iCkcxa8t47m1yF3N7LGLBq8vw+wMHXZt2aUD/xd3ZteFXxj71xxh9+2Gt6DmnE9UerMQddSrnriv9QgYT2s7m6yXfUrBgGOcOnGbRiFU0a1eboBALuzb/RumyccQfuYAAikSHY8twYs/2oHslBDIYKUFKERi6QWIUYNUkQUJHA4z+4xD0DMJY8rLnI0QQIrQ/Qr9ImBZJmMGECR8CH78HvMvvx+n14fB4cTi9ZNvcuB1eIq1WIkKsZKTasaBhMhnQHW7SUrLp+OqjHP/1LF8s/Z7mXR6m4t2lch/T4/Iwsd1cdL/OuI8HEhR6aeOxFWNX8+74NTTp9Ah93+qau5e+ccX3DGs6kehCkczdNjl3Bi9F+bMbsSdfCzgupTwppfQAK4EWN2C9ym0iOCyIMWsH8mSvpnw0ex3jWr+GM2c+1Sad69P3rZfZ8VVg6Ob3M3KS41P5cvGm3OX2bz3E0McmU7hMQbpMbsfMzaMoVb0Mn76xnuw0G6061uGnb36jXr0KJJxNo0LpgqQnZuPIciP8YEZgkhqariF1gS5zPhBIqSPQMWsaZjQEOsLyFxcEWR4EzARpYRjwAn4kEq+u4/T58folui4QfoHwg/RK3A4fmal2Ii0WwkIsHNp1hnr1K7N+5S/c/UA5KlUvwfxhqwiLDuG5wX/085FSMrvbIo7tPsnAJT0u6dEvpeSDyR+zYtxqGnd8mH6LuuXuwa+c8jHTXniDO+pWYvbW8ao9gXJVNyLkiwLxf/r+XM7PLiGE6CqE2CmE2JmSknIDHla5lRhNRnrO7UyP2Z3Y9vlOBjwyhuycS/sfe7FBzhj9HiZ3mIvX4yWueCyPvdgAo8lI/JEE1s75kqf6Pk6XiW2JKRyFEIInOtZDM5tYMuETnupUl5iC4ezYeJAixaJwpztxZbkJt1qItFqwYMTo1zDpGmbdCLoRn67h0w340QLT/Ekv2u/DIOLqsykJYQYtFIEfKXX8UuLWdZx+Ha+u4dcF0i8w6BpmqRGkGTCj4XfrJJ1Oo3hcBD6fn1CDRmaaned6N+Kzxd/y246TdB3bmpDwP053/HLxJr55bwvPj36aB5rXvKSO1TM+Y8nw92nQvi793noZTdPw+/3M77eMt4e9zyNtH2Til8NV73flL92IkL9St6TLrm+XUi6UUtaUUtYsUKDADXhY5VbUss9jjPl4EKd+PcOAR8aQHB+YLPuxlxrSdfrz/LD6ZyZ3mIuu68QWjcZkNpJwIonoQlHUfuKPA5Fb1v7CexM/ok6Lmmz9fDe//niU9j0acHT/OWrcU5Lzp1IpHBeORRfY0py47V58Lj8+l47HLfF6BT6/GY9uxKMb8UojPinw/d56Qb9w1ecgdQfo6Xj1bNy6B4fuw+knsB7dgPQbMfgMaB6BdEukS8eiaxQIC8aiGThzMIk69Sry9YfbefiJ6sQVCmfFtC+o2aAq9VvXyn2cgz8e4c2+y6jZuDpth7a8pIbPF2xg0eB3eeiZ+xm0vBdGkxG3083k9nNYO2cdT/ZuyuB3el+x342i/NmNCPlzwJ+bWxcDLp+dWfnPuL9ZTcZ9NoSk08n0f3g08UfOA/B0/2a8PON5tqzZxtwef7RA2L/lEEazgbCoUACWjlrFu+M/ov/ibvSf/yIlKxZm8gvz8GdmUbxsHLs3/UaRopFIm5esC3ZMukawMGCRGma/hubV8Hs1nF4Nm9eMw2/FqVtwSSNedCRGpOfnqz8BzxZA4pQO7LqOSzfi0s04/WbcXhNuj4bfI5AeCR7Q3RJ3tgd7soPiBSPR/X50mwshBJ36Pcq8ISvx+/x0n/hM7gHVjJQsJrafQ4ESsQx5p1duLxuAze9vYW6PRdz3+D0MWt4bTdOwZzkY1Gg8P6zexktTO9BzTmd1Bo1yXW5EyO8AygshSgshzEAb4LMbsF7lNlazcXWmbRyFy+aiz/3DObz9GACtX21Gm8FPsm7hRt4e8h5SSp7s1ZSfPtvJWwNXMKjxePZuPsCo1a9Sq8ndmC0mXp37PC6Hm3mvLOWZjg9w4Vw6FcvEkZaQSVxUKKHCgN/uR3foSJdEuAGPAa/PiNNnIstnJstvxS5D8BCEExM41yF122V1SymR9mX4CcYhJQ5pwaYHk+0LIttjweU1IL0GDF4Ns0/DogvCNBMRFgt4dM4dvsB995Vl24aDPPn8g+z/8Sg/r9/H80OaU6R0HABej4+J7eaQmZrN8PdeuWTO1W9X/sjU51/njnqVGfnhq5gtJrIuZjPssUkc2X6cEav68cxAdchLuX7/OOSllD6gF/A1cAj4UEp58K/vpfwXVLy3HHO3TSIsOpQhj07gwI+HAeg8qR3NujXmwxmfsXz0KmKLRjNmzQBqP3EPj7/UkDlbx1O03B+TapiMBu5+uAp+NFLPJFO3yR1s33iAcuUL4st04browuABs09g9ILBI5BuDY/LgN1tweYJIstnJd1n5qJPI813ESmzkJmDkfJ/2i8414B3F2m6DZtuId0fRLoviEyvFYfXjMdjQndp4ALplEiHxJ3pxnXRSeGoMILMRhKPJBJXOJInnq3FotFrqFSjNC1f/qO1wPLRq9j33UH6zn+J8veUzv35gR8PM73TPKo+WImJ64ZhCbKQkZLJ4MbjObbrJMM/6Eu91qqDpPL/c0POk5dSfimlrCClLCulnHgj1qnkD4VLF2TG5tFExkUw9NEJHNh6CCEEvd7oQpNOgTbGXy3ZTNnqJan+UNVL2uB6XB5+/eE3vvvwJwoWiaRw+SJ8seIn2nWvjwAsPj+2iw5iw4MJ0gXCqSOcoDlBcwmkx4jTYyLNYyXVE8JFXzjpeghZ0kKaLsC9EZnRC+nZhfQeQc+egcwajpsgMnUvqX4Laf5Q0r0hpLuDyXab8LqNGNwGTB6B2Sew6oIwowmL1Eg+mcqddxTj3PFkXhzUlAXDP8Se5aT39Ha5Z8Vs/WQ7H874nMdfakij5+rlPtcjO44ztMkE4krEMvqjAQSFWEk6nUzv2sM4e+gcYz8eSN2nat/sl0/JB9QVr8q/Lq5EAWb9MI6YotEMajSenz7dgaZpvLKgK/c0vIOZLy1g03tbgEAPmvVvbwZg+/q9bFm7Hd2v8+zA5gxZ/DJpSZl8tnATnV5twtHdZ7jn7hLYErLwZXoCoesTmLwCg0sgHQZcLiNZTitp7jCSvaEkesNJ9FlI012k6Qak+3tkWlvkxWZgX4idCM75Ukj0WUj2hXHBE0aKO5RMjxWXy4zuMCCcAuEAHBLd5sed6iJUGCgQG8aBH45w9/3l8Dvd/PzVPjoOf5IyVYsBcP5YIlNfmEfFe8vSfebzudsn8dQFRjafQmSBcF77biwRseGBSbYbjcOWbmfGt2O5t8ndN/11U/IHFfLKTRFVMJI5P06gbPWSTHh2Jju+2oPRZGTsJ4O586EqvNblzZy9fFg+djVjn36N37YdpWDJArQf3orCZQpS+d4yPNWzEV++vZms+AsULRJO/P5zWI1GYkOC0Ow60u5Hc4DBKRAOge4043BZSHUFk+QKJ8kTQYIvgrPeIFJ0yWlfBol6EKkylLM+J4m+VBJ9IST7I0n0RpHkjiDZGUKW04rHZcLgMmJ0CoweMHshBCORVgu2CzYKhlsRwAt9GzF/+IeUu7M4T3YNDNM47S7GPTMLs8XI6NX9c5uKZaZmMbD+WHweHxPWDSOmcBSpCWn0rTOCjOQsJq0fTuX7yufhK6fc7lTIKzdNRGw4k9YPp2TV4oxsPpUdX+3BGmxh1Jr+FCxVgOGPT+bkr2eYtmEEv/5wCEuQmdb9HscabMHr8XHqQDwZZy9g1H2smbOe9AMnyLqQQckiEWSey8Tih2BpCOzNu8Ds0hAODZfTTIbDSpIzjPPuSM65o4n3RnHKayHBF0K6LkjRBcl6OKe8QST6ooj3xHDeHcUFVzgZrhAcDjPSbgi8gbjA6ASDU+JLc+NLd1GxQiGO7DxNm+71WTHpU1x2NwPe6ITBoCGlZOZLCzjzWzxD3+1DbNFoAOyZdoY9Nom0pIzAdqlcjHNHE+j74Aiy02xM3zxaBbzyj6mQV26qsKhQZmweTalqxRnTajo/frKd8Ogwpn0zmuDwIEY1n4o12MxLU9pzYt+Z3CkFj+48wbsT1hAcFsTM78aiG81EFI6mQqlIju08TZlSsZhcOnqaB5ntx+AAzQGaXUO3GbA7rKTZQ0lwRBHviuasO5aT7hjOeKM57Q3htDeYeG8kpz2xnPEU4KwzivOOSJLtYWQ5rPgdZjS7Ac0OBgcYHDomlyTaasEsBRcOJ1C+WlHCgwzs+vY3Oo9sScmKgatQ1y3axPert/HCuGep2bg6AH6/n3FPv8aJvacZtbo/lWqV52JiOoMbj8ftcDN98xgq1CibZ6+Tkn+okFduupCIEKZ8PYIy1Usx/pmZ7N60nwLFYpjwxVAc2U4GNhxHnZa16PdWV4QQHPzpCLO6LeSOOpXpPvMFKtxdiravPMr5Y4kULV2AUuUKknHmItLlI8piIcgT2Ms2OsBgB81mwG83k+2wkmIP5Zw9ijPOGE65YjnpjuOEO4bjrmhOuGM57YrltDOas84oEh1hZDiseOxmNLsRs0PD4hJYPRDs17B6wZ6YTfG4MNwODx26P8JbI9Zw5wMVaNb5IQAObD3Mm32XUvPR6jwzoDkQmDhlbvdF7P5mP33mvUjtJ2qQmpDGgEdGk3Uxm0nrh1Oxpgp45cZQIa/kicgCEUz5ajglKhdlTMtp7P32AGWrl2LC50NIOpXM2NYzcht1bf9qL827P8qTvZoAgflhi5SKIaZwFFu/3M9zPR/BbXNRMi4cR2I2ZjcE+wRWN1gcArNDQ7MZ8dosZNqDuGAPI94WxWl7LCftsZxwxHLSGcsJRwynHLGcsUeTZAsnzRaM02ZBZhsxZIPRDkYHaHaJnu5F2HyUKh7DiT1nad+jPu9P+wKz1cig+Z3QNI2si9lM7fgGBYrHMOzdPrkXPL034SO+XLyJtkNb0vTFBqRfyGD4Y5NIPZ/G5PXDKX9PmTx7XZT8R4W8kmdCIkKYtH44cSViGdNqOod+OUa1OpUZsKQH+749yOT2c/D7/Zw+GE9IzhR5p/af5etl37Nh+ff0eb0T1lAz8/oup/Z9JTm18wSlSsSgZXuR6V4MdonJCSYHGG0CkW3Ek20hMzuIZFsoCbYIztqiOZIVzsnsKE5mR3LOFklSdhjptiBc2RbIMmHM1jDaBAY7GOw6RrufKJOZIE0j5VgSd91XBntyJkf2nKbnlLbEFIrE7/Mzvs1s0hIzGLqiT+4EHt+v/pkVY1fToENdOk1oiyPLwYD6Y0g4nsTYjwdR7U9dORXlRrh88khFuYlii0Qzcd0wBjYYy+BG45jx7RgaPfcQWanZLOi/nPcnrqV1v8d5o/dSvl35I1EFI9A0jUFLe3Ji32lkto2L6U5O7TqONTuTxF/PElayAEYBmVl+DGgYhYYUIBCBHjZ+sPk0fF4T2QYzmiYxCIlfCjw+A36PEd1lBpvAaAsEvNEOZieE+DSC0XAk2SgRE0qWQaNR87uY0X0JTTrUoV6LGkBghqd93x1kwOJuVKpVDoDfth1lSoc5VH2wIq/M74rL4Wb8s7M4dzSRqRtGctcj1fLypVDyKfH7ga2bqWbNmnLnzp03/XGVW1dqQhqvPDAcgBnfjqFQqTimvvA6m97dQr+3XqbcPWXwOD3EFoumUKk4Ppj6CWvnfMmQd3qx96fjrJm3kXLlC3D+QjbBJQuT7PJiKRRKpubDGQzeUA1fGPhDwBss8Vt1jEHgM3gRWmBeVikFwmdAdwuMbiPYJEa7wOwAixOCHGDI8CPS3ZQrFs3ZvWd5ZeyTrJjwMaGRwcz9egiWIDM/frqDsa1fo2H7ugxc2gMhBIknL9C37kgsQWbm7ZhCcFgQI5pNZvfGX+m3qHvuhOaK8leEELuklDWvveQf1HCNckuILRLNmLUDcWY7ebXeKC4mpjNwaU9qNK7OG73fxu1wU61OJQqViuOz+RvYtWEfi/bNoEbDO+k4/Emq16nAiQPxPPVyAzKTsyhWIAxnYjYhfo0Qj4bJpmPMBpMNzNkCY6aGIcuIyDQjs8zomWZkpglDlglDpgFzlsCULTDbweoSWJ0g07yY7D5KFonizN6ztHjuATa+8wP2LCdD5nfGEmTmwpkUZnSZT8WaZXklZ4Ynp83JyOZT8Lq9jP1kEMHhQczutpCdX+/jlfldVcAr/yoV8soto/w9ZZi2aTS2DDujmk/BkeVk2HuvUKh0HCOemMy5Y4kAJJ9Noe3QlkQWCPSENxg0IkNNCClZv2ILnfs3IeFwIuVLFcB7wY7FJQn1GjBlB4LeYs8J8GwNc5aGMfOPD4tNw5wtMNkCbwjBLg2LQyIvegmTBiKDLKQeSaTq3SUxer38tuMkfWd2oHTVYjiynYxuNQOAoe/2xhJkxuf1MbHtbOIPn2fEyn6UrlaC5aNW8dWSzbQf/hSPvdQwz7a38t+gQl65pZS7qzQjP+zP6QPxjGoxFWuolSlfj8BoNjKy2WSSz6aw77vf8HsDjcWSz6bS+4HhXExIZ9qGEWSlO9ixbjdPtLmPU7tPU6ZEDL5kJ2anTpjPgDFLJ8ipYcqGYEfgsylLYMoSmLMEZpvAYhc5bwYSqxvkRS9RmhGTV0dkOYmKDuWRR6vy0ZsbadKhDg+3uhcpJW/0XsLpA2cZ9l4fipQNNFh7e+j7/LJuN73eeJF7Gt7Jqmmf8sHkj2nSuT4dx7fJy02t/EeokFduObWa3s2ApT05sPUw41rPIDIugrEfDyT5bCoT282hw8inWDBwBXN6LGZSh7lUf6gK0zeOpOp95egxqQ17vj+ENz2LBxtV5eyes5QrFYs/2UGIRxBkB2O2nyAHWJ0a5iwwZwc+TNlgzJaEp9fx/gAAEetJREFUuQ0YMv1E6Eb8F9zEGExIu5dwTeBxeOjS/1EWj15NtfvL03NKIKg/eX0937y3hXbDW3Hvo3cBsG7hRtbM/Jxm3RrTrFtjtn2xi8VD3uXhNg/Sd0HXvNzEyn+IOrtGuSXVb1sHR5aTOd0XsqD/O/SZ9yIDlvRkUrvZfLtyK8Pe7U1IeDAel5eSVQINwKSUNOnwIMnnLvLBrPV0GPgEjtpl2bf9JJVrleXQ0UTiSkZxIcVOTMEQ7BlejF6d3Lk3JODXCQ0NQtglRk0nCAPC6aNAkJnU+DT6T3qKBUM+IDw6lGELX8RoMvDrD7+xYMAKHnzyXjqMeAqAgz8dYW6PRdR67G56zOnEro37GP/Ma5SpXpKBS3pgMKoJP5SbQ4W8cst64uVGnDuawEezviA8OpSO49uQcCKJZSNXUrRsYZ4b/fQly/8+69Jzg5uReCaVd6d/wcsTAsvs2XacmvUrs3PPWeIKhuFzSvRMNyYjFCwcgWbQcDu9ZKTZcGbaKVskmvjfEilQKBKR7eLiuTQGTnmaFZM+xu30MHnNIKLiwkk4kcT4NrMpUq4gg5b2RNM04o+cZ0yr6cSViGX4B/1IS0xnwrOzKFahCFM3jMxtTqYoN4MKeeWW9tK0DtgzHbw38SOKli9Mu2GtiD9ynhXjVlOwVAEav/DwZfcRQtB/7gs4sl0sHLmGbpOeITwymB++2k/12mW5aHdz7kw6JUpGE104giy7C3SICgmhRGgopw8nce5gIvfVLsuJX07i9fgYMac9KyZ9QtKZVCas7E3JioVx2l2MfHIaul9n/CeDCAq14nK4Gdl8KgCT1g/HaXMxuPF4dL/O6I8GEFkg4iZvQeW/ToW8ckszGAz0XdCVxJMXmPXyW4THhNJ3wcukJ2XwWpc3KVymIHfU/eMq0aw0G+ePJVL5vvIMX/wSE7ssZP7QVXQe1ZLyA5uydObXhEcG0aRxNY6fSuHQtlOXPF5wsJl77imJVZdsWbePYqUL0HdcSxYOX8XJA/GMXNqNOx+ogM/rY3KH1zl3JJEpXw+nWIUi+Lw+pj7/OuePJTLtm1EULlOQwY3HkxJ/kSlfj8g9GKsoN5O6GEq5LWSlZTOk8XjOH0ti1pbxFCodR7e7B+J2epixeTTFKxYFYPRTM/j1+9+Y+tVwKtQsi9fjY0bPZfzw2S5avPQIDz9dmwWTvuDIr/HExIVzZ+2yFCkdi47A6/aScCKZ7d8dRtclT7S9j8Yt7mbyS4tJPp/GsIUvUrtJoIvknB6LWbfoG3rN7Uzz7o0BWNB/OR/N+oIeszvRolcTpj7/Opvf38qQFX1o0L5unm07Jf/4OxdDqZBXbhvJ8am88sBwdF0y8/uxeN0+BtYfgzXUyrwdUwiPDiP5bCoDGo7DkeVk5rdjKFG5KH6/zqLRa/h00bdUqlGa/q+/QPzpi2z8eBd7t53AaXfnPkZMwXDur1+Flh3rcHj7Cd4cuhKjycDIZd2oWivQGXLNrHUsHLSCZwc2p8ukdkBgAu5J7WbTrFtj+rz5EiunfsLbQ9+j4/g2tB/+VJ5sLyX/+TshH5id/iZ/1KhRQyrK33HqwFnZMqaj7HpXf5l5MUv+tu2obGppI199eJR02JxSSinPH0+UzxTtKtuW7C7jj5zPve/3n+yUrcr2lU8U6yUXj1srk8+nSa/HJ88cvyBPHEqQSefSpK7r8tDOk3Jkuzdkk7husm/TqTLpTGruOr58e5NsZHxWjnt2pvT5/FJKKfdvPSQbG5+R/R4aKZ12l/x62beyoWgtJ7SZKXVdv7kbSMnXgJ3y/5m3ak9eue3s+GoPo5+cRvHKRZm9dQI/f7aTqc/N5f4W9zLyw1cxGAyc2n+WQY9OwGQxMuXLQEtjgItJGSyb9CnfrNqGEIJSlYtQvU5FgkKtZF7MZv9Px4g/lkRoZDBt+jblya71c1sEb3p/K9M6zqNGozsZ+/FATGYjSaeT6VtnBJZgC2/umELCiQv0rTOCKg9UZOK6YZgtprzcVEo+o4ZrlP+MX77czchmU2j4XD0GLOnBJ3PXM//VZTTpXJ9XF3VDCMHJX88wpOkk/F4fY9YO4I4/tfFNOJ3CD5/sZNd3hzi65zQel5fgMCtV7i1DzQbVaNz2foJCrLnLf73sO2Z2fYs7H6rC+E8HYQ22YMuw07v2UDJTsnjtu7EITePVeiMJCgtizk8TiS0SnRebRsnH1HCN8p+yfPQq2VC0lq+9OF/qui6XjvhANhSt5fLRq3KXSTiRJDtV6SsfC24vv17+3RXX4/f7rzqs4vP55ZKRK2Uj47NycJMJ0ml3SSmldDlcckCDMbKJ+Vm57/uD0mFzym73DJStC3aRCSeTbvyTVRT594Zr1CmUym3rudFP4/X4WDnlY0pVLc4L454l+VwqK8atRvfrdBzfhsJlCjJ7y3gmtJnFjC7z2ffdQXrM7khIeHDuejTtyt09Us+nMePF+ez+Zj+PdnyYPvNexGQ2out6YF3fHmTgsp5UqlWOIY9O4OS+04xeO5DCpQverE2gKNf2/31XuBEfak9euVH8fr8c3XKqbChay/Vvb5I+n0/O6PKmbChay9WvfZa7nM/rk8vHfCgfNbeRTxfpKj9fsEG6ne4rrtNpd8mPZq+TLaI6ysdDO8gvFm7M3dP3+XxyasfXZUPRWn4wea2UUso3+y6VDUVruen9Lf/+E1b+01AHXpX/Iq/Hy4hmU9i7+QDD3u9LnVa1mNh2NlvWbKPzxHa0Hdoyd9kjO07w1qAVHNh6mLDoUB5scS9l7ixJSEQwTpuLE3tP89NnO8hMzeaehnfQ540uuRcxSSl5a8A7fDTrC54f8wwdRrZm5ZRPWDL8fVr0bEKv17vk1SZQ/iPUgVflP8uR7WTYYxM59PNRpmwYyZ31qjC90zw2vbeFFj2b0H12Rww5nciklOz99gBfLt7E7m/2k51uz11PWHQod9arTKtXHqPag5Vy++FIKVnw6nLWzlmXG+hb1v7CuNYzqN+uDoOW98pdv6L8W276gVfgaeAgoAM1r/d+arhG+TfYsxyyc5VXZLOwDnLXxn3S5/PJ+a8ukw1Fazn0sYnSlmG77D66rsvUhDR5/niiTD1/8YoHYF0Ol5zQZqZsKFrLN/q8Lf1+v9y5Ya9samkje98/9KrDPopyo/E3hmv+aT/5A0Ar4Id/uB5F+ceCw4KYumEkhUrHMbL5FPb/cIhur71A3wVd2b3xV3rWGsqx3ScvuY8QgpjCURQpW4iYItG5e+6/Sz6bQr96o/j+w595cUoHeszuxOFfjjH+mZkUq1iECV8MVV0llVvaPwp5KeUhKeWRG1WMovxTsUVjcpuDjWw2hR1f7eHxro2Y9s0o3A43fe4fxrKRK3FkO/9yPbqus2H5d/SoOZjzxxIZ+8kgnh3UgmO7TzK48XgiCoQz4YuhhEeH3aRnpih/j5oZSsl3IgtEMH3zGIpVLMKYVtP5btWP3FmvCm/tnUG9p+/nvYkf8UL53rwz5kMST174fegRCIztb35/C71rD2N6p3kUKVeIN36ZzP3NanJ8zylGNptCeEwYM78fR1zx2Dx8lopyfa554FUI8Q1wpR6pw6WUn+Ys8x0wQEp51aOpQoiuQFeAEiVK1Dhz5szfrVlRrktmahajW07j4I9HaDesFc+PfQaDwcDh7cd4Z+xqdqzfA0BoZAhRBSPwuLxcOJMCQNHyhWk7tCWNnn8ITdPYuSEws1NIRDCTvxpBycrF8vKpKf9ReXZ2zfWE/J+ps2uUm8Xr8TK3+yK+Wvotd9WvxuDlvYgtGgNAwokkdm3Yx4l9Z7Bn2hGaoGSV4lSuXYG7HqmKpmlIKVkz8wsWD15BiSrFmLhumNqDV/KMCnlFuYqvlmxmXp8lmKwmeszuRIP2dS87yPq/MlIymddnCd+t+om6T93HwKU9CQoNukkVK8rlbnrICyFaAq8DBYAMYK+U8tFr3U+FvJIX4o+cZ9oLb3B4+3Eq3VeedsNacd/j91zW1iA73caXizbxweS1OG0uOo1vw7ODn7zmm4Ki/NvUxVCKcg2Bs2a+550xq0iJv0iBYjFUuLcsUXERSAmJJ5P47aejuBxu7m1yFy+/9oIaf1duGSrkFeU6+bw+flj9Mz99toMTe09jS7cjpaRQ6TjK31OGx7s2otzdpfO6TEW5hAp5RVGUfOzvhLw6T15RFCUfUyGvKIqSj6mQVxRFycdUyCuKouRjKuQVRVHyMRXyiqIo+ZgKeUVRlHxMhbyiKEo+pkJeURQlH1MhryiKko+pkFcURcnHVMgriqLkYyrkFUVR8jEV8oqiKPmYCnlFUZR8TIW8oihKPqZCXlEUJR9TIa8oipKPqZBXFEXJx1TIK4qi5GMq5BVFUfIxFfKKoij5mAp5RVGUfEyFvKIoSj6mQl5RFCUfUyGvKIqSj/2jkBdCTBdCHBZC/CqE+FgIEXmjClMURVH+uX+6J78RqCalvBM4Cgz95yUpiqIoN8o/Cnkp5QYppS/n221AsX9ekqIoinKjGG/gujoDq652oxCiK9A151ubEOLIDXzs38UCqf/Cev9NquZ/3+1WL9x+Nd9u9cLtWXPF/+8dhJTyrxcQ4hug0BVuGi6l/DRnmeFATaCVvNYK/0VCiJ1Sypp59fh/h6r533e71Qu3X823W73w36n5mnvyUsqG13jQF4AngAZ5GfCKoijK5f7RcI0QogkwGHhISum4MSUpiqIoN8o/PbvmDSAM2CiE2CuEWHADavonFubx4/8dquZ/3+1WL9x+Nd9u9cJ/pOZrjskriqIoty91xauiKEo+pkJeURQlH8tXIS+EMAgh9gghvsjrWq5FCBEphFiT0xbikBDi/ryu6VqEEP2EEAeFEAeEEB8IIax5XdP/EkIsEUIkCyEO/Oln0UKIjUKIYzmfo/Kyxv91lZpv2ZYhV6r3T7cNEEJIIURsXtR2NVerWQjRWwhxJOf3elpe1XclV/m9uEsIsS3nGOhOIUSta60nX4U88ApwKK+LuE5zgK+klJWA6tzidQshigJ9gJpSymqAAWiTt1Vd0TKgyf/8bAiwSUpZHtiU8/2tZBmX13wrtwxZxuX1IoQoDjQCzt7sgq7DMv6nZiHEI0AL4E4pZVVgRh7U9VeWcfl2ngaMlVLeBYzK+f4v5ZuQF0IUAx4HFud1LdcihAgH6gFvA0gpPVLKjLyt6roYgSAhhBEIBhLyuJ7LSCl/ANL+58ctgOU5Xy8HnrypRV3DlWq+lVuGXGUbA8wCBgG33NkcV6m5OzBFSunOWSb5phf2F65SswTCc76O4Dr+BvNNyAOzCfyC6XldyHUoA6QAS3OGlxYLIULyuqi/IqU8T2BP5yyQCGRKKTfkbVXXraCUMhEg53NcHtfz/9UZWJ/XRfwVIURz4LyUcl9e1/L/UAGoK4T4RQjxvRDi3rwu6Dr0BaYLIeIJ/D1e8z+8fBHyQogngGQp5a68ruU6GYF7gPlSyrsBO7feEMIlcsaxWwClgSJAiBCiQ95Wlf/ltAzxAe/ldS1XI4QIBoYTGD64nRiBKKA2MBD4UAgh8raka+oO9JNSFgf6kTMa8FfyRcgDDwLNhRCngZVAfSHEu3lb0l86B5yTUv6S8/0aAqF/K2sInJJSpkgpvcBa4IE8rul6XRBCFAbI+XxL/Vt+NX9qGdL+Fm8ZUpbAm/++nL/BYsBuIcSVel7dSs4Ba2XAdgKjALfUAeMreIHA3x7AauC/ceBVSjlUSllMSlmKwMHAzVLKW3YvU0qZBMQLIX7vKNcA+C0PS7oeZ4HaQojgnL2dBtziB4v/5DMCfxzkfP40D2u5Ln9qGdL8Vm8ZIqXcL6WMk1KWyvkbPAfck/N7fiv7BKgPIISoAJi59btSJgAP5XxdHzh2zXtIKfPVB/Aw8EVe13Eddd4F7AR+JfDLFpXXNV1HzWOBw8ABYAVgyeuarlDjBwSOGXgJhE0XIIbAWTXHcj5H53Wd11HzcSAe2JvzsSCv6/yrev/n9tNAbF7XeR3b2Ay8m/P7vBuon9d1XkfNdYBdwD7gF6DGtdaj2hooiqLkY/liuEZRFEW5MhXyiqIo+ZgKeUVRlHxMhbyiKEo+pkJeURQlH1MhryiKko+pkFcURcnH/g8Prv21ECo93gAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"u=np.linspace(3,18,200)\n",
"v=np.linspace(-2,4,200)\n",
"X,Y=np.meshgrid(u,np.exp(v))\n",
"#print(np.exp(Y))\n",
"Z1=L1(X,Y,[10,10,11,9,12])\n",
"lev=[.0001,.001,.01]+[x for x in np.arange(.05,.95,.05)]\n",
"CS=plt.contour(X,np.log(Y),np.exp(Z1-np.max(Z1)),levels=lev)\n",
"junk=plt.clabel(CS,[.0001,.001,.01],fmt=\"%1.4f\")\n"
]
},
{
"cell_type": "code",
<<<<<<< HEAD
"execution_count": 7,
=======
"execution_count": 11,
>>>>>>> f8632e30fe8b89896ed0fbf258653883205028de
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/home/jet08013/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:15: RuntimeWarning: divide by zero encountered in log\n",
" from ipykernel import kernelapp as app\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD8CAYAAACSCdTiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xd4FFUXwOHfZDebTe+NdHpvofdeBEQRBQUpiiCooCgIIqJ0gU9EQOkKAipNmvTeWyChBwhppPe+2TbfHxuCSIdAINz3efJks2Xm7iY5c+fOvedIsiwjCIIglExmxd0AQRAE4ekRQV4QBKEEE0FeEAShBBNBXhAEoQQTQV4QBKEEE0FeEAShBCuyIC9JkkKSpDOSJG0uqm0KgiAIT6Yoe/LDgEtFuD1BEAThCRVJkJckyRvoBCwqiu0JgiAIRUNZRNv5ERgJ2N7rCZIkDQQGAlhbWwdWrFixiHYtCILwcggKCkqWZdn1UV7zxEFekqTOQKIsy0GSJLW41/NkWV4ALACoU6eOfOrUqSfdtSAIwktFkqTIR31NUQzXNAZelSQpAvgTaCVJ0vIi2K4gCILwhJ44yMuyPFqWZW9Zlv2BnsAeWZZ7P3HLBEEQhCcm5skLgiCUYEV14RUAWZb3AfuKcpuCIAjC4xM9eUEQhBJMBHlBEIQSTAR5QRCEEkwEeUEQhBJMBHlBEIQSTAR5QRCEEkwEeUEQhBJMBHlBEIQSTAR5QRCEEkwEeUEQhBJMBHlBEIQSTAR5QRCEEkwEeUEQhBJMBHlBEIQSTAR5QRCEEkwEeUEQhBJMBHlBEIQSTAR5QRCEEkwEeUEQhBJMBHlBEIQSTAR5QRCEEkwEeUEQhBJMBHlBEIQSTAR5QRCEEkwEeUEQhBJMBHlBEIQSTAR5QRCEEuyJg7wkSWpJkk5IkhQiSdIFSZK+K4qGCYIgCE9OWQTbyAdaybKcLUmSOXBIkqStsiwfK4JtC4IgCE/giYO8LMsykF3wo3nBl/yk2xUEQRCeXJGMyUuSpJAkKRhIBHbKsny8KLYrCIIgPJkiCfKyLBtkWa4JeAP1JEmq+t/nSJI0UJKkU5IknUpKSiqK3QqCIAgPUKSza2RZTgf2AR3u8tgCWZbryLJcx9XVtSh3KwiCINxDUcyucZUkyaHgtiXQBrj8pNsVBEEQnlxRzK7xBJZKkqTAdNBYJcvy5iLYriAIgvCEimJ2zVmgVhG0RRAEQShiYsWrIAhCCSaCvCAIQgkmgrwgCEIJJoK8IAhCCSaCvCAIQgkmgrwgCEIJJoK8IAhCCSaCvCAIQgkmgrwgCEIJJoK8IAhCCSaCvCAIQgkmgrwgCEIJJoK8IAhCCSaCvCAIQgkmgrwgCEIJJoK8IAhCCSaCvCAIQglWFOX/hJeYXqcnP0+LVqNDpzF9N31p0eXr0eXr0OXrCu/X5evQaU3367V603O0ptt6rR69zmD60urR6/UY9EYMegMGncH0XW/AaDBi0BsLvhswGmWMBtPPRqMR2SgjG+XC20ajjCzLIJtuAyDL939jkoQkSTdvIpmZYWYmgSRhZiYhmUmYKcwwMzMr+F7wc8GXQqlAobz5/T+3zRUolQqUKiVKcyVKlRJzlen7zdvmFuaYW5i+q9TmBT/feZ+FpQqV2hyVWoXKUoW5SlnYbkEAEeRfGrIso8nRkJOZR25mHnlZeQW3c8nNzCM3K4+8LA25WXlosjXkZWvIy9GgydGgyclHk5NPfq7pS5OrNd3O02I0GIukfeYW5ijNFSjNFSjMTYFOYW4KikpzBWYKM9NjSkVhEDVTmKFUKW8FVoVZYfCVJFMgVijMCgI2hffdCt53D4ayXHBQAJBNPxuNMrLRiCyDbDQWHlhuHkwKDzIGIwaD0fTZ6A0Ybh6I9Eb0OtNBS6/TY/jXwUxXcICTH3TgeQiSJGFhqcLCSoWFlQUWVhaoC26rrS2wtFGjtjZ9WdmqsbSxxNJWjZWtJVZ2lljZWWFlZ4m1nSWWtpZY2ZoeVygUT9w2oXiIIP+C0eTmk5mSRWZyFhnJmWSmZJOZkkVWarbpKz2b7LQcstNzyEnPJSfD9JWbmXurF3sfZgozLG1M//RqGzVqa1NwsHexRW3tgspShdpKjdrK1HO0sLQw9SQtVVhYqky9TEsVqn/1OJUqJSq1CnOLW9+VN3urKiUKpUL0PgGD3vCfsxwd2oKzocKzoH/dLjw70mgLz6a0eabb+XlaNLkatHnawoN0Vmo2SdEpaHLyTQfxbA26fN1Dtc3SRo21vRU2DtZY2Vth42CFraMNNg7W2DrZFH7Zu9hh52yDvasdds62WNlait9tMRNBvpjJskx2eg4psWmkxKaSGp9OWnw6aQkZpCdlkJ6YQXpipimgJ2eRn6e957YsbdTYOFpj62iDtYMV7v6u2DhYY21nhbW9qYdmaXurl2Ztb1XYgzP12tSo1CrxT1lMbg7nqK0sntk+9Tr9rbO4zFxyszSFnYLcm/dl5hXel5OZS1ZaDmnx6URdiiE7LYecjNx7noUozRXYOdti52KLg6sdDm72OLja4+Buj5OHA04eDjiXcsK5lCP2rnaYmYnLhEVNBPmnSJZlslKziY9IJCEymaSoZBKjk0mOSSE5JpWU2DRS49LQau7sTVlYqnB0t8fe1Q7nUo4EVPfFwcXUO7J1MvWUbvaabt6nNBe/TuHRKM2V2DnZYudk+9jbMBgM5KTnms4wU7LISDZ9Zd32s6mjcvV0OGkJ6eRm5t2xHYVSgaO7Pc6lHHHxcsLFyxkXb2fcfF1w93PF3c8FJ09HcSB4RCIqPCGj0UhiVDLRobHEXosn9lo8ceEJxIcnkhCRRG7W7X/MFpYqXLydcfFyonLD8jh5OOJcyrGwN+Pk6YiThwOWNmrRoxZeCApFQW/d+eEPFFqNltT4dFLj0ws7O8kxqaTGp5ESm0bM1XhC9l0kOz3ntteZW5jj7ueCu78bngFulCrjQamyHniX96RUWQ/MVeZF/fZeeCLIPySj0UhsWALh56KIOB9F1KUbRF2O4UZo7G09cbWVBR6l3fAIcKNG8yp4lnbHzc8FD3833HxdsHWyEcFbeOmp1Co8/N3w8He77/Nys/JIik4mPiKJxMgk4sMTiS/4fuXkNbLSbh0EzBRmeJZ2x7eSF74VvfCv6ktANV98K3m91MFfBPm7kGWZG1diuXziGldOhnH1zHXCgiPQ5OQDphkM7v6u+Fbyonbr6vhUKIV3hVJ4lfPEycNBBHFBKCJWtpb4VfbBr7LPXR/PSssm9lo80aGx3AiNJeryDaIvx3Jy6xn0OgNgui7gV8WHcrUCKBdYhgr1ylKmht9LM7wpFcW0rUdVp04d+dSpU898v/diMBi4HhJJ8J7znD1wkQtHQslKzQZMPfPSNf0pVyuAMjX9KV3dD9/K3lhaq4u51YIg3Itepyfmahzh56IIC47g6plwrp2+TkZyFgAqtTkV65ejWtNK1GhRhSqNKqBSq4q51Q8mSVKQLMt1Huk1L2uQz8nI4fiWMxzbfIqgHWfJTDH98r3KeVK1cUWqNK5Axfrl8K3kJeYIC0IJIMsyCZFJhJ64xoUjoZw/fJmwM+EYjTIqtTnVm1emfqdAGnWti5uPS3E3965EkH8AnVbHkQ2n2L3iQOHpnIObPXU71CSwbQ1qtqqKs6fjM2+XIAjFIyczl/MHLxG08ywnt53hxpU4ACo3LE/Lt5vQ6p0mTzTzqKgVS5CXJMkHWAZ4AEZggSzLs+73mmcd5NMS0lk/eytbFu4iPSkT51KOtOjRmGbdG1CxfjkxJUsQBACiQ2M4uPY4+1cf4XpIJCq1Oa3ebsIbw7vgX+Xu1wWepeIK8p6ApyzLpyVJsgWCgNdkWb54r9c8qyCfl53HX99vYO3MzeTnaWnQJZDOg9oR2K66GIIRBOG+wkIi2PTLDnYvP0B+npbWvZvSf0JP3Hxdi61Nz8VwjSRJG4A5sizvvNdznkWQP3/4Mt+/+xPxEUm06NmYvt++hXf5Uk91n4IglDyZKVn8NW0D62dvQWmuZMis/rTr26JYZtEVe5CXJMkfOABUlWU58z+PDQQGAvj6+gZGRkYW2X7/a++fh/m+z2zcfF0YufRjqjau+NT2JQjCyyE+IpHp/edydv9F3vy8Cx9Me/eZB/rHCfJFNhgtSZINsBb49L8BHkCW5QWyLNeRZbmOq+vTO905e+AiU3vPonKj8vxyepoI8IIgFAkPfzem7fqGV4e0Z/X/NrFq+sbibtJDKZIgL0mSOaYAv0KW5XVFsc3HYTQa+eGDebj7uzFx02is7ayKqymCIJRACoWCj2e/T9PuDfht7B8kRCYVd5Me6ImXfEmm85XFwCVZln948iY9vvOHLhNzNY7RK4ZhZWtZnE0RHpFOq/tXWuSCnPc389pnm/La5xfksdfk5qO9mVo3X1uYerew2EhBjvabedz/XWzE+O+CIjJ3Fg8pyD1vZlaQm/5fxUD+XfRDoVTcVujjZmplc7U5KgtT2mVTWmYLLKxUqK3VWNqYvtQ2aqztTNk/re1NGUJF9s8XhyRJDJreh0Nrj7Nr+QF6jXmjuJt0X0Wxrrcx8C5wTpKk4IL7vpJleUsRbPuRRF6IBqB688rPetfCv+Tn5RemTE5PzCQ9KZPM5EwykrPITDXlvs9OzyE7Lcd0Oy0HTW7+Q23bzEwqzF2vUqsKguqtiklKldIURAuKjdwsOHKzyIhkJqEwu1lI5FYxEaAwXe7NA4GMfHsxEL3hVtEPvRG91lQVKzczt7AK1r/zumty8x+6qIq5SmlKE+1kytF+87atoylHu72LLfYFqXodC9L0WtlZiQNDMXH3c8UjwI3Ii9HF3ZQHeuIgL8vyIeC5+EsztzAlIdLeJ+e68GTysvOID08kMSqZxKhkkm6kkBxrSpucEpNKckwqORm5d32thaWqMC2yjaM1nqXdKR9Y5lYOfPtbee+t7CwLC5dYFhYvUb9w5e30Or2paEduPppsU+WtvGzNrepcGblkp+eSk5FTUPTFdPC7ma89MyXrrml5wfR5OpdyxNnLCRcvJ5w9nXD1dsbVxxk3P1c8/F2xc7Z9oT6vF4lWo0Wpev7z3zz/LXwE5euUASBoRwilBnsUc2teXHk5Gm6ExhIdGkvMlThir8cTczWOuLAE0pNuv6ZupjDD2dOULtmnohc1W1bFuZQTjh4OOLqbep0OrnbYudg902IYzwuluRIbByU2DtaPvQ2dVkdmSjYZSZmkJ2aQlpBRkKI3lZQ408E19MQ1kmNS76hNYGVriUdpN7zKeeJV1hOvsh54VyiFb0WvR0oNLNwu/HwUKbFplA8sU9xNeaASFeQDqvlSuoYfa2ZupsP7rV7q9KIPQ5Zl4sMTuXYmnLCQCK6fjST8XBTx4Ym3Pc/V2xmvch406loXz9LueAS44ebnipuPM44eDmJh2VNmrjI3HUgfkHLjZpGam2dZ8eGJxIUnEBsWz/WQSI6sP4lBbyh8vqO7Pf5VfSld3Y8yNfwpU9Mf30peL012xiexcvI61FYWNH+rUXE35YFK1G9TkiTen9yLMZ0ms/SbvxgwtXdxN+m5kpej4dKxq1w4fJmLR0MJPXErH7eZmYR3hVJUqFuG9v1a4lfZu7AQg4Xly9cDfxFJklRYvKNsrYA7HjfoDSREJhEdGkvUpRgiL0QTfj6KjT9vL6z1qlKbU7Z2aao0LE+VgkR9Dq72z/qtPNeObQ5i35+H6fX1Gzi6Pf+fTYlMUPbjoPn8s3AXY/74lBY9Gj+1/TzvjEYjV06FcXJrMEG7Qrh8/BoGvQFJkvCr4k2l+uUpX6cM5WoH4F/VRwTzl5RBb+DG1TjCzoQTejKMyyeucjXoOjqtHgD/qj7UalWNuh1qUqNFlRciJe/TEnEhmmGNx+AR4MbsY1NQWTzb0YJiX/H6sJ52kNfm6xjZ5jtCT1zjmzVf0LDLI30mLzSDwUDI3gscWH2UwxtOkp6YgSRJlAssTe3W1ajWrDKVG5Z/ojFioeTT5uu4GnSdcwcucmbvec4fvIRWo0NtZUFg+xo0696Qhq/WeanqKsRci+PzFuOQjTKzj00ulhw2Isj/S3Z6DqPaTyAsOIKRSz+hZc+S3aNPjEpiy8LdbPt1DymxaVjaqKnfqTYNOtehTvsa2LvYFXcThRdYfl4+Z/df5OjGUxzecJLUONPfWIsejek0qC0V6jz/FyCfRFhIBKM7TMRoMDJ99zgCqvkVSztEkP+P7PQcvnnte84duES/8T15Z0y3EjedLOJCNH9O/Zu9fx5GNsrU7ViTdn1b0qBzbTH8IjwVRqORcwcvsXPpfvavOoImN59qTSvxzpg3CGxbvcT9jx3ddIqpvX/C2t6KKdvG3LMU4bMggvxdaDVafvhgHrtXHKRJt/p8sWRIiUh3kJmSxZKvVrJl0W4srFR0GtiW14e+grtf8aVBFV4+OZm5bP91L2t+2ERSdAo1W1Xl49nv41fJu7ib9sQMBgMrJ67j9/GrKVs7gO/+Homrt3OxtkkE+XuQZZk1P2xm0ajleJZ2Z/SKYS/06eWpHSF8/+5PZKXl0PWjDvT6+g0x51koVtp8HVsW7mLZuL/Iy9bw/pRevPFZ5xe2V58Sl8a0fnM4vfMsbd5txrBfBj4X6zxEkH+AswcuMvXdn0iLT6fPtz14a8SrKJQv1hzv9XO28vOwX/Gr4s1XK4YV29igINxNWmIGswYv4PDfJ2jduykjlnz0wv2PHV5/gpkD56HJyWfIj/3pOKD1c3OwEkH+IWSmZvHjhws4uOYYlRuW5/PFQ/Ct6FUsbXlUWxfv5ocP5tGoa11Grxj2XPQsBOG/ZFlm5aR1/PbNn7Tt25wRSz56boLk/WSmZvHLZ7+x6/cDlK0VwKjlQ5+7YafHCfLIsvzMvwIDA+XitnvFAfl1p75yR/Xb8h9T/5Z1Wl1xN+m+ws9HyR3Vb8tftp8ga/O1xd0cQXig3775U24jdZd3LN1X3E25L6PRKB9Ye0x+y3OA3E75lvzr2D+e23gAnJIfMd6+dD35f0tLSOenjxZxaN1xSlf347MFg6hYr1xxN+uuvu02jeC9F/jtyk9iBaLwQjAYDHzWdCzJN1JZem32c5lmJDE6mblDl3Bkw0nK1PTni8VD7rpa+HlRrJWhXkSO7g6MW/MF364bQUZyJkMbjuGnIQvJSssu7qbdJistm6Obgug0sK0I8MILQ6FQ8PbobiTdSCF4z/nibs5t9Do9q/+3ifcrf0rQjhA++L43c09Mfa4D/OMqUblrHlfj1+pRs1VVln7zFxvmbOXAmqMMmNqbdv1aYGZW/MfBq0HXMRqM1Glfo7ib8lzT5uvITsshOz2HnIzcgmIj+aaCI3ladBod2nwduoLCIjfzxP/7bPbfxUFMxUDMUanNCwuAqK0tsLRWY2VvhbW9JbaONqitLV6IMefiENi2OmZmEheOhFK3Q63ibg4AIfsvMOfjxURciKZB50A++uk9PPzdirtZT40I8gWs7awY8mN/2vdvyU8fLeJ/A35h8/wdfDJnABXqli3WtmWmZAGmM4+XkU6rJyk6mcToFFOGxehkkm+Y0uymJ2SQnpRBelIW+Q9ZeOTfzBRmmJndLBrCHUH/YZirlNi52OLgequgh7OXE24+zrj5uuDm44Kbn8tLlQLgJpVahbWDNZnJWcXdFJJupLBgxDL2/XUEdz9Xvls/koZd6pT4A7QI8v9RpoY/Mw+MZ/eKgyz6cjkf1x9N277NeX9yrwemen1arAoWb2U/Z8NIRS0vW0PEeVNmxIgL0URdjiHmahxJ0SkYjbcHXntXO5w9HXHycMCnQinsXe2wcy6oquRgjZWdJZa2llhamwqOWFgWVJFSm5sqSJkr7nmWZjQaMeiN6PJNZQULKz0VFP/IzcojL0tjqm6Vbir2kZ5kqoCVnpBBxIVoUuPT76gK5ehuj1dZD3wqeOFX2ZuA6r4EVPXFwbXkppzQ6/TkZeVhbV+8CxB3LT/ArA8XYDQa6T22Oz2+fO2lmZ0mgvxdmJmZ0fbd5jTqWpc/Jq9j3Y//cHDNMXqOep3uwzs/83QB/lVNy6ivnLpO1SaVnum+nxaDwUjEuSjOH77MxWNXuXr6OjeuxBU+rra2wKeCF5UblqdUbw88C3LYu3o54erj/FQzIZqZmWGmMsNcpYTHrBVs0BtIiUsjKSqFxOhk4iMSiQ1LIOZqHIc3nGTrkj2Fz3Uu5Uj5wNJUrFeWqk0qUqFOmRKT6TEsOAK9zkBA9ae7nkOW5fv2yAOq+VK/c20GTO1doodm7ualnl3zsGKuxbFo1AoOrTuOq48zA6b0okXPxs90vP79Kp9i52zLzAMTntk+i1rSjRSObQ4iaNc5QvZdKCwT6FzK0ZTyuFaAqfBLdT/c/Vyei+shT0taQjrXz0URfi6aa2fCuRIUVniQM7cwp3KDctRuW536HWsRUM33hR1SWDjyd9bM3Mzq+EVFvir7nwU72TB3G9N2ffPSTEgQi6GespB9F5j3+VKunQmnYr2yDPpfX6o2rvhM9r1q+gYWfrmcOSemvlApGZJjUtnz52EOrjlG6KkwANz9XKjVqhrVm1WiWtNKuPm6vLBBrChlpmRx4UgoZw9cInjvBcJCIgDwCHCjabf6tOnVlIBqvsXbyEeQl51H74CPqNa0It+uG1mk2z666RTbluzhxpVYGnSuwwff935gb74kEEH+GTAajexctp9fv/6DlNg0mr3ZkPcnv0OpMk+3pmxOZi59ynxM6eq+TNs17rn+Y5ZlmTN7zvP37K2c3HoGo1GmfGBpmnSrT6NX6+BTodRz3f7nRUpcGsf/Oc2h9Sc4s/s8Br2BCnXK0OXDtrTo2fiZF6x4VMu+XcXv41cz68gkKjco/8Tby0rLRqFUYGVrSXxEIhZWFiDLDG34FZO3jsGnwouxcv1JiCD/DOXlaFgzYxOrpm9Ar9PT9eOOvDOmG3ZOTy9R2MaftzP740WMXjGMVm83eWr7eVyyLHNi6xl++2YVYSEROLjZ0/G9lrTr2wKvsqKw+pNIT8pk7x+H+GfRbqIuxeDk6UjPkV3pNLCN6drBcybqcgwf1vyCxt3qM2blp0+0LZ1Wx8xB84m5Eoebrwv9J76NZ2n3wo7C/C+WkRyTwpg/PsNoNJboYT6R1qAYJMemyjPemyu3NXtTft2pr7xqxkY5X/N00g7o9Xr54/qj5Ned+srxEYlPZR+P68bVOHlkuwlyW2UPuW+FofLWJXue2ufwMjMajXLQrrPy562+ldsqe8j9Kg2TT+8+W9zNuo02X2v6O3XuJ6fEpT7x9k5sOyNPfHumLMuyvGjUcnnusCXywXXHCh/Pz8uX+5T7WD61I/i21xmNxife9/OGx0hrUHIPec+Is6cjny8ewvzg6VSsX44FI5bxfqVh7Fl5EKPR+OANPAKFQsGo5UMxGIxMfudHdFpdkW7/ce1Yuo8Pa4/kStB1Pp7Vn0Xn/keH/i2f++GEF5EkSdRuXY3pu75h4sYvAfiy/SR++Xwpep2+mFtnMv/zZVw+cY1P5w3EyePxph3nZecVrlc4u+8CVjamNQZvfvEq/lV9ObU9hLSEdMA0F/+tL15l3ax/iA6N4c/v16PT6sSQYAER5ItIQDU/Jm8Zw/c7xmLtYM2U3j/xcb1RBO8t2uXcXmU9Gb5wMBePXuHHDxc88sKdoiTLMovH/MGMAfOo1KAcC0Nm8OqQ9ijNn7/hg5JGkiTqdazFvKBpdB3Snr9/2srXr35PXo6mWNu1/be9bJi7jTc+60yz7g0f+fVhIREMDhzJzEHzWfTlcgAav16P9KRMUuPTsHO2pWK9sljbWXJ8y5nC17Xt24KTW4MZ3uwb1NYWz2WenGLzqF3/ovgqScM1d2MwGOSdv++X3/H7UG4jdZfHdJ4sR1yIKtJ9LB33l9xG6i4vn7imSLf7SG34dpXcVtlD/nHwQlmv0xdbOwRZ3vbrXrm9qqc8quOkYsugeGLbGbm9eQ95RJvvHuvvITcrV57y7ix56+LdcnZ6tvxRvS/lfxbuki+fvCYvGbNSXjfrH1mWZVmTq5FXTd8gr5+zVZZlWY4LT5An9vxB/t+AX+Tc7LwifU/PG8RwzfPBzMyMNr2b8evlWXzwfW/OH7rMwOqf8+Og+aTGpxXJPt4d9yatezflt7F/sumX7UWyzUdx7J8glk9cS/u+LRg69/0XrjBESdO+XwuG/fwBQTvPsuy71c98/5dPXGV89xn4VfFm3NovHvrv4d9DmpY2liREJOFbyQtre2sGTO3NpWNXSIlNpWytAEJPXePq6etYWFqgtlETfjYSAFcfZ4b+/AHDF374UqaOeJAiOa+WJGkJ0BlIlGW5alFssyRQqVW8NaIr7fu3ZPmENWz6ZQd7/jhEj5Gv8cbwzk+0rFqSJL5YPISc9Fxmf7wYtY2atu82L8LW35tWo2X2J0sIqOrLJ3Pffy7HPvU6Aylx6aQkpJOWlEV6UiaZKdlkZ+SSm21KT6DL16HTGgoDjZkkoTBXYK5SoraywNLGAht7K+ycbLB3tsHRzQ5nDwecPRxQWz1/K1I7vt+Ki8eusGr6Rtr0bvbMiuFcOxPOqPYTcXCzZ/KWMQ9dQ3nTvB1s/Hkbzd9sRLnA0tR/pTaV6pcjKToFGkLNllU5u/8ikRdu0OzNBsRdT2DmoPn0G9+TvX8eomWPxsiyjEKhwNbR5im/yxdXkUyhlCSpGZANLHuYIF8SplA+jhtX41g0ajmH/z6Bi5cT/Sb0pM27zVAoHr8XrNVoGdN5Cmf3XeCLJR/Rts/TD/Rbl+xh5qAFfL99DLVaVXvq+7sfnVZPZGgcYWejCL8UQ1RoHDHhiSTHpN2R7wbAwtIcK1tL1FYWqCyUKFVKzBSmE1rZaESvM+Ws0eRpyc3SkJd99zFuR1c7PANc8S3ngX+lUpSp5kOZaj7F3pNMT8qkT7lPaP5mQz5f+OFT39/V09f5st0E1NYWzDww4aELyZ/edZal4/5i8Mx+JEQmM//zpcwLns6eFYdIT8ygbd/meJX1JPxcJLM/XsyEjV9ibW/NloW7uHziGqVr+PHaxx2f8rt7/hTrPHlJkvyBzSLIP9i5g5dYMMI0A6F0dT+zX7IzAAAgAElEQVQGTn+XwLaPn0ZYk5vPN12/J2TveeYFzyCg6tNdFTn6lckkRCax+PwPz7wXr9cZuHgijDMHLnH2yFWunY1CqzHNMlJbWeBT3gOv0m54+rvi5uWIs6cDjm72OLraYuto/cgzfgx6A1npuaQnZZKWlEVqQgaJN1JJiEom5noiUVfiyEzNAcDMTMK3QimqNSxLjaYVqdWsYuGskGdpxvu/cGj9CdYmLHqqw2ihJ68xqv1ErO2tmL57HJ6l3R/6tce3nCZoRwhDfuwPwJxPFmM0GHlzxKv8MflvqjWtRPMejVBZmPNl+wl0/6zzc5OquDg9TpB/ZtMgJEkaCAwE8PV9cZZmPw3VmlZi1pFJ7F91lCVjVjKq/UTqdqjJB9PefawArbayYPyGLzm49thTD/CyLHPp+FVa9Wz8zAK8Xmfg9L6L7F13khM7z5GbpcFMYUb5mn507teM8rX8KVvdF0//os93o1AqcHCxxcHFFv+75IaTZZnUhAyunY3makgkF09eZ8cfR9m0ZD9KcwU1mlSgaZfaNHm1NtaPmezsUQW2rc6OZfuJuHCDMjWeTmKwkH0X+Kbr99i72jF997j79uANBsMdZ6s56Tnk5+aj0+owV5nz/tReDAkcSZfB7WjduylbF+8mLCQCr3Ke5GbmUaam/1N5Hy8D0ZMvZtp8HRvmbGPlpLXkZubS4b1W9JvQ87nNHZ+blcdrTv0ZMLUXb33e5anuKy0xk02/7mfb8kOkJWZi62hNww41qNeuGjWbVnhmQfNR3TzbOLHrHIf/CSY+MhkLSxUtXq9Dt8Ft8C3v+VT3H3oyjE8ajeG7dV/QsMujLY58GEc2nmRSz5l4lnZn6vavcfFyvudzl4xZSVZqNrXbVKfpGw0K78/L0fBZ07EMntmPGs2rALBy8jpCT17ju79HkhiVxIY520iKSaH78C6UD3xx8jU9Tc91T164O5WFOW9+3oX2/VqwfMIaNv68nX1/HeGdMW/w+rBXinxBUV6OhsvHr6LV6Kj/Su1Hfr1ea1pwY27+9IYBcjLz+GPmFjYt2Y8uX0/dNlXp2Lsxga2qPJdL+P9Laa6geuPyVG9cnve/6Ubo6Qi2rzzC3rXH2b7yCM1eDaTvmK6U8n+48etHpVKb/mZ02qJfHLVl0W5mfTifcoGlmbh59H2zPy77dhXh56Po+H5r1vywiaQbKbTr2wIbB2ssrdW069OCNT9swq+yNw6u9jR+vR6JUcnkZefh5utK/0lvizUXRUB8gs8JO2dbhvzYny6D27Fg5O8sGrWcfxbsZOD0d2n8Wr0nHhoxGo1cPn6V+SN+x6usB1dOhXF651kGz+z3SNuxsjP1nrPSc56oPfdycGMQc0f9SWZqDq3fqk+PYR3wLvPwY73/JssyaUlZ3IhIJiE2jZSETNKSs8hMzyU3S4Mmz1QO0FBQ3MPMzJRDXm1pjpW1aWaNg7M1zu72uJdywNPXGXcvRxSKhx8SkiSJioEBVAwMoN9XXVm/YDd/z9/Dka3BvD38FXoMbV/k4+aZqaYqTDYO1kW2TVmWWT5hDcu+XUXdDjUZu/rz+15kNugNhOy/wLBfBuJb0QsrW0uO/3OavX8epsuH7QDo9mknrgWHs2LiWuq0q8H23/bi7ueGpY3pb0wE+KJRVFMo/wBaAC6SJN0AxsmyvLgotv2y8angxYQNozi1I4R5w3/juzdmUKt1NQb/0JeAao8/vhqy7wL7/zpC4651eWtEVwC+6jSZpBspuHrf+3T7v5TmStx8XW4r8FEUdFo9P4/+i23LD1G+ph8T//yEstUf7fpCZloO506GcykkimsXYgi7FEt25u2zY6xt1dg6WGFjqzZNk7S2QKEwQ5IkjEYjOq2BzPRc4qJTycrIIys997ZVxSoLJf7lPShb2YsK1b2pXq80Ht5OD9U+e2cb+o7uSpf3WjB/7Gp+/34TQXsvMvbXQTi4FF1iu8iLMQB4lyuaYSG9Ts+cjxfzz8JdtO3bnOELPnxgAFYoFQRU9WXfn4fp8+1bVG5UgbSEDEL2XSDuekLhRdr3Jr3NxaNX2PDzdgKq+vLepLeLpM3CLUUS5GVZFr+ZIlanXQ3mB89g07wdLBv3Fx/WGkHnD9vRd3yPR850mRqfxpH1J/Es41EY4A+sOYrRYMTF6+EC1L+Vr1OaC0dCiyx/t1ajY3y/eQTtvchbn7Tn3S+7oHzI4aDYyGQObDvH0d0XuXo+BlmWUZorKFPJk2Ydq+NXzh0vfxc8vZ1w8bB/rNk1qUlZJMSkEROZQuTVeK5fjmPfP8Fs+es4AJ6+TjRoWYlmHatTobrPAz8TJ3d7Ri8YQP121Zn1+XI+7zydyauH4e7z8Afb+zl74CLOpRxxLYLt5WTmMrHHD5zaHkLPUa/z3qS3H/p3Xr9TbfatOkLMtTi8ynriX9WHq0FhpCWk41nanbjrCbh4O9Gse0MadA4sMdWwnjci1fALIDMli6Xj/mLzvB1YO1jTb3xPOg1s89Cn+X9M+ZsrQWGMXjEMlYU5sWHxHFx7HHMLJa8PfQWg8B/3YQL3Pwt3MWvIIn459f0Tz96QZZmpgxZzYEMQw/7Xiw69H5xC2WAwcmTXBTYuP8r5U+EAVKjuTd1mFanVqCxlq3ihesDYvUajIydbg0ajQ683DdcozCQs1KahGisr1X0/B6PRSHRYEiEnwji5P5TgY2HodQa8/Jzp1LMB7brXwfohpk9eOnmdse/MwcnDnpn/jMTa7skuJuu0et70/ICm3eo/8Tz5xKgkxnSeQtSlGD6dN5CO77d+pNdnpmax9ofNmCnM6PtdDwBGth3P26Nfx2gwEn05lvbvtSz2tQUvEpFPvoQLPxfJz5/+SvDeC5Su7sdHP71H9WaV7/sag97A9P5zeX1YJyrUKUNCZBIntwVzNSiMToPa3jFrITcrD6sHzFpJS8zgHb8hdP2oPR/O6PNE72nLsoPMHrGS/mNe462h7e/7XFmWObjtHEt/3EFsVAoePk50fLMuLbvUwtXjzguAeblarl2JJ+xqAhHhScREp5IQl05aag4azf0zeCoUZjg5W+Pqbo+XtyO+fi6ULudOuQqeODrdOdadk6Xh8M7zbF9ziotnIrGxt+StAc3p+m6jB549nD18hdFvzqJlt7p8MafffZ/7IIfWn2D8mz8wceOX1Ov4+PPKQ0+F8c2rU8nP0/LNmi+o3frxFr1dCQpj4cjfafNucxp0DmRyr1n0/a4HleqXey5XSj/vRJB/CciyzMG1x5j/xTISo5Jp+XZjBk57977T2GYNXkBejobun3Vhy8JdKMwV1O1Q67YgcGTDSaJDYwnZd556r9R+4GrCSb1mcXJbMCvC5z70Mvb/ys7IpX/drylTzZfJq4fed457Ulw6M8es5czRa/iX96DXR61p2LrybRdBZVnmamg8Rw6GEnQinCuXYzEaTH/ftrZqvH2dcfe0x8XVDnt7S6xt1ajV5iiVCiTJdIaQr9GRk6MlMyOXlJRsEuMzuBGVSkpyVuF+vH2dCKxXmgaNy1Ez0B/lf86oQs9Fs2LObk4eCMWntCufTXyDSrXuf8azbOpG/pi5lRmbvqBKvcefLjiy3QRuXInj92uzH/uC7uH1J5jSaxYObvZM3Dwa/yo+j90egDN7zrFr+QEuHgmlbZ8WvPNVtyfa3stMBPmXiCY3nz+n/s2q6RtRmivoPbY7rw975Z4pVqe/N5fczDyqNa1EvVdq33ZRbuPP2zmy8SQNu9QhoJovswYv4O3R3WjTu9k9939zLvbAab3p/lnnx3oPa3/ZxaJv1zJ711eUrXbvQBJyPIzJn65Ep9Xz3ucd6Nij/m3BPTMzjy0bzrD9n2BuRKViZiZRsYoXNWv7UbmqN2XKe+DsYsptkpGZR2JSJimpOWRmacjNzUenNyDLpt67laUKWxs1jo5WuLna4eJkg5mZRHaWhuvXErh0IYazZ6IIOR1Bfr4eewcrWratQtfudfH2uf36RtDhK/z0zd8kJ2QyaHQnXu3V6J7vUZOTT/96Y6kYGMC4ZYMf6/O8cDiUz1qMe+w1DLIss2r6RhaPXkH5OqWZsHHUQ63XOP5PEDmZefetVqbX6ZFlWaQAfkIiyL+E4q4n8PNnv3JsUxC+lbz4ZM4Aara8+3o0g95AZkoWkplUOL/52OYgFn+1gs8WfEjZWgGoLMxZO3Mz7v6uNHm9/n33PaLtBKIu3eDXSz8+cIjnbj5pOxmFUsGPW7+853OCDl3h2yHLKOXrzNjZvfEOuDW3PCdbw8qlh9m47hSaPB3VavrS7pXqNGxSHnsHKzKz8gg+G03I+WhCr8YTHplMdk7+I7VRZa7Az9eZcmXcqV7Fm9o1/HB3syM/X0fQ8evs2XmBw/svYzAYadm2Cv0+aIGn161CGTnZGqaPXMXxvZfo/XFren3U5p77Wjx+Hevm7WbV5RmPNTb/ZYdJhJ+LYtnVnx45+Z02X8eswQvY8ds+WvRoxOeLhzxwGwmRSfz86a8c2XCSyo0q8OPBCWII5ikTi6FeQp6l3ZmwYRTHNgcxd9gSRrT+jta9mjJoRp87emEKpYLo0FgizkfTaVAb0uLTWT9nK32+7XFboeWti3fT59u3Cn++27J0gPcm9mRYk7H8MeVv3p/8ziO1Oyczj2tno+k98t5nAeGhcUz4ZDm+Zd2YumQAtg63hoWOHAjlx2lbSE/LoVXbqvR4txEBZdzIzc1n1/5L7Nl/mZDz0RiNMhYWSsqXcadNi8r4eDni5mqH2tIc2UzCKBsxyDISEhKm3rxklNFrDSSnZBMTl0ZYeBKHjl5ly45zAAT4udC6RSXat6pCo2YVSE3JZt1fJ1i/+gQH912m/8AWdCs427C2UTN2dm9+/Hoty+fsxt3LiTav3X0RWt3WVVkzdycXT4ZRt/WjJXM9tTOEM7vPMXDau48c4DOSM/nujRmcO3iJ3mO78+64N+87dKbX6Vnzw2aWj1+NJEm8P6UXb3zWSQT455ToyZcg+Xn5/DHlb1ZN24DKUsWAqb155YPWd/zD5mTkYG1vzbUz4Wz8eTv9J72No5upZz/u9Wk4ezoy9OcPCD8fxeZ5O8jP1fLaJx0pWyvgjn1O6/8z+1cfZdHZGY+UoOpKcCTD2k9l7K+DaPRKzTse1+sMfNxtNlkZucxe+wlOrqZpo0ajzOJf9rBqxVHKlvfgs1GdKF/Rk/SMXP5ce4L1/5whL0+Hr48TzRqVp25gAOZqJeeuxnLxejxh0cnEJKajyb//alCFwgxPFzv8vZyoXNqDauVKYae2IPhsNIeOXuXshRsozCRaNqtI7x4NCfBzITkpizn/28bhA6EE1ivN2ElvYG1tCrgGvYHR7y3myrkbzN/8Ge5ed5bFy0jJpmflEQwc353XBz38TBatRsug2iORZVgQPP2RpolGXY5hTKfJpMalMeLXj2jRo/F9n3/+0CVmDV5IxIVoGr9WlyE/9sfN9+ms3BXuJHryLzkLSwv6je9J615NmTV4IbMGL2Dn7/v5bP6g2y6eWdubZock3UghLTG9MMBPf28uSpWSfhN68uf360mMSsbazpJKDcoz7vVpzDww/o5/6Pcm9uTw+hPMHrqESZtGPXRvLivNtGLW/h6LgLatPkHktQTG/dynMMDLsszMqZvZtjmELt0CGTysHUqlGRu3BjNv8T5y87S0alaJN7oGolCZsXH/eb7+5R+SC1bnujvbUtbHlZqVvFFZFKQYlsAoy0iYppFKEsh6U8rhpJQsrt9I4fCZ68gyWFqY06hmAD3ers9I9w5s2hrCxq3B7Dlwmdc712ZAnyaMm9KdrRvP8NOMbYz4ZDnTf+qFtY0ahVLBiO/f4oNOP/DbzO18OaPnHe/Z1tF0ppKdkftQn+FNq/63iZir8UzZ+tUjBfignSFM7DETpUrJjL3fUal+uXs+Nystm4Ujl7N18W7c/VwZv+HLp5IXRyh6IsiXQD4VvJi+exw7lu5j/hfLGFx7BD1Hvc7bX3W7LQhUb16ZNT9s4qtOk3HzcSH6cgwz9n7HloW7SI1Lo0m3+oVT545tPkVmavYdQd7Fy4n3JvZk7qe/se3XvXR8r9VDtVEyuzUv/79kWWb970eoUN2H+i0qFt6/8rdDbNscQq9+Teg3sAU5ufl8PXETx05ep1Z1Xz4d0oaU7Dx+WnWAkNAYLFRKmtQqTaWyHuQa9VyNTeZiZCJ7r0c8VBsdbSyp4ONK0yblsbOwIOZGKvtOXWP38St4udnT59V6rFw8kKUrjrBuUxAngq4zeVw3XulaG2cXW8aNWs2kb/5m4oyemJlJuHo60KlnA9YvO8x7X3S867RP4JGGPSIv3uCPKetp/mYDAttUf+jXbZ6/kzmfLMa3khffrR+JZ8Ddz8JkWebA6qPMHbaEjOQsug/vQp/v3hJz218gIsiXUJIk0b5fS+p3qs284UtZPmENB9YcZfjCwVRpVAEAazsr/rf3O7b/thffSt74VChFZkoW0ZdjqNexFtWbmXLrHlx7jOjLsYWVhm6uYLypy+B2HFx3nIVfrqBu+5oPtYrWztF0NpGRkn3HYxFX4omJSGbod68XBrwrl+NYtvgArdpVoe8HzUnPyOXzMasIj0hi6Ietadm8EtN/3cXek1dxdbTh43eaIasktp0KZeumawCUcrLF38uZKuU8MEgyWoMBjV6PzmAAJJRmEhbmSlRmCszNFGg1OmITMliy/QSyDJ5Odrz6Sg287OzYsPssUxbtZK1fMF8P7EDLZhUYN3kjQ4avYOaUHtRvXI4hn7Zj9v+2sfnvIF59w9Tr7fhmXdb9epCjuy/cMdsms+CzsHV8uJwzep2eqX3mYGVnyZCHzEFk0BuYN3wp6+dspW7HWoz549N7ToFNupHC7I8XcXTjKcrXKcPkLWPuOmQnPN9EjdcSzsHVnlG/D2Xylq/Q5OTzWdOxzB26hLzsvMLntO/Xkkr1y2HjYE3wnvNkJGdSsX45lOZKrp0J5+jmU7w1oisqtYrE6GR+H7+aiT1/KHy9mZkZn/7yATqtnu/7zS1M+HU/Hn4uAMSEJdzx2PlTEQAENrl1MXjxL3uws7Nk6BcdMRiMfD3hb6KiU5j6XXeqVPOmz1fLOHTmOoPebEyvN+qx9MBpZq47iMFgpGvzqjRrUJZstYH9ERFsuXiFbReucuZGHNEZGaRocknR5BKblcX52AT2XA5jy/lQdl27Trgmgzq1/Xi1RVV83OxZtPUEU9bupUYtH8Z+2IHUjFwGfLuSuIxs5s3sjaWlOSO/WUNiUiZdugVSs7Yfvy8+QF6eFgDvAFdcPe25eDryjvcdG5Fk+mx8XR7iNwsrp/xNWEgEw34e8FBTHXMychjTeQrr52zljU87maot3SXAG41G/lmwkwFVP+P0zrMMmtGHn45OEgH+BSWC/EuibodaLDz3A68Oac/6OVsZWP1zTu8+d8fz4sMTqVCnLPYudty4EsuGOVvxq+RDlcam3r+Dmz3vT+7F6Z1n+WvahsLXeZcvxUc/9iNk3wXWzfrnge2xsbfCw8+Fy0HhdzwWE5mMpbUFrp6m4YyoiGROnwznjbcbYG2jZvlfxzh3MYZRw1/B2kHNR5NWoVQqmDHiNY5cj2baqn34ujnQs30t0hRaVp8+z7GwKPw8HKlTxYeA0s4oHZUkSXmEazK4kp3GlexUwvLSiZdz0dtKlPJ1oFZlb6qV8eRaYiprgs4TkpJIl5ZVaVmrLEt3BjF321HGDGlPtXKl+O6XrRy/GMX0CW+Rl6dl6sytAPQZ0Jz09Fz27rxQ+P58y7gRG5Vyx/u+GmwK/KWrej/w8wved4GVk9bR+p0mNHmt3gOfnxiVxPDm4wjec57PFw3mwx/63XXGVNz1BL5sN4EfP1xA+TplWHjuB7oP7/JEJSqF4iWC/EvEytaSj2e/z8wD41GYK/my7Xh+/HABuVm3evU1WlThr2nr+X38ar7uMhXPMh40faM+Hv5ugCn//Z6Vh/Cv5surH92ehqB9vxY0fq0uv379J5eOX31ge2o0Ls/Zw1cw6A233Z+dkYedg1XhUM2RA6EAtOlQjeSULFasPk6r5hWpUd2HEf9bj4OtJV992I5xK3ZyKSqB9zrXI1XSsuzIGeys1TSu4Y/eVuJYfAxBcbHoFEb8vR0p4++Mt489rqVscC1lg6e3HaX9nSnr64zKSsHF5EQOxUSRIOdQq4o3lf3cWBd0noNRkQx4rQEKSWL4/E10bV+DRjUDmP7rLlJzchnUvzlBwZEcPRFG1Ro+ePk4cWj/5cL3Z+tgRXZmHv915sBlPHydcS1158ybf8tKy2Za/7mUKuvB0LkDHvg5Xz19nU8afEV8RCKT/hlNh7tcNzEajWz6ZTsDa3zOlZNhfDpvINN2fvNIM6aE55MI8i+hqk0qMT94Ot2Hm9IcDKz+OWf2mHr11ZpW4vud3+BV1oOB09/lna+6UaqMR+Frt/26l71/HuKLxUPuuPgmSRLDFwzC2cuJCT1nkpV253j7vwW2rExOZh6XTl2/47F/X5C9cP4GPn7OuLjasnFrCDqdngF9mjLz973k5GoZ9UE7xi7djtFopE+nuiw+FERmnoZWgWW5nJ3CocgofN0d8PNzJMdST2huCpeykkgx5GKwMKKyMcPCxgzJUiYTDVdzUriQlUSaeT4unjZU8nfjclISh2KjCazqg5ezHT/vPkadmr6U83Ll69+20a19TUq5OTBx/nbata6Ch7s9azYEIUkSNQP9uXjuxq33JFMwK/+WnKw8zhy4TL22988RYzAYmdpnDukJGXz528dYPiAJ2qG/j/NZ07EozBX8dGTSXWsJJ0YnM7rDRH76aBGVG1Vg4bn/0WlgWzHvvYQQQf4lZWFpwaAZfZh5cAJKlZKRbcYz++NFaHLzKVPDn1bvNKXRq3Vve83xf4JYO3MTg2f2uy3wn9wezIa529i1/AC2jjZ8vfJT0uIzmP7eLxiN9x6fD2xVBZXanAMbTt92v429JVkZeYVBMfZGGn4FK133HQylVnVfcvJ17D5+hd6d67Bk50myNVrebFuLn3YeobKPG3aulmy7cpVKvm5Yuag4kxZHDlpKe9lj4SihtdSSQAaR2mSitMlE6ZKJ0CYTY0gj10KDZG/Ep5QtltYKTiXHkq3WEVjei6AbMURps2hTsyyrT5yndGkXvFzsGb9yFx+93ZTYpAy2HLpI+9ZVOB0SSWZWHr5+LuRk55OZYeq9Z2XkYm13e3A+ujUErUZH89fvPy1x1YyNnNwWzOAf+lKh7r1z3MiyzJ9T/2Z89/8RUN2PuSem4lfZ547n7Px9Px9UG87Fo1cY9stApm77Wsx7L2FEkH/JVWlUgXlnpvP60FfY+PN2hgSOJPRUWOHjJ7cHM7XPT1w+cZWVk9fx5uevFtbkBJjY8wd2LN1HVmo2637czIpJa6lQtwyDpr/Lsc1BrJzy9z33bWWjpmGHGuxddwJNrrbwfvdSjuTl5JORaprfnp2twdZWTWZWHpHRKdSp5c/6vWexUCkp5e3I8ctRvNOmFr/sOUbNAE9i9dlEpqdTrawHJ1NisLOxwNVdTYwhjRRjJt7OapwcJVQ2WsyttSgs8zGzzMfcSovKRoutvRF/V0u0ZnmEaZKxcVbi527P0fhoXD1tsbZUsTcygva1y7Pm5HnaNihPamYupyNiqVLGk/W7zxJYww9ZhouXY7EtCOjZ2aYCJvHRqbj9Z0hm96pjuPs4U6lO6Xt+XkG7zrJ03CqadW9A50Ft7/k8g8HAT0MWsvirlbTo2YgZe8bdcWE2MyWL77rPYFrfOZSu7seCs/+j8yDRey+JRJAXUFtZMOTH/kzb9Q2anHyGNRrDH1P+xmAwULd9TaztrBjacAwV65WjXd8Wha+bOXAeiVHJjPr9E3qP7c4H094lKdp0QbHrR+1p/U4Tlo9fw6mdIffcd+f+zchOz2XXX0cL7wuoaJqeea2gwpGEhCzLxMSmA+Dv68zhM9dpWCOAtYfO4eVix7GoG1hbmJNrbiAtL49SXvaEJMVT2deFsPwklCojbs5KcsyyyJLT8bQzw9dBwsFGg72NBntrDQ42GrzsZXwdlBgVOaSSho29ATd7FZezE/DzsiNZk0OGUouLnTUHoiOp5ufB8qPBtKpdljUHz9K8blnCbiSjtjYVwLgRm47RaDojMZMkcnPyiY1Kxb/crbHuqCtxBB8Mpf07je4ZZBOjkpn8ziz8KnkzfMGgez5Pk5vPuNemsXn+Tnp++Rqjlw/DwvL2NAcntwfzQbXhHN8cxAff92bG3m/vOU9eePGJIC8UqtWqGvNDZtCkWz2WjFnJl20nkByTwidzBvD+lF6E7L9AcmwqABvmbiN43wWm7x5XOPMiNS4dtZWq8ELq0LkD8Kvsw6S3ZxFz9e7lAqvUL0uFWv6sX7incGinYnUflOYKQo6Zxupt7dRkpOeRVdATlhRmJKRkUcbPheCwWOpU9uVMRCy1K/hwIS6Rcn6uXExOxLeUHRcz4injZkOyMQ1HS3C3MyCZZ2KmSMPTRkd5Bwiw0xJgp6WsvREfWwNq80y0Zqk42ejwd7DghjYRD2dzkrVZmNtJGDBgsJLRGvRY2JmTmZePtb0FGq0evZkpoMckZQCm3ntWwUVWG1s1l85EIssylWrdKm24es4OVGpzOvZpetfPKC9Hw/i3fsCgN/LN6uH3TAaXGp/GiFbfcmLLGYbONf3O/n0w0Gq0zB22hK86TsLO2ZbZx6fw1oiuYuZMCSeCvHAbW0cbxvzxGV8sGULoyWsMrPEFRzaepMfIroz49SMMOgPZ6TkE7z3PiCVDCnuJKXFp/Db2D6o0qYRCqUCSJCxt1Hz39xd4lnYn5y6zScB0sfa1Qa2ICUvk6FZTj19tpaJKbT9OHboCgGcpR2JupBa+JrVg2X+e3pR/Jk2rQa1ScjI2hgpeLgQlxFLJ14XwrBR8XNQkaf6b/ygAACAASURBVNMJcFSSYUzCzxY8rbOxME/BWplCKctcStvkU9omHx+rXBxUaZgrEnG1yqacvUS2MQV3OxkzhRalpQ4tWiwdFERmplPRz42jEdFU8nHjdFQcapWSmLRMABILimkjQWxMmqkwuK2a4/suo7JQUqW2PwAJUSnsWfP/9u47vqmqj+P452Sn6V4UWvYUmYIoWxkiCog4ABUH+uiDgLKplA1lK0PcC9yC4ECRIU5AkCEqyIYCLd1turJzz/NHahUZ8ghaqOf9evXVkeTml9vk29uTc3/ne26+r8MZ13mVUjL/kRc5uPMoY5cOJr5O3GnXATj2ywkebzuelD0nmLRiFD0HnTrz6eju4wy9dhwfPv0ZvYd0Z/H3M6nTTM17/zdQIa+c5tezZZ/dPpu4GjFM6j2HZx5/laoN4qlUPQaT1YTJYqRaw8B8bp/XR+IN07iub1va9zm1PXFcjVie2TqDei3OPtbcvudVJNSpxOuzV5X9F3B1xwakHMgg/UQedepV4sSxHIylPeTzCwNj9fklgT8cBzNzqF0lmtwSB069n+hQK3vsGdSMtZHnKSQy2IdOV0JCsBtEDg1CJfFWO6HGDCIM2VQy2alkKiDKmIdNn0ac1U6DEA9mfR5hlgLirODXFxBqFlhtGmklBdStHMnuvCx0AoKCTRzOzCU+OoyM/EC4F5ROSw22WTh8IJNadWLR/Bob1/7M1R3qY7EGhnOWPb0WnU7QZ9CZWxCvXLiar97bzP1T7zxrr5jdG/fyeNvxuB1u5m6YRNvfzZuXUrLq+XUMveYJ8jLsTF+VyOBFA08bwlEqLhXyyllVrR/Pgk3J3PrYTXz49GcMazee1IPpaH4Ne1YBrz7xNuuWfsWU2+dR7+raPDjzbuD0fjR/9mae3qDn/qTeHN+fzvp3A2Pz7W4ItNr9ctUPNLgyHinBXtrUzG4PHMkXOd2YjXpS8woQRtAbdBy252G16YmymUl15lElTE+Q0Y/ZWEB1myTKnEuooYB6QU7izflEGFIJ06USrjtBuP4Ylc12alsLqWJ2E2zIoUEIeMmiRogRzVBMsc9BregQ0twFOHxe4mPCyHYE6jGa9RQ5Av3qnc7A8oIRYVYO7DtJoyZV+XHrEfJzirmuR6Dr5vED6ax5axPdB7QjuvLpZ6xuXb2TF8e8SbtbW9F3zC1n3HffrdrOEzcmE1EpjKe3zKRBq9+ajBXmFTH1jidZ9OhLNGofmDZ7zc0tzvm7UCoeFfLKOZnMRh5d8ACTV44m40gmg1uOZesnO5i1dgJSSvIzC+hwW2tGvzoYOL+FwM+kTfem1G9eg7efWo3b6aFSfATNW9fhs+XbaNCwCnq9joN7ThIaYiE7O3C07PL4kEj8mqTY6yUyNAgpJMdL7FitUDnYRJGvgCBjITVtAk2mU9+mEaw7SVWzl/qWIqqZcokzZlHJmEWCKZ+65jxqW8FIKnVtLgS51A8Ft8zFZvRTLdxMltdOic+NzWpEZ9KRWRw4H8Dj18paOvw6W6gwrwSfT6PlNbVYt3I7waEWru4QOHv4pckrCAo2c9eIm07bHyf2n2T2fc9Qq0l1xi4dcsb+7muXfMnkPnOp1jCBp76eSqXqv019/OW7/fy3+Wi++3g7D88ZwIzV44iMO/dJVkrFpEJeOS9te7fi+V3zqNYwgen95rPo0ZcYvGggfcfcQtd7OwJ/PeABvl2xlWu7NCTrRC7vLVwDwC0D2pCTUcCObw/Q9KrqbP72AA0bVOHQgUxMRj1ulxe3LxCqTp8XnUGgNwikkGS5CxF6NzVCTATpPUiyqWeTGMmlnsWFRZ6gqslGPbOOOsZ86hjzqWtyUcsSjk2mU9OcT4TeRVVLMT4tm+o2HZEWH3leOz7pxWSCkCATDr+XYo8HCTjcHnSljz8vp5jQEAv7fkolJNRKpdhQNq7bTZfeLTBbjOz6dh/bN+yh7+M3Eh4Tesq+KMwrJqnnLAwmAxOXj8BcOrTzKyklb0xdzryBz9L0uiuZu+G3KZJSSlbM/4QRHSehN+hZuDmZO0b1OuciIErFpn7zynmLrRrN/G+mcsfInnzywnqGth5H2qHfZs38PuB9Xh8lhefui65pGns27WfRkFd4d/aHFGTaMfncLFu0lpMp2VzdsT5VqkXx0Rubua5LQ9JO5FG1cgRp6XaqVoog3+6A0pEhv1/DLzXMJj1Cr6HXSZx+Fy7NTpzVT4IFpMymqslOmF5S01iMRTtMlLE6CearSTC3opKpIWYtnSqGLOJNYYTp0gjRe6kVpEeQi1srxKT3YzVCqMWIJjRc/t8WH8krKAFNotcJDu3PoMmVCXy38SBtO9Rj5ZKNAPS+ty0el5dnEt8lrloUvR68/pR94vP6SL5rIblpeUxZOYrKNWNP22fPPv4ar09eRtd7OzJj9biy2TbOYiezBizi+ZFLubZnC57bMYf6Lf/6ouBKxaBCXvm/GIwGHp57L9NXJZJ9PIdHW45l4wdbT7ve04NfZkTHieSm559xO1JKdqz/iQ1vf0tUlQgWbZ7OI3MH0LpHC4Tfx4sTAkvL3XpfW/b/dIKIEAtBNjPZpY29wixm0k7aEUCIxYyU4PH58aMhhESn0zDpfRh0XjyaHbOwU9NiJEhIokQaYcaa1DDGEebfj1k7iVlmEuzfRzW9iThTc6zyAJWNwcQaCvBrWZh1boINPsJMgmCTAYTEJzWkJrEYDAjA7fZTVOCkaqUIcvNKCA+y4HR4aNO+HutWbKdTr+ZUio9gxbPrST2UyZA5d2GyGE/ZJ08PfZUfNvzMY888dNoiHn6fn6ceej7QRXJ4D0a9+igGY6BbePqRTB5vO56v3t3E/dP6Men9UQSHn1/LYqViUyGv/CXX3NyC53bOpVqDeKbcNo+3pq84pYVBu9uu5eShDEZdP6lsbv3v5aXns27pV1RvmMDd4/pgMBo4uvsEqQfS6X5fR7au+5kvV3zPDbe1JKZyGG8v3kDX7o3Ztukw9evGkX2yAM2nEWwxYdUbcLu9lLg9uL2B2Tk6ITHoNMw6H1adF71wIbU0Khs8hBjiiJInEfrKiMhliNit6GK/Q0R/DpZO2LSDxJiuIph0hCwkVO8l0qgj3KjHrAcI/Nfg07TAilEGAxaDASSkZxRgFXosZiO/7DxGzdqx7PnuMB63j9seaE96SjbvLVpD25ub0eL6hqfskw8Xr+GzV76gf2Jvut1/3SmXeVweJveZy9olX3LvpDv575P3lQ3BbFu7i0EtxpB9Ipfpn47j7qTb1JmrShkV8spfVql6DE9+NYXO97RnycR3md73KVyls0uu7taMmZ8lkZOWx2Otx5H6u5OhpJS8Mu4dIuMiuGXwjUDgSPTnb36hVfdmPDj1dq5oWYsXJ76P2+nh3sdu4OCeNBLiwvF6/UTZLGSmFRBiNRFmNlNS5Kaw2I3X4y8dvgkEnAD0QsOs8xOkc2MRPpB5ROlsCH01ROQ7CFMzhAi8DIShGiJsPpivJ0Q7hg4vEXoz4QYDFp2GTmhINHxSw+3XcHp8uFxe/F4Nm9FIqMUMfkg9mkPjBlVITcnlhu6N+fit7+jUqxlVa8Xw1OOvo9PreWT6nb/flWz5dAfPj3ydNrdczX1TTr2sKL+YxG7T2fLJDoYufogBk+4AAkM3b0xdTtJNM6hUPYbnds7h6m6nr5er/LtdlJAXQtwohNgvhDgkhEi8GNtULg8mi4mxS4fyyLx72bjyex5rPY7MY4HFLxq1u4L530zD4/QwvP2Esp44ml9Dp9dxfb/AotH7vj/Epo+2c/inYzTpcAVmi4nBs/tRZHfwwvjlXHdzExJqxvDpm5tp3a4uv2xPwWY1EWm1kptVjNvpQ/hAaIBfIDWBJkXpcL1EJyRmIbAKHQYEei0NEXQ3Qnf6ghlCCETwcARuQvWxBOl16PEBfjTpx+PXcPr8uL1+vF4Nr8dPSbGHwgIXJimoHBmCx+GjMLOI6JgQDu5IQacTPDDiRta8uYndWw7x3+l3nNJO+PCPx5g1YDG1m1bniTeGnvImaUlBCWNvmMbeLQcY9/Ywej0aOMmpKL+YSb3n8PrkZXS+pz0LNk0vawetKL93wSEvhNADzwDdgYZAfyFEw3PfSqlIhBDcPqInyavHkXU8h8FXjy1bkKRO85o89c00LEEmxnSewv5th9Ab9FzZtj5zBj7LnPufYcXCT8nLyKfbfdfRvFOg1W7tRlW5fVBnPl+2hdVLv+Gh0d1JPZpNlehgnCUeasVHkpGSh9/tx6zTYTMYsAgDJvToNB2apsOv6dCkDk2ClH50SIylR+0Yz/EUNdQHEYRJGNFJHxI/Pqnh8mu4/RK/JtA0AX4R+OPiBzySglwHnkIPVatEcGRfBtd1uoJvP/uZWwa0we1w88rUlTRtW4+u/VqX3VVehp2knrMICrMy5YPRp8ykyc8qYMR1kzj60zEmrxxd9kfx1/H37Wt3MXjhQMYsGaLWXFXO6mIcybcCDkkpj0gpPcC7wJnP3FAqtKu7NePpLTOIqBTOEzdOZ+XCT5FSUq1BPPO/nUZodAhjb5jGz9/upfvATgxZ+ADX9mzBwOn9uWf8bTRs/dtyf1tX7yQnJZNgEzw/fAkFqdm07tyQde99T+t29UjZk45Jp6dSqA2jT+Bx+PA6/PhdGnh14Dfg0wx4ND1eqceLDh8+yk7Tkp4zPoYAP0gvmnTilT7cmobDr+HSBF5Nh8+vQ+czYPLrsUgDwTojoSYTQQYDxdkO/IVuEhIi2PXFL0TGhHDHfzoyb/ASdDrB8IX3lo2XO4tdTL5tHiX2EpI/TiQmIaqsgsxj2Yy6fhInD2Yw9ePEspOYdm/ax+BWidgz7cxaN4HeQ7ur8XflnC5GyMcDJ373fWrpz04hhHhYCLFdCLE9Ozv7ItytcimqWj+ehZuTad2zBc8NX8KCR17A7/cTHR/Fk19OJrJyBEk3z2DH+h+5qnNj2ve5hoM7jpzSdOuTF9bz2oT3aHVjM0a++AgiKIinH3uVPgNa4/X60QqdaF4/VaJCKDhZhLfIi84Dei8IL2ge8HgETq8Bp9+CSzPi0oy4pR6PLF2Fynv2zph4fwC8OLVCSjQXDo3ANvxGPD4Tfq8evwf8Lg2f04+72IvD7sLsE8RXCiM71U6jhlVIOZDJoPG9WPfWJvZuP8J/k/tSqWogyDVNY8bdCzmw4wiJbwylZuPfGpZlpGQxrN148tLtJK8eVzbOvuGtbxnbdSqhUSEs2jLzlJbPinI2FyPkz3QYIU/7gZQvSilbSilbxsSoRQkqsqAQKxPfH0W/xFtZ/fIGJveZi7PYSWy1GOZ8PpG4mrEk3TyT71ZtRwiB1+Pj6O7AccJnr3zB5299y2OLH6TjHa1pfVNzBk7qg8ercXRPKnc92pntX+3jmmtrc3J/Fha9njCzGaumx+gTGDwCg1eH9OrxePUUew0UeK0U+YNwyCB8mHATjHSuQErvGeuXJa8hMeGQPko0E8WalWKflUKPhRKPAZ9Hj86jw+TTYdIEQUKPFR3OXBeOzGLq1q3E5k9/5Kq2dYmKDOLVaR/QuntTOt3+W0+ZJRPfY+vqHxj05L20veW3xVnSDqUzpstU3A43T341hSYdGqJpGq9PXsasAYu44tp6LNg4jYS6lf/eX6JSYVyMkE8Ffr/kTAJw8iJsV7mM6XQ6HpxxF4MXDeT71T8wpus0ck7mEV0lkvnfTKVO8xpM7/sUW1fvpPNd7ajeMB6vx8fhH1O4O6nPKUM3+zbtJSwymJenrOSa9nWp2SCOg98fISoymDCDEWeWA1+xF5wSnRtwgubS4XYbKPGYKfJasPss5PnM5Pn15PmzwX8CWTgZKU9duUo6Pwb359g1LyXSgl2zkecNIt9rocRrwusxIF16dC6BdEpkicRb6EWUSCpHBuMqcmH0+EDCf0Z356lhrxMRG8bwBQPKhlU+eWE9787+iJse6lz2RipA6sFAwJcUOJjxWRK1mlTH7XQzvd983pi6nK73dmTmmiTCY8L+kd+hUjFcjJDfBtQVQtQUQpiAfsDHF2G7SgXQe0h3Ji4fScru4wy95gmO/HQMW5iN5E/HUb1hApNvncPWT3eg0+nIOJpF2qEMGrVrUHb75LsWUJRXzLNbZ2ALsTLlnsVESA/2nGISYoLJSy0g3GohRBgwecHgAr0L9C49mttIictErstKtieYXF8ouf5gCjQjdqkD53Jk7m1I1wakeyNa4WRkwSjcmMnVnGT6LGT7gsnzBpPnslHkMuN2GRBuHUa3wOwRWDRBsDCgd2vkpeTRvFk1Du46zoOjbmT5ojWkHsxkxIIBhJSemPTDFz/zzLAltOrenKGLHywL/qM/H2NY2yTcDjez102gQau6FOYWMabLVDau2MrDc+9l9GuDMZqMZ9zPinI2FxzyUkofMARYC+wFlkkp91zodpWKo23vVizYOB0pJY+3TWLrpzsIiw5lzueTqNm4GhN7z+GLt78lsnI4hbnFrH5pAz9v3Muglono9HomrxhFdOVwRj59H5kn8ti39QDVQnX88t1hrmxYBU92CZ5cNwY3mLwCo1ugd4B06HE6jRS4gshxh5DhCSXNG0aaz0y230WWpqH5DiPtg5D5A5GOdyiUZk748jnps5DhCyXDE0amO4R8dxBOlwnNYUDn0CEcgEMii/x4c93Y0BEdFcLeTQdp0qoWZgFfvP89d4++meYdrwAI9HrvM4+q9auQ+PoQ9KWtk9MOpTPuphkYTAYWbEqmTvOaZKfmMrrzFA79cJQJy0Zwx8ie6g1W5S8Rf2wL+09o2bKl3L59+z9+v0r5yk3PZ+ItsznyYwrj3hlO+z7X4ChyMvGW2ezeuI8pH44htloMb0xbTlz1GEKiQug/tnfZ7X1eH7P/8xLfrNhCi+sbUqi3kJFegD/UigwxkSs9uEJ1eEJ1eELBHa6hRfiwhruJCS2iWkg+tWzZVDfnEG+0E2fwYBMlWIUeKSUeNDwEke7Tc9IbznFPJIcdMaQURZFZGIqjwIrIN2K0C4wFAksR2JwQ6tXjTCukRkwIhZmFJM65kykDnqVu02rMWjkcvV5H7sk8Hm8/EZ/Hx+ItM4iOjwQCQzQjOkzA79OYu2EStZpU59gvJxh7wzScRS4mrRjFVV2alNevTLnECCF2SCnPvdr7H6gzXpV/TFTlCGavm0DdFrWYevs8Vj2/jqCQwPzwmo2rMan3HDKOZpL09jAenHHXKQFfUuhg88fbCbYZqVwjlj270hgwpDN+j49Ii5Hi9CJsGLD6dJhcYHCCySHQlRhwlxjJcwRx0hHCcUckx1zRpHiiOeoNI8UbxjGfmeN+C8d9IRz2hpDiiSbFHc0xZxQnS8LJLbHhLDEhSwzoHQKjU2Bxg8ULotBHSXoRdWvGkLo/g4fG3MTTo94iOMzKEy8+hF6vw+VwM77XHIryipn64ZiygD/687GygH/q6ynUalKdbWt+4LHWSWh+jQWbpquAVy6YCnnlHxUcbmPeF5NpdVNzFj36Ei+NeYOg0CDmfTGJ2s1qMPX2eWz+aBt6Q2Dd0V87WW764Ht+2LCbGg0TeHJ9ErZQK8+OfYeHE2/m5MFM6taMwZvtQNh9mF1gdoHRITAUC2SREUeRmZziYE6URHDEEc0RZwwHnLEc9lTioDuGg+4YDnliOeCqxGFXLEcc0RwvjiSzOITiYgv+YiP6Ih2GYoGhBESxBvlezB5J1bhwjuxIocdd1/Ldh9vIPJFL4gsPElkpDK/Hx9Q7n+LIT8cY99ZjZStkZR7LJrHbdHR6HU99M5XqDauy8YOtTLxlNpVrV+KZbbOp2ajaWfejopwvNVyjlAuf18ezj7/GqufXcePATgx7/mGcxS6Sesxk39aDzFg9joKcIp4ZtoR7kvpw8kgWtZtW58YHAq159+9MYWSPOUSEmQlPiObA/iyir6pNjtNFiRk84Xq8pcM2nlCJL0xDhLqxBbuJtDmItpQQbi4hzODGogtMpfRJHQU+KwVeK1nOYPIcVoqLrfgKTegK9JgKBOYigalQw1wkCffrocCDsdhFQtVIrmpRjWUL1/Df5Du55aFAnQsGvcTqlzcw/PmH6f5gJyAwBj+60xQcRU4WbJxOjSursv6Nr5l7/zPUb1WHmZ8lqQ6Syhn9leEaw99VjKKci8FoYOgzDxEaFcJbySsoKXSQ9M4wZqwex/D2E5jUew4T3x/FTQ92Yunk5Yx9fQitewTO+nQWuyjIzCcq1ETGiVzi61bBYLdTeCANc3wUBmGkoNiH1IGm0yGFQAodfmmmpLTdgcNjIs8chMXgJcjgQ6Dh04wUeY04vSYKnWacDhP+IiO6Qj3GIoGxWGAokVjcAptf4M4qIdpiAJOBrj2a8Fziu3Tt17os4N+Y9j6rX95A3zG3lAV8+tHMwDx4p4cnv5pCjSur8t6cj3jlibdoev2VTPs4EUuQWn9VuXjUcI1SboQQ3D+tHw/PGcC3729hUu856HSC2esnUrVBPFNum8tVXRrTpGNDdn2xu+x229bu4qNn13L3uFu5fWQvftx6hNrNa2DS/FhcPnx5ToKlHosTjMUSYxEYCwWGAh0UGHHareQXBJFZGMLJgjBS7BEctUeRYg8nozCEbLsNR6EVrcCMvkCPsVBgKhKYHRDkElhKJL5sBwnRIRRlFTHg0et5aeL7NLq2DkNm9wdgzWtf8sbU97nh3o48MK0vADlpuYzsOAlnkZPZ6yZQq0l1Xhr7Ji8nvknHO1szfZUKeOXiU8M1yiVh1XNrWTz0FRq2qc/Uj8bi9/kZ3WkK6UcySf50HNEJUcTXiWPb2l08N2Ip/5l9D617tMDv8zO291P8/OXP3JnYm08//JGg6BAyil0YYoMoMYPLKvGE6PDZwBcMXquGtEqERQODH6Ev7VepAX49mkuH3qVH7wB96Ri82QFWJ5gK/ZDnJiEqmMz9GTw4ohsfLPoMs9XEorWJhETY+G7Vdqb1nU/jDg2Z8UkieoOerOPZJN6YTE5qLk99M5XaTWuwcNBLfPrieno80pUhix9Er9eX7y9BueSp2TXKZavnoG6Me3sY+78/xIgOE9HpdMz9YhKVasQwuc9c7FkFAOz6cg/9E2+ldY8WaJpGbno+jVpUJyI2jFWvfM0dD7TFfjKfuFArviwHNjdYHGAq1DAVgbkQjHaBwa5D5OrR2c1oeSa0fBPSbkbkGjDk6zHmg6lQYCoEc5HE6gBjacBXjQ0lc38G/R++jk9eWI/m15jy5qOERNjYs3k/M+5eRO2mNZj43nD0Bj15GfmM7jyFvPR8pn/yBDWurFoW8P3G9uaxZ/+jAl7526iQVy4ZHe9sQ/LqcaQdyuCJ7sn4fRozVicRFhNK4g3T2L15H1nHs3EWuwDY+91BPly8hsyULOZ/NRGrSfBG0tsk2CB7617iI4PxZTkI8eowF2kYCzQsRQJzgcBSIDDZdZjzBBa7niC7gaB8PaZ8HeaC0usUBf44WB1gLtIQeW7io0PI+OUkPfpdw9aPt1GQU8zUt4dQrV5lDu48yvhes4lOiGTqh6OxhQWRl5HP2K7TyMuwM+OzJBpcU5dZAxaVBfzAGXepk5yUv5UarlEuOVs+2UFy//nEJEQxa+14jGYjw9qNpyivmKHPPMR7c1cRFhNKTHwksdWjGTDhdl6b8C5fLfsOr86IX5OEhRg5kVZITPM6ZBS5MEUHUah5EeFGnHoNgvR4dBKpB6EHg16Hz6eBDwwa6JwSs1dAoY8QnQFvVgkJMaFk7EunZ/9r2PftLxzbn86kpYNocX1D0g6mM7LzVAxGPU99OZnYatHYswsY1z2ZE/tPMu3jRK64ti7je85i1xe7eXjOAO4Y1au8d7VymfkrwzUq5JVL0u6Ne0nqMROLzcKczydishgZ02UqjkInY18fSlTlCOJqxmILC2LRkFfYu+Ug87+ZQm5GAWP7zKcoI4+ruzVl+88ZmEOsOI163ELiDzLg0mkYw02U+H1IPRiMehDg9Wro/BKT1KF3S0SJH5vQ4bO7qRJpI/1ABnc+2J6tH20j/VgO4195mFZdG5N2MJ0R109G0yRz1k2gZqOq5GfaGdFxIlnHc5i0YjRNr2tIcv8FbFm1g9GvDabrvR3LexcrlyE1Jq9UGI3aXRHod6NpJN00A4BZa8djshhJ7j8ft9ONLSyITR9tIzc9n8VbZ2AJMhNfK5bp7wxB83rZsuZH7h10PWa9QMsuJCEyBHKcREkjlkINc74fS4GGzSEIdghsRWDO92PI8xEljejtHsINRsIF5Kbk8NDIbmxcvoXstHymvjWYVl0bk340i8TuyWh+jSc3TKJmo6oU5hUxvucssk/kMnvdBBq1a8C4m2awZdUOBi8aqAJe+UepkFcuWTUbVWPaqicoKXDwWOsk/D6NhZuTCY8NY3zPWRzdfZyc1Fwat7uirNlXflYByf3m065nC2o3q8nLE5Zzw81NuKp1HdJ+OkFCeBBhBgPe9BKsxX4qG61ECxPRmIgzmgl2gSHfjb7ES72qkeTuSyc02MIjo7qxbN4nlBQ6SF72GM3aNyD9SCajOk+hpMBJ8idPUO2KeIryixnTZSpHfzrG+PdGEF+3MiM6TmT3xn0kvjG0bOFyRfmnqOEa5ZJ3bG8qo66bhM/rZ8qHY4hJiGJY+wm4HW66DezEkZ+O8+j8+8k4mslzI1+nw23X8NDMu3E5PMwbuoRNn/xAq66NuKZHC95+7ktysgpJqFuJsLgwSlxeiovdAFitRkKDLTjyS0jZnYbZbOCWe9vgLyxh5XMbqFYvjomvD6JKjRgO7DjC+F6z8Xt9zFk3gdrNapB1IocnbpxO+uFMJq0cTaN2DRjTZSrHf0llwvKRtOrevJz3pHK5U2PySoWVfjSTpJtnls0zDw63Mbz9BBxFTpp1aow12IrFZqZR2/p0uacDAL8+tz9++UtenvoBliAT/YZ3xxRqY8372ziyj1oDxAAADOxJREFUN/2M9xVbOZzOvZtTOS6M5U+v4cSBDG7o34ZByXdisZn54cvdTO4zj9Co4MARfIN4CnIKGd5hIrkn85j60ViqN0xgdOcpHN+bxuSVo2nd8/96XSrKGamQVyq0nLRcHmudREmhg+mrnqBS9WiSeszk5KEMZq4ZT4NWdTBZTEAg4H8/NfHYvpM8P34Zu77dT1RcGN0HtKPhtXVxezUyTuQhBISEBRFTKZSju4/z+bKtHNx1jMo1Ynh0Zl9adgqsp7pu6VfM/+9LJNSrzMzV44iOjyQnLZdRnaaQdTyHWWvHEx4bxoSeM8lOzWP6qkTVSVK5aFTIKxVe1vFA98bsE7mMWTqExh2uYOR1k8g4mkXSu8Np0+vq027j9/nRGwI943d8+QsfvvgFO778BQCDyUC1unEIncCeXUhuRuCkqxpXxNNzYEdu6N8Gg1GPx+XhhdFvsur5dTTv3JgJ7w4jONzG4R9TmNBrFiV2B8mlof94myQ0TTLlg9E0bF3/H90/SsWmQl75V8jLyGfKbfPYu+UgiW8+xlVdGjO+x0wO/ZDC4EUD6fnfG8qum7LnBBNumc3A6f25rm+bsqP73Aw7O7/aS8reNFIPZwEQEhFE3abVadq2HjWuiC/bxsGdR5n34HMc3X2c24ffzMDk/hiMhrJpnrbQIKZ+NBa9UU/STTNwFrtYsHEa1RtWRVEupr8S8kgp//GPFi1aSEW5EC6HSw7vOEF2M/aVn764Xhbbi+W4m5NlF3G7fHbYa9Ln80kppTyxP00ObZMkuxr6ysTuyXL/tkPnfR85J/Pk4sdeld1M/WTfhEfk1tU7yy77atlmeZO1v3ygwWMy83i2/PHrPfL22IGyb/x/5MEfjlz0x6soUkoJbJf/Z96qI3nlslVS6CC533y2rdnF/dP60Xdsb14YuZQPn/6Mq29sRtI7w7CF2fD7NT5+Zg1vJq+kKK+Yq7o0pscjXWl+fSNsYUGnbNPn9bFn836+eGcTn7/xDX6fn57/vYH7ptxJcLgNTdN4Z8YHLJ30Hg3b1GPyytH89PUvzLx7IZVrVWLqR2NJqFelnPaIUtGp4RrlX8fj9jJv4DN8+c4mOt/TnhEvDWL90q94esgrVKldifHvjaBWk+pA4I/Cx8+tY9Vz68hJy0MIQXy9ysRVjwEB9qxCUnYfx+f1Yw4y06lfW+4c3Yv4OnEA2LMLmPvAM3y/+gc6392e4S8+wrolX7F46Ctc0boe01c9oRb7UP5WKuSVfyUpJW/PWMmSCe9yZdv6TPlgDCl7TjCj/wKK8ksYmNyfWx+/qazTo9/n56dvfmHP5gMc+ekYWcdzAAiJDKZ20+rUv7oOLbs1xWqzlN3Hpg+/Z+GgFym2O3hk3r3c/EgX3pi8nLdnrOTaHi0Y986wU66vKH8HFfLKv9rXyzYz5/7FVKoew6SVowmLDuGp/zzPdx9vp3azGgxeOJDG7a/4v7Z58nAGL4x6nc0fbaNW0+qMXTqU6IRIZt/7NN+v/qFs6cJf16RVlL+TCnnlX2/3xr1MuW0ebqeHkS8PosMdrfnm/S08P2IJOWl5NGrXgNuG9+DaHi0wGM+8+qWUkkM/HGXFgk/4+r3NGM1G7hrXh9tH9gysP3v3QvIz7Aya/wC9Hu32Dz9C5d9MhbyiEDhpauodT7J3y0G6DOjAowsewGAysOaVL3j/qVVkHc8hPCaUhm3qU7NxNeJqxIIQ2LMKSNl9nD2b95NxNAtrsIVuD1xP37G9CYmw8caU5Syb+zFxNWNJemcY9a+uU94PVfmXUSGvKKV8Xh9vTV/B2zNWEh4bxkOz7qbLPR3Q/Brb1uziq/c2cWD7YdIOpqNpv70GYqpGUbtpDVr3bEm7264hNDKEnZ//xNNDXib1QDo3DuzEoPn3ExRiLcdHp/xbqZBXlD84sOMwix59if3bDlOvZW1uG96Djne2LnsT1uf1kZ2aC0BoZDC2sN9mx/z87V6Wzf2ILZ/soHKtSgx7/mHVokApVyrkFeUMNE1j/etf8+6sD0g9kE5k5Qiu79uGBtfUpV7L2sRUjUKn12HPKuTQziPs+/4QW1fv5OCOIwSH2+iXeCu3Pn4TJrOxvB+K8i/3j4e8EOIOYDJwBdBKSnleya1CXikPmqbx3cfbWbf0K75fvROf1w+AEILfvw50OkH9VnXo1L893R/qhNlqLq+SFeUUfyXkzzy94PztBvoAL1zgdhTlb6fT6WjbuxVte7fC4/aSsvs4R348dZ58neY1qdm4mjqpSakwLijkpZR7AbXavHLZMZmN1GtRm3otapd3KYryt1LL/ymKolRgf3okL4T4HIg7w0VJUsqPzveOhBAPAw8DVKtW7bwLVBRFUf66Pw15KWWXi3FHUsoXgRch8MbrxdimoiiKcm5quEZRFKUCu6CQF0LcKoRIBVoDnwoh1l6cshRFUZSL4UJn13wAfHCRalEURVEuMjVcoyiKUoGpkFcURanAVMgriqJUYCrkFUVRKjAV8oqiKBWYCnlFUZQKTIW8oihKBaZCXlEUpQJTIa8oilKBqZBXFEWpwFTIK4qiVGAq5BVFUSowFfKKoigVmAp5RVGUCkyFvKIoSgWmQl5RFKUCUyGvKIpSgamQVxRFqcBUyCuKolRgKuQVRVEqMBXyiqIoFZgKeUVRlApMhbyiKEoFpkJeURSlAlMhryiKUoGpkFcURanAVMgriqJUYCrkFUVRKrALCnkhxFwhxD4hxE9CiA+EEOEXqzBFURTlwl3okfx6oJGUsglwAHjiwktSFEVRLpYLCnkp5Toppa/02y1AwoWXpCiKolwshou4rYHAe2e7UAjxMPBw6bfFQoj9F/G+fxUN5PwN2/07qZr/fpdbvXD51Xy51QuXZ831/98bCCnlua8gxOdA3BkuSpJSflR6nSSgJdBH/tkG/0ZCiO1Sypbldf9/har573e51QuXX82XW73w76n5T4/kpZRd/uRO7wN6AJ3LM+AVRVGU013QcI0Q4kZgLNBRSum4OCUpiqIoF8uFzq5ZDIQA64UQu4QQz1+Emi7Ei+V8/3+Fqvnvd7nVC5dfzZdbvfAvqflPx+QVRVGUy5c641VRFKUCUyGvKIpSgVWokBdC6IUQPwghPinvWv6MECJcCPF+aVuIvUKI1uVd058RQgwXQuwRQuwWQrwjhLCUd01/JIR4VQiRJYTY/bufRQoh1gshDpZ+jijPGv/oLDVfsi1DzlTv7y4bJYSQQojo8qjtbM5WsxBiqBBif+nzek551XcmZ3leNBNCbCl9D3S7EKLVn22nQoU88Diwt7yLOE8LgTVSygZAUy7xuoUQ8cBjQEspZSNAD/Qr36rOaAlw4x9+lghskFLWBTaUfn8pWcLpNV/KLUOWcHq9CCGqAl2B4/90QedhCX+oWQhxPXAL0ERKeSUwrxzqOpclnL6f5wBTpJTNgIml359ThQl5IUQCcDPwcnnX8meEEKFAB+AVACmlR0ppL9+qzosBsAohDEAQcLKc6zmNlPIbIO8PP74FWFr69VKg9z9a1J84U82XcsuQs+xjgPnAGOCSm81xlpoHAbOklO7S62T944Wdw1lqlkBo6ddhnMdrsMKEPLCAwBNMK+9CzkMtIBt4rXR46WUhhK28izoXKWUagSOd40A6UCClXFe+VZ23SlLKdIDSz7HlXM//ayDwWXkXcS5CiF5AmpTyx/Ku5f9QD2gvhNgqhPhaCHF1eRd0HoYBc4UQJwi8Hv/0P7wKEfJCiB5AlpRyR3nXcp4MwFXAc1LK5kAJl94QwilKx7FvAWoCVQCbEOKe8q2q4ittGeID3irvWs5GCBEEJBEYPricGIAI4FpgNLBMCCHKt6Q/NQgYLqWsCgyndDTgXCpEyANtgV5CiBTgXaCTEOLN8i3pnFKBVCnl1tLv3ycQ+peyLsBRKWW2lNILrATalHNN5ytTCFEZoPTzJfVv+dn8rmXI3Zd4y5DaBP74/1j6GkwAdgohztTz6lKSCqyUAd8TGAW4pN4wPoP7CLz2AJYD/443XqWUT0gpE6SUNQi8GfiFlPKSPcqUUmYAJ4QQv3aU6wz8Uo4lnY/jwLVCiKDSo53OXOJvFv/OxwReHJR+/qgcazkvv2sZ0utSbxkipfxZShkrpaxR+hpMBa4qfZ5fyj4EOgEIIeoBJi79rpQngY6lX3cCDv7pLaSUFeoDuA74pLzrOI86mwHbgZ8IPNkiyrum86h5CrAP2A28AZjLu6Yz1PgOgfcMvATC5kEgisCsmoOlnyPLu87zqPkQcALYVfrxfHnXea56/3B5ChBd3nWexz42AW+WPp93Ap3Ku87zqLkdsAP4EdgKtPiz7ai2BoqiKBVYhRiuURRFUc5MhbyiKEoFpkJeURSlAlMhryiKUoGpkFcURanAVMgriqJUYCrkFUVRKrD/ATGBe64PqlByAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"u=np.linspace(3,18,200)\n",
"v=np.linspace(-2,4,200)\n",
"X,Y=np.meshgrid(u,np.exp(v))\n",
"Z2=L2(X,Y,[10,10,11,9,12])\n",
"CS=plt.contour(X,np.log(Y),np.exp(Z2-np.max(Z2)),levels=lev)\n",
"junk=plt.clabel(CS,[.0001,.001,.01],fmt=\"%1.4f\")\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"What is the posterior mode for sigma? It follows Inv-chi^2(n-1,s^2) which has mode (n-1)/(n+1)s^2, or in our case (2/3)s^2 or about .866."
]
},
{
"cell_type": "code",
<<<<<<< HEAD
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"u=list(map(sum,Z2))\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 16,
=======
"execution_count": 58,
>>>>>>> f8632e30fe8b89896ed0fbf258653883205028de
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
<<<<<<< HEAD
"1.437043700373185 [0.68945124 0.99000008 1.26005551 1.65286899 3.30679651]\n"
=======
"With Rounding\n",
"mu: 10.40618592964824 0.49106811542637807 [ 9.10552764 10.01005025 10.38693467 10.7638191 11.74371859]\n",
"sigma: 1.3555419607805625 0.5063610652735683 [0.61111977 0.90438284 1.18631786 1.55614414 3.30679651]\n"
>>>>>>> f8632e30fe8b89896ed0fbf258653883205028de
]
}
],
"source": [
<<<<<<< HEAD
"Psigma=list(map(sum,np.exp(Z2-np.max(Z2))))\n",
"p=Psigma/sum(Psigma)\n",
"s=np.random.choice(np.exp(np.linspace(-2,4,200)),p=p,size=10000)\n",
"print(np.mean(s),np.percentile(s,q=[2.5,25,50,75,97.5]))"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"1.371512604377929 [0.61111977 0.9320659 1.18631786 1.60377754 3.20858216]\n"
]
}
],
"source": [
"Psigma=list(map(sum,np.exp(Z1-np.max(Z1))))\n",
"p=Psigma/sum(Psigma)\n",
"s=np.random.choice(np.exp(np.linspace(-2,4,200)),p=p,size=10000)\n",
"print(np.mean(s),np.percentile(s,q=[2.5,25,50,75,97.5]))"
]
},
{
"cell_type": "code",
"execution_count": 180,
"metadata": {},
"outputs": [],
"source": [
"p=[1/3,1/3]\n",
"p.append([1/200.0]*198)\n"
=======
"print('With Rounding')\n",
"Pmu=list(map(np.sum,np.exp(Z1-np.max(Z1)).T))\n",
"p=Pmu/np.sum(Pmu)\n",
"mu_samples=np.random.choice(np.linspace(3,18,200),p=p,size=5000)\n",
"print('mu:',np.mean(mu_samples),np.var(mu_samples),np.percentile(mu_samples,q=[2.5,25,50,75,97.5]))\n",
"Psigma=list(map(np.sum,np.exp(Z1-np.max(Z1))))\n",
"p=Psigma/np.sum(Psigma)\n",
"sigma_samples=np.random.choice(np.exp(np.linspace(-2,4,200)),p=p,size=5000)\n",
"print('sigma:',np.mean(sigma_samples),np.var(sigma_samples),np.percentile(sigma_samples,q=[2.5,25,50,75,97.5]))\n"
>>>>>>> f8632e30fe8b89896ed0fbf258653883205028de
]
},
{
"cell_type": "code",
"execution_count": 55,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Ignoring Rounding\n",
"mu: 10.3954824120603 0.49760868662407504 [ 9.03015075 10.01005025 10.38693467 10.7638191 11.81909548]\n",
"sigma: 1.4352129437474488 0.5752559640752803 [0.68945124 0.99000008 1.26005551 1.65286899 3.30679651]\n"
]
}
],
"source": [
"print('Ignoring Rounding')\n",
"Pmu=list(map(np.sum,np.exp(Z2-np.max(Z2)).T))\n",
"p=Pmu/np.sum(Pmu)\n",
"s=np.random.choice(np.linspace(3,18,200),p=p,size=5000)\n",
"print('mu:',np.mean(s),np.var(s),np.percentile(s,q=[2.5,25,50,75,97.5]))\n",
"Psigma=list(map(np.sum,np.exp(Z2-np.max(Z2))))\n",
"p=Psigma/np.sum(Psigma)\n",
"s=np.random.choice(np.exp(np.linspace(-2,4,200)),p=p,size=5000)\n",
"print('sigma:',np.mean(s),np.var(s),np.percentile(s,q=[2.5,25,50,75,97.5]))"
]
}
],
"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.5"
}
},
"nbformat": 4,
"nbformat_minor": 2
}