{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Python Problem - the lighthouse\n", "\n", "(from D. Sivia's book, \"Data Analysis - A Bayesian Tutorial\"):\n", "\n", "\n", "\n", "A lighthouse is somewhere off a piece of straight coastline at a position $\\alpha$ along the shore and a distance $\\beta$ out at sea. It emits a series of short highly collimated flashes at random intervals and hence at random azimuths. These pulses are intercepted on the coast by photo-detectors that record only the fact that a flash has occurred, but not the angle from which it came. \n", "$N$ flashes have been recorded so far at positions $\\{x_k\\}$.\n", "\n", "Suppose $\\beta$ is given. Where is the lighthouse?\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Guided solution\n", "\n", "We need to estimate the parameter $\\alpha$. \n", "Let us start by writing the likelihood for this problem; since the flashes are thrown at random azimuths, we know that:\n", "$$P(\\theta_k | \\alpha, \\beta) = \\frac{1}{\\pi}.$$\n", "\n", "Moreover,\n", "$$\\beta \\tanh(\\theta_k) = x_k - \\alpha,$$\n", "and by changing variables we get\n", "$$P(x_k | \\alpha, \\beta) = \\frac{\\beta}{\\pi \\big[ \\beta^2 + (x_k - \\alpha)^2 \\big]}.$$\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Scientific computing and plotting packages\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", " \n", "# Likelihood definition\n", "def likelihood(x, alpha, beta):\n", " return beta / (np.pi * (beta ** 2 + (x - alpha) ** 2))\n", "\n", "# Parameters\n", "alpha = 30.0 # alpha appears here, only for simmulations purposes, we want to find the value of this parameter\n", "beta = 10.0 # beta is given\n", " \n", "#Compute and plot the likelihood\n", "x = np.linspace(-90, 90, 1001)\n", "plt.plot(x, likelihood(x, alpha, beta))\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The above likelihood is the a Cauchy or Lorentz distribution. \n", "We will sample from it so that we can have some synthetic data to work with.\n", "\n", "## Generate some synthetic data" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "from scipy.stats import cauchy\n", "samples = cauchy.rvs(loc = alpha, scale = beta, size = 1000)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In order to write the posterior we use Bayes theorem\n", "$$P(\\alpha | \\{x_k\\}, \\beta) \\propto \\prod_{k = 1}^N P(\\{x_k\\} | \\alpha, \\beta)$$" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD8CAYAAAB3u9PLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xt4XHd95/H3d0YXS9bNF1m+W3IsX+TgOIlxCOGaBOIQigubFAe2zXazm9KGlm7ZdpOF0kKfbJd2t0Af4NnmIRSasIQ0LdQLKQngAE1I4gtxEluyYvkq2bIlW7ZlW9Z1vvvHnDFCSNbIHumcmfm8nujRnDO/c+Y7zkgfnd/vd84xd0dERCQWdgEiIhINCgQREQEUCCIiElAgiIgIoEAQEZGAAkFERAAFgoiIBBQIIiICKBBERCRQEHYBEzF79myvra2d+IY7diS/X399RusREYm6HTt2nHD36nTaZlUg1NbWsn379olvaJb8fjnbiohkMTM7lG5bdRmJiAigQBARkYACQUREAAWCiIgEFAgiIgIoEEREJKBAEBERIMvOQxCJikTC2bKng+bjZ7l6QSVvq5+Npc53EclSCgSRCTpzYYDfeXQ7L+7vurju1lVz+OKHrmNaYTzEykSujLqMRCZgcCjB7z62gx2HTvGXH3gDr/35u/nkHav4YVMHH//HV3D3sEsUuWxpBYKZbTCzZjNrMbMHRnm+2My+FTz/kpnVDnvuwWB9s5ndNmz9QTN7zcx2mpmuKSFZ4e+fP8jP9p3kf7z/Ddy9fjHl0wr5T29dyh/ftoLvvdrO5leOhl2iyGUbNxDMLA58CbgdaADuNrOGEc3uBU65+zLgc8Bng20bgE3AamAD8OVgfynvdPe17r7uit+JyCTr6O7l8z98nZtXzuGudYt+6bmPvP0qrl5QwV99v5m+waGQKhS5MukcIawHWtx9v7v3A48DG0e02Qh8PXj8JHCLJUfYNgKPu3ufux8AWoL9iWSdv/vpfnoHE/zpe0f+PQTxmPEnt63kyOkLPLGtNYTqRK5cOoGwABj+CW8L1o3axt0HgTPArHG2deAZM9thZveN9eJmdp+ZbTez7Z2dnWmUK5J5p3v6+ebWw7zvmvnUzZ4+apu31s/m2sVVfPX5gyQSGkuQ7JNOIIw2l27kp32sNpfa9iZ3v45kV9T9Zva20V7c3R9293Xuvq66Oq1Leotk3Le2tdLTP8TvvH3pmG3MjN+6cQkHTpznZ/tOTmF1IpmRTiC0AcM7TBcCI0fOLrYxswKgEui61LbunvreAXwbdSVJRLk7T2xvZd2SGaycW3HJtrdfPY8ZpYV8c9vhKapOJHPSCYRtQL2Z1ZlZEclB4s0j2mwG7gke3wls8eT8u83ApmAWUh1QD2w1s+lmVg5gZtOBdwO7rvztiGTey62n2dd5nrvWLRy37bTCOO95wzy2NHXQ0z84BdWJZM64gRCMCXwUeBpoAp5w991m9hkze1/Q7BFglpm1AH8EPBBsuxt4AmgEvg/c7+5DQA3wnJm9AmwFvufu38/sWxPJjCd3tFFSGOeONfPTan/HmnlcGBji2T0a85LsktaZyu7+FPDUiHWfGva4F7hrjG0fAh4asW4/cM1EixWZakMJ5/u7jvGuhhrKitM7sf+GulnMLivme68d5Y418ya5QpHM0ZnKIpew49Apus73c9vquWlvE48Zt189ly17Ougd0DkJkj0UCCKX8MzuYxTFY7x9xcRmuN28cg69Awm2Hewav7FIRCgQRMbg7jzTeJybls1Ku7so5YalMykqiPGTZo0jSPZQIIiMofn4WQ539fDuCXQXpZQWFXBD3Ux+8roCQbKHAkFkDD8O/rq/eeWcy9r+7cur2dtxjiOnL2SyLJFJo0AQGcPzLSdYXlNGTcW0y9r+rfXJcYcXdNayZAkFgsgoegeG2Hawi5uWzb7sfdTPKaOqtJCtBxQIkh0UCCKj+PnhU/QOJHjLFQRCLGa8sXYm2w6eymBlIpNHgSAyiudbThCPGTcsnXVF+1lfO5MDJ87T0d2bocpEJo8CQWQUz7Wc5NpFVROebjrS+rqZAGzV+QiSBRQIIiOcuTDAa22nr2j8IGX1/ApKi+JsPaBAkOhTIIiMsONQFwmHN11hdxFAQTzG9UtmKBAkKygQREbYdvAUBTFj7aKqjOzv2sUzeP34WV0OWyJPgSAywo6Dp1i9oJKSonhG9rd2USUJh11HujOyP5HJokAQGaZ/MMErbadZt2RGxva5ZmHySOPVttMZ26fIZFAgiAyz6+gZ+gYTGQ2E2WXFLKgqYWerAkGiTYEgMsyO4CSy62szFwgAaxdV8YqOECTiFAgiw2w/1MXimaXMKb+86xeN5ZpFlbR2XeDkub6M7lckkxQIIgF3Z8ehUxntLkq55uI4wpmM71skUxQIIoFDJ3s4ca4/491FAFcvqCRmqNtIIk2BIBJ4uTUYP5iEI4TpxQXUzynXEYJEmgJBJPBK6xlKCuMsqy6blP2vnl9B41GdiyDRpUAQCbzadpqrF1RQEJ+cH4uG+RUc6+7VwLJElgJBBBgYSrD7aPfFk8gmw6p5FQA0tZ+dtNcQuRIKBBHg9eNn6RtMsGZh5aS9RioQGts1jiDRpEAQ4RfTQa+ZxCOEmdOLmFc5TeMIElkKBBGS4weVJYUsmVU6qa/TMK9CXUYSWQoEEZIzjNYsrMTMJvV1GuZX0NJ5jt6BoUl9HZHLoUCQvNc7METz8bOTOn6Q0jCvgqGEs/f4uUl/LZGJSisQzGyDmTWbWYuZPTDK88Vm9q3g+ZfMrHbYcw8G65vN7LYR28XN7GUz++6VvhGRy7X7aDdDCZ/UGUYpDfM1sCzRNW4gmFkc+BJwO9AA3G1mDSOa3QuccvdlwOeAzwbbNgCbgNXABuDLwf5SPgY0XembELkSqfsUTOaAcsqiGaWUFRdoYFkiKZ0jhPVAi7vvd/d+4HFg44g2G4GvB4+fBG6xZGfsRuBxd+9z9wNAS7A/zGwhcAfwlSt/GyKX77UjZ5hTXszcysxe4XQ0sZixcm45je0KBImedAJhAdA6bLktWDdqG3cfBM4As8bZ9vPAnwCJCVctkkGNR7tZHXTlTIUVc8tpPnYWd5+y1xRJRzqBMNq0i5Gf5LHajLrezN4LdLj7jnFf3Ow+M9tuZts7OzvHr1ZkAvoGh2jpOHexb38qrJhbTnfvIMe7dQkLiZZ0AqENWDRseSFwdKw2ZlYAVAJdl9j2JuB9ZnaQZBfUzWb22Ggv7u4Pu/s6d19XXV2dRrki6WvpOMdgwi+eRTwVlteUA9B8XOcjSLSkEwjbgHozqzOzIpKDxJtHtNkM3BM8vhPY4snj4c3ApmAWUh1QD2x19wfdfaG71wb72+Lu/z4D70dkQlKDu2EEwuvHFAgSLQXjNXD3QTP7KPA0EAe+6u67zewzwHZ33ww8AjxqZi0kjww2BdvuNrMngEZgELjf3XVGjkRGU/tZSgrj1M6aPmWvOXN6EdXlxTpCkMgZNxAA3P0p4KkR6z417HEvcNcY2z4EPHSJff8Y+HE6dYhkWmP7GVbMLScem9wzlEdaUVPO6woEiRidqSx5y91paj87pQPKKcuDQEgkNNNIokOBIHmr/UwvZy4MTOn4QcqKuWX0DiRoPdUz5a8tMhYFguSt1IByw7zyKX/tizONNLAsEaJAkLzV1N6NGayYO/VHCPWpmUYaR5AIUSBI3mps72bJzOS1haZaWXEBC2eU0KyrnkqEKBAkbzW1d4cyfpCyvKacvTpCkAhRIEheOtc3yMGTPTSEHAj7Os8xMKTLeUk0KBAkLzUfm/ozlEdaMbeMgSHn4InzodUgMpwCQfJSY3Bf4zDOQUhJzTTao5lGEhEKBMlLjUe7qSwpZN4U3ANhLFdVlxEzNI4gkaFAkLyUHFAuJ3kfp3BMK4yzZNZ0XtdMI4kIBYLknaGEs+dYNw3zKsMuhfo5ZbzeoSMEiQYFguSdgyfP0zuQYFUIZyiPtLymnEMne+gb1EWAJXwKBMk7Fy9ZEeKAckp9TRlDCeeAZhpJBCgQJO80tXdTEDOWzSkLu5Rf3CxH4wgSAQoEyTtN7d0sm1NGcUE87FKomz1dM40kMhQIknca27tDPUN5uGnB3dp0kTuJAgWC5JWT5/o43t0X6hnKI9XXlLFXXUYSAQoEyStNEThDeaTlNeUcPHleM40kdAoEyStN7eFfw2ik+ppyEg77OzXTSMKlQJC80tjeTU1FMTOnF4VdykXLa5KznTSOIGFTIEheaYrQgHJK3ezpxGOmcQQJnQJB8kbf4BAtHeci1V0EUFwQZ8msUh0hSOgUCJI39h4/x2DCIzWgnLJ8Tjl7O3SEIOFSIEjeiOKAcsrymjIOnTxP74BmGkl4FAiSNxrbu5lWGKN21vSwS/kVmmkkUaBAkLzR1N7NyrkVxGPh3QNhLKlrGu3VpbAlRAoEyQvuTuPR7kh2F0FyplFBzDSwLKFSIEheOHqml+7ewUgOKAMUFcSona27p0m4FAiSF5pS90CIwE1xxlI/p0xXPZVQpRUIZrbBzJrNrMXMHhjl+WIz+1bw/EtmVjvsuQeD9c1mdluwbpqZbTWzV8xst5l9OlNvSGQ0jcEMoxVzo3mEAMmB5UNdPZppJKEZNxDMLA58CbgdaADuNrOGEc3uBU65+zLgc8Bng20bgE3AamAD8OVgf33Aze5+DbAW2GBmb8rMWxL5VU3t3dTOKqWsuCDsUsa0vKYMd9jXqW4jCUc6RwjrgRZ33+/u/cDjwMYRbTYCXw8ePwncYmYWrH/c3fvc/QDQAqz3pNSnvjD48it8LyJjamyP7oByysWZRhpHkJCkEwgLgNZhy23BulHbuPsgcAaYdaltzSxuZjuBDuAH7v7SaC9uZveZ2XYz297Z2ZlGuSK/7FzfIIdO9kTuGkYj1c7STCMJVzqBMNqk7ZF/zY/VZsxt3X3I3dcCC4H1Znb1aC/u7g+7+zp3X1ddXZ1GuSK/bE+Ez1AerqggRp1mGkmI0gmENmDRsOWFwNGx2phZAVAJdKWzrbufBn5McoxBJOMuXrIiolNOh6uvKdPJaRKadAJhG1BvZnVmVkRykHjziDabgXuCx3cCW9zdg/WbgllIdUA9sNXMqs2sCsDMSoBbgT1X/nZEflVj+1kqSwqZXzkt7FLGVT+nnMNdPVzo10wjmXrjTrlw90Ez+yjwNBAHvuruu83sM8B2d98MPAI8amYtJI8MNgXb7jazJ4BGYBC4392HzGwe8PVgxlEMeMLdvzsZb1AkOaBcTnKeQ7Qtrym/ONPo6gWVYZcjeSatOXju/hTw1Ih1nxr2uBe4a4xtHwIeGrHuVeDaiRYrMlFDCaf5WDcfWr8k7FLSkrp72t6OswoEmXI6U1ly2oET5+kdSLAqwmcoD1c7ezqFcdPAsoRCgSA5Lcr3QBhNYTw500iXsJAwKBAkpzW2d1MYt4snfWWD+ppyHSFIKBQIktMaj3azbE45RQXZ81Gvn1NG6ynNNJKplz0/JSKXobG9O/JnKI+UmmnUonssyxRTIEjO6jjbS+fZvsjeA2Esw2caiUwlBYLkrKb25C/UbDtCWDJLM40kHAoEyVmNF2+Kk12BUBiPsXS2bpYjU0+BIDmrsb2bBVUlVJYWhl3KhNXXlPG6uoxkiikQJGc1Hj2TdeMHKctrymntukBP/2DYpUgeUSBITurpH2T/ifNZ112UUj8nObCsmUYylRQIkpOaj53Fnaw9QqjX3dMkBAoEyUmN7dk5oJxSO6uUonhM4wgypRQIkpMaj3ZTPq2AhTNKwi7lshTEYyytnq4jBJlSCgTJSakzlLPhHghjSV7TSEcIMnUUCJJzhhLOnvazWTt+kLJ8Thltpy5wvk8zjWRqKBAk5xw8eZ4LA0NZO36QkhpY1kwjmSoKBMk5F89QzvIjhPrgmkbqNpKpokCQnJO6B0L9nOy5B8JolsxMzjTSEYJMFQWC5JxsvAfCaFIzjXSEIFMlu39iREaRjfdAGMty3T1NppACQXJKR3fyHgirs3z8IGV5TRlHTmumkUwNBYLklNeOnAHgDQsrQ64kMy5ewkLjCDIFFAiSU15tO4NZ9l6yYqTURe40jiBTQYEgOWXXkTNcVV3G9OKCsEvJiCWzplNUoJlGMjUUCJJTXjtyhjcsyI3uIoB4zLiqukxHCDIlFAiSMzq6e+k425dTgQDJgWVd5E6mggJBckauDSinLK8p58jpC5zTTCOZZAoEyRm5NqCcsjyYadR8TN1GMrnSCgQz22BmzWbWYmYPjPJ8sZl9K3j+JTOrHfbcg8H6ZjO7LVi3yMyeNbMmM9ttZh/L1BuS/JVrA8opqWsyNR49E3IlkuvGDQQziwNfAm4HGoC7zaxhRLN7gVPuvgz4HPDZYNsGYBOwGtgAfDnY3yDwcXdfBbwJuH+UfYpMSK4NKKfMr5xGVWkhu4OL9olMlnSOENYDLe6+3937gceBjSPabAS+Hjx+ErjFkncm2Qg87u597n4AaAHWu3u7u/8cwN3PAk3Agit/O5KvcnVAGcDMWD2/QoEgky6dQFgAtA5bbuNXf3lfbOPug8AZYFY62wbdS9cCL6Vftsgvy9UB5ZTV8ytpPnaWgaFE2KVIDksnEEa7B6Gn2eaS25pZGfBPwB+6+6h//pjZfWa23cy2d3Z2plGu5KNcHVBOWT2/gv6hhKafyqRKJxDagEXDlhcCR8dqY2YFQCXQdaltzayQZBh8w93/eawXd/eH3X2du6+rrq5Oo1zJR7k6oJySuljfbg0syyRKJxC2AfVmVmdmRSQHiTePaLMZuCd4fCewxd09WL8pmIVUB9QDW4PxhUeAJnf/m0y8Eclf7p6zA8opdbPLKCmMaxxBJtW4f065+6CZfRR4GogDX3X33Wb2GWC7u28m+cv9UTNrIXlksCnYdreZPQE0kpxZdL+7D5nZW4DfBF4zs53BS/13d38q029Qcl/7meSA8tpFVWGXMmniMWPVvPKLtwcVmQxpHV8Hv6ifGrHuU8Me9wJ3jbHtQ8BDI9Y9x+jjCyIT9vLh0wA5HQiQHFj+9stHSCScWEw/PpJ5OlNZst7O1lMUFcRYlaMDyimr51dwrm+Qw109YZciOUqBIFnv5cOnuXp+RdbfQ3k8q+cnx0g0jiCTJbd/giTnDQwleO3IGa5dPCPsUibd8rllFMSMXZppJJNEgSBZbU/7WfoGEzk/fgBQXBBn5bxyXm07HXYpkqMUCJLVdraeAuDaxbkfCADXLKzi1dYzJBIjzw0VuXIKBMlqLx8+zeyyYhZUlYRdypRYu6iKs32D7D+hM5Yl8xQIktV2tp7m2sVVJM91zH2prrGdrRpHkMxTIEjWOt3Tz/4T5/Ni/CBlaXUZZcUFF7vKRDJJgSBZa2drcnA1X8YPIHnG8pqFlbyiIwSZBAoEyVovHz6NGaxZmD+BAHDNoiqa2rvpHRgKuxTJMQoEyVrbD3Wxcm4FZTl6hdOxrF1UxWDCdYKaZJwCQbLSwFCClw+fZn1t7p+QNlJqzOSVVp2PIJmlQJCs1Hi0m57+Id5YNzPsUqZcTcU05lZMuziGIpIpCgTJStsOdgGwvjb/AgGSRwkKBMk0BYJkpa0Hulgyq5Q5FdPCLiUU1y2p4nBXDx1ne8MuRXKIAkGyjruz/dAp3pinRwfAxfe+7YDOR5DMUSBI1tnXeY6u8/15210EcPWCSkoK4xe7zkQyQYEgWWdr8FdxPg4opxTGY1y7uIqtBxQIkjkKBMk62w52MbusmNpZpWGXEqo31s6k6Vg33b0DYZciOUKBIFln64Eu1tfNyJsL2o1lfd1M3GHHIY0jSGYoECSrHD7Zw5HTF3jT0llhlxK6axdXURAztqnbSDJEgSBZ5fl9JwB481WzQ64kfKVFBaxeUKmBZckYBYJkledaTlBTUcxV1dPDLiUSbqibySutZ3ShO8kIBYJkjUTCeWHfSW5aNjvvxw9SbrxqFv1DCR0lSEYoECRr7Dl2lq7z/dyk7qKLbqibSVE8xnN7T4RdiuQABYJkjZ8F4wc3LVMgpJQWFXDdkir+TYEgGaBAkKzxXMsJllZPZ25lfl6/aCxvra+msb2bE+f6wi5FspwCQbJC/2CCrQe61F00ircER0zPt+goQa6MAkGywraDXfT0D/G25dVhlxI5Vy+opLKkUOMIcsUUCJIVtuzpoKggxk3LdELaSPGYcdOyWTzXcgJ3D7scyWJpBYKZbTCzZjNrMbMHRnm+2My+FTz/kpnVDnvuwWB9s5ndNmz9V82sw8x2ZeKNSG57dk8Hb1o6i9Ki/Lp/crreWl9N+5le9nacC7sUyWLjBoKZxYEvAbcDDcDdZtYwotm9wCl3XwZ8DvhssG0DsAlYDWwAvhzsD+BrwTqRSzp44jz7T5zn5hXqLhrLO1fMAeCHTcdDrkSyWTpHCOuBFnff7+79wOPAxhFtNgJfDx4/CdxiyTOHNgKPu3ufux8AWoL94e4/BXQ2jYxry54OAG5eWRNyJdE1t3IaaxZW8oNGBYJcvnQCYQHQOmy5LVg3aht3HwTOALPS3PaSzOw+M9tuZts7OzsnsqnkiGebO7iqejqL8/xy1+O5dVUNO1tP67aactnSCYTRrhEwcuRqrDbpbHtJ7v6wu69z93XV1eoyyDfn+wZ5aX8XN6+cE3Ypkfeuhhrck+MtIpcjnUBoAxYNW14IHB2rjZkVAJUku4PS2VZkTM82d9A/lOCWVeouGs/KueUsqCpRt5FctnQCYRtQb2Z1ZlZEcpB484g2m4F7gsd3Als8Of9tM7ApmIVUB9QDWzNTuuSDf33tGLPLii7eVF7GZma8q6GGf9t7gnN9g2GXI1lo3EAIxgQ+CjwNNAFPuPtuM/uMmb0vaPYIMMvMWoA/Ah4Itt0NPAE0At8H7nf3IQAz+ybwArDCzNrM7N7MvjXJdhf6h9iyp4PbVs8lHtPVTdNxx5p59A0m+EHjsbBLkSyU1qRud38KeGrEuk8Ne9wL3DXGtg8BD42y/u4JVSp55yevd3BhYIj3vGFe2KVkjesXz2B+5TQ27zzK+69dGHY5kmV0prJE1lOvHWNGaSE31Km7KF2xmPFr18zn3/ae4NT5/rDLkSyjQJBIutA/xI+ajnPb6rkUxPUxnYhfu2Y+gwnnX3ep20gmRj9pEknPNB7jfP8QG9dO6LQVAVbPr2Bp9XS+s/NI2KVIllEgSCT908+PsKCqRN1Fl8HM+HfXLWTrgS72d+raRpI+BYJEzvHuXp7b28kHrltATLOLLstd6xZSEDMe39Y6fmORgAJBIuc7Lx8h4fCB6zRL5nLNKZ/GratqeHJHG32DQ2GXI1lCgSCR4u784442rltcRd3s6WGXk9XuvmExXef7eWa3zlyW9CgQJFJe2HeSlo5zfOiGJWGXkvXeumw2i2aW8A8vHAy7FMkSCgSJlH944RAzSgt57xqdjHalYjHjP95Ux7aDp9hx6FTY5UgWUCBIZBw5fYFnGo/xwTcuZlphfPwNZFy/sW4RlSWFPPzTfWGXIllAgSCR8Y0XDwHw4RsWh1xJ7pheXMBv3biEZxqPawqqjEuBIJHQ3TvAYy8e4t0Nc1k0UzfCyaR73lxLUTzGF7e0hF2KRJwCQSLh0RcO0d07yP3vXBZ2KTlndlkx/+HNtXx75xGa2rvDLkciTIEgoevpH+SR5w7wjhXVvGFhZdjl5KTffcdVlBcX8NdPN4ddikSYAkFC99iLh+g638/v36yjg8lSVVrER95xFVv2dPDCvpNhlyMRpUCQUJ3u6eeLW1p4a/1srl+i6xZNpt9+cx0LZ5Twye+8prOXZVQKBAnV3/6ohXN9g3zijlVhl5LzSori/MWvX82+zvP83U/2h12ORJACQUKz9/hZHn3xIL+xbhEr51aEXU5eeOeKObx3zTy+uKWF14+fDbsciRgFgoQikXAe+OfXmF5cwH+9bUXY5eSVP/u11VSUFHD/N37OhX51HckvKBAkFI+9dIgdh07xp3c0MLusOOxy8kp1eTF/8xtr2dtxjk//v91hlyMRokCQKdd87CwPfa+Jty2v5gPX6Y5oYXjb8mp+7x1X8fi2Vv7++QNhlyMRURB2AZJfzvcNcv///TkVJYX877uuwUw3wAnLx9+9gn2d5/jMdxuZV1nChqvnhl2ShExHCDJlhhLOH3zzZfZ3nuPzH1xLdbm6isIUjxmf/+C1XLOwij/45sv8sFH3Tch3CgSZEu7On23exY/2dPDp963mpmWzwy5JSE5F/fpvr2fV/Ao+8tgO/mXnkbBLkhApEGTSJRLOJ7+zi8dePMzvvH0pv3ljbdglyTCVpYU8du96rlsyg489vpO/fnoPQwkPuywJgQJBJtWF/iH+yxM7+cZLh/nI26/igQ0rwy5JRlE+rZBH713P3esX8aVn9/Hhr7zI4ZM9YZclU0yBIJPm0Mnz3Pl/fsbmV47yx7et4L9tWKFB5AgrLojzlx9Yw1/duYbdR7q57fM/5eGf7qN3QOcq5Atzz55Dw3Xr1vn27dsnvmHql1AWvddsNjiU4Gs/O8j/eqaZwniMv910Le9cOSfssmQCjp6+wCe/s4stezpYUFXC79+8jF+/doHuZJeFzGyHu69Lq60CQTJlKOFsfuUIX/jhXg6e7OGWlXP4i1+/mvlVJWGXJpfp+ZYTfPb7e3i17QxVpYV88I2L2HjNAlbNK9fRXpbIeCCY2QbgC0Ac+Iq7/88RzxcD/wBcD5wEPujuB4PnHgTuBYaAP3D3p9PZ52gUCNHj7uzrPM/mnUd4Ynsbx7p7WTWvgj9613JuXTVHvzRygLvz0oEuvvb8QZ5pPEbCoW72dN7VUMONS2exrnYG5dMKwy5TxpDRQDCzOPA68C6gDdgG3O3ujcPa/B6wxt0/YmabgPe7+wfNrAH4JrAemA/8EFgebHbJfY5GgRC+RMLZf+Icr7SeYWfraX66t5NDJ3swg7cvr+ZD6xdz66oaYjEFQS46ca6Pp3cf46nX2tl6oIuBISdmsLymnJVzy1k5r4LlNWUsmlHKghkllBbp3NewTSQQ0vm/tR5ocff9wc4fBzYCw395bwT+PHj8JPBFS/5puBF43N37gANm1hJQ1Q+nAAAHZElEQVTsjzT2KZPI3RkYcnoHh+gbSNA3OETvQILegSG6ewc4dX6AUz39nDrfz8nz/Rzu6uFwVw+tXT30DSYAmF4U54als/jPb13KLavmMK9SXUO5bnZZMR++YQkfvmEJF/qHePnwKV7cf5JdR7vZeqCL7+w8+kvtZ5QWMr+qhJnTi5hRWsSM0kKqgu9l0wopKYxTUhSjpLCAkqI4JYVxSoviFBXEiMeMgpgRjxmF8V8s66hz8qQTCAuA1mHLbcANY7Vx90EzOwPMCta/OGLb1MVrxttnXnn3537CYMIh+R/uHnwHx0kkfwf/ynof3n60bUdZn3CnbzCR9gFTWXEBC2eUcFX1dN65oprlNeWsXVTF0uoy4joSyFslRXHevGw2bx52kuGZngFaOs/SduoCR05f4MipCxw9fYGungEOd/Vw6nw/3b2DV/S68SAkUmERM8MMDDCz4DvA8PVgwXIsCBSzX14/fHsm+LG+nJ+CiQTbzNIinvjIjZfxKhOTTiCMVvXIXyVjtRlr/WjTXUf99WRm9wH3ASxevHjsKrPc8pryi/9gwz/URvAB/pUP7rBlG2M9w5/7xbpYzJhWEKO4ME5xQSz5dfFxnIqSgot/0VWVFlJcoJklkp7K0kKuXzKT65eM3WZwKMHpCwP09A1xYSD51dM/SO/AEBf6E/T0DzIw5AwlEsF3ZzAx9vLYfyQBqeURzyV8rD++kssTcVkd0RPcqHza1HS9pfMqbcCiYcsLgaNjtGkzswKgEugaZ9vx9gmAuz8MPAzJMYQ06s1KX/zQdWGXIDIlCuKx5CXPy8KuREZK58S0bUC9mdWZWRGwCdg8os1m4J7g8Z3AFk/G7GZgk5kVm1kdUA9sTXOfIiIyhcY9QgjGBD4KPE1yiuhX3X23mX0G2O7um4FHgEeDQeMukr/gCdo9QXKweBC4392HAEbbZ+bfnoiIpEsnpomI5LCJTDvVtYxERARQIIiISECBICIigAJBREQCCgQREQGybJaRmXUChyawyWzgxCSVMxlU7+TKpnqzqVZQvZPtSupd4u7V6TTMqkCYKDPbnu50qyhQvZMrm+rNplpB9U62qapXXUYiIgIoEEREJJDrgfBw2AVMkOqdXNlUbzbVCqp3sk1JvTk9hiAiIunL9SMEERFJU84FgpndZWa7zSxhZutGPPegmbWYWbOZ3RZWjSOZ2YagphYzeyDsekYys6+aWYeZ7Rq2bqaZ/cDM9gbfZ4RZ43BmtsjMnjWzpuCz8LFgfSRrNrNpZrbVzF4J6v10sL7OzF4K6v1WcKn4yDCzuJm9bGbfDZYjW6+ZHTSz18xsp5ltD9ZF8vMAYGZVZvakme0JPsc3TkW9ORcIwC7gA8BPh680swaSl+VeDWwAvmxmod8KLKjhS8DtQANwd1BrlHyN5L/ZcA8AP3L3euBHwXJUDAIfd/dVwJuA+4N/06jW3Afc7O7XAGuBDWb2JuCzwOeCek8B94ZY42g+BjQNW456ve9097XDpm9G9fMA8AXg++6+EriG5L/z5NebvOdu7n0BPwbWDVt+EHhw2PLTwI0RqPNG4Omx6ozKF1AL7Bq23AzMCx7PA5rDrvEStf8L8K5sqBkoBX5O8h7jJ4CC0T4nYX+RvMvhj4Cbge+SvDtrlOs9CMwesS6SnwegAjhAMMY7lfXm4hHCWBYArcOW24J1YYtqXeOpcfd2gOD7nJDrGZWZ1QLXAi8R4ZqD7pedQAfwA2AfcNrdU3ekj9rn4vPAnwCJYHkW0a7XgWfMbEdwn3aI7udhKdAJ/H3QJfcVM5vOFNQ7NXduzjAz+yEwd5SnPuHu/zLWZqOsi8IUq6jWlfXMrAz4J+AP3b3bbLR/6mjw5J0E15pZFfBtYNVozaa2qtGZ2XuBDnffYWbvSK0epWkk6g3c5O5HzWwO8AMz2xN2QZdQAFwH/L67v2RmX2CKurOyMhDc/dbL2KwNWDRseSFwNDMVXZGo1jWe42Y2z93bzWweyb9sI8PMCkmGwTfc/Z+D1ZGuGcDdT5vZj0mOfVSZWUHwV3eUPhc3Ae8zs/cA00h2cXye6NaLux8NvneY2beB9UT389AGtLn7S8HykyQDYdLrzacuo83AJjMrNrM6oB7YGnJNANuA+mCGRhHJge/NIdeUjs3APcHje0j200eCJQ8FHgGa3P1vhj0VyZrNrDo4MsDMSoBbSQ4iPgvcGTSLTL3u/qC7L3T3WpKf1y3u/mEiWq+ZTTez8tRj4N0kJ59E8vPg7seAVjNbEay6heR96Se/3rAHUCZhQOb9JBO2DzjOLw/YfoJk32wzcHvYtQ6r6z3A60Ftnwi7nlHq+ybQDgwE/7b3kuwz/hGwN/g+M+w6h9X7FpLdFa8CO4Ov90S1ZmAN8HJQ7y7gU8H6pST/aGkB/hEoDrvWUWp/B/DdKNcb1PVK8LU79TMW1c9DUNtaYHvwmfgOMGMq6tWZyiIiAuRXl5GIiFyCAkFERAAFgoiIBBQIIiICKBBERCSgQBAREUCBICIiAQWCiIgA8P8Bh11vqj2OC2EAAAAASUVORK5CYII=\n", "text/plain": [ "