{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# 11. Sampling Methods" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "%matplotlib inline\n", "\n", "from prml.rv import Gaussian, Uniform\n", "from prml.sampling import metropolis, metropolis_hastings, rejection_sampling, sir\n", "\n", "np.random.seed(1234)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 11.1.2 Rejection sampling" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAD4CAYAAADFAawfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAgAElEQVR4nO3dd3xc1Z3//9e5M6Pee7XlXuUChmCb4gAhhCRAQqiBkIQNhG/yDaTtfnezqbsp+9gUkmx+7LKBhU1IDIYQCMUYY4NtjJvAsmy5yZIlW1bvXTP3nt8fIwnLluyRNDP3Sv48H495WB7NnPuWPfPRnXNPUVprhBBCOJdhdwAhhBDnJoVaCCEcTgq1EEI4nBRqIYRwOCnUQgjhcO5QNJqWlqYLCgpC0bQQQkxJRUVFjVrr9JG+F5JCXVBQwJ49e0LRtBBCTElKqcrRviddH0II4XBSqIUQwuGkUAshhMNJoRZCCIeTQi2EEA4nhVoIIRwuJMPzhBCTQ3t7O/X19Xi9XrujTGkej4eMjAwSEhLG9Xwp1GJyq9oB1e9B9lLIWQ4RMXYnmjTa29upq6sjNzeX6OholFJ2R5qStNb09PRQXV0NMK5iLYVaTE7VRbDpx3DszaG7tDLw5lyC5/YnUQk5NoabHOrr68nNzSUmRn65hZJSipiYGHJzczl16pQUanGBeOP78M4j9BqxbDPWUGLNIZMGpqs6VlTvofs/riTmwY2o5AK7kzqa1+slOjra7hgXjOjo6HF3MUmhFpPL+3+Edx5hr1HIer2GPjygoJ14juqZlFLA3f0v0PO7q4i6fwNGxjy7EzuadHeEz0T+rWXUh5g8TuxGv/x1Ko0CXrKu8RfpM5xS2TzJrWhfP33/dQ2+1mobggoRXFKoxeTQXoN+5m46iGMdH0er0V+69Sqdp/gMHrObuj/cH8aQQoSGFGoxObz4FczuFp5RN9OlI8/78AaVxnZWkNu0jZb3/xaGgEKEjhRq4XwVW+HYm7ylVnHKTAr4aVu5lBYS4ZVvoH19IQwonMiyLLsjBI0UauFsWsOmf6HLlchOc/GYnupTHl7jwyT76qn5y3dDFFA4xec+9znmzJnDQw89xPe+9z3uu+++Yd/XWrNs2TKeeuqpgNv86le/elY7dgi4UCulXEqp95VSL4cykBDDHN0AJ3byNpfhU2dfPDzv09VMDjGLtNIn6KmvCEFA4QQtLS3s2bOHkpISqqur+d3vfscXv/jFYY959tlnaW5u5q677gq43W9961s8/fTTlJWVBTvymIzljPoh4GCogghxFstCv/kj2lypFJnzx93MBq7Cg5eaF38YxHDCSZKTkyktLSUqKornnnuOpqYmrrjiimGP+c1vfsM999yDxxP4L/yCggIuv/xyHn300WBHHpOACrVSKg/4OPD70MYR4jSlL6Dq9rOZVVjKNe5mWlQSh5hFVvV6vN1tQQwonCIlJYXvfve7fPaznyU1NZWUlBQeeOAB+vr81ybKysrYvn07n/nMZ4Y975lnnkEpNeKtvb0dgFtuuYWnn37a1j7vQM+oHwH+Hhg1qVLqfqXUHqXUnoaGhqCEExcwrdFbfk6TK4Nic/aEm9vJRcTQQ/VrjwQhnHCSqqoqWlpa+NWvfkVKSgrr1q3j4Ycf5rHHHuN3v/sdAG+++SaxsbEsXbp02HOvvPJK3n333aHbK6+8QkpKCjfccMPQVO9Vq1ZRV1dHSUlJ2H+2QeedmaiU+gRQr7UuUkqtGe1xWuvHgMcAVqxYoYOWUFyYTu5G1Zey0309BGH2XCV51JBO/IE/oD/1XZQh19GniuLiYsB/4e9nP/sZAFdffTVbtmxh/fr1fOMb36CoqIgFCxZgnPH/np2dTXZ2NuBfpOqaa65h8eLFrFu3bugxixYtwuVysWvXrrMKfbgEMoV8NXCjUuoGIApIUEr9UWt9d2ijiQta0VN4VSTFvtkQjFnOSrFTX8TN1utUv/Mncq+Ql+9Ifvi3A5Searfl2AtzEvj+JxeN+XnFxcXExMTw7W9/e9j9c+fOZdeuXQDU1taSlpY2ahtdXV3ccMMNGIbByy+/PGyhKrfbTVJSErW1tWPOFiznPa3QWv+j1jpPa10A3AFskiItQqq3Db3/OUqYR7+KCFqz+5lHJzFY7/xH0NoU9isuLubKK68kNTV12P21tbXk5PhXUezt7SUycuSJUr29vdx44410dXWxfv164uPjz3pMZGQkvb29wQ8fIFmUSThPyTqUr5e9rqXnuCoydqZys0cvZU3vuzQf3UXKnEuD1/gUMZ4zWrsVFxefNcKjo6ODTZs28cMf+kf6pKSkjHhG7PV6ueWWW6ipqeHtt98mOTl5xGO0traSkpIS/PABGlNHndb6La31J0IVRgi0Rhc9Sb2RxQlz9I+q41VEIRaKxk3/X9DbFuHX1dXFsWPHOHbs2LD7f/Ob3wBw7733AjBv3jwqKoaPozdNkzvvvJPDhw+zceNG0tPTRzxGQ0MD3d3dzJ07NwQ/QWDkjFo4y6n3UbUlvOe+DnTwl+DsVHEc1/mk127GMk0M1/iH/Qn7lZSUYFkWx44d40c/+hFXXHEFGzZs4Be/+AV//OMfSUryLzmwevVqfvSjH9HQ0DBUkB988EFeffVVnnzySaqqqqiqqgIgNjaWwsLCoWPs2bMHpRSrVq0K/w84QC59C2cpehKfimCvb07IDlHCfJJ1K6f2yCTbya64uJiEhAReeeUVXnzxRa6//npeeeUV1q1bx2233Tb0uDVr1pCSksL69esB/3TytWvX0tPTw+23387KlSuHbt/5zneGHWP9+vVcddVVZ/WBh5OcUQvn8PWhD7xAKXPoU+dfIW+8DjKbj/MmPbv+AB+6KWTHEaFXXFxMYWEhS5YsoaioaNTHRUREcPfdd7N27VruueeeYRNazsU0TZ5//vmhYX92kTNq4Rzlb6P62jnkWhjSw/SpKI4yg+ym7fj6ZVW9yay4uDjgsc3f/va32bx5M0eOHAm4/XXr1hEdHc0dd9wx3ohBIYVaOEfpi/Qb0RzxZYf8UCXMJ44uqrc/E/JjidDQWlNSUsKSJUsCenxeXh5PPPEENTU1YzrG448/jtttb+eDdH0IZzC96EMvc0jPwFShf1keZQZ9ROAt+hOs+VzIjyeCL9Dui9ON9cz4zjvvHNPjQ0XOqIUzVGxB9bZy2LUgLIfzKQ8HmU1ex3v0dspCTcLZpFALZyh9Ea+K5LAvN2yH3M98oujj1NY/hO2YQoyHFGphP9OHPvQKRwhPt8egCvLpJRLf/pfCdkwhxkMKtbBf1XZUdyOHXePfHGA8LOXiKAXkdJXgldEfwsGkUAv7lb6IT0Vw0Jcf9kMfYSZxdHNqt0x+Ec4lhVrYy7LQB//GUQrGtSfiRJUxAwtFT/Ffwn5sIQIlhVrYq3YfqrOOY257FrzpVVFUkUtKwy60lv0uhDNJoRb2KnsDgEPe8I32ONMRZpKh66k7vNu2DEKcixRqYSt9dCO1RjZdKta2DEeYCUDb7mdtyyDEuUihFvbpboaTuyhjhq0xmlQKTSQRVbXZ1hxCjEYKtbBP+WaUtihT9hZqgCPMItdbQVtDtd1RhDiLFGphn6Mb6TViqDJH3lkjnA4zEzcmDTuk+2MyWrt2LZGRkXi93gm1o7Vm2bJlPPXUUwE9/qtf/Sr33XffhI4ZCCnUwh6WhT76BmV6OlrZ/zKsIpdeIjAPv253FDEOxcXFLFy4EI9nYkM8n332WZqbm7nrrrsCevy3vvUtnn76acrKyiZ03POx/x0iLky1xajuBircs+xOAoBWBsfJJ6OzFNM07Y4jxqi4uJjly5dPuJ3f/OY33HPPPQEX/IKCAi6//HIeffTRCR/7XKRQC3sc3YhGccibZ3eSIeVMJ5k2akvftTuKGKO9e/eybNmyob9v2rSJ1NRUHn74YUzTpK2tjQceeICUlBSSkpL4xje+wU9/+lNmz5499JyysjK2b9/OZz7zmWFtP/PMMyilRry1t7dzyy238PTTT2NZVsh+PlmPWthCH91ArZFNt46xO8qQY0wHoH3vi+QWXm5zGhGohoYGampqhgr1b3/7W/7hH/6BX//613zpS1/C6/Vy/fXX09XVxX/+538SGxvLP//zP1NXV8fKlSuH2nnzzTeJjY09a8eYK6+8knff/eCXd3NzM/fccw+XXXYZCQkJrFq1irq6OkpKSgLebWaspFCL8Otpheo9HFOX2Z1kmGaSaCWeiBPb7I5in9f+H9SW2HPsrEL42Nj3JiwuLkYpxcKFC/nSl77EX/7yF1599VXWrFkDwM9//nPKyso4fPgwKSkpACQnJ7N69ephu8MUFRWxYMECDGN4R0N2djbZ2f5dh9rb27nmmmtYvHgx69atA2DRokW4XC527dolhVpMIce3obRFuTEdnDRrWynK9XQW9h+lp6uD6Nh4uxOJAOzdu5eUlBQ+/elP09DQwM6dO4e6NCzL4pFHHuErX/nKUJEGmDnTP8np9MJaW1tLWlraqMfp6urihhtuwDAMXn75ZWJi/J8G3W43SUlJ1NbWhuLH8x8jZC0LMZqKt/GqCCrNDFB2hxmunOlcxH4q3nuNGVfcZnec8BvHGa3diouL0VqzdetWNm7cOKzfef/+/dTX13PttdcOe87JkyeB4YW6t7d3qPieqbe3lxtvvJGuri42bdpEfPzwX+KRkZH09vYG60c6i1xMFGGny9+milws5bI7ylnKmYYGeg+8ZncUEaDi4mK+9rWvcffdd/OFL3yBhoaGoe9VV/snMOXlDb9ovXHjRhISEigoKBi6LyUlhdbW1rPa93q93HLLLdTU1LBhwwaSk5PPekxra+uwM/Zgk0Itwqu9BtV4mONGgd1JRtSjoqkhk7j6PXZHEQHo6+vj0KFDLFmyhMcee4y0tDRuvfVWfD4fAKmpqQBUVFQMPaexsZFf/vKXFBYWotQHH+nmzZs37HEApmly5513cvjwYTZu3Eh6+tmTsxoaGuju7mbu3NCtACmFWoRXxRYAygn/JgGBKmcaOdZJmmur7I4izqO0tBSv10thYSHR0dG88MILHDhwgK9//euAv2sjJyeHhx9+mFdeeYVnnnmG6667DtM0z7rwt3r1aqqqqoadkT/44IO8+uqr/OQnP6GqqoodO3awY8cOSko+uOC6Z88elFKsWrUqZD+nFGoRXhVv02vEcMpMtTvJqI4xHRcWTUV/tTuKOI/i4mJiY2OZNcs/cWr69Ok888wzPProozz55JNERkby3HPPAXDrrbfyq1/9ikceeQTTNIeN+ABYs2YNKSkprF+/HvBPJ1+7di09PT3cfvvtrFy5cuj2ne98Z+h569ev56qrrho6ew8FKdQifLRGl79Fuc4D5bCriKc5QQ5e3PiOvGl3FHEen//85+ns7BzWhXH11Vfj8/n4/Oc/D8DKlSspLi6mu7ubHTt2kJ+fT1tb21ln1BEREdx9992sXbsWYGhCi9b6rNtLL/k3RDZNk+effz7k631IoRbh01yOaq+mymX/annnYio3J8ghpf2A7PoyBe3duxfDMCgsLDzre9/+9rfZvHkzR44cCaitdevWER0dzR133BHsmMNIoRbhU/4WAGWmfbu5BOo4eWTqBuorD9kdRQRZcXExM2fOJDb27M0q8vLyeOKJJ6ipqQmoLa01jz/+OG53aEc6yzhqET4Vb9NhJNJkJTpu/PSZjg9c7Gzd+wqZBQtsTiOC6Qc/+AE/+MEPRv3+WM6O77zzziAkOj85oxbhYVno8i2U63xH908PqiaLftxY5W/bHUUIKdQiTOoPoHpbqHIX2J0kIJZycYJcUjtKQ7oqmhCBkEItwuO4f6GjY75sm4ME7jj5ZOhGGo4ftDuKuMBJoRbhcXwbbUYybUyehY6O45923Lr3ZZuThI6MagmfifxbS6EWoWdZ6MrtVOgcu5OMySky6ceDNTCbcqrxeDz09PTYHeOC0dPTM+6tws5bqJVSUUqpXUqpYqXUAaXUD8d1JHHhajiI6mnm5CTpnx5kKReV5JLWcXBK9lNnZGRQXV1Nd3e3nFmHkNaa7u5uqqurycjIGFcbgQzP6wOu1lp3KqU8wDal1Gta6x3jOqK48Bx/B5hc/dODKsljDtuoLd9P1uwl53/CJJKQkADAqVOnJrx7tzg3j8dDZmbm0L/5WJ23UGv/r9rOweMN3OTXrwjc8a20G0m06vG9SMPtlBnPe95cFrjryXH5x1O3Fb885Qo1+Iv1eIuHCJ+AJrwopVxAETAb+J3WeucIj7kfuB9g2rRpwcwoJjOtJ03/dI92s8ubT7mZikKz3TudTKODPhWBrthqdzxxAQvoYqLW2tRaLwPygEuVUotHeMxjWusVWusVI63ZKi5QDYdR3Y1UuwrsTnJOfdrFX3sXcdxM5qLIOn59XSpul8Hb3tlUkUNK55Ep2U8tJocxjfrQWrcCm4HrQxNHTDnH/WeiZaaz+6ePman04uFzeS3877dv48arV/LjTxVSb8Wy05xPBo3UH5d1P4Q9Ahn1ka6UShr4Ohr4CCCvWBGYynfoMBJpcXD/tNZw1JdGmqub+276MHFxcQB8+uJpXDMnib/6LgWgrUS25xL2COSMOhvYrJTaB+wG3tBaT90ZACJ4tEZXbOO4znH0+h5NOoZmHUNhTMewvfWUUvzizks45p5Fn/ZgHpua46mF8wUy6mMfsDwMWcRU01SG6m7gpHsFmHaHGd1RXxouLO5cPWfYAvQASTERXFNYwHsls5neXorW+qzHCBFqMjNRhE7lwPhpM8vmIKPzaUW5mUKBu5VVK0Y+H/nIoix26gVk6QaaqitGfIwQoSSFWoRO5Xa6jTiadJLdSUZVaSbTj5src91DfdNnumJOGkV6AYbSNO+TfmoRflKoRcjoyu2O758+aqYRr/q4/eqLRn1MTISbiGmX4NUu+o+8Fb5wQgyQQi1Co/UEqu0E1a7pdicZVaflocZKYF5kK3PnzDnnYz+8dCb79ExiW0plXQwRdlKoRWhUvQtAheXc8dO1ln/J1esWZmIY534rXLMgg13WfPKooaX+VDjiCTFECrUIjcp36FNR1FopdicZVZ0VTwQmqxfPPO9jsxOjqY5fikeZNBW/HoZ0QnxACrUICV25nSpy0Mq5L7E6K45MVyczCgLrnslecjWmVjQWywVFEV7OfReJyauzAdV4hJOGcxfn6tVu2nQ0+dH9xMTEBPScKwtnc1BPJ6GzLMTphBhOCrUIvoH+6UpybQ4yujrLPxTv4mmBDx1cnJvAPmMBszhBe0tjqKIJcRYp1CL4qt7FpzycMNPsTjKqOisOA4urCmcE/BylFO0ZlxClvNS+L/3UInykUIug05XvcJJsLOWyO8qo6s040l3dzJ4xtuGD0XOuAKC26JVQxBJiRFKoRXD1tkNtCSdUvt1JRuXVBo06hhxPD4mJiWN67tyZszhmZZPYeTRE6YQ4mxRqEVwndqG0xQkj7/yPtUmjFYvGYElWzJgXWCrMS2S3NZ+ZnKCnuytECYUYTgq1CK7Kd7AwOO4b327L4eC/kKj5cOHYZ03GRbqpiFlMvOqhtuStoGcTYiRSqEVQ6ap3qTWy8CqP3VFGVWfFkWL0smB2wbier6evBqB93/ogphJidFKoRfB4e6G6yNHD8iwN9VYcWa5OMjLGd9afP3MBp3QK1sndQU4nxMikUIvgqS5Cmf2OnujSrGPw4WJ+asS4NwBYlp/Mbms+BZykv68vyAmFOJsUahE8VdsBKHdw/3Sz5Z+FeNnc8S8WNS8rnvf0PJJVB3WHdgYrmhCjkkItgqdyOw1GBr0q2u4ko2qyYvBgsnT2+LtnItwGjSkXA9BxYEOwogkxKinUIjhMH7pqJ5U6x+4k59Sso0lx9ZCTPbHlVzNmLaNFx2FVbAtSMiFGJ4VaBEftPpS3i2q3czcK0Nrf9ZHm6h11261ALZ2Wwm5rHpn9xzFNB+/cK6YEKdQiOAYWYjrmdW7/dIeOxIeLgiT3hHcSX5KXyE5rPumqhbqy4iAlFGJkUqhFcFRup81IpkPF251kVM3afyFx6bTUCbdVkBrLftdCANpLZDy1CC0p1GLitEZXvkuFw/unm6xoFJqLZ088p2EoXDlL6dKR6OPvBCGdEKOTQi0mrvEIqqeJUw7eyBb8/dPJRi/TcoOzj+P83FSKrLkkdx7CsqygtCnESKRQi4mr9J9RlptZNgc5t2YrhmTVTWrqxLs+AOZnx7PTWkAWjTRUHQlKm0KMRAq1mLjj79BlxNOkx7ZkaDj1ajfdRJAXx3l3HA/U/Kx4dlnzAWiVdT9ECEmhFhOjNfr4Nip0LkxwJEUoNVv+STgLsxOC1uacjHhKmEm/dmOWbwlau0KcSQq1mJiWClRnLdUO759uGpg6vmJ2cPqnAaIjXGQkxVOqC0hqK0VrHbS2hTidFGoxMQMjHo6ZwSuAodCsY4hV/cydHtyRKQtykthpzSdL19JcUxnUtoUYJIVaTEzldnqMWBp0st1JzqnZiiHF6CEzMzOo7S7MTmSrtRgDTVOx9FOL0JBCLSZEH9/GcZ3j6P5pn1a06SiyIr1EREQEte15WfG8Z83F1Apv2VtBbVuIQVKoxfi1nkC1VTm+f7pVR6NRzE4L/qp+C7Lj6SaKCnJJaNkv/dQiJKRQi/GrHFh/2nL2jMTBER8XzQj+OiT5yTFEuRXFzCPbOkVrQ03QjyGEFGoxfpXb6FPR1Fopdic5p1YrGjcWhTODf8HTMBRzMuJ4x1yAG1P6qUVISKEW46aPb6eSHLRy9suoRUeTbPSSFeQLiYMW5ibxllmIhaL/yJshOYa4sDn7HSacq6MW1VzGCQfvjzioxYomyeghMTE0MycXZMXTrOOpJov4JlnyVATfeQu1UipfKbVZKVWqlDqglHooHMGEww2s7+HkHccBerWLHiLIiWXCa1CPZl6Wf7ZjqZ7h76durA3JccSFK5Azah/wTa31QuAy4CtKqYWhjSUcr2Ir/SqSk2aa3UnOqXXgQuLczNCtkz0/y9/2e8zHjUnj+6+G7FjiwnTeQq21rtFavzfwdQdwEBx+GiVCTh/fShW5k6J/GmB5CEZ8DEqOjSAt1s1OayEWir4jm0J2LHFhGtO7TClVACwHdo7wvfuVUnuUUnsaGhqCk044U3sNqqmMSsPZ46fBf0YdoUzmTQvtEqzzsxKoMpOoIYME6acWQRZwoVZKxQHPAw9rrdvP/L7W+jGt9Qqt9Yr09PRgZhROc3wrABXk2xzk/Fp0NMkq+FPHz7QgJ5E2K4rj5JFtVUs/tQiqgAq1UsqDv0g/rbX+S2gjCcer2EKfiuaU6ezx01r7R3ykuHqJjY0N6bHmZMZjYlCqZ0o/tQi6QEZ9KOBx4KDW+pehjyScTh/fyvFJ0D/dg4d+3OTFu0I24mPQ4MXKvXqO9FOLoAvknbYauAe4Wim1d+B2Q4hzCadqPYFqOc4JV4HdSc6rZWDEx/wgbhYwmjkZcQDUqVTppxZB5z7fA7TW2wDnLo0mwmugf/qY5fyBP4OFevnM0O/lGBvpJjshgpaeaCpdeVxq7aW1oYakdGev0y0mB2d/dhXOU7GVXiPG8et7gH/VvBjlZXZ+eDbdnZ+VQLMvknKm48akoeilsBxXTH1SqEXgtEZXbKFC5zl6/elBLVY0SaqHcI1CmpedSJv2j/wwMfAdfiMsxxVTnxRqEbiW46j2k5Oif1praNVRpHn6iY4O/jrUI5mbGYeFoknHcYJskltLZH1qERRSqEXgBvqny0xnrz8N0Kkj8OEiP/G8l2GCZnDkR4uOpoJpZOpaGk8cCdvxxdQlhVoE7thmOo0Ex++PCB9MHV+UG76sszPiMBS0qzjKmY4CmoteDNvxxdQlhVoExrLQ5W9TrvMnSf90DADLZoR2RuLpojwuchIjaTEjOUUmfXjQ5W+F7fhi6pJCLQJTuw/V00Sle6bdSQLSoqOJN/opyAvv8LgF2Yk0m5FYykUleaR3lGJZVlgziKlHCrUITPlmAA57nd8/Df4RH8mqO2wjPgbNy0qgXUdhakUF00ilhdrDe8KaQUw9UqhFYI5tpsHIoEuFds2MYDC1ok1Hkh7hw+PxhPXYczLj0CjadBTl+He/ad8r46nFxEihFufn7UFX7eCYdv5qeQBtOgqNwYyUyLAfe97AJgKtVjT1pNFJDK7KbWHPIaYWKdTi/KreRZl9HHfNsDtJQAanji+ZFv7dZ2akxfpHfhixoBQV5JPVe4T+vr6wZxFThxRqcX7HNmMqN+W+ybFuRYuOxsBiaRjW+DhTpNvF9OQoWswoAI5RQDxd1BRvDHsWMXVIoRbnpcs3c5IcvCq8/b3j5d91vI/c7PAXaoD5OYk0mf5ul2P4d8HpLpZ+ajF+UqjFuXU2oGpLqDAK7E4SsBbtX+MjJcWehaPmZyXSoSPxaoNOFUct6cTV7rAli5gapFCLc6t4G4BjAyMYnK5fu+jSkeTEaAzDnpe3/4KionVgdmQZBeSYVbTUnbAlj5j8pFCLcyvbSK8RQ7UZ/gtz49Fi+fuGZ6fH2JZhQfbAyA/tz1BGAS4smnbLLnZifKRQi9FZFvroG5Tp6Y7fdmvQ4Bofy2Zk2JYhPzmGKLeiFf+Y8xPk0IcH8/DrtmUSk9vkePcJe9S8j+pupNw9x+4kAWuxYvBgsqjAvhEqhqGYnR47dEHRUi7/anod+zF9PttyiclLCrUY3ZENaBSHvM7fdmtQi44m2eglMzN8izGNZFFeMs1WNIPLUR+jgCTaqDkgk1/E2EmhFqPSRzdwysihR4Vn4f2J0to/IzDF1UN8fLytWeZnxdOn3fTgH9JYNjBMr/P9F+yMJSYpKdRiZJ31qFPvcVRNjtXyAHrw0IebvDgDZfNSrPOz/DufD86SbFVJNJJM1Ek5oxZjJ4VajKzMP5PuKJNj2jh8UBTnZyfanMR/Rg3QdtoiVscoINd3nI7mertiiUlKCrUY2dENdBnxnDJT7U4SsOaBQr18pr390wDJsRGkxriHNjAAOMJMPPio3/GMjcnEZCSFWpzN9KHL3uSILpgUu7kMatYxxKp+5kxzxpokC7IThkZ+ABwnj14iMEtftjGVmIykUIuznTOgRA0AABqgSURBVNiJ6mun3DXb7iRj0mTFkGp0k5VlzxofZ1qYm0SrFYU1MPLDUi7KmEFO5z68fb32hhOTihRqcbajr2Ph4ohvcuzmAuDVBm06iuwoHxEREXbHAWBeZjwmBu06aui+w8wijm6qd//NxmRispFCLYbTGn3wZSpVPv0q/Avvj5e/f1oxL8O+qeNnmj8wlXzwIifAUQowMejdK9PJReCkUIvhGg6hmo9x2DXP7iRj0jSwrsYls53RPw0wOyMOQ0GbETd0X5+KopI80pp2yqa3ImBSqMVwB/+GRlFqFtidZEyarViilI9FM5zTXRPpdjE9JYpmM2rY/YeZSZpuoq50u03JxGQjhVoMow/+jWojlw6cv4nt6fwXErvIznbOGTXAwpwkms4q1LMAaNu91o5IYhKSQi0+0HIcVbuPw2ryLMIE/l3HW3QUmRFeoqOdNd19UW4inTqSPu0auq9NJVJLGnEnt9iYTEwmUqjFBw76x/eW6lk2BxmbFh2NxmBOqvMufhbm+mdJNlnDL3IeZhY5ZhXNJ4/YEUtMMlKoxQcOvUy9kUmzTrA7yZgMFsEVs5wxfvp0i3P8hbqF4YtEHWQuBpqmrU/akEpMNlKohV9HHbpqB4cmWbcH+At1BCZLZznnQuKg5NgIshMiaNTDz6jrSPMv0nTsVZuSiclECrXwO/wKCs1BJtdsRBi4kOjqJifHeYUaYGl+Mg2+M/rOleIA88j1Hae56pA9wcSkIYVa+JW+SKuRSq2ZbHeSMbE0tOgY0t19xMXFnf8JNliSn0THGRcUAQ4Mdn9sedymZGKykEItoKMWXbGFEuZOqkWYANp0FCYGs5I9dkcZ1WgXFBtUGvWkEl0heymKc5NCLWD/8yhtUaIW2p1kzJos/3jvi2ak25xkdIOFupmzd505wDxyzUoay/eFO5aYRM5bqJVSTyil6pVS+8MRSISf3vcstUY2DZb9C+6PVaMVgxuLi+bk2R1lVEkxEeQkRAxNcz/dfuahgOZtT4Q/mJg0AjmjfhK4PsQ5hF0aj6Jq9lKiFtidZFzqrDgyXF1My3duoQZYMtIFRaBZJVNDOnGVb6AHd8IV4gznLdRa6y1AcxiyCDvsexaNosQ6e1ieVxuU+VJ435tNoxWD0+pIvzZo0THkRvTavpnt+Yx2QRH83R855kkaj+6xIZmYDNzBakgpdT9wP8C0adOC1awIJa3RJc9SaUynQ3+wtkeTFc0+bzYnrCTMgd/le325xKteZrmaWOquxVD2V+0GKw6NYmmuM0d7nO70C4o5ro5h3ythPtewjda3HyV97iV2xBMOF7SLiVrrx7TWK7TWK9LTnXthR5ymugjVcpz9xgcXEeutWF7rm0eNTmCep4kH53Tz7B0FfGlpNKkRFnt9ubzVPxNT2z86pM6KQ6FZs3i63VHO61wXFNtVAseYTkb1RkyfN9zRxCQQtDNqMQntexZTudnvmwEK6sxY3uifS7Ty8bVCzb23fhaPxz/s7dJli/gO8OvXivnV27C5fxZrIo7htvHMus6KI9XoYf5s5++UPnRBsWfkjQ32sojP8Con3llL/lX3hDmdcDoZnneh8vag9z3DYWbRpyKpN2PZ0D+XGMPH15cp7rvz00NF+nQPfWwp//iRAk5YSWzqn23bmbWpFQ1WLFnuLlJTJ8dO6aNdUAQ4xGx6iMS7+8nwhhKTQiDD8/4MvAvMU0qdVErdF/pYIuRKX0L1tvK+azn92sXm/lnEGD6+sczg87fdjDrHxJcHrlnE966fSbWVyF6fPes/N+sYTFwsTI88Z1YnWTYtmQ4dSY8++4OsqdyUMJ/8zr10NlbbkE44WSCjPu7UWmdrrT1a6zyttcx3nQqK/odWI5UyXza7vHn04OHWvE4+d+tNARW+L65ZwMcXplLiy6bBCv8mA3Wm/wLimkX5YT/2eF1SkAL4u2xG8j6L8eCjbuPvwhlLTALS9XEhqj8IVe9SpAqpthI5aqazNKKer302sCI96Ke3XUxylMFW7wx8Ye4CqbPiSDD6WDK3IKzHnYjC3EQiXYoGRp5YVEsGtaQRe/QFGVMthpFCfSEqehJTudlpLmKbt4Ako5e//3ghiYljm5mYEOXh13etoM2KosgXvgknWvsLdabR6bitt84lwm2wLD+RGt8on0CUYi+LyTJPUb/vzfCGE44mhfpC4+1BF/+ZQ8zmbe8cerSHT+d2sfLSFeNq7oq5Gdy2PJNSXyZ1Zni6QNp0FH14mJOkcLnOnkDiZCtnp9NkRY848QWgmAV4cdP91iNhTiacTAr1hebAC6jeNraoSzliprPQ08jXPnvjhC7I/eBTy0iOUuz2TQvL7MX6gT7e1XOdt6PL+Vw6IwVQQz/DmXpVNMUsJL9lO10NVeENJxxLCvWFRGvY/TgtRhp/6v0QHizuX5VHUlLShJqNiXDzT59YTIMVS7mZEqSwo6ux4olWXj600Pnjp8+0PD8Zt8Go/dQAO1mOG5O6V34WxmTCyaRQX0iq3oXqPbyhP8QJK5nlUQ1cf/UVQWn6lovymZsWSZEvD68O3cvK0oqTZiL5rvZJuVRBdISLRdnxo/dTA40qlTKmk378RXx93WFMJ5xKCvWFZNsj9Bqx/Lr3Y8Sqfv7v9YVERgZn527DUPzrLcvp0hEcMEPXJVFjxdOPm4sy3URERITsOKG0anY6jVbMOX+h7eAi4umk+g0ZqiekUF846krh6Ou8Zl1GtU5hZVwTKy+5OKiHuHRGKtfOTWafN4suHZodV6rMJNyY3HzZ/JC0Hw6XzkzFQp1z/PkxCvyb377/P2jLCmM64URSqC8U23+LT0Xwy76bSDZ6ePjmVSEZMfH9m5eBUrwfguF6WvsLdZ6rg8KF84LefrhcPD0Zg3P3U6MUO7mITLOaut1/DVs24UxSqC8EbSfRJc+ySa/ghE5jTUoHixaGZqOA/JQYPnfZNI76Umi2Rl7XYrwadSzdRLA42SQ2NvyzIYMlIcrD3IwYasxzL8+6l4V0EoO1+adhSiacSgr1hWDHo2jL4hd9nyLN6OYrN18R0vUxHv7IAmI9iiJfcC/2VZlJKDQ3XTIrqO3aYdWcDOrNmHMuauVTHt7hEnJ6j9Cw58UwphNOI4V6quuoQ+95nB0s4bDO5cNpXcyaFdpClxjj4aFr53HSjKfaTAhau5VmEtmuTlYsmXyb8J7p8jlpmBjUWOfemWYPS+gkBt8bPwpTMuFEUqinurf/De3r51/77yDD1cWDN18ZltXm7l09g8xYF3t8+VhBmATTZkXSpqOZF9tLcnLyxBu02apZaUS5FdWce4lWn/KwjUvJ7iujftfzYUonnEYK9VTWWIYuepI3rUsptaZxTUYvM2fODMuhI90u/vmThTRb0Rw10ybcXqXpL843LHX2JraBivK4uHJOGse9CeedzVlEIR3EYr75Yxy3caUICynUU9mmf8Gn3PxL/+1kurp44KYrw3r4TyzNYUl2DEW+PHpHWIM5UJaGo2YaGa4uLr9oURAT2utjS3Lo1h4a9bkvjPqUh61cSnbfMep3PBOmdMJJpFBPVSeLoPSvvGheQZXO4KNZvRQUFIQ1glKKn9+xAq92TWh1vUormXYdxYcSO8nKmnzre4zm6nmZGApO6PNPu3+PQppJxPPmd9G+/jCkE04ihXoq0ho2fp8eI46f9N1KnquDBz59jS1R5mbGc+9l+RzxpQ0t9j8WWsM+bzZJRi9f/sTKSbObSyASYzysmJZIpe/8y8uays3rrCHZV0/Niz8IfTjhKFKop6J9z8LxrfzBdx3NxPPJAk1ubq5tcb59wyLSog3e9U3HGuMGA9VWAs06hkviWlm4YPLORhzNDUtyabWiaLPOP5X/CDM5wgxSS/6bviZZWe9CIoV6qulqQr/+j1Qbefyi/yYK3K184SZ7zqYHxUS4+fEty2ixose0x6LWUOzLJk718+ANK6bU2fSgjyzyd+WcsAJYdVAp1rMGFyaNf3owxMmEk0ihnmpe/yd0Tys/9d5FH25umRNBRkaG3am4blEWNyxIodiXQ6UZ2LKqdVYc9VY8F8U0s2xJYYgT2iM3KZp5GdFUWYH9m7SoZLazgtymbbQWvxLidMIppFBPJWVvwr61bGYlr3mXMMfdwt03Xmt3KsB/YfGXd13KnBQPW/pnnHd6uVcb7PbmE628PHDdEgxj6r5Ub1iSR50ZE/BCVtu4lBYSMV76ClZPa4jTCSeYuq/+C01vG/rlr9PiSudnfbeg0NyzLMlRk0OiPC7++OUrifEYbOqfTe8o21GZWrGpfxZNOoarExr40MXLw5w0vG5enotCcdQM7JOPV3n4Cx8jzmyl4akvhDidcAIp1FOB1vDiV6H1BP9j3chRK5PlkQ3c+nFnnE2fLjMhiv+5byXdOoIX+xad1Q1iaXirfyanrEQ+HHuKHz94K273+MdgTwbTU2NZPTOJQ760c679cbqTKoctfIjM2k00b/nvECcUdpNCPRXseBQOvsRG40r+t3cl8aqPhz5WSHR0cFevC5aLC1L54xcvJt6j2NQ/m019syj1pbOzP59X++ZTZSWzOvoU//5/bpnwNmGTxZfWzKVHezhuBv4JaAuXcYJsYjZ9h/6G8hCmE3aTQj3ZVe1Ev/Fdjhhz+b33Olp1DGuSWlg1zl3Fw2Xl3Gy2fu+T3LsskZNWIju90zlipmMqF5dHn+KXX76R1NRzr4MxlVwxO428BA+HrMyAn6OVwQt8DIVF1+M3o/s6Q5hQ2Glqf6ac6tqq0evupV0l8GfzBt735pHj6uCbt18zKS6+Rbhd/PCOy3ng2jYOHTnKtIwkUlJSSExMDMmmBk5mGIr7rprDD/9WSoMVQ7oR2F6JLSqJ5/XHuaP3RVp+fwspD74KxoX1b3chcP67WYysqwn9h0/h62phnbqRDf3zMDG4a64R9qniE5WTlsjVq1Ywe/ZsUlJSLrgiPegzF+cR5YLD1timyR9VM1nPh0lp2EHz2v8TonTCTlKoJ6O+Dnj6Fqymcp4xPsWbffOothJZFVPHfbd90u50YpziozzccnEex7xJdI9xz8ndahk7uIiUI2tpW/+TECUUdpFCPdn0d8Gf70TXFPO860be8xawy5tHnqud79+1xrEXEEVgvnTlbJRS7PGOfRGrDVzJAeaQuOPf6Hz1B7Ik6hQihXoy6ayHJz+BPv4OLxk3sN9XwJb+GbjRfHl5HLNnz7Y7oZiggrRY/m71dI6ZqdSY59795UxaGTzPxylmAXG7fkX7C9+UYj1FSKGeLBqOoH9/LWbtfp5z38z75ly2eQto1HFck1jP7Tddb3dCESRf/+hCMmP9i1gFOq56kFYGf+V69rCEhH2P0/Gnz4O3NzRBRdhIoZ4MDr2Kfvwj9HY087/GHZSaBRT5cik3U7k0qpYffelTeDxj69MUzhXlcfGzWy+izYqixDeO9beV4hWuYQsfIv7oX+n67Wp0c0Xwg4qwkULtZH0d/hmHa++kwRvF4+pOqsw0DvrSKfFls9DTyL9/4VrS0ia+1ZVwlg/Pz+TqOUns82XTakWNvQGl2KxW8yduxmg/gfc/VuItXiddIZOUFGon0hqOvI7+z8vR7/+Rd4zL+C/zdhp9cRzwZbDDO43p7jb+7fYVTJ8+3e60IkR+euvFxHgMNvTPpcOKGFcbR9VMHuNuGnUinhf+jp7HPwHNMotxspFC7TSn9sJTn4Q/3UZbewd/cN/FRr0KL262e6ezyzuNAncb/3rDLAoXT539A8XZMhOiWPvl1ZjKzYb+eWMesjeoVSXye30Hr3MVxsldmL+9hP7134WupiAnFqGidAg+Cq1YsULv2bMn6O1OWVrDsU2w87/g6Ov0GrFsUSvZaS7CUi46LQ9bvTOotRJYHlHHz+6+nHlz59qdWoTJrvIGPvvfO4hTfXw04jDRyjfutuJ0Jx9VW1msD2IakeiL7sV9xUOQODV2d5/MlFJFWusR136QQm2n5nIofRG990+oxiP0GHHsZgnv6ovoJQKfVuz3ZbHPlwUoro6r4cf33+yIjQBEeG0uPcXf/W8Rbiwu81Qxw9XMRDa8SdNNXK6KKNQHUIB32uVErPgczP84RJx7V3QRGlKonaK/G07ugoqt6KOvo2pLADhl5FJkLKfYNxNTuenVbo6ZKZT6MunUkcxwt/K5JXHc8cnrZELLBWxfZQNfeWo7J7rdTDNaWOE5SaLRN6E2E3U7F6sDFOoDJNHuP8uevgr3vI/CzA9D2lyYBOvGTAUTLtRKqeuBXwMu4Pda65+d6/EXfKG2LOiogZYKqD8ItSXo2hKo3Y+y+rEwOGXkcEjN5YA1i1bi6bQiqLXiOGElUWUmYWGQbnTxkcxuvnb79WRljWOYlphyTEvzk+d38GRRIyYGmUYHc12N5LtaiVTm+BvWmumcZJFRzkyrnFRa/MfzxEHOclz5KyB9vr9wp82ByLFNxhHnN6FCrZRyAUeAjwAngd3AnVrr0tGeM+kLtdZgesHs99+8PeDrBW+3fwp3X4f/1tsK3c3Q0+KfNdhRg26vgbYTKPODM51eFUM16Zwgh0PMosTMp1kn0K4jabeiaNbRdGn/LtRRysdsdwsfm5/EzWsuIScnZ0pu6iomprK+hd/+bRdvHOsa2sE8UfWQbnSRpHqJVf3Eqn6ilJcIZRKBiUsF/uk5SbcxQ50g32ggy6wmg0ZcWEPfNyMSIDEPIykfFZcOsekQkwZRiQO3BIiIA08MRMSAOxrcEeCOAleErPA3gnMV6kCWOb0UKNNalw80tha4CRi1UI9X2b9chNs6+6OcYvQX2PASpoceq4aep4e+VlgowMBCaY2BxsDChTn0p3vgNhY9OoJGEmnQSdTrZE7p2RzXWVTqDMqsXKpJOyspgAuLBKOPDKObGfGdrJyVxuWF85g1cyYREeMbjiUuDNMzkvn5fR/FNE3+9u4BNuyt4FCDycm+RMoYeVy9wsKFxoWFgUYp/H+e9p7x/3n2+82Nj3xVzwxVy3RVR7bZRFZPM5l1paSoDpLpwDOGM3oLhQ8X5tDNGEhlYGKg1eC7VQ0k/CClPi2hHva++uDr038CPcJ7L3Bje26XO4FF39k+geONLJBCnQucOO3vJ4EPnfkgpdT9wP0A06ZNG1eY5pgZdLWOPGQo0H/ss/8T/X+3hn1tYA2VbmPoReIbeMGYGHhxD936tId+PPThoYcouomkW0fRQQwdxODFg0LjUgqXy8BQ4FbgNiDfUMx2txPtMYiLcJEaH0lmUizZSbHMzE4hPS3tglx/WQSHy+Xi5suXcPPlSwDo6+ujvKqaY7UtVDa0U9PcSXtPPx19Pnq9Fv2mxmtpfJb/g6OFf/sz8P99eIEbrp1UikmlmA+GhWoLtGUBmjh6iKOHWOX/M5o+olQ/UfQRgW/g5vWfDCn/CZG/TFu4Bt6RH5Rr67QSzdCfZ56MDRrtZO5cJ3nnM54Pst1mFLN6e4mKGsckpXMI2sYBWuvHgMfA3/UxnjYu/ebzwYojxAUpMjKSBXNmsmCO3UlEMAVyObcayD/t73kD9wkhhAiDQAr1bmCOUmqGUioCuAN4KbSxhBBCDDpv14fW2qeU+irwOv7heU9orQ+EPJkQQgggwD5qrfWrwKshziKEEGIEMuVICCEcTgq1EEI4nBRqIYRwOCnUQgjhcCFZPU8p1QBUBr3h8UkDGu0OcR5Oz+j0fOD8jE7PB5IxGCaSb7rWOn2kb4SkUDuJUmrPaAudOIXTMzo9Hzg/o9PzgWQMhlDlk64PIYRwOCnUQgjhcBdCoX7M7gABcHpGp+cD52d0ej6QjMEQknxTvo9aCCEmuwvhjFoIISY1KdRCCOFwF1ShVkp9UymllVIj71VkE6XUvyulDiml9imlXlBKJdmdaZBS6nql1GGlVJlS6v/Zned0Sql8pdRmpVSpUuqAUuohuzONRinlUkq9r5R62e4sI1FKJSmlnht4HR5USq20O9PplFJfH/g/3q+U+rNSKrhbqIwv0xNKqXql1P7T7ktRSr2hlDo68GdyMI51wRRqpVQ+cB1QZXeWEbwBLNZaL8G/kfA/2pwHGNrY+HfAx4CFwJ1KqYX2phrGB3xTa70QuAz4isPyne4h4KDdIc7h18B6rfV8YCkOyqqUygW+BqzQWi/Gv9zyHfamAuBJ4Poz7vt/wJta6znAmwN/n7ALplADvwL+nrO3grOd1nqD1to38Ncd+HfRcYKhjY211v3A4MbGjqC1rtFavzfwdQf+4pJrb6qzKaXygI8Dv7c7y0iUUonAlcDjAFrrfq11q72pzuIGopVSbiAGOGVzHrTWW4DmM+6+CXhq4OungJuDcawLolArpW4CqrXWxXZnCcAXgdfsDjFgpI2NHVcIAZRSBcByYKe9SUb0CP6TBMvuIKOYATQA/zPQPfN7pVSs3aEGaa2rgZ/j/zRcA7RprTfYm2pUmVrrmoGva4HMYDQ6ZQq1UmrjQP/VmbebgH8CvufgfIOP+Q7+j/NP25d08lFKxQHPAw9rrdvtznM6pdQngHqtdZHdWc7BDVwEPKq1Xg50EaSP7MEw0M97E/5fKDlArFLqbntTnZ/2j30Oyif4oO1Cbjet9bUj3a+UKsT/H1ys/Pu/5wHvKaUu1VrX2p1vkFLq88AngGu0cwa3O35jY6WUB3+Rflpr/Re784xgNXCjUuoGIApIUEr9UWvtpEJzEjiptR78NPIcDirUwLVAhda6AUAp9RdgFfBHW1ONrE4pla21rlFKZQP1wWh0ypxRj0ZrXaK1ztBaF2itC/C/KC8KZ5E+H6XU9fg/Gt+ote62O89pHL2xsfL/5n0cOKi1/qXdeUaitf5HrXXewGvvDmCTw4o0A++FE0qpeQN3XQOU2hjpTFXAZUqpmIH/82tw0MXOM7wE3Dvw9b3Ai8FodMqcUU9y/wFEAm8MnPXv0Fp/2d5Ik2Jj49XAPUCJUmrvwH3/NLDHpxib/ws8PfALuRz4gs15hmitdyqlngPew981+D4OmEqulPozsAZIU0qdBL4P/Ax4Vil1H/6lnm8LyrGc8ylbCCHESKZ814cQQkx2UqiFEMLhpFALIYTDSaEWQgiHk0IthBAOJ4VaCCEcTgq1EEI43P8P65RFUaIQ0OYAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "def func(x):\n", " return np.exp(-x ** 2) + 3 * np.exp(-(x - 3) ** 2)\n", "x = np.linspace(-5, 10, 100)\n", "rv = Gaussian(mu=np.array([2.]), var=np.array([2.]))\n", "plt.plot(x, func(x), label=r\"$\\tilde{p}(z)$\")\n", "plt.plot(x, 15 * rv.pdf(x), label=r\"$kq(z)$\")\n", "plt.fill_between(x, func(x), 15 * rv.pdf(x), color=\"gray\")\n", "plt.legend(fontsize=15)\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAgAElEQVR4nO3dd3xc1bXo8d+aot6bLUuy5F7AXTGmmHrJM9UQSDAQhwCJL4Qk5NKSkFwI5IWQvLybF178SGjB9+KLwbSQAA7NhGqwZFuuuFu9WZLV68x+f8yMLMkqI2lGM5LW9/PRx6OZM+csaeQ1e9ZZex8xxqCUUmr0swQ6AKWUUr6hCV0ppcYITehKKTVGaEJXSqkxQhO6UkqNEbZAHTgpKclkZWUF6vBKKTUq5ebmHjfGJPf2WMASelZWFjk5OYE6vFJKjUoikt/XY1pyUUqpMUITulJKjRGa0JVSaozQhK6UUmOEJnSllBojBkzoIhImIl+ISJ6I7BGRh3rZJlREXhCRQyLyuYhk+SNYpZRSffOmbbEVuNAY0yAiduBjEXnLGLOlyza3AjXGmOkisgr4DXCdH+JVSvWirq6OiooK2tvbAx2KGga73U5KSgoxMTFDev6ACd241tdt8BzP/dVzzd2VwC/ct18C/igiYnRtXjWGldY288mhKr62KA2LRQIWR11dHeXl5aSlpREeHo5I4GJRQ2eMobm5meLiYoAhJXWvJhaJiBXIBaYDa40xn/fYJA0odAfVISK1QCJwvMd+1gBrACZPnjzoYJUKFjuLTnDruhwq61uJCLFy6bzUgMVSUVFBWloaERERAYtBDZ+IEBERQVpaGiUlJUNK6F6dFDXGOIwxC4F0YKmInD7oI7n284QxJtsYk52c3OvMVaWC3j/2lPGNP39GiNVCZmIEf3j3IE5n4D6Mtre3Ex4eHrDjK98KDw8fculsUF0uxpgTwGZgRY+HioEMABGxAbFA1ZAiUiqIvbevnNuey2XWxBheveMs7rp4JvvL63lzd2lA49Iyy9gxnNfSmy6XZBGJc98OBy4Gvuyx2evATe7b1wLva/1cjUUv5RYxITqMDd9dRkp0GJfPn8T0lCj+8O5BHAEcpSsF3o3QU4HNIrIT2Aq8Y4z5u4g8LCJXurd5GkgUkUPAXcBP/BOuUoHT4XDy8aHjnDczmfAQKwBWi3DnRTM4WNHAG7sCO0pXypsul53Aol7uf6DL7Rbg674NTangsqPwBPUtHZw3q/v5n8vmpfLYewd57L2DXDYvFWsAO17U+KYzRZXy0ocHKrEInD0tqdv9FotwxwXTOVTRwLaCmgBFN3Y5nc5AhzBqaEJXykv/PFDJwow4YiPspzx2zgxXkt+uCX3YvvWtbzFjxgzuvPNOHnjgAW699dZujxtjWLhwIevWrfN6n9///vdP2c9YpAldKS9UN7axs7iW82am9Pp4UlQokxMi2F5wYoQjG1tqamrIyclh165dFBcXs3btWm655ZZu27z44otUV1dzww03eL3fe+65h/Xr13Po0CFfhxxUNKEr5YWPDlZiDJw7M6nPbRZNjtOEPkzx8fHs3buXsLAwXnrpJaqqqli+fHm3bR577DFWr16N3X7qJ6W+ZGVlcc455/D444/7OuSgogldKS98eOA4cRF25qfH9bnNoow4yupaKK1tHsHIxpaEhAT+/d//nRtvvJHExEQSEhL413/9V1pbWwE4dOgQn376Kddee223573wwguISK9fdXV1AFxzzTWsX79+TNfkNaErNQBjDB8erOSc6Un9drAsmhwPoKP0ISooKKCmpobf//73JCQksHHjRn70ox/xxBNPsHbtWgDee+89IiMjWbBgQbfnnnvuuXz22WedX2+88QYJCQlceumlnVPozzrrLMrLy9m1a9eI/2wjJWAXiVZqtNhXWk9lfSvnzux/uYo5qTGE2CxsL6gJ6Nouo1VeXh7gOoH56KOPAnDhhRfy4YcfsmnTJu666y5yc3OZM2cOFkv3sWhqaiqpqa7feV1dHRdddBGnn346Gzdu7NzmtNNOw2q18sUXX5zyhjBWaEJXagAfHqwE4LwBEnqIzcK8tNigGKE/9Lc97C2pC8ix506K4cErThv08/Ly8oiIiODee+/tdv/MmTP54osvACgrKyMpqe/zGI2NjVx66aVYLBb+/ve/d1uwzGazERcXR1lZ2aBjGy205KLUAHYX15IeH86EmLABt12UEceu4lraOsZundZf8vLyOPfcc0lMTOx2f1lZGZMmTQKgpaWF0NDQXp/f0tLClVdeSWNjI5s2bSI6OvqUbUJDQ2lpafF98EFCR+hKDeBAeT2zJ56aHHqzaHI8T318lC/L6vo9gepvQxkhB1peXt4pHS319fW8//77PPSQ60JpCQkJvY6w29vbueaaaygtLeWf//wn8fHxvR7jxIkTJCQk+D74IKEjdKX60dbh5EhlIzMneJvQXUk8GMouo0ljYyOHDx/m8OHD3e5/7LHHALjpJtfaf7NmzeLo0aPdtnE4HFx//fXs37+fd999l76W5q6srKSpqYmZM2f64ScIDprQlerH0eONdDgNs7wcoafGhjEhJlRnjA7Srl27cDqdHD58mIcffpjNmzfz05/+lIceeognnniCuDjXG+XZZ59NQUEBlZWVnc+9/fbbefPNN3nkkUcoKChgy5YtbNmy5ZRulpycHESEs846a0R/tpGkJRel+rG/vB7A6xG6iLAoI57thTpCH4y8vDxiYmJ44403uPnmm/nVr37FrFmz2LhxIytXruzc7vzzzychIYFNmzaxevVqjDFs2LCB5uZmrruu+2WMr7jiCl5//fXO7zdt2sR55513So1+LNGErlQ/DpTVY7UIU5MjvX7OoslxbNpTRlVDK4lRvZ/AU93l5eUxb9485s+fT25ubp/bhYSE8M1vfpMNGzawevXqbhOH+uNwOHj55Zc72yHHKi25KNWP/eX1TEmKJNRm9fo5CzJc5YFdxbX+CmvMycvL87o3/N5772Xz5s0cOHDA6/1v3LiR8PBwVq1aNdQQRwVN6Er140B5PbO8LLd4TE+JAuBwZaM/QhpzjDHs2rWL+fPne7V9eno6zzzzDKWl3l9QxBjD008/jc02tosSY/unU2oYmtscFFQ38bVF6YN6XmJkCDFhNo5UNvgpsrHF27JJV4MdaV9//fWD2n600hG6Un04VNGAMTBzQtSgniciTE2O4oiO0NUI04SuVB86O1y8bFnsampyJEeO6whdjSxN6Er14UB5PSE2C5kJEQNv3MO05CjK61ppaO3wQ2RK9U4TulJ92F9Wz/TkKGzWwf83mZrkanM8qmUXNYI0oSvVhwPl9V7PEO1parKr7q5lFzWSNKEr1Yva5nZKa1u8niHaU2ZiBCLauqhG1oAJXUQyRGSziOwVkT0icmcv25wvIrUissP99YB/wlVqZBx0nxCdNXFwHS4eYXYr6fHh2rqoRpQ3fegdwN3GmG0iEg3kisg7xpi9Pbb7yBhzue9DVGrkDXYNl95M09ZFNcIGHKEbY0qNMdvct+uBfUCavwNTKpAOlNUTGWIlLS58yPuYmhTF0eONOJ3Gh5Ep1bdB1dBFJAtYBHzey8NnikieiLwlIr2uri8ia0QkR0Ryui5/qVSwOVrVxJTkSET6vij0QKYmR9Lc7qCsbuxeIWc8+eCDDxARdu/eHehQ+uR1QheRKOBl4EfGmJ7zdLcBmcaYBcD/BV7rbR/GmCeMMdnGmOy+FqFXKhgUVDWSmeD9Cou98azQeFjr6GqEeJXQRcSOK5mvN8a80vNxY0ydMabBfftNwC4ifV/JVakg5nAaimqamZw4+AlFXU3ztC5qHV2NEG+6XAR4GthnjPmPPraZ6N4OEVnq3m+VLwNVaqSUnGimw2mGNEO0q5ToUCJDrNrp4qU9e/awYsUKEhISiIyMZM6cOaxduxaAN954g4svvpiUlBRiYmJYtmwZb7/9drfn/+IXvyApKYnPP/+c7OxswsPDOeecczh69CgVFRVcddVVREVFMWfOHN5///1uz83KyuKee+7hl7/8JRMnTiQqKoobb7yR2tr+l0B2Op08+uijTJ8+ndDQUGbOnMm6deu6bfPxxx+zfPlyYmJiiImJYeHChWzcuNEHv7FTedPlcjawGtglIjvc990PTAYwxvwJuBa4XUQ6gGZglTFGzwSpUamguglg2CP0zkW6jusI3RtXXHEFc+bM4bnnniM0NJT9+/d3rsJ49OhRrrjiCu655x4sFgtvvfUWl1xyCR9++CFnn3125z6amppYs2YN9913H5GRkfzwhz9k9erVhIaGcskll/C9732P3/72t3z961+nsLCQiIiTr/Hzzz/P9OnTefLJJyktLeW+++7jO9/5Tr/J9wc/+AHr1q3jgQceYPHixbzzzjvccsstJCYmcvnll1NXV8fll1/OypUreeCBBzqXCj5xwk9XtDLGBORryZIlRqlgtH5Lvsn88d9NYXXjsPf1w+e3mbN+/Z4Pourb3r17fbo/h8NpKupajNPp9Ol++1NZWWkAs3PnzgG3dTgcpr293Xz1q181N998c+f9Dz74oAHMBx980Hnf2rVrDWAeeuihzvv27NljAPPmm2923peZmWni4+NNfX19533PPfecEZHO3+/mzZsNYHbt2mWMMebgwYNGRMyzzz7bLb7Vq1eb7OxsY4wxW7duNYCpq6sbzK+j39cUyDF95FWdKapUD/nVjditQmrs0FsWPaYmRVF8opnmNocPIvM/p9Nw/ZNbOPPX77HqiS0j1nKZkJBARkYGt912Gy+88AIVFRXdHi8qKuKmm24iLS0Nm82G3W7n7bffPuWqRSEhISxfvrzz++nTpwNw4YUXnnJfcXFxt+defPHFREWdnEh29dVXY4xh69atvcb83nvvYbFYuPrqq+no6Oj8uuiii9ixYwcOh4Np06YRFRXFDTfcwF//+lf/jczdNKEr1UNhdRMZ8RFYLUNvWfTwdLocHSVll6rGNnLza+hwGnLza6hqbBuR41osFt5++20mTpzILbfcwsSJE1m+fDnbt2/H6XRy5ZVX8umnn/Lwww+zefNmtm7dyiWXXEJLS/eW0OjoaCyWk2ktJCQEgLi4uFPu6/nclJSUbt9HREQQFRXV55WRjh8/jsPhIDY2Frvd3vn17W9/m46ODkpLS4mPj+edd96hvb2db3zjGyQnJ3PZZZdx5MiRof+y+qFXLFKqh/yqpmHXzz26JvS5k2J8sk9/SooKYUlmPLn5NSzJjCcpKmTEjj179mxefvll2tvb+eijj/jxj3/MZZddxgcffMD27dt56623WLFiRef2zc3NPj1+z08FTU1NNDQ0kJqa2uv2CQkJ2Gw2Pvnkk25vIh6eN4hly5axadMmmpubeffdd7nrrru44YYb2LJli0/jB03oSnVjjKGgqonszHif7G+yu1PGc6I12IkIz393GVWNbSRFhQxrYtVQ2e12Lrzwws7E5xkhh4aGdm6Tn5/PJ5984vV1SL3xzjvv0NDQ0Fl2efXVVxERsrOze93+wgsvxOFwUFtby8UXXzzg/sPDw7niiivYvXs3v/71r30Wd1ea0JXqoqapnfrWDiYnDm9SkUd0mJ24CDuFNaMjoQNYLEJydOjAG/rQzp07ueeee7juuuuYOnUqNTU1/OY3v2HBggUsW7aM9PR07r77bn75y19SX1/Pgw8+SFqab1cgCQ8P57LLLuPee++ltLSUe++9l6uvvpq5c+f2uv2sWbO47bbbWLVqFffddx/Z2dm0tLSwZ88eDhw4wFNPPcUbb7zBM888w1VXXcXkyZMpLi7mz3/+c7eavi9pQleqi/wqV617uD3oXWXER1A4SkbogTJx4kQmTJjAr371K0pKSoiLi+OCCy7gN7/5DaGhobzyyivccccdXHvttaSnp/Ozn/2MDz74wKfT8FetWkV0dDS33norDQ0NXHnllTz++OP9Pmft2rXMnDmTJ598kgceeICYmBjmzp3LrbfeCrhOwIoI999/PxUVFSQnJ3P55ZfzyCOP+CzursQEqF08Ozvb5OTkBOTYSvXlrzuKuXPDDt7+t3OHtdJiV99bn8u+0no233O+T/bX0759+5gzZ45f9j1eZGVlce211/K73/0u0KEA/b+mIpJrjOm1DqRdLkp1UVDlnlTk4xF6cU2zrrqo/E4TulJd5Fc3MSEmlDC71Wf7TE+IoM3hpLxeV11U/qU1dKW6KKhqGvYqiz1lxLsmKBVWN/tkspLyvWPHjgU6BJ/QEbpSXeRXN/qsB90jw12+0ROjyt80oSvl1tLuoLyu1acdLkDnVY/82boYqOYG5XvDeS01oSvlVuijVRZ7CrNbmRATSlGNb2c2etjtdp/PmlSB09zcjN1uH9JzNaEr5Zbvhw4XD3/2oqekpFBcXExTU5OO1EcxYwxNTU0UFxefsq6Mt/SkqFJu+e6Em+mjWaJdZSRE8MXRap/vFyAmxrVGTElJCe3t7X45hhoZdrudCRMmdL6mg6UJXSm3gqpGokNtxEcM7eNufzLiw/nrjmbaHU7sVt9/MPZcDUeNb1pyUcqtoLqJjIQIvyxIlZ4QgdO4Lm+nlL9oQlfKraimmYwE//SJZ8R7Whc1oSv/0YSuFK4TUsUnmkmP9/0JUaDzjWI0rbqoRh9N6ErhWja3qc1Berx/RuipseHYLKKTi5RfaUJXCihyj5z9NUK3WoRJceEU+qkXXSnQhK4UQOekH3+N0MFVdtERuvInTehKcXKEnubPhB4f0XkcpfxhwIQuIhkisllE9orIHhG5s5dtREQeE5FDIrJTRBb7J1yl/KOoppmYMBsxYb7vQffISIjgeEMbTW0dfjuGGt+8GaF3AHcbY+YCy4A7RKTnRfYuAWa4v9YA/V+3SakgU1Tjvw4XD085x19ruig1YEI3xpQaY7a5b9cD+4CeV2ddCfyncdkCxIlIqs+jVcpPimqa/Fo/B11GV/nfoGroIpIFLAI+7/FQGlDY5fsiTk36iMgaEckRkZzKysrBRaqUnxhjRmSEfnJykSZ05R9eJ3QRiQJeBn5kjKkbysGMMU8YY7KNMdnJyclD2YVSPufvHnSPpKgQQm0WinX6v/ITrxK6iNhxJfP1xphXetmkGMjo8n26+z6lgl7xCLQsAogIafHhmtCV33jT5SLA08A+Y8x/9LHZ68C33N0uy4BaY0ypD+NUym/8Pamoq7S4cD0pqvzGm+VzzwZWA7tEZIf7vvuByQDGmD8BbwKXAoeAJuBm34eqlH94Eqw/e9A90uMj2FtS5vfjqPFpwIRujPkY6Hc9UeO6TModvgpKqZFUVNNEdJiN2HD/9aB7pMeHU9Xo6kWPCNHLESjf0pmiatwbiQ4XD0+dXtdFV/6gCV2Ne66E7v9yC7hq6IAu0qX8QhO6GtdcPej+n1Tk4fkkUKwJXfmBJnQ1rp1oaqexzTFiJZeU6FDsVtFOF+UXmtDVuDYSy+Z2ZXGvi6696MofNKGrce1kD/rIJHTw9KLr9H/le5rQ1bh2coQ+MiUX17HCtYau/EITuhrXik80Ex06Mj3oHmlxEVTUt9LS7hixY6rxQRO6GteKappGZIZoV57yTmlty4geV419mtDVuFZYPXKTijw8byBadlG+pgldjVueHvSMhJEdoXsmF+mJUeVrmtDVuFUzwj3oHqmxYVgtoq2Lyuc0oatxyzNCzhjhGrrNamFiTJhOLlI+pwldjVuF1SPfsuiRpq2Lyg80oatxq3NS0QjX0AHSdXKR8gNN6GrcKqxpIjbcTkzYyPWge6THh1NW10K7wznix1ZjlyZ0NW4V1TSPeIeLR1p8OE4DZdqLrnxIE7oatwqrm0iPG/n6OZys2+uJUeVLmtDVuOTqQQ/gCF170ZUfaEJX49LxhjZaO5wB6XABSI0LQ0RH6Mq3NKGrcanQ04MeoBF6qM3KhGjtRVe+pQldjUuBWDa3p4yE8M43FqV8QRO6GpcKq0f+whY9ZcRHUFStCV35jiZ0NS4V1TSTGBlCRIgtYDGkJ0RQWtdCW4f2oivfGDChi8gzIlIhIrv7ePx8EakVkR3urwd8H6ZSvlVU00R6QuDKLeBaQ8YYKNFFupSPeDNCfxZYMcA2HxljFrq/Hh5+WEr5V1FNc0DLLQAZ7jcUraMrXxkwoRtjPgSqRyAWpUaE02kormkmI4AnRKFLQq/WEbryDV/V0M8UkTwReUtETutrIxFZIyI5IpJTWVnpo0MrNTgV9a20OZwBH6FPjAnDZhEdoSuf8UVC3wZkGmMWAP8XeK2vDY0xTxhjso0x2cnJyT44tFKDd7IHPbAjdKtFmBQX3tlxo9RwDTuhG2PqjDEN7ttvAnYRSRp2ZEr5SeeyuQEeoYOnF11LLso3hp3QRWSiiIj79lL3PquGu1+l/MVTs/aspxJI2ouufGnAJlwReR44H0gSkSLgQcAOYIz5E3AtcLuIdADNwCpjjPFbxEoNU1FNEynRoYTZrYEOhYyECKoa22hs7SAyNHA98WpsGPAvyBhz/QCP/xH4o88iUsrPCqubA14/9/CUfYpPNDNzQnSAo1Gjnc4UVeNOYU1TUNTPoWvropZd1PBpQlfjSmuHg5ITzWQlRgY6FIDOXnhN6MoXNKGrcaWophmngayk4Ci5JEWFEG63aqeL8glN6Gpcya9qBCAzSEboIkJ6vPaiK9/QhK7GlWPHXYkzWEou4Kqj6whd+YImdDWuHKtqJDrMRnyEPdChdEqPD6eougnt9lXDpQldjSvHqprISozEPRcuKGTER1Df2kFtc3ugQ1GjnCZ0Na7kVzWSmRgcJ0Q9PNc11VUX1XBpQlfjRrvDSVFN8LQseniua6qrLqrh0oSuxo3immYcThOEI3TtRVe+oQldjRvH3C2LWUnBNUKPDbcTG26nQBO6GiZN6GrcyK9yJcxgG6EDZCVGdL7hKDVUmtDVuHGsqpHIECvJUaGBDuUUU5IiO3vklRoqTehq3Dh2vJHMIGtZ9MhKiqSktpmWdkegQ1GjmCZ0NW7kVzUFzRouPU1JisQYtI6uhkUTuhoXOhxOCmuagmYNl548rZRHj2sdXQ2dJnQ1LpTWttDuMGQF4QlRONl5owldDYcmdDUuHAuyVRZ7ig23kxgZwjFN6GoYNKGrceFYVfCtsthTVlKkjtDVsGhCV+NC/vFGwuwWUqKDr2XRIysxUnvR1bBoQlfjwrGqJjITIrFYgq9l0WNKUgTlda00tXUEOhQ1SmlCV+NCMK6y2JPnxKhOMFJDpQldjXntDif5VU1MSQ7e+jm4etFBO13U0A2Y0EXkGRGpEJHdfTwuIvKYiBwSkZ0istj3YSo1dPlVjbQ5nMyaEB3oUPrlOWGrdXQ1VN6M0J8FVvTz+CXADPfXGuDx4YellO/sL2sAYGaQJ/TIUBsp0aE6QldDZhtoA2PMhyKS1c8mK4H/NK4LIm4RkTgRSTXGlPooRhXEnE7D23vL+Msnx0iIDOG7505l8eT4QIfVzf7yeiwC01OiAh3KgLKSIrUXXQ3ZgAndC2lAYZfvi9z3nZLQRWQNrlE8kydP9sGhVSBt2l3G/357PwcrGpicEMGXZfW8tbuMM6Yk8NNL57AwIy7QIQJwsLyezMRIwuzWQIcyoCmJkbz3ZXmgw1Cj1IieFDXGPGGMyTbGZCcnJ4/koZWPfXzwOLevz0UE/rBqIe/ffR6f/uRCfn7ZHI5VNXLTM19QfCI4rpG5v7yemROCf3QOMCU5kuMNbdS16AWj1eD5IqEXAxldvk9336fGqMr6Vv7txR1MS47ir3ecw8qFadisFiJDbXxn+VReWHMmDqfhB/+9jXaHM6CxtrQ7OHa8MehPiHp0nhjVsosaAl8k9NeBb7m7XZYBtVo/H7ucTsPdG/Ooa27njzcsIjzk1DJGVlIkv/7aPLYVnOB3/9gfgChPOlLZiNPAjFGS0LV1UQ3HgDV0EXkeOB9IEpEi4EHADmCM+RPwJnApcAhoAm72V7Aq8J786AgfHqjkf151OrMnxvS53RULJrHlSBV//vAIZ0xN4MLZE0YwypMOlNcDMGvi6EjonslPOrlIDYU3XS7XD/C4Ae7wWUQqaBXVNPG7t/ez4rSJ3HjGwCe1//3yueTm1/CzV3fzwb1JhNpG/qTk/vJ67FYJ6kW5ugqzW0mLC+fI8YZAh6JGIZ0pqrz2538eAeCBK+Z6dRm3MLuV+y+dQ2ltCxtzivwdXq8OltczNSmKENvo+VOfMSGK/WX1gQ5DjUKj569cBVRFXQsv5BRyzeJ0JsWFe/285TOSWDQ5jsc/OExbx8ifIN1fXs+MUdLh4jFrYjSHKxsCfkJZjT6a0JVXnvzoCB0OJ7efP21QzxMR7rxoBsUnmnl528iO0htbOyisbh41HS4esydG0+4wemJUDZomdDWg6sY21n9ewJULJg3pij/nzUxmQUYcazcfGtFR56EK95T/UXJC1GPWBNfJ5i+17KIGSRO6GtBfPjlKU5uDOy6YPqTniwg/umgGRTXNvDKCo/T9ng6XUTZCn5YSidUi7C+rC3QoapTRhK761dDawbOfHmPFaROH1ct9/qxk5qfH8vgHh3E6jQ8j7NuBsnpCbRYyEoJ7HfSeQm1WpiZF6olRNWia0FW/Xt9RQn1LB989d+qw9iMi3HrOFI5VNfHRoeM+iq5/ByoamDEhCmsQX6WoL7MmRmvJRQ2aJnTVrw1bC5g1IZrFk4e/0NaK0yeSFBXCc1vyfRDZwA6U1Qf9krl9mT0xmqKaZhpa9XJ0ynua0FWf9pTUsrOollVLM7zqOx9IqM3KN7IzeG9fud8X7qqob6GsroW5qX3PZg1ms9yzcLXsogZDE7rq04YvCgmxWbh6UZrP9nnDGZMxwIYvCny2z95sy68BYHFmcK3N7q3Z7s4cTehqMDShq141tzl4bXsxl81LJS4ixGf7TY+P4KLZKTz/RaFfJxrl5tcQYrNw+qRYvx3Dn9LiwokMsWqnixoUTeiqV3/fWUJ9awervpIx8MaDdOOyTI43tPKPPWU+37dHbn4NC9JjR9WU/64sFmGmnhhVgzQ6/9qV323YWsjU5EiWTknw+b7Pm5FMRkI4/+Wnk6Mt7QKh/W4AABezSURBVA52F9eN2nKLx+yJ0ewvr8e1/p1SA9OErk5xqKKe3PwaVn3FNydDe7JYhBvPyOSLo9UcqvD9CHR3cS1tDidLguzapoM1a0I0J5raqahvDXQoapTQhK5O8VJuMVaLcPWidL8d49ol6ditwvNfFA688SDluk+ILhnlI3TPkgVadlHe0oSuunE4Da9tL+a8mckkR4f67ThJUaF89bSJvLytiJZ2h0/3nZNfw5SkSBKj/Bf/SJjd2bqoJ0aVdzShq24+PXycsroWrlnsv9G5xw1LJ3OiqZ1Nu313ctQYw7b8GhaP8nILQEJkCMnRoTpCV17ThK66eTm3iJgwGxfNSfH7sc6cmkhWYgT//bnvetLzq5qoamwb9eUWj9MnxbCrqDbQYahRQhO66tTQ2sGmPWVcvmASYXb/Xy7OYhFWLZ3MF8d8d3I0x10/z84aGwl9SWY8BysaqG1qD3QoahTQhK46vbmrlJZ254iUWzx8fXI0N7+G6DAb05NH11WK+uJpvdxWWBPgSNRooAlddXplWxFTkiJ9shCXtzwnR1/KLaK5bfgnRz31c8soXGGxNwvS47BapHMpA6X6owldAVBY3cSWI9V8bVGaX3rP+3PTmVnUNrfzyvbhXfyiqqGVAxX1ZI+R+jlAZKiNOanRna2YSvVHE7oC4OVtRYjA1Yt9txCXt76SFc+8tFie/vjosC5+8ebuMoyBf5k7wYfRBd6SyfHsKDxBh140Wg3Aq4QuIitEZL+IHBKRn/Ty+LdFpFJEdri/vuP7UJW/OJ2Gl3KLOGtaIunxI391HxHhO8uncKSykQ8OVAx5P3/bUcKMlKjOlQr9rmS791/DsDgznqY2h7YvqgENmNBFxAqsBS4B5gLXi8jcXjZ9wRiz0P31lI/jVH605WgVRTXNfH2J7xfi8tal81KZGBPG0x8fHdLzS2ub+eJYNVcsmDTiJSN/87RgbivQsovqnzcj9KXAIWPMEWNMG7ABWOnfsNRIeimniOgwGytOnxiwGOxWCzedlcUnh6rYWzL4mZF/zysF4MoFk3wdWsClxYUzISZUT4yqAXmT0NOArj1lRe77erpGRHaKyEsi0utQT0TWiEiOiORUVlYOIVzla/Ut7by5u5QrRqj3vD83LJ1MuN06pFH663klzE+PJSsp0g+RBZaIsCQznlwdoasB+Oqk6N+ALGPMfOAdYF1vGxljnjDGZBtjspOTk310aDUcb+x09Z5/fcnI9Z73JTbCzjey03k9r5jC6iavn3f0eCO7imvH5OjcY/HkeAqrm6moawl0KCqIeZPQi4GuI+50932djDFVxhjPGp9PAUt8E57ytxdzCpmeEsXCjJHrPe/PbedPw2618Mu/7/X6OX/LK0EELpuf6sfIAmux1tGVF7xJ6FuBGSIyRURCgFXA6103EJGu/5OuBPb5LkTlL4crG9hWcIJvZKcHzYnE1Nhw7rhgOm/vLefDAwOX5YwxvJ5XwleyEkiNDR+BCAPjtEkxhNgs2o+u+jVgQjfGdADfB/6BK1G/aIzZIyIPi8iV7s1+KCJ7RCQP+CHwbX8FrHxn/ZYC7Fb/rns+FN9ZPoXMxAh+8bc9A1539N19FRyqaGDlwrFbbgEItVmZnxbLF0erAx2KCmJe1dCNMW8aY2YaY6YZY37lvu8BY8zr7ts/NcacZoxZYIy5wBjzpT+DVsPX2NrBxtxCLp2X6td1z4ci1GblgcvncqSykXWfHutzu6qGVn76yk7mpMYEtOVypJw/K5m8olqto6s+6UzRceq1HcXUt3TwrTOzAh1Kry6aM4ELZiXzh/cO9rp8rDGGn726m7rmDn5/3YJRezHowbh4rqut9N19Q598pca2sf+/QJ3CGMO6T49xelrMiC7ENVgPrzyd2HA71z3xGe9/Wd7tsdd2FLNpTxl3fXVm55V9xrqZE6KYnBDB23t9d0EQNbbYAh2AGnlbjlRzoLyB3147P2hOhvYmIyGCV+84i1ufzeE763K4b8VsosNs7Cmp4287SsjOjOe7y6cGOswRIyJcPHcC//VZPg2tHUSF6n9f1Z2O0Meh//zsGPER9lHRt50SHcaGNcu4YFYKj771JT97dTd/zyth4eQ4fn/dQqxjZJlcb108dwJtDqdXHUBq/NG3+HGm5EQzb+8t57vLpwZ8Zqi3IkNt/Hn1EnYUnmBCTBjp8eFB/cnCn7Iz44mLsPPO3nIunTd2++7V0GhCH2c80+pvPGNygCMZHJvVQnZWgn8PMsxVEUeCzWrhwtkpvLevgnaHE7tVP2Srk/SvYRypqGvhuS35XL0ojYyEkV8mV/nGV+dOoLa5na3HtCdddacJfRx5/J+H6XAafnDh9ECHooZh+YxkQmwW3tlbPvDGalzRhD5OlNW2sP7zAq5dnE5m4thbkXA8iQy1sXx6Em/tKqNdr2KkutAa+jjx+AeHcDoN3x9Po/NRUBMfqhvOmMyt63J4c1cpKxeO/GUDVXDSEfo4UFrbzPNfFPL17HStnY8RF8xKYWpSJE99dBRjhn4dVjW2aEIfB37z1pcYDHdcMI5G52OcxSLccs4UdhXX6oJdqpMm9DFu8/4KXttRwu3nTw/IBaCV/1yzOJ34CDtPDfE6rGrs0Rr6GNbY2sHPX93NtORI7rhgWqDDGd8GW8+ftGjATcJDrHxzWSZ/3HyIo8cbmTIGL7+nBkdH6GPY//rHfkpqm/nNNfMJtY2OWaFqcFafmYndYuEZHaUrNKGPWbn5Naz77Birl2X6f4alCpiU6DCuWjSJF3IKOVheH+hwVIBpQh+DSk40c/tzuUyKDefe/zEr0OEoP7v3f8wmKtTGXS/maV/6OKcJfYypb2nnlme30tzm4OlvZxMdZg90SMrPkqND+Z9Xnc6u4lrWbj4U6HBUAGlCH0PaHU6+t34bhyoaePybS8bNhR8UXDovlZULJ/HH9w/1eoUnNT5oQh8jmto6+OHz2/no4HEeuXoe58xICnRIo4bTCTVN7fhifo4v9zVYD195OolRIfzohe0cb2gd+QBUwGnb4hhQWN3Emv/K5cuyOn5+2Ry+8RUfXTDZD612wcbphPtf28W+0jrmpMbwyFXzsAxxmNPhMPz4lZ0cqmgY9r6GIjbCzu+vW8gtz27l63/6jP+8ZanODB5nNKGPYsYYNu+v4O4X83A4DX/59lc4f1ZKr9s6nYaqxjaSokL8d3GIIb4BjEhsvXA64VhVA3tLanEa2FdaR21LO/ER9lO2q21pJy7cjie8nvc5nfCTV3ayv8zVadLXvvztrGlJrP/OGdz8l61c8/inrLtlKXNStfQ2XmhCH6W2F9Tw2037+exIFTNSonjiW9l9TixxOg3XP7mF3PwalmTG8/x3l2EZxKXbektovnpez9jW33oGNc3tvSZ3XyZ+pxN++uou9pScrDfPnhhDXLi9W+wxYTbuf3U3X5bVMTs1hl9fNQ84OaqfnhLFb6+ZT11LR2cyBwixWogJs51yzKH8HgdrSWYCG287i2898znXPv4pt58/jVvOmUJEiP53H+u8eoVFZAXwB8AKPGWMebTH46HAfwJLgCrgOmPMMd+Gqmqb2nl7bxmv55Xw0cHjJEWF8NCVp3H90smE2Pr+bF/V2EZufg2zzWFaCoTaI1avR44DlSQ8SSoqxErRiWYyEyKwWKTb86anRPHo1fNoaHN0T2Yl26ltaqelIJfZxtCSDz9eu7/XkoXTCde/0TLkN6WeP1N+dSP7Sk8mcwtw71dncqK5exKfnhzJ/vIGAPYU13KsqoG4iBD2ldbhcBr2l9Vz14s7+Pmlcwi1WWjpcLUNNrc7OHa8kdgIOxYRYsPt/Oy13af8Hv2V5GdNjOaV753NL17fw+/ePsC6z/L5wYXTuWL+JOIjQ3x3IBVUZKCV2kTEChwALgaKgK3A9caYvV22+R4w3xhzm4isAq42xlzX336zs7NNTk7OcOMfs5rbHJTUNvNlaT07i0+wq6iWrceqaXcY0uPDuS47g5vPmeLVld+NMax6YgstBbnMSY3h11fP6ywT9JVMOhyGwpomYsJs3LIuB4fTYLUIz968tPPNoGvSxhgcxrVW9/pbz6C+tYNv/+ULHE7X31eY3UJbh5O5k2K7JWqH0/Djl3dysLyeGROiOVjR0OuxapraO/fX8zFvdTgMBdWN/PnDI3xZVk+I1UJzuwOAGcmRWCzCoYoGpqdEs7+fSTozU6IoqG7qTN7emJYcxdHKBpyAVeD/rFrE5IQI7n91F1+Wukf/V88fUs3d6YSq2Ll9fnLJOVbNo299SU5+DRaB7MwELpidwtxJMUxPiWJSbNi4vUbraCQiucaY7N4e82aEvhQ4ZIw54t7ZBmAlsLfLNiuBX7hvvwT8UUTE+GFdz0MVDQG9Uouh9x+p50/q+dGNAacBpzEYY3AYQ4fT0OEwtHU4ae1w0NrhpKGlg7qWduqaOyivb+FEU3vnvuxWYfbEGG4+ewqXzUtlfnrsoP4DigjPf3cZtUes3Wq+P311V7dSgieZdDgMNz79OY2tHUTYrYTarDS1dRBmt3YrI9S2tHeOVD0aWzsorGkiKzGS6SlRnWWIlnZX8ttXUttZW3Y64Wev7eZgRQPTJ0Tz6Nfm8fO/7ukcxXrKHwBx4XbmpMawr7SO2RNjMMZgDF6Parv+TB7NTkfn7YOVjZ23+0vmAAcqGrw7aBeHKxsIs1toaXcSYrdy54btzEiJOjn6L6mjqqkVC4IIxEeEePWzeV7Hl0pLWZIZz4ZePrlkZyWw8bYzySuq5b195byzt5zfbPqy8/Fwu5Xk6FASIkNIiAwhMtRGuN1CuN2K3WrBahXsFgsWi2ARsIp0xtbz77BnzIK+UfRmYUYcZ05L9Pl+vUnoaUBhl++LgDP62sYY0yEitUAicLzrRiKyBlgDMHny0C5SvL+svtsf42hjswhWi2CzCKF2K6E2CyE2C1GhNmLC7GQmRvCVKfGkxoYzKS6M6cnRzJwYNey1WCwWIX760s7vK2tbeKG4BEgkrwjuip7DhNgwAA6V1vF5S6lrwzYQwADWVqiOPY3k6FAA4owhbLKDnPyazhF6dJiNzNPPQiwWfvv9BVz1/z5hV3Fd53EXZ8QRN+0rIEJVfSsvlZTS4UxkT7nwk7jTeeR7i3utkwvwyPcWUlnfyg+e38bSv1QNqvTS7WcawGmTotlTMrxp9KdPimF3SV33O9sg3G6hudX15razDODkSew173SwvdBVBlqaFc+GNWcO+LOdfB0NXxytprK+tfN17EpEWJgRx8KMOO7+6iyqGlo5WNHAoYoGjlQ2cryhlZqmNspqW2hud9Dc5qC53UG7w+kegDhx6rLrPnPbedMCltB9xhjzBPAEuEouQ9nHitMn8uUvV/g0Ll/pOjoRTo5iLOIa2QTTx9pTRlJdvp85IYroMBv1LR1Eh9qYnRrN9oITLMmMJykqpMtzXCP/qsY24sJsHD7eyMwJUVjcQ32r1cJr3zuba//8GTuLapmfHstLt53Z+XtIigphSWZ8Z13ck8Q9bxg9WSyCxSJsKzhBh9OQm19DVWNbn9t3PYna9WfyvEH1JdxmYVFGLHlFtSxIjyWvsJbBTKiPtFuw9PGM5va+97Sz6OQbwEA/m0d/r2N/EqNCSYwKZdnUwSUVp9P1KRNcnz67fmLV62x4zzrE8z8D8SahFwNdG5vT3ff1tk2RiNiAWFwnR33OahGsFl05cLiSo0NZmnUymXZNHBaLhe0/v5hDlQ3MnBAFSJ/dJRbLyQQ8u5f2OKvVwsu3ndX7qLvLG4K3nSu9vQn0prfOHs/PND05kqqmdpxOJ9WNbdz30k72lJ4ckecU1GK1CAvS47CK8SqZd36KEWHj7Wdx2WMf97pddJiNxtaOXke7ESFW6t0lof5+tq76ex39wWIRLFpGCVrenBS14TopehGuxL0VuMEYs6fLNncA87qcFP2aMeYb/e1XT4oGXqD6v4fLm7gr61tZ9si7OIzrJOSW+/+lW7LrmvAXT46jrcPBjqLuJRKrRTpLSR4Wca1D3tjq6DbSjw6z0dTmIDsznue/ewbXPP5pZ/nEY1FGLC+uOZOa5na+/9/byM2vYU5qDHtL63C64/zbD84hKSqU5OhQr1+T0fo6qqEZ1klRd038+8A/cLUtPmOM2SMiDwM5xpjXgaeB/xKRQ0A1sMp34St/6Tq6Hk28iTshwk5EqKvEEhFqI6FHR4ynlbPDadhWcIJPfnyhe5ht+MGGHWzLr2FxZjwYV2knIvRkwn7ulqUcPt5IXLiNs3/7AQ6noam1gzd+uJxZE6MREV667Syu+dOn7Cyq7RyN7yqu40RLBykxYWxYcyZVjW0kRtq5/snPO0fYc1JjBp2UR+vrqHzPqxq6MeZN4M0e9z3Q5XYL8HXfhqbU0FU3tdPU5upiaWpzUN3U3i3p9SzdpMScHBFv6FIGMsaV/BMi7FQ3nZzwNDvV1WmT3WUfnmQOrlLTK7efzfGGVr7/39vY1uMcRNckPNiyk1J9GbDk4i9aclH+5Om99yTbDWuW+WXmqTf70JKI8qXh9qErNep4c8LVF6UKb/ahJRE1UjShqzFLE6kab3Q9dKWUGiM0oSul1BihCV0ppcYITehKKTVGaEJXSqkxQhO6UkqNEQGbWCQilUB+QA7euyR6LPcbZII9Pgj+GIM9PtAYfSHY44PhxZhpjEnu7YGAJfRgIyI5fc2+CgbBHh8Ef4zBHh9ojL4Q7PGB/2LUkotSSo0RmtCVUmqM0IR+0hOBDmAAwR4fBH+MwR4faIy+EOzxgZ9i1Bq6UkqNETpCV0qpMUITulJKjRGa0HshIneLiBGRpEDH0pWI/C8R+VJEdorIqyISF+iYAERkhYjsF5FDIvKTQMfTk4hkiMhmEdkrIntE5M5Ax9QbEbGKyHYR+XugY+mNiMSJyEvuv8F9InJmoGPqSUT+zf0a7xaR50UkLAhiekZEKkRkd5f7EkTkHRE56P433hfH0oTeg4hkAF8FCgIdSy/eAU43xszHdeHunwY4HkTECqwFLgHmAteLyNzARnWKDuBuY8xcYBlwRxDGCHAnsC/QQfTjD8AmY8xsYAFBFquIpAE/BLKNMafjugZyMFzf+FlgRY/7fgK8Z4yZAbzn/n7YNKGf6vfAfZy8oHvQMMa8bYzpcH+7BUgPZDxuS4FDxpgjxpg2YAOwMsAxdWOMKTXGbHPfrseViNICG1V3IpIOXAY8FehYeiMiscC5uC4IjzGmzRhzIrBR9coGhIuIDYgASgIcD8aYD4HqHnevBNa5b68DrvLFsTShdyEiK4FiY0xeoGPxwi3AW4EOAldiLOzyfRFBliy7EpEsYBHweWAjOcX/wTWQcAY6kD5MASqBv7jLQk+JSGSgg+rKGFMM/A7Xp+tSoNYY83Zgo+rTBGNMqft2GTDBFzsddwldRN5119d6fq0E7gceCOL4PNv8DFcZYX3gIh19RCQKeBn4kTGmLtDxeIjI5UCFMSY30LH0wwYsBh43xiwCGvFRmcBX3HXolbjefCYBkSLyzcBGNTDj6h33SUVg3F1T1BjzL73dLyLzcP0h5LkvKJwObBORpcaYskDH5yEi3wYuBy4ywTGJoBjI6PJ9uvu+oCIidlzJfL0x5pVAx9PD2cCVInIpEAbEiMhzxphgSkZFQJExxvPJ5iWCLKED/wIcNcZUAojIK8BZwHMBjap35SKSaowpFZFUoMIXOx13I/S+GGN2GWNSjDFZxpgsXH/Ai0cymQ9ERFbg+lh+pTGmKdDxuG0FZojIFBEJwXUS6vUAx9SNuN6hnwb2GWP+I9Dx9GSM+akxJt39d7cKeD/Ikjnu/weFIjLLfddFwN4AhtSbAmCZiES4X/OLCLITt128Dtzkvn0T8Fdf7HTcjdBHuT8CocA77k8RW4wxtwUyIGNMh4h8H/gHrq6CZ4wxewIZUy/OBlYDu0Rkh/u++40xbwYwptHoB8B69xv3EeDmAMfTjTHmcxF5CdiGqyS5nSBYBkBEngfOB5JEpAh4EHgUeFFEbsW1jPg3fHKs4PjUrpRSari05KKUUmOEJnSllBojNKErpdQYoQldKaXGCE3oSik1RmhCV0qpMUITulJKjRH/H9/VkNm3y2NOAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "samples = rejection_sampling(func, rv, k=15, n=100)\n", "plt.plot(x, func(x), label=r\"$\\tilde{p}(z)$\")\n", "plt.hist(samples, density=True, alpha=0.2)\n", "plt.scatter(samples, np.random.normal(scale=.03, size=(100, 1)), s=5, label=\"samples\")\n", "plt.legend(fontsize=15)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 11.1.5 Sampling-importance-resampling" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAgAElEQVR4nO3dd3xc1Znw8d8zM9Jo1KtlWdXGnWJsCzCmxoQspoYXQ0xxADshJJCyoSRhN5CQTRZ28yZvCN4ktMAu3UASQm9mCcXYso1csWxsq1mWZPVeZs77x8zIkqwykmY0o9Hz/Xzm49HMufc+Kn7umXOfe44YY1BKKTXxWYIdgFJKKf/QhK6UUmFCE7pSSoUJTehKKRUmNKErpVSYsAXrwKmpqSYvLy9Yh1dKqQlp8+bNR4wxaQO9F7SEnpeXR0FBQbAOr5RSE5KIFA/2ng65KKVUmNCErpRSYUITulJKhQlN6EopFSY0oSulVJgYNqGLSJSIbBSRQhHZKSI/H6CNXUSeE5F9IvKpiOQFIlillFKD86VssQNYZoxpFpEI4EMRed0Ys6FXmzVAnTFmpoisBO4HvhaAeJVSA2hsbKSqqoqurq5gh6LGICIigilTphAfHz+q7YdN6MY9v26z93ieR/85dy8DfuZ5/gLwoIiI0bl5VRiraGjj4301XL4wE4tFghZHY2MjlZWVZGZm4nA4EAleLGr0jDG0tbVRXl4OMKqk7tONRSJiBTYDM4G1xphP+zXJBEo9QXWLSAOQAhzpt5+bgJsAcnJyRhysUqFiR3kDqx/fRFVTB45IKxeemBG0WKqqqsjMzCQ6OjpoMaixExGio6PJzMzk0KFDo0roPl0UNcY4jTEnA1nAqSJywoiP5N7PQ8aYfGNMflragHeuKhXy3tlVyZV//IQIq4XclGh+985eXK7gfRjt6urC4XAE7fjKvxwOx6iHzkZU5WKMqQfWAxf0e6scyAYQERuQANSMKiKlQtj6z6v45v8UMCs9lr/cspQfnj+bPZVNvLHzcFDj0mGW8DGW36UvVS5pIpLoee4Azgc+79fsZeB6z/MVwHs6fq7C0brNpaTF2nn2piVMiYvi4pOmcVxaTNB76UqBbz30DGC9iGwDNgFvG2NeEZF7ReRST5tHgRQR2Qf8EPhxYMJVKni6nS4+3HuEc2anER3pvvxktQjfO28WeyqbeH1HcHvpSvlS5bINWDjA63f3et4OXOnf0JQKLYVlDTS2d3P27L7Xfy4+aRoPvLuX371bxPITpga14kVNbnqnqFI++qCoGhE4c2Zqn9etFuE7586kqLKZraV1QYoufLlcrmCHMGFoQlfKRx/srWZBViJJMZHHvHfWbHeS31pSP95hhZ2vf/3rzJo1i+9///vcfffdrFmzps/7xhhOPvlknnjiCZ/3eeuttx6zn3CkCV0pH9S3dlJYWn/McIvXlLgospIcmtDHqK6ujoKCArZv3055eTlr165l9erVfdo8//zz1NbWcs011/i839tvv52nnnqKffv2+TvkkKIJXSkffLjvCC4D58xOHbTNydmJbC3RIZexSEpKYteuXURFRfHCCy9QU1PDWWed1afNAw88wKpVq4iIiPB5v3l5eZx55pn84Q9/8HfIIUUTulI++KComrgoGwuyEgdtszAniUMN7VQ2to9jZOElOTmZn/70p1x77bWkpKSQnJzMt771LTo6OgDYt28fH3/8MStWrOiz3XPPPYeIDPhobGwE4IorruCpp54K6zF5TehKDcMYwwdFRzhzZio26+D/ZRbmuJO9DruMTklJCXV1dfz2t78lOTmZdevW8YMf/ICHHnqItWvXAvDuu+8SExPDggUL+mx79tln88knn/Q8Xn31VZKTk7nwwgt7bqFfunQplZWVbN++fdy/t/EStEWilZoo9lY1c7ixnXMGGT/3On5aPJFWC1tL67jghKnjFF34KCwsBNwXMO+77z4Ali1bxgcffMAbb7zBD3/4QzZv3sy8efOwWPqeWDMyMsjIcM+n09jYyHnnnccJJ5zAunXretocf/zxWK1WNm7ceMwJIVxoQldqGB8UVQMMekHUy26zMm9afEj00H/+953sOtQYlGPPnxbPPZccP+LtCgsLiY6O5o477ujz+uzZs9m4cSMAhw8fJjV18OsYLS0tXHjhhVgsFl555ZU+E5bZbDYSExM5fDh8bwDTIRelhrG9vIFpCVFMSxx+AqyF2YlsL2ug2xm+47SBUlhYyNlnn01KSkqf1w8fPsy0adMAaG9vx263D7h9e3s7l156KS0tLbzxxhvExcUd08Zut9PeHr7XOLSHrtQw9hxuYs7UY5PDQBbmJPL4xwfZU9nE8dMSAhzZ4EbTQw62wsLCYypampqaeO+99/j5z90LpSUnJw/Yw+7q6uKKK66goqKC//3f/yUpKWnAY9TX15OcnOz/4EOE9tCVGkKX08X+6hZm+5jQF+W4E0koDLtMJC0tLXzxxRd88cUXfV5/4IEHALj+evfcf3PmzOHAgQN92jidTq6++mr27NnDO++8w2BTc1dXV9Pa2srs2bMD8B2EBk3oSg2huKaFTqeLOem+JfSsJAepsZGa0Edo+/btuFwuvvjiC+69917Wr1/PT37yE37+85/z0EMPkZjoriA644wzKCkpobq6umfbb3/727z22mv86le/oqSkhA0bNrBhw4ZjqlkKCgoQEZYuXTqu39t40iEXpYZQVOlefXG2jwldRNw3GOmcLiNSWFhIfHw8r776KjfeeCO//OUvmTNnDuvWreOyyy7raXfuueeSnJzMG2+8wapVqzDG8Oyzz9LW1sbXvtZ3GeNLLrmEl19+uefrN954g3POOeeYMfpwogldqSHsOdyECMycEuvzNgtzknhndxUNrV0kRPt+N+NkVlhYyIknnshJJ53E5s2bB20XGRnJddddx7PPPsuqVav63Dg0FKfTyYsvvthTDhmudMhFqSEUVTaRlxJDVITV5228d5NuL28IVFhhp7Cw0Ofa8DvuuIP169dTVFTk8/7XrVuHw+Fg5cqVow1xQtCErtQQiiqbmDWC3jnArHR3+y+qmwMRUtgxxrB9+3ZOOukkn9pnZWXx2GOPUVFRMaJjPProo9hs4T0oEd7fnVJj0N7l5GBNKxeemDGi7abE2YmJtLJfE7pPfB026W2kPe2rr756RO0nKu2hKzWI/dUtOF3G5wuiXiLCjLRY9h9pCVBkSg1ME7pSgyiqbALw+aai3makxbC/WhO6Gl+a0JUaRFFlEzaLkJcSM+JtZ6TGUl7fRlunMwCRKTUwTehKDaKosokZaTFE2kb+32RGmvskcLBGe+lq/GhCV2oQeyqbRjx+7uVN6DrsosaTJnSlBtDS0U1pbduoE/r0VG9C10oXNX6GTegiki0i60Vkl4jsFJHvD9DmXBFpEJHPPI+7AxOuUuNjX9XIbvnvLzrSRkZClFa6qHHlSx16N3CbMWaLiMQBm0XkbWPMrn7t/mGMudj/ISo1/vaMocLFy13poj10NX6G7aEbYyqMMVs8z5uA3UBmoANTKpiKDjdht1nISY4evvEgZqTGsr+6BWOMHyNTanAjGkMXkTxgIfDpAG+fLiKFIvK6iAw4u76I3CQiBSJS0Hv6S6VCzYEjLUxPjcFqkVHvY0ZaDE0d3VQ3d/gxMhUs77//PiLCjh07gh3KoHxO6CISC7wI/MAY0/8+3S1ArjFmAfB74K8D7cMY85AxJt8Ykz/YJPRKhYKS2tYx9c4BZqS553TRShc1XnxK6CISgTuZP2WMean/+8aYRmNMs+f5a0CEiAy+kqtSIczlMpTUtpKbMsaEnqqli2p8+VLlIsCjwG5jzG8GaTPV0w4ROdWz3xp/BqrUeKlq6qCj2zXmHnpmogO7zaIXRn20c+dOLrjgApKTk4mJiWHevHmsXbsWgFdffZXzzz+fKVOmEB8fz5IlS3jrrbf6bP+zn/2M1NRUPv30U/Lz83E4HJx55pkcOHCAqqoqvvrVrxIbG8u8efN47733+mybl5fH7bffzi9+8QumTp1KbGws1157LQ0NQ0+B7HK5uO+++5g5cyZ2u53Zs2fzxBNP9Gnz4YcfctZZZxEfH098fDwnn3wy69at88NP7Fi+VLmcAawCtovIZ57X7gJyAIwxfwRWAN8WkW6gDVhp9EqQmqBKalsByBnFLf+9WSzC9NQYDmjpok8uueQS5s2bx5NPPondbmfPnj09szAeOHCASy65hNtvvx2LxcLrr7/O8uXL+eCDDzjjjDN69tHa2spNN93EnXfeSUxMDN/73vdYtWoVdrud5cuX853vfIf/+I//4Morr6S0tJTo6KMn7WeeeYaZM2fy8MMPU1FRwZ133sk3vvGNIZPvd7/7XZ544gnuvvtuFi1axNtvv83q1atJSUnh4osvprGxkYsvvpjLLruMu+++u2eq4Pr6AC1RaIwJymPx4sVGqVD0/KYSk/ujV8yB6uYx7+vbTxaYc/9z/diDGsKuXbv8uj+n02WqGtuNy+Xy636HUl1dbQCzbdu2Yds6nU7T1dVlvvKVr5gbb7yx5/V77rnHAOb999/veW3t2rUGMD//+c97Xtu5c6cBzGuvvdbzWm5urklKSjJNTU09rz355JNGRHp+vuvXrzeA2b59uzHGmL179xoRMY8//nif+FatWmXy8/ONMcZs2rTJAKaxsXEkP44hf6dAgRkkr+qdokr1U1LbikVgWqJjzPuanhpDSW0rnd0uP0QWeC6X4eqHN3D6v7/Lyoc24HKNzwft5ORksrOzufnmm3nuueeoqqrq835ZWRnXX389mZmZ2Gw2IiIieOutt45ZtSgyMpKzzjqr5+uZM2cCsGzZsmNeKy8v77Pt+eefT2zs0cVMLr/8cowxbNq0acCY3333XSwWC5dffjnd3d09j/POO4/PPvsMp9PJcccdR2xsLNdccw1/+9vfAtcz99CErlQ/JbWtTEt0jGpSrv5mpMbi9FxknQhqWjrZXFxHt8uwubiOmpbOcTmuxWLhrbfeYurUqaxevZqpU6dy1llnsXXrVlwuF5deeikff/wx9957L+vXr2fTpk0sX76c9vb2PvuJi4vDYjn6e4uMjAQgMTHxmNf6bztlypQ+X0dHRxMbGzvoykhHjhzB6XSSkJBAREREz+OGG26gu7ubiooKkpKSePvtt+nq6uKqq64iLS2Niy66iP3794/+hzUEXbFIqX6Ka8Ze4eJ1dJKu5hEtNB0sqbGRLM5NYnNxHYtzk0iNjRy3Y8+dO5cXX3yRrq4u/vGPf/CjH/2Iiy66iPfff5+tW7fy+uuvc8EFF/S0b2tr8+vx+38qaG1tpbm5mYyMgVesSk5Oxmaz8dFHH/U5iXh5TxBLlizhjTfeoK2tjXfeeYcf/vCHXHPNNWzYsMGv8YMmdKWOUVLbyj8dP9Uv+8r1XFidKD10EeGZby6hpqWT1NhIPMVr4yoiIoJly5b1JD5vD9lut/e0KS4u5qOPPvJ5HVJfvP322zQ3N/cMu/zlL39BRMjPzx+w/bJly3A6nTQ0NHD++ecPu3+Hw8Ell1zCjh07+Pd//3e/xd2bJnSlemlq76K2pXPMJYteSdERxERaKavzb28ykCwWIS3OPnxDP9q2bRu33347X/va15gxYwZ1dXXcf//9LFiwgCVLlpCVlcVtt93GL37xC5qamrjnnnvIzPTvDCQOh4OLLrqIO+64g4qKCu644w4uv/xy5s+fP2D7OXPmcPPNN7Ny5UruvPNO8vPzaW9vZ+fOnRQVFfHII4/w6quv8thjj/HVr36VnJwcysvL+dOf/tRnTN+fNKEr1Yu3J+2vIRcRITs5mrK6idFDD5apU6eSnp7OL3/5Sw4dOkRiYiJf+tKXuP/++7Hb7bz00kvccsstrFixgqysLP7lX/6F999/36+34a9cuZK4uDjWrFlDc3Mzl156KX/4wx+G3Gbt2rXMnj2bhx9+mLvvvpv4+Hjmz5/PmjVrAPcFWBHhrrvuoqqqirS0NC6++GJ+9atf+S3u3sQEqVw8Pz/fFBQUBOXYSg3mjR0V3PzkFl757pmckJngl31+44kCSmtbefOfz/bL/vrbvXs38+bNC8i+J4u8vDxWrFjBr3/962CHAgz9OxWRzcaYAceBtMpFqV6Ka7w3Ffmnhw6QleSgtK5VZ11UAacJXaleimtbSYqOID4qwm/7zE6OprXTSe04lQCqyUvH0JXqpbS2dcy3/PeXneS+Qam0ro2U2PG92Kh8c/DgwWCH4BfaQ1eql+KasU+b21+2Z396YVQFmiZ0pTy6nC7K69vIDVBCL60NXOmijs+Hj7H8LjWhK+VRUd+O02X8ekEUINZuIyk6gtIA9dAjIiL8ftekCp62tjYiIkZ3DUcTulIexbXuaW793UMHyEqKpjRAd4tOmTKF8vJyWlu1kmYiM8bQ2tpKeXn5MfPK+EoviirlEYiSRa/sZAe7K5r8vl+A+Ph4AA4dOkRXV1dAjqHGR0REBOnp6T2/05HShK6UR2ltK5E2C+lxUX7fd3ZSNO/sqsLlMljGsPD0YLyr4ajJTYdclPIormklO8kRkISblRxNp9NFVVOH3/etlJcmdKU8yupbeypS/O1oLbqWLqrA0YSulEd5XRtZSWNfpWggR0sXNaGrwNGErhTQ3NFNXWsXWUmB6aFnepazC2QtulKa0JXC3TuHo4nX36IirEyJs+uQiwooTehKcfS2/EANuQA6L7oKOE3oSgHl9e4eeqCGXMB9YVSHXFQgaUJXCiira8NuswR0UeTs5GgqGtrocroCdgw1uQ2b0EUkW0TWi8guEdkpIt8foI2IyAMisk9EtonIosCEq1RglNW1kpnkCOiiyNlJ0biMe84YpQLBlx56N3CbMWY+sAS4RUT6r5q6HJjledwEDL0Qn1IhpqyuLaDDLQBZyVqLrgJr2IRujKkwxmzxPG8CdgP9l9u+DPhv47YBSBSRDL9Hq1SAlAWwBt0rO0lr0VVgjWgMXUTygIXAp/3eygRKe31dxrFJHxG5SUQKRKSgurp6ZJEqFSCtnd3UtnQGPKFPTYjCIkcvwCrlbz4ndBGJBV4EfmCMaRzNwYwxDxlj8o0x+WlpaaPZhVJ+F+gadK8Iq4WMBEfP8ZTyN58SuohE4E7mTxljXhqgSTmQ3evrLM9rSoW8srrAlyx6ZSY6eo6nlL/5UuUiwKPAbmPMbwZp9jLwdU+1yxKgwRhT4cc4lQoY780+2QEecgHITHLokIsKGF/mQz8DWAVsF5HPPK/dBeQAGGP+CLwGXAjsA1qBG/0fqlKBUVbfRqTNQmqsPeDHykpy8HJhO91OFzar3gai/GvYhG6M+RAYsjjXuNe9usVfQSk1nsrq2shMDMw86P1lJjpwugyHG9vHZYhHTS7aRVCT3niULHpleo6j4+gqEDShq0mvvK513BK6t1eulS4qEDShq0mtvcvJkebOgJcsemUkuNcr1QujKhA0oatJbTxLFsE9L3panF2n0VUBoQldTWrjMQ96f1lauqgCRBO6mtTGu4cO7koXHUNXgaAJXU1qZXVtRFiFKXGBr0H3ykxycKi+HZfLjNsx1eSgCV1NauX1bUwbpxp0r6ykaDqdLqqbO8btmGpy0ISuJrWyutZxq3DxykrUWnQVGJrQ1aRWWtvWM0/5eDl6c5FWuij/0oSuJi13DXoH2cnj20P3fiLQShflb5rQ1aQVjAoXgBi7jaToCK10UX6nCV1NWqVBqEH3ykzSedGV/2lCV5OWN6FmJ4//rIeZiXpzkfI/Tehq0iqrayXSaiFtHOZB7y8rKZryujbcM08r5R+a0NWkVVbXRmbS+Nage2UmOmjrclLb0jnux1bhSxO6mrTKasdv2tz+vKWLOuyi/EkTupq03AtbBGfVIO+JRCtdlD9pQleTUmtnNzUtnUHroWcluk8kWumi/EkTupqUglnhAhDvsBFnt+ndosqvNKGrSSkY86D3JiJkJUdTqj105Uea0NWkdPQu0eAkdIDsJAeltdpDV/6jCV1NSqW1rdhtwalB98pKiqZMa9GVHw2b0EXkMRGpEpEdg7x/rog0iMhnnsfd/g9TKf9yV7g4EBn/GnSv7GR3LXqN1qIrP/Glh/44cMEwbf5hjDnZ87h37GEpFVjBLFn08k7bq8Muyl+GTejGmA+A2nGIRalxU1rXOu7T5vbnrbDRC6PKX/w1hn66iBSKyOsicvxgjUTkJhEpEJGC6upqPx1aqZFpau+ivrUr6D30LF3oQvmZPxL6FiDXGLMA+D3w18EaGmMeMsbkG2Py09LS/HBopUbOe7v9eK9U1F+M3UZyTCSltdpDV/4x5oRujGk0xjR7nr8GRIhI6pgjUypAvAk0mCWLXtlJDu2hK78Zc0IXkaniKRUQkVM9+6wZ636VCpRg31TUW1ZytF4UVX5jG66BiDwDnAukikgZcA8QAWCM+SOwAvi2iHQDbcBKo4W1KoSV1bURHWklOSYy2KGQleTgrZ2HcboM1iBM46vCy7AJ3Rhz9TDvPwg86LeIlAqwUs+0ucGsQffKToqmy2moamonIyH4nxjUxKZ3iqpJJxRq0L16Shf1wqjyA03oalIxxlBS20pOkGZZ7C/bM46v4+jKHzShq0mlpqWT5o5u8lJCI6FPS/QkdK10UX6gCV1NKsU1LQDkpsYEORK3qAgr6fF2XehC+YUmdDWpHDzi7gnnpYRGQgf3hVEdclH+oAldTSrFNS1YLUJmYuhUlGQnR2sPXfmFJnQ1qRysaSUz0UGkLXT+9LOSHFQ0tNHldAU7FDXBhc5ftVLjoLimhdwQuSDqlZ0UjcvAoXrtpaux0YSuJpWDNa0hl9Czkr2zLmpCV2OjCV1NGvWtnTS0dYXUBVHQhS6U/2hCV5PGwRp3wswNsYSekRCF1SJai67GTBO6mjS8NeihclORl81qITPR0XPCUWq0NKGrSePgkVZEjs6fEkryUmM4eKQl2GGoCU4Tupo0imtayIiPIirCGuxQjjE9JZqDR1rQmafVWGhCV5PGwZqWkBs/98pLjaGl00l1c0ewQ1ETmCZ0NWkU17SSlxp6wy0A0z1zy3inJlBqNDShq0mhqb2LmpbOkO2hH03oOo6uRk8TupoUir0liyF4QRQgM9GBzSIcqNGErkZPE7qaFIpDtAbdy2a1kJMcrT10NSaa0NWkcNA7D3qI1aD3lpcawwFN6GoMNKGrSaG4poW0ODsx9mHXRQ+avJQYimtatXRRjZomdDUpHKxpDbk7RPubnhpNW5eTykYtXVSjowldTQoHj7SQkxya4+deeZ5KFx12UaM1bEIXkcdEpEpEdgzyvojIAyKyT0S2icgi/4ep1Og1tHZR1dTBrPTYYIcyJO8skAe10kWNki899MeBC4Z4fzkwy/O4CfjD2MNSyn+KqpoAmJMeF+RIhjYt0UGk1aKVLmrUhk3oxpgPgNohmlwG/Ldx2wAkikiGvwJUE8Oew00hu+LOnsPuhB7qPXSrRchJidYhFzVq/rjknwmU9vq6zPNaRf+GInIT7l48OTk5fji0CiZjDB9/UcOD7+3jk/01WC3CJSdlcPO5xzF3anyww+uxt7KJmEhrSC0MPZi8lBgdclGjNq41XMaYh4CHAPLz87U2awLr6HbyjScK+MfeI0yJs3PXhXOpauzg6Y0l/PWzQ6w5czo/vXh+sMMEoKiymVnpcYhIsEMZ1vTUaP6xtxqXy2CxhH68KrT4I6GXA9m9vs7yvKbC2H2vf84/9h7hXy+ax3VLcnumpL112Uzuf+NzHv3wACdlJXDZyZlBjhSKKpv48rz0YIfhk7zUGDq6XVQ0tk+ITxQqtPijbPFl4OueapclQIMx5pjhFhU+3t1dyZ8/OsgNS/P4xlkz+swvnhgdyS8uO4FT8pK466Xt7K9uDmKkcKS5g5qWzpAfP/eanqKTdKnR86Vs8RngE2COiJSJyBoRuVlEbvY0eQ3YD+wDHga+E7BoVdBVNrZzxwvbmJcRz4+Xzx2wjc1q4YGrFxJps3DL01tp73KOc5RHFVV6KlymhnaFi5fWoquxGHbIxRhz9TDvG+AWv0WkQpYxhn9+7jPaOp38/uqFQ678k5Hg4P9etYDVjxdw3+uf87NLjx/HSI/aW+n+hBDqJYteU+OjsNu0dFGNjt4pqnz2/p5qPv6ihrsunMvMKcMPYSybm861p+Xw5IZiSmuDs3DDnsomEhwRpMXZg3L8kbJYhBlpsewL8lCVmpg0oSufGGN4cP0+MhMdrDzV95LTW5fNxCLCf73/RQCjG9zeyibmTJAKF6856bEUeWrnlRoJTejKJxv217K5uI5vnTODCKvvfzYZCQ6uOiWLFzaXUj7ONx4ZY9hzuGnCXBD1mjM1nkMN7TS0dQU7FDXBaEJXPlm7fh+psXauys8evnE/3z53JgB/HOdeelVTB43t3RPmgqjXnKnuE5D3gq5SvtKErob1WWk9H+47wjfPmj7khdDBZCY6WLE4m+c2lVLRMH699J5b/qdMtITuvsv2cx12USOkCV0N68H39pHgiODaJbmj3sd3zj0OlzHj2kv39nBnT7Ahl2kJUcRF2XQcXY2YJnQ1pC+qm3lndyU3LM0jdgyr/WQnR/PVhZk8X1A2bmPDRZVNpMZGkhI7MSpcvESEOelxPZ8wlPKVJnQ1pOc2lWKzCNcuGftkatefnkdbl5OXtpT5IbLhFVU2M3uC1J/3N3tqHJ8fbtTl6NSIaEJXg+rsdvHi5jLOmzeFKXFRY97fiVkJLMhK4KlPSwKeqFwuw97Kpgmb0OdOjaOxvZvDje3BDkVNIJrQ1aDe2V1JTUvniOrOh3Pdklz2VTWzYf9QU+yP3d6qZlo6nZyUlRDQ4wSK985WvTCqRkITuhrUMxtLmJYQxdmz0vy2z0sWTCPBEcGTnxb7bZ8D2VxcB8Di3KSAHidQvKWWemFUjYQmdDWg0tpWPtx3hCvzs7H6cV7uqAgrKxZn8eaOw1Q1BW44YXNxHamxkeQkRwfsGIGUGB1JerxdL4yqEdGErga0rsC9CNVVp4z8RqLhXHtaDt0uw3MbS4dvPEpbSupYlJM0oW7572/O1HgdclEjogldHaPb6eL5gjLOng38jm0AABigSURBVJUWkEUWZqTFcubMVJ7ZWILT5f+LozXNHRw40jJhh1u85k6NY191M91OV7BDUROEJnR1jA/3HeFwYzsrA9A797r2tBwONbTzQVG13/e9paQemLjj515z0uPo7HZxsCY4M1WqiUcTujrGS1vKSYyOYNm8KQE7xpfnp5Maa+fpjSV+3/fm4joirMIJmROzwsXLe2FUx9GVrzShqz4a27t4c+dhLjlpGnbbyOdt8VWE1cKV+Vm893kVhxv8e3F0S3EdJ2QmjGremVAyc0osFoE9hxuDHYqaIDShqz5e315BR7eLKxZnBfxYK0/JxukyPRdg/aGz20VhWT2Lcyb2cAu4K4LyUmPYVaE9dOUbTeiqjxe3lDMjLYYF43BDTm5KDGfOTOXZTaV+uzi6q6KRjm7XhB8/91qQlchnpfU6BYDyiSZ01aO0tpWNB2q5YlHWuJX7XX1qDuX1bfxjr38ujnpvKFoUJgl9UW4SR5o7KKsb38VB1MSkCV31+MvWcgC+ujBz3I55/vx0UmIiecZPF0e3FNeRleQgPX7sc8+EAu/QkfdEpdRQNKErwL1c20tbyjh9RkpAas8HE2mzsGJxFu/sruLQGJeoM8ZQUFwbNsMt4K50iYm0sqVEE7oaniZ0BbjvrDxY08rli8avd+616vRcjDE88cnBMe1nX1UzlY0dnJKX7Je4QoHVIpyck6g9dOUTTegKgHUFZURHWrnoxIxxP3ZWUjTLT8jgmU9LaOnoHvV+/l54CIvAPx0/1Y/RBd/inCR2VzSO6WejJgefErqIXCAie0Rkn4j8eID3bxCRahH5zPP4hv9DVYHS2tnNK9sqWH5CBjFjWJVoLFafOZ3G9m5e2Dy6xS+MMbxceIilx6WSFjexVigazsLcJFwGCsvqgx2KCnHDJnQRsQJrgeXAfOBqEZk/QNPnjDEnex6P+DlOFUBv7jxMc0c3V+YHvvZ8MItzkzg5O5E/f3QA1yhKGHeUN3KwppVLFoz/J4xAW5TtviawRYdd1DB86aGfCuwzxuw3xnQCzwKXBTYsNZ7WFZSRkxzNqUEee15z5nQO1rTy7udVI97279sOEWEVLjg+/BJ6QnQEs6bE6ji6GpYvCT0T6H0rX5nntf6uEJFtIvKCiAw4q5OI3CQiBSJSUF3t/0mZ1MiV1rby8Rc1XLEoC4sf5z0fjeUnTGVaQhSPfrh/RNu5XIa/Fx7inNlpJERHBCi64FqUk8TW0vpRfXpRk4e/Bkz/DjxjjOkQkW8BTwDL+jcyxjwEPASQn5+vf5kh4KUt5YjAFYvHv7qlP5vVwo1nTOeXr+3m431HWDoz1aftNpfUUdHQzo+Xzx3ZAQ9tHVn7aQtH1t6PFucm8VxBKfuPtDBzSmzQ4lChzZceejnQu8ed5XmthzGmxhjT4fnyEWCxf8JTgeRyGV7YUsrS41LISgqNlX1WnZ5LVpKDn/19p8/zgP+98BBRERa+PC89wNEFj/fOV61HV0PxJaFvAmaJyHQRiQRWAi/3biAivQcuLwV2+y9EFSgbDtRQWtvGinGYiMtXURFW/vWi+RRVNvPkhuHXHe12unhtewXnzUsPWoXOeJiRGkOCI0IvjKohDZvQjTHdwK3Am7gT9fPGmJ0icq+IXOpp9j0R2SkihcD3gBsCFbDynyc3FJPgiAi5C4n/dHw6Z85M5TdvF1HT3DFk2z99sJ8jzZ0hdVIKBItFyM9N4uMvanSiLjUon+rQjTGvGWNmG2OOM8b80vPa3caYlz3Pf2KMOd4Ys8AY8yVjzOeBDFqN3eGGdt7cWcnXTsnGERla84aLCPdcMp/WTie/fqto0Ha7DjXy/94p4uKTMvjSnMAtxhEqls2bQkltK3urmoMdigpReqfoJPX0p8W4jOG603KDHcqAZqXHcf3SPJ7dVMLTnx47cVdHt5MfPv8ZidGR/OKyE4IQ4fjzXiN4e1dlkCNRoUoT+iTU2e3i6Y2lLJszhZyU0LgYOpDbvjKbc2ancddftvOfb37eZ6jhd+/s5fPDTdx/xYkkxUQGMcrxkx4fxYKsBN7ShK4GEb5XkdSgXt9RwZHmDr6+NC/YoQwpOtLGI1/P51//uoO1679gz+EmYu02dh5qZG9VM1/Lz2bZ3PCtbBnI+fPT+fVbRVQ2tofNFMHKf7SHPgn99yfF5KVEc5aPdd7BZLNa+Pf/cyJ3/NMc3t9TzYb9teQkR3Pb+bO559KBZqAIb+fPd0889u7ukd9Nq8Kf9tAnmR3lDWwuruOnF88P+p2hvhIRbvnSTG46ewYR1iD2QUZyI1KAbkKanR5LTnI0b+86zDWn5QTkGGri0h76JPNf7+8j1m6bkGV+QU3mIUJE+PK8dD76okan01XH0P8hk8juikZe236YG8/II8ERnnOeTAbnz0+ns9vFB0U6H5LqSxP6JPLAu3uJs9tYc+b0YIeixuCUvCQSoyO0fFEdQxP6JLHrUCOv73D3zhOjJ0eZX7iyWd3z1ry1q5Km9q5gh6NCiCb0SeJ37xYRF2VjzZkzgh2K8oOvn55Lc0c3z20qHb6xmjQ0oU8COw818ObOSlafMT1s5wufbE7KSuTUvGT+/NFBn2elVOFPE3qYc7kMP3t5JwmOCFbr2HlYWXPWdMrr2/TOUdVDE3qYe3pjCZsO1vGvF83TypYw8+V56eSmRPPIP0a2wpMKX5rQw9jhhnbuf/1zzpiZMiHrztXQrBbhxqV5bCmp1/VGFaAJPazd/bcddDpd/OryExGZGHeFqpG5Mj+buCjbiNdhVeFJb/0PU69sO8Rbuyr5yfK55KbEBDucyWec1iuNsdu4/vQ8Hly/b0TrsKrwpAk9DG0rq+eOdds4OTsxsDcRTaBFlsPZLV+ayWvbK7jjhW288YOziIvSayWTlSb0MFNW18rqxwtIiY3k4a/nYxvJ/CcjTdDKf8ZwcnREWvn1VQtY8YeP+bdXdnP/ipP8HJyaKHQMPYw0tnex+vFNdHQ7+fMNp5AWZw92SGqcLMpJ4uZzjuO5glLe+1zLGCcrTehhoqKhjese+ZT91S386brFzEqPC3ZIapx9/8uzmDs1jjtf2EZRZVOww1FBoEMuYaDgYC03P7mFts5u/njd4tC9MDYBxtxdLmho7yLREYEx7ufxUTYa27tJdEQQysVCdpuVB69ZyDUPf8qVf/yEx244hcW5ScEOS40jTegTWEe3k//+uJj/ePNzMhMdPP3N05gd4J5574Tn7+Q25n0Pc8IYbv8uF9z11+3srmhk7tR4wPB5RSP2CCsdXS7mZsRx5wVzSY6ODNnEPnNKHC9+eymrHv2Uax/ZwB+uW8yX5kwJdlhqnEjvhXfHU35+vikoKAjKsSc6l8vwt8Jy/u9bRZTVtbFs7hR+e9XJfeZpcbkMNS2dpMZG+l6D7kmIgyW+3glvXkY8v/rqiVgsfd8fbUIebt++7mOw4w+2/97b1Ld1ccOfN+J0GawWAWNw9vvvIcDxmQkDbi8C3U5DaV0rucnRWCwy5pPUoNsP8+mluqmDG/68kV0VjfyfhVn88CuzyUx0jDwAFXJEZLMxJn+g93zqoYvIBcDvACvwiDHmvn7v24H/BhYDNcDXjDEHxxK06ssYw+6KJl7dfohXtlVQXNPK/Ix4nlh9ImfPSu2TtF0uw9UPb2BzcR2Lc5N45ptLfF5ubqjE2tDexe6KRpwuw+6KRhrau0jynET6927vvGBOn57sQEMZvZPUUPv2RbfT8KOXtrGvqnnAE8JA+4+z2/jRi4UUVTYzOz2WHy+fy9ypcXxe0cjcqXG4nC52Vzb3/T1Az/YJURF9vufbvjKLW5/+jNbObmIibTy5+lR++vedPp+k+ifvoU5CNU0dQ56s0+LsPPet0/n9u3v588cH+fu2Q1x3Wi6XL8zkhMx4vdEsTA2b0EXECqwFzgfKgE0i8rIxZlevZmuAOmPMTBFZCdwPfC0QAYc7YwyNbd0camijoqGN3RVNbC9rYHt5A+X1bVgtwtLjUrjtK3O4+MSMARN1TUsnm4vrmGu+oL1EaNhv9Tk5DpVYEx0RzMuI70kwiY6IniTkPuG4t9t5qIEbH9vIrPQ4fnLhXJKiI/mXv+7oO5RxuMl9Ae+f5pIcEzngvsG3Xr/LBT96cRt7PBcCBzoh9N9/fJSNO18oZI8nYe+pbObGxwuwWwVE6Ha6KK5tPeZYAsydGtfTo+/9Pa9+/OgnzpbObnZUNPp8khooeQ/0u/CeRF44VDHsyTrWbuMnF87j+qV5/PbtIh7/+ACPfXSAqfFRfGluGvMz4jluSiwz02JJibW7P5WoCc2XHvqpwD5jzH4AEXkWuAzondAvA37mef4C8KCIiAnAeM62snqe2Vji792OWe/v1BgwGM+/4DLu5y5jcLrcjy6nodPpoqPLSUe3i+aObhrbumhs76K9q+90qLkp0SzMSeQ7XzqOC46fSkrs0OWIqbGRLM5Nor1E+iTH/gZKlr0T38wpsSRE2XC5oK6tEwF+cenxlNW3kZscjTFHk9DMtBh37/ZwE06XwQXsqWzihj9vYk56LPuqW3C6DJ8fbuwZyth5qJHVT2zqSWDeJBYfZaO+rYvYSCs/+ct29lY2MW9aQp8ebu/Y61o7e5I5wMy02GO+ZxF69p/oiKC2pbMnmffW4Rlj6f+e3WrhuLSYntedLoMxhrlT49hd0YhrgL/0WLul52c5d2o8xvP7b2zvJjbS2vNztFiE2pYOdh9qxGmOJu+BTnLek0i3K4XNxXXUtHQOW546LdHBf165gB8vn8v6PdW8u7uSVworeGbj0bnULQKJ0ZEkx0QSY7fhiLDgiLASYbVgswo2iwWrRRABiwji+ZkK0udEqx1/35wzO40LTsjw+359SeiZQO9Z9MuA0wZrY4zpFpEGIAU40ruRiNwE3ASQkzO6FcurGjt4d3fVqLYNtD5/2J4/dPcfvmCxuP8jWC2CzSJYRLBHWLHbLMRF2ZiWGEWcPYJ4h40pcVFkJEaRkeBgZlrsiOcwFxGe+eYSaloWDfqxfLBhGQH+7eYFXPXQJ/y1tJ7SV9oxBjZ5Jn+Ks9to7ewmP8/KAysXsu5QBU5XCoWHYWF2Ag/duZhv/c9mPitr6DnWzsNwck4iO8saWJSTBMZQcLAOF4ATdh4SbkuYT1qcnQRPXAUHa3FEWmnuSAPS2FFOT5v+sf9u5cnsMEf/1B69cRmScOx4sQVIcBmOtHTiinX12WZY3bC5AiCNwnLY8VITuyqaibGn0+wauKoofsap3Hv8mWwqqeU3b35O/qNHiLHbaO1y9pzU4qJa2PTj8zj3T+/Q7HLf1XtKVhJdU06CODu/+s7Jfa6FJBpDVI4Tm+d7T431ffWplFg7KxZnsWJxFsYYDje2s6+qmf3VLRxp7qC2pZPalk5aO520dTmpaemks9tFt8vQ7XTh8nRKXC6D4WjHxStIl+MmpEBdzxjXKhdjzEPAQ+C+KDqafXx5fjpfnp/u17jCkcUiQ/bcvMMy3S5zTE+vrq2LbWUNOA1sLq6j9y+qybPSfMHBWsBwUlYCW0vqAdhW1oDFYuG/rl3E0vvX92yzICuBF751OrWtXaTGRuJ0Gq780ydsLXVvtygnsScxeeNyGmjucB7dR/axbbyxg/tE09TRTZzdRuogn2B6nwgWZicQa7f2OUZ/UTZo7x74vZ0V7p76YNvHRFpJjIpg8a/eoanXTrw/v56v27tZX1TdZz9tXd0sve+9nhNt79+jiPDUmtPYV+0e9x/tWLiIkJHgICPBwVmz0ka1DxV6fKkjKAeye32d5XltwDYiYgMScF8cVSHKOyxjs8gxPb3U2EgW5SRiFXey7V3LHGu3AuA08N1nP+P5by5hYXYCFujZT//pBv60ajFWq4W0ODsiQl1bF9vL3T14q0V48JpFPYnJG5dVIC7KhtUiLMxJ5IWbTz+mjTd2i4i71wu0djmpbR14nc3eJ4ItJfW0eYa2rBZhdlp0T7tT85J4/Xtn0jv3Hp8xsnLQ1k4n7+2p7pPMvawWwerJw3FRNs6fn0ZclLtvFWO3sutQU58TLbhPRtVNHTidLq599FMu/v2HXP3wp7gGGutRk5YvPfRNwCwRmY47ca8ErunX5mXgeuATYAXwXiDGz5X/HB2WOba00f2bc18cFIuFp9ecRk2rewzdGMPS+9fjdBm2FNdR29pFhNWKWNztjXFXWJyal9QzJDIlPqrPsb0J2ft+/x6oN67k6IieXn3v+PrHDpDfa3+DDUP0Pu6inEQQYYtnm6e/cRpHWtzfo3dYJzbKRlN7N7F2K3/9zlKufXQTm0vqcERYenrUMZEWWjrdJ4ZTcpPocrnYVtZAjN3Grc9uxSocU/r4yq1nMDMtli+OtDA7PRaLxcLWfz2ffdXNzJoSwzWPbOzzvfT+ZHFSVgKFZQ04B/hkpdSwCd0zJn4r8CbussXHjDE7ReReoMAY8zLwKPA/IrIPqMWd9FWIG2xYpqalky0ldT1Ju66ti3RPUjbG9EmeIvRp600wz950+qB18EOdTPrHNViy6h/7UPsb7LjG0Geb9F4nntpW91gyQFuXi/p2J8/edPREc6Slk9rmDi7+/YcAWAXWXruI1Fg7RZVNXPT7D3vq2Z9ecwq/ebuIrSX1LM5NYm6Gu2xwbkZ8z/FsNkvP1/2/lyPNHT2fLApL6zkpO5HtZQ0jHkNX4c+nMXRjzGvAa/1eu7vX83bgSv+GpoKlfw+6d9IYqHc8UNvhxvCHe3+kfN1f73Yig58wUmMjj+n1ixzdNj0+iilxdvLzkvt80hAR5kyN67Pt6cel8vyMVJ9v9Or/vfT/fTz9jdMG/OSilN4pqgY0kjtNR3VX6gTgy/c1WBt//0zC9WesRm7Md4qqyWckPWh/97ZDhS/f12BtgvUJRE1uOn2uUkqFCU3oSikVJjShK6VUmNCErpRSYUITulJKhQlN6EopFSaCVocuItVAcVAOPrBU+s0OGWJCPT4I/RhDPT7QGP0h1OODscWYa4wZcEa1oCX0UCMiBYMV64eCUI8PQj/GUI8PNEZ/CPX4IHAx6pCLUkqFCU3oSikVJjShH/VQsAMYRqjHB6EfY6jHBxqjP4R6fBCgGHUMXSmlwoT20JVSKkxoQldKqTChCX0AInKbiBgRGXg59yARkf8Ukc9FZJuI/EVEEoMdE4CIXCAie0Rkn4j8ONjx9Cci2SKyXkR2ichOEfl+sGMaiIhYRWSriLwS7FgGIiKJIvKC529wt4icHuyY+hORf/b8jneIyDMiEjX8VgGP6TERqRKRHb1eSxaRt0Vkr+ffpKH24StN6P2ISDbwFaAk2LEM4G3gBGPMSUAR8JMgx4OIWIG1wHJgPnC1iMwPblTH6AZuM8bMB5YAt4RgjADfB3YHO4gh/A54wxgzF1hAiMUqIpnA94B8Y8wJuJfMDIXlMB8HLuj32o+Bd40xs4B3PV+PmSb0Y/0WuBMIuavFxpi3jDHeZeQ3AFnBjMfjVGCfMWa/MaYTeBa4LMgx9WGMqTDGbPE8b8KdiDKDG1VfIpIFXAQ8EuxYBiIiCcDZuNcPxhjTaYypD25UA7IBDhGxAdHAoSDHgzHmA9xrLfd2GfCE5/kTwFf9cSxN6L2IyGVAuTGmMNix+GA18Hqwg8CdGEt7fV1GiCXL3kQkD1gIfBrcSI7x/3B3JFzBDmQQ04Fq4M+eYaFHRCQm2EH1ZowpB36N+9N1BdBgjHkruFENKt0YU+F5fhhI98dOJ11CF5F3PONr/R+XAXcBdw+3jyDG523zL7iHEZ4KXqQTj4jEAi8CPzDGNAY7Hi8RuRioMsZsDnYsQ7ABi4A/GGMWAi34aZjAXzzj0JfhPvlMA2JE5LrgRjU8464d98uIwKRbU9QY8+WBXheRE3H/IRR6FuHNAraIyKnGmMPBjs9LRG4ALgbOM6FxE0E5kN3r6yzPayFFRCJwJ/OnjDEvBTuefs4ALhWRC4EoIF5EnjTGhFIyKgPKjDHeTzYvEGIJHfgycMAYUw0gIi8BS4EngxrVwCpFJMMYUyEiGUCVP3Y66XrogzHGbDfGTDHG5Blj8nD/AS8az2Q+HBG5APfH8kuNMa3BjsdjEzBLRKaLSCTui1AvBzmmPsR9hn4U2G2M+U2w4+nPGPMTY0yW5+9uJfBeiCVzPP8PSkVkjuel84BdQQxpICXAEhGJ9vzOzyPELtz28jJwvef59cDf/LHTSddDn+AeBOzA255PERuMMTcHMyBjTLeI3Aq8ibuq4DFjzM5gxjSAM4BVwHYR+czz2l3GmNeCGNNE9F3gKc+Jez9wY5Dj6cMY86mIvABswT0kuZUQmAZARJ4BzgVSRaQMuAe4D3heRNbgnkb8Kr8cKzQ+tSullBorHXJRSqkwoQldKaXChCZ0pZQKE5rQlVIqTGhCV0qpMKEJXSmlwoQmdKWUChP/HzLNmEVpGiZwAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "samples = sir(func, rv, n=100)\n", "plt.plot(x, func(x), label=r\"$\\tilde{p}(z)$\")\n", "plt.hist(samples, density=True, alpha=0.2)\n", "plt.scatter(samples, np.random.normal(scale=.03, size=(100, 1)), s=5, label=\"samples\")\n", "plt.legend(fontsize=15)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 11.2 Markov Chain Monte Carlo" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAgAElEQVR4nO3dd3xc1Z338c9vpFHvxbas6l4AV+GYYuqSNdUQmjFx6ISShCwtgTwLgWwIsHk2T9h4SUwJJHgppoVQHJoJ1diWbbnigm31ZklWrzPn+WNmZPU6oxnN/N6vl14e3blz708j+XvPnHvuuWKMQSml1Nhn8XYBSiml3EMDXSml/IQGulJK+QkNdKWU8hMa6Eop5SeCvbXjpKQkk5WV5a3dK6XUmJSTk3PEGJPc23NeC/SsrCw2b97srd0rpdSYJCJ5fT2nXS5KKeUnNNCVUspPaKArpZSf0EBXSik/oYGulFJ+YsBAF5EwEdkoIrkisktEHuplnVAReVlEDojI1yKS5YlilVJK9W0wwxZbgLOMMfUiYgU+F5H3jDEbOq1zA1BtjJkqIsuBx4ArPVCvUqoXtbW1lJeX09bW5u1S1AhYrVbGjRtHTEzMsF4/YKAbx/y69a79Ob+6z7m7DPil8/GrwB9ERIzOzav8WElNE18cqOR781OxWMRrddTW1lJWVkZqairh4eGIeK8WNXzGGJqamigqKgIYVqgP6sIiEQkCcoCpwCpjzNfdVkkFCpxFtYtIDZAIHOm2nZuBmwEyMjKGXKxSvmJ74VFueH4zFXUtRIQEcd4JKV6rpby8nNTUVCIiIrxWgxo5ESEiIoLU1FSKi4uHFeiDOilqjLEZY+YBacAiETl+yHtybGe1MSbbGJOdnNzrlatK+bx/7Crlij99RUiQhczECH7/4X7sdu99GG1rayM8PNxr+1fuFR4ePuyusyGNcjHGHAXWA0u7PVUEpAOISDAQC1QOqyKlfNhHe8q45YUcZkyI4Y3bT+bOc6azt6yOd3eWeLUu7WbxHyP5XQ5mlEuyiMQ5H4cD5wDfdFvtLeAa5+PLgI+1/1z5o1dzChkfHcZLNy1mXHQYF8yZyNRxUfz+w/3YvNhKVwoG10JPAdaLyHZgE/CBMeZtEXlYRC5yrvMMkCgiB4A7gZ97plylvKfdZufzA0c4fXoy4SFBAARZhDvOnsb+8nre2eHdVrpSgxnlsh2Y38vyBzo9bgYud29pSvmWbQVHqWtu5/QZXc//nH9CCk98tJ8nPtrP+SekEOTFES8qsOmVokoN0qf7KrAInDIlqctyi0W4/cypHCivZ0t+tZeq8192u93bJYwZGuhKDdI/91UwLz2O2Ahrj+dOneYI+a0a6CP2gx/8gGnTpnHHHXfwwAMPcMMNN3R53hjDvHnzeP755we9zR/96Ec9tuOPNNCVGoSqhla2F9Vw+vRxvT6fFBVKRkIEW/OPjnJl/qW6uprNmzezY8cOioqKWLVqFddff32XdV555RWqqqpYsWLFoLd79913s2bNGg4cOODukn2KBrpSg/DZ/gqMgdOmJ/W5zvyMOA30EYqPj2f37t2EhYXx6quvUllZyZIlS7qs88QTT7By5Uqs1p6flPqSlZXFqaeeypNPPunukn2KBrpSg/DpviPERViZkxbX5zrz0+MorW2mpKZpFCvzLwkJCfz7v/87V199NYmJiSQkJPDDH/6QlpYWAA4cOMCXX37JZZdd1uV1L7/8MiLS61dtbS0Al156KWvWrPHrPnkNdKUGYIzh0/0VnDo1qd8RLPMz4gG0lT5M+fn5VFdX87vf/Y6EhATWrl3LT3/6U1avXs2qVasA+Oijj4iMjGTu3LldXnvaaafx1VdfdXy98847JCQkcN5553VcQn/yySdTVlbGjh07Rv1nGy1eu0m0UmPFnpI6KupaOG16/9NVzEqJISTYwtb8aq/O7TJW5ebmAo4TmI8++igAZ511Fp9++inr1q3jzjvvJCcnh1mzZmGxdG2LpqSkkJLieM9ra2s5++yzOf7441m7dm3HOscddxxBQUFs3LixxwHBX2igKzWAT/dXAHD6AIEeEmzhhNRYn2ihP/T3XewurvXKvmdPjOHBC48b8utyc3OJiIjgnnvu6bJ8+vTpbNy4EYDS0lKSkvo+j9HQ0MB5552HxWLh7bff7jJhWXBwMHFxcZSWlg65trFCu1yUGsDOohrS4sMZHxM24Lrz0+PYUVRDa7v/9tN6Sm5uLqeddhqJiYldlpeWljJx4kQAmpubCQ0N7fX1zc3NXHTRRTQ0NLBu3Tqio6N7rBMaGkpzc7P7i/cR2kJXagD7yuqYOaFnOPRmfkY8T39+iG9Ka/s9geppw2khe1tubm6PES11dXV8/PHHPPSQ40ZpCQkJvbaw29rauPTSSykpKeGf//wn8fHxve7j6NGjJCQkuL94H6EtdKX60dpu52BFA9PHDzbQHSHuC90uY0lDQwPffvst3377bZflTzzxBADXXOOY+2/GjBkcOnSoyzo2m42rrrqKvXv38uGHH9LX1NwVFRU0NjYyffp0D/wEvkEDXal+HDrSQLvdMGOQLfSU2DDGx4TqFaNDtGPHDux2O99++y0PP/ww69ev57777uOhhx5i9erVxMU5DpSnnHIK+fn5VFRUdLz21ltv5d133+WRRx4hPz+fDRs2sGHDhh6jWTZv3oyIcPLJJ4/qzzaatMtFqX7sLasDGHQLXUSYnx7P1gJtoQ9Fbm4uMTExvPPOO1x33XX8+te/ZsaMGaxdu5Zly5Z1rHfGGWeQkJDAunXrWLlyJcYYXnrpJZqamrjyyq63Mb7wwgt56623Or5ft24dp59+eo8+en+iga5UP/aV1hFkESYnRw76NfMz4li3q5TK+hYSo3o/gae6ys3N5YQTTmDOnDnk5OT0uV5ISAjf//73eemll1i5cmWXC4f6Y7PZeO211zqGQ/or7XJRqh97y+qYlBRJaHDQoF8zN93RPbCjqMZTZfmd3NzcQY8Nv+eee1i/fj379u0b9PbXrl1LeHg4y5cvH26JY4IGulL92FdWx4xBdre4TB0XBcC3FQ2eKMnvGGPYsWMHc+bMGdT6aWlpPPvss5SUDP6GIsYYnnnmGYKD/btTwr9/OqVGoKnVRn5VI9+bnzak1yVGhhATFszBinoPVeZfBttt0tlQW9pXXXXVkNYfq7SFrlQfDpTXYwxMHx81pNeJCJOTozioLXQ1yjTQlepDxwiXQQ5Z7GxyciQHj2gLXY0uDXSl+rCvrI6QYAuZCREDr9zNlOQoympbqG9p90BlSvVOA12pPuwtrWNqchTBQUP/bzI5yTHM8ZB2u6hRpIGuVB/2ldUN+grR7iYnO/rdtdtFjSYNdKV6UdPURklN86CvEO0uMzECER26qEbXgIEuIukisl5EdovILhG5o5d1zhCRGhHZ5vx6wDPlKjU69jtPiM6YMLQRLi5h1iDS4sN16KIaVYMZh94O3GWM2SIi0UCOiHxgjNndbb3PjDEXuL9EpUbfUOdw6c0UHbqoRtmALXRjTIkxZovzcR2wB0j1dGFKedO+0joiQ4JIjQsf9jYmJ0Vx6EgDdrtxY2VK9W1IfegikgXMB77u5emTRCRXRN4TkV5n1xeRm0Vks4hs7jz9pVK+5lBlI5OSIxHp+6bQA5mcHElTm43SWv+9Q04g+eSTTxARdu7c6e1S+jToQBeRKOA14KfGmO7X6W4BMo0xc4H/Bt7sbRvGmNXGmGxjTHZfk9Ar5QvyKxvITBj8DIu9cc3Q+K32o6tRMqhAFxErjjBfY4x5vfvzxphaY0y98/G7gFVE+r6Tq1I+zGY3FFY3kZE49AuKOpviGrqo/ehqlAxmlIsAzwB7jDH/1cc6E5zrISKLnNutdGehSo2W4qNNtNvNsK4Q7WxcdCiRIUE60mWQdu3axdKlS0lISCAyMpJZs2axatUqAN555x3OOeccxo0bR0xMDIsXL+b999/v8vpf/vKXJCUl8fXXX5OdnU14eDinnnoqhw4dory8nIsvvpioqChmzZrFxx9/3OW1WVlZ3H333fzqV79iwoQJREVFcfXVV1NT0/8UyHa7nUcffZSpU6cSGhrK9OnTef7557us8/nnn7NkyRJiYmKIiYlh3rx5rF271g3vWE+DGeVyCrAS2CEi25zL7gcyAIwxfwQuA24VkXagCVhujNEzQWpMyq9qBBhxC71jkq4j2kIfjAsvvJBZs2bxwgsvEBoayt69eztmYTx06BAXXnghd999NxaLhffee49zzz2XTz/9lFNOOaVjG42Njdx8883ce++9REZG8pOf/ISVK1cSGhrKueeey2233cbjjz/O5ZdfTkFBARERx37HL774IlOnTuWpp56ipKSEe++9lxtvvLHf8P3xj3/M888/zwMPPMCCBQv44IMPuP7660lMTOSCCy6gtraWCy64gGXLlvHAAw90TBV89KiH7mhljPHK18KFC41SvmjNhjyT+bO3TUFVw4i39ZMXt5iTf/ORG6rq2+7du926PZvNbsprm43dbnfrdvtTUVFhALN9+/YB17XZbKatrc1897vfNdddd13H8gcffNAA5pNPPulYtmrVKgOYhx56qGPZrl27DGDefffdjmWZmZkmPj7e1NXVdSx74YUXjIh0vL/r1683gNmxY4cxxpj9+/cbETHPPfdcl/pWrlxpsrOzjTHGbNq0yQCmtrZ2KG9Hv79TYLPpI1f1SlGlusmrasAaJKTEDn/IosvkpCiKjjbR1GpzQ2WeZ7cbrnpqAyf95iOWr94wakMuExISSE9P55ZbbuHll1+mvLy8y/OFhYVcc801pKamEhwcjNVq5f333+9x16KQkBCWLFnS8f3UqVMBOOuss3osKyoq6vLac845h6ioYxeSXXLJJRhj2LRpU681f/TRR1gsFi655BLa29s7vs4++2y2bduGzWZjypQpREVFsWLFCv72t795rmXupIGuVDcFVY2kx0cQZBn+kEUX10iXQ2Ok26WyoZWcvGra7YacvGoqG1pHZb8Wi4X333+fCRMmcP311zNhwgSWLFnC1q1bsdvtXHTRRXz55Zc8/PDDrF+/nk2bNnHuuefS3Nx1SGh0dDQWy7FYCwkJASAuLq7Hsu6vHTduXJfvIyIiiIqK6vPOSEeOHMFmsxEbG4vVau34uvbaa2lvb6ekpIT4+Hg++OAD2trauOKKK0hOTub888/n4MGDw3+z+qF3LFKqm7zKxhH3n7t0DvTZE2Pcsk1PSooKYWFmPDl51SzMjCcpKmTU9j1z5kxee+012tra+Oyzz/jZz37G+eefzyeffMLWrVt57733WLp0acf6TU1Nbt1/908FjY2N1NfXk5KS0uv6CQkJBAcH88UXX3Q5iLi4DhCLFy9m3bp1NDU18eGHH3LnnXeyYsUKNmzY4Nb6QQNdqS6MMeRXNpKdGe+W7WU4R8q4TrT6OhHhxZsWU9nQSlJUyIgurBouq9XKWWed1RF8rhZyaGhoxzp5eXl88cUXg74P6WB88MEH1NfXd3S7vPHGG4gI2dnZva5/1llnYbPZqKmp4Zxzzhlw++Hh4Vx44YXs3LmT3/zmN26ruzMNdKU6qW5so66lnYzEkV1U5BIdZiUuwkpB9dgIdACLRUiODh14RTfavn07d999N1deeSWTJ0+murqaxx57jLlz57J48WLS0tK46667+NWvfkVdXR0PPvggqanunYEkPDyc888/n3vuuYeSkhLuueceLrnkEmbPnt3r+jNmzOCWW25h+fLl3HvvvWRnZ9Pc3MyuXbvYt28fTz/9NO+88w7PPvssF198MRkZGRQVFfGnP/2pS5++O2mgK9VJXqWjr3ukY9A7S4+PoGCMtNC9ZcKECYwfP55f//rXFBcXExcXx5lnnsljjz1GaGgor7/+OrfffjuXXXYZaWlp/OIXv+CTTz5x62X4y5cvJzo6mhtuuIH6+nouuuginnzyyX5fs2rVKqZPn85TTz3FAw88QExMDLNnz+aGG24AHCdgRYT777+f8vJykpOTueCCC3jkkUfcVndnYrw0XDw7O9ts3rzZK/tWqi9/21bEHS9t4/1/O21EMy12dtuaHPaU1LH+7jPcsr3u9uzZw6xZszyy7UCRlZXFZZddxm9/+1tvlwL0/zsVkRxjTK/9QDrKRalO8iudFxW5uYVeVN2ksy4qj9NAV6qTvKpGxseEEmYNcts20xIiaLXZKavTWReVZ2kfulKd5Fc2jniWxe7S4x0XKBVUNbnlYiXlfocPH/Z2CW6hLXSlOsmranDbGHSXdGf3jZ4YVZ6mga6UU3ObjbLaFreOcAE67nrkyaGL3hrcoNxvJL9LDXSlnArcNMtid2HWIMbHhFJY7d4rG12sVqvbr5pU3tPU1ITVah3WazXQlXLK88AIFxdPjkUfN24cRUVFNDY2akt9DDPG0NjYSFFRUY95ZQZLT4oq5ZTnDNxMN10l2ll6QgQbD1W5fbsAMTGOOWKKi4tpa2vzyD7U6LBarYwfP77jdzpUGuhKOeVXNhAdGkx8xPA+7vYnPT6cv21ros1mxxrk/g/GrrvhqMCmXS5KOeVXNZKeEOGRCanSEiKwG8ft7ZTyFA10pZwKq5tIT/DMOPH0eNfQRQ105Tka6ErhOCFVdLSJtHj3nxAFOg4UY2nWRTX2aKArhWPa3MZWG2nxnmmhp8SGE2wRvbhIeZQGulJAobPl7KkWepBFmBgXToGHxqIrBRroSgF0XPTjqRY6OLpdtIWuPEkDXSmOtdBTPRno8REd+1HKEwYMdBFJF5H1IrJbRHaJyB29rCMi8oSIHBCR7SKywDPlKuUZhdVNxIQFExPm/jHoLukJERypb6Wxtd1j+1CBbTAt9HbgLmPMbGAxcLuIdL/J3rnANOfXzUD/921SyscUVntuhIuLqzvHU3O6KDVgoBtjSowxW5yP64A9QPe7sy4D/mIcNgBxIpLi9mqV8pDC6kaP9p+DTqOrPG9IfegikgXMB77u9lQqUNDp+0J6hj4icrOIbBaRzRUVFUOrVCkPMcaMSgv92MVFGujKMwYd6CISBbwG/NQYUzucnRljVhtjso0x2cnJycPZhFJu5+kx6C5JUSGEBlso0sv/lYcMKtBFxIojzNcYY17vZZUiIL3T92nOZUr5vKJRGLIIICKkxodroCuPGcwoFwGeAfYYY/6rj9XeAn7gHO2yGKgxxpS4sU6lPMbTFxV1lhoXridFlccMZvrcU4CVwA4R2eZcdj+QAWCM+SPwLnAecABoBK5zf6lKeYYrYD05Bt0lLT6C3cWlHt+PCkwDBrox5nOg3/lEjeM2Kbe7qyilRlNhdSPRYcHEhntuDLpLWnw4lQ2OsegRIXo7AuVeeqWoCnijMcLFxdVPr/OiK0/QQFcBzxHonu9uAUcfOqCTdCmP0EBXAc0xBt3zFxW5uD4JFGmgKw/QQFcB7WhjGw2ttlHrchkXHYo1SHSki/IIDXQV0EZj2tzOLM550XUsuvIEDXQV0I6NQR+dQAfXWHS9/F+5nwa6CmjHWuij0+Xi2Fe49qErj9BAVwGt6GgT0aGjMwbdJTUugvK6FprbbKO2TxUYNNBVQCusbhyVK0Q7c3XvlNQ0j+p+lf/TQFcBraBq9C4qcnEdQLTbRbmbBroKWK4x6OkJo9tCd11cpCdGlbtpoKuAVT3KY9BdUmLDCLKIDl1UbqeBrgKWq4WcPsp96MFBFibEhOnFRcrtNNBVwCqoGv0hiy6pOnRReYAGugpYHRcVjXIfOkCaXlykPEADXQWsgupGYsOtxISN3hh0l7T4cEprm2mz2Ud938p/aaCrgFVY3TTqI1xcUuPDsRso1bHoyo000FXAKqhqJC1u9PvP4Vi/vZ4YVe6kga4CkmMMuhdb6DoWXXmABroKSEfqW2lpt3tlhAtASlwYItpCV+6lga4CUoFrDLqXWuihwUGMj9ax6Mq9NNBVQPLGtLndpSeEdxxYlHIHDXQVkAqqRv/GFt2lx0dQWKWBrtxHA10FpMLqJhIjQ4gICfZaDWkJEZTUNtParmPRlXsMGOgi8qyIlIvIzj6eP0NEakRkm/PrAfeXqZR7FVY3kpbgve4WcMwhYwwU6yRdyk0G00J/Dlg6wDqfGWPmOb8eHnlZSnlWYXWTV7tbANKdBxTtR1fuMmCgG2M+BapGoRalRoXdbiiqbiLdiydEoVOgV2kLXbmHu/rQTxKRXBF5T0SO62slEblZRDaLyOaKigo37VqpoSmva6HVZvd6C31CTBjBFtEWunIbdwT6FiDTGDMX+G/gzb5WNMasNsZkG2Oyk5OT3bBrpYbu2Bh077bQgyzCxLjwjhE3So3UiAPdGFNrjKl3Pn4XsIpI0ogrU8pDOqbN9XILHVxj0bXLRbnHiANdRCaIiDgfL3Jus3Kk21XKU1x91q75VLxJx6IrdxpwEK6IvAicASSJSCHwIGAFMMb8EbgMuFVE2oEmYLkxxnisYqVGqLC6kXHRoYRZg7xdCukJEVQ2tNLQ0k5kqPfGxCv/MOBfkDHmqgGe/wPwB7dVpJSHFVQ1eb3/3MXV7VN0tInp46O9XI0a6/RKURVwCqobfaL/HDoPXdRuFzVyGugqoLS02yg+2kRWYqS3SwHoGAuvga7cQQNdBZTC6ibsBrKSfKPLJSkqhHBrkI50UW6hga4CSl5lAwCZPtJCFxHS4nUsunIPDXQVUA4fcQSnr3S5gKMfXVvoyh000FVAOVzZQHRYMPERVm+X0iEtPpzCqkZ0tK8aKQ10FVAOVzaSlRiJ81o4n5AeH0FdSzs1TW3eLkWNcRroKqDkVTaQmegbJ0RdXPc11VkX1UhpoKuA0WazU1jtO0MWXVz3NdVZF9VIaaCrgFFU3YTNbnywha5j0ZV7aKCrgHHYOWQxK8m3Wuix4VZiw63ka6CrEdJAVwEjr9IRmL7WQgfISozoOOAoNVwa6CpgHK5sIDIkiOSoUG+X0sOkpMiOMfJKDZcGugoYh480kOljQxZdspIiKa5pornN5u1S1Bimga4CRl5lo8/M4dLdpKRIjEH70dWIaKCrgNBus1NQ3egzc7h05xpKeeiI9qOr4dNAVwGhpKaZNpshywdPiMKxkTca6GokNNBVQDjsY7MsdhcbbiUxMoTDGuhqBDTQVUA4XOl7syx2l5UUqS10NSIa6Cog5B1pIMxqYVy07w1ZdMlKjNSx6GpENNBVQDhc2UhmQiQWi+8NWXSZlBRBWW0Lja3t3i5FjVEa6Cog+OIsi925TozqBUZquDTQld9rs9nJq2xkUrLv9p+DYyw66EgXNXwDBrqIPCsi5SKys4/nRUSeEJEDIrJdRBa4v0ylhi+vsoFWm50Z46O9XUq/XCdstR9dDddgWujPAUv7ef5cYJrz62bgyZGXpZT77C2tB2C6jwd6ZGgw46JDtYWuhi14oBWMMZ+KSFY/qywD/mIcN0TcICJxIpJijClxU43Kh9nthvd3l/LnLw6TEBnCTadNZkFGvLfL6mJvWR0WganjorxdyoCykiJ1LLoatgEDfRBSgYJO3xc6l/UIdBG5GUcrnoyMDDfsWnnTup2l/N/397K/vJ70hHC+Ka3jvZ2lfGdSAvedN4t56XHeLhGA/WV1ZCZGEmYN8nYpA5qUGMlH35R5uww1Ro3qSVFjzGpjTLYxJjs5OXk0d63c7PP9R7h1TQ4i8Pvl81h/1xl8+fOz+D/nz+JwZQPXPLuRoqO+cY/MvWV1TB/v+61zgEnJkRypb6W2WW8YrYbOHYFeBKR3+j7NuUz5qYq6Fv7tlW1MSY7ib7efyrJ5qQQHWYgMDebGJZN5+eaTsNkNP/7fLbTZ7F6ttbnNxuEjDT5/QtSl48SodruoYXBHoL8F/MA52mUxUKP95/7LbjfctTaXmqY2/rBiPuEhPbsxspIi+c33TmBL/lF++4+9XqjymIMVDdgNTBsjga5DF9VIDNiHLiIvAmcASSJSCDwIWAGMMX8E3gXOAw4AjcB1nipWed9Tnx3k030V/MfFxzNzQkyf6104dyIbDlbyp08P8p3JCZw1c/woVnnMvrI6AGZMGBuB7rr4SS8uUsMxmFEuVw3wvAFud1tFymcVVjfy2/f3svS4CVz9nYFPav/7BbPJyavmF2/s5JN7kggNHv2TknvL6rAGiU9PytVZmDWI1LhwDh6p93YpagzSK0XVoP3pnwcBeODC2YO6jVuYNYj7z5tFSU0zazcXerq8Xu0vq2NyUhQhwWPnT33a+Cj2ltZ5uww1Bo2dv3LlVeW1zby8uYBLF6QxMS580K9bMi2J+RlxPPnJt7S2j/4J0r1ldUwbIyNcXGZMiObbinqvn1BWY48GuhqUpz47SLvNzq1nTBnS60SEO86eRtHRJl7bMrqt9IaWdgqqmsbMCBeXmROiabMZPTGqhkwDXQ2oqqGVNV/nc9HcicO648/p05OZmx7HqvUHRrXVeaDcecn/GDkh6jJjvONk8zfa7aKGSANdDejPXxyisdXG7WdOHdbrRYSfnj2NwuomXh/FVvpe1wiXMdZCnzIukiCLsLe01tulqDFGA131q76lnee+PMzS4yaMaCz3GTOSmZMWy5OffIvdbtxYYd/2ldYRGmwhPcG350HvLjQ4iMlJkXpiVA2ZBrrq11vbiqlrbuem0yaPaDsiwg2nTuJwZSOfHTjipur6t6+8nmnjowjy4bsU9WXGhGjtclFDpoGu+vXSpnxmjI9mQcbIJ9paevwEkqJCeGFDnhsqG9i+0jqfnzK3LzMnRFNY3UR9i96OTg2eBrrq067iGrYX1rB8Ufqgxp0PJDQ4iCuy0/loT5nHJ+4qr2umtLaZ2Sl9X83qy2Y4r8LVbhc1FBroqk8vbSwgJNjCJfNT3bbNFd/JwAAvbcx32zZ7syWvGoAFmb41N/tgzXSOzNFAV0Ohga561dRq482tRZx/QgpxESFu225afARnzxzHixsLPHqhUU5eNSHBFo6fGOuxfXhSalw4kSFBOtJFDYkGuurV29uLqWtpZ/mJ6QOvPERXL87kSH0L/9hV6vZtu+TkVTM3LXZMXfLfmcUiTNcTo2qIxuZfu/K4lzYVMDk5kkWTEty+7dOnJZOeEM5fPXsuTzsAABhLSURBVHRytLnNxs6i2jHb3eIyc0I0e8vqcMx/p9TANNBVDwfK68jJq2b5ie45GdqdxSJc/Z1MNh6q4kC5+1ugO4tqaLXZWehj9zYdqhnjozna2EZ5XYu3S1FjhAa66uHVnCKCLMIl89M8to/LFqZhDRJe3Fgw8MpDlOM8IbpwjLfQXVMWaLeLGiwNdNWFzW54c2sRp09PJjk61GP7SYoK5bvHTeC1LYU0t9ncuu3NedVMSookMcpz9Y8G1w1E9mmgq0HSQFddfPntEUprm7l0geda5y4rFmVwtLGNdTvdd3LUGMOWvGoWjPHuFoCEyBCSo0PZoyNd1CBpoKsuXsspJCYsmLNnjfP4vk6anEhWYgT/+7X7xqTnVTZS2dA65rtbXI6fGMOOwhpvl6HGCA101aG+pZ11u0q5YO5Ewqyev12cxSIsX5TBxsPuOzm62dl/np3lH4G+MDOe/eX11DS2ebsUNQZooKsO7+4oobnNzqUL3Hdl6EDcfXI0J6+a6LBgpiaPrbsU9cU19HJLQbWXK1FjgQa66vD6lkKyEiNGtf/ZdXL01ZxCmlpHfnLU1X9uGYMzLPZmblocQRbpmMpAqf5ooCsACqoa2XCwiu8tSPPI2PP+XHNSFjVNbby+dWQ3v6isb2FfeZ3f9J8DRIYGMyslumMoplL90UBXALy2pRAR3DoR12CdmBXPCamxPPP5oRHd/OLdnaUYA+fMHu/G6rxvYUY82wqO0q43jVYDGFSgi8hSEdkrIgdE5Oe9PH+tiFSIyDbn143uL1V5it1ueDWnkJOnJHrl7j4iwo1LJnGwooFP9pUPezt/31bMtHFRHTMV+osFmfE0ttr0AiM1oAEDXUSCgFXAucBs4CoRmd3Lqi8bY+Y5v552c53KgzYcqqSwuonLF7p/Iq7BOu+EFCbEhPHM54eG9fqSmiY2Hq7iwrkTR73LyNNcXUhb87XbRfVvMC30RcABY8xBY0wr8BKwzLNlqdH06uZCokOD+dfjJnitBmuQhWtOzuKLA5XsLh76hTRv55YAcNHcie4uzetS48IZHxOq/ehqQIMJ9FSg85iyQuey7i4Vke0i8qqI9NrUE5GbRWSziGyuqKgYRrnK3eqa23h3ZwkXzptIeIjnx573Z8WiDMKtQcNqpb+VW8yctFiykiI9UJl3iQgLM+PJ0Ra6GoC7Tor+HcgyxswBPgCe720lY8xqY0y2MSY7OTnZTbtWI/HOdsfY88sXev5S/4HERli5IjuNt3KLKKhqHPTrDh1pYEdRjV+2zl0WZMRTUNVEeW2zt0tRPmwwgV4EdG5xpzmXdTDGVBpjXHN8Pg0sdE95ytNe2VzA1HFRzEsf+U2g3eGWM6ZgDbLwq7d3D/o1f88tRgTOn5Piwcq8q+MCI22lq34ED2KdTcA0EZmEI8iXAys6ryAiKcaYEue3FwF73Fql8ohvK+rZkn+U+8+b6TMnElNiw7n9zKn85z/28um+Ck6b3v8nOWMMb+UWc2JWAimx4YPbSfHWoRc2cf7QX+NGx02MISTYQk5eNUuP998DlxqZAVvoxph24EfAP3AE9SvGmF0i8rCIXORc7ScisktEcoGfANd6qmDlPms25GMN8uy858Nx45JJZCZG8Mu/7xrwvqMf7innQHk9y+b5b3cLQGhwEHNSY9l4qMrbpSgfNqg+dGPMu8aY6caYKcaYXzuXPWCMecv5+D5jzHHGmLnGmDONMd94smg1cg0t7azNKeC8E1I8Ou/5cIQGB/HABbM5WNHA818e7nO9yvoW7nt9O7NSYrw65HK0nDEjmdzCGu1HV33SK0UD1JvbiqhrbucHJ2V5u5RenT1rPGfOSOb3H+3vdfpYYwy/eGMntU3t/O7KuWP2ZtBDcc5sx7DSD/cM/+Ir5d/8/3+B6sEYw/NfHub41BgWZPjGydDePLzseGLDrVy5+is+/qasy3Nvbiti3a5S7vzu9I47+/i76eOjyEiI4P3d7rshiPIvgzkpqvzMhoNV7Cur5/HL5vjMydDepCdE8MZtJ3P985u48fnN3Lt0JtFhwewqruXv24rJzoznpiWTvV3mqBERzpk9nr9+lUd9SztRofrfV3WlfxEB6C9fHSY+wjomxm2Piwnj5ZtP4scvbuXR9xynZmLCgpmXEccjl5xA0GhNkzvUkTEeGhVzzuzxPPP5IT7dV8F5J+hoF9WVBnqAKT7axPu7y7hpyeRRuSuRO0SGBrN65UK2FhxlQkwYafHhPv3JwpOyM+OJi7Dywe4yDXTVgwZ6gHFdVn/1dzK8XMnQBAdZODErwdtleF1wkIWzZo7joz3ltNnsWIP0NJg6Rv8aAkh5bTMvbMjjkvmpXpkmV7nHd2ePp6apjU2HdUy66koDPYA8+c9vabcbfnzWVG+XokZgybRkQoItfLC7bOCVVUDRQA8QpTXNrPk6n0sXpJKZ6H8zEgaSyNBglkxN4r0dpbTpXYxUJxroAeLJTw5gtxt+fNY0b5ei3GDFdzIorW3m3R0lA6+sAoYGegAoqWnixY0FXJ6dpn3nfuLMGeOYnBTJM58fwpjh34dV+RcN9ADw+Lq9GAy3n6l95/7CYhGuP3US2wtrdMIu1UED3c99srecN7YWcesZU0mL19a5P7l0QRrxEVaeHuZ9WJX/0UD3Yw0t7fzijZ1MSY7k9jOneLsc5WbhIUF8f3EmH+4p49CRBm+Xo3yAXljkx377/l6Ka5pY+8OTCA0eG1eFBrRh3Hhj5Umz+NM/D/LnLw7x8LLjPVCUGku0he6ncvKqee7Lw6xcnEm2XmHpt8ZFh3Hx/Im8tKmA/WV13i5HeZkGuh8qPtrEbWtymBgbzj3/OsPb5SgPu+dfZxIVGsydr+TquPQAp10ufqauuY3rn9tEY4uNtbcuIjrM6u2SAtNw7ls6TMnRofzHxcdz25ot/M/6b7njX/Rag0ClLXQ/0mazc9uaLRwor+fJ7y8MmBs/KDjvhBSWzZvIf3/c+x2eVGDQFrqfaGxt565Xcvls/xEev3QOp05L8nZJo2MUW8K+7uGLjmfDwUp++vJWXv7hSSRF+da9YpXnaaD7gYKqRm7+aw7flNbyf86fxRUn+sANk4cTtAPcFMJuN1Q2tJIUFeLx+dDtdqhpbiMu3MpYmXo9NsLK766cx/XPbeLyP37FX65fpFcGBxgN9DHMGMP6veXc9UouNrvhz9eeyBkzxg38wiGEbUewTTnRqzeVsNsNVz21gZy8ahZmxvPiTYuxOO9W5K7wdW0nJiyYX7y5kz0ltcxKieE/lh1PXUv7gNsfSh2eOmAsnpTIH1Ys4M6Xt3Hpk1/y/PWLmJWiXW+BQgN9jNqaX83j6/by1cFKpo2LYvUPspmU5L5ZFO12qG5q5fF/7OWbklrCMmxdQnQ02e2GfWV1bM6rxmY35ORVU9nQSnJ0KHY73P/mDvaU1DJ1XBSPXzpn0Lel6xyqxhzbzrRxUewvr8dmN+wpqeXnr29nf3k9s1JieOTiE7B0O/PUbjPkVzXwp08P8k1p3YAHgc4197XNvursEf6dDs52O9z3xna+KanlpKR4tjUkctn/fM5t2VFcNy+SCKsFux0qY2ePyqccNfoGFegishT4PRAEPG2MebTb86HAX4CFQCVwpTHmsHtLVTWNbby/u5S3cov5bP8REiND+OWFs1nxnUxCgkd2frvXcCuuweac96k5P4eag0HER3hm1IzdDpV1LT2CxtUy33y4iojQYBpbbSzMjCcpKgRw1LynpBab3bC3tI6fvbadxy+d2yMgXQcoAeIjQjDGEX67i2uZnBzJv58/q2M7+8rqmJocxYHyeiYnRbCvtA47sKeklqNNrcSFh3S8V+3tdq5+9msaW20d+9pVVMMdL2+hoKqJaeOjuxxk7HY4XFnP7qKajm3WNLd1eV/7OtAMFP7VTa3sKq4FoKiimmeumMzvNzXxn1/V8VxuAz86MZpt+w7x97KSHp9ylH8YMNBFJAhYBZwDFAKbROQtY8zuTqvdAFQbY6aKyHLgMeBKTxQcKJpabRTXNPFNSR3bi46yo7CGTYeraLMZ0uLDufOc6Vx/6iS33Pm9e4vxnn+d4Qg3Z5gHWYRZKTHEhVt7vK5z8AzUhdDX+sdC64seoVXT2EZzfg6zMAS1wv9bPp+sxDCkZBsAceFWpo6LYm+p46Ka/WV1vQbkfW9s7wi741JiuPG0yR3ff1vRwA//msOM8VHsLqnDbuBART124GBFA66R3Ta74bH39gDCN6V1zEyJprnV1iXMAQyQV9kE0OUgA3Df67nsKjl2AdDM8VEYYzAGRHr+Lu46Z/qx8C+u6fGzddb9bZ8YFcRTFySwubiFR7+o48F/1gDxgGHT4Sre3l7Cwqx4JsaGaWvdTwwmDRYBB4wxBwFE5CVgGdA50JcBv3Q+fhX4g4iI8cC8ngfK633mTi2GYz9e95/U9aMbA3YDdmMwxmAzhna7od1maG2309Juo6XdTn1zO7XNbdQ2tVNW18zRxraObVmDhBkTorn25CwumDOROWmxw/4PaLcbahq7Bm/nVu6ekloEmJUSw57iGqaOi+L+82aREBnasX737hjH8EjT0d3QWyvSEao7+Ka0lpkTogE61u84gDj33zm04sKtjlqcAZeVGNnlgCECj186h5+9tp39ZXXMmhjbceBxHUCMMXxTUtvxml0ltfzby9u61NdiM9R1es/tzt+nrdvvtXMYuw4ILqEWaLPTcQBwcR1k7MZ0eX2QCO0GrntuU0c3TUF1Y8cno91FNfzm3T0d2wsJDiI6tO8pHOIjQjhuYqzjPU6JIT7C8Skme2Ioay8LYVtpG79Yl8/uugjsBn7ykqO7JtwaRHJ0KAmRISRGhhAZGkyY1UK4NQhrkIWgIMFqsWCxCBZx1O36HXT/O+z+Zyk9DjMKYF56HCdNSXT7dgcT6KlAQafvC4Hv9LWOMaZdRGqAROBI55VE5GbgZoCMjOHdpHhvaR2PrftmWK/1BcEWIcgiBFuEUGsQocEWQoItRIUGExNmJTMxghMnxZMSG87EuDCmJkczfUKUW+Zi6evEYpwxhGXYHMsz4omfuoj/mHQil6/+ijcLayhcb3jxpnmIRbp0gdhMIpDIzhIBY7CZJHYVC3fFziY5uuuQuYqaZl4uKgYSyS0Ci4Dduf6d4+cSlmHv2H/clBM7kkGAR26b1+/oFrEb7rl2NiKOi2xEpMvPuiAzHmvafHLz+x+fvfMoRIRYaGw9FsmOOo+tExoELbZeXgzMGR/FtxWNNLR1jfRFmQnETTmRiroWdppj/yWOnxjDm8W12EwyO4tg79+b2F5YS7h1GvXOnWwvA3Ce6G6FK/7ezKu3nNxrV4kAv7m99/dKgPmp8Pb8RVQ2tAKGA+UNHKio51BFA0fqW6hubKWkppmmNhtNrTaa2my02ezOBoi9y/ugRuaW06d4LdDdxhizGlgNkJ2dPaw/j6XHT+CbXy11a13u0qXlyLFWjEUcLRtvf6ytbGglJ6+a9m4nFkWEF29a3CUIqpta2VFY0+MkpGsbnbtjFmbGgzFsyT/apX/bxW43VDW0dFk2Jy2WnUW1LMyMJzk6tMf+O7NYpMcBovO2ux+kRHr+rHNSY7GIY78hQRY2Hq7udXudw1xwhHm41UKTM6T7CnOA7cX1vS5/4qp5iDh+hkVZCWzOq2Juehyv/nAxK57e6KgvLZZc5/td389OtuUf7fhddH8fXO9fX+8VdH0vk6PDhhwqdrvjUyY4Pn329ylV9W2wJ+6HajCBXgR0Htic5lzW2zqFIhIMxOI4Oep2QRYhyKIzBw5HUlQICzPjO8Kvc/B2D82+1u28fEFGHH9YsYDk6FCMoUsguwImIcLKiqe/ZnNeNVGhQTS22MjOcgRvVWNbx/qu1vVQ9XWQ6lznCWmxbC84it3AzqJavvj5WRypb+H8Jz4/9vPjCPttna6ydOVTU5udiBALTa12+suseWkxbCus7bHc4vq0IcJLN3c9cLkOZPHhwVy+egPbnaHe2fPXLuS2F7fR0GIjMiyYhIju5zL6HtLpbhaLYNFuFJ81mEDfBEwTkUk4gns5sKLbOm8B1wBfAZcBH3ui/1yNTG8t8aGu2/fyY4HcOWBcYWozjmB8944lzJgQ3dFiHam+Djyd60yMtHLVU193rDMuOpSkyBAiQyw0tNqJsArr7z6T5OhQvvc/X/Qays1tdt79yak8+NauLq17i8BfbljEtOQokqNDWb76a3Lyqwm3BtHY2k52VkKXn7P7gdNiERIjQ7jqqQ1sLzjKnLRYgi3CpryjAESGWJiZEkuz8xNCY6uNqsa2Ltvo66CmAs+Age7sE/8R8A8cwxafNcbsEpGHgc3GmLeAZ4C/isgBoApH6Csf1F/3xWDXHWgbnQNme2ENc9Pj2F5Yw8LM+I4wd5f+DlKd6+y+zpGGFhqc3SuNbY4hJhaLhdduPYVL//ilo96MuI6fITsrgZkpMbx080mU1zZz65otbC+qISIkiGue2Uh2VgIv3rS4owWeEGHt8gmkP527sXY4P0EYu6G6sbXj/cru45MV9P/JSwUW8VZDOjs722zevNkr+1aeZYxh+epjXQD/e+N3Bh1uo6W8tplFj3zU8f3G+89mXEwY0LU/untXkovrYqfzn/gMm3Gc7P7qvrOH1TLu/n69dPPiHu/TQNMejOa0CMq7RCTHGJPd23N6pahyu95azb7WBeA6QZmTX91xYtalc8u+r759i8UxlDQ7K2HELePBdIUN9KloKJ+8lP/SFroKWO5o1WrLWI02baEr1Qt3tGq1Zax8id7gQiml/IQGulJK+QkNdKWU8hMa6Eop5Sc00JVSyk9ooCullJ/w2jh0EakA8ryy894l0W26Xx/j6/WB79fo6/WB1ugOvl4fjKzGTGNMcm9PeC3QfY2IbO5rsL4v8PX6wPdr9PX6QGt0B1+vDzxXo3a5KKWUn9BAV0opP6GBfsxqbxcwAF+vD3y/Rl+vD7RGd/D1+sBDNWofulJK+QltoSullJ/QQFdKKT+hgd4LEblLRIyIJHm7ls5E5D9F5BsR2S4ib4hInLdrAhCRpSKyV0QOiMjPvV1PdyKSLiLrRWS3iOwSkTu8XVNvRCRIRLaKyNverqU3IhInIq86/wb3iMhJ3q6pOxH5N+fveKeIvCgiYT5Q07MiUi4iOzstSxCRD0Rkv/PfeHfsSwO9GxFJB74L5Hu7ll58ABxvjJkD7APu83I9iEgQsAo4F5gNXCUis71bVQ/twF3GmNnAYuB2H6wR4A5gj7eL6MfvgXXGmJnAXHysVhFJBX4CZBtjjsdxD2RfuL/xc8DSbst+DnxkjJkGfOT8fsQ00Hv6HXAv4HNni40x7xtj2p3fbgDSvFmP0yLggDHmoDGmFXgJWOblmrowxpQYY7Y4H9fhCKJU71bVlYikAecDT3u7lt6ISCxwGo4bwmOMaTXGHPVuVb0KBsJFJBiIAIq9XA/GmE+Bqm6LlwHPOx8/D1zsjn1poHciIsuAImNMrrdrGYTrgfe8XQSOYCzo9H0hPhaWnYlIFjAf+Nq7lfTw/3A0JOzeLqQPk4AK4M/ObqGnRSTS20V1ZowpAn6L49N1CVBjjHnfu1X1abwxpsT5uBQY746NBlygi8iHzv617l/LgPuBB3y4Ptc6v8DRjbDGe5WOPSISBbwG/NQYU+vtelxE5AKg3BiT4+1a+hEMLACeNMbMBxpwUzeBuzj7oZfhOPhMBCJF5PverWpgxjF23C09AgF3T1FjzL/0tlxETsDxh5DrvNlvGrBFRBYZY0q9XZ+LiFwLXACcbXzjIoIiIL3T92nOZT5FRKw4wnyNMeZ1b9fTzSnARSJyHhAGxIjIC8YYXwqjQqDQGOP6ZPMqPhbowL8Ah4wxFQAi8jpwMvCCV6vqXZmIpBhjSkQkBSh3x0YDroXeF2PMDmPMOGNMljEmC8cf8ILRDPOBiMhSHB/LLzLGNHq7HqdNwDQRmSQiIThOQr3l5Zq6EMcR+hlgjzHmv7xdT3fGmPuMMWnOv7vlwMc+FuY4/x8UiMgM56Kzgd1eLKk3+cBiEYlw/s7PxsdO3HbyFnCN8/E1wN/csdGAa6GPcX8AQoEPnJ8iNhhjbvFmQcaYdhH5EfAPHKMKnjXG7PJmTb04BVgJ7BCRbc5l9xtj3vViTWPRj4E1zgP3QeA6L9fThTHmaxF5FdiCo0tyKz4wDYCIvAicASSJSCHwIPAo8IqI3IBjGvEr3LIv3/jUrpRSaqS0y0UppfyEBrpSSvkJDXSllPITGuhKKeUnNNCVUspPaKArpZSf0EBXSik/8f8BDOx6KALVEqgAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "samples = metropolis(func, Gaussian(mu=np.zeros(1), var=np.ones(1)), n=100, downsample=10)\n", "plt.plot(x, func(x), label=r\"$\\tilde{p}(z)$\")\n", "plt.hist(samples, density=True, alpha=0.2)\n", "plt.scatter(samples, np.random.normal(scale=.03, size=(100, 1)), s=5, label=\"samples\")\n", "plt.legend(fontsize=15)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "### 11.2.2 The Metropolis-Hastings algorithm" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAgAElEQVR4nO3dd3ib1dn48e8tW947XolHnJANZBpIWC0BWsIKFCiBkjJLKbTQt4wW+hYK/ZXSvh1vR0rfMAotFGgYLWWkBEiaUggkzt7be8VxvJes8/tDkmM7si3bkiVL9+e6fEV6dPQ8t6z41tF57uccMcaglFJq9LP4OwCllFLeoQldKaWChCZ0pZQKEprQlVIqSGhCV0qpIBHurwOnpqaavLw8fx1eKaVGpYKCgiPGmDR3j/ktoefl5bFhwwZ/HV4ppUYlESns6zEdclFKqSChCV0ppYKEJnSllAoSmtCVUipIaEJXSqkgMWBCF5EoEflMRLaIyA4RedRNm0gReUVE9ovIpyKS54tglVJK9c2TssU2YKExplFErMBHIvKuMWZdtza3ArXGmEkisgT4KXCtD+JVSrlRX19PVVUVHR0d/g5FDYPVaiU9PZ2EhIQhPX/AhG4c8+s2uo7n/Ok95+5i4IfO268CvxMRMTo3rwpi5XUtfLy/hivnZGGxiN/iqK+vp7KykqysLKKjoxHxXyxq6IwxtLS0UFpaCjCkpO7RhUUiEgYUAJOAZcaYT3s1yQKKnUHZRKQOGAMc6bWf24HbAXJzcwcdrFKBYntpHbc8t56qhjaiI8K4+NSxfoulqqqKrKwsYmJi/BaDGj4RISYmhqysLMrKyoaU0D06KWqM6TTGzAaygdNF5JRBH8mxn+XGmHxjTH5amtsrV5UKeO/vrOSaP3yCNcxCbkoMv35/H3a7/76MdnR0EB0d7bfjK++Kjo4e8tDZoKpcjDHHgNXARb0eKgVyAEQkHEgEaoYUkVIBbPXuKm7/8wYmZ8Txxp1ncu8XprCnsoGVOyr8GpcOswSP4byXnlS5pIlIkvN2NHAhsLtXszeBG523rwY+1PFzFYxWFBSTGhfJy7fPJz0hiktnjmNiWqzfe+lKgWc99LHAahHZCqwHVhlj3hKRx0TkcmebZ4AxIrIf+A7wPd+Eq5T/2DrtfLTvCJ+bkkZMhOP0U5hFuOf8yQHRS1fKkyqXrcAcN9sf7na7FbjGu6EpFVi2lNRR32rj3Ck9z/9cOnMcv/lgH79+fx8XnZzp14oXFdr0SlGlPLR2bzUicPak1B7bwyzCnZ+fxJ7KBjYV1/opuuBlt9v9HcKooQldKQ+t3VfNzOwkkmMjTnjsnCmOJL+p6NhIhxV0vvrVrzJ58mTuueceHn74YW699dYejxtjmD17Ns8//7zH+/zmN795wn6CkSZ0pTxQ19zBluJjfG5yqtvH0+OjyEqKZlOxJvThqK2tZcOGDWzbto3S0lKWLVvGLbfc0qPNX//6V44ePcr111/v8X7vu+8+XnzxRfbv3+/tkAOKJnSlPPDR/iPYDSeMn3c3JzeJzdpDH5bk5GR27txJVFQUr776KjU1NZxzzjk92vzmN79h6dKlWK1Wj/ebl5fH2WefzZNPPuntkAOKJnSlPLB2bzXxUeHMzknqs82c3GRKj7VQWd86gpEFl5SUFH7wgx/wla98hTFjxpCSksLXv/512traANi/fz8ff/wxV199dY/nvfLKK4iI25/6+noArrrqKl588cWgHpPXhK7UAIwxrN1XzdmTUgkP6/tPxpXsdRx9aIqKiqitreVXv/oVKSkprFixgm9/+9ssX76cZcuWAfDBBx8QGxvLrFmzejz33HPP5ZNPPun6efvtt0lJSeHiiy/uuoT+zDPPpLKykm3bto34axspflskWqnRYl9VI+V1rdx9fv/TVZw8LgFrmLC5+BgXnZI5QtEFjy1btgCOE5hPPPEEAAsXLmTt2rWsXLmS73znOxQUFDB9+nQslp4frGPHjmXsWMd8OvX19Zx//vmccsoprFixoqvNySefTFhYGJ999tkJHwjBQhO6UgNYu7ca6H/8HCDKGsaMcYlsKvJ/6eKj/9jBzrJ6vxx7xrgEHrns5EE/b8uWLcTExHD//ff32D5lyhQ+++wzACoqKkhNdX9iGqCpqYmLL74Yi8XCW2+91WPCsvDwcJKSkqioCN4LwHTIRakBbCutY1yio4plIHNykthaUoetM3jHaX1ly5YtnHvuuYwZM6bH9oqKCsaNGwdAa2srkZGRbp/f2trK5ZdfTlNTEytXriQ+Pv6ENpGRkbS2Bu85Du2hKzWAPRUNTM08MTm4Myc3iec+PsyeygZOHpfo48j6NpQesr9t2bLlhIqWhoYGPvzwQx591LFQWkpKitsedkdHB1dddRXl5eX861//Ijk52e0xjh07RkpKiveDDxDaQ1eqH7ZOOwerm5iS4WFCz3Ekks1ajz4oTU1NHDhwgAMHDvTY/pvf/AaAG290zP03depUDh061KNNZ2cn1113HXv27OH999+nr6m5q6uraW5uZsqUKT54BYFBE7pS/Thc00x7p93jhJ6TEs2Y2AitdBmkbdu2YbfbOXDgAI899hirV6/mwQcf5NFHH2X58uUkJTkqiM466yyKioqorq7ueu43vvEN3nnnHR5//HGKiopYt24d69atO6GaZcOGDYgIZ5555oi+tpGkQy5K9WNvZQOAx0MuIsLsnKSAODE6mmzZsoWEhATefvttbr75Zn784x8zdepUVqxYweLFi7vaff7znyclJYWVK1eydOlSjDG8/PLLtLS0cO21PZcxvuyyy3jzzTe77q9cuZLPfe5zJ4zRBxNN6Er1Y29lAyJwUlqcx8+Zk5vEB7urqGvuIDHG86sZQ9mWLVs49dRTmTlzJgUFBX22i4iI4IYbbuDll19m6dKlPS4c6k9nZyevvfZaVzlksNIhF6X6sbeygfEpMURHhHn8nNnOcfTtZXW+CivobNmyxePa8Pvvv5/Vq1ezd+9ej/e/YsUKoqOjWbJkyVBDHBU0oSvVjz0VDR6Pn7tMznD05g9UN/oipKBjjGHbtm3MnDnTo/bZ2dk8++yzlJeXD+oYzzzzDOHhwT0oEdyvTqlhaLN1crimmUWnjB3U89LjI4mNCONgdZOPIgsung6bdDfYnvZ11103qPajlfbQlerDweomOu2GKR6eEHURESakxWoPXY04TehK9aGrwmWQQy4AE1PjtIeuRpwmdKX6sLeygXCLMCE1dtDPnZgWS1ldC60dnT6ITCn3NKEr1Yc9FY1MSI0lInzwfyYT0+IwBg7XaC9djRxN6Er1YW9lw6DHz10mOnv1OuyiRpImdKXcaG63UVzbzJT0oSX0CV0JXU+MqpEzYEIXkRwRWS0iO0Vkh4jc46bN50WkTkQ2O38e9k24So2M/VWNGANTMz2/QrS72MhwMhOitIeuRpQndeg24F5jzEYRiQcKRGSVMWZnr3b/NsZc6v0QlRp5eyocFS6Dvaiou4lpsRw4ogldjZwBe+jGmHJjzEbn7QZgF5Dl68CU8qd9VY1EhFsYP2bwFS4uE9NiOVjdiDHGi5Ep1bdBjaGLSB4wB/jUzcMLRGSLiLwrIm5n1xeR20Vkg4hs6D79pVKB5mB1IxPGxBJmkSHvY2JqHA2tNo40tnsxMuUva9asQUTYvn27v0Ppk8cJXUTigNeAbxtjel+nuxEYb4yZBfwW+Ju7fRhjlhtj8o0x+X1NQq9UICg62kzumJiBG/ZjYpqjd39Ih13UCPEooYuIFUcyf9EY83rvx40x9caYRuftdwCriPS9kqtSAcwY40joKcNL6K4pd7XSRY0UT6pcBHgG2GWM+WUfbTKd7RCR0537rfFmoEqNlOqGNlo77IwfZg99XFI0EeEWDmoP3SM7duzgoosuIiUlhdjYWKZPn86yZcsAePvtt7nwwgtJT08nISGB+fPn89577/V4/g9/+ENSU1P59NNPyc/PJzo6mrPPPptDhw5RVVXFFVdcQVxcHNOnT+fDDz/s8dy8vDzuu+8+fvSjH5GZmUlcXBxf+cpXqKvrfwpku93OE088waRJk4iMjGTKlCk8//zzPdp89NFHnHPOOSQkJJCQkMDs2bNZsWKFF35jJ/KkyuUsYCmwTUQ2O7c9BOQCGGP+AFwNfENEbEALsMTomSA1ShUebQYYdg89zCLkjYnRHrqHLrvsMqZPn84LL7xAZGQke/bs6ZqF8dChQ1x22WXcd999WCwW3n33XRYtWsTatWs566yzuvbR3NzM7bffzgMPPEBsbCx33303S5cuJTIykkWLFnHnnXfys5/9jGuuuYbi4mJiYo6/xy+99BKTJk3iqaeeory8nAceeIDbbrut3+T7rW99i+eff56HH36YuXPnsmrVKm655RbGjBnDpZdeSn19PZdeeimLFy/m4Ycf7poq+NgxHy1RaIzxy8+8efOMUoFoxYZiM/67b5kDVQ3D3tfX/7TBnPc/q4cfVD927tzp1f11dtpNVX2rsdvtXt1vf6qrqw1gtm7dOmDbzs5O09HRYb7whS+Ym2++uWv7I488YgCzZs2arm3Lli0zgHn00Ue7tu3YscMA5p133unaNn78eJOcnGwaGo6/5y+88IIRka7f7+rVqw1gtm3bZowxZt++fUZEzHPPPdcjvqVLl5r8/HxjjDHr1683gKmvrx/Mr6Pf9xTYYPrIq3qlqFK9FB1txiKQnTy8Hjo4TowWHW2mo9Puhch8z243XPfUOhb85AOWLF+H3T4yX7RTUlLIycnhjjvu4JVXXqGqqqrH4yUlJdx4441kZWURHh6O1WrlvffeO2HVooiICM4555yu+5MmTQJg4cKFJ2wrLS3t8dwLL7yQuLjjF5JdeeWVGGNYv36925g/+OADLBYLV155JTabrevn/PPPZ/PmzXR2dnLSSScRFxfH9ddfz9///nff9cydNKEr1UtRTRNjE6OHNClXbxPT4rDZHSdZR4OapnYKCmux2Q0FhbXUNI1MyaXFYuG9994jMzOTW265hczMTM455xw2bdqE3W7n8ssv5+OPP+axxx5j9erVrF+/nkWLFtHa2tpjP/Hx8Vgsx9+3iIgIAJKSkk7Y1vu56enpPe7HxMQQFxfX58pIR44cobOzk8TERKxWa9fPTTfdhM1mo7y8nOTkZFatWkVHRwdf/vKXSUtL45JLLuHgwYND/2X1Q1csUqqXQi9UuLi4ShcPVjcNaqFpf0mNi2De+GQKCmuZNz6Z1LiIETv2tGnTeO211+jo6ODf//433/3ud7nkkktYs2YNmzZt4t133+Wiiy7qat/S0uLV4/f+VtDc3ExjYyNjx7pfsSolJYXw8HD+85//9PgQcXF9QMyfP5+VK1fS0tLC+++/z3e+8x2uv/561q1b59X4QRO6UicoPtrMBdMzvLKv8c4PhuJR0kMXEV762nxqmtpJjYvAWbw2oqxWKwsXLuxKfK4ecmRkZFebwsJC/vOf/3i8DqknVq1aRWNjY9ewyxtvvIGIkJ+f77b9woUL6ezspK6ujgsvvHDA/UdHR3PZZZexfft2fvKTn3gt7u40oSvVTWOb48rO4V5U5JISG0FMRBjFtaMjoQNYLEJafOTADb1o69at3HfffVx77bVMnDiR2tpafvrTnzJr1izmz59PdnY29957Lz/60Y9oaGjgkUceISvLuzOQREdHc8kll3D//fdTXl7O/fffz5VXXsmMGTPctp86dSp33HEHS5Ys4YEHHiA/P5/W1lZ27NjB3r17efrpp3n77bd59tlnueKKK8jNzaW0tJT/+7//6zGm702a0JXqpthLJYsuIkJ2cjTFR707PBBsMjMzycjI4Mc//jFlZWUkJSVx3nnn8dOf/pTIyEhef/117rrrLq6++mqys7P5/ve/z5o1a7x6Gf6SJUuIj4/n1ltvpbGxkcsvv5wnn3yy3+csW7aMKVOm8NRTT/Hwww+TkJDAjBkzuPXWWwHHCVgR4aGHHqKqqoq0tDQuvfRSHn/8ca/F3Z0YP5WL5+fnmw0bNvjl2Er1ZeX2Cu54oYB/fPNsTs1O9Mo+b31uPaXHWlj57XO9sr/edu3axfTp032y71CRl5fH1Vdfzc9//nN/hwL0/56KSIExxu04kFa5KNVN0VHHVZ3eGnIByEmJoaS2RWddVD6nCV2pboqONpMUYyUx2uq1fWYnR9PYZuNYc4fX9qmUOzqGrlQ3hTXeK1l0yXFVutQ2kxw7cmWAynOHDx/2dwheoT10pbop9mINukt2crRz33piVPmWJnSlnGyddkpqW4Y9y2Jv3XvovqLj88FjOO+lJnSlnMrrWrHZjdd76AlRjjH5Eh8ldKvV6vWrJpX/tLS0YLUO7RyOJnSlnIq6atCHvo5oX3JSfFeLnp6eTmlpKc3NzdpTH8WMMTQ3N1NaWnrCvDKe0pOiSjkV1jgSureHXABykmPYU9ng9f0CJCQkAFBWVkZHh1bSjGZWq5WMjIyu93SwNKEr5VR4tImIMAsZCVFe33dOSgwf7K7CbjdYhrHwdF9cq+Go0KZDLko5FR9tJjs5mjAfJNzs5GjabXaqG9u8vm+lXDShK+VUUtvSVZHibTnOxTJ8dWJUKdCErlSXktoWspw1496Wk6K16Mr3NKErBTS32zja1N51EZC3uZazGy3zoqvRSRO6UkBpraPn7I11RN2JsoaRFh85quZFV6OPJnSlcAy3AGQl+aaHDui86MrnNKErxfGTlTk+GnJx7DuGkmPaQ1e+owldKaDkWAsR4RZS43y39FpOSjRlx1qxddp9dgwV2gZM6CKSIyKrRWSniOwQkXvctBER+Y2I7BeRrSIy1zfhKuUbJbUtZCdF++SiH5ec5Bg67YbyulafHUOFNk966DbgXmPMDGA+cJeI9F41dREw2flzO9D/QnxKBRhfliy6jMSsiyq0DZjQjTHlxpiNztsNwC6g93Lbi4E/GYd1QJKIjPV6tEr5SGlts89KFl26Li7SE6PKRwY1hi4iecAc4NNeD2UBxd3ul3Bi0kdEbheRDSKyobq6enCRKuUjrR2dHGls91nJoktmYhQijvF6pXzB44QuInHAa8C3jTH1QzmYMWa5MSbfGJOflpY2lF0o5XUjUbIIEBFuITMhqqvmXSlv8yihi4gVRzJ/0RjzupsmpUBOt/vZzm1KBTxXyaKvh1zA8aGh87koX/GkykWAZ4Bdxphf9tHsTeCrzmqX+UCdMabci3Eq5TMlPr5KtLus5GhKdchF+Ygn86GfBSwFtonIZue2h4BcAGPMH4B3gIuB/UAzcLP3Q1XKN0qPtWANE9LjfVeD7pKdHM3bW8vptBufTNOrQtuACd0Y8xHQ7/8841j36i5vBaXUSCqpbWGcj2vQXbKSYrDZDZX1rYzz8Zi9Cj16pagKeSUjULLo4qp1L9ETo8oHNKGrkOe4StT34+dwvJKmVOd0UT6gCV2FtNaOTqob2nx+laiL65uAli4qX9CErkJa2TFXhcvIJPQoaxipcRE65KJ8QhO6CmkjWbLokpWkpYvKNzShq5DWdZXoCPXQwfHhoUMuyhc0oauQVnqsmXCLkDECNegurouLHNW+SnmPJnQV0kpqWxibFEV42Mj9KWQlRdNms1Pd2DZix1ShQRO6CmkltS0+n5SrN610Ub6iCV2FtJLa5q55ykeKa7xeT4wqb9OErkJWm62Tyvq2Ea1wgeMXF2npovI2TegqZJXWjmwNukt8lJWEqHAdclFepwldhSxXD9m11udIyk6O0SEX5XWa0FXIKvFTDx0c4+i60IXyNk3oKmQV1zZjDRMyEqJG/NhZSdGU1motuvIuTegqZLnmQffHQhPZydE0tXdS19Ix4sdWwUsTugpZIzkPem/ZOi+68gFN6CpkjeQ86L1lOY+rCV15kyZ0FZJc86DnpPi3h66VLsqbNKGrkOSPaXO7S4qxEhsRppUuyqs0oauQ5Eqk/hpDFxGyk2MoPqo9dOU9mtBVSCr240VFLjkpWouuvEsTugpJJbXNRIRZSIsbuXnQe8tOjqFEa9GVFw2Y0EXkWRGpEpHtfTz+eRGpE5HNzp+HvR+mUt5VUttCVnI0Fj/UoLtkJ0fT2GbjWLPWoivv8KSH/hxw0QBt/m2Mme38eWz4YSnlWyVH/VeD7uIa7inWYRflJQMmdGPMWuDoCMSi1IgpqW3xW4WLi2sedj0xqrzFW2PoC0Rki4i8KyIn99VIRG4XkQ0isqG6utpLh1ZqcJrbbdQ0tfu9h56d4rpaVHvoyju8kdA3AuONMbOA3wJ/66uhMWa5MSbfGJOflpbmhUMrNXilAVDhApAQZSUx2qpDLsprhp3QjTH1xphG5+13AKuIpA47MqV8pNjPNejd5aRE65CL8pphJ3QRyRQRcd4+3bnPmuHuVylf8ec86L3lJMdoD115TfhADUTkJeDzQKqIlACPAFYAY8wfgKuBb4iIDWgBlhgtrFUBrKS2hchw/9agu2QnR/Ph7iqMMTj7RUoN2YAJ3Rhz3QCP/w74ndciUsrHip0li4GQQHNSYmiz2aluaCPdDwttqOCiV4qqkBMIJYsuXaWLOuyivEATugopxhgO1zSR6+cKFxfX9L16YlR5gyZ0FVKONXfQ0Gpj/JjASOjHF7rQHroaPk3oKqQcrmkCIG9MrJ8jcYiOCCM1LlJ76MorNKGrkFJY4+gJ56UGRg8dnLXo2kNXXqAJXYWUwppmRPy3UpE72VqLrrxEE7oKKYU1TYxLjCbKGubvULrkJEdTfqwVW6fd36GoUU4Tugoph2uaAuaEqEtOSgw2u6GivtXfoahRThO6CimFNc2MD5AToi46ja7yFk3oKmTUt3ZQ09QecD1015wyOo6uhksTugoZRa4KlwBL6OOSohFxrKKk1HBoQlchw1WDHmhDLhHhFsYlRlOoCV0NkyZ0FTJcNeiBNuQCjrr4w0ea/B2GGuU0oauQcfhIE+nxkcREDDjJ6IjLGxPLoSNN6MzTajg0oauQUVjTHDCX/Pc2ITWW+lYbtc0d/g5FjWKa0FXICMQadBfXB80hHXZRw6AJXYWE5nYbVQ1t5KUGZg/dFZeOo6vh0ISuQoLrhGigzIPeW25KDBY5Xomj1FBoQlchoTDAps3tLSLcQlZytA65qGHRhK5CQlcPPUDH0MHxYaM9dDUcmtBVSDhc00xKbASJ0VZ/h9KnCamxHD7SrKWLasg0oauQUBjAFS4ueWNiaWyzcaSx3d+hqFFKE7oKCYePNDE+QE+IukxwVbrosIsaogETuog8KyJVIrK9j8dFRH4jIvtFZKuIzPV+mEoNXWObjbK6Vialx/k7lH5p6aIaLk966M8BF/Xz+CJgsvPnduDJ4YellPfsq2wAYEpGvJ8j6V92cjRhFtEeuhqyARO6MWYtcLSfJouBPxmHdUCSiIz1VoBqdKhtaqe53ebvMNzaV9kIBH5Ct4ZZyEmO5vARnXVRDY03ZinKAoq73S9xbiv3wr5VgDtQ3cjvVx/gb5tLiY0IY+mC8dx05gTS4iP9HVqXPZUNRFkt5AT4GDo4hl20Fl0N1YhOOycit+MYliE3N3ckD628zG43PPj6Nv5aUExkuIUbzsilqqGN3685wNP/PsT9X5zKbedM9HeYAOytbGBSehxhFvF3KAPKGxPLZ4eOYoxBJPDjVYHFGwm9FMjpdj/bue0ExpjlwHKA/Px8LbYdxZ781wFe2VDMTWfm8c2Fk0iNc/TID1Q38vjbu/jxO7uYkhHPuVPS/BypI6GfNSnV32F4ZEJqLM3tnVQ3tJGeEOXvcNQo442yxTeBrzqrXeYDdcYYHW4JYhuLavnlqr1cOnMsj1w2oyuZA5yUFsfvrp/L5PQ4/uuVzVT5eSX7uuYOKuvbAn783MVV6aLDLmooPClbfAn4BJgqIiUicquI3CEidzibvAMcBPYDTwF3+ixa5Xd1LR3c/dImxiZG8fiXTnU7LBAdEcay6+fS1G7jnpc302n335exvVWOCpepoyShTxijtehq6AYccjHGXDfA4wa4y2sRqYD2/Te2UV7Xyoo7FpAQ1fdl9JMz4nls8Sk88OpWfr96P986f/IIRnncXmfJ4uSMwK5BdxmXFIU1TDiklS5qCPRKUeWxDYeP8tbWcu5eOJm5uckDtr9mXjaXnDqWZWv2c6SxbQQiPNHeigZiI8LISor2y/EHKzzMwvgxsRyobvR3KGoU0oSuPPa71ftJiY3ga+dO8Ki9iPCdL0yh3WbnqX8f9HF07u2tbGRyRvyoqhiZmhHPnooGf4ehRiFN6Moj20vrWLOnmlvPnjCoRZZPSovjslnj+PMnhRxtGvlJp/ZVNYya8XOXqZnxFB1tDtgLtVTg0oSuPPL7NfuJjwxn6YLxg37uN8+bREtHJ898NLK99JrGNo40to+a8XOXqZmOD6C9lTrsogZHE7oa0P6qBt7dXsFXzxzf74nQvkzOiOfiU8fy/MeFHGseuV66KyG6EuRo4fpGsaei3s+RqNFGE7oa0O/XHCAqPIxbzvJs7Nydby2cRGObjWc/OuTFyPq3d5RMytVbbkoM0dYwdus4uhokTeiqX+V1Lfx9cxnXnZ7LmLihz88yLTOBL56cwfOfFNLa0enFCPu2t7KBxGgr6QE0r4wnLBZhSkZc1weSUp7ShK76tWJDCZ12w01n5g17XzeemUddSwf/2FI2/MA8sLeygSkZcaOqwsVlaqZWuqjB04Su+mS3G15ZX8zZk1K9srjygoljOCktlhc+LfJCdP0zxrC3snHUDbe4TMmI50hju9/q99XopAld9emj/UcoPdbCktNzBm7sARHhhvnj2VJ8jO2ldV7ZZ18Ka5qpa+ng5HGJPj2Or0zLTADQXroaFE3oqk8vry8iOcbKhTMyvLbPL83NJtoaxgvrCr22T3cKCmsByM8b+IrWQOSqzNGErgZDE7py60hjG6t2VnLV3Gwiw8O8tt/EaCuXzxrH3zeXUd/a4bX99lZQVEt8VDiT0kZXDbpLalwEKbERmtDVoGhCV269vrGEjk7jteGW7m6YP56Wjk5eLyjx+r5dNhbWMjc3GcsoWNTCHRFhakY8u7XSRQ2CJnR1AmMML68vJn98Msa0Uj8AABiISURBVJPSvX9S8dTsRGblJPHCp0U4Juv0rvrWDvZUNjBv/OgcbnGZmhnPvsoG7H6cfliNLprQ1QkKCms5WN3Etad5v3fu8pUzctlf1cgG51i3N20uOoYxjPqEPi0znub2TkpqW/wdiholNKGrE7y2sZRoaxgXnzrWZ8e4dOZY4iPDeckHJYwFhbVYBGblJHl93yNpivPE6G6dAkB5SBO66qG1o5O3tpax6JRMYiN9t4Z4TEQ4i+eM461t5V6f32VjUS3TMhOI82H8I2FKhla6qMHRhK56eH9XJQ2tNr40N9vnx7ru9FzabXbe2OR2TfEh6bQbNhUdG/XDLQBxkeHkpETriVHlMU3oqofXN5YyNjGKBSeN8fmxTh6XyKzsRF76zHsnR/dWNtDYZguKhA5walYiW4qP+TsMNUpoQlddqhva+Nfeaq6Yk0XYCJX7XXd6LnsrG9lY5J2k5bqgKFgS+tzcZEpqW6isb/V3KGoU0ISuury5pYxOu+FLc7JG7JiXzRpHbEQYL33mnZOjGwtrSYuPJDt5dKwhOpC5zg+mjT6oBlLBRxO66vL6xhJmZicyeQQntIqNDGfxnCze2lpGrReWqCsoqmVebvKonGHRnZPHJRARbmFjkSZ0NTBN6ApwlMbtKKsf0d65y1cXjKe1w85fhtlLL6ltprCmedTO3+JOZHgYM7MSu4aSlOqPJnQFwKsbSgi3CJfNGjfix56WmcDZk1L50yeHabfZh7yft7aWA/DFkzO9FFlgmDs+me2l9bTZRmZhEDV6eZTQReQiEdkjIvtF5HtuHr9JRKpFZLPz5zbvh6p8paPTzt82l3L+9PRhrUo0HLeePYHK+jbe2VY+5H38Y0sZs3OSyEkZ/tztgWRubjLtnXa2l+oFRqp/AyZ0EQkDlgGLgBnAdSIyw03TV4wxs50/T3s5TuVDa/ZUc6SxnWvm+e5S/4F8bkoaE9NieeajQ0MqYTxQ3ciOsnq/fMPwtbnjHVe86olRNRBPeuinA/uNMQeNMe3Ay8Bi34alRtKKDcWkxkXwualpfovBYhFuOWsC20rrWH948InrrS3liDimFAg26fFR5KRE64lRNSBPEnoWUNztfolzW29XichWEXlVRNx29UTkdhHZICIbqqurhxCu8raaxjY+3F3FlXOysIb595TKVXOzSYqx8sxHBwf1PGMMb24p5YwJKWQkRPkoOv+al5tMQWGtT2anVMHDW3/B/wDyjDEzgVXA8+4aGWOWG2PyjTH5aWn+6w2q4/62uQyb3XC1H4dbXKIjwvjKGbm8t7OSXeWejxfvKm/gQHVTUA63uMwdn0xVQxulx3TmRdU3TxJ6KdD9rz3bua2LMabGGONazfZpYJ53wlO+9mqBo/bcteSZv33tnIkkRVv54Zs7PO6N/mNrGWEWYdEpwTfc4jI311GKqeWLqj+eJPT1wGQRmSAiEcAS4M3uDUSk+1/S5cAu74WofGV7aR27yuu5Zp7vJ+LyVFJMBPd9cSqfHjrK2x5UvBhj+MeWMs6elEpKbMQIROgf0zLjiYkIY5OXpkhQwWnAhG6MsQHfBP6JI1H/1RizQ0QeE5HLnc3uFpEdIrIFuBu4yVcBK+958dNCoqwWLp818hcT9WfJabnMGJvA42/vornd1m/b1zaWUlLbwpV+uCBqJIWHWZidk8S6gzX+DkUFMI/G0I0x7xhjphhjTjLG/Ni57WFjzJvO2w8aY042xswyxpxnjNnty6DV8NU1d/DGplKumJ1FYozV3+H0EGYRHl18MmV1rfxhzYE+25Uea+HRN3dw+oQULg/i8XOX86ams7uigeKjzf4ORQUovVI0RK0oKKa1w87SBeP9HYpbp+WlsHj2OP6w9iAf7TtywuN2u+G7r26l0xh+cc2sUbsY9GBcMCMDcMxZr5Q7mtBDkN1u+PO6QvLHJ3PyuER/h9On/75kBhPGxHLTHz/j1YKSHo+98GkhH+0/wn9fMiPorgzty4TUWCalx2lCV30a3Wt0qSH5175qCmuaufcLU/0dSr/S4iNZ8Y0FfOOFAu5bsYW9lQ3ER4azo6ye1Xuq+NyUNK473f/lliPpwhkZPLX2IHUtHSRGB9ZQmfI/7aGHoD9/UkhqXCQXjYJJrBKirPzxptP50twslq89yC9W7WV3RT2LTsnkf66ZGTTT5HrqgukZ2OyGNXuq/B2KCkDaQw8xRTXNrN5TxbcWTiYi3Auf52WbBv+ccXMG1Twi3MIvrpnFXedNIj0+kvio0O2ZzslJIjUukvd2VrJ4dnBX9qjB0x56iFn+7wOEW4TrT8/1dyiDIiKclBYX0skcHHPeXDA9nX/tqR7WVMMqOGlCDyGlx1p4ZX0xX87PITMxOOc8CQUXTM+gsc2mNenqBDrkEkJ+v3o/AHeeN8nPkQSwERhCGq6zJ6cSbQ1j1c5Kzp2icyKp47SHHiJKapv564Zirj0th6yk4FhAOVRFWcM4d0oq726v0FWMVA+a0EPE79ccQBDu/Lz2zoPB9WeM50hjG29tGfoKTyr4aEIPASW1zaxw9s7Hae88KJw7OZXJ6XE8PcQVnlRw0oQeAn7y7m5EhG98/iR/h6K8RES47ZwJ7Cqv5xM9Oaqc9KRokHt/ZyVvby3n3gunBE7vfKROPA7lOKPI4tlZ/GzlHp759yHOPCnV3+GoAKA99CDW0NrBD/6+nakZ8Xz9c9o7DzZR1jBumD+eD3ZXcaC60d/hqACgCT2I/fyfe6iob+WJq071zlWhKuDcMH88EeEWnv3okL9DUQFA/8qDVEHhUf60rpAbF+Qxx7l8mQo+afGRfGlOFis2lLCnosHf4Sg/0zH0IFR2rIVvvLCRrKRo7vviIGZUDPIx52B13xensmpnJfeu2Mwbd56FNUz7aaFKE3qQaWjt4Jbn1tPS3skLt51BXKS+xT7n56tLU+Mi+fGVp3DHCxv5/eoD3HPBZK/tW40u+lEeRDo67dz1l03sr2rkyRvmMSUj3t8hqRFy0SljuWL2OH774T62l9b5OxzlJ5rQg0RLeyfffmUza/dW8/iVp3L2ZC1jCzWPXn4KKbER/Ncrmzna1O7vcJQfaEIPAqXHWrj6Dx/zzrZyHrp4Gl8+bfSs4mO3Q21zB3qx4/Alxlj532tnU3S0mav/8DGlx1r8HZIaYZrQR7m1e6u5/LcfUVTTzDM35nP7ub6vNx9MEu6vrd0OD/1tGzf98TMefGMbtk5zQtvjzw/MjB9oH0hnTkrlz7eeQXVDG1c/+TH7KrXyJZSIv/5Q8vPzzYYNG/xy7GCwvbSOn/1zD2v3VjMxNZblX81nUnrc8Hbqwck9VxLeVV7P9LEJPH7FqQDUNrcjAskxEbhWhXPX1mJxbK9rdSTpm59bT6fdEGYRJqfHsa+qkUnpcfzsqpkI4vb5/cVW19pBUrQVY47f7r5Knd0OtS3tCD1j7b2f2uZ2DAaLCInRVupbbW731d/r692+Bx9PubuzrJ4b//gZrR2d3L1wMksXjCfKGubTY6qRISIFxph8d495VAIhIhcBvwbCgKeNMU/0ejwS+BMwD6gBrjXGHB5O0OpETW02PthdxZuby3h/VyWJ0VYeungaX12Qd+If6xBLEAdKRnWtHewqr6fTbthVXs/hmkb+b+1BdpTVA3BSagy/+PJsLGKh8GgTu8rq6DSwq7yeutYOEqOsPPjGVnaX1zMtM55pmQnsrqhnUnoceyoaMMCeiga++9pWHlw0vcex6lo7SI5xv2JR9+Q6LTMBMOyuaGD62AT+3+JTaGizkRAVzkNvbOuKdWpGPD+7eiZhFumxnwff2MaOsuMnFqOtYbTZOpkxLrHHh0rv34Xr9bn7wDvhd+rjypgZ4xJ4/Rtn8t9/286P39nFH/9ziG9fMIWLZ47VyqcgNuA7KyJhwDLgQqAEWC8ibxpjdnZrditQa4yZJCJLgJ8C1/oi4FBh67RT1dDG7op6tu7aw7bKdv5T0k6rzZAea+Gu/DhunxdHYmQ9VG/1rFfo1Ffbvnqc3R83xjAtM57dFQ1EhYdxz8ub6f4d78CRZpYsX8dJ6XGONtZwWjpsTBubQEJUOIdrGrsS6o7yBv732tmkxEaQEBXOdU9/Sku7Y37vvRUNiNCV8KeNTSCpn1XuuyfX3RX1YEzXB8m9KzZzqLqJyRlx7Ks8fon8nsoG/uuVTdx89kROGRtPY7sdY5zP76alwxHTzrI6Dh5pxCIQFxlOY5uNaRlx7K5sZFpmPMYYjrW0H0/yZfUcbWrl56v2sausjskZ8fz0qp4fIL6UkxLD87eczsf7j/DEyt088NpW/vtv25l/0hgWTk3j5KxEJqXFkRwbMSLxKN8bcMhFRBYAPzTGfNF5/0EAY8xPurX5p7PNJyISDlQAaaafnQ91yOXwkSb+tbd60M/zhe4vz/TYfvy+MQZjwG4Mdue/tk6DzW6n3WanzWanzdZJQ6uN+lYb9S0dVNW3UtnQRqfdsRcBTkoOZ0F2BJdOiSZ/bMQJvUpPhyZ6t3X1XpOirRxr6eCmP37WNQTy3M2nkxhlpbalHXunnZ+s3M2+ykaiwi202ez0t6KlBboetwDTMuPoNMLeyoYevysBpmTG8+Ciadz6/Iau1zwt05H8Hnpjm6M3PzaBn1w5s+t19f5QMga+9/o2R/LPdJRr7q5owBomtHYcj3RqRhx7Kt3Pe2IRmJweR3iYpetDp7vIcKHN1vO/dExEOL+7fja/WLWvK07sdnaUN3Q9p6PT8d67jv8/V8/ud+jIHbsdahJnkBoXgQz0ie2GMYbPDh3l/V2VrNpZyeGa5q7HUmIjSI+PJCU2gpTYCOIiw4myhhFlDSMiTAgPsxAeJoSJYBFBhK4YBHp0CnpHNpRYQ8EpWQnMG58ypOcOd8glCyjudr8EOKOvNsYYm4jUAWOAI70CuR24HSA3d2iLFO8oq+eRN3cM6bmBRAQiwy1EhocREW4hPjKc+GgrCVHhTEwbw7jEaDITonjxs0L2VjQwJimBR6+bj8VN766moY1Xy8qx2cewo0y4N3EGafGRbo/bve32Utjzjxa2ldQxb3wyf7ntDKJyOykorGVebjIJE/JZ8tQ6Pjtc63x2uuOnY+DXFxMRRnP78dV0tnatw3DikmnbyuHgqk4ic+ZSUHSMmdmJ/OyOBdQ0dfBaeQU2eyrby+DehOmkJ0Rhtxuue2qdI87xybz0tfkAHLK2sN1eS6Q1mb/cNp/91Y0s+vW/e3yAPH3TediBC365tkd8ABjYWgGzsxN44s5TsWPnoTd2sKuigelj4ykoPTHJ0wYHwibzWlmt43daJvzxpnm88qyzs9IBU9Jj2FvlSKDbK+CBhOlkJHi+puvx1/tB1+t19/+gPyLCGRPHcMbEMTx08XRKj7Wwr6qR/ZWNHDzSSHVDO7XN7ewoq6e53UZLeyetHXbaO3Uhal+443MnDTmh92dEB9OMMcuB5eDooQ9lHxfMSGfjDy70alzD0f3PqmdPRboetAhYnL2bMIsQbhGP/iCrG9r44T920GmgoLCWmqZ2t4k6NS6CeeOTuxJcalzfX6G7t52ZnciWkjo67YYNhbUcbe7gpa/Np6apndS4CI40tlNQWNvnvgBOGRfPrvIG5uYmUd/awZ7KJoAeyVLo+Q3Gna2ldXz83YVYLNLVC02Ni2BubhKfHa6l08A3X9rEy874CgprsdlN1+8FYGPRMTqN49/alg6mZsZzWl5ytw8kuPvlzfz2+rm0dvS9dNvmknouXfafrvuzshJO7Ho6xUeGc/qEZObmJvPZ4aN02g2/fG9vjza/+vJsLvndxz1+H4NR09TOhsLarvepr/8HnhIRspNjyE6O4byp6f22NcbQaTfY7KbrW6brWxQGTLd3tvf38QAp/AlIUVbfFBh6ktBLge6FzdnObe7alDiHXBJxnBz1usjwMCLDQ+NsvaeJWkR6JOL+vuZ2b5scHc6c//c+Da02YiLCSImxYrFIV7JwHb97QuzaD3BanqO3eLS5A7vdzhk/+dDtMZ+/OZ+7XtpMQ6uN2IgwpmXEUlDcs7c7b3wy6QmRXbHb7YaapnZ+e90czvzpajrtho3OZNbX76X3NhHh5dsXsLuinkt/+xF2Z7J3xJ7ChsJaThkbzwOLpnFabhLXPvUpm4pPvMpyS6+e+Qu3nc6ElBjqW21MzYzHYrHw2+vncOZPPqDTwLbSeubkJLG1pI5545OYPi6R0/NSKChyxNY7Gbtea1/vXUqMlZiIsB7v00gREcLDhBD5kxv1PEno64HJIjIBR+JeAlzfq82bwI3AJ8DVwIf9jZ8rzwwmUXdPxANxta1uaKO5zQZAc5uNo80dPfbhSojVjW0Yu52vv1DAttJ65uUmsewr80iLdyTgtPhIqupb+zxeanxU13FaOzr5/Q2ngRw/B2Fx7qN7MncNqczNTWJebjIbi3omane/F3fbLBZh+tgETstL6Ur2afGRbtv+4YZ8FjzxQdd4tztzc5M466RURISsbtvT4yPJ73aMv9x2BkebO7p9sLh/H90NH/X+9na0uaPf90kplwETunNM/JvAP3GULT5rjNkhIo8BG4wxbwLPAH8Wkf3AURxJX3nBYBL1YKXGRfRIQu6+AVgs0jXe+8adZ/f54ZIWH8npeY4e8uycRFo7OtlZ3shpeclMy4zvcZzuPXF3ug+pbCw6xn++txCLyAmJuvfvpa/flbsPABFOaJueEOnouR8+yszsJH5//RwsFuFbL2/uGqJ69Y4FbmN3d4zu++8rNnfDR73befI+KQV6YVHIG+jr/lD3ZQw99juY4xhjWLL8eK/15dvnj1i1hLs4vfk76s3T1+rLGNTo0l+ViyZ0FZBCKYGF0mtVwzfsK0WVGmm+HGoKNKH0WpVv6eRcSikVJDShK6VUkNCErpRSQUITulJKBQlN6EopFSQ0oSulVJDwWx26iFQDhX45uHup9JodMsAEenwQ+DEGenygMXpDoMcHw4txvDHmxGlL8WNCDzQisqGvYv1AEOjxQeDHGOjxgcboDYEeH/guRh1yUUqpIKEJXSmlgoQm9OOW+zuAAQR6fBD4MQZ6fKAxekOgxwc+ilHH0JVSKkhoD10ppYKEJnSllAoSmtDdEJF7RcSISKq/Y+lORP5HRHaLyFYReUNEkvwdE4CIXCQie0Rkv4h8z9/x9CYiOSKyWkR2isgOEbnH3zG5IyJhIrJJRN7ydyzuiEiSiLzq/D+4S0QW+Dum3kTkv5zv8XYReUlEogIgpmdFpEpEtnfbliIiq0Rkn/PfZG8cSxN6LyKSA3wBKPJ3LG6sAk4xxswE9gIP+jkeRCQMWAYsAmYA14nIDP9GdQIbcK8xZgYwH7grAGMEuAfY5e8g+vFrYKUxZhowiwCLVUSygLuBfGPMKTiWzAyE5TCfAy7qte17wAfGmMnAB877w6YJ/US/Ah4AAu5ssTHmPWOMzXl3HZDtz3icTgf2G2MOGmPagZeBxX6OqQdjTLkxZqPzdgOORJTV/7NGlohkA5cAT/s7FndEJBE4F8f6wRhj2o0xx/wblVvhQLSIhAMxQJmf48EYsxbHWsvdLQaed95+HrjCG8fShN6NiCwGSo0xW/wdiwduAd71dxA4EmNxt/slBFiy7E5E8oA5wKf+jeQE/4ujI2H3dyB9mABUA390Dgs9LSKx/g6qO2NMKfBzHN+uy4E6Y8x7/o2qTxnGmHLn7Qogwxs7DbmELiLvO8fXev8sBh4CHg7g+Fxtvo9jGOFF/0U6+ohIHPAa8G1jTL2/43ERkUuBKmNMgb9j6Uc4MBd40hgzB2jCS8ME3uIch16M48NnHBArIjf4N6qBGUftuFdGBEJuTVFjzAXutovIqTj+I2xxLtSbDWwUkdONMRX+js9FRG4CLgXON4FxEUEpkNPtfrZzW0ARESuOZP6iMeZ1f8fTy1nA5SJyMRAFJIjIC8aYQEpGJUCJMcb1zeZVAiyhAxcAh4wx1QAi8jpwJvCCX6Nyr1JExhpjykVkLFDljZ2GXA+9L8aYbcaYdGNMnjEmD8d/4LkjmcwHIiIX4fhafrkxptnf8TitByaLyAQRicBxEupNP8fUgzg+oZ8BdhljfunveHozxjxojMl2/r9bAnwYYMkc599BsYhMdW46H9jpx5DcKQLmi0iM8z0/nwA7cdvNm8CNzts3An/3xk5Droc+yv0OiARWOb9FrDPG3OHPgIwxNhH5JvBPHFUFzxpjdvgzJjfOApYC20Rks3PbQ8aYd/wY02j0LeBF5wf3QeBmP8fTgzHmUxF5FdiIY0hyEwEwDYCIvAR8HkgVkRLgEeAJ4K8iciuOacS/7JVjBca3dqWUUsOlQy5KKRUkNKErpVSQ0ISulFJBQhO6UkoFCU3oSikVJDShK6VUkNCErpRSQeL/A7voUv9bYbB+AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "samples = metropolis_hastings(func, Gaussian(mu=np.ones(1), var=np.ones(1)), n=100, downsample=10)\n", "plt.plot(x, func(x), label=r\"$\\tilde{p}(z)$\")\n", "plt.hist(samples, density=True, alpha=0.2)\n", "plt.scatter(samples, np.random.normal(scale=.03, size=(100, 1)), s=5, label=\"samples\")\n", "plt.legend(fontsize=15)\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.3" } }, "nbformat": 4, "nbformat_minor": 1 }